地震作用下土坡穩(wěn)定性的變分極限平衡分析

中圖分類號(hào)：TU435 文獻(xiàn)標(biāo)志碼：A DOI：10.16385/j.cnki.issn.1004-4523.202306032

Abstract：Basedonthebasicprincipleofvariationalmethodandlimitequilibriummethod，thestabiltyofsoillopeunderearth quakeactionisaalyzedacurately.Combiningthelimitequlibriummethodofsopestabilityaalysisandthepseudostaticmeth od，theauxilaryfunctionalunderonstraintconditionsisonstructedbyintroducingLagrangeultiplritothqilibimequtionofslidingsoil.Thefrstorderordinarydiferentialequationswiththebasicunknowsofpotentialslidingsurface，oalstress of slidingsurface，forceofslidingbodysafetyfactorandLagrange muliplierareobtainedbyusing Eulerequation.Thecoupled nonlineardiferentialequationsaresolvednumericallybyusingtheshoting method，andanaccuratesolutionforslopestability analysis under seismic action is obtained.The efectiveness ofthe modeland method is verified by numerical examples.

Keywords： variational limit equilibrium method ;homogeneous soil slope;earthquake action;shooting method

考慮地震作用下的邊坡穩(wěn)定性分析是巖土工程的一個(gè)重要課題。擬靜力法將作用在邊坡上的地震作用等效為靜態(tài)條件下的慣性力，將動(dòng)力學(xué)問(wèn)題簡(jiǎn)化為靜力學(xué)問(wèn)題，由于其物理概念明確，計(jì)算效率高而被廣泛應(yīng)用于地震作用分析中[1-3]。極限平衡方法、強(qiáng)度折減法、極限分析法是目前邊坡穩(wěn)定性分析的常用方法[45]。其中，基于極限平衡理論對(duì)邊坡進(jìn)行分析時(shí)，由于靜力平衡原理分析滑動(dòng)體所得靜力平衡方程絕大部分為超靜定問(wèn)題，所需求解的未知變量往往大于所建立方程的數(shù)量，所以在解決此類問(wèn)題時(shí)都會(huì)引人某種假定，不同簡(jiǎn)化方法采取了不同的假定條件。極限分析法將復(fù)雜巖土材料本構(gòu)關(guān)系簡(jiǎn)化成理想剛塑性狀態(tài)下的應(yīng)力-應(yīng)變關(guān)系，從而得到極限狀態(tài)下的剛塑性體極限破壞荷載，一般包括上限定理和下限定理。解析極限分析法最主要的核心問(wèn)題是需要構(gòu)造一個(gè)符合條件的運(yùn)動(dòng)許可速度場(chǎng)。強(qiáng)度折減法不需要提前假設(shè)邊坡的滑動(dòng)面，可以滿足靜力許可、土體的應(yīng)力-應(yīng)變之間的本構(gòu)關(guān)系等條件，但是難以確定邊坡的初始應(yīng)力狀態(tài)和邊坡接近破壞時(shí)的彈塑性本構(gòu)關(guān)系。由于極限平衡法概念明確，計(jì)算方法簡(jiǎn)便，因此得到廣泛應(yīng)用。然而，傳統(tǒng)的極限平衡法存在條間力假設(shè)、滑動(dòng)面試算迭代過(guò)程繁雜等不足，使得計(jì)算精度和效率受限。為此，KOP'ASCY提出了基于變分原理的邊坡穩(wěn)定性極限平衡分析方法，該法通過(guò)變分原理將安全系數(shù)的求解著作是滑動(dòng)面和應(yīng)力泛函的極值問(wèn)題，具有嚴(yán)謹(jǐn)?shù)臄?shù)學(xué)基礎(chǔ)從而被廣泛關(guān)注[7-9]。

