FLAC3D中岩石能量耗散模型的开发与应用

2021-09-01宋子枫郑冬杰神文龙勾攀峰韦四江

煤炭学报 2021年8期

王猛，宋子枫，郑冬杰，神文龙，勾攀峰,韦四江,2

(1.河南理工大学能源科学与工程学院，河南焦作 454000; 2.中国矿业大学深部煤炭资源开采教育部重点实验室，徐州江苏 221116; 3.煤炭安全生产河南省协同创新中心，河南焦作 454000)

准确描述岩石强度与变形是评价工程岩体安全和稳定性的基础，现有应力强度理论与破坏准则很难有效描述岩石复杂的强度变化与整体破坏行为[1]。岩石破坏是能量驱动下的一种状态失稳现象[2]，受载岩石变形过程伴随能量演化，峰前以能量积聚为主，峰后主要表现为能量耗散与释放，岩石变形破坏可看成是不同形式能量相互转换的结果。因此，从能量角度研究岩石变形和失稳，更接近其破坏本质[3-5]。

基于能量守恒理论，国内外学者对岩石加卸载过程能量演化规律开展了大量研究，取得了丰富的研究成果。具有代表性的有:谢和平等[6-8]阐述了岩石变形破坏过程能量耗散与损伤、能量释放与整体破坏的关系，定义了基于能量耗散与释放原理的岩石强度失效准则和岩体整体破坏准则。张志镇等[9-11]研究了不同岩石类型能量演化的异同点，建立了岩石能量演化随加载条件的自我抑制模型，揭示了岩石能量耗散的围压效应。宫凤强、李淼等[12-13]分别研究了动静态劈裂过程岩石能量构成与耗散特征，揭示了动静荷载作用下岩石能量耗散机制及各向异性特征。与此同时，许江等[14]、杨磊等[15]、尹光志和苏国韶[16-17]、李杨杨等[18]分别研究了含孔隙水砂岩、煤岩组合体、不同加载速率以及循环加卸载等条件下岩石能量积聚与耗散特征，丰富了岩石能量理论的研究成果。

现有研究常聚焦于室内试验，多是针对特定加卸载条件、岩样类型以及其他条件下(含水率、温度[19]等)岩石能量演化特征的描述，研究对象着重关注岩石屈服、峰值等特征点，对于贯穿岩石变形破坏全过程的能量实时演化及分配规律认识不足[9]。岩石能量积聚、耗散及其状态间的转化决定了其内部裂隙扩展和变形破坏，虽然从能量角度反演岩石破坏路径对实验结果意义不大，但对于评估工程岩体区域稳定、精准定位加固位置至关重要。目前，实验室常用岩石破坏监测手段有声发射定位[20],CT扫描[21],SEM电镜扫描[22],CCD相机[23]等，其中，仅有声发射可以实时定位受载岩样变形破坏路径，但该方法多用于单轴压缩试验，且存在操作复杂、误差大等缺点，很难推广应用于工程岩体的监测。分析工程岩体稳定性时多是基于围岩应力场、位移场及支护场的耦合作用关系进行综合评判。

为增加能量理论模型分析工程问题时的适用性，笔者基于能量守恒和有限差分理论，推导岩石能量耗散有限差分方程式，采用FISH语言将其应用于FLAC3D应变软化本构模型。以室内试验结果验证采用耗散能演化描述受载岩石变形破坏的合理性，进一步讨论能量耗散模型在深部高应力巷道中的应用，揭示工程岩体峰后变形破坏全过程耗散能演化特征与分配规律，为从能量耗散角度评估工程岩体安全和稳定性提供理论依据，弥补FLAC3D软件采用塑性区定性判断地下工程围岩破坏时的不足。

1 岩石能量耗散的室内实验

1.1 岩石单、三轴压缩试验

