一种新式突变判据用于三维边坡稳定性分析

谢旭奎, 张 杰

(湖南致力工程科技有限公司，湖南 长沙 410208)

有限元强度折减法由Zienkiewice[1]于1975年提出，国内外许多学者在这方面做了大量的工作[2-12]。目前对边坡稳定性的有限元分析大多都集中于二维模型，而二维模型分析的结果推广到实际边坡应用时，往往存在较大的误差，因此当前多采用三维模型分析边坡稳定性。目前国内外采用的边坡有限元分析软件主要为FlAC 3D、hyperworks、algor、nastran、abaqus等，虽然ANSYS软件操作复杂，岩土方面运用较少，但ANSYS在模型建立和网格划分优势明显，因此本文采用ANSYS进行边坡稳定性分析。基于强度折减法的边坡临界破坏状态主要有以下传统判据[7-10]:①边坡破坏的一个重要依据是有限元计算是否收敛，当计算收敛时，边坡稳定;计算不收敛，则边坡破坏。②塑性区是否贯通是边坡失稳的另外一个标志，在有限元计算中可以得到不同折减强度下的边坡塑性应变图，当塑性区贯通时，认为边坡破坏。③位移的急剧增加。在计算中当发现边坡的滑动面位移急剧增加时，认为边坡破坏。许多学者对坡体塑性区贯通状态、有限元计算收敛性、位移值突变进行研究，而对坡体的塑性应力应变值变化进行较少研究，且对位移值突变研究仅限于定性的考虑位移突变为无穷大的可能性，但实际模拟中，随着折减系数的增大，边坡位移突变到某定值后会存在屈服下降，针对此种情况，本文提出一种新式突变判据，即在位移和塑性应力应变值两类指标上，采用新的曲线拟合方式，从理论上分析得出边坡极限平衡状态下的安全系数。另外对边坡的等效塑性区变化进行全过程模拟，在以往边坡有限元分析中也比较少见。

1 有限元强度折减法

1.1 基本原理

有限元强度折减法的基本原理是通过不断的调整岩土体的抗剪强度指标，即将粘聚力和内摩角、泊松比、弹性模量等用一个折减系数Ki进行折减，得到折减后的新抗剪强度指标，然后通过弹塑形有限元数值计算确定边坡内的应力场、应变场或位移场，并且基于对应变或者位移的某些分布特征以及有限元中的某些数学特征进行分析判断，不断增加折减系数Ki的值，直至边坡失稳，将此极限状态下所对应的折减系数定义为边坡的稳定安全系数。

1.2 参数调整

根据整体强度折减，内摩擦角和黏聚力调整如下：

(1)

弹性参数(弹性模量和泊松比)调整强度折减时，

对于泊松比调整如下：

(2)

对于弹性模量调整如下：

(3)

1.3 屈服准则

采用理想弹塑性本构模型的Drucker-Prager屈服准则。对于理想的弹塑性模型，土体进入屈服阶段就意味着破坏，将Mohr-Coulomb准则作为屈服准则变换得到Mohr-Coulomb屈服条件，在三维主应力空间中的屈服函数的表达式如(4)所示，Drucker-Prager屈服准则是对Mohr-Coulomb准则不收敛的改进，以此来修正Von Mises准则，ANSYS自带DP1屈服模块(Mohr-Coulomb准则外接点外接圆)，表达式(4)中的b和sy为DP1准则修正后的材料常数。

(4)

1.4 边坡失稳判据

在通过对X方向塑性应变、XY方向的剪切塑性应变、塑性应变强度，等效塑性应变，Von Mises塑性应变5个塑性指标和空间位移、x方向位移、y方向位移3个位移指标的分析的基础上(将坡体整体位移值与塑性应力应变值中各自的最大变化值作为研究值)，采用Gauss、Lorentz、Vigot曲线进行非线性拟合，定量获得突变的极大值，当最大变化值达到极大值时，边坡发生破坏，为此当位移或塑性指标无限靠近极大值时，边坡极限平衡状态所对应的折减系数，即为安全系数。

2 边坡地质条件

实例边坡位于湖南省怀化市麻阳苗族自治县某集镇。集镇分为侵蚀堆积河谷阶地地貌和剥蚀构造丘陵地貌两种地貌区。集镇沿锦江两岸的河谷阶地地貌区地段地势较为平坦，高程约为190～200 m，河谷地段宽度约为850～1 000 m。沿锦江两岸河谷地段之外为低山丘陵区，高程一般在230～260 m，山坡坡形呈弧形状，自然坡度45°～55°，山间沟谷长且宽缓。

图1 边坡地理位置图