考慮地震作用下邊坡穩(wěn)定性的變分極限平衡分析方面也有一系列的研究成果[10-14]。其中，結(jié)合擬靜力法，LESHCHINSKY等[]考慮對(duì)數(shù)螺旋形破壞機(jī)制的二維土坡，利用變分極限平衡法得到水平地震作用下邊坡的穩(wěn)定性設(shè)計(jì)圖表。BAKER等[13采用變分極限平衡法，考慮滿足邊坡穩(wěn)定性的極值條件，得到了地震作用下問(wèn)題的理論解。基于變分原理，欒茂田等[14]研究了地震作用下的層狀非均質(zhì)土坡的穩(wěn)定性，利用數(shù)值方法討論了土坡的臨界破壞機(jī)制及其極限抗震性能。

由于邊坡穩(wěn)定性分析涉及可動(dòng)的邊界條件，在確定滑面的過(guò)程中通常采用最優(yōu)化搜索算法進(jìn)行數(shù)值求解[15-16]，最優(yōu)化方法適用于復(fù)雜地質(zhì)和地形條件的邊坡，但在進(jìn)行滑動(dòng)面搜索時(shí)當(dāng)問(wèn)題包含較多自由度時(shí)，計(jì)算量較大，計(jì)算有時(shí)也不容易收斂到安全系數(shù)的整體極值上面。本文以水平地震作用下均質(zhì)土坡的穩(wěn)定性為研究對(duì)象，結(jié)合擬靜力法，將變分極限平衡法得到的基本方程轉(zhuǎn)化為具有固定邊界的兩點(diǎn)邊值問(wèn)題，采用打靶法求解該非線性微分方程組，可以直接得到邊坡臨界滑動(dòng)面、滑動(dòng)面上的法向應(yīng)力和條間力分布的精確解，并與已有文獻(xiàn)的分析結(jié)果進(jìn)行對(duì)比，驗(yàn)證模型和方法的有效性。

1變分極限平衡法的數(shù)學(xué)模型

1. 1 基本方程

考慮如圖1所示的在水平地震作用下的均質(zhì)土坡，其潛在滑坡體以邊坡坡面函數(shù) y_p（x） 和潛在滑動(dòng)面函數(shù) y（x） 為邊界。在Mohr-Coulomb強(qiáng)度準(zhǔn)則下引人安全系數(shù) F_s^[6] ，則有：

式中， τ 和 σ 分別為滑動(dòng)面上的切向和法向應(yīng)力； c 為土體的黏聚力； φ 表示土體的內(nèi)摩擦角。

采用擬靜力法，將水平地震作用力等效為靜力條件下的慣性力，則滑動(dòng)體的靜力平衡條件為：

（2b） （2c）式中， k_h 表示水平地震作用系數(shù)； γ_m 表示土體重度；γ_m（y_p-y） 表示土體單位重量； α 為滑動(dòng)面與水平面的夾角；滑動(dòng)面弧長(zhǎng)； x₀ 和 x_n 分別為滑動(dòng)面函數(shù)y（x） 的兩端點(diǎn)。利用幾何關(guān)系 cosα=dx/dl ，sin ，由于靜力平衡滿足坡面上各個(gè)點(diǎn)，故而滿足： H（x₀）=H（x_n）=0 V（x₀）=V（x_n）=0 和 M（x₀）=M（x_n）=0 。通過(guò)引入滑動(dòng)體上水平力 H_? 豎向力 V 以及彎矩 M ，極限平衡方程（2）可表示為：

方程（3）表明，滑動(dòng)面函數(shù) y（x），σ（x） 與參數(shù)F_s 相關(guān)。邊坡的穩(wěn)定分析即在水平地震作用系數(shù) k_h 下尋找2個(gè)未知函數(shù) y（x） 和 σ（x） ，在滿足平衡方程（3）條件下實(shí)現(xiàn)安全系數(shù) F_s 最小，即

由此可知，危險(xiǎn)滑動(dòng)面的最小安全系數(shù)的實(shí)質(zhì)就是求某個(gè)泛函的極值問(wèn)題。為了方便后面的分析和計(jì)算，令：

式（5a）和（5b）即為傳統(tǒng)條分法中的條間力，在FELLENIUS、BISHOP、PRICE和JANBU等的極限平衡分析中，對(duì)于條間力均有不同的假定[17]。

1. 2 變分分析

方程（4）實(shí)際上定義了一個(gè)變分問(wèn)題。根據(jù)變分法的基本原理，考慮3個(gè)極限平衡方程中的1個(gè)方程為目標(biāo)泛函，同時(shí)另2個(gè)方程就必須滿足所需的積分約束，這樣將問(wèn)題可轉(zhuǎn)換成標(biāo)準(zhǔn)的等周問(wèn)題。通過(guò)構(gòu)造含拉格朗日乘子 λ₁ 和 λ₂ 的輔助泛函 F ，可將有積分約束的等周問(wèn)題轉(zhuǎn)換成無(wú)約束的變分極值問(wèn)題，即

F=F₀+λ₁F₁+λ₂F₂

式中，

F₁=c+σtanφ-F_sσy^′-k_c（y_p-y）

該等周問(wèn)題的極值必須滿足歐拉方程[6]，即

由方程（8a）可得滑動(dòng)面函數(shù)的微分方程：

由方程（8b）可得法向應(yīng)力函數(shù)的微分方程：

1.3邊界條件與橫截條件

到目前為止有了約束方程（5）滑動(dòng)面方程（9）和應(yīng)力方程（10），為獲得問(wèn)題的特定解，還需要滿足下列條件。

邊界條件：

橫截條件：

式（6）結(jié)合式（7）（9）～（11），并代入方程（12），經(jīng)過(guò)整理可以得到：

至此，基于變分法的基本原理和極限平衡法得到了在水平地震作用下均質(zhì)土坡穩(wěn)定性問(wèn)題的控制方程及其相應(yīng)的邊界條件。

1.4狀態(tài)方程及打靶法求解

式（9）（10）為滑動(dòng)面函數(shù)與應(yīng)力函數(shù)相互耦合的非線性常微分方程組，一般很難得到其解析解。打靶法將兩點(diǎn)邊值問(wèn)題的數(shù)值解轉(zhuǎn)化為一階方程組初值問(wèn)題的數(shù)值解，通過(guò)尋求滿足原邊界條件的初值問(wèn)題解從而得到能夠滿足任意精度要求的邊值問(wèn)題的精確解，為此類常微分方程組的求解提供有效途徑[18-19]。通過(guò)打靶法將方程（5）、（8）、（9）轉(zhuǎn)化為一階常微分方程組的形式，建立相應(yīng)的狀態(tài)方程。

引入未知變量 W₀ 和 W_n

由于 x₀ 和 X_n 的大小與自變量 x 無(wú)關(guān)，則有：

再引入新自變量 Ψ_t. ，并做下列變量變換：

x=W₀+t（W_n-W₀），t∈[0，1]

以 x 為自變量，那么微分方程組的自由邊界問(wèn)題就可以轉(zhuǎn)換為固定邊界問(wèn)題。引入新的自變量 t 和表示積分弧長(zhǎng)區(qū)間大小的變量 W₀ 和 W_n 后，問(wèn)題就轉(zhuǎn)變?yōu)?t 在 [0，1] 上的兩點(diǎn)邊值問(wèn)題。選取 （y，σ，H，V，M F_s，λ₁，λ₂ W₀ ， W_n ）作為未知量，結(jié)合方程（5）（9）、（10），則問(wèn)題相應(yīng)的狀態(tài)方程可表示為：

