0 引言

岩石蠕变理论是深入解释和分析隧道、高边坡等工程岩体流变现象的重要依据，随着中国基础设施建设的加快，岩石蠕变理论的研究范围愈发宽广^[1-2]。岩石受水、冻融、高温、风化等外界多种因素的共同影响，蠕变行为显著，威胁岩体工程的长期稳定性^[3-4]。

目前关于岩石在不同外界因素共同影响下的蠕变特性研究已取得较多成果，王艳春^[5]以唐口煤矿页岩为研究对象，开展考虑温度-应力-化学耦合效应的三轴压缩蠕变试验，探究不同耦合因素对页岩蠕变特性的影响规律，构建可反映pH值和温度影响的蠕变模型；Qiang等^[6]进行粉砂岩在应力-化学-水流耦合作用下的蠕变试验，分析了多场耦合条件下粉砂岩蠕变作用机理；陈卫忠等^[7]选择黏土岩作为试验对象，设计完成了不同围压下的增、降温排水蠕变试验，基于Perzyna应力理论，构建热-水-力耦合蠕变损伤模型；张明珠等^[8]开展隧洞砂岩在不同温度和围压下的蠕变试验，建立温度-围压耦合作用下的一维本构模型；冯学志等^[9]完成砂岩在水化学-冻融耦合下的蠕变力学试验，探索不同pH值和冻融条件对蠕变速率及长期强度的影响。目前对不同外界因素共同影响的蠕变模型已有较多研究，但是对工程建设背景下冻融和应力耦合的岩石蠕变特性及损伤模型研究相对较少。

鉴于此，本研究以汶川-马尔康高速公路卓克基隧道炭质板岩为研究对象，开展不同冻融循环次数和应力加载下的三轴压缩蠕变试验，分析蠕变应变特征及长期强度。基于试验成果，考虑冻融及应力的耦合效应，构建冻融、受荷总损伤变量，从而建立一个新的考虑冻融、受荷损伤发展的蠕变损伤模型，拓展到三维应力状态，并分析岩石损伤发展规律。采用本研究新建损伤模型识别炭质板岩蠕变数据，验证新建损伤模型的可行性。

1 蠕变试验 1.1 研究背景

研究背景为汶川-马尔康高速公路卓克基隧道，该隧道位于马尔康县境内，长1 917 m，海拔4 564 m，最大埋深约为172 m，隧道穿越三叠系上统新都桥组(T₃x)地层，岩性以炭质板岩为主，围岩等级为Ⅳ-Ⅴ级，围岩蠕变效应显著。

1.2 试验材料及设备

本研究以炭质板岩为研究对象，采用RLW-2000型压缩蠕变仪开展三轴压缩蠕变试验，试验中设置不同冻融循环次数和逐级加载应力。为还原隧道围岩实际应力状态，依据取样深度实测地应力，将试验围压设置为20 MPa。隧道内取样后制备成Φ50 mm×100 mm的圆柱样，岩样基本物理参数如表 1所示。参考文献[10]，将冻融循环次数设为0，5，15，30次，冻结-融化温度设为― 20~20 ℃，冻结和融解时间均为12 h，完成1次冻结和融解作为1个循环。

1.3 试验方案设计

在进行蠕变试验前，开展常规三轴压缩试验，确定炭质板岩破坏偏应力。以冻融循环0次为试验工况，围压设为20 MPa，另设两组平行样，得到常规三轴试验结果，如图 1所示。

破坏偏应力q_f取3组岩样的平均值，则q_f为142.49 MPa。以Δq=Kq_f/c(K为强度削弱系数，c为蠕变荷载等级数)作为蠕变试验每级轴向荷载增量Δq的设定依据^[11]，参考文献[12]，取c为6，K为0.7，则每级加载偏应力分别为16.62，33.24，49.86，66.48，83.1 MPa和99.72 MPa，每一级加载历时设为100 h。

1.4 蠕变试验结果

