高应力区岩石统计损伤本构模型研究

贾 逸，魏良帅，黄安邦，和 铭，黄细超，蓝康文

(1. 中国地质科学院探矿工艺研究所，四川 成都 611734；2. 中国地质调查局地质灾害防治技术中心，四川 成都 611734；3. 成都理工大学地质灾害防治与地质环境保护国家重点实验室，四川 成都 610059；4. 西南交通大学地球科学与环境工程学院，四川 成都 610031)

《水利水电工程地质勘察规范》[1]定义高应力区是指最大主应力超过20 MPa且强度应力比介于2～4的高地应力区域，深部巷道、隧洞多处于高应力区[2]。高应力区岩体具有成岩作用复杂、非均质特性明显和所处应力水平高等特点，往往处于较为复杂的应力状态，为大型地下工程、水电工程的安全施工与建设带来严峻挑战[2]。自从COOK[3]于1965年提出岩石应力-应变关系的概念，人们对于岩石破裂有了更深认识，国内外学者对岩石本构关系进行深入研究，构建大量岩石本构模型，但绝大多数难以阐述高应力区岩体变形力学特性。深部岩体受上覆岩层荷载和高应力水平影响，其节理、裂隙等初始缺陷更是复杂，岩体损伤理论将岩体内部天然缺陷理想地简化为岩体的损伤，这与岩体特性较吻合[4-5]。KRAJCINOVIC等[6]基于岩石材料内部缺陷随机分布，有机地结合了连续损伤理论与统计强度理论，据此提出统计损伤的概念。

目前，关于岩石统计损伤本构模型已有很多文献进行讨论分析。曹文贵等[7-10]、刘成学等[11]、韦立德等[12]、蒋维等[13]、石崇等[14]、ZHU等[15]基于LEMAITREL[16-17]应变等价性假说，假定岩石微元强度服从不同概率密度函数(如WEIBULL分布[8-10,12,14,17]、幂函数分布[11]、正态分布[7]等)，并根据是否考虑损伤阈值、节理裂隙尺寸和力学参数变化等情况构建了大量统计损伤本构模型。这类模型对于一般工程岩体变形力学特性虽有较好辨识效果，但对于复杂应力状态下的高应力区岩体未必适用，究其原因，目前基于传统强度准则所建立的岩石微元强度度量方法对于高应力区岩体可能存在一定局限。曹文贵等[7-11]对岩石微元强度度量方法进行一些改进，弥补了岩石微元强度的度量方法的一些不足，但现阶段针对高应力区岩石统计损伤模型的研究较少，未有得到广泛认可的系统性统计损伤本构模型。

已有研究中，确定微元强度的度量方式多是通过传统MOHR-COULOMB(M-C)准则和DRUCKER-PRAGER(D-P)准则等而建立。基于D-P准则确定微元强度的合理性可能受该准则自身保守性的限制[14]；传统M-C准则的MOHR强度包络线近似直线，这样的近似方法在应力变化范围较小的浅部低应力区产生的误差不大，得到广大水利水电、岩土行业学者认可。然而在高应力条件下，应力变化范围很大，岩石在高应力状态下其MOHR强度包络线将远远偏离直线，将MOHR强度包络线近似为直线可能会造成较大误差甚至偏于危险[18]。故作者引入一种双参数抛物线型MOHR强度准则[19]建立微元强度度量方式，假定由该度量方式得到的微元强度服从幂函数概率密度分布，由此构建能反映高应力区岩石变形破坏全过程的统计损伤本构模型，此外对高应力区岩体损伤累积的扩展过程进行详细分析。

1 岩石强度准则

传统M-C准则是目前最常用的岩石强度准则之一，但传统M-C准则的MOHR强度包络线为线性，而高应力水平下的岩体处于复杂的非线性应力状态下，应力变化范围很大，传统的近线性MOHR强度包络线不适用于高应力区岩体。已有研究中，MOHR强度包络线的形式有双曲线、抛物线和摆线等形式，本文引入一种双参数二次抛物线型MOHR破坏准则(图1)，考虑到剪应力的对称性，公式如下[19]：

