王 鑫，張昌鎖

(太原理工大學(xué) 礦業(yè)工程學(xué)院，太原 030024)

在巖石爆破中，巖石的破壞主要是由于炸藥的爆轟作用形成的應(yīng)力波破壞及爆生氣體的準(zhǔn)靜態(tài)作用產(chǎn)生破壞。其中大部分巖石破壞是因?yàn)楸鷼怏w滲入裂縫產(chǎn)生的“氣楔”作用[1]。因此在數(shù)值模擬研究中，較為準(zhǔn)確地描述爆轟產(chǎn)物的膨脹做工過程是為巖土爆破設(shè)計(jì)及優(yōu)化提供理論指導(dǎo)的關(guān)鍵。

JWL狀態(tài)方程能夠比較地準(zhǔn)確描述爆轟產(chǎn)物的膨脹過程中狀態(tài)參數(shù)的變化，不含化學(xué)反應(yīng)，在工程設(shè)計(jì)中應(yīng)用較廣[2]。JWL狀態(tài)方程是一種通過試驗(yàn)與計(jì)算相結(jié)合的半經(jīng)驗(yàn)狀態(tài)方程[3]，未知參數(shù)由標(biāo)準(zhǔn)圓筒試驗(yàn)法確定[4]。JWL狀態(tài)方程參數(shù)的獲得較為繁瑣，需要進(jìn)行不斷調(diào)整，并且圓筒試驗(yàn)中圓筒很快就會(huì)破裂，對(duì)于中后期的氣體膨脹過程描述不夠準(zhǔn)確[5]。因此很多學(xué)者常利用數(shù)值計(jì)算與數(shù)值模擬結(jié)合的方法來確定未知參數(shù)。王成利用基因遺傳算法根據(jù)CJ條件和Hugoniot關(guān)系式對(duì)炸藥爆轟產(chǎn)物γ律狀態(tài)方程對(duì)P-V曲線進(jìn)行擬合[6]，得到JWL狀態(tài)方程參數(shù)，并利用數(shù)值模擬證明了方法的合理性。王言金利用貝葉斯分析方法研究了炸藥的不確定參數(shù)[7]，并對(duì)狀態(tài)方程進(jìn)行了貝葉斯參數(shù)標(biāo)定，利用標(biāo)定參數(shù)對(duì)平面爆轟驅(qū)動(dòng)飛片試驗(yàn)進(jìn)行了數(shù)值模擬并取得了良好的結(jié)果。沈飛等利用圓筒中實(shí)驗(yàn)中能量守恒關(guān)系及圓筒壁特征位置膨脹速度經(jīng)驗(yàn)計(jì)算公式[8]，建立了一種確定JWL狀態(tài)方程參數(shù)的簡(jiǎn)便算法，確定了四種常用炸藥的JWL狀態(tài)方程參數(shù)，具有良好的擬合效果。李科斌利用水下爆炸實(shí)驗(yàn)法結(jié)合實(shí)驗(yàn)結(jié)果[9]，利用Autodyna程序測(cè)定了銨油炸藥的JWL狀態(tài)方程參數(shù)。目前，狀態(tài)方程參數(shù)標(biāo)定只能通過少數(shù)實(shí)驗(yàn)：圓筒膨脹實(shí)驗(yàn)、水下爆炸實(shí)驗(yàn)、爆轟驅(qū)動(dòng)平板實(shí)驗(yàn)及地面掛高爆轟實(shí)驗(yàn)，且一種實(shí)驗(yàn)標(biāo)定的狀態(tài)方程不能用另一種實(shí)驗(yàn)來驗(yàn)證。很少文獻(xiàn)通過兩種不同的爆轟實(shí)驗(yàn)來驗(yàn)證參數(shù)的準(zhǔn)確性，對(duì)于JWL狀態(tài)方程參數(shù)調(diào)整沒有依據(jù)。工程爆破中對(duì)于炸藥要求具有一定的猛度、爆力，但在大多數(shù)爆破設(shè)計(jì)數(shù)值模擬過程中，缺少對(duì)于炸藥材料爆力、猛度指標(biāo)的驗(yàn)證。

