馬泗洲，劉科偉，楊家彩，李旭東，郭騰飛

（中南大學(xué)資源與安全工程學(xué)院，湖南 長沙 410083）

隨著地表淺層礦產(chǎn)資源的逐漸枯竭，通過深部開采獲取資源將逐漸常態(tài)化[1-2]。與傳統(tǒng)淺部開采相比，深部開采需要更加先進的技術(shù)和設(shè)備，從而保證資源的有效利用與合理開發(fā)。其中，鉆爆法仍是深部硬巖開挖的主要手段，廣泛應(yīng)用于巷道開挖、井筒掘進和采場爆破等工程領(lǐng)域。淺部開采所受地應(yīng)力較小，地壓對巖體爆破碎裂的影響可忽略不計，而深部開采所受地應(yīng)力較大，千米深井實測地應(yīng)力結(jié)果可達30～50 MPa[3]，深部巖體在高地應(yīng)力下爆破時會表現(xiàn)出與淺部不同的力學(xué)行為。

關(guān)于地應(yīng)力下巖體爆破，國內(nèi)外學(xué)者已開展諸多試驗、理論及數(shù)值計算方面的研究。其中，最早的研究可追溯至1966 年由Nicholls 等[4]開展的地應(yīng)力下預(yù)裂爆破試驗，首次提出爆破徑向裂紋更傾向于沿最大主應(yīng)力方向擴展?；谠撛囼?，Kutter 等[5]通過不同材料小型爆破試驗發(fā)現(xiàn)單向壓力對爆破裂紋的萌生和擴展具有導(dǎo)向作用，從而印證了Nicholls 等的結(jié)論。國內(nèi)方面，肖正學(xué)等[6]基于前人的研究，結(jié)合室內(nèi)與現(xiàn)場爆破試驗數(shù)據(jù)，發(fā)現(xiàn)了類似的爆破裂紋擴展行為規(guī)律，并指出初始應(yīng)力場會改變爆轟波的傳播規(guī)律。此后，劉殿書等[7]基于動光彈實驗記錄的動態(tài)等差條紋圖，分析了不同初始應(yīng)力下爆破應(yīng)力波的變化過程，驗證了肖正學(xué)等[6]的觀點，并將該結(jié)論進一步具體化。在爆破試驗測量儀器方面，楊仁樹等[8-9]將動態(tài)焦散線方法與高速攝影技術(shù)相結(jié)合，建立了一套適用于爆炸、侵徹等超動態(tài)問題研究的數(shù)字激光動態(tài)焦散線試驗系統(tǒng)，通過該系統(tǒng)解決了巖體爆破領(lǐng)域難以實時觀測爆破過程的一些問題。岳中文等[10]通過該系統(tǒng)分析了單向初始壓力作用下切縫藥包爆破裂紋的擴展行為，指出切縫藥包可以有效控制爆炸能量分布，使能量沿切縫方向集中釋放，進而實現(xiàn)控制爆破中材料的定向斷裂。隨后，Yang 等[11]進一步開展了切縫藥包控制爆破相關(guān)研究，基于彈性與斷裂力學(xué)理論，分析了初始應(yīng)力對切縫藥包爆破裂紋擴展行為的影響，討論了不同切縫角度下的爆破損傷特征，揭示了深部巖體切縫藥包爆破的破裂機理。

此外，數(shù)值模擬以其簡便/普適性，被廣泛應(yīng)用于巖體爆破過程研究。常用的數(shù)值方法主要包括有限元、離散元兩類，其中使用最為廣泛的軟件當(dāng)屬LS-DYNA。Ma 等[12]基于該軟件中的Johnson-Holmquist-2(JH-2) 模型模擬了單孔爆破裂紋擴展過程，討論了應(yīng)力加載速率、自由面距離、節(jié)理位置和初始壓力等因素對爆破裂紋擴展行為的影響，并用于控制爆破工程指導(dǎo)。Xie 等[13-14]基于該軟件中的拉壓損傷Riedel-Hiermaier-Thoma (RHT) 模型模擬了地應(yīng)力下巖體掏槽爆破裂紋擴展過程，發(fā)現(xiàn)了不同壓力條件下掏槽爆破效果不佳的原因，即靜水壓力抑制裂紋擴展及非靜水壓力影響損傷分布，進而提出改變掏槽孔距離及布置方位的優(yōu)化設(shè)計方法。Li 等[15-16]通過數(shù)值模擬與圖像處理相結(jié)合的方法研究了深部巖體爆破的塊度特征，分析了地應(yīng)力對巖體爆破的塊度尺寸及形狀影響，研究結(jié)果表明，隨著地應(yīng)力的增加，巖體爆破塊度會更大更圓，且塊度粒徑分布的范圍會更廣。

上述研究多集中在初始應(yīng)力作用下巖體爆破動態(tài)力學(xué)響應(yīng)現(xiàn)象分析，有必要深入了解初始應(yīng)力對巖體爆破破裂機制的影響。此外，如何定量分析巖體爆破裂紋擴展行為和損傷特征與初始應(yīng)力的關(guān)系也有待深入挖掘。本文基于彈性力學(xué)建立初始應(yīng)力作用下單孔爆破動靜組合理論模型，分析靜載模型應(yīng)力分布及動載模型應(yīng)力演化過程，討論初始應(yīng)力對巖體爆破損傷的影響，探索破裂機理；采用顯式動力學(xué)分析軟件LS-DYNA，介紹RHT 模型中巖石材料參數(shù)的標定及修正方法，并通過爆破實驗進行模型驗證；使用驗證的模型分析不同初始應(yīng)力作用下巖體爆破應(yīng)力波演化過程及其對裂紋擴展行為的影響，并通過分形維數(shù)對爆破損傷量化分析，為深部巖體爆破工程提供理論指導(dǎo)。

