0 引言

化学驱已成为国内外大多油田提高采收率的主要方法，其中聚合物驱技术是应用最广泛的方法^[1]，以大庆油田为例，已有50余个区块成功开展了工业化规模的聚合物驱油。相较于水驱，聚合物驱采收率平均提高了约13%，正在实施的二类油层聚合物驱提效项目，也有望在水驱基础上提高原油采收率10%以上^[2-5]。注聚过程中为了提高井注入能力，往往会提高其注入压力，当注入压力高于地层破裂压力时，将在近井地带形成诱导裂缝^[6]，给油藏进一步挖潜带来困难。因此，明确诱导裂缝特征及参数成为国内外学者研究的热点。

关于诱导裂缝对压力动态的影响研究人员已做过很多研究。Hagoort^[7]首先给出注水诱导裂缝的概念，并提出可采用与井筒储集效应相似的变流量方法求解关井压力测试时裂缝闭合导致的流量变化；Koning^[8]在Hagoort的基础上利用拉氏变换等方法求解出二维瞬态椭圆方程的解析解，给出了在裂缝存在条件下水驱的压力动态；Hoek等^[9-10]研究了水驱在考虑动态裂缝闭合不同边界条件下的井底压力解，在处理裂缝闭合时，将裂缝储集系数分为常数及随压力变化2种情况，并给出无限导流裂缝闭合的解析解及半解析解。国内学者对诱导裂缝成因机理及其延伸规律以及注水诱导裂缝关井闭合的压力动态曲线也进行了相关研究^[11-13]。

在实际水驱开发后的注聚阶段，由于聚合物具有剪切变稀、浓度扩散、渗透率下降等特征^[14]，考虑诱导裂缝闭合条件下的数学模型的数值求解就显得尤为方便，PEBI变尺度混合网格及其符合物理渗流的网格划分在数值计算中凸显优势^[15-16]。为此，本文采用数值试井方法，基于PEBI网格剖分，综合考虑聚合物溶液黏浓关系、剪切效应、渗透率下降以及吸附机理，建立考虑诱导裂缝的聚合物驱垂直裂缝井数学模型，并对注聚井压力动态曲线特征及影响因素进行分析，形成考虑诱导裂缝闭合的注聚井压力降落动态研究方法。该研究结果可进一步拓展注聚井压力动态分析方法，以期为准确认识诱导裂缝规模及参数提供重要的理论支持。

1 物理模型的建立 1.1 混合PEBI网格剖分方法

网格剖分是油藏空间离散的基础，在渗流过程中PEBI网格剖分可分为圆模块、圆角矩形模块、平面模块几种形式^[15]。PEBI网格的优势在于可利用变尺度缩小网格节点数目，减少计算量，且网格剖分更符合物理渗流过程，描述流动更准确。

1.2 物理模型

假设无限大地层中有一口注聚井，一条与井筒共面的垂直诱导裂缝，且在垂直方向上贯穿整个油层，以定流量Q注入，建立考虑诱导裂缝的均质油藏聚合物驱垂直裂缝井物理模型（图 1）。对含垂直裂缝井的无限大地层进行空间离散，在裂缝周围采用复合线性流动的圆角矩形模块进行网格剖分；在远井地带，采用复合椭圆及径向流特征的变尺度网格剖分（图 2）。

假设条件：外边界无限大，地层上下为不渗透边界，诱导裂缝关于井筒对称，裂缝末端封闭；诱导裂缝内压力相同；单相微可压缩流体流动；考虑井筒储存和表皮效应。

2 聚合物渗流特征数学模型 2.1 浓度分布模型

聚合物溶液在多孔介质内流动，会发生浓度扩散与对流作用，综合Fick定律、对流通量变化及物质守恒原理，得到聚合物溶液在地层中浓度扩散方程的散度形式为

$ \operatorname{div}\left(D \cdot \nabla C_{\mathrm{p}}\right)-\frac{Q}{2 {\rm{ \mathsf{ π} }} r h \varphi} \cdot \operatorname{div}\left(C_{\mathrm{p}}\right)=\frac{\partial C_{\mathrm{p}}}{\partial t} $

(1)

