基于重力增加法的邊坡失穩(wěn)破壞全過(guò)程模擬

，，， ，

(武漢大學(xué) a.水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室；b.水工巖石力學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，武漢 430072)

1 研究背景

邊坡失穩(wěn)破壞[1-2]是一個(gè)漸進(jìn)變化的過(guò)程。在外力作用下，坡體內(nèi)部應(yīng)力狀態(tài)不斷發(fā)生改變，巖體部分范圍出現(xiàn)應(yīng)力集中，產(chǎn)生局部滑裂面。隨著局部裂縫延伸，形成一條貫通坡體的裂縫，最終裂縫上部的巖體破碎滑落堆積成穩(wěn)定狀態(tài)。傳統(tǒng)剛體極限平衡法由于模型簡(jiǎn)易、計(jì)算高效而廣泛運(yùn)用于邊坡穩(wěn)定分析領(lǐng)域[3-7]，但需要人為預(yù)先假定滑裂面的位置和形狀，并只能模擬邊坡臨界失穩(wěn)狀態(tài)。

近年來(lái)，隨著計(jì)算機(jī)軟件技術(shù)的不斷成熟，有限元法、離散元法等在邊坡失穩(wěn)分析領(lǐng)域中迅速發(fā)展。重力增加法作為極限平衡分析的一種工具通常被用來(lái)結(jié)合有限元法、離散元法評(píng)估邊坡的穩(wěn)定性?？祦喢鞯萚8]將重力增加法運(yùn)用于有限元法中，對(duì)邊坡破壞過(guò)程中的位移場(chǎng)進(jìn)行分析，得到了邊坡失穩(wěn)破壞的最危險(xiǎn)滑動(dòng)面；王述紅等[9]通過(guò)流形元程序采用重力增加法模擬了塊狀巖體邊坡的大變形破壞過(guò)程；Li等[10]將重力增加法應(yīng)用于RFPA程序中，開(kāi)發(fā)了RFPA-GIM程序，對(duì)實(shí)際工程中的邊坡破壞進(jìn)行了數(shù)值模擬；Utili等[11]將重力增加法應(yīng)用于離散元法中，模擬了理想化巖質(zhì)邊坡的破壞過(guò)程。有限元法雖然能反映巖體在變形過(guò)程中的應(yīng)力-應(yīng)變關(guān)系，但無(wú)法模擬巖體在破壞過(guò)程中復(fù)雜的裂縫演化過(guò)程。而離散元法雖然能夠模擬邊坡巖體的大變形，呈現(xiàn)邊坡破壞的動(dòng)態(tài)演化過(guò)程，卻無(wú)法反映邊坡巖體破壞之前的連續(xù)狀態(tài)。

為了模擬邊坡巖體從連續(xù)狀態(tài)向非連續(xù)狀態(tài)轉(zhuǎn)化的漸進(jìn)破壞全過(guò)程，部分學(xué)者開(kāi)始采用連續(xù)-離散耦合分析方法[12-13]。Intrieri 等[14]運(yùn)用連續(xù)-離散耦合分析方法研究了Torgiovannetto di Assisi滑坡的誘發(fā)機(jī)制和演化過(guò)程；Elmo等[15]和Havaej等[16]采用連續(xù)-離散耦合程序ELFEN研究了不同巖質(zhì)邊坡的破壞過(guò)程；Mahabadi等[17]開(kāi)發(fā)了程序Y-Geo，對(duì)比分析了海蝕作用下均質(zhì)巖體和節(jié)理巖體的漸進(jìn)破壞過(guò)程。

本文嘗試將重力增加法(Gravity Increase Method，GIM)應(yīng)用于連續(xù)-離散耦合分析方法(Combined Finite-Discrete Element Method，F(xiàn)DEM)中，得到基于重力增加法的連續(xù)-離散耦合分析方法(FDEM-GIM)，即在有限元法中引入斷裂力學(xué)的內(nèi)聚力模型，將界面單元插入邊坡模型的表層巖體中，建立連續(xù)-離散耦合的邊坡計(jì)算模型，采用重力增加法對(duì)邊坡臨界破壞狀態(tài)進(jìn)行分析，依據(jù)邊坡系統(tǒng)總動(dòng)能突增判斷邊坡的極限狀態(tài)，計(jì)算得到邊坡的安全系數(shù)和滑面的位置和形狀。以紅石巖堰塞體邊坡工程為實(shí)例，對(duì)比了FDEM-GIM與剛體極限平衡法的邊坡穩(wěn)定分析結(jié)果，驗(yàn)證了基于重力增加法的連續(xù)-離散耦合分析方法進(jìn)行邊坡穩(wěn)定分析的合理性；通過(guò)繼續(xù)增加重力加速度，模擬了邊坡達(dá)到臨界失穩(wěn)狀態(tài)之后的后續(xù)破壞過(guò)程，得到了邊坡的最終滑落堆積方量。