1 動靜組合理論模型

初始應(yīng)力作用下巖體爆破問題可簡化為平面應(yīng)變問題，如圖1 所示。假定含圓孔無限大平板為均勻連續(xù)且各向同性的彈性介質(zhì)，外側(cè)分別受水平初始應(yīng)力σx和豎直初始應(yīng)力σy作用，圓孔內(nèi)部受動態(tài)載荷p(t) 作用，初始應(yīng)力下巖體單孔爆破模型可看作靜態(tài)荷載和動態(tài)荷載共同作用下的組合模型。

圖1 理論分析幾何模型Fig. 1 Geometrical model for theoretical analysis

1.1 靜態(tài)應(yīng)力分布特征

根據(jù)彈性力學(xué)理論，可分別計算不同初始應(yīng)力條件下炮孔近場(r=a)與遠場(r=10a)的巖體徑向應(yīng)力、環(huán)向應(yīng)力和剪切應(yīng)力[17]

式中：a和r分別為炮孔半徑和距炮孔中心的距離；定義豎直方向應(yīng)力σy為壓力p，取值20 MPa；K定義為側(cè)壓系數(shù)，K=σx/σy。

炮孔近場和遠場的應(yīng)力分布特征如圖2 所示，近場為孔壁自由面，徑向應(yīng)力與剪切應(yīng)力均為0。由圖2(a)和2(b)可知：水平方向，炮孔近場環(huán)向應(yīng)力隨側(cè)壓系數(shù)的增大呈線性減小的趨勢，而炮孔遠場幾乎不受側(cè)壓系數(shù)影響；豎直方向，炮孔近場與遠場環(huán)向應(yīng)力均隨側(cè)壓系數(shù)的增大呈線性增長趨勢，且近場的增長速率更加顯著。此外，側(cè)壓系數(shù)K=0 時，孔壁在水平與豎直方向分別出現(xiàn)了壓應(yīng)力和拉應(yīng)力集中，考慮到巖體抗拉/壓強度特性，初始拉應(yīng)力存在處巖體破壞將更為嚴重，根據(jù)環(huán)向應(yīng)力變化規(guī)律可知，拉應(yīng)力轉(zhuǎn)化為壓應(yīng)力的臨界側(cè)壓系數(shù)K=1/3。炮孔遠場徑向應(yīng)力與環(huán)向應(yīng)力分布規(guī)律相似，僅在方向上發(fā)生90°偏轉(zhuǎn)，如圖2(c)所示，即豎直方向上徑向應(yīng)力的大小與側(cè)壓系數(shù)無關(guān)，而水平方向上隨側(cè)壓系數(shù)的增大徑向應(yīng)力也隨之增加。剪切應(yīng)力在水平與豎直方向上均為0，其他方向隨側(cè)壓系數(shù)增大呈先減后增的規(guī)律，其方向會在K=1.0 前后發(fā)生變化，如圖2(d)所示。

圖2 炮孔周邊靜態(tài)應(yīng)力分布Fig. 2 Stress distribution around the borehole under static load

距炮孔中心不同位置各應(yīng)力變化過程如圖3 所示，水平與豎直方向上各應(yīng)力隨距離的增加均趨于平緩變化，不同的是，環(huán)向與徑向應(yīng)力主導(dǎo)方向會在側(cè)壓系數(shù)K=1.0 前后發(fā)生變化。K＞1.0 時，水平方向以環(huán)向應(yīng)力為主導(dǎo)，如圖3(a)所示，K＜1.0 時，豎直方向以環(huán)向應(yīng)力為主導(dǎo)，如圖3(b)所示。通常認為環(huán)向應(yīng)力是產(chǎn)生徑向裂紋的主要原因，因此，各向異性壓力條件下徑向裂紋在不同方向上的分布有所差異，且隨著水平、豎直方向初始應(yīng)力差的增加，徑向裂紋擴展各向異性也會更加明顯。

圖3 不同位置處巖體靜態(tài)應(yīng)力的變化Fig. 3 Stress variation of rock mass at different distances under static load

1.2 動態(tài)應(yīng)力演化規(guī)律

將炸藥產(chǎn)生的爆炸荷載等效為動態(tài)荷載p(t)作用于炮孔壁處，p(t)=pVN(eγ/n)ntne-γt，其中，pVN為炸藥起爆后炮孔壁處峰值壓力，n為控制應(yīng)力波演化的模型參數(shù)，通常取整數(shù)， γ 為壓力衰減指數(shù)。根據(jù)上述基本假設(shè)，極坐標系下巖體中爆炸應(yīng)力波的控制方程可表示為[18]：

式中： φ(r,t) 為位移勢函數(shù)，t為應(yīng)力波傳播時間，vp為巖體的縱波波速，(r,t) 為爆炸荷載下巖體的徑向應(yīng)力，其動態(tài)徑向應(yīng)力和環(huán)向應(yīng)力隨時間變化的關(guān)系可由下式求得：

