有限元極限平衡法在黃土邊坡穩(wěn)定性分析中的應用探討

王　月，葉朝良

(石家莊鐵道大學土木工程學院，石家莊　050043)

?

有限元極限平衡法在黃土邊坡穩(wěn)定性分析中的應用探討

王月，葉朝良

(石家莊鐵道大學土木工程學院，石家莊050043)

摘要：對目前邊坡穩(wěn)定分析中的極限平衡和有限元兩類方法進行分析闡述，在此基礎(chǔ)上，結(jié)合黃土邊坡工程算例，利用geostudio軟件將有限元與極限平衡法相結(jié)合進行黃土邊坡穩(wěn)定性分析，并對其計算結(jié)果進行比較。算例分析結(jié)果表明：極限平衡法在求解過程中作出各種假定來解決超靜定問題；有限元法克服了極限平衡法的不足，考慮土體的本構(gòu)關(guān)系，但不能確定可能破裂面，對邊坡失穩(wěn)破壞的判斷沒有一個定量的標準，將有限元法和極限平衡法相結(jié)合，取長補短，既能反映邊坡的穩(wěn)定和變形之間的關(guān)系，準確描繪出破壞面的特征，又能得到更精準可靠的安全系數(shù)；有限元極限平衡法可自動搜索出坡肩帶有垂直裂縫的組合滑動面，比極限平衡法搜索出的單一圓弧滑動面更貼近實際情況，提高了邊坡穩(wěn)定性分析的合理性，為類似邊坡穩(wěn)定性分析方法提供參考。

關(guān)鍵詞：黃土邊坡；極限平衡法；有限元；穩(wěn)定性分析；裂縫

黃土是地球上分布廣泛且性質(zhì)特殊的一種沉積物，其結(jié)構(gòu)疏松，具有垂直節(jié)理。近年來，隨著黃土地區(qū)的工程建設(shè)規(guī)模不斷擴大，經(jīng)常會遇到黃土邊坡的防護和治理等問題。黃土邊坡在失穩(wěn)之前，常常在坡頂出現(xiàn)張拉裂縫，其失穩(wěn)破壞情況比一般邊坡更加復雜，因此，找出一種適用于黃土邊坡的邊坡穩(wěn)定性分析方法是十分必要的。

目前，邊坡穩(wěn)定性計算方法主要有極限平衡法和有限元法， 但各自在應用中存在一些問題。極限平衡法可以滿足力和力矩的平衡，摩爾-庫侖破壞準則和應力邊界條件，其方法簡單，概念清楚，因此在工程中被廣泛應用。但其未考慮土體本身的應力-應變關(guān)系和實際的工作狀態(tài)，所得的安全系數(shù)只是假定滑裂面上的平均安全度，求出的條間力和滑條底部反力并不是真實數(shù)值，沒有考慮邊坡失穩(wěn)的破壞機理。有限單元法是目前在邊坡穩(wěn)定評價中應用最廣泛的數(shù)值方法。但是，這類數(shù)值方法沒有給出一個明確的穩(wěn)定性安全系數(shù)以及可能的破壞面，限制了其在工程中的應用。

近年來許多學者嘗試將兩種方法相結(jié)合并取得了一定的成果，曾亞武，田偉明[1]將有限元分析計算的應力結(jié)果通過應力張量變換的方法，求得滑動面上的應力分布，然后通過積分或求和的方法來求解對應滑面穩(wěn)定性安全系數(shù)。楊輝[2]在有限元法分析應力的基礎(chǔ)上將穩(wěn)定性分析問題轉(zhuǎn)化為帶有約束條件的數(shù)學規(guī)劃問題，可以描述為：設(shè)定安全系數(shù)Fs為目標函數(shù)，利用有限數(shù)目的坐標節(jié)點(xi，yi)將曲線v離散，在已知應力場內(nèi)，給定一組坐標點的x坐標，變化尋找相應的y坐標，使其所確定的曲線v對應的安全系數(shù)為最小。鞏留杰[3]根據(jù)模型同一高度塑性應變特點尋找滑裂面，輸出滑裂面上節(jié)點主應力σx、σy、τxy，然后計算滑裂面上總的正應力、切應力，進而計算安全系數(shù)，本文所采用的有限元極限平衡法通過限元法分析計算巖土體內(nèi)較為真實的應力分布，采用插值方法得到已給定滑動面上的應力值，然后結(jié)合極限平衡法的概念得出邊坡的穩(wěn)定性安全系數(shù)，結(jié)合黃土邊坡的具體實例進行分析，對不同方法的優(yōu)缺點進行了分析研究。

