分段黏结非连续变形分析方法及其在砂岩破裂分析中的应用

2023-10-18邓庆龙王培瑞焦玉勇

煤炭学报 2023年9期

郑飞，邓庆龙，李芷，王培瑞，靳陆，焦玉勇

(中国地质大学(武汉) 工程学院,湖北武汉 430074)

深部煤炭开采中，高地温、高地压、高岩溶水压及开采扰动影响更为显著，岩石力学问题也更为复杂[1]。冲击地压、煤与瓦斯突出灾害与岩石的多尺度断裂行为密切相关，在实验室尺度研究煤层顶底板岩层代表性岩石在动静荷载作用下的变形、断裂、强度等特性及演化规律，对于认识上述灾害的发生规律具有重要的意义[2-5]。

大量试验及相关研究表明，岩石的变形及强度特性与细观尺度的颗粒、微裂纹等特征具有较强的相关关系。精细刻画岩石细观结构，考虑颗粒及微裂纹等岩石细观特征，进行宏观变形及强度特性的分析是重要的研究课题。在岩石的室内实验方面，声发射[6]、CT 扫描[7]、3D 打印[8]等先进的技术得到运用，推动了对岩石强度特性、动力学特性的研究。

数值模拟技术由于其灵活性、定量性、过程可视性等特点，也是岩石多尺度断裂分析研究的重要手段。分析岩石断裂过程的数值方法，通常可划分为基于连续介质假定的方法、基于非连续介质假定的方法、以及混合连续与非连续框架进行求解的方法。基于连续介质假定的方法包括扩展有限元(XFEM)[9]、岩石破裂过程分析方法(RFPA)[10]、边界单元方法[11]等；非连续类数值模拟方法主要包括基于显式接触力求解的离散单元法(Discrete Element Method,DEM)和基于隐式位移求解的非连续变形分析方法(Discontinuous Deformation Analysis，DDA)[12]；混合连续与非连续框架求解的方法包括FDEM[13-14]、CDEM[15]、拉格朗日元-离散元耦合法[16]等方法。

非连续类数值计算方法可以考虑较大尺度断层、节理、层理等不连续面的张开、闭合、滑移，也可以模拟较小尺度的岩石断裂问题，其中DDA 方法在接触的数学计算、块体的积分、计算收敛性方面具有严密的数学基础。不同于采用显式时间积分计算的DEM，DDA 可视为基于Newmark 隐式积分和接触开闭迭代求解的隐式非连续计算方法，在离散块体系统大变形、大位移计算中可以保持时间步取值的无条件稳定。DDA 方法最初用于解决裂隙切割的块体系统的静态稳定性及动态响应问题。KE[17]提出采用子块体划分DDA 方法模拟裂隙岩体变形问题，焦玉勇等[18]则通过引入节理约束的三角形子块体剖分及黏结单元来模拟节理岩体的渐进破裂过程，甯尤军等[19]和ZHENG 等[20]相继对黏结单元的破坏准则、单元积分模型进行修正和优化。面对深部开采中与岩层破断、非连续大变形相关的岩石力学问题，DDA 方法具有天然的优势。许多学者利用DDA 方法研究了深部沿空掘巷过程[21]、煤层开采导致的地表沉降[22]、煤矿岩体大变形和开挖过程[23]、矿山陷落柱[24]等问题。

非连续类数值模拟方法已成为分析岩石断裂问题的重要手段，传统子块体DDA 方法基于边-边接触关系和法向、切向强度参数的修正模拟断裂过程，对复杂岩石材料变形破坏过程的适用性和计算精度有所欠缺。基于离散表征黏结单元力学破坏过程的朴素思想，笔者提出了一种分段式黏结模型，提升传统黏结模型在模拟复杂本构关系时的适用性，更准确地模拟岩石动静荷载下的断裂问题。

1 基于分段式黏结模型的DDA 方法

1.1 子块体剖分DDA 方法的基本原理

岩体可以看作是由节理、裂隙等不连续面切割形成的断续介质系统。考虑不连续面的在几何方面的有限延展性和在力学方面的非连续变形特性，采用考虑非连续界面约束的子块体几何剖分和黏结/接触单元进行力学近似求解，是常用的非连续计算方法[17]。子块体剖分后，黏结单元使用由岩石性质等效的虚拟节理参数，接触单元则根据真实节理性质进行赋值[18]。

1.1.1数学原理与控制方程

