基于最小勢能法的二維邊坡整體滑動方向研究

孫 加 平，詹 鳳 程，陳 延 可，康 念 坤

(江西工業工程職業技術學院，江西 萍鄉 337000)

邊坡的穩定性評價是巖土領域的研究熱點問題[1-3]。最小勢能法是一種新的評價方式，由于無需迭代，得到了一些學者的青睞[4-6]。最小勢能法總體思路是先確定系統勢能最小時的虛位移，然后根據胡克定律、莫爾-庫倫強度準則確定滑面的法向力以及剪切力，最后依據穩定系數的計算結果評價邊坡的穩定性。從目前的研究成果來看，穩定系數的定義方式主要分為兩類：① 通過抗滑力矩與下滑力矩的比值確定[7-9]；② 將力視為矢量，通過選取合適的投影軸，將作用力在滑體上投影得到抗滑力與下滑力，二者的比值確定穩定系數[10-13]。其中，第一種定義方式僅適用于圓弧形滑面，而對任意形狀滑面適用性較差。并且，該種定義方式將力的運算視為標量，忽略了力的矢量特性。因此，第二種定義方式應用較為廣泛，該種定義方式考慮了力為矢量這一要素，但對于投影軸(即滑體的整體運動方向)的合理選取研究不夠完善。同時，剪切勢能的計算模型忽略了微元體間作用力的影響，理論上缺乏嚴謹性。

因此，針對利用最小勢能法分析二維邊坡穩定性的不足，本文對剪切勢能計算模型、滑體整體運動趨勢開展研究，以期尋找最為合理的滑體運動方向，完善最小勢能理論在邊坡穩定性分析中的應用。為驗證文中研究成果的可靠性，計算分析了兩個邊坡的穩定系數以及滑動方向，并將分析結果與其他算法進行對比。

1 二維邊坡最小勢能分析方法

1.1 滑移面彈性勢能

如圖1所示的任意形狀滑動面邊坡，其方程為y=fx。建立的模型采用如下假定[14-15]：① 滑面法向的彈性變形可用剛度系數為k的彈簧模擬，且k=mdl；② 邊坡受到一個外力R=R1，R2產生微小許可的位移d=d1，d2，該位移稱為虛位移，且其使得滑體系統勢能最小。

圖1 任意形狀滑動面Fig.1 Slip surface of arbitrary shape

通過引入彈性地基梁模型，可得滑面上儲存的彈性勢能Ve為

(1)

式中：m為與土體性質有關的常數，kN/m3；η為滑面外法線方向矢量；dl為滑面微段弧長，m。

1.2 滑移面剪切勢能

土體受剪切變形影響的深度為h，則儲存的剪切勢能為[7]

(2)

則式(2)可以轉化為

(3)

式中：d′為滑床與滑面產生的相對位移，m。

土體的剪切變形示意圖如圖2所示，易知：

d′=-t·d

(4)

圖2 土體剪切變形示意Fig.2 Schematic diagram of soil shear deformation

通過式(2)～(3)可知，剪切勢能計算的關鍵在于滑面與滑床相互滑移時產生的剪應力τ′。從現有的研究來看，對于切應力τ′的計算有兩種方式：文獻[7]認為τ′沿著滑動面均勻分布，可通過滑體力矩平衡方程確定τ′；文獻[13，15]以滑體任意微元體為研究對象，通過虛位移方向的靜力平衡方程確定τ′。一般而言，滑面上任意一點的剪應力并不相同，且文獻[7]的計算方式僅僅適用于圓弧形滑面，應用具有局限性。同時，滑體中單元體的受力情況較為復雜，任何關于條間力的假定都會給計算結果帶來影響，而文獻[13，15]忽略了條間力的影響，理論上略顯不嚴謹。

鑒于此，本文提出了一種新的滑面剪應力計算方法，通過滑體整體沿著虛位移方向的靜力平衡關系確定滑面上的剪應力τ′，這樣處理避免了條間力的假定，且適用于任意形狀的滑移面，即

(5)

(6)

根據式(3)、(4)以及(6)可得：

(7)

1.3 滑體外力勢能與總勢能

滑體系統儲存的總勢能V為

(8)

根據最小勢能原理可知，系統勢能最小虛位移d應滿足

(9)

1.4 穩定系數的計算

