压裂井组非线性渗流模型求解

黄迎松

（中国石化胜利油田分公司勘探开发研究院，山东东营 257015）

低渗透油藏中流体的流动不再符合线性达西渗流规律，只有当压力梯度大于启动压力梯度后流体才会流动。随着驱替压力梯度增大，流体在低渗透地层中的渗流过程经历非线性和拟线性2个渗流阶段。目前，在研究低渗透油藏渗流时，大都将非线性和拟线性渗流阶段简化为超过拟启动压力梯度后的拟线性渗流段，忽略了非线性渗流阶段对流体流动的影响[1-3]。尹芝林等人[4]基于动态渗透率的概念，采用统一的运动形式描述非线性和拟线性渗流阶段的渗流规律，认为与采用拟线性渗流规律相比，采用非线性渗流规律计算出的压力变化更为平缓。综合考虑前人的研究成果，笔者提出分段描述非线性、拟线性渗流规律，同时考虑压裂裂缝中高速流动特征的求解方法，以精细描述低渗透油藏的流动规律，精确刻画地层中的压力分布，提高低渗透压裂注采井组生产指标的计算精度。

开发低渗透、特低渗透油藏时一般采用水力压裂，学者们对此进行了很多研究[5-7]，有部分学者将裂缝和地层看成统一系统，但这种处理会造成模拟不准确；压裂裂缝中流体的流速快，可能会出现高速非达西渗流，需要根据雷诺数判断流动形态[8]。此外，裂缝的导流能力沿缝长及随裂缝中压力的变化而变化[9-10]。基于以上特征，笔者考虑非线性和拟线性渗流段的渗流特征，建立了低渗透三维油水两相达西渗流和高速非达西渗流耦合的数学模型，采用有限差分法进行了求解，并分析了计算结果。

1 压裂注采井组耦合数学模型的建立

1.1 基本假设

建立模型前，进行以下假设：1）油藏内流体的流动为等温流动，油藏外边界封闭；2）地层岩石和流体微可压缩；3）三维地层中有油水两相参与渗流，分油藏和压裂裂缝2个区域分别建立渗流方程；4）油藏区域考虑非线性和拟线性2种渗流的特征，认为油相和水相的启动压力梯度为常数；5）裂缝系统考虑达西和高速非达西流动，且考虑裂缝导流能力的时变特征；6）考虑重力和毛细管力的影响。

1.2 油藏系统的渗流方程

图1为低渗透地层中流体流速与压力梯度的关系曲线（图中，GA为流体的启动压力梯度，GB为线性流动阶段的拟启动压力梯度，GC为开始呈现线性流动时的压力梯度，单位均为10-1MPa/cm）。

当驱替压力梯度不小于GA且不大于GC时，流体流速-压力梯度曲线呈现下凹的非线性段，可用二次函数近似描述非线性段：

图1 低渗透地层中流体流速与压力梯度的关系Fig. 1 Relationship between fluid flow velocity and pressure gradient in low permeability formation

当驱替压力梯度大于GC时，流体流速-压力梯度关系曲线为一条直线，可描述为：

式中：p为压力，10-1MPa；v为流体流速，cm/s；K为渗透率，D；μ为流体黏度，mPa·s；a，b和c为二次函数方程的系数。

考虑三维油水两相流动，且油藏区域与裂缝之间存在交互流动项，根据式（1）和式（2），可得到油藏系统的渗流方程。

式中：D为油藏深度，cm；Kr为相对渗透率；ρ为流体密度，g/cm3；φ 为油藏孔隙度；S 为饱和度；τlmf为油藏和压裂裂缝系统交互流动项，g/（cm3·s）；qlm为单位时间单位体积的产量项，g/（cm3·s）；下标 l=o，w；o，w分别代表油相和水相；下标m代表油藏；下标f代表裂缝；为压力梯度，10-1MPa/cm；t为时间，s；g 为重力加速度，m/s2。

1.3 裂缝系统的渗流方程

因为裂缝宽度较小，可以忽略流体在宽度y方向的流动，建立坐标系 Ox′z′，x′轴沿裂缝延伸方向，z′轴与油藏坐标系的z轴相同，当其与油藏系统坐标系不产生混淆的情况下也可将裂缝坐标系记为（x，z）。裂缝中流体的流动形态用卡迪雷夫雷诺数来判断。

