程燦宇，羅富榮，戚承志，王 霆

（1. 北京建筑工程學院 工程結(jié)構(gòu)與新材料北京市高校工程研究中心，北京 100044；2. 北京市軌道交通管理公司，北京 100037）

1 引 言

我國是一個地質(zhì)災害多發(fā)的國家，據(jù)不完全統(tǒng)計，2010年僅四川省就因地質(zhì)災害造成直接經(jīng)濟損失90多億元[1]，而且滑坡泥石流等地質(zhì)災害每年都給人民群眾的生命和財產(chǎn)安全帶來巨大的損失。我國西高東低的地勢導致我國西部地區(qū)地質(zhì)災害問題更加嚴重，伴隨著西部大開發(fā)的推進和大型工程項目的快速上馬，在我國西南和黃土高原地區(qū)以及西北部新疆地區(qū)地質(zhì)災害問題變得更加嚴重，涌現(xiàn)出一系列亟待解決的邊坡穩(wěn)定問題。因此，在我國邊坡穩(wěn)定分析現(xiàn)階段仍是一個很大的課題。

在邊坡穩(wěn)定分析領域一種常用的方法是有限元強度折減法，早在1975年該方法就被Zienkiewice[2]用來求解邊坡穩(wěn)定問題，只是當時由于需要花費大量的機時而在具體應用領域受到限制。Wong[3]給出了用有限元方法分析邊坡穩(wěn)定時誤差產(chǎn)生的原因?，F(xiàn)在，隨著計算機硬件技術(shù)和有限元軟件技術(shù)的飛速發(fā)展，運用有限元強度折減法分析邊坡穩(wěn)定已經(jīng)成為新的趨勢[4]。Ugai[5]、Matsui和 San[6]、Ugai和Leshchinsky[7]、Griffiths 和 Lane[8－9]、Manzari 和Nour[10]等都對此做了進一步的研究，研究結(jié)果均表明：有限元強度折減法能得到比較穩(wěn)定的與極限平衡法接近的安全系數(shù)和臨界滑動面。應用于巖土工程的分析軟件也都可以采用強度折減法計算邊坡穩(wěn)定問題，但目前還沒人研究針對具體的工況條件各軟件計算邊坡穩(wěn)定時的偏差以及各軟件計算特點問題。

本文研究了FLAC3D、MIDAS/GTS和ANSYS在分析邊坡穩(wěn)定時各自的特點；配合D-P屈服準則和M-C屈服準則，計算了軟黏土、弱膨脹土、硬黏土3種工況下不同坡度邊坡的安全系數(shù)；研究了3種工況條件下M-C屈服準則和D-P屈服準則計算結(jié)果的偏差；分析了3種工況下滑動面的規(guī)律以及各軟件得到的滑動面特征；總結(jié)了軟件計算時的注意事項。本文的研究可以為進行邊坡穩(wěn)定分析的研究人員提供參考。

2 有限元強度折減法計算原理

Duncan[11]指出，邊坡安全系數(shù)可以定義為使邊坡剛好達到臨界破壞狀態(tài)時對土的剪切強度進行折減的程度。通過逐步減小抗剪強度指標，將黏聚力c、內(nèi)摩擦角φ值同時除以折減系數(shù)K，從而得到一組新的強度指標c′、φ′，反復計算直至邊坡達到臨界破壞狀態(tài)，此時的強度折減系數(shù)即為該邊坡的安全系數(shù)。

用這種方法分析邊坡穩(wěn)定性時，在彈塑性有限元數(shù)值分析中，先將折減系數(shù)K取的小一些，以保證開始時是一個近乎彈性問題，然后不斷增加折減系數(shù)K，反復分析邊坡直至K增加至某一值時邊坡達到臨界狀態(tài)。有限元強度折減法的優(yōu)點是安全系數(shù)可以直接求出，滑裂面的形狀和位置都不需要事先求出，而且從計算過程中可以看到土坡逐步破壞的過程。

3 強度準則及失穩(wěn)判據(jù)

巖土材料具有復雜的本構(gòu)特性，而邊坡的穩(wěn)定性分析主要關(guān)心的是力和強度問題，而不是位移和變形問題，因此，可在有限元強度折減法中采用理想彈塑性本構(gòu)模型[12]。

對于巖石、混凝土和土壤等顆粒狀材料其受壓屈服強度遠大于受拉屈服強度，且材料受剪時顆粒會膨脹，常用的Mises屈服條件不適用于這類材料，巖土工程中常采用M-C屈服準則和D-P準則[13]。其表達式如下：

式（4）中的α、k與（3）中的常數(shù)c、φ關(guān)系為

