基于Phase2的软基边坡有限元强度折减法分析

王海建

(中水珠江规划勘测设计有限公司，广东 广州 510610)

0 引言

对于边坡稳定分析方法，极限平衡法研究历史最为悠久，至今已成为最为通用和被认可的方法，如其中的瑞典圆弧法、毕肖普法等。但极限平衡法很难获得边坡尤其是支护结构的应力应变情况，并需要事先假定屈服面，而有限元方法优点有①不需要进行滑动面假设，且可以得到应力应变分布场；②能够模拟边坡的失稳过程及塑性变形区分布情况；③可以考虑各种支挡结构与岩土材料的共同作用，为边坡支护设计提供科学依据。20世纪70年代Zienkiewize等[1-3]就利用有限元方法进行边坡稳定分析，但是受计算机条件的限制，此法没有得到广泛推广。90年代以来，随着计算机算力的飞速发展，以及岩土材料本构模型、屈服准则等的深入研究，有限元强度折减法在边坡稳定分析中的应用越来越广泛[4-5]。但对于软土基础边坡的研究和应用较少，如郑颖人、赵尚毅等[5-7]的研究，大多基于岩石边坡或均质土坡，其潜在滑裂面一般在坡脚区域，右边界(坡顶方向)范围R的影响较左边界L(坡脚方向)的影响更大[7-11]，且认为网格疏密对计算精度的影响甚至大于单元类型的影响，每10平方米不少于3个节点时，计算精度较为理想。刘金龙等[12]的研究认为网格尺寸在0.5～1m比较合适。但对于软基边坡，滑裂面一般穿透软土层，同时应力应变分布与岩石或均质边坡不同，其精度影响参数与均质边坡差异较大，需对此进一步深入的研究，以更好地分析边坡稳定和采取合适的支护措施。

1 基本原理

有限元强度折减法，就是在有限元计算边坡稳定时，将边坡岩土体抗剪切强度参数c、φ值等逐渐按比例折减降低，直到其达到破坏状态为止，程序可以根据弹塑性计算结果得到破坏滑动面(塑性应变和位移突变的地带)，此时的折减系数即为边坡的稳定安全系数。

在进行边坡稳定分析时，岩土材料模型一般假定为弹塑性材料，常用的Mohr-Coulomb强度准则(简称M-C准则)表达式为：

τn=c+σntanφ

(1)

式中，c—土的黏聚力；φ—土的内摩擦角；σn、τn—滑移面上的正应力与切应力。

假设强度折减系数为SRF，对于(1)式有：

(2)

则(2)式变为：

(3)

式中，SRF—折减系数，也就是边坡稳定安全系数；c′、φ′—折减后的黏聚力和内摩擦角。

有限元强度折减法初始条件SRF=1，对于稳定的边坡，程序收敛，此时逐渐增大SRF，直至程序不收敛；若初始边坡处于不稳定状态，此时减小SRF，即增大c′、φ′值去进行试算，直到边坡位移变形出现不收敛的拐点，同时塑性区刚好处于贯通与未贯通这样的一个分界点，此时对应的折减系数，即为求得的边坡稳定系数。

用折减系数法求解实际边坡稳定，目前常用的有限元计算软件如ANSYS、MARC、MASTRAN等，所用的岩土材料屈服准则均采用广义Mises准则(又称D-P准则)，表达式为：

(4)

式中，αφ、k—与材料性质相关的参数；J2—应力偏量第二不变量。

为便于有限元计算，将M-C准则与D-P准则表达统一，即对于三轴应力状态下土体，其剪切破坏的强度准则表达为：

(5)

与D-P准则一致考虑应力不变量时，M-C准则下的屈服面表达式为：

(6)

式中，θ—应力罗德角；I1—应力张量第一不变量。

用有限元强度折减法求边坡稳定安全系数，国内郑颖人、赵尚毅等[5-7]自20世纪90年代开始做了大量的研究，包括屈服准则、流动法则、有限元模型本身以及参数选取等对计算精度的影响，重点研究本构模型和屈服准则的选择，如M-C屈服面外接圆、内接圆、等面积圆等以及不同屈服准则的转换，其中采用等面积圆计算结果与M-C准则计算结果最为接近，该准则求得的稳定安全系数与传统Spencer法的误差在5%左右，且与M-C屈服准则计算的结果十分接近。

材料参数的差异也是影响计算精度的主要因素之一。由于软土层具有高压缩性、抗剪强度低、透水性低、触变性、流变性、不均匀性等特点，软基边坡发生破坏时，滑动面通常都位于软土层与相对持力层交界处，导致采用有限元分析软土边坡时，网格划分和边界条件需要专门研究。

本文采用Rocscience公司岩土有限元软件Phase2，其内置了强大的本构模型和屈服准则，来研究软基边坡分析中的几个关键问题。

2 失稳判据

有限元强度折减法分析边坡稳定性的一个关键问题是如何根据有限元计算结果来判别边坡是否处于破坏状态[6-8]。目前的失稳判据主要有：

(1)采用力和位移的不收敛作为边坡失稳的标志。

(2)以广义塑性应变或者等效塑性应变从坡脚到坡顶贯通作为边坡破坏的标志[9-12]。