式(4) 經(jīng)過拉普拉斯逆變換后并不能求得其具體的解析表達式，因此需要結(jié)合數(shù)值反演的方法進行求解[19]。目前關(guān)于數(shù)值反演的算法有很多，如Durbin 算法、Crump 算法和Stehfest 算法等，其中Durbin 算法穩(wěn)定可靠性較好[20-21]，通過該算法可求得其爆炸荷載下巖體中應(yīng)力波隨時間的變化規(guī)律：

式中：α 可取任意實數(shù)，0≤α≤Re(m)；T為求解時間間隔，0≤t≤T/2。

距炮孔中心不同位置處巖體動態(tài)應(yīng)力波演化規(guī)律如圖4 所示。爆炸荷載下動態(tài)徑向應(yīng)力先是急劇增加，增長至峰值應(yīng)力后迅速衰減，且孔壁處應(yīng)力最大，如圖4(a)所示。較大的徑向壓應(yīng)力會使得炮孔近區(qū)巖體破壞程度加劇，最終形成壓縮粉碎區(qū)。相較于孔壁處而言，其他位置的徑向應(yīng)力會出現(xiàn)由壓縮應(yīng)力轉(zhuǎn)變?yōu)槔鞈?yīng)力的現(xiàn)象。需要說明的是，徑向拉應(yīng)力的存在會進一步誘導(dǎo)環(huán)向裂紋形成，隨著距離的增加，峰值徑向壓應(yīng)力逐漸減小，且衰減速度相對較快，而峰值徑向拉應(yīng)力幾乎沒有發(fā)生顯著變化，整體呈先增后減，最后趨于平緩的變化趨勢，如圖4(b)所示。動態(tài)環(huán)向應(yīng)力與徑向應(yīng)力的變化規(guī)律有所區(qū)別，如圖4(c)所示，其應(yīng)力時程曲線出現(xiàn)三處應(yīng)力峰值：爆炸壓應(yīng)力達到峰值后迅速衰減轉(zhuǎn)變?yōu)槔瓚?yīng)力，達到拉應(yīng)力峰值后又繼續(xù)衰減轉(zhuǎn)變?yōu)閴簯?yīng)力，相較于壓應(yīng)力峰值而言，拉應(yīng)力峰值相對較大，這也是徑向裂紋較為發(fā)育的原因之一。峰值環(huán)向壓應(yīng)力與拉應(yīng)力隨距離的變化趨勢則剛好與峰值徑向應(yīng)力相反，其拉應(yīng)力峰值相較于壓應(yīng)力較高，壓應(yīng)力隨距離的增加無明顯的變化，且距炮孔中心5 倍的炮孔半徑之外，拉伸與壓縮應(yīng)力變化趨勢幾乎相同，如圖4(d)所示。

圖4 不同距離處巖體的動態(tài)應(yīng)力波演化Fig. 4 Dynamic stress waves evolution in rock mass at various distances

2 數(shù)值計算模型及驗證

理論模型將巖石視為彈性介質(zhì)是為了便于獲取數(shù)學(xué)解析表達式，而在實際工程中，巖體爆破時的力學(xué)行為是復(fù)雜且非線性的。因此，在彈性理論模型基礎(chǔ)上，通過非線性數(shù)值模型可彌補理論模型的不足并論證理論模型的合理性。數(shù)值模擬計算結(jié)果的可靠性很大程度取決于材料本構(gòu)方程，LS-DYNA 軟件中常用于模擬巖石材料的有RHT 模型、JH-2 模型、CSCM 模型等。通過數(shù)值計算結(jié)果發(fā)現(xiàn)，RHT 材料模型可以更準確地描述深部巖體受到爆炸荷載所表現(xiàn)的力學(xué)響應(yīng)和破壞行為[22]。

2.1 RHT 模型參數(shù)標定

RHT 模型參數(shù)標定需要巖石的基本物理力學(xué)參數(shù)，通過試驗測得花崗巖的力學(xué)參數(shù)包括：密度（2 620 kg·m-3）、單軸抗壓強度（162 MPa）、泊松比（0.18）、剪切模量（21.9 GPa）、縱波波速（4 600 m/s）、橫波波速（2 820 m/s）。巖石的抗壓/拉強度具有明顯的應(yīng)變率效應(yīng)，在RHT 巖石材料模型中，材料的應(yīng)變強化效應(yīng)由應(yīng)變率增強因子Fr為:

A和N為材料的破壞面參數(shù)，用于描述RHT 模型中材料的破壞面強度。通過歸一化壓力滿足3＞Fr，得到：

RHT 材料模型相較于其他模型的優(yōu)點之一是引入了偏應(yīng)力張量第三不變量及羅德角，其中羅德角因子可以描述失效面子午線失效壓縮強度的折減，折減值Q與相對壓力的函數(shù)關(guān)系可表示為：

式中：Q0為拉伸與壓縮荷載下材料子午線半徑之比，B為巖石材料羅德角相關(guān)系數(shù)，參考文獻[23]中的試驗數(shù)據(jù)，進行線性回歸求得：Q0=0.68，B=0.05，如圖5(b)所示。

表1 不同圍壓下花崗巖力學(xué)參數(shù)Table 1 Mechanical parameters of granite sample under various confining pressure

圖5 RHT 模型破壞面的相關(guān)參數(shù)擬合Fig. 5 Fitting curve of failure surface parameters for RHT model

