基于極限分析上限法雙層土坡穩定性分析

王 迪， 王宏權， 王曉飛， 徐惠民， 張 飛

(1.巖土力學與堤壩工程教育部重點實驗室， 江蘇 南京 210098；2.河海大學 土木與交通學院， 江蘇 南京 210098； 3.廣東省南粵交通大豐華高速公路管理中心， 廣東 梅州 514300； 4.安徽省交通規劃設計研究總院股份有限公司， 安徽 合肥 230088； 5.南京市水利規劃設計院股份有限公司， 江蘇 南京 210022)

邊坡穩定問題是土力學三大經典問題之一，是巖土力學與工程研究中的重點領域。關于邊坡穩定性問題，國內外學者提出了各種分析方法(極限平衡法、極限分析法、數值模擬法)，并被廣泛應用于實際工程。邊坡穩定性分析最早是基于Coulomb的擋土墻土壓力計算方法和Rankine的主動土壓力與被動土壓力計算方法，建立了最初的極限平衡法。隨后，條分理論的提出對邊坡穩定性極限平衡分析法的發展有重要貢獻，基于土體條分的各種假設建立了諸多方法[1-7]。由于土體塑性理論的建立，極限定理的思想被用于邊坡穩定性分析中，提出了極限分析法。伴隨著計算機技術和數值模擬方法的快速發展，各種數值分析方法被提出并引入到邊坡穩定性分析中。但極限平衡法[8-9]由于將土條簡化為剛體，忽略了巖土體的應力-應變關系，沒有考慮土體發生破壞時的變形協調條件，計算得到的安全系數僅是假定滑動面上的平均安全度，結果并不十分理想。數值模擬法雖然具有簡易的操作性和很強的可視化效果，但在本構模型、屈服準則以及失穩判據的選取方面缺少統一性，不同的模擬軟件本身也對計算精度有著不同程度的影響。極限分析法的優勢在于，它可以不必一步步地對各個階段的彈塑性狀態進行分析而直接求得破壞荷載。從理論上，極限分析法相比極限平衡法更為嚴格，極限分析法分為上限解法和下限解法，上限解需要建立運動許可的速度場(即破壞面)，下限解需要建立靜力許可的應力場。若上限解與下限解相等時，則就認為獲得了“真解”。

目前，邊坡穩定極限分析法的研究主要針對均質土坡，但實際工程中的邊坡問題常是非均質的，關于非均質邊坡的穩定性研究仍需要進一步的改進與完善。Chen等[10]采用對數螺旋線的轉動破壞機制，建立能量平衡方程得到均質邊坡穩定性的上限解，隨后該方法又被Chen等[11]拓展到非均質各向異性的邊坡分析中。Chen[12]對分層位置在坡腳處的雙層邊坡進行分析，給出了限制深度下的穩定系數。然而，Chen[12]研究中只考慮了分層位置在坡腳處的情況，自然界中成層邊坡的分層位置不一定都在坡腳處，對于任意分層位置雙層土坡，目前尚未開展相關研究。此外，邊坡失穩破壞時其破壞面的分布，一直是工程界比較關心的問題，但是針對不同參數變化對邊坡滑動面分布規律的研究相對較少，尤其是對于常見的雙層邊坡。

針對以上研究的不足和需要，本文基于Chen[12]提出的雙層邊坡穩定性的極限分析上限解法，推導了任意分層位置的雙層土坡在不同破壞機制下臨界高度的計算公式，依托MATLAB軟件平臺，通過編制不同破壞機制下的穩定系數的計算程序，分析了臨界高度與分層深度的關系，并得到不同條件下雙層土坡臨界高度隨分層深度改變的線性變化規律。同時探究了雙層土坡臨界滑動面位置的分布規律，建立雙層土坡的穩定性分析方法，為工程設計提供一定的參考依據。

1 任意分層位置雙層土坡穩定性極限分析方法

極限分析法[13-14]是基于理想彈塑性體的基本定理-上限定理和下限定理求解物體極限應力狀態的分析方法。對于巖土工程問題，通過一些簡單的假定和邊界條件的簡化[15-16]，極限分析法提供了一種尋求合理真解范圍的有效途徑。

1.1 基于坡底破壞機制雙層土坡臨界高度計算公式

圖1為雙層土坡坡底破壞機制示意圖，假設雙層土坡破壞機制為旋轉機構，破壞面為對數螺旋面，通過建立對數螺旋面上的能量耗損率與旋轉土體重力所做外功率的等量關系，得到雙層土坡極限分析上限解，進而用于雙層土坡的穩定性評價。

圖1 坡底破壞機制示意圖

