基于最大拉应变准则的冻融岩石损伤模型研究

候 超 ，靳晓光 ,2,3，何 杰 ，张 驰

（1.重庆大学土木工程学院，重庆 400044；2.重庆大学山地城镇建设与新技术教育部重点实验室，重庆 400044；3.重庆大学煤矿灾害动力学与控制国家重点实验室，重庆 400044）

据相关数据统计，永久性冻土区和季节性冻土区约占我国国土面积的3/4，当环境温度低于0 ℃时寒区岩体内水分冻结，温度高于0 ℃时则消融.随着季节变换及昼夜更替环境温度在0 ℃上下震荡，岩体内水分发生周期性的冻结和消融，即冻融循环[1-2].随着“一带一路”倡议的实施，大量的工程建设在西部寒区兴起，越来越多的岩土工程遭遇到冻融问题，其中岩石在冻融循环作用下的损伤模型一直是研究者关注的问题，对寒区岩石的损伤劣化机理及岩土工程稳定控制具有重要的理论意义[3-6].

目前，对冻融循环引起的岩石损伤劣化机理已经有了较为一致的认识：岩石矿物颗粒在低温作用下体积收缩，而孔隙和裂隙内的水冻结成冰产生约9.08%的体积膨胀，由于不同矿物颗粒缩胀率有所差异，导致颗粒间的缩胀变形不协调，矿物颗粒与微孔隙之间产生巨大的冻胀力，对岩石造成一定的损伤[7]；温度升高时，岩石内孔隙、裂隙中的冰体融化，随着冻胀力的释放及水分的迁移，进一步促进岩石的损伤；在长期的冻融循环作用下岩石内部微裂纹不断萌生和发展，造成岩石承载力下降[8-9].

迄今为止，国内外学者从宏观、微观及细观等多方面对冻融岩石的损伤劣化机理进行了大量研究：在冻融岩石损伤模型方面，张慧梅等[10]基于岩石内部缺陷分布的随机性，运用损伤力学理论建立了温度-荷载耦合作用下岩石的损伤模型，并通过试验验证模型的合理性；Huang 等[11]假设岩石微元体强度服从Weibull 分布，以弹性模型作为损伤变量结合最大拉应变准则建立了冻融岩石的损伤模型，对模型进行了验证，并应用于寒区隧道的稳定性分析；Jia 等[12]将冻融作用视为疲劳损伤，建立了冻融岩石的疲劳损伤模型，揭示了反复冻融下岩石力学强度衰减机理；为实现对冻融岩石变形全过程的预测，袁超等[13]选用Drucker-Prager 准则，建立了能反映冻融岩石各变形阶段特征的损伤模型；张峰瑞等[14]对经历不同冻融循环次数的花岗岩开展细观特征分析和剪切蠕变试验，基于试验结果提出了冻融岩石损伤黏性元件，构建了花岗岩冻融剪切蠕变本构模型.

总体看来，对岩石的冻融损伤模型研究较少，模型多基于莫尔-库伦准则及Drucker-Prager 准则，不能准确反映岩石的张拉破坏特性和峰后强度特征.众所周知，岩石的抗拉强度远低于抗压强度，岩石微孔隙内水结成冰形成的冻胀力可视为拉应力，岩石在冻融和受荷作用下易发生张拉破坏.因此，本文基于最大拉应变破坏准则和统计损伤理论建立了考虑冻融和荷载耦合作用的岩石损伤模型，推导了模型参数的理论解，利用已有研究结果对本文建立的损伤模型进行了验证，探讨了不同冻融次数下岩石的总损伤曲线演化规律，对模型参数进行了分析，并采用数值模拟方法计算了冻融作用对隧道工程稳定性的影响规律.

1 损伤本构模型建立

1.1 损伤变量及本构方程

基于Lemaitre 应变等价准则[15]，岩石在荷载作用下的有效应力为

式中： σi为名义应力；D为受荷损伤变量.

天然岩石中含有大量微孔隙和裂隙，为了表征岩石微元体的强度引入Weibull 分布函数，则岩石微元体的强度概率密度函数可表示为