本文基于ANSYS/LS-DYNA有限元軟件，針對(duì)工程中常用的CHNO型炸藥，使用一種簡(jiǎn)便算法得出JWL狀態(tài)方程的未知參數(shù)，將該參數(shù)作為原始參數(shù)。以工業(yè)炸藥的爆力、猛度性能指標(biāo)為標(biāo)準(zhǔn)，對(duì)鉛柱壓縮試驗(yàn)、鉛鑄擴(kuò)孔試驗(yàn)進(jìn)行數(shù)值模擬，通過兩種不同的實(shí)驗(yàn)來檢驗(yàn)計(jì)算參數(shù)的準(zhǔn)確性并根據(jù)數(shù)值實(shí)驗(yàn)結(jié)果對(duì)于計(jì)算參數(shù)進(jìn)行調(diào)整修正。為巖石爆破數(shù)值模擬研究提供一種快速確定工業(yè)炸藥參數(shù)的方法。

1 JWL狀態(tài)方程參數(shù)計(jì)算

JWL狀態(tài)方程的標(biāo)準(zhǔn)形式及等熵形式為[10]

(1)

PS=Ae-R1V+Be-R2V+CV-(ω+1)

(2)

式中：P爆轟產(chǎn)物壓力；V為爆轟產(chǎn)物的相對(duì)體積；E為爆轟產(chǎn)物單位體積的內(nèi)能。式中共有A、B、C、R1、R2、ω六個(gè)未知參數(shù)。

炸藥在爆轟的過程中滿足CJ爆轟條件，可以得到一組方程式，代表CJ爆轟條件對(duì)于JWL狀態(tài)方程中未知參數(shù)的約束

(3)

(4)

(5)

式中：VCJ是CJ點(diǎn)處爆轟產(chǎn)物的相對(duì)體積；ρ0是炸藥的初始密度；D是炸藥的爆速；PCJ是CJ點(diǎn)處的爆轟產(chǎn)物壓力；ECJ為CJ點(diǎn)的爆轟產(chǎn)物的內(nèi)能。

對(duì)于凝聚態(tài)炸藥，在CJ點(diǎn)處有

(6)

(7)

(8)

式中：ρJ為CJ點(diǎn)處爆轟產(chǎn)物的密度；r是爆轟產(chǎn)物的多方指數(shù)，與炸藥性質(zhì)有關(guān)，對(duì)于高密度的凝聚態(tài)炸藥r通常取3[11]。文獻(xiàn)[12]提出一種適用于混合炸藥的經(jīng)驗(yàn)公式，具體形式如下

(9)

將式(6)、(7)、(8)分別帶入到式(3)、(4)、(5)中，便可以得到未知參數(shù)A、B、C、R1、R2、ω同炸藥密度ρ0和爆速D之間的關(guān)系。

根據(jù)大量實(shí)驗(yàn)證明，R1、R2、ω并非是相互獨(dú)立的變量，R1、R2、ω的取值范圍分別為[4,7]、[0.8,2]、[0,1][13]。此外根據(jù)文獻(xiàn)[14]對(duì)于CHNO型炸藥，ω一般取0.33，且R1與R2有類線性關(guān)系，R2=0.27R1，所以僅需要在[4,7]的范圍內(nèi)選取R1的值便可以對(duì)A、B、C進(jìn)行求解。

得到一組未知參數(shù)A、B、C的解之后就需要對(duì)于解的正確性進(jìn)行判斷。選擇標(biāo)準(zhǔn)圓筒試驗(yàn)中爆轟產(chǎn)物比容為2.4和7.0的函數(shù)值是否接近于零來判斷解是否符合要求。

(10)

式中：u為圓筒特征處的速度值；ρm為圓筒材料的密度；ρm=8.93 g/cm3；ri0和re0分別為圓筒的初始內(nèi)外徑，ri0=12.7 mm,re0=15.2 mm;E0是能量計(jì)算中的一個(gè)基準(zhǔn)值，對(duì)于大多數(shù)的凝聚態(tài)炸藥根據(jù)文獻(xiàn)[15]，有

