FEM/DEM法在典型柱狀危巖體破壞過程數(shù)值分析中的應(yīng)用

陳小婷，黃波林(防災(zāi)減災(zāi)湖北省重點實驗室/三峽大學(xué)，湖北 宜昌 443002)

目前開展斜坡變形破壞數(shù)值模擬的方法有基于有限單元的方法(FEM)、基于離散單元的方法(DEM)和連續(xù)離散單元方法(FEM/DEM耦合方法)。在處理復(fù)雜材料小變形機(jī)制方面，F(xiàn)EM有著明顯的優(yōu)勢[1]；而DEM在處理離散介質(zhì)的大變形問題上，具有非常重大的優(yōu)勢[2]。FEM/DEM耦合分析方法結(jié)合了連續(xù)和非連續(xù)數(shù)值方法的優(yōu)點，既可計算連續(xù)介質(zhì)的小變形，又可在裂隙間開展大變形，同時在采用斷裂力學(xué)的相關(guān)準(zhǔn)則后，可以對裂紋的產(chǎn)生和擴(kuò)展進(jìn)行模擬[3]。本文選用該方法對典型柱狀危巖體變形破壞進(jìn)行數(shù)值模擬分析。

1 FEM/DEM方法研究現(xiàn)狀

自1989年Munjiza[4]提出FEM/DEM耦合方法以來，這一方法被應(yīng)用在了廣泛的領(lǐng)域，其有效性和靈活性也得了極大的提高。Munjiza等、Owen等[5]、Xiang J等[6]和Latham等[7]推動了在FEM/DEM模型框架中的裂隙擴(kuò)展算法。在這個算法中，裂隙以有限單元網(wǎng)格的形式被植入。在有限元法中引入損傷力學(xué)和斷裂力學(xué)理論，當(dāng)單元網(wǎng)格受力超過強(qiáng)度準(zhǔn)則后，連續(xù)單元變成離散單元。離散單元的開裂發(fā)生在有限單元的邊界上，并采用離散元法處理開裂后塊體間的接觸。基于FEM/DEM方法，Morris等開發(fā)了Livermore distinct element code (LDEC)來進(jìn)行地質(zhì)體的數(shù)值模擬[8]，Rockfield公司則推出了ELFEN。國內(nèi)LI S H等[9]、魏懷鵬等[10]、馮春等[11]發(fā)展了連續(xù)—非連續(xù)離散單元法(CDEM)。Eberhardt等[12]利用EFLFEN對1991年Swiss的Randa巖質(zhì)崩滑體進(jìn)行了數(shù)值分析，研究了巖橋的破壞、崩滑體的運動及堆積過程。Vyazmensky等[13]使用基于FEM/DEM方法的ELFEN軟件研究了巖體開挖引起地面沉降的機(jī)理，指出該法可以加深對開挖過程中巖體的力學(xué)行為理解，并改進(jìn)沉降預(yù)測的方法，使用FEM/DEM方法的ELFEN軟件對露天采礦引發(fā)的邊坡破壞及尋找合適開挖方式進(jìn)行了研究。Mahabadi等[14]使用FEM/DEM方法構(gòu)建了Y-GUI程序，并對動力條件下的巴西圓盤試驗進(jìn)行數(shù)值模擬，模擬結(jié)果與試驗相吻合。Chang Xiaolin等[15]利用FEM/DEM方法和物理試驗同時開展了巴西圓盤試驗、單軸抗壓試驗、單邊切口梁試驗、平面滑移和撞擊試驗，模擬結(jié)果與試驗結(jié)果吻合程度較高。

本次研究選用了基于MunjizaY程序的Y-GUI程序開展箭穿洞危巖體的變形破壞過程研究。Y-GUI程序包括Mahabadi等建立的FEM/DEM前處理模塊和FEM/DEM方法的Y程序計算模塊。該模型包括摩爾庫倫模型、巖石斷裂模型和摩擦模型，可以用于模擬分析巖質(zhì)斜坡的變形破壞過程。

2 典型柱狀危巖體概況

箭穿洞危巖體位于三峽庫區(qū)重慶市巫山縣長江左岸，上距巫山縣城約15 km，位于巫峽著名的神女峰西側(cè)坡腳處，距下游神女溪口和下游青石居民聚集區(qū)1.3 km。在構(gòu)造上，箭穿洞斜坡位于神女峰背斜的核部。神女峰背斜核部出露的最老地層為三疊系大冶組三段(T1d3)，然后依次出露大冶組四段(T1d4)和三疊系嘉陵江組1-4段(T1j1-4)。

