層狀巖體邊坡整體安全系數數值求解方法對比分析

汪良峰 吳小剛 徐鼎平,3,4 周玉新

(1.中鋼集團馬鞍山礦山研究院有限公司；2.金屬礦山安全與健康國家重點實驗室；3.中國科學院武漢巖土力學研究所；4.巖土力學與工程國家重點實驗室)

巖土工程中的邊坡穩定性評價是相對“古老”而又十分重要的研究課題。隨著工程建設持續而快速地開展，越來越多的邊坡穩定性問題受到學者們的廣泛關注和研究。邊坡穩定性分析需要研究2個主要問題：一是潛在主滑動面，二是整體安全系數。層狀巖體作為人類工程活動最密切的地質體之一，在地殼淺部占比超過70%，因此，層狀巖體邊坡穩定性分析也同樣面臨這2個急需解決的問題。層狀巖質邊坡的一個突出特點是巖體的強度既受層面和層間巖石力學性質的協同控制，也受到邊坡臨空面與層面產狀組合關系的制約。倘若忽略層狀巖體的這一突出特點，而將其假定或等效為各向同性巖體進行處理，很可能錯估邊坡的潛在破壞模式，給出與實際相差較大的邊坡整體安全系數，進而誤導邊坡工程的開挖、加固設計與失穩治理。

隨著計算機技術的發展和彈塑性力學理論的逐步完善，各種反映層狀巖體橫觀各向同性力學特征的力學模型和模擬手段被提出并應用于層狀巖體邊坡的穩定性分析中。就數值分析方法的經驗積累和工程實用性而言，連續介質分析方法(有限元法、有限差分法等)仍是當前邊坡穩定數值分析最主要的手段，因而采用它們進行層狀巖體邊坡穩定性數值分析方法也日益增多。趙尚毅等[1]采用低強度彈塑性薄層單元模擬巖體中的軟弱夾層，采用接觸面單元模擬層面，進而對層狀巖體邊坡采用強度折減法進行穩定性分析。林杭等[2]、張社榮等[3]采用低強度彈塑性薄層單元來模擬軟硬互層巖體中的軟弱夾層，然后采用強度折減法研究了軟弱夾層產狀、邊坡角與邊坡整體安全系數、潛在破壞模式間的關系。吳順川等[4]則基于遍布節理模型，通過對薄層狀巖質邊坡層間巖石和層面折減相同系數進行了穩定性問題研究?？傮w而言，當前這些層狀巖體邊坡穩定性數值模擬方法是在客觀認識層狀巖體結構組成對其力學性質影響，將層狀巖體視為層間巖石和層面的復合材料的前提下，采用不同模擬手段來反映其力學性質的層面依賴性特征及其對邊坡穩定的影響。在這些研究的基礎上，進而形成了4種層狀巖體邊坡穩定性的安全系數數值求解方法，即基于層面接觸單元的強度折減法、基于層面薄層單元的強度折減法、基于隱式節理等效連續模型的強度折減法和基于層狀巖體點安全系數的面積加權平均法。其中，基于層狀巖體點安全系數的面積加權平均法是作者在層狀巖體橫觀各向同性力學模型基礎上推導得到其點安全度表達式之后，借鑒楊濤等[5]的各向同性巖體邊坡整體安全系數面積加權法思想在本文中首次提出。這些研究對加深層狀巖體邊坡變形破壞的認識和促進層狀邊坡穩定定量評價的發展做出了重大貢獻。但是，當前對層狀邊坡穩定性分析往往是孤立地采用上述某一方法具體開展，計算結果的合理性也多采用層狀巖體各向同性化后的經典條分法計算結果來驗證。事實上，邊坡穩定性數值分析方法由于考慮了巖土體的本構關系及變形對計算結果的影響，它實質上是與基于剛體極限平衡理論的經典條分法并行的一種分析方法。而且，層狀巖體邊坡的經典條分法在固有的條間力假設情況下又附加了巖體力學性質各向同性這一假設。顯然，對于同一分析對象，這兩者計算結果間的相互校核不僅在理論上缺乏嚴密性，在工程上也缺乏說服力。因此，對現有順層巖質邊坡穩定安全系數數值求解方法進行系統、全面地橫向對比研究，對于客觀認知和準確評價層狀巖體邊坡穩定十分必要。