E0=(0.204-0.0734ρ0)ρ0D2

(11)

對(duì)于特征處的速度值，文獻(xiàn)[8]提出一種較為準(zhǔn)確的經(jīng)驗(yàn)計(jì)算公式

(12)

(13)

計(jì)算程序圖如圖1所示，∈為極小正實(shí)數(shù)。

圖1 計(jì)算程序圖Fig. 1 Computational Diagram

2 數(shù)值模擬

2.1 材料參數(shù)選取

采用關(guān)鍵字*MAT_HIGH_EXPLOSIVE_BURN和關(guān)鍵字*EOS-JWL來描述炸藥材料及爆轟產(chǎn)物的狀態(tài)方程。對(duì)于2#巖石乳化炸藥，產(chǎn)品說明書主要性能指標(biāo)如表1所示。采用上文中的計(jì)算方法，利用MATLAB軟件編寫計(jì)算程序，計(jì)算JWL參數(shù)，得到的材料參數(shù)如表2所示。

表1 2#巖石乳化炸藥主要性能指標(biāo)Table 1 Main performance index of 2 # rock emulsion explosive

表2 2#巖石乳化炸藥材料計(jì)算參數(shù)Table 2 Material calculation parameters of 2 # rock emulsion explosive

鉛采用Johnson-cook本構(gòu)模型及Gruneisen狀態(tài)方程描述，其具體形式如下

(14)

(γ0+αμ)E

(15)

式中：μ=ρ/ρ0-1；α為γ0的體積修正系數(shù)；C0、Si、γ0為狀態(tài)方程中的常數(shù)。具體參數(shù)值見表3所示。

表3 鉛材料參數(shù)Table 3 Lead material parameter

2.2 數(shù)值模型建立

2.2.1 鉛柱壓縮實(shí)驗(yàn)

鉛柱壓縮法，是測(cè)定炸藥猛度的一種方法。在鋼板中央放置φ40×60 mm的鉛柱，鉛柱上放置φ41×10 mm圓鋼板，鋼板上放50 g炸藥，包裝成直徑d=40 mm紙筒，放入雷管，把藥筒與鋼板固定好。引爆后，鉛柱壓縮的高度，用mm表示的值稱為該炸藥的猛度[16]。為了保證計(jì)算的準(zhǔn)確性，采用同試驗(yàn)參數(shù)相同的設(shè)置。利用ANSYS軟件前處理建立三維計(jì)算模型，采用1/4建模，模型尺寸為：鉛柱φ40×60 mm，鋼板φ41×10 mm，炸藥φ40×30 mm，裝藥量為50 g。采用多物質(zhì)流固耦合(ALE)算法進(jìn)行模擬，空氣部分φ60×100 mm，劃分網(wǎng)格采用均勻網(wǎng)格劃分，網(wǎng)格大小為1 mm，單元總數(shù)為262080。起爆位置為裝藥中心，數(shù)值模型如圖2所示。

圖2 鉛柱壓縮計(jì)算模型圖Fig. 2 Calculation model diagram of lead column compression

2.2.2 鉛鑄擴(kuò)孔實(shí)驗(yàn)

鉛鑄擴(kuò)孔實(shí)驗(yàn)是測(cè)定炸藥做功能力的一種方法。系將10 g炸藥置于圓柱形鉛鑄中央的孔中(鉛鑄直徑及高均為200 mm，孔徑25 mm，深125 mm，用細(xì)石英砂填充到孔口。引爆炸藥后，爆轟產(chǎn)物將孔擴(kuò)張為梨形，測(cè)量孔的擴(kuò)張?bào)w積，以此衡量炸藥的做功能力[17]。為了保證計(jì)算的準(zhǔn)確性，采用同試驗(yàn)參數(shù)相同的設(shè)置。

