工程地质体在冻融、荷载等多种环境因素下的力学特性研究已成为冻融灾变防治与冻融岩体稳定性急需解决的关键性科学问题，这对于指导工程实践具有重要的理论意义^[1-2].天然岩体是由许多矿物经长期的地质作用而形成的集合体，内部存在天然形成的结构面以及经过地应力作用形成的后期结构裂隙等，这些宏微观缺陷在冻融作用下势必会导致岩体的损伤破坏等劣化现象^[3].

目前，国内外学者在岩石冻融损伤方面的研究已经取得一些成果，如Papk等^[4]通过对玄武岩、闪长岩进行饱水冻融循环试验，并结合CT与SEM等技术观察岩石在冻融过程中内部微缺陷的结构变化，来研究冻融条件下岩石的物理力学性能；Chen等^[5-6]研究了水的饱和度对岩石在冻融作用下劣化性能的影响；Yamabe等^[7-8]研究了温度在岩石冻融循环试验中的影响；Khanlari等^[9-10]通过冻融循环和将液氮作为冻融条件研究砂岩和煤岩的物理性质；Matsuoka^[11]通过试验方法研究了冻结作用下岩石破裂力学机制；Nicholson等^[12]对含有宏观裂隙的沉积岩在冻融风化作用下岩石物理性能的劣化进行了研究；Prick^[13]通过冻融循环作用对多孔石灰岩样品的体积膨胀行为进行了研究；Ruiz等^[14]通过冻融循环试验利用X射线扫描对白云岩孔隙结构的演变特征进行了研究；Bayram^[15]通过冻融循环与单轴抗压试验对土耳其寒区的天然石灰岩进行了力学强度损失研究，建立了预测强度劣化统计模型；Huang等^[16]考虑热应力和水在椭圆腔(裂缝)中的流动影响，提出了一种预测岩石椭圆腔中冻胀压力的数学模型.徐光苗等^[17-18]通过对岩石进行冻融试验，分析了岩石损伤劣化破坏的机理及其引起岩石内部微缺陷发生变化的影响因素；文献[19-21]采用CT扫描试验，建立以CT数为变量的损伤函数，分析岩石的细观损伤扩展规律，以及建立CT数、岩石强度与冻融循环次数的联系；张慧梅等^[22-23]通过冻融单轴抗压试验，建立了单向应力状态下冻融岩石损伤模型；刘红岩等^[24-25]以相似材料预制的节理岩体进行了冻融试验，结果发现，冻融对节理岩体力学性能及其损伤破坏程度有很大影响；陈有亮等^[26-27]通过在不同环境条件下进行冻融单轴压缩试验，分析了岩石的初始损伤对冻融破坏的影响；李新平等^[28]对寒区裂隙岩体在冻融单轴荷载作用下建立了损伤劣化模型并对其损伤劣化机理进行了讨论.

从已有研究成果来看，以往主要是针对完整岩石的冻融试验研究，而对于具有不同几何特征的裂隙岩体冻融试验的研究甚少，已有研究主要停留在试验层面上，且主要是以单向受力状态为主^[22-28].对于在不同冻融循环次数下的三向受力状态的岩体冻融损伤模型少有报道，且建立冻融荷载下的损伤模型主要是完整岩石微细观损伤模型，未能将工程实际中岩体内部宏观尺度的缺陷所引起的损伤考虑进来.后来刘红岩等^[24-25]只是在试验的层面对预制的节理岩体做了冻融试验研究，发现这种宏、微观尺度的缺陷对岩体的力学性能影响较大.鉴于此，为更好地反映工程实际中岩体的受力情况，在前人研究的基础上，根据Lemaitre应变等效性假设，将宏观节理引起的损伤考虑进来，建立宏细观缺陷的冻融受荷节理岩体三向应力状态下的复合损伤模型，并通过现有的试验资料对模型进行合理性验证.

1 节理裂隙岩体冻融受荷损伤演化方程 1.1 岩石受荷作用下细观损伤

