基于PLAXIS和OPTUM的路堤加筋数值模拟*

左建忠 蒋 鑫 付用国 陈欣妮 黄岑嶷 邱延峻

(1.西南交通大学土木工程学院 成都 610031； 2.西南交通大学道路工程四川省重点实验室 成都 610031；3.西南交通大学高速铁路线路工程教育部重点实验室 成都 610031)

随着电算技术的迅猛发展，基于筋土分离技术[1-2]以真实模拟筋-土之间相互作用的有限元法日益受到重视，且目前多倾向采用岩土工程专用程序开展土工格栅[3-4]等加筋土工程的分析。Abdi等[5]利用PLAXIS中的线弹性Geogrid单元模拟土工格栅，探究了拉拔试验模型箱和土工格栅尺寸的合理性。唐晓松等[6]使用PLAXIS建立了土工格栅界面特性研究模型，界面采用弹塑性模型，对比分析了室内直剪试验与数值方法。蒋鑫等[7]则使用Phase2程序，将Liner-Geosynthetic、Joint单元分别用于模拟格栅和界面，探究了土-结构物界面对路基稳定性的动态影响。

显然众多可用于加筋土结构的数值分析程序模拟技术各有差异，若处理不当则可能造成计算结果失真。另一方面，路堤作用下筋-土作用规律的相关研究尚欠深入。本文选取由荷兰开发的PLAXIS[8]与由丹麦开发的OPTUM[9]这2款具有代表性的岩土有限元程序，从筋材模拟、筋-土界面模拟两处出发，比较其异同点，讨论这2款程序中界面参数的等效转换方法，进而探究筋材内力等关键力学响应随路堤填筑的全过程动态演变规律。

1 筋材模拟

PLAXIS与OPTUM中均内置了Geogrid结构单元用于模拟土工格栅，Geogrid单元节点构成见图1，注意各节点之间仅传递轴力，各节点均具有2个平移自由度。其中PLAXIS中的Geogrid单元实质上为多个节点构成的线单元，当土体采用15节点、6节点实体单元离散时，Geogrid单元分别由5节点、3节点所定义；OPTUM中则对应于15节点和6节点实体单元，Geogrid单元均为2节点的桁架单元。各个Geogrid单元前后相连时，上个单元的末尾节点与下个单元的首个节点空间坐标相同。

图1 Geogrid单元节点组成

针对Geogrid结构性单元，PLAXIS中有4种本构模型，分别为线弹性、弹塑性、弹塑性(N-ε)、黏弹性；OPTUM则只有弹塑性本构模型。2款程序关于筋材的本构模型参数见表1。

表1 筋材本构模型参数

相对而言，PLAXIS内嵌加筋材料的本构模型更为丰富，若需考虑加筋材料的黏弹性行为，则只能使用PLAXIS。

2 筋-土界面模拟

PLAXIS和OPTUM这2款程序均采用设置零厚度[10-11]的界面单元以刻画筋-土之间应力的传递，下面详细予以阐述。

2.1 PLAXIS筋-土界面模拟

PLAXIS中的筋-土界面由服从Mohr-Coulomb破坏准则的弹塑性模型模拟，界面抗剪强度为

τmax=ci+σntanφi

(1)

式中：ci、φi、σn、τmax分别为界面的黏聚力、内摩擦角、界面所受的有效正应力、最大剪应力。除Mohr-Coulomb剪应力破坏准则，拉伸截断准则也可以被考虑在界面中，相应的屈服包络线见图2a)，拉伸强度指定为σt。

图2 考虑拉伸截断的屈服包络线

界面的强度、刚度等性质与筋材相邻土体有关，PLAXIS中界面性质的定义方法有2种，第一种方法是为相邻土体指定一个小于1.0的折减因子Ri(Ri=1.0时表示筋材与相邻土体假定为完全黏合状态)，此时界面的强度参数(用下标i予以标识)由式 (2)～(5) 定义。