离散块体系统可通过连接或接触关系形成统一的分析对象。块体中任意一点p(x,y)的位移向量u={u,v}T可表示为u=Td，其中，T为位置矩阵；d为表征块体刚体位移和应变自由度变量的向量[12]：

式中，u0、v0为块体沿x轴和y轴的平动位移量；r0为块体在xoy平面的旋转角(逆时针为正)；εx、εy、γxy分别为x轴和y轴方向的应变和剪切应变。

p点的位置矩阵T可由式(2)计算[12]：

式中，x0、y0为块体的形心点坐标分量。

块体之间的接触/黏结关系可由多种方式进行求解，例如罚函数法、拉格朗日乘子法等。由于不增加额外的自由度，以及对于隐式求解中矩阵对阵正定性质的满足，罚函数法使用更为广泛。考虑所有块体之间的接触/黏结关系以后构造总势能泛函，基于最小势能原理可获得块体系统的动力平衡方程[20]：

式中，M为质量矩阵；C为阻尼矩阵；KE为弹性刚度矩阵；KB为黏结单元刚度矩阵；KC为块体接触矩阵；F为系统外力向量；D为块体的广义位移自由度向量，其上1 个点表示对时间的一阶导数，2 个点表示对时间的二阶导数。

1.1.2时间积分方案与隐式求解模式

求解动力学方程，Newmark 积分法是常用的一种积分方案。DDA 方法中原始的常加速度积分方案可视为Newmark 积分方法的一个隐式求解方案的特例[25]。从tn时刻到tn+1时刻，间隔为Δt的求解步中，常加速度积分方案可表示为

将式(5)代入动力平衡方程，可获得每一时步的位移增量Dn+1的求解公式：

由此可见，每一个计算步中需要求解一个整体矩阵，可以证明该矩阵满足对称正定特征[12]。

1.2 分段式块体黏结模型

基于可变形块体及黏结/接触的方法可实现岩石材料的连续变形及非连续破坏过程分析。常用的子块体DDA 方法中，每个子块体均为常应变弹性单元，采用考虑拉伸破坏的修正摩尔-库伦准则，破裂发生在子块体边界上，不考虑块体内部发生破坏。黏结模型可视为一种特殊的接触模型，可模拟黏结单元两侧块体的连续变形、黏结单元的损伤及失效。基于罚函数的处理方法，在子块体间引入黏结模型，可以同时模拟分析域的连续变形、渐进破裂、以及破裂后的动力学过程。随着对岩体及岩石材料精细力学模拟的发展，传统的黏结模型需要更新以更好的还原岩体结构/岩石材料的原始几何形态及物理状态，尤其是对于岩石裂纹面的非线性破裂过程，因此本研究提出分段式黏结单元模型。

1.2.1分段式黏结模型概述

基于离散化表征的思想，笔者提出了一种分段式单元黏结模型，以更准确的表征块体间黏结的逐步损伤失效过程。如图1 所示，黏结模型在几何上从属于2 个离散可变形块体的2 条重叠边，在力学上通过本构关系定义2 个边的相对位移与作用力的关系。几何方面，长度为lb的一对黏结边可离散化为长度均等的nb个线段，在每对线段的中心点沿着边的法向和切向可施加一对罚弹簧进行黏结对边的位移约束。分段式黏结模型主要包括以下参数：黏结分段数量nb，黏结长度lb，法向黏结刚度切向黏结刚度和可定义为

其中，Ec、h分别为块体的特征弹性模量和特征长度；αb为黏结刚度比，表征法向黏结刚度与块体弹性模量和特征长度的比例关系；βtn为切向刚度比，表征切向黏结刚度与法向黏结刚度的比值。基于分段策略，黏结单元内非线性应力分布、非线性变形-应力关系可以等效为分段线性的模式进行表征和求解，而每个分段的应力、变形由设置在分段形心点处的黏结弹簧模拟。

1.2.2黏结本构关系及分段式表征

边-边黏结单元的平均正应力σm可由式(9)获得：

边-边黏结单元的平均剪应力τm可由式(10)获得：

基于摩尔库伦强度准则判别黏结单元的法向拉伸破坏和切向剪切破坏。若黏结单元平均正应力超过单元抗拉强度值，黏结发生拉伸断裂，黏结单元失效，两侧块体由黏结状态变为纯接触状态。

式中，Tb为黏结单元的抗拉强度。

若黏结单元平均剪应力超过剪切极限值，则黏结单元发生剪切断裂，黏结单元失效，同时两侧块体由黏结状态变为纯接触状态。

式中，cb为边-边黏结单元的黏聚力；φb为边-边黏结单元的内摩擦角。