试验岩样取自徐矿集团三河尖煤矿吴庄区运输大巷，选取岩性均匀、结构完整的取芯岩柱加工岩样，采用GCTS-RTX3000岩石力学试验系统进行单轴和常规三轴压缩试验，测试岩石基本力学参数，围压等级设置为0,5,15和25 MPa，测试结果分别如图1和表1所示。

表1 常规单、三轴试验结果Table 1 Results of uniaxial and triaxial tests

图1 岩石应力-应变曲线Fig.1 Stress-strain curves of the rock samples

1.2 岩石能量耗散特征

岩样在试验机外力作用下产生变形，假设该物理过程与外界没有热交换，即封闭系统，外力功所产生的总输入能量为Wz，根据热力学第一定律[6]得到

Wz=We+Wd

(1)

式中，We为岩石可释放弹性应变能;Wd为耗散能。

由式(1)可知，岩石总输入能量Wz一部分作为可释放的弹性能储存在其内部，另一部分则伴随岩石变形破坏耗散和释放，此过程不可逆[8-9]。岩石损伤程度越高，耗散能越大，残余弹性能越少。

对于单轴压缩试验，试样总输入能量Wz等于轴向输入能量Wa;而常规三轴压缩试验，轴向输入能量Wa应包含试样总输入能量Wz以及环向膨胀对液压油做功所释放能量Wh两部分，为此，常三轴压缩条件下岩石总输入能量[6]可表示为

(2)

式中，ε1，ε3分别为轴向应变和环向应变;σ1，σ3分别为轴向应力和围压。

根据弹性力学相关理论，岩石可释放弹性能We[6]可由下式计算:

(3)

则岩石耗散能Wd可由下式求得

Wd=Wz-We

(4)

采用式(2)～(4)分别计算岩样总输入能量Wz、弹性能We和耗散能Wd，如图2所示。由图2可知，岩样总输入能量随轴向应变增加逐渐增大，应力达到峰值前主要以弹性能形式储存在岩样内部，耗散能变化不大。当岩样加载至峰后阶段，岩样总输入能量依然持续增加，但主要伴随峰后破坏而耗散，此外，岩样内部储存的大部分弹性能也将转化成耗散能释放，能量转化程度与岩样变形破坏程度正相关。

图2 室内试验岩样能量演化规律Fig.2 Evolution law of energy of rock samples in laboratory test

2 FLAC3D能量耗散模型的开发

σ1<σ2<σ3

(5)

相应的主应变增量Δεi[25]可分解为

(6)

单元体弹性应变增量和应力增量之间的关系[25]可表示为

(7)

(8)

式中，α1,α2为切变模量G和体积模量K控制的岩石材料常数，其中，α1=K+4G/3，α2=K-2G/3。

基于Mohr-Coulomb的应变软化模型在峰后软化过程控制强度参数描述材料的峰后软化行为，对于剪切破坏，非关联流动法则[25]改写为

(9)

式中，λs为塑乘因子;gs为剪切势函数，其表达式为

(10)

式中，ψ为剪胀角。

因此，联立式(6)，(7)和(9)，由总应变表示的应力增量Δσi表达式[25]为

(12)

(13)

式中，φ为内摩擦角。

根据FLAC3D应变软化模型屈服准则，岩体材料剪切破坏面上应力[25]满足:

(14)

式中，c为黏聚力。

拉伸破坏面上应力[25]满足:

ft=σ3-σt=0

(15)

式中，σt为抗拉强度。

式(14)和(15)分别为FLAC3D应变软化模型判定材料压剪破坏和张拉破坏的失稳准则。

图3 单元能量计算简图Fig.3 Schematic diagram for the energy calculation

(16)

假设t时刻对应第n循环，t+Δt时刻对应第n+1循环，则第n+1循环内总能量增量可表示为

(17)

第n循环至第n+1循环运算的物理意义表示岩体受外力作用时，某一运算时步内单元网格应力-应变行为的差分运算。t+Δt时刻单元总能量等于n+1循环总能量增量累加，即

(18)

不考虑峰后弹性模量衰减，单元可释放弹性应变能可由下式计算:

