基于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) 隨著路堤填筑進程深入，筋材軸向拉力整體作用逐漸發揮，而筋土界面抗剪作用則主要發揮在靠近坡腳一側。設計時需注意提高坡腳附近局部位置的抗剪強度，而在考慮筋材軸向拉斷破壞時，應重點關注靠近路堤中心或危險滑動面穿過的位置。