边坡稳定分析中有限元法影响因素探讨

王 静

(中国电建集团西北勘测设计研究院有限公司，西安 710065)

0 前 言

在中国水电工程建设中，边坡问题尤为突出，可能成为工程建设的制约性因素[1]。目前在工程设计中普遍以现行规范中规定的刚体极限平衡法分析边坡稳定为主，而因数值分析法本身的特点加之边坡自身的难点，工程界专家和学者虽做了大量的研究工作，一直未达成统一的认识。数值分析中有限元法应用最广，得到工程界认可，但在边坡稳定性分析中同样存在上述问题。

本文研究了数值分析法中有限元法在边坡稳定性分析中的关键技术。在总结前人已有研究成果的基础上，通过典型工程案例分析，系统地研究了边坡稳定性分析中有限元法的模型范围、网格尺寸、边界条件、收敛容差、收敛准则及判别准则等问题，得出了统一的认识及一般规律性的结论，具有建设性的指导意见，可供边坡稳定性分析中参考、借鉴。

1 基本原理

1.1 强度折减法简介

有限元法得到的结果是边坡岩土体中各点的变形、应力、应变以及各点应力状态下的塑性屈服状态，这对了解边坡结构的实际运行状态有很大的帮助。但目前中国现行规范中给出的边坡稳定性安全评价标准为较直观的刚体极限平衡法下的安全系数，这在数值分析中很难做到，因此出现了后来在有限元法基础上的强度折减法。

所谓的有限元强度折减法[2-4]是通过降低强度，使系统达到不稳定状态，有限元计算不收敛，此时的折减系数就是边坡稳定安全系数。它是由有限元和强度折减相结合的方法对边坡的稳定进行分析。有限元强度折减法原理[5-6]:不断降低岩土c、φ值，直到破坏。至破坏时，c、tanφ降低倍数就是安全系数。

(1)

(2)

(1) 边坡有限元强度折减法与传统边坡刚体极限平衡法相比之优势[7]

1) 有限元强度折减法充分发挥了用数值计算方法计算水电水利工程边坡稳定性的优势，较以往的计算方法更加快捷、方便，为岩土边坡从业人员提供了更好的技术支持；

2) 有限元强度折减法能够考虑岩土体的非线性本构关系、开挖和支护结构的作用效应等；

3) 有限元强度折减法不需要做太多假定；

4) 通过分析可以直观地反映坡体的实际滑动面；

5) 强度折减法可以获得边坡的应力、变形分布等信息，适用于复杂地质条件下边坡稳定性分析[5]。

(2) 有限元强度折减法自身存在的问题

1)c值和φ值等比例折减问题。

尚无定论，本文按同时等比例折减进行考虑。

2) 多个滑动面的边坡，折减范围问题。

本文按研究对象(单一滑动面)进行折减。

3) 折减至计算不收敛还是发散问题。

本文按计算至不收敛考虑，折减程度即为安全系数。

1.2 有限元法的影响因素

有限元法在边坡稳定性分析中的关键技术包括：① 模型范围;② 网格尺寸;③ 边界条件;④ 收敛容差;⑤ 收敛准则;⑥ 判别准则。这6种因素将对边坡稳定性分析成果造成一定影响，甚至影响到边坡稳定性的判断。

本文通过某一典型工程作为案例进行分析，得出一般规律性的结论，可降低或消除有限元法分析过程中人为因素的影响。

2 工程实例

2.1 工程概况

某水电站位于澜沧江河段，坝址两岸岸坡自然坡度一般为35°～55°，局部60°～75°。两岸冲沟发育，地表植被较丰茂。右岸坝肩高程1 900.00 m以下为陡壁，基岩裸露，高程1 970.00～2 095.00 m有第四系坡积碎石土层，顺河向宽度约130 m，表层经后期改造为耕植土，此处覆盖层最大厚度约32 m，结构松散，因此有必要研究右岸坝肩覆盖层滑坡体和倾倒体滑坡体的稳定性问题，确保工程的长期运行安全性。

本文以右岸坝肩的覆盖层边坡为工程实例，采用有限元强度折减法对其进行稳定性分析。有限元计算主要针对右岸坝肩覆盖层滑坡体和倾倒体滑坡体，具体部位如图1。

