1、有限元強(qiáng)度折減法1 背景1974年，Smith & Hobbs1使用有限元方法分析了u=0條件下的邊坡穩(wěn)定性并與Taylar2的結(jié)果進(jìn)行比照，得到了很好的一致性；1975年，Zienkiewicz等3考慮c、進(jìn)行有限元邊坡穩(wěn)定性分析，其結(jié)果與圓弧滑面解有較好吻合；1980年Griffiths4驗(yàn)證了一系列具有不同材料特性和形狀的邊坡穩(wěn)定性并通過與Bishop& Morgenstern5的結(jié)果進(jìn)行了比照確定了數(shù)據(jù)的可靠性；此后也有研究證實(shí)了利用有限元方法進(jìn)行邊坡穩(wěn)定性分析的可靠性6,7,8,9；在文獻(xiàn)9中，引入一些案例證明了有限元強(qiáng)度折減法的準(zhǔn)確性，并證明了有限元強(qiáng)度折減法在分析

2、非均質(zhì)邊坡時(shí)相對于傳統(tǒng)方法的優(yōu)越性。2001年，鄭穎人等10把有限元強(qiáng)度折減法引入國內(nèi)，并對此進(jìn)行了后續(xù)研究11,12,13,14。相較于一些傳統(tǒng)的邊坡穩(wěn)定型分析方法，有限元強(qiáng)度折減法有以下幾個(gè)優(yōu)點(diǎn)9：(1) 不必假設(shè)滑面的位置和形狀，當(dāng)土體自身強(qiáng)度缺乏以抵抗剪應(yīng)力時(shí)土體失穩(wěn)會(huì)自然發(fā)生。(2) 由于有限元強(qiáng)度折減法中沒有條分的概念，因此也不必假設(shè)條間力，在整體失穩(wěn)之前土體都處于整體穩(wěn)定狀態(tài)。(3) 使用有限元方法能夠查看破壞過程。2 有限元強(qiáng)度系數(shù)折減法1.模型參數(shù)邊坡模型主要包括六個(gè)參數(shù)，分別是：膨脹角、內(nèi)摩擦角、黏聚力c、彈性模量E、泊松比、重度。膨脹角影響土體屈服后的體積變形，假設(shè)<

3、;0，那么土體屈服后體積減小，假設(shè)>0那么體積增大，=0那么體積不變。=的情況被稱之為關(guān)聯(lián)流動(dòng)法那么，但是此時(shí)值通常高于實(shí)驗(yàn)觀測值，特別是在側(cè)限條件下會(huì)提高土的承載力預(yù)測值。邊坡穩(wěn)定型問題通常是處于無側(cè)限條件下，此時(shí)膨脹角的選取不再重要9，因此文獻(xiàn)9選取=0條件下的非關(guān)聯(lián)流動(dòng)法那么，并且通過案例分析可以得出此膨脹角的選取可以得出準(zhǔn)確的平安系數(shù)以及滑動(dòng)面。c和指Mohr-Coulomb準(zhǔn)那么中邊坡土體的有效黏聚力和內(nèi)摩擦角；E和是土體材料的彈性參數(shù)，這兩個(gè)參數(shù)對土體穩(wěn)定性分析的影響較小；是土體的重度。應(yīng)用有限元方法進(jìn)行邊坡穩(wěn)定性分析中最重要的三個(gè)參數(shù)是c、 、和。2.屈服條件(1)Mohr

4、-Coulomb準(zhǔn)那么Mohr-Coulomb準(zhǔn)那么用大小主應(yīng)力表示如式(1)所示：1'-3'2=1'+3'2sin'-c cos' (1)其中， 1'、3'分別指土中一點(diǎn)的大小主應(yīng)力。在主應(yīng)力空間中，如果不考慮1、2、3之間的大小關(guān)系，屈服面是一個(gè)不等角六棱錐，在平面上是一個(gè)等邊不等角六邊形。(2)D-P準(zhǔn)那么D-P準(zhǔn)那么可以寫成式(2)形式：-I1+J2=kf (2)其中I1為第一應(yīng)力不變量、J2為第二偏應(yīng)力不變量，和kf為試驗(yàn)常數(shù)。在主應(yīng)力空間中其屈服面為一個(gè)圓錐，在平面上是一個(gè)圓形。3)D-P準(zhǔn)那么轉(zhuǎn)換為Mohr-Coul

