基于上限定理的含裂縫邊坡穩定性圖表分析方法

高明輝， 張雅賢， 王奇智， 趙建青， 孫 超， 周記名， 張娜娜

(1.河北省地質資源環境監測與保護重點實驗室， 石家莊 050021； 2.中央軍委后勤保障部安置住房保障中心， 北京 100036； 3.河北科技大學建筑工程學院， 石家莊 050018； 4.中國土木工程集團有限公司， 北京 100038)

大量道路邊坡現場勘測發現，邊坡表面分布著大量裂縫，而裂縫的存在對其整體安全性的影響顯著，因此有眾多學者對裂縫邊坡展開了深入的研究[1-5]。并對裂縫的擴裂范圍進行了限定，但常常忽略了裂縫位置存在的影響[6]，導致得到的臨界裂縫破壞狀態不是最不利裂縫誘發的。因此，近年來，不斷有學者對裂縫擴展的形態和出現的區域進行理論上的創新和方法的改進。Utili[7]使用Chen[8]提出的極限分析上限法的同時，考慮了裂縫可能存在的不同深度和位置；Michalowski[9]則將裂縫的擴展與邊坡的破壞機制相結合進行研究分析，而后以此研究為基礎，進行了深入的探討與分析。

導致裂縫邊坡失穩和滑坡的因素很多，典型的如降雨[10]、地震[11]、交通負荷等或這些因素的共同作用。在邊坡危險性評估分析時，穩定性系數Ns、安全系數Fs及動力荷載下的邊坡位移大小，均可作為評估其是否失穩的指標，而安全系數Fs一直是預測坡體失穩的常用方法。強度折減法是一種簡單、有效、系統的確定安全系數取值的技術，學者們在巖土工程設計計算中廣泛采用了這一技術[12-14]。當使用強度折減法評估邊坡穩定性，需要假想坡體兩個重要的參數：黏聚力c和內摩擦角φ，同時折減到相同的安全系數。

對于常見的均質邊坡設計，Baker等提出了基于極限平衡法的靜力和動力條件下邊坡穩定性分析的完整設計圖；Michalowski等[15-16]提供了一系列基于極限分析方法考慮孔隙水壓力和地震力的設計圖表；Steward等[17]提出了計算5種不同破壞類型下土質邊坡安全系數的兩種設計圖表。前人的研究表明，各類方法，如極限平衡法(limit equilibrium method)、有限元法(finite element method)或有限差分法(finite difference method)、極限分析法(limit analysis method)等，已被廣泛使用以準確得到邊坡的穩定性參數。以往研究中，較為常見的是傳統的極限平衡法[18]，特別是生成需要大量安全系數數據設計圖表時，與有限元方法相比，它的時間消耗更少。然而，傳統的極限平衡法并未考慮到土體的本構關系。而極限分析上限法考慮了能量的守恒，有效解決了土體應力應變關系的問題，并且給出了精確的臨界上限值，使得當前極限分析上限法成為分析動靜態邊坡穩定時經常采用的方法。

現基于經典的極限分析上限法，結合強度折減法對裂縫邊坡的安全系數進行綜合性分析，針對常見地震荷載條件下、裂縫中的靜水壓力以及坡頂部荷載作用等進行研究分析，將臨界裂縫的形態及其對應的深度、位置展示，并與有限元軟件Optum-G2的結果進行對比，以期為邊坡安全性初步快速評估提供科學依據。

1 基本破壞機制以及公式推導

R0為摩爾圓半徑；δ為剪脹角； f′c為單軸抗壓強度；σ為應力； τ為剪切力；速度矢量δw； (δu, δv)為速度分矢量圖1 邊坡破壞強度包絡線Fig.1 Slope failure strength envelope

理想的塑性邊坡往往表現為達到極限應力時開始發生屈服。經典Mohr-Coulomb(M-C)屈服準則往往用于評估邊坡安全性時首選的破壞準則，即屈服準則中采用線性的強度包絡線。采用M-C屈服方法時，抗剪強度由土體的黏聚力c和內摩擦角φ來預測，如圖1所示。而對土體在實測的壓縮狀態往拉伸方向外推得到其抗拉強度，所以并不是實際實驗而得到的抗拉強度。而針對裂縫邊坡，裂縫的存在可以有效地解決邊坡抗拉強度的問題，因此用經典的M-C屈服準則對分析是合理的。而使用強度折減法進行邊坡穩定性分析時候，需要同時將土體的黏聚力c和內摩擦角φ折減到一定的程度Fs,可表示為

(1)

式(1)中：cf、φf分別為折減后的黏聚力和內摩擦角。

極限分析上限法基于能量的守恒，破壞機制可以包含連續變形的區域，也可以由剛性塊構成：平移(或轉動)，由中間脆弱的剪切帶隔開，通常被視為運動場中的速度不連續性。當用線性屈服條件描述土體的強度并用相關聯流法則約束變形時，速度矢量δw與強度屈服面呈垂直狀態，可將其分解成兩部分：水平分量δu與豎直分量δv。通過應力矢量(σ,τ)和速度矢量(δu,δv)的乘積來表示非連續面上每單位面積的能量耗散率dc，可表示為

