線性與非線性強度的土石壩壩坡穩(wěn)定分析下限法

周建烽，王均星，陳 煒，羅貝爾

（1.武漢大學 水資源與水電工程科學國家重點實驗室，湖北 武漢 430072；2.福建省水利規(guī)劃院，福建 福州 350001）

1 引 言

目前，土石壩穩(wěn)定分析在工程上廣泛采用的是剛體極限平衡法[1－4]。剛體極限平衡法概念清晰、計算簡便，工程實踐經(jīng)驗比較豐富。但它假定失穩(wěn)塊體為整體移動的剛體，只考慮滑移體上力的平衡，不考慮結(jié)構(gòu)的屈服，不能確定相應(yīng)的應(yīng)力和位移分布。為求改進，國內(nèi)外眾多學者都在尋求一種更有效的穩(wěn)定分析方法，其中最受關(guān)注的是有限元法和塑性極限分析方法。在有限元法中，一般采用有限元強度儲備系數(shù)法，它可以不作任何假定，能夠考慮土體的本構(gòu)關(guān)系，能夠模擬邊坡的破壞過程及其滑移面形狀，并且可以很方便地考慮土體和錨固、支護的共同作用。采用有限元強度儲備系數(shù)法進行邊坡穩(wěn)定分析是近年來的熱點，國內(nèi)外許多學者都進行了有益的探索[5－9]，尤其是可以利用它來進行邊坡失穩(wěn)的仿真計算。

塑性極限分析方法避開彈塑性分析的全過程，直接研究結(jié)構(gòu)的極限狀態(tài)，求解結(jié)構(gòu)的極限荷載，為評價結(jié)構(gòu)的安全度提供必要的數(shù)據(jù)，其理論嚴謹，在土工問題極限承載能力上是一種十分有效的方法[10]。自引入有限元法解決了構(gòu)造靜力場和機動場的難題之后，結(jié)合數(shù)學規(guī)劃手段，已成功地用于計算地基承載力及邊坡穩(wěn)定分析等方面。

本文采用塑性極限分析方法進行土石壩的壩坡穩(wěn)定分析，在Sloan[11]的工作基礎(chǔ)上，考慮了滲流作用和地震荷載，推導(dǎo)了基于非線性規(guī)劃的土石壩有限元塑性極限分析下限法模型，并采用迭代算法對堆石、砂土等無黏聚土進行非線性強度指標的壩坡穩(wěn)定計算。最后，對幾個典型土坡和具體的土石壩工程進行了計算，并與多種方法的分析結(jié)果相比較。

2 塑性極限分析的下限原理

在物體及其邊界上，滿足下列條件的應(yīng)力場稱為靜力容許應(yīng)力場：

（3）在力邊界上，荷載與真實荷載相同。

靜力容許應(yīng)力場有多個，每個靜力容許應(yīng)力場對應(yīng)著一個外荷載。因為結(jié)構(gòu)并未破壞，因而這個荷載一定比極限荷載小。下限定理可以表述為：在所有的與靜力容許應(yīng)力場對應(yīng)的荷載中，最大的荷載為極限荷載，與最大荷載相應(yīng)的靜力容許應(yīng)力場為結(jié)構(gòu)的極限狀態(tài)應(yīng)力場。

3 土石壩壩坡穩(wěn)定的有限元極限分析下限法模型

土石壩壩坡穩(wěn)定的有限元塑性極限分析下限法非線性規(guī)劃數(shù)學模型由三部分組成：①非線性約束條件，即屈服條件；②線性約束條件：包括平衡條件、應(yīng)力間斷條件和邊界條件；③目標函數(shù)，即安全系數(shù)。本文采用線性三角形單元離散結(jié)構(gòu)體，每個單元都需滿足上述的約束條件與目標函數(shù)。

3.1 結(jié)構(gòu)的離散

