爆破動力荷載下邊坡穩(wěn)定性的時程分析

陳鵬輝 樓曉明,2 周文海(.福州大學紫金礦業(yè)學院，福建 福州 3506；2.福州大學爆炸技術(shù)研究所，福建 福州 3506)

爆破動力荷載下邊坡穩(wěn)定性的時程分析

陳鵬輝1樓曉明1,2周文海1(1.福州大學紫金礦業(yè)學院，福建 福州 350116；2.福州大學爆炸技術(shù)研究所，福建 福州 350116)

針對露天礦開采過程中，爆破開挖作業(yè)對最終邊坡穩(wěn)定的影響問題，采用動力有限元法分析爆破過程中邊坡的穩(wěn)定性隨時間變化情況。首先建立礦山邊坡的二維有限元模型，利用有限元強度折減法求解，得到邊坡的最危險滑動面及折減安全系數(shù)；再根據(jù)得到的最危險滑動面，重新建立邊坡和炸藥的三維實體模型，布置爆破孔網(wǎng)參數(shù)并設(shè)置各個炮孔的起爆時間，通過LS-DYNA大型動力有限元程序，求解邊坡在炸藥爆炸后30 s時間內(nèi)的動力響應(yīng)過程。結(jié)合極限平衡法和動力有限元法計算結(jié)果，計算邊坡在同時起爆方式下每一時刻的安全系數(shù)。提出應(yīng)根據(jù)邊坡時程安全系數(shù)及邊坡滑面單元是否同時發(fā)生整體破壞等條件綜合評價邊坡穩(wěn)定性。數(shù)值模擬結(jié)果能夠動態(tài)顯示爆破過程中邊坡的應(yīng)力、應(yīng)變狀態(tài)，計算結(jié)果與極限平衡法計算結(jié)果接近，可以為露天礦爆破設(shè)計提供參考。

動力穩(wěn)定 安全系數(shù) 時程分析 LS-DYNA 數(shù)值模擬

爆破荷載下巖質(zhì)邊坡的動力響應(yīng)及控制是礦山開發(fā)和工程建設(shè)中面臨的一個亟待解決的課題。邊坡穩(wěn)定是露天礦安全生產(chǎn)的前提，隨著露天礦不斷進入深凹開采，爆破地震效應(yīng)對邊坡穩(wěn)定的影響問題日益突出。炸藥爆炸過程是一個短暫且復雜的過程，一般認為，爆破動力荷載對邊坡穩(wěn)定性的影響主要是因為反復的爆破振動引起邊坡穩(wěn)定性減弱，巖體結(jié)構(gòu)抗剪強度指標降低而產(chǎn)生的衰減失穩(wěn)以及爆破振動自身產(chǎn)生的地震慣性力導致坡體整體的下滑力增大，最終引起整個坡體的慣性失穩(wěn)[1]。

衰減失穩(wěn)通常采用變形破壞分析法和流動破壞分析法；而慣性失穩(wěn)主要的分析方法有Newmark滑塊分析法、擬靜力法、Makidisi Seed法以及有限元法等?，F(xiàn)階段巖質(zhì)邊坡的動力穩(wěn)定數(shù)值模擬中，爆破動力荷載的加載方式大多采用以下幾種方法：采用實測或人工合成的爆破振動速度、加速度時程曲線等轉(zhuǎn)換為爆破動力荷載[2-5]，或者采用理論計算的等效爆轟壓力施加在孔壁或炮孔連心線所在的平面上等方法進行動力響應(yīng)時程分析[6-7]。

1 有限元法+極限平衡法