从统计损伤的角度考虑，因岩石内部各种缺陷的随机性，考虑到岩石材料在加载过程中的损伤是一个连续的过程，且在缺陷尺度很小的情况下，假设^[29-30]：1)岩石材料在宏细观上为各向同性；2)岩石微元破坏前服从虎克定律, 即微元具有线弹性性质.因此，采用Weibull概率分布，密度函数为

$ P\left( \varepsilon \right) = \frac{m}{{{F_0}}}{\left( {\frac{\varepsilon }{{{F_0}}}} \right)^{m - 1}}\exp \left[ { - {{\left( {\frac{\varepsilon }{{{F_0}}}} \right)}^m}} \right]. $

(1)

式中：P(ε)为岩石在微应变下对应的强度分布概率，ε为岩石在力作用下的微应变，m和F₀为Weibull分布参数.其中m反映的是岩石的脆性参数，F₀代表了岩石的宏观平均强度^[29].

随荷载不断增加，当加载到一定的应变水平时，岩石内部微元数目的破坏数量为

$ {N_{\rm{s}}} = \int_0^\varepsilon {NP\left( x \right){\rm{d}}x} = N\left[ {1 - {{\rm{e}}^{ - \left( {\frac{\varepsilon }{{{F_0}}}} \right)m}}} \right]. $

(2)

定义岩石损伤变量D为岩石在荷载作用下岩石微元的破坏数n与无损伤时总微元数目N的比值，其取值范围为0~1，即

$ D = \frac{{{N_{\rm{s}}}}}{N}. $

(3)

将式(2)代入式(3)即可得到岩石在受力状态下的损伤演化方程：

$ D = 1 - {{\rm{e}}^{ - \left( {\frac{\varepsilon }{{{F_0}}}} \right)m}}. $

(4)

1.2 物理参数的求解

依据岩石的的本构关系曲线，采用具有普遍明确物理意义的几何关键点求导法.首先，在峰荷前(包括峰荷点)建立本构曲线关系图上的几何约束条件：

$ {\varepsilon _1} = 0,{\sigma _1} = 0; $

$ {\varepsilon _1} = 0,\frac{{{\rm{d}}{\sigma _1}}}{{{\rm{d}}{\varepsilon _1}}} = E; $

$ {\varepsilon _1} = {\varepsilon _{\rm{f}}},{\sigma _1} = {\sigma _{\rm{f}}}; $

$ {\varepsilon _1} = {\varepsilon _{\rm{f}}},\frac{{{\rm{d}}{\sigma _1}}}{{{\rm{d}}{\varepsilon _1}}} = 0. $

ε_f，σ_f分别为完整岩石的峰值应变、峰值应力.具体推导过程见文献[31].两参数的表达式如下：

