有限元強(qiáng)度折減系數(shù)法計(jì)算土坡穩(wěn)定安全系數(shù)的精度研究

摘要：有限元強(qiáng)度折減系數(shù)法在邊坡穩(wěn)定分析中的應(yīng)用正逐漸受到人們的重視。本文較為全面地分析了土體屈服準(zhǔn)那么的種類、有限元法自身計(jì)算精度以及H(坡高)、β(坡角)、C(粘聚力)、Φ(摩擦角)對(duì)折減系數(shù)法計(jì)算精度的影響，并給出了提高計(jì)算精度的具體措施。通過(guò)對(duì)106個(gè)算例的比擬分析，說(shuō)明折減系數(shù)法所得穩(wěn)定平安系數(shù)比簡(jiǎn)化Bishop法平均高出約5.7%，且離散度極小，這不僅驗(yàn)證了文中所提措施的有效性，也說(shuō)明了將折減系數(shù)法用于分析土質(zhì)邊坡穩(wěn)定問(wèn)題是可行的。關(guān)鍵詞：強(qiáng)度折減系數(shù)邊坡穩(wěn)定屈服準(zhǔn)那么誤差分析自弗倫紐期于1927年提出圓弧滑動(dòng)法以來(lái)，至今已出現(xiàn)數(shù)十種土坡穩(wěn)定分析方法，有極限平衡法、極限分析法、有限元法等。不少研究說(shuō)明，各種方法所得穩(wěn)定平安系數(shù)都比擬接近，可以說(shuō)，這些方法已經(jīng)到達(dá)了相當(dāng)高的精度。近年來(lái)，由于計(jì)算機(jī)技術(shù)的長(zhǎng)足開(kāi)展，基于有限元的折減系數(shù)法在邊坡穩(wěn)定分析中的應(yīng)用備受重視。與極限平衡法相比，它不需要任何假設(shè)，便能夠自動(dòng)地求得任意形狀的臨界滑移面以及對(duì)應(yīng)的最小平安系數(shù)，同時(shí)它還可以真實(shí)的反映坡體失穩(wěn)及塑性區(qū)的開(kāi)展過(guò)程。到目前為止，已有很多學(xué)者對(duì)折減系數(shù)法進(jìn)行了較為深入的研究[1,2,3]，并在一些算例中得到了與極限平衡法十分接近的結(jié)果。但總體說(shuō)來(lái)，此法仍未在工程界得到確認(rèn)和推廣，究其原因在于影響該法計(jì)算精度的因素很多，除了有限元法引入的誤差外，還依賴于所選用的屈服準(zhǔn)那么。

此論文');">論文的目的有兩點(diǎn)：(1)力圖全面分析屈服條件和有限元法本身對(duì)折減系數(shù)法計(jì)算精度的影響，并提出應(yīng)選用何種屈服準(zhǔn)那么以及提高有限元法計(jì)算精度的具體措施；(2)結(jié)合工程實(shí)例，分析對(duì)邊坡穩(wěn)定平安系數(shù)影響最大的4個(gè)主要參數(shù)(H坡高、β坡角、C粘聚力、Φ摩擦角)對(duì)折減系數(shù)法計(jì)算精度的影響。從以往的計(jì)算結(jié)果來(lái)看，嚴(yán)格法(Spencer)所得穩(wěn)定平安系數(shù)比簡(jiǎn)化Bishop法平均高出約2%～3%，而通過(guò)106個(gè)算例的比擬分析，說(shuō)明：折減系數(shù)法所得穩(wěn)定平安系數(shù)比簡(jiǎn)化Bishop法平均高出約5.7%，且誤差離散度極小，可以認(rèn)為是正確的解答[4]。這有力地說(shuō)明了將有限元折減系數(shù)法用于分析土坡穩(wěn)定問(wèn)題是可行的，但必須合理地選用屈服條件以及嚴(yán)格地控制有限元法的計(jì)算精度，同時(shí)也說(shuō)明：有限元折減系數(shù)法所得平安系數(shù)稍微偏高，其原因有待進(jìn)一步研究。

1折減系數(shù)法的根本原理