1常用邊坡穩(wěn)定性分析方法

1.1極限平衡法

極限平衡法假設(shè)土條為剛體，依據(jù)純粹的靜力學原理，即力的平衡方程和力矩的平衡以及摩爾-庫侖強度準則來判斷邊坡的穩(wěn)定性。極限平衡法不涉及應變和位移，沒有考慮位移的兼容性，因而不可能獲得真實的應力分布，這一缺點造成的該方法的局限性。

極限平衡條分法采用了簡化假定的方法，減少了未知量的數(shù)目，表1總結(jié)了常用的幾種極限平衡條分法的條間力假定及主要特征。極限平衡法的詳細內(nèi)容見文獻[4]。

表1　典型極限平衡法的條間力假設(shè)及主要特征

1.2有限元法

有限元法將土體離散化成有限個通過節(jié)點相互聯(lián)系的單元，按照彈塑性理論，得出土體的應力應變關(guān)系，可以解決各種加荷歷史和地質(zhì)條件復雜的邊坡問題。近年來國內(nèi)外很多學者運用有限元法對邊坡進行穩(wěn)定性分析[5-12]，發(fā)現(xiàn)盡管該方法有其自身的優(yōu)越性但也存在一定的弊端。該方法常常把塑性區(qū)的貫通作為判別邊坡失穩(wěn)的判據(jù)，但很多工程實例證明塑性區(qū)出現(xiàn)貫通邊坡并不一定失穩(wěn)。而有限元強度折減法通常把計算不收斂作為邊坡失穩(wěn)的判據(jù)，這種情況下迭代方法和收斂容差的選取對計算結(jié)果產(chǎn)生了很大的影響。有限元法對邊坡失穩(wěn)并沒有一個明確的判斷標準，是塑性區(qū)的貫通還是計算不收斂，將有限元法和極限平衡法相結(jié)合，以安全系數(shù)作為邊坡失穩(wěn)的判據(jù)，很好地彌補了有限元法的缺點。

1.3有限元極限平衡法

邊坡失穩(wěn)從根本上說是由于各種因素引起了坡體內(nèi)部應力和強度的變化，所以邊坡的穩(wěn)定性與其應力狀態(tài)密切相關(guān)。本文采用的有限元極限平衡法通過有限元計算輸出模型區(qū)域內(nèi)真實應力場分布，采用插值方法得到假定滑動面上的應力值，依據(jù)下述公式計算沿滑動面的安全系數(shù)，進而在真實的應力場中得到確定的最危險滑動面，不需要對條間作用力做任何假定，通過有限元法計算出應力場，不存在迭代和收斂問題，易于工程中的推廣和應用。以傳統(tǒng)的極限平衡法為基礎(chǔ)使該方法更具有說服力，極限平衡法雖然并不精確，但其具有悠久的歷史，經(jīng)過了多年的使用和觀察校準，這種可信度是一種完全新的方法所不能提供的。

有限元極限平衡法依據(jù) Mohr-Coulomb 強度準則確定土體中任意一點的土體抗剪強度，滑動面安全系數(shù)Fs定義為沿滑動面土體抗剪強度與實際剪應力的比值

通過彈塑性有限元分析，計算出對于每個單元高斯積分點的σx、σy、τxy，進而利用差值函數(shù)求出各土條中點的σx、σy、τxy，用摩爾圓來計算土條底部的正應力和剪切力，按照上述公式計算Fs值，重復上述處理直到第n個土條。

2工程應用

2.1工程概況