為構(gòu)造靜力許可應(yīng)力場，結(jié)構(gòu)物采用不共節(jié)點的三角形線性單元離散，以節(jié)點應(yīng)力為未知量?？紤]滲流作用后，每個節(jié)點都增加一個水頭變量Hi。

圖1 下限分析法單元離散圖Fig.1 Element discrete of lower bound limit analysis

在下限法中，單元內(nèi)任意一點的應(yīng)力分量、水頭可表示為節(jié)點應(yīng)力分量、節(jié)點水頭的線性函數(shù)，即

式中：A為三角形單元的面積，ξi、ηi、ζi的具體表達式參見一般有限元著作。

3.2 平衡約束條件

土石壩穩(wěn)定計算要考慮多種荷載的影響，其中最主要的外荷載為孔隙水壓力和地震荷載。

對于孔隙水壓力的考慮，采用有效應(yīng)力原理，根據(jù)水土分算原則，即把土骨架當作有限元極限分析法中的研究對象，把孔隙水壓力當作是一種外力作用于土骨架上。本文在考慮滲流作用時，將其視做體積力考慮并表示為

式中：γw為水的重度；J為水力坡降。

可以把Fp分解為兩個分量：

式中：Fpx、Fpy分別為滲流作用力在x、y 方向的分量； Jx、 Jy分別為水力坡降在x、y 方向的分量。

土石壩用動力法作抗震計算時，由于壩料的非線性特性使得計算仍然十分復(fù)雜。迄今為止，對于土石壩的地震穩(wěn)定性和永久變形的動力方法和準則在工程界仍然沒有共同的看法。我國的水工建筑物抗震規(guī)范[12]規(guī)定土石壩抗震穩(wěn)定計算可采用擬靜力法。擬靜力法將地震荷載作為靜荷載來檢驗結(jié)構(gòu)的強度，在靜力法的形式中，納入了動力分析法的內(nèi)容，它具有簡便、實用的優(yōu)點。因此，本文在有限元極限分析下限法中，采用了擬靜力法來考慮地震荷載。

下限法平衡方程表示為

式中：kh、kv分別為水平和垂直地震加速度代表值；γ′為土體重度（浸潤線以下取浮重度，浸潤線以上取天然重度）；γ為包含孔隙水在內(nèi)的土體重度，應(yīng)力拉正、壓負。

平衡條件寫成矩陣形式如下：

式中：上標e表示單元號；且有

對于其他線性約束條件即應(yīng)力間斷條件和邊界條件的建立，可以參考文獻[13]，本文不再贅述。

3.3 屈服約束條件

對于一般受力下的巖土材料，通常采用Mohr-Coulomb屈服準則。在有限元極限分析下限法模型中，對于平面應(yīng)變問題，每個單元節(jié)點應(yīng)滿足屈服條件，可表示為

式中：i=1，2，3，e=1，2，…，NE，NE為總單元數(shù)，分別為第e 個單元的3個應(yīng)力分量、材料的黏聚力和內(nèi)摩擦角。

3.4 目標函數(shù)及非線性數(shù)學規(guī)劃模型

本文采用強度儲備系數(shù)來定義安全系數(shù)，即

式中：φ0、c0為進行強度折減以后的強度參數(shù)。

式（13）可簡寫為

將強度儲備系數(shù)設(shè)為目標函數(shù)的下限法數(shù)學規(guī)劃模型為

對上述形成的土石壩壩坡穩(wěn)定分析下限法的非線性數(shù)學規(guī)劃模型求解，本文采用專業(yè)優(yōu)化軟件Lingo8.0，通過非線性規(guī)劃求解最終可以獲得土石壩壩坡抗滑安全系數(shù)值及相應(yīng)的靜力許可應(yīng)力場。

4 非線性強度指標的壩坡穩(wěn)定計算

高土石壩一般多采用堆石料，對此類材料在高土石壩的高應(yīng)力水平條件下采用線性強度指標已不能滿足要求，有必要采用非線性的強度指標。