该边坡属于人工开挖非均质边坡，边坡位于锦江的转弯处的右岸(图1)，边坡体纵长约135～136 m，高约41 m，坡脚高程约215 m，自然坡度约45°～50°，在坡脚有人工修建道路和开挖山体所形成的宽约51 m，高约15 m高陡边坡，坡度为85°。边坡地层分为4层，从上往下，第1层为新生界第4系全新统Qh粉质粘土，为紫红色，硬塑状，以黏粒及粉粒为主，含少量高岭土，稍有光泽，层理特征不明显，无摇震反应，干强度低，韧性差。该层曾发生土质滑坡，滑坡后缘有约15～20 mm拉张裂隙，滑体两侧有剪切裂隙。第2～第4层全为中生界白垩系上统K1L的泥质粉砂岩，软弱～极软薄厚层状，为紫红色，岩石硬度低，岩石抗风化能力差，节理极发育，倾角近乎水平。第5层为基岩。该边坡粉质黏土层曾发生土质滑坡，现基本稳定。具体建模所需地质参数见表1。

表1 边坡地层参数表Table1 Formationparametertableofslop地层岩性平均层厚/m弹性模量/Pa泊松比密度/(g·m-3)粘聚力/Pa内摩擦角/(°)1a3.01.00E+070.317002640019.92b10.01.00E+070.318502860023.93c13.01.00E+070.319333710029.54d15.02.00E+070.2820007500036.15e20.01.00E+080.25224300注:a为粉质黏土,b为全风化泥质粉砂岩,c为中-强风化泥质粉砂岩,d为强风化泥质粉砂岩,e为中风化泥质粉砂岩;地层参数是均处于天然状态下的。

3 ANSYS三维边坡模型建立

取人工开挖地段的宽度为模型边坡宽度，将模型简化成整体高61 m，宽41 m的几何实体，具体参数见图2，各地层厚度按表1设置，定义单元为solid45，采用ANSYS的Drucker-Prager模块，网格大小为1.5 m，网格扫掠(sweep)成六面体单元，生成40 018个单元(Element)，共有44 030个节点(NODE)，模型底面自由度全约束，侧面法向约束，表面自由，模型法向施加重力荷载。定义材料属性：所计算的模型为5层结构，1～4层假定为理想弹塑性材料，第5层(中风化泥质粉砂岩层)假定为理想弹性材料，其中对1～4层进行强度折减，折减系数分别为0.7、0.8、0.9、0.95、1，1.01、1.02、1.03、1.04、1.05、1.06、1.07、1.08、1.09、1.1、1.11、1.12、1.13、1.14、1.15。(对折减系数在1～1.15之间重点分析，得出边坡每一步的内在变形机制)根据式(1)～式(3)，对材料的初始抗剪强度参数和弹性参数进行折减，将折减的数值代入计算。求解器设置：求解打开自动步长，时间设为1，采用PCG求解器进行求解，求解子步数为20，打开Nonline框进行非线性搜索，在Analysis Options打开Large Displacement Static大变形。

图2 三维边坡模型图

4 计算结果分析

4.1 传统判据

4.1.1计算收敛性判据

从表2可以看出，当折减系数达到1.1时，计算仍收敛，边坡处于平衡状态，当折减系数达到1.11时，计算不收敛，边坡已经失稳破坏，此时的安全系数介于1.1～1.11。

表2 ANSYS计算收敛性结果表Table2 ANSYSConvergenceResultsTable折减系数0.70.80.911.011.021.031.041.051.06收敛情况收敛收敛收敛收敛收敛收敛收敛收敛收敛收敛折减系数1.071.081.091.11.111.121.131.141.15收敛情况收敛收敛收敛收敛不收敛不收敛不收敛不收敛不收敛

4.1.2塑性区贯通性判据

在岩土工程中，常用等效塑性应变来描述岩土体的剪应变。根据等效塑性应变云图(见图3)，由于边坡曾经失稳，在第一层粉质黏土中发生过土质滑坡。当折减系数为0.7～0.8时，岩土体材料参数被加强，一定程度上“还原”土质滑坡滑移状态，最大塑性应变出现在坡面上，坡脚出现塑性应变区，滑坡后缘位置与现场情况符合；折减系数为0.9～1.0时，土质滑坡已经滑移，现基本稳定，与现场情况符合，最大塑性应变区转移至坡脚上端1.5 m处；折减系数为1.02时等效塑性应变区向坡体内部延伸，因坡体内部为泥质粉砂岩，延伸区成线性条带状，折减系数为1.04～1.06时，坡体内部出现塑形应变区，与坡脚向坡体内部延伸的塑形区连接，此时计算收敛，边坡处于蠕变阶段；折减系数为1.08～1.10时，塑形区向坡体延伸，形成贯通，贯通面尾部未到自由表面，此时计算收敛，坡体处于蠕变阶段，折减系数为1.11～1.13时，塑形区形成贯通，贯通面到达自由表面，计算不收敛，坡体处于破坏变形阶段，安全系数介于1.1～1.11，为岩质滑坡，滑裂面为如图3所示。结合等效塑性应变矢量图(见图4)，边坡先是发生土质滑坡，随着折减系数的增大然后发生岩质滑坡，岩质滑坡体方量约15 000～16 000 m3，传统判据所得安全系数介于1.1～1.11。