在傳統(tǒng)的邊坡穩(wěn)定性分析中，極限平衡法和有限元法一般都是分開分析的。極限平衡法因其原理簡單、計算簡便，在工程上具有廣泛的應(yīng)用，積累了豐富的實踐經(jīng)驗。而有限元法不需要預(yù)先假設(shè)滑動面，同時考慮了巖體應(yīng)力-應(yīng)變關(guān)系，能夠真實反映巖體的實際情況，更適合分析大型、復雜的地質(zhì)邊坡。利用有限元法得到邊坡的最危險滑面，通過有限元軟件對滑面進行單元劃分。將動力分析后得到的滑動面單元應(yīng)力結(jié)果σxi、σyi、σxyi，結(jié)合極限平衡法計算邊坡的安全系數(shù)[8-9]，如圖1及公式(1)。

圖1 有限元法+極限平衡法

(1)

式中，σi為作用在弧段上的法向應(yīng)力(有限元計算以拉為正，壓為負)；τi為作用在弧段上的剪應(yīng)力；σxi，σyi，σxyi為有限元計算得到的應(yīng)力;θi為弧段Δli與x軸的夾角;τfi為抗剪強度;li為弧段長度。

2 工程實例

在礦山正常生產(chǎn)作業(yè)過程中，爆破作業(yè)從設(shè)計到施工是一個系統(tǒng)工程。在利用有限元軟件進行建模分析時，受到模型材料參數(shù)、地形地質(zhì)條件等因素的制約。特別是在建立三維模型的時候，劃分單元數(shù)量會急劇增加，對計算機的硬件配置提出更高的要求。以青海威斯特銅礦為基礎(chǔ)，建立一個臺階的計算模型。巖石密度2.5g/cm3，彈性模量22.5GPa，泊松比0.22，黏聚力0.22MPa，內(nèi)摩擦角27.4°。首先采用有限元強度折減法[10-13]計算二維平面應(yīng)變狀態(tài)下邊坡的最危險滑動面，再按照現(xiàn)場爆破參數(shù)建立三維實體模型及炸藥模型，運用LS-DYNA程序求解爆破過程中邊坡滑面上單元的應(yīng)力應(yīng)變狀態(tài)，分析邊坡在爆破動力作用下穩(wěn)定性隨時間變化情況。

2.1 數(shù)值建模

模型采用cm-g-us單位制，臺階高度24m(并段后)，坡面角70°，臺階尺寸如圖2。設(shè)計孔網(wǎng)參數(shù)為孔距6m，排距5m，孔徑150mm，孔深14m，裝藥長度9m，填塞長度5m，采用耦合裝藥方式，單孔藥量為192kg，共布置3排10個炮孔。臺階模型尺寸根據(jù)鄭穎人等[10]的研究經(jīng)驗，坡腳至左邊界延伸坡高的1.5倍，坡頂線至右邊界延伸坡高的至少2.5倍，上下邊界距離至少為2倍的坡高建模。以臺階橫向為X軸，豎向為Y軸，縱向為Z軸，三維狀態(tài)下Z軸延伸50 m，即臺階爆破區(qū)域長50 m，寬36 m。臺階巖體采用塑性隨動強化模型，炸藥采用2#巖石乳化炸藥，炸藥密度0.95 g/cm3，爆速3 600 m/s，填塞材料和臺階巖體材料一致。

圖2 臺階平面模型

在爆炸場的數(shù)值模擬中，采用美國勞倫斯·利弗摩爾實驗室的Lee等提出JWL狀態(tài)方程[14]：

(2)

式中,p為狀態(tài)方程決定的壓力；A，B，R1，R2，w為描述JWL方程獨立的5個物理常數(shù);V是相對體積；E0為初始內(nèi)能。炸藥及狀態(tài)方程具體取值如表1所示。

表1 炸藥及其狀態(tài)方程參數(shù)Table 1 Explosive and its parameters in state equation

2.2 網(wǎng)格劃分及加載約束

