基于非連續變形分析方法的爆破荷載下不同類別巖體的動力響應

江 巍，程圣國，馮 強

(三峽大學三峽庫區地質災害教育部重點實驗室，湖北 宜昌 443002)

爆破是工程巖體開挖的常見手段，20世紀60年代以來爆破技術已經得到了長足的發展，但爆破帶來的不利影響問題(如洞室安全[1]、邊坡穩定[2]、近區民居破壞[3]等)也愈發明顯。隨著爆破理論和計算機技術的發展，數值模擬已經成為分析爆破產生不利影響的重要手段，動力有限元[4]、離散元[5]、有限差分法[2]等多種計算方法被引入該領域并取得了一定的研究成果。

采用連續數值計算方法分析爆炸荷載的效應時將工程巖體視為均一介質，這與巖體中廣泛存在裂隙和節理的實際情況并不相符[6]，而且對于爆炸荷載作用過程中，巖體結構面的張開拉裂、塊體飛散等動態過程，采用不連續數值計算方法更為合適。非連續變形分析(discontinuous deformation analysis，DDA)方法是石根華先生提出的用于統一分析塊體系統靜力和動力行為的不連續數值方法[7]，關于其有效性和正確性Mac Laughlin等[8]已經對1986—2006年的相關工作做了全面的綜述。該方法在巖土工程領域已經得到了廣泛的研究與應用，但是在模擬工程巖體中的爆破問題方面則剛剛起步，如朱傳云等[9]進行了臺階爆破的仿真模擬，劉紅巖等[10]進行了爆破漏斗形成過程模擬等。

對于工程巖體爆破振動引起的不利影響，主要通過考察工程區巖體的動力響應進行衡量。工程巖體根據巖體質量常劃分為不同類別，本文以非連續變形分析方法為工具，對不同類別的巖體在爆破荷載作用下的動力響應進行數值研究，分析巖體振動的相關參數(起振時間、峰值振速)隨距離的變化特征，為爆破的控制提供技術依據。

1 爆破荷載的等效

用數值方法模擬爆破后巖體的動力響應，首要前提就是確定爆破荷載。目前爆破荷載模型可以粗略劃分為經驗公式法的爆破荷載模型和考慮爆生氣體狀態的爆破荷載模型兩類。

經驗公式法建立的爆破荷載模型將作用在爆孔壁上的爆破荷載壓力簡化為隨時間變化的函數，基本形式為

式中：Pb為壓力峰值;f(t)通常為指數型的時間滯后函數。考慮藥包在爆孔中爆炸時，爆孔周圍巖體在瞬間高壓作用下變成似流體介質，形成粉碎區，該類方法直接將壓力作用于爆孔壁上，則如何定義粉碎區的巖石力學特性將成為一個難題。

考慮爆生氣體狀態的爆破荷載模型將爆孔內的爆生氣體壓力P 與爆腔體積V 對應起來：

式中的爆腔體積V 包含了爆孔體積和孔壁裂隙擴展產生的體積。李寧等[11-12]均對這一模型做了優秀的研究工作，并取得了豐富的成果，但該方法要求將爆孔周圍網格劃分得十分細小，在面對大型的計算模型時，會產生過多的計算單元和龐大的計算量。

鑒于上述兩類方法的局限性，尋求一種適當的爆破荷載簡化等效方法是有必要的，本文采用徐全軍等[13]提出的粉碎區邊緣施加爆炸壓力方法，該方法通過分析水孔法測量得到的靠近炮孔巖體中的壓力波形，獲取粉碎區邊緣的壓力曲線，其表達式如下：

式中：Pe為作用在粉碎區邊緣的壓力峰值，MPa;α為衰減系數，與壓力持續時間有關。對于完全耦合裝藥，作用在粉碎區邊緣的壓力峰值為

式中：γ為絕熱指數;ρ0為炸藥初始密度，kg/m3;Qv為單位質量爆熱，J/kg;D為爆速，m/s;d0為爆孔直徑，m;de為粉碎區直徑，m。爆炸荷載的持續時間為

式中：R1為離爆源的距離，m;C為爆生氣體膨脹后的音速，取370m/s。以銨油炸藥為對象[13]，可得Pe=5.87MPa，T=0.00216 s，α=4500，粉碎區邊緣壓力曲線如圖1 所示。

圖1 粉碎區邊緣壓力曲線

