基于强度折减法和冻融界面理论的块石护坡铁路路基稳定性分析

2020-10-28侯彦东叶泽国周凤玺

水利与建筑工程学报 2020年5期

侯彦东,叶泽国,周凤玺

(兰州理工大学土木工程学院，甘肃兰州 730050)

多年冻土区道路工程往往具有重要的经济、政治及战略意义，其所涉及的冻土问题一直是国内外学者研究的热点和难点。青藏铁路自建成以来，已历经16～17个冻融循环[1]，其下伏土层的温度场已经发生了质的改变。大量研究表明，块石结构能有效降低路基下伏土层地温，抬升多年冻土上限。然而，块石结构也会由于其自重大，改变路基附加应力的作用范围，进而对下伏土层的承载性提出考验。

单纯利用极限平衡方法计算冻土边坡或冻土路基的稳定性结果往往并不十分理想。武鹤等[2]利用楔形极限平衡法，考虑滑动面端摩阻力的影响推导出了边坡安全系数的公式，同时又讨论了边坡融化深度、土体黏聚力、内摩擦角、饱和重度等与安全系数之间的关系。但其并没有考虑融化层土体力学性质的改变等影响因素。李永亮等[3]通过有限元强度折减法和传统极限平衡法的对比，利用数值模拟的方式，探究了三种失稳判断标准的适用性。但其只是考虑了在判别是否失稳时的标准，并未从根本上探究失稳的原因。还有学者采用了强度折减法，通过不断对土体的黏聚力和内摩擦角进行强度折减或者通过不断对土体尝试加载以使得边坡达到极限状态[4]。Pei等[5]对普通路基和块石护坡路基在四个不同季节的应力应变状态进行了计算分析，同时给出了两种路基在不同季节的安全系数计算结果。何彬彬[6]利用双强度折减法分析了多年冻土路基的安全稳定性，探讨了路基高度与融化盘位置的关系，研究了安全系数的变化规律，其中，融化盘位置是根据温度场数值计算结果估测的。

近年来，冻融界面这一概念越来越被国内外学者所重视，对于冻融界面的土体参数以及对土体稳定性的影响程度也成为各位学者一直探究的问题。当温度升高时，冻土内部部分土体融化，融化土与冻土之间会形成一个固液交互面，此界面即为“冻融界面”，高樯等[7]在对冻土斜坡失稳的研究过程中针对冻融界面进行了直剪试验，发现冻结区和融化区物理参数相差甚远，而冻融界面的性质成为了影响冻土边坡失稳的关键。彭慧等[8]分别探讨了圆弧法和折线法对青藏高原多年冻土公路路基稳定性评价的适用性，虽提及冻融界面，但未赋予其相应力学界面参数。褚志成等[9]在研究多年冻土边坡在冻融作用下的稳定性时考虑了土体冻结和融化情况下土体力学性质的不同，但其并没有设置单独的冻融界面参数。

本文在计算块石护坡路基稳定性时，参考了以上学者的研究方法的优点及不足，在计算路基稳定性时，采用了折减系数法与冻融界面理论相结合的方式，并且单独设置冻融界面参数，使得稳定系数在计算过程中更加准确，能够更加清晰地描述多年冻土区铁路路基在运行过程中的稳定性情况。

1 基本理论

1.1 路基的应力应变关系

假设路基土体为各向同性的线弹性材料，并且不考虑水分、温度等因素引起路基应变。因此，静力平衡方程可表示为：

0=▽·S+Fv

(1)

式中：Fv为自重体力；S为总应力；S总应力可表示为：

S=Sad+C·ε

(2)

式中：Sad为初始应力；C为弹性刚度矩阵，C=C(E，υ)；ε为附加应力引起的应变，可表示为：

(3)

式中：u为路基产生的位移。

1.2 数值模型的屈服及强度准则

岩土力学理论分析常使用Mohr-Columb准则作为强度准则，但因在三向应力状态下，摩尔-库仑破坏面具有角隅性质，即应力落在破坏面的尖角处，导致导数的方向不定[10]。针对这一问题，1952年Drucker与Prager对三向应力状态下的Mohr-Columb准则函数进行了改进，提出了D-P准则，弥补了摩尔-库仑准则的局限性。

