邊坡穩定極限分析方法（三維邊坡在對數螺旋破壞機制下的穩定性分析）

2023-04-13 11:08:10 2

摘要：

為研究三維邊坡在對數螺旋破壞機制下的穩定性，基於極限分析上限法提出了非相關聯流動法則下的速度場及安全係數的求解方法，可有效的解決相關聯法則中摩擦角過大可能導致速度場不收斂的問題。基於Michalowski構建了考慮體積應變的均質三維邊坡的對數螺旋線破壞機制，並利用MATLAB編制計算程序計算不同內摩擦角φ、不同邊坡傾角β工況下邊坡的穩定係數，分析內摩擦角、邊坡傾角對穩定係數和破壞模式的影響規律，引入強度折減技術計算邊坡安全係數，通過與Michalowski方法對比發現計算結果接近。研究結果表明：在三維邊坡對數螺旋破壞機制下相對於一定寬度的邊坡適用，且經過計算穩定係數符合規律。研究成果為邊坡的設計計算提供參考。

關鍵詞：

邊坡工程; 極限分析; 三維邊坡; 非相關流動法則;

作者簡介：

董鳴鏑(1970—),男，博士研究生，主要研究方向為工程項目管理。E-mail:[email protected];

基金：

國家重點研發計劃(2018YFC0407102);

引用：

董鳴鏑． 基於對數螺旋破壞機制的三維邊坡穩定性分析［J］. 水利水電技術( 中英文) ，2021，52( 3) : 169-176.

DONG Mingdi． Three dimensional slope stability analysis based on logarithmic spiral failure mechanism［J］. Water Ｒesources and Hydropower Engineering，2021，52( 3) : 169-176.

圖1 垂直條分法滑動土體及條柱受力示意

邊坡穩定性研究較常採用的方法主要有三種：極限平衡法、極限分析法和有限元法。極限平衡法，即條分法，是最早被提出的邊坡穩定性分析方法。1936年，極限平衡法被正式用於邊坡穩定性分析，這一方法被命名為瑞典條分法。瑞典條分法假設邊坡滑裂面為圓弧，在邊坡體表面劃分直線，使之成為一系列豎直條塊，忽略條塊間的相互作用力，對各條塊進行靜力分析，建立滑坡體的靜力、力矩平衡方程，解方程求得邊坡安全係數。而後極限平衡法研究分支不斷發展，1954年，JANBU 將邊坡滑裂面優化為非圓弧線，但該方法並不滿足力矩平衡方程。1955年，BISHOP 研究時引入條塊間的法向力，以土體抗剪強度與滑裂面剪應力的比值作為安全係數，將平衡方程不斷迭代求得結果。隨著計算機技術的飛速發展，又產生了MONGENSTERN-PRICE法 (1965)和SPENCER法 (1967)。此後，隨著有限元理的建立和發展，ZIENKIEWICZ將邊坡的穩定性分析帶入到有限元時代，後面發展成為目前應用最廣的數值分析法 。由於過程簡明，維計算方法結果往往偏保守，誤差甚至達15% ,而三維邊坡穩定性分析能更加貼近真實邊坡破壞形態。儘管如此，相對於二維計算方法，三維計算方法尚未成熟，國內外諸多學者進行了這方面的研究。

隨著理論體系的成熟，能較好的滿足工程要求，二維邊坡穩定性計算方法目前被廣泛採用 。但是二邊坡三維穩定性分析方法主要有三維極限平衡法和三維極限分析法。與二維計算方法原理類似，三維極限平衡法將邊坡劃分為系列垂直條塊，如圖1所示，利用靜力平衡方程求解邊坡的安全係數 。

在三維極限平衡分析中，垂直條塊間的力的大小、方向、作用點均未知，且單個垂直條塊前後左右均有作用力，同二維邊坡相比，未知變量個數遠大於平衡方程個數。根據Lma&Fredlund, 由於求解三維邊坡安全係數存在許多獨立變量，需要引入大量假設，以有n行m列垂直條塊的邊坡為例，需要引入8 mn個假設 。也有其他學者採用不同的未知量假設。由於垂直條塊數量眾多，方程數量龐大，算法的收斂性是個需要解決的問題 ,另外搜索最優三維滑動面也是一個巨大的挑戰。