選取位于隴西地區(qū)的天然黃土邊坡為例，考慮不同分析方法，對其進行邊坡穩(wěn)定性評價。工程所在地區(qū)屬于中溫帶半干旱大陸性氣候，干旱少雨，黃土邊坡坡高20 m，坡面傾角56.4°，坡面平整，土質(zhì)均勻，為砂質(zhì)黃土，可將該斜坡簡化為如下計算模型：土體采用莫爾-庫侖理想彈塑性模型，模型底部采用剛性邊界，材料參數(shù)：容重γ=18.5 kN/m3，彈性模量E=45 MPa，v=0. 3，黏聚力c=16.0 kPa，內(nèi)摩擦角27°，該邊坡為老黃土單一結(jié)構(gòu)邊坡，分析過程不考慮地下水的影響。

2.2計算結(jié)果對比分析

采用geostudio程序，通過有限元SIGMA/W模塊，計算出應力狀態(tài)，將結(jié)果導入邊坡SLOPE/W模塊，計算得出滑坡的潛在滑動面及安全系數(shù)，實現(xiàn)有限元法和極限平衡法的結(jié)合。

2.2.1最危險滑裂面形狀

有限元的計算模型如圖1所示，限定模型兩側(cè)的水平位移和模型底部兩個方向的位移，采用四邊形和三角形單元劃分網(wǎng)格，該有限元模型共劃分301個單元，包含337個節(jié)點。

為了便于分析邊坡坡頂處的變形情況，分別提取坡頂處(截面aa)、19 m處(截面bb)、18 m處(截面cc)、17 m處(截面dd)、16 m處(截面ee)截面上各點的水平方向應力值、應變值和各截面處的剪應力值，見圖2～圖4。由圖2可以看出，坡頂處和19 m處截面水平向應力的變化趨勢明顯區(qū)別于其他三截面處，應力值先逐漸增大，在水平位置為18 m附近，水平應力值達到峰值，而后又迅速減小。由此可推斷，在水平位置為18 m左右的坡肩附近有可能產(chǎn)生裂縫，而18 m處截面的水平應力雖然也在水平位置18 m附近達到峰值，但其曲線較為平緩，應力值沒有發(fā)生較大的變化。17、16 m截面處水平應力曲線并沒有較為明顯的峰值出現(xiàn)，由此可見，裂縫的開展在邊坡坡頂截面(aa)和18 m處截面(cc)之間，即坡肩裂縫的深度在2 m之內(nèi)。

圖2　各截面水平應力分布情況

圖3　各截面水平應變分布情況

圖4　各截面的剪應力分布情況

圖3中坡頂截面和19 m截面處的水平向應變在水平位置18 m附近出現(xiàn)峰值，說明該位置處的土體發(fā)生了較大的拉伸變形，進而可推斷水平位置18 m附近是坡頂裂縫最可能出現(xiàn)的地帶。同時從圖中可以發(fā)現(xiàn)18、17、16 m截面處的應變曲線較為平緩，并沒有明顯的峰值出現(xiàn)，且全部為負值，說明坡頂豎向裂縫的開展不會到達18 m截面處，與前面的論證相符。

各截面的剪應力分布如圖4所示，從該圖中可以得到如下兩個推斷：(1)所取5個截面的剪應力均在18 m附近發(fā)生正負的改變，剪應力的方向發(fā)生變化，其位置與邊坡上部的最大拉應力區(qū)一致。說明邊坡在該區(qū)域內(nèi)易發(fā)生剪切破壞；(2)坡頂截面處和19 m截面處的剪應力曲線與x軸的交點基本上是重合的，即坡頂截面和19 m截面處的剪應力是在同一水平位置達到零值，發(fā)生方向的改變，據(jù)此可以推斷，至少在距離坡頂1 m內(nèi)，斷裂是垂直開展的。這與前面的分析一致，說明了該邊坡的破壞會在坡肩處形成垂直裂縫，且裂縫的水平位置在18 m附近，而裂縫的豎向開展范圍在坡頂截面(aa)與18 m處截面(cc)之間。

