流固耦合作用下页岩破裂过程的数值模拟

娄义黎，邬忠虎，王安礼，左宇军，刘 镐，孙文吉斌

(1.贵州大学土木工程学院，贵州 贵阳 550025；2.贵州省质安交通工程监控检测中心有限责任公司，贵州 贵阳 550000；3.贵州大学矿业学院，贵州 贵阳 550025)

近年来，随着石油、天然气等常规油气能源的消耗，清洁高效的页岩气成为勘探开发的焦点，并有望成为常规油气资源的理想替代品。勘探发现，我国页岩气储量巨大，在油气资源开发中具有很高战略地位[1-5]。由于页岩气储层物理特性差、渗透性极低，给页岩气的开发带来了极大的困难[6-7]。目前，主要利用压裂技术对页岩进行压裂改造、增透，实现页岩气的开发[8-9]。因此，开展页岩力学特性、破裂过程的研究是迫切需要的，这将为页岩气开发过程中水平井设计、压裂效果评估等提供理论依据。

目前，关于页岩力学特性和破坏模式方面，国内外学者已做了大量的研究，并取得了显著的成果。H.Niandou 等[10]分析了Tournemire 页岩的变形特性、力学特性和破裂模式，发现页岩的塑性变形和破裂模式均表现出明显的各向异性。F.Valès 等[11]对Tournemire 页岩进行了单轴和三轴压缩实验，研究了Tournemire 页岩的力学特性与流体饱和度之间的关系，得出页岩的力学特性对流体饱和度较为敏感。Chen Ping 等[12]对四川盆地龙马溪组页岩样品进行了显微压痕实验，分析页岩力学性质，发现页岩细观弹性模量和压痕硬度分布不均匀，且细观弹性模量平均值与宏观模量值近似。Qiao Lyu 等[13]对含CO2的页岩进行了单轴抗压实验，发现CO2吸附降低了页岩的抗压强度，提高了页岩的延展性。杨志鹏等[14]通过巴西圆盘劈裂实验对含层理页岩进行了研究，发现页岩的抗拉强度随着层理倾角的增加逐渐减小，且不同层理方向的页岩表现出3 种破坏模式。候振坤等[15]对页岩进行单轴压缩实验，结果表明，页岩的抗压强度受层理影响显著，弹性模量随层理倾角的增加而降低，且页岩表现出张拉破坏、滑移破坏及剪切破坏的3 种破坏模式。衡帅等[16]通过真三轴实验研究了页岩层理方向对水力裂缝延伸的影响，结果表明，页岩断裂韧性和水力裂缝的扩展受层理的影响显著。魏元龙等[17]对含裂缝页岩开展了单轴实验研究，发现裂隙的存在使页岩屈服应力、破裂压力和峰值强度降低，页岩主要表现出拉剪贯通和拉贯通两种破坏模式。随着研究的深入，研究者在页岩气储层物性和渗透性等方面做了大量的研究。Chen Dong 等[18]分析页岩裂缝渗透率与储层压力的关系，对不同渗透率数据进行拟合，发现页岩渗透率随着模型系数、裂缝压缩性的增大而减小，裂缝渗透率与储层压力呈正相关关系。Wu Zhonghu 等[19]通过岩石力学实验研究黔北地区下寒武统页岩孔隙率和渗透率，得出页岩的孔隙率和渗透率随有效应力的增大而减小。张宏学等[20]基于裂隙法向刚度的概念，提出了与有效应力相关的页岩储层渗透率模型，并对页岩渗透率进行了测试，发现页岩的渗透率随着有效应力的增大而降低。由以上研究可知，目前，页岩气储层的研究主要集中于力学特性、单轴或三轴加载条件下的破坏模式和渗透特性等方面，而关于页岩在流固耦合作用下破裂过程及声发射特征的研究鲜见报道，因此，开展这方面的研究将对页岩气水力压裂开采提供重要的理论支持。

笔者以黔北地区下寒武统牛蹄塘组页岩为研究对象，通过RFPA2D-Flow 建立7 组不同层理倾角页岩的数值模型，并对该模型进行渗流–应力耦合数值模拟，分析页岩试样的破裂过程，并研究该过程中声发射信号演化规律。

1 实验方法