RHT 模型中材料的破壞積累量由損傷參數(shù)Dr確定，其損傷變量和累積塑性應(yīng)變表達式如下：

RHT 模型中剪切應(yīng)力分量與壓應(yīng)力分量相耦合，壓力與應(yīng)力的關(guān)系由Mie-Grüneisen 形式的多項式曲線連同p-α 孔隙壓實曲線關(guān)系描述，其狀態(tài)方程為：

根據(jù)上述經(jīng)驗公式可確定RHT 模型大部分參數(shù)，其余難以確定且對數(shù)值計算結(jié)果比較敏感參數(shù)，如屈服面參數(shù)和，殘余面參數(shù)Af和Nf等，可在RHT 模型原始文獻[24]推薦取值的基礎(chǔ)上，結(jié)合分離式霍普金森壓桿(SHPB)動態(tài)力學(xué)試驗進行修正和優(yōu)化[25]。SHPB 系統(tǒng)中入射桿和透射桿的長度分別為2 000 和1 500 mm，巴西圓盤巖樣的直徑和厚度分別為50 和25 mm，通過試錯法不斷調(diào)整RHT 模型中的參數(shù)，直至數(shù)值模擬結(jié)果與試驗結(jié)果匹配[26]，其比較結(jié)果如圖6 和圖7 所示。

圖6 數(shù)值模擬與試驗中巖體沖擊破壞過程比較Fig. 6 Comparison of failure process between numerical simulation and test

圖7 數(shù)值模擬與試驗中巖體應(yīng)力時程曲線的對比Fig. 7 Comparison of stress with time between numerical simulation and test

圓盤試樣在沖擊荷載下的破壞過程及特征如圖6 所示，試樣與桿件接觸端裂紋萌生起裂，此后裂紋沿弱面不斷擴展延伸，進而兩端的裂紋相遇聚結(jié)，最終巖樣被裂紋貫穿產(chǎn)生宏觀大面積裂隙。此外，在接觸面會產(chǎn)生較大的破碎區(qū)，試樣整體為拉伸劈裂破壞。在對比巖樣動態(tài)破壞過程及特征的基礎(chǔ)上，進一步比較了其應(yīng)力時程曲線，如圖7 所示。首先進行動態(tài)應(yīng)力平衡驗證，如圖7(a)所示，試樣兩端處于較好的受力平衡狀態(tài)。然后基于三波法求得應(yīng)力時程曲線，如圖7(b)所示，由圖可知，巖樣的應(yīng)力時程曲線變化趨勢基本一致，且峰值強度相同。值得注意的是，試驗中的應(yīng)力時程曲線相對粗糙，且峰后表現(xiàn)的延性更為顯著，這可能與試樣的均質(zhì)性有關(guān)。為說明破壞過程與應(yīng)力時程曲線的對應(yīng)關(guān)系，將破壞過程中的特征點對應(yīng)地標注在應(yīng)力時程曲線上（其中ti、tg、tc和tp分別代表裂紋的起裂、增長、凝聚和貫穿的時間）。總體而言，數(shù)值模擬與實驗結(jié)果整體吻合較好，該模型能夠較準確地描述巖樣損傷區(qū)域范圍及動態(tài)破壞過程，具體的模型參數(shù)如表2 所示。

表2 巖石RHT 模型材料參數(shù)Table 2 RHT model parameters for rock mass

2.2 炸藥及空氣材料

LS-DYNA 軟件采用高能炸藥模型*MAT_HIGH_EXPLOSIVE_BURN 模擬試驗中的炸藥材料，并結(jié)合Jones-Wilkens-Lee (JWL)狀態(tài)方程來描述其爆炸壓力、體積與能量之間的關(guān)系：

式中：pe為爆轟壓力，V為相對體積，為炸藥單位體積內(nèi)能，Ae、Be、R1、R2、ω 為材料常數(shù)。具體炸藥模型的材料參數(shù)如表3 所示[27]，其中vd和pCJ分別表示爆轟速度和初始壓力。

表3 炸藥模型材料參數(shù)Table 3 Parameters for the explosive material

爆破時常采用不耦合裝藥，LS-DYNA 中空氣介質(zhì)使用典型材料模型*MAT_NULL，并結(jié)合線性多項式狀態(tài)方程*EOS_LINEAR_POLYNOMIAL 定義氣體壓力、密度和內(nèi)部能量之間的關(guān)系：

式中：pa為氣體壓力，為單位氣體體積的內(nèi)部能量，u為動態(tài)黏度系數(shù)，u=(ρ/ρ0)-1 ，ρ 和ρ0分別是材料的密度和初始密度，C4=C5=γ-1 ，γ 為比熱比，通常C0=C1=C2=C3=C6=0。具體空氣模型的材料參數(shù)如表4 所示[27]。

表4 空氣模型材料參數(shù)Table 4 Mateiral parameters for the air

2.3 數(shù)值模型驗證

數(shù)值模型參數(shù)確定后，可通過相關(guān)試驗[27]驗證選取的巖石、炸藥及空氣材料的合理性。試驗中圓柱巖樣的高度和直徑分別為150 和144 mm，炮孔和藥卷直徑分別為6.45 和1.70 mm，藥卷周圍包裹聚乙烯，并在其外側(cè)嵌套銅管，確保試樣爆裂后碎塊不會飛濺，便于后期掃描完整斷面的裂紋形態(tài)。在ANSYS 軟件中建立相同尺寸的幾何模型并劃分網(wǎng)格，模型共包含約332 萬個單元和338 萬個節(jié)點，其詳細結(jié)構(gòu)尺寸及橫截面局部網(wǎng)格如圖8 所示。