1.2.3分段式黏结模型的力学求解

黏结单元的本构模型采用罚函数方法计算，边-边黏结单元的法向约束和切向约束的能量泛函可表示为

式中，p0、q0为时步初始时刻分段si中线点的空间位置向量；ni和ti为当前分段的法向和切向单位向量；up和uq为当前时步p点和q点的位移增量向量，可由式u=Td求得。

以第si分段为例，根据能量泛函最小化原理，法向约束能量泛函对未知量D的二阶导数项和一阶导数项分别以子矩阵和子向量的形式加入总刚度矩阵和总广义力向量：

同理，切向约束能量泛函对未知量D的二阶导数项和一阶导数项也分别加入总刚度矩阵和总广义力向量：

黏结单元各分段按上述公式求解并叠加可得到整个边-边黏结单元模型在DDA 计算框架中的罚函数求解公式。

2 模型连续变形与应力校验

准确模拟材料在连续变形阶段的变形与应力是模拟岩石变形和破坏的基础。对于分段黏结模型，需要确定其几何参数(分段数目、黏结单元尺寸)、变形参数(黏结刚度比、切向刚度比)、强度参数(抗拉强度、黏聚力、内摩擦角)。

采用简支梁中心受压、单轴压缩和巴西劈裂这3个实验，通过对比模型的解析解与模拟结果，对分段式黏结模型的几何参数、变形参数进行校验并提出合适的取值范围或推荐值。

2.1 简支梁中心受压变形验证

本文选取简支梁模型如图2 所示，长度l=1 m，高度H=0.1 m，宽度取单位值b=1 m。简支梁模型的密度ρ=2 500 kg/m3，弹性模量E=50 GPa，泊松比μ=0.2，抗拉强度T=12 MPa，黏聚力c=63 MPa，内摩擦角φ=42°。在梁的左端约束其水平位移和垂直位移，右端只约束其垂直位移，分别模拟简支梁中的铰支座和滚动支座。在梁的中间处施加Fp=20 kN 的集中荷载，梁中轴线中心点处垂直变形w的解析解为

图2 不同块体尺寸的简支梁模型Fig.2 Illustrations of simplely supported beam models with different block sizes

其中，I为梁的截面惯矩，计算公式为

可求得梁中心点处垂直位移的理论值为-0.1 mm。

在本算例模拟中，考虑了4 种块体(黏结单元)尺寸(15.0、10.0、7.5、5.0 mm)以及6 种黏结刚度比(5、10、25、50、75、100)，黏结单元分段数目取2，切向黏结刚度比取1.0，设置梁的中轴线中心点为监测点。监测点得到的垂直位移如图3(a)所示，模拟结果与理论结果的相对误差如图3(b)所示，可以看出，随着黏结刚度比增加，监测点的位移逐渐接近理论值，且随着划分块体尺寸的增大，监测点的位移越来越接近理论值。当黏结刚度比大于25 时，所有模型的相对误差都达到了10%以下，达到可以接受的精度。另外，控制其他变量相同，取黏结单元分段数为2、3、5 时的计算结果表明结果相对偏差远小于1%，取分段数为2 即可满足要求。

图3 简支梁监测点的垂直位移与相对误差Fig.3 Simulated vertical displacement and relative error of the track point in simply supported beam example

并且当黏结刚度比比大于50 后，相对误差的变化趋于稳定，且削弱了块体尺寸的影响。故建议黏结刚度比应在[50,100]为宜。选取黏结刚度比为50，块体尺寸为7.5 mm 的简支梁模型，在梁中心处分别施加1Fp至5Fp的力，在线弹性阶段模型的相对误差保持在4.09%，不发生显著变化，表明未发生破坏阶段变形计算的稳定性。

2.2 试样单轴压缩应力验证

岩石单轴压缩试验是测定岩石强度的基本试验，其测定的单轴抗压强度在地下工程中应用广泛，具有重要的工程意义。笔者建立长100 mm、宽50 mm 的矩形作为标准试样单轴压缩二维模型，模型的物理力学参数与简支梁一致。如图4 所示，在模型上边界施加0.1 MPa 的均布载荷，下边界固定，模型中心处设置一监测点。对于均匀试样，在轴向应力应在垂直方向上均匀分布。

图4 不同块体尺寸单轴压缩模型Fig.4 Illustrations of uniaxial compression models with different block sizes