赵尚毅等[11]对基于D-P准则的边坡安全系数定义及其转化进行了深入研究和对比，建立了基于D-P准则的边坡稳定安全系数与传统Mohr-Columb条件下安全系数的关系表达式。在COMSOL中将路基变形视为在一定范围内弹性变形，但大部分为塑性变形。因此，本文选用D-P准则匹配Mohr-Columb准则作为路基强度准则。

Mohr-Columb强度准则表达式为[10]：

(4)

D-P强度准则表达式为[11]：

(5)

式中:α、k为与岩土材料内摩擦角φ和黏聚力c有关的常数。

根据岩土力学中Mises屈服面与塑性流动法则，在二维平面应变的条件下，D-P准则中的材料参数可以由材料的黏聚力和内摩擦角得出，因此：

(6)

(7)

1.3 冻融界面性质

多年冻土区活动层含水丰富，冻结期普遍存在丰富的地下冰晶体，融化期时部分冰晶融化为孔隙水[7]。因为冻结锋面的下伏冻土不透水，融化的土石混合体活动层通常为高含水状态。在冻土边坡活动层开始融化时，会形成一个特殊的冻-融土交界面，土石混合体作为一种特殊的工程地质体[12]。在冻结区，土体的力学性质类似于岩石，强度极大，而在正融或正冻状态的剧烈相变区内冻土的力学强度因冰水相变作用而发生较大的变化，即冻融交界面为冻结岩土体与融化或正融岩土体力学性质变化的分界[13]，往往是物理力学特性发生突变的位置。这个特殊交界面的位置、形态、附近土体温度及其参数往往是影响多年冻土区土石混合体斜坡稳定性的关键。因此，本文在计算路基稳定性时，采用了强度折减法，将冻融界面考虑到稳定性计算过程当中，试图提高多年冻土路基边坡安全系数计算结果的准确性。

1.4 强度折减法原理

强度折减有限元法是在1975年提出的，之后强度折减法对路基的安全稳定性计算在处理边坡问题中被广泛应用。强度折减法是在寻找土体的极限状态时，通过工况的实际情况不断地尝试对土体的黏聚力和摩擦角进行折减，使土体达到破坏状态[14]。具体过程为：将土体黏聚力和内摩擦角进行折减，再将它们导入模型中进行计算。经过多次计算直至计算结果不收敛，最终求得的该极限状态下的折减系数即为此时工况的安全系数。

强度折减法的具体公式为：

(8)

式中：c为强度折减前土体的黏聚力；c′为强度折减后土体的黏聚力；F为路基的安全系数。

(9)

式中：φ为强度折减前土体的内摩擦角；φ′为强度折减后土体的内摩擦角；F为路基的安全系数。

2 几何模型确定过程

由于本文计算安全稳定系数时考虑了冻融界面理论，所以除了基本的建模过程外，还需要进行多软件联合操作，来实现强度折减法和冻融界面理论相结合的计算过程。具体方法如下：

(1) 将FLUENT计算完成的路基温度场进行后处理，以文件的形式导入到Surfer绘图软件中，根据温度场数据建立网格，画出温度等值线图，准确定位0℃等温线，即冻融交界面。

(2) 将在Surfer中找到的0℃等温线保存为dxf格式，以便于在AutoCAD中打开，同时，删减掉多余的等温线及模型边框，将0℃等温线，也就是冻融界面，将其竖直上下偏移一定距离(本文取30 cm)便能得到的冻融交界带。

(3) 将建立好冻融带的AutoCAD文件另存为dxf格式，之后在COMSOL软件中将另存的dxf文件打开，与提前建立好的二维路基模型形成联合体；在COMSOL软件中，建立不同变量，分别设置各个土层和各个土层冻融界面参数，赋予土层参数和处于不同土层中冻融界面的参数。

3 几何模型及参数设置

本文将Hou等[15]对年平均地温为-0.5℃的块石护坡路基在RCP4.5排放情景下运营至第5年、20年和50年10月15日(此时下部的土层达到最大融化深度)温度场计算结果为块石护坡路基力学参数的计算基础，如图1所示。从-0.5℃等温线的演化过程可以看出，块石护坡路基从5 a运营至20 a期间，块石护坡结构可有效消散路基布设期间引起的热扰动，并降低其正下方多年冻土温度，20 a时路基热稳定状态明显优于5 a时的情形。而对比20 a和50 a的温度场可以发现，路基运营至50 a时天然上限和路基下部的人为上限均由于气候变暖的影响而明显下降，下伏冻土也增温明显，热稳定状态也明显弱于20 a时的情形。