5、omb準(zhǔn)那么首先引入?yún)?shù)b，如式(3)所示：b=2-31-3 (3)(4)那么，I1和J2分別可轉(zhuǎn)化為式(4)：I1=3(1+3)2+(b-12)(1-3)J2=3(1-3)其中=1-b+b2將其帶入(2)，得式(5)：1-3=1.52-b+31+3+kf2-b+3 (5)(6)與式(1)比照可知兩個(gè)準(zhǔn)那么之間的轉(zhuǎn)換關(guān)系如式(6)所示：sinb=1.52-b+32cbcosb=kf2-b+3因此，當(dāng)b=0時(shí)，即外角點(diǎn)外接DP圓的兩個(gè)試驗(yàn)常數(shù)分別如式(7)所示，當(dāng)b=1時(shí)，即內(nèi)角點(diǎn)外接DP圓的兩個(gè)試驗(yàn)常數(shù)分別如式(8)所示。=2sin3(3-sin)，kf=6ccos3(3-sin) (7)=2

6、sin3(3+sin)，kf=6ccos3(3+sin) (8)3.平安系數(shù)的定義(1)Mohr-Coulomb準(zhǔn)那么中的平安系數(shù)1955年，Bishop15首先在邊坡穩(wěn)定性分析中提出了抗剪強(qiáng)度折減的概念，在有限元強(qiáng)度折減法中通過將坡體的強(qiáng)度參數(shù)：黏聚力c和內(nèi)摩擦角同時(shí)除一個(gè)折減系數(shù)Ft，得到一組新的c和值，作為一個(gè)新的強(qiáng)度參數(shù)輸入進(jìn)行試算，當(dāng)計(jì)算不收斂時(shí)，對應(yīng)的Ft即為所求的平安系數(shù)，此時(shí)坡體到達(dá)極限狀態(tài)，發(fā)生剪切破壞。c=c/Ft=arctan(tan/Ft)(2)D-P(Drucker-Prager)準(zhǔn)那么中的平安系數(shù)取Ft為D-P準(zhǔn)那么中的強(qiáng)度折減系數(shù)，那么D-P準(zhǔn)那么可以表示為式(9

7、)，-FtI1+J2=kfFt (9)(3)不同屈服條件下平安系數(shù)轉(zhuǎn)換13(10)首先引入Mohr-Coulomb等面積圓屈服準(zhǔn)那么，在平面上，其屈服面是一個(gè)圓，并且面積與Mohr-Coulomb準(zhǔn)那么的不等角六邊形相等，Mohr-Coulomb等面積圓屈服準(zhǔn)那么中的試驗(yàn)參數(shù)如式(10)所示：-=sin3(3cos-sinsin)kf=3ccos3cos-sinsin式中=arcsin-23Asin+49A2sin2-4sin23+1A23-1122sin23-1，A=9-sin263簡稱外接圓屈服準(zhǔn)那么為DP1準(zhǔn)那么，其試驗(yàn)常數(shù)分別為1，kf1；Mohr-Coulomb等面積圓屈服準(zhǔn)那么為D

8、P2準(zhǔn)那么，其試驗(yàn)常數(shù)分別為2，kf2。把DP1準(zhǔn)那么表示為f1=J2=1I1+kf1，DP2準(zhǔn)那么可表示為f2=J2=2I1+kf2。令=12=kf1kf2=f()，f1=1I1+kf1=2I1+kf2，所以f1f2=2I1+kf22I1+kf2=f()。由此可知，是的函數(shù)，當(dāng)取不同值時(shí)可以得到不同的值如表1所列：表 1 不同內(nèi)摩擦角時(shí)的值4.失穩(wěn)判據(jù)目前兩個(gè)比較主流的失穩(wěn)判據(jù)分別是有限元計(jì)算中力不平衡和位移的不收斂以及廣義塑性應(yīng)變或者等效塑性應(yīng)變從坡腳到坡頂貫穿。Griffiths9和鄭穎人11,12,13,14都使用計(jì)算不收斂作為失穩(wěn)判據(jù)。Griffiths9提出，當(dāng)在用戶定義的最大迭代

