硬脆性页岩断裂韧性二维数值模拟研究

张景轩，范 晓，陈 波，刘向君，梁利喜，熊 健

(1.西南石油大学油气藏地质及开发工程国家重点实验室，四川 成都 610500；2.东方地球物理公司物探技术研究中心成都分中心，四川 成都 610000)

页岩储层作为页岩气赋存的场所，具有自生自储的特点，储层特征表现为极低孔、极低渗。页岩气资源的高效开发需要水平井钻井、分段压裂改造等一系列关键技术[1]。压裂过程中裂缝起裂和延伸受页岩断裂韧性的影响，当裂缝尖端应力强度因子超过页岩断裂韧性值时，裂缝开始扩展[2]。因此，研究页岩断裂韧性及其影响因素对评价页岩储层可压性有重要意义。

目前，针对岩石断裂韧性研究，主要采用室内试验测试和数值模拟方法[3-5]。衡帅[6]开展切口与层理呈不同方位的圆柱形试样三点弯曲试验，研究页岩断裂韧性的各向异性特征，并揭示其断裂机制的各向异性。Atkinson等[7]通过直缝巴西圆盘实验法(SNBD)测量断裂韧性，并给出了I型和II型断裂韧性计算公式。崔振东等[8]采用人字形切槽巴西圆盘法(CCNBD)，研究泥质砂岩的断裂力学特性，讨论了该试样类型的有效尺寸和断裂机制。张盛等[9]分别以5种不同形状的圆盘试件测定了大理岩断裂韧度，提出了用每种圆盘确定断裂韧度的公式。由此可见，室内测试岩石断裂韧性的方法较多，但众多方法都有各自的适应性。国际岩石力学学会(ISRM)推荐人字形切槽巴西圆盘法为岩石断裂韧性的测试方法[10-11]。因为该方法的制样简单，测试结果离散度较低，结果可靠。同时，ISRM也提出了岩石断裂韧性测试样品的最小有效试样直径为75 mm[12]。研究岩石断裂韧性的数值方法，包括位移不连续法(DDM)、有限单元法(FEM)、扩展有限元法(XFEM)、离散元法(DEM)和边界元法(BEM)等[13]。满轲[14]在LS-DYNA程序的基础上，进行了玄武岩动态断裂韧性测试，得到了动态断裂韧性随着加载率增大而增大的规律。程远方等[15]用拟三维模型描述煤层气井压裂裂缝扩展行为时提出，断裂韧性是缝高的重要影响因素。由此可见，岩石断裂韧性特征的研究已经取得了一定的成果。针对硬脆性页岩断裂韧性研究，陈建国等[16]基于SNBD法研究了渝东南酉阳地区的龙马溪组页岩断裂韧性特性，并基于地球物理测井信息，建立了页岩断裂韧性的预测模型；梁利喜等[17]研究了四川长宁地区龙马溪组页岩断裂韧性的预测模型，并讨论了矿物组成、TOC含量对其影响。赵金洲等[18]基于页岩储层断裂韧性、脆性等特征，提出了表征页岩气可压性的评价方法；张明明[19]开展了三点弯曲半圆盘法测试页岩Ⅰ-Ⅱ复合型断裂韧性的实验和模拟，建立了裂纹扩展角的预测方法。然而，这些研究并未考虑层理角度、层理密度及层理面力学性质等对页岩断裂韧性的影响，以及这些因素对页岩裂缝前缘形态的影响。

在前人的研究基础上，以龙马溪组页岩为研究对象，运用以统计损伤的本构关系与有限单元法结合的RFPA2D岩石破裂过程数值模拟软件，构建人字形切槽巴西圆盘法数值模拟过程，从层理角度和层理面密度以及层理面力学性质等角度，研究页岩断裂韧性的变化规律；在此基础上，讨论页岩裂缝前缘形态的影响因素。同时，也研究了页岩断裂韧性与脆性指数间的关系。

1 模型参数选取与验证