利用geostudio軟件分別用有限元極限平衡法和傳統(tǒng)的極限平衡法(Morgenstern-Price法、Bishop法、Ordinary條分法和Bishop)對該邊坡進行穩(wěn)定性分析，發(fā)現(xiàn)用有限元極限平衡法得到的滑裂面在坡肩的地方呈垂直劃裂狀，其滑面為垂直裂隙面與圓弧面組合形成的滑裂面，如圖5所示，裂縫的水平坐標為18.1，且裂縫的開展范圍在坡頂截面(aa)與18 m處截面(cc)之間，與之前的推斷相符。而4種極限平衡法所找到的最危險滑動面的形狀和位置基本一致，呈簡單圓弧狀。

圖5　滑動面形狀對比

極限平衡法假設(shè)土條為剛體，用分析剛體的辦法來分析土體，且對條間力的作用做了各種假設(shè)，其計算出的滑動面上的應力狀態(tài)顯然是不真實的。因此在某些情況下極限平衡法得出的最危險滑動面的位置有所偏差，而有限元法極限平衡法將土坡當成變形體，搜索出的滑動面形狀更接近真實狀態(tài)。基于有限元的邊坡穩(wěn)定性分析方法考慮了土體的本構(gòu)關(guān)系以及變形對應力的影響，基于土體內(nèi)真實的應力狀態(tài)搜索最危險滑裂面，得到的滑裂面在坡肩呈垂直滑裂狀態(tài)，坡肩裂縫的形成減少了劃弧的長度，裂縫不斷加深與圓弧滑動面相交，形成由圓弧面和垂直裂隙面共同組成的滑動面。裂縫的形成是一個需要從地質(zhì)和力學的角度綜合考慮的復雜問題，極限平衡法只是簡單地將土條視為剛體，通過靜力學平衡方程求解，只得到簡單的圓弧狀的滑裂面，而無法得到坡肩部分由張拉破壞產(chǎn)生的直立的裂縫，顯然與實際情況存在出入。

常規(guī)的極限平衡法不能得到真實的滑裂面形狀，針對這種情況，部分學者對常規(guī)的極限平衡法進行修正，得到了考慮裂縫的極限平衡法，該方法認為，若在受力分析中出現(xiàn)了拉應力則認為是裂縫，進而搜索最危險的裂縫，在最危險裂縫處與圓弧滑裂面相交。該方法存在如下兩個問題：(1)在搜索滑裂面之前，是否考慮裂縫需要人為的分析并進行選擇，而在有限元極限平衡法中，不需要分析是否考慮裂縫，基于土體內(nèi)真實應力狀態(tài)的分析搜索最危險滑裂面，這是一個自動化的過程，所判定的最危險滑裂面的形狀位置真實可靠；(2)極限平衡法本身對條間力引入了各種假定，求出的條間力和滑條底部反力并不是真實數(shù)值，由土條的受力出現(xiàn)拉應力來判斷裂縫位置，是不準確的。由常規(guī)極限平衡法、考慮裂縫的極限平衡法和有限元極限平衡法這3種方法得到的滑裂面的形狀位置如圖6所示。從該圖可以看出，考慮裂縫的極限平衡法搜索出的滑動面為垂直裂縫加圓弧形，其中圓弧段與未考慮裂縫的極限平衡法所搜索出的圓弧滑動面基本上是重合的。而采用有限元極限平衡法所得到的裂縫位置與采用考慮裂縫的極限平衡法時十分接近，說明采用有限元極限平衡法得到的裂縫位置是可靠的。

圖6　不同方法的滑裂面位置對比

2.2.2穩(wěn)定性安全系數(shù)

極限平衡方法之間的區(qū)別在于條塊間相互作用力的大小和方向的假定不同，其中普通條分法假定條塊間的作用力為零，因此計算誤差較大。Morgenstern-Price法和Bishop法所得的結(jié)果比較相近，且相比于其他條分法結(jié)果更為準確。

Janbu法硬性規(guī)定了土條側(cè)向力作用點的絕對位置，迭代過程不能靈活調(diào)整，當條塊數(shù)量太多、條塊寬度太小時，計算結(jié)果往往不收斂。普通條分法不考慮條塊間作用力，引起了較大的誤差 ， 使結(jié)果偏小。在實際工程中不建議采取這兩種方法。

