单轴压缩下页岩类脆性材料细观损伤演化的离散元分析

钟寄宸， 王 蕉*， 李翔宇,2， 楚锡华

(1.西南交通大学 四川省应用力学与结构安全重点实验室,力学与航空航天学院,成都 610031；2.西南交通大学 牵引动力国家重点实验室应用力学研究所,成都 610031；3.武汉大学 工程力学系，武汉 430072)

1 引 言

在宏观尺度上，页岩是典型的横观各向同性材料，其破坏通常表现为脆性破坏，难以监测。研究者有必要对页岩的力学特性进行测试，以确定页岩的各向异性对力学性能、断裂模式和抗压能力的影响因素。对页岩或类似页岩的脆性材料的储层结构和力学性能的清晰认识将为工程实践提供启示。

层状岩石的单轴压缩试验表明，岩石的抗压强度随着层理角的增加而下降[1-3]。层理面对沉积岩力学行为影响的巴西圆盘试验研究表明各向异性导致了复杂的断裂力学行为[4，5]。单轴压缩试验和巴西试验结合起来研究层状岩石的力学行为,结果表明页岩的单轴抗压强度UCS、弹性模量和破坏模式受层理面影响显著。衡帅等[8]总结了龙马溪黑页岩的四种典型破坏模式;Yang等[9]也总结了单轴压缩条件下，岩石类脆性材料的四种典型破坏模式。此外，也有研究探讨了页岩各向异性对井壁断裂的影响[10,11],结果显示，样品的抗压强度曲线呈U型。页岩的力学行为很大程度上取决于层理面的粗糙度。为了进一步研究页岩的各向异性，Ge等[12]基于分形方法研究了岩石节理粗糙度的测量和量化。以上学者通过物理实验对页岩的力学性能进行了研究，并提供了直观的物理现象和可靠的实验数据。但是物理实验受实验条件和实验设备的限制较多，且成本巨大。相比于传统物理实验，数值方法研究成本较低且能收集到更多数据，所以很多学者采用数值方法对页岩进行研究。

数值模拟是岩石力学和岩土工程的一个重要研究方法。在岩土力学中，广泛使用Cundall等[13]提出的离散元法DEM为基础的离散元模拟。采用嵌入光滑节理接触的粘结颗粒离散元模型模拟横观各向同性岩石的力学行为[14-17]。在此基础上，Xia等[18]研究了光滑节理的参数和校准方法，本文采用了同样的方法。Xu等[19]提出了一个新的模型，用以模拟各向异性岩石的力学行为。其中，岩石基体和层理面分别使用平节理接触和光滑节理接触进行模拟。Lin等[20]用离散元法研究了层理面参数对双孔层状岩石力学性能的影响。实验结果表明岩石的抗压强度和弹性模量随层理面的变化呈中间低两端高的形态。杨宇江等[21]引入多重分形理论，计算了单轴加载过程中峰前某一阶段单元体积应变的多重分形谱。除了压缩试验，也有少量关于页岩拉剪行为的研究[23]，发现试样在混合形式的荷载条件下，裂缝通常以拉剪I/II型混合模式起裂，而后以I型为控制方式稳定扩展。上述学者通过数值方法研究了岩石的力学特性，极大地丰富了各向异性岩石力学响应的数据库。但是所做的贡献主要集中在完整岩石的破坏表征和力学性能等方面，对于细观损伤的研究尚不充分。

为了更好地理解和表征岩石内部的细观损伤，本文采用颗粒流模型分析不同层理角的页岩样品在单轴压缩下的破坏情况，引入裂纹张量和裂纹密度，研究微裂纹随加载过程的演化以及层理角对裂纹密度的影响。此外，使用平均配位数，分析配位数与试样应力状态及细观损伤的联系。本文在细观尺度上定义页岩类脆性材料的内部损伤,并将其与宏观应力应变相结合,对横观各向同性脆性岩石的破坏过程和内在机理的理解提供参考。