此時D-P準則屈服面在π平面上的投影為M-C屈服準則在π平面上投影的六邊形的外角點外接圓[14]。

邊坡失穩(wěn)判據(jù)目前有3種比較主流的觀點：①有限元數(shù)值迭代不收斂判據(jù)[15－17]；②特征部位位移突變判據(jù)[18－20]；③等效塑性應變區(qū)貫通判據(jù)[21－22]。前人大量的研究成果表明，這3種有限元強度折減法邊坡失穩(wěn)判據(jù)具有統(tǒng)一性[23]，因此，本文根據(jù)實際情況綜合考慮判據(jù)1和判據(jù)3作為本文計算中判定邊坡失穩(wěn)的依據(jù)。

4 不同軟件計算結(jié)果的對比分析

FLAC3D、MIDAS/GTS和ANSYS這3款軟件中除了ANSYS都既提供了D-P模型又提供了M-C模型，所以在對比不同軟件的計算結(jié)果的同時，本文同時對比了不同軟件不同模型下計算結(jié)果的差異。嚴格地說，這4款軟件中FLAC3D應用的是有限差分程序，不是嚴格意義上的有限元程序，在FLAC3D中對控制偏微分方程進行的是直接逼近來推導方程[24]，這與有限元法不同。因此，通過本文對FLAC3D計算結(jié)果和其他軟件計算結(jié)果的對比，在一定程度上也可以看到在邊坡分析中有限元法和有限差分法的差別。

本文計算中選用3種材料參數(shù)，分別為：工況1、軟黏土；工況2、硬黏土；工況3、弱膨脹土[25]。材料參數(shù)如表1所示。

表1 材料參數(shù)Table 1 Material parameters

4.1 算例

針對研究的問題：不同軟件計算結(jié)果的對比分析，本算例采用在許多論文中都得到印證的一個算例[8,23－25]作為對比的參照：坡高為20 m，邊坡傾角為 45°。

4.1.1 M-C屈服準則下計算結(jié)果的對比

3種工況下利用FLAC3D和MIDAS計算，利用其中內(nèi)置的強度折減計算模塊調(diào)整相同的收斂參數(shù)得到的安全系數(shù)如表2所示。

表2 不同工況下求得的穩(wěn)定安全系數(shù)Table 2 Safety factors by different conditions

從表2可以看出，3種工況下MIDAS的計算結(jié)果都要比FLAC3D的計算結(jié)果大，其中硬黏土工況下二者偏差最大，達到10.81% 。

為了得到進一步的規(guī)律，調(diào)整邊坡坡度利用MIDAS和FLAC3D分別進行3種工況下邊坡安全系數(shù)的計算，計算發(fā)現(xiàn)MIDAS和FLAC3D的計算結(jié)果差異不僅與工況有關(guān)而且還與邊坡的坡角有關(guān)，計算結(jié)果如圖1所示。

從圖1中可以看出，當邊坡土體為軟黏土時，F(xiàn)LAC3D和MIDAS的計算結(jié)果吻合得比較好，經(jīng)計算結(jié)果最大偏差不超過 3%，而且隨著坡度的增加二者計算結(jié)果的差別在 45°坡時趨于最小。與此相對，當邊坡土體為硬黏土或弱膨脹土時，F(xiàn)LAC3D和MIDAS的計算結(jié)果偏差相對較大，最大偏差發(fā)生在45°的硬黏土坡，最大偏差為10.81%。當坡度較小時FLAC3D的計算結(jié)果要明顯大于MIDAS的計算結(jié)果，當坡度達到 30°時二者計算結(jié)果趨于吻合，且MIDAS的計算結(jié)果略大于FLAC3D的計算結(jié)果。而當坡度繼續(xù)增加時二者的計算結(jié)果偏差也在增大，而且這時MIDAS的計算結(jié)果大于FLAC3D的計算結(jié)果，與軟黏土工況不同的是：在工況2、3下邊坡坡度達到 45°時二者偏差達到最大，隨著坡度的繼續(xù)增加，相對偏差趨于減小。

圖1 3種工況下MIDAS和FLAC計算的安全系數(shù)對比Fig.1 Comparison of safety factors for three conditions calculated by MIDAS and FLAC

4.1.2 D-P屈服準則下計算結(jié)果的對比

選擇同樣的3種工況在ANSYS和MIDAS中選擇D-P準則進行手動折減，分別計算15°、30°、45°、60°時邊坡的穩(wěn)定安全系數(shù)。手動折減時設置和上面M-C準則下相同的收斂容差值。得到典型的45°坡時的安全系數(shù)如圖2所示。

