滑坡运动二维黏壶分析模型

雷先顺，李 漪

(1.湖北省电力勘测设计院有限公司，湖北 武汉 430040；2.合肥工业大学，安徽 合肥 230009)

0 引言

在我国，每年发生的各种地质灾害中，滑坡灾害发生的频率最高，产生的危害最大[1]。为了对滑坡灾害进行有效的预测和防治，需要对其发生机理、运动过程和灾害评估等问题进行深入研究。

关于滑坡的静力学分析和灾害评估，已经有大量的研究成果[2-3]。而滑坡运动过程分析是滑坡静力学分析的后续步骤，并且关系到最终的滑坡灾害评估结果，具有重要的理论研究意义和实际用途。在滑坡运动过程理论研究方面，主要有质点能量法、流体理论、条分法和离散元法等理论，但是各有其局限性。质点能量法[4]将整个滑体当做一个整体，运用能量转化关系求得最终结果，但是没有考虑滑道的性质与滑体各位置的变形和速度不同。流体理论[5]运用流体力学的观点，将滑体看作流动性的连续介质，但是不能反映滑坡过程中的破碎、相互作用等现象，且不适用于一些整体性较好的滑坡。条分法[6-9]是一种运用较广泛的方法，它把滑体分成若干个竖向条块，可以和滑坡静力学分析相结合，但是没有考虑滑体在运动过程的变形和内部耗能。离散元理论[10]基于微观颗粒间的相互作用原理，用一个个颗粒组成的堆积体来模拟实际的滑体，可以十分方便地分析滑体内部颗粒的力学行为，是一种较准确的方法。

此外，基于灾变链式理论，肖盛燮[11]等提出一种考虑孔隙水和空气阻尼耦合作用的滑坡力学模型。江巍[12]等对DDA分析方法进行了改进，可以很好地模拟滑坡运动过程。

滑体在下滑过程中，是其重力势能转化成滑体动能和热能，并最终全部转化成热能、声能等其他能量的过程。其中热能包括：滑体和滑道摩擦产生的热能和滑体自身变形、内部相互摩擦碰撞产生的热能。特别是对于高速岩石滑坡，在其运动过程中，内部岩块往往会发生崩解，相互碰撞作用[13]十分频繁。所以，在对滑体运动过程进行分析时，有必要考虑这一部分能量损耗。本文基于条分法的特点，考虑滑体的内部耗能[14]，采用黏壶连接代替弹簧连接，再根据条块运动方程，得到滑坡最终的各项评估指标。

1 滑坡运动分析模型

1.1 分析模型的提出

在条分法中，将滑坡体离散成垂直条，且假定条块为刚体，或者在滑动过程中滑块仍然保持垂直状态。但事实上，这些垂直块体在运动过程中，会变成曲线状，并且颗粒组成越松散，滑速越大，滑体变形越普遍，甚至会出现碎屑流现象，这与传统条分法的假定显然有很大的不同。

根据滑体在运动过程中能量转化的规律，对于消耗的热能，其中滑体和滑道摩擦消耗的能量可以通过滑块底部摩擦力做功求得。而滑体由于自身变形，内部颗粒间互相作用消耗的能量则无法定量求出，吴越[14]等通过试验发现，在滑坡过程中，滑体内部耗能所占比例较大，是不能忽略的，且内部变形越大，内部耗能所占比例越大。这里，基于滑坡过程中内部耗能的特点，并结合传统条分法的一些优点，提出一种滑坡运动过程黏壶分析模型。

将滑坡运动考虑为平面问题，即不考虑滑体的横向扩展问题，将滑体分成等质量的若干个滑块，滑块间用黏壶连接。相邻前后滑块间速度不同时，会导致黏壶变形，利用黏壶变形做功来模拟滑坡过程中的内部能耗。

另外，在滑体运动轨迹方面，也做出了相应的简化处理，将运动轨迹简化成滑坡启动区、滑坡加速区和滑坡减速堆积区。滑坡运动模型如图1所示。

图1 滑坡运动模型

由滑体剪切试验可知，滑体剪切段的剪应力[15]为：

式中：G为滑体剪切模量；γ为滑体某处剪应变，μ为滑体内某处剪切位移；μm为滑体剪切破坏时的临界位移；τm为滑体的残余剪切强度。

只考虑滑体发生剪切破坏后的过程，即滑体内部的剪切力达到其残余剪切强度。另外，当相邻滑块速度不同时，两者会出现一定的相对剪切滑动，产生剪切力；反之，则没有剪切力。从而，可以得到滑体内崩解剪切面上的剪切力：

式中：τ为崩解面上的剪切力；V1和V2分别为相邻两滑块的速度。

而对崩解剪切面进行分析，由库伦—摩尔准则：

式中：σ为剪切面上的正压力；θ为每个剪切面上的残余内摩擦角。可以由现场勘察资料得到。