2 PFC2D数值模型

2.1 接触模型概述

由于实验手段的限制，在物理实验过程中，很难直接观察材料内部的细观运动。为进一步研究岩石破坏过程中材料的内部损伤演化过程，本文采用基于DEM的颗粒流模型，每个样本使用三种接触模型，其中，线性平行粘结模型用于模拟岩石基体，光滑节理模型用于模拟层理面。当颗粒间的力超过粘结强度时，粘结发生断裂，颗粒之间的接触行为由线性接触模型描述。此外，摩擦力在这三种接触模型中均存在，三种模型描述如下。

线性模型中接触力的计算如下。

法向力Fn和切向力Fs为

Fn=Fn+knAΔδn,Fs=Fs-ksAΔδs

(1,2)

式中A为等效面积(与接触颗粒半径相关)，kn和ks分别为线性模型法向和切向刚度，Δδn和Δδs分别为相对法向和相对切向位移增量。

平行粘结模型中接触力的计算如下。

(3，4)

光滑节理模型中接触力的计算如下。

(5，6)

2.2 页岩数值模型描述

将层理角β从0调整到90°(β定义为层理面与x轴的夹角)，建立7个层理面厚度恒定的数值模型(层理面厚度db=1.4×rh i)，如图4所示。考虑到真实岩石不同层的厚度不尽相同，因此页岩各层的厚度设置略有不同。

h1=16.5 mm,h2=21 mm,h3=19.5 mm
h4=21 mm,h5=18 mm,h6=21 mm
h7=24 mm

对层理角为0°～90°的矩形试样(75 mm×150 mm)进行单轴压缩测试。颗粒与壁面的摩擦系数设为0.1。模拟中的每个数值模型大约由40300个颗粒组成。颗粒半径遵循0.21 mm～0.35 mm 的均匀分布。在模拟过程中，上壁用于加载，下壁固定。根据Zhang等[24]的研究，加载速度设定为0.05 m/s(≤0.08 m/s)，时间步长设定为2.5×10-8s，以确保整个加载过程是准静态的。

2.3 页岩材料参数校正

表1 层理面的材料参数

表2 岩石基质的材料参数

图1 线性模型

2.3.1 破坏模式对比

将数值模拟结果与文献[8]的物理实验的破坏模式进行对比，如图2所示。图2(a1～d1)显示了页岩在单轴压缩时的破坏模式，图2(a2～d2)对应为数值模拟加载至模型破坏时的应变分布(应变云图计算过程参见文献[25])，两种方法得到的破坏模式具有相似性。图2(a)的β=0；图2(b)的β=30°，岩石中产生多个贯穿层理面和沿层理面的破坏面；图2(c)的β=60°，岩石中产生了沿层理面的破坏面，其中，贯穿性破坏沿着60°的层理面，导致岩石试样发生明显的剪切滑移破坏；图2(d)的β=90°，岩石中产生多个贯穿层理面和沿层理的破坏面，以后者为主。

图2 平行粘结模型

除物理实验对应的4组层理角外，对层理角为15°，45°和75°的情况进行了数值模拟，破坏模式如 图3 所示。结合图5，对破坏模式进行分析，随着层理角从0°增加到90°，岩石的破坏模式完成了从基质和层理面的混合破坏到沿层理面的单一破坏再到基质和层理面的混合破坏的演化，如图6所示。

图3 光滑节理模型

图4 不同层理角度页岩的数值模型(β =0°和β =30°)

图5 实验和数值模拟的破坏模式对照

图6 不同层理角页岩的破坏模式

2.3.2 宏观力学性能对照

将物理实验和数值模拟所得的抗压强度和弹性模量进行比较，结果如图7所示。考虑到文献[8，26]对于页岩抗压强度描述的差异，本文将文献[8]的物理实验测得的抗压强度和数值模拟的抗压强度做归一化处理以对比抗压强度的变化规律。具体数值列入表3和表4。