图1 高温冻土区块石护坡路基温度场[10]

将温度场计算结果导入COMSOL有限元软件，之后把温度与参数的关系式以变量的形式编写入COMSOL有限元软件计算模型中，使得模型能够自动根据各个土层中点的温度选取参数进行计算。参数的设置主要分为两个部分：土体参数，冻融界面参数。其中路基下部土层从地表处往下依次为砂土(厚度2 m)、黏土(含土冰层，厚度4 m)和强风化泥岩(未见底)，温度与冻土力学参数的关系式参考文献[16]，如下：

E=a1+b1|T|m

(10)

v=a2+b2|T|

(11)

c=a3+b3|T|

(12)

φ=a4+b4|T|

(13)

式中：m为弹性模量非线性指数取0.6；E为土体的弹性模量；v为土体的泊松比；c为土体的黏聚力；φ为土体的摩擦角；a1—a4，b1—b4为试验常数，土体处于融化状态时b1—b4=0。土层的力学参数如表1所示。主要通过公式(10)—公式(13)和表1的参数设置路基填土、砂土、黏土和强风化泥岩等土层的力学参数。

基于高樯[7]的研究成果确定了冻融界面的参数设置，其中内摩擦角和黏聚力如表2所示。根据表格中的参数，砂土层选择含水率13%时的界面参数，粉土层选取含水率为25%时的冻融界面参数，黏土层选取含水率为25%时的冻融界面参数进行计算。

表1 各层土体的试验参数

表2 冻融界面物理参数

4 荷载及边界条件设置

荷载和边界条件的设置可以引起应力应变计算结果的变化，荷载和边界的类型也有很多种。本文对路基稳定性的计算中，不考虑温度场和水分场的变化，故边界条件都设置为力学边界。在已有研究中，其他的专家学者在研究块石护坡路基稳定性时，往往只考虑块石护坡施加后对路基温度场的影响，而忽略了块石护坡自身重力对路基所产成的力学效应。本文不仅考虑块石护坡对路基产生的热学影响，同时考虑块石护坡因自身重力产生的对路基和下部土层产生的力学影响，进而探究块石护坡对铁路路基的综合效果。本文首先依据花岗岩的密度及多孔介质的孔隙率计算了块石结构的等效重度，在块石护坡铺设的位置按照不同位置处块石结构自重应力荷载大小不同而分段设置了线性荷载，用来等效块石护坡对路基和下部土层的荷载。坡脚处荷载长度为1.5 m，大小从0 N/m到23 520 N/m，边坡上第一段为均布荷载，其大小为23 520 N/m，之后逐渐减小为0 N/m，路基上部为竖直向下的48 031.08 N/m的铁路荷载，如图2所示。

图2 块石护坡路基荷载分布

5 结果分析与讨论

5.1 块石护坡路基位移场对比分析

根据之前学者对冻土边坡稳定性的研究，笔者以文献[15]中所计算的路基温度场为基础，计算了高温冻土区块石护坡路基运营第5年、第20年和第50年的安全系数，结果如表3所示。可以看出，考虑冻融界面情况下路基边坡安全系数随着铁路的运营时间的增长有先增加后减小的规律，而不考虑冻融界面情况下路基边坡安全系数似乎与路基的运营时间关系不大，均维持在2.8左右的水平。此外，还将考虑冻融界面理论和不考虑冻融界面情况下的路基的总位移及X方向位移分量进行了对比分析，如图3—图8所示。

表3 路基运行5 a、20 a和50 a的安全系数

对比图3和图4中总位移在X方向的分量可得，在不考虑冻融界面时，强度折减至极限状态下，路基的内部出现一个从路基中心向左延伸至坡脚的潜在滑移面，即潜在滑移面出现在了路基内部，而未出现在冻土上限附近冻融界面处，路基潜在破坏模式为整体破坏。对比图4中考虑冻融界面后X方向位移分量和总位移，可以发现，在考虑冻融界面后，强度折减至极限状态下，最大的位移差的位置出现在路基下部黏土层(含土冰层)中冻融界面处，路基潜在破坏模式为局部破坏，这与前人对冻土边坡失稳机理的诠释[17]是一致的。同时，考虑和不考虑冻融界面情况下路基的安全系数计算结果也差异较大，分别为1.05和2.81。由此可见，冻融界面的存在对路基的安全稳定性分析结果有着至关重要的影响。

