基于点稳定系数法的非饱和河岸边坡稳定性分析

，,,,

(中国地质大学(武汉) 教育部长江三峡库区地质灾害研究中心，武汉 430074)

基于点稳定系数法的非饱和河岸边坡稳定性分析

王艳昆，熊承仁,王菁莪,邹宗兴,韩天骅

(中国地质大学(武汉) 教育部长江三峡库区地质灾害研究中心，武汉 430074)

河岸边坡的整体稳定性计算，无法定量反映水位变化对边坡各部位局部稳定性的影响。基于点稳定系数法与非饱和土力学理论，运用非饱和流固耦合模型进行河岸边坡的稳定性数值分析。以某大型模型试验为例，通过绘制边坡在水位升降作用下各时段所对应的点稳定系数分布云图，揭示边坡的破坏原因与内部稳定性变化过程，并与实际观测结果进行对比。研究表明：点稳定系数法能合理地描述河岸边坡各位置的稳定性及其动态变化，与实际吻合良好；边坡临空面稳定性较低，水位缓慢上升将导致有效应力降低及非饱和区基质吸力丧失，易诱发浅层失稳；水位下降形成的动水压力，是导致边坡稳定性下降的主因，坡脚处最易失稳，致使边坡发生牵引破坏；边坡内部点稳定系数随水力梯度的增大而降低，非稳定区分布面积随水力梯度的增大而增大，并主要分布于浸润线以下、坡脚以上。研究成果为此类边坡的稳定性评价及治理提供了一定的科学依据。

河岸边坡；水位升降；非饱和流固耦合；数值分析;点稳定系数法

1 研究背景

河岸边坡广泛存在于自然界中，水位升降作用下河岸边坡的稳定性是人们关注的焦点。

近年来，国内外学者一般采用非饱和土力学研究河岸边坡稳定问题，朱朋等[1]依据非饱和土渗流理论和极限平衡理论对三峡库区某滑坡在库水位波动及降雨条件下的渗流及稳定性进行了分析；周永强等[2]运用非饱和土渗流理论和抗剪强度理论，以付家坪子高陡滑坡为例，分析了库水位变化与降雨联合作用工况下滑坡的稳定性；沈细中等[3]利用非饱和流固耦合模型，定义材料参数为饱和度的函数，对小浪底库岸1#滑坡在不同工况下的安全性进行了评价；倪卫达等[4-5]采用饱和-非饱和渗流有限元算法获取三峡库区黄土坡临江1#滑坡在水库运营期间地下水渗流场的动态变化过程，并以此为基础研究滑坡的变形和稳定性演化规律；周建烽等[6]运用有限元塑性极限分析下限法求解非饱和非稳定渗流作用下边坡稳定问题；张旭等[7]对库水作用下的三峡库区巫峡段龚家坊堆积层滑坡的稳定性等问题进行饱和-非饱和数值计算。运用非饱和流固耦合模型并考虑岩土体材料的湿化过程的计算结果更加符合实际，但以上关于河岸边坡在水位升降作用下的稳定性研究多是基于整体的方法，通过给出潜在滑面的稳定性系数对边坡进行定量评价，对于边坡各部位的稳定性状况及其变化无法定量展现。

点稳定系数法可用来衡量斜坡体内单个质点稳定程度[8]，它可揭露边坡的初始破坏位置以及潜在不稳定区域随孔隙水压力变化的发展趋势[9]。因此，为研究水位升降作用下非饱和土边坡各部位的稳定程度及演化机制，运用点稳定系数法定量评估边坡内部稳定性对孔隙水压力变化的响应情况。以某大型模型试验为例，通过非饱和流固耦合数值计算,得到边坡各时段的内部应力及孔隙水压力值；利用点稳定系数的定义计算得到边坡各位置的点稳定系数并绘制了分布云图，以数值分析的手段揭示了河岸边坡在水位升降条件下的失稳原因，为此类边坡的稳定性评价及治理提供了一定的科学依据。

2 饱和-非饱和流固耦合理论

2.1 本构方程

二维空间下，由非饱和流固耦合理论[10-11]，非饱和土应力应变的增量关系可写为

{Δσ}=[D]{Δε}-

[D]{mH}Δ(ua-uw)+{Δua} 。

(1)

式中：σ为法向应力;[D]为本构关系矩阵;ε为法向应变;{mH}为与基质吸力(ua-uw)有关的非饱和土结构模量矩阵;ua为孔隙气压力;uw为孔隙水压力。

2.2 水流连续方程

单个土体单元的孔隙水二维流动可由达西等式给出,即

(2)

