基于非連續(xù)變形分析方法的爆破荷載下不同類別巖體的動力響應(yīng)

江 巍，程圣國，馮 強(qiáng)

(三峽大學(xué)三峽庫區(qū)地質(zhì)災(zāi)害教育部重點(diǎn)實(shí)驗(yàn)室，湖北 宜昌 443002)

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

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

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

1 爆破荷載的等效

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

經(jīng)驗(yàn)公式法建立的爆破荷載模型將作用在爆孔壁上的爆破荷載壓力簡化為隨時間變化的函數(shù)，基本形式為

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

考慮爆生氣體狀態(tài)的爆破荷載模型將爆孔內(nèi)的爆生氣體壓力P 與爆腔體積V 對應(yīng)起來：

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

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

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

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

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

圖1 粉碎區(qū)邊緣壓力曲線

該方法的基本思想為：當(dāng)藥包爆炸后，以爆孔為中心周圍的巖體按照由近至遠(yuǎn)形成粉碎區(qū)、破碎區(qū)和振動區(qū)，將爆破荷載直接施加在粉碎區(qū)的邊緣，則可以避開粉碎區(qū)的巖石力學(xué)特性定義難題?；诖朔椒ㄐ烊姷萚13]采用非線性有限元程序?qū)Ρ普駝有?yīng)進(jìn)行了分析，取得了一定的成果，但在計算中為使用有限元技術(shù)將破碎區(qū)和振動區(qū)的巖體均近似為均勻連續(xù)介質(zhì)。本文以非連續(xù)變形分析方法為手段，視破碎區(qū)和振動區(qū)的巖體由線彈性體的巖塊和節(jié)理組成，破碎區(qū)巖體中除巖塊發(fā)生彈性變形外更主要的是節(jié)理新生、擴(kuò)張和延伸等導(dǎo)致巖體發(fā)生破碎行為，振動區(qū)的巖體中節(jié)理受力過小變形細(xì)微，兩側(cè)巖塊之間力的傳遞良好，其力學(xué)行為表現(xiàn)接近連續(xù)介質(zhì)。

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

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

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

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

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

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

施加邊界條件和塊體系統(tǒng)的運(yùn)動學(xué)條件，則可對方程(8)進(jìn)行求解，得到每一塊體的位移和變形狀態(tài)。

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

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

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

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

3 不同類別巖體在爆破荷載下的動力響應(yīng)

在實(shí)際工程中，通常根據(jù)一定的標(biāo)準(zhǔn)將巖體進(jìn)行質(zhì)量分級，GB 50218—1994《工程巖體分級標(biāo)準(zhǔn)》[14]在我國應(yīng)用廣泛，該標(biāo)準(zhǔn)先考慮巖石的堅硬程度和巖體的完整程度進(jìn)行基本質(zhì)量分級：

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

3.1 堅硬程度差異導(dǎo)致的不同類別巖石的動力響應(yīng)

采用同樣結(jié)構(gòu)面分布的厚層狀大理巖[15]和厚層狀細(xì)粒石英砂巖[16]作為樣本，兩種巖體的相關(guān)數(shù)據(jù)如表1 所示，建立的計算模型如圖2 所示。

表1 堅硬程度不同的巖石樣本質(zhì)量分級參數(shù)

圖2 計算模型

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

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

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

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

3.2 完整程度差異導(dǎo)致的不同類別巖石的動力響應(yīng)

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

表2 完整程度不同的巖石樣本質(zhì)量分級參數(shù)

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

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

根據(jù)計算結(jié)果，隨著距離的增加，塊體的振速峰值變小和起始振動時間變長這一基本趨勢沒有改變，但是結(jié)構(gòu)面的增加對于爆破荷載的傳遞產(chǎn)生很大的影響。同樣距離的情況下，大理巖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 顯示，在同樣堅硬的巖石組成情況下，巖石的完整程度也會影響爆破荷載的傳遞，同樣距離時巖石愈完整，振速峰值愈大，起始振動時間愈短。值得關(guān)注的是，結(jié)構(gòu)面的增多對于爆破荷載的傳遞時間影響比巖石堅硬程度更加明顯，距離20 m 處，大理巖2的起始振動時間比大理巖1 早3 ms。將計算結(jié)果與文獻(xiàn)[6]的結(jié)果進(jìn)行比較，爆破荷載的傳遞受完整程度影響的變化趨勢也基本相似，而且在數(shù)值上完整程度高的大理巖2的計算結(jié)果更接近于采用動力有限元計算得到的結(jié)果，這說明在完整程度較優(yōu)的情況下，采用均質(zhì)連續(xù)介質(zhì)近似是可以接受的。

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

4 結(jié) 論

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

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

[1]高渙渙.爆破荷載對領(lǐng)近洞室的及支護(hù)結(jié)構(gòu)的的影響機(jī)理分析[D].西安：西安理工大學(xué)，2008.

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

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

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

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

[6]張秀麗，焦玉勇，劉泉聲，等.節(jié)理對爆炸波影響的數(shù)值研究[J].巖土力學(xué)，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 方法在臺階爆破仿真模擬中的應(yīng)用[J].巖石力學(xué)與工程學(xué)報，2002，21(增刊2)：2461-2464.

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

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

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

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

[14]GB50218—1994 工程巖體分級標(biāo)準(zhǔn)[S].

[15]田昌貴，陳世華.工程巖體分級標(biāo)準(zhǔn)在地下采礦工程中的應(yīng)用[J].采礦技術(shù)，2005，5(4)：89-93.

[16]蔡斌.國標(biāo)《工程巖體分級標(biāo)準(zhǔn)》應(yīng)用中的幾個問題[J].巖土力學(xué)，2003，24(增刊1)：74-76.

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