爆破動力荷載下邊坡穩定性的時程分析

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

爆破動力荷載下邊坡穩定性的時程分析

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

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

動力穩定 安全系數 時程分析 LS-DYNA 數值模擬

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

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

1 有限元法+極限平衡法

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

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

(1)

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

2 工程實例

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

2.1 數值建模

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

圖2 臺階平面模型

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

(2)

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

表1 炸藥及其狀態方程參數Table 1 Explosive and its parameters in state equation

2.2 網格劃分及加載約束

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

圖3 臺階二維和三維網格劃分

3 結果分析

3.1 靜力求解結果

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

圖4 塑性區域貫通

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

圖5 邊坡最危險滑面

表2 極限平衡法和靜態有限元法計算結果Table 2 The results of limit equilibrium method and FEM

3.2 動力有限元求解結果

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

(2)邊坡位移響應分析。由計算結果得知，坡頂位置出現較大的位移波動，波動范圍為-0.35～0.25 cm，坡腳及坡面單元基本集中在-0.1～0.1 cm。在測點的Y軸、Z軸位移監測中，最大值均出現在坡頂及其附近單元。邊坡內滑面單元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。雖然最值發生的時間不一致，但是最大值均發生在靠近坡腳的地方。

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

圖7 測點Y軸振動速度

3.3 安全系數時程分析

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

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

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

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

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

4 總 結

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

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

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

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

[1] 劉立平,雷尊宇,周富春.地震邊坡穩定分析方法綜述[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] 鄭穎人,葉海林,黃潤秋,等.邊坡地震穩定性分析探討[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] 劉漢龍,費 康,高玉峰.邊坡地震穩定性時程分析方法[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] 許名標,彭德紅.邊坡爆破振動測試及響應規律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] 許紅濤,盧文波,周創兵,等.基于時程分析的巖質高邊坡開挖爆破動力穩定性計算方法[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] 陳祖煜.土質邊披穩定分析——原理·方法·程序[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] 鄭穎人,趙尚毅.有限元強度折減法在土坡與巖坡中的應用[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] 趙尚毅,鄭穎人,張玉芳.有限元強度折減法中邊坡失穩的判據探討[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] 馬建勛,賴志生,蔡慶娥,等.基于強度折減法的邊坡穩定性三維有限元分析[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] 鄭穎人,趙尚毅.巖土工程極限分析有限元法及其應用[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