式中：kx与ky分别为x，y方向的渗透系数;vw为渗流速度;γw为水的重度;θw为体积含水量函数。

非饱和状态下，土体体积含水率函数及渗透系数函数可用VG模型表示：

(3)

(4)

m=1-1/n。

(5)

式中：θr,θs分别为残余含水率与饱和含水率;α为进气值的倒数;h为吸力水头;ψ为吸力压力(kPa);ks为饱和渗透系数;m,n为拟合参数。

3 点稳定系数定义

摩尔-库伦强度准则下点稳定系数可定义为边坡内一点抗剪强度与其当前状态所受剪应力的比值,即

(6)

式中：LFS为点稳定系数;τf为抗剪强度;τ*为剪应力。

当土体处于非饱和状态时，非饱和区域的抗剪强度可用Fredlund[10]双应力变量表示,即

τf=c′+(σ-ua)tanφ′+(ua-uw)tanφb。

(7)

式中：c′为有效黏聚力;φ′为有效内摩擦角;(σ-ua)为净法向应力;tanφb为与基质吸力(ua-uw)有关的参数。

1996年Vanapalli等[12]将tanφb简化修正为

(8)

此时黏聚力为

(9)

式中cb为表观黏聚力。

将式(9)代入式(7)可得

τf=C+(σ-ua)tanφ′。

(10)

图1为边坡内点稳定系数定义示意图，非饱和状态下，非饱和区点稳定系数LFS由图中几何关系可得

(11)

经变量代换得

(12)

简化起见，将孔隙气压力ua设定为0，则

(13)

式中σ1,σ3分别为最大主应力与最小主应力。

图1 点稳定系数定义示意图Fig.1 Illustration of the concept of local factor of safety

当土体处于饱和状态时，同样由几何关系，点稳定系数可用下式表示为

(14)

4 算例分析

贾官伟等[13]通过大型模型试验研究了水位骤降引致河岸边坡滑坡的原因及失稳模式，并指出动水压力是滑坡产生的主要原因，松散填土边坡的失稳模式为有多重滑面的牵引破坏模式。模型水位缓慢上升最高水位5.6 m时，填土边坡坡面逐渐发生了局部崩塌，随着水位的上升,坍塌的范围逐渐扩大，最终边坡坡肩退到原位置后方1 m左右；水位下降时，以1 m/h速度下降至3 m水位，在这期间随着水位的下降，坡顶出现张拉裂缝且深度和宽度越来越大，并分别在第56,100,100 min稍后发生滑块1,2,3的滑动。该试验结果现象显著，试验模型形态简洁经典，因此以贾官伟等研究中的模型试验为研究对象，通过数值模拟的手段分析水位升降对河岸边坡稳定性的影响。试验二维概化模型如图2所示，模型水平距离为15 m，高度6 m，坡高4 m，初始坡度为45°，P-4,P-5,P-6为埋设于边坡内部不同位置处的孔隙水压力计。假设模型材料各项同性，水位变动条件下边坡的应力应变场与孔隙水压力场由GeoStudio中SEEP/W模块与SIGMA/W模块计算，点稳定系数分布云图由MatLab与Surfer软件联合生成。

图2 模型试验示意图Fig.2 Sketch of model test

4.1 模型本构与材料参数

水位升降数值分析中,模型本构均选用线弹性模型，岩土体类型为粉土，其有效黏聚力和内摩擦角分别为c′=1 kPa和φ′=30°，土颗粒相对密度Gs=2.69，饱和渗透系数ks=5.3×10-6m/s，天然体积含水量为0.18。变形模量与泊松比依据工程地质手册[14]取粉土经验值E=18 MPa,υ=0.3，令饱和体积含水量值等于液限值θs=wL=0.317，则粉土的天然重度经换算得γ=20.2 kN/m3。粉土体积含水量函数经工程类比法确定，其表达式为

(15)

当土体处于非饱和状态时，其体积含水量可通过式(15)得到。

4.2 水位上升

大型模型试验中，水位上升条件由可升降水箱与模型底部带出水点的管网控制，通过提升水箱，地下水经底部管网系统通过反滤层直接渗入土体，坡内外水位同时缓慢上升至5.6 m。由于水位上升条件为坡内外水位同时缓慢上升，因此在水位上升过程中浸润线始终保持水平。考虑模型初始未蓄水与最高水位2种情况，有限元网格剖分如图3所示，网格划分类型为三角形与四边形的混合类型,边界条件为底部全约束，两侧法向约束，当初始水位高程为0 m时，模型底部施加0 m总水头；当达5.6 m最高水位时，模型底部施加5.6 m总水头，并在坡面施加随水深变化的水压力边界。

