有限元強(qiáng)折減法

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

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

4、準(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)系，屈服面是一個不等角六棱錐，在平面上是一個等邊不等角六邊形。(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)力空間中其屈服面為一個圓錐，在平面上是一個圓形。3)D-P準(zhǔn)則轉(zhuǎn)換為Mohr-Coulomb準(zhǔn)則首先引入?yún)?shù)b，

5、如式(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)對比可知兩個準(zhǔn)則之間的轉(zhuǎn)換關(guān)系如式(6)所示：sinb=1.52-b+32cbcosb=kf2-b+3因此，當(dāng)b=0時，即外角點(diǎn)外接DP圓的兩個試驗(yàn)常數(shù)分別如式(7)所示，當(dāng)b=1時，即內(nèi)角點(diǎn)外接DP圓的兩個試驗(yàn)常數(shù)分別如式(8)所示。=2sin3(3-sin)，kf=6ccos3(3-sin) (7)=2sin3(3+sin)，kf=6

6、ccos3(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ù)Ft，得到一組新的c和值，作為一個新的強(qiáng)度參數(shù)輸入進(jìn)行試算，當(dāng)計(jì)算不收斂時，對應(yīng)的Ft即為所求的安全系數(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)，-FtI1+J2=kfFt (9)(3

7、)不同屈服條件下安全系數(shù)轉(zhuǎn)換13(10)首先引入Mohr-Coulomb等面積圓屈服準(zhǔn)則，在平面上，其屈服面是一個圓，并且面積與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)則為DP2準(zhǔn)則，其試驗(yàn)常數(shù)分別為2，kf2。把DP1準(zhǔn)則表示

8、為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)取不同值時可以得到不同的值如表1所列：表 1 不同內(nèi)摩擦角時的值4.失穩(wěn)判據(jù)目前兩個比較主流的失穩(wěn)判據(jù)分別是有限元計(jì)算中力不平衡和位移的不收斂以及廣義塑性應(yīng)變或者等效塑性應(yīng)變從坡腳到坡頂貫通。Griffiths9和鄭穎人11,12,13,14都使用計(jì)算不收斂作為失穩(wěn)判據(jù)。Griffiths9提出，當(dāng)在用戶定義的最大迭代數(shù)目下計(jì)算仍不收斂時，則沒有任何一種應(yīng)力分布方式可以同時滿足

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

10、效準(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ì)邊坡每一個安全系數(shù)對應(yīng)的迭代次數(shù)如表2所列，當(dāng)真正的安全系數(shù)接近時需要更多的迭代次數(shù)。表 2 例一計(jì)算結(jié)果當(dāng)安全系數(shù)為1.4時，無量綱位移Emax/H2突變，并且此時計(jì)算無法收斂，在此情況下有限元計(jì)算結(jié)果與Bishop & Morgenstern5給出的結(jié)果

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

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

13、軟弱層強(qiáng)度只有其周圍土體60%( Cu2/Cu1=0.6)時的潛在滑面，此時圓弧滑面和沿軟弱層的三段線式滑面都有可能發(fā)展，至少存在兩種明顯的滑動機(jī)制。圖 4有軟弱層的不排水黏性邊坡圖 5 不同軟弱層強(qiáng)度時的安全系數(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°,求坡角分別為30°,35°,40°,45

14、76;,50°時邊坡的安全系數(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所示，巖體中存在兩組方向不同的軟弱結(jié)構(gòu)面，貫通率100%，第一組軟弱結(jié)構(gòu)面傾角為30°，

15、平均間距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)度折減，求得坡體破壞時的運(yùn)動矢量如圖8所示，滑動面如圖9(a)所示，它是最先貫通的塑性區(qū)，塑性區(qū)貫通并不等于破壞，當(dāng)塑性區(qū)貫通后繼續(xù)發(fā)展到一定程度，巖體發(fā)生整體破壞，同時出現(xiàn)第二條貫通的塑性面，如圖9(b)所示

16、。求得的穩(wěn)定安全系數(shù)如表4所列，其中，極限平衡方法計(jì)算結(jié)果是根據(jù)最先貫通的那一條滑動面求得的。圖 7 巖質(zhì)邊坡節(jié)理 圖 8 極限狀態(tài)時的塑性區(qū)圖 9坡體破壞時的運(yùn)動矢量圖表 4 例四計(jì)算結(jié)果例五，存在接觸問題的邊坡穩(wěn)定性分析12當(dāng)邊坡中存在如圖10所示的硬性結(jié)構(gòu)面時，不能按照例四中軟弱結(jié)構(gòu)面的方法進(jìn)行處理，可以采用接觸單元來模擬硬性接觸面的不連續(xù)性。按照Mohr-Coulomb定律來定義接觸面上的摩擦行為，如式11所示，則其接觸面上的安全系數(shù)定義如式11所示。=c+tan，0 Ft=c/c=tan/tan (11)圖 10 無充填的硬性結(jié)構(gòu)面圖11所示為兩個直線滑面組成的折線型滑體ABMCD。

17、巖體重度=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ì)邊坡在滑動面AB，BM上布置接觸單元，坡體達(dá)到極限狀態(tài)后的破壞滑動如圖12所示，并把有限元計(jì)算結(jié)果，與傳統(tǒng)極限平衡方法Spencer法進(jìn)行對比，接觸單元的相關(guān)力學(xué)參數(shù)以及兩種計(jì)算結(jié)果對比如表5所列。圖 12 坡體達(dá)到極限狀態(tài)后的破壞滑動表

18、 5例五計(jì)算結(jié)果另外，在單元劃分的過程中，在兩個滑動面的交匯處形成了尖角，在尖角處形成較大的應(yīng)力集中，求解時會產(chǎn)生病態(tài)方程。為了避免這些建模問題，需要在實(shí)體模型上，使用線的倒角來使尖角光滑化，或者在曲率突然變化的區(qū)域使用更細(xì)的網(wǎng)格。例六，泥巖層上粉質(zhì)粘土邊坡計(jì)算分析坡體材料力學(xué)參數(shù)為彈性模量40MPa、泊松比0.28、重度19kN/m3、黏聚力20kPa、內(nèi)摩擦角25°，地基材料力學(xué)參數(shù)為彈性模量400MPa、泊松比0.23、重度24kN/m3、黏聚力400kPa、內(nèi)摩擦角32°。模型幾何尺寸及邊界條件如圖13所示：圖 13 模型幾何尺寸及邊界條件通過ABAQUS有限元軟件

19、計(jì)算，當(dāng)邊坡的安全系數(shù)為1.3時，計(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 centrifuged and built-up slopes. Ge Âotechnique 24, No. 4, 531-559:2Taylor, D. W. (1937). Stability of earth slope

20、s. 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. Geotechnique 25,671-689:4Griffiths, D. V. (1980). Finite element analyses of walls, footings and slopes. PhD thesis, Uni

21、versity 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 element analysis