数值模拟技术被喻为第三种科学研究方法，在科学研究领域应用广泛。RFPA2D-Flow 是基于有限元理论，充分考虑岩石破裂过程中伴随的非线性、非均匀性和各向异性等特点而提出的可用于分析岩石在渗流、应力条件下的数值分析系统。在通常的物理岩石力学实验中，很难取得足够的不同方位、不同角度的岩心，但在RFPA2D-Flow 中，可以通过设置单元的力学参数来建立各种不同力学参数的数值模型，并利用该模型来模拟岩石在渗流–应力耦合作用下的破裂全过程[21-22]。图1 是通过RFPA2D-Flow建立的数值模型，模型中单元颜色的灰度表征材料弹性模量的大小，亮灰色的单元弹性模量较大，图1 中亮灰色的单元是页岩基质；灰黑色的基元弹性模量较小，图1 中灰黑色的单元是层理面。

图1 页岩加载示意Fig.1 Schematic diagram of shale loading

本文建立了7 组不同层理倾角的页岩数值模型，分别为0°、15°、30°、45°、60°、75°和90°。每组模型均为圆柱体(直径50 mm，高100 mm)，加载方式采用位移加载 ΔS=0.000 2 mm，固定围压p1=10 MPa，渗透压差Δp=p3–p2，渗透压取8 MPa。在RFPA2D-Flow 中以水为流体，通过在数值模型中设置渗流边界条件，使模型上下边界形成渗透压差，进而控制流体在储层中的分布。

假设在该模型中页岩弹性模量和强度服从Weibull 随机分布[21]，即：

式中s表示细观元素的属性，如弹性模量或强度；s0与细观元素参数的平均值有关；m表示材料的均质度，均质度m越大，材料越均匀。

试样的初始力学性质参数，如均质度m，弹性模量E，单轴抗压强度σc，孔隙水压力系数等具体值见表1[23]。

2 结果及讨论

2.1 页岩力学特性及破裂过程的分析

图2 是不同层理倾角页岩的应力–应变关系曲线。由图2 可知，层理倾角不同，页岩的应力–应变关系曲线有较大差异，但仍表现出相同的力学作用阶段，大体可以分为弹性、屈服和破坏3 个阶段：初始加载阶段未出现明显的压密(实验采用位移加载，且设置了初始位移)；弹性阶段，应力与应变呈线性关系，且样品表现出较好的弹性特征；屈服特征阶段，其应力–应变关系曲线斜率减小。受围压和渗透压的作用，页岩的脆性指数降低，相较于单轴压缩，在流固耦合作用下页岩的变形过程表现出较明显的屈服阶段。峰值过后，继续施加载荷，应力几乎竖直跌落，虽然在围压和渗透压的作用下页岩脆性指数有所降低，但围压和渗透压较低，页岩仍显现出较强的脆性特征，这与左宇军等[24]对黔北页岩的研究结果基本一致，如图3 所示，但在数值上有所差异，这是因为页岩储层所处地质条件不同，所受构造应力环境不同，造成页岩储层的力学性质不同。本文所取页岩参数较低，且伴有层理弱面，故实验结果数值较低。

表1 页岩数值模型参数[23]Table 1 Numerical model parameters of shale

图2 不同层理倾角页岩的应力–应变关系Fig.2 Stress-strain relationship of shale with different bedding angle

表2 为模拟所得的不同层理倾角页岩的弹性模量和抗压强度，由表中可以发现，在渗流–应力耦合作用下，页岩的抗压强度和弹性模量较表1 中的初始值显著降低，这是由于页岩基质与层理弱面的弱胶结作用及页岩基质固有各向异性所致；另一方面，在流体的作用下，流体对页岩基质具有一定的软化作用，并对页岩基质和层理弱面的弱胶结作用具有一定削弱，进而对页岩的力学性质产生重要影响。

图3 页岩的应力–应变关系[24]Fig.3 Stress-strain relationship of shale

将表2 中弹性模量和抗压强度分布绘制成图4。由图4 中可以看出，不同层理倾角的页岩，其抗压强度表现出显著的各向异性。当α=0°，即层理面与加载方向垂直时，页岩的抗压强度最大，为40.14 MPa；当α=45°时，页岩的抗压强度最低，为34.63 MPa，且为对应0°试样的0.863 倍。倾角在α=0°～15°和α=75°～90°时，抗压强度变化不大，但α=15°～45°时，随着倾角的增加，抗压强度急剧降低，在α=45°～75°内，随着倾角的增加，抗压强度迅速增加，最终形成V 型变化趋势。层理倾角不同的页岩，其弹性模量表现出显著的各向异性。当α=45°时，页岩的弹性模量最大，为41.32 GPa；当α=15°时，弹性模量最低，为38.56 GPa，为对应45°试样的0.93 倍。倾角在α=0°～30°和α=45°～60°内，随着倾角增加，弹性模量变化不大；而在α=30°～45°内，弹性模量随着倾角的增大而迅速增大；在α=60°～90°内，弹性模量随着角度的增加先降低后增加。