利用打靶法可以將邊界方程（18）（19）轉(zhuǎn)化為初值問(wèn)題，從而利用數(shù)值分析方法得到方程組（17）中所有未知量的數(shù)值解，即可以得到滑動(dòng)面位置y（x） 的數(shù)值解。有關(guān)打靶法計(jì)算的具體過(guò)程可參考文獻(xiàn)[19]，這里不再贅述。

2 數(shù)值算例

為了驗(yàn)證本文所提出模型和計(jì)算方法的有效性，選取澳大利亞計(jì)算機(jī)應(yīng)用協(xié)會(huì)（AssociationofComputerAidedDesign，ACADS）中的算例[i6]，當(dāng)水平地震作用系數(shù) k_h=0 時(shí)，本文的計(jì)算模型與該算例相同，即不考慮地震作用的情形。邊坡的計(jì)算模型如圖2所示。取土體的黏聚力 c=3kPa ，內(nèi)摩擦角φ=19.6^° ，土體重度 γ_m=20kN/m³ 。

計(jì)算結(jié)果與不同形式的條分法、強(qiáng)度折減法以及使用搜索算法的變分極限平衡法[1進(jìn)行了比較，結(jié)果如表1所示?？梢钥闯觯?dāng)水平地震作用系數(shù)k_h=0 時(shí)所得結(jié)果0.991與變分極限平衡法和Morgenstern-Price法的安全系數(shù)相等。

