李 斌，許夢國，王文杰

(1.武漢科技大學，湖北 武漢 430081；2.冶金礦產(chǎn)資源高效利用與造塊湖北省重點實驗室，湖北 武漢 430081)

無底柱分段崩落法是20世紀50年代以來隨著新的采礦工藝和設備的出現(xiàn)而逐步發(fā)展起來的的一種采礦方法，由于其高效率、高度機械化、工藝簡單、生產(chǎn)安全和成本低廉的特點，無底柱分段崩落法在金屬礦山得到迅速推廣，特別是鐵礦山更為廣泛。無底柱分段崩落法的采礦機理主要是礦石和圍巖的爆破崩落[1]。然而，無底柱分段崩落法在實際爆破落礦過程中會產(chǎn)生較多的大塊，不僅增加了成本、降低了采礦效率，還影響了采礦的安全生產(chǎn)。

爆破產(chǎn)生不同塊度的巖塊可以認為是爆破產(chǎn)生各種裂隙的結(jié)果，研究爆破裂隙的產(chǎn)生和傳播規(guī)律可以幫助更好的了解大塊的產(chǎn)生機理。爆破的化學反應過程是炸藥從具有高能量的密致物質(zhì)轉(zhuǎn)化為高壓、高溫的氣體產(chǎn)物[2-3]。近年來通過對巖石爆破的不斷研究，國內(nèi)外專家學者一致認為在巖石爆破過程中，有兩類荷載作用在巖石上：爆炸應力波作用和爆生氣體膨脹作用[4]。然而，很缺乏對爆破過程中巖石裂隙產(chǎn)生和發(fā)展規(guī)律的認識，主要原因是爆破過程非常迅速，一般在幾個毫秒就結(jié)束了，很難對爆破過程進行深入的了解。

為了對無底柱分段崩落法爆破過程進行比較深入的研究，采用非線性顯示動力學軟件AUTODYN進行爆破模擬，來分析無底柱分段崩落法爆破落礦過程中爆破裂隙的產(chǎn)生和傳播規(guī)律。

1 數(shù)值模擬計算模型

在本文中，采用AUTODYN 2D平面對稱建立計算模型，假定巖體是連續(xù)、各向同性的。本次模擬的中深孔布置及參數(shù)如圖1所示，炮孔孔徑φ為80mm。

圖1 中深孔布置

1.1 材料模型

在本文模擬中用到的材料有兩種：巖石(鐵礦石)和炸藥(ANFO)。

1.1.1 巖石動力學模型

本次模擬某鐵礦無底柱分段崩落法中深孔爆破落礦，鐵礦石可以認為是脆性巖石，這類脆性材料具有高抗壓強度和低抗拉強度的性質(zhì)，當受到荷載作用時因為微裂隙的擴展而表現(xiàn)出漸進式破壞。針對這種連續(xù)介質(zhì)破壞模型，Johnson和Holmquist[5-7]提出了一種適用于脆性材料的本構(gòu)模型(J-H模型)。J-H模型可以用來研究脆性材料(如巖石、混凝土和陶瓷等)沖擊和爆炸問題。J-H本構(gòu)模型包括有強度模型和失效模型，在使用J-H本構(gòu)模型時，必須同時使用其強度模型和失效模型[8]。

具體而言，J-H模型有兩種形式：JH-1[9]模型和JH-2[7]模型。JH-1模型主要是用來計算大變形而并沒有考慮由于不斷增加的變形導致的漸進式破壞。JH-2模型通過添加多節(jié)段應力-應變曲線來彌補JH-1模型的不足。

1.1.1.1 多項式狀態(tài)方程

本文采用多項式狀態(tài)方程來描述巖石，通常情況下多項式狀態(tài)方程是和J-H模型一起使用的[10]。巖石在失效前，多項式狀態(tài)方程可以表示為式(1)。

P=K1μ+K2μ2+K3μ3+ΔP

(1)

式中：P為壓力；μ為壓縮量，μ=V0/V-1=ρ/ρ0-1，V和ρ代表巖石當前的體積和密度，V0和ρ0代表巖石初始的體積和密度；K1，K2，K3為常數(shù)(K1通常是體積模量)。