如圖3，模型網(wǎng)格二維狀態(tài)下嚴格按照四邊形，三維狀態(tài)六面體進行網(wǎng)格劃分。二維分析模型總共劃分4 521個單元，三維模型總共劃分170 958個單元。二維模型底部和左右邊界施加位移約束，邊坡面積坡頂為自由邊界；三維模型除臺階面和臺階坡面為自由邊界以外，其余邊界均施為透射邊界(無反射邊界)條件。二維模型施加重力荷載，三維模型利用動力松弛加載預(yù)應(yīng)力，在動力求解過程中保持重力加載。

圖3 臺階二維和三維網(wǎng)格劃分

3 結(jié)果分析

3.1 靜力求解結(jié)果

二維靜力計算時，采用二分法不斷地折減內(nèi)摩擦系數(shù)和黏聚力，根據(jù)臺階坡面廣義塑性應(yīng)變或者等效塑性應(yīng)變從坡腳到坡頂貫通、計算不收斂、邊坡單元或節(jié)點位移發(fā)生突變等邊坡失穩(wěn)判據(jù)[11]判斷邊坡是否失穩(wěn)。得到在折減系數(shù)取2.750時，邊坡塑性區(qū)域已經(jīng)貫通，但邊坡并沒有破壞，處于臨界破壞狀態(tài)。當折減系數(shù)繼續(xù)增大時，邊坡塑性應(yīng)變區(qū)域貫通(見圖4)，計算機計算不收斂，節(jié)點位移發(fā)生突變，同時發(fā)生邊坡失穩(wěn)，發(fā)生破壞。

圖4 塑性區(qū)域貫通

利用有限元強度折減法得到靜力狀態(tài)下邊坡的最危險滑面(見圖5)，結(jié)合傳統(tǒng)的極限平衡法，計算邊坡靜態(tài)下的安全系數(shù)，具體結(jié)果見表2。傳統(tǒng)的極限平衡法計算結(jié)果和有限元計算結(jié)果比較接近，相對誤差在3%左右。

圖5 邊坡最危險滑面

表2 極限平衡法和靜態(tài)有限元法計算結(jié)果Table 2 The results of limit equilibrium method and FEM

3.2 動力有限元求解結(jié)果

為便于分析，本次模擬第1排3個炮孔同時起爆，同段起爆藥量為576 kg，計算時間分別為1 s、5 s、30 s 3種時長。選取邊坡坡腳(H2)、坡面(H5738)、坡頂(H5381)及上臺階面單元(H28944)進行分析，各點距離爆心的水平距離分別為10.80、15.53、19.80、52.58 m，如圖6所示。

圖6 單元位置

(1)邊坡應(yīng)力響應(yīng)分析。在炸藥爆炸過程中，對炮孔中間位置每隔10 ms進行應(yīng)力云圖切片觀察。炸藥起爆后，孔壁壓力迅速上升，炮孔壁上單元的最大應(yīng)力達49 MPa，達到峰值后快速衰減，在60 ms左右下降到0.2 MPa。而邊坡表面的測點有效應(yīng)力明顯小于炮孔周圍的單元，坡腳單元有效應(yīng)力大于其他測點，其最大值達1 MPa。另外，邊坡測點的最大主應(yīng)力，最小主應(yīng)力以及最大剪應(yīng)力都在坡腳附近單元出現(xiàn)應(yīng)力集中現(xiàn)象，但出現(xiàn)的時間都很短暫，在爆炸結(jié)束后，基本趨于穩(wěn)定值。

(2)邊坡位移響應(yīng)分析。由計算結(jié)果得知，坡頂位置出現(xiàn)較大的位移波動，波動范圍為-0.35～0.25 cm，坡腳及坡面單元基本集中在-0.1～0.1 cm。在測點的Y軸、Z軸位移監(jiān)測中，最大值均出現(xiàn)在坡頂及其附近單元。邊坡內(nèi)滑面單元X軸位移最大值為5.4 s時的0.019 8 cm，最小值為0.3 s時的-0.029 8 cm；Y軸位移最大值為6.3 s時的0.031 4 cm，最小值為0.6 s時的-0.079 cm；Z軸位移最大值為0.3 s時的0.000 916 cm，最小值為0.9 s時的-0.002 14 cm。雖然最值發(fā)生的時間不一致，但是最大值均發(fā)生在靠近坡腳的地方。

