利用PLAXIS軟件計算考慮降雨的邊坡穩(wěn)定性

孔郁斐，周夢佳，宋二祥，楊 軍，張龍英，施洪剛，劉 劍

(1.清華大學 土木工程安全與耐久教育部重點實驗室，北京 100084；2.中國建筑股份有限公司，北京100000）

利用PLAXIS軟件計算考慮降雨的邊坡穩(wěn)定性

孔郁斐1，周夢佳1，宋二祥1，楊 軍1，張龍英2，施洪剛2，劉 劍2

(1.清華大學 土木工程安全與耐久教育部重點實驗室，北京 100084；2.中國建筑股份有限公司，北京100000）

降雨條件下的邊坡穩(wěn)定性分析，需要同時考慮水的滲流與邊坡內(nèi)力進行耦合計算。非飽和土的特殊性質(zhì)增加了計算難度，其滲透性、強度都會隨含水量的變化而變化。介紹了非飽和土的有效應力原理，并比較了PLAXIS和Geo-Studio兩款軟件使用的有效應力原理和破壞準則的異同。PLAXIS軟件中用有限元法計算非飽和土的滲流問題，利用簡化的Bishop有效應力公式進行固結、變形及邊坡穩(wěn)定性計算，土體本構關系根據(jù)需要選用；Geo-Studio軟件中也使用有限元法計算非飽和滲流，但其邊坡安全系數(shù)計算方法為極限平衡法，強度準則為Fredlund雙變量理論。最后結合算例詳細介紹了PLAXIS軟件中進行非飽和土邊坡計算的建模方法。算例顯示，PLAXIS軟件中建立的計算模型可以準確反映降雨過程中邊坡安全系數(shù)的變化規(guī)律。今后的工程設計中，可以考慮使用這一軟件進行非飽和土邊坡設計計算。

非飽和土；降雨滲流；土坡穩(wěn)定；耦合計算；PLAXIS軟件；Geo-Studio軟件

降雨型滑坡是常見的自然災害，非飽和土邊坡在降雨過程中安全系數(shù)逐漸降低。降雨條件下的邊坡穩(wěn)定性計算是滲流分析和受力分析的耦合問題。土的滲透系數(shù)、強度都會隨其含水量變化而變化，因此瞬態(tài)滲流計算和安全系數(shù)計算中都必須考慮非飽和土的特殊性質(zhì)。經(jīng)過20世紀的發(fā)展，非飽和土理論已比較完善［1］，并向?qū)嵱没~進。巖土工程相關的數(shù)值計算軟件也正逐漸開發(fā)與非飽和土有關的功能，比較知名的PLAXIS、Geo-Studio等軟件，都能進行非飽和土邊坡穩(wěn)定性計算。

Geo-Studio是加拿大開發(fā)的用于巖土工程和巖土環(huán)境模擬計算的軟件，較早涉足非飽和土計算。它包括8個模塊，其中SEEP/W模塊可進行飽和/非飽和的二維穩(wěn)態(tài)/瞬態(tài)滲流計算，得出孔隙水壓分布［2］；然后將孔壓計算結果導入SLOPE/W模塊，利用極限平衡法計算邊坡安全系數(shù)［3］。對均質(zhì)土坡而言，降雨時表層土壤首先浸潤，深層土壤受影響較晚，故水分分布不均勻，使用“圓弧滑移面”的假設不一定合適；另外很多坡體本身就不均勻，真實滑移面形狀很復雜［4］，基于“圓弧滑移面”假設的極限平衡法不適用，SLOPE/W雖然也能搜索任意形狀滑移面，但有時算法不穩(wěn)定，結果可靠性有待驗證。

荷蘭開發(fā)的巖土工程有限元設計計算軟件PLAXIS使用“強度折減有限元法”計算邊坡穩(wěn)定性［5］，計算出的滑移面更接近真實滑移面，因此對非均勻的土坡使用這一方法結果更可靠。PLAXIS可以進行2D/3D的彈塑性計算、滲流計算、固結分析、安全分析、動力計算等，近幾年增加了非飽和土的計算功能，其2011版可以將非飽和土滲流計算和安全系數(shù)的強度折減有限元計算耦合［6］，值得探索和嘗試。

本文將首先介紹非飽和土有效應力原理，這是數(shù)值計算的重要理論基礎，并對PLAXIS與Geo-Studio進行了比較，接著詳細介紹PLAXIS中的建模方法，并結合算例討論PLAXIS的計算結果。

1 非飽和土有效應力原理

如前所述，邊坡的滲流計算、安全系數(shù)計算都要考慮非飽和土的特殊性質(zhì)，問題復雜，影響因素眾多，吳宏偉等［7-8］曾就相關的難點、計算方法作過詳細的闡釋，本文不多介紹。兩個相關商業(yè)軟件的主要區(qū)別在于使用的有效應力原理和強度準則不同，本節(jié)將重點介紹其區(qū)別。

飽和土的有效應力原理由太沙基提出［9-10］，該理論認為：①飽和土體內(nèi)任意一平面上受到的總應力可分為有效應力和孔隙水壓力兩部分，在同一點處水壓向各方向都是相等的；②土的變形(壓縮）與強度的變化都只取決于有效應力變化，其表達式為：

有學者將其擴展到孔隙中有多種流體(水、油、氣等）的情況［11］：

該式又叫“廣義有效應力公式”。式中：uβ為第β種流體的分壓，αβ為其分項系數(shù)，表征其在總應力中的貢獻。不難理解，有多種流體共同存在時，每種流體分壓的分項系數(shù)αβ≤1。土是由土骨架、水、空氣組成的三相體：飽和土中僅有水，飽和度S=1，αw=1，αs=0，式(2）退化為式(1）；干燥土僅有空氣，S=0，αa=1，αw= 0；非飽和土的2個系數(shù)均在0～1之間變化。

圖1 飽和度較大的非飽和土Fig.1 Sketch of the nearly-saturated soil

考察a-a面以上的土骨架，由土骨架受力平衡可得：

故

此時仍然有αw=1，αa=0，與飽和土類似。

但若土中空氣較多，內(nèi)部空氣不再只以氣泡的形式存在、能與大氣連通，則土中氣體的平衡方程中不僅包含浮力、土骨架作用力，還包含外部空氣壓力以及表面張力，有效應力就更加復雜了。1960年畢肖普等提出了修正的有效應力公式(Bishop有效應力公式）：

式中：χ是與飽和度有關的系數(shù)，飽和土χ=1，干土χ=0。這一公式是式(2）的簡化形式，在式(2）中令αw=χ，αa=1-χ可得到式(5）。式中σij-uaδij被稱作凈法向應力張量，ua-uw被稱作基質(zhì)吸力。土飽和時，ua=uw，凈法向應力等于飽和土的有效應力，基質(zhì)吸力為零；土不飽和時，基質(zhì)吸力為正。

將Bishop有效應力公式代入摩爾-庫倫準則(τ=c′+σ′tanφ′）可得基于Bishop有效應力的摩爾-庫倫準則：

式中：φ′為有效內(nèi)摩擦角，χtanφ′體現(xiàn)了基質(zhì)吸力對抗剪強度的貢獻。

從數(shù)學上看這一簡化非常簡單，但在工程上它非常有意義，我們可以在不同飽和度下進行土的三軸試驗，確定函數(shù)χ(S），就可以得到非飽和土的有效應力，再根據(jù)“土的變形及強度變化都只取決于有效應力的變化”的假設就可以將原有的飽和土應力應變分析的理論應用于非飽和土。Zerhouni總結了幾種不同土的χ-S關系曲線(見圖2），這些數(shù)據(jù)均由試驗測得［12］，可見χ(S）函數(shù)沒有特別的規(guī)律可言。1和4兩種土的飽和度達到90%左右時，χ=1，與圖1的情形有相似之處。

圖2 幾種不同土的χ-S關系曲線Fig.2 Curves for several different soils

不過隨著研究的深入，人們也發(fā)現(xiàn)Bishop有效應力公式存在一定的局限性。例如，它無法解釋一些土在浸濕過程中的體縮現(xiàn)象：保持土的總應力不變，使土浸濕，則土中有效應力減小，應變本應持續(xù)增加，但某些土會先濕脹、再濕陷。所以有人提出疑問：在非飽和土中直接照搬飽和土的有效應力的概念似乎有些牽強，當水和空氣共同存在時，針對飽和土所作的“土的變形及強度變化都只取決于有效應力的變化”假設還成立么?

不一定成立。在飽和土中水的體積應變很小，體變的內(nèi)力功可以忽略，土進入塑性以后土骨架之間的摩擦使能量耗散；非飽和土中空氣可以壓縮，空氣-水的界面可以變形，空氣的壓縮和表面張力做功都會消耗能量，故耗能能力受到流體性質(zhì)的影響，有效應力不再是唯一的決定因素。

D.G.Fredlund等［1］提出，可以用凈法向應力和基質(zhì)吸力兩個獨立的應力張量來描述土的狀態(tài)，克服單一有效應力變量的缺陷，在Fredlund雙變量理論中，土的應力應變關系如下［1］：

似乎有了這個公式之后，對非飽和土的計算就能得心應手了，其實不然，要獲得完整的參數(shù)著實不易，兩個應力張量讓破壞面更難描述。比較常用的還是更簡單直觀的假設：線彈性、理想彈塑性或者剛塑性。

摩爾-庫倫準則也被擴展到雙變量理論中(擴展的摩爾-庫倫準則）：

式中：tanφb表征基質(zhì)吸力對抗剪強度的貢獻，且φb并非常數(shù)，隨飽和度變化。

沈珠江院士曾指出，既然φb并非常數(shù)，在強度準則里式(8）就不見得比式(6）優(yōu)越［13］。很多學者提出了更為復雜的強度公式，但現(xiàn)在仍然是擴展的摩爾-庫倫準則使用最廣泛。

PLAXIS和Geo-Studio分別使用了以上兩種有效應力理論，并且都作了簡化。

PLAXIS中使用了Bishop有效應力公式，由于χ(S）函數(shù)很復雜，直接令χ=S；也就是說在圖2中取對角線，這與實際情況存在一些差距，最接近圖2中的第7種土(White Clay），用于其他土時需要對參數(shù)進行評估。這一有效應力可以用于多種本構模型，當然最簡單的是理想彈塑性模型及摩爾-庫倫準則。雖然Bishop有效應力公式有局限性，但是處理非特殊土時，這一理論是能夠勝任的。

Geo-Studio中使用了Fredlund雙變量理論，并假設φb為常數(shù)，計算時需根據(jù)滑移面穿過區(qū)域的含水量變化范圍，選取合適的φb，否則極易低估或高估基質(zhì)吸力對強度的貢獻。雖然Fredlund雙變量理論比Bishop公式有優(yōu)勢，但使用極限平衡法計算邊坡穩(wěn)定性時不涉及土的應力-應變關系，優(yōu)勢便無從體現(xiàn)了。

在此也介紹一下非飽和土的土水特征曲線(SWCC）和滲透性函數(shù)。非飽和土的體積含水率θw，基質(zhì)吸力ψ(=ua-uw），滲透系數(shù)k都會隨著飽和度S變化。ψ(θw）或ψ(S）函數(shù)關系即為“土水特征曲線”，k(S）或k(ψ）函數(shù)關系即為“滲透性函數(shù)”。這兩種函數(shù)也較為復雜，研究者們曾提出很多種函數(shù)形式，比較知名的有Gardner方程、van Genuchten方程［14］、Fredlund-Xing方程等［15］。

無論如何，用兩款軟件進行計算時都應謹慎選用數(shù)據(jù)。除了有效應力原理、強度準則和安全系數(shù)計算方法，在其他方面這兩款軟件有很多相似之處，都能進行多種邊界條件設定、內(nèi)置了較多土水特征曲線和滲透性函數(shù)的方程。值得一提的是，PLAXIS還內(nèi)置了數(shù)個和非飽和土數(shù)據(jù)庫，如USDA(美國農(nóng)業(yè)部）數(shù)據(jù)庫、Hypres數(shù)據(jù)庫等，方便用戶在試驗數(shù)據(jù)有限的情況下選取合適的參數(shù)快速建模。

2 PLAXIS建模與計算

在PLAXIS 2D(2011）中建立如圖3所示的均質(zhì)邊坡模型，模型大小為30 m×10 m，邊坡高5 m，長8 m，地下水位位于坡腳下方3 m處，網(wǎng)格劃分疏密程度為“很細”；模型為平面應變模型，單元為15節(jié)點單元。土體類型選用USDA數(shù)據(jù)庫中的“沃土”：粒徑小于2 μm的質(zhì)量占總質(zhì)量的比例p0.0002=20%，2～50 μm的占40%，50 μm～2 mm的占40%，按國內(nèi)分類屬于細粒土。

土的飽和滲透系數(shù)ksat=0.247 9 m/d，土的滲透性函數(shù)和土水特征曲線使用van Genuchten方程擬合［16］，結果如圖4所示，圖中兩條曲線由軟件自動生成，其中kr=k/ksat表示相對滲透系數(shù)，為無量綱量，S為飽和度。容易看出飽和度為0.7時基質(zhì)吸力水頭約為0.5 m，相對滲透系數(shù)為0.01。

圖3 PLAXIS中的邊坡模型及網(wǎng)格Fig.3 The slope model and mesh used in PLAXIS

圖4 模型中土的滲透性函數(shù)(左）和土水特征曲線(右）Fig.4 Permeability function and SWCC used in the model

對于一個地下水位恒定、表面沒有降雨/蒸發(fā)的邊坡而言，其地下水位以上的基質(zhì)吸力分布是線性的，水位以上x處基質(zhì)吸力水頭為x。但天然邊坡的基質(zhì)吸力水頭一般有上限，上限為零點幾米到幾米不等，可以在表面定義一個極小的入滲速率(或稱“虛擬降水”），使得內(nèi)部穩(wěn)態(tài)滲流時基質(zhì)吸力最大值在我們需要的范圍內(nèi)。本算例計算中我們首先進行穩(wěn)態(tài)滲流計算，令上表面有1 mm/d的虛擬降水，這樣邊坡中最大基質(zhì)吸力水頭為0.719。

計算以下兩工況：工況1，降雨速率為50 mm/d，持續(xù)時間分別為0.5，1.0和2.0 d；工況2，降雨速率為10 mm/d，持續(xù)時間分別為2.5，5.0和10.0 d。兩種工況最后的總降雨量均為100 mm。

表1和圖5展示了算例中的邊坡在不同工況下的安全系數(shù)、滑移面深度計算結果。圖6、圖7分別為不同時刻邊坡的基質(zhì)吸力分布和滑移面，計算所得的7組滑移面都比較淺，工況1中滑移面深度從1.82 m降至1.11 m，降雨結束時安全系數(shù)接近1.00，這個理想的邊坡已經(jīng)非常危險了；工況2中滑移面深度從1.82降至1.45 m。

圖5 兩種工況下邊坡安全系數(shù)隨降雨量的變化Fig.5 Variation of safety factor with rainfall in different working conditions

表1 安全系數(shù)、滑移面深度隨降雨量的變化Tab.1 Variations of safety factor and depth of sliding surface with rainfall

圖6 不同工況下的邊坡基質(zhì)吸力分布Fig.6 Matric suction distribution in different working conditions

圖7 不同時刻強度折減法計算得到的邊坡位移云圖Fig.7 Shadings of displacements from φ/c reduction calculation

可以看出降雨過程中該邊坡安全系數(shù)明顯減小，且工況1中安全系數(shù)降幅更大。隨著雨水入滲，坡中土的含水量會自上而下逐漸增加，隨著飽和度的提高，土的基質(zhì)吸力下降、有效應力減小，強度降低，因此邊坡的安全系數(shù)會降低。降低的速率取決于很多因素，在總降雨量相同的情況下，如果降雨速率快、時間短(工況1），雨水會主要集中在土體上部，淺層滑移面經(jīng)過處的土的強度降幅大，更容易發(fā)生淺層滑坡；如果降雨速率慢，雨水有充足的時間下滲至更深處，淺層滑移面經(jīng)過處的土的強度降幅小。這與筆者此前用Geo-Studio計算發(fā)現(xiàn)的規(guī)律相符［17］。

3 結 語

本文總結了非飽和土有效應力原理的不同形式，在此基礎上比較了PLAXIS和Geo-Studio兩款軟件計算原理的差異，并介紹了在PLAXIS中進行考慮降雨的邊坡穩(wěn)定性計算的建模方法。通過本文的分析可以發(fā)現(xiàn)，PLAXIS中使用簡化的Bishop有效應力原理可以反映非飽和土的性質(zhì)，在使用中需要注意合理選取參數(shù)，與Geo-Studio使用的雙變量理論相比，兩款軟件的強度準則沒有明顯的優(yōu)劣之分。算例顯示，PLAXIS中建立的計算模型可以正確反映降雨過程中邊坡安全系數(shù)的變化規(guī)律。在今后的工程設計中，可以考慮使用PLAXIS軟件進行非飽和土邊坡設計計算。

參 考 文 獻：

［1］FREDLUND D G，RAHARDJO H.Soil mechanics for unsaturated soils［M］.New York：John Wiley&Sons，1993.

［2］GEO-SLOPE International Ltd.Seepage modeling with SEEP/W 2007-an engineering methodology［M］.3rd ed.Canada：GEOSLOPE International Ltd，2008.

［3］GEO-SLOPE International Ltd.Stability modeling with SLOPE/W 2007 version-an engineering methodology［M］.3rd ed. Canada：GEO-SLOPE International Ltd，2008.

［4］張建民，王睿，張嘎.降雨條件下含軟弱夾層土坡的離心模型試驗研究［J］.巖土工程學報，2010(10）：1582-1587. (ZHANG Jian-min，WANG Rui，ZHANG Ga.Centrifuge modeling of rainfall-induced deformation of slopes with weak layers［J］.Chinese Journal of Geotechnical Engineering，2010(10）：1582-1587.(in Chinese））

［5］宋二祥.土工結構安全系數(shù)的有限元計算［J］.巖土工程學報，1997，19(2）：1-7.(SONG Er-xiang.Finite element

analysis of safety factor for soil structures［J］.Chinese Journal of Geotechnical Engineering，1997，19(2）：1-7.(in Chinese））

［6］PLAXIS.PLAXIS Scientific Manual 2011［M］.［S.l.］：［s.n.］，2011.

［7］吳宏偉，陳守義，龐宇威.雨水入滲對非飽和土坡穩(wěn)定性影響的參數(shù)研究［J］.巖土力學，1999，20(1）：1-14.(WU Hong-wei，CHEN Shou-yi，PANG Yu-wei.Parametric study of effects of rain infiltration on unsaturated slopes［J］.Rock and Soil Mechanics，1999，20(1）：1-14.(in Chinese））

［8］NG C W W，SHI Q.Influence of rainfall intensity and duration on slope stability in unsaturated soils［J］.Quarterly Journal of Engineering Geology，1998，31(2）：105-113.

［9］TERZAGHI K.Theoretical soil mechanics［M］.New York：John Wiley&Sons，1943.

［10］陳仲頤，周景星，王洪瑾.土力學［M］.北京：清華大學出版社，1994.(CHEN Zhong-yi，ZHOU Jing-xing，WANG Hongjin.Soil mechanics［M］.Beijing：Tsinghua University Press，1994.(in Chinese））

［11］NUTH M，LALOUI L.Effective stress concept in unsaturated soils：Clarification and validation of a unified framework［J］. International Journal for Numerical and Analytical Methods in Geomechanics，2008，32(7）：771-801.

［12］ZERHOUNI M I.R?le de la pression interstitielle négative dans le comportement des sols-Application au calcul des routes［D］.Paris：Ecole Centrale Paris，1991.

［13］沈珠江.當前非飽和土力學研究中的若干問題［C］∥區(qū)域性土的巖土工程問題學術論文集.北京：原子能出版社，1996：1-9.(SHEN Zhu-jiang.Current issues of unsaturated soil mechanics［C］∥Regional Colloquium on Soil Geotechnical Proceedings.Beijing：Atomic Energy Press，1996：1-9.(in Chinese））

［14］VAN GENUCHTEN M T.A closed-form equation for predicting the hydraulic conductivity of unsaturated soils［J］.Soil Science Society of America Journal，1980，44(5）：892-898.

［15］ LEONG E C，RAHARDJO H.Review of soil-water characteristic curve equations［J］.Journal of Geotechnical and Geoenvironmental Engineering，1997，123(12）：1106-1117.

［16］GALAVI V.Groundwater flow，fully coupled flow deformation and undrained analyses in PLAXIS 2D and 3D［R］.Delft：Research Department，PLAXIS，2010.

［17］孔郁斐.降雨對非飽和土邊坡穩(wěn)定性影響研究［D］.北京：清華大學，2012.(KONG Yu-fei.Rainfall′s effect on the

stability of unsaturated slopes［D］.Beijing：Tsinghua University，2012.(in Chinese））

Slope stability analysis in consideration of rainfall influence based on PLAXIS software

KONG Yu-fei1，ZHOU Meng-jia1，SONG Er-xiang1，YANG Jun1，ZHANG Long-ying2，SHI Hong-gang2，LIU Jian2
(1.Key Laboratory of Ciυil Engineering Safety and Durability of China Ministry of Education，Tsinghua Uniυersity，Beijing 100084，China；2.China State Construction Engineering Co.，Ltd.，Beijing 100000，China）

In order to analyze the slope stability in consideration of the rainfall influences，it is necessary to couple the calculation of the rainfall infiltration and the slope internal force.Special characteristics of the unsaturated soils make the calculation more complex，because the soil′s hydraulic permeability and strength would change with moisture content.In the paper the principle of the effective stress of the unsaturated soils is introduced.The similarities and differences of the principle of effective stress and failure criteria used in the PLAXIS and Geo-Studio softwares are compared.In PLAXIS，the finite element method is used for unsaturated seepage calculation and a simplified Bishop effective stress formula is applied to the calculation of consolidation，deformation and slope stability.The constitutive models for the saturated soils can be used on that basis.The Geo-Studio uses the finite element method for the unsaturated seepage calculation and the limit equilibrium method for the calculation of the slope safety factor，in consideration of Fredlund′s bilinear strength criterion.Finally，details of the method for modeling the unsaturated slopes in PLAXIS are presented with calculation examples.The calculation examples show that the models in PLAXIS can accurately reflect rainfall-induced slope safety factor variation.Therefore，this software can be adopted for the unsaturated soil slope design and calculation in the future.

unsaturated soil；rainfall infiltration；slope stability；coupling calculation；PLAXIS software；Geo-Studio software

TU413.6+2 文獻標心碼：A

1009-640X(2014）03-0070-07

2013-11-12

孔郁斐(1991-），男，江蘇泰興人，博士研究生，主要從事巖土工程方面的研究。E-mail：kyf12@mails.tinghua.edu.cn