引言

水力压裂技术是页岩气藏效益开发的关键技术手段，在压裂过程中，裂缝的起裂扩展规律及导流特性是决定储层的改造体积及有效压裂周期的两大关键因素^[1-5]。受地应力、岩石强度以及天然裂缝发育程度等因素的影响，压裂会形成包括充填裂缝、剪切滑移缝和原位闭合缝在内的多种复杂裂缝^[6-7]。随着缝网复杂程度的增大，人工裂缝中支撑剂的平均支撑浓度小于0.5 kg/m$^2$或者支撑面积小于10%，大部分次生裂缝没有支撑剂充填^[8]，但由于地层的非均质性、偏差应力和射孔相位角的不匹配等因素，造成自支撑压裂裂缝壁面剪切滑移^[9]，而剪切裂缝面的凸起能够有效支撑裂缝^[10]，在闭合应力作用下仍然能保持残余缝宽，从而形成具有一定导流能力的自支撑裂缝。国外学者利用Barnett页岩露头进行了实验研究^[11]，结果显示，无支撑剂的滑移裂缝可以在28 MPa的地层闭合压力下提供15 mD·cm的导流能力，比未滑移的自支撑缝高1$\sim$2个数量级，能够为气体流动提供有效通道。同时，自支撑裂缝导流能力受杨氏模量和闭合压力共同影响^[9]，由于矿物组分和岩石力学性质的差异，导致滑移后的自支撑裂缝样品的导流能力存在较大的差异^[11]，准确评价和预测自支撑裂缝导流能力，对页岩气藏压裂后的产能预测有重要的影响。

1 导流能力预测模型

自支撑裂缝导流能力受裂缝滑移程度、表面粗糙度、岩石力学性质以及闭合应力等多种因素的共同影响。目前，对于裂缝导流能力的数学模型主要利用欧式几何参数建立，包括单一裂缝面几何参数和对应面叠合参数。前者的表征参数应用最广泛的是表面粗糙度^[12-14]，采用10条标准轮廓线来表征节理表面轮廓线的粗糙度，并得到国际岩石力学协会推荐使用，但是由于样本量过少，无法客观地进行定量对比，因而在使用中的误差较大^[15]。近年来，有学者提出采用曲折比^[16]、偏度和峰度^[17]等参数来表征裂缝表面的粗糙程度；此外，也有学者提出了采用吻合度^[18]、接触比^[19]和裂缝开度标准方差^[20]等裂缝面叠合参数来描述裂缝开度，进而建立裂缝导流能力数学模型。然而，这些统计学参数所包含的信息量有限，难以全面体现粗糙表面的随机行为和细节特征^[21]，且受测量仪器误差的影响相对较大，而分形理论具有解决复杂非线性问题的优势，分形维数可以克服测量仪器分辨率和取样长度等因素的影响，较为全面地反映裂缝面的形态特征^[22-25]。本文采用分形理论来表征裂缝面叠合后流动空间的分布特征，建立考虑闭合应力作用下的裂缝导流能力预测模型。

1.1 物理模型

自支撑裂缝由裂缝表面微凸体相互接触形成流动空间，将流体在裂缝内的流动进行简化处理(图 1)，假设流体只沿裂缝长度方向流动，且为达西流；裂缝在闭合应力作用下只发生弹性变形；将各个流线等效为大量毛细管束，并引入表征毛管直径大小分布的分形维数$D_{\rm{f}}$来定量分析毛细管的数量，同时引入表征毛细管迂曲程度的分形维数$D_{\rm{t}}$来定量分析毛细管的实际长度$L_{\rm{e}}$，依此即可定量表征自支撑裂缝的流动空间。

1.2 数学模型

多孔介质中的渗流遵循分形标度定律^[26]，则直径大于$r$的毛细管数量为

$ N(R > r) = {({r_{\max }}/r)^{{D_{\rm{f}}}}} $

(1)

式中：