图1 右坝肩滑坡体具体部位图

2.2 计算参数

计算参数如表1。

表1 边坡岩体的物理力学参数取值表

3 计算结果

3.1 模型范围敏感性分析

传统极限平衡法对边界范围只需要将破坏面包括在内，计算结果就不会有影响，但有限元分析中边界范围的大小会造成计算结果发生变化。本节计算主要对模型计算范围进行敏感性分析，对于网格尺寸暂不做讨论，为提高计算效率，网格尺寸初选4 m。

3.1.1 模型简介

(1) 模型1简介

模型原始边界，模型Z向长405 m，Y向高533 m；网格尺寸均取4 m；侧面采用法向约束，底面采用全约束，模型见图2。

图2 有限元模型1

(2) 模型2简介

模型左、右边界各增加0.5倍的倾倒体滑体水平宽度(200 m)，竖直方向增加1倍的倾倒体滑体水平宽度(400 m)；模型Z向长805 m，Y向高933 m；网格尺寸均取4 m；侧面采用法向约束，底面采用全约束，模型见图3。

(3) 模型3简介

模型左、右边界各增加1倍的倾倒体滑体水平宽度(400 m)，竖直方向增加2倍的倾倒体滑体水平宽度(800 m)；模型Z向长1 205 m，Y向高1 333 m；网格尺寸均取4 m；侧面采用法向约束，底面采用全约束，模型见图4。

(4) 模型4简介

模型左、右边界各增加1.5倍的倾倒体滑体水平宽度(600 m)，竖直方向增加3倍的倾倒体滑体水平宽度(1 200 m)；模型Z向1 605 m，Y向1 733 m；网格尺寸均取4 m；侧面采用法向约束，底面采用全约束，模型见图5。

图3 有限元模型2

图4 有限元模型3

图5 有限元模型4

3.1.2 敏感性分析

(1) 位移分析

采用程序默认的最大不平衡比率收敛值(1×10-5)进行计算分析，各模型位移变化对比见表2。各模型范围对Y、Z向位移的影响见图6～7。

表2 各模型位移分析表

图6 模型范围对Y向位移的影响

图7 模型范围对Z向位移的影响

由位移变化结果可知，模型1(原始边界)的Y向位移误差精度控制在6%，Z向位移误差精度控制在3%；模型2(左、右边界各增加200 m，竖直方向增加400 m)的Y、Z向位移误差精度均控制在2%左右；模型3(左、右边界各增加400 m，竖直方向增加800 m)的Y、Z向位移误差精度均控制在1%左右。

(2) 屈服状态分析

由计算结果可知，模型1～4的覆盖层滑坡体部位全部呈现屈服状态，倾倒体滑坡体基岩分界面中上部出现塑性屈服，下部无屈服区域出现，边界范围变化对塑性区屈服状态的变化影响不大。

结合模型1～4的位移分析和塑性屈服状态分析可知，模型范围变化对塑性区屈服状态的变化影响不大，但对滑坡体变位影响较大，故采用位移分析结果。

位移分析结果可知，模型3的Y、Z向位移误差精度均控制在1%左右，故选取模型3是合理的。

3.2 网格尺寸敏感性分析

网格尺寸是建立有限元模型的重要环节之一，会对计算精度和计算规模产生显著的影响。结合上一节分析可知，模型3边界范围是合理的，故以模型3为基础进行网格敏感性分析，即左、右边界各增加1倍400 m，竖直方向增加800 m。

3.2.1 滑坡体网格尺寸敏感性分析

本节计算整个模型的网格尺寸分别取2、4及8 m。变位数值见表3，网格尺寸对Y、Z向位移的影响见图8～9。

表3 各模型位移分析表

图8 滑坡体的网格尺寸对Y向位移的影响图

图9 滑坡体的网格尺寸对Z向位移的影响图

由位移变化结果可知，网格尺寸为2 m时，Y向位移误差精度控制在0.23%，Z向位移误差精度控制在0.08%；网格尺寸为4 m时，Y向位移误差精度控制在0.24%，Z向位移误差精度控制在0.16%。为保证计算对象的计算精度，计算对象的网格尺寸取2 m，网格尺寸(2 m)为滑体水平宽度(400 m)的0.5%。

