爆炸沖擊波數(shù)值模擬及超壓計算公式的修正

吳彥捷，高軒能

（華僑大學(xué) 土木工程學(xué)院，福建 廈門361021）

炸藥爆炸的瞬時能釋放巨大的能量并產(chǎn)生各種效應(yīng)，但對建筑物破壞力最強(qiáng)、影響區(qū)域最大的是爆炸沖擊波，而沖擊波超壓是評價炸藥爆炸效應(yīng)的重要參數(shù)之一.關(guān)于爆炸空氣沖擊波超壓的計算，早期常用的經(jīng)驗(yàn)公式有 GB 6722-2003《爆破安全規(guī)程》［1］、Henrych公式［2］、Sadovskyi公式［2-3］、Brode公式［4］、Aliansov公式［5］、Baker公式［6］等.20世紀(jì)末，不少學(xué)者也提出相應(yīng)的計算公式［7-11］，但這些經(jīng)驗(yàn)公式主要基于試驗(yàn)數(shù)據(jù)和理論分析結(jié)果得到，因爆炸作用的時程極短，爆炸沖擊波迅速衰減，使實(shí)驗(yàn)結(jié)果的準(zhǔn)確性受到影響.近年來，隨著計算機(jī)技術(shù)的快速發(fā)展，數(shù)值模擬方法也愈來愈多地成為爆炸效應(yīng)研究的重要手段，眾多學(xué)者通過爆炸試驗(yàn)與數(shù)值模擬相結(jié)合的方法對沖擊波超壓進(jìn)行一系列的研究，并提出相應(yīng)的計算公式［12-21］.由前人的研究可知，爆炸數(shù)值模擬結(jié)果與理論值之間存在一定誤差，且基于有限試驗(yàn)得到的若干爆炸沖擊波經(jīng)驗(yàn)公式適用范圍和標(biāo)準(zhǔn)不一致，有必要對TNT炸藥爆炸沖擊波超壓計算公式進(jìn)行進(jìn)一步修正.本文通過數(shù)值模擬結(jié)果與經(jīng)驗(yàn)公式的對比分析，驗(yàn)證了計算模型及參數(shù)取值的可信性，提出爆炸沖擊波修正計算方法.

1 爆炸沖擊波超壓ΔPf的計算

分別選取常用比例距離區(qū)段ˉR＝0.3～1.0 m·kg-1／3及ˉR＝1.0～15.0 m·kg-1／3，收集文獻(xiàn)［1-21］中提到的沖擊超壓經(jīng)驗(yàn)公式，對其進(jìn)行對比分析，研究不同經(jīng)驗(yàn)公式之間的誤差.其中，文獻(xiàn)［6-11］中的公式適用于高爆炸藥的沖擊波超壓計算，文獻(xiàn)［16-21］中的公式則是針對國防及核爆炸提出的，其余文獻(xiàn)公式均適用于計算TNT炸藥在無限空氣中爆炸的沖擊波超壓.為便于參考比較，將各經(jīng)驗(yàn)公式計算結(jié)果進(jìn)行加權(quán)平均，設(shè)為平均值，并將平均值作為沖擊波超壓的中位值進(jìn)行考慮，對比結(jié)果如圖1所示.

沖擊波超壓平均值函數(shù)化表達(dá)式為

圖1 經(jīng)驗(yàn)公式的 ΔPf-ˉR 曲線匯總Fig.1 ΔPf-ˉR curves of all empirical formulae

由圖1（a）可知：當(dāng)比例距離為0.3～1.0 m·kg-1／3時，各經(jīng)驗(yàn)公式總體變化趨勢基本保持一致，隨著比例距離ˉR的增大，沖擊波超壓由于能量耗散而逐漸衰減.文獻(xiàn)［6-11］提到的高爆炸藥沖擊波計算結(jié)果普遍大于平均值，而文獻(xiàn)［16-21］針對國防及核爆炸公式計算結(jié)果則普遍低于平均值，這說明炸藥性能及爆炸類型使沖擊波超壓計算結(jié)果之間存在偏差.當(dāng)比例距離ˉR＝0.3 m·kg-1／3時，各經(jīng)驗(yàn)公式的最大值和最小值分別為90.5，8.69 MPa，之間相差近10倍，且與平均值35.3 MPa的相對誤差達(dá)到-61.0%和75.4%，存在較大波動.在選取不同經(jīng)驗(yàn)公式計算時，可能造成結(jié)果偏差過大，但Sadovskyi公式［2］計算結(jié)果與平均值較為吻合.因此，在比例距離ˉR小于1.0 m·kg-1／3時，Sadovskyi公式［2］可以用于估算小范圍內(nèi)的爆炸沖擊波超壓.

