何志亮

（中铁水利水电规划设计集团有限公司，江西南昌，330029）

由于暴雨或地震影响，山区通常会发生临时或永久性河流阻塞。由这一过程形成的滑坡堆积坝可能会因漫顶、坝体突然坍塌或渐进性破坏引起的侵蚀而失效。大坝的滑坡是一个非稳定渗流过程，边坡稳定分析与非稳定渗流分析相结合有助于滑坡堆积坝的稳定性研究。

许多学者采用不同的边坡稳定性分析方法，重点是结合水文模型与渗流稳定分析。肖耀廷等人[1]采用了VIC模型，对滑坡体稳定性进行了有限元评估与分析；鲜丹等人[2]以甘肃某村的滑坡体为研究对象，研究了初始含水率对边坡稳定性的影响原理；辛长虹等人[3]基于非饱和渗流分析理论，结合有限元数值方法，分析了大坝位移和应力规律；付建新等[4]基于交替隐式有限差分法，将不同含水率作为初始条件进行二维入渗数值计算与分析；陈玺等人[5]对溃坝时的坝坡稳定进行分析；周建云[6]分析了降雨入渗对水库坝坡失稳的影响，结合2D模型对水库溃坝进行演算；孙蔚等人[7]考虑了渗流和岩土体的耦合作用，进行了滑坡堆积坝的稳定性分析。但这些研究大多局限于二维分析，仅适用于三维长边坡。然而，天然斜坡和滑坡堆积坝在狭窄的U形或V形山谷中的破坏是三维的，因此三维方法更适合于分析此类稳定性问题。

1 数值模型

采用极限平衡法对边坡稳定性进行了评价。根据坝体内部水的入渗情况，计算安全系数，寻找安全系数最低的临界滑动面。

由两个模型组成，瞬态渗流模型计算坝内孔隙水压力和含水量的变化。图1为边坡瞬态稳定分析耦合模型通用流程。由于上游水库水位逐渐升高，边坡稳定模型根据坝体孔隙水压力和水分运动的变化计算边坡的安全系数和临界滑动面的几何形状。上游水库水位（WL）在每个时间步由水量平衡方程确定，考虑流入流量（Qin）、坝体内入渗率（Qsh）和水槽几何形状（宽度和坡度），耦合模型的总体轮廓如图1所示。

图1 边坡暂态稳定分析耦合模型通用流程Fig.1 General flow chart of coupling model for transient stabili-ty analysis of slope

1.1 渗流模型

在狭窄的山谷中形成的滑坡堆积坝渗流可能产生三维效应。为考虑非饱和滑坡堆积坝三维渗流问题，建立三维渗流模型。为了评价饱和土的孔隙水压力变化，采用了基于压力的修正方程：

式中：h为水压头；K x(h)、K y(h)、K z(h)为x、y和z方向的导水率；C(h)=∂θ/∂h，为比湿容量，θ为土壤体积含水量；t为时间；S w为饱和比；S s为比库容；x、y为水平空间坐标；z为垂直空间坐标，取向上为正。公式（1）表示非饱和域和饱和域中的流动。对于饱和域，K x(h)=K y(h)=K z(h)=K s，θ=θs，C(h)=0，S w=1，其中K s和θs分别为饱和导水率和饱和含水量。S s取决于固体基质和流体的可压缩性，因此在非饱和无侧压多孔介质中，S s接近于0。

为求解式（1）中的方程，结合式（2）和式（3）：

式中：α和n v为方程参数；S e为有效饱和度；θs和θr分别为饱和含水量和残余含水量。

1.2 边坡稳定性模型

用极限平衡分析法评价滑坡堆积坝的瞬态边坡稳定性，包括计算安全系数和寻找安全系数最低的临界滑动面。基于对任意滑移柱的极限平衡分析，本研究中考虑的边坡是平缓的（小于45°）。一般情况下，安全系数只能由一个公式而不是联立方程确定，因此迭代程序被简化。边坡稳定性模型结合了基于动态规划的最小化方法和随机数生成方法，将该模型与渗流模型耦合，进行边坡稳定分析。