在本算例模拟中，考虑了4 种块体尺寸(10.0、7.5、5.0、2.5 mm)以及6 种黏结刚度比(5、10、25、50、75、100)，黏结单元分段数目取2，切向黏结刚度比取1.0，模型中心的应力如图5(a)所示，模拟结果与理论结果的相对误差如图5(b)所示。由图5 可知，所有模型的相对误差均低于2%，表明分段式黏结模型在连续变形阶段的应力分布精度足够高。且随着黏结刚度比的增加，中心点的轴向应力越接近理论值，当黏结刚度比大于25 时，相对误差低于0.5%。控制其他变量相同，取黏结单元分段数为2、3、5 时的计算结果表明结果相对偏差远小于1%，取分段数为2 即可满足要求。选取黏结刚度比为50，块体尺寸为5 mm 的单轴压缩模型，在模型上边界施加0.1～1.0 MPa 的分布力，在线弹性阶段模型的相对误差保持在0.13%，不发生显著变化，表明未发生破坏阶段应力计算的稳定性。

图5 单轴压缩监测点的轴向应力与相对误差Fig.5 Simulated stress and relative error of the track point in uniaxial compression test

结合简支梁中心受压与单轴压缩的模拟结果进行分析可知，在几何参数方面，块体尺寸、黏结单元数目对于算例中连续变形、及应力计算精度的影响不显著；在变形参数方面，当黏结刚度比≥50 时，可以保证变形与应力的相对误差都满足精度要求，同时减小了块体尺寸的影响，因此在后续模拟中将采用黏结刚度比50 进行实验；而剪切刚度比值设为1.0，可满足算例中对于计算精度的要求。对比原有黏结模型[18-19]，分段式黏结模型中黏结刚度比可以克服对于不同块体尺寸下取不同刚度的复杂取值问题，根据经验在[50,100]内可获得合适的精度要求。

2.3 巴西圆盘应力分布验证

巴西劈裂试验是一种广泛使用的间接测量岩石抗拉强度的方法。本例中采用平台巴西劈裂试验，巴西圆盘模型如图6 所示，半径R=0.025 m，厚度取单位值w=1 m，物理力学参数设置与简支梁一致，黏结刚度比取50，块体尺寸取1～2 mm。模拟时将圆盘的下边界固定，在上边界加载平台上假定试验机施加了一集中力F=10 kN，平台加载角2α=12°，则加载平台上的均布荷载P为

图6 巴西劈裂模型Fig.6 An illustration of the Brazilian splitting model

式中，α为弧度。

平台巴西劈裂中圆盘各个点的x方向应力与y方向应力的解析解参考黄耀光等[26]的计算公式。为验证应力计算的准确性，在圆盘中心沿垂直方向均匀设置9 个监测点。监测点x方向应力与y方向应力的模拟结果与解析解的对比如图7(a)所示，相对误差如图7(b)所示。由图7 可以看出，监测点应力的模拟结果与解析解基本吻合。y方向应力的相对误差均低于10%，x方向应力的相对误差表现为越远离端面精度越高，在靠近端面时会表现出较大的误差，由于数值求解时模型与边界条件与解析求解时并不完全一致，其结果符合圣维南原理。算例进一步验证了分段式黏结模型的准确性。

图7 巴西劈裂监测点应力与相对误差Fig.7 Stress measurement of track point in Brazilian splitting test and relative errors

3 砂岩准静态变形及强度特性

在天然的岩体中，总是存在着许多裂隙、节理和断层等不连续面，这些缺陷很大程度的影响了岩体的强度特性以及断裂模式，所以探究岩石在含预制缺陷时的强度特性与裂纹扩展规律对地下工程围岩破坏问题有重要价值。笔者选取煤矿开采等地下工程中最为常见的岩石种类之一——砂岩为研究对象，设置了2 组模拟实验，探究了缺陷对砂岩强度特性和变形的影响，表1 总结了2 组实验所使用到的参数。

表1 模拟参数Table 1 Simulation parameters

3.1 单轴压缩强度实验模拟

选择对杨圣奇等[27]的含孔洞的砂岩单轴压缩实验进行模拟，探究孔洞对单轴压缩强度特性以及裂纹萌生的影响。该砂岩试样致密均匀，平均密度为2 620 kg/m3，高度为120 mm，宽度为60 mm，厚度为30 mm。选择长度为120 mm 宽度为60 mm 的矩形截面作为单轴压缩模拟二维模型，分别建立孔洞直径d=0(完整岩样)、5、10、15 mm 的模型，在孔洞周围划分更小的块体以更好的捕捉裂纹萌生过程，如图8 所示。单轴压缩模拟过程中，模型的下边界固定，上边界以0.5 mm/s 的速度垂直向下加载。