ci=Ricsoil

(2)

tanφi=Ritanφsoil

(3)

σt,i=Riσsoil

(4)

(5)

式中：csoil、φsoil、ψsoil分别为相邻土体的黏聚力、内摩擦角、剪胀角。

而界面的刚度参数如剪切模量Gi、弹性模量Ei、泊松比υ、压缩模量Eoed则由式 (6)～(9) 定义。

(6)

Ei=2Gi(1+2υi)

(7)

υi=0.45

(8)

(9)

式中：Gsoil为相邻土体的剪切模量，PLAXIS界面的泊松比υi为定值0.45。PLAXIS还可采用kn/ks模式定义界面刚度，此时界面法向刚度kn,切向刚度ks与上述模量符合式 (10) 的换算关系。

(10)

式中：ti为界面虚拟厚度，由虚拟厚度因子和单元尺寸定义。当界面法向、切向刚度可从试验得到时，此模式直接定义界面刚度参数相对更为方便。

第2种定义界面性质的方法则是为界面单独设置土体材料。此种情况下，土体材料施加于界面仅用于计算，界面的实际厚度依旧视作为0，界面的各种参数仍然可以通过上述各式计算。

2.2 OPTUM筋-土界面模拟

OPTUM中筋-土界面同样基于Mohr-Coulomb破坏准则的弹塑性模型[9]。界面拉伸截断准则对应的屈服包络线见图2b)，与PLAXIS不同之处在于除拉伸强度σt，还考虑拉伸截断角φt。若需考虑较为复杂的界面拉伸截断设置时需使用OPTUM。

界面的刚度参数设置则可采用E-υ(弹性模量Ei，泊松比υi)或K-G(体积模量Ki，剪切模量Gi)方式，Ki与Ei的关系满足式 (11)。

(11)

在OPTUM程序中，筋-土界面是Geogrid结构单元的一部分。在定义Geogrid单元时可为筋-土界面指定单独设置的土体材料(零厚度)；亦可为界面指定一个强度折减因子Ri，然后按照c-φ或c-tanφ模式予以折减，分别满足式 (12)、(13)。

ci=Ricsoil，φi=Riφsoil

(12)

ci=Ricsoil, tanφi=Ritanφsoil

(13)

2.3 2款程序筋-土界面性质的等效转换

前述表明，PLAXIS中的折减模式实则对应于OPTUM中c-tanφ的模式。因程序自身设定，PLAXIS中界面泊松比恒定为0.45，故如考虑PLAXIS与OPTUM刚度参数等效，则当在PLAXIS中为筋材相邻土体指定折减因子时，可利用式 (6)～(8), (11) 中关系在OPTUM中为界面单独设置材料使得其刚度参数满足式 (14)。

(14)

在强度参数等效方面，PLAXIS与OPTUM在设置拉伸截断时有所区别。由图2可知，对应于PLAXIS，实际上φi固定为90°，故计算时在OPTUM中设置φt=90°；剪胀角ψi均设置为0。其余强度参数2款程序则保持一致。

可利用此等效方法建立计算模型，保证2款程序在同一问题的界面模型一致，从而使得2款程序计算结果具有平行可比及规律映证的前提。

3 加筋土路堤数值模型及结果分析

3.1 数值模型建立

3.1.1几何模型

模型示意见图3，该路堤高6 m、宽24 m，边坡坡度为1∶1.5，在路堤底部通长铺设一层土工格栅；假定为平面应变问题，并考虑到对称性，取右半结构计算。模型宽度取65 m，约为3.1倍路堤底部宽度，地基计算深度取20 m，以尽量减小边界条件的影响。地下水位线在地基表面下1 m。

图3 加筋土路堤有限元计算模型(单位：m)

3.1.2本构模型与材料参数

路堤与地基土材料均假定为Mohr-Coulomb模型，具体参数见表2。为便于平行比较，采用弹塑性本构模型定义土工格栅，轴向抗拉刚度为3 000 kN/m，极限抗拉强度为100 kN/m。