圖3給出了水平地震作用系數(shù) k_h=0 時(shí)臨界滑動(dòng)面上法向應(yīng)力分布情況。由圖3可知臨界滑動(dòng)面坡角處的法向應(yīng)力值為 11.23kPa ，而坡頂出現(xiàn)的法向應(yīng)力為一 ?3.13kPa ，表明在坡頂存在張拉應(yīng)力。圖4分別繪出了采用Morgenstern-Price法、文獻(xiàn)［16」的搜索算法以及本文方法得到的ACADS算例滑動(dòng)面位置?？梢钥闯?，不同方法所得到的滑動(dòng)面形狀相似，但其位置和形狀并不相同，并且與采用搜索算法的變分極限平衡法得到的結(jié)果基本吻合，但Morgenstern-Price方法由于利用了傳統(tǒng)的極限平衡法中對(duì)條間力的假設(shè)等原因?qū)е禄瑒?dòng)面范圍更大。

圖5（a）和（b）分別繪出了本文方法得到的法向和切向條間力分布曲線，并與Morgenstern-Price法得到的結(jié)果進(jìn)行了比較?？梢钥闯?，無(wú)論是法向還是切向，條間力峰值基本一致，但其出現(xiàn)的位置不同，Morgenstern-Price法得到的結(jié)果均要靠近坡頂，并且兩種方法得到的切向條間力分布存在較大的差別，這是由于基于傳統(tǒng)條分法的Morgenstern-Price法雖然考慮了土條間法向力和條間剪切力，且力系平衡條件滿足靜力平衡條件，但同時(shí)人為引入了王條間的剪切力和法向力之間假定的關(guān)系式，即條間力函數(shù)，使得計(jì)算結(jié)果出現(xiàn)偏差。由圖5還可以看出，坡頂張力的出現(xiàn)導(dǎo)致條間力在坡頂處出現(xiàn)轉(zhuǎn)折，這與龐峰等[20]的研究結(jié)果是一致的。

為驗(yàn)證本文提出的方法在分析地震作用下邊坡穩(wěn)定性的有效性，參照文獻(xiàn)22]，水平地震作用系數(shù)分別取為 k_h=0.1 和 k_h=0.2 ，進(jìn)行土坡算例驗(yàn)證分析。

首先考慮地震作用系數(shù) k_h=0.1 時(shí)土坡穩(wěn)定性情況。從圖6中能夠得到最小安全系數(shù)為1.619，與算例2基本一致，并得到坡腳處的法向應(yīng)力值為12.94kPa ，在坡頂處的法向應(yīng)力值為一 3.01kPa 。正如方程（4）所預(yù)期的那樣，最小安全系數(shù)值是對(duì)滑動(dòng)面函數(shù)和應(yīng)力函數(shù)極值問(wèn)題的求解。從圖中可以看到在坡頂處出現(xiàn)的拉應(yīng)力與陳建功等[16]的結(jié)果一樣，這種一致性說(shuō)明邊坡發(fā)生滑動(dòng)時(shí)在坡頂將會(huì)產(chǎn)生張拉裂縫。相較于應(yīng)力不參與變分分析的優(yōu)化方式[22]，本文的應(yīng)力和滑動(dòng)面是兩個(gè)獨(dú)立的函數(shù)，最后求解的最小安全系數(shù)值通過(guò)各狀態(tài)方程聯(lián)立得到，這使得最終結(jié)果不僅有滑動(dòng)面及其應(yīng)力分布對(duì)應(yīng)的最小安全系數(shù)，還包含滑動(dòng)土體的條間力分布情況，具體的函數(shù)曲線如圖7所示。

考慮水平地震系數(shù) k_h=0.2 時(shí)得到的最小安全系數(shù)為1.391，這與文獻(xiàn)[21]中的算例結(jié)果是一樣的。圖8繪出了滑動(dòng)面及其法向應(yīng)力分布圖，可以看出坡腳的法向應(yīng)力為 54.67kPa ，而在坡頂?shù)姆ㄏ驊?yīng)力值為 -11.04kPa 。通過(guò)與圖6比較可以發(fā)現(xiàn)，水平地震系數(shù)的變化顯著影響了邊坡潛在滑動(dòng)面的位置及其法向應(yīng)力的大小和分布。圖9同樣是該邊坡臨界滑動(dòng)面的條間力分布情況，能夠看到在坡頂處出現(xiàn)的轉(zhuǎn)折更為顯著。