$ \left\{ \begin{array}{l} m = \frac{{{\varepsilon _{\rm{f}}} - \left( {1 - 2\mu } \right){\sigma _3}/E}}{{{\varepsilon _{\rm{f}}}\ln \frac{{E{\varepsilon _{\rm{f}}}}}{{{\sigma _{\rm{f}}} - 2\mu {\sigma _3}}}}},\\ {F_0} = \left[ {{\varepsilon _{\rm{f}}} - \left( {1 - 2\mu } \right){\sigma _3}/E} \right]{\left[ {\frac{{m{\varepsilon _{\rm{f}}}}}{{{\varepsilon _{\rm{f}}} - \frac{{\left( {1 - 2\mu } \right){\sigma _3}}}{E}}}} \right]^{1/m}}. \end{array} \right. $

(5)

1.3 岩石冻融作用下细观损伤

岩石在冻融循环作用下，其内部结构会受到损伤劣化，最终导致岩石力学性能弱化.岩石在冻融循环过程中内部微细观缺陷发生着非常复杂的损伤演化过程，损伤劣化的程度可以通过岩石外在的宏观物理量来反映，因此，可以采用宏观唯象损伤力学方法^[28].以岩石的弹性模量作为内部损伤度量的基准来描述岩石的冻融损伤变量，即

$ {D_n} = 1 - \frac{{{E_n}}}{{{E_0}}}. $

(6)

式中：E₀和E_n分别为冻融前和n次冻融循环后的岩石弹性模量.

1.4 岩石冻融荷载下细观损伤

根据Lemaitre的应变等价性假设，文献[28, 32]将岩石的初始损伤状态作为岩石的基准损伤状态，然后提出了推广后的应变等价原理.将岩石冻融后的损伤状态视为第一种损伤状态，冻融受荷造成的损伤视为第二种损伤状态，如下为两种损伤状态的本构关系：

$ {\sigma _n} = {E_0}\left( {1 - {D_n}} \right){\varepsilon _n}, $

(7)

$ \sigma = {E_n}\left( {1 - D} \right)\varepsilon . $

(8)

式中E_n为冻融后的弹性模量.

将式(6)代入式(8)可得到如下本构关系：

$ \sigma = {E_0}\left( {1 - {D_{\rm{t}}}} \right)\varepsilon , $

(9)

其中

$ {D_{\rm{t}}} = {D_n} + D - {D_n}D. $

(10)

式中：D_t为岩石冻融荷载耦合后的损伤，D为岩石受荷后引起的细观损伤，D_nD为两种损伤的耦合项.式(10)表明，岩石在冻融与荷载耦合作用下加剧了损伤的进行，表现出非线性特征的性质.分析认为，岩石在冻融过程中使得岩石内部的微观缺陷增强，相对密度提高，因此，岩石在受力的作用下，内部晶粒间的滑移与错动受到阻碍.相反，在岩石受力过程中，内部微缺陷闭合压密，导致冻融时液体的存在方式发生了变化，对液体产生的冻胀力和渗透力起到了一定的限制作用.由两种方式引起的损伤相互作用，相互影响，最终导致岩石的力学性质发生了变化.因此，两种作用形式下引起的损伤不能简单地理解为叠加，耦合项前的“-”号代表由于耦合作用使得总损伤有所弱化^{[28, 32]}.

由式(4)、(6)及(10)可得由于岩石的冻融受荷引起的损伤演化关系为

$ {D_{\rm{t}}} = 1 - \frac{{{E_n}}}{{{E_0}}}{{\rm{e}}^{ - \left( {\frac{\varepsilon }{{{F_0}}}} \right)m}}. $

(11)

由式(11)考虑两种极端情况：当岩石仅受到荷载损伤时，E_n=E₀，此时D_t=D；当岩石仅受到冻融作用时，ε=0，此时D_t=D_n.

1.5 冻融荷载下节理岩体宏细观缺陷耦合损伤

基于文献[33-34]中两种损伤变量耦合的基本条件假设：1)采用人的视力能不能看见宏细观缺陷作为区分的标准，假设宏细观两种损伤具有均匀一致的各向同性.2)采用不同的描述方法来分别表述宏细观缺陷.3)依据损伤力学中Lemaitre提出的基本假设，利用应变等效原理进行耦合，而不是两种损伤变量简单的叠加.如图 1所示，假设图 1(a)~(d)分别代表同时含有宏细观缺陷岩体、仅含宏观缺陷岩体、仅含细观缺陷岩体以及完全不含缺陷岩体，对应的弹模分别为E₁₂、E₁、E₂、E₀，在外载作用下对应的应变量分别为ε₁₂、ε₁、ε₂、ε₀，基于Lemaitre应变等效假设，则有

$ {\varepsilon _{12}} = {\varepsilon _1} + {\varepsilon _2} - {\varepsilon _0}. $

(12)

假设在外载作用方向上造成的宏、细观损伤变量分别为D₁和D₂，耦合后的损伤变量为D₁₂，则有

$ \frac{\sigma }{{{E_{12}}}} = \frac{\sigma }{{{E_1}}} + \frac{\sigma }{{{E_2}}} - \frac{\sigma }{{{E_0}}}. $

(13)

根据Lemaitre假设可得