三維邊坡穩定性分析，除了極限平衡法，另外一種較常用的方法是極限分析法。1975年，GIGER將極限分析法引入到三維豎直切坡穩定性分析中去，開創了三維邊坡極限分析的先例，但其研究方法僅適用於三維豎直切坡，並不能分析一般邊坡。基於Mohr-Coulomb準則和相關聯流動法則，MICHALOWSKI 將三維邊坡劃分為一系列滑塊，相鄰滑塊的交界面垂直於滑體的對稱軸，並採用二維邊坡的方法構建速度場，以此分析三維邊坡穩定性。FARZANEH 改進了上述方法，將其用於求解三維非均質邊坡的穩定性問題，並討論了用迭代算法求無約束和有約束問題的最小上限解。SUOBAR 用MICHALOWSKI的方法計算三維擋土牆的被動土壓力，由於此方法將三維速度場簡化成二維，其應用性有限。CHEN 引入中性面概念，將中性面劃分為一系列傾斜條塊，並採用一定方法對中性面進行空間擴展，將邊坡離散成為一系列三維傾斜條塊，應用Mohr-Coulomb準則和相關聯流動法則構建速度場，最終求出三維邊坡的安全係數，並在此基礎上提出了一種應用單純形自動搜索臨界滑動面的方法。CHEN在構建三維速度場時引入兩個假設，破壞了解的嚴謹性，但其進行的大量算例表明兩個假設引起的誤差不大。孫平 對上述方法中由於相關聯法則中摩擦角過大可能導致速度場不收斂的問題進行了研究，並提出了非相關聯流動法則下的速度場及安全係數的求解方法。MICHALOWSKI 構建了考慮體積應變的均質三維邊坡的對數螺旋線破壞機制，並提出相應研究思路。

1975年，CHEN 在其著作「Limit Analysis and Soil Plasticity」中提出一種二維邊坡旋轉破壞機制。該機構將二維邊坡視作剛體，假設其滑裂面為對數螺旋線，利用極限分析上限法推導出穩定係數的解析解，並編制計算程序求得其穩定係數。此後，對於二維均質邊坡，諸多學者採用對數螺旋線破壞模式進行穩定性分析。但實際工程中，邊坡都呈現三維形態，採用二維破壞模式低估了實際邊坡的複雜性，且採用二維破壞模式所得安全係數小於實際情況,影響實際加固和設計方案的選取。

為解決上述問題，眾多學者對三維邊坡破壞模式進行了研究，其中Michalowski引入對數螺旋線圓錐破壞模式，對二維對數螺旋線模式進行空間擴展，並考慮邊坡體的體積應變，推導三維邊坡穩定係數的解析公式，編制相應計算程序。本文將對此破壞模式進行詳細介紹。

1.1 三維邊坡極限分析理論

極限分析上限定理認為，在一個假設的滿足速度邊界條件和應變與速度相容條件的速度場中，由外荷載功率等於所消耗的內力功率所確定的荷載一定不小於實際破壞荷載。滿足上述條件的速度場叫做運動許可速度場。因此，如果能找到一個運動許可速度場，則土體內的塑性流動必將發生或早已發生。根據虛功率原理推導出上限定理為：在所有的機動容許的塑性變形速度場相對應的荷載中，極限荷載為最小。即

式中，σi*j 為由εi*j 按塑性變形法則求出的應力；T 、F 分別為作用物體上的面力和體力；Δv為速度間斷線兩側切向速度變化量；u 、εi*j 分別為速度場中的速度和應變率；c、φ為抗剪強度指標中的粘聚力和摩擦角；dA, dV, dL為積分變量，分別對應積分面元，積分體積元和積分線元。

