冻融作用下多年冻土边坡局部稳定性研究

2020-01-08褚志成雷胜友原喜忠郭长志王汉文陈鹏辉陈虹宇

水利与建筑工程学报 2019年6期

褚志成，雷胜友，原喜忠，郭长志，王汉文，陈鹏辉，陈虹宇

(长安大学公路学院, 陕西西安 710064)

随着我国经济的蓬勃发展，为满足人民日益增长对美好生活的需求，我国在寒区基建的投入也逐年增加。在全球气候变暖的大环境下，我国多年冻土分布整体呈现出南界北移的趋势，寒区边坡的失稳滑塌问题也逐渐增多，给人们的生产生活带来了极大的隐患。究其缘由，周期性气温冷暖交替引起的反复冻融是造成多年冻土边坡失稳破坏的重要原因。在冻融过程中，冷季随着气温的下降，季节活动层土体水分在冻结锋面积聚，而到了暖季，春融期积雪融化以及原冻结土层融化使得土体水分处于过饱和状态[1-3]。就多年冻土边坡而言，暖季气温的升高直接影响到坡体温度场的分布，边坡处于融化期，活动层水分的增加降低了潜在滑动面的抗剪强度，造成边坡土体承载力急剧下降；当颗粒空隙完全充满水时，应力的传递方式发生了变化，由起初的颗粒传递变为孔隙水传递，导致土体的粘结作用降低，从而降低了土体的强度。

目前，冻土边坡冻融滑塌稳定性分析已有大量的研究成果，武鹤等[4]考虑了滑坡面的端摩阻力作用，并利用基线平衡原理推导出了寒区土质边坡安全系数的计算公式。刘义高[5]利用有限差分法从土体强度的角度分析了西藏地区边坡的稳定性，并揭示了该地区边坡的破坏机理。葛琪等[6]考虑了水分迁移对季冻区边坡稳定性的影响，结合有限差分法，分析了季冻区正融土坡的失稳状态。高樯等[7]利用有效应力极限平衡法在不考虑冻土融化过程中的孔隙水压力增大的影响的情况下，推导出了斜坡稳定性系数的计算公式。曾韬睿等[8]考虑了融化土体的渗透力，基于修正的传递系数法推导了冻融边坡稳定性的计算公式，并表明融化深度是冻土边坡稳定性的主要影响因素。孙国栋等[9]考虑冻融循环以及渗流对寒区高陡边坡的作用，提出了基于冻融折线型滑移面的概念，并结合框架锚杆支护边坡得到了寒区高陡边坡稳定性评价方法。可看出，对于寒区土质边坡冻融稳定性分析中，大多采用静力平衡法、传递系数法、强度折减等整体分析方法，而且大多未考虑冻融循环对土体的损伤作用以及冻土边坡的局部冻融的特殊性，所以上述方法在多年冻土边坡稳定性分析中不能完全适用。基于此，本文以东北大兴安岭地区某一级公路路堑土质边坡为研究对象，考虑冻融循环对边坡土体的损伤作用以及考虑多年冻土边坡的局部冻融，建立有限元边坡稳定性分析方法，为多年冻土地区边坡的稳定性评价以及边坡防护治理提供参考。

1 计算方法

1.1 带有相变的非稳态平衡方程

根据热传导和质量迁移原理，仅考虑对流换热和土体在冻融过程中的相变问题，则二维非稳态热传导方程表示为:

(1)

式中：ρ为土体密度；T表示土体温度；q为内热源强度；在相变区为[TbTp]时，土体的导热系数λ及容积热容量C的取值如式(2)、式(3)所示：

(2)

(3)

式中：Cf、Cu分别表示冻土和融土的比热；λf、λu分别表示冻土和融土的导热系数；L为水在结晶或融化时释放的潜热，故相变热Q=Lρd(w-wu)，其中wu为未冻水含量，ρd为土体的干密度。在土体融化后，未冻水含量认为是定值，如式(4)所示[10]：

(4)

1.2 土体破坏准则的选取

