0 引言

缝洞型油藏的孔洞、基质及裂缝系统相互影响，为了描述其储层复杂的孔隙结构，出现了三重介质模型的概念^[1-2]。描述缝洞型碳酸盐岩油藏的渗流特征，须充分考虑裂缝系统的分形特征、致密基质块的非线性渗流特征，以及储层的应力敏感性^[3-5]。

国内外学者针对三重介质模型进行了深入的研究。姚军等^[6]创新变井筒储存三重介质油藏试井理论，对介质间窜流、弹性储容比及变井储参数对压力响应的影响进行了研究。Camacho-Velazquez等^[7]对所建立的三重孔隙模型考虑了溶洞系统与井筒连通和不连通2种情况下的解，并对样板曲线进行了详细的研究。李成勇等^[8]建立了三重介质油藏的水平井各向异性渗流数学模型，并对水平井的压力动态进行了系统的研究。张冬丽等^[9]建立了三重介质数值试井模型，很好地反映出了缝洞型储层的非均质性特征，该数值试井模型可进一步考虑缝洞型油藏中具体的裂缝和溶洞，进而推广至两相流试井。Gomez等^[10]利用全局最优化解法对缝洞型储层试井曲线解释方法进行研究，提高了计算效率和试井解释精度。Li等^[11]对不同边界条件下缝洞型储层的动态特征进行了研究，形成了一套系统的刻画缝洞型储层动态特征的方法。不仅如此，分形理论在储层渗流领域的应用也取得了较大的进展。Chang等^[12]提出了分形裂缝网络嵌入基质的表征公式，并通过对扩散方程系数的适当修正来描述分形系统中的单相流动。Acuna等^[13]将分形理论应用于裂缝性储层试井中，并基于试井解释明确了分形参数的物理意义。Razminia等^[14]利用分形几何和分数阶微积分的概念分析了两区复合储层的数学模型，并推导了径向复合系统中井底瞬态压力解。张本健等^[15]针对裂缝性气藏试井渗流问题，使用分形网络模拟气藏裂缝系统并将分形网络嵌入到气藏中，计算出试井模型的流动动态特征典型曲线，结果表明分形网络能够很好地模拟裂缝性气藏的渗流特征。

虽然目前对三重介质模型的研究较多，但仍存在无法精确表征储层非均质性特征的问题，而分形理论的应用使得非均质储层，尤其是裂缝的表征更灵活且更贴近实际，但该理论在缝洞型储层中的应用研究较少且未成体系。此外，对缝洞型储层致密基质块非线性渗流特征的研究也不可忽视。研究低渗储层最常用的渗流模型主要有拟启动压力梯度模型、分段模型及连续模型，这3种模型均由不同数学函数拟合实验数据所得，缺乏物理背景，表征不够灵活。

针对缝洞型碳酸盐岩储层渗流模型目前存在的问题，综合考虑裂缝分形特征和基质系统非线性渗流机理，完善非线性渗流理论，建立三重介质油藏水平井分形非线性渗流模型；基于有限元原理求解水平井动态压力，并对比各模型的压力动态曲线，总结渗流规律。在分析非线性参数、分形指数及水平井长度等相关参数的敏感性之后，将所建立模型进行实际井动态压力拟合并解释相关渗流参数，以期验证模型的准确性及实用价值。

1 缝洞型碳酸盐岩油藏渗流机理 1.1 基质非线性渗流

将低渗储层微观渗流机理、边界层流动理论和毛细管渗流模型相结合，简化并改进非线性渗流模型。考虑边界层厚度与流体屈服应力值，在HagenPoisseuille定律的基础上进行修正^[16]可得

$ v = \frac{K}{\mu }{\left( {1 - \frac{\delta }{{{r_0}}}} \right)^4}\left[ {1 - \frac{{8{\tau _0}}}{{3{r_0}\left( {1 - \frac{\delta }{{{r_0}}}} \right)\nabla p}}} \right]\nabla p $

(1)

式中：v为通过岩心的流速，cm/s；K为岩心渗透率，mD；μ为流体黏度，mPa·s；r₀为毛管半径，mm；δ为边界层厚度，mm；τ₀为流体屈服应力，MPa；▽_p为压力梯度，MPa/m。

对于相同毛管，δ/r₀与压力梯度成反比^[17]，故假设δ/r₀= a₁/▽_p；对于同一种流体，屈服应力值保持不变，假设8τ₀/3r₀=a₂。将式（1）变形并简化可得

$ v = \frac{K}{\mu }\left[ {1 - \frac{{4{a_1} + {a_2}}}{{\nabla p}} + \frac{{6a_1^2 + 3{a_1}{a_2}}}{{\nabla p\left( {\nabla p - {a_1}} \right)}}} \right]\nabla p $

(2)

式中：a₁，a₂均为实验拟合得到的常数值。

将表征多项式c₁=4a₁+a₂, c₂=a₁, c₃=6a₁²+3a₁a₂和条件"▽_p=0时, v=0"代入式(2), 并约去c₃, 得到简化后的模型为

$ \begin{array}{l} v = \frac{K}{\mu }\left[ {1 - \frac{{{c_1}}}{{\nabla p}} - \frac{{{c_1}{c_2}}}{{\nabla p\left( {\nabla p - {c_2}} \right)}}} \right]\nabla p = \\ \;\;\;\;\;\frac{K}{\mu }\left( {1 - \frac{{{c_1}}}{{\nabla p - {c_2}}}} \right)\nabla p \end{array} $