图3 水位上升工况模型有限元网格剖分Fig.3 Finite element meshes of model in the condition of water level rising

首先通过SEEP/W模块中稳态流分析类型，计算边坡初始孔隙水压力分布，然后利用SIGMA/W模块中初始应力分析类型调用孔隙水压力并计算边坡应力场分布，将模型各节点的坐标信息及其对应的最大主应力、最小主应力、孔隙水压力提取出来，运用式(13)、式(14)中点稳定系数的定义,在MatLab中编程计算各节点所对应的点稳定系数值，最终插值生成坡体内点稳定系数分布云图(见图4)。

图4 水位上升前后点稳定系数分布云图Fig.4 Contours of local factor of safety before and after water level rising

从图4中可以看出边坡内部点稳定系数分布的差异性。当蓄水位为0 m时，由图4(a)可知点稳定系数随深度的增大而减小，在垂直坡面方向约0.5 m深度范围内，点稳定系数随深度的减小而增大，表明浅层坡面易发生失稳；临空面坡脚靠上部区域由于应力集中作用其稳定系数最小，安全系数位于[1.2，1.3]之间，说明坡脚处为最易破坏区域；坡顶与坡底区域点稳定系数均大于3，安全储备充足。整体上看，边坡内部点稳定系数均大于1，表明边坡处于稳定状态。

图4(b)为最高水位5.6 m时的点稳定系数分布图，此时边坡内部出现非稳定区，该区域点稳定系数小于1，其分布位置主要位于浸润线以下、坡脚以上部分；非稳定区在临空面约0.5 m深度范围内安全系数最小，为[0.8，0.9]，由临空面向坡体中部扩展的非稳定区安全系数为[0.9，1.0]，临空面处安全系数最小，表明此处将首先发生滑动破坏。水位上涨后，模型试验观测到的现象为“在水位上升过程中，填土边坡坡面逐渐发生局部崩塌，随着水位的上升坍塌的范围逐渐扩大，最终边坡坡肩退到原位置后方1 m左右”，这说明数值模拟结果与实际相符。

水位上升一方面使得边坡原始非饱和区逐渐饱和，基质吸力逐渐丧失，致使该区域的表观黏聚力C值降低。在图1中表现为非饱和摩尔-库伦强度包络线下移，削弱了边坡的抗剪强度。另一方面随着水位的上升，坡体内孔隙水压力逐渐增大，导致边坡内有效应力的降低，在图1中表现为莫尔应力圆左移，局部稳定性降低。由水位上升前后的点稳定系数云图对比可知,水位上升将导致边坡内部稳定性降低，尤以浸水区域降低幅度最为明显；临空面浅层范围内稳定系数较低，水位上升易诱发该区域浅层失稳。

4.3 水位下降

大型模型试验中，水位下降的初始水位为5.6 m，坡外水位下降速度为1 m/h，经156 min，坡外水位从5.6 m降至3.0 m。水位的上升使得原始边坡临空面浅层范围处于非稳定区域，导致坡面发生局部崩塌，并随水位上升坍塌范围扩大，填筑边坡的坡度由初始状态的1∶1(45°)逐步变为约1∶1.5(约33°)，因此水位下降工况下边坡的初始坡度为33°。重新建立数值模型并剖分网格如图5所示，有限元网格划分类型为三角形与四边形的混合类型。为便于观测边坡内部各参数的变化情况，在边坡中部设置监测线A，其顶点坐标为(7,4)，深度为2 m。模型边界条件为底部全约束，两侧法向约束，坡面同时施加水位变化边界以及水压力变化边界。

图5 水位下降工况模型有限元网格剖分Fig.5 Finite element meshes of model in the condition of water level drawdown

模型初始应力计算方法与第4.2节相同。水位下降过程模拟选用SIGMA/W模块中的流固耦合分析类型，定义模型的体积含水率函数及渗透系数函数，调用初始应力场与初始孔压场，以1 min为步长，计算整个水位骤降过程的应力应变场与孔隙水压力场。

监测点处孔隙水压力计算值与实测值的对比如图6所示，总体上看计算值与实测值变化趋势较为吻合，说明模型计算结果具有一定合理性。随着坡外水位的下降，孔隙水压力的实测值稍大于计算值，其原因在于边坡在不断的破坏解体过程中，阻断了原有的渗流通道，使得水流不畅，从而导致孔隙水压力实测值大于计算值。

图6 孔隙水压力实测值与计算值对比Fig.6 Comparison of pore water pressure between measurement and calculation