圖2所示為壓力P和壓縮量μ的關系圖，式(1)中的ΔP就是圖中的壓力增量，ΔP取決于材料的破壞程度，它的值變化范圍由D=0時的ΔP=0到D=1時的ΔPmax。D是破壞程度，取值范圍為0～1。而對于拉伸應力(μ<0)，多項式狀態(tài)方程表示為式(2)。

P=K1μ

(2)

圖2 壓力P和壓縮量μ關系

1.1.1.2 JH-2本構(gòu)模型

JH-2模型假定材料最初是彈性和各向同性的，滿足式(3)。

σij=-p(εkk)δij+2Gξij

(3)

Hugoniot彈性極限時的應力可以有如式(2)所示的關系式。

(4)

式中：pHEL為Hugoniot彈性極限的壓力；σHEL為Hugoniot彈性極限的應力。

JH-2模型的強度是服從塑性力學中的米塞斯條件的[10]。即

(5)

式中σ1,σ2,σ3為主應力。

有效應力通?？梢杂脕碜鳛槊枋鰪椥詷O限破壞的標準。對于服從米塞斯條件的模型，有效應力可以表示為式(6)。

(6)

完整材料的破壞強度表達式為式(7)。

(7)

(8)

材料完全破壞時的強度表達式為式(9)。

(9)

式中B，C，M為材料常數(shù)。

材料在破壞過程中的強度表達式為式(10)。

(10)

式中D為當前的破壞程度。

破環(huán)程度隨時間變化的關系式為式(11)。

(11)

式中等效塑性應變εf的表達式為式(12)。

εf=D1(p*+T*)D2

(12)

式中D1和D2為材料的破壞常數(shù)。

1.1.2 確定JH-2模型材料參數(shù)

巖石材料的強度一般是應用Koek-Brown失效標準來評估，尤其是像鐵礦石這類堅硬巖石[11]。作為一個經(jīng)驗公式，Koek-Brown失效標準的表達式[12]為式(13)。

(13)

式中：σ1和σ3為三軸試驗中的軸向應力和側(cè)向應力；σc為單軸抗壓強度；m和s為與材料有關的常數(shù)，本文模擬的是完整巖石材料，s取1.0，根據(jù)文獻[12]，m取32。

在研究動力學問題時通常需要采用動態(tài)強度，動態(tài)強度可以超過相對應靜態(tài)強度一個甚至數(shù)個數(shù)量級[13-14]，在本次模擬中選取動態(tài)強度是靜態(tài)強度的10倍?？梢缘玫綆r石的基本參數(shù)如表1所示。

表1 巖石的基本參數(shù)

巖石材料常數(shù)A和N根據(jù)文獻[11]按照式(7)進行曲線擬合，得到的數(shù)據(jù)分別為0.86和0.63。考慮到難測取巖石完全破壞時的數(shù)據(jù)，這里假定殘余強度為完整巖石強度的30%[15]，得到B和M的值分別為0.29和0.63。參照其他脆性材料的破壞參數(shù)，D1和D2選取為0.006和0.8。應變率常數(shù)C是隨研究材料而變化的，Johnson GR研究的玻璃的C值為0.003[16]，Holmquist TJ研究的碳化硅的C值為0.009[17]，由于缺乏鐵礦石這方面的研究數(shù)據(jù)，參照以上數(shù)據(jù)將本次研究的C值選取為0.006。

1.1.3 炸藥模型

本次模擬采用的炸藥為ANFO。描述爆炸產(chǎn)物的狀態(tài)方程很多，然而應用最普遍的是JWL(Jones-Wilkins-Lee) EOS[18]，JWL EOS是描述高能炸藥爆轟產(chǎn)物壓力-體積關系的狀態(tài)方程，JWL狀態(tài)方程的P-V關系如式(14)所示。

(14)

式中：P為壓力；A、B、R1、R2、ω為炸藥常數(shù)；V是爆炸產(chǎn)物的相對體積；E0為初始內(nèi)能密度。

本次研究采用的ANFO炸藥各參數(shù)如表2所示。

表2 ANFO炸藥參數(shù)

1.2 計算模型

ANSYS AUTODYN在研究巖石破壞方面有了成功的應用案例[19-20]。在本次研究中，建立的計算模型如圖3所示。

