吴伯建,朱珍德,顾祖军

(1.河海大学 岩土工程科学研究所,江苏 南京210098;2.河海大学 岩土力学与堤坝工程教育部重点实验室,江苏 南京 210098;3.通州建总集团有限公司,江苏 南通 226300)

0 引 言

边坡广泛存在于各种工程建设中,对其稳定性研究历来是岩土工程领域的一个研究热点。采用有限元分析边坡稳定性满足了土体本身的应力应变关系,且可以处理复杂的边界条件和材料的各向异性,但不同于极限平衡法,由于无法直接获取安全系数指标,有限元分析在实际工程应用时存在局限性,强度折减法的引入改变了这一格局,国内众多学者对两者结合的有限元强度折减法展开了大量研究。郑颖人等[1]详细分析了有限元强度折减法的计算精度和影响因素。王栋等[2]通过开发具有自动搜索安全系数功能的边坡稳定分析有限元模型探讨了单元阶次对安全系数的影响。唐芬等[3]对强度指标采用了不同的折减系数进行了定量计算。杨光华等[4]基于Duncan-Chang非线性本构模型提出在弹性阶段对弹性模量也进行折减的变模量弹塑性模型强度折减法,获得了更符合实际的变形场。

本文采用通用有限元软件ABAQUS结合强度折减法对土坡进行分析研究,以了解由于有限元软件

本身网格划分和土坡边界范围选取对计算结果精度造成的影响,最后利用得到的结论对实际土坡进行数值分析,将所得结果同土坡实际滑动情况及极限平衡法求得的安全系数展开比较。

1 有限元强度折减法

1.1 原理

强度折减理论最初出现于1975年,Zienkiewize在土工弹塑性有限元数值分析中将其提出,运用该理论时对抗剪强度参数不断折减,直至达到极限状态即实际发挥的抗剪强度与外荷载产生的实际剪应力相等,此时的折减系数就是安全系数Fs。

有限元软件ABAQUS中实现强度折减法首先定义一个场变量FV,取为强度折减系数Fr,然后对土体强度指标粘聚力 c和内摩擦角φ按下式进行调整[5]:

开始分析时Fr取值小于1,以放大强度,避免对模型施加重力荷载时发生破坏,分析过程中线性增加Fr,直至计算中止。

1.2 失稳判据

目前存在的土坡失稳判据主要有三种:有限元数值计算收敛与否(判据1)、特征部位位移的拐点(判据2)和塑性区是否形成连续的贯通区(判据3)。

上述三种判据中,判据1的影响因素众多,单元类型、网格密度、迭代次数、收敛限定值等都会导致数值计算收敛发生变化,判据2和判据3的物理概念明确,但同样存在上述计算精度方面的问题。

一些研究人员对边坡失稳判据的选取进行了研究:裴利剑等[6]认为三类失稳判据在理论上具有统一性和一致性,判据1最为方便可靠;刘金龙等[7]建议联合采用判据2、3作为失稳判据,李红等[8]认为判据1、2具有普遍适用性。本文进行有限元计算时分别采用三种判据来获得安全系数值,对三种判据进行评价。

2 有限元精度分析

采用一均质土坡[9]作为分析算例,该算例被众多学者采用进行了研究。土坡参数见表1,初始选择的坡体几何尺寸见图1。

采用传统极限平衡条分法计算得到的安全系数Fs值见表2,在对有限元强度折减法获取的安全系数进行分析评价时,以表2中Morgenstern-price法的Fs值作为标准。

表1 土坡参数

图1 土坡几何尺寸(单位:m)

表2 不同条分法所得Fs值

与传统极限平衡法相比,分析时还需采用的参数包括弹性模量E、泊松比 ν、剪胀角 ψ,研究[10]表明三者对边坡分析结果均有影响,本文不再进行探讨,统一取 E=100 MPa,泊松比 ν=0.35,ψ=0。

2.1 网格划分

划分网格是建立有限元模型的重要环节之一,会对计算精度和计算规模产生显著的影响。ABAQUS划分网格由三部分内容组成:确定单元类型、选择划分技术及布设种子以控制网格密度。分析时网格单元类型分别采用一阶三角形单元CPE3、二阶三角形单元CPE6、一阶四边形单元CPE4、二阶四边形单元CPE8;网格划分技术选择最为灵活的自由网格划分选项;布设种子尺寸控制为全局尺寸1,在CPE8的基础上对坡体局部加密,将坡面及坡顶区域单元尺寸设为0.5(CPE8′),按上述设定进行数值模拟,采用三种判据得到的安全系数值列表如表3。