裂缝中流体流动的卡迪雷夫雷诺数计算式为：

裂缝中流体的运动方程为：

β为非达西因子，由介质参数孔隙度和渗透率决定，可表示为：

裂缝系统的连续性方程为：

τlmf为油藏和裂缝系统的交互流动项，可表示为：

σ是基质块形状因子，由基质岩块形状的维数及其特征长度决定：

式中：d为裂缝面的维数；L为基质部分的特征长度，m。

根据双重介质关于形状因子的计算方法，利用达西公式推导油藏和裂缝系统的交互流动项的计算公式：

式中：Re为流体雷诺数；Kf为裂缝渗透率，D；φf为裂缝孔隙度；wf为裂缝的宽度，m；Df为裂缝穿过油藏网格的长度，m；Dx为油藏网格x方向的步长，m；Dy为油藏网格y方向的步长，m。

M. Y. Soliman[9]研究发现，裂缝的导流能力随着缝长增长呈线性降低或指数降低。水力裂缝随裂缝中压力的变化张开或闭合，该过程中裂缝的导流能力也会发生改变。因此，裂缝的导流能力是缝长和压力的函数，函数形式与裂缝性质有关。

式中：i为生产过程中的时刻序号，i0为初始生产时刻。

1.4 辅助方程

考虑油水两相流动，必须满足以下辅助方程。

油藏系统：

1.5 定解条件

油藏系统的初始条件为：

裂缝系统的初始条件为：

油藏外边界为封闭条件：

井底边界采用定压条件：

式中：Q为油井地面产量，cm3/s；B为流体体积系数。

2 数学模型的求解

采用有限差分法求解压裂注采井组耦合数学模型，采用顺序解法求解压力和饱和度。为提高计算结果的精度，计算时要减小时间步长，以达到工程计算精度的要求。非线性方程组先进行线性化处理，得到线性差分方程组，避免使用迭代法求解非线性差分方程组，以减少计算量和提高计算速度，并给出了油藏系统和裂缝系统的具体差分格式。

2.1 油藏系统的差分格式

式中：Cp和Cl为岩石和流体的压缩系数，10-1MPa-1。

以第一项为例说明差分处理方法：

渗流方程后2项的差分格式为：

y和z方向的差分处理方式同上，不再赘述。

2.2 裂缝系统的差分格式

裂缝系统在x方向上采用非等距网格，根据裂缝的长度和方位及压裂井的坐标，可求得裂缝系统通过的油藏网格和所穿过的长度，每个油藏网格对应裂缝系统的一个网格，油藏和裂缝之间的流动交互项只存在于相应的网格，通过交互项建立耦合方程组。

从而，裂缝系统渗流方程x方向的差分格式为：

同理可以写出z方向的差分格式，结合窜流、产量项以及右端项的差分可以得到裂缝高速非达西渗流的差分格式。油藏系统和裂缝系统相交的网格上才会有窜流交互项，窜流项的差分按照窜流公式离散即可。

由于主要考虑的是油水两相流动，在求解饱和度时可以只求解油相或者水相的饱和度，然后根据辅助方程求出另外一相的饱和度。水相的饱和度可以采用下式求解：

式中：ΔAwn为采用显示处理方法得到的n+1时刻压力项前的系数；GWWT为式（27）中的所有n时刻项的组合；Vp为网格块的孔隙体积，Vp= ΔxΔyΔzφ。

2.3 求解框图

根据上述差分模型，可以将之转化为计算机模型进行求解分析，求解框图见图2。

3 算例分析

3.1 算例计算结果

图2 求解过程框图Fig. 2 Solution process block diagram

选择五点法井网中的1个压裂注采井组，只压裂中心注水井，不压裂采油井，假设裂缝为双翼对称裂缝，裂缝延伸方向与地层最大水平主应力方向一致。井组的基本参数为：储层厚度8.50 m，孔隙度0.21，地层深度1 000.00 m，原始地层压力10.0 MPa，储层平均渗透率4.0 mD，原油黏度3. 5 mPa·s，地层水黏度0.45 mPa·s，井距200.00 m，裂缝半长 80.00 m。根据上文建立的压裂注采井组耦合数学模型编制求解程序，计算该注采井组的动态。该注采井组的网格系统平面如图3所示。