图 2为炭质板岩经历不同冻融作用后的蠕变试验结果，将各级加载偏应力标于曲线下方。通过Boltzmann线性叠加^[12]的方式处理如图 2所示，绘制分别加载蠕变曲线如图 3所示，以冻融循环5次为例。在图 3中，择取0~90 h时刻的蠕变数据，绘制等时偏应力-轴向应变曲线，如图 4所示。以冻融循环5次为例，由于炭质板岩在偏应力83.10 MPa的加载下于78.33 h时刻结束了加速蠕变，考虑到等时偏应力-轴向应变曲线的完整性，将加载时间78.33 h的应变近似地作为90 h时刻数据点。

由图 2可看出，炭质板岩的蠕变试验曲线呈台阶形，逐级抬升，轴向应变随着时间增长而递增。冻融循环0次下在第6级加载发生蠕变破坏，而冻融循环5，15，30次下的炭质板岩在第5级加载下蠕变破坏，这说明冻融循环作用于岩石可能会降低蠕变破坏的难度。

由图 3可看出，岩石受偏应力作用，瞬时弹性应变瞬间产生，接着蠕变应变随时间增长，增长速率递减，当蠕变应变趋稳时，完成了从衰减到稳定的阶段过渡。稳定蠕变阶段占据了整个蠕变时程的大部分时间比例，当应力超过长期强度，岩石还发生加速蠕变。

由图 4可看出，等时偏应力-轴向应变曲线由0 h段和15~90 h曲线簇组成，其中0 h段为线性段，15~90 h曲线簇具有明显的非线性特征，表现出逐渐向轴向应变横轴倾斜的趋势。采用取拐点^[13]的方法，确定炭质板岩在冻融循环5次下的长期强度为58.84 MPa，通过这样的方法，得到冻融0，15，30次下的炭质板岩长期强度分别为76.02，55.63，51.49 MPa，长期强度随着冻融作用的增强而递减。

2 岩石蠕变损伤模型 2.1 损伤变量的定义

岩石是一种非线性非均质的复杂材料，经历风化、冻融、浸水等诸多因素作用，经受压力会导致岩石性能劣化，损伤即为致使岩石性能劣化的结构微观变化^[14]。炭质板岩经历不同程度的冻融作用，在外界应力作用下，岩石发生蠕变变形，内部损伤持续累积，在岩石蠕变模型中，通过损伤变量考虑材料内部损伤，从某种程度上可体现岩石内部结构损伤劣化过程。结合本研究炭质板岩的实际赋存环境及试验条件，认为岩石损伤受冻融和应力控制，岩石损伤变量应考虑冻融循环及应力的耦合作用，由此还需构建能反映冻融循环与应力耦合作用造成岩石损伤的演化方程。

假设岩石由M个微单元组成，其中无损微单元数量为M₁，受外界荷载和冻融循环作用导致损伤的微单元数量为M₂，M₂中包含冻融损伤微单元M_d和受荷损伤微单元M_s，由此定义冻融、受荷总损伤变量D如式(1)所示：

(1)

从考虑岩石内部受损微单元所占比例出发，将冻融损伤变量D_d定义为M_d与M的比值，受荷损伤变量D_s定义为受荷损伤微单元M_s与微单元总数量减去冻融损伤微单元后的比值，则有：

(2)

联立式(1)~ (2)可得：

(3)

式(3)即为考虑冻融循环及受荷的损伤变量，从式(3)可看出，D并不是由D_d和D_s单一相加而得，这从某种意义上表征了冻融循环和应力荷载致使岩石发生损伤的过程具有耦合关联，两者相互影响，从量值大小来看，式(3)计算得到D的量值小于式(1)。

冻融与荷载作用于岩石时，由于两者作用机理不同，使岩石材料损伤劣化的效能也不同。冻融、荷载作用过程之间发生耦合关联，促使岩石材料冻融损伤及受荷损伤区域之间不断相互反馈，致使岩石材料力学性能降低。

根据唯象损伤力学理论^[15]，将D_d定义为：

(4)

式中，E₀为冻融0次的弹性模量；E_d为岩石冻融后的弹性模量，在本研究中，认为E₀和E_d分别为通过常规三轴压缩试验确定的冻融前后的弹性模量。