Morgenstern-Price法得到的邊坡安全系數(shù)較大，收斂較好。其原因是它在靜力平衡要求、滑裂面的形狀以及多余未知數(shù)的迭代各方面均不做任何假定，能較好反映出最危險滑面各土條間相互作用力的情況，Morgenstern-Price法是計算任意形狀滑動面安全系數(shù)最為嚴謹?shù)囊环N條分法。

利用有限元應力結(jié)果計算得的安全系數(shù)為1.149，比各種條分法算得安全系數(shù)都有所提高，最大相差14%(普通條分法)，最小相差10%(Bishop法)，比考慮裂縫的極限平衡法最大相差17%，最小相差10%。在實際工程中，利用極限平衡法得到安全系數(shù)低于規(guī)范時，邊坡往往還處于穩(wěn)定狀態(tài)，計算結(jié)果和實際情況存在偏差。安全系數(shù)產(chǎn)生差異的原因是沿滑裂面上應力的分布不同。圖7為條間力為一恒定函數(shù)的Morgenstern-Price法、有限元法、考慮裂隙的極限平衡法得到的正應力分布對比圖，由該圖可以得出如下結(jié)論。(1)水平位置17～18 m處采用極限平衡法時滑動面上存在拉應力，而考慮裂隙的極限平衡法則認為在此段產(chǎn)生裂縫，進而搜索出最危險裂縫在水平位置18 m處，邊坡由此開裂，此處正應力為零。而在18～30 m處，采用考慮裂隙的極限平衡法時滑動面上的正應力與采用極限平衡法是基本相同，故采用這兩種方法搜索到的圓弧滑裂段的形狀位置是相同的。(2)坡趾區(qū)有限元極限平衡法得到的正應力大于極限平衡法，在以有限元為基礎(chǔ)的分析方法中很好地體現(xiàn)出了邊坡坡趾的剪應力集中現(xiàn)象，而極限 平衡法中正應力源自于土條的重力，無法反映出應力集中現(xiàn)象。有限元極限平衡法考慮了土體的本構(gòu)關(guān)系，以及變形對應力的影響，得到的滑動面形狀和滑弧位置更貼和實際情況。

圖7　滑裂面上的正應力分布

另外，從表2中可看出，采用Bishop法進行計算時，裂縫的存在對安全系數(shù)并無影響，而采用Janbu法時，裂縫的存在使安全系數(shù)降低了4.07%。進一步計算分析表明，采用Bishop法進行邊坡穩(wěn)定性計算時，裂縫對安全系數(shù)的影響是及其微小的，這顯然與實際情況不相符，而采用Janbu法考慮裂縫時，安全系數(shù)降低的最為顯著，故Bishop法不適用于黃土邊坡的穩(wěn)定性分析。

由圖8可看出，邊坡坡趾存在應力集中現(xiàn)象，剪應力沿潛在滑動面方向逐漸增大，尤其是坡腳部位剪應力較大，當其剪切應力超過了該點的抗剪強度，則該點產(chǎn)生剪切破壞，隨著滑坡體滑動的逐步發(fā)展，潛在滑動面上發(fā)生應力重分配，破壞點逐步增多進而形成連續(xù)破壞面。垂直方向上的應力分布中值得注意的一點為，應力等值線的輪廓并不是一個恒定的厚度，如圖9所示，應力為80 kPa的等值線在坡趾處發(fā)生突變，一部分輪廓厚度明顯減小，這意味著垂直方向上的應力不僅受重力的影響，也受到剪切應力的影響。

表2　安全系數(shù)計算結(jié)果

圖8　剪應力分布

圖9　垂直向應力分布

3結(jié)論與建議

通過對邊坡穩(wěn)定性分析方法的分析總結(jié)，并結(jié)合工程算例的計算分析，得到如下結(jié)論和建議。

(1)針對具有特殊破壞面形狀的黃土邊坡，建議有限元極限平衡法將更具可靠性。對于節(jié)理發(fā)育的黃土邊坡，傳統(tǒng)的極限平衡法無法真實地反映坡肩裂縫的形成，而有限元極限平衡法自動搜索出了由垂直裂隙面和圓弧面共同組成的最危險滑裂面，更貼近現(xiàn)實中黃土邊坡滑動面的形狀。