表2 页岩弹性模量和抗压强度的模拟结果Table 2 Simulation results of elastic modulus and compressive strength of shale

图4 不同层理倾角页岩的抗压强度和弹性模量Fig.4 Compressive strength and modulus of elasticity of shale with different bedding angle

为表示层理特性对岩石力学特性的影响，文献[25]提出了用层理效应系数S(α)来表征层理结构效应，即不同层理方向页岩的抗压强度与层理倾角α=0°时抗压强度的比值。

式中χ(α)为不同层理方向下页岩的抗压强度；χ(0°) 为层理倾角α=0°时的抗压强度。

将实验数据代入式(2)中计算，得到不同层理倾角页岩的层理效应系数(表3)。由表3 可以发现，当α=45°时，抗压强度的层理效应系数最大，为0.14，表明45°试样的层理结构效应最为显著。且层理效应系数随着层理倾角的增大，大致表现出先增加后降低的变化趋势。因此，可知层理方向对页岩抗压强度有着显著的影响。

表3 页岩抗压强度的层理效应系数Table 3 The bedding effect coefficient of shale compressive strength

图5a 是不同层理倾角页岩的破裂过程和声发射图片，图5b 是文献[24]通过物理实验所得页岩的破坏模式。从图5 中可以发现，30°试样数值模拟结果与物理实验的破坏模式相同，去除层理弱面对页岩破坏模式的影响，两种方法得到的页岩试样破坏模式基本一致，其结果表明RFPA2D-Flow 能很好地对层理页岩进行渗流–应力耦合模拟。图5a 中可以看出，页岩的破坏模式受层理方向影响显著。基于不同层理倾角试样的裂纹分布，最终可将其分为3种破坏模式。

a.斜I 型(15°、60°、75°、90°) 当α=15°时，首先在试样右边基质中形成初始裂纹，而后裂纹上下延伸，并在试样左下方也萌生出新的裂纹，最终形成两条斜I 型破裂裂纹。当α=60°时，首先在试样右边层理面上出现初始裂纹，继续施加荷载，裂纹沿层理面延伸，并贯通整个层理面，最终形成斜I型破坏模式。当α=75°时，初始裂纹首先出现在试样左上方，而后沿试样左右两边延伸，左边沿层理面扩展，右边扩展方向与层理面大致成45°角，最后形成斜I 型的破坏模式。当α=90°时，试样初始裂纹出现在页岩右边基质中，然后沿上下方向扩展，直至形成斜I 型的破坏模式，尽管90°试样表现出斜I型破坏模式，但可以明显观察到页岩有崩坏的趋势。出现斜I 型破坏模式，主要是由于低角度试样(15°)和高角度试样(60°、75°、90°)均受到拉剪应力，当拉剪强度达到其极限强度时，页岩发生破坏，但在这种情况下，围压的作用较为显著，抑制多裂纹的产生，从而形成斜I 型破坏模式。

b.V 型(0°) 当α=0°时，层理面对页岩整体强度有较强的增韧作用，且压应力垂直于层理面，层理面不会发生剪切破坏。随着载荷的增加，页岩基质中拉应力达到峰值并产生微裂纹，载荷继续增加，裂纹向上下方向延伸，与层理面大致成45°角，且在试样左下方萌生新的裂纹，直至形成V 型破坏模式。

c.火焰型(30°、45°) 当α=30°时，在加载过程中页岩受到拉应力和剪切力的作用，当拉剪力超过其拉剪强度时，在试样中下部出现初始裂纹，随着载荷的增加，裂纹向各个方向延伸，尽管围压可以抑制多裂纹的产生，但试样的层理效应较为显著，从而形成较复杂的火焰型破坏模式。当α=45°时，试件的破坏过程与30°试样相似，但其破坏程度不及30°试样。

图5 声发射(AE)图片中不同颜色表示不同的破坏特征，其中红色表示单元发生拉伸破坏，白色表示单元发生剪切破坏，黑色表示已破坏。从声发射图片可以看出，不同层理倾角试样其主要破坏特征为拉伸破坏，并伴随着剪切破坏。这是因为页岩抗压不抗拉，加载过程中页岩出现拉应力集中，随着载荷的增加，拉应力率先达到抗拉强度，从而使页岩发生破坏。

图5 页岩破坏模式Fig.5 The failure mode of shale

2.2 页岩破坏过程中声发射分布演化特征