表2 土层材料参数

采用前述2.3节中筋-土界面参数等效转换方法定义筋-土界面，具体阐述如下。在PLAXIS中，定义土工格栅相邻土体时为其赋值Ri=0.7，上、下界面分别对应路堤填土、淤泥质土；在OPTUM中分别为上、下界面单独设置2种材料，其参数同PLAXIS中Ri=0.7时等效。上、下界面具体参数见表3，其中刚度与强度参数分别由式 (2)～(9)、(14) 计算得到，界面虚拟厚度2款程序均取为6.7 cm。

表3 上、下界面参数

3.1.3单元剖分

为尽量减少网格密度对计算结果的影响，2款程序均采用高精度15节点三角形单元离散地基与路堤，加筋位置处适当加密网格；而土工格栅方面PLAXIS、OPTUM分别采用5节点Geogrid单元、2节点Geogrid单元模拟，格栅上下侧均设置界面单元模拟土-筋的相互作用。网格剖分具体见图3。

3.1.4边界条件及施工模拟

模型底部双向位移皆约束，左右两侧竖直向自由、水平向位移约束，其他位置自由。格栅最左端随左边界一同设置为水平向位移固定，竖直向自由。在地基初始地应力平衡后分3层加载激活路堤以模拟路堤分步逐层填筑，每层高度2 m。

3.2 计算结果与分析

3.2.1筋材轴向拉力

路堤逐层填筑过程中2款程序计算所获筋材轴向力分布情况见图4。

图4 土工格栅轴向拉力分布

随着路堤的填筑，2款程序所获筋材轴向拉力N均相应增大，且增大幅度加剧，OPTUM计算结果略大于PLAXIS，且随填筑层数的增加，2款程序之间结果差值也逐渐增大。从轴向力的横向分布规律看，2款程序基本一致，靠近路堤中心一侧轴向力较大，靠近路堤坡脚(x=21 m)轴向力逐渐减小至0，随填筑层数的增加，坡脚处筋材末端轴向力趋近于0的长度逐渐增长；轴向力分布曲线的整体重心越发靠近路堤中心位置；填筑1、2层后轴向力曲线的峰值出现在路堤中心处，但在其他位置也出现盆状峰(分别位于x=16.5 m、x=12 m附近)；填筑第3层后，轴向力曲线则呈“勺状”分布，从路堤中心至坡脚，轴向力先增大后逐渐减小至0，OPTUM的峰值位置相对更靠左侧。

3.2.2界面剪切应力

尽管上、下界面的材料属性不同，但经分析发现模型中土工格栅上下界面的剪切应力τ分布规律大致呈对称分布，故仅绘制下界面剪应力分布见图5，并定义相对位移朝向路堤坡脚时(即下侧地基土相对格栅出现向右侧的相对剪切位移)剪应力为正。

图5 下界面剪应力分布

由图5可见，宏观上看，依次填筑3层路堤时，从路堤中心至坡脚位置，剪应力分布大致上都先出现负剪应力后逐渐变化为正剪应力，这是由于不同位置界面的相对位移方向不同，且正负剪应力转换点位置2款程序较为接近，整体上转换点位置随着路堤填筑逐渐向路堤中心线一侧移动，同前述3.2.1节中格栅轴向力重心移动方向一致，这应是随填筑层数的增加，路堤潜在滑动面逐渐向左侧移动所致。2款程序的剪应力计算差值最大处基本上都位于剪应力峰值位置附近。填筑第2、3层后，路堤填筑造成坡脚处的侧向推挤效应明显，导致界面相对位移较大，剪应力在路堤坡脚处急剧上升。

3.2.3路堤稳定性

2款程序利用强度折减计算得到的最危险滑动面及稳定安全系数Fs随路堤填筑而动态演变的规律见图6。