(3)边坡特征点的位移发生突变。当岩土材料强度不断降低时，边坡内部位移场中特征点的位移也会随着发生变化；当强度降低至临界点时，位移变化的斜率由平滑曲线发生突变，此时可以认为边坡已经失稳破坏。

根据郑颖人等[10-14]的研究表明，以上判据得到的分析结果通常相差不大，但在某些特殊土层和边界条件下，如果采用单一判据，可能会出现较大的差异，所以一般采用两种或两种以上来作为边坡破坏的标志。后续案例分析可以较好地验证这一点。

3 模型网格划分和边界

3.1 网格划分

Phase2中有限元网格划分有3节点三角形、4节点四边形、6节点三角形和8节点四边形等几种。一般来说节点越多，单元划分越密，计算结果越精准。但也会大大增加计算时间，对于复杂边坡也可能导致迭代次数不够或出错。所以需要研究网格划分的影响，在保证计算精度的同时也有利于加快计算。

张鲁渝、郑颖人等[7-11]在均质边坡中研究认为，网格疏密对计算精度的影响甚至大于单元类型的影响，每10平方米不少于3个节点时，计算精度较为理想。刘金龙等[12]的研究认为网格尺寸在0.5～1m比较合适。为便于讨论，这里选用刘金龙等人所采用的算例边坡作为研究对象，如图1所示。坡高H=20m，坡角45°，土的重度γ=20kN/m3。土体采用M-C准则的理想弹塑性本构模型，土的粘聚力c=42kPa，内摩擦角φ=17°。弹性模量E=5MPa，泊松比v=0.3。

图1 有限元计算模型

不同网格划分有限元计算结果与常用的极限平衡法(毕肖普法、Spencer)进行对比，见表1。

表1 不同网格单元划分计算结果对比表

从计算结果可以得出以下结论：

(1)总体上，网格疏密的要求本质上是对于节点的要求。节点越多、网格单元越密，计算结果越精准，同时计算用时越长。

(2)节点数对于计算精度的影响比单元数更为密切，如单元数1000左右的8节点四边形单元计算结果与单元数1600左右的6节点三角形单元计算精度已经与极限平衡法近似，其节点数也类似，分别为3036和3298。

(3)根据本模型的尺寸，算出节点间距在0.5m左右时(本例中节点数3000左右)，计算结果已经满足工程要求，与极限平衡法偏差在1%左右。

(4)节点数类似时，8节点四边形单元和6节点三角形单元计算结果明显更精准。

因此，为便于工程应用，建议有限元模型计算时，单元划分节点间距控制在0.5～1m，采用8节点四边形单元或6节点三角形单元。再加密节点对于计算结果影响已经不大。

3.2 模型边界

模型边界范围的大小在有限元法中对计算结果的影响比在传统极限平衡法中表现的更为敏感。其中传统极限平衡法要求滑移面在边界之内即可，但有限元法中边界的范围直接影响到应力-应变的分布[10-12]，从而影响稳定计算结果。如图2所示边坡高度为H示例，张鲁渝、郑颖人等[10]针对均质边坡的研究表明，右边界(坡顶方向)范围R的影响较左边界L(坡脚方向)的影响更大，当L>1.5H，B>H，R>2.5H时，边界再扩大对于计算结果的影响小于1%。

其中，关于底边界，该结论未考虑到软基的情况，在基础中有软土层分布时，滑动面往往从软土层底面滑出，模型下边界须在软基底面以下，并深入相对坚实土层一定深度。

在此基础上，分析模型边界的影响。采取边坡模型如图3所示，软土层重度γ=17kN/m3，粘聚力c=5kPa，内摩擦角φ=10°，弹性模量E=3MPa，泊松比v=0.35。

图3 软基边界范围示意(修正后H′)

关于模型边界约束条件，底部采用完全固定边界，对于两侧的约束条件，建议根据实际判断是否会产生位移，相应采用固定边界或水平约束条件。

为了得到能使计算结果趋于稳定的边界范围，分别对左端、右端、底端3个边界范围的取值大小进行了分析。计算时，令3个边距中的1个变化，其余2个不变，初始条件相对边距比值均为1。计算结果见表2，如图4所示。

表2 不同边距比计算结果对比表

图4 不同边距比计算结果对比

由表2结果可以看出，由于软土层的压缩模量低、变形较大，单元应力应变分布不同，使计算结果的精度偏差较大，但基本在5%以内。

经过对比分析，塑性破坏区基本位于软土层，基于此，认为可以近似把边坡高度中的H算至软土层的底边，模型各边界范围边距比均以修正后的H′来调整。上述例子中根据该方式修正后，即L>1.5H′，B>1.5H′，R>2.5H′时，边界扩大对于计算结果的影响变动在1%以内。修正H后的计算结果如图5所示。

图5 按修正后的H′计算结果

需要指出的是，与张鲁渝等[10]人研究结论有所不同的是，对于软土基础边坡，底边界范围B对于计算结果的影响比右边界R更为敏感，这是需要注意的地方。

4 案例分析