圖2 45°坡不同工況下安全系數(shù)對比（D-P準則）Fig.2 Safety factors for 45 degrees slope under different conditions (D-P criterion)

從圖2中可以看出當坡體材料為軟黏土和硬黏土時，MIDAS計算的安全系數(shù)比ANSYS計算得到的結(jié)果大，最大相差5.73%，邊坡土體為弱膨脹土時ANSYS計算得到的安全系數(shù)比MIDAS計算得到的結(jié)果大1.69%。可見45°坡時MIDAS和ANSYS采用D-P準則計算時得到的安全系數(shù)具有較高的吻合度。

表3 改變坡角時3種工況下邊坡的安全系數(shù)Table3 Slope safety factors under different slope angles and conditions

在研究坡角對軟件計算偏差的影響時，為了減少采用塑性區(qū)貫通作為判別準則產(chǎn)生的人為誤差，計算中設定相同的收斂容差后統(tǒng)一以計算不收斂作為邊坡破壞的依據(jù)。從表3中可以得到比較一致的結(jié)論：采用D-P模型時，MIDAS計算的安全系數(shù)大于ANSYS計算的安全系數(shù)。對硬黏土進行分析時偏差最明顯，最大相差14.58%，最小相差1.83%，且呈現(xiàn)出偏差隨邊坡坡度的增加而增大的趨勢。在對軟黏土邊坡進行分析時也呈現(xiàn)出這種趨勢。但對弱膨脹土邊坡進行分析時，偏差隨邊坡坡度變化的規(guī)律沒有硬黏土時那么明顯，大約都在5%左右。

4.1.3 不同軟件在D-P、M-C準則下計算結(jié)果的對比

在MIDAS和ANSYS中調(diào)用的D-P屈服準則的屈服面在π平面上的投影為 M-C屈服準則的屈服面在π平面上投影的外接圓（見圖3），因此，從理論上說采用 M-C屈服準則得到的邊坡穩(wěn)定安全系數(shù)小于采用D-P準則得到的邊坡穩(wěn)定安全系數(shù)，本例通過計算印證了這一推導，并且綜合各軟件計算結(jié)果得到了針對具體問題時，兩個準則計算結(jié)果的偏差大?。ㄒ妶D4）。

圖3 各屈服準則在π 平面的投影Fig.3 Yield surfaces on the π plane

圖4 3種工況下4種計算方式的結(jié)果對比Fig.4 Comparison of results by using four calculation methods for three conditions

從圖中可以看到利用MIDAS（D-P）的計算結(jié)果和MIDAS（M-C）的計算結(jié)果有較好的一致性，同時ANSYS（D-P）的計算結(jié)果和FLAC（M-C）的計算結(jié)果有較好的一致性。因此，我們進行如下的偏差統(tǒng)計，如表4所示。

表4 4種計算模式的相對偏差統(tǒng)計Table 4 Statistics of relative deviation of four calculated modes

從表4可以看出，D-P準則和M-C準則同時用于軟黏土工況和弱膨脹土工況時計算結(jié)果的偏差一致且相對較小，而當D-P準則和M-C準則分別被用來計算硬黏土邊坡時，計算結(jié)果的偏差明顯變大。由于D-P準則計算得到的安全系數(shù)偏大，因此，當D-P準則被用于硬黏土邊坡分析時應該特別小心。

4.2 計算分析

4.2.1 關(guān)于折減系數(shù)的選取

采用手動折減安全系數(shù)進行邊坡穩(wěn)定分析時，最初一定要選取較小的折減系數(shù)。比如，采用ANSYS進行分析時，當邊坡坡角較小時，計算中先出現(xiàn)塑性區(qū)貫通，然后出現(xiàn)計算不收斂，而隨著邊坡坡度的增大，計算不收斂和塑性區(qū)貫通可能會同時發(fā)生，甚至計算不收斂時塑性區(qū)仍未貫通，而且在計算不收斂的情況下會出現(xiàn)塑性區(qū)貫通后又不貫通的情況。圖 5～7反映了 ANSYS（D-P）分析45°軟黏土邊坡時塑性區(qū)隨安全系數(shù)變化的關(guān)系。

圖5 K=1.43（計算不收斂，塑形區(qū)恰好貫通）Fig.5 K =1.43(without convergence of calculation, the plastic zone just run through)

圖6 K =1.45（計算不收斂，塑性區(qū)未貫通）Fig.6 K =1.45(without convergence of calculation,and without breakthrough of plastic zone)

圖7 K =1.46（計算不收斂，塑性區(qū)貫通）Fig.7 K =1.46 (without convergence of calculation, but with breakthrough of plastic zone)