根據力與虛位移的關系，便可求得滑面上的法向力N以及剪切力T，即

(10)

式中：c為土體的凝聚力，kPa；φ為土體的內摩擦角，°；σ為滑面上的法向應力，kPa。

本文將穩定系數K定義為滑體上的作用力在整體運動方向上投影得到的抗滑力與下滑力的比值，即

(11)

式中：Fanti為抗滑力，kN；Fslid為下滑力，kN；s為滑體整體運動方向的單位矢量。

2 3種整體下滑方向模型

從式(11) 可知，運用最小勢能法進行邊坡穩定性分析，不同的整體下滑方向得到的穩定系數是不同的。通過分析計算，文中對3種整體下滑方向進行了探究。

(1) 方法1。由于系統產生一個微小位移之后，處于穩定的平衡狀態，因此將滑體的虛位移方向視為邊坡的整體運動方向，即

s=d

(12)

(2) 方法2。由于邊坡會沿著抵抗力最小的方向發生運動，并且坡體內部的應力會自動調整以發揮最大的抗滑能力[16]。顯然，對于二維邊坡的穩定性分析而言，滑面上一點處的抗滑力應為該點的切線方向。由于滑動方向與抗滑方向相反，則滑體整體運動方向為各微面抗滑方向矢量合成的相反方向。因此，滑面上的抗滑力在坐標軸上的分量為

(13)

則滑體整體運動方向的單位矢量s為

(14)

(3) 方法3。將滑體的整體運動方向視為剪入點C到剪出點A的方向，即

s=xA-xC，yA-yC

(15)

式中：xA、xC為A、C的橫坐標；yA、yC為A、C的縱坐標。

3 算例驗證

3.1 算例1

邊坡的幾何尺寸如圖3所示，坡高為20 m，滑面為圓弧形，邊坡的幾何尺寸信息已在圖中示出。已知：土體c=15 kPa，φ=30°，容重γ=18.4 kN/m3。運用本文方法、極限平衡法得到的計算結果如表1所示。

圖3 邊坡算例1示意圖(尺寸單位：m)Fig.3 Schematic diagram of slope calculation example 1

從表1的計算結果可以看出：最小勢能法得到的結果與極限平衡法較為接近，相對誤差(與Bishop法對比)分別為1.0%、0.2%以及4.7%，滿足工程許可的范圍。并且，基于最小勢能法的3種計算方式得到的穩定系數相對誤差亦在工程允許范圍內，表明文中的3種計算方式在該算例中均是適用的。表1中方法1、方法2給出的下滑方向角度數值相同，但穩定系數不同。經過計算可知，方法1、方法2得到的角度結果在小數點第13位不同，而在MATLAB程序中參與運算的位數為30位，故而出現了上述情況。由于Bishop法的計算精度得到了業界的認可，多個行業的規范均采用該種計算方法，因此從計算結果來看，方法2的整體下滑運動方向模型更為合理。

表1 穩定系數計算結果Tab.1 Stability factor calculation results

為了進一步驗證方法2是否最為合理，下文將通過兩種方式進行驗證。

3.1.1固定滑面方式

在上述算例的基礎上，保持邊坡的幾何形狀以及滑動面不變，通過改變土體的c、φ值來計算不同方法下的邊坡穩定系數，并將得到的結果與Bishop法進行對比，如圖4所示。

圖4 穩定系數K與c、φ關系Fig.4 Relationship between stability factor K and c、φ

圖4表示了邊坡的穩定系數與土體凝聚力c、內摩擦角φ之間的關系。從圖4(a)中可以看出：在保持γ=18.4 kN/m3、φ=30°不變的的情況下，穩定系數隨著c的增大而呈現線性增大的趨勢。從方法1～3的計算結果來看，方法2的計算結果與Bishop法幾乎是一致的，相對誤差不到1%，而方法2，3與Bishop法的相對誤差分別在6.7%、10%范圍內。從而可知，基于最小勢能法的第二種計算方式最為合理。并且以基于方法2得到的計算結果穩定性較好，而方法1，方法3的計算結果隨著c值的減小，其結果與Bishop法的相對誤差越來越大，但是當c降至0時，方法1，3與Bishop法的相對誤差仍在工程許可范圍內。