(a)Ki=0.7

(i)Ki=1.1

4.2 新式突变判据

4.2.1位移指标判据

从图5中可以看出，4种位移都存在突变与屈服阶段，3种曲线均有明显的极值点，其中Z方向位移非常小，变化幅度小，可忽略不计。空间位移、X方向位移、Y方向位移曲线趋势相近，其中空间位移与Y方向位移曲线数值基本相等，当折减系数K增大到极值点之前均发生位移突变，继续增大折减系数，位移曲线均发生了屈服，因此可以定性地以这3种位移方式下的位移曲线作为研究对象，采用最小二乘法拟合来定量确定极值点处的安全系数。假定位移与折减系数关系满足Gauss、Lorentz、Vigot曲线方程，采用最小二乘法对该曲线进行拟合，Gauss曲线方程为：

(4)

Lorentz曲线方程为：

(5)

Vigot曲线方程为：

Vigot曲线为Gauss曲线型卷积Lorentz曲线形型，与Gauss、Lorentz具有相同的单调性，所以当Ki=Kic时，此时的U→极大值，即认为边坡处于极限平衡状态。采用最小二乘法拟合，拟合结果的可靠度用拟合相关系数表示，相关系数越接近1，残差平方和越接近0，拟合效果越好。

图6～图8给出了不同位移方式下曲线拟合情况，并根据式Ki=Kic求解对应的安全系数，拟合结果和安全系数列于表3中。分析表中的拟合结果发现：Gauss、Lorentz曲线对X方向位移指标拟合不收敛，拟合效果差，Vigot曲线对3种位移方式拟合效果都较好，特别是空间总位移，3种位移方式下得到安全系数介于1.123 6～1.133 9，相差较小。根据最佳拟合结果得出的安全系数为1.123 8，将安全系数1.123 8代入ANSYS计算所得结果不收敛，边坡发生破坏。从Ki=0.7到Ki=1.123 8的X方向位移云图可以看出，大位移变形区由坡体内部过渡到坡体表面，在坡体内部形成，明显滑移面，位置与等效塑性应变云图(如图3)贯通区吻合。

图5 位移与折减系数的关系曲线

图6 空间总位移曲线拟合

图7 X方向位移曲线拟合

图8 Y方向位移曲线拟合

(a)Ki=0.7

表3 位移曲线拟合结果Table3 Curvefitresultofdisplacement拟合对象拟合方式拟合参数U0KicAw相关系数残差平方和安全系数Gauss拟合3.91321.12450.09060.06970.84900.03691.1245空间位移Lorentz拟合3.79831.12360.15300.07950.90710.02271.1236Vigot拟合3.79071.12380.15806.3184E-9[1]0.0826[2]0.90700.02431.1238Gauss拟合不收敛X方向位移Lorentz拟合不收敛Vigot拟合0.51941.13390.21391.5855E-06[1]0.0970[2]0.92410.02771.1339Gauss拟合3.91491.12440.08980.06900.84890.03691.1244Y方向位移Lorentz拟合3.79871.12360.15270.07930.90700.02271.1236Vigot拟合3.67981.12870.23461.7480E-06[1]0.1258[2]0.88960.02891.1287注:[1]代表系数WG;[2]代表系数WL

4.2.2塑性区指标判据

从图11中可以看出，9种塑性区指标都存在突变阶段与屈服阶段，由于9种曲线均有明显的极值点，其中Y、Z方向塑性应变、XZ、YZ剪切应变非常小，变化幅度小，可以忽略不计。对其余5种塑性指标与折减系数关系的分析依然采用Gauss、Lorentz、Vigot曲线拟合，图12～图16给出了不同塑性区指标下曲线拟合情况，并根据式Ki=Kic求解对应安全系数，拟合结果和安全系数分析见表4，分析表中的拟合结果发现：Vigot曲线对5种塑性指标都拟合，且拟合效果最优，Gauss、Lorentz曲线只对塑性应变强度不拟合，其余4种拟合效果较好，5种塑性指标拟合得到的安全系数介于1.127 8～1.129 6，相差非常小，可以忽略不计，根据拟合效果，将X方向塑性应变拟合所得结果1.129 6作为安全系数。由X方向塑性应变云图可以看出，随着折减系数增大，塑性区发生贯通。