我国东南沿海地区尤其是珠三角等广泛分布有软土地层，边坡稳定分析和支护设计是关键点和难点。下面以某一堤岸边坡为例，采用有限元强度折减法进行分析。该堤岸边坡高约10m，下部分布有软弱淤泥质土层厚度约8m，本次计算工况为低水位时的边坡稳定。

4.1 模型建立

根据前述的分析，建立外部模型边界范围。模型网格主要采用8节点四边形单元。模型网格共969个单元，3024个节点，如图6所示。

图6 计算案例有限元模型

4.2 岩土参数

该段堤防地层分布主要有①素填土层/新填土、②-1(含泥)砂层、②-2粉质粘土、②-3淤泥质土、③-1残积粉质粘土。其中②-3淤泥质土呈流塑状，含水量高、孔隙比大、抗剪强度低，是影响岸坡稳定的主要软土层。岩土参数见表3。

表3 岩土参数表

4.3 结果分析

4.3.1计算结果

采用有限元强度折减法计算该边坡最大剪切应变云图如图7所示。可以看出其最大剪切应变发生在软土层底部与相对坚实土层的交界面区域。通过云图可以看出，在折减系数为0.992时，塑性破坏区即将贯穿边坡。

图7 最大剪切应变云图

最大位移云图如图8所示。折减系数与最大位移关系曲线如图9所示。通过折减系数与最大位移关系曲线可以看出，在折减系数为0.982时，边坡最大位移0.295m，此时最大位移产生突变，但程序仍可以收敛，塑性区未贯通；当折减系数为0.992时，边坡最大位移0.415m，此时继续增加折减系数，程序开始不收敛，且塑性区贯通。因此可以认为，0.992即为折减系数的临界值，且同时满足了边坡失稳的3个判据。

图8 最大位移云图

图9 折减系数与最大位移关系曲线

4.3.2计算结果与极限平衡法对比

同一个模型采用极限平衡法计算，得出稳定安全系数为0.998～1.01。不同计算方法得到的稳定安全系数对比见表4，同时将极限平衡法计算求得的圆弧滑面与有限元强度折减法计算结果进行叠加，如图10所示。通过对比可以看出，该滑面与有限元强度折减法计算所得的最大剪切应变区基本吻合，同样位于软土层底部。

表4 稳定安全系数对比表

图10 有限元强度折减法与极限平衡法滑动面(红色虚线)

4.3.3支护后计算结果

从稳定计算结果可知，边坡基本处于临界状态，拟采用在坡脚增设支护桩增加其稳定性[15]。支护桩采用刚性桩，连续密排布置，在Phase2中采用混凝土梁单元模型，厚度采用0.2m，泊松比0.2，杨氏模量取3×107kPa。桩底按穿透软土层1m(方案一)、2m(方案二)两种情况分别进行计算。

最大剪切应变分布如图11—12所示。可以看出增设支护桩后，稳定系数增加，桩越长，稳定系数越大。其中方案一桩底下部仍有剪切应变区，桩身有脱离土层失效的风险，临界折减系数1.112时的最大位移0.952m。方案二桩底嵌入深延长后，桩底以下基本无剪切应变分布，变形情况明显改善。临界折减系数1.128时的最大位移0.881m。

图11 最大剪切应变分布云图(桩底穿透1m)

图12 最大剪切应变分布云图(桩底穿透2m)

桩身弯矩分布如图13—14所示。可以看出不同强度折减系数(SRF)下桩身弯矩分布情况，与极限平衡法得出的桩身弯矩分布情况基本一致。对于方案一桩底嵌入持力层深度偏小，最大弯矩582kNm，方案二延长桩底进入持力层长度后，桩身最大弯矩525kNm。桩底锚入持力层加深后发挥锚固作用更加明显。可以推断桩长继续增加则稳定系数增加，但造价也相应增加较多。从结果来看，采用方案二桩底进入持力层2m已满足稳定需要和桩身承载能力，不需再继续加长。

图13 桩身弯矩分布图(桩底穿透1m)

图14 桩身弯矩分布图(桩底穿透2m)

5 结论

(1)有限元强度折减法中边坡失稳的3个判据，在实际应用时宜采用两种或两种以上综合确定，可以使计算结果更为精确。

(2)网格划分疏密的要求本质上是对于节点的要求，宜达到节点间距0.5～1m较为合适。软土基础边坡破坏时滑动面基本沿软土层底部边界，模型三个边界条件中，坡高H宜采用软土层高度修正后的H′来调整，当L>1.5H′，B>1.5H′，R>2.5H′，计算精度满足工程需要，其中底边界范围B对于计算结果的影响更为敏感。但是当软土层过于深厚时，边界范围应根据试算结果进一步调整。

(3)采用上述网格划分和边界条件，有限元强度折减法计算得出的结果与极限平衡法计算得出的稳定系数和滑动面均较为一致，结果偏差也在1%以内。

(4)有限元强度折减法计算可以得出软土边坡应力应变分布，为边坡支护提供依据，可以据此选出合适支护桩桩型。本文支护桩采用的刚性桩为弹性材料模型，采用柔性桩支护时或者边坡破坏时桩土作用方式和机理还有待进一步研究。