圖3中綠色區(qū)域代表巖石，藍色條狀代表ANFO炸藥，紅色標記點為炸藥起爆點(孔底起爆)，11個炮孔同時起爆，數(shù)字標記點為高斯點，邊界條件為無反射邊界，邊界材料為巖石。模型尺寸為24000mm×42000mm，考慮到炮孔直徑為80mm，以及計算機性能，按40mm×40mm劃分網(wǎng)格，總共有630000個單元。由于模型是關于炮孔6對稱的，故只在炮孔6～11之間布置高斯點，這些高斯點設置在孔底、孔口和中間位置，分別用來研究孔底、孔口和中間部位的爆破過程。為便于后面的數(shù)據(jù)記錄和分析，每個部位的高斯點縱坐標一樣，為盡量獲取較多的數(shù)據(jù)，橫坐標間距40mm。

2 模擬結(jié)果

基于以上計算模型、狀態(tài)方程、JH-2本構(gòu)模型和相關參數(shù)，在AUTODYN中進行計算，計算到3475次時爆破已到孔口位置，可以認為爆破過程基本完成了，得到的材料狀態(tài)如圖4所示。

從圖4中可以明顯觀察到：孔口位置產(chǎn)生非常密集的爆破裂隙，中間位置爆破裂隙也比較密集，而在靠近孔底位置爆破裂隙顯得很稀疏，特別是中間5個炮孔之間。可以認為巖石在孔口位置破壞程度比較大，而在孔底位置破壞程度較小。

在應力波作用下巖石的破碎過程主要由兩方面決定：①應力波產(chǎn)生的最大應力；②應力波的比沖能和比能[21]。這里主要從最大爆破應力來分析研究扇形中深孔的爆破過程。

從圖3可以看出,炮孔布置是以炮孔6對稱的，考慮到炮孔6、7、炮孔7、8、炮孔8、9和炮孔10、11相鄰孔布置形式一樣，就選取炮孔7、8之間的巖石作為研究對象，以及布孔形式有差異的炮孔9、10。

2.1 孔底爆破過程分析

統(tǒng)計孔底位置高斯點的最大爆破應力，如圖5所示。

圖5 孔底最大爆破應力統(tǒng)計

圖5中(a)、(b)分別是炮孔7和炮孔8之間孔底，炮孔9和炮孔10之間孔底最大爆破力統(tǒng)計圖。圖5中可以觀察到，(a)曲線呈現(xiàn)兩邊低中間高的現(xiàn)象，平均最大爆破應力值為：3.2385GPa，而(b)曲線卻不同，呈現(xiàn)兩邊高中間低的現(xiàn)象，平均最大爆破應力為1.2062GPa。

炮孔8孔底到炮孔7的水平距離為2200mm，布置高斯點57個。圖6中3號高斯點靠近炮孔7，56號高斯點靠近炮孔8，28號高斯點位于這段研究距離的中間位置。從56號高斯點應力曲線圖可以看出炮孔8在0ms時爆破到這個位置，爆破應力急劇上升并達到最大值；從3號高斯點應力曲線圖炮孔7在0.5ms時爆破到這個位置，爆破應力急劇上升并達到最大值；從28號高斯點應力曲線可以看出在最大應力出現(xiàn)前有一段平穩(wěn)應力值較小的曲線，這是炮孔8爆破應力波到達這個位置造成的，隨后炮孔7爆破應力波達到，而后應力急劇上升并到達最大值，該值將近是3號和56號高斯點最大爆破應力的3倍。從而出現(xiàn)了炮孔中間位置的最大爆破應力比靠近炮孔位置高的現(xiàn)象。

圖6 爆破應力曲線(3、28和56號高斯點)

炮孔10孔底到炮孔9的水平距離為1860mm，布置高斯點48個。圖7中2號高斯點靠近炮孔9，47號高斯點靠近炮孔10，24號高斯點位于這段研究距離的中間位置。炮孔9和炮孔10都在0ms時爆破到這個位置，2號、47號高斯點爆破應力急劇上升并達到最大值；24號高斯點處炮孔9和爆破產(chǎn)生的應力波幾乎同時到達，產(chǎn)生一系列的波峰和波谷，但是疊加后的爆破應力卻不到2號和47號高斯點最大爆破應力的一半。從而出現(xiàn)了中間曲線值較小的現(xiàn)象。

