初始應(yīng)力下巖體爆破損傷特性及破裂機(jī)理*

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

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

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

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

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

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

1 動(dòng)靜組合理論模型

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

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

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

根據(jù)彈性力學(xué)理論，可分別計(jì)算不同初始應(yīng)力條件下炮孔近場(chǎng)(r=a)與遠(yuǎn)場(chǎng)(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。

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

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

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

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

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

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

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

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

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

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

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

2 數(shù)值計(jì)算模型及驗(yàn)證

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

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

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

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

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

式中：Q0為拉伸與壓縮荷載下材料子午線半徑之比，B為巖石材料羅德角相關(guān)系數(shù)，參考文獻(xiàn)[23]中的試驗(yàn)數(shù)據(jù)，進(jìn)行線性回歸求得：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)變表達(dá)式如下：

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

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

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

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

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

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

2.2 炸藥及空氣材料

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

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

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

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

式中：pa為氣體壓力，為單位氣體體積的內(nèi)部能量，u為動(dòng)態(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ù)值模型驗(yàn)證

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

根據(jù)監(jiān)測(cè)點(diǎn)記錄的數(shù)據(jù)，繪制相應(yīng)的環(huán)向應(yīng)力演化曲線，如圖22 所示。各向同性壓力下水平與豎直方向監(jiān)測(cè)點(diǎn)是中心對(duì)稱的，因此其環(huán)向應(yīng)力時(shí)程曲線基本重合，如圖22(a)所示。而各向異性壓力下水平與豎直方向上的環(huán)向應(yīng)力時(shí)程曲線有所區(qū)別，主要表現(xiàn)在拉壓應(yīng)力峰值的差異性，最大主應(yīng)力方向上的壓應(yīng)力削弱，而拉應(yīng)力增強(qiáng)，如圖22(b)所示。這可以很好地解釋裂紋傾向于最大主應(yīng)力方向擴(kuò)展的原因。值得注意的是，不同初始?jí)毫ο?，其環(huán)向應(yīng)力的變化規(guī)律具有相似性，即爆炸初期因其強(qiáng)烈的沖擊作用產(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)力時(shí)程曲線Fig. 22 Time histories of hoop stress under various pressure conditions

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

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

3.3 裂紋擴(kuò)展行為

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

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

徑向裂紋的擴(kuò)展過(guò)程如圖24(c)所示。在0.1 ms 內(nèi)，不同壓力下裂紋擴(kuò)展速度幾乎相同，這是因?yàn)楸ǔ跗诘谋Z壓力遠(yuǎn)大于初始應(yīng)力，初始應(yīng)力對(duì)巖體爆破裂紋擴(kuò)展的抑制作用被弱化。而0.1 ms 后，裂紋的擴(kuò)展速度發(fā)生了明顯的變化，隨初始?jí)毫Φ脑黾悠鋽U(kuò)展速度明顯降低，且隨著時(shí)間的推移，這種趨勢(shì)愈加明顯。初始?jí)簯?yīng)力對(duì)裂紋的擴(kuò)展具有明顯的抑制作用，深部巖體在爆破時(shí)可能會(huì)因?yàn)楦叩貞?yīng)力擠壓作用導(dǎo)致破碎困難，此時(shí)建議結(jié)合預(yù)裂爆破等手段進(jìn)行卸壓以削弱地應(yīng)力的抑制作用，從而提高爆破破巖效率。

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

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

3.4 損傷分形維數(shù)

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

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

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

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

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

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

不同初始應(yīng)力作用下巖體爆破裂紋分形維數(shù)擬合線如圖27 和圖28 所示，由圖可知，所有的分形維數(shù)數(shù)據(jù)都可以用一條直線很好地進(jìn)行擬合。各向同性壓力下，隨著圍壓的不斷增加，爆破裂紋的分形維數(shù)不斷減小，初始?jí)毫? 時(shí)分形維數(shù)最大：Db=1.403，當(dāng)圍壓增長(zhǎng)至10、20 和30 MPa 時(shí)，對(duì)應(yīng)的爆破裂紋分形維數(shù)分別為1.341、1.214 和1.087，如圖27 所示。各向異性壓力下，側(cè)壓系數(shù)K為0.25、0.5、1.5 和2.0 時(shí)，對(duì)應(yīng)的爆破裂紋分形維數(shù)分別為1.384、1.247、1.131 和1.078，隨著側(cè)壓系數(shù)的增加，分形維數(shù)呈現(xiàn)降低的變化趨勢(shì)，如圖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)力下單孔爆破動(dòng)靜組合理論模型，研究了靜態(tài)應(yīng)力分布及動(dòng)態(tài)應(yīng)力演化特征，揭示了巖體爆破破裂機(jī)制。通過(guò)相應(yīng)的數(shù)值模擬，討論了初始應(yīng)力對(duì)巖體爆破壓力演化及裂紋擴(kuò)展行為的影響，并引入分形維數(shù)定量表征巖體爆破損傷特性，結(jié)論如下：

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

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

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

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