由圖1（b）可知：當(dāng)比例距離ˉR為1.0～15.0 m·kg-1／3時，早期常用經(jīng)驗(yàn)公式的沖擊波超壓普遍低于平均值，處于經(jīng)驗(yàn)公式的中位值或下位值.如果采用常用的經(jīng)驗(yàn)公式對沖擊波超壓進(jìn)行估算時，結(jié)果容易產(chǎn)生偏差，低估了爆炸沖擊波的威力，導(dǎo)致計算結(jié)果安全性降低.當(dāng)比例距離ˉR＝1.0 m·kg-1／3時，各經(jīng)驗(yàn)公式的平均值為1.29 MPa，而最大值和最小值分別為2.81，0.51 MPa，之間相對誤差減小到-54.1%和60.5%.當(dāng)比例距離ˉR從0.3 m·kg-1／3變化到1.0 m·kg-1／3，相對誤差呈現(xiàn)明顯減小的趨勢.隨著比例距離ˉR的增大，各經(jīng)驗(yàn)公式計算結(jié)果逐漸趨于相近，相互之間的誤差也逐漸減小.

通過以上分析可知：各沖擊波超壓經(jīng)驗(yàn)公式之間均存在一定偏差，且適用范圍不一致.這是因?yàn)樵缙诘慕?jīng)驗(yàn)公式多數(shù)是通過爆炸試驗(yàn)測得，而在爆炸試驗(yàn)過程中，由于試驗(yàn)條件、測試范圍的限制，以及地面或其他剛性物體產(chǎn)生的反射波效應(yīng)，可能使沖擊波超壓得到增強(qiáng)，也使試驗(yàn)實(shí)測結(jié)果離散性較大.

2 TNT炸藥爆炸數(shù)值模擬

圖2 空氣與炸藥的有限元模型Fig.2 Finite element model of explosive and air

2.1 材料參數(shù)取值

應(yīng)用LS-DYNA有限元軟件建立自由空爆模型.空氣尺寸為20 m×20 m×20 m，炸藥為立方體0.2 m×0.2 m×0.2 m，網(wǎng)格尺寸按0.2 m×0.2 m×0.2 m均勻劃分.以炸藥為中心取1／8模型進(jìn)行計算，單元類型均為8節(jié)點(diǎn)3D-SOLID164，采用ALE（Arbitrary Lagrange-Euler）算法.在XOY，XOZ，YOZ平面采用對稱約束，其他面采用透射邊界以模擬無限空氣域，空氣與炸藥有限元模型如圖2所示.炸藥和空氣按均勻連續(xù)介質(zhì)考慮，炸藥采用MAT＿HIGH＿EXPLOSIVE＿BURN材料模型和JWL（Jones-Wilkins-Lee）狀態(tài)控制方程，爆炸沖擊壓力［22］為

式（2）中：A，B，R1，R2，ω為輸入?yún)?shù)；V為相對體積；E0為初始內(nèi)能.TNT炸藥的材料參數(shù)取值如下：炸藥密度ρ為1 630 kg·m-3；炸藥爆速D為6 713 m·s-1；爆壓PCJ為18.5 GPa；A，B分別為540.9，9.4 GPa；R1，R2，ω分別為4.5，1.1，0.35；E0為8 GJ·m-3；初始相對體積V0為1.0.

空氣采用MAT＿NULL空材料模型和線性多項(xiàng)式方程EOS＿LINEAR＿POLYNOMIAL，即

式（3）中：μ＝ρ／ρ0-1，ρ為空氣質(zhì)量密度，ρ0為參考質(zhì)量密度；E為單位體積內(nèi)能.線性多項(xiàng)式狀態(tài)方程遵守 Gamma定律，c0，c1，c2，c3，c4，c5，c6為實(shí)常數(shù).空氣的材料輸入?yún)?shù)如下：ρ為1.290 kg·m-3；c0～c3都為0；c4，c5都為0.4；c6為0；E0為0.25 MJ·m-3；V0為1.0.

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

由爆炸數(shù)值計算模型生成相應(yīng)的K文件，將全部信息輸入K文件，通過LS-DYNA求解器進(jìn)行計算，運(yùn)用常用后處理軟件LS-PREPOST4.0對以上求解結(jié)果進(jìn)行后處理；提取比例距離ˉR＝0.6～1.0 m·kg-1／3及ˉR＝1.0～6.0 m·kg-1／3范圍內(nèi)有限測點(diǎn)的沖擊波超壓.數(shù)值模擬結(jié)果，如表1所示.

表1 沖擊波超壓數(shù)值模擬結(jié)果Tab.1 Numerical simulation results of shock wave overpressure

結(jié)合表1數(shù)值模擬結(jié)果，擬合得到?jīng)_擊波超壓計算公式，具體表達(dá)式為