表3 数值模型力学性能

表4 不同层理角度页岩单轴压缩试验结果

图7 实验和数值模拟的对照

图7(a)显示了物理实验和数值模拟中不同层理角页岩的抗压强度。数值模拟结果显示为虚线，文献[8]物理实验的结果显示为实线。很明显，从0°～60°，页岩的抗压强度随着层理角的增加而降低，这种下降在30°～90°区间内更明显。这一现象与层理面的弱化效应有关。随着层理角的增加，层理面两侧的基质更容易沿着层理面产生相对运动。由于层理面内颗粒的粘结力和抵抗变形的能力相对较弱，这些区域更容易受到破坏。因此，页岩的抗压强度下降。在60°～90°范围内，页岩的抗压强度随着层理角的增加而明显增加。

当层理角为75°时，岩石的抗压强度增加。因此，在60°～75°区间可能存在某个角度阈值。当层理角超过这个阈值时，层间弱化效应对岩石抗压强度的影响减少，抗压强度增加，页岩的抗压强度在层理面的影响下显示出明显的各向异性。阈值大小还与层理面粗糙度和颗粒排列方式等密切相关，其具体作用机理还需要进一步研究。

图7(b)显示了不同层理角页岩的弹性模量。物理实验和数值结果都表明，页岩的弹性模量随着层理角的增加而增加，在75°～90°区间内变化更为明显。与抗压强度类似，弹性模量在层理面的影响下也显示出明显的各向异性，数值模拟结果与文献[8]物理实验有良好的一致性。

3 细观分析

3.1 裂纹张量定义

(7)

(8)

式中d为裂纹两侧颗粒的平均直径，ni和nj分别为裂纹法线的i和j方向的分量，N为裂纹总数量，A为颗粒的总面积。

3.2 层理面对裂纹密度的影响

图8显示了不同层理角的页岩样本中裂纹密度随轴向应变的关系曲线，可以看出，当层理角不大于60°时，裂纹密度随着层理角的增加而增加，这与层理面的弱化效应有关。随着层理角的增加，层理面两侧的基体更容易沿着层理面产生相对运动。因此，在相同的荷载条件下，裂纹密度随着层理角的增加而增加。

图8 不同层理角页岩的裂纹密度

层理角等于60°时，当应变从0变化到1.33×10-3，裂纹密度迅速增加。继续加载，当应变从 1.33×10-3变化到2.16×10-3，裂纹密度曲线的斜率没有明显变化。可以观察到应力集中的现象大多发生在层理面附近，大量裂纹沿着层理面产生，此时，基体中的应力较小，所以在基体中产生的微裂纹数量较少。裂纹密度的增长主要由沿层理面的粘结破坏引起，当应变达到2.16×10-3时，岩石完全破坏，失去原有的承载能力。

层理角大于60°时，当应变从0到1.33×10-3变化时，裂纹密度增长缓慢，当应变从1.33×10-3到2.16×10-3变化时，裂纹密度迅速增加。这一现象表明，岩石的破坏机制发生了改变，在60°～75°之间存在某个层理角阈值。当层理角超过阈值时，破坏模式将从主要由层理面控制变为层理面和岩石基体共同控制，层理面弱化效应的影响减弱。因此，在初始加载阶段，β=75°的页岩的抗破坏能力比β=60°的强。当层理角增加到90°时，岩石在初始加载阶段抵抗变形的能力进一步增强。

3.3 平均配位数分析

对15°的页岩样本展开进一步分析。样本的应力和裂纹数量曲线如图9所示。根据裂纹数量的变化，整个加载过程可以分为三个阶段，第I阶段，裂纹数量缓慢上升；第II阶段，裂纹数量增长加快；第III段，裂纹的数量迅速上升，试样迅速发生破坏。对照试样的应力曲线，在I和II阶段一直保持上升状态，仅通过应力-应变曲线，无法反映出试样内部的细观损伤。

