基于有限元極限平衡法的巖土結構穩(wěn)定分析

邵龍?zhí)叮鹾迫唬獫珊?/p>

(1.大連理工大學工業(yè)裝備結構分析國家重點實驗室，遼寧 大連 116024；2.中山大學土木工程學院，廣東 珠海 519082；3.大連理工大學白俄羅斯國立大學聯(lián)合學院，遼寧 大連 116024)

在土力學中，有三大經(jīng)典問題，即邊坡穩(wěn)定、土壓力和地基承載力。有學者指出，這三個問題都可表征為邊坡的穩(wěn)定分析問題。土壓力問題是一個具有垂直表面的邊坡穩(wěn)定問題，地基承載力問題是一個具有水平表面的邊坡穩(wěn)定問題[1]，而邊坡穩(wěn)定最常見的地質災害是滑坡，滑坡的發(fā)生是從一點的破壞開始的[2]。我國是一個地質災害十分嚴重的國家，據(jù)統(tǒng)計，我國每年由地質災害所造成的經(jīng)濟損失達200億～500億人民幣[3]。

邊坡穩(wěn)定性分析是坡體結構設計、失穩(wěn)防護的基礎。目前，邊坡穩(wěn)定性分析方法主要有極限平衡法、極限分析法和以有限元為主的數(shù)值分析方法等。在分析土壓力問題時，鄭穎人等[4]采用極限平衡法的思想假定擋土墻后土體處于極限平衡狀態(tài)計算土體主動與被動土壓力。也有學者將邊坡劃分為若干豎條，通過分析作用在土條上力和力矩的平衡關系對邊坡的穩(wěn)定性能進行評價。極限分析方法是應用剛塑性材料處于極限狀態(tài)的上下限定理求解極限載荷的一種解析方法。

采用有限元法分析邊坡穩(wěn)定性又可分為有限元強度折減法與有限元極限平衡法。有限元強度折減法是將土體強度參數(shù)進行折減，將土體破壞時的折減系數(shù)視為安全系數(shù)，該方法存在的問題是破壞判別標準不統(tǒng)一、需要反復迭代計算、不同的材料取相同的折減系數(shù)不合理以及強度折減可能會影響材料的應力應變關系參數(shù)從而導致有可能違背材料自身性質等[5-7]。有限元極限平衡法依據(jù)極限平衡理論定義安全系數(shù)，根據(jù)有限元計算得到的應力分布確定臨界滑動面及其相應的安全系數(shù)，進而進行巖土結構穩(wěn)定性分析。大量實踐結果表明，有限元極限平衡法可用于各種巖土結構的穩(wěn)定分析[8-11]。

1 滑動穩(wěn)定安全系數(shù)定義

應用莫爾-庫倫強度準則，邵龍?zhí)兜萚11]將土體內一點的極限平衡條件推廣到曲面，得到任意曲面上任意形狀土體的極限平衡條件。對于并未處于極限平衡狀態(tài)的實際巖土結構，采用折減抗剪強度或者增加荷載的方法令其達到沿曲面的極限平衡狀態(tài)，由此可以得到定義安全系數(shù)，并且發(fā)現(xiàn)滑動面上抗剪力與剪切力代數(shù)和之比的安全系數(shù)定義的物理意義，既是曲面上任意形狀的土體達到極限平衡狀態(tài)時抗剪強度強度折減系數(shù)的中值，也是使曲面上土體沿著曲面達到極限平衡的超載系數(shù)的中值[12]。

1.1 土體內一點的極限平衡條件

土體抵抗剪切破壞的能力稱為土的抗剪強度。在常見壓力范圍內，土的抗剪強度主要由顆粒之間的聯(lián)結性質決定。莫爾-庫倫通過試驗發(fā)現(xiàn)土在滑動面上的抗剪強度與法向應力成正比，見式(1)。

τf=σntanφ+c

(1)

式中：τf為土的抗剪強度；c為黏聚力；φ為內摩擦角；σn為滑動面上的正應力。

當土工結構中任意一點的土體在某一方向上的剪應力τ等于土的抗剪強度τf時，就稱該點的土處于極限平衡狀態(tài)，一點土體的極限平衡條件見式(2)。

τ=τf

(2)

當一點土體處于極限平衡狀態(tài)時，該點對應的土體微元體也在對應方向上達到極限平衡狀態(tài)[12]。

1.2 滑動面上土體的極限平衡條件

將一點的極限平衡狀態(tài)擴展到土體沿某一曲線(面)的極限平衡狀態(tài)。假設l為土體區(qū)域內的任意形狀連續(xù)曲線，土體沿l達到極限平衡狀態(tài)是指在曲線任意微元長度上，沿曲線切線方向土體的滑動力(剪應力)與阻滑力(抗剪力)相等。