圖8 數(shù)值模型及局部網(wǎng)格Fig. 8 Configuration of numerical model and local mesh

炸藥起爆后產(chǎn)生大量沖擊波作用于孔壁，其應(yīng)力峰值遠高于巖體的動態(tài)抗壓強度，在沖擊波作用下，炮孔周圍形成壓剪粉碎區(qū)（Zone Ⅰ）。沖擊波穿過耦合介質(zhì)消耗許多能量進而衰減轉(zhuǎn)變成應(yīng)力波，其應(yīng)力峰值低于巖體抗壓強度，因而壓剪粉碎區(qū)的范圍不再繼續(xù)擴展。然而巖石動態(tài)抗拉強度相比抗壓強度較低，因此在應(yīng)力波作用下形成了拉剪破碎區(qū)（Zone Ⅱ）。應(yīng)力波在向自由面?zhèn)鞑サ耐瑫r，徑向裂紋沿裂隙尖端繼續(xù)擴展，最終形成邊界處的裂隙發(fā)育區(qū)（Zone Ⅲ）。除徑向裂紋外，在拉剪破碎區(qū)還有明顯的環(huán)向裂紋產(chǎn)生，這與實驗和理論也是相符的，其比較結(jié)果如圖9 所示。從裂紋區(qū)域分布及擴展特征分析，數(shù)值模擬與試驗和理論結(jié)果具有高度一致性。

圖9 巖體爆破裂紋形態(tài)理論、試驗與模擬結(jié)果比較Fig. 9 Comparison of rock blasting crack patterns among theoretical, experimental and simulated results

為進一步驗證數(shù)值模型的可靠性及合理性，在比較爆破裂紋擴展過程及分布特征的基礎(chǔ)上，進一步對比分析了爆炸壓力分布。其結(jié)果數(shù)據(jù)主要來源于理論計算、試驗測量及數(shù)值模擬，其中，理論計算結(jié)果可由經(jīng)驗公式[22]求得：

式中：pmax為炸藥的起爆峰值壓力，ρs和ρe分別為巖體及炸藥的密度，vd為炸藥的爆速，vp為巖體的波速，γe為爆轟產(chǎn)物的絕熱膨脹系數(shù)。壓力衰減規(guī)律可由下式求得：

式中：p(r)為爆炸壓力，ζ 為爆炸壓力衰減指數(shù)。

巖體不同位置處的理論峰值壓力可由式(13)與式(14)聯(lián)合計算求得，而數(shù)值模擬結(jié)果則是通過沿徑向布置若干單元測點輸出的峰值壓力，試驗結(jié)果為Banadaki 室內(nèi)試驗測量的壓力，如圖10 所示。由圖可知，理論計算與數(shù)值中峰值壓力的衰減趨勢基本一致，且試驗中測得的壓力值也基本分布在其壓力衰減曲線上。通過裂紋分布特征及壓力衰減規(guī)律來看，數(shù)值模擬與試驗和理論結(jié)果較為吻合，故該模型適用于本文的研究工作。

圖10 壓力衰減曲線比較Fig. 10 Comparison of attenuation curve of peak pressure

3 單孔爆破數(shù)值模擬

在數(shù)值模型驗證的基礎(chǔ)上，進行單孔爆破數(shù)值模擬計算，分析不同初始壓力下巖體靜態(tài)應(yīng)力初始化分布及壓力演化過程對裂紋擴展行為及損傷特征的影響。模擬主要分為2 組進行：各向同性初始壓力與各向異性初始壓力條件下的單孔巖體爆破，不同工況下的應(yīng)力加載條件如表5 所示。

表5 初始應(yīng)力加載條件Table 5 Initial stress conditions in numerical simulation

為驗證理論模型的合理性，數(shù)值模型采用平面單層網(wǎng)格，約束垂直平面位移確保計算在平面應(yīng)變條件下進行。模型幾何尺寸的邊長為4 m，其中炮孔半徑為21 mm，不耦合系數(shù)為1.3。模型中單元數(shù)量共計約108 萬個，單元的平均尺寸為4 mm。為保證計算效率，且要解決爆炸大變形可能引起的網(wǎng)格畸變問題，數(shù)值模型中巖體單元采用拉格朗日算法，炸藥與空氣單元則采用任意拉格朗日歐拉算法，通過流固耦合關(guān)鍵字與多物質(zhì)組關(guān)鍵字實現(xiàn)爆炸荷載的有效傳遞。在模型四周施加圍壓來模擬初始應(yīng)力，并設(shè)置無反射邊界條件以模擬無限區(qū)域巖體，在距炮孔中心不同位置處布置若干測點，記錄壓力變化過程，數(shù)值模型及壓力測點布置如圖11 所示。

3.1 靜態(tài)應(yīng)力初始化