為了分析地震作用對(duì)邊坡穩(wěn)定性的影響，選取前面的ACADS算例中的邊坡模型，取水平地震系數(shù)分別為 k_h=0，0.10，0.2，0.3 圖 10～12 分別繪出了不同水平地震系數(shù)下邊坡潛在滑動(dòng)面及其條間力分布情況。

從圖10可以看出，隨著地震作用的增大，邊坡潛在滑動(dòng)面在坡頂?shù)奈恢醚?x 軸正向移動(dòng)，也就是滑動(dòng)體的范圍在不斷增加，穩(wěn)定安全系數(shù)逐漸變小，邊坡穩(wěn)定性越差。從圖11、12可以看出，法向和切向條間力的變化趨勢(shì)均呈現(xiàn)從零增加到最大值后逐漸減小再增長(zhǎng)為零的一個(gè)過(guò)程。隨著地震作用的增加，條間力的也隨之增大。

3結(jié)論

結(jié)合擬靜力法，建立了土坡在水平地震作用下穩(wěn)定性問(wèn)題的變分極限平衡方程，結(jié)合數(shù)值算例，采用打靶對(duì)問(wèn)題進(jìn)行了求解并驗(yàn)證了模型和方法的有效性。

（1）變分極限平衡法較傳統(tǒng)的極限平衡法在數(shù)學(xué)上更為嚴(yán)謹(jǐn)，分析結(jié)果更加可靠。針對(duì)地震作用下二維邊坡的穩(wěn)定性問(wèn)題，基于打靶法能夠準(zhǔn)確快速地獲得穩(wěn)定系數(shù)、潛在滑動(dòng)面、滑動(dòng)面上法向應(yīng)力以及條件作用力的分布。

（2）隨著水平地震作用系數(shù)的增大，滑動(dòng)體的范圍將逐漸增大，邊坡的穩(wěn)定性逐漸降低，滑動(dòng)體的條間力幅值也隨之增大。

參考文獻(xiàn)：