4.2.2 工況對滑動面的影響

上面算例定義了3種典型工況：軟黏土、硬黏土、弱膨脹土。以45°坡為例，3種工況下ANSYS（D-P）計算的滑動面如圖8所示。

圖8 45°坡角的各邊坡滑動面Fig.8 Sliding surfaces of slopes with 45 degrees

從圖中可以看出，硬黏土的滑動面最淺，軟黏土和弱膨脹土的滑動面都較深，而且模擬結(jié)果顯示邊坡破壞時軟黏土的位移量最大，弱膨脹土次之，硬黏土最小。改變坡角可以得到同樣的規(guī)律。

不同軟件得到的關(guān)于工況影響的規(guī)律是相同的，但是不同軟件計算得到的滑動面的位置卻是有差異的，以30°軟黏土邊坡為例，ANSYS和MIDAS計算的安全系數(shù)分別為1.70和1.78。而且ANSYS計算得到的滑動面較深切，變換坡度和工況后仍能發(fā)現(xiàn)這一規(guī)律。

圖9 30°軟黏土坡的屈服面Fig.9 The yield faces of 30 degrees slope

5 結(jié) 論

（1）采用 ANSYS(D-P)、MIDAS(D-P)、MIDAS(M-C)、FLAC(M-C)共 4種分析方法，在軟黏土、硬黏土和弱膨脹土 3種工況下的計算表明：FLAC(M-C)最保守，MIDAS(D-P)的結(jié)果最大。

（2）采用M-C屈服準則時當邊坡坡度較?。?5°左右）時MIDAS/GTS的計算結(jié)果小于FLAC，計算中顯示30°及以上邊坡FLAC得到的安全系數(shù)均小于由MIDAS/GTS得到的。

（3）邊坡的滑動面位置與土體材料和坡角大小有關(guān)，軟黏土邊坡的滑動面較深，硬黏土滑動面較淺。

（4）計算軟件的選取也會影響對滑動面深淺的判斷，MIDAS計算得到的滑動面比ANSYS計算得到的滑動面淺。

（5）對軟黏土邊坡分析時，D-P準則的計算結(jié)果和M-C準則的計算結(jié)果偏差較小（12%～25%）；對硬黏土邊坡進行分析時，D-P準則的計算結(jié)果與M-C準則的計算結(jié)果偏差較大（12%～42%）；弱膨脹土工況的情況和軟黏土工況類似。工況條件對分析結(jié)果的影響規(guī)律通過不同軟件分析印證是一致的。

綜上所述，當邊坡土體為軟黏土或弱膨脹土時采用M-C屈服準則和D-P屈服準則分析均可；當邊坡土體為硬黏土時D-P準則計算結(jié)果和M-C準則計算結(jié)果偏大較大，從安全考慮本文推薦硬黏土邊坡盡量采用M-C屈服準則進行分析；當邊坡支護時，如果參照 MIDAS的分析結(jié)果，設計者應考慮MIDAS計算得到的滑動面較其他軟件得到的淺。

[1]四川省國土資源廳. 四川省地質(zhì)環(huán)境公報[J]. 資源與人居環(huán)境, 2011, 8: 53－55.

[2]ZIENKIEWICZ O C, HUMPHESON C, LEWIS R W.Associated and non-associated visco-plasticity in soil mechanics[J]. Geotechnique, 1975, 25(4): 671－689.

[3]WONG F S. Uncertainties in FE modeling of slope stability[J]. Computers & Structures, 1984, 19: 777－791.

[4]連鎮(zhèn)營, 韓國城, 孔憲京. 強度折減有限元法研究開挖邊坡的穩(wěn)定性[J]. 巖土工程學報, 2001, 23(4): 407－411.LIAN Zhen-ying, HAN Guo-cheng, KONG Xian-jing.Stability analysis of excavation by strength reduction FEM[J]. Chinese Journal of Geotechnical Engineering,2001, 23(4): 407－411.

[5]UGAI K. A method of calculation of total factor of safety of slopes by elasto-plastic FEM[J]. Soils and Foundations, 1989, 29(2): 190－195.

[6]MATSUI T, SAN K C. Finite element slope stability analysis by shear strength reduction technique[J]. Soils and Foundations, 1992, 32(1): 59－70.

[7]UGAI K, LESHCHINSHY D. Three-dimensional limit equilibrium and finite element analysis: A comparison of results[J]. Soils and Foundations, 1995, 35(4): 1－7.

[8]GRIFFITHS D V, LANE P A. Slope stability analysis by finite elements[J]. Geotechnique, 1999, 49(3): 387－403.