利用ANSYS軟件前處理建立三維計(jì)算模型，采用1/4建模，模型尺寸為：鉛鑄φ200×200 mm，孔徑為25 mm，孔深為125 mm，炸藥裝藥高度為20 mm，裝藥量為10 g。采用多物質(zhì)流固耦合(ALE)算法進(jìn)行模擬，空氣單元設(shè)置底部大小與模型等大，高度為模型高度的1.5倍，劃分網(wǎng)格采用均勻網(wǎng)格劃分，網(wǎng)格大小為2 mm，單元總數(shù)為136540。起爆位置為裝藥中心，數(shù)值模型如圖3所示。

圖3 鉛鑄擴(kuò)孔模型圖Fig. 3 Explosion power test model

2.3 數(shù)值模擬結(jié)果

2.3.1 鉛柱壓縮實(shí)驗(yàn)?zāi)M結(jié)果

利用后處理軟件LS-PrePost對(duì)于數(shù)值模型進(jìn)行軸對(duì)稱處理，鉛柱壓縮過程的變形云圖，如圖4所示。在炸藥起爆后，爆轟波從裝藥中心向四周擴(kuò)散。沖擊波壓力通過鋼板傳遞到鉛柱的上表面，鉛柱開始?jí)嚎s。當(dāng)壓力達(dá)到鉛柱材料的屈服強(qiáng)度后，鉛柱塑性變形急劇增長(zhǎng)，由于剛性地面的作用，壓力無法傳遞，鉛柱上表面出現(xiàn)擴(kuò)散狀變形。隨著壓縮的進(jìn)行，鉛柱上表面逐漸出現(xiàn)“蘑菇狀”卷曲。

圖4 鉛柱壓縮過程變形云圖Fig. 4 Deformation cloud image of lead column during compression

在整個(gè)鉛柱壓縮過程中，鋼板的變形量極小，因此在鋼板的不同位置取點(diǎn)，記錄鋼板z方向上的位移量作為鉛柱壓縮量，繪制鉛柱壓縮量-時(shí)間曲線如圖5所示。對(duì)應(yīng)測(cè)點(diǎn)的z方向上測(cè)點(diǎn)的速度-時(shí)間曲線如圖6所示。

圖5 壓縮量-時(shí)間曲線圖Fig. 5 Compression-time graph

圖6 測(cè)點(diǎn)速度-時(shí)間曲線圖Fig. 6 Velocity-time curve of measuring point

炸藥起爆后，鉛柱在40 μs左右達(dá)到了最大壓縮速度，之后隨著時(shí)間的增長(zhǎng)壓縮速度逐漸減小，并在800 μs時(shí)接近于0，800 μs之后鉛柱的壓縮量基本不再變化，認(rèn)為鉛柱壓縮過程完成。

2.3.2 鉛鑄擴(kuò)孔試驗(yàn)?zāi)M結(jié)果

利用后處理軟件LS-PrePost對(duì)于數(shù)值模型進(jìn)行軸對(duì)稱處理，觀察鉛鑄擴(kuò)孔過程，其變形云圖如圖7所示。隨著炸藥的起爆，爆生氣體迅速膨脹，以近似球面的形式作用于鉛鑄孔壁以及堵塞材料上，早期爆腔的形狀基本為球形。隨著爆生氣體的膨脹，爆生氣體的壓力也在不斷減小，爆腔的擴(kuò)張速度也逐漸減慢。鉛鑄底部為剛性地面，且石英砂堵塞裝填密度較低，鉛鑄內(nèi)孔逐漸被擴(kuò)張為梨形,上表面內(nèi)孔邊緣出現(xiàn)了一定尺寸的凸起。

圖7 鉛鑄擴(kuò)孔過程變形云圖Fig. 7 Explosion power test deformation nephogram