考虑到冻融损伤变量D_d作用于岩石材料后，弹性模量的变化是一种衰减型过程，损伤变量逐渐累积，最终趋于稳定，E_d最终会趋于一个稳定值，故采用一个指数型函数修正式(4)，于是有

(5)

式中，α为与损伤相关的指数；t为时间。不妨取E₀为5 GPa，E_d为1 GPa，绘制不同α取值条件下D_d随时间t的变化曲线，如图 5所示。

由图 5可看出，在不同α取值条件下，D_d均随时间t的增长而递增，当t增长到一定程度时，D_d值增长至恒定不变。α值越大，D_d值增至恒定值的速率越快。α值可控制冻融损伤累积速率，使式(5)的结构形式更为灵活，由此看来，采用指数型函数修正式(4)的方法是可行的，式(5)即为冻融损伤变量D_d的演化方程。

岩石材料由于其非均质非线性的特点，当长期受外界荷载影响时，材料内部颗粒位移、微缺陷发育甚至贯通，其内部微结构的变化具有随机性，本研究将这种随机性理解为一种时效概率，且受荷损伤随时间而连续变化，由此致使受荷损伤区域随着时间推移在岩石材料内部随机分布，基于此，本研究假设岩石材料受荷载作用下的时效损伤服从某种概率密度分布函数φ(t)，则受荷损伤变量D_s如式(6)所示：

(6)

Weibull分布在材料寿命检验方面应用广泛，本研究假设岩石在荷载作用下的时效损伤服从Weibull分布，则φ (t)为：

(7)

式中β与m为材料力学性能相关的参数。

将式(7)代入式(6)可得：

(8)

式(8)即为D_s的演化方程。

将式(5)和式(11)同时代入式(3)变形可得：

(9)

式(9)即为D的演化方程。

2.2 蠕变模型

考虑到炭质板岩在受荷瞬间的应变响应是弹性的，蠕变模型中应包含一个弹簧体(H体)。炭质板岩在蠕变过程中，衰减和稳定蠕变行为以黏弹性应变为主，采用由一个H体和牛顿体(N体)并联而成的Kelvin模型，与H体串联。当岩石外界应力超过长期强度后，炭质板岩黏塑性应变累积，还会表现出加速蠕变行为，此时引入一个具有开关功能的黏塑性体，几者组装后得到如下模型，如图 6所示，其模型结构为H-H|N-N|S。

图 6模型结构相应的状态方程为：

(10)

式中，σ为应力；ε为应变；η为黏滞系数；σ_s为长期强度；上标圆点为变量对t的一阶导，下标(1)~ (3)用于区分模型中不同元件的应力及应变，式(10)的第5式仅当应力超过长期强度后生效。

求解式(10)可得：

(11)

式(11)即为图 6模型的一维本构方程，σ-σ_s为开关，仅当σ≥σ_s时生效。

根据损伤力学理论^[15]，若以各向同性描述岩石某一力学参数P的损伤劣化规律，则任意时刻经历时效损伤的参数P (t，α，β，m)如式(12)所示：

(12)

将式(12)代入式(11)可得：

(13)

式(13)即为本研究所建反映冻融、受荷损伤的岩石一维蠕变方程。

在隧道工程中，岩体处于3向应力状态，本研究蠕变试验亦是三轴压缩的形式，因此还需将一维蠕变模型拓展到三维空间。假设岩石各向同性，依据广义胡克定律，三维应力本构关系如式(14)所示：

(14)

式中，σ_m为球应力张量；ε_m为体积应变张量；S_ij为偏应力张量；e_ij为偏应变张量；K为体积模量; G为剪切模量。将岩石内部应力张量分解，则有

(15)

式中，σ_ij为应力张量；δ_ij为单位张量。

本研究不考虑岩石材料的体积蠕变，将体积蠕变忽略不计，仅考虑岩石受荷瞬时的体积变形，引入屈服函数F，将式(13)拓展到三维空间，如式(16)所示：

(16)

式中，(S_ij)₀为初始偏应力张量；(S_ij)_s为岩石长期强度所对应的偏应力张量；F₀为F的初始值；Q为塑性势函数。

F可取如下形式

(17)

式中，J₂为应力偏量不变量，下标2为第2不变量。