假設滑動破壞區A-B-C-D-A繞旋轉中心O做剛體旋轉，而對數螺旋面A-E′-D以下土體保持靜止，因此，A-E′-D面是一個薄層的速度間斷面。設邊坡的高度為H，雙層土分界面E-E′以上土體高度為αH，基準線OA、OD和OE′的傾角分別為θ0、θh和θm，基準線長度分別為r0、rh和rm，邊坡的坡角為β，∠BDC為β′，AB的長度為L。考慮土體重度變化不大，假定雙層土坡上下兩層土體的重度均為γ；為便于構建統一滑動面，上下兩層土體的內摩擦角都為φ；上下兩層土體的黏聚力不同，分別為c1和c2。

對數螺旋面的方程為：

r(θ)=r0e(θ-θ0)tanφ

(1)

因此，基準線OD的長度為：

rh=r0e(θh-θ0)tanφ

(2)

在所考慮的滑動破壞區域A-B-C-D-A內，土體重力所做的外功率W，詳見Chen在《Limit Analysis and Soil Plasticity》[12]一書，在此僅給出其表達式。

(3)

內部能量耗損主要發生在速度間斷面A-E′-D上。能量耗散率的微分，可以由該面的微分面積rdθ/cosφ與滑動面上的黏聚力ci以及該速度間斷面上的切向間斷速度Vcosφ的乘積計算求得，因此，雙層土坡的總的內部能量耗散率可以沿著整個速度間斷面積分A-E′-D求得，即：

(4)

令外功率W與內能耗散率D相等，即式(3)與式(4)相等，可以得到：

(5)

化簡式(5)，即可得到雙層土坡的無量綱的臨界高度表達式：

f(θ0，θh，β′，H/r0,L/r0)

(6)

其中：f為目標函數；θ0，θh，β′，H/r0,L/r0為待搜索優化的變量。

1.2 基于坡面破壞機制雙層土坡穩定性極限分析方法