箭穿洞斜坡在地貌上呈現(xiàn)陡崖與中緩坡交替的臺階狀地貌，共有三級臺階。第一級臺階下發(fā)育箭穿洞危巖體，其臺階面為T1d4。由于巖層位于背斜核部，地層較平緩，產(chǎn)狀為260°∠6°。因此，箭穿洞危巖體斜坡總體為近水平層狀岸坡，巖性組成為高程155 m以下的中層狀含泥質(zhì)條帶灰?guī)r(T1d3)和高程155 m以上的灰色厚層礫屑、砂屑灰?guī)r(T1d4)。

箭穿洞危巖體的上游側(cè)邊界為沖溝，下游側(cè)邊界為一上陡下緩的大型裂隙，其上部產(chǎn)狀為324°∠65°，中段局部近直立轉(zhuǎn)緩傾，下部產(chǎn)狀324°∠45°。危巖體從臨空面陡崖至后緣共發(fā)現(xiàn)有卸荷裂縫，產(chǎn)狀為276°～260°∠75°～85°。后緣高程為278～305 m，基座高程為155 m，平均高差為135 m。平均橫寬約55 m，平均厚度約50 m，體積約36×104m3，主崩方向為260°(圖1)，其典型剖面圖見圖2。危巖體基座較破碎，但主柱體完整。

圖1 箭穿洞危巖體平面圖Fig.1 Geological engineering map of the Jianchuandong rockmass1—裂隙；2—巖體邊界；3—地層界線；4—巖層和結(jié)構(gòu)面產(chǎn)狀

泥質(zhì)條帶灰?guī)r位于箭穿洞危巖體的基座部位，相對厚層灰?guī)r較為軟弱?；鶐r體長期受壓，較破碎。2008年三峽庫區(qū)175 m蓄水后，泥質(zhì)條帶灰?guī)r處于庫水波動帶。水位周期性波動造成該基座巖體處于浸泡-風(fēng)干周期性作用下，造成基座巖體加速劣化[10]。基座的加速變形破壞，導(dǎo)致危巖體穩(wěn)定性降低，發(fā)生基座壓裂式崩塌。

圖2 箭穿洞危巖體工程地質(zhì)剖面圖Fig.2 Geo-engineering section map of the Jianchuandong dangerous rockmass

柱狀危巖體的變形破壞過程取決于危巖體的應(yīng)力應(yīng)變場。同時，由于危巖體入江會產(chǎn)生涌浪，入江方式不同會造成涌浪的大小差異。因此，需要對該柱狀危巖體的變形破壞過程進(jìn)行分析。

3 危巖體FEM/DEM模型構(gòu)建

野外調(diào)查表明，箭穿洞危巖體存在緩傾層面和較少的近垂向裂隙，除此以外沒有其他裂隙發(fā)育。危巖體整體可劃分為完整厚-塊狀巖體結(jié)構(gòu)，要模擬這樣一種巖體結(jié)構(gòu)的變形破壞過程，需要能實現(xiàn)大變形和形成新裂隙的數(shù)值工具。本文采用FEM/DEM模型的Y-GUI程序，對箭穿洞危巖體進(jìn)行數(shù)值模型構(gòu)建、參數(shù)設(shè)置和破壞全過程計算。

圖3為箭穿洞危巖體的FEM/DEM模型，模型長384 m，高222 m。模型總單元數(shù)為700個，其中危巖體及基座巖體單元數(shù)為250個，它們是可離散單元；斜坡巖體單元數(shù)為450個，它們是不可離散單元。模型中巖體和結(jié)構(gòu)面所用參數(shù)見表1，模型采用彈性Mohr-Coulomb準(zhǔn)則和拉張破壞準(zhǔn)則。

4 危巖體破壞過程分析

從初期的危巖體受力情況來看，基座和基座附近巖體明顯有應(yīng)力集中現(xiàn)象(圖4)。危巖體基座的拉應(yīng)力主要來源于壓力，即壓致拉裂效應(yīng)。從此時基座巖體的運動矢量圖來看(圖5)，基座巖體受壓向兩側(cè)和向下運動明顯。

圖3 箭穿洞危巖體的FEM/DEM模型Fig.3 FEM/DEM model of the Jianchuandong dangerous rockmass

表1 模型中巖體和結(jié)構(gòu)面參數(shù)Table 1 Rockmass and structure face parameters for the numerical model

圖4 基座巖體的應(yīng)力集中現(xiàn)象Fig.4 Stress concentration in the basic rockmass注：紅色為應(yīng)力陰影

圖5 基座巖體的運動方向Fig.5 Moving direction of the basic rockmass注：紅色為向x-方向，藍(lán)色向x+方向