式中：Q为流量，m³/s；r为网格到井中心的径向距离，m；h为地层厚度，m；φ为孔隙度，%；D为扩散系数，m²/s；t为时间，s；C_p为聚合物质量浓度，kg/m³。

2.2 黏度模型

聚合物溶液的黏度与其在多孔介质内的浓度及剪切速率有关。聚合物溶液黏度μ_p、剪切速率$ \dot \gamma $的关系可表达为^[17]

$ \mu_{p}=\mu_{\infty}+\frac{\mu_{0}-\mu_{\infty}}{1+\left(\frac{\dot{\gamma}}{\dot{\gamma}_{1 / 2}}\right)^{P_{\alpha}-1}} $

(2)

式中：μ_p为聚合物流动时的有效黏度，Pa·s；μ_∞为剪切率趋于无穷大时聚合物溶液黏度，近似等于水的黏度，Pa·s；μ₀为聚合物零剪切黏度，Pa·s；$ \dot \gamma $为剪切速率，s^-1；$ {\dot \gamma _{1/2}} $为μ_p = (μ₀ + μ_∞) /2μ时所对应的剪切速率，s^-1；P_α为无因次常数。

零剪切黏度与浓度的关系可表示为

$ \mu_{0}=\mu_{\mathrm{w}}\left[1+\left(A_{1} C_{\mathrm{p}}+A_{2} C_{\mathrm{p}}^{2}+A_{3} C_{\mathrm{p}}^{3}\right)\right] $

(3)

式中：A₁，A₂，A₃为与聚合物溶液有关的常数，其单位分别为（kg/m³）^-1，（kg/m³）^-2和（kg/m³）^-3；μ_w为水的黏度，Pa·s。

2.3 渗透率模型

聚合物溶液在多孔介质流动过程中，往往伴随着吸附滞留现象，增大附加阻力从而导致地层渗透率下降，通过引入渗透率下降系数R_k来修正渗透率模型，建立聚合物溶液修正渗透率模型^[18]。

$ R_{\mathrm{k}}=1+\frac{\left(R_{\mathrm{k} \max }-1\right) b_{\mathrm{p}} C_{\mathrm{p}}}{\left(1+b_{\mathrm{p}} C_{\mathrm{p}}\right)} $

(4)

式中：b_p为实验确定常数，10^-4；R_{k max}^[19]为最大渗透率下降系数，1.621。

渗透率模型为

$ K_{\mathrm{p}}=\frac{K}{R_{\mathrm{k}}} $

(5)

式中：K为地层渗透率，mD；K_p考虑聚合物吸附作用后的渗透率，mD。

3 考虑诱导裂缝聚合物驱垂直裂缝井数学模型 3.1 模型的建立

考虑诱导裂缝的均质油藏聚合物驱垂直裂缝井数值试井模型主要包括基本渗流微分方程，井以定产量生产、外边界无限大的边界条件以及初始条件3个部分，由此建立考虑诱导裂缝的聚合物驱垂直井数值试井数学模型。

运动方程为

$ v=-\frac{1 \times 10^{-15} K_{\mathrm{p}}}{\mu_{\mathrm{p}}} \nabla P $

(6)

式中：v为渗流速度，m/s；$ \nabla P $为压力梯度，Pa/m。

液体状态方程为

$ B=\frac{B_{0}}{1+C_{\mathrm{L}}\left(P-P_{a}\right)} $

(7)

式中：B为体积系数；B₀为原始体积系数；C_L为液体压缩系数，Pa^-1；P为网格压力，Pa；P_a为原始大气压，Pa。

岩石状态方程为

$ \varphi=\varphi_{0}\left[1+C_{\mathrm{f}}\left(P-P_{\mathrm{a}}\right)\right] $

(8)

式中：φ₀为原始孔隙度，%；C_f为岩石压缩系数，Pa^-1。

连续性方程为

$ \operatorname{div}(\rho v)+\frac{\partial(\rho \varphi)}{\partial t}=0 $

(9)

式中：ρ为流体密度，kg/m³。

将式（6）、式（7）、式（8）同时代入式（9）中，得到考虑诱导裂缝的均质油藏单相流数值试井模型微分方程

$ \operatorname{div}\left(\frac{1 \times 10^{-15} K_{\mathrm{p}}}{\mu_{\mathrm{p}} B} \nabla P\right)=\frac{\partial}{\partial t}\left(\frac{\varphi}{B}\right) $

