2020, Vol. 32
碳酸盐岩储集空间类型多样、孔隙结构复杂，储层裂缝、孔洞发育，具有很强的非均质性[1]。郑松青等[2]针对碳酸盐岩缝洞型储层呈离散介质的特征，提出缝洞网络模型，并在此基础上建立了油水两相流数学模型，简化了储层描述参数，计算成本低，可以进行大规模的工程应用；贾永禄等[3]基于基质、微裂缝、溶蚀孔洞和大裂缝组成的四重介质，建立了缝洞型四重介质的碳酸盐岩试井模型；熊钰等[4]针对大尺度的碳酸盐岩缝洞型油气藏，考虑缝洞的多尺度特征，建立了井打在大尺度裂缝上的裂缝线性流试井解释模型；程飞[5]针对缝洞型碳酸盐岩储层的非均质性和各向异性强的特征，建立了碳酸盐岩油藏储层类型识别模板，并划分了3种储层类型；姜瑞忠等[6]针对低渗油藏存在的压敏效应和启动压力梯度现象，建立了双重介质低渗油藏的斜井数学模型，并绘制了斜井试井曲线；杜鑫等[7]针对缝洞型油气藏，建立了井筒与溶洞相连通情况下的流动方程并给出了井筒和溶洞之间的压力差表达式；刘建峰等[8]建立了孔洞和裂缝组成渗流空间的典型双重介质模型；王蓓等[9]针对磨溪龙王庙组缝洞型碳酸盐岩储层，建立了多尺度离散裂缝模型，明确了该气藏高、低渗分布区域；姜瑞忠等[10]建立了缝洞型碳酸盐岩油藏水平井分形非线性渗流模型，并结合测井数据将新模型应用于矿场渗流参数的解释，但这些有关碳酸盐岩储层的渗流模型均未考虑二次压力梯度非线性项的影响。
聂仁仕等[11-13]建立了考虑井储和表皮效应的不稳定渗流数学模型，同时，该模型考虑了二次压力梯度非线性项的影响；Nie等[14]针对底水油藏和厚油层油藏，建立了不同边界条件下考虑二次压力梯度非线性项的球面渗流方程；张林等[15]建立了不同边界条件下考虑二次压力梯度非线性项的三区复合渗流方程；张强等[16]建立了考虑启动压力梯度和二次梯度非线性项影响的渗流模型；Wang等[17]针对多层复合储层，建立了考虑二次压力梯度非线性项的多区复合非线性模型；Guo等[18]针对地下地层中气井生产，建立了考虑二次压力梯度非线性项的直井模型；王美楠等[19]建立了考虑启动压力梯度、动边界、二次梯度非线性项的低渗透变形介质油藏渗流模型；Lu等[20]建立了考虑二次压力梯度非线性项的多孔介质水平井渗流模型。目前，针对碳酸盐岩缝洞型三孔介质储层，尚无有关二次压力梯度三孔渗流模型及非线性渗流特征的研究。
流体在地下储层中的渗流具有二次压力梯度的非线性渗流特性[18-20]，缝洞型碳酸盐岩储层也不例外，采用常规的三孔介质线性渗流模型无法精准描述流体在缝洞型碳酸盐岩储层中的渗流特征。笔者首次建立二次压力梯度三孔渗流模型，研究缝洞型碳酸盐岩储层受二次压力梯度控制的非线性渗流模型的缺失问题，以期实现对缝洞型碳酸盐岩储层渗流特征的精准描述。
1 渗流模型 1.1 物理模型平面径向流物理模型假设如下：①裂缝孔洞型三孔介质水平等厚地层中心一口井以定产量q生产；②地层厚度为h，地层渗透率为常数且与压力无关，原始地层压力为pi；③井半径为rw，井底压力为pwf；④井到地层外边界的距离为re，外边界可以是无穷大地层、封闭地层或恒压边界；⑤在弹性驱动方式下，黏度为μ的单相液体渗流符合等温达西定律；⑥考虑天然裂缝与井筒沟通，溶洞系统和基质系统中流体均向天然裂缝系统发生拟稳态窜流，因溶洞渗透率大于基质渗透率，故溶洞系统中流体优先向天然裂缝系统窜流；⑦考虑井筒储集效应和表皮效应的影响。
1.2 数学模型及其求解 1.2.1 数学模型一个完整的数学模型由渗流控制微分方程、初始条件、内外边界条件构成。对缝洞型碳酸盐岩储层，其渗流控制微分方程包括裂缝系统方程、基质系统方程和溶洞系统方程。
(1) 渗流控制微分方程
裂缝系统方程：
| $ \frac{{{\partial ^2}{p_{{\rm{fD}}}}}}{{\partial {\kern 1pt} r_{\rm{D}}^2}} + \frac{1}{{{r_{\rm{D}}}}}\frac{{\partial {\kern 1pt} {p_{{\rm{fD}}}}}}{{\partial {\kern 1pt} {r_{\rm{D}}}}} - \beta {\left( {\frac{{\partial {\kern 1pt} {p_{{\rm{fD}}}}}}{{\partial {\kern 1pt} {r_{\rm{D}}}}}} \right)^2} + {\lambda _{{\rm{mf}}}}{{\rm{e}}^{ - 2S}}({p_{{\rm{mD}}}} - {p_{{\rm{fD}}}}) + {\lambda _{{\rm{vf}}}}{{\rm{e}}^{ - 2S}}({p_{{\rm{vD}}}} - {p_{{\rm{fD}}}}) = {\omega _{\rm{f}}}{e^{ - 2S}}\frac{{\partial {p_{{\rm{fD}}}}}}{{\partial {t_{\rm{D}}}}} $ | (1) |
式中：pfD为无因次裂缝系统压力；pmD为无因次基质系统压力；pvD为无因次溶洞系统压力；rD为无因次径向距离；β为无因次二次压力梯度项系数；λmf为基质系统向裂缝系统的窜流系数；λvf为溶洞系统向裂缝系统的窜流系数；ωf为裂缝系统的弹性储容比；S为井的表皮系数；tD为无因次时间。下标f，m，v，D分别代表裂缝(fracture)、基质(matrix)、溶洞(vug)和无因次(Dimensionless)。
基质系统方程：
| $ {\lambda _{{\rm{mf}}}}{{\rm{e}}^{ - 2S}}({p_{{\rm{fD}}}} - {p_{{\rm{mD}}}}) = {\omega _{\rm{m}}}{{\rm{e}}^{ - 2S}}\frac{{\partial {\kern 1pt} {p_{{\rm{mD}}}}}}{{\partial {\kern 1pt} {t_{\rm{D}}}}} $ | (2) |
式中：ωm为基质系统的弹性储容比。
溶洞系统方程：
| $ {\lambda _{{\rm{vf}}}}{{\rm{e}}^{ - 2S}}({p_{{\rm{fD}}}} - {p_{{\rm{vD}}}}) = {\omega _{\rm{v}}}{{\rm{e}}^{ - 2s}}\frac{{\partial {\kern 1pt} {p_{{\rm{vD}}}}}}{{\partial {\kern 1pt} {t_{\rm{D}}}}} $ | (3) |
式中：ωv为溶洞系统的弹性储容比。
(2) 初始条件
| $ {p_{{\rm{fD}}}}{|_{{t_{\rm{D}}} = 0}} = {p_{{\rm{mD}}}}{|_{{t_{\rm{D}}} = 0}} = {p_{{\rm{vD}}}}{|_{{t_{\rm{D}}} = 0}} = 0 $ | (4) |
(3) 内边界条件
| $ {C_{\rm{D}}}{\left. {\frac{{\partial {p_{{\rm{vD}}}}}}{{\partial {t_{\rm{D}}}}}} \right|_{{r_{\rm{D}}} = 1}} - {\left. {\left( {{r_{\rm{D}}}\frac{{\partial {p_{{\rm{fD}}}}}}{{\partial {r_{\rm{D}}}}}} \right)} \right|_{{r_{\rm{D}}} = 1}} = 1 $ | (5) |
式中：CD为无因次井筒储集系数。
(4) 外边界条件
无限大外边界
| $ {\mathop {{\rm{lim}}}\limits_{{r_{\rm{D}}} \to \infty } {p_{{\rm{fD}}}} = 0} $ | (6) |
恒压边界
| $ {{p_{{\rm{fD}}}}{|_{{r_{\rm{D}}} = {r_{{\rm{eD}}}}}} = 0} $ | (7) |
式中：reD为无因次外边界距离。
封闭边界
| $ {{{\left. {\frac{{\partial {\kern 1pt} {p_{{\rm{fD}}}}}}{{\partial {\kern 1pt} {r_{\rm{D}}}}}} \right|}_{{r_{\rm{D}}} = {r_{{\rm{eD}}}}}} = 0} $ | (8) |
上述数学模型的无因次参数定义如下：无因次压力
| $ {p_{j{\rm{D}}}} = \frac{{{k_{\rm{f}}}h}}{{1.842 \times {{10}^{ - 3}}\mu {\kern 1pt} q{\kern 1pt} B}}({p_{\rm{i}}} - {p_j})\quad (j = {\rm{f}},{\rm{m}},{\rm{v}}) $ | (9) |
式中：pf为裂缝系统压力，MPa；pm为基质系统压力，MPa；pv为溶洞系统压力，MPa；pi为初始油藏压力，MPa；q为流量，m3/d；μ为黏度，mPa·s；B为地层流体的体积系数，无因次；kf为裂缝渗透率，D；h为储层厚度，m。
基于有效井径定义的无因次距离
| $ {r_{\rm{D}}} = \frac{r}{{{r_{\rm{w}}}{{\rm{e}}^{ - S}}}} $ | (10) |
式中：r为径向距离，m；rw为井半径，m。
无因次时间
| $ {t_{\rm{D}}} = \frac{{3.6{k_{\rm{f}}}t}}{{\mu {\kern 1pt} r_{\rm{w}}^2({\varphi _{\rm{f}}}{C_{{\rm{ft}}}} + {\varphi _{\rm{v}}}{C_{{\rm{vt}}}} + {\varphi _{\rm{m}}}{C_{{\rm{mt}}}})}} $ | (11) |
式中：t为时间，h；φf为裂缝孔隙度；φv为溶洞孔隙度；φm为基质孔隙度；Cft为裂缝弹性综合压缩系数，MPa-1；Cvt为溶洞弹性综合压缩系数，MPa-1；Cmt为基质弹性综合压缩系数，MPa-1。
无因次井筒储集系数
| $ {C_{\rm{D}}} = \frac{{{C_{\rm{s}}}}}{{6.283{\kern 1pt} {\kern 1pt} 2{\kern 1pt} {\kern 1pt} h{\kern 1pt} {\kern 1pt} r_{\rm{w}}^2({\varphi _{\rm{f}}}{\kern 1pt} {C_{{\rm{ft}}}} + {\varphi _{\rm{m}}}{\kern 1pt} {C_{{\rm{mt}}}} + {\varphi _{\rm{v}}}{\kern 1pt} {C_{{\rm{vt}}}})}} $ | (12) |
式中：Cs为井的弹性储容系数，m3/MPa。
弹性储容比
| $ {\omega _j} = \frac{{{\varphi _j}{C_{j{\rm{t}}}}}}{{{\varphi _{\rm{f}}}{C_{{\rm{ft}}}} + {\varphi _{\rm{m}}}{C_{{\rm{mt}}}} + {\varphi _{\rm{v}}}{C_{{\rm{vt}}}}}}\quad (j = {\rm{f}},{\rm{m}}{\kern 1pt} {\kern 1pt} {\rm{v}}) $ | (13) |
无因次二次压力梯度项系数
| $ \beta = \frac{{1.842 \times {{10}^{ - 3}}q{\kern 1pt} B{\kern 1pt} \mu {\kern 1pt} {C_{\rm{ \mathsf{ ρ} }}}}}{{{k_{\rm{f}}}{\kern 1pt} h}} $ | (14) |
式中：Cρ为地层流体的压缩系数，MPa-1。
1.2.2 积分变换求解模型引入基于tD的拉氏变换
| $ \begin{array}{*{20}{c}} {L[{p_{j{\rm{D}}}}({r_{\rm{D}}},{t_{\rm{D}}})] = {{\bar p}_{j{\rm{D}}}}({r_{\rm{D}}},\mu ) = \int\limits_0^\infty {{p_{j{\rm{D}}}}} ({r_{\rm{D}}},{t_{\rm{D}}}){e^{ - u{\kern 1pt} {t_{\rm{D}}}}}{\rm{d}}{\kern 1pt} {t_{\rm{D}}}}\\ {(j = {\rm{f}},{\rm{m}}{\kern 1pt} {\kern 1pt} {\rm{v}})} \end{array} $ | (15) |
式中：L[]为Laplace算子；pD为变量pD对应的La‐ place空间变量；u为Laplace变量。
对式(1)-(3)经基于tD的Laplace变换后，得
| $ \frac{{{\partial ^2}{\kern 1pt} {{\bar p}_{{\rm{fD}}}}}}{{\partial {\kern 1pt} r_{\rm{D}}^2}} + \frac{1}{{{r_{\rm{D}}}}}\frac{{\partial {\kern 1pt} {{\bar p}_{{\rm{fD}}}}}}{{\partial {\kern 1pt} {r_{\rm{D}}}}} - \beta {\left( {\frac{{\partial {\kern 1pt} {{\bar p}_{{\rm{fD}}}}}}{{\partial {\kern 1pt} {r_{\rm{D}}}}}} \right)^2} + {\lambda _{{\rm{mf}}}}{{\rm{e}}^{ - 2S}}({\bar p_{{\rm{mD}}}} - {\bar p_{{\rm{fD}}}}) + {\lambda _{{\rm{vf}}}}{{\rm{e}}^{ - 2S}}({\bar p_{{\rm{vD}}}} - {\bar p_{{\rm{fD}}}}) = u{\kern 1pt} {\omega _{\rm{f}}}{{\rm{e}}^{ - 2S}} $ | (16) |
| $ {{\lambda _{{\rm{mf}}}}{{\rm{e}}^{ - 2S}}({{\bar p}_{{\rm{fD}}}} - {{\bar p}_{{\rm{mD}}}}) = u{\kern 1pt} {\omega _{\rm{m}}}{{\rm{e}}^{ - 2S}}{{\bar p}_{{\rm{mD}}}}} $ | (17) |
| $ {{\lambda _{{\rm{vf}}}}{{\rm{e}}^{ - 2S}}({{\bar p}_{{\rm{fD}}}} - {{\bar p}_{{\rm{vD}}}}) = u{\kern 1pt} {\omega _{\rm{v}}}{{\rm{e}}^{ - 2S}}{{\bar p}_{{\rm{vD}}}}} $ | (18) |
将式(17)变形，得
| $ {\bar p_{{\rm{mD}}}} = \frac{{{\lambda _{{\rm{mf}}}}}}{{u{\kern 1pt} {\omega _{\rm{m}}} + {\lambda _{{\rm{mf}}}}}}{\bar p_{{\rm{fD}}}} $ | (19) |
将式(18)变形，得
| $ {\bar p_{{\rm{vD}}}} = \frac{{{\lambda _{{\rm{vf}}}}}}{{u{\kern 1pt} {\omega _{\rm{v}}} + {\lambda _{{\rm{vf}}}}}}{\bar p_{{\rm{fD}}}} $ | (20) |
将式(19)和式(20)带入式(16)，得
| $ {\frac{{{\partial ^2}{{\bar p}_{{\rm{fD}}}}}}{{\partial {\kern 1pt} r_{\rm{D}}^2}} + \frac{1}{{{r_{\rm{D}}}}}\frac{{\partial {\kern 1pt} {{\bar p}_{{\rm{fD}}}}}}{{\partial {\kern 1pt} {r_{\rm{D}}}}} - \beta {{\left( {\frac{{\partial {\kern 1pt} {{\bar p}_{{\rm{fD}}}}}}{{\partial {\kern 1pt} {r_{\rm{D}}}}}} \right)}^2} - f(u){{\bar p}_{{\rm{fD}}}} = 0} $ | (21) |
| $ {f(u) = \left( {\frac{{{\omega _{\rm{m}}}{\kern 1pt} {\lambda _{{\rm{mf}}}}}}{{u{\kern 1pt} {\omega _{\rm{m}}} + {\lambda _{{\rm{mf}}}}}} + \frac{{{\omega _v}{\lambda _{vf}}}}{{u{\kern 1pt} {\omega _v} + {\lambda _{vf}}}} + {\omega _{\rm{f}}}} \right)u{\kern 1pt} {\kern 1pt} {{\rm{e}}^{ - 2S}}} $ | (22) |
当ωf = 1时，非线性模型则简化为常规单一介质模型，式(22)可以简化为
| $ {f(u) = u{\kern 1pt} {\kern 1pt} {{\rm{e}}^{ - 2S}}} $ | (23) |
换句话讲，在相同的油井生产和地层边界条件下，Laplace空间下的三孔介质模型和均质模型的差异仅仅表现在f(u)的函数表达式不同。因此，如果首先在Laplace空间中得到均质模型的解，就可以直接通过改变f(u)来得到Laplace空间下三孔介质模型的解。
因为式(16)是关于pfD的非线性微分方程，该方程包括pfD对rD的一阶、二阶偏微分以及Laplace空间中一阶导数的二次幂，将该方程线性化十分困难。
先令ωf = 1，则考虑单一介质情形，渗流控制方程中就不存在窜流项，则式(1)变为
| $ \frac{{{\partial ^2}{{\bar p}_{{\rm{fD}}}}}}{{\partial {\kern 1pt} r_{\rm{D}}^2}} + \frac{1}{{{r_{\rm{D}}}}}\frac{{\partial {\kern 1pt} {p_{{\rm{fD}}}}}}{{\partial {\kern 1pt} {r_{\rm{D}}}}} - \beta {\left( {\frac{{\partial {\kern 1pt} {p_{{\rm{fD}}}}}}{{\partial {\kern 1pt} {r_{\rm{D}}}}}} \right)^2} = {{\rm{e}}^{ - 2S}}\frac{{\partial {\kern 1pt} {p_{{\rm{fD}}}}}}{{\partial {\kern 1pt} {t_{\rm{D}}}}} $ | (24) |
作如下变量代换，将式(24)线性化
| $ {p_{{\rm{fD}}}} = - \frac{1}{\beta }{\rm{ln}}(x + 1) $ | (25) |
式中：x为无因次裂缝压力的中间代换变量。
再经基于tD的Laplace变换后，可得到Laplace空间中的数学模型
| $ \frac{{{\partial ^2}\bar x}}{{\partial {\kern 1pt} r_{\rm{D}}^2}} + \frac{1}{{{r_{\rm{D}}}}}\frac{{\partial {\kern 1pt} \bar x}}{{\partial {\kern 1pt} {r_{\rm{D}}}}} = f(u)\bar x $ | (26) |
式中：x为代换变量x对应的Laplace空间变量。
| $ f(u) = u{\kern 1pt} {{\rm{e}}^{ - 2S}} $ | (27) |
内边界条件
| $ {\left. {\frac{{d{\kern 1pt} \bar x}}{{d{\kern 1pt} {r_{\rm{D}}}}}} \right|_{{r_{\rm{D}}} = 1}} - (u{\kern 1pt} {C_{\rm{D}}} - \beta ){\bar x_{\rm{w}}} = \frac{\beta }{u} $ | (28) |
式中：xw为变量x在井筒壁面位置处的值。
无限大外边界
| $ {\mathop {{\rm{lim}}}\limits_{{r_{\rm{D}}} \to \infty } \bar x = 0} $ | (29) |
恒压边界
| $ {\bar x{|_{{r_{\rm{D}}} = {r_{{\rm{eD}}}}}} = 0} $ | (30) |
封闭边界
| $ {{{\left. {\frac{{\partial \bar x}}{{\partial {r_{\rm{D}}}}}} \right|}_{{r_{\rm{D}}} = {r_{{\rm{eD}}}}}} = 0} $ | (31) |
用式(22)替换式(27)，就可得到经变量代换线性化后的裂缝孔洞型三孔介质渗流偏微分控制方程
| $ {\frac{{{\partial ^2}\bar x}}{{\partial {\kern 1pt} r_{\rm{D}}^2}} + \frac{1}{{{r_{\rm{D}}}}}\frac{{\partial {\kern 1pt} \bar x}}{{\partial {\kern 1pt} {r_{\rm{D}}}}} = f(u)\bar x} $ | (32) |
| $ {f(u) = \left( {\frac{{{\omega _{\rm{m}}}{\lambda _{{\rm{mf}}}}}}{{u{\kern 1pt} {\omega _{\rm{m}}} + {\lambda _{{\rm{mf}}}}}} + \frac{{{\omega _{\rm{v}}}{\lambda _{{\rm{vf}}}}}}{{u{\kern 1pt} {\omega _{\rm{v}}} + {\lambda _{{\rm{vf}}}}}} + {\omega _{\rm{f}}}} \right)u{\kern 1pt} {\kern 1pt} {{\rm{e}}^{ - 2S}}} $ | (33) |
式(32)的通解为
| $ \bar x = A{I_0}[{r_{\rm{D}}}\sqrt {f(u)} ] + B{\kern 1pt} {K_0}[{r_{\rm{D}}}\sqrt {f(u)} ] $ | (34) |
在井壁处，rD = 1，x = xw，有
| $ {I_0}[\sqrt {f(u)} ]A + {K_0}[\sqrt {f(u)} ]B - {\bar x_{\rm{w}}} = 0 $ | (35) |
将式(34)代入式(28)，得
| $ \begin{array}{*{20}{c}} {\sqrt {f(u)} {I_1}[\sqrt {f(u)} ]A - \sqrt {f(u)} {K_1}[\sqrt {f(u)} ]B - }\\ {({C_{\rm{D}}}u + \beta ){{\bar x}_{\rm{w}}} = \frac{\beta }{u}} \end{array} $ | (36) |
(1) 无穷大外边界模型的解
故无穷大外边界条件模型在Laplace空间的通解为
| $ \bar x = - \frac{\beta }{u}\frac{{{K_0}[{r_{\rm{D}}}\sqrt {f(u)} ]}}{{\{ ({C_{\rm{D}}}u + \beta ){K_0}[\sqrt {f(u)} ] + \sqrt {f(u)} {K_1}[\sqrt {f(u)} ]\} }} $ | (37) |
井底rD = 1，当p = pwf时，pD = pwD，则x = xw，故无穷大地层模型井底动态压力在Laplace空间为
| $ {\bar x_{\rm{w}}} = {\left. {\bar x} \right|_{{r_{\rm{D}}} = 1}} = - \frac{\beta }{u}\frac{{{K_0}(\sqrt \xi )}}{{\{ ({C_{\rm{D}}}u + \beta ){K_0}[\sqrt {f(u)} ] + \sqrt {f(u)} {K_1}[\sqrt {f(u)} ]\} }} $ | (38) |
(2) 恒压外边界模型的解
将式(34)代入恒压外边界条件式(30)，有
| $ {I_0}[{r_{{\rm{eD}}}}\sqrt {f(u)} ]A + {K_0}[{r_{{\rm{eD}}}}\sqrt {f(u)} ]B = 0 $ | (39) |
联解式(35)，(36)，(39)，将其写为矩阵形式，利用克莱姆法则或高斯迭代消元法，即可解出系数A和B，从而得出外边界条件为恒压边界时井底压力动态在Laplace空间的解。
(3) 封闭外边界模型的解
将式(34)代入封闭外边界条件式(31)，有
| $ {I_1}[{r_{{\rm{eD}}}}\sqrt {f(u)} ]A - {K_1}[{r_{{\rm{eD}}}}\sqrt {f(u)} ]B = 0 $ | (40) |
联解式(35)，(36)，(40)，将其写为矩阵形式，利用克莱姆法则或高斯迭代消元法，即可解出系数A和B，从而得出外边界条件为封闭边界时井底压力动态在Laplace空间的解。
对上述三类外边界条件下Laplace空间的解，利用Stehfest数值算法进行数值反演，得到实空间的解xw，再根据压力pD与x的变换关系式，可求出裂缝孔洞型三孔介质非线性渗流模型实空间井底压力动态响应的数值解，进而编程绘制pwD~tD/CD与p'wDtD/CD~tD/CD的无因次双对数渗流特征曲线。
2 非线性渗流特征分析 2.1 非线性渗流过程分析图 1为受非线性项系数影响的无穷大裂缝孔洞型三孔介质地层拟稳态窜流非线性渗流模型的典型曲线。该曲线主要受无因次二次梯度非线性项系数β控制。在其他参数一定时(CD = 10，S = 10，ωf = 0.001，ωm = 0.989，ωv = 0.01，λmf = 1×10-10，λvf = 1×10-8)，通过分别设定β值为0，0.01和0.05，来分析非线性渗流特征曲线与线性渗流特征曲线的差异。当β = 0时，该曲线则蜕化为线性模型的渗流典型曲线(图 1曲线①)。从图 1可以看出，非线性渗流特征曲线向下偏离线性渗流特征曲线，随β增大，偏离量越大。图 1主要有如下7个流动阶段：
|
下载原图 图 1 受非线性项系数(β)影响的三孔介质拟稳态窜流模型典型曲线 Fig. 1 Type curves of pressure transients controlled by varying nonlinear coefficient(β) |
第Ⅰ阶段：纯井筒储集阶段。此流动阶段的曲线特征不受非线性二次压力梯度项的影响。
第Ⅱ阶段：井筒储集和表皮效应阶段。在此阶段，非线性二次压力梯度项的影响开始出现，非线性渗流特征曲线也开始偏离线性渗流特征曲线。
第Ⅲ阶段：早期裂缝径向流阶段。裂缝系统中的流体径向流入井筒，而溶洞系统和基质系统中的流体还未开始流动。与第Ⅱ阶段相比，线性模型和非线性模型典型曲线的差异更为明显。
第Ⅳ阶段：溶洞系统向裂缝系统窜流阶段。压力导数曲线呈“V形”，反映了溶洞系统向裂缝系统的窜流。线性模型和非线性模型典型曲线的差异进一步增加。
第Ⅴ阶段：裂缝系统和溶洞系统整体径向流阶段。非线性二次压力梯度项的影响越来越明显，其压力导数曲线不符合“0.5线”定律，而是位于“0.5线”下方。
第Ⅵ阶段：基质系统向裂缝系统窜流阶段。压力导数曲线呈“V形”，表明基质系统开始向裂缝系统窜流。线性模型和非线性模型特征曲线的差异越来越大。
第Ⅶ阶段：裂缝系统、溶洞系统和基质系统整体径向流阶段。非线性二次压力梯度项的影响越来越明显，其压力导数曲线不符合“0.5线”定律，而是位于“0.5线”下方，随着时间的推移，压力导数曲线逐渐偏离“0.5线”，如图 1中曲线②和曲线③所示。随着时间的推移，非线性模型的压力导数曲线逐渐偏离线性模型的压力导数曲线。
非线性三孔介质曲线的外边界反映特征与常规线性模型曲线的外边界反映特征类似，此处不再赘述。
该曲线除了受β控制，同时还受别的因素影响。图 2为受溶洞向裂缝窜流的窜流系数(λvf)的非线性渗流模型典型曲线。在其他参数一定时(CD= 10，S = 1，ωf = 0.001，ωm = 0.989，ωv = 0.01，β = 0.05)，分别设定λvf值为1×10-7，1×10-8和1×10-9。随λvf增加，溶洞向裂缝窜流时间推迟，在图 2中表现为压力导数曲线的“V形”向左移动。
|
下载原图 图 2 受溶洞向裂缝窜流的窜流系数(λvf)影响的典型曲线 Fig. 2 Type curves of pressure transients controlled by varying inter-porosity flow factor(λvf) |
图 3为受基质向裂缝窜流的窜流系数(λmf)的非线性渗流模型典型曲线。在其他参数一定时(CD= 10，S = 1，ωf = 0.001，ωm = 0.989，ωv = 0.01，β = 0.05)，分别设定λmf值为1×10-9，1×10-10和1×10-11。溶洞向裂缝窜流时间推迟，在图 3中表现为压力导数曲线的“V形”向右移动。
|
下载原图 图 3 受基质向裂缝窜流的窜流系数(λmf)影响的典型曲线 Fig. 3 Type curves of pressure transients controlled by varying inter-porosity flow factor(λmf) |
图 4为受裂缝系统弹性储容比(ωf)、基质系统弹性储容比(ωm)影响的非线性渗流模型典型曲线。其他参数一定时(CD = 10，S = 1，ωv = 0.1，λvf = 1× 10-10，λmf = 1× 10-8，β = 0.05)，分别设定ωf值为0.001，0.01和0.08，相应的ωm值为0.889，0.89和0.82。随ωf的增加，在图 4中表现为压力导数曲线的“V形”变浅、变窄。
|
下载原图 图 4 受裂缝系统弹性储容比(ωf)影响的典型曲线 Fig. 4 Type curves of pressure transients controlled by fluid capacitance coefficient(ωf) |
图 5为受溶洞系统弹性储容比(ωv)、基质系统弹性储容比(ωm)影响的非线性渗流模型典型曲线。其他参数一定时(CD = 10，S = 1，ωf = 0.001，λvf = 1× 10-10，λmf= 1×10-8，β = 0.05)，分别设定ωv值为0.01，0.01和0.5，相应的ωm值为0.989，0.889和0.499。随ωv的增加，在图 5中表现为压力导数曲线的第1个“V形”变得更深、更宽，第2个“V形”变得更窄、更浅。
|
下载原图 图 5 受溶洞系统弹性储容比(ωv)影响的典型曲线 Fig. 5 Type curves of pressure transients controlled by fluid capacitance coefficient(ωv) |
为了定量描述非线性项系数的影响，定义绝对偏差(DV)与相对偏差(RDV)：
| $ DV = \left| {{\rm{常规线性模型的值}} - {\rm{非线性模型的值}}} \right| $ | (41) |
| $ RDV = \left| {\frac{{DV}}{{{\rm{常规线性模型的值}}}}} \right| \times 100\% $ | (42) |
其他参数一定时(CD= 10，S = 1，ωf = 0.001，ωm = 0.989，ωv = 0.01，λmf = 1×10-10，λvf = 1×10-8)，分别设定β值为0.01和0.05，分析非线性项系数对渗流典型曲线的影响(表 1-2)。
|
下载CSV 表 1 非线性与线性渗流典型曲线的定量偏差数据(β=0.01) Table 1 Theoretical offset of type curves between the linear and nonlinear models(β=0.01) |
|
|
下载CSV 表 2 非线性与线性渗流典型曲线的定量偏差数据(β=0.05) Table 2 Theoretical offset of type curves between the linear and nonlinear models(β=0.05) |
当β=0.01时，由表 1中数据知：当tD/CD=103时，无因次压力曲线的绝对偏差和相对偏差分别为0.579 8和8.83%，无因次压力导数曲线的绝对偏差和相对偏差分别为0.085 4和17.37%；当tD/CD=107时，无因次压力曲线的绝对偏差和相对偏差分别为1.4和13.43%，无因次压力导数曲线的绝对偏差和相对偏差分别为0.089 7和25.40%。
当β=0.05，由表 2中数据知：当tD/CD=103时，无因次压力曲线的绝对偏差和相对偏差分别为2.135 0和32.52%，无因次压力导数曲线的绝对偏差和相对偏差分别为0.272 4和55.40%；当tD/CD= 107时，无因次压力曲线的绝对偏差和相对偏差分别为4.527 2和43.43%，无因次压力导数曲线的绝对偏差和相对偏差分别为0.241 1和68.28%。
绝对偏差和相对偏差均随着无因次时间的增加而增加。在同一时刻，压力导数的相对偏差值大于压力的相对偏差值，例如，当tD/CD=103，β=0.01时，压力导数的相对偏差值就(17.37%)大于压力的相对偏差(8.83%)。绝对偏差和相对偏差值均随β的增加而增加，例如，当tD/CD=103时，β=0.01的压力相对偏差值为8.83%，而β=0.05的压力相对偏差值为32.52%。
3 结论(1) 针对裂缝孔洞型三孔介质储层建立的二次压力梯度非线性渗流模型比不考虑二次压力梯度的常规线性渗流模型能更精准地描述储层的渗流特征。
(2) 随油井生产时间的增加，裂缝孔洞型三孔介质储层非线性渗流特征曲线逐渐偏离常规线性渗流特征曲线，且低于常规线性渗流特征曲线。
(3) 非线性与线性渗流特征曲线之间的偏离量可用“绝对偏差”和“相对偏差”等2个参数来定量描述，其偏差随着二次压力梯度非线性系数的增加而增加。当非线性系数为0时，非线性三孔渗流模型可蜕化为常规线性三孔渗流模型。
| [1] |
田瀚, 沈安江, 张建勇, 等. 一种缝洞型碳酸盐岩储层胶结指数m计算新方法. 地球物理学报, 2019, 62(6): 2276-2285. TIAN H, SHEN A J, ZHANG J Y, et al. New calculation method of cementation exponent m for crack-vuggy carbonate reservoirs. Chinese Journal of Geophysics, 2019, 62(6): 2276-2285. |
| [2] |
郑松青, 李阳, 张宏方. 碳酸盐岩缝洞型油藏网络模型. 中国石油大学学报(自然科学版), 2010, 34(3): 72-75. ZHENG S Q, LI Y, ZHANG H F. Fracture-cavity network model for fracture-cavity carbonate reservoir. Journal of China University of Petroleum, 2010, 34(3): 72-75. |
| [3] |
贾永禄, 孙高飞, 聂仁仕, 等. 四重介质油藏渗流模型与试井曲线. 岩性油气藏, 2016, 28(1): 123-127. JIA Y L, SUN G F, NIE R S, et al. Flow model and well test curves for quadruple-media reservoirs. Lithologic Reservoirs, 2016, 28(1): 123-127. |
| [4] |
熊钰, 叶海峰, 蔡明金, 等. 井打在大尺度裂缝上的缝洞型油气藏试井模型. 科学技术与工程, 2017, 17(34): 49-54. XIONG Y, YE H F, CAI M J, et al. A well test model for wells drilled in large-scale fracture of fractured-vuggy carbonate reservoirs. Science Technology and Engineering, 2017, 17(34): 49-54. |
| [5] |
程飞. 缝洞型碳酸盐岩油藏储层类型动静态识别方法:以塔里木盆地奥陶系为例. 岩性油气藏, 2017, 29(3): 76-82. CHENG F. Integrated dynamic and static identification method of fractured-vuggy carbonate reservoirs:a case from the Ordovician in Tarim Basin. Lithologic Reservoirs, 2017, 29(3): 76-82. |
| [6] |
姜瑞忠, 沈泽阳, 崔永正, 等. 双重介质低渗油藏斜井压力动态特征分析. 岩性油气藏, 2018, 30(6): 131-137. JIANG R J, SHEN Z Y, CUI Y Z, et al. Dynamical characteristics of inclined well in dual medium low permeability reservoir. Lithologic Reservoirs, 2018, 30(6): 131-137. |
| [7] |
杜鑫, 李冬梅, 徐燕东, 等. 井洞相连的缝洞型油藏试井新模型. 水动力学研究与进展A辑, 2018, 33(5): 12-21. DU X, LI D M, XU Y D, et al. A new analytical well test model for fractured vuggy reservoirs with vug connecting with wellbore. Chinese Journal of Hydrodynamics, 2018, 33(5): 12-21. |
| [8] |
刘建峰, 郑开兵. 一种典型缝洞型油藏双重介质试井特征. 辽宁化工, 2019, 48(2): 56-58. LIU J F, ZHENG K B. Well testing characteristics of a typical fractured-vuggy reservoir with dual medium. Liaoning Chemical Industry, 2019, 48(2): 56-58. |
| [9] |
王蓓, 刘向君, 司马立强, 等. 磨溪龙王庙组碳酸盐岩储层多尺度离散裂缝建模技术及其应用. 岩性油气藏, 2019, 31(2): 124-133. WANG B, LIU X J, SIMA L Q, et al. Multi-scale discrete fracture modeling technology for carbonate reservoir of Longwangmiao Formation in Moxi area and its application. Lithologic Reservoirs, 2019, 31(2): 124-133. |
| [10] |
姜瑞忠, 张春光, 郜益华, 等. 缝洞型碳酸盐岩油藏水平井分形非线性渗流. 岩性油气藏, 2019, 31(6): 118-126. JIANG R Z, ZHANG C G, GAO Y H, et al. Fractal nonlinear seepage model of horizontal wells in fractured-vuggy carbonate reservoirs. Lithologic Reservoirs, 2019, 31(6): 118-126. |
| [11] |
聂仁仕, 贾永禄, 石国新, 等. 二次压力梯度非线性渗流理论与应用. 大庆石油地质与开发, 2010, 29(5): 96-103. NIE R S, JIA Y L, SHI G X, et al. Theory and its application for nonlinear percolation with quadratic pressure gradient. Petroleum Geology and Oilfield Development in Daqing, 2010, 29(5): 96-103. |
| [12] |
贾永禄, 聂仁仕, 王永恒, 等. 二次梯度非线性渗流问题求解. 石油天然气学报, 2008, 30(4): 119-123. JIA Y L, NIE R S, WANG Y H, et al. Solution for the problem of quadratic gradient nonlinear percolation. Journal of Oil and Gas Technology, 2008, 30(4): 119-123. |
| [13] |
贾永禄, 聂仁仕, 陈珂, 等. 考虑二次梯度影响的渗流模型与试井样版曲线. 西南石油大学学报, 2007, 29(5): 69-71. JIA Y L, NIE R S, CHEN K, et al. The flow model considering the quadratic pressure gradient and the well test type curves. Journal of Southwest Petroleum University, 2007, 29(5): 69-71. |
| [14] |
NIE R S, DING Y. Research on the nonlinear spherical percolation model with quadratic pressure gradient and its percolation characteristics. Natural Science, 2010, 2(2): 98-105. |
| [15] |
张林, 贾永禄, 张福祥, 等. 考虑二次压力梯度非线性渗流三区复合油藏模型分析. 大庆石油学院学报, 2011, 35(5): 54-59. ZHANG L, JIA Y L, ZHANG F X, et al. Analysis of considering quadratic pressure gradient nonlinear percolation in 3-zone composite reservoir. Journal of Daqing Petroleum Institute, 2011, 35(5): 54-59. |
| [16] |
张强, 王永清, 杨玲智, 等. 考虑启动压力和二次梯度的应力敏感油藏稳定渗流模型. 大庆石油学院学报, 2011, 35(5): 22-24. ZHANG Q, WANG Y Q, YANG L Z, et al. Stable flow models for stress sensitive reservoir taking into consideration the start-up pressure gradient and the quadratic gradient. Journal of Daqing Petroleum Institute, 2011, 35(5): 22-24. |
| [17] |
WANG X L, FAN X Y, HE Y M, et al. A nonlinear model for fluid flow in a multiple-zone composite reservoir including the quadratic gradient term. Journal of Geophysics and Engineering, 2013, 10(4): 1-11. |
| [18] |
GUO J C, NIE R S. Nonlinear flow model for well production in an underground formation. Nonlinear Processes in Geophysics, 2013, 20(3): 311-327. |
| [19] |
王美楠, 尹洪军, 钟会影. 考虑二次梯度项及动边界的低渗透变形介质油藏渗流规律. 大庆石油地质与开发, 2014, 33(1): 75-80. WANG M N, YIN H J, ZHONG H Y. Analysis on the seepage flow rules for low-permeability oil reservoir with deformed media considering quadratic gradient term and dynamic boundary. Petroleum Geology and Oilfield Development in Daqing, 2014, 33(1): 75-80. |
| [20] |
LU C, NIE R S, GUO J C, et al. Modeling the nonlinear flow behavior for horizontal well production in an underground porousmedia formation. Journal of Porous Media, 2015, 18(1): 43-56. |