初始應(yīng)力下巖體爆破數(shù)值模擬分兩步計算：第一步在邊界施加初始應(yīng)力，待應(yīng)力初始化穩(wěn)定后再進行爆破計算。應(yīng)力初始化方法有很多，如動力松弛法、顯-隱式轉(zhuǎn)化法和dynain 文件法等，綜合考慮模擬效果和計算效率，采用dynain 文件法更穩(wěn)定，該方法可用顯/隱式求解功能實現(xiàn)應(yīng)力初始化[28]。使用關(guān)鍵字*INTERFACE_SPRINGBACK_LSDYNA 將應(yīng)力初始化結(jié)果輸出，然后通過關(guān)鍵字*INCLUDE 將第一步的計算結(jié)果文件包含進來，再進行第二步計算，具體的計算流程如圖12 所示。

圖12 初始應(yīng)力下巖體爆破計算流程Fig. 12 Flow chart for the rock blast under initial stress

為說明施加于巖體上的初始應(yīng)力處于穩(wěn)定狀態(tài)，以工況E-3 為例進行分析。選取炮孔周邊的四個測點，輸出其壓力時程曲線，如圖13 所示。初始應(yīng)力時程曲線主要包括兩個階段——應(yīng)力增長階段與穩(wěn)定階段，加載初期應(yīng)力以一定速度穩(wěn)定增長，增長至預(yù)設(shè)的壓力目標值時便趨于穩(wěn)定狀態(tài)。此外，不同測點的應(yīng)力時程曲線基本重合，因此施加于巖體上的初始應(yīng)力是穩(wěn)定可靠的。

圖13 應(yīng)力初始化驗證Fig. 13 Verification of stress initialization

在后處理軟件LS-PREPOST 中建立局部柱坐標系，繪制巖體徑向、環(huán)向及剪切應(yīng)力的分布云圖，如圖14～圖18 所示。云圖可直觀地顯示炮孔周圍巖體的應(yīng)力分布，如側(cè)壓系數(shù)K=0 時，豎直方向上徑向與環(huán)向應(yīng)力出現(xiàn)局部拉應(yīng)力集中，剪切應(yīng)力在一、三象限為壓應(yīng)力，二、四象限為拉應(yīng)力，如圖14 所示；側(cè)壓系數(shù)K=1.0 時，徑向與環(huán)向應(yīng)力在炮孔壁處成同心圓狀分布，且不存在剪切應(yīng)力，如圖16 所示。上述的模擬結(jié)果與1.1 節(jié)中理論模型分析的結(jié)果高度一致，再次驗證了數(shù)值模型的合理性及應(yīng)力初始化方法的可靠性。

圖14 應(yīng)力初始化模擬結(jié)果 (K=0)Fig. 14 Numerical results of the stress initialization (K=0)

圖15 應(yīng)力初始化模擬結(jié)果 (K=0.5)Fig. 15 Numerical results of the stress initialization (K=0.5)

圖16 應(yīng)力初始化模擬結(jié)果 (K=1.0)Fig. 16 Numerical results of the stress initialization (K=1.0)

圖17 應(yīng)力初始化模擬結(jié)果 (K=1.5)Fig. 17 Numerical results of the stress initialization (K=1.5)

圖18 應(yīng)力初始化模擬結(jié)果 (K=2.0)Fig. 18 Numerical results of the stress initialization (K=2.0)

3.2 壓力演化過程

由理論分析可知，環(huán)向應(yīng)力對徑向主裂紋的擴展行為具有決定性作用，為進一步探究裂紋的擴展機制，本節(jié)主要選取E-1（無初始壓力）、E-4（各向同性壓力）及A-4（各向異性壓力）三種工況，分析環(huán)向應(yīng)力的演化過程，如圖19～圖21 所示。工況E-1 及E-4 下環(huán)向應(yīng)力演化規(guī)律相似，如圖19 與圖20 所示。爆炸初期，炮孔周邊產(chǎn)生較強的爆轟波，環(huán)向應(yīng)力呈圓形輻射狀傳播，接著沖擊波能量不斷衰減，環(huán)向應(yīng)力逐漸降低并趨于穩(wěn)定狀態(tài)。此外，相較于無初始壓力而言，各向同性壓力作用下環(huán)向拉應(yīng)力波的擴展范圍及擴展時間受到明顯的抑制作用：無初始壓力時，其擴展時間停留在0.5 ms 前后，40 MPa 初始壓力時，其擴展時間僅停留在0.3 ms 前后，且其擴展范圍也大幅度縮小。相較于各向同性壓力條件，各向異性壓力下環(huán)向應(yīng)力演化過程有顯著差別，尤其在0.2 ms 后，環(huán)向拉應(yīng)力波在水平與豎直方向分布出現(xiàn)差異，差異性隨時間增加愈加明顯，其傳播方向逐漸集中于最大主應(yīng)力方向，如圖21 所示。

圖20 各向同性壓力下環(huán)向應(yīng)力的演化過程 (E-4)Fig. 20 Evolution of hoop stress under equibiaxial pressure (E-4)

圖21 各向異性壓力下環(huán)向應(yīng)力的演化過程 (A-4)Fig. 21 Evolution of hoop stress under anisotropic pressure (A-4)