Bishop等將土坡穩(wěn)定平安系數(shù)F定義為沿整個(gè)滑移面的抗剪強(qiáng)度與實(shí)際抗剪強(qiáng)度之比，工程中廣為采用的各種極限平衡條分法便是以此來(lái)定義坡體穩(wěn)定平安系數(shù)。有限元強(qiáng)度折減系數(shù)法的根本思想與此一致，兩者均可稱之為強(qiáng)度儲(chǔ)藏平安度。因后者無(wú)法直接用公式計(jì)算平安系數(shù)，而需根據(jù)某種破壞判據(jù)來(lái)判定系統(tǒng)是否進(jìn)入極限平衡狀態(tài)，這樣不可防止地會(huì)帶來(lái)一定的人為誤差。盡管如此，仍開(kāi)展了一些切實(shí)可行的平衡判據(jù)，如：限定求解迭代次數(shù)，當(dāng)超過(guò)限值仍未收斂那么認(rèn)為破壞發(fā)生；或限定節(jié)點(diǎn)不平衡力與外荷載的比值大??；或利用可視化技術(shù)，當(dāng)廣義剪應(yīng)變等值線自坡角與坡頂貫穿那么定義坡體破壞[3]。文中平衡判據(jù)?。寒?dāng)節(jié)點(diǎn)不平衡力與外荷載的比值大于10-3時(shí)便認(rèn)為坡體破壞。

有限元折減系數(shù)法的根本原理是將土體參數(shù)C、Φ值同時(shí)除以一個(gè)折減系數(shù)Ftrial，得到一組新的C′、Φ′值，然后作為新的材料參數(shù)帶入有限元進(jìn)行試算，當(dāng)計(jì)算正好收斂時(shí)，也即Ftrial再稍大一些(數(shù)量級(jí)一般為10-3)，計(jì)算便不收斂，對(duì)應(yīng)的Ftrial被稱為坡體的最小平安系數(shù)，此時(shí)土體到達(dá)臨界狀態(tài)，發(fā)生剪切破壞，具體計(jì)算步驟可參考文獻(xiàn)[2]，文中如無(wú)特別說(shuō)明，計(jì)算結(jié)果均指到達(dá)臨界狀態(tài)時(shí)的折減系數(shù)。(1)(2)2屈服準(zhǔn)那么的影響

用折減系數(shù)法求解實(shí)際邊坡穩(wěn)定問(wèn)題時(shí)，通常將土體假設(shè)成理想彈塑性體，其中本構(gòu)模型常選用摩爾-庫(kù)侖準(zhǔn)那么(M-C)、Drucker-Prager準(zhǔn)那么以及摩爾-庫(kù)侖等面積圓[5]準(zhǔn)那么。

摩爾-庫(kù)侖準(zhǔn)那么可用不變量I1,J2,θσ表述成如下形式：(3)Drucker-prager準(zhǔn)那么：(4)式中：I1為應(yīng)力張量第一不變量；J2為應(yīng)力偏量第二不變量；θσ是應(yīng)力洛德角。

圖1各屈服準(zhǔn)那么在π平面上的曲線M-C準(zhǔn)那么較為可靠，它的缺點(diǎn)在于三維應(yīng)力空間中的屈服面存在尖頂和棱角的不連續(xù)點(diǎn)，導(dǎo)致數(shù)值計(jì)算不收斂，所以有時(shí)也采用抹圓了的M-C修正準(zhǔn)那么[6]，它是用光滑連續(xù)曲線來(lái)逼進(jìn)摩爾-庫(kù)侖準(zhǔn)那么，此法雖然方便了數(shù)值計(jì)算，但不可防止地會(huì)引入一定的誤差；而D-P準(zhǔn)那么在偏平面上是一個(gè)圓，更適合數(shù)值計(jì)算。通常取M-C準(zhǔn)那么的外角點(diǎn)外接圓、內(nèi)角點(diǎn)外接圓或其內(nèi)切圓作為屈服準(zhǔn)那么，以利數(shù)值計(jì)算。各準(zhǔn)那么的參數(shù)換算關(guān)系見(jiàn)表1。