2 基本理論與計(jì)算方法

2.1 連續(xù)-離散耦合分析方法

本文將內(nèi)聚力模型[18-19]用于連續(xù)-離散耦合分析方法中，模擬巖石材料的開(kāi)裂過(guò)程。假設(shè)巖石為膠凝材料[20]，在巖石破壞過(guò)程中實(shí)體單元只發(fā)生彈性變形，損傷和斷裂發(fā)生在界面單元上。在各實(shí)體單元之間插入無(wú)厚度四節(jié)點(diǎn)界面單元，實(shí)現(xiàn)短暫的“連續(xù)”效果，為了便于示意，界面單元被拉開(kāi),如圖1所示。采用考慮單元?jiǎng)偠韧嘶膽?yīng)力-分離準(zhǔn)則來(lái)模擬裂縫的產(chǎn)生和擴(kuò)展。

圖1 界面單元與實(shí)體單元的共節(jié)點(diǎn)方式Fig.1 Mode of common nodes between cohesive elements and elastic elements

考慮到巖石等準(zhǔn)脆性材料的破壞大多是由拉斷和剪斷導(dǎo)致，因此本文采用二次應(yīng)力準(zhǔn)則定義裂縫的起裂，即

(1)

巖石材料抗剪強(qiáng)度的計(jì)算采用帶拉斷的Mohr-Coulomb準(zhǔn)則，即

(2)

式中：c為材料的內(nèi)聚力；φ為材料內(nèi)摩擦角。

當(dāng)界面單元開(kāi)始出現(xiàn)損傷后，加載仍然存在，直至界面單元完全失效，產(chǎn)生裂縫。在裂縫不斷擴(kuò)展時(shí)，界面單元的本構(gòu)關(guān)系為

(3)

式中：kn，ks分別為界面單元的法向剛度、切向剛度；δn，δs分別為界面單元的法向應(yīng)變、切向應(yīng)變；D為無(wú)量綱的損傷因子，當(dāng)D=1時(shí)，界面單元失去承載能力，可將界面單元從巖石材料中剔除。

采用基于線性軟化的Benzeggagh-Kenane準(zhǔn)則，界面單元的應(yīng)力-分離曲線如圖2，對(duì)應(yīng)公式為

(4)

圖2 界面單元應(yīng)力-分離曲線Fig.2 Stress-separation curves of cohesive elements

2.2 FDEM-GIM邊坡失穩(wěn)破壞分析方法

本文采用大型有限元軟件ABAQUS進(jìn)行邊坡失穩(wěn)破壞過(guò)程的數(shù)值模擬。為真實(shí)反映邊坡模型的初始狀態(tài)，對(duì)邊坡進(jìn)行地應(yīng)力平衡計(jì)算。施加重力作用，計(jì)算時(shí)忽略上部巖體及節(jié)理裂隙的開(kāi)裂行為。

重力增加法常被用于計(jì)算邊坡安全系數(shù)，其基本原理為在控制巖體的抗剪強(qiáng)度參數(shù)c和φ不變的前提下，通過(guò)逐漸增加重力加速度(即增大自重力作用)，使得邊坡處于臨界失穩(wěn)狀態(tài)，此時(shí)對(duì)應(yīng)的重力加速度(即臨界重力加速度glimit)與實(shí)際重力加速度g0的比值即為邊坡的安全系數(shù)，即

(5)

式中：glimit為邊坡處于臨界破壞狀態(tài)時(shí)的重力加速度；g0為邊坡實(shí)際重力加速度，通常取g0=9.8 m/s2；Fs為邊坡的安全系數(shù)。

本文邊坡臨界失穩(wěn)狀態(tài)的判別標(biāo)準(zhǔn)是邊坡的系統(tǒng)總動(dòng)能發(fā)生突增，其計(jì)算公式為

(6)

式中：EK為邊坡的系統(tǒng)總動(dòng)能；ρ，v，V分別為積分點(diǎn)處的密度、速度和體積。

文中將連續(xù)-離散耦合分析方法與重力增加法相結(jié)合，采用FDEM-GIM分析方法模擬邊坡失穩(wěn)破壞的全過(guò)程，具體流程如圖3所示。

圖3 邊坡失穩(wěn)破壞模擬流程Fig.3 Flow chart of simulation of slope instability

3 工程實(shí)例分析

3.1 紅石巖堰塞體概況