該方法的基本思想為：當藥包爆炸后，以爆孔為中心周圍的巖體按照由近至遠形成粉碎區、破碎區和振動區，將爆破荷載直接施加在粉碎區的邊緣，則可以避開粉碎區的巖石力學特性定義難題。基于此方法徐全軍等[13]采用非線性有限元程序對爆破振動效應進行了分析，取得了一定的成果，但在計算中為使用有限元技術將破碎區和振動區的巖體均近似為均勻連續介質。本文以非連續變形分析方法為手段，視破碎區和振動區的巖體由線彈性體的巖塊和節理組成，破碎區巖體中除巖塊發生彈性變形外更主要的是節理新生、擴張和延伸等導致巖體發生破碎行為，振動區的巖體中節理受力過小變形細微，兩側巖塊之間力的傳遞良好，其力學行為表現接近連續介質。

2 非連續變形分析基本理論與爆破荷載的施加

2.1 非連續變形分析方法基本理論簡介

非連續變形分析方法中的基本單元是塊體。塊體的位移和變形是由分時步的位移和變形累加而成的，在每一時步所有單元的位移都是小量。考慮每個塊體為常應力和常應變，塊體中一點(x，y)的位移(u，v)可用位移不變量表示為

式中：(u0，v0)為塊體重心(x0，y0)的剛體平動;r0為塊體繞點(x0，y0)的轉動角;εx，εy，γxy為該塊體的法向和切向應變。塊體單元任意點的變形均可由上述位移不變量乘形函數矩陣得到。

塊體系統方程組的建立根據系統勢能最小化原理獲得，系統的總勢能包括塊體單元的應變能、初始應力的勢能、點荷載和線荷載的勢能、體荷載勢能、錨桿作用力的勢能、慣性力勢能以及塊體相互接觸力的勢能。塊體系統總勢能Π的一般表達形式

式中：D，K，F 分別為整個塊體系統的位移不變量、總剛度矩陣和總荷載向量。 Π取極值的條件為其對D的每一分量求偏導均為零，于是可得支配方程：

施加邊界條件和塊體系統的運動學條件，則可對方程(8)進行求解，得到每一塊體的位移和變形狀態。

2.2 爆破荷載的施加及邊界條件

在第1 節中爆破荷載已經被等效為隨時間變化的作用于粉碎區邊緣的壓力，因此在非連續變形分析計算模型中該荷載為一分布力，考慮非連續變形分析塊體單元幾何形狀一般為多邊形，則該分布力體現為非連續變形分析方法中的線荷載形式。因此，爆破荷載的作用能夠非常簡便地在程序中得到體現，在非連續變形分析計算過程中的每一時步開始時，將該時刻對應的粉碎區邊緣壓力施加于邊緣塊體即可。

在模擬爆破荷載在巖體中的動力響應時，必然需要采用有限區域代替無限區域或者半無限區域，這樣建立的模型邊界需要進行一定的處理。經典非連續變形分析中對于邊界通常設置為固定邊界，其不發生任何變形和移動，這樣的邊界對于應力波起著完全反射的作用，因此必須對其做一定的修改。借鑒FLAC3D動力分析中的黏滯邊界來實現應力波的無反射，該邊界可以吸收有限域向外界傳遞的能量，具體操作為在邊界塊體的沿邊界側兩個端點分別施加獨立的法線方向黏滯力f n和切線方向黏滯力fs：

式中：ρ為塊體密度，kg/m3;cp和cs分別為質點的縱波和橫波波速，m/s;vn和vs分別為該端點的法向速度和切向速度，m/s。該黏滯力體現為非連續變形分析方法中的點荷載形式，很容易添加至方程(8)中。

3 不同類別巖體在爆破荷載下的動力響應

在實際工程中，通常根據一定的標準將巖體進行質量分級，GB 50218—1994《工程巖體分級標準》[14]在我國應用廣泛，該標準先考慮巖石的堅硬程度和巖體的完整程度進行基本質量分級：

式中：BQ為巖體質量基本指標;Rc為巖塊飽和單軸抗壓強度，MPa;Kv為巖體完整性系數。然后再根據地下水、巖體軟弱結構面和地應力的情況進行修正，巖體按照計算得到的BQ 值共分為5 級：大于550時為Ⅰ級;450和550 之間為Ⅱ級;350和450 之間為Ⅲ級;250和350 之間為Ⅳ級;小于250 時為Ⅴ級。本文將該問題進行適當簡化，僅考慮由于巖石堅硬程度和完整程度變化而導致的不同級別巖體。由于實際工程中，對于動力響應主要關注點為爆炸荷載在巖體中的傳遞時間和造成的不利振動效應，因此在計算中著重考察距離爆孔一定距離處巖石的起始振動時間及其峰值振速并進行分析。

3.1 堅硬程度差異導致的不同類別巖石的動力響應

采用同樣結構面分布的厚層狀大理巖[15]和厚層狀細粒石英砂巖[16]作為樣本，兩種巖體的相關數據如表1 所示，建立的計算模型如圖2 所示。