表3 不同单元类型所得Fs值对比

为了直观了解不同网格划分得到的Fs值与Morgenstern-price法Fs值的偏差,作相对误差柱形图:

图2 Fs值相对误差柱形图

结合图2、表3进行分析:二阶单元的精度明显优于一阶单元,四边形单元的精度优于三角形单元,在同等网格密度划分时,选用二阶四边形单元效果最好;在不同网格密度情况下,网格过于稀疏会导致误差值增大,本例中通过对局部重点位置的加密得到了较好的结果,图2中网格局部加密的误差值是最小的,仅为0.93%。在对边坡进行网格划分时,建议优先采用二阶四边形单元,在坡面及坡顶滑动面经过的区域进行网格加密措施,以获取较精确的结果。

对比三种判据得到的Fs值,判据1得到的结果波动较大,有限元计算的收敛性与坡体的极限平衡状态没有表现出直接对应关系,若将其作为唯一的失稳判据是不合理的;而坡体失稳时,特征点位移突变和塑性区贯通则几乎是同时发生的,5组数据中有3组的结果是一致的,剩余2组中CPE3单元误差较大,选用CPE8单元数据分析,判据2、3对应的FV值分别为 0.9828、1.0145,当 FV=0.9828时,塑性区未贯通,见图3(a),FV=1.0145塑性区贯通,见图3(b),而实际发生塑性区贯通的临界点是介于0.9828～1.0145之间的,当精度足够时位移突变点也会发生在该区间内,这样的分析结果同Morgensternprice法Fs值是吻合的。建议采用判据2、3结合使用进行评价。

图3 塑性区云图

2.2 边界范围

传统极限平衡条分法在计算安全系数时,边界范围只需将破坏面包括在内,计算结果就不会有影响,但有限元分析中边界范围的大小会造成计算结果发生变化,为了了解土坡受边界范围变化所产生的计算偏差,对不同边界范围的土坡展开分析。

定义L1为坡脚到左端边界的距离,L2为坡顶到右端边界的距离,L3为上下边界总高减去坡高后的距离。土坡的边界范围主要受这三个参数控制,在分析网格精度中采用的边界范围已经将坡体的破坏区包括在内,在此基础上将研究对象参数设为变量值,保持其他两参数值不变,对研究的对象参数进行放大,即扩大边界范围,得到的安全系数见图4。

图4 Fs值随单一边界参数变量变化的折线图

分析中选取的特征点为坡脚点A和坡顶点B,位置见图1,将对应FV=0.9828时的点A处水平位移U1值及点B处竖向位移U2值列表进行误差分析见表4。

表4 位移相对误差对比

通过图4可知任意一个边界参数的变化都会影响安全系数Fs的值,边界扩大到原来两倍时,精度值有所上升,但扩大至3倍精度反倒减低,说明一味的追求边界范围扩大并不可取。在L1、L2、L3中,边界L1对安全系数的影响相对较小。针对本文实例,安全系数的变化范围仅在0.004内,边界范围的扩大对分析结果有影响,但影响甚微,建议在一般分析中,要求边界范围稍大于坡体破坏区即可,不应盲目扩大边界范围,由于追求外界边界所产生的计算量增加及收敛问题,有时往往得不偿失。

表4反应了特征点 A、B位移对边界范围的敏感性。随着L1值的不断扩大,点A水平位移和点B竖向位移均增加,点 A水平位移受L1影响更加明显;L2对两特征点的位移影响均衡;L3主要影响了点B的竖向位移,竖向位移增加趋势明显。由于土坡滑动方向向左,左边界L1主要控制水平位移,下边界L3主要控制竖向位移这一结论与实际也较符合。在计算时适当的扩大L1和L3可以获得更大的精度。

由于网格划分时采用了二阶四边形单元,在此基础上进行了对边界范围的影响分析,可见得到的计算结果偏差不大,在有限元自身影响因素中,网格划分的优劣比边界范围的大小重要的多。

3 实例分析

为了充分说明所得结论的实际应用效果,采用一实际发生滑动的土坡进行分析,该土坡在强降雨条件下发生了滑动,滑动面的大致形态由工程地质测绘推测得到,位置如图5所示。坡体内土体主要为粉质粘土,局部有孤石分布,分析时视为均质土坡,土体参数为:γ=20.0 kN/m3,c=12.5 kPa,φ=14.5°。Morgenstern-price法计算得到的Fs=1.276,滑动面位置见图5。