$N(R > r$—直径大于$r$的毛细管数量，无因次；

$R$—毛细管直径的测量尺度，mm；

$r$—毛细管直径，mm；

$r_{\max }$—最大毛细管直径，mm；

$D_{\rm{f}}$—毛管直径大小分布分形维数，无因次。

毛细管实际长度$L_{\rm{e}}$大于裂缝的表观长度$L$，可由下式获得^[27]

$ {L_{\rm{e}}} = {r^{1 - {D_{\rm{t}}}}}{L^{{D_{\rm{t}}}}} $

(2)

$L_{\rm{e}}$—毛细管实际长度，mm；

$L$—裂缝的表观长度，mm；

$D_{\rm{t}}$—毛管迂曲程度分形维数，无因次；

$r_{\max }$与最大裂缝开度$\lambda_{\max }$呈正相关关系，即$r_{\max }$=$m\lambda_{\max }$，0 < $m\leqslant$1。

令$m = {({\lambda _{\max }}/r)^n}$，则有$- 1 < n \leqslant0$。

由此，可将式(1)改写为

$ N(R > r) = {({\lambda _{\max }}/r)^{{D_{\rm{f}}}{'}}} $

(3)

式中：

$D_{\rm{f}}{'}$—裂缝开度分形维数，${D_{\rm{f}}}{'} = {D_{\rm{f}}} + n$。

由于1 < $D_{\rm{f}}$ < 2，则$D_{\rm{f}}{'}$的取值在0$\sim$2；裂缝开度可通过对裂缝表面进行数字化处理获得^[18]。

假设裂缝空间由$N$个毛细管束(图 1)组成，则直径介于$r$和$r$+d$r$的毛细管数量

$ - {\rm{d}}N(r) = {D_{\rm{f}}}{'}{\lambda _{\max }}^{{D_{\rm{f}}}{'}}{r^{ - ({D_{\rm{f}}}{'} + 1)}}{\rm{d}}r $

(4)

因此，裂缝流动空间可以通过毛细管的数量和截面积来表示

$ V = \int_{{r_{\min }}}^{{r_{\max }}} {\dfrac{{{\rm{\pi }} {r^2}}}{4}{L_{\rm{e}}}} {\rm{d}}N\left( r \right) $

(5)

在闭合应力作用下，裂缝表面微凸体相互接触的变形机理复杂，会同时发生弹性变形和塑性变形甚至破碎，为简化数学模型，假设微凸体仅发生弹性变形。因此，毛细管直径$r$将随闭合应力的增大而线性减小，毛细管长度$L_{\rm{e}}$也将随闭合应力的增大而线性增大

$ r = \left( {1 - \sigma /E} \right){r_0} $

(6)

$ {L_{\rm{e}}} = \left( {1 + \sigma \upsilon /E} \right)L $

(7)

式中：

$\sigma$—闭合应力，MPa；

$E$—杨氏模量，GPa；

$r_0$—未加载闭合应力时的毛细管直径，mm；

$\upsilon$—泊松比，无因次。

将式(4)、式(6)和式(7)代入式(5)，即可获得随闭合应力变化的裂缝空间体积

$ V = \dfrac{{{\rm{\pi }} {L^{{D_{\rm{t}}}}}{D_{\rm{f}}}{'}\lambda _{\max 0}^{3 - {D_{\rm{t}}}}{{\left( {1 - \sigma /E} \right)}^{3 - {D_{\rm{t}}}}}}}{{4\left( {3 - {D_{\rm{f}}}{'} - {D_{\rm{t}}}} \right)}}\left[ {1 - {{\left( {\dfrac{{{r_{\min 0}}}}{{{\lambda _{\max 0}}}}} \right)}^{3 - {D_{\rm{f}}}{'} - {D_{\rm{t}}}}}} \right] $

(8)

式中：

$V$—裂缝空间体积，mm$^3$；

$r_{\min 0}$—未加载闭合应力时毛细管最小直径，mm；

$\lambda _{\max 0}$—未加载闭合应力时裂缝最大开度，mm。

在初始条件下，即不考虑闭合应力时的初始裂缝空间体积为

$ {V_0} = \dfrac{{{\rm{\pi }} {L^{{D_{\rm{t}}}}}{D_{\rm{f}}}{'}\lambda _{\max 0}^{3 - {D_{\rm{t}}}}}}{{4\left( {3 - {D_{\rm{f}}}{'} - {D_{\rm{t}}}} \right)}}\left[ {1 - {{\left( {\dfrac{{{r_{\min 0}}}}{{{\lambda _{\max 0}}}}} \right)}^{3 - {D_{\rm{f}}}{'} - {D_{\rm{t}}}}}} \right] $

(9)

流体在单个毛细管中的流动符合Navier-Stokes方程^[26]

$ q = \dfrac{{{\rm{\pi }} {r^4}}}{{128\mu }}\dfrac{{{\rm{d}}p}}{{{\rm{d}}{L_{\rm{e}}}}} $

(10)

式中：

$q$—单个毛细管的流量，m$^3$/s；

$p$—流动压力，MPa；

$\mu$—流体黏度，mPa·s。

整个裂缝的流量即为$N$个毛细管流量的总和

$ Q = \int_{{r_{\min }}}^{{r_{\max }}} q {\rm{d}}N\left( r \right) $

(11)

将式(4)和式(10)代入式(11)，即可得到随闭合应力变化的裂缝流量

$ Q = \dfrac{{{\rm{\pi }} \Delta p{D_{\rm{f}}}{'}\lambda _{\max 0}^{3 + {D_{\rm{t}}}}{{(1 - \sigma /E)}^{3 + {D_{\rm{t}}}}}}}{{128\mu {L^{{D_{\rm{t}}}}}\left( {3 + {D_{\rm{t}}} - {D_{\rm{f}}}{'}} \right)\left( {1 + \sigma \upsilon/E} \right)}}\left[ {1 - {{\left( {\dfrac{{{r_{\min 0}}}}{{{\lambda _{\max 0}}}}} \right)}^{3 + {D_{\rm{t}}} - {D_{\rm{f}}}{'}}}} \right] $

(12)

在初始条件下，即不考虑闭合应力时的流量为

$ {Q_0}\!=\!\dfrac{{{\rm{\pi }} \Delta p{D_{\rm{f}}}{'}\lambda _{\max 0}^{3 + {D_{\rm{t}}}}}}{{128\mu {L^{{D_{\rm{t}}}}}\left( {3\!+\!{D_{\rm{t}}} - {D_{\rm{f}}}{'}} \right)}}\left[ {1\!-\!{{\left( {\dfrac{{{r_{\min 0}}}}{{{\lambda _{\max 0}}}}} \right)}^{3 + {D_{\rm{t}}} - {D_{\rm{f}}}{'}}}} \right] $

(13)

式中：

$\Delta p$—流体入口端与出口端的压差，MPa；

$Q_0$—不考虑闭合应力时的流量，m$^3$/s。

通过数字化手段将两个对应裂缝面进行叠加即可获得初始裂缝流动空间$V_0$，通过导流能力实验可获得初始裂缝流量$Q_0$，在此基础上，联立上述两个参数的表达式，即可计算得到裂缝开度分形维数${{D_{\rm{f}}}{'}}$和迂曲度分形维数$D_{\rm{t}}$。

裂缝导流能力为裂缝渗透率和裂缝开度的乘积，其中，裂缝渗透率$K_{\rm{f}}$可由达西定律给出

$ {K_{\rm{f}}} = \dfrac{{Q\mu L}}{{A\Delta p}} $

(14)

式中：

$K_{\rm{f}}$—裂缝渗透率，mD；

$A$—裂缝截面积，mm$^2$，可由所有毛细管的总面积进行表征。

另外，裂缝开度$w$可由毛细管截面积获得

$ w = \dfrac{A}{{N(r)}} $

(15)

因此，自支撑裂缝导流能力为

$ {C_{\rm{f}}} = K_{\rm{f}}w = \dfrac{{Q\mu L}}{{\Delta pN(r)}} $

(16)

式中：$C_{\rm{f}}$—裂缝导流能力，D·cm。

将式(3)和式(12)代入式(16)，得到了利用裂缝开度分形维数和迂曲度分形维数进行表征，同时考虑闭合应力作用的自支撑裂缝导流能力数学模型

$ {C_{\rm{f}}}\!\!=\!\!\dfrac{{{\rm{\pi }} \!{D_{\rm{f}}}{'}r_{\min 0}^{3\!+\! {D_{\rm{t}}}}{{(1\!-\!\sigma /E)}^{3\!+\!{D_{\rm{t}}}}}{L^{1\!-\! {D_{\rm{t}}}}}}}{{128\left(\!{3\!+\!{D_{\rm{t}}}\!-\!{D_{\rm{f}}}{'}} \!\right)\left(\!{1\!+\!\sigma v/E}\!\right)}}\left[ {{{\left(\!{\dfrac{{\lambda _{\max\!0}^{}}}{{{r_{\min\!0}}}}}\!\right)}^{3\!+\!{D_{\rm{t}}}\!-\! {D_{\rm{f}}}{'}}}\!-\!1} \right] $

(17)

上述数学模型从流动机理出发，引入两个分形维数(裂缝开度分形维数和迂曲度分形维数)来表征裂缝空间，同时基于弹性力学理论，考虑了裂缝空间在闭合应力作用下的变化过程，从而可以获得不同闭合应力下裂缝导流能力，为验证理论模型预测结果的准确性，开展了不同闭合应力下的裂缝导流能力实验。

2 实验测试 2.1 样品准备

实验所用岩样取自威远地区龙马溪组一段3号小层，杨氏模量为23.5 GPa，泊松比为0.15，密度2.57 g/cm$^3$。由于页岩含有大量极小弱面，难以有效劈裂岩板，所以采用柱塞岩芯进行剖缝后测定导流能力。制备3块1寸标准岩芯并开展基础物性测试，测试平均孔隙度4.28%，通过巴西劈裂法获取人工裂缝(图 2)。

2.2 导流能力测试

参照Freed^[28]等人的自支撑裂缝导流能力测试实验，确定实验中的滑移量为2.5 mm，通过在相反的断面处贴铜箔垫片的方法来产生错位。

实验采用精密气驱实验装置(图 3)，选择氮气为驱替介质，用造泡流量计测量气体流量。具体实验步骤为：(1)将打磨好的柱塞岩芯放入岩芯夹持器中，根据设定的滑移量放入金属垫片；(2)通入氮气，待流体通过稳定后，记录出口端压力，和对应时间段的流量；(3)改变围压至1 MPa，待流体通过稳定后，再次记录出口端压力和对应的流量；(4)继续依次改变围压(5$\sim$40 MPa)，分别记录各围压下的出口端压力和对应的流量。计算出柱塞岩芯气测导流能力大小，实验结果如图 3所示。

实验结果显示，裂缝滑移后的自支撑裂缝具备一定的导流能力，初期可达30$\sim$50 D·cm；3个自支撑裂缝样品导流能力的变化规律基本相同，在闭合应力增大到10$\sim$20 MPa之前，导流能力下降速度较快，随着压力的进一步增大，导流能力变化相对平缓；由于3个样品取自同一块大岩芯，其岩石力学性质和物性均相同，分析认为裂缝导流能力数值上的差异主要受裂缝面粗糙程度的影响，3#样品的裂缝滑移后微凸体接触产生的叠加高度最大，从而形成较大的缝宽，造成裂缝导流能力的整体增大。

3 结果讨论

利用本文建立的裂缝导流能力数学模型，对上述3个样品的裂缝导流能力进行理论预测，并将其与实验结果进行对比(图 5)，结果显示，模型计算的预测曲线变化趋势与实测结果基本一致，且数值基本吻合。通过模型计算得到1#、2#和3#共3个样品的裂缝开度分形维${{D_{\rm{f}}}{'}}$分别为0.80、0.90和0.95，迂曲度分形维数$D_{\rm{t}}$分别为1.8、1.6和1.5。由此可以看出，由于裂缝面形态的差异，造成3个样品裂缝闭合后的裂缝开度分形维数${{D_{\rm{f}}}{'}}$和迂曲度分形维数$D_{\rm{t}}$不同。

从前述物理模型中可知，${{D_{\rm{f}}}{'}}$表征裂缝开度的分布，$D_{\rm{t}}$表征裂缝流动方向的迂曲程度，通过计算不同分形维数对应的导流能力分析其敏感性。

裂缝开度分形维数${{D_{\rm{f}}}{'}}$取值在0$\sim$2，分别计算${{D_{\rm{f}}}{'}}$为0.4、0.6、0.8、1.0和1.2时对应裂缝导流能力(图 5)，结果显示，裂缝开度分形维数${{D_{\rm{f}}}{'}}$与导流能力呈正相关关系；另外，随着分形维数的逐渐增大，不同闭合应力下裂缝导流能力整体增大，但增大的幅度逐渐减小。分析认为，这是由于较大的裂缝开度分形维数表征裂缝开度越大且相对更加均匀，另外，模型中假设裂缝内流动为达西流，没有考虑毛管力以及闭合应力增大过程中毛管数量的减少。

迂曲度分形维数$D_{\rm{t}}$取值在1.0$\sim$2.0，分别计算$D_{\rm{t}}$为1.1、1.3、1.5、1.7和1.9时对应裂缝导流能力(图 7)，结果显示，迂曲度分形维数$D_{\rm{t}}$与导流能力呈负相关关系，裂缝导流能力随着迂曲度分形维数的增大而减小，且在对应闭合压力下的裂缝导流能力下降幅度相同，分析认为，这是由于随着迂曲度分形维数的增大，裂缝流动方向的迂曲程度增大，造成渗流阻力逐渐增大。

根据上述敏感性分析结果可以对实验测试结果进行分析。1#样品的裂缝开度分形维数为0.8，表明其裂缝开度整体小于另外两块样品；迂曲度分形维数为1.8，即其裂缝空间迂曲程度均高于另外两块样品，表明流体的渗流阻力相对较大，两种因素共同作用导致其导流能力整体偏低。对于3#样品而言，其裂缝开度分形维数为0.95，迂曲度分形维数为1.5，表明其裂缝开度整体大于前两个样品，同时其裂缝空间的迂曲程度也最小，因而，在相同的闭合应力下，其裂缝导流能力均高于另外两个样品。

实验结果与模型预测结果的对比分析表明，将裂缝流动通道用具有分形特征的毛细管束来描述的方法是准确可行的，通过计算裂缝开度分形维数${{D_{\rm{f}}}{'}}$和迂曲度分形维数$D_{\rm{t}}$可以有效预测闭合应力作用下自支撑裂缝的导流能力，基于分形理论建立的自支撑裂缝导流能力预测模型可大幅减少室内实验的工作量，包括岩芯样品的制备和实验耗时；另外，采用分形方法进行预测，其结果具有尺度不变性，可以消除室内实验结果与大尺度人工裂缝之间的尺度效应。虽然上述预测结果存在一定的误差，但其变化趋势及最终的长期导流能力是吻合的，根据某一区域或储层岩石力学特征，可建立相应的经验参数，从而使预测结果更为准确。

4 结论

(1) 利用迂曲度分形维数$D_{\rm{t}}$和裂缝开度分形维数${{D_{\rm{f}}}{'}}$来描述闭合应力下的自支撑裂缝空间，建立了考虑闭合应力作用的自支撑裂缝导流能力预测数学模型。模型预测结果与实验结果基本吻合，预测方法准确可靠，研究区自支撑裂缝的裂缝开度分形维数在0.9左右，迂曲度分形维数在1.5$\sim$1.8。

(2) ${{D_{\rm{f}}}{'}}$表征裂缝开度的分布，与导流能力呈正相关关系随着分形维数的逐渐增大，不同闭合应力下裂缝导流能力整体增大，但增大的幅度逐渐减小；$D_{\rm{t}}$表征裂缝流动方向的迂曲程度，与导流能力呈负相关关系，随$D_{\rm{t}}$数值的变化，裂缝导流能力的变化程度基本相当。

(3) 基于分形理论建立的自支撑裂缝导流能力预测模型可大幅减少室内实验的工作量，同时预测结果具有尺度不变性，可有效表征水力压裂形成的自支撑裂缝渗流特征，为气井的产能评价及预测提供支撑。