表1 堅硬程度不同的巖石樣本質量分級參數

圖2 計算模型

模型長、高均為25m，上下左右四側邊界均為無反射的黏滯邊界。模型中的節理分為天然節理和虛擬節理兩種，設置虛擬節理后巖石在爆破荷載下的斷裂行為能夠被模擬，另外由于非連續變形分析塊體單元均為常應力常應變單元，虛擬節理的添加可以更加精確地分析塊體的應力場和速度場。天然節理不考慮其抗拉強度，內摩擦角取45°，黏聚力取0.05MPa。虛擬節理在計算過程中達到破壞準則后轉變為真實節理，其破壞有拉伸破壞和剪切破壞兩種情況。當虛擬節理面上的拉應力達到節理面抗拉強度值時，拉伸破壞發生;當虛擬節理面上的剪應力達到摩爾-庫倫準則容許的最大剪應力時，剪切破壞發生。虛擬節理的內摩擦角取45°，黏聚力取2.00MPa，抗拉強度取0.5MPa。

計算時采用的兩種巖石塊體材料參數為：大理巖密度2750kg/m3，彈性模量70GPa，泊松比0.10;石英砂巖密度2650 kg/m3，彈性模量55 GPa，泊松比0.20。將第1 節獲得的壓力曲線施加于粉碎區邊緣，分別對距離爆孔中心3m，5m，10 m和20m 處的塊體振動速度進行分析，計算得到的不同堅硬程度的兩種巖體各距離處的塊體振速隨時間變化曲線如圖3 所示，為了便于觀察，僅取了速度曲線的首個半振幅。

圖3 不同堅硬程度巖石振速隨時間變化曲線

觀察圖3可知，隨著距離的增加，塊體的振速峰值變小，起始振動時間變長，但是兩種堅硬程度不同的巖石確存在一定差異。同樣距離的情況下，石英砂巖的起始振動時間晚于大理巖，振速峰值小于大理巖。距離爆孔3 m 處，石英砂巖的起始振動時間比大理巖晚0.3 ms，振速峰值比大理巖小0.2 m/s。隨著距離的增長，這種趨勢愈發明顯，距離爆孔5m處，石英砂巖的起始振動時間比大理巖晚0.5 ms，振速峰值比大理巖小0.4m/s。至距離爆孔較遠處時，由于此時兩種巖體的振速均已較小，所以差別縮小，但起始振動時間的差別仍在繼續增大，距離20 m 時石英砂巖的起始振動時間比大理巖晚2 ms。這說明，在同樣的結構面情況下，巖石的堅硬程度對爆破荷載的傳遞影響主要體現為振速峰值與起始振動時間，同樣距離情況下，巖石愈堅硬，振速峰值愈大，起始振動時間愈短。將計算結果與文獻[16]的結果進行比較，發現堅硬程度對爆破荷載的傳遞影響趨勢基本一致，但是本文計算的起始振動時間與振速峰值均偏低，其原因在于文獻[16]中巖體為均質連續介質，而本文中節理的存在在一定程度上阻礙了爆炸荷載的迅速響應。

3.2 完整程度差異導致的不同類別巖石的動力響應

為分析完整程度差異導致的不同類別巖石動力響應，將前面所做的厚層狀大理巖的結構面進行改變，以獲得一種新的大理巖樣本，并將其與原樣本的計算結果進行比較分析，兩種樣本的相關數據見表2，計算模型的區域尺寸和邊界條件也保持不變，區別在于大理巖2的結構面分布更加稀疏，間距增大。

表2 完整程度不同的巖石樣本質量分級參數

為便于分析比較，仍然采用3.1 節中的大理巖材料參數和節理參數，對距離爆孔中心3m，5m，10m和20 m 處的塊體振速進行計算分析，得到的不同完整程度的兩種巖體各距離處的振速隨時間變化曲線的首個半振幅如圖4 所示。

圖4 不同完整程度巖石振速隨時間變化曲線