圖3為沖擊波超壓數(shù)值計算結(jié)果與經(jīng)驗(yàn)公式的比較，并給出了擬合曲線式（4）.由圖3可知：數(shù)值模擬結(jié)果總體上與經(jīng)驗(yàn)公式符合較好，且與Henrych公式最為接近，但均小于各經(jīng)驗(yàn)公式計算結(jié)果.這是因?yàn)閿?shù)值模擬計算基于理論狀態(tài)方程建模，是理想狀態(tài)的爆炸效應(yīng)，得出的沖擊波超壓較試驗(yàn)數(shù)據(jù)偏小也是合理的.因此，數(shù)值模擬結(jié)果可以作為爆炸沖擊波超壓的下限值.隨著比例距離ˉR的增大，數(shù)值模擬結(jié)果與經(jīng)驗(yàn)公式的誤差逐漸減小.總體而言，數(shù)值模擬結(jié)果與經(jīng)驗(yàn)公式的總體變化趨勢一致且吻合較好，驗(yàn)證了計算模型與參數(shù)取值的可信性.

3 沖擊波超壓計算公式修正

已有的沖擊波超壓經(jīng)驗(yàn)公式適用范圍和標(biāo)準(zhǔn)不一致，且相互之間誤差較大，通過數(shù)值模擬及經(jīng)驗(yàn)公式計算結(jié)果進(jìn)行對比分析，提出數(shù)值模擬的系數(shù)修正方法，并擬合出統(tǒng)一的爆炸沖擊波超壓計算公式.

3.1 沖擊波超壓修正計算方法

由圖3可知：常用經(jīng)驗(yàn)公式計算結(jié)果普遍處于沖擊波超壓的中位值或下位值，而數(shù)值模擬結(jié)果則處于下限值.與經(jīng)驗(yàn)公式相比，數(shù)值模擬結(jié)果偏低，若抗爆設(shè)計過程中以數(shù)值模擬結(jié)果為依據(jù)進(jìn)行設(shè)計時，有存在低估爆炸沖擊波超壓的危險.為了確?？贡O(shè)計的安全性，在數(shù)值模擬結(jié)果ΔPf的基礎(chǔ)上乘以一個修正系數(shù)ηp，得到修正后的沖擊波ΔP′f，具體表達(dá)式為

隨著比例距離ˉR的增大，數(shù)值模擬結(jié)果更為精確，即修正系數(shù)ηp與比例距離ˉR成反比關(guān)系，且趨勢逐漸平緩.因此，選用冪函數(shù)η＝aˉRb＋c形式進(jìn)行擬合能夠準(zhǔn)確地描述修正函數(shù).

選取常用比例距離區(qū)段ˉR為0.6～1.0及1.0～6.0 m·kg-1／3，對圖1中的各經(jīng)驗(yàn)公式進(jìn)行加權(quán)平均，得到平均值ΔP′f，e，如式（1）所示.將ΔP′f，e作為參考，通過ΔP′f及ΔP′f，e求得ηp，利用修正系數(shù)對數(shù)值模擬計算結(jié)果進(jìn)行修正，以求得到?jīng)_擊波超壓的中位值，使得修正后的ΔP′f接近ΔP′f，e.擬合得到修正函數(shù)如圖4所示.

修正函數(shù)ηp具體表達(dá)式為

3.2 誤差分析

作為參考，通過驗(yàn)證該修正計算方法的適用性，對修正函數(shù)進(jìn)行誤差分析，結(jié)果如表2所示.表2中：δ為ΔP′f和ΔP′f，e的相對誤差百分比.從表2可知：相對誤差δ基本控制在±1%范圍內(nèi)，只有少數(shù)超過±1%，擬合程度較高.說明通過上述修正計算方法能夠有效地修正數(shù)值模擬計算結(jié)果，并得到?jīng)_擊波超壓的中位值，使得數(shù)值模擬結(jié)果能夠有效地運(yùn)用于抗爆設(shè)計中.

3.3 修正后的超壓計算公式

通過以上提出的超壓修正計算方法，得到統(tǒng)一的沖擊波超壓修正計算公式，具體表達(dá)式為

表2 修正函數(shù)ηp誤差分析Tab.2 Error analysis of modified functionηp

4 結(jié)論

通過上述對TNT炸藥爆炸數(shù)值模擬及沖擊波超壓計算公式的修正，可以得出以下3點(diǎn)結(jié)論.

1）早期常用的經(jīng)驗(yàn)公式普遍處于經(jīng)驗(yàn)公式的中位值或下位值，存在了低估爆炸沖擊波的危險.在比例距離ˉR小于1.0 m·kg-1／3時，Sadovskyi公式［2］可以用于估算小范圍內(nèi)的爆炸沖擊波超壓.