9、數(shù)目下計(jì)算仍不收斂時(shí)，那么沒有任何一種應(yīng)力分布方式可以同時(shí)滿足Mohr-Coulomb準(zhǔn)那么以及整體穩(wěn)定，這種情況可看做邊坡失穩(wěn)判據(jù)。邊坡失穩(wěn)與數(shù)值計(jì)算不收斂同時(shí)發(fā)生，并伴隨著極大的節(jié)點(diǎn)位移，并以1000作為最大的迭代步數(shù)。鄭穎人14提出，有限元的計(jì)算迭代過程就是尋找外力和內(nèi)力到達(dá)平衡狀態(tài)的過程，整個(gè)迭代過程直到一個(gè)適宜的收斂標(biāo)準(zhǔn)得到滿足才停止。可見，如果邊坡失穩(wěn)破壞，滑面上將產(chǎn)生沒有限制的塑性變形，有限元程序無法從有限元方程組中找到一個(gè)既能滿足靜力平衡又能滿足應(yīng)力-應(yīng)變關(guān)系和強(qiáng)度準(zhǔn)那么的解，此時(shí)不管是從力的收斂標(biāo)準(zhǔn)，還是從位移的收斂標(biāo)準(zhǔn)來判斷有限元計(jì)算都不收斂。3 案例分析例一，不含地基的均

10、質(zhì)邊坡9該邊坡如圖1所示，有限元程序采用Mohr-Coulomb失效準(zhǔn)那么，建立平面應(yīng)變條件下八節(jié)點(diǎn)四邊形單元減縮積分計(jì)算模型，其強(qiáng)度參數(shù)為=20°，c/h=0.05。邊坡坡度為26.57°(2:1)，坡底水平，其邊界條件為坡底約束豎直方向位移與水平方向位移，左側(cè)約束水平方向位移，其余面為自由面。施加重力荷載后使平安系數(shù)從0.8到1.4逐步變化直至計(jì)算不收斂圖 1 不含地基的均質(zhì)邊坡每一個(gè)平安系數(shù)對應(yīng)的迭代次數(shù)如表2所列，當(dāng)真正的平安系數(shù)接近時(shí)需要更多的迭代次數(shù)。表 2 例一計(jì)算結(jié)果當(dāng)平安系數(shù)為1.4時(shí)，無量綱位移Emax/H2突變，并且此時(shí)計(jì)算無法收斂，在此情況下有限元計(jì)

11、算結(jié)果與Bishop & Morgenstern5給出的結(jié)果吻合良好，如圖2所示。圖 2 平安系數(shù)與無量綱位移邊坡失穩(wěn)時(shí)(FOS=1.4)節(jié)點(diǎn)位移矢量和網(wǎng)格變形如圖3(a)和圖3(b)所示，由此可得到邊坡的潛在滑動(dòng)面。圖 3平安系數(shù)為1.4計(jì)算不收斂時(shí)邊坡變形(a)節(jié)點(diǎn)位移矢量(b)網(wǎng)格變形例二，有軟弱層的不排水黏性土邊坡在本案例中，使用Tresca準(zhǔn)那么(u=0)進(jìn)行總應(yīng)力分析。邊坡幾何形狀如圖4所示，地基厚度與邊坡高度相同，該邊坡有一個(gè)軟弱層，在有限元計(jì)算中，令其抗剪強(qiáng)度(Cu2)在一定范圍內(nèi)變化但其周圍土體抗剪強(qiáng)度保持Cu1/H=0.25不變。利用有限元方法計(jì)算該邊坡的平安系數(shù)結(jié)

