基于未冻水膜压力判据的冻土水热力耦合冻胀模型研究

王自力; 李金峰; 滕继东; 张升; 盛岱超

doi:10.11779/CJGE20220227

基于未冻水膜压力判据的冻土水热力耦合冻胀模型研究 English Version

王自力^1,,
李金峰¹,
滕继东^{1, 2, ,},
张升^{1, 2},
盛岱超^{1, 2, 3}

1.
中南大学土木工程学院, 湖南长沙 410075
2.
高速铁路建造技术国家工程研究中心, 湖南长沙 410075
3.
悉尼科技大学土木与环境工程学院, 悉尼 2007

基金项目:

国家自然科学基金项目 52178376

国家自然科学基金项目 51878665

国家自然科学基金项目 U1834206

湖南省湖湘青年英才项目 2020RC3008

湖南省研究生科研创新项目 CX20220110

中南大学研究生自主探索创新项目 2022ZZTS0019

详细信息

作者简介:
王自力(1998—)，女，博士研究生，主要从事土体冻胀机理及数值计算方面的研究工作。E-mail: wzl1126@csu.edu.cn

通讯作者:
滕继东, E-mail: jdteng@csu.edu.cn

中图分类号: TU445
计量
- 文章访问数: 386
- HTML全文浏览量: 87
- PDF下载量: 119
出版历程
- 收稿日期: 2022-03-03
- 网络出版日期: 2023-05-18
- 刊出日期: 2023-04-30

THM coupled model for simulating frost heave based on a new water film pressure criterion

WANG Zili^1,,
LI Jinfeng¹,
TENG Jidong^{1, 2, ,},
ZHANG Sheng^{1, 2},
SHENG Daichao^{1, 2, 3}

1.
College of Civil Engineering, Central South University, Changsha 410075, China
2.
National Engineering Research Center of High-Speed Railway Construction Technology, Changsha 410075, China
3.
School of Civil and Environmental Engineering, University of Technology Sydney, Sydney 2007, Australia

摘要

摘要: 冻胀融沉是多年冻土区的主要病害，其受水分场、温度场、应力场的复杂耦合作用影响。基于水膜理论提出了冻土未冻水膜压力作为冰透镜体生成的判据，并重新对水分迁移驱动作用进行描述，建立了以温度、土体孔隙比为变量的全耦合模型。通过考虑已冻区冰基质的影响，推导了涵盖原位冻胀与冰分凝两部分的冻胀量计算公式。基于Matlab和COMSOL Multiphysics的联立平台，提出了模型冰透镜体实时分布的数值求解方法，实现了冻土温度、水分、应力、冰透镜体分布的全耦合数值求解。通过与室内土柱冻结试验及现有水热力模型（Thermal-Hydraulic-Mechanical，又称THM模型）冻胀量结果进行对比分析，证明了温度、含水率与冻胀计算上的可靠稳定。最后通过探讨温度梯度、上覆压力、渗透系数与压缩模量对土柱冻结的影响发现，温度梯度能显著增加土体冻胀量，上覆压力会导致更多水分向冻结锋面迁移，但对冻胀量起着抑制作用，渗透系数与压缩模量均对冻胀量产生正影响。为冻胀理论研究与数值实现提供了新思路。
- 冻土 /
- 水热力耦合 /
- 冰透镜体 /
- 水膜理论 /
- 冻胀融沉
Abstract: The frost heave and thaw settlement are the main frost damage in cold areas, which are the complex coupling process of water, temperature and stress fields. In this study, a coupled thermal-hydraulic-mechanical model is developed based on the water film theory, in which the temperature and void ratio of soils are the input variables. The novelty of this model is that the frozen water film pressure is used as the criterion for the generation of ice lens. The driving force of water migration is newly defined, and the frost heave includes the pristine frost heave and the amount of ice segregation. The fully coupled model is numerically solved based on the Matlab and COMSOL Multiphysics, generating the results of soil temperature, moisture, stress and the layered ice lens. The simulated results are then compared with those of the laboratory freezing tests, which shows that they match quite well and verify the validity of the proposed model. The simulation indicates that temperature gradient can promote the frost heave, and the overburden pressure can attract more water to the freezing front but decrease the amount of the frost heave. In addition, both the hydraulic conductivity and the compressive modulus have positive effects on the frost heave. The proposed model provides a new approach to understand the frost heave.
- frozen soil /
- thermal-hydraulic-mechanical coupling /
- ice lens /
- water film theory /
- frost heave and thaw settlement

HTML全文

0. 引言

尽管对于大型滑坡的启动和运动机理的认识仍然存在争议^[1]，但是滑动面的摩擦强度衰减是大型滑坡发生高速运动的主要诱因已成为广泛共识^[2-3]。滑动面的摩擦强度在滑坡发展过程中是一个动态演化的参量，呈现出显著的速率相关性，这种相关性决定了滑坡的滑移模式、滑动速度和移动距离^[4-5]。大型滑坡的勘探资料和反演分析表明，滑体高速滑移时滑动面呈现出远低于其正常摩擦强度的超低摩擦性^[6-7]，典型的如日本的Nikawa滑坡^[8]和中国的西藏易贡滑坡^[9]。滑带土的残余强度决定了滑动面的摩擦强度，因此，室内环剪试验和旋剪试验能定量地给出滑动面摩擦强度随速率的衰减程度，揭示高速滑坡的触发机制^[10-15]。试验研究表明粗粒滑带土的残余强度随剪切速率的增加呈现减小的规律^[10-13]；且随着剪切变形量和剪切速率的增加，粗粒滑带土内部颗粒破碎、细化，孔隙水压力逐渐累积，并可能触发剪切面液化现象，导致残余强度显著下降^[14-15]。这种剪切速度增加和滑动面上强度衰减的相互促进机制，促使滑体不断加速，并最终发展成为高速远程滑坡。