$ \left\{ \begin{array}{l} {E_{12}} = {E_0}\left( {1 - {D_{12}}} \right),\\ {E_1} = {E_0}\left( {1 - {D_1}} \right),\\ {E_2} = {E_0}\left( {1 - {D_2}} \right), \end{array} \right. $

(14)

将式(14)代入式(13)可得

$ {D_{12}} = 1 - \frac{{\left( {1 - {D_1}} \right)\left( {1 - {D_2}} \right)}}{{1 - {D_1}{D_2}}}. $

(15)

对于式(15)中两种极端情况为：1)当D₂=0，将其代入式(15)可得：D₁₂=D₁，此时岩体的耦合损伤变量就等于岩体的宏观损伤变量；2)当D₁=0，将其代入式(15)得：D₁₂=D₂，即此时岩体的耦合损伤变量等于岩体的细观损伤变量.

根据上面推导出的宏细观损伤耦合公式(15)，将完整岩石初始损伤状态视为岩石的基准状态，将预制节理引起的损伤视为宏观损伤D₀，岩石在冻融与荷载耦合作用引起的损伤视为细观损伤D_t，将预制节理岩体在冻融与荷载共同作用下引起的总损伤视为D_m，则式(15)变为

$ {D_{\rm{m}}} = 1 - \frac{{\left( {1 - {D_0}} \right)\left( {1 - {D_{\rm{t}}}} \right)}}{{1 - {D_0}{D_{\rm{t}}}}}, $

(16)

其中

$ {D_0} = 1 - \frac{{{E_{\rm{j}}}}}{{{E_0}}}. $

(17)

式中：E_j为节理岩体的弹性模量，E₀为完整岩石在基准损伤状态下的弹性模量.D_t的含义及计算同前.

2 模型实例验证与分析

为进一步验证模型的有效性及合理性，引用文献[35]中的试验数据，文献中试验采用三峡库区石英白砂岩的三轴压缩，围压恒定设为10 MPa，预制节理倾角设定为30°，泊松比μ=0.11，将各贯通率的试验组岩样分别进行0，40，80，120，160次的循环冻融试验，结果表明，冻融循环次数越大，含水率越高，贯通率越大，试件损伤破坏越严重.根据文献[35]中试验数据，如表 1所示，依据本文所建立的损伤模型，分别讨论岩样贯通率(K=0(完整)，0.25，0.5，0.75)与冻融循环荷载次数(0，40，80，120，160次)在模型应用中的效果.岩样受力后经历裂隙压密、弹性、屈服、破坏、残余应力等阶段，在加载过程中，裂隙压密及线弹性阶段微缺陷损伤不会发生，进入非线性弹性阶段微缺陷开始萌生，数量开始增多，到达屈服时，内部微裂隙发展迅速加快，直到承载力峰值点.此后随着变形的不断增大，承载力下降，表现出应变软化的特征，此时微裂纹逐渐贯通，形成宏观节理.在工程实际中最关心的是岩石的最大承载力即应力峰值点，因此，可以将峰值点处的损伤值作为衡量岩石破坏损伤程度的一个基准值.利用式(4)计算完整岩石受荷作用下的细观损伤D，这里应变ε取峰值应变，峰值应变可通过文献[35]中试验得出的已知峰值强度与弹性模量计算得出，当ε=ε_f峰值应变时，天然状态下的损伤D=0.066，饱和状态下的损伤D=0.079.岩石冻融损伤(D_n)由式(6)计算得出，由式(10)计算得出冻融与荷载耦合作用下的损伤D_t，由式(16)、(17)计算考虑宏细观损伤耦合作用引起的损伤D_m.

图 2为天然与饱和状态下总损伤随冻融循环次数的变化, 可以看出，各贯通率的损伤曲线随冻融次数的增加呈递增的变化规律.饱和状态下各贯通率的损伤严重程度明显大于天然状态，在贯通率K=0的情况下，冻融循环次数在80次内，两种状态下的损伤变化差异并不明显.80次之后，两者出现明显的差别，损伤值大概最大相差0.4.随贯通率不断增大，两种状态下的损伤随冻融次数的变化曲线出现明显差异，且冻融次数越多，差异性越明显.在冻融120次之前(包括120次)损伤曲线增长速率较快，超过120次后，曲线变化速率放缓，且两种状态下损伤最为严重的曲线随着冻融次数的增加均未趋向于1，两种状态下K=0.75对应的曲线损伤最为严重，极限值分别趋向于0.55与0.91.分析表明，由于水的冻融作用，饱和试验组的损伤随冻融次数增加变化效果较为敏感，且含水率越高，损伤越严重，刚度弱化越大，抗冻性能就弱.而损伤未趋向于1，表明岩体的承载能力未完全丧失，因围压的存在，对提高岩体的强度有一定的优化作用.