σ=aτ2+b(1)

式中：σ，τ——破坏面上的正应力和剪应力；

a，b——抛物线待定参数。

图1 二次抛物线型Mohr破坏包络线Fig.1 Curves of the parabolic Mohr failure envelope

图1中(σ1+σ3)/2和(σ1-σ3)/2分别代表岩石在任意应力状态下的MOHR圆心坐标和半径。假设岩石在某一应力状态(σ1，σ3)达到屈服，由此可得：

(2)

对MOHR包络线方程即式(1)求偏微分可得：

(3)

对式(2)求偏微分可得：

(4)

由于包络线和MOHR应力圆外切，故在外切点曲线斜率相等，联立式(3)、(4)得：

(5)

式(1)可变形为：

(6)

将式(5)(6)代入式(2)可得：

(7)

式(7)即本文所建高应力区岩石强度准则表达式，该式具备一定非线性特征，明显不同于传统M-C准则所体现的线性形式。

2 岩石统计损伤本构模型

2.1 基本假定及本构关系

LEMAITREL[18]应变等价性假说的核心假定是宏观应力σ引起损伤材料上的应变ε与有效应力σ′引起无损材料的应变ε′相等，即：

ε=ε′ (8)

考虑三向应力状态，则有：

(9)

式中：D——损伤变量。

假设无损岩石服从广义HOOKE定律，即：

(10)

式中：E——弹性模量；

μ——泊松比。

由于岩石屈服后存在一个残余强度，故引入一个损伤修正系数q[20]，即：

(11)

式中：σr，σ1p——残余强度和峰值强度。

岩石常规三轴试验中，围压σ3=σ2，再结合式(8)～(11)可得损伤本构模型，即：

σ1=Eε1(1-qD)+2μσ3(12)

2.2 统计损伤本构模型

实际上，岩石强度取决于岩石材料中未损部分，本文认为岩石材料一旦受损便失去承载能力，参考曹文贵等[7-11]微元强度度量方法，基于式(7)建立本文高应力区岩石微元强度度量方式：

(13)

联立式(9)、(12)和(13)可得：

(14)

假设岩石微元强度F=f(σ′)[7-14]服从某种概率分布P(F)，则其损伤变量D为：

统计分布主要有Weibull分布[8-10,12,14,17]、幂函数分布[11]、正态分布[7]等，幂函数形式简单，参数较少，其概率密度分布函数P(F)为：

(16)

将式(16)代入式(15)可得：

(17)

结合式(12)、(14)和(17)可得：

(18)

式(18)即本文所建模型。

3 模型参数确定

模型参数的确定是本文模型建立的关键，目前主要有两大类模型参数确定方法。其一是数据拟合法[9]，虽然该方法较为简单，但拟合求解法不能严格满足各项求解条件，参数本身物理意义也不明确。其二是多元函数求极值的方法[11]，利用应力-应变曲线中峰值点的极值性，由此确定的模型参数不仅物理意义较为明确，更可得到适用于不同应力状态的参数表达式，该方法相较前者更具优越性。将式(14)代入(18)再求偏微分得到：

设岩石全应力-应变曲线的峰值点为p，由于其峰值点的极值性，故其峰值条件为：

(20)

将式(20)代入式(19)可得：

式(18)可变形为：

(22)

联立式(21)和式(22)可得：

(23)

式(19)～(25)即本文模型参数m和F0的确定方法。而参数a,b的取值基于图1所示的拟合结果，3个围压条件下的参数a,b值相同。

4 模型验证

4.1 试验方案