土体的本构模型是边坡冻融稳定性分析的重要前提，本文选取D-P(Drucker-Prager)准则匹配Mohr-Coulomb条件。如式(5)和式(6)所示，在平面应变情况下，当采用关联流动法则时，此准则成为DP4准则，较好的克服了Mohr-Coulomb屈服面存在棱角不好收敛的问题。

(5)

(6)

式中：α、k是土体内摩擦角与黏聚力有关的常数；c、φ是土体黏聚力和内摩擦角。

1.3 基于有限元的热-力耦合多年冻土边坡稳定性计算方法

冻融作用下多年冻土边坡稳定性是一个复杂的动态平衡问题，极端的气候条件使得边坡土体性质具有明显的时效性，并且该地区边坡有着整体稳定性较好，坡体浅层易受极端气候反复冻融的影响，而永冻层相对较稳定的特点。基于此，笔者考虑土体冻融过程中，环境温度及冻融循环对边坡土体的综合作用，利用多物理场有限元软件强大的函数编辑功能，编写了热-力耦合动态if函数，即利用温度来界定冻融交界面，以此解决在数值计算中难以考虑局部冻融的问题。

强度折减法有着不需假设滑动面以及能够适应各种复杂情况的优点。就边坡而言，其本质就是在有限元计算中，使得岩土体的抗剪强度参数黏聚力以及内摩擦角不断折减同样的倍数，直至最终岩土体达到临界的破坏状态，此时的折减系数即为边坡的稳定性系数[11-12]，即：

cF=c/F
φF=tan-1(tanφ/F)

(7)

式中：cF、φF分别为折减后的黏聚力和内摩擦角；F为折减系数即边坡的稳定性系数。

值得注意的是，上述强度折减法是针对边坡整体稳定性而言的，对仅受到局部冻融的多年冻土边坡而言不能完全适用。本文是在热-力耦合动态if函数的基础上，利用强度折减原理仅对发生冻融作用的浅层土体不断折减其强度，以坡顶位置处位移发生突变作为边坡失稳破坏的判断标准，以此分析边坡稳定性以及计算其安全系数。

2 模型建立过程

2.1 工程概况

某一级公路位于我国东北大兴安岭内蒙古境内，整体呈东西走线。根据气象数据得到该区年平均气温为-3℃，最大冻结深度为3.1 m，根据中国季节性冻土标准冻深线，该地区属于多年冻土区。根据课题组在该段钻孔实测数据显示地表大多覆盖0.5 m～2.0 m的草炭土，以下多为粉质黏土、碎石土以及砂砾。人工边坡坡高度在8 m～15 m之间，坡率大多集中在1∶1.1～1∶3之间，土质边坡为黏性土，这里选取坡率为1∶1.2、坡高为15 m的土质边坡进行分析。

2.2 冻融循环对土体强度的影响

(1) 对抗剪强度的影响。抗剪强度指标在边坡稳定性分析中的基本参数。目前大量的研究成果表明，土体在长期的冻融作用下，其强度都有不同程度的损伤[13-14]。文献[14]给出了土体在冻融作用下，土体抗剪强度指标损伤系数与冻融循环次数的关系：

(8)

式中：kc、kφ分别是黏聚力和内摩擦角的冻融损伤系数；n为冻融循环次数，这里指时间周期，单位为a。则得到经过冻融循环后的抗剪强度参数为：

c′=ckc,φ′=φkφ

(9)

式中：c、φ分别是土体初始黏聚力和内摩擦角；c′、φ′分别为土体损伤后的黏聚力和内摩擦角。

(2) 对弹性模量的影响。弹性模量是岩土工程中变形和稳定性分析中的一个重要参数。文献[15]以轴向应变为2.0%时所对应的应力作为土体弹性模量受冻融作用影响的标准，得到土体的弹性模量随冻融循环次数先降低而后增大规律。文献[16]针对东北地区典型黏性土通过不固结不排水三轴试验，得到了黏性土应力-应变关系。本文选取轴向应变为2.0%时所对应的的偏差应力增量与轴向应变增量的比值作为土体的弹性模量，即：