12、果如圖5所示，對于均質(zhì)邊坡情況，Cu2/Cu1=1，有限元計(jì)算結(jié)果與Taylor2的結(jié)論很接近，隨著軟弱層的強(qiáng)度逐漸減小，在Cu2/Cu10.6時(shí)，結(jié)果發(fā)生了明顯的變化。分別假定圓弧滑面和穿過軟弱面的三段線滑面并利用Janbu法計(jì)算平安系數(shù)，可見在Cu2/Cu10.6處也發(fā)生了滑動(dòng)機(jī)制的轉(zhuǎn)換，當(dāng)Cu2/Cu1>0.6時(shí)，潛在滑面形狀為圓弧，當(dāng)Cu2/Cu1<0.6時(shí)，潛在滑面為結(jié)構(gòu)軟弱面。圖6更加清晰的展示了這一現(xiàn)象，圖6(a)為均質(zhì)邊坡(Cu2/Cu1=1)時(shí)的潛在滑面，可見此時(shí)的滑面形狀為圓弧滑面，與Taylor2的預(yù)測相同；圖6(c)為軟弱層強(qiáng)度只有其周圍土體20%( Cu2

13、/Cu1=0.2)時(shí)的潛在滑面，此時(shí)潛在滑面沿軟弱層開展；圖6(b)為軟弱層強(qiáng)度只有其周圍土體60%( Cu2/Cu1=0.6)時(shí)的潛在滑面，此時(shí)圓弧滑面和沿軟弱層的三段線式滑面都有可能開展，至少存在兩種明顯的滑動(dòng)機(jī)制。圖 4有軟弱層的不排水黏性邊坡圖 5 不同軟弱層強(qiáng)度時(shí)的平安系數(shù)圖 6 不同軟弱層強(qiáng)度下的網(wǎng)格變形 (a) Cu2/Cu1=1.0 (b) Cu2/Cu1=0.6 (c) Cu2/Cu1=0.2例三，不同坡度邊坡平安系數(shù)計(jì)算13，驗(yàn)證Mohr-Coulomb等面積圓屈服準(zhǔn)那么均質(zhì)邊坡,坡高H=20m,土容重=25kN/m3,黏聚力c=42kPa,內(nèi)摩擦角=17°,求坡

14、角分別為30°,35°,40°,45°,50°時(shí)邊坡的平安系數(shù)。計(jì)算結(jié)果如表3所列。表 3 平安系數(shù)計(jì)算結(jié)果從表中計(jì)算結(jié)果可以看出,采用外接圓屈服準(zhǔn)那么計(jì)算的平安系數(shù)比傳統(tǒng)的方法大許多,采用莫爾-庫侖等面積圓屈服準(zhǔn)那么計(jì)算的結(jié)果與傳統(tǒng)極限平衡方法(Spencer法)計(jì)算的結(jié)果十分接近,說明采用莫爾-庫侖等面積圓屈服準(zhǔn)那么來代替莫爾-庫侖不等角六邊形屈服準(zhǔn)那么是可行的,這樣使計(jì)算大為方便。而采用外接圓屈服準(zhǔn)那么計(jì)算的平安系數(shù)要比莫爾-庫侖等面積圓屈服準(zhǔn)那么計(jì)算的結(jié)果大(1.21)倍。例四，存在兩組節(jié)理面的巖質(zhì)邊坡穩(wěn)定性分析12如圖7所示，巖體中存

15、在兩組方向不同的軟弱結(jié)構(gòu)面，貫穿率100%，第一組軟弱結(jié)構(gòu)面傾角為30°，平均間距10m；第二組軟弱結(jié)構(gòu)面傾角75°，平均間距10m。巖體重度為25kN/m3，彈性模量1×1010Pa，泊松比0.2，黏聚力1MPa，內(nèi)摩擦角38°，兩組節(jié)理參數(shù)相同，重度為17kN/m3，彈性模量1×107Pa，泊松比0.3，黏聚力0.12MPa，內(nèi)摩擦角24°。按照二維平面應(yīng)變問題建立有限元模型，按照連續(xù)介質(zhì)處理。通過有限元強(qiáng)度折減，求得坡體破壞時(shí)的運(yùn)動(dòng)矢量如圖8所示，滑動(dòng)面如圖9(a)所示，它是最先貫穿的塑性區(qū)，塑性區(qū)貫穿并不等于破壞，當(dāng)塑性區(qū)貫穿

16、后繼續(xù)開展到一定程度，巖體發(fā)生整體破壞，同時(shí)出現(xiàn)第二條貫穿的塑性面，如圖9(b)所示。求得的穩(wěn)定平安系數(shù)如表4所列，其中，極限平衡方法計(jì)算結(jié)果是根據(jù)最先貫穿的那一條滑動(dòng)面求得的。圖 Error! Main Document Only. 巖質(zhì)邊坡節(jié)理 圖 Error! Main Document Only. 極限狀態(tài)時(shí)的塑性區(qū)圖 Error! Main Document Only.坡體破壞時(shí)的運(yùn)動(dòng)矢量圖表 4 例四計(jì)算結(jié)果例五，存在接觸問題的邊坡穩(wěn)定性分析12當(dāng)邊坡中存在如圖10所示的硬性結(jié)構(gòu)面時(shí)，不能按照例四中軟弱結(jié)構(gòu)面的方法進(jìn)行處理，可以采用接觸單元來模擬硬性接觸面的不連續(xù)性。按照Mohr-