图 3为天然与饱和状态下各冻融次数的损伤随贯通率大小变化规律.可以看出，各冻融循环次数的损伤随贯通率的不断增大呈递增的趋势.天然状态下，在较低的冻融循环次数下，如40次，总损伤随贯通率的增加变化率较快，曲线变化呈下凹型增长，当冻融循环达80次之后，曲线随贯通率的增大变化率变缓，曲线变化呈上凹型增长.饱和状态下，各冻融次数下的总损伤随贯通率的增加一直处于迅速增长的状态，较天然状态更为敏感.对比两种状态下的损伤曲线发现，饱和状态下，各冻融次数间的总损伤变化曲线幅差值较天然状态大，尤其在80次冻融循环之后，两种状态下总损伤随贯通率的变化差异更大，但极限值都未趋向于1.在120次冻融循环后，两种状态下的试件在K=0.5和0.75情况下，天然状态下损伤变量超过0.5，饱和状态下损伤变量超过0.8，饱和状态的损伤程度更严重，损伤也未达到1.说明在饱和状态下岩体经过数次的冻融循环完全可以丧失承载能力，因围压的存在导致存在部分残余应力.天然试样组的损伤变量随贯通率的增加变化不太明显，而饱和试验组的损伤变量随贯通率的增加变化较大，损伤严重，说明贯通率的增加对饱和试验组中损伤变量的影响较天然试验组强烈.

图 4为两种状态下各贯通率试验峰值强度与理论峰值强度随冻融循环次数变化.由式(9)及式(16)可知，考虑宏细观耦合作用引起总损伤岩样的峰值强度计算公式为：σ_p=E₀(1－D_m)ε_f+2μσ₃，其中，ε_f为峰值应变，σ_p为峰值应力，μ=0.11为泊松比，天然状态E₀=26.8 GPa，饱和状态E₀=22.67 GPa，其他参数意义同前.由图 4可知，随冻融循环次数增加，理论曲线与试验曲线都近似呈线性化逐渐降低，理论值与试验值总体上吻合良好.图 4各贯通率下的岩样，经过160次冻融循环后，天然状态下的峰值强度分别为σ_p=116.62, 102.87, 83.92 MPa，与未冻融前相比，强度分别降低18.5%, 20%, 22.9%；饱和状态下的峰值强度σ_p=84.93, 65.3, 63.65 MPa，与未冻融前相比，强度分别降低29.7%, 39.1%, 36.7%.其中饱和状态下贯通率K=0.75时，理论值与试验值较前两种贯通率下的误差大些.造成这种误差大的内因可能是贯通率越来越大，岩体的冻融损伤占据主导作用，大大弱化了节理引起的损伤作用，外因可能是含水率较高受外界温度等影响较敏感.从以上数据分析可知，岩体在长期循环冻融下，随冻融次数增加，饱和试验组的强度衰减较大，抗冻融破坏能力下降较快.

对比图 4(a)、(b)、(c)发现，各贯通率下的曲线变化趋势基本一致.从以上数据分析可知，两种理论峰值强度与试验值峰值强度经过160次冻融循环后分别在降低，与实际相符，但各贯通率下峰值强度降低的程度不同，贯通率越大，峰值强度衰减越严重.分析认为，贯通率越大，水的贯通性比较好，集水能力较强，在冻融往复循环中，水的冻胀力与渗透力前后应力差大，造成强度衰减程度就大.饱和状态下的理论峰值强度与试验值都低于天然状态下的峰值强度，在前120次冻融循环中天然试验组的峰值强度变化曲线斜率较小，峰值强度衰减较慢，而饱和试验组峰值强度曲线斜率较陡，峰值强度衰减较快.分析表明，贯通率、冻融循环次数以及含水率的高低对岩体强度劣化有一定影响.