本次试验以三叠系中统竹卡组(T2z)英安岩为研究对象，采用YSJ-01-00岩石三轴流变试验机进行常规三轴压缩试验，试样取自某水电站坝肩边坡平洞洞底，取样点实测最大主应力22 MPa，强度应力比2.7，平洞开挖过程中未见岩爆现象，洞壁岩体有较多新生裂隙，存在剥离和掉块现象，参考《水利水电工程地质勘察规范》[1]认为该区域属于高应力区。取样后密封运回实验室制备成Φ50 mm ×100 mm的圆柱样，围压设置为10，20，30 MPa。图2为岩样典型破坏形态，以张剪组合破坏和剪切破坏为主。

图2 岩样破坏形态Fig.2 Failure pattern of rock specimen

4.2 模型验证

现通过高应力区英安岩变形破坏全过程曲线对本文模型进行验证。引用文献[9]中基于线性M-C准则所建立的适用于浅部地层低应力水平岩石的统计损伤本构模型，与本文所建适用于高应力区岩石的统计损伤本构进行对比验证。表1为本文模型参数，图3为两种模型的试验值与预测值对比曲线。由图3可看出，本文所建模型的理论曲线与试验曲线吻合较好，平均相关性系数为0.988 2，误差较小，而文献[9]模型的理论曲线难以描述高应力区英安岩的峰后阶段，高应力区英安岩有较大的应力变化范围，文献[9]模型对于高应力区英安岩的峰值强度和残余强度的识别误差较大，平均相关性系数仅有0.795 3，实际工程应用可能会出现较大误差甚至偏于危险。本文所建模型能够较好地辨识高应力区英安岩应力-应变全过程曲线。为了验证本文模型对高应力区其它岩石的适用性，引用文献[21]、[22]中高应力区砂岩和高应力区花岗岩的三轴压缩试验数据，模型验证结果如图4、图5所示。

表1 统计损伤本构模型参数Table 1 Parameters of the statistical damage model

由图4、图5可看出，本文所建模型对高应力区砂岩和花岗岩的辨识效果较好，证明本文模型反映高应力区岩石变形破坏全过程的适用性。

图3 模型预测对比Fig.3 Model prediction contrast

图4 高应力区砂岩拟合对比曲线Fig.4 Fitting contrast curves for sandstone in high stress areas

图5 高应力区花岗岩拟合对比曲线Fig.5 Fitting contrast curves for granite in high stress areas

4.3 累积损伤的扩展过程

岩石变形破坏过程实际上是累积损伤的扩展过程。图6为依据式(17)绘制的英安岩损伤累积曲线。图6曲线近似“S”型曲线，该曲线初始阶段损伤无明显累积，对应岩石初始加载时应变量较小。接着是迅速攀升阶段，对应随着岩石逐渐加载，内部微裂纹相互摩擦挤压，损伤不断累积汇聚。最后曲线渐缓，表明岩石内部结构完全破坏。对比图3、图6发现，在不同围压下，岩石损伤皆在峰值点附近开始迅速累积，在相同的损伤量下，高围压的轴向应变总是高于低围压，说明低围压的损伤累积更快，围压可能会抑制损伤累积扩展。

图6 损伤累积曲线Fig.6 Cumulative damage curves

5 结论

(1)通过研究英安岩在高应力状态下的变形力学特性，引入适用于高应力区岩体的双参数抛物线型Mohr强度准则，建立一种适用于高应力区岩石的微元强度度量方式。假定本文微元强度度量方式服从幂函数概率密度分布，结合连续损伤理论，考虑损伤变量修正，构建高应力区岩石统计损伤本构模型。

(2)通过引用基于线性Mohr强度准则而建适用于浅部低应力水平岩体的相关统计损伤模型，依据英安岩试验数据与本文所建模型进行对比验证，证明本文引入的双参数抛物线型Mohr强度准则应用于高应力区岩体的优越性，验证了本文模型的可行性和正确性。并采用相关文献中高应力区砂岩及花岗岩试验数据进行辨识，证明本文模型描述高应力区岩石变形破坏全过程的合理性和适用性。

(3)分析英安岩累积损伤的扩展过程，岩石损伤曲线近“S”型，峰值点附近损伤迅速累积，最后无限趋近于1，围压对损伤累积扩展具有抑制作用。