根據計算結果，隨著距離的增加，塊體的振速峰值變小和起始振動時間變長這一基本趨勢沒有改變，但是結構面的增加對于爆破荷載的傳遞產生很大的影響。同樣距離的情況下，大理巖2的起始振動時間早于大理巖1，塊體振速峰值大于大理巖1。距離爆孔3 m 處，大理巖2的起始振動時間比大理巖1 早0.2ms，振速峰值比大理巖1 高0.3 m/s。這種差異也有隨距離增長而增大的傾向，距離爆孔5m處，大理巖2的起始振動時間比大理巖1 早0.5 ms，振速峰值比大理巖高0.4 m/s。圖4 顯示，在同樣堅硬的巖石組成情況下，巖石的完整程度也會影響爆破荷載的傳遞，同樣距離時巖石愈完整，振速峰值愈大，起始振動時間愈短。值得關注的是，結構面的增多對于爆破荷載的傳遞時間影響比巖石堅硬程度更加明顯，距離20 m 處，大理巖2的起始振動時間比大理巖1 早3 ms。將計算結果與文獻[6]的結果進行比較，爆破荷載的傳遞受完整程度影響的變化趨勢也基本相似，而且在數值上完整程度高的大理巖2的計算結果更接近于采用動力有限元計算得到的結果，這說明在完整程度較優的情況下，采用均質連續介質近似是可以接受的。

值得注意的是，本文采用的節理仍然是一種非常理想化的力學節理，真實的巖體節理是非常復雜的，其特性及分布狀態同樣也會對巖體中爆破荷載的動力響應造成一定的影響，如朱哲明等[17]針對節理中充填水或者不同材料等情況下應力波傳遞過程進行的分析等。目前非連續變形分析方法仍處于發展之中，對于節理面的不同充填材料和不同充填度等的處理手段仍然不夠成熟，因此對這個問題進行進一步深入的研究是很有價值的。

4 結 論

a.在分析爆破荷載引起的工程巖體中的動力響應問題時非連續變形分析方法能夠較好處理巖體的不連續變形行為，但是在處理巖體中多種復雜的節理特性方面還值得進一步研究，如節理中的充填物質性質和充填程度如何反映等。

b.不同堅硬程度和完整程度的巖體對于爆破荷載的動力響應存在明顯差異。巖石的堅硬程度愈高，完整程度愈好，巖體對爆破荷載的動力響應愈迅速，能量損耗愈小，振速峰值愈大，完整程度對于爆破荷載傳遞時間的影響更為明顯。

[1]高渙渙.爆破荷載對領近洞室的及支護結構的的影響機理分析[D].西安：西安理工大學，2008.

[2]陳占軍，朱傳云，周小恒.爆破荷載作用下巖石邊坡動態響應的FLAC3D模擬研究[J].爆破，2005，22(4)：8-13.

[3]史秀志，周健，崔松，等.露天采礦爆破振動對民房危害預測的DDA 模型及應用[J].中南大學學報：自然科學版，2011，42(2)：441-448.

[4]許紅濤，盧文波，周小恒.爆破振動場動力有限元模擬中爆破荷載的等效施加方法[J].武漢大學學報：工學版，2008，41(1)：67-71.

[5]劉亞群，李海波，李俊如，等.爆破荷載作用下黃麥嶺磷礦巖質邊坡動態響應的UDEC 模擬研究[J].巖石力學與工程學報，2004，23(21)：3659-3663.

[6]張秀麗，焦玉勇，劉泉聲，等.節理對爆炸波影響的數值研究[J].巖土力學，2008，29(3)：712-721.

[7]SHI Gen-hua.Discontinuous deformation analysis：a new numerical model for the statics and dynamics of block systems[D].Berkeley，California：Department of Civil Engineering University of California，1988.

[8]Mac LAUGHLIN M M，DOOLIN D M.Review of validation of the discontinuous deformation analysis(DDA)method[J].International Journal for Numerical and AnalyticalMethodsin Geomechanics，2006，30(4)：271-305.

[9]朱傳云，戴晨，姜清輝.DDA 方法在臺階爆破仿真模擬中的應用[J].巖石力學與工程學報，2002，21(增刊2)：2461-2464.

[10]劉紅巖，楊軍，陳鵬萬.爆破漏斗形成過程的DDA 模擬分析[J].工程爆破，2004，10(2)：17-20.

[11]李寧，SWOBODA G.爆破荷載的數值模擬與應用[J].巖石力學與工程學報，1994，13(4)：357-364.

[12]甯尤軍，楊軍，陳鵬萬.節理巖體爆破的DDA 方法模擬[J].巖土力學，2010，31(7)：2259-2263.

[13]徐全軍，龍源，張慶明，等.微差爆破震動疊加起始位置數值模擬[J].力學與實踐，2000，22(5)：45-48.

[14]GB50218—1994 工程巖體分級標準[S].

[15]田昌貴，陳世華.工程巖體分級標準在地下采礦工程中的應用[J].采礦技術，2005，5(4)：89-93.

[16]蔡斌.國標《工程巖體分級標準》應用中的幾個問題[J].巖土力學，2003，24(增刊1)：74-76.

[17]朱哲明，李元信，周志榮，等.爆炸荷載下缺陷巖體的動力響應[J].巖石力學與工程學報，2011，30(6)：1157-1167.