(2)有限元極限平衡法與考慮裂隙的極限平衡法所得到的裂縫位置十分接近，且其計算的安全系數(shù)提高了大約10%，說明采用有限元極限平衡法進行黃土邊坡穩(wěn)定性分析的合理性和精確性。

參考文獻：

[1]曾亞武，田偉明.邊坡穩(wěn)定性分析的有限元法與極限平衡法的結(jié)合[J].巖石力學與工程學報，2005(S2):5355-5359.

[2]楊輝.邊坡穩(wěn)定分析的有限元極限平衡法原理及程序?qū)崿F(xiàn)[D].北京：北京交通大學，2014.

[3]鞏留杰.基于有限元計算的邊坡穩(wěn)定極限平衡法研究[D].長沙：湖南大學，2012.

[4]陳祖煜.土質(zhì)邊坡穩(wěn)定分析：原理·方法·程序[M].北京：中國水利水電出版社，2003.

[5]孫聰，李春光，鄭宏，等.同時考慮張拉及剪切破壞的邊坡上限原理有限元法[J].巖石力學與工程學報，2015(S1):2783-2791.

[6]唐亞明，馮衛(wèi)，李政國.黃土滑塌研究進展[J].地球科學進展，2015(1):26-36.

[7]湯連生，桑海濤，羅珍貴，等.土體抗拉張力學特性研究進展[J].地球科學進展，2015(3):297-309.

[8]遲世春，關(guān)立軍.基于強度折減的拉格朗日差分方法分析土坡穩(wěn)定性[J].巖土工程學報，2004(1):42-46.

[9]唐亞明.陜北黃土滑坡風險評價及監(jiān)測預警技術(shù)方法研究[D].北京：中國地質(zhì)大學，2012.

[10]段釗.黃土滑坡觸發(fā)機理研究[D].西安：長安大學，2013.

[11]王振文.黃土深路塹邊坡穩(wěn)定性分析及邊坡坡型的合理設(shè)計[J].鐵道標準設(shè)計，2013(4):17-20.

[12]李守升.某黃土邊坡滑塌機理分析及治理方案研究[J].鐵道標準設(shè)計，2014(S1):70-72.

收稿日期：2015-10-13； 修回日期：2015-12-21

基金項目：國家自然基金項目(50978172)；鐵道部科技研究開發(fā)計劃項目(2010G018-B-4)

作者簡介：王月(1991—)，女,碩士研究生，E-mail：1296739730@qq.com。

文章編號：1004-2954(2016)07-0038-05

中圖分類號：U213.1+3

文獻標識碼：A

DOI:10.13238/j.issn.1004-2954.2016.07.009

Application of Finite Element Limit Equilibrium Method in Stability Analysis of Loess Slope

WANG Yue， YE Chao-liang

(College of Civil Engineering， Shijiazhuang Tiedao University， Shijiazhuang 050043， China)

Abstract：The limit equilibrium and finite element methods for slope stability analysis are analyzed. On the basic of it， the limit equilibrium method based on the finite element stress field is analyzed by using GeoStudio software with reference to the engineering calculation of loess slope. The example analysis results show that the limit equilibrium method in solving process has made various assumptions to solve the statically indeterminate problem; the finite element method can overcome the deficiency of the limit equilibrium method in consideration of the constitutive relationship， but can not determine potential rupture surface. There is no quantitative criterion for judging the failure of slope. The combination of the finite element method with the limit equilibrium method can not only reflect the relationship between slope stability and deformation but also accurately depict the failure surface features with more accurate and reliable safety factor. The limit equilibrium method based on the finite element stress field can automatically search for the combined sliding surface with vertical cracks on the slope shoulder， which is much closer to the actual situation than the single arc sliding surface obtained with the limit equilibrium method. Thus the rationality of the slope stability analysis is improved to provide references for similar slope stability analysis．

Key words：Loess slope; Limit equilibrium method; Finite element method; Stability analysis; Cracks