(3)

式中：c₁和c₂是通过对实验数据进行非线性拟合得到的参数，c₁反映了流体的屈服应力以及边界层对渗流的影响，c₂反映了边界层对渗流的影响。

当c₁=0时，式（3）简化为达西模型；当c₂=0时，式（3）简化为拟启动压力梯度模型。c₁/(▽p-c₂)相当于启动压力梯度项。当v > 0时，▽p > c₁ + c₂，最小启动压力梯度为▽p_min = c₁ + c₂。启动压力梯度的存在与否以及取值大小均取决于c₁，c₂的值，即取决于矿场实际^[18]。

1.2 裂缝分形特征

利用分形维数及异常扩散系数来表征分形裂缝中的渗透率和孔隙度，能够精确地反映地下储层复杂的裂缝形态。设分形体内流体聚集在某节点处，且节点体积处处相等，则渗流节点的数量可结合节点密度求出，即

$ \alpha = \frac{{{K_{{\rm{fw}}}}{{\left( {{D_{\rm{f}}}/{d_{\rm{s}}}} \right)}^2}}}{{{\varphi _{{\rm{fw}}}}\mu {C_{{\rm{ft}}}}r_{\rm{w}}^2}} $

(4)

式中：D_f为分形维数；α为节点的数量；d_s为谱维数；r_w为井筒半径，m；K_fw为生产井处的渗透率，mD；φ_fw为生产井处的孔隙度；C_ft为裂缝综合压缩系数，MPa^-1。

考虑维数为d的欧氏空间内嵌入分形体，可得到用分形维数表征的孔隙度表达式为

$ {\varphi _{\rm{f}}} = \frac{{{V_{\rm{f}}}}}{V} = \frac{{{V_{\rm{s}}}\alpha {r^{{D_{\rm{f}}} - 1}}\Delta r}}{{B\alpha {r^{d - 1}}\Delta r}} = \frac{{\alpha {V_{\rm{s}}}}}{B}{r^{{D_{\rm{f}}} - d}} $

(5)

式中：B为描述对称性的几何常量；V为储层总体积，m³；V_f为储层孔隙体积，m³；V_s为储集流体的节点体积，m³；r为距生产井的距离，m。

裂缝孔隙度主要受孔隙空间聚集方式（用分形维数表征）的影响，且与有效上覆压力呈指数关系，所以考虑应力敏感性的裂缝分形孔隙度表达式^[19-20]为

$ {\varphi _{\rm{f}}}(r) = {\varphi _{{\rm{fw}}}}\left( {\frac{r}{{{r_{\rm{w}}}}}} \right){r^{{D_{\rm{f}}} - d}}\exp \left[ { - {C_{\rm{f}}}\left( {{p_{\rm{i}}} - {p_{\rm{f}}}} \right)} \right] $

(6)

式中：C_f为裂缝孔隙压缩系数，MPa^-1；p_i为原始地层压力，MPa；p_f为裂缝系统压力，MPa。

裂缝渗透率不仅取决于孔隙空间的聚集方式，还与裂缝之间的连通性（用异常扩散系数表征）密不可分。储层压力变化会引起裂缝张合状态的变化，从而影响储层渗透率，在对应力敏感地层不稳定试井的研究中发现，渗透率与基岩上覆压力之间也成指数关系。在分形孔隙度定义的基础上引入异常扩散系数便可得到考虑应力敏感性的分形渗透率^[19-20]，即

$\begin{array}{l} {K_{\rm{f}}}(r) = {K_{{\rm{fw}}}}\left( {\frac{r}{{{r_{\rm{w}}}}}} \right){r^{{D_{\rm{f}}} - d - \theta }}\\ \exp \left[ { - {\alpha _{\rm{K}}}\left( {{p_{\rm{i}}} - {p_{\rm{f}}}} \right)} \right] \end{array} $

(7)

式中：α_K为渗透率模数，MPa^-1；θ为异常扩散系数。

2 缝洞型碳酸盐岩分形油藏水平井非线性渗流模型 2.1 物理模型

建立圆形三重介质碳酸盐岩分形油藏水平井物理模型^[21]（图 1），并设该模型符合以下假定条件：①模型由基质、裂缝和溶洞系统组成，基质块向裂缝和溶洞系统发生拟稳态窜流，溶洞系统向裂缝发生拟稳态窜流；②外边界为封闭或定压条件，储层厚度为h，原始压力为p_i，水平井半长为L，水平井纵向深度为z_w，油藏半径为r_e；③考虑储层的表皮系数、井筒储集系数以及裂缝渗透率的各向异性，裂缝系统水平渗透率为K_fh，垂直渗透率为K_fv，基质系统渗透率为K_m，溶洞系统渗透率为K_v；④水平井位于油藏平面中心，产量为Q）；⑤忽略重力与毛管压力的影响；⑥考虑裂缝系统的应力敏感性及裂缝的分形特征，基质向裂缝系统的窜流为非线性渗流。

2.2 数学模型

在笛卡尔坐标系中，将流体流动的质量守恒方程和状态方程与非线性运动方程[式（3）]相结合，可得到缝洞型碳酸盐岩油藏水平井分形非线性渗流模型^[22]，即