3.2.2 基岩网格尺寸敏感性分析

由3.2.1节分析可知，滑坡体网格尺寸取2 m。为提高计算效率，本节计算基岩网格尺寸分别取倾倒体滑体水平宽度(400 m)的1/10、1/20及1/200，即40 m、20 m及2 m。变位数值见表4，网格尺寸对Y、Z向位移的影响见图10～11。

由位移变化结果可知，网格尺寸为2 m时，Y向位移误差精度控制在0.93%，Z向位移误差精度控制在0.24%；网格尺寸为20 m时，Y向位移误差精度控制在4.7%，Z向位移误差精度控制在1.2%。为提高计算速度，非计算对象的网格尺寸的计算误差可控制在5%。所以基岩的网格尺寸宜选定20 m，不超过计算对象网格尺寸的10倍。

表4 各模型位移分析表

图10 基岩的网格尺寸对Y向位移的影响图

图11 基岩的网格尺寸对Z向位移的影响图

为保证计算精度和计算速度，滑坡体的网格尺寸取2 m，网格尺寸(2 m)为滑体水平宽度(400 m)的0.5%；基岩的网格尺寸宜选定20 m，不超过计算对象网格尺寸的10倍。

3.3 边界条件敏感性分析

结合上述分析，模型3是合理的。滑坡体网格尺寸取2 m，基岩网格尺寸取20 m，在此基础上分析基底全约束和法向约束对计算结果的影响，对应的位移分析见表5。

表5 位移分析表

由表5可知，2种约束模型所得计算结果相近,计算误差控制在0.3%，所以认为约束条件对计算结果影响很小。

3.4 收敛容差敏感性分析

综合前面各节分析结果，在选定的模型范围、网格尺寸、模型边界条件下，现对收敛容差的影响做敏感性计算分析。收敛容差除影响计算精度外，还对计算效率影响很大，有限元程序默认的收敛准则为能量收敛准则(见表6)。

表6 验证模型位移对比表(能量收敛准则)

通过5种收敛容差拟定方案分析成果可知：

(1) 从收敛容差与模型变位关系曲线(见图12)可知，当收敛容差低于1.0×10-3时，模型最大变位变化幅度不大，因此要保证计算精度，收敛容差不得高于1.0×10-3。

图12 能量收敛准则下模型最大变位与收敛容差关系曲线图

(2) 本文推荐边坡稳定分析中，能量法收敛容差设置为1.0×10-5。

3.5 收敛准则敏感性分析

数值分析计算中，大多采用能量法收敛准则。边坡结构为弹塑性体，数值模型为弹塑性模型，容易发生塑性流动，属于材料软化方面的求解。位移收敛准则的模型位移对比见表7。

表7 验证模型位移对比表(位移收敛准则)

通过5种收敛容差拟定方案分析成果可知：

(1) 位移收敛准则与能量收敛准则下边坡最大变形差别不大。

(2) 不同量级收敛容差下的位移收敛较能量法稳定性要更好、曲线更平稳,见图13。

(3) 综合分析，采用位移收敛准则，当收敛容差低于1.0×10-3时，模型最大变位结果一致。

图13 位移收敛准则下模型最大变位与收敛容差的关系曲线图

3.6 不同判断准则安全系数分析

采用有限元强度折减法分析边坡稳定性过程中，安全系数的选取取决于失稳标准的选取，但目前在边坡稳定分析中并没有统一标准，常用失稳判据标准主要有计算不收敛标准、塑性贯通标准及位移突变标准[8-9]。

本节主要利用强度折减法，分别采用上述3种失稳判据标准对坝肩倾倒体的安全系数进行分析，讨论其合理性与实用性。

模型范围：左、右边界各增加1倍的倾倒体滑体水平宽度400 m，竖直方向增加2倍的倾倒体滑体水平宽度800 m；滑坡体的网格尺寸取2 m，基岩的网格尺寸宜选定20 m；底部采用法向约束；收敛准则采用位移法，收敛容差采用1.0×10-5。

3.6.1 覆盖层滑坡体

(1) 计算不收敛标准

天然工况下，当折减系数为1.078时，计算不收敛，故覆盖层滑坡体安全系数为1.077。