1.1 数值模型介绍

RFPA2D(Rock Failure Process Analysis)用弹性损伤力学的本构关系来描述岩石细观单元的力学性质，在考虑到细观非均匀性的前提下，将岩石的非均匀性参数引入到基元中去，细观基元的性质是简单的弹塑性或者弹脆性，当基元应力达到破坏准则时则发生破坏，采用刚度退化处理，达到用连续介质力学方法处理物理非连续性的问题，从而实现岩石宏观非线性[20]。

按照应变等价原理，应力作用在受损材料上引起的应变与有效应力作用在受损材料引起的应变等价，受损材料单元的本构关系可通过无损材料中的名义应力得到[21]，即：

(1)

式中：E和E0分别为损伤后单元的弹性模量和初始弹性模量，MPa；D为损伤变量，D=0时，表示材料无损伤，而当D=1时，表示材料完全丧失承载能力。

1.2 参数选取

为研究层理面对断裂韧性以及裂缝前缘形态的影响，进一步探讨页岩断裂韧性及诱导裂缝前缘形态影响因素，利用RFPA2D构建了层理性页岩的巴西圆盘法数值模拟。该软件用统计损伤的本构关系来描述岩石材料的非均匀性和缺陷分布的随机性，再将细观单元性质的统计分布与有限单元法结合，以弹性力学为应力分析工具，以弹性损伤理论及修正后的Coulomb破坏准则对单元进行变形及破坏处理，从而实现对非均匀材料破坏过程的模拟，设计了如图1所示的技术路线。

图1 参数确定流程

由于RFPA2D系统中所输入的模型参数不等同于岩石力学实验得到的岩石宏观物性参数，而是细观单元参数的平均值。岩石的孔隙度、渗透率、密度三个参数对裂缝的形态影响很小，在这里假设模型中这几个参数均匀分布，故本研究以龙马溪组页岩基础物性和力学参数值作为参考，主要确定细观弹性模量、抗压强度、泊松比三个力学参数。

首先，以物理实验获得的力学参数作为参考，对模型输入的各细观力学参数平均值，以及影响该力学参数大小的相关参数进行调节，初步确定力学模型参数。然后，对模型的参数进行验证，以该套参数建立数值模型，并开展单轴压缩实验模拟，得到数值模型的宏观力学参数，将获得的力学参数与物理实验力学参数进行对比。按照上述方法，结合测得的龙马溪组岩心基本物性参数与力学参数的平均值，对软件输入的各参数进行调节，主要是力学参数的调节。最终确定模型参数见表1所示。

以表1中模型参数建立模型，模型尺寸为50 mm×25 mm的标准岩心，取层理面密度的平均值0.45条/mm，将模型共划分网格400×200个。计算时采用位移加载，单步增量为0.005 mm。其模拟过程及计算结果如图2所示。

表1 模拟参数

图2 计算前后模型

1.3 模型验证

根据应力-应变数据，绘制出单轴压缩下数值模型的应力-应变曲线，取其中极大值为抗压强度。取应力应变曲线中弹性阶段进行拟合，通过拟合直线的斜率可得到弹性模量、泊松比的数据，如图3所示。通过计算，单轴压缩情况下，抗压强度为98.872 MPa，弹性模量为243.77 MPa，泊松比为0.202 3。

图3 力学参数计算

对比表2中数据，物理试验与模拟计算所得结果较吻合，与第一组试验的抗压强度误差为1.9%，弹性模量误差为1.6%，泊松比误差为1.8%；与第二组试验的抗压强度误差为2.7%，弹性模量误差为0.2%，泊松比误差为1.2%，抗压强度、弹性模量、泊松比误差均在5%以内。因此，该模拟参数能较好地反应岩心的真实情况，适用于该龙马溪页岩的模拟计算。

表2 物理试验力学参数与RFPA模拟宏观力学参数对比

2 实验结果讨论与分析