受云南省魯?shù)榭h發(fā)生在2014年8月3日的地震影響，在魯?shù)榭h火德紅鄉(xiāng)李家山村和巧家縣包谷垴鄉(xiāng)紅石巖村交界的牛欄江干流上，兩岸山體大規(guī)模塌方形成堰塞湖，震后地形地貌如圖4。

圖4 地震后地形地貌Fig.4 Topographic features after earthquake

堰塞湖集水面積11 832 km2，堰塞體位于紅石巖水電站取水壩下游600 m處，在取水壩與廠房之間。堰塞體頂部高程1 222 m，估算堰塞體總方量約1 200萬(wàn)m3，兩岸邊坡坡高在800～1 100 m之間。由于震后邊坡巖性強(qiáng)度較為薄弱，在外力作用下，可能導(dǎo)致兩岸邊坡繼續(xù)失穩(wěn)滑塌，對(duì)坡腳堆積體和工程安全產(chǎn)生不利影響。因此,有必要對(duì)兩岸邊坡的穩(wěn)定性進(jìn)行分析，研究邊坡失穩(wěn)破壞的全過(guò)程。

圖5 剖面g-g邊坡計(jì)算模型和邊坡材料分區(qū)Fig.5 Computational model and material partition of section g-g

3.2 計(jì)算模型與材料參數(shù)

計(jì)算選取紅石巖堰塞體右岸邊坡的一個(gè)典型剖面g-g斷面，剖面g-g的右岸邊坡坡高為1 020 m。進(jìn)行二維FDEM-GIM邊坡失穩(wěn)分析，計(jì)算模型如圖5(a)所示。在滑動(dòng)部分的各實(shí)體單元之間插入厚度為0的四節(jié)點(diǎn)界面單元。為了提高模擬的精度并減少計(jì)算量，僅對(duì)上部滑動(dòng)部分巖體的網(wǎng)格進(jìn)行了加密。上部巖體與卸荷裂隙網(wǎng)格尺寸為2 m，下部巖體網(wǎng)格尺寸為20 m。剖面g-g右岸邊坡模型的實(shí)體單元數(shù)量為26 376個(gè)，界面單元數(shù)量為32 775個(gè)。根據(jù)相應(yīng)工程資料，材料分區(qū)如圖5(b)所示，計(jì)算模型參數(shù)見(jiàn)表1。結(jié)合實(shí)測(cè)資料并通過(guò)實(shí)驗(yàn)室單軸壓縮數(shù)值試驗(yàn)反復(fù)調(diào)整[21-22]，獲得表層巖體和節(jié)理裂隙的界面單元力學(xué)特性參數(shù)，見(jiàn)表2。

表1 計(jì)算模型參數(shù)Table 1 Parameters of computational model

表2 表層巖體和節(jié)理裂隙的界面單元力學(xué)參數(shù)Table 2 Mechanical parameters of cohesive elements of surface rock mass and joint fractures

3.3 邊坡臨界失穩(wěn)狀態(tài)分析

圖6 剖面g-g系統(tǒng)總動(dòng)能歷時(shí)曲線Fig.6 Curve of time duration of systematic total kinetic energy of section g-g

根據(jù)FDEM-GIM計(jì)算結(jié)果提取得到了剖面g-g邊坡模型的系統(tǒng)總動(dòng)能歷時(shí)曲線，如圖6所示。當(dāng)t=16.8 s時(shí)，剖面g-g的系統(tǒng)總動(dòng)能發(fā)生突增，邊坡處于臨界失穩(wěn)狀態(tài)，對(duì)應(yīng)臨界重力加速度glimit為1.680g0(文中取g0=9.8 m/s2)，邊坡安全系數(shù)Fs=1.680。

為了驗(yàn)證計(jì)算結(jié)果的可靠性，采用剛體極限平衡法對(duì)剖面g-g進(jìn)行邊坡穩(wěn)定分析，將分析結(jié)果與FDEM-GIM的計(jì)算結(jié)果進(jìn)行對(duì)比。表3列出了FDEM-GIM與剛體極限平衡法求得的邊坡安全系數(shù)。

表3 剖面g-g右岸邊坡模型安全系數(shù)Table 3 Factor of safety of slope

由表3知，通過(guò)剛體極限平衡法計(jì)算得到的邊坡安全系數(shù)為1.616，與FDEM-GIM所得到的邊坡安全系數(shù)基本一致。圖7為剛體極限平衡法與FDEM-GIM模擬得到的邊坡滑裂面對(duì)比，當(dāng)邊坡處于臨界失穩(wěn)狀態(tài)，剛體極限平衡法與FDEM-GIM搜索得到的最危險(xiǎn)滑裂面的位置以及形狀基本吻合，因此采用基于重力增加法的連續(xù)-離散耦合分析方法分析邊坡的安全系數(shù)是可行的。