當滑動面通過雙層土坡坡腳以上時，土坡將發生坡面破壞，見圖2。此時，通過引入高度變量h(h

圖2 雙層土坡坡面破壞機制

這里，引入變量n，并假設：

(7)

此時，由幾何關系得：

(8)

A-B-D滑動破壞區域土重所做的外功率為：

(9)

內部能量耗損主要發生在速度間斷面A-E′-D上。能量耗散率的微分，可以由該面的微分面積rdθ/cosφ與滑動面上的黏聚力ci以及該速度間斷面上的切向間斷速度Vcosφ的乘積計算求得，因此，雙層土坡的總的內部能量耗散率可以沿著整個速度間斷面積分A-E′-D求得，即：

(10)

令外功率W與內能耗散率D相等，即式(9)與式(10)相等，則有：

化簡這個等式，則有：

(11)

其中：f為目標函數；θ0，θh，β′，h/r0，L/r0為待搜索優化的變量。

1.3 雙層土坡上限解優化計算過程

通過上面基于雙層土坡的坡底破壞和坡面破壞機制建立的能量平衡方程，得到雙層土坡的無量綱的臨界高度γH/c1的表達式：

(12)

其中：f是目標函數。公式(12)給出了雙層土坡臨界高度γH/c1的一個上限。為了得到γH/c1的上限解，需要對目標函數f進行優化計算，求其最小值，從而獲得γH/c1的最小上限解。

目標函數f有五個待優化的參數，為了提高計算效率，本文將采用陳祖煜[17]院士提出的隨機搜索方法。對于給定的邊坡參數值(坡高H，坡角β，土體黏聚力比值c2/c1，以及上層土的厚度αH)，通過對優化參數的初始值，初始搜索步長以及精度進行設定，進行迭代計算求解，從而得到目標函數f的最小值。優化過程中，角度(θ0，θh，β′ )的精度為0.01°，比值(H/r0，L/r0，h/r0)的精度為0.001。

2 雙層土坡臨界高度與分層深度關系

依托MATLAB數學計算平臺，針對不同破壞機制編制相應的邊坡穩定程序進行優化計算，綜合考慮各個影響因素的相互影響，通過整理計算結果，繪制了不同因素下雙層土坡的臨界高度與分層深度的關系圖(見圖3)。

圖3 雙層土坡臨界高度圖(φ = 20°)

從圖3可以看出，對于均質土坡(α=1)，當已知邊坡的幾何參數和土性參數時，其臨界高度為一定值，這與前人研究結果一致。對于雙層土坡，當邊坡的坡度β、坡內土體的內摩擦角φ以及上下兩層土的黏聚力比值c2/c1一定時，雙層土坡的臨界高度γH/c1隨分層位置α的變化呈現不同的規律性。

3 不同因素對雙層土坡最危險滑動面的影響

土坡發生失穩破壞時存在三種不同的破壞機制，不同破壞機制下其滑動面的位置具有明顯的差異性，而土坡的破壞機制取決于其幾何參數和坡內土體的力學參數。因此，針對雙層土坡在不同的幾何參數和力學參數影響下，研究其臨界滑動面位置的分布規律，旨在為工程實踐提供一定的參考。圖4和圖5為內摩擦角φ=5°時雙層土坡在不同的坡角β、上下土層不同分界面位置α及其黏聚力比值c2/c1下，其臨界滑動面的位置分布圖。

圖4 α = 0.50時雙層土坡滑動面位置分布圖(φ= 5°)

從圖4中可以看出，隨著黏聚力比值c2/c1的增大，雙層土坡臨界滑動面的位置越來越接近于坡面，且對于極緩坡度(β不大于15°)的邊坡，主要發生坡底破壞。隨著坡角β和黏聚力比值c2/c1的增大(即下層土較上層土堅硬程度的增加)，邊坡由坡底破壞逐漸發展為坡腳破壞，這在圖4(b)中表現的尤為明顯。隨著坡角β的繼續增加，邊坡主要發生坡腳破壞。然而當坡角增大到一定程度(坡角β不小于60°)時，隨著黏聚力比值c2/c1的增加，邊坡將在兩層土的分界面處發生坡面破壞，由圖4(c)中可以清晰的看出，當c2/c1≥ 5時，這些邊坡發生了坡面破壞。

圖5為α=1.0時雙層土坡滑動面位置分布圖。坡底破壞主要出現在較緩邊坡，較陡邊坡主要發生坡腳破壞，且隨著下層土堅硬程度的增加，雙層土坡的破壞機制發生改變；對于上硬下軟型土坡，隨著坡角β的增加，對于較陡邊坡也可發生坡底破壞。

圖5 α=1.0時雙層土坡滑動面位置分布圖(φ= 5°)

對比圖4可以看出，當分層深度α發生變化時，相同坡角β和黏聚力比值c2/c1下的雙層土坡的破壞形式發生了改變，這在圖4(b)和圖5(b)中表現的較為明顯，同時也說明了分層深度α的變化改變了雙層土坡的臨界失穩機制。對于該邊坡，隨著分層深度α的增加，土坡由坡腳破壞轉變為坡底破壞。

同理，通過研究內摩擦角φ=10° 時雙層土坡在不同的坡角β、上下土層不同分界面位置α及其黏聚力比值c2/c1下，其臨界滑動面的位置分布，當內摩擦角φ增大到10°時，對于很緩的邊坡(如β<15°)，其破壞形式基本上均為坡腳破壞，且隨著黏聚力比值的增加，其滑動面的位置基本保持不變，但對于垂直邊坡，當黏聚力比值很大(c2/c1≥5)時，雙層土坡將發生坡面破壞。

隨著分層深度的變化，當α為0.5時，當黏聚力比值很大(c2/c1≥5)時，除坡度很緩的邊坡(β<15°)，雙層土坡均可發生坡面破壞，其他規律與之前具有很好的一致性。

當分層深度在坡腳處(α=1.0)時，與之前規律不同的是，對于很陡(β≥75°)的邊坡，無論是上軟下硬型雙層土坡還是上硬下軟型雙層土坡，均發生坡腳破壞，其滑動面位置基本上不隨黏聚力比值的變化而發生改變。對于其他坡度的上硬下軟型雙層土坡，則有會發生坡底破壞。

4 結 論

本文推導了雙層土坡不同破壞機制下臨界高度的計算公式，編制了相關程序進行雙層土坡穩定性計算；在不同幾何參數和力學參數條件下，對雙層土坡滑動面的分布規律進行了探究，得出以下結論：

(1) 對于上硬下軟型土坡(c2/c1<1)，雙層土坡的臨界高度γH/c1隨分層深度α的增加而線性增加，也就是說臨界高度γH/c1隨著上層較硬土層厚度的增加而線性增加。對于上軟下硬型土坡(c2/c1>1)，雙層土坡的臨界高度γH/c1隨分層深度α的增加而減小，即臨界高度γH/c1隨著上層較軟土層厚度的增加而減小。

(2) 坡底破壞主要出現于內摩擦角較小(φ<10°)的較緩(β<45°)邊坡，但對于上硬下軟型雙層土坡，在內摩擦角較大且分層深度在坡腳位置處時，也可發生坡底破壞；坡面破壞主要發生在分層深度在坡面中部及中部以下位置處的上軟下硬型土坡；坡腳破壞在不同分層深度和不同力學參數下的雙層土坡中均有出現。

(3) 雙層土坡臨界滑動面隨著黏聚力比值c2/c1的增加越接近與坡面，隨著分層深度的增加和坡角β的增加越來越容易出現坡面破壞，且臨界滑動面隨內摩擦角φ和坡角β的增加，不同黏聚力比之下的滑動面的位置愈加接近。