本文在考慮材料強度的非線性時，采用Duncan對數(shù)模式為

式中：φ為土體滑動摩擦角（°）；φ0為一個大氣壓下的摩擦角； Δφ為增加一個對數(shù)周期下φ 的減小值；σ1為土體的大主應(yīng)力值（本文以拉應(yīng)力為正）；Pa為大氣壓力。

該模式在我國有十分豐富的應(yīng)用經(jīng)驗，碾壓式土石壩設(shè)計規(guī)范規(guī)定對于粗粒料非線性抗剪強度指標采用該模式。

在下限法模型中，本文采用了迭代算法進行非線性強度指標的壩坡穩(wěn)定計算，具體步驟如下：

（1）擬定一個初始內(nèi)摩擦角 φ0，一般取φ0。

（2）用本文第3節(jié)中的方法進行計算，得到第j 次迭代的安全系數(shù) kj及應(yīng)力場{σj}，并獲得第一主應(yīng)力

（3）由式（16）計算出各單元內(nèi)摩擦角φe。

（4）利用步驟（3）的結(jié)果重復(fù)步驟（2），得到第j+1次迭代的安全系數(shù)kj+1及應(yīng)力場{σj+1}，驗算前、后兩次安全系數(shù)是否滿足，ε為收斂精度。

（5）如果滿足收斂精度，則停止迭代，得到的安全系數(shù)kj+1即為所求的安全系數(shù)。如不滿足則重復(fù)步驟（2）～（4）。

5 算例分析

5.1 含軟弱夾層的算例

圖2所示為一含有軟弱夾層的土坡，F(xiàn)redlund[13]、Kim[14]等分別采用剛體極限平衡法及基于線性規(guī)劃的極限分析法對此算例進行過分析。該土坡頂寬L=18.3 m，高H=12.2 m，水頭高度HW=6.1 m，坡比為1:2，均質(zhì)土體材料參數(shù)為c′=28.76 kPa，γ=18.89 kN/m3，軟弱夾層材料參數(shù)c′=0 kPa，φ′=10°。

本文采用3種不同密度的網(wǎng)格（74個單元、156個單元、238個單元，如圖3所示）進行計算，計算結(jié)果見表1。由本文方法還可以得到極限狀態(tài)下的應(yīng)力場。由于篇幅原因，此處只給出密網(wǎng)格工況3、4的第2主應(yīng)力等值線圖，詳見圖4、5。由主應(yīng)力等值線圖可看出，在孔隙水壓力和土體自重作用下，主應(yīng)力大小隨著深度加深而線性增加，符合一般自重應(yīng)力場的分布規(guī)律。

從表1的計算結(jié)果可知，安全系數(shù)隨著網(wǎng)格的加密而增加，并且當單元數(shù)達一定程度時，增加網(wǎng)格數(shù)量提高的精度并不明顯。本文方法在網(wǎng)格密到一定程度時計算的結(jié)果比瑞典法高2%～5%左右，比簡化畢肖普法和摩根斯頓-普賴斯法小3%～5%左右。從塑性極限分析的觀點來看，剛體極限平衡法考慮了力的平衡條件，在滑動面上滿足屈服條件，其解答可歸為下限法的范疇，但并不能判斷其他部位是否滿足屈服條件，而且瑞典法就連平衡條件也沒有完全滿足，因而不是嚴格的下限解；另外，它假設(shè)了滑動面，即假設(shè)了結(jié)構(gòu)破壞的一種機構(gòu)，并且從一系列破壞機構(gòu)中尋找安全系數(shù)的最小值作為最終解，其解答也可歸為上限法的范疇，但它并沒有考慮變形協(xié)調(diào)關(guān)系，事實上也不是一種完全的上限解，因而由極限平衡法得到的結(jié)果既非下限解也非上限解。

