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

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

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

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

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

本文從這4種層狀巖體邊坡整體穩(wěn)定性數(shù)值分析方法的基本原理出發(fā)，同時，基于相同的邊界條件、計算收斂標(biāo)準(zhǔn)和相同的力學(xué)參數(shù)，采用這4種方法對已有文獻(xiàn)中的經(jīng)典算例進行精細(xì)分析對比驗證，系統(tǒng)總結(jié)各種求解方法的優(yōu)勢、局限性和適用范圍，以期為使用者求解層狀巖體邊坡整體穩(wěn)定性提供借鑒和參考。

1 整體穩(wěn)定性數(shù)值分析方法原理

1.1 方法一

(1)

(2)

1.2 方法二

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

圖1 層狀巖體力學(xué)模型坐標(biāo)系

1.3 方法三

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

1.4 方法四

1.4.1 整體安全系數(shù)定義

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

(3)

式中，F(xiàn)si為指定網(wǎng)格單元的點安全度;Ai為指定網(wǎng)格單元的點面積。

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

1.4.2 點安全度

1.4.2.1 層間巖石點安全度

層間巖石應(yīng)力狀態(tài)由全局坐標(biāo)系下的應(yīng)力狀態(tài)決定，其點安全度Fsm依據(jù)空間應(yīng)力狀態(tài)和Mohr-Coulomb剪切包絡(luò)線的位置關(guān)系(圖2)采用最小距離法[6]獲得，即

Fsm=AC/AB=

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

(4)

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

1.4.2.2 層面點安全度

層面應(yīng)力狀態(tài)由層面產(chǎn)狀決定的局部坐標(biāo)系上的應(yīng)力狀態(tài)所決定，其安全度Fsb基層面應(yīng)力狀態(tài)由層面產(chǎn)狀決定的局部坐標(biāo)系上的應(yīng)力狀態(tài)所決定。層面的抗剪安全度基于強度儲備法定義，即

圖2 全局坐標(biāo)系下的應(yīng)力狀態(tài)與Mohr-Coulomb剪切強度包絡(luò)線關(guān)系

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

(5)

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

1.4.2.3 層狀巖體點安全度

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

Fs=min(Fsm,Fsb) .

(6)

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

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

2.1 算例1

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

圖3 算例1幾何模型

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

圖4 算例1數(shù)值計算模型

方法容重/(kN/m3)層間巖石層面彈性模量層間巖石/GPa層面泊松比層間巖石層面黏聚力/MPa層間巖石層面內(nèi)摩擦角/(°)層間巖石層面一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種求解方法給出的整體安全系數(shù)和潛在滑動面(方法一以位移云圖給出，方法二、三以等效塑性應(yīng)變云圖給出，方法四以點安全度云圖給出)?？梢钥闯觯?種方法得到的整體安全系數(shù)大小存在差異，方法一、二和四與文獻(xiàn)[1，7]給出的原解答1.05十分接近，方法三則與該解答有較大差異。此外，這4種方法給出的潛在滑動面非常接近。從滑動面也可以看出，該邊坡表現(xiàn)為沿第3組層面滑動，其尾部后緣存在一定的拉破壞并往上延至坡頂，使得邊坡的潛在滑動面為非圓弧形滑動。這與巖體各向同性處理后獲得的圓弧形深層滑動模式(圖6)有顯著差別。

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

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

2.2 算例2

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

圖7 算例2幾何模型

圖8 算例2數(shù)值計算模型

方法容重/(kN/m3)層間巖石層面彈性模量層間巖石/GPa層面泊松比層間巖石層面黏聚力/MPa層間巖石層面內(nèi)摩擦角/(°)層間巖石層面一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種求解方法給出的整體安全系數(shù)和潛在滑動面(方法一以位移云圖給出，方法二、三以等效塑性應(yīng)變云圖給出，方法四以點安全度云圖給出)?？梢钥闯觯@4種方法得到的整體安全系數(shù)大小存在差異，方法一、二和四與原解答1.22十分接近，方法三則與該解答有較大差異。此外，這4種方法給出的潛在滑動面非常接近。方法一、二給出的滑動面表明，該順層邊坡表現(xiàn)為沿第一個層面滑動，而后在下部穿切層間巖石，從坡腳滑出的“順層-圓弧”復(fù)合滑動模式。基于隱式節(jié)理等效連續(xù)模型的強度折減法也完全反映了這一復(fù)合破壞模式。這與各向同性處理后的邊坡圓弧形滑動(圖10)差異明顯 。

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

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

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

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

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

3 結(jié) 論

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

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

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

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