图9 应力与裂纹数量曲线

为了更好地表征岩石试样的细观损伤，引入平均配位数，分析配位数为3～7的颗粒数量和平均配位数的变化。平均配位数可表示为

(9)

式中Cave为平均配位数，N为颗粒总数量，Ck是id为k的颗粒的配位数。

本文使用的数值模型平均配位数约为 4.6，大于平均配位数的配位数定义为大配位数(5,6,7)，反之为小配位数(3,4)。从加载初时刻起，每间隔1.25 ms统计一次配位数，共统计12次，不同配位数颗粒数量和平均配位数曲线经过归一化处理。0代表12组取样点的最小值，1代表12组取样点的最大值。从图10可以看出，岩石试样的小配位数颗粒数量先下降后上升，大配位数颗粒数量先上升后下降，平均配位数先上升然后伴随较小幅度的下降，然后急剧下降。小配位数颗粒数量曲线及大配位数颗粒数量曲线的转折点都在临界状态附近(黑色竖线代表临界状态出现位置，即应力峰值点附近)。

图10 不同配位数颗粒数量及平均配位数的变化曲线

图11显示了加载过程中平均配位数与裂纹数量的对应关系，可以看出，平均配位数和裂纹数量变化趋势有很好的一致性。从初始状态到临界状态，平均配位数和裂纹数量都缓慢上升。超过临界状态后，裂纹数量迅速上升，平均配位数迅速下降。在样本发生破坏之前(即I和II阶段)，平均配位数和裂纹数量的变化可以相互佐证，I阶段，平均配位数稳步上升，对应裂纹数量从无到有，并线性增加；II阶段，平均配位数有小幅下降，对应裂纹数量继续增加，但增速明显大于I阶段。最后样本发生破坏，平均配位数急剧下降，对应裂纹数量急剧上升。因此，配位数能够在一定程度上反映试样内部细观结构的变化。

图11 平均配位数与裂纹数量曲线

单轴压缩过程中，岩石颗粒受到挤压。在外力的作用下，部分小配位数颗粒的稳定结构破坏，颗粒之间变得紧密以维持内部的稳定状态，配位数上升。此时颗粒之间的相互作用更加紧密，接触力上升；当接触力进一步增大，试样达到临界状态，大配位数颗粒的接触结构发生脆性破坏，迅速分裂成多个小配位数颗粒的接触结构(对应着大配位数颗粒的减少和小配位数颗粒的增多)，裂纹大量增加，试样发生破坏。因此，颗粒集合体在单轴压缩条件下的应力状态及内部细观损伤程度与平均配位数的变化关系是对应的。

4 结 论

本文通过离散元模拟再现了不同层理角的页岩试样在单轴压缩下的破坏模式。在此基础上，系统地研究了层理角对岩石破坏过程的影响。为描述岩石内部细观损伤，引入裂纹张量和平均配位数，研究了单轴压缩条件下，层理角对岩石内部细观损伤的影响以及平均配位数与岩石应力状态和细观损伤的联系。结论如下，随着层理角从0°增加到90°，岩石破坏模式从基质和层理面的混合破坏到沿层理面的单一破坏再到基质和层理面的混合破坏。同等条件下，页岩的裂纹密度会随着层理角的增加而增加。相较于应力应变曲线，平均配位数曲线能够更准确地反映岩石的细观损伤；应力应变曲线峰值滞后于裂纹数量曲线和平均配位数曲线的转折点。尽管在第III阶段初，试样内部已经出现大量裂纹，但裂纹尚未聚结产生贯穿性裂纹，所以试样仍有一定强度，继续加载，裂纹贯通，试样达到临界状态，迅速发生破坏。小配位数颗粒的数量先下降后上升，大配位数颗粒的颗粒数量先上升后下降；颗粒集合体在单轴压缩条件下的应力应变及裂纹数量曲线与平均配位数曲线有良好的一致性。