由此可以說明，本文方法得到的結(jié)果為嚴格的下限解，是偏于安全的。本文結(jié)果比Kim等[15]計算的下限解略高，這是因為Kim在對屈服函數(shù)的線性化時采用內(nèi)接正多邊形對屈服圓進行逼近，根據(jù)下限定理得到的解必然比本文方法小。

圖2 含有軟弱夾層的土坡示意圖Fig.2 Sketch of slope with a weak interlayer

圖3 有限元網(wǎng)格劃分Fig.3 Meshes of finite elements

表1 多種方法的安全系數(shù)Table 1 Factor of safety for slopes obtained with different methods

圖4 工況3第2主應(yīng)力圖（單位:kPa）Fig.4 Contours of second principal stress in slope at condition 3 (unit:kPa)

圖5 工況4第2主應(yīng)力圖（單位：kPa）Fig.5 Contours of second principal stress in slope at condition 4 (unit:kPa)

5.2 線性和非線性強度指標的土壩壩坡穩(wěn)定計算

某土壩壩坡坡角為α，坡比為1:2，壩料的線性強度參數(shù)c=0，φ=45°，非線性強度參數(shù)為c=0，φ0=55°，Δφ=10°。本文計算不同 γ H值下的安全系數(shù)，以獲得安全系數(shù)隨 γ H變化規(guī)律，計算結(jié)果見表2。圖6、7分別為 γ H=4 000 kPa情況下的第1和第2主應(yīng)力圖（本文以拉應(yīng)力為正）。由圖可看出，主應(yīng)力（有效應(yīng)力）等值線分布良好，符合一般規(guī)律。

由表2中結(jié)果可知：

（1）當采用線性強度指標計算時，本文方法與瑞典圓弧法和簡化畢肖普法計算出的安全系數(shù)極為接近，皆近似等于tan φ/tanα，由此可以看出，本文方法計算結(jié)果是可靠的。

（2）采用非線性強度指標的計算結(jié)果兩者相差較大，達15%左右，并且在 γ H較大時，本文算出的結(jié)果比后者要??；在 γ H較小時，本文算出的結(jié)果反而大。究其原因，剛體極限平衡法無論是進行線性還是非線性計算，僅能計算滑動面上的應(yīng)力，只滿足滑動面上的屈服條件。其計算結(jié)果既不是下限解也不是上限解。

表2 線性和非線性強度指標下的壩坡安全系數(shù)Table 2 Dam slope factors of safety with linear and nonlinear strength indexes

圖6 第一主應(yīng)力圖（單位：kPa）Fig.6 Contours of the first principal stress in slope (unit:kPa)

圖7 第二主應(yīng)力圖（單位：kPa）Fig.7 Contours of the second principal stress in slope (unit:kPa)

5.3 土石壩穩(wěn)定分析算例

土石壩規(guī)范[16]規(guī)定土石壩穩(wěn)定計算應(yīng)分別對施工期（包括竣工期）、穩(wěn)定滲流期、水庫水位降落期和正常運用遇地震4種工況進行計算。由于本文研究目的在于驗證塑性極限分析這一方法在土石壩穩(wěn)定分析中的可行性，所以只分別對心墻壩的穩(wěn)定滲流期遇地震工況、斜心墻壩的正常蓄水位遇地震工況進行計算。至于其他工況，不同之處只是在坡外水水頭、材料參數(shù)和孔隙水壓力取值上的差別。

5.3.1 心墻壩算例

某心墻壩壩體斷面圖和材料分區(qū)如圖8所示，壩體力學參數(shù)見表3，壩高為32.5 m，上游水位為27.5 m，下游無水，水的重度取為9.81 kN/m3，地震烈度為Ⅷ度。

表3 土體力學參數(shù)Table 3 Mechanical parameters of soils

