2 上限分析
在極限狀態(tài)下，假設均質邊坡的潛在滑動面y（x）［極坐標系下為r（θ）］通過坡腳，坡面函數(shù)為y1（x），α和β為邊坡傾角.邊坡坡頂破壞點（即滑入點）至坡腳的垂直高度為H，水平距離為L；坡頂點與坡腳的垂直高度為H1，水平距離為L2；坡頂點與滑入點的距離為L1，建立如圖2所示的均質邊坡旋轉破壞機制.滑動面以上巖土體可視為剛性體，以角速度沿滑動面繞點O轉動，且邊坡失穩(wěn)時瞬時變形忽略不計.假定均質邊坡巖土體材料服從三參數(shù)非線性破壞準則和相關聯(lián)流動法則，根據(jù)極限分析上限定理，塑性速度矢量δw與滑動面的夾角為φt.
圖2 均質邊坡旋轉破壞機制Fig.2 Rotational failure mechanism of homogeneous slope
2.1 能耗分析
非線性破壞準則下，滑動面上單位面積的內能耗散功率為［8，27］：
將式（3）代入式（2）并積分可得臨界狀態(tài)下滑動面上總的內能耗散功率
式中：r表示滑動面上一點的旋轉半徑；θ為極徑與水平方向的夾角，順時針為正.
在本文中，外力包含重力和地震力，此處采用水平和垂直地震加速度系數(shù)kx、ky表征地震對邊坡穩(wěn)定性的影響，由于都是體積力，因此可一并計算.值得注意的是，當kx=0、ky=1.0 時與自然條件下相同，即只有重力作用.基于圖2 的邊坡滑動破壞機制，在體積力作用下，外力功率可計算如下
式中：We為滑動體外力功率；W為坐標軸x、x=xn與滑動面所圍區(qū)域的外力功率；WⅠ為區(qū)域Ⅰ的外力功率；WⅡ為區(qū)域Ⅱ的外力功率；WⅢ為區(qū)域Ⅲ的外力功率.
式中：γ為巖土體材料容重，單位kN/m3.
根據(jù)虛功率原理，有
將式（6）代入式（7）中并整理可得
在極限狀態(tài)下，滑動體的力學平衡條件和材料屈服條件必須同時滿足，且認為當能量平衡和力學平衡條件同時滿足時可以獲得最小的上限解答［15］.因此，在地震擬靜力荷載和自重作用下，滑動體的力學平衡方程如下
式中：Q、V分別為滑動體水平和豎直方向的合力；M為滑動體對坡頂滑入點的矩.
因此，滑動體的總虛功可表示為
式中：δu、δv和δΩ分別為水平、豎直和轉動的虛位移.
將式（10）代入式（11）可得總虛功為：
當I的一階變分為零時，邊坡處于臨界狀態(tài).此時，泛函I滿足如下歐拉方程：
式中：G表示臨界狀態(tài)下邊坡潛在滑動面上正應力和切應力滿足的相關關系，在本文中即為三參數(shù)破壞準則.
對于Baker提出的三參數(shù)非線性屈服準則，有：
于是將式（15）（22）代入式（21）可得：
由式（23）可知，線性準則下，即φt為常量時，均質邊坡潛在滑移面為對數(shù)螺旋線；反之，在非線性準則下，邊坡潛在滑移面并非單一對數(shù)螺旋線.
同理，將式（16）（18）代入式（17c）中可得：
通過變分原理的歐拉方程，得到了滑動面及其正應力分布的一階微分方程組，求解該常微分方程組，還需要邊界條件.由變分極值原理可知，當極值曲線的端點不固定時，泛函取極值的必要條件為滿足變分橫截條件［28］.由于坡頂破壞點在坡頂表面上移動，因此，泛函取極值的變分橫截條件為：
需要注意的是，在滑動面端點處，有：
故含有kx和ky的兩項消失，再結合式（15）和（22）可得：
由于滑動體的瞬時變形忽略不計，可采用變形前的尺寸進行分析.因此，根據(jù)圖2 的幾何關系，有以下等式成立：
2.2 求解步驟
根據(jù)2.1節(jié)能耗分析，可以獲得邊坡給定參數(shù)條件下的臨界高度，求解流程如圖3 所示.詳細計算步驟如下：
圖3 邊坡臨界高度求解流程圖Fig.3 Flow chart for solving critical height of slope
1）已知邊坡坡頂傾角及坡面傾角、巖土體強度參數(shù)和地震系數(shù)；
2）假定一組輸入?yún)?shù)（r0，θ0），根據(jù)變分橫截條件［（式（30）］求解坡頂破壞點處的正應力σ0；
3）選擇迭代步長Δθ，使得θi+1=θi+Δθ，根據(jù)Runge-Kutta 法在區(qū)間［θi，θi+1］上求解式（23）和式（26），得到滑動面上點（i+1）的坐標（ri+1，θi+1）和應力值σi+1、τi+1；
4）計算點（i+1）處累計的內能耗散功率和外力功率，再根據(jù)虛功率平衡原理，求解式（8）得到坡頂破壞長度L1；
5）若L1大于零，則根據(jù)式（31）計算相應的邊坡坡度和高度，反之則返回第3步；
6）判斷計算邊坡坡度是否和已知邊坡坡度相等，若不相等，則返回第2步；
7）輸出（r0，θ0，H），得到均質邊坡臨界高度最優(yōu)上限解及其對應的潛在滑動面.
其中，通過第2 步～第6 步可以得到已知坡度下邊坡的臨界高度，進而利用二維免疫算法得到最優(yōu)上限解.需要說明的是，由于本文構建的邊坡潛在滑動面通過坡腳，因此上述推導公式只適用于坡腳破壞機制，對于坡面及坡腳以下這兩類破壞機制，該程序尚不能獲得其臨界高度.
2.3 免疫算法
免疫算法是一種受生物免疫系統(tǒng)的啟示而形成的智能搜索算法.此算法適用于求解連續(xù)函數(shù)的極值問題，對于非連續(xù)函數(shù)的極值問題，也具有很強的全局搜索能力.基本原理是將優(yōu)化問題中待優(yōu)化的問題、可行解、可行解質量與免疫系統(tǒng)中抗原、抗體、免疫細胞與抗原的親和度一一對應.因此可以將生物免疫應答中的進化過程轉換成數(shù)學上的尋優(yōu)過程，免疫算法采用群體搜索策略，通過迭代計算，最終以較大的概率得到問題的最優(yōu)解.主要步驟如下.
1）假定免疫個體維數(shù)為D1，免疫個體數(shù)目為N，則初始種群的位置表示為：
式中：rand（0，1）表示在［0，1］上服從均勻分布的隨機數(shù)；xij（0）表示初代種群中第j個粒子的第i個分量的位置.
2）計算個體的親和度、抗體濃度以及激勵度，并按激勵度大小進行排序.
式中：J表示激勵度；U為該個體的親和度值，即函數(shù)值；K為抗體濃度；p、q稱為激勵度系數(shù).
3）取激勵度前N/2 個個體進行免疫操作（克隆、變異和克隆抑制），免疫后的種群重新計算激勵度.
4）隨機生成N/2個個體的新種群，計算個體的親和度、抗體濃度以及激勵度；免疫種群與隨機種群合并，按激勵度大小排序，進行免疫迭代.
5）反復執(zhí)行步驟2）～4），直至達到最大免疫迭代次數(shù)或所要求的收斂精度結束.
3 對比分析
3.1 參數(shù)轉換
為了驗證本文非線性上限變分分析結果的準確性和有效性，將本文方法與極限分析有限元方法（OPTUM G2 2021）計算結果進行對比.對比分析中采用非線性Hoek?Brown 準則，由于OPTUM G2 中的Hoek-Brown 破壞準則是由地質參數(shù)D、σci、GSI 和mi表示，因此在對比驗證之前應進行參數(shù)的轉換.Hu等［29］提出了一種將Hoek?Brown 參數(shù)與非線性Mohr?Coulomb破壞準則強度參數(shù)的轉換計算方法.由于非線性Mohr?Coulomb 準則與本文的三參數(shù)破壞準則表達形式一致，故也可采用此方法得到Hoek?Brown準則參數(shù)與三參數(shù)破壞準則參數(shù)之間的轉換關系，詳細步驟如下.
Hoek 等［30］提出了一種方法實現(xiàn)兩種形式的Hoek?Brown 破壞準則之間的參數(shù)轉換.Hoek-Brown破壞準則的主應力形式為：
式中：σ1和σ3分別為最大和最小主應力；a、mb和s是與D、σci、GSI 和mi相關的無量綱參數(shù)，D為擾動因子，σci為單軸抗壓強度，GSI為地質強度參數(shù)，mi為材料常數(shù)［31］.
此外，以切向應力和法向應力表示的Hoek-Brown破壞準則為：
式中：τ和σ分別為切向和法向應力；A2和B2為反映材料特性的無量綱參數(shù)；σtm為單軸抗拉強度.
當給定參數(shù)D、σci、GSI和mi時，可以通過線性回歸分析得到參數(shù)A2、B2和.另外，對比式（1）和式（37）可知，切應力形式的Hoek-Brown 破壞準則與三參數(shù)破壞準則表達式一致，故可以獲得參數(shù)（A2、B2、σci和σtm）和（n、T和A）之間的轉換關系，具體關系如下：
3.2 算例驗證
算例邊坡的形狀尺寸和材料參數(shù)均取自孫超偉等［32］的研究，具體參數(shù)為：H=10 m，α=0°，β=45°，γ=25 kN/m3，kx=0，ky=1.0，彈性模量E=5 000 MPa，泊松比μ=0.3，σci=25 MPa，mi=2，GSI=5，D=0.
采用OPTUM G2 建立邊坡模型，土體材料服從Hoek-Brown 破壞準則和相關聯(lián)流動法則，模型設置標準邊界條件，即左右邊界設置法向約束、底部邊界設置固定約束.結合強度折減和網(wǎng)格自適應技術，采用三節(jié)點三角形單元，單元數(shù)量為3 000個.
由OPTUM G2 的有限元極限分析法（FELA）強度折減分析，可以獲得均質邊坡的安全系數(shù)、剪切耗散帶以及滑動帶上的主應力分布，一般情況下可認為該剪切耗散帶即為邊坡的潛在滑動面.此外，通過本文方法亦可獲得邊坡的安全系數(shù)、臨界狀態(tài)下的潛在滑動面及其正應力和切應力分布，由式（39）可將該正應力和切應力轉化為主應力，從而與極限分析上限有限元滑動帶上的主應力結果進行對比，以此來驗證本文方法的準確性和有效性，如圖4所示.
圖4 對比分析結果Fig.4 Comparative analysis results
圖4（a）為計算滑動面與極限分析上限有限元剪切耗散帶的對比.可以看出，兩者形狀基本吻合.本文方法計算的安全系數(shù)Fs=1.073，OPTUM G2 計算的安全系數(shù)Fs=1.125，相對誤差為4.6%.由于切應力形式的Hoek?Brown 準則與三參數(shù)破壞準則之間的參數(shù)轉換遵循一一對等原則，因此相對誤差的主要來源在于兩種形式的Hoek?Brown 破壞準則的參數(shù)轉換.
圖4（b）為兩種方法計算所得滑動面上σ1和σ3的對比.不難看出，兩者都是先增大后減小的變化趨勢；兩者出現(xiàn)極大值的位置基本一致，并且出現(xiàn)最大主應力的位置不在坡腳處.因此，以上對比可以證明本研究的準確性和有效性，能夠為邊坡的加固設計提供理論支持和合理參考.
4 參數(shù)分析
4.1 穩(wěn)定系數(shù)分析
根據(jù)前述計算方法，求解了給定邊坡的臨界高度.為了方便研究，采用了Baker［33］定義的穩(wěn)定系數(shù)Fn，計算公式如式（40）：
式中：Fs為邊坡安全系數(shù)，本文求解邊坡臨界高度，即Fs=1.
為探究巖土體非線性強度參數(shù)、地震荷載、邊坡傾角等因素對邊坡穩(wěn)定性的影響規(guī)律，采用上述極限分析上限定理，計算邊坡的臨界高度，再由式（40）計算邊坡的穩(wěn)定系數(shù)，以Fn為縱坐標繪制如圖5 所示的影響規(guī)律曲線.各圖例中非線性強度參數(shù)（T、A、n）、邊坡傾角（β）及水平地震加速度系數(shù)（kx）的取值如表1所示，此外ky=1.0，γ=25 kN/m3，α=0°.
表1 不同圖例下各參數(shù)取值Tab.1 Values of parameters in different legends
圖5 不同參數(shù)對邊坡穩(wěn)定系數(shù)影響規(guī)律分析Fig.5 The influence of different parameters on slope stability coefficient
圖5（a）研究了無量綱參數(shù)T從10～30 變化以及n從0.5～0.9 變化的條件下對穩(wěn)定系數(shù)Fn的影響.計算結果表明：當無量綱參數(shù)T相同時，隨著n的增大，邊坡穩(wěn)定系數(shù)逐漸減小，當n≥0.6 時，減小的趨勢出現(xiàn)加快，且當T越大時，減小的趨勢越顯著；當n相同時，隨著無量綱參數(shù)T的增大，邊坡穩(wěn)定系數(shù)逐漸增大，且當n越大時，增大的幅值越小.
圖5（b）研究了無量綱參數(shù)A從1.0～3.0 變化以及n從0.5～0.9 變化的條件下對穩(wěn)定系數(shù)Fn的影響.計算結果表明：當無量綱參數(shù)A相同時，穩(wěn)定系數(shù)隨n的增大呈減小的趨勢，且這種趨勢越來越緩慢；當n相同時，無量綱參數(shù)A對穩(wěn)定系數(shù)幾乎沒有影響.
圖5（c）研究了邊坡傾角β從40°～80°變化以及n從0.5～0.9 變化的條件下對穩(wěn)定系數(shù)Fn的影響.計算結果表明：當邊坡傾角β相同時，隨著n的增大，邊坡穩(wěn)定系數(shù)逐漸減小，當n≥0.8 時，減小的趨勢出現(xiàn)明顯地減緩，且當β越大時，減緩趨勢越顯著；當n相同時，穩(wěn)定系數(shù)隨邊坡傾角β的增大而減小，且減小的趨勢越來越緩慢.
圖5（d）研究了地震系數(shù)kx從0.10～0.30 變化以及n從0.5～0.9 變化的條件下對穩(wěn)定系數(shù)Fn的影響.計算結果表明：當?shù)卣鹣禂?shù)kx相同時，穩(wěn)定系數(shù)隨n的增大而減小，當n≥0.6 時，減小的趨勢出現(xiàn)明顯地加快，且當kx越小時，減小的趨勢越顯著；當n相同時，隨著kx的增大，邊坡穩(wěn)定系數(shù)逐漸減小，且當n越大時，減小的幅值越小.
4.2 滑動面及其應力分析
計算邊坡穩(wěn)定系數(shù)是為了給邊坡穩(wěn)定性評估提供量化指標，但是當邊坡穩(wěn)定性不足時，則需要對其采取加固防護措施，而加固措施的選擇取決于最危險滑動面位置.因此，確定邊坡潛在滑動面位置顯得尤為重要.本節(jié)研究了非線性強度參數(shù)、邊坡傾角、地震系數(shù)等因素對潛在滑動面及其應力分布的影響規(guī)律.
由圖6（a）可知，L1（坡頂點距坡頂破壞點的距離）隨著無量綱參數(shù)T的增大而增大，H（坡頂破壞點至坡腳的垂直高度）隨著T的增大同樣增大；由圖6（b）可知，坡頂破壞點（y=0 m 處）的正應力隨著T的增大而減小，滑動面上的最大正應力隨著T的增大而增大，且出現(xiàn)最大正應力的位置不在坡腳處.
圖6 不同T時滑動面及其應力分布Fig.6 Sliding surface and its stress distribution with different T
由圖7（a）可知，隨著無量綱參數(shù)A的增大，L1逐漸增大，并且H隨著A的增大同樣增大；由圖7（b）可知，坡頂破壞點的正應力為0，這是由于當T=0 時最小正應力為0，不存在負值，且根據(jù)邊界條件式（30）計算也可以得到此結果，同時由圖可知滑動面上的最大正應力隨著A的增大而增大，且坡腳處出現(xiàn)正應力最大值.
圖7 不同A時滑動面及其應力分布Fig.7 Sliding surfaces and its stress distribution with different A
由圖8（a）可知，隨著邊坡傾角β的增大，L1逐漸減小，并且H隨著β的增大同樣減?。挥蓤D8（b）可知，滑動面上的最大正應力隨著β的增大而減小，當β≤60°時，出現(xiàn)最大正應力的位置不在坡腳處.
圖8 不同β時滑動面及其應力分布Fig.8 Sliding surface and its stress distribution with different β
由圖9（a）可知，L1隨著地震系數(shù)kx的增大而增大，但是H隨著kx的增大而減小；由圖9（b）可知，坡頂破壞點的正應力隨著kx的增大而減小，并且滑動面上的最大正應力隨著kx的增大同樣減小，出現(xiàn)最大正應力的位置不在坡腳處.
圖9 不同kx時滑動面及其應力分布Fig.9 Sliding surface and its stress distribution with different kx
5 結論
基于三參數(shù)非線性破壞準則和極限分析上限定理，建立均質邊坡旋轉失穩(wěn)機制，考慮地震荷載的作用，根據(jù)力學平衡方程結合變分原理推導了邊坡潛在滑動面及其應力分布函數(shù)的微分方程組，采用Runge-Kutta 方法進行求解，并根據(jù)虛功率原理結合免疫算法求解了邊坡的最小臨界高度.本文無需假定邊坡潛在滑動面形式，且避免采用“外切線技術”引入瞬時抗剪強度指標，保證了滑動面上抗剪強度參數(shù)與應力的非線性關系，可真實反映巖土體材料的非線性特征，并且適用于任何形式的非線性破壞準則，相較于目前嵌入非線性破壞準則尚不完善的極限分析有限元方法，本文方法在處理非線性破壞準則方面具有一定的優(yōu)越性.該方法進一步分析了非線性強度參數(shù)、地震荷載、邊坡傾角等因素對邊坡穩(wěn)定性的影響規(guī)律，主要結論如下：
1）非線性強度參數(shù)n和邊坡傾角β對邊坡穩(wěn)定系數(shù)的影響顯著，無量綱參數(shù)T和地震系數(shù)kx隨著n的增大對穩(wěn)定系數(shù)的影響越來越小，而無量綱參數(shù)A對穩(wěn)定系數(shù)幾乎沒有影響.
2）L1隨著無量綱參數(shù)T、A和地震系數(shù)kx的增大而增大，隨著邊坡傾角β的增大而減??；H隨著T、A的增大而增大，隨著β、kx的增大而減小.
3）潛在滑動面上最大正應力隨著無量綱參數(shù)T、A增大而增大，隨著邊坡傾角β和地震系數(shù)kx增大而減小，與各參數(shù)影響邊坡臨界高度H的結果相似，并且當T=0時，坡頂破壞點（y=0 m處）正應力為0.