(10)

3.2 模型的求解

有限体积法具有处理复杂边界渗流问题的优越性，本文采用有限体积法进行离散求解。PEBI网格单元与其邻接网格关系如图 3所示。

首先对浓度方程在控制体内进行积分

$ \left.\int\limits_{t_{n}}^{t_{n}+1} \int\limits_{V} \operatorname{div}\left(D \cdot \nabla C_{\mathrm{p}}\right)-\frac{Q}{2 {\rm{ \mathsf{ π} }} r h \varphi} \cdot \operatorname{div}\left(C_{\mathrm{p}}\right)\right] \mathrm{d} V \mathrm{d} t=\int\limits_{t_{n}}^{t_{n}+1} \int\limits_{V} \frac{\partial C_{\mathrm{p}}}{\partial t} \mathrm{d} V \mathrm{d} t $

(11)

根据高斯定理，将体积积分[式（11）]转化为面积积分，再根据PEBI网格特点，得到质量浓度方程离散形式为

$ \begin{aligned} \sum\limits_{j} D \cdot &G_{i j}\left[\left(C_{\mathrm{p}}\right)_{j}-\left(C_{\mathrm{p}}\right)_{i}\right]^{n+1}-\frac{Q}{2 {\rm{ \mathsf{ π} }} r h \varphi}\left(C_{\mathrm{p}}\right)_{j}^{n} V_{j}= \\ &\frac{V_{i}}{\Delta t}\left[\left(C_{\mathrm{p}}\right)_{i}^{n+1}-\left(C_{\mathrm{p}}\right)_{i}^{n}\right] \end{aligned} $

(12)

对式（10）在控制体积和时间段上积分，得到

$ \int\limits_{t_{n}}^{t_{n+1}} \int\limits_{V} \operatorname{div}\left(\frac{1 \times 10^{-15} K_{\mathrm{p}}}{\mu_{\mathrm{p}} B} \nabla P\right) \mathrm{d} V \mathrm{d} t \int\limits_{t_{n}}^{t_{n+1}} \int\limits_{V} \frac{\partial}{\partial t}\left(\frac{\varphi}{B}\right) \mathrm{d} V \mathrm{d} t $

(13)

根据高斯定理，将式（13）左边的体积分简化为网格单元界面周围的面积分

$ \int\limits_{t_{n}}^{t_{n+1}} \int\limits_{V} \operatorname{div}\left(\frac{1 \times 10^{-15} K_{\mathrm{p}}}{\mu_{\mathrm{p}} B} \nabla P\right) \mathrm{d} V \mathrm{d} t=\int\limits_{t_{n}}^{t_{n}+1} \int\limits_{V} \operatorname{div}\left(\frac{1 \times 10^{-15} K_{\mathrm{p}}}{\mu_{\mathrm{p}} B} \nabla P\right) \mathrm{d} S \mathrm{d} t=\int\limits_{t_{n}}^{t_{n+1}} \sum\limits_{j} \frac{1 \times 10^{-15} K_{\mathrm{p}}}{\mu_{\mathrm{p}} B} \frac{\left(P_{j}-P_{i}\right)^{n}}{d_{i j}} \omega_{i j} \mathrm{d} t $

(14)

由于PEBI网格局部正交特性，最终离散得到

$ \sum\limits_{j} T_{i j}\left(P_{j}-P_{i}\right)^{n+1}=\frac{V_{i}}{\Delta t}\left(P_{i}^{n+1}-P_{i}^{n}\right)\left[\frac{\varphi_{0} C_{\mathrm{f}}+\varphi_{0} C_{\mathrm{L}}}{B_{0}}\right] $

(15)

式中：P_i为网格i的压力，Pa；P_j为与网格i邻接的网格j的压力，Pa；T_ij为传导系数，为PEBI网格任意2个相邻网格中心点之间的流动系数λ_ij与其几何因子G_ij的乘积，m³/（Pa·s）；λ_ij为任意2个相邻网格中心点之间的流动系数，m³/（Pa·s）；G_ij为网格i和j之间的几何因子，m；ω_ij为相邻网格i，j之间流体流动的截面面积，m²；$ \nabla {C_p} $为聚合物质量浓度梯度，kg/m⁴；S为表皮系数；Δt为关井时间，s；V_i为网格i的体积，m³；V_j为网格j的体积，m³。