圖1 滑動面上一微元體沿切線方向的力和力矩的平衡Fig.1 The equilibrium condition in the sliding surfacealong the tangent direction

如圖1所示，若土體沿滑動面任意一點都處于極限平衡狀態(tài)，見式(3)。

τi-τfi=0

(3)

由式(3)可知，滑動面上每一點對應的微元長度上土體的滑動力與阻滑力相平衡，得式(4)。

(4)

此時，整個滑動面上土體滑動力的合力與阻滑力的合力也相平衡(式(5))，對于滑動面外任意一點，滑動力矩與阻滑力矩的合力矩也平衡，見式(6)。

(5)

(6)

任意一個曲面上，如果任意一點的土體都處于極限平衡狀態(tài)，則以該曲面為底的任意形狀的土體都處于極限平衡狀態(tài)，稱為土體沿著該曲面達到極限平衡狀態(tài)。邵龍?zhí)兜萚11]證明，土體沿著曲面達到極限平衡的充分必要條件見式(7)，證明過程見圖2。

(7)

圖2 曲面上土體極限平衡的充分必要條件證明Fig.2 Proof of the necessary and sufficient conditions for the soil’s limit equilibrium on the curved surface

1.3 沿曲面滑動穩(wěn)定安全系數(shù)的定義

對于一個穩(wěn)定的土體結構，人為地減小土體的抗剪強度使其達到極限平衡狀態(tài)。假設R(l)為沿曲面(線)l使土體各點均達到極限平衡狀態(tài)的強度折減系數(shù)函數(shù)，l為土體結構內的任意曲面(線)，那么土體沿著曲面(線)l達到極限平衡的充要條件見式(8)。

(8)

設定l為閉區(qū)間，在剪應力方向不變的條件下可以應用積分中值定理，得式(9)。

(9)

土體沿該曲面的安全系數(shù)可以表示為式(10)。

(10)

同樣地，假設P(l)為沿曲面(線)l使土體各點均達到極限平衡狀態(tài)的超載系數(shù)函數(shù)，那么土體沿著曲面(線)l達到極限平衡的充要條件為式(11)。

(11)

此時中值定理可以表示為式(12)。

(12)

由此可知，安全系數(shù)K不僅是土體在曲面上達到極限平衡的抗剪強度強度折減系數(shù)的中值，也是使曲面上土體沿著曲面達到極限平衡的超載系數(shù)的中值。曲面(線)l既可以是貫穿巖土結構的整體曲面，也可以是在巖土結構內部的局部曲面，說明用此安全系數(shù)作為評價標準，既可以分析土體結構整體的滑動穩(wěn)定性，也可以分析土體的局部破壞區(qū)域[12-13]。

2 安全系數(shù)和最危險滑動面的求解

根據(jù)上述理論，在實際計算中，土體沿某一滑動面的安全系數(shù)可以表示為式(13)。

(13)

將滑動面離散成m個微小單元i，安全系數(shù)可以表示為式(14)。

(14)

采用一維線性單元平面坐標轉換與二階高斯積分，安全系數(shù)可以進一步表示為式(15)。

(15)

式中，|J|為Jacobian矩陣行列式，其表達式為式(16)。

(16)

對于任意土體微元，在確定了高斯積分點坐標值以后，首先需要判斷該高斯點在哪一個有限元單元網(wǎng)格內，然后再根據(jù)距離倒數(shù)加權平均法得到該點應力，并結合土體微元的方向得到高斯積分點沿土體微元方向的剪應力與抗剪強度，根據(jù)式(15)可計算出土體沿該滑動面的安全系數(shù)。

3 應用算例

3.1 均質土坡穩(wěn)定性分析

采用費康等[14]的算例，某均質邊坡高H=10 m，邊坡傾角為45°(圖3)。土體容重為20 kN/m3，內摩擦角φ=20°，黏聚力c=12.38 kPa，彈性模量E=100 MPa，泊松比v=0.35。有限元計算利用商業(yè)軟件Abaqus6.14，采用理想彈塑性本構模型進行計算。模型底部為固定約束，左右邊界為水平約束。

圖4為有限元極限平衡法分析結果，圖4中邊坡區(qū)域內線段表示安全系數(shù)為1.00～1.01，1.01～1.02，1.02～1.03，1.03～1.04，1.04～1.05的局部滑動面線段。從圖4中可以看出，安全系數(shù)臨近1.0的線段主要集中在x=7到x=15之間，說明該范圍內的土體安全系數(shù)較低。