山中景色随四时而变。出自北宋著名的画家兼山水画理论家郭熙写的《山水训》：“真山之烟岚，四时不同。春山淡冶而如笑，夏山苍翠而欲滴，秋山明净而如妆，冬山惨淡而如睡。”正写出了山水随季节变化而迥然不同的物象特征，但这特征正如描写的美人的姿态一样“淡冶而如笑，苍翠而欲滴，明净而如妆，惨淡而如睡”，以此来比喻山水的春夏秋冬。

图2显示了斜坡内的三维滑动体和垂直分割柱。图3显示了作用在图2中典型柱上的力：Wij为柱的重量；Pij为作用在柱顶的竖向外力；T ij和Nij为作用在柱脚上的剪力和总法向力；Q ij为作用在柱侧面的所有柱间力的合力。将Q ij分解为平行于x轴的Q1和以一定角度β=tan-1(ηtanαxzij)倾斜的分量Q2，其中η为常数。此外，图3中的αxzij和αyzij为柱基在xz和yz平面上与水平方向的倾角。

图2 三维斜坡的阶段状态系统和划分方案Fig.2 System of stage state and classification scheme of 3D slope

图3 作用在典型柱上的力Fig.3 The force acting on the typical column

在整个滑动体的水平和垂直方向上应用平衡条件，引入两个安全系数Fsh和Fsv，这两个系数都取决于η。对于较小的η值，通常满足方程F=Fsh=Fsv，计算出的Fsh值受η变化的影响不大。可以由η=0得出Fsh并基于此近似值计算3D安全系数，安全系数F可通过式（4）表示：

式中：c ij和φ为土的有效强度参数；u ij为柱基孔 隙 水 压 力 ；和为作用于柱上下游表面的水压力合力。

莫尔-库仑方程可用于计算非饱和土的抗剪强度，从而得出公式（4）。然而，在这项研究中，由于负孔隙水压力，剪切强度的增加被忽略了。

在公式（4）中，W ij、P ij和uij（无黏性土的cij=0）是时间相关的，可以在每个时间步长通过与3D瞬态渗流模型耦合确定。渗流模型的输出（每个单元的孔隙水压力和含水量）和边界条件（地表水的深度）可用于确定这些物理量。

2 试验研究

测量水分分布和观察破坏面的试验总结如表1所示。使用平均直径为1 mm的混合硅砂在水槽中制成三角坝。方程（2）和（3）中的参数（包括θr）通过对pF仪试验获得的土壤持水量数据进行非线性回归分析进行估算。

表1 试验总结Table 1 Test results

利用含水率反射计测量了渗流过程中含水率随时间的变化。为了测量坝坡在滑动过程中的运动，在槽壁面坝体内放置红色的泥沙条带。在渡槽侧边放置数码摄像机，拍摄渡槽突然滑动时的滑动面形状。通过对水槽两侧的视频分析，得出了坝体滑动过程中的滑动面形状。

使用长500 cm、宽30 cm、高50 cm的水槽。设定水槽的坡度为20°。坝体形状及尺寸见图4和图5。含水率反射计的排列如图5所示。在水槽的B侧开孔，并从该侧将12个仪器插入坝体内。进行了两组试验来测量坝体不同位置的水分剖面。

图4 坝体形状及坝顶断面Fig.4 Dam shape and dam crest section

图5 含水率反射计（1-12）布置图和B侧视图Fig.5 Layout diagram of moisture content reflectometry(1-12)and view of B side

3 结果分析

库水位逐渐上升导致水渗入坝体，增加了孔隙水压力和坝体重量。当坝体自重增大引起的剪应力大于孔隙水压力增大引起的抗剪应力时，坝体就会发生滑动。