本文從這4種層狀巖體邊坡整體穩定性數值分析方法的基本原理出發，同時，基于相同的邊界條件、計算收斂標準和相同的力學參數，采用這4種方法對已有文獻中的經典算例進行精細分析對比驗證，系統總結各種求解方法的優勢、局限性和適用范圍，以期為使用者求解層狀巖體邊坡整體穩定性提供借鑒和參考。

1 整體穩定性數值分析方法原理

1.1 方法一

(1)

(2)

1.2 方法二

基于層面薄層單元的強度折減法將層間巖石處理為各向同性巖體，層面采用薄層實體單元模擬，層間巖石和薄層均采用各向同性的理想彈塑性力學模型。求解安全系數時，同時對層間巖石和薄層抗剪強度參數折減相同的系數并計算至最終狀態，不斷試算直至計算不收斂，以此時的折減系數作為邊坡整體安全系數，以此時等效塑性云圖貫通所確定的滑動面為邊坡潛在滑動面。層狀巖體力學模型中，局部坐標系由層面傾向d、走向s和法向n構成；全局坐標系由正東向x,正北向y和豎直向z組成，α和β分別為走向角和傾向角。見圖1。

圖1 層狀巖體力學模型坐標系

1.3 方法三

基于隱式節理等效連續模型的強度折減法在將層狀巖體視為層間巖石和層面的復合材料的基礎上，分別基于全局坐標系和局部坐標系(圖1)建立基于線性Mohr-Coulomb準則和最大拉應力準則的層間巖石和層面力學模型，而后合成得到層狀巖體的橫觀各向同性理想彈塑性模型，即Ubiquitous-Joint模型[6]；以此模型為基礎，再分別對層間巖石和層面抗剪強度參數折減相同的系數，不斷試算直至不收斂，以此時的折減系數作為邊坡整體安全系數，以此時等效塑性云圖貫通所確定的滑動面為邊坡潛在滑動面。Ubiquitous-Joint模型本質上是考慮層狀巖體強度各向異性的一類隱式節理等效連續模型，該模型通過特定的幾何參數(層面傾向、傾角)和強度參數(層面黏聚力和摩擦角)來反映定向節理/層面的力學性質，因而無需在網格模型里給出特定的接觸面單元或實體單元來模擬層面。

1.4 方法四

1.4.1 整體安全系數定義

基于層狀巖體點安全系數的面積加權平均法同樣是基于Ubiquitous-Joint模型，不同的是基于該模型在未進行強度參數折減的情況下進行彈塑性計算獲得邊坡體內的真實應力場，而后插值獲得臨界滑動面每一個網格質心上的應力分量，再基于其上的應力分量，求得每一個網格單元的點安全系數，最后基于點安全系數的面積加權法求得應力-應變分析獲得的邊坡整體安全系數，即

(3)

式中，Fsi為指定網格單元的點安全度;Ai為指定網格單元的點面積。

基于式(3)這一定義，要確定邊坡的整體安全系數，勢必要先求得各網格單元的點安全度。

1.4.2 點安全度

1.4.2.1 層間巖石點安全度

層間巖石應力狀態由全局坐標系下的應力狀態決定，其點安全度Fsm依據空間應力狀態和Mohr-Coulomb剪切包絡線的位置關系(圖2)采用最小距離法[6]獲得，即

Fsm=AC/AB=

[2Cm+(σ1+σ3)tanφm]cosφm/(σ1-σ3) ,

(4)

式中，cm為層間巖石黏聚力,MPa；φm為層間巖石內摩擦角,(°)；σ1和σ3分別為指定網格單元的最大主應力和最小主應力,MPa。

1.4.2.2 層面點安全度

層面應力狀態由層面產狀決定的局部坐標系上的應力狀態所決定，其安全度Fsb基層面應力狀態由層面產狀決定的局部坐標系上的應力狀態所決定。層面的抗剪安全度基于強度儲備法定義，即

圖2 全局坐標系下的應力狀態與Mohr-Coulomb剪切強度包絡線關系

Fsb=(cb+σbntanφb)/τb,

