有限元強度折減法

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

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

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

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

5、omb準那么首先引入?yún)?shù)b，如式(3)所示：b=2-31-3 (3)(4)那么，I1和J2分別可轉化為式(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)比照可知兩個準那么之間的轉換關系如式(6)所示：sinb=1.52-b+32cbcosb=kf2-b+3因此，當b=0時，即外角點外接DP圓的兩個試驗常數(shù)分別如式(7)所示，當b=1時，即內(nèi)角點外接DP圓的兩個試驗常數(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準那么中的平安系數(shù)1955年，Bishop15首先在邊坡穩(wěn)定性分析中提出了抗剪強度折減的概念，在有限元強度折減法中通過將坡體的強度參數(shù)：黏聚力c和內(nèi)摩擦角同時除一個折減系數(shù)Ft，得到一組新的c和值，作為一個新的強度參數(shù)輸入進行試算，當計算不收斂時，對應的Ft即為所求的平安系數(shù)，此時坡體到達極限狀態(tài)，發(fā)生剪切破壞。c=c/Ft=arctan(tan/Ft)(2)D-P(Drucker-Prager)準那么中的平安系數(shù)取Ft為D-P準那么中的強度折減系數(shù)，那么D-P準那么可以表示為式(9

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

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

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

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

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

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

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

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

15、在兩組方向不同的軟弱結構面，貫穿率100%，第一組軟弱結構面傾角為30°，平均間距10m；第二組軟弱結構面傾角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°。按照二維平面應變問題建立有限元模型，按照連續(xù)介質(zhì)處理。通過有限元強度折減，求得坡體破壞時的運動矢量如圖8所示，滑動面如圖9(a)所示，它是最先貫穿的塑性區(qū)，塑性區(qū)貫穿并不等于破壞，當塑性區(qū)貫穿

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

17、Coulomb定律來定義接觸面上的摩擦行為，如式11所示，那么其接觸面上的平安系數(shù)定義如式11所示。=c+tan，0 Ft=c/c=tan/tan (11)圖 10 無充填的硬性結構面圖11所示為兩個直線滑面組成的折線型滑體ABMCD。巖體重度=20kN/m3，彈性模量E=109Pa?；瑝KABCD面積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 折線型平面滑動巖質(zhì)邊

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

19、角25°，地基材料力學參數(shù)為彈性模量400MPa、泊松比0.23、重度24kN/m3、黏聚力400kPa、內(nèi)摩擦角32°。模型幾何尺寸及邊界條件如圖13所示：圖 13 模型幾何尺寸及邊界條件通過ABAQUS有限元軟件計算，當邊坡的平安系數(shù)為1.3時，計算不收斂，通過Slide軟件利用瑞典條分法計算得到的邊坡平安穩(wěn)定系數(shù)為1.295，瑞典條分法計算的滑面和有限元計算的塑性區(qū)如圖14所示，可見兩種方法計算出的平安系數(shù)和滑面吻合性較好。5 參考文獻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