鉛鑄擴(kuò)孔數(shù)值模擬的結(jié)果，由于后處理軟件LS-PrePost無法對(duì)于擴(kuò)孔空腔進(jìn)行直接的體積測(cè)量，所以需要對(duì)數(shù)值實(shí)驗(yàn)結(jié)果進(jìn)行二次處理。利用LS-PrePost對(duì)于鉛鑄擴(kuò)孔結(jié)果模型沿著不同徑向進(jìn)行切分，輸出擴(kuò)孔空腔邊緣各節(jié)點(diǎn)坐標(biāo)，利用ANSYS前處理軟件，對(duì)于空腔進(jìn)行重建，生成K文件，利用LS-PrePost測(cè)量重建空腔模型的體積。測(cè)量不同時(shí)間的空腔體積繪制空腔體積—時(shí)間曲線如圖8所示。取空腔邊緣具有代表性的測(cè)點(diǎn)，繪制速度-時(shí)間曲線圖如圖9所示?？涨惑w積的增長(zhǎng)速度在隨著時(shí)間的推移逐漸減小。在1000 μs左右，空腔體積基本上不再變化，擴(kuò)孔過程基本完成。

圖8 空腔體積-時(shí)間曲線Fig. 8 Cavity volume-time curve

圖9 測(cè)點(diǎn)速度-時(shí)間曲線圖Fig. 9 Velocity-time curve of measuring point

3 數(shù)值模擬結(jié)果分析

對(duì)于不同方向上測(cè)得的空腔體積求平均值，得到鉛鑄擴(kuò)孔數(shù)值實(shí)驗(yàn)的結(jié)果，測(cè)得的鉛鑄擴(kuò)孔體積為335.053 mL。鉛鑄擴(kuò)孔實(shí)驗(yàn)測(cè)得的炸藥爆力值按下式計(jì)算

X=(V2-V1)(1+K)-22

(16)

式中：X為炸藥的爆力值，mL；V1為起爆前鉛鑄內(nèi)孔的體積，mL；V2為起爆后鉛鑄擴(kuò)孔空腔的體積，mL；K為溫度修正系數(shù)，實(shí)驗(yàn)室溫為20℃，K=-2%。

按照2#巖石乳化炸藥說明書爆力值280 mL，代入式(16)中計(jì)算得，起爆后鉛鑄擴(kuò)孔空腔的體積為369.51 mL，數(shù)值模擬結(jié)果的誤差為9.3%，滿足要求。

對(duì)于鉛柱壓縮數(shù)值模擬結(jié)果，鉛柱壓縮量可以直接通過測(cè)量鋼板位移量來確定。經(jīng)過測(cè)量，鉛柱壓縮量為13.08 mm。按照2#巖石乳化炸藥說明書壓縮量16 mm計(jì)算，數(shù)值模擬實(shí)驗(yàn)結(jié)果的誤差為18.25%，超過了15%，誤差較大，所以需要對(duì)參數(shù)進(jìn)行修正。

PS=Ae-R1V+Be-R2V+CV-(ω+1)

(17)

在JWL狀態(tài)方程中，第一項(xiàng)控制高壓段，第二項(xiàng)控制中壓段，R1、A為高壓段控制參數(shù)，R2、B為中壓段控制參數(shù)。以計(jì)算的得到的2#巖石乳化炸藥JWL狀態(tài)方程參數(shù)為原始參數(shù)集，控制其他參數(shù)不變，將A、B、R1、R2作為因素，設(shè)置三個(gè)水平，進(jìn)行正交試驗(yàn)，對(duì)這四個(gè)參數(shù)進(jìn)行敏感性分析并確定相應(yīng)的參數(shù)調(diào)整方向，正交試驗(yàn)表及結(jié)果如表4所示。

表4 正交試驗(yàn)表及結(jié)果Table 4 Orthogonal numerical test scheme and results

對(duì)正交試驗(yàn)結(jié)果進(jìn)行方差分析(F檢驗(yàn))，以起爆后空腔體積表征爆力，鉛柱壓縮量表征猛度。試驗(yàn)因素及誤差的自由度分別為2和6，當(dāng)顯著性水平分別為α=0.05、α=0.01時(shí)，查表可知F臨界值為5.140、10.90；計(jì)算得出試驗(yàn)因素的F比如表5所示。

表5 試驗(yàn)因素F值表Table 5 Test factor F value