上式所確定的力系T 及F 大於或等於真實的極限荷載。也就是說由上述方程確定的力系T 、F 是真實極限荷載的一個上限解，我們最終需要通過程序優化出來一個最小上限解。按照極限分析理論，能量耗散發生在土體塑性區內。然而，為了簡化問題，對邊坡進行穩定性分析時，經常假設邊坡體由系列剛性塊體組成，此時能量耗散只發生在速度間斷線上，而速度間斷線是各相鄰剛性塊體之間的窄過渡層。本文仍然採用剛性塊體假設，不考慮邊坡的體積應變。

1.2 破壞機制的構建

如圖2所構建的穩定性分析機構示意圖，假設一均質邊坡高度為H,寬度為B,傾角為β,土體黏聚力為c,內摩擦角為φ,重度為γ。當邊坡失穩時，以角速度ω圍繞旋轉中心O旋轉，曲線AC是式(1)確定的對數螺旋線，為坡體的滑裂面，該曲線通過坡趾C。數螺旋線A′C′確定了三維邊坡滑坡體的上界，曲線AC與曲線A′C′相交於一點。以O為原點向邊坡體做射線OA,交曲線A′C′於A′,以AA′為直徑沿垂直於紙面方向做圓，該圓以O為旋轉中心進行旋轉，直徑漸變，為上述射線與曲線A′C′和曲線AC兩交點確定的線段。圓圍繞O點旋轉與邊坡體相交形成的曲面為三維邊坡的滑裂面，整個旋轉機構為螺旋錐體，按照相關聯流動法則，其頂角夾角為2φ,整個滑坡體的上邊界A′C′和下邊界AC由式(1)和式(2)確定，至此，邊坡的三維旋轉機構構建完成

式中，ro和r′o分別為圖2所示穩定性分析機構中線段OA、OA′的長度；θo為線段OA與水平面的夾角；θ滑面AC上任意一點於機構旋轉中心連線於水平線所成夾角；φ為邊坡土體的摩擦角。

圖2 三維邊坡旋轉破壞機制

滑面AC上任意一點於機構旋轉中心連線於水平線所成夾角記為θ,每一個θ都有一個其對應的旋轉斷面的圓(如圖2中陰影部分的圓所示),則其半徑為用R表示(圓心點與旋轉中心O的連線的長度),可由下式計算

所有旋轉斷面的圓心所構成的曲線的極坐標方程可表示為

以每個斷面圓圓心為原點，做如上圖2所示局部坐標系，即以圓心點與旋轉中心的連線為y軸，y軸的垂直方向為x軸，則邊坡體內任意一點的速度為

式中，vi為滑體中任意一點在繞中心O旋轉過程中的速度；rmi為任意一點所對應斷面圓圓心旋轉中心的連線的長度；yi為任意一點在其所在局部坐標系中的縱坐標值；ω為滑體繞O點旋轉的的角速度。該模型可用於有寬度限制的邊坡穩定係數的計算。對於無寬度限制的邊坡，該破壞機制不太適用 。

如圖3所示，旋轉體體積微元可由式(7)計算得到，其中dV表示體積微元的大小，r 微元體所對應斷面圓圓心旋轉中心的連線的長度，y表示微元體在其所在局部坐標系中的縱坐標值。

圖3 旋轉體體積微元

則根據力做功功率的計算公式P=F·V,可求得滑坡體微元重力功率可表示為

式中，γ為邊坡土體的容重；vi為微元體運動的速度其數值，由公式(6)計算可得；θ為體積微元與機構旋轉中心連線O於水平線所成夾角。

通過滑坡體微元重力做功進行積分，求得邊坡滑坡體重力功率可表示為

式中，參數θc, θb, θhx1* ,x2* , y , a, d為積分計算公式中的積分邊界，其數值可由下式計算求得

在求解公式(9)時，先對y進行積分，由圖3可知，y的積分範圍在區間(θ ,θ )內為(a,y ),在區間(θ ,θ )內為(d,y ),其中x1* 、x2* 分別表示坡頂、坡面某一旋轉斷面的最大寬度。