图8 含孔洞岩样数值模型Fig.8 Numerical models of rock samples with holes

进行数值模拟实验前，需对模型的子块体弹性参数(弹性模量、泊松比)、黏结刚度参数(黏结刚度比，经验取值50；切向与法向刚度比，取值1.0)、黏结强度参数(抗拉强度、黏聚力、内摩擦角)进行校准，采用“试错法”对以上参数进行标定。

笔者通过进行大量的标准岩样单轴压缩和巴西劈裂试验模拟寻找参数规律，确定的标定步骤为：①通过模型的单轴压缩试验来校准弹性模量、泊松比，根据岩样的真实参数进行赋值并进行±5%范围内调整；②通过模型的巴西劈裂试验来校准模型的抗拉强度，赋予模型真实岩样的抗拉强度并进行±5%范围内调整；③通过单轴抗压强度来校验黏聚力与内摩擦角，对比强度和破坏模式确定最终取值。结合该实验完整岩样单轴压缩曲线，标定结果见表1，单轴压缩模拟的应力应变响应曲线与实验曲线的对比如图9 所示，可以看出2 条曲线峰值强度与总体趋势基本吻合。

图9 完整岩样实验与模拟应力应变曲线对比Fig.9 Comparisons of experimental and simulated stress-strain curves for uniaxial compression of intact rock sample

使用同一套参数进行含孔洞砂岩的模拟实验，得到的应力应变响应曲线与裂纹萌生状态如图10 所示。由图10 可知，在含孔洞模型中，拉伸裂纹最先在孔洞中心附近上下同时萌生，同时，强度特性表现为随孔洞直径的增加，单轴压缩峰值强度与峰值应变均出现明显的下降，这些现象均与实验室实验观察到的一致。对比图9、图10 可以看出，使用分段式黏结模型模拟可以很好的捕捉到试样的抗压强度与裂纹萌生点，最大的差别在于加载的初始阶段，这是由于实验室实验所使用的试样内部存在一些微裂隙，在位移加载初期表现为压实阶段，而模拟时是一个完整的连续体，并没有这个阶段。

图10 应力应变曲线对比Fig.10 Comparisons of stress-strain curves

3.2 巴西圆盘劈裂实验模拟

本算例探究砂岩巴西圆盘在单裂纹、多裂纹条件下的裂纹扩展模式，对ZHAO 等[28]的实验进行了数值模拟，使用的参数取值见表1。

单裂纹巴西圆盘的模型如图11(a)所示，圆盘的直径为62 mm，裂隙的长度为20 mm，宽度为1 mm。裂隙与水平方向分别呈0°、30°、60°、90°，块体尺寸为1～3 mm，在有裂纹的地方划分更小的子块体以更好地捕捉裂纹萌生和裂纹扩展的过程。将模型的下边界固定，在上边界中以0.5 mm/s 的速度垂直向下进行位移加载模拟试样巴西劈裂过程。由图11(b)、(c)可知，当裂纹倾角为0°时，裂纹从预制裂纹的中部上下同时萌生，并沿着垂直方向扩展，直至贯穿。对于其他3 种裂纹倾角模型，裂纹都是从预制裂纹的尖端开始萌生，并沿着加载方向继续扩展，直至完全贯穿。

图11 单裂纹巴西圆盘模型及其裂纹萌生和破坏模式Fig.11 Brazilian disk models with single crack and their failure mode

含有多条裂纹的巴西圆盘模型如图12(a)所示，圆盘的直径、裂纹的大小、子块体剖分策略、加载方式均与单裂纹模型一致，多裂纹构成的裂纹系统有如图①、②、③、④四种形式。由图12(b)、(c)可以看出，对于模型①和模型②，裂纹从预制裂纹L1 和L2 的尖端同时萌生，靠近试样端面的2 个尖端从裂纹萌生开始，不断沿着加载方向扩展直至贯穿，另外2 个尖端裂纹则不断沿着靠近模型中轴线的方向扩展。对于模型③，裂纹从预制裂纹L2 的上下尖端同时萌生，并沿着加载方向扩展直至贯通破坏，裂纹扩展模式与单裂纹巴西劈裂类似。对于模型④，裂纹最先从预制裂纹L1 和L3 的上下中部萌生，紧接着预制裂纹L2 的上下尖端也同时出现裂纹，并且逐渐向预制裂纹L1和L3 合并，最后一起上下贯穿。