图11 9种塑性区指标与折减系数关系曲线

4.4 失稳判据的对比分析

通过对比传统判据和新式判据的安全系数，可以发现：

a.传统的计算收敛性与塑性区贯通性判据：2种判据原理简单，而且在ANSYA有限元计算过程中基本同步，折减系数都介于1～1.1。其值明显比新式判据所得安全系数小，这是因为计算结果与网格划分、收敛准则的参数设置有关，具有一定的人为自由性，计算收敛性与模型参数设置也有关，如果模型参数取值不当，也可能导致结果提前不收敛，则所得安全系数偏小。对于理想弹塑性单元，当其应力达到屈服状态后，该单元周围依然存在弹塑性约束，那么这些约束会限制这个单元的塑性应变发展，坡体依然是稳定状态，只有继续增大折减系数，坡体内单元的应变逐渐增大，塑性单元逐渐增多，才形成整体滑坡，若仅以贯通作为判断标准，则所得安全系数偏小。

图12 塑性应变强度曲线拟合

图13 Von Mises塑性应变曲线拟合

图14 等效塑性应变曲线拟合

图15 X方向塑性应变

图16 XY剪切塑性应变

b.位移指标突变与塑性指标突变判据得到的安全系数十分接近，两者都可以通过最小二乘法对指标与折减系数关系进行非线性拟合，且假定它们之间的关系满足Gauss、Lorentz、Vigot曲线方程，从而确定边坡位移指标与塑性指标与折减系数的定量关系，能够表征边坡在折减系数达到一定程度后产生多种位移与塑性应以应变的极大值，然后出现屈服阶段，用来确定安全系数具有明确的物理意义，建议将两者结合分析。

(a)Ki=1

表4 塑性指标曲线拟合结果表Table4 Tableofplasticitycurvefittingresults拟合对象拟合方式拟合参数U0KicAw相关系数残差平方和安全系数Gauss拟合不收敛塑性应变强度Lorentz拟合不收敛Vigot拟合0.04361.12930.02510.0430[1]0.0085[2]0.91260.00331.1293Gauss拟合0.02921.12930.01320.04040.91190.00101.1293vonMises塑性应变Lorentz拟合0.01171.12920.02150.04840.89820.00121.1292Vigot拟合0.02581.12930.01470.0423[1]0.0100[1]0.91260.00111.1293Gauss拟合0.02921.12930.01330.04030.91180.00111.1293等效塑性应变Lorentz拟合0.01161.12920.02170.04830.89800.00121.1292Vigot拟合0.02591.12930.01470.0423[1]0.0097[2]0.91250.00111.1293Gauss拟合0.01801.12960.01120.04070.91570.00071.1296X方向塑性应变Lorentz拟合0.00341.12960.01820.04880.89810.00081.1296Vigot拟合0.01661.12960.01180.0452[1]0.0052[2]0.91590.00081.1296Gauss拟合0.02441.12790.00670.03910.89770.00031.1279XY剪切塑性应变Lorentz拟合0.01541.12780.01090.04630.89690.00031.1278Vigot拟合0.01991.12780.00860.0322[1]0.0237[2]0.90200.00031.1278注:[1]代表系数WG;[2]代表系数WL

c.采用传统极限平衡法进行计算验证，其中摩根斯顿-普莱斯法所得安全系数为1.132，Janbu法所得安全系数为1.119，安全系数十分接近，滑面位置如图18所示，与新式突变判据所得滑面位置吻合。

5 结论

a.采用有限元强度折减法，以等效塑性应变为基础，在一定程度上模拟了边坡从土质滑坡到岩质滑坡等效塑性应变全过程；采用传统判据分析三维边坡稳定性，由计算收敛性所得安全系数为1.1～1.11，由塑性区贯通性判据所得安全系数为1.1～1.11，两者结果偏于保守。

(a)摩根斯顿-普莱斯法

b.采用新式突变判据，3种位移值和5种塑性应力应变值都存在突变与屈服阶段，采用最小二乘法对边坡位移值和塑性应力应变值与折减系数关系进行非线性拟合，假定它们之间的关系满足Gauss、Lorentz、Vigot曲线方程，效果达到预期，其中Vigot曲线拟合效果最好，且应用到此边坡不存在不收敛情况，因此Vigot曲线方程拟合为该新式突变判据理想的拟合方式。

c.新式突变判据所得安全系数分别为：位移指标突变判据的安全系数取自空间总位移曲线Vigot拟合结果，为1.123 8；塑性指标突变判据的安全系数取自X方向塑性应变曲线Vigot拟合结果为1.129 6。所得安全系数十分接近，且均达到极限平衡状态采用极限平衡法对比验证，结果可行，因此建议将两者结合分析。