(10)

式中：E表示土体的弹性模量；σ、ε表示在三轴试验时的应力与应变，脚标表示该应变时的应力值。

根据上式计算得到黏性土在经历冻融循环后的弹性模量的变化趋势，计算结果如图1所示。可看出黏性土在经历冻融循环13次后弹性模量都发生了不同程度的下降，前两次的冻融对土体弹性模量损伤幅度最大，而后随着冻融循环次数的增加，弹性模量表现出先增大后减小的趋势。总体而言，弹性模量的平均降低幅值为未冻融土弹性模量的38%。

图1 黏性土弹性模量随冻融循环次数变化曲线

2.3 模型建立

(1) 几何模型。通过对该段某处边坡实际形态的简化处理建立边坡分析模型，边坡高度为15 m，坡率1∶1.1，重度17.8 kN/m3，初始黏聚力39.3 kPa,初始内摩擦角20.9°，初始弹性模量20 MPa，泊松比为0.3，几何模型如图2所示。

图2 边坡几何模型(单位:m)

(2) 参数选取。冻融作用对土体强度的影响如2.2节所述，土体冻结时力学参数抗剪强度、泊松比以及弹性模量与土温的关系如式(11)所示[17],其热力学参数参见表1和表2。

E=a1+b1|T|m
v=a2+b2|T|
c=a3+b3|T|
φ=a4+b4|T|

(11)

式中:ai、bi是试验参数，各个参数的取值如表1所示：在土体温度T>0℃时，bi均为0；m取0.6。

表1 土体力学试验参数[18]

表2 土体热物理参数

(3) 边界条件。气候环境对边坡土体的作用主要取决于边坡表面与大气的对流换热过程。根据近60年东北地区气候变化趋势特征，预测东北地区年平均气温未来升温速率为0.05℃/a，在考虑大气逐渐升温的情况下，气温和地表温度随时间变化如公式(12)所示[19]：

(12)

式中:t为时间，单位s；Ta为气温；Ts为地表温度；Tm为年平均气温；A为年平均气温较差；g(t)为考虑全球变暖升温项。

(4) 初始条件。初始温度场对整个坡体状态的影响较为重要，将年平均地温作为整个坡体初始温度场计算的初始值，把g(t)=0时的天然地表温度作为整个坡体表面的温度边界条件，最后得到热-力耦合阶段的初始温度场。

为消除自重对边坡位移的影响，应首先平衡地应力，排除由土体自重对坡体应力及位移的影响，把此时的应力状态作为热-力耦合阶段的初始应力场。

2.4 天然温度场验证

根据实测钻孔数据，该段天然地层从上依次为0.3 m草炭土、2.2 m粉质黏土、2.2 m碎石土以下为砂砾，按照实际地层情况建立80 m×40 m的矩形几何模型，以此模拟天然地表地温变化情况。限于篇幅，几何模型不再展示。

图3所示为天然地表钻孔实测温度与数值模拟值地温对比曲线，由图3可看出，模拟值与实测值随深度的变化曲线吻合度较高，在7月30日和11月18日两个时间点除了在接近地表位置有所差异(可能是接近地表位置测温钢管内的温度受到外界空气的影响较大所致)，在冻土上限以下，基本一致，证明了本模型的适用性及合理性。也得到该段年平均地温为-1℃，年变化深度在地表以下9 m处。

3 多年冻土边坡稳定性分析

3.1 边坡冻融演变过程

(1) 温度场演变过程。随着气温对土体周期性的作用，坡面季节活动层随之经历着反复的融化和冻结，坡体温度场演变过程中，融化深度呈现出先增加后减少的趋势。在5月初出现浅层融化现象直至9月份达到最大融化深度，而后进入冷季，气温的下降，冷量向坡体传递，原来的融化层逐渐开始缩小，直至次年1月末融化层完全消失，之后坡体一直处于冻结状态。如果将坡体融化深度增长的过程认为边坡处于融化期，坡体融化深度减小的过程认为边坡处于冻结期，则边坡融化期为5月—9月份，冻结期为10月—次年4月份。通过“0.00℃”等值线位置可以明显看出坡体浅层融化夹层的存在，在10月份出现到次年2月份坡体中部会形成保持长达5个月的融化夹层，融化夹层在10月份出现并逐渐发展直至次年2月份消失。

