引言

水力压裂施工过程中，在裂缝净压力作用下，压裂液向地层滤失，压裂液的滤失对裂缝几何尺寸、支撑剂的分布都有较大的影响，特别是在裂缝性储层中，由于天然裂缝的存在，压裂液的滤失量会更大。

对于裂缝性储层压裂液滤失，国内外做了大量的研究，Penny G等^[1]对比了静态和动态条件下压裂液的滤失速度变化，对基岩和裂缝中的滤失控制机理进行了实验评价，指出裂缝对压裂液滤失起控制作用；Barree R D等^[2]提出一种随压力变化的压裂液滤失模型，分析对天然裂缝储层进行压裂后，压裂液滤失对裂缝延伸的影响。文献[3]提出了一种考虑多参数的压裂液在双重介质中的滤失模型，考虑了影响滤失量的多种实验参数，压裂液在实际地层中滤失的3个区域，并推导了压裂液在油藏区为双重介质渗流机理控制下滤失速度方程。李勇明等考虑基岩和天然裂缝对压裂液滤失的影响，建立了裂缝性低渗透储层压裂液滤失模型^[4-6]。以上对裂缝性储层压裂液滤失模拟研究均对天然裂缝进行了简化处理，未考虑天然裂缝实际分布情况，不能准确描述裂缝性储层压裂液滤失特征。

描述裂缝的参数主要有裂缝的位置、长度、方位、开度、密度等^[7-8]。对油气藏内所有的裂缝都进行测量和模拟是不可能的，一般通过测井、录井、地震等方法获得裂缝分布数据，利用概率统计方法得出裂缝几何参数概率分布密度函数，最后运用Monte-Carlo法生成每条裂缝的几何参数：中心点位置、产状、半径和开度^[9-11]。

在参考上述方法的基础上，建立储层二维离散裂缝网络模型，借助等效渗透率张量原理将裂缝性介质转化为均质各向异性连续介质，建立了基于渗透率张量的裂缝性储层压裂液滤失数学模型。

1 裂缝性地层压裂液滤失的数学模型 1.1 裂缝性介质整体渗透率等效原理

在相同压力梯度和外边界条件下，若同种流体通过相同尺寸的非均质裂缝型介质体和均质各向异性连续介质体的流量相等，则裂缝型介质体可由均质各向异性介质体等效替代，裂缝型介质体的整体渗透率可由均质各向异性介质体的渗透率表示，此即为流量等效原理。对所研究的非均质裂缝型介质区域(图 1a，Γ₁，Γ₂，Γ₃，Γ₄-网格边界；n₁，n₂，n₃，n₄-网格边界Γ₁，Γ₂，Γ₃，Γ₄的单位外法向量)进行网格划分，计算每个网格块的等效渗透率张量，用等效渗透率张量表征的均质各向异性连续介质体代替原裂缝介质体，形成新的等效研究区域(图 1b)，即可利用全张量渗透率连续介质渗流理论分析非均质体渗流。

1.2 裂缝性储层等效渗透率张量的计算

Lee S H等提出对油藏中存在的多级规模的裂缝进行分级处理^[12-13]，根据裂缝长度l_f与离散网格单元尺寸l_g的比值，将裂缝划分为三类：短裂缝$({l_{\rm{f}}}/{l_{\rm{g}}} \ll 1)$，中裂缝${l_{\rm{f}}}/{l_{\rm{g}}} \approx 1)$和长裂缝$({l_{\rm{f}}}/{l_{\rm{g}}} \gg 1)$，利用网格边界对裂缝进行截断处理。截断处理之后，网格内的裂缝包括两种情况：裂缝存在于网格内和裂缝跨网格分布。分别采用不同的边界条件，对分布在各自网格的裂缝进行处理。此外还假设：裂缝为具有一定开度的两平行板构成；基岩为均质各向同性多孔介质；渗流属于单相稳定渗流；不考虑流体压缩性和多孔介质变形对渗流的影响。

基岩和裂缝中的渗流满足达西定律和质量守恒定律

$ {v_{\rm{m}}} = - \frac{{{K_{\rm{m}}}}}{\mu }\nabla {p_{\rm{m}}} $

(1)

$ {v_{{\rm{f}},i}} = - \frac{{{K_{{\rm{f}},i}}}}{\mu }\nabla {p_{{\rm{f}},i}} $

(2)

$ - \nabla {v_{\rm{m}}} + {q_{{\rm{fm}}}} = 0 $

(3)

$ - \nabla {v_{{\rm{f}},i}} + {q_{{\rm{fm}},i}} = 0 $

(4)

式中：

$v_{\rm{m}}$-基岩中的渗流速度，cm/s；

$K_{\rm{m}}$-基岩的渗透率，D；