据文献[26]，RFPA2D-Flow 数值模型中每个破坏单元可以看作是一个微破裂，而微破裂的连通则形成宏观裂纹，因此，破坏单元数和位置特征反映了岩石试样的损伤程度。假设模型中一个单元微破裂代表一个声发射事件，而每一步中破坏单元总数可作为当前时刻的声发射率，加载过程中累计破坏单元数即为累计声发射数，如果能够统计出加载过程中试样破裂的声发射事件，就能得到试样在破裂损伤过程中的声发射演化规律。不同层理倾角下，页岩流固耦合实验的应力、声发射计数、累计声发射与步数的关系如图6 所示。从图6 中可以发现，页岩试样的应力–步数曲线与声发射–步数曲线有很好的相关性。在加载过程中，除了45°试样的声发射变化分为4 个阶段，其他不同层理倾角试样的声发射变化可分为3 个阶段，且低倾角(0°～30°)试样的声发射演化规律不同于高倾角(60°～90°)试样。

a.当α=0°、15°、30°，即低倾角试样时，累计AE 曲线呈现出“平缓线性—台阶—平缓”的变化规律。在应力加载初期，试样处于压密阶段，没有声发射信号出现，继续施加载荷，试样出现微裂纹，声发射信号随之出现，当载荷达到峰值应力的90%左右，微裂纹延伸扩展，进而出现声发射第一次小幅激增，在峰值应力后，页岩中微裂纹急剧扩展，直至贯通形成宏观裂纹，而声发射计数出现大幅激增达到峰值，在这之后，页岩已完全破坏，AE 计数趋于平缓。

b.当α=60°、75°、90°时，累计AE 曲线呈现出“平缓—线性—陡增”的变化规律。在加载初期，与低倾角试样相同，未出现声发射信号，继续施加载荷，页岩中出现微裂纹，累计AE 呈线性增加，当载荷达到峰值应力时，AE 计数跳跃式增加达到峰值，累计AE 激增至最大，峰值过后，声发射的变化情况和与低倾角试样相同。

图6 不同层理角度下页岩的加载应力、声发射计数、累计声发射计数与步数的关系Fig.6 Relationship among loading stress,acoustic emission count,cumulative acoustic emission count and number of steps in shale under different bedding angles

c.当α=45°时，声发射的演化规律较高倾角试样和低倾角试样有所不同。加载过程中，AE 计数出现两个峰值，累计AE 曲线表现为“平缓—线性—激增—平缓—激增”的变化规律。在峰值应力前，累计AE 变化规律与低倾角试样相同，但在峰值应力后，AE 计数出现第二个峰值，累计AE 出现大幅激增。这是因为载荷达到峰值应力后，裂纹并未完全贯通，页岩没有完全破坏，而加载仍在继续，裂纹完全贯通，页岩完全破坏，进而AE 计数突然增大达到第二个峰值。

总体来看，不同倾角试样的声发射演化规律不同。从累计AE 峰值来看，层理倾角α=90°试样的累计AE 最大，虽然90°试样破坏模式为斜I 型，但其伴有崩坏现象，所以其累计AE 最大是合理的；层理倾角α=30°和45°试样次之，层理倾角α=15°、60°、75°试样的累计AE 较为接近，且最小。这也与上述页岩破坏模式分析结论一致，破坏模式越复杂累计AE 越大，因此，声发射信号可以较好地反映页岩的破坏规律。

3 结论

a.流固耦合下页岩的破裂过程可以分为弹性、屈服和破坏3 个阶段，且在破坏过程中表现出明显的脆性特征。

b.页岩的抗压强度和弹性模量受层理方向的影响显著，均表现出明显的各向异性。随着层理倾角的增加，抗压强度峰值曲线呈V 型。其中，0°试样抗压强度最大，45°试样抗压强度最低，是最大抗压强度的0.863 倍。层理效应在45°时最为明显。层理倾角逐渐增加，弹性模量表现出减小–增加–减小–增大的变化趋势。其中45°试样弹性模量最高，0°试样最低，是最大弹性模量的0.93 倍。

c.流固耦合下页岩的破坏模式可分为斜I 型(15°、60°、75°、90°)、V(0°)型和火焰型(30°、45°)3种破坏模式。

d.不同层理倾角下，页岩破裂过程中声发射信号呈现出不同的演化规律。当α=45°时，累计AE 曲线表现为“平缓–线性–激增–平缓–激增”的变化规律，最为复杂；低层理倾角(0°、15°、30°)试样的累计AE 曲线变化规律较复杂，表现为“平缓—线性—台阶—平缓”的变化规律；高层理倾角(60°、75°、90°)试样的累计AE 曲线变化规律最简单，呈现出“平缓—线性—陡增”的变化规律。

e.黔北地区牛蹄塘组页岩的力学特性、破裂过程及声发射能量演化规律，对黔北页岩气开采过程中水力压裂和水平井设计提供重要的理论支撑。