裂隙最開始新生在基座泥質(zhì)條帶灰?guī)r以及與泥質(zhì)條帶灰?guī)r接觸的灰?guī)r處。基座泥質(zhì)條帶灰?guī)r區(qū)由于壓力作用裂隙雖然出現(xiàn)破壞但沒有分離，因此顯示為局部出現(xiàn)空洞，而不是裂隙(紅色區(qū)域)，基座裂隙最先出現(xiàn)在兩端(圖6)。上部危巖體的裂隙大多是從底部開始往上延伸，其形成機(jī)制是典型的壓致拉裂，且裂隙最開始出現(xiàn)的位置靠近內(nèi)側(cè)(圖7)。

圖6 基座巖體最初破壞位置圖Fig.6 Original failure location map of the basic rockmass注：黑色圈為破壞位置，藍(lán)色為壓破壞區(qū)域，紅色為拉張破壞區(qū)域

圖7 危巖體初始破壞圖Fig.7 Original failure of the main rockmass注：紅色線為巖體出現(xiàn)張開的裂隙

當(dāng)基座巖體和危巖體巖體開始出現(xiàn)新裂隙后，裂隙網(wǎng)絡(luò)迅速發(fā)展，最終將危巖體分成3個相互關(guān)聯(lián)的巖柱。裂隙網(wǎng)絡(luò)的發(fā)展和基座的破壞導(dǎo)致了單個巖柱出現(xiàn)破壞。破壞從臨空面開始，基座開始垮塌，然后上部的巖柱或滑動或垮落或傾倒。當(dāng)臨空側(cè)的巖柱倒塌開始后，剩余的巖柱也加速運動，最后整個危巖體解體失穩(wěn)。整個破壞過程，有滑移、墜落、傾倒等狀態(tài)，與坐落過程類似。圖8展示了裂隙網(wǎng)絡(luò)的發(fā)展、危巖體的失穩(wěn)過程。

FEM/DEM方法預(yù)演了箭穿洞危巖體可能的變形破壞過程。從上述分析來看，箭穿洞危巖體今后裂隙的主要形成機(jī)制為壓致拉裂，這與定性的變形失穩(wěn)模式分析一致。壓致拉裂主要是基座巖體被壓縮導(dǎo)致的拉張破壞。破壞的基座巖體開始向臨空面運動，直接導(dǎo)致了危巖體開始垮落。值得注意的是，基座巖體雖然壓裂了，并未全部擠出，只有部分基座巖體離開了原位。圖8顯示這部分運動離開的巖體約占全部基座巖體的50%。

5 FEM/DEM討論

箭穿洞危巖體案例分析充分展示了FEM/DEM方法的獨特之處，它通過明確地考慮斷裂和碎裂過程來模擬從連續(xù)到不連續(xù)的轉(zhuǎn)變過程。小變形和應(yīng)力應(yīng)變由有限元進(jìn)行分析；顯式時間推進(jìn)方案被用來整合離散單元牛頓運動方程，因此可以完全動態(tài)地模擬破壞過程。FEM/DEM方法對幫助分析大小變形混合的斜坡破壞與運動十分有利，適用于巖質(zhì)危巖體或滑坡的破壞分析。

圖8 箭穿洞危巖體變形破壞全過程數(shù)值模擬圖(S為時步)Fig.8 Failure process map of the Jianchuandong dangerous rockmass

同時，需要注意的是，當(dāng)前的FEM/DEM分析框架僅考慮了自重作用，并沒有考慮水及振動等因素。對水因素的疊加，在FEM/DEM框架中可采用增加外部壓力和增減重度的方式進(jìn)行；對振動因素的增加，可采用脈沖式加速度來進(jìn)行。增加考慮水和振動工況后，將使得FEM/DEM的適用范圍更為廣泛，其分析結(jié)果將更接近工程實際。

6 結(jié)論與建議

(1)柱狀危巖體的基座處于應(yīng)力集中狀態(tài)，既有拉應(yīng)力集中區(qū)也有壓應(yīng)力集中區(qū)，應(yīng)力集中造成柱狀危巖體的基座附近邊界巖體易破壞。

(2)柱狀危巖體裂隙的延展大部分是從下而上發(fā)展的，也有一部分是從上向下發(fā)展的，裂隙的延展導(dǎo)致了危巖體的解體和最終的破壞。

(3)FEM/DEM方法預(yù)演了巖體壓致拉裂的裂隙網(wǎng)絡(luò)形成和破碎巖體運動過程。數(shù)值模擬驗證了箭穿洞危巖體會發(fā)生基座壓裂型崩塌這一定性分析?；谠撐r體未來發(fā)展的認(rèn)識，建議對該柱狀危巖體基座進(jìn)行補強(qiáng)防護(hù)。

(4)FEM/DEM耦合了有限元法和離散元法的優(yōu)點，有利于分析巖體破裂及運動過程，是未來數(shù)值模擬的發(fā)展方向之一。