图6 填筑各层后安全系数(Fs)及危险滑动面

由图6可见，由PLAXIS、OPTUM分别确定的最危险滑动面有所不同，相对而言OPTUM得到的滑动面位置更靠深处(左下侧)，填筑第3层后，2款程序的滑动面最为相似，均从路堤中心附近向坡脚外侧发展。填筑第1、2层后滑动面差异明显，当填筑第1层后，PLAXIS、OPTUM滑动面分别出现在边坡坡面、坡面以下位置；当填筑第2层后，PLAXIS中滑动面更为靠近边坡且未穿过加筋位置，而OPTUM中的滑动面则已经穿过加筋位置朝更深处发展。滑动面位置出现差异的原因可能在于：①2款程序表征滑动面的方式不同，PLAXIS、OPTUM分别通过增量位移、塑性剪切耗散的云图来表征；②2款程序开展强度折减分析时涉及非线性收敛控制的相关参数，如荷载增量步、收敛精度等存在一定差异。但安全系数FS计算值二者相差不大，均随路堤填筑减小。

4 筋材受拉与界面抗剪的动态发展

显然加筋路堤的稳定性受到土体、筋材、筋-土界面的共同影响。不妨定义相对轴向力系数nN=N/NP。其中：NP为筋材极限抗拉强度(本文为100 kN/m)；定义界面相对剪切应力系数nτ=|τ/τmax|，τmax为界面上各点的抗剪强度，由式(1)计算得到。这2个系数数值大小愈趋近于1，表明筋材受拉或界面抗剪的作用发挥程度越大，越接近于破坏。图4、图5结果已表明2款程序计算得到的轴向力与界面剪切应力分布规律类似，而PLAXIS软件应力结果后处理更为方便，故仅提取PLAXIS计算结果，绘制填筑路堤各层后沿筋材方向相对力系数的变化规律见图7。

图7 相对轴向力与相对剪切应力系数分布

由图7可见，随路堤填筑，相对轴向力系数nN基本上在整个筋材长度上均有所提高，这表明筋材的整体拉力作用随着路堤的填筑而逐渐发挥。而筋土界面的抗剪强度与界面处的正应力密切相关，界面的正应力随路堤填筑而增大，界面抗剪强度提高。需要注意的是，图中nτ=0时代表此处界面相对位移为0，填筑第1、2层后相对剪切应力nτ在各个位置均未超过0.2，界面的抗剪强度发挥程度还较低，填筑第3层后靠近坡脚处位置nτ急剧上升至0.81，其余位置nτ仍较小。筋土界面抗剪作用集中发挥在靠近坡脚一侧。

5 结论

1) PLAXIS程序采用由5节点、3节点所定义的Geogrid单元模拟筋材，并拥有线弹性、弹塑性、弹塑性(N-ε)、黏弹性等4种筋材本构模型；OPTUM程序则采用2节点的Geogrid单元模拟筋材，只有弹塑性本构模型；OPTUM中考虑界面拉伸截断准则更为全面，定义界面强度折减的方式相较于PLAXIS更多样。

2) 本文所提出的界面属性等效变换方法较为有效，利用此方法可以保证2款程序中界面模型的一致性。2款程序计算结果差异主要受Geogrid单元类型、网格剖分方式等因素影响，计算时应尽量细分网格、注意网格疏密合理布设、提高网格剖分质量等以减少其影响。PLAXIS计算得到的轴向拉力小于OPTUM，2款程序所获下界面剪切应力差值在最大剪应力附近位置相较其他位置更大，OPTUM中危险滑动面位置更靠左下侧。

3) 随着路堤填筑进程深入，筋材轴向拉力整体作用逐渐发挥，而筋土界面抗剪作用则主要发挥在靠近坡脚一侧。设计时需注意提高坡脚附近局部位置的抗剪强度，而在考虑筋材轴向拉断破坏时，应重点关注靠近路堤中心或危险滑动面穿过的位置。