含初始损伤冻结砂岩动态破坏机理及损伤本构关系

王 磊，王远鹏，秦 越，苏宏明

(西安科技大学 建筑与土木工程学院，西安 710054)

1 研究背景

西部矿区侏罗系煤层上覆巨厚白垩系基岩(200～400 m，主要由泥岩、砂岩和泥砂互层构成)，具有弱胶结、低强度、易水解等特性[1]，在多应力场耦合作用下内部随机分布多尺度损伤缺陷，包括节理、裂隙等宏观缺陷和孔隙、裂纹等细观缺陷，属于复合损伤材料[2-3]。为达到止水、提高软弱围岩强度的目的，矿井建设常采用冻结法施工，冻结壁在井筒掘砌中会承受爆破、机械凿岩等动载作用，引起富水冻结砂岩初始宏细观缺陷与新生裂隙不断耦合累积直至岩石失效破坏，造成立井渗水、涌水等事故频发[4]。因此，研究白垩系富水冻结砂岩动态压缩破坏机理对西部地区矿井建设及寒区支护工程优化具有重要意义。

外载作用下，相对于岩石整体损伤的动态演变，其内部按尺度划分的宏细观损伤是相互联结的统一体，在多应力场环境中必然存在耦合现象。现有研究成果忽略两者的内在联系或将其分割开单独分析，无法知悉岩石宏细观损伤的耦合演变机制。赵光明等[5]基于软岩的应变硬化与塑形流动特性，以损伤体取代 ZWT 模型中弹簧元件。袁璞等[6]考虑干湿循环和动载耦合作用对砂岩的劣化，得到相应损伤演化方程。Wang等[7]改进现有黏弹性本构模型，考虑材料的应变率效应和温度效应。上述本构模型均未考虑岩石的初始损伤，仅研究冲击过程产生的损伤。张力民等[8]考虑节理岩体的宏细观缺陷，基于TCK模型和复合损伤变量建立岩石动态损伤本构模型，对单轴压缩下节理岩体力学特性变化规律进行描述。刘红岩等[9]构建单轴加载下非贯通节理岩体动态损伤本构模型，其复合损伤变量考虑岩石的初始宏细观缺陷。赵航等[10]分别从弹性波波幅和统计强度理论定义裂隙岩体的宏细观损伤变量，通过连续损伤理论建立了宏细观裂隙耦合的裂隙岩体损伤本构模型。李克钢等[11]优化裂隙岩体损伤本构模型，定义损伤变量包含初始宏细观耦合缺陷，表示岩石损伤场-应力场耦合作用机制。陈蕴生等[12]用奇异损伤变量和分布损伤变量表征轴向加载下孔隙和微裂隙两种损伤，运用CT技术研究细观损伤演化规律。孙传猛等[13]通过数字图像处理技术Monte-Carlo理论计算煤岩初始损伤值，构建模型反映煤岩单轴受压破坏全过程。上述模型虽考虑岩石初始损伤，但研究方向集中在静力学领域，对于初始损伤在动载作用下变化规律的研究相对较少。崔新壮等[14]、王倩倩等[15]考虑水泥砂浆材料含有的初始损伤，建立了动载作用下的损伤演化模型，但未考虑初始损伤与动态损伤的耦合作用，仅视为简单的累加关系，无法反映其真实特性。

基于上述研究，本文首先对冻结砂岩破裂损伤演化机理和强度与变形的应变率效应进行分析。其次，根据岩石损伤Lemaitre假设，对单轴冲击下考虑初始宏细观缺陷的损伤变量演变机制进行详细阐述，并引入影响因子δ表征岩样完整程度，再基于应变强度准则和双参数Weibull分布，定义动态损伤变量。最后，根据单轴压缩下损伤演化的正交异性，构建考虑初始宏细观损伤的冻结砂岩动态损伤本构模型，并基于实测数据确定本构模型参数。

2 动态冲击试验

2.1 试样制备

选取甘肃省五举煤矿冻结凿井穿越的白垩系富水弱胶结砂岩为研究对象，岩样呈红褐色，表面岩性均匀。原岩经钻、切、磨等工序制备成标准试件(D×H=50 mm×25 mm)，控制端面及轴向不平行度(±0.02 mm)，借助烘箱和超声波测速仪筛选试件以降低物理力学性质的离散性，其基本物理力学参数如表1所示。