圖7 爆破應力曲線(2、24和47號高斯點)

通常情況下，孔間距越小爆破效果越好，而這里模擬的結(jié)果發(fā)現(xiàn)2200mm孔間距的平均最大爆破應力達到1860mm孔間距的2.5倍，造成這種情況的原因就是起爆時間的間隔。炮孔9和炮孔10幾乎同時起爆，其孔底位置有很好的應力波疊加或碰撞，然而最大應力卻偏??；炮孔8提前炮孔7約0.5ms爆破，炮孔8爆破應力波優(yōu)先在孔底位置傳播，其作用降低了這部分巖石的強度和完整性，當隨后炮孔7爆破產(chǎn)生的應力波到達時，巖石得到很大程度的破壞。Liqing Liu和P.D.Katsabanis的研究[22]提出：相鄰炮孔同時起爆時應力波的疊加和碰撞會加強對巖石的破壞能力，當起爆間隔超過50μs時,最大爆破應力相比降低了46%。本文的研究結(jié)果卻截然相反，在這里對后者提出質(zhì)疑。

2.2 中間位置爆破過程分析

統(tǒng)計炮孔中間位置高斯點的最大爆破應力，如圖8所示。

圖8中(a)、(b)分別是炮孔7和炮孔8之間中間位置、炮孔9和炮孔10之間中間位置最大爆破力統(tǒng)計圖?？梢杂^察到，圖8中(a)曲線與圖5中(b)曲線非常相似，但是中間曲線更平穩(wěn)，平均最大爆破應力值為2.5700GPa；而圖8中(b)曲線呈現(xiàn)兩邊向中間遞減的現(xiàn)象，平均最大爆破應力為2.8741GPa。

炮孔7和炮孔8中間某部位的水平距離為1440mm，共計37個高斯點。圖9中58號高斯點靠近炮孔7，91號高斯點靠近排孔8，75號高斯點位于這段研究距離的中間位置。炮孔8提前炮孔7約0.4ms起爆，從91號高斯點應力曲線可以觀察到，隨著炮孔8的爆破，靠近炮孔8的地方爆破應力急劇上升并達到最大值；從58號高斯點曲線可以看出，當炮孔8爆破應力波達到后炮孔7才開始爆破，然后應力急劇上升達到最大值；從75號高斯點曲線可以看出，在0.8ms左右炮孔8產(chǎn)生的應力波到達這位置，應力峰值達到1.8GPa，在1.1ms左右炮孔7產(chǎn)生的應力波到達，達到最大應力2.2GPa，增加量很小。對比圖6可以發(fā)現(xiàn)，隨著炮孔間距的減小，最大爆破應力卻呈現(xiàn)降低的現(xiàn)象，這是因為炮孔間距的降低，當炮孔8產(chǎn)生的爆破應力波作用巖石時，產(chǎn)生較大的應力使巖石產(chǎn)生一定程度的破壞(包括爆破裂隙和自由面)，從而降低了炮孔7應力波的作用。

圖8 孔底最大爆破應力統(tǒng)計

圖9 爆破應力曲線(58、75和91號高斯點)

炮孔9和炮孔10中間某部位的水平距離為1270mm，共計33個高斯點。圖10中80號高斯點靠近炮孔9，78號高斯點靠近炮孔10，64號高斯點位于這段研究距離的中間位置。由于炮孔傾斜角度的影響在這個水平產(chǎn)生了大約0.1ms的爆破間隔，可以發(fā)現(xiàn)圖10 與圖9中的曲線非常相似，主要不同是圖10中64號高斯點的最大爆破應力(2.3068GPa)出現(xiàn)在第一個波峰,接著第二個應力波到達造成第二個波峰(爆破應力約為1.9GPa)。

圖10 爆破應力曲線(50、64和78號高斯點)

通過以上分析可以得出：當相鄰炮孔的間距減小到某一程度時，由于起爆之間的間隔，第一個應力波在巖石中傳播時就會對巖石產(chǎn)生一定程度的破壞，降低了第二個應力波對巖石的作用，并且隨著距離的減小，第一個應力波產(chǎn)生的爆破應力越大，第二個應力波產(chǎn)生的應力波就越小。