由於不考慮體積應變，內能耗散僅發生在速度間斷面上，不考慮邊坡的體積應變，內能耗散發生在空間曲面AC上，Chen 在其著作「Limit Analysis and Soil Plasticity」中指出內部耗損率可以由該面的微分面積dS與粘聚力c以及與跨該面的切向間斷速度(vcosφ)的連乘積沿整個間斷面積分，即得到總的內部能量耗損率，其表達式為

其中參數θo, θB,θh為同上所述，積分變量dα為積微元體在局部坐標系中與原點連線同x軸之間的夾角，α 和α為積分變量dα對應的積分邊界，其數值可由公式(18)和公式(19)計算公式求得，其中參數α 、α 可由公式(10)和公式(11)計算求得

2.1 穩定係數計算

根據極限分析上限定理，外力做功功率與內能耗散功率相等即

整理式(20)後可知，穩定係數γH/c可表示成θo, θh, r′0/r0的函數，求邊坡的最小穩定係數γH/c本質上是非線性規劃問題，使用MATLAB編制計算程序，能方便的求解穩定係數。若給出土體的c,γ,還可通過優化程序求得邊坡的臨界高度。

如圖4所示，旋轉體內插入一寬度為b的塊體，當插入塊的寬度b不斷增大，三維邊坡破壞機制寬度將不斷增大。當b→∞時，有限寬度邊坡向無限寬度邊坡轉換，最終變為二維平面應變模式。插入塊的重力功率和內能耗散功率可由二維情況的重力功率和內能耗散功率乘以b獲得。

圖4 三維邊坡插入塊體示意

求解邊坡的穩定係數公式比較複雜，本文採用MATLAB軟體編製程序進行計算。求解邊坡的最小穩定係數屬於非線性約束二次規劃問題，MATLAB的fmincon函數是求解含約束優化問題的比較高效的函數之一。本文採用fmincon函數進行計算。該函數在給定的約束條件內對穩定係數γH/c進行優化，最終求得γH/c的最小值。除了用fmincon函數求解最小穩定係數，還可以採用枚舉法，先根據約束條件劃分比較粗糙的求解區間，並搜索該區間內的最小穩定係數，確定最小安全係數所在的目標區間後，細化該區間進行第二次循環，便可找到較優的穩定係數。枚舉法雖然不能非常精確的求得最小穩定係數，經過比較，其精度完全可以滿足計算要求，但是由於枚舉法需要計算每個步長對應的解，即使增大步長，工作量依然巨大，導致計算時間較長，故本文採用fmincon函數求解穩定係數。

2.1.1 穩定係數對比

為驗證結果的正確性，需要與二維邊坡穩定係數 進行對比，結果如表1所列。

表1 普通均質邊坡二維與三維極限分析穩定係數

由表1可知，當B/H=10時，三維邊坡穩定係數略大於二維情況穩定係數，這是由於當b不斷增大三維邊坡逐漸向二維邊坡模式轉變，三維穩定係數不斷趨近於二維穩定係數，同時也可表明，二維平面應變破壞機制比三維破壞機制計算結果危險 。

另將本文計算結果與Michalowski計算結果對比，列於表2至表4中。可以發現，本文計算結果與Michalowski計算結果相差無幾，且絕大部分小於後者，而極限分析上限法所求為上限解，其值越小越具有意義，從與二維計算結果及Michalowski計算結果對比可知，本文計算可行。

2.1.2 規律分析

圖4展示了不同內摩擦角及不同邊坡傾角下γH/c的變化規律。以圖4(a)為例，當內摩擦角φ和邊坡傾角β均為定值時，邊坡的穩定係數γH/c將隨著邊坡寬高比B/H的增大不斷減小，當B/H5時，γH/c減小幅度比較有限，將趨於穩定，這是由於隨著B/H的增大，三維模型逐漸趨近於二維平面應變模型。

當內摩擦角φ和邊坡寬高比B/H為定值時，γH/c隨著邊坡傾角β的增大而減小，這表明當邊坡寬度一定時，隨著邊坡傾角的增大，邊坡失穩的臨界高度不斷減小，即邊坡傾角越大越容易失穩。

表2 本文方法與Michalowski三維極限分析穩定係數對比分析(φ=0°)