受荷损伤变量D为失效单元数Nf与总单元数N的比值，如式（3）.

则可得到受荷损伤变量为

参照张慧梅等[16]提出的冻融损伤模型，将冻融作用下承受荷载的岩石划分为冻融损伤部分、受荷损伤部分和未损伤部分，则

式中： σr为损伤部分承受的残余应力；A为受荷面积；A1为未损伤部分面积；An为冻融损伤部分面积；A2为受荷损伤部分面积.

岩石的冻融损伤变量、受荷损伤变量和总损伤变量如式（8）～（10）.

联立式（8）～（10）可得

联立式（1）和式（6）、（7）、（11）可得

假定岩石只在主应力方向发生损伤，且未损伤部分服从广义胡克定律，结合岩石变形协调关系得到冻融作用下承受三轴压缩荷载岩石在主应力方向的损伤模型为

式中：En为岩石冻融循环n次后的弹性模量； µn为岩石冻融循环n次后的泊松比； σ3为围压； σ1为主应力方向应力；ε1为主应力方向应变.

由宏观损伤力学可知，岩石的冻融损伤变量为

式中：E0为未受冻融作用岩石的弹性模量.

将式（5）、（14）代入式（11）、（13）得

假设岩石微元体强度服从最大拉应变准则，如式（16）.

将式（16）代入式（15）得

1.2 模型参数确定

参数f0、m可用峰值应力法确定[11]，峰值点处应力、应变为

式中： σ1为轴向应力.

峰值点处应力对应变的偏微分等于0，则

将式（18）、（19）代入式（17）得

将式（21）代入式（20）得

由式（22）可得

则

将式（24）代入式（21）得

由式（25）及Weibull 分布函数参数定义可得

由式（24）得

最后，引入修正系数k对式（17）进行修正得

k的取值为

式（28）及式（29）即为本文所建立的基于最大拉应变破坏准则并考虑冻融和荷载耦合作用的寒区冻融岩石损伤模型.

2 模型验证

2.1 实例一：低孔隙率硬质花岗岩

Tan 等[17]对采自西藏嘎隆拉山高速公路隧道中的花岗岩开展了饱和状态下，温度变化范围为-40～+ 40 ℃，循环次数高达150 次的冻融循环试验.文献[17]所采用的花岗岩孔隙率为0.67%，平均单轴抗压强度为135.73 MPa，属于致密的硬质岩石.为验证本文所建立损伤本构模型的适用范围及合理性，利用文献中围压为10 MPa，冻融循环次数为0、50、100、150 次的试验曲线与本文损伤模型计算的理论曲线进行对比分析.同时，基于Huang 等[11]建立的损伤本构模型计算了低孔隙率硬质花岗岩的应力-应变曲线，将其作为对比曲线之一.

表1 给出了不同冻融次数下低孔隙率硬质花岗岩的力学参数.表1 中，N为冻融次数.由于冻融作用对硬质岩石的泊松比影响程度较低，因此，本文在计算过程中，假设不同冻融次数下低孔隙率硬质花岗岩的泊松比数值不变，均为0.15.试验与理论模型的对比结果如图1 所示，图中：理论模型1 代表本文所建立的损伤本构模型；理论模型2 代表文献[11]中的损伤本构模型.

表1 低孔隙率硬质花岗岩力学参数Tab.1 Mechanical parameters of hard granite with low porosity

从图1 中可以看出：当花岗岩冻经历0、50、100 次冻融循环作用后，试验数据曲线与理论模型1、理论模型2 在峰值强度前均具有较好的对应关系；不同的是，在峰后阶段，理论模型1 可较好反映出岩石的残余强度，而理论模型2 的应力在峰值强度后跌落速率较快，无法体现岩石的残余强度特征；当冻融循环次数达到150 次时，理论模型2 在峰值强度前与试验数据吻合度较理论模型1 高，此时，虽然理论模型1 的峰值强度与试验数据曲线之间存在一定的差距，但曲线整体趋势较为一致.可见，本文所建立的冻融岩石损伤本构模型可以适用于低孔隙率硬质岩石，且能较好地体现岩石的残余强度特征.