大量滑坡实例和室内试验发现滑体运动过程中，其自身发生了很大的变形，必然存在颗粒间的相互摩擦，滑体内部会损失部分能量，这里将这部分内部耗能转化成滑块间黏壶轴向变形做功耗能。

则黏壶的黏滞力T等于崩解面上的总剪力Q，黏壶轴向变形量等于上下相邻滑块的相对运动距离。可以得到：

将式(3)带入式(4)中，可得：

这里，将滑体分成n个相互平行的等质量块体，从下到上依次编号为1至n号，假设滑体在运动过程中相邻块体间均发生了剪切破坏，即滑体内出现n-1个内部崩解剪切面，则有n-1个黏壶。并且，假设滑体中位于上层的滑块先滑出，滑块由上到下依次运动，则崩解剪切面上的正压力σ即为块体的自身重力在垂直作用在剪切面上的分量，然后推导得到黏壶模型中的黏滞力为：

岩土体的残余强度大小与剪切速率快慢关系密切，式(6)中k为与剪切速率有关的修正系数，针对不同的滑坡形式，可以选取不同的值。

从而可以得到黏壶黏滞力的表达式：

1.2 建立运动方程

假设滑块所受到的黏壶黏滞力均是平行于滑面方向，且取沿滑面向下方向为正方向。对于相邻的两个滑块，若V1＞V2则滑块1所受的黏壶黏滞力为抗滑力，滑块2所受的黏壶黏滞力则为下滑力；反之，若V1＜V2，则结果相反。从而，可以建立每个滑块的运动学方程，沿滑面方向，滑块所受的合力为：

则有：

在具体计算时，可以假设所有滑块的初速度均为零，而且滑坡开始运动时，所有黏壶均处于拉伸状态。那么，可以得出初始时，滑块的加速度。然后，再取一个很小的时间间隔Δt，可认为在(t，t+Δt)时段内滑块的加速度不变，即滑块做匀加速直线运动，则变形块在t+Δt时刻的速度为：

在t+Δt时刻的位置为：

当1号滑块的位移达到最大值时，停止计算，即可以得到最终滑坡的冲程和运动时间t总。相应地，也可以得到滑坡过程中的内部耗能和滑体自身的变形量。其中，内部耗能可以通过黏壶轴向力做功求出，第i个黏壶耗能：

则整个滑体的内部耗能为：

滑体自身的变形量即为1号滑块和n号滑块运动距离的差值：

1.3 计算步骤

通过C语言编程，可以非常方便地进行计算。且不需要往复迭代，最终可以得到滑坡的冲程、任意时刻的滑体变形量、最大速度、耗能情况等参数。这里，将第一块滑块速度为0时作为程序终止的判定条件。具体计算步骤如图2所示。

图2 计算步骤流程图

2 计算与分析

根据前文，该滑坡分析模型的具体推导过程，现运用该黏壶分析模型对前期的室内模型试验进行计算分析，并和离散元软件(PFC 2D)的计算结果进行对比。

2.1 模型计算结果

模型试验装置如图3所示。其中，滑坡启动区、滑坡加速区和减速堆积区的坡度分别为30°、30°和0°，三段长度分别为0.4 m、3.6 m和10 m。试验所用颗粒材料是碎石子，其粒径范围为0～20 mm，堆积体重度为13.5 kN/m3。室内试验的结果滑体的最终冲程为0.92 m，堆积长度为1.1 m，总共耗时4.5 s。在其运动过程中，出现了明显的颗粒跳动、相互碰撞等现象，滑体自身发生了很大的变形。如果用常规的条分法，显然是不符合实际的。

图3 室内模型试验图

首先，运用本文提出的分析模型进行计算。因为模型试验中滑体的运动路径均为直线段，在分析模型中，也只考虑3个直线段。在该计算模型中，参数的取值(包括各位置残余内摩擦角θ、修正系数k和滑面的摩擦系数)对最终结果的影响较大。滑体和滑面的参数取值需要经过试验和反分析相结合的方法得到。滑体的体积为40 L，滑体的质量为55.1 kg。这里，需要确定的系数包括滑体各位置残余内摩擦角、修正系数k和滑面摩擦系数。由于滑体的组成比较均匀，所以可以认为滑体各位置处颗粒的残余内摩擦角θ相同，即在确定黏壶黏滞力时，可以将每个崩解面上的tanθ取相同的数值。三段运动路径中，启动区和加速区的滑面材料相同，均为光滑板，而减速区滑面采用粗糙板。三段路径的长度和坡度均和模型试验保持一致，这里三段路径的长度S1=0.4 m，S2=3.6 m，S3=10 m；三段路径坡度分别为α1=α2=30°，α3=0°；三段路径的摩擦系数分别为0.51、0.36和0.62。碎石子的直剪摩擦角为40.1°。由于该模型试验中滑体全部由无黏性颗粒组成，其运动形式类似于碎屑流，内部崩解速率相对较快。综合参考文献[14]中的取值，这里将修正系数K设定为0.15。程序计算中，将计算时步取为0.05 s。

