基于有限元重度增加法的涉水岸坡稳定性分析研究

刘中才

(西藏众兴工程设计咨询有限公司，西藏 山南 856100)

当前，已建和在建的水利工程数量数以万计，我国水库数量逐渐增加[1]。水库建成蓄水后水位逐渐上升，造成山体长期浸泡在水中，水位反复升降直接削弱了山岩土体的各类力学强度，非常容易造成滑坡、泥石流等地质灾害，这直接影响到水利枢纽工程的运行安全[2-4]。

因此，评估不同水位时涉水岸坡岩体的各项非线性特征显得尤为重要。极限平衡法包含瑞典圆弧法、费伦纽斯(Fellenius)法和简布(Janbu)法等[5-7]。而数值分析法包含了有限元法和离散元法等[8-9]。刚体极限平衡法由于仅能计算出稳定安全系数和滑弧，难以分析边坡内部的应力应变特征而受到一定的限制。而数值分析法，按相同比例折减岩体内摩擦角和黏聚力[10]，理论上存在缺陷。本文结合了数值有限元分析法，提出采用重度增加法计算分析某涉水岸坡的稳定性，旨在为该边坡工程的安全评估提供参考和借鉴。

1 研究方法

1.1 渗流计算原理

有限元法渗流计算基本机理：根据渗流作用力概念，将作用在滑动面上和划分土块表面的水压力转换为等价的体积力。把各节点水头值换算成各单元渗透力，不需要考虑各单元体接触面边界上的孔隙水压力和边坡外水压力。滑坡体中的地下水位随水库、河谷水位升降而变化，因而在滑坡体中就会形成饱和区和非饱和区。非饱和区土壤水与饱和区地下水运动密切相关，因此，由于水流运动导致边坡岩土体状态改变即为饱和与非饱和渗流问题，其渗透各向异性的饱和与非饱和渗流控制方程如式(1)、式(2)所示[11-13]：

(1)

(2)

式中：x、y表示位置坐标；Kx、Ky分别为岩土体的水平和垂直方向的渗透系数；ρw为水的密度；mw为释水系数；Q为边界流量；θw为体积含水率，ua为孔隙气压力，uw为孔隙水压力，h为压力水头，t为时间。

1.2 稳定计算原理

计算工具采用岩土有限元软件PLAXIS，对该涉水岸坡的应力应变和稳定性进行定量计算和分析。该软件包含丰富的岩土体本构模型，可以应用于各类计算工况[14]。

鉴于当前边坡稳定性方法均存在一定的理论缺陷问题，本文中边坡稳定性分析时结合有限元数值分析法，采用有限元重度增加法[15-16]来分析边坡的稳定性。由于有限元重度增加法仅需要梯度增加重度大小数值，然后不断循环迭代计算，直到有限元计算不收敛或者节点位移突变作为判别条件，此时重度数值增加的倍数即定义为边坡的安全系数，相比于强度折减法更为可靠，重度增加法计算式如式(3)所示：

(3)

式中：Fs为重度增加法定义的安全系数，γ为自然条件下边坡岩土体的重度，γ′为根据梯级倍数增加后达到临界状态时边坡岩土体的重度。

2 结果与分析

2.1 计算模型前处理

本次研究的边坡对象位于西藏江达县境内，由于含泥量较高，属于强风化山体，山体的结构主要包含土、碎石和块石。该边坡顶部高程约3650 m，坡脚高程约2930 m，边坡坡度平均33°。岸坡山体表层为碎石土，最大厚度超过50 m，边坡基础部分为弱风化岩体。该涉水边坡剖面及网格划分情况见图1。参照岩土试验结果，获得了该岸坡碎石土和弱风化基岩的有限元计算力学参数，见下表1所示。

表1 涉水岸坡模型有限元计算力学参数

采用有限元重度增加法计算前，设置了不同水位的计算工况，分析了三种不同库水位，高程分别为2950 m、2930 m、2910 m。本次计算主要考虑稳定渗流工况，不计算非稳定渗流，即认为岩土体内部水位随着外面水位变化，与之保持一致，不存在滞后的现象。采用重度增加法计算边坡稳定性时同时涉及表层碎石土和弱风化基岩的重度参数改变。岩土体材料模型选用摩尔库伦模型，从而可以考虑涉水岸坡的非线性特征。

2.2 岸坡应力应变分析