稳定的水流量来自水槽上游，测量出坝体中水分含量。测得的水分含量是距离水槽侧面18 cm处的平均值（见图5）。由于水在坝体内运动的三维效应，水分运动不是恒定的。尽管两个试验中使用的排放量几乎相等（见图6），但3D-1和3D-2测量的水分分布略有不同，这可能由两个试验中饱和导水率的变化所致。图7显示了不同水分剖面的模拟和试验结果比较，具有良好的一致性。

图6 水分剖面比较（示例3D-1和3D-2）Fig.6 Comparison of water profiles(experiment 3D-1 and 3D-2)

图7 水分分布比较（模拟和试验）Fig.7 Comparison of water distribution(simulation and experi-ment)

对比两个试验3D-3和3D-4以观察边坡的破坏情况。试验3D-3（图8）中，在930 s处观察到斜坡失效。而在试验3D-4的情况下，在1 030 s处观察到边坡的破坏，滑动面也比3D-3更深（图9），这可能由两个试验在压实和饱和导水率方面的不均匀性不同而导致。

图8 930 s时的滑动面（试验：3D-3）Fig.8 Sliding surface at 930 seconds(experiment:3D-3)

图9 1 030 s时的滑动面（试验：3D-4）Fig.9 Sliding surface at 1 030 seconds(experiment:3D-4)

在三维边坡稳定性初步分析的基础上，在每个阶段平面上生成了上千个状态。渗流模型中用于模拟的其他水力条件/参数和网格系统为Qin=29.8 cm3/s，K s=0.000 3 m/s，θs=0.287，θr=0.045，α=5.5 m-1，n v=3.2，Δt=0.01 s，块体尺寸为10 mm。边坡稳定模型采用的单元大小Δx=5 cm、Δy=3 cm。

在t=0 s和t=770 s时，模拟的临界滑移面如图10所示。模拟安全系数（F）一开始时大于1。770 s时，模拟的安全系数小于1，但试验中观察到的破坏时间约为930 s。

模拟过程中，考虑蒸发作用是以减少29 cm3/s的流量实现的，考虑水力传导率的不确定性是通过减少0.000 28 m/s的饱和传导率来进行的。这两种情况下，大坝都在790 s时失效。Ks=0.000 25 m/s时，模拟失效时间为830 s。因此，失效时间也取决于饱和导水率。初步分析（图11）表明，使用此3D方法计算的安全系数，当滑坡堆积坝破坏时为1.05~1.06。因此，如果使用更严格的边坡稳定性分析方法，破坏时间将延长并接近观察到的时间。此外，在计算中也忽略了水槽侧壁的摩擦力。水槽两面破坏面（图10）与试验吻合较好。

图10 模拟临界滑动面（3D）Fig.10 Simulation of criticalsliding surface(3D)

图11 用3D方法模拟不同Ks和安全系数下的破坏时间Fig.11 Simulation of failure time by 3D method with different Ks and safety factors

滑动面的形状取决于许多因素，如土体物理力学性质、植被分布、孔隙水压力等，可进一步完善模型，采用更严格的边坡稳定分析方法，是三维边坡稳定分析实际应用的关键。

4 结语

通过水槽试验和数值模拟，对滑坡堆积坝的三维瞬态边坡稳定性进行了分析，得出以下结论：

（1）在上游水库稳定泄流试验水槽中，对滑坡堆积坝的突然破坏进行了研究。库水位逐渐上升使水渗入坝体，增大了动剪应力，当动剪应力大于抗剪应力时，坝体突然垮塌。

（2）采用极限平衡法建立了边坡稳定与非稳定渗流耦合模型。通过数值模拟和水槽试验，研究了滑坡溃坝机理。结果表明，数值模拟和试验结果在坝体水分运动和预测临界滑动面方面非常接近。

（3）然而，在模拟中，坝体发生破坏的时间较早。今后应通过引入更严格的边坡稳定性分析方法对模型进行进一步的改进，使其既可用于滑坡堆积坝的边坡稳定性分析，也可用于天然边坡的稳定性分析。