圖7 剛體極限平衡法與FDEM-GIM邊坡滑裂面對(duì)比Fig.7 Comparison of the slope slip surface obtained respectively from rigid body limit equilibrium method and FDEM-GIM

圖8 FDEM-GIM模擬邊坡開(kāi)裂-滑落-堆積全過(guò)程Fig.8 Whole process of fracture, sliding, and deposition of slope simulated by FDEM-GIM

3.4 邊坡破壞全過(guò)程模擬

圖8為FDEM-GIM模擬得到的剖面g-g右岸邊坡的失穩(wěn)破壞動(dòng)態(tài)演化全過(guò)程。當(dāng)t=10.0 s時(shí)，邊坡在g0下能夠維持其自身穩(wěn)定。隨著重力加速度的增大，表層巖體首先從其底部逐漸開(kāi)裂，當(dāng)t=16.8 s時(shí)，底部出現(xiàn)一條局部貫通裂縫，邊坡處于臨界失穩(wěn)狀態(tài)，此貫通裂縫為坡體的第一滑裂面，是邊坡潛在最危險(xiǎn)滑裂面。當(dāng)繼續(xù)增加重力加速度，使達(dá)到2g0之時(shí)，底部巖體不斷相互碰撞并破碎成若干小塊，邊坡開(kāi)始產(chǎn)生一條沿卸荷裂隙延伸的裂縫。當(dāng)t=32.0 s時(shí)，所加重力加速度仍然為2g0，邊坡逐漸產(chǎn)生一條貫通坡體的裂縫，邊坡巖體在向下運(yùn)動(dòng)的過(guò)程中，不斷撞擊破碎成若干巖塊，呈現(xiàn)出明顯的坍塌破壞模式，最終堆積成穩(wěn)定狀態(tài)。

由于邊坡失穩(wěn)破壞是一個(gè)漸進(jìn)破壞過(guò)程，當(dāng)邊坡達(dá)到極限狀態(tài)時(shí)，繼續(xù)增加重力加速度會(huì)導(dǎo)致上部巖體隨之開(kāi)裂。邊坡表層巖體的底部滑落使得上部失去了承托，在后續(xù)破壞階段，邊坡沿著第一滑裂面自下而上延伸出一條貫通坡體的裂縫。從圖8中邊坡滑落堆積的情況可以看出，邊坡在出現(xiàn)第一次滑落堆積之后仍然有后續(xù)滑落過(guò)程。因此FDEM-GIM能夠模擬邊坡在達(dá)到臨界失穩(wěn)狀態(tài)之后的開(kāi)裂行為，邊坡最終滑落方量為貫通裂縫上部體積，而剛體極限平衡法僅能用來(lái)分析邊坡臨界破壞狀態(tài)。

4 結(jié) 論

本文基于重力增加法，采用連續(xù)-離散耦合分析方法對(duì)邊坡巖體進(jìn)行失穩(wěn)破壞的全過(guò)程模擬，以紅石巖堰塞體邊坡工程為實(shí)例，將剛體極限平衡法得到的邊坡安全系數(shù)與FDEM-GIM得到的邊坡安全系數(shù)進(jìn)行對(duì)比，驗(yàn)證了基于重力增加法的連續(xù)-離散耦合分析方法運(yùn)用于邊坡失穩(wěn)破壞問(wèn)題的合理性，并模擬了邊坡巖體從未發(fā)生破壞到裂縫產(chǎn)生、擴(kuò)展、破壞直至滑落堆積的全過(guò)程,結(jié)論如下：

(1)將FDEM-GIM計(jì)算結(jié)果與傳統(tǒng)剛體極限平衡法的結(jié)果進(jìn)行對(duì)比，結(jié)果顯示2種方法得到的邊坡安全系數(shù)基本相同，最危險(xiǎn)滑裂面位置與形狀基本吻合，驗(yàn)證了基于重力增加法的連續(xù)-離散耦合分析方法進(jìn)行邊坡失穩(wěn)破壞分析的合理性。

(2)邊坡加固措施必須依賴(lài)于潛在滑動(dòng)面位置的確定，基于重力增加法的連續(xù)-離散耦合分析方法能自動(dòng)搜索獲得邊坡臨界破壞滑面，克服傳統(tǒng)剛體極限平衡法需要預(yù)先人為假定滑裂面的缺陷。剛體極限平衡法僅能用于分析邊坡臨界失穩(wěn)狀態(tài)，不能預(yù)測(cè)后續(xù)滑塊的形成。而基于重力增加法的連續(xù)-離散耦合分析方法能模擬邊坡破壞的全過(guò)程，并且可以顯示模擬后續(xù)滑坡體的形成，為邊坡防治處理提供可靠依據(jù)。