图5 土坡滑动面位置

进行有限元分析前,选取的边界范围如图5中坡体边界所示,左右边界为150 m,上下边界为50 m,L1=14 m,L2=30 m,L3=16 m。网格采用二阶四边形单元,除了坡面位置的网格尺寸加密为2,其余三面设为4。计算得到的场变量FV随特征点A处U1及特征点 B处U2的变化曲线图见图6,特征点位置见图5,图6中数值计算不收敛对应的 FV=1.280,位移突变处对应的FV=1.260,土坡塑性区等值线图见图7,在FV=1.265时塑性区发生贯通现象。故三种判据得到的的 Fs值分别为1.280、1.260、1.265。

观察图7可以看出,在场变量FV=0.926时,坡体前缘部位塑性区贯通,坡体前缘部位现场实际形态是由于坡脚处修筑临时道路形成的,在强降雨前已有位移。在强降雨影响下,整个坡体内的土体强度下降,造成了土坡形成图5中的裂缝及滑动面。

图6 FV随位移的变化曲线

图7 塑性区随FV的变化过程

对比图5及图7(d),条分法和有限元强度折减法求得的滑动面位置基本一致,同实际推测的滑动面比较仅滑动面后端延伸有差别,这是由于现场实际条件造成的,在此处存在孤石,影响了滑动面的延伸。有限元强度折减法在分析滑动面形成过程优势明显,直观的反映了坡体滑动面形成的全过程,且未遗漏对土坡的前缘部位的分析。

由三类判据得到的安全系数同条分法的误差控制在1.3%内,达到了良好精度。

4 结 论

(1)网格划分即对单元类型、网格疏密的设定影响了计算精度。同等条件下,二阶单元的精度优于一阶单元,四边形单元的精度优于三角形单元,对坡体滑动部位的网格进行局部加密取得的效果最为良好,既使精度得到了一定提高,也避免了计算量不合理的增加。实际应用时,建议采用二阶四边形单元,局部加密网格对研究对象进行分析。

(2)由三种判据得到的安全系数可知:收敛性判据的波动性太大,建议实际中结合位移和塑性区判据进行评价。

(3)当采用二阶四边形单元时,土坡边界范围的选择对分析结果的影响甚微,建议边界范围稍大于坡体破坏区即可。边界范围中边界L1主要影响坡体的水平位移,边界L3主要影响坡体的竖向位移。

(4)利用本文结论对一实际土坡进行数值模拟,得到的结果与条分法的结果接近,滑动面位置与实际滑动位置吻合,达到了用于工程实际分析的精度,为类似均质土坡稳定性分析提供参考。

[1]郑颖人,赵尚毅.有限元强度折减法在土坡与岩坡中的应用[J].岩石力学与工程学报,2004,23(19):3381-3388.

[2]王 栋,年廷凯,陈煜淼.边坡稳定有限元分析中的三个问题[J].岩土力学,2007,28(11):2309-2313,2318.

[3]唐 芬,郑颖人.边坡稳定安全储备的双折减系数推导[J].重庆交通大学学报(自然科学版),2007,26(4):95-100.

[4]杨光华,张玉成,张有祥.变模量弹塑性强度折减法及其在边坡稳定分析中的应用[J].岩石力学与工程学报,2009,28(7):1506-1512.

[5]费 康,张建伟.ABAQUS在岩土工程中的应用[M].北京:中国水利水电出版社,2010.

[6]裴利剑,屈本宁,钱闪光.有限元强度折减法边坡失稳判据的统一性[J].岩土力学,2010,31(10):3337-3341.

[7]刘金龙,栾茂田,赵少飞,等.关于强度折减有限元方法中边坡失稳判据的讨论[J].岩土力学,2005,26(8):1345-1348.

[8]李 红,宫必宁,陈 琰.有限元强度折减法边坡失稳判据[J].水利与建筑工程学报,2007,5(1):79-82.

[9]Dawson E M,Koth W H,Drescher A.Slope stability analysis by strength reduction[J].Gé otechnique,1999,49(6):835-840.

[10]凌平平,杨雪强,向胜华,等.弹性模量、泊松比和剪胀角对边坡安全系数的影响[J].岩土工程界,2007,10(6):57-59.