初始應力下巖體爆破損傷特性及破裂機理*

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

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

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

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

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

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

1 動靜組合理論模型

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

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

1.1 靜態應力分布特征

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

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

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

圖2 炮孔周邊靜態應力分布Fig. 2 Stress distribution around the borehole under static load

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

圖3 不同位置處巖體靜態應力的變化Fig. 3 Stress variation of rock mass at different distances under static load

1.2 動態應力演化規律

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

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

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

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

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

圖4 不同距離處巖體的動態應力波演化Fig. 4 Dynamic stress waves evolution in rock mass at various distances

2 數值計算模型及驗證

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

2.1 RHT 模型參數標定

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

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

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

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

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

圖5 RHT 模型破壞面的相關參數擬合Fig. 5 Fitting curve of failure surface parameters for RHT model

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

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

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

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

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

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

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

2.2 炸藥及空氣材料

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

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

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

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

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

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

2.3 數值模型驗證

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

圖8 數值模型及局部網格Fig. 8 Configuration of numerical model and local mesh

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

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

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

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

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

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

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

3 單孔爆破數值模擬

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

表5 初始應力加載條件Table 5 Initial stress conditions in numerical simulation

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

3.1 靜態應力初始化

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

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

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

圖13 應力初始化驗證Fig. 13 Verification of stress initialization

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

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

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

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

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

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

3.2 壓力演化過程

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

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

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

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

圖22 不同壓力條件下環向應力時程曲線Fig. 22 Time histories of hoop stress under various pressure conditions

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

圖23 不同距離峰值環向應力變化Fig. 23 Peak hoop stress variation at various distances

3.3 裂紋擴展行為

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

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

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

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

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

3.4 損傷分形維數

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

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

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

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

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

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

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

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

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

4 結 論

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

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

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

(3) 各向同性壓力下爆破裂紋帶近似呈圓形，裂紋擴展半徑及數量隨壓力的增加顯著減??；各向異性壓力下爆破裂紋帶近似呈橢圓形，裂紋帶的軸長比隨側壓系數的增加顯著增大；高地應力下巖體爆破破碎效果不佳時，建議結合預裂技術先卸壓再爆破；

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