$ \begin{array}{*{20}{c}} {\frac{\partial }{{\partial {x_{\rm{D}}}}}\left( {r_{\rm{D}}^{{D_{\rm{f}}} - d - \theta }{{\rm{e}}^{ - {\alpha _{{\rm{KD}}}}{p_{{\rm{fD}}}}}}\frac{{\partial {p_{{\rm{fD}}}}}}{{\partial {x_{\rm{D}}}}}} \right) + \frac{\partial }{{\partial {y_{\rm{D}}}}}\left( {r_{\rm{D}}^{{D_{\rm{f}}} - d - \theta }{{\rm{e}}^{ - {\alpha _{{\rm{KD}}}}{p_{{\rm{fD}}}}}}\frac{{\partial {p_{{\rm{fD}}}}}}{{\partial {y_{\rm{D}}}}}} \right) + \frac{\partial }{{\partial {z_{\rm{D}}}}}\left( {L_{\rm{D}}^2{{\rm{e}}^{ - {\alpha _{{\rm{KD}}}}{p_{{\rm{fD}}}}}}\frac{{\partial {p_{{\rm{fD}}}}}}{{\partial {z_{\rm{D}}}}}} \right) = }\\ {{\omega _{\rm{f}}}{{\left( {{h_{\rm{D}}}{L_{\rm{D}}}} \right)}^2}\frac{{\partial {p_{{\rm{fD}}}}}}{{\partial {t_{\rm{D}}}}} + {\lambda _{{\rm{mf}}}}{\delta _{{\rm{LVDf}}}}\left( {{p_{{\rm{fD}}}} - {p_{{\rm{mD}}}}} \right) + {\lambda _{{\rm{vf}}}}\left( {{p_{{\rm{fD}}}} - {p_{{\rm{vD}}}}} \right)} \end{array} $

(8)

$ {\lambda _{{\rm{mf}}}}{\delta _{{\rm{LVDf}}}}\left( {{p_{{\rm{fD}}}} - {p_{{\rm{mD}}}}} \right) + {\lambda _{{\rm{mv}}}}{\delta _{{\rm{LVDv}}}}\left( {{p_{{\rm{vD}}}} - {p_{{\rm{mD}}}}} \right) = {\left( {{h_{\rm{D}}}{L_{\rm{D}}}} \right)^2}\left( {1 - {\omega _{\rm{f}}} - {\omega _{\rm{v}}}} \right)\frac{{\partial {p_{{\rm{mD}}}}}}{{\partial {t_{\rm{D}}}}} $

(9)

$ {\lambda _{{\rm{vf}}}}\left( {{p_{{\rm{fD}}}} - {p_{{\rm{vD}}}}} \right) + {\lambda _{{\rm{mv}}}}{\delta _{{\rm{LVDv}}}}\left( {{p_{{\rm{vD}}}} - {p_{{\rm{mD}}}}} \right) = {\left( {{h_{\rm{D}}}{L_{\rm{D}}}} \right)^2}{\omega _{\rm{v}}}\frac{{\partial {p_{{\rm{vD}}}}}}{{\partial {t_{\rm{D}}}}} $

(10)

式中：r_D为无因次渗流半径；L_D为无因次水平井半长；h_D为无因次储层厚度；t_D为无因次时间；p_fD为无因次裂缝压力；p_mD为无因次基质压力；p_vD为无因次溶洞压力；δ_LVDf为裂缝的无因次中间变量；δ_LVDv为溶洞的无因次中间变量；λ_mf为基质向裂缝窜流系数；λ_vf为溶洞向裂缝窜流系数；λ_mv为基质向溶洞窜流系数；ω_f为裂缝弹性储能比；ω_v为溶洞弹性储能比；下标D代表无因次；下标m，f，v分别代表基质、裂缝和溶洞系统。

若$\left({{p_{{\rm{AD}}}} - {p_{{\rm{mD}}}}} \right) < \frac{{{c_{{\rm{ID}}}} + {c_{2{\rm{D}}}}}}{2}$, ${\delta _{{\rm{LVDA}}}} = 0$, 否则${\delta _{{\rm{LVDA}}}} = \frac{{1 - {c_{1{\rm{D}}}}}}{{2\left({{p_{{\rm{AD}}}} - {p_{{\rm{mD}}}}} \right) - {c_{2{\rm{D}}}}}}({\rm{A}} = {\rm{f}}, {\rm{v}})$。

初始条件：

$ {p_{1{\rm{D}}}}\left( {{x_{\rm{D}}},{y_{\rm{D}}},{z_{\rm{D}}},{t_{\rm{D}}} = 0} \right) = 0\;\;\;\;(1 = {\rm{m}},{\rm{f}},{\rm{v}}) $

(11)

封闭外边界条件：

$ {\left. {\frac{{\partial {p_{1{\rm{D}}}}}}{{\partial {r_{\rm{D}}}}}} \right|_{{r_D} = {r_{{\rm{eD}}}}}} = {\left. {\frac{{\partial {p_{1{\rm{D}}}}}}{{\partial {z_{\rm{D}}}}}} \right|_{{z_D} = 0,1}} = 0 $

(12)

定压边界条件：

$ {p_{1{\rm{D}}}}\left( {{r_{\rm{D}}} = {r_{{\rm{eD}}}}} \right) = {p_{1{\rm{D}}}}\left( {{z_{\rm{D}}} = 0,1} \right) = 0 $

(13)

内边界条件：