在实际模型试验中第40 min与第100 min分别产生了滑块1,2的滑动，因此重点研究0,40,100,156 min时刻所对应的稳定性情况，并绘制各时刻点稳定系数分布云图。为便于对比点稳定系数与整体稳定系数变化情况，将模型试验中观测到的3个滑面即滑面1,滑面2，滑面3在SLOPE/W中定义，计算各时段每个滑面的整体稳定性系数大小，计算方法选用滑面应力法、Morgenstern-Price(M-P)条分法以及Janbu条分法。点稳定性系数分布云图与整体稳定性计算结果见图7。

从图7中可看出，随着水位的下降，边坡内部点稳定系数降低，且非稳定区域逐步扩大。

(1) 初始时刻时，由点稳定系数分布云图中可知，浸润线以下边坡临空面0.2 m深度范围内以及边坡中后部区域处于非稳定状态。3条潜在滑面均未穿过非稳定区，不同计算方法计算的各滑面整体稳定性系数均在1.20以上，表明此时边坡整体处于稳定状态。

(2) 在第40 min时，水位下降至4.9 m，浸润线向坡外弯曲，说明边坡内部水位下降速度滞后于坡外，这将产生一定的动水压力，加大边坡下滑力，降低边坡稳定性。由点稳定系数云图可知非稳定区域分布面积较初始时刻明显扩大，坡脚靠上部区域稳定性最低，初始时边坡中后部与坡面处的非稳定区域在此时已形成贯通。各滑面部分穿过非稳定区，整体稳定性系数与初始时刻相比降低0.2左右，稳定性下降较为明显。

(3) 当试验进行至100 min时，浸润线向坡外弯曲程度增大，动水压力增大；非稳定区域分布面积进一步扩大，其分布位置主要位于浸润线以下坡脚以上；3条滑面大部分处于非稳定区，整体安全系数进一步降低，安全系数均<1，表明整体处于失稳状态。

第156 min时，水位降至3 m并停止下降，浸润线弯曲程度达到最大，非稳定区域分布面积与各滑面穿过非稳定区长度也达到最大，整体稳定性降至最低；非稳定区分布仍主要位于浸润线以下坡脚以上区域，而该区域为动水压力的主作用区，说明动水压力是导致边坡稳定性下降的主因。坡脚处的点稳定系数依旧为全局最低区域，说明在水位下降过程中，坡脚处为边坡最易发生破坏部位，易引致边坡发生牵引破坏。

图8为监测线A各时段所对应的水力梯度与点稳定系数。因监测剖面顶点高度为4 m，则在0～100 min之内，A均位于浸润线以下。

图8 不同时间下监测线A处水力梯度与点稳定系数曲线Fig.8 Profiles of hydraulic gradient and local factor of safety at different instances at monitoring line A

从图8中可看出，在此时间段内，随着水位的下降，水力梯度逐渐增大，点稳定系数逐渐减小。相同时刻下水力梯度随深度的增大而减小。当时间超过100 min时，监测剖面逐渐露出水面，由于露出水面部分由饱和状态转为非饱和状态，抗剪强度增大，水力梯度减小，点稳定系数有所回升。总体而言，水力梯度与点稳定系数呈正相关，水力梯度的变化对点稳定系数的改变具有直接影响。

5 结 论

(1) 以模型试验为基准，采用非饱和流固耦合理论与点稳定系数相结合的方法,分析了水位升降作用下非饱和河岸边坡各部位的稳定性分布及其动态变化过程，数值分析结果与实际较为吻合，对河岸边坡的稳定性评价及治理具有借鉴意义。

(2) 边坡临空面浅层范围内，点稳定系数随深度减小而逐渐降低，水位上升将导致边坡内有效应力降低，非饱和区基质吸力丧失，抗剪强度减弱，易诱发边坡浅层失稳。

(3) 水位骤降产生的动水压力是边坡失稳的主因，随着水位下降，浸润线向坡外弯曲程度逐渐增大，整体稳定性下降，非稳定区域分布面积逐渐增大，其分布位置主要位于浸润线以下、坡脚以上区域；坡脚区稳定系数最低，易导致边坡发生牵引破坏；边坡浸润线以下部分的点稳定系数随水力梯度不断增大而逐渐减小，浸润线以上部分由于饱和到非饱和状态的过渡，水力梯度减小，点稳定系数回升。

需要指出的是，文中的研究对象为松散堆积体均质边坡，结构较为简单，对具有复杂工程地质条件河岸边坡的失稳模式及稳定性演化规律仍需进一步研究。另外文中的耦合分析所选用的线弹性模型并不能考虑边坡塑性破坏后应力重分布的情况，因此对于边坡破坏过程的模拟存在不足。