2.2 实例二：高孔隙率软质红砂岩

张慧梅等[16]对未受冻融循环作用及受冻融循环后的红砂岩开展了三轴力学试验，分析了不同冻融循环次数和围压下红砂岩的力学特性.文献[16]所采用的红砂岩孔隙度为14.26%，单轴抗压强度约为4 MPa，属于高孔隙率软岩[18-19].为验证本文所建立的损伤本构模型在高孔隙率软岩中的适用性，利用文献[16]中的试验数据曲线与本文计算的理论曲线进行对比分析.类似的，此处也加入了与文献[11]中的损伤本构模型的对比.表2 给出了不同冻融次数和围压下高孔隙率软质红砂岩的力学参数.试验与理论模型的对比结果如图2 所示.

图2 高孔隙率软质红砂岩试验与理论曲线对比Fig.2 Comparison of experimental and theoretical curves of soft red sandstone with high porosity

表2 高孔隙率软质红砂岩力学参数Tab.2 Mechanical parameters of soft red sandstone with high porosity

由图2 可知：对于经历不同冻融作用的高孔隙率软质红砂岩，在峰值应力前，试验数据与理论模型1、理论模型2 的吻合性均较好；在峰值强度后阶段，理论模型2 应力快速跌落，与试验数据中岩石具有一定的残余强度形变规律不符；理论模型1 则可较好地反映高孔隙率软质红砂岩的峰后形变特征.以上现象进一步证明了理论模型1 的合理性，也证明了本文所建立的损伤本构可适用于孔隙率较高的软岩.

由图2、3 可知：本文所建立的冻融岩石损伤本构模型具有较广泛的适用性，不仅适用于低孔隙率硬岩，也适用于高孔隙率软岩；同时，本文所建立的模型还可较好地反映出不同类别岩石的峰后形变特征.

3 损伤变量演化规律

由式（5）、（11）、（14）及文献[16-17]中的试验数据计算得到两类岩石在冻融循环下的总损伤演变曲线，如图3 和图4 所示.

图3 低孔隙率硬质花岗岩总损伤演变曲线Fig.3 Total damage evolution curves of hard granite with low-porosity

图4 高孔隙率软质红砂岩总损伤演变曲线Fig.4 Total damage evolution curves of soft red sandstone with high porosity

从图3 中可以看出：不同冻融循环次数下花岗岩的总损伤演变曲线均呈“S”形；未受冻融影响的花岗岩初始阶段总损伤近似为0，受冻融循环作用的花岗岩在初始阶段总损伤为冻融损伤量，且初始阶段的总损伤量随冻融次数的增加而增加；随着外力的加载，岩石内部微裂纹和缺陷不断发展扩张，受荷损伤开始快速增加，花岗岩的总损伤曲线也呈现出加速损伤的趋势，岩石逐渐进入屈服阶段；当岩石强度逐渐丧失后，总损伤曲线斜率开始变缓，并最终趋于稳定值1，此时岩石已完全破坏，总损伤达到最大值.由图4 可知：红砂岩的总损伤曲线也呈现出近似“S”形的趋势，其演变规律与花岗岩类似，不再展开描述.

基于以上分析，将低孔隙率硬质花岗岩和高孔隙率软质红砂岩总损伤演变曲线分为初始损伤阶段、加速阶段和完全损伤阶段（图5）：在初始损伤阶段，岩石内的总损伤量较低，且随着应变的增加，总损伤量近乎不变；在加速损伤阶段，总损伤量随应变的增加而急剧增长，岩石内产生大量的裂纹；在完全损伤阶段，总损伤的增长趋势逐渐变缓并最终达到稳定，岩石在此阶段发生宏观破坏.

图5 冻融受荷岩石总损伤演变规律Fig.5 The evolution of total damage of rock under loading and freeze-thaw action

4 模型参数分析

4.1 参数m 和 f0 对应力-应变曲线的影响

为确定本文所建立的损伤模型中参数m和f0的物理意义及对应力-应变曲线的影响，以冻融循环50 次，围压为10 MPa 的低孔隙率硬质花岗岩以及冻融循环5 次，围压为4 MPa 的高孔隙率软质红砂岩为例进行参数分析，如图6、7 所示.

