水力压裂(hydraulic fracture，又称为水力劈裂、水力裂解技术)是开采页岩气时所用的方法，也是油气井增产的主要措施^[1].水力压裂的数值计算方法及模型有多种，如弥散裂缝模型^[2]、双孔隙度模型^[3]、分离裂缝模型^[3]、粘结模型^{[4, 5]}、扩展有限元法^{[6, 7]}等.基于断裂力学理论的模型有PKN，PKN-C，KGD，KGD-C，P-3D和P-3D-C等^{[1, 8]}.弥散裂缝模型是通过修改破坏区域的材料特性(如：刚度、渗透系数等)实现裂缝发展过程的模拟，无需引入真正的裂缝和重构网格，适用于岩石材料的拉伸破坏和土体材料的剪切破坏过程的模拟^[2,3].但当采用此法进行数值计算时，Bazant^[9,10]等指出网格的形状会影响计算结果，即网格敏感性.本文通过对比不同本构关系和有限元网格的单元高宽比(AR)对水力压裂问题中垂直裂缝开裂区域的数值计算结果，分析并探讨了如何回避网格敏感性问题.

本文所采用的弥散裂缝模型是基于完全流-固耦合的弹塑性理论，其主要控制方程如下：

ABAQUS软件中，渗透系数又称水力传导系数(hydraulic conductivity)，在各向同性介质中，它被定义为单位水力梯度下的单位流量，表示流体通过孔隙骨架的难易程度，其单位为m/s.本文通过修改渗透系数的数值模拟破坏岩石中流体的流动.结合文献^{[11, 12]}中关于渗透系数的修改法，本文所采取的渗透系数k的修改方式是通过式(7)实现的，其中系数h表示有效平均应力σ′_m的双曲正切函数如式(8)所示，材料的抗拉强度P_t的单位是MPa.

图 1(Figure 1)

Fig. 1

图1 不同ξ取值下的渗透系数变化曲线 Fig.1 Modification of permeability with different ξ

水力压裂的数值分析软件的国内代表有基于渗流应力耦合法的RFPA^2D软件^[14].而本文采用ABAQUS软件的用户自定义子程序实现水力压裂的数值计算功能.在该软件中材料参数可以设置为场变量(field variable)函数.在求解之前可以通过子程序SDVINI设计材料分布情况并初始化状态变量(state variable)，而后通过子程序USDFLD在计算过程中根据状态变量更新场变量的数值实现材料参数随着积分点位置变化而变化.程序中共考虑8个场变量，依次为：弹性模量、泊松比、抗拉程序、渗透系数、密度、孔隙比、饱和度和混合黏滞系数，用于描述材料特性和渗透系数的计算.

文献^[3]指出，在弥散裂缝模型中存在网格敏感性和收敛问题.这个问题产生主要是由不准确的破坏模型的零能耗引起的^[9].目前的主要解决方法是添加一个位置限制器^[15]，常用的有：非局部连续体模型^[9]、裂缝带模型^{[10, 16]}、特征值长度法^{[17, 18]}等.

考虑本文所使用的方法是在求解过程中通过提高岩石材料的渗透系数实现流体在破坏岩石中的流动，属于一种弥散裂缝模型，因此同样存在网格敏感性问题.在不增加数值计算量的条件下，本文探讨单元高宽比(AR)与本构关系的选取对计算结果的影响，其中材料的弹性部分采用线弹性(LE)和多孔弹性(PE)，材料的屈服和破坏采用Drucker-Prager (DP)和Mohr-Coulomb (MC)两种模型.通过三组算例进行验证：1) PE-DP；2) LE-DP；3) LE-MC.每种组合方式的AR取值范围是0.1~19.

为了简化计算，算例中的材料参数均为常量，各参数取值如表1所示.

表 1

Table 1

表 1(Table 1)

表1 材料参数 Table 1 Material parameters 参数名称 取值 参数名称 取值 杨氏模量/MPa 4063 抗拉强度/MPa 13.8 对数体积模量 0.005 渗透系数/(mm·s-1) 5×10-7 泊松比 0.33 孔隙率 0.06 摩擦角/(°) 30 通量率/(mm·s-1) 1×10-3 膨胀角/(°) 30 注水时间/s 3600

表1 材料参数 Table 1 Material parameters

算例根据如图2所示力学模型建立有限元模型及网格划分.注水点位于模型中间A点处，宽度为1mm，注水时间为1h.模型长度和宽度分别为101mm和101mm，周围分别受水压力P_w、垂 直压应力P_y和水平压应力P_x作用，数值分别为10.34，29.13和29.13MPa.所采用的单元的多种AR取值以及单个有限元网格形状如图3所示，其中4种典型AR取值下的有限元网格图列于图4中.

图 2(Figure 2)

Fig. 2

图2 力学模型 Fig.2 Mechanical model

图 3(Figure 3)

Fig. 3

图3 相同注水宽度不同AR取值下的单元形状 Fig.3 Element shape under different ARs with the same width of injection line

图 4(Figure 4)

Fig. 4