從圖4(b)中可以看出：在保持γ=18.4 kN/m3、φ=30°不變的的情況下，從定性的角度來看，各種方法下的穩定系數的計算結果均隨著內摩擦角的增大而逐步增加，表明文中計算方法的合理性。從定量的角度來看，方法1的計算結果與Bishop法幾乎是一致的，而方法1，3的計算結果只有在內摩擦角φ處于25°～30°的范圍時，其計算結果與Bishop法較為一致，與Bishop法相對誤差為1%～9.2%，滿足工程上的精度要求；但當內摩擦角φ處于5°～20°時，與Bishop法相對誤差范圍為22.3%～118.4%，相對誤差較大。可見，內摩擦角對方法1，3的穩定系數的影響較為敏感，導致了計算結果穩定性差。作者認為這主要與邊坡整體下滑模型有關。同時，從圖4(b)中可以看出，方法2的計算結果與Bishop法幾乎是一致的。因此，將邊坡的整體運動方向視為各微面抗滑方向矢量合成的相反方向是合理且可行的。

3.1.2非固定滑面方式

上文在滑面固定的情況下，通過改變土體的抗剪強度指標，計算各種方法的穩定系數，以此來探究二維邊坡的整體運動最合理方式。下面將使滑體的滑入點、滑出點位置不變，通過設置不同的滑面來探究上述問題，需要注意的是這些變化的滑面的圓心均在滑入點與滑出點組成線段的垂直平分線上，選取的圓心半徑依次為：37.24，40.39，48.01，63.19，77.01，97.15，107.52 m，各種方法的穩定系數如圖5所示。

從圖5可以看出：基于最小勢能原理的3種分析方式與Bishop的結果較為接近，方法1與其相對誤差為1.1%～17.3%，方法2與其相對誤差為5%內，方法3與其相對誤差為10%以內。從定量的角度來看，方法2的計算方式更為合理。其次，從穩定系數與滑面半徑的變化關系來看，只有方法2與Bishop法保持一致，即隨著滑面深度的增加，穩定系數逐漸提高。當然，Bishop法不能作為評價一種新的邊坡穩定性分析方法的絕對標準，但其被土木、水利等行業寫進了規范。因此，將本文方法同其進行對比，具有一定的說服力。

圖5 非固定滑面下的各種方法的穩定系數Fig.5 Stability factors of various methods under non fixed sliding surface

綜上，基于最小勢能原理的方法2計算模型，即將邊坡的整體運動方向視為各微面抗滑方向矢量合成的相反方向是更加合理且有效的，可以用于邊坡的穩定性分析。

3.2 算例2

為驗證本文算法對其余滑動面的適用性，將其應用于非圓弧形滑動面。如圖6所示的滑動面為折線的均質邊坡[14]。已知：土體凝聚力c=20 kPa，φ=15°，容重γ=30 kN/m3。

圖6 邊坡算例2示意(尺寸單位：m)Fig.6 Schematic diagram of slope calculation example 2

運用文中確定的邊坡整體滑動方向模型計算可得穩定系數為0.683，下滑方向s與x軸所夾銳角為41.032°，而文獻[14]得到的穩定系數為0.657，下滑方向為40.421°。采用Janbu法得到的穩定性系數為0.726。從上述計算結果可知：運用本文計算模型得到的邊坡整體滑動方向具有一定的合理性，所得結果與其余算法的相對誤差在10%以內，對于非圓弧形滑動面的邊坡穩定性分析同樣具有適用性，可以將其應用于邊坡的穩定性分析。

4 結 論

(1) 本文從滑體整體平衡的角度建立平衡方程，求解滑面上的剪應力，進而建立剪切勢能計算模型。相比以往的研究，避免了微元體條間力的假設，理論分析上更嚴謹。

(2) 本文提出了3種滑體運動方向，并以固定、非固定滑面為研究對象，從定性以及定量的角度進行對比驗證，發現方法2的整體下滑方式(邊坡的整體運動方向視為各微面抗滑方向矢量合成的相反方向)與Bishop法的結果最為接近。

(3) 本文構建的滑體運動模型不僅適用于圓弧形滑面，從整體滑動方向、穩定系數計算表結果來看，對于折線形滑面同樣具有適用性。此外，僅僅針對固定滑面的邊坡穩定性進行了研究與探討，對于任意形狀滑動面的搜索仍可以進一步研究。