对考虑诱导裂缝的聚合物驱垂直裂缝井物理模型进行简化，通过计算得到地层总流量为

$ \begin{aligned} Q=&Q_{i}^{n+1}+Q_{\mathrm{c}}^{n+1}=\sum\limits_{i}\left[W I_{i} \cdot \lambda\left(P_{\mathrm{wf}}^{n+1}-P_{i}^{n+1}\right)\right]+\\ & \frac{C}{\Delta t}\left(P_{\mathrm{wf}}^{n+1}-P_{\mathrm{wf}}^{n}\right) \end{aligned} $

(16)

式中：Q_i井筒与相邻网格i间的流量，m³/s；Q_c为井筒储集效应的流量，m³/s；P_wf为井底压力，Pa；C为井筒储集系数，m³/Pa；λ为网格i与井筒间的传导率，m²/（Pa·s）；WI_i为生产指数，m，WI_i = θh/[ln(r_i/r_w)+ S]，其中θ为网格i相对井的张开弧度，rad；r_i为网格i的径向距离，m；r_w为井筒半径，m；S为表皮系数。

结合式（12）、式（15）、式（16）及边界条件得到考虑诱导裂缝的聚合物驱垂直裂缝井数值离散方程，求解得到井底压力动态变化，并对计算结果进行无量纲化，无量纲化的计算公式为

$ P_{\mathrm{wD}}=\frac{2 \times 10^{-15} {\rm{ \mathsf{ π} }} K h}{Q \mu_{\mathrm{p}} B_{0}}\left(P_{\mathrm{wf}}-p_{i}\right) $

(17)

$ t_{\mathrm{D}}=\frac{1 \times 10^{-15} K}{\varphi_{\mathrm{p}} \mu_{\mathrm{p}} C_{\mathrm{t}} r_{\mathrm{w}}^{2}} t $

(18)

$ C_{\mathrm{D}}=\frac{C}{2 \pi h \varphi_{\mathrm{p}} C_{\mathrm{t}} r_{\mathrm{w}}^{2}} $

(19)

$ F_{\mathrm{CD}}=\frac{K_{\mathrm{f}} w_{\mathrm{f}}}{K L_{\mathrm{f}}} $

(20)

式中：P_wD为无量纲井底压力；t_D为无量纲时间；C_D为无量纲井筒储集系数；F_CD为无量纲导流系数；C_t为流体和裂缝的总压缩系数，Pa^-1；p_i为原始地层压力，Pa；r_w为井筒半径，m；φ_p为考虑聚合物吸附作用后的孔隙度，%；K_f为裂缝渗透率，mD。

3.3 诱导裂缝闭合模型

注聚井关井测压，井底流压P_wD逐渐降低，当其降低至低于裂缝开始时的闭合压力P_wD⁰时，诱导裂缝开始闭合。Dam等^[20]的实验结果表明，关井测压时诱导裂缝首先从尖端开始闭合，之后逐渐向井筒靠拢。为此本文考虑裂缝横向闭合（图 4），裂缝半长变化满足式（21）^[21]。

$ \begin{array}{ll} L_{\mathrm{f}}\left(P_{\mathrm{wD}}\right)=L_{\mathrm{f} 0} & P_{\mathrm{wD}}>P_{\mathrm{wD}}^{0} \\ L_{\mathrm{f}}\left(P_{\mathrm{wD}}\right)=L_{\mathrm{f} 0} \times\left(1-\frac{P_{\mathrm{wD}}^{0}-P_{\mathrm{wD}}}{ { delpat }}\right) & P_{\mathrm{wD}} \leqslant P_{\mathrm{wD}}^{0} \end{array} $

(21)

式中：L_f为裂缝半长，m；L_f0为初始裂缝半长，m；P_wD为井底压力；P_wD⁰为关井后裂缝开始闭合时的压力；delpat为裂缝闭合速度因子。

诱导裂缝闭合会导致其渗透率发生变化，张阳等^[22]的研究发现，即使裂缝完全闭合，其渗透率仍比地层渗透率高950.4倍左右。说明即使压裂过程中没有支撑剂，裂缝的存在也可显著提高地层渗流条件。