由徐干成、鄭穎人(1990)提出的摩爾庫(kù)侖等效面積圓準(zhǔn)那么[5]實(shí)際上是將M-C準(zhǔn)那么轉(zhuǎn)化成近似等效的D-P準(zhǔn)那么形式。該準(zhǔn)那么要求偏平面上的摩爾-庫(kù)侖不等邊六角形與D-P圓面積相等。計(jì)算說(shuō)明它與摩爾-庫(kù)侖準(zhǔn)那么十分接近。

見(jiàn)圖1，r1為外角外接圓半徑；r2為內(nèi)角外接圓半徑；r3為內(nèi)切圓半徑；摩爾-庫(kù)侖準(zhǔn)那么構(gòu)成的六角形面積為(5)對(duì)半徑為r的圓面積S=πr2，令S=Smorl得(6)(7)式(7)與式(4)對(duì)應(yīng)項(xiàng)相等，可得(8)

表1各準(zhǔn)那么參數(shù)換算編號(hào)準(zhǔn)那么種類αφkDP1外角點(diǎn)外接D-P圓DP2內(nèi)角點(diǎn)外接D-P圓DP3內(nèi)切D-P圓DP4等面積D-P圓注：表中αφ、k是與D-P有關(guān)的材料參數(shù)。表2不同屈服準(zhǔn)那么所得最小平安系數(shù)

φ/°

0.110253545DP10.5251.0441.7692.2543.051DP20.5250.9301.3321.5301.887DP30.4540.8481.2791.4991.870DP40.4770.8961.3961.6892.182簡(jiǎn)化Bishop法0.4940.8461.3161.6232.073(DP1-Bishop)/Bishop0.0630.2340.3440.3550.472(DP2-Bishop)/Bishop0.0630.0990.012-0.080-0.090(DP3-Bishop)/Bishop-0.0810.002-0.028-0.099-0.098(DP4-Bishop)/Bishop-0.0340.0590.0610.0410.053注：H=20mm；β=45°；C=42kPa。算例分析說(shuō)明(表2、圖2)：DP4準(zhǔn)那么與簡(jiǎn)化Bishop法所得穩(wěn)定平安系數(shù)最為接近。對(duì)有效算例(Φ≠0)的誤差進(jìn)行統(tǒng)計(jì)分析可知，中選用DP4準(zhǔn)那么時(shí)，誤差的平均值為5.7%，且離散度很小(圖3)。而DP1的平均誤差為29.5%，同時(shí)采用DP2、DP3準(zhǔn)那么所得計(jì)算結(jié)果的離散度非常大，均不可用。因此在數(shù)值分析中可用DP4準(zhǔn)那么代替摩爾-庫(kù)侖準(zhǔn)那么。

圖2Φ～折減系數(shù)曲線

圖3DP4準(zhǔn)那么的計(jì)算誤差3不同流動(dòng)法那么的影響

有限元計(jì)算中，采用關(guān)聯(lián)還是非關(guān)聯(lián)流動(dòng)法那么，取決于ψ值(剪脹角)：ψ=φ，為關(guān)聯(lián)流動(dòng)法那么；ψ≠0，為非關(guān)聯(lián)流動(dòng)法那么?？傮w說(shuō)來(lái)，采用非關(guān)聯(lián)流動(dòng)法那么所得破壞荷載比同一類型材料而采用關(guān)聯(lián)流動(dòng)法那么所得破壞荷載小，如忽略剪脹角(ψ=0)，將會(huì)得到較為保守的結(jié)果。值得注意的是：當(dāng)ψ=0時(shí)，正好與鄭穎人等提出的廣義塑性力學(xué)理論相符[7]，這時(shí)對(duì)應(yīng)的塑性勢(shì)面與q軸垂直。表3不同流動(dòng)法那么的影響

φ=10°φ=17°φ=25°非關(guān)聯(lián)0.8711.1051.363關(guān)聯(lián)0.8871.1371.425相對(duì)誤差0.0180.0290.045

β=45°；C=40kPa；H=20m；DP4準(zhǔn)那么。表4網(wǎng)格疏密對(duì)計(jì)算結(jié)果的影響

節(jié)點(diǎn)數(shù)

57711112250DP40.6610.6180.593簡(jiǎn)化Dishop法0.5830.5830.583(DP4-Bishop)/Bishop0.1340.0600.017