图 5为两种状态下各冻融次数的峰值强度理论与试验值随贯通率变化规律.可以看出，各冻融循环次数下，理论与试验峰值强度都随着贯通率的增大呈近似线性递减的趋势，各冻融次数下，天然与饱和状态下的峰值强度曲线近似平行，斜率相当，反映出两种状态下峰值强度衰减的速率基本相同.对比图 5中各图发现，随着冻融循环次数的增加，天然与饱和状态峰值强度曲线间的幅值差在增大，冻融循环达到一定次数(120次)时这种间距差变化不再明显，说明贯通率的大小对峰值强度衰减的速率影响较小，而冻融循环次数对峰值强度的衰减影响较大.

图 6为两种状态下各贯通率的弹性模量理论与试验值随冻融次数变化规律.弹性模量的理论计算公式可表示为E=E₀(1－D_m)，天然状态弹性模量E₀=26.8 GPa，饱和状态弹性模量E₀=22.67 GPa.从图 6可以看出，理论值与试验值随冻融次数的变化趋势基本一致，且理论曲线与试验曲线吻合良好.天然与饱和状态的弹性模量随着冻融次数的增加呈现递减的趋势，说明在发生冻融循环时，由于冻胀力与渗透力的往复循环作用，导致内部微缺陷结构发生变化，加之节理裂隙等宏观缺陷的耦合作用，最终导致弹性模量弱化，岩体抗变形性能得到劣化，饱和试验组的衰减程度更大.随着贯通率的不断增大，天然与饱和状态曲线间的幅值差在随着冻融循环次数的增大而增大，说明不但贯通率与冻融次数对弹性模量有劣化影响，且含水率越高越加速这种劣化效果.

图 7为两种状态下各冻融次数的弹性模量理论与试验值随贯通率变化规律.可以看出，理论弹性模量与试验值随贯通率的变化基本一致，理论曲线与试验曲线吻合良好.各冻融次数下弹性模量随贯通率增加呈递减的变化趋势，说明节理贯通率大小对岩体物理力学性能有很大影响，且饱和组试样较天然组变化衰减较为敏感.随冻融循环次数的增加，两种状态下的弹性模量曲线间的幅值差在不断增大，天然组试件弹性模量曲线的斜率变化不明显，饱和组试件的弹性模量曲线斜率变化比较显著，说明饱和试验组弹性模量弱化比较严重，抗冻能力差.当达到一定冻融次数后，弹性模量衰减变化的程度也在变缓.

3 结论

1) 天然与饱和状态下各贯通率的损伤演化曲线随冻融次数的增加而增大，各冻融循环次数的损伤演化曲线随贯通率的增大而增加，且饱和试验组受冻融效应较天然组增加敏感.冻融次数较多，贯通率较大，损伤越严重，岩体刚度衰减越大，抗冻性能越差.

2) 天然与饱和状态下各贯通率的试验峰值强度与理论峰值强度随冻融循环次数增加而降低，各冻融次数的试验峰值强度与理论峰值强度随贯通率增加而降低.随冻融次数增加，饱和试验组的峰值强度曲线斜率较大，衰减较快；随贯通率增加，两种状态下的理论与试验峰值强度曲线斜率近似平行，衰减的速率基本相同，说明峰值强度衰减受冻融次数的影响较大.

3) 天然与饱和状态下各贯通率的弹性模量理论值与试验值随冻融次数增加而递减，各冻融次数的弹性模量理论值与试验值随贯通率增加而递减.无论是贯通率还是冻融循环次数的增大，天然与饱和状态曲线间的幅值差都在增大，且饱和试验组反映更加敏感，说明贯通率与冻融次数对弹性模量具有一定的劣化影响，

4) 所建模型计算结果的理论曲线与试验曲线吻合良好，验证了所建模型的可行性与有效性.