图3 温度实测值与模拟值对比曲线

以坡顶位置作为特征点，提取坡顶特征点以下的融化层厚度，得到边坡在一个温度周期的融化层厚度随时间的变化曲线见图4。由图4可看出，边坡融化层厚度在温度周期内和环境温度变化呈现出很强的相关性，并随着环境温度的变化而变化，但与环境温度在时间上存在一定的滞后性。环境温度在4月份到7月份是升温阶段，坡体表层从5月份开始逐渐融化直到10月份达到最大融化深度2.9 m；而后环境温度下降，坡体出现融化夹层，坡体内部融化层厚度慢慢减小，随之进入冷季，到次年2月份融化夹层完全消失；在此期间坡体融化层厚度经历先增大后减小最终完全消失的变化趋势。总体上，坡体融化层厚度与环境温度在时间上存在约1个月的滞后期。

图4 边坡融化深度与环境温度随时间变化曲线

(2) 应力及位移场分析。冻融作用下，坡体应力场分布表现出明显的时效性。在整个温度周期内最大压应力都出现在边坡底部；在边坡融化期，拉应力区集中在边坡浅层，并且随着温度的升高，拉应力区不断增加；在冻结期，边坡拉应力区随着温度的降低不断减少。而冻融循环对边坡浅层土体强度的损伤作用，是边坡在融化期发生失稳破坏不可忽视的原因。

图5为坡顶特征点竖向位移和水平位移随时间发展的变化曲线。从位移变化情况看，在整个循环周期内，特征点竖向位移和水平位移随时间变化趋势基本一致，整条曲线呈“台阶”状，可分为快速发展阶段和稳定发展阶段。具体的，在温度周期循环下，随着气温的逐渐升高，季节活动层开始融化，土体融化使得含水量逐渐增加，导致土体强度的急剧下降，此时位移出现第一个快速发展阶段；当边坡处于冻结期，活动层逐渐冻结，土体在冻结状态时的强度远远大于融化状态时的强度，并且随着温度的降低在不断增大，在此期间特征点位移进入稳定发展阶段并且竖向位移出现轻微的回升。如此反复特征点位移在不断增加，以位移发生突变作为判定边坡发生破坏的标准，从图中可看出，在经历8个周期循环后位移发生突变，此时数值计算也出现了不收敛的情况，可以判定此时边坡达到了破坏的临界状态。

图5 坡顶位移随时间变化过程

3.2 安全系数的讨论

依据第1节所述多年冻土边坡稳定性分析方法以及3.1(2)所述，可认为在第9年即经历8次冻融循环后边坡达到失稳的临界状态。即可根据强度折减原理得到该边坡安全系数为1.494。

(1) 滑动面的确定。如图6所示为边坡安全系数F=1.494时的位移等值线分布图，可以明显看出有位移突变点的存在，图中粗实线为位移等值线沿着突变点位移增大的方向形成的实线，该方向上的特征节点全部都发生了突变，此时坡体已经是失稳状态，故此可以以该方向作为边坡的最危险滑动面。

图6 F=1.494时的位移分布图

(2) 对比极限平衡法。极限平衡法根据不同的假设条件有如表3所示的几种方法。本节利用GEO-SLOPE极限平衡法方法来验证有限元方法的合理性，计算模型如图7所示。因冻土边坡的失稳破坏常常发生在暖季，所以可仅考虑边坡处于融化期的状态。考虑冻融影响，将有限元计算得到的边坡处于临界状态时的土体强度参数作为边坡融化层的参数，融化层厚度按照边坡最大融化深度设置；永冻层的力学参数按照冻结状态设置。计算结果如表3所示，可看出极限平衡法的5种方法与本文所得安全系数相差不大，证明了本文方法的准确性。