$ \begin{array}{*{20}{c}} {{{\left. {\mathop {\lim }\limits_{{\varepsilon _{{\rm{D}} \to 0}}} \int_{{z_{{\rm{w}}0}} - \frac{{{\varepsilon _{\rm{D}}}}}{2}}^{{z_{{\rm{wD}}}} + \frac{{{\varepsilon _{\rm{D}}}}}{2}} {\left( {r_{\rm{D}}^{{D_{\rm{f}}} - d - \theta }{{\rm{e}}^{ - {\alpha _{{\rm{KD}}}}{p_{{\rm{fD}}}}}}\frac{{\partial {p_{{\rm{fD}}}}}}{{\partial {r_{\rm{D}}}}}} \right)} {\rm{d}}{z_{{\rm{wD}}}}} \right|}_{{r_{\rm{D}}} = {\varepsilon _{\rm{D}}}}} = - \frac{1}{2},}\\ {\left| {{z_{\rm{D}}} - {z_{{\rm{wD}}}}} \right| \le \frac{{{\varepsilon _{\rm{D}}}}}{2}} \end{array} $

(14)

各无因次量定义如下：${\lambda _{{\rm{mf}}}} = {\sigma _{\rm{m}}}\frac{{{K_{\rm{m}}}}}{{{K_{{\rm{fhi}}}}}}{L^2}, {\lambda _{{\rm{mv}}}} = {\sigma _{\rm{m}}}\frac{{{K_{\rm{m}}}}}{{{K_{\rm{v}}}}}{L^2}$, ${\lambda _{{\rm{vf}}}} = {\sigma _{\rm{v}}}\frac{{{K_{\rm{v}}}}}{{{K_{{\rm{fhi}}}}}}{L^2}, {L_{\rm{D}}} = \frac{L}{h}\sqrt {\frac{{{K_{{\rm{fvi}}}}}}{{{K_{{\rm{fhi}}}}}}} $, ${x_{\rm{D}}} = \frac{x}{L}$, ${y_{\rm{D}}} = \frac{y}{L}, {z_{\rm{D}}} = \frac{z}{h}$, ${z_{{\rm{wD}}}} = \frac{{{z_{\rm{w}}}}}{h}, {r_{\rm{D}}} = \frac{r}{L}$, ${h_{\rm{D}}} = \frac{h}{{{r_{\rm{w}}}}}\sqrt {\frac{{{K_{{\rm{fhi}}}}}}{{{K_{{\rm{fvi}}}}}}} $, ${\alpha _{{\rm{KD}}}} = \frac{{1.842 \times {{10}^{ - 3}}Q\mu B}}{{{K_{{\rm{fhi}}}}h}}{\alpha _{\rm{K}}}$, ${\omega _{\rm{f}}} = \frac{{{\varphi _{\rm{f}}}{C_{{\rm{ft}}}}}}{{{\varphi _{\rm{f}}}{C_{{\rm{ft}}}} + {\varphi _m}{C_{{\rm{mt}}}} + {\varphi _{\rm{v}}}{C_{{\rm{vt}}}}}}$, ${\omega _{\rm{v}}} = \frac{{{\varphi _{\rm{v}}}{C_{{\rm{vt}}}}}}{{{\varphi _{\rm{f}}}{C_{{\rm{ft}}}} + {\varphi _{\rm{m}}}{C_{{\rm{mt}}}} + {\varphi _{\rm{v}}}{C_{{\rm{vt}}}}}}$, ${c_D} = \frac{{{K_{{\rm{fhi}}}}h{l_{\rm{c}}}}}{{1.842 \times {{10}^{ - 3}}Q\mu B}}c$, ${\delta _{{\rm{LVDA}}}} = 1 - \frac{{{c_{\rm{D}}}}}{{2\left({{p_{{\rm{AD}}}} - {p_{{\rm{mD}}}}} \right) + {c_{\rm{D}}}}}\quad ({\rm{A}} = {\rm{f}}, {\rm{v}})$, ${C_{\rm{D}}} = \frac{C}{{2{\rm{ \mathit{ π} }}h\left({{\varphi _{{\rm{fw}}}}{C_{{\rm{tf}}}} + {\varphi _{\rm{m}}}{C_{{\rm{tm}}}}} \right)r_{\rm{w}}^2}}$, ${p_{1{\rm{D}}}} = \frac{{{K_{{\rm{fhi}}}}h\left({{p_{\rm{i}}} - {p_1}} \right)}}{{1.842 \times {{10}^{ - 3}}Q\mu B}}(1 = {\rm{m}}, {\rm{f}}, {\rm{v}})$, ${t_{\rm{D}}} = \frac{{{K_{{\rm{fhi}}}}t}}{{\left({{\varphi _{\rm{f}}}{C_{{\rm{ft}}}} + {\varphi _{\rm{m}}}{C_{{\rm{mt}}}} + {\varphi _{\rm{v}}}{C_{{\rm{nt}}}}} \right)\mu r_{\rm{w}}^2}}$。式中：K_fhi为裂缝初始水平渗透率，mD；K_fvi为裂缝初始垂向渗透率，mD；x_D，y_D，z_D为无因次的方向坐标；C_ft，C_mt分别为裂缝与基质系统的综合压缩系数，MPa^-1；φ_f，φ_m分别为裂缝与基质系统的孔隙度；α_KD为无因次渗透率模数；C为井筒储集系数，m³/MPa；δ_LVDA为无因次中间变量；c_1D，c_2D为无因次非线性渗流参数；l_c为特征长度，m；z_wD，ε_D，r_eD分别为无因次水平井纵向深度、微变量和泄流半径。