针对静态或准静态I型断裂韧性测试，国内外学者作了大量研究，测试方法包括：V型切口的三点弯曲圆棒试样(CB)、短棒试样(SR)以及人字形切槽巴西圆盘试样(CCNBD)等。CCNBD试样在很多方面优于CB和SR试样[8]：①试样体积较小；②加工工艺较简单，免除了加工尖锐直裂纹的困难，也能承受较高的临界载荷；③不仅适合测试I型裂纹的断裂韧度值，而且很容易实现纯II型或I-II复合型裂纹断裂模式。故选取人字形切槽巴西圆盘试样进行断裂韧性模拟，几何参数如图4所示。

图4 CCNBD试样示意

采用式(2)计算CCNBD试样断裂韧性值[17]：

(2)

由以上CCNBD测试岩石断裂韧性值的原理可知，断裂韧性与材料本身的物理力学参数无直接相关，仅与施加于试样上的荷载以及试样的几何参数(D,B,a,a0,a1)有关。又由于模拟研究的各模型几何参数保持一致，故可用RFPA2D程序计算的岩石破坏应力反映断裂韧性的变化规律。首先基于岩石三轴实验模拟，确定出页岩的三轴力学参数，结合巴西圆盘的几何参数，分别从层理的角度(即与加载方向的夹角，以10°为间隔依次从0°变化至90°)和层理的密度(以0.2条/mm为间隔依次从0.2条/mm变化至1.0条/mm)以及层理面力学性质(0.25～0.75倍基质强度折减)建立人字形切槽巴西圆盘模型，采用位移加载，单步增量取0.002 mm，模型运算条件如表1所示。

2.1 层理面对页岩裂缝前缘形态影响

层理面是影响裂缝前缘形态的一个重要因素，由于层理面弱面的影响，裂缝在延伸至层理面时，可能出现转向、分叉或穿过现象。在模型中划分不同层理面，分别从层理面角度、层理面线密度、层理面力学性质三方面分析层理面对页岩裂缝前缘形态的影响。

固定层理的面密度，以0.4条/mm为例，当页岩层理与加载方向呈不同角度时，其裂缝前缘形态的变化如图5所示。图中可见，各个页岩试样均在预制裂缝处开始宏观断裂，但裂纹的扩展路径却差别较大：当层理夹角为0°时，在加载过程中裂缝沿着预制裂缝尖端起裂和延伸，未发生转向，这可能是因为裂缝和层理都作为弱结构面且都相互平行；当层理角度小于40°时，在加载过程中，裂缝起裂时沿着预制裂缝延伸，当遇到层理后，裂缝并没有穿过层理而是发生转向，继而沿着层理方向延伸，此时，裂缝转向角近似等于层理角度，裂缝延伸过程中层理起到较重要的作用；当层理夹角在40～70°之间时，除了在延伸时前缘随层理出现转向外，圆盘端部出现多条沿层理方向的次生裂纹，这可能是由于在峰值载荷后轴向压应力骤降，使卸荷后的拉张应力和剪切应力集中于试样的端部区域，而层理面的胶结强度较小并且力学性质相对薄弱，最终导致圆盘端部沿层理面出现多条次生裂纹；当层理夹角大于70°时，层理在加载方向上的间隔变小，圆盘在轴向应力与拉张应力的作用下，层理面的剪切应力作用减弱，裂缝先沿预置裂缝方向扩展，后穿过层理发生分叉，并且出现次生裂纹的数量明显减少。

裂纹自尖端起裂后，弹性应变能急剧释放，表现有峰后载荷骤降、起裂偏离角度、次生裂纹产生等[8]，转向后的裂缝沿层理面继续延伸扩展，且在沿层理面继续延伸的过程中会进一步沟通水力裂缝或天然裂缝而形成相对复杂的裂缝网络[6]。以层理角度30°、60°、90°为例，分析密度0.2条/mm、0.4条/mm、0.6条/mm的裂缝前缘形态(见图6)，层理密度对裂缝前缘形态影响规律较复杂：①层理密度较小时，裂缝在预置裂缝方向开裂后，受层理面影响转向层理面方向，最终贯穿圆盘；若角度增大，伴随有次生裂纹出现，增大至80°以上时，裂缝穿过层理面，会出现转向和分叉；②随着层理密度增大，裂缝不均匀扩展现象增多，圆盘端部伴随产生的次生裂纹数量也逐渐增加，在角度较小时最明显；在角度增大后，裂缝穿过层理面，并出现裂缝分叉。由此，进一步说明层理面对于裂缝前缘形态有着重要影响。