2.3 孔口爆破過程分析

統(tǒng)計炮孔孔口位置高斯點的最大爆破應力，如圖11所示。

圖11 孔口最大爆破應力統(tǒng)計圖

圖11中(a)、(b)分別是炮孔7和炮孔8之間孔口、炮孔9和炮孔10之間孔口最大爆破力統(tǒng)計圖，平均最大應力分別為8.5367GPa和7.3526GPa，兩條曲線非常相似，呈現(xiàn)U形。

炮孔8孔口到炮孔7的水平距離為660mm，共計17個高斯點。圖12中95、97號高斯點靠近炮孔7，107、109號高斯點靠近炮孔8，從圖12中可以看出，107、109號高斯點的最大爆破應力是由炮孔8爆破應力波造成的，95、97號高斯點的最大爆破應力是由炮孔7爆破應力波造成的，并且里炮孔越近最大爆破應力越大，從而出現(xiàn)了圖11中的U形曲線。炮孔9和炮孔10之間的孔口位置情況與此類似，這里就不再贅述。

圖12 爆破應力曲線(95、97、107和109號高斯點)

3 討論

在孔底位置，由于炮孔間距比較大，單個應力波作用巖石產(chǎn)生的最大爆破應力較小，若相鄰炮孔同時起爆，兩個應力波疊加和碰撞作用可以使中間部分巖石的最大爆破應力達到靠近炮孔巖石最大爆破應力的一半左右；而當相鄰炮孔間隔起爆時，由于前一個炮孔爆破應力波(小于巖石的靜態(tài)單軸抗壓強度)作用，降低了中間巖石的強度及完整性，當?shù)诙€應力波達到時，巖石受到的爆破應力大幅度上升并到達最大值，其值可以達到靠近炮孔巖石最大爆破應力的3倍。在前期的現(xiàn)場分析中[23]，同條進路相鄰5排炮孔，相同的爆破工藝，前3排巖石的完整性比較好，爆破產(chǎn)生較多的大塊，而后2排巖石受到前排爆破的影響，完整性較差、頂板比較破碎，爆破效果很好，基本上沒有大塊的出現(xiàn)。這說明完整性很好的巖石受到爆破應力波作用降低強度后，當再次有爆破應力波作用時可以很好的破碎巖石。在文獻[24]中，大塊主要出現(xiàn)在孔底位置，可見改善相鄰炮孔孔底之間的最大爆破應力非常重要。

當相鄰炮孔間距減小時，第一個應力波作用巖石的應力較大(超過巖石的靜態(tài)抗壓強度)時，巖石會產(chǎn)生一定程度的破壞，降低了第二個應力波的作用，并且當?shù)谝粋€應力波產(chǎn)生的應力越大時，第二個應力波的作用越小。

在孔口位置，由于相鄰炮孔間距很小，巖石的最大爆破應力主要是由距離近的炮孔爆破應力波作用造成，距離炮孔越近，最大爆破應力越大。

4 結(jié)論與展望

本文應用AUTODYN進行扇形孔爆破模擬，主要通過分析炮孔之間的最大爆破應力產(chǎn)生機理來解釋巖石爆破裂隙的分布狀態(tài)，主要有以下結(jié)論。

1) 炮孔越是密集的位置，爆破應力越大，爆破裂隙也越多，從圖4中可以清晰看到這一點，這就意味著孔底位置可能有大塊的出現(xiàn)。

2) 相鄰炮孔間距較大時，尤其是孔底位置，應避免同時(或間隔時間很短)起爆，而應該確定合理的間隔時間，進行間隔起爆。

3) 相鄰炮孔間距較小時，間隔起爆的作用就不明顯了，巖石受到的最大爆破應力主要受距離近的炮孔爆破應力波的影響。

為改善扇形孔的爆破質(zhì)量，下一步的工作就是確定合理的炮孔間距范圍，以及對應不同孔間距的起爆時間間隔。

[1] D.Brunton,S.J.Fraser.Parameters influencing full scale sublevel caving material recovery at the Ridgeway gold mine[J].International Journal of Rock and Mining Sciences,Volume 47,Issue 4,June 2010：647-656.

