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

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

1 引 言

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

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

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

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

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

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

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

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

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

3.1 結構的離散

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

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

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

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

3.2 平衡約束條件

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

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

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

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

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

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

下限法平衡方程表示為

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

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

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

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

3.3 屈服約束條件

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

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

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

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

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

式（13）可簡寫為

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

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

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

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

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

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

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

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

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

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

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

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

（5）如果滿足收斂精度，則停止迭代，得到的安全系數(shù)kj+1即為所求的安全系數(shù)。如不滿足則重復步驟（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所示）進行計算，計算結果見表1。由本文方法還可以得到極限狀態(tài)下的應力場。由于篇幅原因，此處只給出密網(wǎng)格工況3、4的第2主應力等值線圖，詳見圖4、5。由主應力等值線圖可看出，在孔隙水壓力和土體自重作用下，主應力大小隨著深度加深而線性增加，符合一般自重應力場的分布規(guī)律。

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

由此可以說明，本文方法得到的結果為嚴格的下限解，是偏于安全的。本文結果比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主應力圖（單位:kPa）Fig.4 Contours of second principal stress in slope at condition 3 (unit:kPa)

圖5 工況4第2主應力圖（單位：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ī)律，計算結果見表2。圖6、7分別為 γ H=4 000 kPa情況下的第1和第2主應力圖（本文以拉應力為正）。由圖可看出，主應力（有效應力）等值線分布良好，符合一般規(guī)律。

由表2中結果可知：

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

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

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

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

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

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

土石壩規(guī)范[16]規(guī)定土石壩穩(wěn)定計算應分別對施工期（包括竣工期）、穩(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（本文方法考慮的是土石壩上下游整體，不是某個方向的滑動），單元主應力矢量見圖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 心墻壩主應力矢量圖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，單元主應力矢量見圖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 斜心墻壩主應力矢量圖Fig.13 Vectogram of principal stresses of inclined core wall dam

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

6 結 論

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

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

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

（4）目前，塑性極限分析在水工領域應用還比較少。為發(fā)揮其自身優(yōu)勢使其今后能夠在水工方面得到推廣，今后還需做的工作是：①通過對大量的實際工程進行計算分析，與現(xiàn)行的各種計算方法進行比較后，確定與規(guī)范相對應的下限法壩坡抗滑安全系數(shù)的取值范圍；②由于求解是結合優(yōu)化方法進行的，尋求一種有效的求解大規(guī)模稀疏矩陣的算法和求解器是必要的；③結合上限法構造一套力學概念清晰、計算方法先進、計算結果有效全面的土壩壩坡穩(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]鄭穎人，趙尚毅.有限元強度折減法在土坡與巖坡中的應用[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]李小春，袁維，白冰，等.基于局部強度折減法的邊坡多滑面分析方法及應用研究[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水工建筑物抗震設計規(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碾壓式土石壩設計規(guī)范[S].北京:中國水利水電出版社，2002.