图3 不考虑冻融界面时路基运行5 a时极限状态下

图4 考虑冻融界面后路基运行5 a时极限状态下

同样的，从图5、图6可得，在路基运行20 a时，采用强度折减法折减至极限状态下，路基及下部土层中X方向位移分量和总位移规律与运行至5 a时基本相同，路基及地基土的潜在破坏模式也基本相同。然而，考虑冻融界面情况下，路基下部土体潜在破坏时，X方向的位移大小和局部破坏的位置(最大位移差发生的位置)也相比于5 a时的情形也发生了变化。从表3可以看出，20 a时路基安全系数也从5 a时的1.05增至1.31，这与块石护坡在5 a～20 a期间对路基下伏土层温度场的改善(上限位置的调整及上限附近土体温度的降低)有很大的关系。

图5 不考虑冻融界面时路基运行20 a时极限状态下

图6 考虑冻融界面后路基运行20 a时极限状态下

由图1可知，当路基运行至第50年，天然场地冻融界面下降至-5 m，路基中心位置下降至-3.2 m左右。由图7来看，在不考虑冻融界面时，潜在滑移面仍出现在路基内部，从路基中心向左延伸至路基坡角处，最大位移差能够达到0.16 m。而从图8来看，考虑冻融界面情况下，强度折减到极限状态时，X方向位移分量在路基坡脚正下方-4 m以上范围内，位移差有一个明显的突变过程，最大值达0.04 m，同时，冻融界面附近，多处都出现了数值不同的位移突变点。纵观图4、图6和图8中考虑冻融界面时总位移云图，在运行的第50年，考虑冻融界面情形下，潜在局部破坏位置由单点转变为多点且存在逐渐向路基中心迁移的趋势，同时，路基边坡安全系数也由20 a时的1.31降低至0.85。这是由于当路基运营至50 a时，路基下部冻融界面所处土层已经由沙土层为主转变为黏性土层为主。由表2、表3及文献[7]可知，当黏性土为含土冰层时，冻融界面处强度参数相比于冻结区衰减非常显著，而当砂类土为含土冰层时，冻融界面处强度参数几乎不变。此外，路基运营至50 a时，冻融界面的规模、形态的发展变化也对该计算结果产生了显著影响。

图7 不考虑冻融界面时路基运行50 a时极限状态下

图8 考虑冻融界面后路基运行50 a时极限状态下

5.2 安全系数计算结果可靠性讨论

本文分别对高温冻土区块石护坡路基考虑和不考虑冻融界面两种情况下的安全系数进行了计算，结果表明，考虑和不考虑冻融界面情形下路基边坡的安全系数差异较大，破坏模式则分别呈破脚处上限位置黏土层(含土冰层)中的局部破坏和路堤填土整体破坏；考虑冻融界面情形下，路基边坡安全系数和人为上限附近土层温度、所处位置、形态及路基下部温度场关系密切，而不考虑冻融界面情形下路基边坡安全系数和路基下部温度场似乎关系不大。

由于冻融界面准确定位、不同位置土层强度参数准确划分，运行至第5年时，由于路基铺设时对地基温度场产生了热扰动，上限附近温度较高，安全系数仅为1.05；运行至20 a时由于上限附近热稳定性改善，安全系数增至1.32；运行至50 a时，由于上限下降，原上限附近冻土部分已经融化，说明此时路基已处于热不稳定状态，但在该状态下不考虑冻融界面情形下安全系数计算结果仍高达2.78，反观考虑冻融界面的情形，此时潜在破坏面全部出现在黏土层(含土冰层)中，安全系数降低至0.84，说明考虑冻融界面情形下，强度折减的极限状态的位移云图能很好反映路基热稳定性降低对路基边坡稳定性造成的不利影响。

总体来看，考虑冻融界面情形下，路基边坡安全系数与路基热稳定性呈正相关，计算结果较为合理。此外，过去多项研究表明，冻融斜坡滑动面多位于黏土与上限附近的冰层的交界面[18]，这也与本文考虑冻融界面情形的研究结果相契合。因此，在作多年冻土路基边坡评估时(特别是当冻融界面出现在粘性土层中时)应明确冻融界面力学参数、准确定位人为上限所处位置、形态、及其附近冻土温度。