17、Coulomb定律來定義接觸面上的摩擦行為，如式11所示，那么其接觸面上的平安系數(shù)定義如式11所示。=c+tan，0 Ft=c/c=tan/tan (11)圖 10 無充填的硬性結(jié)構(gòu)面圖11所示為兩個(gè)直線滑面組成的折線型滑體ABMCD。巖體重度=20kN/m3，彈性模量E=109Pa。滑塊ABCD面積433m2，滑面AB=20m，傾角為15°，AD=25m，DC=19.32m，BC=19.82m；滑塊BCM面積196.5m2，滑面BM=28.03m，傾角為45°，CM=19.82m。CM面上施加有線性變化的面荷載，PM=400kPa，PC=0。圖 11 折線型平面滑動(dòng)巖質(zhì)邊

18、坡在滑動(dòng)面AB，BM上布置接觸單元，坡體到達(dá)極限狀態(tài)后的破壞滑動(dòng)如圖12所示，并把有限元計(jì)算結(jié)果，與傳統(tǒng)極限平衡方法Spencer法進(jìn)行比照，接觸單元的相關(guān)力學(xué)參數(shù)以及兩種計(jì)算結(jié)果比照方表5所列。圖 12 坡體到達(dá)極限狀態(tài)后的破壞滑動(dòng)表 5例五計(jì)算結(jié)果另外，在單元?jiǎng)澐值倪^程中，在兩個(gè)滑動(dòng)面的交匯處形成了尖角，在尖角處形成較大的應(yīng)力集中，求解時(shí)會(huì)產(chǎn)生病態(tài)方程。為了防止這些建模問題，需要在實(shí)體模型上，使用線的倒角來使尖角光滑化，或者在曲率突然變化的區(qū)域使用更細(xì)的網(wǎng)格。例六，泥巖層上粉質(zhì)粘土邊坡計(jì)算分析坡體材料力學(xué)參數(shù)為彈性模量40MPa、泊松比0.28、重度19kN/m3、黏聚力20kPa、內(nèi)摩擦

19、角25°，地基材料力學(xué)參數(shù)為彈性模量400MPa、泊松比0.23、重度24kN/m3、黏聚力400kPa、內(nèi)摩擦角32°。模型幾何尺寸及邊界條件如圖13所示：圖 13 模型幾何尺寸及邊界條件通過ABAQUS有限元軟件計(jì)算，當(dāng)邊坡的平安系數(shù)為1.3時(shí)，計(jì)算不收斂，通過Slide軟件利用瑞典條分法計(jì)算得到的邊坡平安穩(wěn)定系數(shù)為1.295，瑞典條分法計(jì)算的滑面和有限元計(jì)算的塑性區(qū)如圖14所示，可見兩種方法計(jì)算出的平安系數(shù)和滑面吻合性較好。5 參考文獻(xiàn)1 Smith, I. M. & Hobbs, R. (1974). Finite element analysis of c

20、entrifuged and built-up slopes. Ge Âotechnique 24, No. 4, 531-559:2Taylor, D. W. (1937). Stability of earth slopes. J. Boston Soc. Civ. Eng. 24, 197-246:3Zienkiewicz, O. C., Humpheson, C. & Lewis, R. W.(1975). Associated and non-associated viscoplasticityand plasticity in soil mechanics. Ge

21、otechnique 25,671-689:4Griffiths, D. V. (1980). Finite element analyses of walls, footings and slopes. PhD thesis, University of Manchester.5Bishop, A. W. & Morgenstern, N. R. (1960). Stability coefficients for earth slopes. Geotechnique 10, 129-150:6 Griffiths, D. V. (1989). Computation of collapse loads in geomechanics by finite elements. Ing Arch 59, 237-244:7 Potts, D. M., Dounias, G. T. & Vaughan, P. R. (1990). Finite eleme