根據(jù)監(jiān)測點記錄的數(shù)據(jù)，繪制相應(yīng)的環(huán)向應(yīng)力演化曲線，如圖22 所示。各向同性壓力下水平與豎直方向監(jiān)測點是中心對稱的，因此其環(huán)向應(yīng)力時程曲線基本重合，如圖22(a)所示。而各向異性壓力下水平與豎直方向上的環(huán)向應(yīng)力時程曲線有所區(qū)別，主要表現(xiàn)在拉壓應(yīng)力峰值的差異性，最大主應(yīng)力方向上的壓應(yīng)力削弱，而拉應(yīng)力增強，如圖22(b)所示。這可以很好地解釋裂紋傾向于最大主應(yīng)力方向擴展的原因。值得注意的是，不同初始壓力下，其環(huán)向應(yīng)力的變化規(guī)律具有相似性，即爆炸初期因其強烈的沖擊作用產(chǎn)生較高的壓應(yīng)力，巖體發(fā)生壓縮破壞，隨后迅速衰減轉(zhuǎn)變?yōu)槔瓚?yīng)力，巖體發(fā)生拉伸破壞，最后趨于穩(wěn)定狀態(tài)，其演化規(guī)律與1.2 節(jié)中理論模型分析結(jié)果一致。

圖22 不同壓力條件下環(huán)向應(yīng)力時程曲線Fig. 22 Time histories of hoop stress under various pressure conditions

不同距離處峰值環(huán)向應(yīng)力的變化特征如圖23 所示。各向同性壓力下，峰值環(huán)向壓應(yīng)力的變化趨勢一致，即隨著距離的增加而降低，降低速率呈先增后減的規(guī)律。同樣地，峰值環(huán)向拉應(yīng)力隨距離的增加也顯著降低，且隨著初始壓力的增加，其降低速率變得更加明顯，如圖23(a)所示。各向異性壓力下，峰值環(huán)向拉應(yīng)力在水平與豎直方向有所差異，最大主應(yīng)力方向峰值環(huán)向拉應(yīng)力相對較大，且隨著距離的增加，差異性更加明顯，尤其在0.3 m 外，如圖23(b)所示。

圖23 不同距離峰值環(huán)向應(yīng)力變化Fig. 23 Peak hoop stress variation at various distances

3.3 裂紋擴展行為

各向同性壓力下巖體爆破裂紋擴展特征如圖24 所示。徑向主裂紋在各方向分布勻稱，大致呈圓形裂紋帶狀分布，如圖24(a)所示。其擴展半徑隨壓力的增加而不斷減小，且減小速率逐漸降低，當(dāng)圍壓增加至40 MPa 左右，裂隙擴展半徑幾乎不再發(fā)生明顯變化，此外，徑向主裂紋的數(shù)量隨壓力增大也會顯著減少，與裂紋擴展半徑有相似的變化規(guī)律，如圖24(b)所示。經(jīng)過指數(shù)函數(shù)擬合，可獲得圓形裂紋帶擴展半徑Lc(m)與圍壓p(MPa)之間的關(guān)系：

圖24 各向同性壓力下爆破裂紋擴展特征Fig. 24 Characteristics of blasting crack propagation under different equibiaxial pressures

徑向裂紋的擴展過程如圖24(c)所示。在0.1 ms 內(nèi)，不同壓力下裂紋擴展速度幾乎相同，這是因為爆炸初期的爆轟壓力遠大于初始應(yīng)力，初始應(yīng)力對巖體爆破裂紋擴展的抑制作用被弱化。而0.1 ms 后，裂紋的擴展速度發(fā)生了明顯的變化，隨初始壓力的增加其擴展速度明顯降低，且隨著時間的推移，這種趨勢愈加明顯。初始壓應(yīng)力對裂紋的擴展具有明顯的抑制作用，深部巖體在爆破時可能會因為高地應(yīng)力擠壓作用導(dǎo)致破碎困難，此時建議結(jié)合預(yù)裂爆破等手段進行卸壓以削弱地應(yīng)力的抑制作用，從而提高爆破破巖效率。

各向異性壓力條件下巖體爆破裂紋擴展特征如圖25 所示。爆破裂紋在水平及豎直方向大致對稱分布，徑向主裂紋呈橢圓形裂紋帶擴展，如圖25(a)所示。橢圓裂紋帶橫軸長度Lx與豎軸長度Ly的比值定義為橢圓裂紋帶軸長比c，即c=Lx/Ly。隨側(cè)壓系數(shù)的增大，橫軸裂紋長度無明顯變化，豎軸裂紋長度顯著減小，不同方向主裂紋長度的變化規(guī)律與初始應(yīng)力具有良好的一致性。此外，橢圓裂紋帶軸長比會隨側(cè)壓系數(shù)的增加而增大，且增長的速率也逐漸提高，如圖25(b)所示，也就是說，當(dāng)初始水平與豎直應(yīng)力差增大時，主裂紋帶分布的各向異性會更加顯著。經(jīng)過指數(shù)函數(shù)擬合，可獲得橢圓裂紋帶的軸長比c與側(cè)壓系數(shù)K的關(guān)系：

徑向裂紋在水平與豎直方向的擴展過程如圖25(c)所示，其擴展規(guī)律與各向同性壓力下類似，在此不再贅述。需要說明的是，不管側(cè)壓系數(shù)如何變化，徑向主裂紋總是會傾向于最大主應(yīng)力方向擴展。初始應(yīng)力對裂紋擴展長度抑制的同時，對其擴展的方向也有誘導(dǎo)作用，因此，在實際工程中，若要對巖體進行爆破卸壓，炮孔應(yīng)盡可能沿最大主應(yīng)力方向布置，以便于促進炮孔間裂隙的貫通。

3.4 損傷分形維數(shù)