图6 参数 m 对应力-应变曲线的影响Fig.6 Influence of parameter m on the stress-strain curve

从图6 中可以看出：随着m的增加，低孔隙率硬质花岗岩和高孔隙率软质红砂岩的峰值强度均有所增加；在峰值应力后，应力跌落速率随m的增大而增大，说明岩石的脆性特征逐渐增强.从图7 中可以看出：随着f0的增加，两类岩石的应力-应变曲线峰值点升高且向右移动，达到峰值应力所需的应变量也随之增加；岩石的塑性变形特征随着f0的增加而增加.由此可以推断出，参数m和f0的物理意义分别为岩石的脆性变形特征和塑性变形特征.

图7 参数f0对应力-应变曲线的影响Fig.7 Influenceofparameterf0on the stress-strain curve

4.2 参数 m 和 f0 对总损伤曲线的影响

由式（5）、（11）和式（14）计算得到总损伤量，并分析参数m和f0对总损伤量的影响，如图8、9 所示.由图8 可知：随着参数m的增加，2 类岩石的初始损伤阶段皆有所延长，即，岩石进入加速损伤阶段的临界应变值变大；加速损伤阶段的斜率随m的增加而增加，且岩石进入完全损伤阶段的临界应变值随m的增加而减小，说明了岩石提前进入完全损伤阶段.由图9 可知：随着参数f0的增加，2 类岩石的初始损伤阶段皆有所延长，这与参数m对总损伤曲线的影响规律一致；不同的是，随着f0的增加，2 类岩石加速损伤阶段的斜率均逐渐减小，进入完全损伤阶段的临界应变值增大，即岩石将逐渐滞后进入完全损伤阶段.以上现象进一步表明参数m可反映岩石的脆性形变特征，而参数f0则代表岩石的塑性形变特征.

图8 参数 m 对总损伤量的影响Fig.8 Influence of parameter m on the total damage

图9 参数 f0 对总损伤量的影响Fig.9 Influence of parameter f0 on the total damage

5 工程算例分析

为分析冻融循环对大型工程稳定性的影响，本节将以隧道工程为算例进行数值模拟分析.图10 给出了工程算例的计算模型.如图所示，该隧道断面为三心拱形，拱顶埋深24.5 m，隧道围岩为软质红砂岩.该隧道采用全断面开挖法，开挖循环进尺为2 m，计算开挖总长度为20 m.隧道支护方式为锚喷支护，喷射混凝土等级为C20，锚喷衬砌厚度0.15 m，锚杆截面面积为490.63 mm2，长度为4 m，锚杆间弧长1.8 m.

图10 隧道工程算例计算模型（单位：m）Fig.10 Calculation model of tunnel engineering example （unit: m）

基于文献[16]中所给出不同冻融次数下的红砂岩力学参数，在大型岩土工程有限元软件MIDAS GTS 中分别计算了冻融循环次数为0、5、10 次时，隧道开挖后的位移场和应力场分布.不同工况下隧道围岩的力学参数见表3，在计算过程中假设围岩的容重和内摩擦不随冻融次数而改变.计算得到的围岩竖向位移场、围岩最大主应力矢量场以及衬砌结构最大主应力场分别如图11～13 所示.

图11 不同冻融次数下隧道围岩竖向位移Fig.11 Vertical displacement of surrounding rock for tunnel under different freeze-thaw cycles

表3 不同冻融次数下隧道围岩力学参数Tab.3 Mechanical parameters of surrounding rock for tunnel under different freeze-thaw cycles

从图11 中可以看出：隧道开挖后拱顶和底部的围岩垂直位移量较大；当围岩经历0、5、10 次冻融循环后，隧道拱顶的最大下沉量分别为3.75 、3.99 、4.42 mm，经历10 次冻融循环后，隧道拱顶的围岩下沉量增加了17.87%；类似地，当围岩经历0、5、10 次冻融循环后，隧道底部的最大隆起量分别为4.73 、5.05 、5.64 mm，经历10 次冻融循环后，隧道底部的围岩隆起量增加了19.24%.