图4所示为该井组注水开发90 d后的压力和饱和度计算结果。从图4可以看出：压力和饱和度围绕压裂裂缝形成等值线，沿着裂缝两翼呈现对称分布；裂缝附近的等压线密集，距离裂缝越远，等压线越稀疏；靠近油井井底时，等压线又逐渐变密。同样，注入水也是沿着压裂裂缝逐渐向外扩散。计算结果显示的压力和饱和度的分布与地层实际情况相符，说明建立的模型是正确的。

3.2 压裂裂缝对注采井组动态的影响

图3 五点法井网网格系统平面示意Fig. 3 Schematic diagram of the five-spot well pattern system

图4 压裂注采井组压力和饱和度的计算结果Fig. 4 Calculation results of the pressure and saturation of fractured injection-production well group

为了研究压裂裂缝对注采井组生产动态的影响，首先计算了考虑裂缝与不考虑裂缝时的注采动态，结果见图5。从图5可以看出：注水井未压裂时，油井的产油量迅速降低，维持在很低的水平；注水井压裂后，产油量先降低后升高，与注水井未压裂时相比产油量明显增高，表明压裂注水井具有明显的增产效果，其主要原因是压裂裂缝提高了注水井的吸水能力，油藏能量得到补充，保持地层压力不降低，生产压差稳定。

图5 注水井压裂对注采井组产量的影响Fig. 5 Effect of fracturing in water injection well on the production of injection-production well group

3.3 压裂裂缝导流能力对注采井组动态的影响

由文献[10]可知，压裂裂缝导流能力是地层闭合压力的函数。目前，压裂井进行数值模拟时多认为裂缝导流能力随时间变化[11-16]，但实际上压裂裂缝导流能力的变化主要还是由闭合压力变化引起的，因此，笔者将导流能力设计成裂缝中压力的函数。在此基础上对比文献算法和本文算法的计算结果，结果见图6。

图6 裂缝导流能力对压裂注采井组产量的影响Fig. 6 Effect of fracture flow conductivity on the production of fractured injection-production well group

由图6可以看出，2种算法计算的2条产量-生产时间曲线具有相同的趋势，但本文算法计算出的产量要略低于文献算法。文献算法中，注水井压裂裂缝的导流能力随生产时间增长而逐渐降低；而本文算法中，注水井压裂裂缝导流能力的变化取决于压裂裂缝中压力的变化：裂缝中的压力大于闭合压力时，裂缝导流能力随裂缝中压力降低而降低；裂缝中压力小于等于闭合压力时，此时裂缝的导流能力趋近于零。注水井压裂裂缝的导流能力决定了注采井组补充地层能量的能力，因此也影响了整个注采井组的产量。这也是本文计算结果与文献计算结果存在一定的差异的原因。

4 结 论

1）根据低渗透油藏非线性渗流特征和压裂裂缝导流能力的变化规律，建立了低渗透油藏压裂井组油藏和裂缝耦合的非线性数学模型。该模型综合考虑了低渗透油藏非线性和拟线性渗流规律，以及压裂裂缝出现达西和非达西渗流的情况，模拟结果与地层实际情况相符，说明该模型正确有效。

2）提出了利用Taylor展开将油藏系统非线性渗流方程转化为线性差分方程组的方法，编制了计算机求解程序，计算结果证明该方法有效，能够用来模拟低渗透油藏压裂井注采组的生产动态，与其他非线性模型相比，解法简单。

3）低渗透油藏具有非线性渗流特征，综合考虑非线性段和拟线性段渗流比仅考虑拟线性渗流更加符合地层实际渗流情况，所以不能忽略非线性段渗流的影响。

4）五点法井网压裂注采井组的计算结果表明，注水井压裂可以产生明显的增产效果，压裂裂缝的导流能力随着裂缝中压力的变化而变化，裂缝导流能力降低，造成注水井注水效果变差，导致油井产量降低。