表3 本文與Michalowski三維極限分析穩定係數對比分析(φ=15°)

表4 本文與Michalowski三維極限分析穩定係數對比分析(φ=30°)

當寬高比B/H和邊坡傾角β為定值時，γH/c隨著摩擦角φ的增大而不斷增大。這表明當邊坡傾角、邊坡寬度不變時，邊坡失穩的臨界高度會隨著內摩擦角的增大而增大，即土體內摩擦角越大越不容易失穩。

圖5 不同工況下的穩定係數

2.2 安全係數計算

邊坡穩定性的評定指標有很多種，如臨界高度Hcr、穩定係數N =γH/c,安全係數F ,而工程中最常用的邊坡穩定性評價指標為安全係數。現實邊坡往往具有一定的安全儲備，通過折減土體強度參數，把潛在滑動土體的抗剪強度參數折減Fs倍，邊坡恰好過渡到臨界極限平衡狀態。對滑動土體的抗剪強度進行折減，實際上就是對滑面土體的強度參數進行折減使滑面上的內部耗能減小並使其達到與外力做功相等的臨界狀態。折減的倍數反映了邊坡的安全儲備，稱之為安全係數Fs。自強度折減技術與極限分析法結合以來，眾多學者採用安全係數對二維邊坡進行了穩定性分析，並取得了非常豐富的成果，然而，採用安全係數Fs分析三維邊坡穩定的研究相對較少。本節將推導三維邊坡安全係數公式，利用MATLAB編制求解公式，利用安全係數這一評價指標對三維邊坡的穩定性進行分析

式中，Fs為剪切強度折減係數；c和φ為原始抗剪強度參數；cf和φf為折減後的抗剪強度參數。

採用這一方法時，Fs以隱式出現在求解的方程式中，計算時常需要進行迭代。將式(21)折減後的強度參數cf和φf代替式中(20)的強度參數c和φ進行迭代求解，可得Fs是關於位置參數θo, θh, r′0, r0的隱函數

為驗證該三維邊坡安全係數計算結果的正確性，選取典型邊坡計算，對比FLAC3D計算結果，列於表5。本文三維計算案例中選取固定b/H值，由表5對比結果可知，本文計算結果與有限差分軟體FLAC3D計算結果相差不大，證明了本文計算方法的有效性。

表5 普通邊坡FLAC3D與三維極限分析安全係數對比

本章基於極限分析上限法，構建了三維邊坡破壞機制，推導了穩定係數、安全係數解析公式。通過計算不同工況下邊坡的穩定係數、安全係數，主要得到如下結論：

(1)三維邊坡對數螺旋破壞機制中，邊坡傾角、內摩擦角一定時，穩定係數隨著邊坡寬高比B/H的增大先急劇減小，後趨於穩定；內摩擦角、寬高比一定時，穩定係數隨著邊坡傾角的增大而減小；寬高比、邊坡傾角一定時，穩定係數隨著內摩擦角的增大而增大。

(2)當插入塊的寬度b不斷增大時，穩定係數不斷減小，最終趨近於二維平面應變穩定係數，但仍然比二維結果大，在表明該破壞機制中的邊坡由三維狀態逐漸向二維平面應變狀態轉變的同時，也說明同等條件下二維計算結果更加危險。

(3)同等工況下，三維邊坡安全係數與FLAC3D計算結果相差無幾，表明引入安全係數有效。

水利水電技術

水利部《水利水電技術》雜誌是中國水利水電行業的綜合性技術期刊（月刊），為全國中文核心期刊，面向國內外公開發行。本刊以介紹我國水資源的開發、利用、治理、配置、節約和保護，以及水利水電工程的勘測、設計、施工、運行管理和科學研究等方面的技術經驗為主，同時也報導國外的先進技術。期刊主要欄目有:水文水資源、水工建築、工程施工、工程基礎、水力學、機電技術、泥沙研究、水環境與水生態、運行管理、試驗研究、工程地質、金屬結構、水利經濟、水利規劃、防汛抗旱、建設管理、新能源、城市水利、農村水利、水土保持、水庫移民、水利現代化、國際水利等。