(3)邊坡速度響應(yīng)分析。由于礦山邊坡是一個臨時邊坡，并且邊坡的狀態(tài)隨時可能發(fā)生改變，目前尚無一個能夠判斷邊坡穩(wěn)定性的安全振速統(tǒng)一標準。從圖7時程曲線中可知，坡頂單元的振動速度明顯高于其他部位單元，出現(xiàn)速度放大效應(yīng)。坡頂振動速度最大值達到17 cm/s，坡腳振動速度最大值為8 cm/s。測點合速度中，坡頂速度依然大于其他部位，最大值為0.9 s時的18.9 cm/s。

圖7 測點Y軸振動速度

3.3 安全系數(shù)時程分析

根據(jù)程序計算結(jié)果，在后處理軟件LS-PREPOST中提取出滑面單元在爆炸過程中的X、Y、XY軸應(yīng)力，利用式(1)進行時程安全系數(shù)計算(見圖8)。計算結(jié)果如表3所示。

圖8 邊坡三維有限元計算模型

表3 動力有限元計算結(jié)果Table 3 The results of dynamic FEM

在爆破動力作用下，邊坡內(nèi)部應(yīng)力狀態(tài)時刻在發(fā)生變化，其時程安全系數(shù)在靜態(tài)安全系數(shù)一定范圍內(nèi)上下波動，隨著時間的增加，時程安全系數(shù)又趨于穩(wěn)定。由于3種時長的計算結(jié)果在提取時間的間隔上有差異，結(jié)果也存在一定差異。從側(cè)面反映了炸藥爆炸動力荷載對邊坡巖體的瞬時作用特性。爆炸過程中，時程安全系數(shù)在個別時刻由于受到爆炸沖擊波的作用，個別單元應(yīng)力發(fā)生突變，導致整體安全系數(shù)在瞬間會突然增大或減小，最大值達到12.434，是初始值的4.5倍，最小值為0.907，比初始值降低66.68%，但很快均返回到正常范圍。劉漢龍等[2]采用最小平均安全系數(shù)，即靜力下的安全系數(shù)Fs0與最小動力安全系數(shù)Fs min差值的0.65倍作為安全系數(shù)的平均振幅來反應(yīng)安全系數(shù)因振動作用而偏離的幅度，計算得最小平均安全系數(shù)

雖然滑面單元在個別時刻應(yīng)力發(fā)生突變，但滑面單元的最大剪應(yīng)力并沒有達到破壞值。模擬結(jié)果顯示，邊坡的有效塑性應(yīng)變只發(fā)生在坡腳的部分單元，邊坡的坡面及內(nèi)部并沒有發(fā)生失效破壞，邊坡整體還是穩(wěn)定的。當邊坡的初始安全系數(shù)比較低時，隨著這種單元應(yīng)力突變次數(shù)的增加，邊坡內(nèi)部就很有可能因損傷累積效應(yīng)發(fā)生失穩(wěn)破壞。單元應(yīng)力狀態(tài)的突變是邊坡失穩(wěn)的潛在危險因素。因此，采用時程安全系數(shù)中的最小值或者最小平均安全系數(shù)作為邊坡穩(wěn)定性的判據(jù)是不夠準確的。應(yīng)該同時結(jié)合邊坡滑面單元的塑性應(yīng)變及位移情況，最大剪應(yīng)力是否達到破壞等條件綜合評價邊坡的穩(wěn)定性。

4 總 結(jié)

(1)利用動力有限元法，能直觀顯示在爆破動力荷載下邊坡安全系數(shù)隨時間的變化情況。同時，在計算精度上與傳統(tǒng)極限平衡法相當接近，為邊坡穩(wěn)定性設(shè)計及爆破參數(shù)優(yōu)化提供參考。