注：H=20m；β=45°；φ=45°；c=10000Pa。筆者對(duì)采用不同流動(dòng)法那么的算例進(jìn)行了初步分析，表3的計(jì)算結(jié)果說(shuō)明：對(duì)同一邊坡，不管采用關(guān)聯(lián)流動(dòng)法那么還是非關(guān)聯(lián)流動(dòng)法那么，計(jì)算結(jié)果相差不大。這是因?yàn)樗鼈冎慌c坡體的體積變形有關(guān)，而在邊坡穩(wěn)定分析中，坡體常常為無(wú)約束天然坡體，體積變形對(duì)坡體穩(wěn)定影響并不明顯。然而，從破壞時(shí)位移大小及塑性區(qū)的分布來(lái)看，還是會(huì)有一些差異，有時(shí)并不能簡(jiǎn)單的忽略這種差異[8]。文中所有的算例均取ψ=0，即滿足非關(guān)聯(lián)流動(dòng)法那么，算例結(jié)果顯示出較好的精度。

4有限元法引入的誤差

如前所述，本構(gòu)模型的選擇合理與否會(huì)對(duì)有限元折減系數(shù)法的計(jì)算精度造成較大影響，除此之外，有限元法本身也是誤差的主要來(lái)源之一。

4.1網(wǎng)格的疏密網(wǎng)格疏密對(duì)單元精度的影響甚至大于單元類型的影響，對(duì)于精度較低的單元，可通過(guò)加密網(wǎng)格來(lái)到達(dá)較高的精度。表4列出了不同疏密的網(wǎng)格對(duì)計(jì)算結(jié)果的影響，由表4可知，對(duì)于折減系數(shù)法，有限元網(wǎng)格不能太稀，否那么結(jié)果將不可用。通過(guò)大量算例證實(shí)，對(duì)于4節(jié)點(diǎn)矩形單元，當(dāng)單元密度到達(dá)每10m2不少于3個(gè)節(jié)點(diǎn)時(shí)，計(jì)算精度較為理想，如果再增加節(jié)點(diǎn)，計(jì)算精度應(yīng)還能提高，但此時(shí)消耗的機(jī)時(shí)也將成倍增長(zhǎng)。

4.2邊界范圍邊界范圍的大小在有限元法中對(duì)計(jì)算結(jié)果的影響比在傳統(tǒng)極限平衡法中表現(xiàn)的更為敏感，在極限平衡法中只要所求滑移面在邊界之內(nèi)就不會(huì)對(duì)計(jì)算結(jié)果有影響，平安系數(shù)只與劃分的土條有關(guān)，而與土條外的區(qū)域無(wú)關(guān)，有限元法那么不然，邊界的大小直接影響到應(yīng)力-應(yīng)變的分布。表5邊界條件對(duì)折減系數(shù)的影響

相對(duì)邊距比

00.51.01.52.02.53.0L/H

R/H

B/H1.129

1.097

1.1061.124

1.078

1.1171.124

1.121

1.1201.120

1.122

1.1311.122

1.122

1.1241.121

1.120

1.1321.129

1.123

1.131

注：表中所有模型均選用DP4屈服準(zhǔn)那么。L為坡腳到左端邊界的距離〔左

邊距〕；R為坡頂?shù)接叶诉吔绲木嚯x〔右邊距〕；B為邊坡到底端邊

界的距離〔底邊

距〕；H為坡高。為了得到能使計(jì)算結(jié)果趨于穩(wěn)定的邊界范圍，分別對(duì)左端、右端、底端三條邊界范圍的取值大小進(jìn)行了分析(表5、圖4)，計(jì)算時(shí)，令3個(gè)邊距中的1個(gè)變化，其余2個(gè)不變。由圖4可知：左邊界對(duì)計(jì)算結(jié)果的影響最不敏感，不同的取值相差不到1%，底端邊界次之，最大相差在1%左右，右端邊界對(duì)計(jì)算精度的影響最大，到達(dá)5%。經(jīng)比照分析得：當(dāng)坡角到左端邊界的距離為坡高的1.5倍，坡頂?shù)接叶诉吔绲木嚯x為坡高的2.5倍，且上下邊界總高不低于2倍坡高時(shí)，計(jì)算精度最為理想。