2.3 模型求解

对于缝洞型碳酸盐岩油藏水平井分形非线性渗流模型，首先采用隐式求解裂缝系统，应用Galerkin法，得到裂缝系统的有限元方程为

$ \begin{array}{*{20}{c}} {\iiint\limits_{{\Omega _{\text{e}}}} {{N_i}\left[ {\frac{\partial }{{\partial {x_{\text{D}}}}}\left( {r_{\text{D}}^{{D_{\text{f}}} - d - \theta }{{\text{e}}^{ - {\alpha _{{\text{KD}}}}{p_{{\text{fD}}}}}}\frac{{\partial {p_{{\text{fD}}}}}}{{\partial {x_{\text{D}}}}}} \right) + \frac{\partial }{{\partial {y_{\text{D}}}}}\left( {r_{\text{D}}^{{D_{\text{f}}} - d - \theta }{{\text{e}}^{ - {\alpha _{{\text{KD}}}}{p_{{\text{fD}}}}}}\frac{{\partial {p_{{\text{fD}}}}}}{{\partial {y_{\text{D}}}}}} \right)} \right]{\text{d}}V} + \iiint\limits_{{\Omega _{\text{e}}}} {{N_i}\left[ {\frac{\partial }{{\partial {z_{\text{D}}}}}\left( {{{\text{e}}^{ - {\alpha _{{\text{KD}}}}{p_{{\text{fD}}}}}}L_{\text{D}}^2\frac{{\partial {p_{{\text{fD}}}}}}{{\partial {z_{\text{D}}}}}} \right)} \right]{\text{d}}V} = } \\ {\iiint\limits_{{\Omega _{\text{e}}}} {{N_i}{{\left( {{L_{\text{D}}}{h_{\text{D}}}} \right)}^2}{\omega _{\text{f}}}\frac{{\partial {p_{{\text{fD}}}}}}{{\partial {t_{\text{D}}}}}{\text{d}}V} + \iiint\limits_{{\Omega _{\text{e}}}} {{N_i}{\delta _{{\text{LVDf}}}}{\lambda _{{\text{mf}}}}\left( {{p_{{\text{fD}}}} - {p_{{\text{mD}}}}} \right){\text{d}}V} + \iiint\limits_{{\Omega _{\text{e}}}} {{N_i}{\lambda _{{\text{vf}}}}\left( {{p_{{\text{fD}}}} - {p_{{\text{vD}}}}} \right){\text{d}}V}} \end{array} $

(15)

式中：N_i为形函数的分量，i=1，2，…，n。

以格林公式为基础，通过分部积分可得到内部单元和封闭条件外边界单元的有限元方程为

$ \begin{array}{*{20}{c}} {\iiint\limits_{{\Omega _{\text{e}}}} {r_{\text{D}}^{{D_{\text{f}}} - d - \theta }{{\text{e}}^{ - {\alpha _{{\text{KD}}}}{p_{{\text{fD}}}}}}\left( {\frac{{\partial {N_i}}}{{\partial {x_{\text{D}}}}}\frac{{\partial {p_{{\text{fD}}}}}}{{\partial {x_{\text{D}}}}} + \frac{{\partial {N_i}}}{{\partial {y_{\text{D}}}}}\frac{{\partial {p_{{\text{fD}}}}}}{{\partial {y_{\text{D}}}}}} \right){\text{d}}V} + \iiint\limits_{{\Omega _{\text{e}}}} {{{\text{e}}^{ - {\alpha _{{\text{KD}}}}{p_{{\text{fD}}}}}}\left( {L_{\text{D}}^2\frac{{\partial {N_i}}}{{\partial {z_{\text{D}}}}}\frac{{\partial {p_{{\text{fD}}}}}}{{\partial {z_{\text{D}}}}}} \right){\text{d}}V} + \iiint\limits_{{\Omega _{\text{e}}}} {{N_i}{\delta _{{\text{LVDf}}}}{\lambda _{{\text{mf}}}}{p_{{\text{fD}}}}{\text{d}}V} + } \\ {\iiint\limits_{{\Omega _{\text{e}}}} {{N_i}{\lambda _{{\text{vf}}}}{p_{{\text{fD}}}}{\text{d}}V} + \iiint\limits_{{\Omega _{\text{e}}}} {{N_i}{{\left( {{L_{\text{D}}}{h_{\text{D}}}} \right)}^2}{\omega _{\text{f}}}\frac{{\partial {p_{{\text{fD}}}}}}{{\partial {t_{\text{D}}}}}{\text{d}}V} = \iiint\limits_{{\Omega _{\text{e}}}} {{N_i}{\delta _{{\text{LVDf}}}}{\lambda _{{\text{mf}}}}{p_{{\text{mD}}}}{\text{d}}V} + \iiint\limits_{{\Omega _{\text{e}}}} {{N_i}{\lambda _{{\text{vf}}}}{p_{{\text{vD}}}}{\text{d}}V}} \end{array} $

(16)

由此可得，单元有限方程的矩阵形式为