[9]GRIFFITHS D V, LANE P A. Slope stability analysis by finite elements[J]. Geotechnique, 1999, 49(3): 387－403.

[10]MANZARI M T, NOUR M A. Significance of soil dilatancy in slope stability analysis[J]. Journal of Geotechnical and Geoenvironmental Engineering,ASCE, 2000, 126(1): 75－80.

[11]DUNCAN J M. State of the art: Limit equilibrium and finite-element analysis of slope[J]. Journal of Geotechnical Engineering, ASCE, 1996, 122(7): 577－596.

[12]水利部水利水電規(guī)劃設計總院, 黃河勘測規(guī)劃設計有限公司. SL386－2007水利水電工程邊坡設計規(guī)范實施指南[S]. 北京: 中國水利水電出版社, 2009.

[13]譚長建, 張娟, 董城. ANSYS高級工程應用實例分析與二次開發(fā)[M]. 北京: 電子工業(yè)出版社, 2006.

[14]余天慶, 王勛文, 劉再華. 彈性與塑性力學[M]. 北京:中國建筑工業(yè)出版社, 2009.

[15]DAWSON E M, ROTH W H, DRESCHER A. Slope stability analysis by strength reduction[J]. Geotechnique,1999, 49(6): 835－840.

[16]趙尚毅, 鄭穎人, 孫玉芳. 有限元強度折減法中邊坡失穩(wěn)的判據(jù)探討[J]. 巖土力學, 2005, 26(2): 332－336.ZHAO Shang-yi, ZHENG Ying-ren, SUN Yu-fang. Study of slope failure criterion in strength reduction finite element method[J]. Rock and Soil Mechanics, 2005,26(2): 332－336.

[17]王棟, 年廷凱, 陳煜淼. 邊坡穩(wěn)定有限元分析中的三個問題[J]. 巖土力學, 2007, 28(11): 2310－2313, 2318.WANG Dong, NIAN Ting-kai, CHEN Yu-miao. Three problems in slope stability analyses with finite element method[J]. Rock and Soil Mechanics, 2007, 28(11):2310－2313, 2318.

[18]李紅, 宮必寧, 陳琰. 有限元強度折減法邊坡失穩(wěn)判據(jù)[J]. 水利與建筑工程學報, 2007, 5(1): 79－82.LI Hong, GONG Bi-ning, CHEN Yan. Study of criteria for evaluating stability of slope with FEM based on shear strength reduction methods[J]. Journal of Water Resources and Architectural Engineering, 2007, 5(1):79－82.

[19]MANZARI M T, NOUR M A. Significance of soil dilatancy in slope stability analysis[J]. Journal of Geotechnical and Geoenvironmental Engineering,ACSE, 2000, 123(1): 75－80.

[20]宋二祥. 土工結(jié)構(gòu)安全系數(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.

[21]欒茂田, 武亞軍, 年廷凱. 強度折減有限元中邊坡失穩(wěn)的塑形區(qū)判據(jù)及其應用[J]. 防災減災學報, 2003, 23(3):1－8.LUAN Mao-tian, WU Ya-jun, NIAN Ting-kai. A criterion for evaluating slope stability based on development of plastic zone by shear strength reduction FEM[J]. Journal of Disaster Prevention and Mitigation Engineering,2003, 23(3): 1－8.

[22]周翠英, 劉祚秋, 董立國. 邊坡變形破壞過程的大變形有限元分析[J]. 巖土力學, 2003, 24(4): 644－647, 652.ZHOU Cui-ying, LIU Zuo-qiu, DONG Li-guo. Large deformation FEM analysis of slopes failure[J]. Rock and Soil Mechanics, 2003, 24(4): 644－647, 652.

[23]裴利劍, 屈本寧, 錢閃光. 有限元強度折減法邊坡失穩(wěn)判據(jù)的統(tǒng)一性[J]. 巖土力學, 2010, 31(10): 3337－3341.PEI Li-jian, QU Ben-ning, QIAN Shan-guang.Uniformity of slope instability criteria of strength reduction with FEM[J]. Rock and Soil Mechanics, 2010,31(10): 3337－3341.

[24]劉波, 韓彥輝. FLAC原理、實例與應用指南[M]. 北京:人民交通出版社, 2005.

[25]閆宇, 趙旻. 渠道膨脹土邊坡土體強度參數(shù)取值原則初探[J]. 水利水電工程設計, 2008, 27(2): 37－39, 46.YAN Yu, ZHAO Min. The preliminary discussion on the principle of taking value for strength parameters of canal expansive soil side slope[J]. Design of Water Resources and Hydroelectric Engineering, 2008, 27(2): 37－39,46.