5可行性研究5.1模型的建立及研究方案某一處于施工階段的土質(zhì)邊坡，不考慮孔隙水壓力的影響，土坡天然容重γ=25kN/m3。坡體底邊界為固定約束，左右邊界為水平約束，其它邊界為自由端。計(jì)算程序采用商業(yè)有限元軟件Ansys，計(jì)算單元采用平面4節(jié)點(diǎn)矩形單元，當(dāng)坡高保持不變時(shí)(20m)，節(jié)點(diǎn)1111個(gè)，對(duì)于坡高H與坡度β變化的情況，為保證足夠的精度，必須重新劃分單元。同時(shí)還應(yīng)在臨界滑移面可能區(qū)域?qū)卧W(wǎng)格進(jìn)行局部加密，見(jiàn)圖5。有限元計(jì)算所需參數(shù)共6個(gè)，見(jiàn)表6，其中所有算例均取ψ=0。為防止個(gè)別特例在計(jì)算結(jié)果上存在的巧合，文中對(duì)影響坡體穩(wěn)定平安系數(shù)的4個(gè)主要參數(shù)(坡高H、坡角β、粘聚力C、摩擦角Φ)進(jìn)行了比照分析，計(jì)算方案如下：保持4個(gè)參數(shù)中的3個(gè)為常量，只取1個(gè)參數(shù)為變量，共4組計(jì)算方案，計(jì)算結(jié)果見(jiàn)表2、表7、表8、表9。

圖4邊界距離與坡高比～折減系數(shù)曲線

圖5有限元單元網(wǎng)格劃分

表6有限元計(jì)算參數(shù)Drucker-Prager準(zhǔn)那么彈性模量E/kPa1000帕松比ν0.3土體密度γ/(kN/m3)25磨擦角?/°變量粘聚力C/Pa變量表7c為變量時(shí)的最小平安系數(shù)(節(jié)點(diǎn)數(shù)1111個(gè))

C/kPa

0.120406090DP40.3040.7931.1011.3791.781簡(jiǎn)化Bishop法0.2540.7521.0361.3021.685(DP4-Bishop)/Bishop0.1970.0550.0630.0590.057

注：H=20m；β=45°；?=17°。

表8H為變量時(shí)的最小平安系數(shù)(節(jié)點(diǎn)數(shù)≥1190個(gè))

H/m

1020304050DP41.7331.1280.9230.8200.735簡(jiǎn)化Bishop法1.6121.0640.8670.7640.698(DP4-Bishop)/Bishop0.0750.0600.0650.0730.053

注：β=45°；c=42kPa；?=17°。表9β為變量時(shí)的最小平安系數(shù)(節(jié)點(diǎn)數(shù)≥1210個(gè))

坡角β/°

3035404550DP41.4551.3231.2141.1281.044簡(jiǎn)化Bishop法1.3981.2691.1561.0640.987(DP4-Bishop)/Bishop0.0410.0430.0500.0600.058

注：H=20m；c=42kPa；?=17°。表2、表7、表8、表9列出了各種情況下的最小平安系數(shù)。作為比照，文中還給出了極限平衡法(簡(jiǎn)化Bishop法)的計(jì)算結(jié)果，考慮到算例坡體土質(zhì)均勻，可采用圓弧滑移面，在此情況下簡(jiǎn)化Bishop法已具有足夠的精度[4]。誤差定義詳見(jiàn)各表。

5.2結(jié)果分析屈服準(zhǔn)那么對(duì)計(jì)算精度的影響很明顯。在相同網(wǎng)格密度下，DP4的計(jì)算精度明顯好于DP1、DP2、DP3(圖2)。因?yàn)楦鞣N方案所顯示的規(guī)律相同，所以文中只給出了方案1的詳細(xì)解答(表2)，而其它方案只給出DP4的計(jì)算結(jié)果。

Φ值的大小對(duì)計(jì)算精度的影響是明顯的。如圖2所示，Φ值增大，誤差也呈增大的趨勢(shì)。在文中所給出的4個(gè)屈服準(zhǔn)那么中，DP4精度最高，其與簡(jiǎn)化Bishop法最接近，平均誤差為5.7%。