巖石損傷破裂過程具有明顯的分形特征，因此，巖石的損傷演化可以看作是分形的演化過程[29]。同樣，爆破裂紋的擴展過程也可以用分形維數(shù)來表征。本文采用計盒分形維數(shù)法定量描述初始應(yīng)力下爆破裂紋的分布形態(tài)，分形維數(shù)D的數(shù)學(xué)表達式為[30]

式中：δi為方盒的邊長，N(δi)為覆蓋爆破裂紋的方盒數(shù)量。

將爆破裂紋損傷圖像轉(zhuǎn)換成黑白二值圖像，然后將二值圖像劃分成若干個邊長為δi的方盒，最后再統(tǒng)計覆蓋爆破裂紋的所有盒子數(shù)量N(δi)，如圖26 所示。當(dāng)方盒長度δi趨于0 時，方盒數(shù)量N(δi)與方盒長度δk分別取對數(shù)后的比值將趨于Db，對數(shù)據(jù)點用最小二乘法擬合，其線性方程為：

圖26 爆破裂紋盒子劃分Fig. 26 Box division for blasting cracks

式中：線性方程斜率Db為爆破裂紋的分形維數(shù)，S為線性擬合曲線的縱軸截距。

此外，不同初始應(yīng)力下巖體爆破破裂的分形損傷ω 可用公式ω=(D-D0)/(Dmax-D0)來量化比較，其中，D0為巖石爆破前的初始分形維數(shù)，Dmax為巖石完全損傷時的最大分形維數(shù)。對于完整巖石而言，初始分形維數(shù)D0=0。同時，本文理論模型及數(shù)值模型是基于二維平面應(yīng)變條件下進行的，巖石破裂也是二維平面的破裂，因此最大分形維數(shù)Dmax=2。也就是說，初始應(yīng)力下巖體爆破破裂的分形損傷ω 與分形維數(shù)Db呈線性相關(guān)，因此僅選擇一個分形參量的變化規(guī)律分析即可。

不同初始應(yīng)力作用下巖體爆破裂紋分形維數(shù)擬合線如圖27 和圖28 所示，由圖可知，所有的分形維數(shù)數(shù)據(jù)都可以用一條直線很好地進行擬合。各向同性壓力下，隨著圍壓的不斷增加，爆破裂紋的分形維數(shù)不斷減小，初始壓力為0 時分形維數(shù)最大：Db=1.403，當(dāng)圍壓增長至10、20 和30 MPa 時，對應(yīng)的爆破裂紋分形維數(shù)分別為1.341、1.214 和1.087，如圖27 所示。各向異性壓力下，側(cè)壓系數(shù)K為0.25、0.5、1.5 和2.0 時，對應(yīng)的爆破裂紋分形維數(shù)分別為1.384、1.247、1.131 和1.078，隨著側(cè)壓系數(shù)的增加，分形維數(shù)呈現(xiàn)降低的變化趨勢，如圖28 所示。

圖27 各向同性壓力下爆破裂紋分形維數(shù)擬合線Fig. 27 Fractal dimension of blasting cracks and its fitting lines under equibiaxial pressure

圖28 各向異性壓力下爆破裂紋分形維數(shù)擬合線Fig. 28 Fractal dimension of blasting cracks and its fitting lines under anisotropic pressure

4 結(jié) 論

本文基于彈性力學(xué)建立了初始應(yīng)力下單孔爆破動靜組合理論模型，研究了靜態(tài)應(yīng)力分布及動態(tài)應(yīng)力演化特征，揭示了巖體爆破破裂機制。通過相應(yīng)的數(shù)值模擬，討論了初始應(yīng)力對巖體爆破壓力演化及裂紋擴展行為的影響，并引入分形維數(shù)定量表征巖體爆破損傷特性，結(jié)論如下：

(1) 動靜組合理論模型可以較好地反應(yīng)初始應(yīng)力下巖體爆破時的靜態(tài)應(yīng)力分布及動態(tài)應(yīng)力演化過程，解釋其爆破碎裂機理，并在相應(yīng)數(shù)值模型中得到驗證；

(2) 初始應(yīng)力下巖體爆破裂紋的擴展行為可以通過環(huán)向應(yīng)力分析；各向同性壓力下環(huán)向應(yīng)力分布完全對稱，各向異性壓力下環(huán)向應(yīng)力主要集中于最大主應(yīng)力方向；初始應(yīng)力對裂紋擴展范圍抑制的同時，對裂紋擴展方向也具有導(dǎo)向作用，隨著水平與豎直初始應(yīng)力差的增加，裂紋分布差異性顯著增強；合理調(diào)整環(huán)向應(yīng)力分布，可以改善巖體爆破碎裂效果；

(3) 各向同性壓力下爆破裂紋帶近似呈圓形，裂紋擴展半徑及數(shù)量隨壓力的增加顯著減?。桓飨虍愋詨毫ο卤屏鸭y帶近似呈橢圓形，裂紋帶的軸長比隨側(cè)壓系數(shù)的增加顯著增大；高地應(yīng)力下巖體爆破破碎效果不佳時，建議結(jié)合預(yù)裂技術(shù)先卸壓再爆破；

(4) 初始應(yīng)力下巖體爆破破裂過程具有明顯的分形特征，各向同性壓力下分形維數(shù)隨圍壓的增加而減小，各向異性壓力下呈現(xiàn)相同的變化規(guī)律，分形維數(shù)可以定量表征巖體爆破破裂的損傷特性。