表1 砂岩物理力学参数

为实现充分饱水，先将试件置于高温箱中以105 ℃恒温干燥24 h，然后置于真空饱和装置中以恒定0.1 MPa负压先干抽6 h后湿抽4 h，最后加水持续浸泡48 h以上。结合煤矿立井安全监测冻结壁温度场分布情况，设置试件温度为-30 ℃，饱水试件密涂凡士林并用保鲜膜包裹后放于低温控温箱中，以0.02 ℃/min恒定降温速率缓慢降温，且达到设定温度后继续稳定48 h以上，通过改变冲击气压对试件实现6种加载应变率。

2.2 分离式霍普金森压杆(SHPB)系统及工作原理

试验所用SHPB系统结构组成如图1所示，装置由动力源、杆系、信号采集系统、数据处理系统及缓冲阻尼共五部分构成。

图1 SHPB系统结构组成Fig.1 Structure of SHPB system

冲击试验时，冻结试件夹持在入射杆和透射杆之间，气室释放压力迫使子弹高速撞击入射杆，杆中产生应力脉冲并向前传播，接触试件后，应力脉冲在试件两端面-压杆界面多次反射、透射，分别被应变片1与2记录收集至超动态应变仪，经三波法公式计算可得试样动态力学参数，即

(1)

3 岩石动态损伤力学机理分析

3.1 冻结砂岩破裂损伤演化机理

根据公式(1)计算-30 ℃冻结砂岩不同应变率下动态应力-应变曲线和试件典型破坏形态如图2所示，由图可知，岩石破裂损伤呈4阶段变化：

(1)线弹性变形AB段：初始孔隙闭合完成，应变增幅远大于应力，损伤主要来源于无损部分的应变累积，动弹性模量随加载应变率的上升逐渐变大。

(2)微裂纹演化BC段：裂隙冰与岩石基质相互作用，应变增长迅速而应力增加缓慢，岩石内部微裂隙开始扩展、贯通，并衍生新裂纹，占据试样变形主体，整体塑性变形增强，损伤程度加剧。

(3)裂隙非稳定扩展CD段：位于曲线峰值两侧区域，岩石内部微裂纹非稳定扩展，沿主应力方向桥接、贯通，形成宏观主控裂隙。

(4)应变软化DE段：内部裂纹呈非线性扩展，岩石软化屈服，迅速失去承载能力并失稳破坏。

图2 动态应力-应变曲线Fig.2 Dynamic stress-strain curves

综上所述，冲击荷载作用下岩石破裂损伤由三部分构成：①初始细观损伤的扩张、加剧；②初始宏观缺陷的扩展、贯通；③无损部分裂隙的衍生、发育。由此可知，岩石初始损伤不可忽略，动态冲击荷载下岩石破坏是已有缺陷与新生裂隙之间的相互演变，属于损伤渐进软化的过程。

3.2 应变率效应

为进一步表征动载作用下冻结砂岩的强度和变形特性，引入动态强度增长因子DIF和动态应变增长因子DEIF概念，计算公式如下：

(2)

(3)

式中：σd、εd分别为动态抗压强度及对应应变；σs、εs分别为静态抗压强度和静态峰值应变。

不同应变率下-30 ℃冻结砂岩的DIF和DEIF值变化情况如图3所示。

图3 DIF、DEIF与应变率关系Fig.3 Relations of DIF and DEIF versus strain rate

由图3拟合关系式可知，DIF与DEIF值均随应变率的升高呈线性增长，应变率强化效应显著，其中动态压缩强度低于静态抗压强度，峰值应变率相关性则相反。究其原因，冲击载荷作用下冻结砂岩微裂纹的衍生与汇集、裂纹区的桥接与贯通均伴随能量耗散，且应变率较高时，试样呈多条主裂隙扩展，砂岩破碎程度加剧，承载能力下降而变形水平上升，需要吸收更多的能量。

4 动态损伤演变机制