表1為瑞典圓弧法、BISHOP法、有限元強度折減法、有限元極限平衡法計算得到的安全系數(shù)。 圖5為有限元極限平衡法計算得到的臨界滑動面上土體的剪應力和抗剪強度分布情況。從圖5中可以看出， 坡體剪應力分布符合一般規(guī)律。 在坡體表面(x=16.4處)，剪應力接近于0，隨滑動面土體埋深的增加，土體剪應力與抗剪強度增加。在x=11.7處，滑動面埋深達到最大，同時剪應力與抗剪強度達到峰值，隨后滑動面埋深變化不大，剪應力與抗剪強度變化不大。在坡腳處(x=2處)滑動面埋深雖為0，但其沿臨界滑動面的剪應力并不為0，說明坡腳處承受一定的水平力作用，建議在坡腳處增加抗滑措施。在x=7到x=15之間的滑動面安全系數(shù)較低，這也與局部穩(wěn)定性分析得到的結論吻合。

圖3 均質邊坡示意圖Fig.3 Schematic diagram of homogeneous slope

圖4 計算得到的局部破壞區(qū)與臨界滑動面Fig.4 Calculated local failure area and critical sliding surface

表1 安全系數(shù)比較Table 1 Comparison of safety factors

3.2 土壓力分析

某混凝土重力式擋土墻高H=3 m，彈性模量E=20 GPa，泊松比v=0.2。墻后為干砂，彈性模量E=50 MPa，泊松比v=0.3，內摩擦角φ=37°，黏聚力c=0 kPa，重度為17 kN/m3。假定墻體為理想的剛性體，不會產(chǎn)生破壞。

基于RANKINE土壓力理論得到的主動土壓力為19.125 kPa，滑動面破裂角為63.5°，墻后土體在該主動土壓力作用下處于極限平衡狀態(tài)，即沿該滑動面的安全系數(shù)為1.0。 為了便于展開臨界滑動面的搜索，在新的有限元計算模型中，撤去墻體結構，依據(jù)理論土壓力大小和分布形式對墻后土體施加水平梯度載荷及重力載荷，其中主動土壓力大小和分布依據(jù)RANKINE土壓力理論計算得到，底部為固定約束，右邊界為水平簡支約束，計算模型見圖7。

圖5 臨界滑動面上剪應力和抗剪強度分布Fig.5 Shear stress and shear strength distributionon the critical sliding surface

圖6 重力式擋土墻結構示意圖Fig.6 Schematic diagram of gravity retainingwall structure

圖7 有限元計算模型Fig.7 Finite element calculation model

圖8為局部與整體穩(wěn)定性分析結果，從圖8中可以看出局部破壞區(qū)域形狀接近楔體，并貫穿整個填土，說明土體處于極限平衡狀態(tài)，符合主動土壓力特征。利用有限元極限平衡法計算程序搜索得到的最危險滑動面趨近于直線，安全系數(shù)為1.012，滑動面與水平方向夾角接近63.5°，與RANKINE土壓力理論的滑裂面方向基本一致，但是臨界滑動面端點并非位于墻角處，可能是底部約束影響墻角處的應力場導致的。

圖8 計算得到的局部破壞區(qū)與臨界滑動面Fig.8 Calculated local failure area and critical sliding surface

圖9為有限元極限平衡法計算得到的臨界滑動面上土體的剪應力和抗剪強度分布情況。由圖9可知，在擋土墻墻角附近，土體雖承受較大的水平壓力作用，但其沿臨界滑動面的剪應力并不是最大的。整體而言，土體沿滑動面的剪應力與抗剪強度相差不大，表明土體失穩(wěn)，沿滑動面達到極限平衡狀態(tài)。

圖9 臨界滑動面上剪應力和抗剪強度分布Fig.9 Shear stress and shear strength distribution onthe critical sliding surface

3.3 地基承載力分析

采用劉士乙[15]的算例，某均質地基深15 m，寬30 m，彈性模量E=30 MPa，泊松比v=0.3(圖10)。 土體容重為12 kN/m3，內摩擦角φ=10°，黏聚力c=10 kPa，地基中部作用有2 m寬的均布載荷P分別為40 kPa、60 kPa、75 kPa、85 kPa，計算地基在不同載荷作用下的局部破壞區(qū)。

圖11為有限元極限平衡程序在不同載荷下得到的局部破壞區(qū)。由圖11可以看出，在載荷作用下，地基主動區(qū)下方的土體首先發(fā)生剪切破壞，然后破壞區(qū)逐漸增大并向土體兩側擴展，最終載荷為85 kPa時局部破壞區(qū)貫穿至地基表面。該算例應用PRANDTL理論計算得到的地基承載力為82.8 kPa，與有限元極限平衡法得到的結果差別不大。分別向地基兩側搜索，有限元極限平衡法計算程序得到的兩條臨界滑動面匯總如圖12所示，安全系數(shù)平均為1.017，與劉士乙[15]得到的臨界滑動面基本一致。