(2) 塑性区贯通标准

因坝肩覆盖层滑坡体本身处于临界状态，故在天然状态下，覆盖层滑坡体屈服状态本身就处于整个塑性区全部贯通，因此无法用塑性区贯通作为判断准则。

(3) 位移突变标准

以位移突变作为边坡破坏标准，选取覆盖层滑坡体坡脚处节点作为控制点，可得不同折减系数所对应的位移值，位移数值见表8，位移与安全系数关系曲线见图14。

由图14可知，采用位移突变标准时，当折减系数为0.9时，位移发生突变，即认为安全系数为0.9。

表8 覆盖层滑坡体位移数值表

图14 覆盖层滑坡体位移与安全系数关系曲线图

通过采用计算不收敛标准及位移突变标准对坝肩覆盖层滑坡体的安全系数进行分析，所得安全系数分别为1.077及0.9，即计算不收敛所得安全系数相对较大。同时，因为覆盖层滑坡体本身处于塑性区贯通状态，所以无法用塑性区贯通标准取得安全系数。

3.6.2 倾倒体滑坡体

(1) 计算不收敛标准

天然工况下，当折减系数为1.23时，计算不收敛，故倾倒体滑坡体安全系数为1.22。

(2) 塑性区贯通标准

天然状态下，倾倒滑坡体本身处于塑性区从坡脚到坡顶贯通的状态，因此无法用塑性区贯通作为判断准则。

(3) 位移突变标准

以位移突变作为边坡破坏标准，选取倾倒体滑坡体坡脚处节点作为控制点，可得不同折减系数所对应的位移值，位移数值见表9，位移与安全系数关系曲线见图15。

由图15可知，采用位移突变标准时，当折减系数为1.21时，位移发生突变，即认为安全系数为1.21。

表9 倾倒体滑坡体位移数值表

图15 倾倒体滑坡体位移与安全系数关系曲线图

通过采用计算不收敛标准及位移突变标准对右岸坝肩倾倒体滑坡体的安全系数进行分析，所得安全系数分别为1.22和1.21，即计算不收敛所得安全系数相对较大。同时，因为倾倒体滑坡体本身处于塑性区贯通状态，所以无法用塑性区贯通标准取得安全系数。

计算不收敛：计算不收敛所得安全系数与刚体极限平衡法的误差在5%左右，此计算精度可满足工程需要。

塑性区贯通：本工程数值计算中存在大片的塑性区，但这些进入塑性的单元并没有产生无限制的塑性流动，而是处于塑性极限平衡状态。所以，塑性区贯通并不一定意味着边坡整体破坏，塑性区贯通是破坏的必要条件，但不是充分条件。

4 结 语

本文主要针对目前边坡稳定性分析中有限元强度折减法做了大量的研究工作，总结分析了边坡设计计算中的关键技术，研究成果可供工程设计提供参考和借鉴。主要得出如下结论。

(1) 模型范围：横河向两侧各伸出计算对象1倍尺寸，竖直向延伸至计算对象2倍深度，模型精度相对误差可控制在1%以内。

(2) 网格尺寸：计算对象网格尺寸控制在计算对象尺度的1/100，可满足精度要求，误差在1‰以内。计算对象以外网格尺寸可适度放大些，可取计算对象网格尺寸的3～5倍为宜。

(3) 边界条件：底部边界法向或固结约束对计算结果影响甚微，误差万分之一量级。

(4) 收敛容差：有限元法默认收敛准则为能量收敛准则，在边坡稳定分析中收敛容差不宜低于1.0×10-3，建议采用1.0×10-5。

(5) 收敛准则：位移收敛准则与能量收敛准则下边坡最大变形差别不大，不同量级收敛容差下的位移收敛较能量法稳定性要更好。

(6) 塑性区分布范围人为因素不易控制，对边坡稳定性安全系数难以定量分析判断，因此不建议采用塑性区贯通标准进行边坡稳定性评价。

(7) 采用计算不收敛和位移突变分析边坡的稳定性均是合适的，但计算不收敛法人为因素较多，不易控制，因此建议采用边坡表部特征点位移突变标准做定量分析。具体为：不同折减系数下的特征点位移突变时所对应的折减系数，即为边坡稳定性安全系数。