[1] 朱 朋,卢书强,薛聪聪,等.库水位升降与降雨条件下滑坡的渗流及稳定性分析[J].长江科学院院报,2015,32(11):87-92.

[2] 周永强,盛 谦.库水位变化和降雨作用下付家坪子高陡滑坡稳定性研究[J].长江科学院院报,2014,31(2):57-61.

[3] 沈细中,杨文丽,兰 雁.小浪底库岸1#滑坡体安全评价[J].岩石力学与工程学报,2011,30(3):589-595.

[4] 倪卫达,唐辉明,胡新丽,等.黄土坡临江I号崩滑体变形及稳定性演化规律研究[J].岩土力学,2013,34(10):2961-2970.

[5] WANG Jing-e,XIANG Wei,LU Ning.Landsliding Triggered by Reservoir Operation:A General Conceptual Model with a Case Study at Three Gorges Reservoir[J].Acta Geotech,2014,9(5):771-788.

[6] 周建烽,王均星,陈 炜.非饱和非稳定渗流作用下边坡稳定的有限元塑性极限分析下限法[J].岩土工程学报,2014,36(12):2300-2305.

[7] 张 旭,谭卓英,周春梅.库水位变化下滑坡渗流机制与稳定性分析[J].岩石力学与工程学报,2016,35(4):713-723.

[8] 马建全,王念秦,张新社.基于点稳定系数法的斜坡稳定性分析[J].吉林大学学报(地球科学版),2015,45(2):533-540.

[9] LU Ning,SENER-KAYA B,WAYLLACE A,etal.Analysis of Rainfall-induced Slope Instability Using a Field of Local Factor of Safety[J].Water Resources Research,2012,48(9),doi:10.1029/2012WR011830.

[10] FREDLUND D G,RAHARDJO H.Soil Mechanics for Unsaturated Soils[M].New York:John Wiley and Sons,Inc.,1993.

[11] KRAHN J.Seepage Modeling with SEEP/W 2007:An Engineering Methodology[M].2nd ed.Calgary,Alberta,Canada:Geo-Slope International Ltd.,2010.

[12] VANAPALLI S K,FREDLUND D G,PUFAHL D E,etal.Model for the Prediction of Shear Strength with Respect to Soil Suction[J].Canadian Geotechnical Journal,1996,33(3):379-392.

[13] 贾官伟,詹良通,陈云敏.水位骤降对边坡稳定性影响的模型试验研究[J].岩石力学与工程学报,2009,28(9):1798-1803.

[14] 常士骠.工程地质手册[M].第4版.北京:中国建筑工业出版社,2006.

(编辑：赵卫兵)

Stability of Unsaturated Bank Slope Based onLocal Factor of Safety

WANG Yan-kun，XIONG Cheng-ren，WANG Jing-e，ZOU Zong-xing，HAN Tian-hua

(Three Gorges Research Center for Geo-hazard under Ministry of Education,China University of Geosciences,Wuhan 430074,China)

Overall stability calculation for river bank slopes could not reflect the impact of water level fluctuation on various parts of the slope stability quantitatively.On the basis of the theory of unsaturated soil mechanics and local factor of safety method (LFS),we conducted numerical analysis on the stability of river bank slope by the model of unsaturated liquid-solid coupling.Taking a large-scale model test for example,we plotted the contours of LFS at different instances under water level fluctuation to reveal the causes of destruction and the processes of internal stability.Compared with actual observations,the analysis result demonstrates that 1) the LFS method could describe the river bank slope stability and its dynamic changes at various slope positions reasonably; 2) the free face of slope is of low stability,and water level rise would reduce the effective stress and lead to the loss of matrix suction in the unsaturated zone of slope,which is easy to induce shallow slide; 3) the hydrodynamic pressure formed by rapid drawdown is the main cause of the decline of slope stability and the most vulnerable zone is the toe of slope where retrogressive failures often emerges; 4) along with the increase of hydraulic gradient,the internal local stability coefficient decreases and the instability zone,which distributes mainly between the saturation line and the toe of slope,expands.

river bank slope; water level fluctuation; unsaturated liquid-solid coupling;numerical analysis; local factor of safety

P642.22

A

1001-5485(2017)10-0124-06

2016-06-17;

2016-07-05

国家自然科学基金项目(41272309)

王艳昆(1990-)，男，河南新乡人,硕士研究生，主要从事地质灾害的评价与治理研究，(电话)15927508327(电子信箱)wykyynl@foxmail.com。

熊承仁(1965-)，男，湖北武汉人，副教授，博士，主要从事岩土工程与工程地质研究，(电话)027-67885020(电子信箱)xiongcr@126.com。

10.11988/ckyyb.20160627 2017,34(10):124-129