計算模型中共劃分了738個單元，共有2 214個不共用節(jié)點，1 069條公共邊（見圖9）。由下限法計算得到的壩坡抗滑安全系數(shù)見表4（本文方法考慮的是土石壩上下游整體，不是某個方向的滑動），單元主應(yīng)力矢量見圖10。

圖8 心墻壩斷面示意圖Fig.8 Sketch of core wall dam section

圖9 心墻壩有限元網(wǎng)格劃分Fig.9 Meshes of finite element of core wall dam

表4 心墻壩壩坡抗滑安全系數(shù)Table 4 Factors of safety for core wall dam

圖10 心墻壩主應(yīng)力矢量圖Fig.10 Vectogram of principal stresses of core wall dam

5.3.2 斜心墻壩算例

某斜心墻壩壩體斷面圖和材料分區(qū)如圖11所示，壩體及壩基材料物理力學性質(zhì)指標見表5。水的重度取為9.81 kN/m3，地震烈度為Ⅵ度。

模型共劃分了901個單元，共有2 703個不共用節(jié)點、1 311條公共邊（見圖12）。由下限法計算得到的壩坡抗滑安全系數(shù)見表6，單元主應(yīng)力矢量見圖13。

圖11 斜心墻壩斷面示意圖Fig.11 Sketch of inclined core wall dam section

表5 壩體及壩基物理力學性質(zhì)指標表Table 5 Physico-mechanical parameters of dam body and foundation

圖12 斜心墻壩有限元網(wǎng)格示意圖Fig.12 Meshes of finite element of inclined core wall dam

表6 斜心墻壩壩坡抗滑安全系數(shù)Table 6 Factors of safety for inclined core wall dam

圖13 斜心墻壩主應(yīng)力矢量圖Fig.13 Vectogram of principal stresses of inclined core wall dam

由心墻和斜心墻土石壩計算結(jié)果可知，由下限法算出的壩坡抗滑安全系數(shù)結(jié)果與瑞典圓弧法和簡化畢肖普法基本相符，主應(yīng)力分布規(guī)律與常規(guī)有限元法的結(jié)果基本一致，自重應(yīng)力為主要分布，符合一般的應(yīng)力分布規(guī)律。

6 結(jié) 論

（1）本文在Sloan的工作基礎(chǔ)上，將有限元塑性極限分析下限法用于土石壩壩坡穩(wěn)定分析，理論嚴格且物理意義明確，可求得壩坡抗滑安全系數(shù)的下限解及相應(yīng)的靜力許可應(yīng)力場。通過幾個典型土坡和具體的土石壩工程的分析發(fā)現(xiàn)，本文方法的結(jié)果與極限平衡法相符，并且無論是使用線性強度指標還是非線性強度指標進行計算，都能獲得較好的結(jié)果，證明了本文方法的正確性。

（2）剛體極限平衡法在一定意義上滿足了靜力平衡條件，假定了破壞機制，一定程度上滿足了上下限定理的一些條件，但并不滿足完全屈服條件和變形協(xié)調(diào)條件，其計算結(jié)果既不是下限解也不是上限解。

（3）本文方法具有以下優(yōu)點：①不用假定滑動面，避免了由于滑動面位置和形狀帶來的計算誤差，并在解決了規(guī)劃問題的前提下，可以對三維問題進行分析；②能夠較方便地處理復(fù)雜幾何形狀和非均質(zhì)材料問題；③回避了結(jié)構(gòu)材料復(fù)雜的本構(gòu)關(guān)系；④可以考慮各種外荷載。

（4）目前，塑性極限分析在水工領(lǐng)域應(yīng)用還比較少。為發(fā)揮其自身優(yōu)勢使其今后能夠在水工方面得到推廣，今后還需做的工作是：①通過對大量的實際工程進行計算分析，與現(xiàn)行的各種計算方法進行比較后，確定與規(guī)范相對應(yīng)的下限法壩坡抗滑安全系數(shù)的取值范圍；②由于求解是結(jié)合優(yōu)化方法進行的，尋求一種有效的求解大規(guī)模稀疏矩陣的算法和求解器是必要的；③結(jié)合上限法構(gòu)造一套力學概念清晰、計算方法先進、計算結(jié)果有效全面的土壩壩坡穩(wěn)定分析方法。