(2)通過LS-DYNA程序可以直接模擬三維狀態(tài)下高能炸藥爆炸的整個動力作用過程。根據(jù)工程實際，建立爆破實體模型，設(shè)置任意起爆時間和起爆順序，簡化了爆破動力的加載方式和加載時間。同時，結(jié)合有限元強度折減法分析邊坡在爆破動力下的時程響應(yīng)規(guī)律，簡化了最危險滑動面搜索的工作量。

圖9 極限平衡法與有限元法計算結(jié)果比較(T=30 s)

(3)為便于模擬分析，本研究沒有考慮節(jié)理裂隙、孔隙水等因素對巖質(zhì)邊坡的影響。另外，由于炸藥爆炸過程時間及其短暫，遠小于現(xiàn)實的延期時間和爆破后沖擊波對邊坡的作用和反應(yīng)時間，如果要完整計算長時間下邊坡的動力響應(yīng)，可以根據(jù)實際情況在不影響實驗結(jié)果的前提下對模型參數(shù)進一步優(yōu)化，以減少計算機計算時間和對計算機內(nèi)存空間的占用。

[1] 劉立平,雷尊宇,周富春.地震邊坡穩(wěn)定分析方法綜述[J].重慶交通學院學報,2001,20(3):83-88. Liu Liping,Lei Zunyu,Zhou Fuchun.The evaluation of seismic slope stability analysis methods[J].Journal of Chongqing Jiaotong University,2001,20(3):83-88.

[2] 鄭穎人,葉海林,黃潤秋,等.邊坡地震穩(wěn)定性分析探討[J].地震工程與工程振動,2010(2):173-180. Zheng Yingren,Ye Hailin,Huang Runqiu,et al.Study on the seismic stability analysis of a slope[J].Journal of Earthquake Engineering and Engineering Vibration,2010(2):173-180.

[3] 劉漢龍,費 康,高玉峰.邊坡地震穩(wěn)定性時程分析方法[J].巖土力學,2003(4):553-556. Liu Hanlong,Fei Kang,Gao Yufeng.Time history analysis method of slope seismic stability [J].Rock and Soil Mechanics,2003(4):553-556.

[4] 許名標,彭德紅.邊坡爆破振動測試及響應(yīng)規(guī)律ANSYS時程分析[J].巖石力學與工程學報,2012(S1):2629-2635. Xu Mingbiao,Peng Dehong.Blasting vibration tests and analysis time-history analyses of slope responses[J].Chinese Journal of Rock Mechanics and Engineering,2012(S1):2629-2635.

[5] 陳 明,盧文波,舒大強,等.爆破振動作用下邊坡極限平衡分析的等效加速度計算方法[J].巖石力學與工程學報,2009(4):784-790.

Chen Ming,Lu Wenbo,Shu Daqiang,et al.Calculation method of equivalent acceleration for limit equilibrium analysis of slope under blasting vibration[J].Chinese Journal of Rock Mechanics and Engineering,2009(4)：784-790.

[6] 趙婉婷,盧文波,楊建華,等.深孔臺階爆破振動模擬中的等效荷載施加邊界比較[J].爆破,2012(2):10-14. Zhao Wanting,Lu Wenbo,Yang Jianhua,et al.Comparison of equivalent load in boundaries in deep-hole bench blasting vibration simulation[J].Blasting,2012(2):10-14.

[7] 許紅濤,盧文波,周創(chuàng)兵,等.基于時程分析的巖質(zhì)高邊坡開挖爆破動力穩(wěn)定性計算方法[J].巖石力學與工程學報,2006(11):2213-2219. Xu Hongtao,Lu Wenbo,Zhou Chuangbing,et al.Time history analysis method for evaluating dynamic stability of high rock slope under excavation blasting[J].Chinese Journal of Rock Mechanics and Engineering,2006(11)：2213-2219.