2）基于LS-DYNA有限元程序?qū)崿F(xiàn)TNT炸藥爆炸的數(shù)值模擬計算是可行的，且數(shù)值模擬計算結(jié)果可以作為爆炸沖擊波超壓的下限值.

3）通過沖擊波超壓公式修正計算方法對數(shù)值模擬結(jié)果進(jìn)行修正，能夠有效地提高數(shù)值模擬結(jié)果的可用性，提出統(tǒng)一的計算公式并得到?jīng)_擊波超壓中位值，為抗爆設(shè)計提供參考.

［1］ 國家質(zhì)量監(jiān)督檢驗(yàn)檢疫總局.GB 6722-2003爆破安全規(guī)程［S］.北京：中國標(biāo)準(zhǔn)出版社，2003：47.

［2］ 亨利奇J.爆炸動力學(xué)及其應(yīng)用［M］.熊建國，譯.北京：科學(xué)出版社，1987：127.

［3］ SADOVSKYI M A.Mechanical action of air shock waves of explosion，based on experimental data［M］.Moscow：Izd Akad Nauk SSSR，1952：1-2.

［4］ BRODE H L.Blast wave from a spherical charge［J］.Phys Fluids，1959，2（2）：217.

［5］ 李翼祺，馬素貞.爆炸力學(xué)［M］.北京：科學(xué)出版社，1959：260.

［6］ BAKER W E.Explosions in air［M］.Austin：University of Texas Press，1974：6-10.

［7］ FAHMY A，TOLBA F.Response of FRP-retrofitted reinforced concrete panels to blast loading［M］.Canada：Carleton University，2001：12.

［8］ 葉序雙，戴鎮(zhèn)華，莊兆鈴.非理想剛性地面球形裝藥爆炸沖擊波超壓的試驗(yàn)研究［J］.解放軍理工大學(xué)學(xué)報：自然科學(xué)版，1988（4）：72-82.

［9］ 惠君明，陳天云.炸藥爆炸理論［M］.南京：江蘇科學(xué)技術(shù)出版社，1995：145-146.

［10］ 葉曉華.軍事爆破工程［M］.北京：解放軍出版社，1999：216-219.

［11］ WU Cheng-qing，HAO Hong.Modeling of simultaneous ground shock and air blast pressure on nearby structures from surface explosions［J］.Int J Impact Eng，2005，31（6）：699-717.

［12］ 楊鑫，石少卿，程鵬飛.空氣中 TNT爆炸沖擊波超壓峰值的預(yù)測及數(shù)值模擬［J］.爆破，2008，25（1）：15-19.

［13］ 仲倩，王伯良，黃菊，等.TNT空中爆炸超壓的相似律［J］.火炸藥學(xué)報，2008，33（4）：32-35.

［14］ 董桂旭，杜茂華，黃雪峰.某型炸藥的沖擊波超壓峰值計算公式參數(shù)的修正［J］.海軍航空工程學(xué)院學(xué)報，2010，25（5）：542-544.

［15］ 胡華權(quán)，裴明敬，許學(xué)忠，等.燃料空氣炸藥爆炸威力評價方法研究［C］∥第四屆全國爆炸力學(xué)實(shí)驗(yàn)技術(shù)學(xué)術(shù)會議.合肥：［出版者不詳］，2006：313-318.

［16］ 王儒策，趙國志.彈丸終點(diǎn)效應(yīng)［M］.北京：北京理工大學(xué)出版社，1993：45.

［17］ 喬登江.核爆炸物理概論［M］.北京：原子能出版社，1988：30-70.

［18］ 陶俊林，黎澤朋，陳小偉.核爆聚變電站發(fā)電時的沖擊波超壓傳播規(guī)律［J］.高壓物理學(xué)報，2012，26（2）：221-226.

［19］ 國家技術(shù)監(jiān)督局，中華人民共和國建設(shè)部.GB 50038-2005人民防空地下室設(shè)計規(guī)范［S］.北京：中國建筑工業(yè)出版社，2006：161.

［20］ 王年橋.防護(hù)結(jié)構(gòu)計算原理與設(shè)計［M］.南京：工程兵工程學(xué)院出版社，1998：176-178.

［21］ 樊海青.爆炸沖擊波作用下多層多跨框架結(jié)構(gòu)的塑性極限分析［D］.北京：北方工業(yè)大學(xué)，2009：21.

［22］ 段雷琳，高軒能，江媛.內(nèi)爆炸下球面鋼網(wǎng)殼結(jié)構(gòu)的沖擊波超壓計算分析［J］.華僑大學(xué)學(xué)報：自然科學(xué)版，2013，34（5）：557-562.