[1］陳昌富，楊宇，龔曉南.基于遺傳算法地震荷載作用 下邊坡穩(wěn)定性分析水平條分法[J].巖石力學(xué)與工程學(xué) 報(bào)，2003，22（11）：1919-1923. CHEN Changfu，YANG Yu，GONG Xiaonan. Horizontal slice method of slope stability analysis under seismic load based on genetic algorithm[J]. Chinese Journal of Rock Mechanics and Engineering，20o3，22（11）： 1919-1923.

[2]LI ZY，XIAO SG.Seismic overallstability of embankment slopes retained by multi-step cantilever retaining walls using pseudo-static method[J]. Geotechnical and Geological Engineering，2023，41（3）：1893-1907.

[3］張飛，林立曜，舒爽，等.轉(zhuǎn)彎凸型邊坡靜動(dòng)力穩(wěn)定性 分析及空間效應(yīng)研究[J].巖土工程學(xué)報(bào)，2022，44 （8）：1558-1566. ZHANG Fei，LIN Liyao，SHU Shuang，et al.3D effectson static and dynamic stability ofconvex curved slopes[J].Chinese Journal of Geotechnical Engineering，2022，44（8）：1558-1566.

[4]TANTZ，HUANG MS，SHI ZH.A perturbation method for upper bound analysis of stability of slopes based on rigid translational/rotational moving elements [J].Computers and Geotechnics，2023，161：105595.

[5] 董建華，朱彥鵬，馬巍.錨固邊坡地震動(dòng)力穩(wěn)定性計(jì) 算方法[J].振動(dòng)工程學(xué)報(bào)，2013，26（4）：633-640. DONG Jianhua，ZHU Yanpeng，MA Wei. Dynamic stability calculation method for anchor slope based on stressstate[J]. Journal of Vibration Engineering， 2013，26（4）：633-640.

[6]KOP'ASCY J. Three-dimensional stress distribution and slip surfaces in earth works at rupture[C]//Proceedings of the 4th International Conference on Soil Mechanics and Foundations Engineering.London，UK：1957： 339-342.

[7] BAKER R，GARBER M. Theoretical analysis of the stability of slopes[J].Géotechnique，1978，28（4）： 395-411.

[8]LI X G，LIU W N. Study on the action of the active earth pressure by variational limit equilibrium method [J].International Journal for Numerical and Analytical Methods inGeomechanics，2010，34（10）：991-1008.

[9]周鳳璽，朱順望，梁玉旺，等.變分極限平衡法對(duì)土質(zhì) 邊坡穩(wěn)定性的精確分析[J].巖土工程學(xué)報(bào)，2023，45 （7）： 1341-1346. ZHOU Fengxi， ZHU Shunwang，LIANG Yuwang，et al.Exact analysis of soil slope stability by using variational limit equilibrium method[J].Chinese Journal of Geotechnical Engineering，2023，45（7）：1341-1346.

[10] LING H I， LESHCHINSKY D. Seismic performance of simple slopes[J]. Soils and Foundations，1995， 35（2）：85-94.

[11]SARKAR S，CHAKRABORTY M. Pseudo-static slope stability analysis for cohesive-frictional soil by usingvariational method[M]//Advances in Computer Methods and Geomechanics. Singapore： Springer Singapore，2020：159-171. stability of slopes：design chartslJ]. Journal of Geotechnical Engineering，1994，120（9）：1514-1532.

[13]BAKERR，SHUKHAR，OPERSTEINV，et al. Stability charts for pseudo-static slope stability analysis[J]. Soil Dynamics and Earthquake Engineering，2O06，26 （9）：813-823.

[14］欒茂田，金崇磐，林皋，等.層狀非均質(zhì)土坡抗震穩(wěn)定 性的變分解法[J].地震工程與工程振動(dòng)，1993，13 （4）：73-80. LUAN Maotian， JIN Chongpan， LIN Gao，et al. Variational procedure for aseismic stability of layered soil slopes[J]. Earthquake Engineering and Engineering Vibration，1993，13（4）：73-80.

[15］李亮，鄧東平，趙煉恒.任意滑動(dòng)面搜索新方法應(yīng)用 于三維復(fù)雜土坡穩(wěn)定性分析[J].巖土工程學(xué)報(bào)， 2011，33（4）：544-553. LI Liang，DENG Dongping，ZHAO Lianheng.Application of new method of arbitrary sliding surface search to stability analysis of complex three-dimensional slopes [J].Chinese Journal of Geotechnical Engineering， 2011，33（4）：544-553.

[16］陳建功，李會(huì)，賀自勇.基于變分法的均質(zhì)土坡穩(wěn)定 性分析[J].巖土力學(xué)，2019，40（8）：2931-2937. CHEN Jiangong， LI Hui， HE Ziyong.Homogeneous soil slope stability analysis based on variational method[J]. Rock and Soil Mechanics，2019，40（8）：2931-2937.

[17］陳祖煜.土質(zhì)邊坡穩(wěn)定分析：原理·方法·程序[M].北 京：中國(guó)水利水電出版社，2003. CHEN Zuyu. Soil Slope Stability Analysis：Theory， Methods and Programs[M]. Beijing：China Water amp;. Power Press，2003.

[18]ZHOU FX，MA Q.Exact solution for capillary interactions between two particles with fixed liquid volume[J]. Applied Mathematics and Mechanics，2016，37（12）： 1597-1606.

[19]ZHOUFX，LISR，LAIYM.Three-dimensional analysis for transient coupled thermoelastic response of afunctionally graded rectangular plate[J].Journal of Sound and Vibration，2011，330（16）：3990-4001.

[20］龐峰，李亮.摩根斯坦普萊斯方法條間力函數(shù)影響分 析[J].濟(jì)南大學(xué)學(xué)報(bào)（自然科學(xué)版），2012，26（4）： 412-414. PANG Feng，LI Liang. Discussion on inter-slice func tion inMorgenstern-Pricemethod[J].Journal of University of Jinan（Science and Technology），20l2，26（4）： 412-414.

[21]QIN C B，CHIAN S C.New perspective on seismic slope stability analysis[J]. International Journal of Geomechanics，2018，18（7）：06018013.

[22] HOU C， ZHANG R，LI Y，et al. Comparisons of safe ty factors for slope in nonlinear soils[J].KSCE Journal of Civil Engineering，2021，25（10）：3737-3749.