[8] 陳祖煜.土質(zhì)邊披穩(wěn)定分析——原理·方法·程序[M].北京:中國水利水電出版社,2003:255-256. Chen Zuyu.Soil Slope Stability Analysis：Theory Methods and Programs[M].Beijing:China Water Power Press,2003:255-256.

[9] 王成華.土力學[M].武漢:華中科技大學出版社,2010:346-348. Wang Chenghua.Soil Mechanics[M].Wuhan:Huazhong University of Science & Technology Press,2010:346-348.

[10] 鄭穎人,趙尚毅.有限元強度折減法在土坡與巖坡中的應(yīng)用[J].巖石力學與工程學報,2004,19:3381-3388. Zheng Yingren,Zhao Shangyi.Application of strength reduction fem in soil and rock slope[J].Chinese Journal of Rock Mechanics and Engineering,2004,19:3381-3388.

[11] 趙尚毅,鄭穎人,張玉芳.有限元強度折減法中邊坡失穩(wěn)的判據(jù)探討[J].巖土力學,2005(2)：332-336. Zhao Shangyi,Zheng Yingren,Zhang Yufang.Study on slope failure criterion in strength reduction finite element method[J].Rock and Soil Mechanics,2005(2):332-336.

[12] 馬建勛,賴志生,蔡慶娥,等.基于強度折減法的邊坡穩(wěn)定性三維有限元分析[J].巖石力學與工程學報,2004,16:2690-2693. Ma Jianxun,Lai Zhisheng,Cai Qing′e,et al.3D FEM analysis of slope stability based on strength reduction method[J].Chinese Journal of Rock Mechanics and Engineering,2004,16:2690-2693.

[13] 鄭穎人,趙尚毅.巖土工程極限分析有限元法及其應(yīng)用[J].土木工程學報,2005(1):91-98. Zheng Yingren,Zhao Shangyi.Limit state finite element method for geotechnical engineering analysis and its applications[J].China Civil Engineering Journal,2005(1):91-98

[14] Livemore Software Technology Corporation.LS-DYNA Version 971 Keyword User′s Manual[M].California:Livemore Software Technology Corporation,2007.

(責任編輯 徐志宏)

Time History Analysis of Slope Stability under the Blasting Loading

Chen Penghui1Lou Xiaoming1,2Zhou Wenhai1
(1.CollegeofZijinMining,FuzhouUniversity,Fuzhou350116,China;2.ResearchInstituteofExplosionTechnology,FuzhouUniversity,Fuzhou350116,China)

For the significant impact of blasting on slope stability during open-pit mining process,the dynamic finite element method is used to analyze the variation of slope stability with the time in blasting.Firstly,the finite element model of mine slope is established,then the finite element strength reduction method is adopted to get the most dangerous sliding surface and safety coefficient of the slope.After that,based on the most dangerous sliding surface obtained,a three-dimensional entity model of slope and explosives was created,the blasting parameters are arranged and the initiation time for each hole is set up.Next,this model is calculated after the initiation.The dynamic response process after slope explosive for 30 s is solved by LS-DYNA large dynamic finite element program.Combining with the limit equilibrium method and dynamic finite element method,the safety factors of slope at every moment under unified detonation are calculated.Slope stability should be comprehensively evaluated by the time-history safety coefficient and whether the elements of the sliding surface are destroyed at the same time.The simulation results can dynamically display the security status of slope blasting process,and the result is very close to that by the limit equilibrium method.This research could make a reference for the blasting design.

Dynamic stability,Safety coefficient,Time history analysis,LS-DYNA,Numerical simulation

2014-10-29

福建省自然科學基金項目(編號：2001J01310)。

陳鵬輝(1990—)，男，碩士研究生。通訊作者 樓曉明(1972—)，男，博士，副教授。

TD235

A

1001-1250(2015)-11-029-05