白垩系砂岩是多种造岩矿物经漫长地质演变形成的矿物集合体，处于地层深部复杂多应力场环境时存在无裂损伤区，且损伤区易形成应力集中，在冻结凿井施工中常诱发冻结壁失稳破坏、井壁渗水等事故。由于宏细观损伤的动力学响应存在显著差异，仅用单一损伤变量表示整体损伤演化显然不符合实际情况。因此，有必要基于初始宏细观缺陷对冻结砂岩动态损伤演变机制展开研究。

4.1 考虑初始宏细观损伤变量

西部白垩系弱胶结砂岩内含随机分布的孔隙、裂纹等微缺陷，统计损伤理论认为岩石可视作由大量具有一定强度的规则微元弹性体构成的集合。设岩石含有的总微元数目为N，由3部分构成：初始细观损伤N1、初始宏观损伤N2、无损微元体N3，即N=N1+N2+N3，则岩石微单元结构用图4表示。

图4 岩石微单元结构Fig.4 Micro unit structure of rock

定义初始损伤变量为

(4)

式中D1、D2分别为初始宏观损伤变量和初始细观损伤变量。

无外载作用时，初始损伤构成简单，可近似视作简单数学累加，则初始损伤变量D12为

D12=D1+D2。

(5)

鉴于岩石内部宏细观损伤难以准确定量区分，为简化计算，引入初始损伤影响因子δ以表征岩样完整程度，其值与初始损伤程度成反比，即

δ=1-D12。

(6)

4.2 基于Weibull分布的动态冲击损伤

外载作用下，微元强度服从随机分布规律，当微元体承受应力超过极限时，即发生破坏。已有研究表明[16]，岩石材料在冲击载荷下损伤来源于局部微元体的非线性渐进破坏，微元弹性体强度服从双参数 Weibull 统计分布，相应的概率密度函数为

(7)

式中：P(ε)表示岩石微元强度概率分布函数；ε表示分布变量，基于岩石应变强度理论，ε也视作应变；m、α为反映岩石材料力学性质的 Weibull 分布参数。

进一步，由冲击载荷作用下破坏的微元体数目可定义统计损伤变量，即

(8)

式中：D为统计损伤变量；ΔN为已破坏的微元数目。

在任意应变区间[ε，ε+dε]内产生破坏的微元数目为NP(x)dx，当岩石材料达到某一应变ε时，有

(9)

联立式(8)与式(9)，由微元体破坏概率定义得动态冲击损伤变量演化方程为

(10)

4.3 考虑初始宏细观损伤的动态耦合损伤变量

冻结砂岩在经历冲击载荷作用时，新生动态损伤不能视作在原有损伤的基础上简单累加。一方面，初始宏观裂隙与细观孔洞在原有基础上扩展贯通，损伤进一步加剧；另一方面，新生裂隙与原有缺陷发生交互作用，不同尺度损伤缺陷的演化定量表示为损伤变量的耦合[17]，最终形成动态总损伤。

根据岩石损伤理论和Lemaitre假设可知，具有初始宏细观耦合损伤岩石的弹性模量E12可等效为

E12=E0(1-D12) 。

(11)

式中：E0表示无损完整岩石的动态弹性模量；E12为考虑初始宏细观损伤的岩石初始动态弹性模量。

动态冲击损伤在初始损伤的基础上进一步发展形成，则动态损伤状态下岩石的弹性模量E123为

E123=E12(1-D3) 。

(12)

代入式(11)整理得

E123=E0(1-D12-D3-D12D3) 。

(13)

反推可知岩石在冲击载荷作用下动态损伤变量D123为

D123=D12+D3-D12D3。

(14)

代入式(6)、式(10)可得

(15)

式(15)即为考虑初始损伤的岩石正交异性损伤演化方程。

5 动态损伤本构模型

5.1 正交异性损伤本构方程

假设外载作用下损伤演化各向同性，根据连续损伤力学基本关系式有

σ=Eε(1-D) 。

(16)

由于岩石在单轴冲击压缩作用下裂纹、孔隙扩展与受力方向近似平行，即损伤演化表现出显著的正交异性。因此，基于各向同性方程修正得到砂岩动态正交异性损伤本构关系为

σ=Eε(1-D)2。