图5 层理密度0.4条/mm时的结果

图6 不同层理密度时不同层理倾角的裂缝扩展结果

图7 不同力学性质层理条件下各层理倾角的裂缝扩展结果

取层理面线密度为0.6条/mm进行研究，分别将层理面力学性质按表1中页岩基质的0.25、0.5、0.75倍强度折减。从图7中可看出，层理面力学性质从层内强度的0.25倍增大至0.75倍过程中，力学性质也对裂缝前缘形态有重要的影响；当层理面与页岩基质的力学性质越接近，层理越不易受剪切应力影响，继而产生的次生裂纹数量减少，这可能是由于层理面的力学性质较小时，在试样端部越易形成应力集中，促使圆盘沿层理产生次生裂纹；在角度较小时，随层理面力学性质增强，试样的破坏形态由顺层理破坏转变为裂缝转向、分叉；当层理角度较大时，圆盘的裂缝前缘形态基本不发生变化，都是起裂转向后穿过较多层理，最终贯穿圆盘。

2.2 层理面对页岩断裂韧性影响

断裂韧性对裂缝扩展压力产生影响，进而影响到裂缝的扩展及裂缝前缘的形态。经过模拟测试的页岩力学参数计算得出其对应的断裂韧性变化规律如图8和图9所示。①在层理角度相同条件下，随着层理密度增大，模型中层理数量增多，数值模拟试样的断裂韧性整体呈减小趋势；②在层理密度相同条件下，试样的断裂韧性随角度变化的先减小后增大，与层理角度呈波折变化的趋势；③在角度为20°时，裂缝在预置裂缝方向开裂后，受层理面影响转向层理面方向，此时的断裂韧性值最小；④在角度为90°时，层理面与加载方向垂直，裂缝在延伸的过程中，不易穿切层理面，故此时的断裂韧性值最大。

图8 不同层理密度条件下断裂韧性与层理角度的关系

将层理面力学性质按基质强度的1/4、1/2、3/4倍折减，取层理面线密度为0.6条/mm，层理面角度为0～90°时为例，分析断裂韧性的情况，如图10所示。在层理角度相同的条件下，随着层理面强度由层内岩石强度的1/4增大至3/4,数值模拟试样的断裂韧性逐渐增大；随着层理面强度逐渐减小，岩石呈现的断裂韧性值减小得越快；在层理面线密度为0.6条/mm时，随着层理面强度由层内岩石强度的1/4增大至3/4，同一层理面角度间断裂韧性值差异随层理面角度增大而先增后减。

图9 不同层理角度下断裂韧性与层理密度的关系

图10 不同层理力学性质下断裂韧性与层理角度的关系

2.3 页岩断裂韧性与脆性指数的关系

页岩的脆性也能够显著影响井壁的稳定性，关系到设备效率和安全钻进时泥浆密度的选择；是评价储层力学特性的关键指标，同时还对压裂的效果产生显著影响。对于页岩气藏，页岩的脆性越高，压裂时裂缝形态越复杂，越容易形成缝网，达到体积改造的目的；若地层的脆性较低，在压裂时地层不易被压开，或者是形成对称双翼缝，延伸时裂缝前缘形态较为单一，不利于裂缝缝网的形成[22]。