由此建立动态渗透率模型

$ \begin{aligned} &K_{\mathrm{f}}(i=0)=K_{\mathrm{f} 0} \\ &K_{\mathrm{f}}(i=1)=K_{\mathrm{f} 1} \end{aligned} $

(22)

式中：K_f0为裂缝初始渗透率，mD；K_f1为裂缝闭合渗透率，mD；i为网格节点号。基于上述方法，其计算流程如图 5所示。

4 压力动态曲线特征及参数敏感性

通过计算得到无量纲井底压力及其压力导数，绘制典型压力和压力导数动态曲线，从而研究不同参数对动态曲线的影响（表 1）。

4.1 压力曲线特征分析

考虑诱导裂缝的聚合物驱垂直裂缝井典型试井曲线（图 6）存在5个流动段：纯井筒储存阶段（Ⅰ段），该阶段曲线主要受井筒储集系数影响，在双对数坐标中，无因次压力和压力导数曲线重合，呈现一条斜率为1.00的直线；双线性流阶段（Ⅱ段），为诱导裂缝内及裂缝周围流体的线性流，该阶段压力和压力导数曲线互相平行，压力导数曲线呈斜率为0.25的直线段；线性流阶段（Ⅲ段），即诱导裂缝内流体向地层的流动，压力导数曲线呈斜率为0.50的线性流阶段；椭圆流阶段（Ⅳ段），此阶段压力导数曲线呈斜率为0.36的直线；径向流阶段（Ⅴ段），压力波及范围逐渐扩大，水驱压力导数曲线呈水平直线段而聚合物驱则为一条上翘的曲线，这是由于聚合物在多孔介质内的剪切变稀作用，渗流阻力变化从而曲线上翘。本模型计算的聚合物驱垂直裂缝的压力动态特征曲线与徐有杰等^[23]得出的特征曲线流动阶段及各阶段特征一致，说明了该方法正确可靠。

4.2 参数敏感性

聚合物浓度的变化对试井压力产生影响，聚合物浓度越大其黏度越大，渗流阻力也越大（图 7）。从图中可看出，随着浓度增加，双线性流及非牛顿径向流阶段曲线均靠上；当聚合物注入质量浓度为0时，径向流阶段压力导数曲线呈水平直线段，与水驱现象一致。裂缝长度也会影响试井压力，随着裂缝半长的增加，双线性流持续时间变长（图 8）。

4.3 诱导裂缝闭合压力动态特征

为验证方法的正确性，本次研究将聚合物设置为注入水，计算注水井关井后的裂缝闭合压力，并将压力动态曲线与Hoek^[21]的参数解析解进行对比研究（图 9）。从图上可以看出，基于PEBI网格的考虑裂缝闭合的数值计算结果与Hoek^[21]的参数解析解吻合良好。诱导裂缝闭合时，压力导数曲线出现一个突起的尖峰。

将诱导裂缝闭合过程看作裂缝长度的变化，通过引入裂缝闭合速度因子（delpat）来控制裂缝闭合的快慢。诱导裂缝闭合使无因次压力及压力导数曲线迅速上升，压力导数曲线上升趋势更为显著，形成一个“尖峰”（图 10）。delpat越小，裂缝闭合越快，压力导数“尖峰”也就越明显。相反，delpat越大，裂缝闭合越慢，压力导数“尖峰”越不显著，曲线越靠右。delpat值的变化主要对后期拟径向流出现时间的早晚产生影响，聚合物驱与水驱诱导裂缝闭合的区别主要在裂缝闭合后，由于聚合物扩散及剪切稀化等流变特性使得其径向流阶段压力导数曲线上翘。

5 结论

（1）考虑诱导裂缝的注聚井压力动态曲线划分为5个流动段：纯井筒储存阶段，双线性流阶段，线性流阶段，椭圆流阶段，非牛顿径向流阶段。

（2）聚合物浓度主要影响试井压力变化的后期。曲线上表现为随着浓度增加，曲线上翘明显；裂缝半长越小，线性流阶段持续时间越短。

（3）注聚井关井测压时，裂缝闭合速度因子（delpat）越大，闭合越慢，压力导数曲线峰值越低，曲线越靠右，晚期拟径向流出现时间越晚。