(5)

式中，cb為指定網格單元的層面黏聚力,MPa；φb為指定網格單元的層面內摩擦角,(°)；τb和σbn分別為指定網格單元層面上的剪應力和法向應力,MPa。

1.4.2.3 層狀巖體點安全度

由于巖體的破壞總是沿著最薄弱路徑破壞，因此，層狀巖體點安全度Fs取Fsm和Fsb的最小值，即

Fs=min(Fsm,Fsb) .

(6)

2 典型算例對比研究及分析

基于FLAC3D[6]采用2個經典算例開展上述4種層狀巖體邊坡整體穩定性數值分析方法的對比研究。比較對象為潛在滑動面的形狀、位置及整體安全系數。在強度折減法中，失穩判據為計算不收斂。

2.1 算例1

圖3為文獻[1]提供的一個處于臨界狀態的順層巖質邊坡幾何模型。該邊坡層理面間距為10 m，層面順傾且傾角為40°,邊坡角為65°。

圖3 算例1幾何模型

4種邊坡整體穩定性數值分析方法的計算模型采用單元厚度模型，以等效進行平面應變分析。數值計算模型幾何尺寸為水平方向125 m，豎直方向95 m,見圖4。4種安全系數求解方法所采用的物理力學參數根據文獻[1,7]給定并列于表1。

圖4 算例1數值計算模型

方法容重/(kN/m3)層間巖石層面彈性模量層間巖石/GPa層面泊松比層間巖石層面黏聚力/MPa層間巖石層面內摩擦角/(°)層間巖石層面一25-10kn=10 000 kN/mks=3 846 kN/m0.2-1.00.123824二2517100.01 GPa0.20.31.00.123824三、四25-10-0.2-1.00.123824

注：kn為層面法向剛度，ks為層面剪切剛度。

圖5為4種求解方法給出的整體安全系數和潛在滑動面(方法一以位移云圖給出，方法二、三以等效塑性應變云圖給出，方法四以點安全度云圖給出)?？梢钥闯?，4種方法得到的整體安全系數大小存在差異，方法一、二和四與文獻[1，7]給出的原解答1.05十分接近，方法三則與該解答有較大差異。此外，這4種方法給出的潛在滑動面非常接近。從滑動面也可以看出，該邊坡表現為沿第3組層面滑動，其尾部后緣存在一定的拉破壞并往上延至坡頂，使得邊坡的潛在滑動面為非圓弧形滑動。這與巖體各向同性處理后獲得的圓弧形深層滑動模式(圖6)有顯著差別。

圖5 算例1中4種求解方法的邊坡整體安全系數及滑動面

圖6 各向同性處理后的算例1邊坡滑動面

2.2 算例2

文獻[3]基于Sarma 極限平衡法和有限元強度折減法探討了多個復雜多層軟弱夾層邊坡巖體的破壞機制及穩定性特征。選取該文中的層面傾角和坡角都為45°、層面間距為30 m的順層巖質邊坡(圖7)為研究對象來進一步驗證這4種安全系數數值求解方法給出的邊坡整體安全系數和潛在滑動面的差異。仍采用單元厚度的模型，計算網格模型幾何尺寸為水平方向500 m，豎直方向200 m，見圖8。4種方法所采用的物理力學參數根據文獻[3]給定并匯總于表2。

圖7 算例2幾何模型

圖8 算例2數值計算模型

方法容重/(kN/m3)層間巖石層面彈性模量層間巖石/GPa層面泊松比層間巖石層面黏聚力/MPa層間巖石層面內摩擦角/(°)層間巖石層面一26.5-2kn=555 000 kN/mks=15 000 kN/m0.2-0.30.053019.3二26.51020.26 GPa0.20.450.30.053019.3三、四26.5-2-0.2-0.30.053019.3

注：kn為層面法向剛度，ks為層面剪切剛度。