国内外学者提出了数十种表征岩石脆性的指标，主要有基于强度、硬度、应力应变、贯入实验、破碎程度和应变能等六大类岩石脆性评价方法，各有其局限性，有的并不适合评价页岩[22]。采用基于岩石峰后应力跌落速率及能量比的岩石脆性特征评价方法来评价页岩脆性，能够充分反应岩体本身力学参数，及外在加载条件变化对于脆性破坏特征的影响[23]。通过计算峰后应力降速率BPOST和岩石破坏时所释放弹性能与峰前储存总能量比值BE，对脆性指数Bi进行综合评价[24]，其计算公式如式(3)：

(3)

式中：A为修正系数，取值0.001；σp与为σr峰值应力和残余应力，MPa；εp与εr为峰值应变和残余应变。

采用与之前相同的数值模拟理论及方法，以龙马溪地层岩石物理实验数据为基础，利用RFPA软件对该地层进行数值模拟实验。分别从层理的角度(以10°为间隔依次从0°变化至90°)和层理的密度(0.4条/mm，0.6条/mm)两方面来研究其对脆性指数的影响。根据图3中的应力-应变曲线，结合公式(3)可计算出页岩的脆性指数，其结果如图11所示。在此基础上，在相同层理参数条件下，统计其断裂韧性值与脆性指数间的关系，其关系如图12所示。

图11 脆性指数与层理角度的关系

从图11中可看出，页岩脆性指数有随层理角度先增大再减小的趋势；同层理密度条件下，层理角度20°时岩石的脆性指数最大，原因是此时破坏形态主要为裂缝尖端起裂后随即沿层理面的剪切滑移破坏，破坏形态较为单一；而角度90°时破坏形态的脆性指数最小，此时的层理角度与位移加载方向垂直，相比角度20°时较易沿层理直接贯穿圆盘试样，其破坏模式为穿切层理面的劈裂型层状剪切破坏，故其表现的脆性较小；层理面作为片理状构造的弱结构，随着层理的密度减小，岩石的脆性随之减小。

同时，从图12中可看出，页岩断裂韧性与其脆性指数之间呈较好的负线性关系，即随着岩石的脆性指数增大，表现的断裂韧性减小。说明了页岩的脆性越强，其可压性越高，故其断裂韧性越低。

图12 不同层理密度和角度下脆性指数与断裂韧性的关系

3 结论

以四川盆地长宁地区龙马溪组为研究对象，在已有岩石力学特性的基础上，用RFPA建立人字形切槽巴西圆盘试样模型，模拟计算分析层理面对页岩断裂韧性的影响以及断裂韧性测试过程中的裂缝前缘形态，得到了以下结论：

(1)随着层理密度增大，断裂韧性整体呈减小趋势；断裂韧性随角度的变化先减小后增大，在20°时出现最小值；在层理角度相同的条件下，随着层理面力学性质增强，断裂韧性逐渐增大；在层理面线密度为0.6条/mm时，随着层理面力学性质增强，断裂韧性值差异随层理面角度增大而先增后减。

(2)当层理角度为0°时，试样发生劈裂破坏，裂缝未发生转向；当层理角度较小时，裂缝前缘随层理面发生转向，并沿层理发生脆性剪切破坏；当角度增大时，层理面作用减弱，试样端部出现次生裂纹；当角度增至80°以上时，裂缝在延伸过程中穿过层理会出现分叉现象。

(3)层理密度较小时，裂缝在预置裂缝方向开裂后，受层理面影响转向层理面方向，最终贯穿圆盘，若角度增大，伴随有次生裂纹出现，若增大至80°以上时，裂缝穿过层理面，并出现转向和分叉；层理密度较大时，裂缝不均匀扩展现象增多，圆盘端部伴随产生的次生裂纹数量也逐渐增加，在角度较小时最明显，而当角度增大后，裂缝穿过层理面，并出现裂缝分叉。

(4)脆性指数与断裂韧性之间呈现负相关关系，即随着岩石的脆性指数增大，表现的断裂韧性减小，并且两者之间有较好的相关性；脆性指数有随角度先增大再减小的趋势。