dc=cvcosφ

(2)

式(2)中：v為單位旋轉速度。

在M-C屈服條件下，將邊坡外部做功與內能能量耗散(滑裂面做功)相平衡，如式(3)所示，破壞機制如圖2所示。邊坡右側BE為豎直裂縫，即預設的邊坡后緣的產生拉應力區的裂縫位置，裂縫具有一定的位置x(與破址角的直線距離)及一定的深度y,而裂縫在不同深度、位置的限制下產生不同的安全系數，從而形成最小的安全系數，x、y可由式(5)、式(6)得到。裂縫形態如圖3所示，計算目標為可能出現最不利影響時裂縫所處的位置x和所達到的深度y。

β為坡腳；θ0、θc、θh為旋轉角；khG為水平加速度； kvG為豎直加速度；G為重力，q為上部荷載； r0、rh分別為旋轉點到A點和D點的距離；H為邊坡高度圖2 裂縫邊坡破壞機制Fig.2 Failure mechanism of slope with crack

L2為裂縫處到坡上緣的距離圖3 臨界裂縫位置及深度Fig.3 Critical crack position and depth

(3)

在實際破壞條件下，地震加速度并沒有固定傳遞方向。而計算地震作用引起的外部功率的推導方法與邊坡自身重量引起的外部功率的方法一致。在計算水平外力作用下的外部功的比率時，唯一的區別是速度方向。因此，將極限分析法與強度折減法相結合，可以得到裂縫邊坡在地震條件下的安全系數計算公式為

[(f1-f2-f3-t1+t2+t3)(1+λkhG)+

(4)

(5)

sinθ0}-1+cotβ

(6)

考慮裂縫中靜水壓力εw影響時，引入無量綱系數K表示裂縫中的含水量，K=0時表示裂縫中未出現靜水壓力，而K=1時表示裂縫中充滿了靜水壓力。則加入裂縫中的靜水壓力后的安全系數Fs、fr的計算公式分別為

(7)

(8)

式(8)中：γw為水的容重。

考慮邊坡上部荷載作用時，引入無量綱系數qt，qt=q/γH。整理后的邊坡安全系數Fs和fq的表達式分別為