$\mu$-黏度，mPa.s；

$p_{\rm{m}}$-基岩的孔隙压力，MPa；

$v_{{\rm{f}}, i}$-第i条裂缝的渗流速度，cm/s；

$K_{{\rm{f}}, i}$-第i条裂缝的渗透率，D；

$p_{{\rm{f}}, i}$-第i条裂缝与基岩接触面处的裂缝压力，MPa；

$q_{\rm{fm}}$-基岩和所有裂缝之间的体积流量，cm³/s；

$q_{{\rm{fm}}, i}$-第i条裂缝与周围基岩之间的体积流量，cm³/s。

将运动方程代入连续性方程，得到描述基岩和网格块中第i条裂缝中流体渗流微分方程

$ \frac{{{K_{\rm{m}}}}}{\mu }\frac{{{\partial ^2}{p_{\rm{m}}}}}{{\partial {x^2}}} + \frac{{{K_{\rm{m}}}}}{\mu }\frac{{{\partial ^2}{p_{\rm{m}}}}}{{\partial {y^2}}} + {q_{{\rm{fm}}}} = 0 $

(5)

$ \frac{{{K_{{\rm{f}},i}}}}{\mu }\frac{{{\partial ^2}{p_{{\rm{f}},i}}}}{{\partial {x^2}}} + \frac{{{K_{{\rm{f}},i}}}}{\mu }\frac{{{\partial ^2}{p_{{\rm{f}},i}}}}{{\partial {y^2}}} + {q_{{\rm{fm}},i}} = 0 $

(6)

接触边界条件

$ {p_{{\rm{f}},i}}{|_\mathit\Omega } = {p_{\rm{m}}}{|_\mathit\Omega } $

(7)

$ ({v_{{\rm{f}},i}}{|_\mathit\Omega }) \cdot \boldsymbol{n} = ({v_{\rm{m}}}{|_\mathit\Omega }) \cdot \boldsymbol{n} $

(8)

式中：

$\varOmega$-裂缝与基岩接触面；

$\boldsymbol{n}$-边界的外法向单位向量。

在含有裂缝的正方形基岩网格块中(图 1)，网格块的周期边界(纵向Γ₁，Γ₃；横向Γ₂，Γ₄条件可表示为^[14]

$ \left\{ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {p(x,y = 0) = p(x,y = {l_y}) - {l_y}\frac{{\partial p}}{{\partial y}}}\\ {v(x,y = 0) \cdot {\boldsymbol{n}_{\rm{1}}} = - v(x,y = {l_y}) \cdot {\boldsymbol{n}_{\rm{3}}}} \end{array}}\\ {\begin{array}{*{20}{c}} {p(x = 0,y) = p(x = {l_x},y) - {l_x}\frac{{\partial p}}{{\partial x}}}\\ {v(x = 0,y) \cdot {\boldsymbol{n}_2} = - v(x = {l_x},y) \cdot {\boldsymbol{n}_4}} \end{array}} \end{array}} \right. $

(9)

式中：p-压力，MPa；l_x，l_y-网格x向、y向长度，cm。

通过边界元法对求解区域进行离散，可降低求解的难度，可提高计算精度。因此，可以采用边界元法对上述方程组进行求解。

基岩和裂缝中的渗流方程在本质上来说是Laplace方程

$ K{\nabla ^2}{p_\varsigma } = 0 $

(10)

式中：

K-渗透率，D；

$p_\varsigma ^*$-边界节点处裂缝压力或者基岩压力，MPa。

利用第二格林公式和δ函数的性质，通过极限过程推导出求解区域边界上各点的边界积分方程，其一般的形式为

$ {c^i}{p_\varsigma }^i + \int\limits_\varsigma {v_\varsigma ^*{p_\varsigma }{\rm{d}}\mathit\Gamma } = \int\limits_\varsigma {{v_\varsigma }p_\varsigma ^*{\rm{d}}\mathit\Gamma } $

(11)

其中

$ {c^i} = \frac{{{\alpha _i}}}{{2{\rm{\pi }}}} $

(12)

$ p_\varsigma ^{\rm{*}} = \frac{1}{{2{\rm{\pi }}K}}\ln \frac{{\sqrt K }}{r} $

(13)

$ v_\varsigma ^{\rm{*}} = K\frac{{\partial p_\varsigma ^{\rm{*}}}}{{\partial {n_i}}} $

(14)

$ r = \sqrt {{{\left( {x - {x_{\rm{i}}}} \right)}^2} + {{\left( {y - {y_{\rm{i}}}} \right)}^2}} $

(15)

式中：

c-系数，无因次；

$p_\varsigma ^*$-渗流方程的Laplace基本解，MPa；