(17)

代入式(15)可得考虑初始损伤的岩石受荷全过程应力-应变本构方程为

(18)

同理，考虑初始宏细观损伤岩石的弹性模量为

E12=E123(1-D12)2=Eδ2。

(19)

经历动载作用后岩石的弹性模量为

(20)

5.2 Weibull分布参数解析解

由动态应力-应变曲线原点至峰值点段可得下列几何条件[18]：

(21)

对式(18)关于应变ε求导得

(22)

显然式(21)中条件①、②分别满足式(18)、式(22)，需要通过条件③、④求解分布参数m、α。

由条件③与式(18)联立可得

(23)

两边取对数并整理得

(24)

由条件④与式(22)联立可得

(25)

整理可得

(26)

联立式(24)、式(26)可得：

(27)

联立式(22)与式(5)可得

(28)

令Ed为应力-应变曲线峰值点处的割线模量，则

(29)

将式(29)代入式(27)得

(30)

则式(15)可进一步化简为

(31)

定义应力-应变曲线峰值点对应的损伤值为临界损伤Dcr，其表达式为

(32)

由此可知，临界损伤值与岩石的初始损伤、弹性模量、动态抗压强度和峰值应变有关，是同时具备四个指标特性的综合量。

5.3 冻结砂岩动态损伤本构模型及参数确定

综上所述，考虑初始宏细观损伤的冻结砂岩动态损伤本构关系为

(33)

式中E为冲击载荷下岩石动态弹性模量，取应力-应变曲线初始弹性阶段斜率。

上述考虑初始宏细观损伤的岩石动态损伤本构关系式需要确定E、δ、ε、εd和m共5个参数。根据动态应力-应变曲线可知：E表示砂岩初始弹性模量，可取应力-应变曲线线弹性阶段斜率；参考相关文献[19-20]，选取初始损伤影响因子δ=0.95；应变ε和εd用实测数据代入；Weibull分布参数m通过线性回归法拟合确定。冻结砂岩动力学基本参数如表2所示。

表2 冻结砂岩动力学基本参数

利用冲击试验获取的动态应力-应变曲线拟合确定 Weibull 分布参数m值，根据拟合曲线变化形态和决定系数R2判断拟合效果，不同加载应变率下拟合曲线如图5所示。

图5 试验数据与拟合曲线对比Fig.5 Comparison between test data and fitting curves

由图5可知：曲线整体拟合效果好，决定系数R2最大为0.990，最小为0.915，表明理论结果与试验曲线吻合良好，能较好反映峰前段强度特征。随着应变率升高，冻结砂岩动态应力-应变曲线的线弹性特征越突出，应变率效应显著，本构模型对裂隙非稳定扩展阶段的描述更加符合实际情况，反映冻结砂岩的峰值强度和变形特性，拟合确定的本构模型参数如表3所示。

表3 本构模型参数

6 结 论

(1)冻结砂岩动态损伤经历“线弹性变形、微裂纹演化、裂隙非稳定扩展、应变软化”4个阶段，损伤由初始细观损伤扩张、初始宏观缺陷贯通、无损部分裂隙衍生以及初始宏细观缺陷与新生裂隙的交互演变，动态增长因子DIF和DEIF变化表明砂岩强度和变形的应变率强化效应显著。

(2)初始损伤由表示宏细观缺陷的变量累加构成，初始损伤影响因子δ反映岩样完整程度，不同尺度损伤缺陷耦合形成动态总损伤。由应变强度准则和双参数 Weibull 统计分布确定冲击损伤D3，由岩石损伤理论和 Lemaitre 假设反推得到总损伤D123，临界损伤Dcr是具备初始损伤、动弹性模量、动态抗压强度及峰值应变等指标特性的综合量。

(3)单轴冲击下岩石裂纹、孔隙扩展近似平行受力方向，损伤演化呈正交异性，基于各向同性损伤方程修正得到砂岩动态损伤本构关系，求解Weibull 分布参数解析解得到考虑初始宏细观损伤的冻结砂岩动态损伤本构模型。线性回归法拟合确定 Weibull 分布参数m值，决定系数均>0.910，理论模型能较好反映峰前段强度特征。