$ \begin{gathered} \iiint\limits_{{\Omega _{\text{e}}}} {\left[ {\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{B}}^T}\mathit{\boldsymbol{GB}} + {\lambda _{{\text{mf}}}}\mathit{\boldsymbol{DN}}{\mathit{\boldsymbol{N}}^T} + {\lambda _{{\text{vf}}}}\mathit{\boldsymbol{EN}}{\mathit{\boldsymbol{N}}^T} + {\omega _{\text{f}}}{{\left( {{L_{\text{D}}}{h_{\text{D}}}} \right)}^2}\frac{{\mathit{\boldsymbol{N}}{\mathit{\boldsymbol{N}}^T}}}{{\Delta {t_{\text{D}}}}}} \right]{\text{d}}V\mathit{\boldsymbol{P}}_{\text{f}}^n} = \\ \iiint\limits_{{\Omega _{\text{e}}}} {{\omega _{\text{f}}}{{\left( {{L_{\text{D}}}{h_{\text{D}}}} \right)}^2}\frac{{\mathit{\boldsymbol{N}}{\mathit{\boldsymbol{N}}^T}}}{{\Delta {t_{\text{D}}}}}{\text{d}}V\mathit{\boldsymbol{P}}_{\text{f}}^{n - 1}} + \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\iiint\limits_{{\Omega _{\text{e}}}} {{\lambda _{{\text{mf}}}}\mathit{\boldsymbol{DN}}{\mathit{\boldsymbol{N}}^T}{\text{d}}V\mathit{\boldsymbol{P}}_{\text{m}}^{n - 1}} + \iiint\limits_{{\Omega _{\text{e}}}} {{\lambda _{{\text{vf}}}}\mathit{\boldsymbol{EN}}{\mathit{\boldsymbol{N}}^T}{\text{d}}V\mathit{\boldsymbol{P}}_{\text{v}}^{n - 1}} \hfill \\ \end{gathered} $

(17)

假设E为单位矩阵，其他各矩阵的表达形式为

$ \mathit{\boldsymbol{B}} = \left[ {\begin{array}{*{20}{l}} {\frac{{\partial {N_1}}}{{\partial {x_{\rm{D}}}}}}&{\frac{{\partial {N_2}}}{{\partial {x_{\rm{D}}}}}}& \cdots &{\frac{{\partial {N_n}}}{{\partial {x_{\rm{D}}}}}}\\ {\frac{{\partial {N_1}}}{{\partial {y_{\rm{D}}}}}}&{\frac{{\partial {N_2}}}{{\partial {y_{\rm{D}}}}}}& \cdots &{\frac{{\partial {N_n}}}{{\partial {y_{\rm{D}}}}}}\\ {\frac{{\partial {N_1}}}{{\partial {z_{\rm{D}}}}}}&{\frac{{\partial {N_2}}}{{\partial {z_{\rm{D}}}}}}& \cdots &{\frac{{\partial {N_n}}}{{\partial {z_{\rm{D}}}}}} \end{array}} \right],\mathit{\boldsymbol{N}} = \left[ {\begin{array}{*{20}{c}} {{N_1}}\\ {{N_2}}\\ \vdots \\ {{N_n}} \end{array}} \right],{\mathit{\boldsymbol{P}}_{\rm{f}}} = \left[ {\begin{array}{*{20}{c}} {{p_{{\rm{fD}}1}}}\\ {{p_{{\rm{fD}}2}}}\\ \vdots \\ {{p_{{\rm{fD}}n}}} \end{array}} \right],{\mathit{\boldsymbol{P}}_{\rm{m}}} = \left[ {\begin{array}{*{20}{c}} {{p_{{\rm{mD}}1}}}\\ {{p_{{\rm{mD}}2}}}\\ \vdots \\ {{p_{{\rm{mD}}n}}} \end{array}} \right],{\mathit{\boldsymbol{P}}_{\rm{v}}} = \left[ {\begin{array}{*{20}{c}} {{p_{{\rm{vD1}}}}}\\ {{p_{{\rm{vD}}2}}}\\ \vdots \\ {{p_{{\rm{vD}}n}}} \end{array}} \right], $

$\begin{array}{l} \mathit{\boldsymbol{C}} = \left[ {\begin{array}{*{20}{c}} {{{\rm{e}}^{ - {\alpha _{{\rm{KD}}}}{p_{{\rm{fD}}1}}}}}&{}&{}&{}\\ {}&{{{\rm{e}}^{ - {\alpha _{{\rm{KD}}}}{P_{{\rm{fD}}2}}}}}&{}&{}\\ {}&{}& \ddots &{}\\ {}&{}&{}&{{{\rm{e}}^{ - {\alpha _{{\rm{KD}}}}{p_{{\rm{fD}}n}}}}} \end{array}} \right],\\ {\mathit{\boldsymbol{H}}_{\rm{A}}} = \left[ {\begin{array}{*{20}{c}} {{\delta _{{\rm{LVDA}}1}}}&{}&{}&{}\\ {}&{{\delta _{{\rm{LVDA}}2}}}&{}&{}\\ {}&{}& \ddots &{}\\ {}&{}&{}&{{\delta _{{\rm{LVDA}}n}}} \end{array}} \right]({\rm{A}} = {\rm{F}},{\rm{v}}),\\ \mathit{\boldsymbol{G}} = \left[ {\begin{array}{*{20}{c}} {r_{\rm{D}}^{{D_{\rm{f}}} - d - \theta }}&{}&{}&{}\\ {}&{r_{\rm{D}}^{{D_{\rm{f}}}} - d - \theta }&{}&{}\\ {}&{}&{L_{\rm{D}}^2}&{} \end{array}} \right] \end{array} $

对矩阵方程进一步简化可得

$ {\mathit{\boldsymbol{K}}_{\rm{e}}}\mathit{\boldsymbol{P}}_{\rm{f}}^n = {\mathit{\boldsymbol{F}}_{\rm{e}}} $