${v_\varsigma }$-边界节点处的外法向裂缝渗流速度或基岩渗流速度，cm/s；

$v_\varsigma ^*$-Laplace方程的基本解的外法向导数与渗透率的乘积，cm/s；

$\varGamma $-网格块外边界或基岩与裂缝的交界面。

为了计算方便，采用常单元对求解区域进行剖分，将边界Γ分为N个小单元，以各小单元的中点为节点。对边界积分方程进行分段积分，可得

$ {c^i}{p^i} + \sum\limits_{j = 1}^N {\int\limits_{{\mathit\Gamma _{\rm{j}}}} {{v^*}p{\rm{d}}\mathit\Gamma } } = \sum\limits_{j = 1}^N {\int\limits_{{\mathit\Gamma _{\rm{j}}}} {{p^*}v{\rm{d}}\mathit\Gamma } } $

(16)

根据常单元的特性，以每一个单元上的节点值代替单元上的值，即在每一个单元上p_j和v_j是常数，则有

$ {c^i}{p^i} + \sum\limits_{j = 1}^N {\left( {\int\limits_{{\mathit\Gamma _{\rm{j}}}} {{v^*}{\rm{d}}\mathit\Gamma } } \right)} {p_j} = \sum\limits_{j = 1}^N {\left( {\int\limits_{{\mathit\Gamma _{\rm{j}}}} {{p^*}{\rm{d}}\mathit\Gamma } } \right){v_j}} $

(17)

p^*和v^*是基本解及其法向导数，它们在每一个单元上的积分是已知的。式(17)可化为

$ \begin{array}{*{20}{c}} {\sum\limits_{j = 1}^N {{H_{i,j}}{p_j}} = \sum\limits_{j = 1}^N {{G_{i,j}}{v_j}} } \end{array} $

(18)

式中${H_{i,j}} = \int\limits_{{\mathit\Gamma _j}} {{v^*}{\rm{d}}\mathit\Gamma } + {c^i}{\delta _{i,j}}$；δi, j-Kronecker函数，当i=j时取1，i≠j时取0；${G_{i,j}} = \int\limits_{{\mathit\Gamma _j}} {{p^*}{\rm{d}}\mathit\Gamma } $。

全部节点离散后得到方程组可表示为

$ \boldsymbol{Hp} = \boldsymbol{Gv} $

(19)

式中：

$ \mathit{\pmb{H}}=\left[\begin{matrix} {{H}_{11}} & {{H}_{12}} & \cdots & {{H}_{1N}} \\ {{H}_{21}} & {{H}_{22}} & \cdots & {{H}_{2N}} \\ \cdots & \cdots & \cdots & \cdots \\ {{H}_{N1}} & {{H}_{N2}} & \cdots & {{H}_{NN}} \\ \end{matrix} \right]； $

$ \boldsymbol{G}=\left[\begin{matrix} {{H}_{11}} & {{H}_{12}} & \cdots & {{H}_{1N}} \\ {{H}_{21}} & {{H}_{22}} & \cdots & {{H}_{2N}} \\ \cdots & \cdots & \cdots & \cdots \\ {{H}_{N1}} & {{H}_{N2}} & \cdots & {{H}_{NN}} \\ \end{matrix} \right]； $

$ \boldsymbol{p}={{\left[\begin{matrix} {{p}_{1}} & {{p}_{2}} & \cdots & {{p}_{N}} \\ \end{matrix} \right]}^{\text{T}}}={{\left[\begin{matrix} {{v}_{1}} & {{v}_{2}} & \cdots & {{v}_{N}} \\ \end{matrix} \right]}^{\text{T}}}。 $

求解全部节点离散后得到方程组，通过等效渗透率张量原理计算，得到裂缝型介质的等效渗透率张量。

1.3 压裂液在地层中渗流模型

利用前面计算的渗透率全张量表征储层的渗透性，假设压裂液微可压缩，渗流符合达西定律，并忽略重力影响，二维单相渗流的广义达西公式和连续性方程为

$ v=-\frac{K}{{{\mu }_{\text{e}}}}\cdot \nabla p $

(20)

$ n-\nabla (\rho v)=\frac{\partial (\rho \phi )}{\partial t} $

(21)

式中：

v-渗流速度，cm/s；

$\rho$-密度，g/cm³；

$\phi$-孔隙度，%；

t-时间，s。

渗流速度的分量形式

$ n-\frac{1}{B}\left( \frac{\partial {{v}_{x}}}{\partial x}+\frac{\partial {{v}_{y}}}{\partial x} \right)={{C}_{\text{l}}}\frac{\partial p}{\partial t} $

(22)

式中：

$v_x$，$v_y$-x向，y向渗流速度，cm/s；