[2] Bhandari S.Engineering rock blasting operations[M].Brookfield;A.A.Balkema,1997.

[3] Fourney WL.Mechanisms of rock fragmentation by blasting.Comprehensive rock engineering[M].Vol.4,ed J.A.Hudson,(Pergamon),1993:39-69.

[4] H.K.Kutter,C.Fairhurst.On the fracture process in blasting[J].International Journal of Rock Mechanics and Mining Science,Volume 8,Issue 3,May 1971：181-188.

[5] Holmquist TJ,Templeton DW,Bishnoi KD.Constitutive modeling of aluminum nitride for large strain,high-strain rate,and high-pressure applications[J].Int J Impact Eng,2001,25:211-231.

[6] Johnson GR,Holmquist TJ,Beissel SR.Response of aluminum nitride(including a phase change) to large strains,high strain rates,and high pressures[J].J Appl Phys,2003,94(3):1639-1646.

[7] Johnson GR,Holmquist TJ.An improved computational constitutive model for brittle materials.In:Schmidt SC,Shaner JW,Samara GA,Ross M,editors.High pressure science and technology-1993[M].Woodbury,NY:AIP Press;1994.p.981-4.

[8] Century Dynamics Inc.AUTODYN Theory Manual,Revision 4.3[M].Concord,California,2005.

[9] Johnson GR,Holmquist TJ.A computational constitutive model for brittle materials subjected to large strains,shock-wave and high strain-rate phenomena in materials.Ed.M.A.Meyers,L.E.Murr and K.P.Staudhammer,Marcel Dekker lnc[M].New York,1992,PP.1075-1081.

[10] M.M.Dehghan Banadaki,B.Mohanty.Numerical simulation of stress wave induced fractures in rock[J].International Journal of Impact Engineering,40-41(2012):16-25.

[11] Hoek E,Brown ET.Practical estimates of rock mass strength[J].International Journal of Rock Mechanics and Mining Sciences,1997;34(8):1165-1186.

[12] Hoek E,Brown ET.Underground excavations in rock[M].London:Institute of Mining and Metallurgy,1980.

[13] Rinehart JS.Dynamic fracture strength of rocks.Proceedings of the seventhsymposium on rock mechanics[M].Pennsylvania State University,1965.

[14] Prasad U,Mohanty B,Nemes J.Dynamics strength and fragmentation in selected rocks under impact loading.In:Girard J,et al.,editors.Proceedings of the fourth North American rock mechanics symposium—NARMS 2000[M].Rotterdam:Balkema;2000.p.577-582.

[15] G.W.Ma,X.M.An.Numerical simulation of blasting-induced rock fractures[J].International Journal of Rock Mechanics and Mining Sciences,2008(45):966-975.

[16] Johnson GR,Holmquist TJ.Response of boron carbide subjected to large strains,high strain rates,and high pressures[J].J Appl Phys,1999(85):8060-8063.

[17] Holmquist TJ,Johnson GR.Characterization and evaluation of silicon carbide for high-velocity impact[J].J Appl Phys,2005;97:93502-1—93502-12.

[18] Lee EL,Hornig HC,Kury JW.Adiabatic expansion of high explosive detonation products.Report NO.UCRL50422.California Univ[M].United States:Lawrence Radiation Lab;1968.

[19] Zheming Zhu,Bibhu Mohanty,Heping Xie.Numerical investigation of blasting-induced crack initiation and propagation in rocks[J].International Journal of Rock Mechanics and Mining Sciences,2007(44):412-424.

[20] Chen SG,Zhao J.A study of UDEC modeling for blasting wave propagation in joint rock mass[J].International Journal of Rock Mechanics and Mining Sciences,(35)1998:93-99.

[21] 鈕強.巖石爆破機理[M].沈陽：東北工學院出版社，1990.

[22] Liu Liqing,P.D.Katsabanis.A numerical study of the effects of accurate timing on rock fragmentation[J].International Journal of Rock Mechanics and Mining Sciences,1997,34(5):817-835.

[23] 李斌，許夢國.無底柱分段崩落法大塊產(chǎn)生機理與解決措施[A].王運敏.中國采選技術十年回顧與展望[C].北京：冶金工業(yè)出版社，2012.