本研究炭质板岩蠕变试验中，中间主应力的量值等于围压，则有

(18)

令F₀=1，取屈服函数作为塑性势函数(Q=F)^[16-18]，将式(18)代入式(16)可得:

(19)

式中，σ₁-σ₃为偏应力；(σ₁-σ₃)_s为长期强度。式(19)即为本研究考虑冻融、受荷损伤的岩石三维蠕变损伤模型的轴向蠕变方程。

3 参数求解及模型验证 3.1 参数求解方法

式(19)用于本研究炭质板岩蠕变行为验证，该模型包含G₁，G₂，K，η₂，η₃，α，β，m和(σ₁-σ₃)_s共9个未知参数。其中(σ₁-σ₃)_s为岩石长期强度，通过图 4所示的方法来获取，G₁和K通过广义胡克定律单独求解。另外6个参数则利用数学软件1stOpt，通过擅于处理非线性问题的BFGS算法进行求解。在本研究中，认为E₀和E_d分别为通过常规三轴压缩试验确定的冻融前和冻融后的弹性模量。每一级加载均导致岩石损伤递增累积，而本研究新建模型基于图 3(分别加载蠕变曲线)计算验证得到的参数α，β，m在每一级加载下均有不同值，因此取不同加载等级下α，β，m的平均值作为每种冻融循环条件下的损伤参数，使得每种冻融循环条件下D的累积过程为逐级递增过程。以上即为本研究新建模型的参数求解方法。

3.2 模型验证

利用本研究新建模型，对炭质板岩分别加载蠕变曲线进行模拟，得到模拟结果与试验值对比曲线，如图 7所示，为了比较新建模型模拟效果，引入未考虑损伤的图 6模型，同样拓展为三维形式，与新建模型进行对比。表 2为模型参数表，E₀和E_d不再列出，值得注意的是，当冻融循环次数为0时，E₀=E_d，此时D_d=0，参数α无唯一解，不影响D_s和D的计算。

由图 7和表 2可看出，本研究考虑冻融、受荷耦合效应的蠕变损伤模型对炭质板岩蠕变行为的模拟效果较好，平均R²达0.987 8。而图 6模型对岩石破坏前加载等级的蠕变行为模拟尚可，破坏前加载模拟的平均R²为0.942 5，对最后一级加载下的蠕变曲线模拟较差，尤其是描述岩石加速蠕变行为存在较大误差。

3.3 损伤分析

将表 2中参数代入式(9)，得到D的演化曲线，如图 8所示。

由图 8和表 2可看出，未经历冻融时，炭质板岩冻融、受荷总损伤变量D的累积速度较慢，此时D仅通过外界荷载作用下产生损伤，当岩石经历冻融后，D的累积速度加快。岩石经历冻融循环次数越高，α值越大，D_d值累积越快，对D累积演化的贡献越大，由此D更快地累积到恒定值。总体而言，冻融作用的加强，使得炭质板岩内部结构的损伤发展更快，对岩石蠕变破坏起一定促进作用。

经过模型验证和分析，本研究新建模型对炭质板岩在冻融环境下的蠕变行为模拟具有一定的可行性，但对于描述不同类型岩石在冻融环境下的蠕变特性还有待进一步验证研究。

4 结论

(1) 炭质板岩的蠕变试验曲线呈台阶形，轴向应变随着时间增长而递增。随着冻融作用的增强，长期强度随之递减，冻融循环作用可能会降低岩石蠕变破坏的难度。

(2) 通过概率密度分布的方式定义受荷损伤变量，依据唯象损伤力学理论定义冻融损伤变量，考虑两者之间的耦合作用，构建冻融、受荷总损伤变量。

(3) 针对炭质板岩蠕变行为表现，确定岩石蠕变基础模型结构为H-H N-N S，通过所建冻融、受荷耦合对基础模型进行损伤演化，得到新的蠕变损伤模型，并拓展到三维应力状态。

(4) 结合蠕变模型组建方式，给出模型参数解析方法。采用本研究新建损伤模型模拟炭质板岩蠕变数据，并引入一个传统模型比较验证，对比验证结果，充分显示本研究新建损伤模型的可行性。