Solution of Nonlinear Seepage Model for Fracture Well Groupin Low Permeability Reservoirs
-
摘要:
精细描述低渗透油藏中流体流速与与压力梯度的非线性关系，是准确计算低渗透油藏压裂井组产量的基础。为此，在描述低渗透油藏非线性渗流特征的基础上，建立了低渗透油藏和压裂裂缝耦合的非线性数学模型，该模型根据渗流特征将渗流过程分为非线性渗流阶段和拟线性渗流阶段进行计算。利用Taylor展开对非线性数学模型进行线性化处理，建立了有限差分方程组，并编制了计算机求解程序。算例分析表明：采用非线性数学模型计算出的地层中压力和饱和度的分布符合地层实际情况；五点法井网压裂井组注水井的裂缝导流能力会随着裂缝闭合而降低，注水效果变差，导致油井产量降低。研究结果表明，低渗透油藏和压裂裂缝耦合的非线性数学模型可以较准确地描述低渗透油藏中流体流速与压力梯度的非线性关系，为准确计算低渗透油藏压裂井组产量奠定基础，为低渗透油藏注水开发提供指导。
Abstract:Having a closely detailed description of the nonlinear relationship between flow velocity and pressure gradient in low permeability reservoir is necessary for accurately developing the frac design, and calculating the production of a group (or unit) of wells that have been hydraulically fractured. Therefore, based on the description of the nonlinear seepage characteristics of low permeability reservoir, a nonlinear mathematical model of coupling low permeability reservoir and hydraulic fractures was established, which divided the seepage process into the nonlinear seepage stage and quasi-linearity stage according to the seepage characteristics. The Taylor expansion was used to linearize the nonlinear mathematical model, and established the finite difference equations, and then formed the computer solving model. The results of example analysis showed that the distributions of formation pressure and saturation calculated by the nonlinear mathematical model were in line with the actual situations of the stratum; the fracture flow conductivity of injection well in the fractured five-spot well pattern decreased with the formation closure, which led to poor water injection effect and low oil well production. Thus, the fracture design should be modified in accordance with the study’s results. The study results indicated that the nonlinear mathematical model and hydraulic fracture coupling could accurately describe the nonlinear relationship between flow velocity and pressure gradient in low-permeability reservoir. This breakthrough establishes a foundation to calculate the production of fractured well group in low-permeability reservoir accurately, and provides a guidance for water flooding development of low permeability reservoir.
-
低渗透油藏中流体的流动不再符合线性达西渗流规律，只有当压力梯度大于启动压力梯度后流体才会流动。随着驱替压力梯度增大，流体在低渗透地层中的渗流过程经历非线性和拟线性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–1 MPa/cm）。
![]()
图 1 低渗透地层中流体流速与压力梯度的关系Figure 1. Relationship between fluid flow velocity and pressure gradient in low permeability formation当驱替压力梯度不小于GA且不大于GC时，流体流速–压力梯度曲线呈现下凹的非线性段，可用二次函数近似描述非线性段：
v=a(dpdx)2+bdpdx+c (1) 当驱替压力梯度大于GC时，流体流速–压力梯度关系曲线为一条直线，可描述为：
v=Kμ(dpdx−GB) (2) 式中：p为压力，10–1 MPa；v为流体流速，cm/s；K为渗透率，D；μ为流体黏度，mPa∙s；a，b和c为二次函数方程的系数。
考虑三维油水两相流动，且油藏区域与裂缝之间存在交互流动项，根据式（1）和式（2），可得到油藏系统的渗流方程。
当
|∇pl|⩾ 时，有：\begin{split} & \quad\quad\dfrac{\partial }{{\partial x}}\left[ {\dfrac{{{\rho _l}{K_x}{K_{{\rm{r}}l}}}}{{{\mu _l}}}\left( {\dfrac{{\partial {p_l}}}{{\partial x}} - {G_{{\rm{B}}l}}} \right)} \right] + \\ & \quad\quad \dfrac{\partial }{{\partial y}}\left[ {\dfrac{{{\rho _l}{K_y}{K_{{\rm{r}}l}}}}{{{\mu _l}}}\left( {\dfrac{{\partial {p_l}}}{{\partial y}} - {G_{{\rm{B}}l}}} \right)} \right]+\\ & \dfrac{\partial }{{\partial z}}\left[ {\dfrac{{{\rho _l}{K_z}{K_{rl}}}}{{{\mu _l}}}\left( {\dfrac{{\partial {\rho _l}}}{{\partial z}} + {\rho _l}g\dfrac{{\partial D}}{{\partial z}} - {G_{{\rm{B}}l}}} \right)} \right] - \\ & \qquad\quad {\tau _{l{\rm{mf}}}} + {q_{l{\rm{m}}}} = \dfrac{{\partial \left( {{p_l}\phi {S_l}} \right)}}{{\partial t}} \end{split} (3) 当
{G_{{\rm{C}}l}} > \left| {\nabla {p_l}} \right| \geqslant {G_{{\rm{A}}l}} 时，有：\begin{array}{l} \qquad\qquad \dfrac{\partial }{{\partial x}}\left\{ {{\rho _l}\left[ {{a_l}{{\left( {\dfrac{{\partial {p_l}}}{{\partial x}}} \right)}^2} + {b_l}\dfrac{{\partial {p_l}}}{{\partial x}} + {c_l}} \right]} \right\} +\\ \qquad\qquad \dfrac{\partial }{{\partial y}}\left\{ {{\rho _l}\left[ {{a_l}{{\left( {\dfrac{{\partial {p_l}}}{{\partial y}}} \right)}^2} + {b_l}\dfrac{{\partial {p_l}}}{{\partial y}} + {c_l}} \right]} \right\}+\\ \dfrac{\partial }{{\partial z}}\left\{ {{\rho _l}\left[ {{a_l}{{\left( {\dfrac{{\partial {p_l}}}{{\partial z}} + {\rho _l}g\dfrac{{\partial D}}{{\partial z}}} \right)}^2} + {b_l}\left( {\dfrac{{\partial {p_l}}}{{\partial z}} + {\rho _l}g\dfrac{{\partial D}}{{\partial z}}} \right) + {c_l}} \right]} \right\} - \\ \qquad\qquad\;\qquad {\tau _{l{\rm{mf}}}} + {q_{l{\rm{m}}}} = \dfrac{{\partial \left( {{p_l} \phi {S_l}} \right)}}{{\partial t}} \end{array} (4) 式中：D为油藏深度，cm；Kr为相对渗透率；ρ为流体密度，g/cm3；ϕ为油藏孔隙度；S为饱和度；τlmf为油藏和压裂裂缝系统交互流动项，g/（cm3∙s）；qlm为单位时间单位体积的产量项，g/（cm3∙s）；下标l=o，w；o，w分别代表油相和水相；下标m代表油藏；下标f代表裂缝；
\left| {\nabla p} \right| 为压力梯度，10–1MPa/cm；t为时间，s；g为重力加速度，m/s2。1.3 裂缝系统的渗流方程
因为裂缝宽度较小，可以忽略流体在宽度y方向的流动，建立坐标系Ox′z′，x′轴沿裂缝延伸方向，z′轴与油藏坐标系的z轴相同，当其与油藏系统坐标系不产生混淆的情况下也可将裂缝坐标系记为（x，z）。裂缝中流体的流动形态用卡迪雷夫雷诺数来判断。
裂缝中流体流动的卡迪雷夫雷诺数计算式为：
R{e_l} = \frac{{{v_{l{\rm{f}}}}\sqrt {{K_{\rm{f}}}} {\rho _l}}}{{17.50{\mu _l}\phi_{\rm{f}}^{3/2}}} (5) 裂缝中流体的运动方程为：
\left\{ {\begin{array}{*{20}{l}} {{v_{l{\rm{f}}}} = \dfrac{{{K_{\rm{f}}}{K_{rl}}}}{{{\mu _l}}}\nabla p}&{R{e_l} \leqslant 0.3}\\ { - \nabla {p_l} = \dfrac{{{\mu _l}}}{{{K_l}}}{v_{l{\rm{f}}}} + \beta {\rho _l}v_{l{\rm{f}}}^2}&{R{e_l} > 0.3} \end{array}} \right. (6) β为非达西因子，由介质参数孔隙度和渗透率决定，可表示为：
{\beta = \dfrac{c}{{{K^a}{\phi ^b}}}} (7) 裂缝系统的连续性方程为：
- \nabla \cdot \left( {{{\rm{w}}_{\rm{f}}}{\rho _l}{v_{l{\rm{f}}}}} \right) + {w_{\rm{f}}}{\tau _{l{\rm{mf}}}} + {w_{\rm{f}}}{q_{l{\rm{f}}}} = {w_{\rm{f}}}\frac{\partial \left( {{\rho _l}{\phi_{\rm{f}}}{S_l}} \right)}{{\partial t}} (8) τlmf为油藏和裂缝系统的交互流动项，可表示为：
{\tau _{l{\rm{mf}}}} = \sigma \frac{{{K_{\rm{m}}}{K_{rl}}}}{{{\mu _l}}}\left( {{p_{l,{\rm{m}}}} - {p_{l,{\rm{f}}}}} \right) (9) σ是基质块形状因子，由基质岩块形状的维数及其特征长度决定：
\sigma = \frac{{4d\left( {d + 2} \right)}}{{{L^2}}} (10) 式中：d为裂缝面的维数；L为基质部分的特征长度，m。
根据双重介质关于形状因子的计算方法，利用达西公式推导油藏和裂缝系统的交互流动项的计算公式：
{\tau _{l{\rm{mf}}}} = \frac{{4{D_{\rm{f}}}}}{{{D_x}{D_y}}}\left( {\frac{1}{{{D_x}}} + \frac{1}{{{D_y}}}} \right)\frac{{{K_{\rm{m}}}{K_{{\rm{r}}l}}}}{{{\mu _l}}}\left( {{p_{l,{\rm{m}}}} - {p_{l,{\rm{f}}}}} \right) (11) 式中：Re为流体雷诺数；Kf为裂缝渗透率，D；ϕf为裂缝孔隙度；wf为裂缝的宽度，m；Df为裂缝穿过油藏网格的长度，m；Dx为油藏网格x方向的步长，m；Dy为油藏网格y方向的步长，m。
M. Y. Soliman[9]研究发现，裂缝的导流能力随着缝长增长呈线性降低或指数降低。水力裂缝随裂缝中压力的变化张开或闭合，该过程中裂缝的导流能力也会发生改变。因此，裂缝的导流能力是缝长和压力的函数，函数形式与裂缝性质有关。
{K_{\rm{f}}}\left( i \right) = {K_{\rm{f}}}\left( {{i_0}} \right){{\rm{e}}^{ - \alpha \left( {{p_{{\rm{f}},i}} - {p_{{\rm{f}},{\rm{m}}}}} \right) - 0.75{D_{\rm{f}}}}} (12) 式中：i为生产过程中的时刻序号，i0为初始生产时刻。
1.4 辅助方程
考虑油水两相流动，必须满足以下辅助方程。
油藏系统：
{p_{{\rm{m}},{\rm{cow}}}} = {p_{{\rm{m}},{\rm{o}}}} - {p_{{\rm{m}},{\rm{w}}}} (13) {S_{{\rm{m}},{\rm{o}}}} + {S_{{\rm{m}},{\rm{w}}}} = 1 (14) 裂缝系统：
{p_{{\rm{f}},{\rm{cow}}}} = {p_{{\rm{f}},{\rm{o}}}} - {p_{{\rm{f}},{\rm{w}}}} (15) {S_{{\rm{f}},{\rm{o}}}} + {S_{{\rm{f}},{\rm{w}}}} = 1 (16) 1.5 定解条件
油藏系统的初始条件为：
\left\{ {\begin{array}{*{20}{c}} {{p_{{\rm{m}},{\rm{o}}}}\left( {x,y,z,O} \right) = {p_{{\rm{oi}}}}\left( {x,y,z} \right)}\\ {{S_{{\rm{m}},{\rm{w}}}}\left( {x,y,z,O} \right) = {S_{{\rm{wi}}}}\left( {x,y,z} \right)} \end{array}} \right. (17) 裂缝系统的初始条件为：
\left\{ {\begin{array}{*{20}{c}} {{p_{{\rm{f}},{\rm{o}}}}\left( {x,z,O} \right) = {p_{{\rm{oi}}}}\left( {x,z} \right)}\\ {{S_{{\rm{f}},{\rm{w}}}}\left( {x,z,O} \right) = {S_{{\rm{wi}}}}\left( {x,z} \right)} \end{array}} \right. (18) 油藏外边界为封闭条件：
\dfrac{{\partial {p_{{\rm{m}},{\rm{o}}}}}}{{\partial n}}\Bigg| {}_{{\Gamma}} \Bigg. = 0 (19) 井底边界采用定压条件：
{Q_l} = {K_{\rm{f}}}\left( i \right)\dfrac{{{\lambda _l}}({p_{\rm{f}}} - {p_{{\rm{wf}}}})}{{{B_l}{\mu _l}}} (20) 式中：Q为油井地面产量，cm3/s；B为流体体积系数。
2. 数学模型的求解
采用有限差分法求解压裂注采井组耦合数学模型，采用顺序解法求解压力和饱和度。为提高计算结果的精度，计算时要减小时间步长，以达到工程计算精度的要求。非线性方程组先进行线性化处理，得到线性差分方程组，避免使用迭代法求解非线性差分方程组，以减少计算量和提高计算速度，并给出了油藏系统和裂缝系统的具体差分格式。
2.1 油藏系统的差分格式
1）当
\left| {\nabla {p_l}} \right| \geqslant {G_{{\rm{C}}l}} 时，油藏系统处于拟线性段，其渗流方程的差分格式为：\begin{array}{l} \;\; \frac{{{\lambda _{lxi + \frac{1}{2}}}\left( {\dfrac{{p_{li + 1jk}^{n + 1} - p_{lijk}^{n + 1}}}{{\Delta x}} - {G_{{\rm{B}}l}}} \right) - {\lambda _{lxi - \frac{1}{2}}}\left( {\dfrac{{p_{lijk}^{n + 1} - p_{li - 1jk}^{n + 1}}}{{\Delta x}} - {G_{{\rm{B}}l}}} \right)}}{{\Delta x}} +\\ \;\; \frac{{{\lambda _{lyj + \frac{1}{2}}}\left( {\dfrac{{p_{lij + 1k}^{n + 1} - p_{lijk}^{n + 1}}}{{\Delta y}} - {G_{{\rm{B}}l}}} \right) - {\lambda _{lyj - \frac{1}{2}}}\left( {\dfrac{{p_{lijk}^{n + 1} - p_{lij - 1k}^{n + 1}}}{{\Delta y}} - {G_{{\rm{B}}l}}} \right)}}{{\Delta y}}+\\ \frac{{{\lambda _{lzk \!+\! \frac{1}{2}}}\!\!\left( {\dfrac{{p_{lijk + 1}^{n + 1} - p_{lijk}^{n + 1}}}{{\Delta z}}\! -\! \rho g \!-\! {G_{{\rm{B}}l}}} \right)\! -\! {\lambda _{lzk - \frac{1}{2}}}\left( {\dfrac{{p_{lijk}^{n + 1} \!-\! p_{lijk - 1}^{n + 1}}}{{\Delta z}}\! -\! \rho g\! -\! {G_{{\rm{B}}l}}} \right)}}{{\Delta z}} -\\ \qquad\quad {\left( {{\tau _{l{\rm{mf}}}}} \right)_{ijk}} + {\left( {{q_{l{\rm{m}}}}} \right)_{ijk}} = {\left( {{\beta _l}} \right)_{ijk}}\dfrac{{p_{lijk}^{n + 1} - p_{lijk}^n}}{{\Delta t}} +\\ \qquad\quad\quad\quad\quad\quad{\left( {\phi {\rho _l}} \right)_{ijk}}\dfrac{{S_{lijk}^{n + 1} - S_{lijk}^n}}{{\Delta t}} \end{array} (21) \!{\text{其中}}\qquad\qquad\qquad\quad{\lambda _l} = \frac{{{\rho _l}K{K_{{\rm{r}}l}}}}{{{\mu _l}}}\quad\quad\quad\quad (22) {\beta _l} = {\rho _l}\phi {S_l}\left( {{C_{\rm{p}}} + {C_{\rm{l}}}} \right) (23) 式中：Cp和Cl为岩石和流体的压缩系数，10–1 MPa–1。
2）当
{G_{{\rm{C}}l}} > \left| {\nabla {p_l}} \right| \geqslant {G_{{\rm{A}}l}} 时，油藏系统处于非线性段，其有限差分格式关键在于对\dfrac{\partial }{{\partial x}}\left[ {\rho a{{\left( {\dfrac{{\partial p}}{{\partial x}}} \right)}^2}} \right] 项的处理，笔者采用Taylor展开进行差分求解。\frac{\partial }{{\partial x}}\left[ {\rho a{{\left( {\frac{{\partial p}}{{\partial x}}} \right)}^2}} \right] = \frac{{{{\left( {\rho a} \right)}_{i + \frac{1}{2}}}{{\left( {\dfrac{{\partial p}}{{\partial x}}} \right)}^2}_{i + \frac{1}{2}} - {{\left( {\rho a} \right)}_{i - \frac{1}{2}}}{{\left( {\dfrac{{\partial p}}{{\partial x}}} \right)}^2}_{i - \frac{1}{2}}}}{{\Delta x}} (24) 以第一项为例说明差分处理方法：
\begin{split} & \qquad\qquad \left[ {{{\left( {\dfrac{{\partial p}}{{\partial x}}} \right)}^2}} \right]_{i + \frac{1}{2}}^{n + 1} \approx \left[ {{{\left( {\dfrac{{\partial p}}{{\partial x}}} \right)}^2}} \right]_{i + \frac{1}{2}}^n +\\ & \quad\quad 2\left( {\dfrac{{\partial p}}{{\partial x}}} \right)_{i + \frac{1}{2}}^n\left[ {\left( {\dfrac{{\partial p}}{{\partial x}}} \right)_{i + \frac{1}{2}}^{n + 1} - \left( {\dfrac{{\partial p}}{{\partial x}}} \right)_{i + \frac{1}{2}}^n} \right] = \\ & \quad\;\;\; 2\dfrac{{p_{i + 1}^n - p_i^n}}{{\Delta x}} \dfrac{{p_{i + 1}^{n + 1} - p_i^{n + 1}}}{{\Delta x}} - {\left( {\dfrac{{p_{i + 1}^n - p_i^n}}{{\Delta x}}} \right)^2} \end{split} (25) 渗流方程后2项的差分格式为：
\begin{array}{l} \quad\qquad\qquad\quad\quad \dfrac{\partial }{{\partial x}}\left[ {\rho b\left( {\dfrac{{\partial p}}{{\partial x}}} \right) + \rho c} \right] =\\ \dfrac{{{{\left( {\rho b} \right)}_{i + \frac{1}{2}}}\dfrac{{{p_{i + 1}} - {p_i}}}{{\Delta x}} - {{\left( {\rho b} \right)}_{i - \frac{1}{2}}}\dfrac{{{p_i} - {p_{i - 1}}}}{{\Delta x}} + {{\left( {\rho c} \right)}_{i + \frac{1}{2}}} - {{\left( {\rho c} \right)}_{i - \frac{1}{2}}}}}{{\Delta x}} \end{array} (26) x方向的差分格式乘以
V = \Delta x\Delta y\Delta z 之后有：\begin{array}{*{20}{l}} \begin{array}{l} \quad \left[ {\dfrac{{2{{\left( {\rho a} \right)}_{i + \frac{1}{2}}}\Delta y\Delta z}}{{\Delta {x^2}}}\left( {p_{i + 1}^n - p_i^n} \right) + \dfrac{{{{\left( {\rho b} \right)}_{i + \frac{1}{2}}}\Delta y\Delta z}}{{\Delta x}}} \right]p_{i + 1}^{n + 1} + \\ \quad \left[ {\dfrac{{2{{\left( {\rho a} \right)}_{i - \frac{1}{2}}}\Delta y\Delta z}}{{\Delta {x^2}}}\left( {p_i^n - p_{i - 1}^n} \right) + \dfrac{{{{\left( {\rho b} \right)}_{i - \frac{1}{2}}}\Delta y\Delta z}}{{\Delta x}}} \right]p_{i - 1}^{n + 1} +\\ \end{array}\\ \begin{array}{l} \left[ {\dfrac{{2{{\left( {\rho a} \right)}_{i + \frac{1}{2}}}\Delta y\Delta z}}{{\Delta {x^2}}}\left(\! {p_{i + 1}^n \!-\! p_i^n} \!\right) \!+\! \dfrac{{2{{\left( {\rho a} \right)}_{i - \frac{1}{2}}}\Delta y\Delta z}}{{\Delta {x^2}}}\left( {p_i^n \!-\! p_{i - 1}^n} \right) \!+ } \right.\\ \quad\quad\quad\quad \left. {\dfrac{{{{\left( {\rho b} \right)}_{i + \frac{1}{2}}}\Delta y\Delta z}}{{\Delta x}} + \dfrac{{{{\left( {\rho b} \right)}_{i - \frac{1}{2}}}\Delta y\Delta z}}{{\Delta x}}} \right]p_i^{n + 1} + \\ \end{array}\\ \begin{array}{l} \left[ {\dfrac{{{{\left( {\rho a} \right)}_{i - \frac{1}{2}}}\Delta y\Delta z}}{{\Delta {x^2}}}{{\left( {p_i^n - p_{i - 1}^n} \right)}^2} \!- \!\dfrac{{{{\left( {\rho a} \right)}_{i + \frac{1}{2}}}\Delta y\Delta z}}{{\Delta {x^2}}}{{\left( {p_{i + 1}^n \!-\! p_i^n} \right)}^2} + } \right.\\ \qquad\qquad\qquad \left[ {{{\left( {\rho c} \right)}_{i + \frac{1}{2}}} - {{\left( {\rho c} \right)}_{i - \frac{1}{2}}}} \right]\Delta y\Delta z \end{array} \end{array} (27) y和z方向的差分处理方式同上，不再赘述。
2.2 裂缝系统的差分格式
裂缝系统在x方向上采用非等距网格，根据裂缝的长度和方位及压裂井的坐标，可求得裂缝系统通过的油藏网格和所穿过的长度，每个油藏网格对应裂缝系统的一个网格，油藏和裂缝之间的流动交互项只存在于相应的网格，通过交互项建立耦合方程组。
1）当
R{e_l} \leqslant 0.3 时，裂缝系统渗流方程的差分格式为：\begin{split} \begin{split} &\qquad \frac{{{\gamma _{lxi \!+\! \frac{1}{2}}}\!\dfrac{{2\left( {p_{{{\rm{f}}li} \!+\! 1,j}^{n + 1} \!-\! p_{li,j}^{n \!+\! 1}} \right)}}{{\Delta {x_i} \!+\! \Delta {x_{i + 1}}}} \!-\! {\gamma _{lxi \!-\! \frac{1}{2}}}\!\dfrac{{2\left( {p_{{\rm{f}}li,j}^{n \!+\! 1} \!-\! p_{li - 1,j}^{n + 1}} \right)}}{{\Delta {x_i} \!+\! \Delta {x_{i - 1}}}}}}{{\Delta {x_i}}} \!+\! \\ & \qquad \frac{{{\gamma _{lzj + \frac{1}{2}}}\dfrac{{p_{{\rm{f}}lij + 1}^{n + 1} - p_{{\rm{f}}lij}^{n + 1}}}{{\Delta z}} - {\gamma _{lzj - \frac{1}{2}}}\dfrac{{p_{{\rm{f}}lij}^{n + 1} - p_{{\rm{f}}lij - 1}^{n + 1}}}{{\Delta z}}}}{{\Delta z}} +\\ \end{split}\\ \begin{split} &\quad {\left( {{w_{\rm{f}}}{\tau _{l{\rm{mf}}}}} \right)_{ij}} + {\left( {{w_{\rm{f}}}{q_{l{\rm{m}}}}} \right)_{ij}} = {\left( {{w_{\rm{f}}}{\beta _l}} \right)_{ij}}\dfrac{{p_{{\rm{f}}lij}^{n + 1} - p_{{\rm{f}}lij}^n}}{{\Delta t}} + \\ & \quad\qquad\qquad {\left( {{w_{\rm{f}}}\phi {\rho _l}} \right)_{ij}}\dfrac{{S_{{\rm{f}}lij}^{n + 1} - S_{{\rm{f}}lij}^n}}{{\Delta t}} \end{split} \end{split} (28) \!{\text{其中}}\qquad\qquad\qquad{\gamma _l} = {w_{\rm{f}}}\frac{{{\rho _l}K{K_{rl}}}}{{{\mu _l}}}\qquad\qquad\quad (29) 2）当
R{e_l} > 0.3 时，对式（6）两边关于x求偏导，可得到：\frac{{\partial \left( {{v_{lx}}} \right)}}{{\partial x}} = \frac{1}{{ - \dfrac{{{\mu _l}}}{{{k_l}}} + 2\beta {\rho _l}{v_{lx}}}}\dfrac{{{\partial ^2}{p_{{\rm{f}}l}}}}{{\partial {x^2}}} (30) 从而，裂缝系统渗流方程x方向的差分格式为：
\begin{array}{l} \qquad\qquad\qquad\quad \dfrac{{\partial \left( {{w_{\rm{f}}}{\rho _l}{v_{lx}}} \right)}}{{\partial x}} =\\ \dfrac{{{w_{\rm{f}}}{\rho _l}}}{{ - \dfrac{{{\mu _l}}}{{{K_l}}} + 2\beta {\rho _l}v_{lxi,j}^n}}\dfrac{{\dfrac{{2\left( {p_{{\rm{f}}li + 1,j}^{n + 1} - p_{{\rm{f}}li,j}^{n + 1}} \right)}}{{\Delta {x_i} + \Delta {x_{i + 1}}}} + \dfrac{{2\left( {p_{{\rm{f}}li,j}^{n + 1} - p_{{\rm{f}}li - 1,j}^{n + 1}} \right)}}{{\Delta {x_i} + \Delta {x_{i - 1}}}}}}{{\Delta {x_i}}} + \\ \quad c{w_{\rm{f}}}{v_{lx}}{\rho _l}\dfrac{{2\left( {p_{{\rm{f}}li + 1,j}^{n + 1} - p_{{\rm{f}}li,j}^{n + 1}} \right)}}{{\Delta {x_i} + \Delta {x_{i + 1}}}} + {\rho _l}{v_{lx}}\dfrac{{2\left( {w_{{\rm{f}}i + 1,j}^n - w_{{\rm{f}}i,j}^n} \right)}}{{\Delta {x_i} + \Delta {x_{i + 1}}}}\\ \end{array} (31) 同理可以写出z方向的差分格式，结合窜流、产量项以及右端项的差分可以得到裂缝高速非达西渗流的差分格式。油藏系统和裂缝系统相交的网格上才会有窜流交互项，窜流项的差分按照窜流公式离散即可。
{\left( {{\tau _{l{\rm{mf}}}}} \right)_{ijk}} = {\alpha _{ijk}}\left( {{p_{l{\rm{m}},ijk}} - {p_{{\rm{f}}l,ij}}} \right) (32) \!{\text{其中}}\quad\quad\qquad{\alpha _{ijk}} = \sigma \left(\dfrac{{{\rho _l}{{{K}}_{\rm{m}}}{K_{rl}}}}{{{\mu _l}}}\right)_{ijk}\qquad\quad\quad (33) 由于主要考虑的是油水两相流动，在求解饱和度时可以只求解油相或者水相的饱和度，然后根据辅助方程求出另外一相的饱和度。水相的饱和度可以采用下式求解：
\begin{split} & \quad\quad\quad {\left( {\Delta A_{\rm{w}}^n\Delta p_{\rm{w}}^{n + 1} + GWWT + {Q_{\rm{w}}}} \right)_{ijk}} =\\ & \quad\quad\quad \dfrac{1}{{\Delta t}}\left[\left( {\dfrac{{{V_p}{S_{\rm{w}}}}}{{{B_{\rm{w}}}}}} \right)_{ijk}^{n + 1} - \left( {\dfrac{{{V_p}{S_{\rm{w}}}}}{{{B_{\rm{w}}}}}} \right)_{ijk}^{n1}\right] \end{split} (34) 式中：∆Awn为采用显示处理方法得到的n+1时刻压力项前的系数；GWWT为式（27）中的所有n时刻项的组合；Vp为网格块的孔隙体积，Vp= ∆x∆y∆zϕ。
2.3 求解框图
根据上述差分模型，可以将之转化为计算机模型进行求解分析，求解框图见图2。
3. 算例分析
3.1 算例计算结果
选择五点法井网中的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可以看出：压力和饱和度围绕压裂裂缝形成等值线，沿着裂缝两翼呈现对称分布；裂缝附近的等压线密集，距离裂缝越远，等压线越稀疏；靠近油井井底时，等压线又逐渐变密。同样，注入水也是沿着压裂裂缝逐渐向外扩散。计算结果显示的压力和饱和度的分布与地层实际情况相符，说明建立的模型是正确的。
![]()
图 4 压裂注采井组压力和饱和度的计算结果Figure 4. Calculation results of the pressure and saturation of fractured injection-production well group3.2 压裂裂缝对注采井组动态的影响
为了研究压裂裂缝对注采井组生产动态的影响，首先计算了考虑裂缝与不考虑裂缝时的注采动态，结果见图5。从图5可以看出：注水井未压裂时，油井的产油量迅速降低，维持在很低的水平；注水井压裂后，产油量先降低后升高，与注水井未压裂时相比产油量明显增高，表明压裂注水井具有明显的增产效果，其主要原因是压裂裂缝提高了注水井的吸水能力，油藏能量得到补充，保持地层压力不降低，生产压差稳定。
![]()
图 5 注水井压裂对注采井组产量的影响Figure 5. Effect of fracturing in water injection well on the production of injection-production well group3.3 压裂裂缝导流能力对注采井组动态的影响
由文献[10]可知，压裂裂缝导流能力是地层闭合压力的函数。目前，压裂井进行数值模拟时多认为裂缝导流能力随时间变化[11-16]，但实际上压裂裂缝导流能力的变化主要还是由闭合压力变化引起的，因此，笔者将导流能力设计成裂缝中压力的函数。在此基础上对比文献算法和本文算法的计算结果，结果见图6。
![]()
图 6 裂缝导流能力对压裂注采井组产量的影响Figure 6. Effect of fracture flow conductivity on the production of fractured injection-production well group由图6可以看出，2种算法计算的2条产量–生产时间曲线具有相同的趋势，但本文算法计算出的产量要略低于文献算法。文献算法中，注水井压裂裂缝的导流能力随生产时间增长而逐渐降低；而本文算法中，注水井压裂裂缝导流能力的变化取决于压裂裂缝中压力的变化：裂缝中的压力大于闭合压力时，裂缝导流能力随裂缝中压力降低而降低；裂缝中压力小于等于闭合压力时，此时裂缝的导流能力趋近于零。注水井压裂裂缝的导流能力决定了注采井组补充地层能量的能力，因此也影响了整个注采井组的产量。这也是本文计算结果与文献计算结果存在一定的差异的原因。
4. 结 论
1）根据低渗透油藏非线性渗流特征和压裂裂缝导流能力的变化规律，建立了低渗透油藏压裂井组油藏和裂缝耦合的非线性数学模型。该模型综合考虑了低渗透油藏非线性和拟线性渗流规律，以及压裂裂缝出现达西和非达西渗流的情况，模拟结果与地层实际情况相符，说明该模型正确有效。
2）提出了利用Taylor展开将油藏系统非线性渗流方程转化为线性差分方程组的方法，编制了计算机求解程序，计算结果证明该方法有效，能够用来模拟低渗透油藏压裂井注采组的生产动态，与其他非线性模型相比，解法简单。
3）低渗透油藏具有非线性渗流特征，综合考虑非线性段和拟线性段渗流比仅考虑拟线性渗流更加符合地层实际渗流情况，所以不能忽略非线性段渗流的影响。
4）五点法井网压裂注采井组的计算结果表明，注水井压裂可以产生明显的增产效果，压裂裂缝的导流能力随着裂缝中压力的变化而变化，裂缝导流能力降低，造成注水井注水效果变差，导致油井产量降低。
-
![]()
图 1 低渗透地层中流体流速与压力梯度的关系
Figure 1. Relationship between fluid flow velocity and pressure gradient in low permeability formation
![]()
图 4 压裂注采井组压力和饱和度的计算结果
Figure 4. Calculation results of the pressure and saturation of fractured injection-production well group
![]()
图 5 注水井压裂对注采井组产量的影响
Figure 5. Effect of fracturing in water injection well on the production of injection-production well group
-
[1] 宋付权,刘慈群. 含启动压力梯度油藏的两相渗流分析[J]. 石油大学学报(自然科学版), 1999, 23(3): 47–50. doi: 10.3863/j.issn.1674-5086.1999.03.014 SONG Fuquan, LIU Ciqun. Analysis of two-phase fluid flow in low permeability reservoir with the threshold pressure gradient[J]. Journal of the University of Petroleum, China(Edition of Natural Science), 1999, 23(3): 47–50. doi: 10.3863/j.issn.1674-5086.1999.03.014
[2] 程时清, 陈明卓. 油水两相低速非达西渗流数值模拟[J]. 石油勘探与开发, 1998, 25(1): 41–43. doi: 10.3321/j.issn:1000-0747.1998.01.012 CHENG Shiqing, CHEN Mingzhuo. Numerical simulation of oil-water low-velocity non-Darcy flow[J]. Petroleum Exploration and Development, 1998, 25(1): 41–43. doi: 10.3321/j.issn:1000-0747.1998.01.012
[3] 周涌沂,彭仕宓,李允,等. 低速非达西渗流的全隐式模拟模型[J]. 石油勘探与开发, 2002, 29(2): 90–93. doi: 10.3321/j.issn:1000-0747.2002.02.024 ZHOU Yongyi, PENG Shimi, LI Yun, et al. Fully implicit simulation model for low-velocity non-Darcy flow[J]. Petroleum Exploration and Development, 2002, 29(2): 90–93. doi: 10.3321/j.issn:1000-0747.2002.02.024
[4] 尹芝林,孙文静,姚军. 动态渗透率三维油水两相低渗透油藏数值模拟[J]. 石油学报, 2011, 32(1): 117–121. doi: 10.3969/j.issn.1001-8719.2011.01.020 YIN Zhilin, SUN Wenjing, YAO Jun. Numerical simulation of the 3D oil-water phase dynamic permeability for low-permeability reservoirs[J]. Acta Petrolei Sinica, 2011, 32(1): 117–121. doi: 10.3969/j.issn.1001-8719.2011.01.020
[5] 吕广忠,鞠斌山,栾志安. 油藏水力压裂区域分解模拟算法[J]. 石油大学学报(自然科学版), 1998, 22(5): 61–63. LYU Guangzhong, JU Binshan, LUAN Zhian. Domain decomposition simulation method for hydraulic fracturing area of reservoir[J]. Journal of the University of Petroleum, China (Edition of Natural Science), 1998, 22(5): 61–63.
[6] 苏玉亮,王霞,李涛,等. 人工裂缝对低渗透油田开发的影响研究[J]. 钻采工艺, 2006, 29(4): 33–34. doi: 10.3969/j.issn.1000-7393.2006.04.011 SU Yuliang, WANG Xia, LI Tao, et al. Influence of created fracture on low permeability reservoir development[J]. Drilling & Production Technology, 2006, 29(4): 33–34. doi: 10.3969/j.issn.1000-7393.2006.04.011
[7] 何勇明,孙尚如,徐荣伍, 等. 低渗透油藏污染井压裂增产率预测模型及敏感性分析[J]. 中国石油大学学报(自然科学版), 2010, 34(3): 76–79. doi: 10.3969/j.issn.1673-5005.2010.03.016 HE Yongming, SUN Shangru, XU Rongwu, et al. Prediction model for fracturing incremental recovery of damaged well in low-permeability reservoir and sensitivity analysis[J]. Journal of China University of Petroleum(Edition of Natural Science), 2010, 34(3): 76–79. doi: 10.3969/j.issn.1673-5005.2010.03.016
[8] BELHAJ H A, AGHA K R, NOURI A M, et al. Numerical modeling of Forchheimer’s equation to describe darcy and non-Darcy flow in porous media[R]. SPE 80440, 2003.
[9] SOLIMAN M Y. Numerical model estimates fracture production increase[J]. Oil and Gas, 1986, 84(41): 70–74.
[10] 温庆志,张士诚,王秀宇,等. 支撑裂缝长期导流能力数值计算[J]. 石油钻采工艺, 2005, 27(4): 68–70. doi: 10.3969/j.issn.1000-7393.2005.04.020 WEN Qingzhi, ZHANG Shicheng, WANG Xiuyu, et al. Numerical calculation of long - term conductivity of propping fractures[J]. Oil Drilling & Production Technology, 2005, 27(4): 68–70. doi: 10.3969/j.issn.1000-7393.2005.04.020
[11] 任勇,郭建春,赵金洲,等. 压裂井裂缝导流能力研究[J]. 河南石油, 2005, 19(1): 46–48. doi: 10.3969/j.issn.1673-8217.2005.01.017 REN Yong, GUO Jianchun, ZHAO Jinzhou, et al. A study on flow conductivity of fractures in a fractured well[J]. Henan Petroleum, 2005, 19(1): 46–48. doi: 10.3969/j.issn.1673-8217.2005.01.017
[12] 胥元刚,张琪. 变裂缝导流能力下水力压裂整体优化设计方法[J]. 大庆石油地质与开发, 2000, 19(2): 40–43. doi: 10.3969/j.issn.1000-3754.2000.02.014 XU Yuangang, ZHANG Qi. Overall optimizing designation method for hydraulic fracturing under variable fracture diverting capacity[J]. Petroleum Geology & Oilfield Development in Daqing, 2000, 19(2): 40–43. doi: 10.3969/j.issn.1000-3754.2000.02.014
[13] 孔祥言．高等渗流力学[M]．合肥: 中国科学技术大学出版社, 1999: 76–77. KONG Xiangyan. Advanced seepage mechanics[M]. Hefei: Press of University of Science and Technology of China, 1999: 76–77.
[14] 李淑霞, 谷建伟．油藏数值模拟基础[M]．东营: 中国石油大学出版社, 2009: 97–101． LI Shuxia, GU Jianwei. Fundamentals of numerical reservoir simulation[M]. Dongying: China University of Petroleum Press, 2008: 97–101.
[15] 戴嘉尊, 邱建贤．微分方程数值解法[M]．南京: 东南大学出版社, 2002． DAI Jiazun, QIU Jianxian. Numerical solutions for differential equations[M]. Nanjing: Southeast University Press, 2002.
[16] 张建国, 雷光伦．油气层渗流力学[M]．东营: 石油大学出版社, 1998: 46–47． ZHANG Jianguo, LEI Guanglun. Seepage mechanics of oil and gas reservoir[M]. Dongying: Petroleum University Press, 1998: 46–47.
-
期刊类型引用(3)
1. 张文锐，代玉杰，刘贵满. 环形井网稠油渗流压力场数值模拟. 当代化工. 2022(01): 196-201+205 . 
百度学术
2. 李瑞轩，黄云龙，李源. 油藏渗流模型与数值模拟技术研究进展. 石油化工应用. 2021(07): 11-15+29 . 百度学术
3. 房娜，姜光宏，程奇，李广龙，王双龙. 裂缝性油藏不同见水模式下的注水优化. 断块油气田. 2020(05): 633-637 . 百度学术
其他类型引用(2)
下载:
计量
- 文章访问数: 1212
- HTML全文浏览量: 557
- PDF下载量: 52
- 被引次数: 5