圖9為4種求解方法給出的整體安全系數和潛在滑動面(方法一以位移云圖給出，方法二、三以等效塑性應變云圖給出，方法四以點安全度云圖給出)?？梢钥闯?，這4種方法得到的整體安全系數大小存在差異，方法一、二和四與原解答1.22十分接近，方法三則與該解答有較大差異。此外，這4種方法給出的潛在滑動面非常接近。方法一、二給出的滑動面表明，該順層邊坡表現為沿第一個層面滑動，而后在下部穿切層間巖石，從坡腳滑出的“順層-圓弧”復合滑動模式?；陔[式節理等效連續模型的強度折減法也完全反映了這一復合破壞模式。這與各向同性處理后的邊坡圓弧形滑動(圖10)差異明顯 。

圖9 算例2中4種求解方法的邊坡整體安全系數及滑動面

圖10 各向同性處理后的算例2邊坡滑動面

綜上所述，4種層狀巖體邊坡安全系數的數值求解方法給出的結果都能滿足實際需求，但是每一種方法都存在固有缺陷。例如基于層面接觸單元的強度折減法和基于層面薄層單元的強度折減法由于要將層面在網格模型中直接用幾何形式表現出來，故只能針對層面間距較大、邊坡巖體結構相對比較簡單的順層巖質邊坡，否則數值求解比較復雜的以薄層巖體(圖11)為主的順層巖質邊坡安全系數時，網格模型建立比較困難。

圖11 兩河口水電站廠址區域陡傾薄層巖體

現有算例也表明，基于隱式節理等效連續模型的強度折減法求解的層狀巖體邊坡穩定的整體安全系數與其他方法差異較大，這很可能是由于在層狀巖體破壞過程中，層面強度與層間巖石強度的調動與發揮過程存在明顯不同，在進行強度折減時同時除以相同的折減系數，忽略了這個事實。方法四相對而言適用性更廣，理論也相對嚴密。因為它是基于未折減強度參數獲得的應力場，對于巖體中的層面也不要求在網格模型中予以顯性表達，因此，求解層狀巖體邊坡不存在理論上的明顯缺陷。但是該方法使用前需要確知潛在滑動面的位置和形狀。算例1和2之所以能采用該方法求得安全系數，是因為這2個邊坡都已接近臨界狀態，基于未折減強度參數給出的位移場或等效塑性應變場已經能較為明確地給出潛在滑動面。因此，如何在基于未折減強度參數的應力場條件下獲得非臨界狀態下的邊坡潛在滑動面是該方法求解安全系數的一個關鍵問題。上述4種安全系數數值求解方法的優劣對比見表3。

3 結 論

(1)層狀巖體是一類典型的橫觀各向同性材料，其力學性質由層面和層間巖石協同控制，并受臨空面與層面產狀組合關系的制約，因此,其邊坡潛在滑動面與均質各向同性邊坡的圓弧形滑動面差異較大。在對其進行邊坡穩定性評價時，不能籠統地以各向同性材料對待，即無論是采用何種力學模型或者模擬手段，均應通過不同的方式表現這種材料力學性質的層面依賴性特征。

表3 層狀巖體邊坡4種整體安全性數值分析方法主要特點

(2)將層狀巖體中的層面處理為實體單元和將層面處理為接觸面單元，再基于各向同性力學模型的強度折減法可以獲得相對接近的安全系數計算結果，但這類方法要將層面在網格模型中予以顯性表示，因此，只能針對層面間距較大、巖體結構相對比較簡單的層狀邊坡，對于比較復雜的層狀邊坡(薄層巖體邊坡等)，因網格模型建立比較困難則難以應用?；陔[式節理等效連續模型的強度折減法可以采用均一化網格處理復雜層狀邊坡穩定性分析，但是其獲得的整體安全系數計算結果與其他方法差別較大，表明層狀巖體在破壞過程中，其層間巖石和層面的強度調動過程和機理存在較大差異，求解安全系數時對這二者折減相同系數的做法有待商榷，尚需進一步研究。

(3)在順層巖質邊坡潛在滑動面上，基于接觸面單元的強度折減法、基于層面薄層單元的強度折減法、基于隱式節理等效連續模型的強度折減法均能給出比較一致的形狀。層狀巖體點安全系數的面積加權平均法則取決于邊坡本身的狀態，若處于臨界狀態附近，則破壞模式與基于隱式節理等效連續模型的強度折減法結果相同；若邊坡處于非臨界狀態，則需借助其他方式來確定潛在滑動面。