分析计算了库水位2950 m、2930 m、2910 m时该涉水边坡的应力应变特征。其中，应力为剪切应力，应变为剪切总应变。由于坡体强度主要由抗剪强度进行反应和体现，因此，边坡应力应变分析以剪应力和总应变为主。图2为不同水位时边坡的剪应力分布情况，分析可知，三种情况下山体岩体内部的剪应力分布情况比较一致，表层碎石土整体出现应力松弛现象，剪应力非常小，介于200～-200 kPa之间，正值受拉，负值受压。当水位为2950 m和2930 m时，在高程3430 m附近出现了较大的正值剪应力，大约1900 kPa，此处岩层分界坡度较陡，极易出现断层错位。当水位为2910 m时，在高程3430 m处的剪应力相对减小，仅为1800 kPa。三种水位下剪应力的最大值均位于模型底部居中位置，越往山体内部走剪应力数值越小，这对该岸坡的表层滑动影响很小，基本可以忽略不计。以上分析说明，在不同蓄水位时涉水岸坡的应力状态差别很小，且表层碎石土应力松弛现象较为明显，当在外界自然条件作用下，在坡度较陡的位置容易发生局部垮塌，应引起关注。

图2 不同水位时涉水岸坡剪应力分布

图3为不同水位时边坡的总应变分布情况，阴影部位范围大小反映了总应变大小。分析可知，三种水位情况下山体岩体内部的总应变分布情况差别较大。库水位为2950 m时总应变分布范围最大，主要分布在岩层分界线两侧，剪应变集中高程在3380～3480 m，3380 m以下往坡脚位置总应变范围逐渐减小。当库水位为2930 m时此处的总应变分布量相比2950 m时大大减小，但坡体中部基本没有发生应变。当库水位为2910 m时，在高程3380～3480 m的总应变分布量与2930 m时相差不大，但坡脚位置附近的剪应变基本消失，中部同样也没有发生剪切应变。分析说明，在不同蓄水位时涉水坡体内部的剪应变分布范围区别相差较大，主要体现在涉水边坡的中上部位、岩层分界线坡度较陡的位置。

图3 不同水位时涉水岸坡总应变分布

2.3 岸坡稳定性分析

基于数值有限元分析基本原理，采用重度增加法计算了不同水位时该岸坡的稳定安全系数，见表2。库水位在2950 m时，岸坡的安全系数仅为1.145，稳定性较差。库水位2930 m时，边坡稳定性系数为1.248，库水位为2910 m时，边坡的安全系数最大，为1.259。在稳定渗流条件下，蓄水位越高，涉水边坡的稳定性反而越差，因为蓄水位越高，水体浸泡的边坡岩体也越多，导致岩体抗剪强度降低。同时，根据重度增加法的计算结果，导出了三种水位时涉水岸坡的最不利滑裂面，发现不同水位时的边坡滑裂面完全相同，滑裂面起于边坡临空面顶部，内部与岩层分界线相切，剪出口位于坡面中下部。图4仅给出了蓄水水位2930 m时该涉水岸坡的最不利滑裂面。

表2 涉水岸坡稳定安全系数

图4 岸坡水位2930 m时最不利滑裂面

3 结 论

本文针对西藏某涉水岸坡的稳定性问题，以有限元软件PLAXIS为计算平台，结合数值有限元分析法，提出采用有限元重度增加法分析边坡在不同蓄水位时的岩土体应力应变特征，并且定量计算不同水位工况时的岸坡稳定安全系数，得到以下结论。

(1)不同蓄水位时涉水岸坡的剪应力状态差别很小，仅在岩层分界线处存在局部差异。

(2)不同蓄水位时坡体内部剪应变分布区别相差较大，蓄水位越高，剪应变分布量越大，在最高水位2950 m时剪应变集中在高程3380～3480 m的岩层分界线处。

(3)当山体内部处于稳定渗流状态时，蓄水位越高，边坡的稳定性越差，最高水位和最低水位安全系数相差可达10%。

(4)三种水位时的边坡滑裂面完全相同，起于边坡临空面顶部，内部与岩层分界线相切，剪出口位于坡面中下部。

证明了该岸坡在高水位时稳定性欠佳，遇到暴雨地震等特殊工况时存在滑动的可能，且滑动模式属于高位滑坡，因此该涉水岸坡的治理应重点关注坡体中上部位。