(19)

则单元耗散能计算式如下:

(20)

图4 单元能量计算流程Fig.4 Energy calculation flow for the model cells

3 模型校验与运算

3.1 岩样模型建立

建立标准岩样数值模型，尺寸φ×H=50 mm×100 mm，共划分17 280单元，如图5所示。模型底部固定位移边界，四周施加环向边界应力模拟围压，上边界通过施加轴向位移对岩样加压，加载速率2.5×10-5mm/step，计算时忽略岩石自重影响。

模型采用基于Mohr-Coulomb的应变软化本构模型，以表1给出的岩样力学参数对模型进行初始赋值。岩石峰后软化阶段黏聚力和内摩擦角随应变的变化数据由试验结果计算得到，见表2。

表2 岩石峰后软化参数Table 2 Post-peak parameters of rock samples

依据图5建立的数值模型，表1和2所列岩石力学参数进行数值分析，图6给出了模拟岩样应力应变曲线与试验结果的对比，由图6可知，不同围压下岩样应力应变曲线模拟结果较为匹配试验结果，单轴压缩和低围压条件下岩样脆性特征得到较好描述，表明所选取的岩石参数和数值本构模型合理，可用于后续能量模块的模拟分析。

图5 岩样数值模型Fig.5 Numerical model of rock samples

图6 应力-应变曲线校验Fig.6 Verifications of the stress-strain curves

3.2 岩样能量模型校验

考虑到耗散能是由岩石总能量减去弹性应变能计算得到，限于篇幅，校验能量模型时只校验总能量和耗散能两个指标。如图7(a)所示，不同围压下岩样总能量数值模拟结果与试验结果基本吻合，表明FISH语言编写的能量算法可靠性较强。

但需要指出的是，模拟岩样峰后耗散能演化曲线与试验结果吻合度较好，但峰前存在差异性，如图7(b)所示。峰前阶段，岩样试验耗散能主要呈缓慢线性增加;加载至峰值80%时，耗散能近似指数增加;一旦到达峰后某破坏点，耗散能演化为斜率较大的直线增长，最终进入残余阶段。数值模拟可近似反演耗散能演化趋势，耗散能在峰后增长曲线与试验结果匹配较好，且总量大致相等;但峰前阶段存在显著差异，模拟峰前耗散能近似等于0，主要原因是能量模型峰前假设不发生塑性变形，无法准确描述孔隙、松软岩体初始压密阶段和峰前塑性阶段，导致围压越大，峰前耗散能相差越大。但是，对于地下工程岩体的分析更关注于峰后变形破坏，由于该模型可较好模拟岩石峰后能量耗散特征，可用于岩石峰后变形破坏的模拟。

图7 能量模型校验Fig.7 Verifications of the energy model

3.3 岩样峰后破坏路径反演

图8给出了不同围压下岩石破坏模式对比，其中，数值模拟给出的是岩石残余阶段耗散能密度演化云图，对应试验岩样残余阶段破坏照片。由图8可知，岩样单轴压缩试验呈纵向劈裂破坏，模拟岩样耗散能集中区域位于岩样中线偏上位置，并向周边演化出一“八”字分支集中线，岩样耗散能演化沿中线呈轴对称分布，反映了岩样单轴压缩时的纵向劈裂特征，与试验结果吻合度较高。随着围压增加，试验岩样逐渐向剪切破坏演化，破断角逐渐减小;模拟岩样在5 MPa时耗散能密度集中成“y”型分布，长边集中程度显著大于短边，为岩样主破断面;随围压进一步增加，其破坏模式演化成斜切面，且围压越大，耗散能越大，破断角越小，与试验结果基本吻合，表明采用耗散能密度演化揭示岩石峰后破坏行为是可行的。

图8 岩石破坏模式校验Fig.8 Verifications of the failure modes of rock samples