这里先探究滑块不同分块数n的影响，如图4所示。从图4中可以看出，当滑块数n取不同的数值时，对1号滑块最终的运动距离S有一定的影响。当n=2时，S为5.23 m；当n=7时，S为5.33 m；当n=12时，S为5.27 m；当n继续增大时，S最终结果保持在5.22左右。即当n取不同值时，对最终结果的影响在10%之内。实际计算时，可以将n取一系列的数值，然后再找出计算结果的范围。另外，可以看出滑体的运动距离都是随时间的变化呈先增加后逐渐趋于稳定。当运动到2.7 s时，1号滑块的运动距离几乎保持不变，而后面的滑块继续运动，和模型试验的试验现象比较接近。

图4 滑块分块数对1号滑块运动距离的影响规律

因此，以n=7为例，则模型中黏壶的黏滞力为：

T=σ×tanθ×k= 9.7 N

然后，将各参数带入编制好的程序中，经过计算可以得到，7个滑块的最终运动距离依次为：5.33 m、5.05 m、4.90 m、4.67 m、4.42 m、4.14 m和3.84 m。即滑体最终的冲程是1.33 m，相对于模型试验的结果，误差为44.5%。滑体的堆积长度为1.49 m，相对于模型试验的结果，误差为35.4%；滑体的内部变形为1.09 m，相对于初始堆积长度，变形比例达到172%。表1为计算结果和模型试验结果对比情况。另外，分析滑体堆积长度变化情况，如图5所示。

表1 理论计算与模型试验结果对比 m

图5 滑体长度变化历时图

可以看出，开始阶段滑体前端的滑块先运动，而后面滑块运动较慢，导致前后滑块的距离越来越远。当运动到1.8 s左右时，后方滑块速度大于前端滑块，前后滑块的相对距离逐渐减小，直至到达稳定值。这和试验中观察到的现象相同，即滑体先逐渐在斜面上展开，然后再逐渐在减速区上堆积，该模型可以较好地模拟滑体的内部变形特性。

2.2 计算结果对比分析

离散元是一种比较好的分析方法，可以模拟物体发生的变形，并且可以分析其运动过程。PFC2D的数值模型和室内试验的模型大致相同，如图6所示。主要由装料槽、滑动加速面和减速堆积面组成，且尺寸和模型试验尺寸相同，滑道强度摩擦角为31.8°。堆积体颗粒总共包括800个球，半径在0～2 cm范围内，并按照实际砂子的颗粒级配分布。球间摩擦系数为0.15，球间法向和切向刚度均为1e8，球的密度为2.49×103kg/m3。通过计算得到滑体的最终冲程为1.10 m，和试验结果0.92 m的差距较小。

图6 PFC模型

图7为两种计算方法下，模型试验中滑坡体的最终堆积状态。本文计算的结果和试验现象较为吻合，且计算模型方法较简单，易于编程，但是由于没有考虑滑体的形状特征，不能模拟出滑体运动过程中滑体各位置处高度的变化，从而无法得到滑体最终堆积形态，这是需要改进的地方。而运用PFC分析软件可以较好地模拟室内模型试验的结果，最终的堆积形态十分相似，但由于在离散元方法中，各试验参数取值较困难，需要经过试算，方法不易掌握。

图7 计算结果对比图

而后，分析滑体运动过程中各参数值的变化情况。在滑坡评估中，比较重要的是滑坡最终的冲程，图8为两种计算方法下，1号滑块运动距离随时间变化图。可知，2种方法中滑体运动距离随时间的变化规律大致相同，滑体的运动距离都是随时间的变化呈先增加后逐渐趋于稳定。本文计算模型得出的滑坡运动距离相对于PFC 2D计算结果较大，相对误差为20%。另一方面，2种方法计算的滑坡运动时间不同，本文计算的运动时间为2.7 s，而运用PFC 2D计算的运动时间为3.9 s，实际运动时间为4.5 s。总的来说，本文计算模型的计算结果均使最终的风险偏大，且误差相对较小。

图8 1号滑块运动距离历时图

3 结论

本文基于滑坡运动过程中能量转化规律，考虑滑体运动过程中自身崩解、碰撞等内部耗能，提出一种滑坡运动分析模型，并编制了相关计算程序。

1)采用本文分析模型，可以简单、准确地计算滑坡过程中的最大滑速、滑坡冲程、运动时间等滑坡致灾指标。

2)基于室内模型试验，分别利用该分析模型和离散元分析进行计算，并对比计算结果，证明了该分析模型的适用性。

3)本文计算模型可以较好地模拟出现变形的滑坡现象，相比传统刚体运动模型，更加符合实际情况。

4)由于没有考虑滑体的结构组成，滑道定义过于简化等原因，本文提出的滑坡分析模型存在一定的不足。如不能模拟出滑坡过程中的滑体高度变化情况、不能分析滑体起始破坏情况等。