C、β、H對(duì)計(jì)算精度的影響不明顯，圖6分別給出了不同C、β、H值對(duì)計(jì)算結(jié)果的影響，它們對(duì)結(jié)果的影響約在1%左右。

圖6C、H、β～折減系數(shù)曲線值得注意的是：當(dāng)Φ、C分別為零時(shí)，DP4的誤差較大，這是因D-P類本構(gòu)只適用于摩擦型材料，當(dāng)?=0時(shí)，計(jì)算是不收斂的，在這里筆者代入?=0.1以近似?=0；當(dāng)C為零時(shí)，處理方法同Φ。因計(jì)算誤差較大，所以此時(shí)折減系數(shù)法將不再適用。

需要指出的是，極限平衡法采用的圓弧滑動(dòng)面也會(huì)對(duì)計(jì)算精度造成一定的影響，特別是當(dāng)β很大時(shí)，如果采用任意滑動(dòng)面的Spencer法，那么計(jì)算精度還能提高。

6結(jié)論

(1)通過(guò)4組計(jì)算方案共計(jì)106個(gè)算例的比擬分析，說(shuō)明：折減系數(shù)法所得穩(wěn)定平安系數(shù)比簡(jiǎn)化Bishop法平均高出約5.7%，且誤差離散度極小。這有力地說(shuō)明了將有限元折減系數(shù)法用于分析土坡穩(wěn)定問(wèn)題是可行的，但必須合理地選用屈服條件并嚴(yán)格地控制有限元法的計(jì)算精度。(2)邊界范圍的大小在有限元法中對(duì)計(jì)算結(jié)果的影響比在傳統(tǒng)極限平衡法中表現(xiàn)的更為敏感，當(dāng)坡角到左端邊界的距離為坡高的1.5倍，坡頂?shù)接叶诉吔绲木嚯x為坡高的2.5倍，且上下邊界總高不低于2倍坡高時(shí)，計(jì)算精度最為理想。(3)Φ值的大小對(duì)有限元折減系數(shù)法的計(jì)算精度有較大影響，且隨Φ值增大誤差隨之增大，增大幅度因準(zhǔn)那么類型不同而不同，其中準(zhǔn)那么DP4計(jì)算誤差隨Φ值增大影響相對(duì)較小。C、β、H對(duì)計(jì)算精度的影響不明顯，約在1%左右。(4)有限元法的優(yōu)點(diǎn)不僅僅在于求出折減系數(shù)，如果此法可行，那么對(duì)于具有復(fù)雜地貌、地質(zhì)的邊坡那么可自動(dòng)求出任意形狀的臨界滑移面，并能模擬出土坡失穩(wěn)及施工開(kāi)挖的自然過(guò)程，這是傳統(tǒng)極限平衡法無(wú)法做到的。因此本文對(duì)該法的可行性分析是重要且必要的。目前的工作是個(gè)根底，如何在有限元折減系數(shù)法中考慮多種土層邊坡、孔隙水的影響、非自重外荷以及巖質(zhì)邊坡中存在的大量節(jié)理等仍需做深入研究。參考文獻(xiàn)：

[1]JiangGL，MagnanJP.Stabilityanalysisofembankments:comparisonoflimitanalysiswithmethodsofslices[J].Geotechnique1997,47(4):857-872.

[2]GriffithsDV，LanePA.Slopestabilityanalysisbyfiniteelements[J].Geotechnique,1999,49(3):387-403.

[3]連鎮(zhèn)營(yíng)，韓國(guó)城，孔憲京.強(qiáng)度折減有限元法開(kāi)挖邊坡的穩(wěn)定性[J].巖土工程學(xué)報(bào)，2001，23(4)：407-411.

[4]龔曉南.土工計(jì)算機(jī)分析[M].北京：中國(guó)建筑工業(yè)出版社，2000.

[5]徐干成，鄭穎人.巖土工程中屈服準(zhǔn)那么應(yīng)用的研究[J].巖土工程學(xué)報(bào)，1990，(2)：93-99.

[6]MenétreyPh,WillamK.Atriaxialfailurecriterionforconcreteanditsge