图4 典型AR下的有限元网格 Fig.4 FEM mesh under typical ARs (a)—AR=9.18； (b)—AR=1.00; (c)—AR=0.50； (d)—AR=0.10.

为了描述水力压裂开裂区域的基本特性，本文考虑三种典型裂缝特征：1) 水力压裂裂缝的垂直影响范围，即裂缝高度(H)；2) 水力压裂裂缝的水平影响范围，即裂缝宽度(W)；3) 水力压裂裂缝的形状比率，即裂缝高宽比(H/W).相同注水时间下，三种本构关系组合(PE-DP,LE-DP和LE-MC)情况下，如图3所示的AR取值的典型12种工况下的水力压裂裂缝分布分别如图5~图7所示.从图5可以看出，采用PE-DP组合，当AR<4.8时，水力压裂的分布区域均不同，且开裂区域不稳定，而当单元高宽比AR=0.1(见图5a)时，数值计算的收敛性也受到影响.

采用LE-DP的本构关系组合时，计算结果有了明显改善，如图6所示.当AR≥0.56时，水力压裂的传播高度基本一致，而开裂宽度在AR≤2.46时不尽相同.当AR=0.1时，数值计算不收敛.

采用LE-MC的本构关系组合时，在AR≥1.0时，如图7所示,计算结果中水力压裂的传播高度基本一致.在AR<1.0时，水力压裂的计算区域出现分岔等不稳定情况，虽然AR=0.1时数值计算可以收敛，但其开裂方向却受到网格形状的控制而沿着水平方向(即单元的长边方向)发展.

图 5(Figure 5)

Fig. 5

图5 水力压裂破坏区域(PE-DP) Fig.5 Hydraulic fracture zones (PE-DP)

图 6(Figure 6)

Fig. 6

图6 水力压裂破坏区域(LE-DP) Fig.6 Hydraulic fracture zones (LE-DP)

图 7(Figure 7)

Fig. 7

图7 水力压裂破坏区域(LE-MC) Fig.7 Hydraulic fracture zones (LE-MC)

图8~图10为3种本构关系组合情况下，AR取值在0.1~19共67种工况下的裂缝开裂区域特性对比.可见，当采用PE-DP组合(图8)，AR<2.8时，开裂区域的高度、宽度变化极其不稳定.当AR≥2.8时，裂缝高度在平均值上下波动.当采用LE-DP组合(图9)，AR≥2.6时，裂缝高度维持在h₁上下波动.当0.5h₂上下浮动，且h₂>h₁.而对于图10所示的LE-MC的本构关系组合，当AR≥1.0时，裂缝高度在某一平均值上下波动.在AR<1.0时，计算的开裂区域不稳定，如出现图7所示的分岔等特征.

图 8(Figure 8)

Fig. 8

图8 水力压裂特性(PE-DP) Fig.8 Features of hydraulic fracture zones (PE-DP)

图 9(Figure 9)

Fig. 9

图9 水力压裂特性(LE-DP) Fig.9 Features of hydraulic fracture zones (LE-DP)

图 10(Figure 10)

Fig. 10

图10 水力压裂特性(LE-MC) Fig.10 Features of hydraulic fracture zones(LE-MC)

裂缝高度随着AR值的增加波动范围增大，因为在保持单元宽度相同时，AR值的增加导致单元高度增加.降低了裂缝高度数据的提取精度，进而引起裂缝高度的波动.当采用PE-DP组合，出现不稳定的裂纹发展较早(AR<2.8)，这是由于渗透系数修改的耦合项和多孔介质弹性本构关系的非线性影响，而采用线弹性本构关系则减少了一个耦合因素，因此对于LE-DP/MC组合，当AR<1时才开始出现不稳定.由于渗透系数耦合项的存在网格敏感性是不能回避的.

此外，对于不同的本构关系组合存在着不同的合适AR取值.对于PE-DP的本构关系组合，当AR>2.8时可以得到稳定结果，但裂缝的发展受单元形状影响，对于已知裂缝开裂方向的工程问题可以控制网格的长边方向平行于开裂方向；对于LE-DP/MC的本构关系组合，当AR≥1时即可以得到稳定的裂缝高度.当AR=1时，可以认为网格形状不影响裂缝发展方向，所以对于未知开裂方向的工程问题建议控制AR为1较为合适.对于已知裂缝开裂方向的工程问题,当采用PE-DP组合时，建议AR的取值范围是2.8~10.0；当采用PE-DP/MC组合时，建议AR的取值范围是1.0~8.0.

1) 通过大量数值计算结果的对比分析可知，无论采用何种本构关系的组合，均不能完全回避弥散裂缝模型中的网格敏感性问题.

2) 对于已知水力压裂裂缝传播方向的工程问题可以通过恰当地控制单元高宽比AR的取值得到稳定的开裂区域.当采用PE-DP，LE-DP和LE-MC 的本构关系组合时，建议AR临界值范围依次为2.8~10.0，1.0~8.0和1.0~8.0.综合考虑所有本构模型，建议AR取值在2.8~8.0之间.

3) 对于未知开裂方向的问题，本文建议使用LE-DP或LE-MC准则，且控制单元的AR取值为1.0.