(18)

其中K_e，F_e的表达形式为

$ \begin{gathered} {K_{\text{e}}} = \iiint\limits_{{\Omega _{\text{e}}}} {\left[ {\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{B}}^T}\mathit{\boldsymbol{GB}} + {\lambda _{{\text{mf}}}}{\mathit{\boldsymbol{H}}_{\text{f}}}\mathit{\boldsymbol{N}}{\mathit{\boldsymbol{N}}^T} + {\lambda _{{\text{vf}}}}\mathit{\boldsymbol{EN}}{\mathit{\boldsymbol{N}}^T} + } \right.} \hfill \\ \;\;\;\;\;\;\;\left. {{\omega _{\text{f}}}{{\left( {{L_{\text{D}}}{h_{\text{D}}}} \right)}^2}\frac{{\mathit{\boldsymbol{N}}{\mathit{\boldsymbol{N}}^T}}}{{\Delta {t_{\text{D}}}}}} \right]{\text{d}}V \hfill \\ \end{gathered} $

(19)

$ \begin{gathered} {\mathit{\boldsymbol{F}}_{\text{e}}} = \iiint\limits_{{\Omega _{\text{e}}}} {{\omega _{\text{f}}}{{\left( {{L_{\text{D}}}{h_{\text{D}}}} \right)}^2}\frac{{\mathit{\boldsymbol{N}}{\mathit{\boldsymbol{N}}^T}}}{{\Delta {t_{\text{D}}}}}{\text{d}}V\mathit{\boldsymbol{P}}_{\text{f}}^{n - 1}} + \hfill \\ \;\;\;\;\;\;\;\iiint\limits_{{\Omega _{\text{e}}}} {{\lambda _{{\text{mf}}}}{\mathit{\boldsymbol{H}}_{\text{f}}}\mathit{\boldsymbol{N}}{\mathit{\boldsymbol{N}}^T}{\text{d}}V\mathit{\boldsymbol{P}}_{\text{m}}^{n - 1}} + \hfill \\ \;\;\;\;\;\;\;\iiint\limits_{{\Omega _{\text{e}}}} {{\lambda _{{\text{vf}}}}\mathit{\boldsymbol{EN}}{\mathit{\boldsymbol{N}}^T}{\text{d}}V\mathit{\boldsymbol{P}}_{\text{v}}^{n - 1}} \hfill \\ \end{gathered} $

(20)

式中：K_e为非线性系数矩阵；F_e为单元载荷向量；上标n表示时刻变量。

式（18）为三重介质分形非线性模型裂缝系统的有限元单元平衡方程，同理可采用显式求解方法得到基质和溶洞系统的有限元单元平衡方程，且计算结果表明隐式求解裂缝系统、显式求解基质和溶洞系统的方法符合计算精度。

考虑到无限导流能力模型比均匀流量模型求解困难，通常选取0.7 L处为水平井在这2种模型下的等价压力点，从而可以借助均匀流量模型评价水平井井底压力。通过把水平井考虑成单元源汇项并积分，应用Delta函数将水平井划分成n个节点，并进行无因次化得到单元源汇项有限元方程。

所求压力解经过Laplace变换后，通过Duhamel原理引入井筒储集与表皮系数：

$ {\bar p_{{\text{wD}}}} = \frac{{s{{\bar p}_{\text{D}}} + \frac{S}{{{L_{\text{D}}}}}}}{{s\left[ {1 + {C_{\text{D}}}s\left( {s{{\bar p}_{\text{D}}} + \frac{S}{{{L_{\text{D}}}}}} \right)} \right]}} $

(21)

式中：S为表皮因数；s为Laplace变量；C_D为井筒储集系数。

采用Stehfest数值反演，得到考虑井筒储集效应与表皮因数的水平井无因次井底压力为

$ {p_{{\text{wD}}}} = \frac{{\ln 2}}{{{t_{\text{D}}}}}\sum\limits_{i = 1}^N {{V_i}} {\bar p_{{\text{wD}}}} $

(22)

$ {V_i} = {( - 1)^{\left( {\frac{N}{2} + i} \right)}}\sum\limits_{k = \frac{{i + 1}}{2}}^{\max \left[ {i,\frac{N}{2}} \right]} {\frac{{{k^{\frac{N}{2}}}(2k + 1)!}}{{(k + 1)!k!\left( {\frac{N}{2} - k + 1} \right)!(i - k + 1)!(2k - i + 1)!}}} $

(23)

3 渗流规律

图 2为碳酸盐岩油藏水平井的达西、非线性以及分形非线性模型的压力动态曲线。各基本参数取值如下：C_D= 100，S = 0，ω_f = 0.01，ω_v = 0.1，λ_mf= 0.001，λ_mv = 0.000 5，λ_vf = 0.1，L_D = 5，h_D = 200，α_KD = 0.02。假设非线性参数c_D= c_1D= c_2D，分形指数β = 2 + θ - D_f，则各模型其他参数的取值如下：达西模型c_D= 0，θ = 0，D_f = 2；非线性流模型c_D= 0.3，θ = 0，D_f = 2；分形非线性模型c_D = 0.5，θ = 0.3，D_f = 1.9。在考虑储层应力敏感性的基础上，可以根据实际情况适当调整各模型的α_KD值。