${\mu _{\rm{e}}}$-压裂液的有效黏度，mPa.s；

B-压裂液体积系数，无因次；

$C_{\rm{l}}$-压裂液的综合压缩系数，MPa^-1。

封闭外边界条件

$ \left\{ {\begin{array}{*{20}{l}} {{{\left. {\frac{{\partial p}}{{\partial y}}} \right|}_{y = {W_y}}} = 0}\\ {{{\left. {\frac{{\partial p}}{{\partial x}}} \right|}_{x = {W_x}}} = 0}\\ {{{\left. {\frac{{\partial p}}{{\partial x}}} \right|}_{x = 0}} = 0}\\ {{{\left. {\frac{{\partial p}}{{\partial y}}} \right|}_{({L_{\rm{f}}} < x{W_x},y = 0)}} = 0} \end{array}} \right. $

(23)

式中：

$W_x$、$W_y$-研究区域x方向、y方向的长度，m；

$L_{\rm{f}}$-裂缝长度，m。

内边界条件和初始条件分别为

$ {\left. p \right|_{(x{L_{\rm{f}}},y = 0)}} = {p_{\rm{f}}} $

(24)

$ {\left. p \right|_{t = 0}} = {p_{\rm{i}}} $

(25)

式中：

$p_{\rm{f}}$-人工裂缝压力，MPa；

$p_{\rm{i}}$-初始压力，MPa。

2 实例分析

在200 m×200 m的区域内，生成符合如下分布规律的离散裂缝网络：裂缝位置均匀分布；裂缝倾角符合正态分布，均值为60°，标准差20°；裂缝长度符合正态分布，均值为1.5 m，标准差为0.4 m；裂缝开度符合正态分布，均值为90 µm，标准差20 µm。结合储层基本参数(表 1)，进行模拟计算。

根据上述数据生成储层离散裂缝网络模型，并划分网格计算每个网格的等效渗透率分布，如图 2所示，其中用红线所圈部分为人工裂缝与天然裂缝的相交位置。将网格的裂缝分布情况与等效后的渗透率张量进行对比，可知基岩渗透率、裂缝密度、裂缝开度和裂缝方位等参数对裂缝性介质的渗透率张量都有影响，长裂缝对相应网格的等效渗透率张量起主要控制作用，短裂缝对等效渗透率张量的影响相对较小，不含裂缝网格块的等效渗透率张量就是基岩的渗透率。

在计算渗透率张量的基础上，模拟了压裂施工10 min和60 min时，水力裂缝一侧压力分布，如图 3所示。从图 3中可以看出由于天然裂缝分布的非均匀性，造成压力传播的不均匀。

图 4为压裂施工10 min、60 min时水力裂缝一侧压裂液滤失速度分布，由于水力裂缝一侧天然裂缝分布的不均匀性以及与水力裂缝的位置关系不同，使沿裂缝长度方向压裂液的滤失速度变化较大，

尤其当某处天然裂缝与水力裂缝连通时，该处的压裂液滤失速度明显增大%（图 2中标注了天然裂缝与水力裂缝相交的位置），对比发现这几处的压裂液滤失速度是其他位置的3~4倍。进而研究了这几处在不同时间的滤矢量与总滤失量之间的关系，发现这几处的滤失量占总滤失量的45%以上，如图 5所示。所以在对裂缝性储层进行压裂施工时，要对压裂液的性能进行调整，以适应裂缝性储层的要求。

随着裂缝净压力的增加，裂缝不同位置处的滤失速度均有所增加，但是不同位置处滤失速度增加的程度不同(图 6)。天然裂缝与人工裂缝连通处的滤失速度的增加幅度远大于其他各处，这是由于：(1)裂缝净压力的增加，使滤失压差增加，滤失量增加；(2)裂缝净压力的增加，使天然裂缝开度增加，从而增加了天然裂缝的渗流能力，使滤失速度增加。

3 结语

(1)提出了基于离散裂缝网络的裂缝性储层压裂液滤失模型，综合考虑了油藏边界、压裂液的非牛顿特性以及储层的天然裂缝分布，计算结果更符合实际。

(2)天然裂缝分布及其与水力裂缝的位置关系的不同，使得沿裂缝长度方向，滤失速度变化较大，若天然裂缝与水力裂缝连通，则压裂液滤失主要发生在连通处，当天然裂缝和人工裂缝存在三处连通时，这三处的滤失量占裂缝总滤失量的45%以上。所以，在裂缝性储层压裂施工中，需考虑压裂液的降滤性以及对天然裂缝的封堵性能。

(3)由于该模型在很大程度上依赖于储层复杂的天然裂缝分布，所以在如何获取更为准确的天然裂缝分布参数方面，有待进一步研究。