工程现场，准确掌握围岩变形破坏路径，精准定位加固位置是控制其稳定的前提。为此，反演岩石峰后耗散能演化过程可为揭示其主破坏区演化、评判岩石稳定性提供新途径。图9给出了围压25 MPa时不同峰后目标点岩样耗散能密度演化云图。岩样初始进入峰后阶段(B点)，经历裂隙孕育形成初始主裂隙面，并随着轴向应变增加，主裂隙面耗散能逐渐增加，一旦岩样进入残余阶段(F点)，主裂隙面耗散能密度增幅减小并趋于稳定，表明岩样此时已发生宏观破断，继续加载将导致其结构失稳。

图9 围压25 MPa岩样峰后破坏路径Fig.9 Post-peak destruction path of rock sample with a confining pressure of 25 MPa

4 能量耗散模型的应用

将能量耗散模型应用于分析吴庄区运输大巷，尝试从能量耗散角度揭示巷道变形破坏路径，为巷道加固提供依据。建立试验巷道三维数值模型，尺寸X×Y×Z=60 m×40 m×60 m，水平、底边界限制位移，上边界施加20 MPa载荷模拟巷道800 m埋深，依据地应力测试结果，取侧压系数0.8。巷道断面宽×高=5 m×4 m，模拟巷道锚杆支护参数与现场一致。

模拟采用应变软化模型，由于岩石实验参数无法直接应用于工程岩体，模拟前需校验岩体参数。基于实验获取的岩石强度参数衰减规律，采用RocLab软件进行岩体参数转换，以获取的初始岩体参数进行赋参，以实测巷道变形作为已知值，校验获取的岩体参数见表3,4。同时，为防止单元畸变中断运算，模拟时未设置set large命令。由于巷道周边围岩受力多处于三向不等压状态，计算单元能量时考虑中间主应力，各应力分量增量计算方法见式(12)，再由式(17)～(20)分别计算单元总能量、弹性能和耗散能。

表3 校验岩体参数(峰前)Table 3 Calibrated parameters of the rock masses

图10给出了试验巷道模拟结果与现场对比图。试验巷道开挖后，采用锚网索支护，巷道顶帮变形得到控制，顶帮变形均在200 mm以内;但底臌严重，最大底臌量接近1 000 mm，模拟巷道变形与现场实测数据基本吻合。

表4 校验岩体参数(峰后)Table 4 Calibrated parameters of the rock masses

图10 试验巷道变形破坏对比Fig.10 Comparisons on the roadway deformations

模拟巷道塑性区主要集中于两帮和底板，基于软件拉剪屈服破坏准则，仅可以区分底帮围岩变形破坏方式，但很难量化底帮围岩的破坏程度，也就无法准确定位巷道后期加固区域。而利用开发的能量耗散模型，从巷道耗散能集中区分布可实现对巷道围岩主失稳区域的实时定位，同时，依据围岩耗散能集中程度可以量化岩体变形破坏程度。如图10所示，巷道底板破坏后，底角位置分别存在两个耗散能集中区，均呈条状分布，一条贯通与底角，诱发底角剪切破坏;另一条则向底板中部演化，加剧底臌。

基于3.3节，研究开发的能量耗散模型除可定位围岩主破坏位置，还可以通过控制运算时步反演巷道变形破坏路径。图11给出了不同时步对应的巷道围岩耗散能密度演化云图，运算平衡时对应的耗散能云图如图10所示。由图11可知，巷道开挖后耗散能集中区位于巷道两底角处，表明直角相比于拱形更易产生剪切破坏，这与前人研究结论相吻合。随着时间增加，巷道底角破坏后，耗散能集中区域将向深部围岩转移，并逐渐演化成带状分布，一端向巷道两帮转移，一端朝向巷道底板方向，诱发巷道底臌和两帮收敛。对于破坏后的巷道，如对底帮尤其是底角围岩不加以控制，剧烈底臌将进一步引发巷帮松动破坏，甚至危及顶板稳定，最终诱发巷道整体灾变失稳。