不同模型在不同流动阶段曲线的形态各异，但压力动态曲线仍然呈9个特征流动阶段（图 2）。以达西模型为例，具体流动阶段划分及不同流动阶段压力动态曲线特征^[23-24]如表 1所列。通过分析和对比非线性、分形非线性模型的压力动态曲线与达西模型的压力动态曲线的异同点（图 2）可发现，非线性参数主要影响基质-裂缝窜流强度，而分形指数使得溶洞-裂缝窜流后的流动阶段的压力动态曲线发生上翘。

4 敏感性分析

首先对非线性参数和分形指数进行敏感性分析，各基本参数取值如下：C_D = 100，S = 0，ω_f = 0.1，ω_v = 0.5，λ_m = 0.001，λ_v = 0.1，L_D = 5，h_D = 200，α_KD = 0.02。图 3与图 4分别为非线性参数c_D和分形指数β对压力动态曲线的影响。

非线性参数主要对基质向裂缝的窜流段有显著影响，基本不影响溶洞向裂缝系统的窜流段。随着非线性参数的增大，基质向裂缝的窜流强度减弱，持续时间变短。在压力导数曲线上表现为第2个下凹段延迟出现，下凹程度减弱，但窜流段结束的时间不受影响。

分形指数对压力动态曲线的影响主要体现在溶洞向裂缝窜流后的流动阶段，而对早期径向流和早期线性流阶段只产生微小的影响。从溶洞-裂缝拟径向流阶段开始，压力及其导数曲线随着时间的增加而逐渐上翘，且上翘幅度随着分形指数的增大而加剧。分形指数并不影响各个窜流段出现的时间以及窜流强度，所以压力导数曲线上下凹段的深度和持续时间均保持不变。

图 5与图 6分别为水平井长度L_D和纵向位置z_wD对压力动态曲线的影响，各基本参数取值如下：C_D = 100，c_D= 0，β = 0.1，S = 0，ω_f = 0.01，ω_v = 0.1，λ_m= 0.001，λ_v= 0.1，h_D= 200，α_KD= 0.02。研究表明，水平井越长，压力及其导数曲线越低，早期径向流越明显，从而后续流动阶段出现的时间延迟。当水平井长度小到一定程度时，早期径向流易被井筒储集效应所掩盖。考虑储层分形和应力敏感性特征，不同长度水平井的压力动态曲线并不在后期汇聚在一起，而是呈平行分布。在其他参数不变的情况下，水平井越靠近油层中部，即z_wD的值越接近0.5，早期径向流阶段出现得越早，同时早期线性流持续时间越长，但水平井纵向上位置的变化并不影响其他流动阶段。

5 应用实例分析

为了验证本文模型的正确性，进一步认识具有分形特征的缝洞型碳酸盐岩油藏的渗流特征，对四川CN油田S63井的渗流参数进行了解释。为了便于试井解释，假定裂缝垂向渗透率与水平渗透率相等，即β = 1。S63井的基础参数取值如表 2所列。

各个模型的异同已在渗流规律及敏感性分析部分进行了详细对比，但为方便对比、分析本文分形非线性模型的合理性，同时应用达西模型和本文模型进行试井拟合（图 7）。2个模型的拟合结果均显示压力导数与压力曲线分开。实际压力导数曲线和理论曲线在前期拟合程度相对较低，这是由于实际储层的井筒储集系数发生变化所导致的，在实际试井过程中普遍存在类似的问题。在渗流的中后期，虽然实测数据不完整，压力导数曲线仍明显出现反映孔洞缝三重介质储层渗流特征的2个下凹段，且相对于达西模型，本文模型的拟合结果更加精确，理论计算值和实测值拟合得更好。选取本文的分形非线性模型进行S63井的试井解释，表 3为该模型所解释的渗流参数。

实际生产曲线的基质向裂缝的窜流段明显比达西模型所呈现的下凹程度小，且试井解释下的非线性参数为0.15，基质渗透率为1.52 mD，这都体现了基质确实存在非达西渗流。拟合结果表明，采用本文模型进行渗流参数解释更加精确，这说明储层考虑分形特征和非线性渗流后更符合油井的实际情况，可应用本文的分形非线性模型进行实际试井解释及渗流规律分析。

6 结论

（1）引入分形理论描述裂缝系统的分形特征，同时应用非线性渗流新模型描述致密基质块的渗流特征，在此基础上建立缝洞型油藏水平井分形非线性渗流模型，并基于有限元方法求解水平井动态压力。

（2）通过对比达西模型与分形非线性模型的异同点，总结缝洞型碳酸盐岩油藏水平井的渗流规律，分析各个流动阶段的特点，将渗流过程划分为9个流动阶段，即井筒储存、表皮效应过渡、早期径向流、早期线性流、中期径向流、溶洞向裂缝窜流、溶洞-裂缝拟径向流、基质向裂缝窜流以及系统拟径向流。

（3）非线性参数主要影响基质向裂缝系统的窜流段，随着非线性参数的增大，基质向裂缝系统窜流强度减小，持续时间变短。分形指数不影响窜流出现的时间及强度，它对压力动态曲线的影响主要体现在溶洞向裂缝窜流后的阶段，压力及导数曲线随时间增加而逐渐上翘，且上翘幅度随分形指数的增大而加剧。

（4）利用本文的分形非线性模型进行实际井动态压力拟合并解释相关渗流参数，考虑了储层分形特征和非线性渗流后，模型更贴近矿场实际，可用其进行实际试井解释及渗流规律分析。