,,, ,
(武漢大學(xué) a.水資源與水電工程科學(xué)國家重點實驗室;b.水工巖石力學(xué)教育部重點實驗室,武漢 430072)
邊坡失穩(wěn)破壞[1-2]是一個漸進變化的過程。在外力作用下,坡體內(nèi)部應(yīng)力狀態(tài)不斷發(fā)生改變,巖體部分范圍出現(xiàn)應(yīng)力集中,產(chǎn)生局部滑裂面。隨著局部裂縫延伸,形成一條貫通坡體的裂縫,最終裂縫上部的巖體破碎滑落堆積成穩(wěn)定狀態(tài)。傳統(tǒng)剛體極限平衡法由于模型簡易、計算高效而廣泛運用于邊坡穩(wěn)定分析領(lǐng)域[3-7],但需要人為預(yù)先假定滑裂面的位置和形狀,并只能模擬邊坡臨界失穩(wěn)狀態(tài)。
近年來,隨著計算機軟件技術(shù)的不斷成熟,有限元法、離散元法等在邊坡失穩(wěn)分析領(lǐng)域中迅速發(fā)展。重力增加法作為極限平衡分析的一種工具通常被用來結(jié)合有限元法、離散元法評估邊坡的穩(wěn)定性??祦喢鞯萚8]將重力增加法運用于有限元法中,對邊坡破壞過程中的位移場進行分析,得到了邊坡失穩(wěn)破壞的最危險滑動面;王述紅等[9]通過流形元程序采用重力增加法模擬了塊狀巖體邊坡的大變形破壞過程;Li等[10]將重力增加法應(yīng)用于RFPA程序中,開發(fā)了RFPA-GIM程序,對實際工程中的邊坡破壞進行了數(shù)值模擬;Utili等[11]將重力增加法應(yīng)用于離散元法中,模擬了理想化巖質(zhì)邊坡的破壞過程。有限元法雖然能反映巖體在變形過程中的應(yīng)力-應(yīng)變關(guān)系,但無法模擬巖體在破壞過程中復(fù)雜的裂縫演化過程。而離散元法雖然能夠模擬邊坡巖體的大變形,呈現(xiàn)邊坡破壞的動態(tài)演化過程,卻無法反映邊坡巖體破壞之前的連續(xù)狀態(tài)。
為了模擬邊坡巖體從連續(xù)狀態(tài)向非連續(xù)狀態(tài)轉(zhuǎn)化的漸進破壞全過程,部分學(xué)者開始采用連續(xù)-離散耦合分析方法[12-13]。Intrieri 等[14]運用連續(xù)-離散耦合分析方法研究了Torgiovannetto di Assisi滑坡的誘發(fā)機制和演化過程;Elmo等[15]和Havaej等[16]采用連續(xù)-離散耦合程序ELFEN研究了不同巖質(zhì)邊坡的破壞過程;Mahabadi等[17]開發(fā)了程序Y-Geo,對比分析了海蝕作用下均質(zhì)巖體和節(jié)理巖體的漸進破壞過程。
本文嘗試將重力增加法(Gravity Increase Method,GIM)應(yīng)用于連續(xù)-離散耦合分析方法(Combined Finite-Discrete Element Method,F(xiàn)DEM)中,得到基于重力增加法的連續(xù)-離散耦合分析方法(FDEM-GIM),即在有限元法中引入斷裂力學(xué)的內(nèi)聚力模型,將界面單元插入邊坡模型的表層巖體中,建立連續(xù)-離散耦合的邊坡計算模型,采用重力增加法對邊坡臨界破壞狀態(tài)進行分析,依據(jù)邊坡系統(tǒng)總動能突增判斷邊坡的極限狀態(tài),計算得到邊坡的安全系數(shù)和滑面的位置和形狀。以紅石巖堰塞體邊坡工程為實例,對比了FDEM-GIM與剛體極限平衡法的邊坡穩(wěn)定分析結(jié)果,驗證了基于重力增加法的連續(xù)-離散耦合分析方法進行邊坡穩(wěn)定分析的合理性;通過繼續(xù)增加重力加速度,模擬了邊坡達到臨界失穩(wěn)狀態(tài)之后的后續(xù)破壞過程,得到了邊坡的最終滑落堆積方量。
本文將內(nèi)聚力模型[18-19]用于連續(xù)-離散耦合分析方法中,模擬巖石材料的開裂過程。假設(shè)巖石為膠凝材料[20],在巖石破壞過程中實體單元只發(fā)生彈性變形,損傷和斷裂發(fā)生在界面單元上。在各實體單元之間插入無厚度四節(jié)點界面單元,實現(xiàn)短暫的“連續(xù)”效果,為了便于示意,界面單元被拉開,如圖1所示。采用考慮單元剛度退化的應(yīng)力-分離準(zhǔn)則來模擬裂縫的產(chǎn)生和擴展。
圖1 界面單元與實體單元的共節(jié)點方式Fig.1 Mode of common nodes between cohesive elements and elastic elements
考慮到巖石等準(zhǔn)脆性材料的破壞大多是由拉斷和剪斷導(dǎo)致,因此本文采用二次應(yīng)力準(zhǔn)則定義裂縫的起裂,即
(1)
巖石材料抗剪強度的計算采用帶拉斷的Mohr-Coulomb準(zhǔn)則,即
(2)
式中:c為材料的內(nèi)聚力;φ為材料內(nèi)摩擦角。
當(dāng)界面單元開始出現(xiàn)損傷后,加載仍然存在,直至界面單元完全失效,產(chǎn)生裂縫。在裂縫不斷擴展時,界面單元的本構(gòu)關(guān)系為
(3)
式中:kn,ks分別為界面單元的法向剛度、切向剛度;δn,δs分別為界面單元的法向應(yīng)變、切向應(yīng)變;D為無量綱的損傷因子,當(dāng)D=1時,界面單元失去承載能力,可將界面單元從巖石材料中剔除。
采用基于線性軟化的Benzeggagh-Kenane準(zhǔn)則,界面單元的應(yīng)力-分離曲線如圖2,對應(yīng)公式為
(4)
圖2 界面單元應(yīng)力-分離曲線Fig.2 Stress-separation curves of cohesive elements
本文采用大型有限元軟件ABAQUS進行邊坡失穩(wěn)破壞過程的數(shù)值模擬。為真實反映邊坡模型的初始狀態(tài),對邊坡進行地應(yīng)力平衡計算。施加重力作用,計算時忽略上部巖體及節(jié)理裂隙的開裂行為。
重力增加法常被用于計算邊坡安全系數(shù),其基本原理為在控制巖體的抗剪強度參數(shù)c和φ不變的前提下,通過逐漸增加重力加速度(即增大自重力作用),使得邊坡處于臨界失穩(wěn)狀態(tài),此時對應(yīng)的重力加速度(即臨界重力加速度glimit)與實際重力加速度g0的比值即為邊坡的安全系數(shù),即
(5)
式中:glimit為邊坡處于臨界破壞狀態(tài)時的重力加速度;g0為邊坡實際重力加速度,通常取g0=9.8 m/s2;Fs為邊坡的安全系數(shù)。
本文邊坡臨界失穩(wěn)狀態(tài)的判別標(biāo)準(zhǔn)是邊坡的系統(tǒng)總動能發(fā)生突增,其計算公式為
(6)
式中:EK為邊坡的系統(tǒng)總動能;ρ,v,V分別為積分點處的密度、速度和體積。
文中將連續(xù)-離散耦合分析方法與重力增加法相結(jié)合,采用FDEM-GIM分析方法模擬邊坡失穩(wěn)破壞的全過程,具體流程如圖3所示。
圖3 邊坡失穩(wěn)破壞模擬流程Fig.3 Flow chart of simulation of slope instability
受云南省魯?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萬m3,兩岸邊坡坡高在800~1 100 m之間。由于震后邊坡巖性強度較為薄弱,在外力作用下,可能導(dǎo)致兩岸邊坡繼續(xù)失穩(wěn)滑塌,對坡腳堆積體和工程安全產(chǎn)生不利影響。因此,有必要對兩岸邊坡的穩(wěn)定性進行分析,研究邊坡失穩(wěn)破壞的全過程。
圖5 剖面g-g邊坡計算模型和邊坡材料分區(qū)Fig.5 Computational model and material partition of section g-g
計算選取紅石巖堰塞體右岸邊坡的一個典型剖面g-g斷面,剖面g-g的右岸邊坡坡高為1 020 m。進行二維FDEM-GIM邊坡失穩(wěn)分析,計算模型如圖5(a)所示。在滑動部分的各實體單元之間插入厚度為0的四節(jié)點界面單元。為了提高模擬的精度并減少計算量,僅對上部滑動部分巖體的網(wǎng)格進行了加密。上部巖體與卸荷裂隙網(wǎng)格尺寸為2 m,下部巖體網(wǎng)格尺寸為20 m。剖面g-g右岸邊坡模型的實體單元數(shù)量為26 376個,界面單元數(shù)量為32 775個。根據(jù)相應(yīng)工程資料,材料分區(qū)如圖5(b)所示,計算模型參數(shù)見表1。結(jié)合實測資料并通過實驗室單軸壓縮數(shù)值試驗反復(fù)調(diào)整[21-22],獲得表層巖體和節(jié)理裂隙的界面單元力學(xué)特性參數(shù),見表2。
表1 計算模型參數(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
圖6 剖面g-g系統(tǒng)總動能歷時曲線Fig.6 Curve of time duration of systematic total kinetic energy of section g-g
根據(jù)FDEM-GIM計算結(jié)果提取得到了剖面g-g邊坡模型的系統(tǒng)總動能歷時曲線,如圖6所示。當(dāng)t=16.8 s時,剖面g-g的系統(tǒng)總動能發(fā)生突增,邊坡處于臨界失穩(wěn)狀態(tài),對應(yīng)臨界重力加速度glimit為1.680g0(文中取g0=9.8 m/s2),邊坡安全系數(shù)Fs=1.680。
為了驗證計算結(jié)果的可靠性,采用剛體極限平衡法對剖面g-g進行邊坡穩(wěn)定分析,將分析結(jié)果與FDEM-GIM的計算結(jié)果進行對比。表3列出了FDEM-GIM與剛體極限平衡法求得的邊坡安全系數(shù)。
表3 剖面g-g右岸邊坡模型安全系數(shù)Table 3 Factor of safety of slope
由表3知,通過剛體極限平衡法計算得到的邊坡安全系數(shù)為1.616,與FDEM-GIM所得到的邊坡安全系數(shù)基本一致。圖7為剛體極限平衡法與FDEM-GIM模擬得到的邊坡滑裂面對比,當(dāng)邊坡處于臨界失穩(wěn)狀態(tài),剛體極限平衡法與FDEM-GIM搜索得到的最危險滑裂面的位置以及形狀基本吻合,因此采用基于重力增加法的連續(xù)-離散耦合分析方法分析邊坡的安全系數(shù)是可行的。
圖7 剛體極限平衡法與FDEM-GIM邊坡滑裂面對比Fig.7 Comparison of the slope slip surface obtained respectively from rigid body limit equilibrium method and FDEM-GIM
圖8 FDEM-GIM模擬邊坡開裂-滑落-堆積全過程Fig.8 Whole process of fracture, sliding, and deposition of slope simulated by FDEM-GIM
圖8為FDEM-GIM模擬得到的剖面g-g右岸邊坡的失穩(wěn)破壞動態(tài)演化全過程。當(dāng)t=10.0 s時,邊坡在g0下能夠維持其自身穩(wěn)定。隨著重力加速度的增大,表層巖體首先從其底部逐漸開裂,當(dāng)t=16.8 s時,底部出現(xiàn)一條局部貫通裂縫,邊坡處于臨界失穩(wěn)狀態(tài),此貫通裂縫為坡體的第一滑裂面,是邊坡潛在最危險滑裂面。當(dāng)繼續(xù)增加重力加速度,使達到2g0之時,底部巖體不斷相互碰撞并破碎成若干小塊,邊坡開始產(chǎn)生一條沿卸荷裂隙延伸的裂縫。當(dāng)t=32.0 s時,所加重力加速度仍然為2g0,邊坡逐漸產(chǎn)生一條貫通坡體的裂縫,邊坡巖體在向下運動的過程中,不斷撞擊破碎成若干巖塊,呈現(xiàn)出明顯的坍塌破壞模式,最終堆積成穩(wěn)定狀態(tài)。
由于邊坡失穩(wěn)破壞是一個漸進破壞過程,當(dāng)邊坡達到極限狀態(tài)時,繼續(xù)增加重力加速度會導(dǎo)致上部巖體隨之開裂。邊坡表層巖體的底部滑落使得上部失去了承托,在后續(xù)破壞階段,邊坡沿著第一滑裂面自下而上延伸出一條貫通坡體的裂縫。從圖8中邊坡滑落堆積的情況可以看出,邊坡在出現(xiàn)第一次滑落堆積之后仍然有后續(xù)滑落過程。因此FDEM-GIM能夠模擬邊坡在達到臨界失穩(wěn)狀態(tài)之后的開裂行為,邊坡最終滑落方量為貫通裂縫上部體積,而剛體極限平衡法僅能用來分析邊坡臨界破壞狀態(tài)。
本文基于重力增加法,采用連續(xù)-離散耦合分析方法對邊坡巖體進行失穩(wěn)破壞的全過程模擬,以紅石巖堰塞體邊坡工程為實例,將剛體極限平衡法得到的邊坡安全系數(shù)與FDEM-GIM得到的邊坡安全系數(shù)進行對比,驗證了基于重力增加法的連續(xù)-離散耦合分析方法運用于邊坡失穩(wěn)破壞問題的合理性,并模擬了邊坡巖體從未發(fā)生破壞到裂縫產(chǎn)生、擴展、破壞直至滑落堆積的全過程,結(jié)論如下:
(1)將FDEM-GIM計算結(jié)果與傳統(tǒng)剛體極限平衡法的結(jié)果進行對比,結(jié)果顯示2種方法得到的邊坡安全系數(shù)基本相同,最危險滑裂面位置與形狀基本吻合,驗證了基于重力增加法的連續(xù)-離散耦合分析方法進行邊坡失穩(wěn)破壞分析的合理性。
(2)邊坡加固措施必須依賴于潛在滑動面位置的確定,基于重力增加法的連續(xù)-離散耦合分析方法能自動搜索獲得邊坡臨界破壞滑面,克服傳統(tǒng)剛體極限平衡法需要預(yù)先人為假定滑裂面的缺陷。剛體極限平衡法僅能用于分析邊坡臨界失穩(wěn)狀態(tài),不能預(yù)測后續(xù)滑塊的形成。而基于重力增加法的連續(xù)-離散耦合分析方法能模擬邊坡破壞的全過程,并且可以顯示模擬后續(xù)滑坡體的形成,為邊坡防治處理提供可靠依據(jù)。