當(dāng)試驗(yàn)因素值小于10.90，大于5.140時(shí)，認(rèn)為該因素對(duì)于結(jié)果影響比較顯著；當(dāng)試驗(yàn)因素F值大于10.90時(shí)，認(rèn)為該因素對(duì)于結(jié)果影響十分顯著；并且F值越大表明對(duì)于結(jié)果的影響也就越顯著。通過表4、5分析JWL狀態(tài)方程中的四個(gè)參數(shù)變化對(duì)于結(jié)果的影響：

(1)爆力對(duì)于參數(shù)R1的變化較為敏感，并且擴(kuò)孔空腔體積隨著參數(shù)R1的減小而增大；對(duì)于剩余的三個(gè)參數(shù)敏感性A?R2?B，鉛鑄擴(kuò)孔空腔體積分別隨著參數(shù)A、B、R2的增大而增大。

(2)猛度對(duì)于參數(shù)R1的變化極為敏感，并且壓縮量隨著參數(shù)R1的減小而增大；對(duì)于剩余的三個(gè)參數(shù)敏感性A?R2?B，壓縮量分別隨著參數(shù)A、B、R2的增大而增大。

(3)猛度對(duì)于參數(shù)R1的敏感性明顯強(qiáng)于爆力。

猛度表征炸藥對(duì)于接觸介質(zhì)的破壞能力，反映的是炸藥瞬間爆轟釋放的能量，主要受爆速及爆轟壓力影響。但從正交試驗(yàn)結(jié)果來看，猛度也與爆生氣體作用有關(guān)，體現(xiàn)在猛度對(duì)與JWL狀態(tài)方程高壓段參數(shù)R1、A的敏感性上。相較于猛度，爆力對(duì)JWL狀態(tài)方程高壓段參數(shù)R1、A的敏感性較弱，因?yàn)楸Ψ从车氖潜鷼怏w膨脹做功的能力，受爆生氣體整個(gè)膨脹過程的影響。爆生氣體作用時(shí)間較長(zhǎng)，受JWL狀態(tài)方程參數(shù)A、B、R1、R2的協(xié)同作用影響。

對(duì)于炸藥材料參數(shù)的修正，主要修正JWL狀態(tài)方程中的四個(gè)參數(shù)A、B、R1、R2，其余參數(shù)的計(jì)算公式經(jīng)過前人大量的驗(yàn)證大多誤差不超過3%。對(duì)于A、B、R1、R2的修正，結(jié)合爆力、猛度對(duì)于參數(shù)的變化關(guān)系，先匹配猛度值，率先確定R1、A，再匹配爆力，確定R2、B的值。必要時(shí)可采用插值的方法，將爆力、猛度結(jié)果誤差控制在5%以內(nèi)。

修正后的2#巖石乳化炸藥材料參數(shù)如表6所示。

表6 修正后2#巖石炸藥材料參數(shù)Table 6 Modified material parameters of 2 # rock explosive

4 結(jié)論

(1)以工業(yè)炸藥產(chǎn)品說明書上的爆速、密度作為已知參數(shù)，以2#巖石乳化炸藥為例，采用一種簡(jiǎn)便算法直接計(jì)算得出了JWL狀態(tài)方程未知參數(shù)。

(2)利用計(jì)算得到的2#巖石乳化炸藥的JWL狀態(tài)方程未知參數(shù)，進(jìn)行鉛柱壓縮數(shù)值模擬和鉛鑄擴(kuò)孔數(shù)值模擬，兩種實(shí)驗(yàn)基本模擬出了現(xiàn)場(chǎng)試驗(yàn)的結(jié)果形狀，對(duì)計(jì)算參數(shù)進(jìn)行修正后，數(shù)值模擬結(jié)果的相對(duì)誤差大大減小。

(3)分析了爆力及猛度對(duì)JWL狀態(tài)方程中的四個(gè)主要參數(shù)A、B、R1、R2的敏感性，R1>A>R2>B，確定了工業(yè)炸藥爆力及猛度關(guān)于JWL狀態(tài)方程未知參數(shù)的變化關(guān)系及計(jì)算參數(shù)修正方法。為復(fù)雜組分炸藥、含添加劑炸藥及新型炸藥的JWL狀態(tài)方程參數(shù)的擬合提供了理論依據(jù)及參考指導(dǎo)。