但是需要指出的是，剪切面摩擦强度衰减对大型滑坡的失稳启动的影响研究仍然较少。目前，临界位移是判别滑坡失稳一个常用指标^[16-17]。Jibson^[18]的研究表明大多数土体强度在5~10 cm的剪切变形下可降低至残余值，据此建议了临界位移的参考值。但是在一些地震滑坡现象中发现，滑体位移超过1 m仍未发生破坏（Dong等^[19]、Yang等^[20]）。传统的临界位移没有反映滑动面上的强度速率衰减效应，无法反映大型滑坡的启动机理，更无法考虑地震荷载与强度速率衰减效应对滑坡启动的同步作用。因此，在地震作用下，滑动面摩擦强度是如何衰减并对大型滑坡失稳和启动产生影响是需要进一步研究的^[20]。

数值模拟方法是研究滑坡灾害机理，找出滑坡启动和运动规律的有效手段。离散单元法（DEM）^[21]、非连续变形分析（DDA）^[22]、光滑粒子流体动力学方法（SPH）^[23-25]和物质点法（MPM）^[26-27]等方法已经广泛用于模拟大变形的滑坡问题当中。但是目前的数值模拟研究主要关注大型滑坡在滑体高速滑动时的超低摩擦现象，较少考虑地震荷载与摩擦面强度衰减效应在滑坡失稳启动中的影响规律。为此，本文基于SPH和FEM耦合的方法，发展能够在滑坡启动和滑移过程中考虑滑动面摩擦强度速率衰减效应影响的动力数值方法，通过实际地震诱发滑坡的模拟，再现地震作用下大型滑坡启动和运动过程，探讨滑动面摩擦强度的速率衰减效应对大型地震滑坡启动的影响机制。

1. SPH和FEM耦合的数值计算方法

1.1 数值计算方法流程

有限元法（FEM）采用连续介质力学原理描述，能够适应复杂的几何形状和边界条件，通过适当的边界处理，能够合理考虑地震动的输入，降低地震波在边界上反射产生的误差。但是，FEM方法在处理大变形问题中有一定的局限性，存在网格畸变问题。光滑粒子流体动力学法（SPH）采用拉格朗日方法描述，物体的变形及运动信息携带在粒子上，广泛地应用于模拟流体和固体变形问题，能够较好地模拟大变形现象。但是SPH方法在边界处理上相对复杂，需要在边界上人为设定初始速率、荷载或约束，这种设置难以考虑复杂的边界条件，也很难完全再现真实的边界条件，更适用于滑坡发生后的滑体运动过程的模拟^[23-25]。

本文拟采用SPH和FEM耦合的数值计算方法，建立一种能够描述滑坡由初始稳定状态、启动状态到滑移破坏状态全过程的模拟方法。该方法能够合理考虑复杂边界条件的影响，也能够合理再现滑坡中的滑体的大变形发展过程。所建立的耦合计算方法的流程图如图 1所示。

图 1 SPH和FEM耦合的数值模拟方法流程图

Figure 1. Flow chart of coupled SPH and FEM numerical simulation method

下载: 全尺寸图片幻灯片

（1）确定滑体和滑床间的滑动面位置。一些滑坡可以根据现场地勘直接确定。对于现场难于确定的情况可以由有限元的强度折减法计算确定滑动面的具体位置。

（2）建立SPH和FEM耦合的几何模型。首先，将滑坡的几何模型导入有限元网格进行网格剖分，其中滑动面设置为接触单元。然后，将滑体的有限元网格转化为SPH粒子，并在滑床的有限元网格的外侧增加设置无限元网格，以模拟远场地基对地震波能量的吸收作用以及无穷远处位移为零的边界条件。

（3）设置不同部分的本构模型并输入材料参数。需根据不同部分岩土材料的性质，分别定义滑体和滑床的本构模型以及滑动面的接触模型。

（4）设置计算步骤并分别定义边界条件和荷载条件。根据地震滑坡前后的荷载变化设置两个计算步骤：重力作用下的静力计算步骤和地震作用下的动力计算步骤。在动力计算中需要设置地震荷载的作用。

（5）开展静力计算分析，确定重力作用下滑坡内部的应力和应变，并将计算结果存储。

（6）开展动力计算分析，将静力计算得到的应力应变作为初应力和初应变，进一步开展动力计算分析，模拟地震荷载作用下滑坡启动和滑移的过程。

1.2 滑床和滑体的本构模型

考虑外部降雨、荷载作用、地震及滑体运动的冲击作用，均会导致岩土体发生结构损伤，在本文的静力和动力计算步骤中滑床和滑体均采用强度软化的Mohr-Coulomb模型。假定强度参数会随着塑性变形的增加而降低，采用如下强度参数软化模型：

$c = {c_{\text{r}}} + ({c_{\text{p}}} - {c_{\text{r}}})\exp ( - \eta {\varepsilon _{{\text{eq}}}}) \text{，}$

(1)

$\varphi = {c_{\text{r}}} + ({\varphi _{\text{p}}} - {\varphi _{\text{r}}})\exp ( - \eta {\varepsilon _{{\text{eq}}}}) 。$