3.4 尾礦庫穩(wěn)定性分析

圖13是采用多管放礦上游式筑壩的尾礦庫。 初期壩采用庫區(qū)開采的巖石筑壩，壩頂標高為135.0 m，上游坡與下游坡坡比均為1∶2，壩頂寬6.0 m。 子壩采用尾砂筑壩，外坡為1∶5，尾礦庫最終堆筑標高為250 m。尾礦壩各部分材料參數(shù)見表2。

圖10 地基結構示意圖Fig.10 Schematic diagram of foundation structure

圖11 不同荷載下地基的局部破壞區(qū)Fig.11 Local failure zone of foundation under different loads

圖12 計算得到的臨界滑動面Fig.12 Critical sliding surface calculated byfinite element limit equilibrium method

對圖13所示的斷面，采用有限元程序進行有限元應力應變計算，可以得到壩體計算區(qū)域的應力分布。為便于計算，取尾礦庫尺寸的十分之一進行建模。有限元計算模型如圖14所示。啟用自主開發(fā)的巖土結構穩(wěn)定分析的有限元極限平衡法計算程序，自動讀入有限元計算得到的應力場信息，在操作界面輸入有限元極限平衡法分析所需的參數(shù)。分析結果表明，尾礦庫內未出現(xiàn)局部破壞區(qū)域，計算得到的整體臨界滑動面及其安全系數(shù)如圖15所示。圖16為采用BISHOP法和瑞典圓弧法計算得到的臨界滑動面及其安全系數(shù)，圖17為有限元強度折減法計算終止時的增量位移圖。由圖15～圖17可以看出，基于有限元極限平衡法得到的臨界滑動面與其他方法計算得到的臨界滑動面位置基本一致。

表3給出了各種不同方法計算得到的安全系數(shù)。從表3可以看出，4種計算方法得到的臨界滑動面差別不大，基于有限元應力場計算得到的安全系數(shù)略大于傳統(tǒng)極限平衡法計算的結果，整體安全系數(shù)較高。改變荷載工況，可以分析壩體局部和漸進破壞過程，由此可以預估破壞區(qū)域和滑動危險性。

圖13 尾礦壩示意圖Fig.13 Sketch of the tailings dam

表2 尾礦壩材料參數(shù)Table 2 Material parameters of tailings pond

圖14 有限元計算模型Fig.14 Finite element calculation model

圖15 計算參數(shù)輸入和計算結果顯示Fig.15 Parameter input and calculation results display

圖16 BISHOP法和瑞典圓弧法得到的臨界滑動面及其安全系數(shù)Fig.16 Critical sliding surface and safety factor obtained by BISHOP and Ordinary method

4 結 論

本文對傳統(tǒng)極限平衡理論進行了擴展，在討論土體內一點的極限平衡狀態(tài)與土體沿著任意曲面的極限平衡狀態(tài)的基礎上，闡釋了以曲面為底的任意土體沿著曲面達到極限平衡的充分必要條件，由此定義了土體沿(整體或局部)曲面的安全系數(shù)，將此理論用于土體學三大經(jīng)典問題的分析，并與其他理論得到的結果進行對比，得到以下結論如下所述。

圖17 有限元強度折減法計算終止時的增量位移圖Fig.17 Incremental displacement diagram at the end of the FEM strength reduction

表3 安全系數(shù)比較Table 3 Comparison of safety factors

1) 對于邊坡穩(wěn)定性問題，有限元極限平衡法計算得到的臨界滑動面上土體剪應力與抗剪強度分布符合實際特征，其位置和安全系數(shù)與瑞典圓弧法、簡化BISHOP法、有限元強度折減法得到的結果相差不大。

2) 對于擋土墻問題，將RANKINE理論計算得到的主動土壓力直接施加至墻后土體用于有限元計算，有限元極限平衡法得到的局部破壞區(qū)形狀接近楔形體，并貫穿整個填土區(qū)，搜索得到的臨界滑動面接近于直線，與RANKINE土壓力理論的滑動面方向一致，計算得到的安全系數(shù)接近于1。但是由于底部約束的影響，臨界滑動面端點不在墻角處，而是位于墻角上方不遠處。

3) 對于地基承載力問題，采用逐漸增加地基載荷的方式，得到了地基的漸進破壞過程。結果表明，在極限載荷作用下，地基土體RANKINE主動區(qū)下方土體首先發(fā)生剪切破壞，然后破壞區(qū)逐漸增大,最終貫穿至地基表面，得到的極限載荷與PRANDTL理論計算的結果很接近。