从图12 中可以看出：隧道开挖后拱顶与底部分布拉应力积聚区，隧道两侧围岩分布压应力积聚区；当围岩经历0、5、10 次冻融循环后，隧道围岩的最大拉应力极值分别为0.297 、0.300、0.305 MPa，经历10 次冻融循环后，隧道围岩最大拉应力极值增加了2.7%；当围岩经历0、5、10 次冻融循环后，隧道围岩的最大压应力极值分别为0.448、0.450、0.457 MPa，经历10 次冻融循环后，隧道围岩的最大压应力极值增加了2.01%.隧道围岩的最大主应力随冻融循环变化较小，与文献[11]中得出的结论一致，这是由于围岩内主要是自重应力，自重应力在冻融循环下基本保持不变.因此，围岩的最大主应力变化量很小，而由于岩石力学性质的改变，围岩位移随冻融循环的变化量则较大.

此外，隧道衬砌结构的应力分布也受到冻融作用的影响，从图13 中可以看出：隧道衬砌拱顶和底部主要为拉应力积聚区，衬砌两侧分布压应力积聚区，随着冻融循环次数的增加，隧道衬砌拉、压应力极值也逐渐增大；当围岩经历0、5、10 次冻融循环后，隧道衬砌最大拉应力极值分别为4.88、5.25、5.93 MPa，经历10 次冻融循环后，隧衬砌最大拉应力极值增加了21.52%；当围岩经历0、5、10 次冻融循环后，隧道衬砌最大压应力极值分别为3.75、3.99、4.42 MPa，经历10 次冻融循环后，隧衬砌最大压应力极值增加了17.87%.冻融作用后隧道围岩的力学性质劣化，导致围岩压力传递至衬砌结构，引起衬砌应力的增加.由此可知，虽然冻融作用对隧道围岩应力场影响不大，但随着围岩的软化也可引起隧道结构的破坏，影响隧道的整体稳定性.因此，在寒区进行岩土工程建设时，应特别注意对冻融循环对岩土性质的影响并采取适当的防治措施.

图13 不同冻融次数下隧道衬砌最大主应力Fig.13 Maximum principal stress in tunnel lining under different freeze-thaw cycles

6 结 论

本文基于统计损伤力学理论，假设岩石微元体强度服从Weibull 分布，结合最大拉应变破坏准则建立了考虑冻融和荷载耦合作用的寒区冻融岩石损伤本构模型；基于试验数据及已有的损伤本构模型验证了所建立的损伤本构模型的正确性，并对总损伤演变过程及模型参数进行了分析；最后，在数值模拟软件中分析了冻融循环对隧道工程稳定性的影响，得到的主要结论如下：

1） 以低孔隙率硬质花岗岩和高孔隙率软质红砂岩为验证对象，结合试验数据及已有的损伤本构模型确定所建立模型的合理性和适用范围，结果表明，本文所建立的模型可较好地重现2 类岩石的应力-应变全过程曲线，且能较好地反映岩石的峰后强度.

2） 不同冻融循环次数下低孔隙率硬质花岗岩和高孔隙率软质红砂岩的总损伤演变曲线均呈“S”形，且总损伤演变曲线可分为初始损伤阶段、加速阶段和完全损伤阶段.

3） Weibull 分布参数m和f0分别代表岩石的脆性特征和塑性特征；随着m增加，2 类岩石的应力峰值逐渐增加，峰后应力跌落速率逐渐增大，总损伤曲线加速损伤阶段的斜率也逐渐增加，岩石提前进入完全损伤阶段；随着f0的增加，岩石达到峰值应力所对应的应变逐渐增加，总损伤曲线加速损伤阶段斜率逐渐减小，岩石滞后进入完全损伤阶段.

4） 采用数值模拟方法，计算了冻融循环对隧道围岩位移、应力及隧道衬砌应力的影响，随着冻融次数的增加，隧道围岩的垂直位移、最大主应力及衬砌最大主应力均逐渐增加；冻融循环对隧道围岩的应力分布影响较小，但由于围岩物理力学性质的劣化，围岩压力传递至衬砌可导致隧道结构产生损伤.