(2)

式中：c和φ为土体的黏聚力和内摩擦角；c_r为残余黏聚力； ${\varphi _{\text{r}}}$ 为残余内摩擦角；c_p为峰值黏聚力；φ_p为峰值内摩擦角；η为影响因子；ε_eq为塑性剪应变。当塑性变形很小时，黏聚力和内摩擦角分别为峰值黏聚力和峰值内摩擦角。该模型中，随着塑性变形的增加，黏聚力和内摩擦角逐渐减小至残余值，其中η的具体确定方法及取值参见文献[28]。

1.3 滑动面的摩擦强度速率衰减模型

滑动面采用Coulomb摩擦模型描述，如下：

$\tau = \left\{ \begin{gathered} \mu \sigma {\text{sign}}v{\text{ (}}\left| v \right| > 0{\text{)}} \hfill \\ \min (\left| {{\tau _{{\text{app}}}}} \right|,\mu \sigma {\text{sign}}({\tau _{{\text{app}}}})){\text{ (}}v = 0{\text{)}} \hfill \\ \end{gathered} \right. 。$

(3)

式中：τ，μ分别为滑动面上的剪应力和摩擦系数；v为剪切速度；σ为动面上作用的正应力；τ_app为 $v = 0$ 时，作用于滑动面上的剪应力；sign为符号函数，x > 0时sign(x)=1；x=0时sign(x)=0；x < 0时sign(x)=-1。

本文计算中摩擦系数μ采用Lucas等^[2]建议的考虑剪切速率影响的滑动面摩擦强度衰减模型：

$\mu (v) = \frac{{{\mu _0} - {\mu _{\text{w}}}}}{{1 + {{\left\| v \right\|} \mathord{\left/ {\vphantom {{\left\| v \right\|} {{v_{\text{w}}}}}} \right. } {{v_{\text{w}}}}}}} + {\mu _{\text{w}}} 。$

(4)

式中：μ₀和μ_w分别为静摩擦系数和稳态摩擦系数； ${v_{\text{w}}}$ 为特征速度。Lucas等^[2]广泛地分析了已有研究成果，指出式（4）中相关参数的取值范围为μ₀∈[0.5，0.9]，μ_w∈[0.08，0.16]， ${v_{\text{w}}}$ ∈[0.8 m/s，5.2 m/s]。Lucas等^[2]还进一步研究了太阳系内不同行星上发生的大型滑坡现象，发现当μ₀=0.84，μ_w=0.11， ${v_{\text{w}}}$ =4 m/s，在数值计算中可以较好地再现大型滑坡发生的过程，该组参数具有较广泛的适用性。

基于ABAQUS中的SPH模块进行了二次开发，通过用户子程序定义了相应的滑体和滑床的本构模型以及接触面的摩擦强度衰减模型。SPH算法中采用三次样条函数作为光滑核函数，时间积分方法采用中心差分法对运动方程进行显式的时间积分。SPH-FEM耦合法能充分利用两者的优势，在滑坡大变形区域（滑体）使用SPH方法，在小变形区域（滑床）采用FEM方法，边界采用无限元模拟，提高了计算效率的同时也提高了模拟的准确性。

2. 工程实例分析

2.1 唐家山滑坡工程概况

唐家山滑坡处于龙门山断裂带北川—映秀段西北约2.5 km，属于“V”型高山峡谷地貌。地震前唐家山斜坡上部坡度为39.01°，中部为34.65°下部为33.94°，属于中陡顺向岸坡结构，地震前天然状态下边坡整体稳定。滑坡顶部高程约1300 m，滑坡前后缘高差约540 m，水平距离约800 m，整个滑坡下滑时间约0.5 min，滑体方量约2.037×10⁷ m^3[29-30]，根据现场调研资料推测最大下滑速度约为28 m/s^[31]。现场地勘资料显示唐家山滑坡的地质剖面及滑体最终堆积形态见图 2，滑体结构从上至下分别为：强—中风化岩，厚4~10 m；弱风化灰岩，厚3~10 m；新鲜灰岩夹灰绿色含绿泥石的细粒状砱块岩、灰色含砱泥灰岩，厚约50 m，间夹层间剪切带和泥化夹层^[30]。

图 2 唐家山滑坡地质剖面图^[29]

Figure 2. Geological profile of Tangjiashan landslide^[29]

下载: 全尺寸图片幻灯片

2.2 数值计算模型

本文建立的唐家山数值计算模型如图 3所示。结合场地的实际条件，滑床基岩采用有限元网格以构建近场模型，基岩左右侧和下侧均设置无限元网格以构造远场条件，此外还设置了6个监测点，分别位于滑坡顶面和滑动面上的前段、中段和后段，用于记录监测点状态的变化过程。

图 3 汶川地震什邡八角镇地面加速度记录

Figure 3. Ground acceleration records of Bajiao Town, Shifang City during Wenchuan earthquake

下载: 全尺寸图片幻灯片

滑体和滑床材料强度衰减由式（1），（2）描述，材料参数见表 1，其中滑体和滑床的重度γ、弹性模量E、峰值强度c_p和φ_p根据曹琰波等^[21]的研究成果确定；残余强度c_r和φ_r根据唐家山滑坡（堰塞体）室内物性试验，进行地质工程经验类比分析确定^[31]；η根据Soga等^[28]的建议确定。滑动面采用Lucas等建议的速率衰减模型模拟，参数 ${\mu _0}$ ， ${\mu _{\text{w}}}$ 和 ${v_{\text{w}}}$ 也采用了Lucas等^[2]的建议取值： ${\mu _0}$ =0.84， ${\mu _{\text{w}}}$ =0.11， ${v_{\text{w}}}$ =4.0 m/s。此外，在地震作用下唐家山滑坡的破坏过程模拟中还需要输入地震波。本文采用什邡八角镇地面加速度实测数据（如图 4所示）作为输入地震波，地震水平和竖向加速度均施加在有限元网格的下边界处，如图 4所示。

表 1 唐家山滑坡的材料参数

Table 1. Material parameters of Tangjiashan landslide

滑坡	γ/(kN·m^-3)	E/GPa	$\nu$	c_p/MPa	c_r/MPa	φ_p/(°)	φ_r/(°)	η
滑体	24.0	2.63	0.2	1.28	0.60	40	28	400
滑床	26.5	12	0.1	2.65	2.65	40	40	0.0

下载: 导出CSV

| 显示表格

图 4 唐家山数值计算模型

Figure 4. Numerical model for Tangjiashan landslide

下载: 全尺寸图片幻灯片

3. 模拟结果分析

3.1 静力计算的模拟结果

在重力作用下，斜坡竖向应力值的范围为0.09~16 MPa，水平应力值范围为0.47~4.3 MPa，说明重力作用下斜坡处于均匀受压状态，同时剪应力仅在坡脚出现应力集中，最大剪应力约为2 MPa。斜坡竖向应力和水平应力等值线平滑，从底部到顶部均匀过渡，并无应力集中或突变的现象，计算结果符合经典力学的理论解。此外，进一步采用强度折减法分析了重力作用下边坡的稳定性，经计算重力作用下坡体的安全系数为1.61，边坡处于整体稳定状态。

3.2 地震诱发滑坡的模拟结果

图 5~7给出了唐家山滑坡的数值模拟结果。图 5给出了不同时刻的唐家山滑坡的形态及水平位移云图，图 6给出了各测点的速度模拟结果，图 7给出了滑动面上3个观测点的摩擦系数的变化过程。模拟结果显示，唐家山滑坡运动总的持续时间约为28~30 s；滑体的最大水平位移约332 m，最大竖向位移约为266 m，最大总位移约370 m；各测点中最大水平速度约为24.8 m/s，最大竖向速度约为18.7 m/s，最大合速度约为27.2 m/s；滑动面上各观测点出现的最小摩擦系数均大致约为0.25。

图 5 模拟的不同时刻的滑坡形态和水平位移云图

Figure 5. Simulated landslide morphology and horizontal displacement at different time

下载: 全尺寸图片幻灯片

图 6 测点的速度的模拟结果

Figure 6. Simulated results of velocity at measuring points

下载: 全尺寸图片幻灯片

图 7 滑动面摩擦系数的变化曲线

Figure 7. Curves of frictional coefficient of sliding surface

下载: 全尺寸图片幻灯片

进一步地，根据图 7中模拟的唐家山滑坡的滑动面摩擦强度随着速度的变化过程，将滑动过程分为4个阶段（图 7（a））：①启动阶段；②摩擦衰减阶段；③低摩擦滑移阶段；④逐步稳定阶段。分述如下：①启动阶段：由于地震的发生，滑动面上的摩擦系数出现较大波动和逐渐降低的趋势，滑体开始出现滑动，该阶段时间较短，滑体滑动速度较低。该阶段，滑动面上的摩擦系数由初始值降低至tanθ（θ为滑坡面的平均倾角），此时滑坡的瞬时安全系数为1，标志着滑坡的启动。后面将进一步结合数值模拟的结果分析地震作用和摩擦强度衰减对启动的影响机理。②摩擦衰减阶段。滑坡启动后随着滑体滑动速度的增加，滑动面的摩擦系数发生大幅降低。根据已有基于速度进行滑坡分类的研究成果^[32-33]，将滑坡由启动至滑移速度达到10 m/s的过程定义为摩擦衰减阶段，该阶段中滑移速度逐步增加到极高的水平，滑动面摩擦强度大幅衰减。如图 6所示，在地震发生9 s左右，滑动面上各测点合速度陆续达到10 m/s。③低摩擦滑移阶段中，滑动面上的摩擦系数始终保持在较低的值，滑体保持在超过10 m/s的高速滑移的状态，因此μ < μ|_{v=10 m/s}。地震后10 s左右，滑坡前端的水平速度达到14.6 m/s，水平位移达到了40 m左右；地震后15 s时滑体前缘已运动至对岸，水平位移达140 m左右；25 s时，滑体的三分之一已爬升到对岸山坡上，前缘水平位移达317 m。地震发生14~20 s之间，各测点的速度先后达到最大值，水平最大速度出现在前缘测点^#1和^#2，在15 s左右达到24.5~26.8 m/s；竖向最大速度出现在尾部测点^#5和^#6，在19 s时均达到18.7 m/s左右；最大合速度出现在测点^#2和^#4，达到27.2 m/s左右。④逐步稳定阶段中，地震发生25 s时滑坡由于对岸山坡的阻碍，滑动速度迅速降低至10 m/s以下，滑动面上摩擦系数迅速增大，并促使速度进一步降低，摩擦系数逐步增加至初始值。此后，滑体滑移15 m后，在28~30 s左右迅速停止并稳定，前端最大总位移约370 m，部分滑体堆积至对岸山坡坡面。

对比模拟的停止后的滑体边界与滑坡的实际轮廓线（如图 5），可见计算结果与滑坡实际轮廓线十分吻合^[29]，且模拟得到的滑坡滑移过程的速度及最终堆积形态与已有的离散元模拟结果^{[21, 34]}以及DDA数值模拟结果^[35]均较为吻合。本文的静动力数值模拟结果与现场调研结果^[29-31]也较为一致。此外，在加速阶段滑动面上3个监测点处的摩擦系数从峰值降到较低水平需约30~40 m左右的水平滑动距离，时间约在地震发生9 s后，而从低摩擦恢复到较高水平只需约15 m的水平滑动距离，需要的时间约为2 s左右，中间较长距离均保持低摩擦的运动状态。这种摩擦强度先衰减再恢复，且衰减需要的距离明显大于恢复需要的距离的现象，与Sone等^[36]的变速试验结果相吻合。唐家山滑坡的模拟结果初步证实了本文所采用的SPH-FEM数值模拟方法的有效性以及所采用的滑动面摩擦强度速率衰减模型的合理性。

边界条件的处理一直以来都是SPH应用的难点，目前的方法主要有虚粒子边界条件^[37]、排斥粒子边界条件^[38]、速度无滑移边界条件^[39]、以及动力粒子边界条件^[40]等。李云屹等^[25]通过在计算域外设置边界粒子，并对边界粒子施加地震位移时程，实现了瑞利波作用下缓倾场地流滑大变形的SPH数值模拟。该方法物理概念清晰，对于给定位移边界条件的情况计算效果较好，但是存在时间步长和支撑域半径的取值对计算精度影响显著的问题^[25]。本文通过将FEM与SPH两种方法耦合，利用SPH实现滑体的大变形发展过程模拟，利用FEM可以方便地模拟水平和竖向地震加速度作用的边界条件，两者耦合的方法具有较好的适用性。

3.3 滑动面摩擦强度衰减对滑坡启动的影响

图 8给出了地震发生后滑体上受力的示意图。考虑地震的影响后，滑动面对滑体产生的动摩擦力f_d和作用在滑体上的动下滑力T_d分别为

${f_{\text{d}}} = \mu (v)(mg\cos \theta + m{a_z}\cos \theta - m{a_x}\sin \theta ) \text{，}$

(5)

${T_{\text{d}}} = mg\sin \theta + m{a_z}\sin \theta + m{a_x}\cos \theta 。$

(6)

图 8 滑体受力示意图

Figure 8. Schematic diagram of force on sliding body

下载: 全尺寸图片幻灯片

式中：m为滑体的质量；g为重力加速度为9.80 N/kg； ${a_z}$ 和 ${a_x}$ 分别为地震产生的竖向和水平加速度；θ为滑动面的平均倾角。定义动摩擦力f_d和动下滑力T_d的比值由动态参数R定义，即

$R = \frac{{{f_{\text{d}}}}}{{{T_{\text{d}}}}} 。$

(7)

当R值初次小于1时，出现了整个滑体上作用的下滑力大于摩擦力的情况，笔者认为此时整个滑动面将处于贯通状态，可以作为判断滑坡开始启动的标志。

图 9（a）给出了地震发生前5 s内滑体6个测点平均速度和平均摩擦系数间的关系。地震发生初期，尽管滑体的剪切速率较低，但是滑动带的平均摩擦强度在波动的同时仍出现逐渐降低的趋势。这与Yang等^[3]的高速旋转剪切摩擦试验规律一致，Yang等发现在恒定的剪切速率下，速率由0.0023 m/s增加至0.0124，0.0436 m/s时，摩擦强度出现了明显的衰减。较低剪切速率下摩擦强度的衰减效应将对滑坡的启动产生显著影响。图 9（b）给出了基于该摩擦系数计算的参数R值的变化曲线，由图 9（b）可见在地震第2秒末时首次出现了R < 1的情况，可以判别滑坡启动。如果不考虑摩擦强度的速率衰减效应，即滑动面摩擦系数μ始终为μ₀=0.84，计算表明在地震发生10 s末才可能出现R < 1的情况，这说明了强度的速率衰减效应对大型滑坡启动的影响。

图 9 滑动带摩擦系数变化和滑坡失稳启动

Figure 9. Variation of friction coefficient of sliding surface and failure of landslide

下载: 全尺寸图片幻灯片

由于临界位移判别方法在大型滑坡的失稳启动中存在局限性，Yang等^[20]认为采用tanθ与摩擦系数间关系判别大型滑坡的启动更为合理，即当滑动面摩擦系数小于tanθ时滑坡启动。本文的滑坡的滑动面倾角由上自下逐渐降低，滑动面的角度略大于滑坡的坡面角。为了简化这里将本滑坡分为上中下3段，假定滑动面倾斜角约等于滑坡坡面倾角，上中下坡面倾角分别约为39.01°，34.65°，33.94°。据此分析滑动面上摩擦系数与滑动面倾角间关系。图 9（a）对比了滑动面平均摩擦系数与上中下段的tanθ与摩擦系数间的关系。由图可见滑动面上段在刚发生地震时，摩擦系数就出现首次小于tanθ的情况，说明滑坡上段滑动面后缘裂隙开始产生。上段滑动面最先发生破坏的模拟结果与已有的数值模拟结果是一致的^[21]。在地震发生的1 s末，滑动面中段出现摩擦系数首次小于tanθ，而在地震发生的2 s末滑动面下段产生裂隙。地震后滑动面上的裂隙逐渐贯通，在2 s末处于整体贯通状态，可以判别滑坡失稳启动。由此得到的判别结果与本文所提出的基于R的判别结果是一致的。

进一步地分析模拟的临界位移的大小。地震2 s末时，各测点的临界位移在1~2 m，远大于Jibson^[18]建议的5~10 cm的临界位移值，模拟结果与Dong等^[19]和Yang等^[20]发现的现象一致。如假设唐家山滑坡模拟中地震荷载在1 s末结束，由图 9（a）可知此刻滑动面下段的摩擦系数仍大于tanθ。所以尽管此时各测点的位移为0.45~1.0 m，可能仍无法发生整体滑移破坏。这是Dong等^[19]和Yang等^[20]发现地震后滑坡位移达到1 m仍未发生破坏的一种可能的解释。

4. 结论

本文基于SPH和FEM耦合的数值计算方法，引入了滑动面摩擦强度的速率衰减模型，发展了一种能够描述滑坡由初始稳定状态、启动状态到滑移破坏状态全过程的模拟方法。

（1）开展了唐家山滑坡的数值模拟，模拟结果与现场勘查结果及室内试验规律较为一致，初步说明了该数值方法能合理考虑地震等复杂边界条件的影响，能合理反映滑动面摩擦强度的速率衰减效应对滑体启动和滑移过程的影响，再现大型滑坡的启动和滑移过程。

（2）基于模拟的剪切面上摩擦系数的演化过程，将滑坡过程分为4个阶段：启动阶段、摩擦衰减阶段、低摩擦滑移阶段和逐步稳定阶段。滑坡在启动和运动过程中，随着滑移速度的增加，滑动带上的摩擦系数逐渐降低，速度增加和摩擦强度衰减的相互促进，触发了滑体的高速滑移。

（3）提出采用滑体上作用的动摩擦力f_d和动下滑力T_d的比值R作为判别指标用于判断大型滑坡的启动，认为首次出现R < 1时滑动面整体贯通并出现失稳启动。采用R的判断结果与采用滑动面摩擦系数（考虑了速率衰减效应）的判别结果一致，揭示了滑坡启动中滑动面摩擦强度衰减和渐进贯通过程，揭示了地震作用与滑动面摩擦参数速率衰减效应共同作用触发了大型滑坡的失稳启动的内在机理。

图 1 含冰透镜体的冻结土体示意图

Figure 1. Schematic graph of frozen soils with ice lens

下载: 全尺寸图片幻灯片

图 2 饱和冻土三相图

Figure 2. Three-phase graph of saturated frozen soils

下载: 全尺寸图片幻灯片

图 3 未冻水膜示意图

Figure 3. Schematic graph of unfrozen water film

下载: 全尺寸图片幻灯片

图 4 水热力耦合流程图

Figure 4. Flow chart of THM model

下载: 全尺寸图片幻灯片

图 5 冻结试验中的温度分布图

Figure 5. Distribution of temperature in freezing tests

下载: 全尺寸图片幻灯片

图 6 冻结试验中的含水率分布图

Figure 6. Distribution of water content in freezing tests

下载: 全尺寸图片幻灯片

图 7 不同时刻冰透镜体分布图（白色代表冰透镜体）

Figure 7. Distribution of ice lens at different time (white color indicates ice lens)

下载: 全尺寸图片幻灯片

图 8 72 h冰透镜体分布对比图

Figure 8. Distribution of ice lens at 72 h

下载: 全尺寸图片幻灯片

图 9 不同时刻未冻水膜压力随位置变化图(试验4)

Figure 9. Charge of water film pressure with position at different time (^#4)

下载: 全尺寸图片幻灯片

图 10 试验1模拟结果对比图

Figure 10. Comparison of simulated results in Test 1

下载: 全尺寸图片幻灯片

图 11 试验4模拟结果对比图

Figure 11. Comparison of simulated results in Test 4

下载: 全尺寸图片幻灯片

图 12 不同模拟工况下含水率结果图

Figure 12. Results of water content under different simulation conditions

下载: 全尺寸图片幻灯片

图 13 不同模拟工况下冻胀量与冻深结果图

Figure 13. Results of frost heave and depth under different simulation conditions

下载: 全尺寸图片幻灯片

表 1 冻胀模型主要变量说明

Table 1 Description of main variables in THM model

变量	物理意义	公式	来源
孔隙冰饱和度S_i	孔隙冰与孔隙的体积之比	$T + {T_0}:1 - {[1 - (T - {T_0})]^\alpha }$	Tice等^[19]
平均水压力P	驱动水分迁移的作用力	$[{\sigma _{\text{t}}} - \sigma + (1 - \chi )L\ln (T/{T_0}){\rho _i}]/[\chi + (1 - \chi ){\rho _{\text{i}}}/{\rho _{\text{w}}}]$	Zhou等^[13]
未冻水膜压力P_Lh	未冻水膜上的相对压力	$P + g(h){\rho _{\text{w}}} + {P_{{\text{atm}}}} - {P_{{\text{or}}}}$	Gilpin^[3]
已冻区渗透系数k	单位水力梯度下的水流量	$T + {T_0}:{\zeta _{\text{w}}}{(1 - {S_{\text{i}}})^{2{\text{b}} + 3}}$	Teng等^[25]
土体体积热容C	单位体积土体温度改变1 ℃需要的热量	$C_{\mathrm{s}} \theta_{\mathrm{s}}+C_{\mathrm{w}} \theta_{\mathrm{w}}+C_{\mathrm{i}} \theta_{\mathrm{i}}$	徐学祖等^[26]
导热系数λ	单位时间通过单位面积传递的热量	$\lambda _{\text{s}}^{{\theta _{\text{s}}}} + \lambda _{\text{w}}^{{\theta _{\text{w}}}} + \lambda _{\text{i}}^{{\theta _{\text{i}}}}$	徐学祖等^[26]

下载: 导出CSV

表 2 冻胀模型输入参数

Table 2 Input parameters in THM model

参数	取值	参数	取值	参数	取值
试验1：k₀ /(m·s^-1)	3.0×10^-10	试验3：T₀/℃	-0.3	${\lambda _{\text{w}}}$ /(W·(m·K)^-1)	0.58
试验2：k₀ /(m·s^-1)	2.5×10^-10	试验4：T₀/℃	-0.3	C_s/(kJ·(m³·K) ^-1)	2160
试验3：k₀ /(m·s^-1)	2.5×10^-10	试验1： ${\rho _{\text{s}}}$ /(kg·m^-3)	2720	C_i/(kJ·(m³·K) ^-1)	1874
试验4：k₀ /(m·s^-1)	2.5×10^-10	试验2： ${\rho _{\text{s}}}$ /(kg·m^-3)	2700	C_w/(kJ·(m³·K) ^-1)	4180
试验1：E_s/MPa	2.5	试验3： ${\rho _{\text{s}}}$ /(kg·m^-3)	2700	L/(kJ·kg^-1)	334.56
试验2：E_s/MPa	1.2	试验4： ${\rho _{\text{s}}}$ /(kg·m^-3)	2768	${\xi _{\text{w}}}$	1
试验3：E_s/MPa	0.8	${\rho _{\text{i}}}$ /(kg·m^-3)	917	b	3
试验4：E_s/MPa	12	${\rho _{\text{w}}}$ /(kg·m^-3)	1000	α	-5
试验1：T₀/℃	0	${\lambda _{\text{s}}}$ /(W·(m·K)^-1)	1.2	β	-8
试验2：T₀/℃	-0.5	${\lambda _{\text{i}}}$ /(W·(m·K) ^-1)	2.22	$\sigma _{{\text{tensile}}}^ *$ /kPa	12

下载: 导出CSV

表 3 不同试验的土性与边界条件参数

Table 3 Soil properties and boundary conditions of different tests

试验编号	顶部温度/℃	底部温度/℃	初始温度/℃	底部状态	初始含水率/%	土的种类	干密度/(g·cm^-3)	上覆压力/ kPa	土体高度/cm	冻结时间/h	来源
^#1	-3℃	3℃	3℃	补水	17	青海-西藏黏土	1.70	100	10	72	文献[16]
^#2	-1.7℃	1℃	1℃	封闭	25.71	内蒙古粉质黏土	1.63	0	15	120	文献[26]
^#3	-2.3℃	0.7℃	0.7℃	补水	18.9	黏土	1.70	0	10	72	文献[28]
^#4	1→-1 20 ℃/h -1→-8 0.1 ℃/h	1℃	1℃	补水	21.7	内蒙古粉质黏土	1.70	500	5	72	文献[24]

下载: 导出CSV

表 4 不同模拟工况参数表

Table 4 Parameters of different simulation conditions

工况	顶温/℃	上覆压力/ kPa	渗透系数/ (10^-10 m·s^-1)	压缩模量/ MPa
^#1	-3	100	3	2.5
^#5	-1.5	100	3	2.5
^#6	-4.5	100	3	2.5
^#7	-3	50	3	2.5
^#8	-3	200	3	2.5
^#9	-3	100	30	2.5
^#10	-3	100	1	2.5
^#11	-3	100	3	0.5
^#12	-3	100	3	5.0

下载: 导出CSV

参考文献(28)

[1]	TABER S. Frost heaving[J]. The Journal of Geology, 1929, 37(5): 428-461. doi: 10.1086/623637
[2]	MILLER R D. Freezing and heaving of saturated and unsaturated soils[J]. Highway Research Record, 1972, 393: 1-11
[3]	GILPIN R R. A model for the prediction of ice lensing and frost heave in soils[J]. Water Resources Research, 1980, 16(5): 918-930. doi: 10.1029/WR016i005p00918
[4]	FOWLER A C. Secondary frost heave in freezing soils[J]. SIAM Journal on Applied Mathematics, 1989, 49(4): 991-1008. doi: 10.1137/0149060
[5]	REMPEL A W, WETTLAUFER J S, WORSTER M G. Premelting dynamics in a continuum model of frost heave[J]. Journal of Fluid Mechanics, 2004, 498: 227-244. doi: 10.1017/S0022112003006761
[6]	KONRAD J M, MORGENSTERN N R. A mechanistic theory of ice lens formation in fine-grained soils[J]. Canadian Geotechnical Journal, 1980, 17(4): 473-486. doi: 10.1139/t80-056
[7]	HARLAN R L. Analysis of coupled heat-fluid transport in partially frozen soil[J]. Water Resources Research, 1973, 9(5): 1314-1323. doi: 10.1029/WR009i005p01314
[8]	TAYLOR G S, LUTHIN J N. A model for coupled heat and moisture transfer during soil freezing[J]. Canadian Geotechnical Journal, 1978, 15(4): 548-555. doi: 10.1139/t78-058
[9]	原国红. 季节冻土水分迁移的机理及数值模拟[D]. 长春: 吉林大学, 2006. YUAN Guohong. The Mechanism and Numerical Simulation of Water Transfer in Seasonal Freezing Soil[D]. Changchun: Jilin University, 2006. (in Chinese)
[10]	李智明. 冻土水热力场耦合机理研究与工程应用[D]. 哈尔滨: 哈尔滨工业大学, 2017. LI Zhiming. Study on Mechanisum of Moisture-Heat-Stress Coupling for Frozen Soil and Engineering Application[D]. Harbin: Harbin Institute of Technology, 2017. (in Chinese)
[11]	白青波, 李旭, 田亚护, 等. 冻土水热耦合方程及数值模拟研究[J]. 岩土工程学报, 2015, 37(增刊2): 131-136. http://manu31.magtech.com.cn/Jwk_ytgcxb/CN/abstract/abstract16270.shtml BAI Qingbo, LI Xu, TIAN Yahu, et al. Equations and numerical simulation for coupled water and heat transfer in frozen soil[J]. Chinese Journal of Geotechnical Engineering, 2015, 37(S2): 131-136. (in Chinese) http://manu31.magtech.com.cn/Jwk_ytgcxb/CN/abstract/abstract16270.shtml
[12]	THOMAS H R, CLEALL P, LI Y C, et al. Modelling of cryogenic processes in permafrost and seasonally frozen soils[J]. Géotechnique, 2009, 59(3): 173-184. doi: 10.1680/geot.2009.59.3.173
[13]	ZHOU J Z, LI D. Numerical analysis of coupled water, heat and stress in saturated freezing soil[J]. Cold Regions Science & Technology, 2012, 72: 43-49.
[14]	HUANG X, RUDOLPH D L. Coupled model for water, vapour, heat, stress and strain fields in variably saturated freezing soils[J]. Advances in Water Resources, 2021, 154: 103945. doi: 10.1016/j.advwatres.2021.103945
[15]	LAI Y M, PEI W S, ZHANG M Y, et al. Study on theory model of hydro-thermal–mechanical interaction process in saturated freezing silty soil[J]. International Journal of Heat & Mass Transfer, 2014, 78: 805-819.
[16]	MING F, LI D Q. Experimental and theoretical investigations on frost heave in porous media[J]. Mathematical Problems in Engineering, 2015: 1-9.
[17]	曾桂军, 张明义, 李振萍, 等. 正冻土中冰透镜体形成力学判据的分析讨论[J]. 冰川冻土, 2015, 37(1): 192-201. ZENG Guijun, ZHANG Mingyi, LI Zhenping, et al. Review of mechanical criterion for formation office lens in freezing soil[J]. Journal of Glaciology and Geocryology, 2015, 37(1): 192-201. (in Chinese)
[18]	NIXON, DERICK J. Discrete ice lens theory for frost heave in soils[J]. Canadian Geotechnical Journal, 1991, 28(6): 843-859. doi: 10.1139/t91-102
[19]	TICE A R, ANDERSON D M, BANIN A. The prediction of unfrozen water contents in frozen soils from liquid limit determinations[C]// Symposium on Frost Action on Roads. Oslo, 1973.
[20]	季雨坤. 冰透镜体生长机制及水热力耦合冻胀特性研究[D]. 徐州: 中国矿业大学, 2019. JI Yu-kun. Ice Lens Growth Mechanism and Hydro-Thermal-Mechanical Coupling Research on Frost Heave[D]. Xuzhou: China University of Mining and Technology, 2019. (in Chinese)
[21]	ZHOU Y, ZHOU G Q. Intermittent freezing mode to reduce frost heave in freezing soils-experiments and mechanism analysis[J]. Canadian Geotechnical Journal, 2012, 49(6): 686-693. doi: 10.1139/t2012-028
[22]	SHENG D C, ZHANG S, YU Z, et al. Assessing frost susceptibility of soils using PCHeave[J]. Cold regions science and technology, 2013, 95(11): 27-38.
[23]	TENG J D, LIU J L, ZHANG S, et al. Frost heave in coarse-grained soils: experimental evidence and numerical modelling[J]. Géotechnique, 2022: 1-12
[24]	曹宏章. 饱和颗粒土冻结过程中的多场耦合研究[D]. 北京: 中国科学院研究生院, 2006. CAO Hongzhang. Research on Fields Coupling in Saturated Granular Soil Freezing Process[D]. Beijing: University of Chinese Academy of Sciences, 2006. (in Chinese)
[25]	TENG J D, YAN H, LIANG S, et al. Generalising the kozeny-carman equation to frozen soils[J]. Journal of Hydrology, 2020, 594: 125885.
[26]	徐学祖, 王家澄, 张立新. 冻土物理学[M]. 北京: 科学出版社, 2001. XU Xuezu, WANG Jiacheng, ZHANG Lixin. Permafrost Physics[M]. Beijing: Science Press, 2001. (in Chinese)
[27]	O'NEILL K, MILLER R D. Exploration of a rigid ice model of frost heave[J]. Water Resources Research, 1985, 21(3): 281-296.
[28]	周家作. 土在冻融过程中水、热、力的相互作用研究[D]. 北京: 中国科学院研究生院, 2012. ZHOU Jiazuo. Study on the Interaction of Water, Heat and Force in Soil during Freezing and Thawing[D]. Beijing: Graduate University of Chinese Academy of Sciences, 2012. (in Chinese)