[1]BISHOP A W.The use of the slip circle in stability analysis of slope[J].Geotechnique，1955，5(1):7－17.

[2]MORGENSTERN N R，PRICE V E.The analysis of the stability of general slip surfaces[J].Geotechnique，1965，15(1):79－93.

[3]SPENCER E.A method of analysis of the stability of embankments assuming parallel inter-slice forces[J].Geotechnique，1967，17(1):11－26.

[4]SARMA S K.Stability analysis of embankments and slopes[J].Geotechnique，1973，23(3):423－433.

[5]TAMOTSU，MASTUI，SAN KA CHING.Finite element slope stability analysis by shear strength reduction technique[J].Soils and Foundations，1992，32(1):59－70.

[6]連鎮(zhèn)營，韓國城，孔憲京.強度折減有限元法研究開挖邊坡的穩(wěn)定性[J].巖土工程學報，2001，23(4):407－411.LIAN Zhen-ying，HAN Guo-cheng，KONG Xian-jing.Stability analysis of excavation by strength reduction FEM [J].Chinese Journal of Geotechnical Engineering，2001，23(4):407－411.

[7]鄭穎人，趙尚毅.有限元強度折減法在土坡與巖坡中的應(yīng)用[J].巖石力學與工程學報，2004，23(19):3381－3388.ZHENG Ying-ren，ZHAO Shang-yi.Application of strength reduction FEM to soil and rock slope[J].Chinese Journal of Rock Mechanics and Engineering，2004，23(19):3381－3388.

[8]陳國慶，黃潤秋，周輝，等.邊坡漸進破壞的動態(tài)強度折減法研究[J].巖土力學，2013，34(4):1140－1146.CHEN Guo-qing，HUANG Run-qiu，ZHOU Hui，et al.Research on progressive failure for slope using dynamic strength reduction method[J].Rock and Soil Mechanics，2013，34(4):1140－1146.

[9]李小春，袁維，白冰，等.基于局部強度折減法的邊坡多滑面分析方法及應(yīng)用研究[J].巖土力學，2014，35(3):847－854.LI Xiao-chun，YUAN Wei，BAI Bing，et al.Analytic approach of slope multi-slip surfaces based on local strength reduction method and its application[J].Rock and Soil Mechanics，2014，35(1):847－854.

[10]陳惠發(fā)，詹世斌.極限分析與土體塑性[M].北京:人民交通出版社，1995.

[11]SLOAN S W.Lower bound limit analysis using finite elements and linear programming[J].International Journal for Numerical and Analytical Methods in Geomechanics，1988，12(1):61－77.

[12]中華人民共和國行業(yè)標準編寫組.SL203－97水工建筑物抗震設(shè)計規(guī)范[S].北京:中國水利水電出版社，1997.

[13]王漢輝，王均星，王開治.邊坡穩(wěn)定的有限元塑性極限分析[J].巖土力學，2003，24(5):733－738.WANG Han-hui，WANG Jun-xing，WANG Kai-zhi.Plastic limit analysis of slope stability using finite element[J].Rock and Soil Mechanics，2003，24(5):733－738.

[14]FREDLUND D G，KRAHN J.Comparison of slope stability methods of analysis[J].Canadian Geotechnical Journal，1977，14(3):429－439.

[15]KIM J，SALGADO R，LEE J.Stability analysis of complex soil slopes using limit analysis[J].Journal of Geotechnical and Geoenvironmental Engineering，2002，128(7):546－557.

[16]中華人民共和國行業(yè)標準編寫組.SL274－2001碾壓式土石壩設(shè)計規(guī)范[S].北京:中國水利水電出版社，2002.