[(f1-f2-f3-t1+t2+t3+qtfq)+

(9)

(10)

(11)

2 程序優化與結果分析

先前求解安全系數的計算方法存在顯著的不足，常常存在搜索精度不高、迭代計算復雜以及無法獲取全局最優解的問題。文章為解決這類問題，采用隨機搜索法與序列二次優化法(sequential quadratic programming,SQP)結合的計算方法，即采用MONTE-CARLO算法在一定精度范圍內確定一個初值，再使用SQP方法確定全局最優解，計算效率顯著提高，且精準度更高。具體計算邏輯如圖4所示。

圖4 優化算法邏輯圖Fig.4 Optimization algorithm logic diagram

考慮破壞工況為：坡腳30°～80°、內摩擦角10°～30°、參數組合c/γH為0.05～0.3、地震系數(0～0.2)g(g為重力加速度)、裂縫靜水壓力系數：0～0.3及上部荷載系數：0.1～0.15。為驗證計算的準確性，將結果與經典結果[19]進行對比，對比參數如表1所示，計算結果如表2所示，結果顯示，針對裂縫邊坡，可以明顯看到其相對于完整邊坡(無裂縫)呈現出更加不利的安全系數，因此對于裂縫邊坡的破壞穩定性研究，尤為關鍵。

表1 算例參數設計

表2 不同計算方法結果對比

2.1 考慮地震作用的影響

針對地震的影響，采用擬靜力法進行研究分析，即針對khG=(0～0.2)g進行研究分析。橫坐標表示邊坡坡腳β，縱坐標表示計算求得的安全系數Fs，考慮內摩擦角的變化，以邊坡參數組合c/γH為變量進行研究分析，研究地震力作用時，為僅考慮地震荷載影響的結果，因此并未加入裂縫中的靜水壓力以及坡頂上部荷載的影響。結果顯示，地震荷載對于邊坡安全系數的影響非常顯著，圖5(a)～圖5(c)展示了不同地震荷載下的裂縫邊坡的安全系數，結果顯示考慮地震荷載的作用大大折減了Fs，如針對c/γH=0.05、β=30°、φ=20°的邊坡，不考慮地震荷載時，其安Fs=1.24，坡體整體處于安全狀態，而考慮地震影響時候，如khG=0.2g時，裂縫邊坡的Fs=0.85，相比較靜態降低了31%，且此時坡體處于滑落狀態，因此地震作用力對于邊坡安全性的影響尤為顯著。而豎向地震動對于邊坡穩定性也存在一定程度的影響，圖5(d)展示了不同程度的豎向分析系數λ對于邊坡安全系數的影響，結果顯示，考慮豎向地震荷載的存在對于邊坡的影響顯著，考慮豎向地震系數λ為正時，即向下的豎向加速度，邊坡的安全系數隨之降低。在針對考慮地震力時的帶裂縫邊坡安全分析時，工程師可直接采用圖5中一系列圖表參數進行選取評估。

2.2 考慮邊坡上部荷載作用的影響

將上部荷載的作用加入裂縫邊坡的討論中，考慮上部荷載系數qt(qt=q/γH)為0.05～0.15的工況，計算結果整理到圖6中，與圖5坐標軸一致，區別了地震荷載與上部荷載的作用，結果表明：上部荷載作用下，邊坡的安全系數產生一定的折減，尤其針對c/γH較大時候的工況，而邊坡的Fs同樣隨著上部荷載qt的增加而降低。上部荷載的作用將會對邊坡的穩定性帶來較大程度的影響，如針對c/γH，β=80°、φ=10°的邊坡，當不考慮上部荷載作用時，根據圖6(a)可見，邊坡Fs=1.18，而考慮上部荷載的作用，當qt=0.15時，邊坡Fs=0.95，降低了19.5%，折減程度明顯，因此邊坡穩定性分析時，需將上部荷載的作用進行討論。

圖5 地震荷載作用下裂縫邊坡的安全系數Fig.5 Safety factor of slope with crack under earthquake load

圖6 上部荷載作用下裂縫邊坡的安全系數Fig.6 Safety factor of slope with crack under surcharge load

2.3 裂縫中靜水壓力的影響

將裂縫中的靜水壓力的作用加入裂縫邊坡的討論中，考慮無量綱的靜水壓力系數K為0.2～0.6的條件，結果如圖7所示，在坡度較大的裂縫中，同樣產生較大的折減，同樣針對c/γH=3、β=80°、φ=10°的邊坡，在考慮空隙水壓力K=0.3時，邊坡Fs=1.02，折減程度為14%，同樣顯著。因此對裂縫邊坡設計時，不可忽略裂縫存在的靜水壓力的作用。

圖7 靜水壓力作用下裂縫邊坡的安全系數Fig.7 Safety factor of slope with crack under hydrostatic pressure

3 裂縫位置搜索及圖表應用

3.1 裂縫位置的搜索

基于表1的工況，進行了3種條件下的邊坡破壞臨界裂縫的位置與深度，計算結果如表3所示?？梢悦黠@看出，更低的安全系數將導致距離破址更遠的裂縫位置x/H與距離坡頂面更大的裂縫深度y/H。為更加形象展示計算結果，將滑裂面破壞形態進行展示，即分別考慮khG=0.2g，上部荷載系數qt=0.26以及裂縫靜水壓力系數K=0.2與僅考慮自身重力破壞的邊坡臨界滑裂面進行對比分析，并將其裂縫的位置、以及裂縫深度進行對比分析。結果如圖8所示?？梢悦黠@看出，僅考慮自重條件下，裂縫的位置更加靠近坡頂面，而考慮不同條件的破壞工況時，裂縫出現在較遠處。為驗證裂縫位置、深度的準確性，采用有限元軟件Optum G2計算結果進行了對比，對比兩類工況：即khG=0.2g、qt=0.26的兩類，從圖9可知，計算得到的臨界滑裂面與軟件結果幾乎一致，由此再次驗證本文公式推導、計算結果均具備合理性。

圖8 典型條件下臨界狀態滑裂面形態Fig.8 The morphology of sliding crack surface in critical state under typical conditions

圖9 本文結果與Optum G2的破壞面對比Fig.9 Comparison between the current analytical results and Optum G2 results

3.2 圖表選取運用

圖5～圖7可作為工程設計時安全系數的初步評估。為反映圖表的有效性，針對一個案例進行示范。如假設實測β=45°、φ=20°、c=20 kPa、γ=17 kN/m3、H=12 m的邊坡，在不考慮任何工況條件時，根據圖5(b)顯示，Fs=1.27時邊坡處于安全狀態，而針對khG=0.2g時可讀取Fs=0.95，處于不安全的狀態，因此設計時需要進行更多考慮。

4 結論

邊坡后緣往往分布著裂縫，當前對于裂縫邊坡的研究尚且缺乏，為準確分析裂縫邊坡安全系數，將LAM、強度折減法、擬靜力法綜合研究推導。計算時，用蒙特卡洛隨機搜索法與序列二次優化法結合進行數值計算。得出如下結論。

(1)采用蒙特卡洛與SQP算法相結合可以有效解決先前算法的缺陷，且結果精度更加準確，計算速度更快。

(2)三類典型破壞條件：地震、靜水壓力、上部荷載的存在均會對邊坡的安全系數產生一定程度的影響，而對于坡度較高、內摩擦角較小的邊坡影響更為顯著。

(3)對臨界破壞狀態下裂縫的位置、深度還原，結果表明，相同參數條件下，產生較小的安全系數，而較小的安全系數易產生更深的臨界裂縫深度以及更遠的臨界裂縫距離。

(4)針對不同條件的計算結果與傳統文獻經典解進行了對比分析，其精度優于傳統分析結果，且計算結果與有限元結果基本吻合，因此繪制了大量安全系數的圖表，可作為快速評估裂縫邊坡安全系數的參考。