高軒能,吳彥捷
(華僑大學土木工程學院,福建 廈門361021)
炸藥爆炸瞬時能釋放巨大的能量并產生各種效應,但破壞力最強、影響區域最大的是爆炸沖擊波。早期研究爆炸沖擊波效應主要以實驗為主,但因爆炸作用的時程極短,通常在數十毫秒內,爆炸沖擊波即從最大值變為零,影響試驗結果的準確性。近年來,隨著計算機技術的快速發展,數值模擬方法已成為研究爆炸效應的重要手段。
爆炸空氣沖擊波超壓計算常用方法有Sadovskyi公式、Henrych 公式、Brode公式,Aliansov公式 和TM5-1300表格等[1-4]。葉序雙等[5]通過測量沖擊波傳播速度轉換計算不同測點處的沖擊波壓力,對非理想剛性地面球形炸藥爆炸沖擊波超壓進行實驗研究并得到了大藥量TNT 炸藥在地面爆炸時的沖擊波超壓計算公式;仲倩等[6]通過爆炸試驗測得空氣沖擊波峰值超壓,對經驗計算公式進行了系數修正;劉偉等[7-10]進行了TNT爆炸試驗,并與有限元數值計算結果進行了對比研究,兩者符合較好;楊鑫等[11]將數值計算結果與Henrych等沖擊波超壓經驗公式進行了比較,指出數值計算結果普遍小于經驗公式;李秀地等[12]在坑道爆炸試驗數據的基礎上,運用數值模擬方法研究了T 型坑道爆炸沖擊波的傳播衰減規律;盧紅琴等[13-14]討論了有限元網格密度和空氣方程等參數變化對數值計算結果精度的影響。由此可見,現有研究主要集中在試驗研究或數值計算,對引起爆炸數值模擬與經驗計算公式結果之間誤差程度不一的原因探討較少。
本研究采用LS-DYNA 有限元程序建立TNT爆炸的數值計算模型,研究了空氣沖擊波的傳播特性,結合經驗公式和已有試驗數據,驗證計算模型及參數取值的可信性,分析探討不同參數取值對沖擊波超壓的影響。
應用LS-DYNA 有限元程序建立自由空爆模型。空氣尺寸取12m×12m×12m 的較小空域,以節省運算時間,炸藥尺寸取0.2m×0.2m×0.2m 的立方體,網格尺寸按0.1m×0.1m×0.1m 劃分。空爆以炸藥為中心取1/8模型計算,單元類型取8節點3D-SOLID164,采用ALE(Arbitrary Lagrange-Euler)算法。在XOY、XOZ、YOZ 平面內采用對稱約束,其他面采用透射邊界以模擬無限空氣域。為考慮剛性地面對沖擊波超壓的影響,地面單元類型取為SHELL163,采用MAT_RIGID剛體材料模型。
剛性地面具體參數為:*MAT_RIGID;3 7 830 2.07e30 0.300 0.00 0.00 0.00 0.00;1 7 7;0.00 0.00 0.00 0.00 0.00 0.00。
炸藥和空氣按均勻連續介質考慮。炸藥采用MAT_HIGH_EXPLOSIVE_BURN 材料模型和JWL(Jones-Wilkins-Lee)狀態控制方程,爆炸沖擊波壓力為:

式中:A、B、R1、R2、ω 為輸入參數;V 為相對體積;E0為初始內能。
TNT 的材料參數見表1。

表1 炸藥的材料參數Tab 1 Material parameters of explosive
空氣采用MAT_NULL 空材料模型和線性多項式方程EOS_LINEAR_POLYNOMIAL,即:

式中:μ=ρ/ρ0-1;E 為單位體積內能;ρ 為空氣密度;ρ0為參考密度。
線性多項式狀態方程遵守Gamma定律,空氣的材料輸入參數見表2。

表2 空氣的材料參數Table 2 Material parameters of air
炸藥在無限空氣中發生爆炸時,爆炸波以炸藥為圓心向四周傳播,沖擊波隨著距離的增大能量逐漸耗散,超壓峰值迅速衰減。自由空氣爆炸沖擊波在不同時刻的傳播過程如圖1所示。
1.3.1 數值計算結果的公式表達
基于上述建立的空中爆炸數學模型,為了研究更大范圍內沖擊波超壓與比例距離的關系,模型中空氣尺寸擴展為20m×20m×20m,網格尺寸取為0.2m×0.2m×0.2m,炸藥仍為0.2m×0.2m×0.2m的立方體,同樣取1/8模型計算。經過一系列數值計算,得到比例距離0.6~6.0m/kg1/3內的沖擊波超壓(見圖2),以及1.0~3.0m/kg1/3內的正壓作用時間(見圖3)。經過擬合后,得到函數化表達的沖擊波超壓和正壓作用時間計算公式分別如下:


圖1 空爆沖擊波在不同時刻的傳播圖Fig.1 The spread process of the air explosion shock waves at different times

式中:Δpf為沖擊波超壓,MPa;,R 為 計算點到爆心的距離,m;m 為炸藥藥量,kg。

式中:t+為正壓作用時間,s。
1.3.2 數值計算結果與試驗和經驗公式計算值的比較
常用的空中爆炸沖擊波超壓(105Pa)經驗公式見式(5)~(10)。
我國《爆破安全規程》GB6722-2003公式[15]:

Sadovskyi公式[1]:


Aliansov公式[4]:

常用的正壓作用時間(s)經驗公式有:
Henrych公式[2]:

Sadovskyi公式[1]:

式(10)中:B 取1.35[4]。
將沖擊波超壓數值計算結果與經驗公式進行比較,結果見圖2,圖2中同時給出了擬合曲線(擬合曲線方程見式3)。

圖2 由數值計算和經驗公式得到的Δpf-R 曲線Fig.2 Δpf-curves obtained by the numerical calculation and empirical formulae
從圖2可以看出,數值計算結果與經驗公式計算值總體符合較好,而與Henrych 公式最為接近。隨著比例距離的增大,數值計算結果與經驗公式計算值的誤差逐漸減小。
將正壓作用時間數值計算結果與經驗公式進行比較,結果見圖3,圖3 中同時給出了擬合曲線(擬合曲線方程見式4)。

圖3 正壓作用時間數值計算與經驗公式計算結果的比較Fig.3 Comparison of the results obtained by the positive pressure time numerical simulation with the empirical formula
為進一步驗證數值計算結果的準確性,將相關文獻中TNT 炸藥爆炸試驗的沖擊波超壓實測數據[16-20]、經驗公式和數值計算結果進行比較,結果如圖4所示。因爆心附近的參數測量較困難,爆炸試驗實測的數據均在比例距離=1.0~15.0m/kg1/3,因此,圖4未給出小于1.0的值。

圖4 TNT 爆炸試驗數據、經驗公式及數值計算結果的比較Fig.4 Comparison of the TNT explosion experimental data and the results obtained by empirical formula and numerical calculation
從圖4可以看出,數值計算結果、試驗結果和經驗公式計算結果的變化趨勢一致,隨著比例距離的增大,三者的結果趨于相近,相互之間的誤差越來越小。另一方面,不同試驗結果或不同經驗公式計算值之間也存在較大誤差。例如,當為1.5m/kg1/3時,各經驗公式的Δpf最大、最小和平均值分別為0.9138,0.1768 和0.4625MPa,最大值是最小值的5.2倍,試驗所得Δpf的最大,最小和中間值分別為1.1107、0.0831和0.3279MPa,最大值是最小值的13.4倍,數值計算結果為0.2270MPa;由此可見,爆炸試驗結果給出了沖擊波超壓的上位值(圖中曲線幅寬的較大值)和中位值(平均值),Sadovskyi等經驗公式得到中位值和下位值(圖中曲線幅寬的較小值),數值計算則給出了下限值(最低值)。這是因為,在爆炸試驗過程中,由于試驗條件、測試范圍的限制,以及地面或其他剛性物體產生的反射波效應,可能使沖擊波超壓得到增強,也使實測結果離散性較大。數值計算是基于理論狀態方程建模,得到的是理想狀態的爆炸效應,沖擊波超壓偏小是合理的。
通過TNT 爆炸數值計算,以及與經驗公式計算值和試驗數據的對比分析,驗證了計算模型和參數取值的可信性。與試驗結果相比,數值計算結果可以作為爆炸沖擊波超壓的下限值,而Henrych公式、Sadovskyi公式和我國《爆破安全規程》GB6722-2003公式給出的是中位值和下位值,存在低估爆炸沖擊波超壓的危險。
本研究從文獻中按TNT 的不同材料參數取值,給出具有代表性的3組,如表3所示,分別進行數值計算。炸藥尺寸仍為0.2m×0.2m×0.2m,空氣尺寸為8m×8m×8m,取1/8模型計算。對表3材料代表值進行計算,結果見表4。
從表4可以看出,第3 組與第1 組的計算結果基本相同,在比例距離=0.45~1.65m/kg1/3范圍內,兩者誤差的絕對值不超過4%;第2組計算結果比第1 組和第3 組均小,比例距離小 于1.05m/kg1/3時,相 對 誤 差 的 絕 對 值 超 過10%,但隨著比例距離的增大,差值逐漸減小,比例 距 離大 于1.05m/kg1/3后,第2 組 與 第1 組的超壓相對誤差絕對值不超過8%。表明隨著比例距離的增大,3組材料參數取值的數值計算結果逐漸趨同。

表3 TNT 的不同材料參數Table 3 Different material parameters of TNT

表4 不同材料參數下的沖擊波超壓Table 4 The shock waves overpressure under different material parameters
TNT爆炸的數值計算模型:空氣尺寸為10m×10m×10m,炸藥尺寸分別為0.2m×0.2m×0.2m、0.3m×0.3m×0.3m、0.4m×0.4m×0.4m 和0.5m×0.5m×0.5m,即TNT 藥量分別為13.04、44.01、104.32和203.75kg組,取1/8模型計算,計算結果見表5。從表5可以看出,相同比例距離下,當TNT 藥量從13.04kg增加到203.75kg時,雖然藥量增加了14.625倍,但相應的沖擊波超壓增減沒有超過20%,呈比較平穩增加的狀態。表明在同等條件下,TNT 藥量的增減不會顯著影響TNT 爆炸的數值計算結果,僅會引起沖擊波超壓的小幅度增減。

表5 不同TNT 藥量下沖擊波超壓的數值計算結果Table 5 The numerical calculation results of shock wave overpressure under different TNT dosage
為探討網格劃分密度對數值計算結果的影響,分別按0.1、0.2、0.3和0.4m 的網格尺寸劃分單元。空氣尺寸為12m×12m×12m,炸藥取0.1m×0.1m×0.1m 的立方體,取1/8 模型計算,計算結果如圖5所示。
從圖5可看出,網格劃分密度對數值計算結果的影響程度取決于比例距離。當比例距離小于2.0m/kg1/3時,網格尺寸對計算精度有較大影響,例如,當比例距 離 為1.44m/kg1/3時,Henrych 公 式 的Δpf為0.344MPa,按0.1、0.2、0.3和0.4m 網格密度計算得到的Δpf分別為0.305、0.259、0.220和0.208MPa,與Henrych公式計算結果的誤差分別為-11.4%、-24.8%、-36.1%和-39.6%,誤差增大了28.2%,表明加密單元網格可以有效提高爆炸模擬計算的精度。隨著比例距離的增大,網格劃分密度對計算結果的影響逐漸減小。當比例距離為5.0m/kg1/3時,按不同網格劃分密度計算的結果,相對誤差小于2.5%,網格劃分密度對超壓的影響可以忽略不計。

圖5 不同單元網格劃分下的Δpf- 曲線Fig.5 Δpf-curves under different mesh sizes
考慮兩種模型:空氣尺寸8m×8m×8m、立方體炸藥0.2m×0.2m×0.2m 的1/8模型和空氣尺寸16m×16m×16m、立方體炸藥0.4m×0.4m×0.4m 的整體模型。為考察建模方式對爆炸沖擊波空間分布的影響,分別提取Z=0平面(水平面)和對角線平面(斜平面)上不同比例距離的沖擊波超壓數值計算結果進行對比,見表6。
由表6可知,相同比例距離下,整體模型數值計算的沖擊波超壓一般是1/8模型的1.03~1.25倍,多數計算點誤差在10%內,僅少數計算點誤差超過20%,但與比例距離沒有呈現明顯的規律性。主要原因在于爆炸荷載是動態荷載,取1/8模型進行計算時,有可能無法計算部分反射波的增強效應。

表6 不同建模方式下沖擊波超壓的數值計算結果Table 6 The numerical calculation results shock wave overpressure under different modeling ways
分別建立立方體(邊長12m)、圓柱體(半徑12m、徑高比1∶1)和球體(半徑12m)3 種不同空氣域的1/8模型,TNT 均為0.1m×0.1m×0.1m的立方體,數值計算結果如圖6所示。從圖6可以看出,3種不同空氣域形狀對爆炸沖擊波超壓的影響很小,立方體空氣域計算結果稍大些,也更符合經驗公式。

圖6 不同空氣域形狀下的Δpf- 曲線Fig.6 Δpf-curves under different air domain shapes
分別建立立方體、圓柱體(徑高比1∶1)和球形體TNT 在立方體空氣域內的1/8爆炸模型,TNT藥量分別為104.32、109.00和106.63kg,以便于網格劃分,藥量之差小于5%,其對超壓的影響可以忽略不計,數值計算結果如圖7所示。

圖7 不同TNT 形狀的Δpf-曲線Fig.7 Δpf-curves with different shapes of TNT
從圖7可以看出,與網格劃分密度的影響類似,炸藥形狀對數值計算結果的影響程度取決于比例距離 的 大 小。比 例 距 離小 于1.5m/kg1/3時,TNT 形狀對沖擊波超壓的影響較大,此時,在相同的比例距離下,圓柱體TNT 得到的沖擊波超壓最大,球形體次之,立方體最小。比例距離大于1.5m/kg1/3后,TNT 形狀對計算結果的影響很小,可以忽略不計。
(1)基于LS-DYNA 有限元程序實現TNT 爆炸沖擊波超壓及正壓作用時間的數值計算是可行的。數值計算結果可以作為爆炸沖擊波超壓的下限值。
(2)與試驗數據相比,用Henrych 公式、Sadovskyi公式和我國《爆破安全規程》GB6722-2003公式計算沖擊波超壓給出的是中位值和下位值,存在低估沖擊波超壓的危險。
(3)建模方式和空氣域形狀對TNT 爆炸數值計算結果的影響可以忽略不計。在相同比例距離下,數值計算的沖擊波超壓隨TNT 藥量的增加有小幅度增加。
(4)材料參數取值、單元網格密度和TNT 形狀對數值計算結果的影響與比例距離密切相關。在比例距離較小的范圍內(如小于2.0m/kg1/3或更小)進行數值計算時,不能忽視網格劃分密度、TNT形狀和材料參數取值對計算結果的影響。
[1] Sadovskyi M A.Mechanical Action of Air Shock Waves of Explosion,Based on Experimental Data[M].Moscow:Izd Akad Nauk SSSR,1952.
[2] 亨利奇.J.爆炸動力學及其應用[M].熊建國,譯.北京:科學出版社,1987.
[3] Brode H L.Blast wave from a spherical charge[J].The Physics of Fluids,1959,2(2):217-229.
[4] 李翼祺,馬素貞.爆炸力學[M].北京:科學出版社,1992.
[5] 葉序雙,戴鎮華,莊兆鈴.非理想剛性地面球形裝藥爆炸沖擊波超壓的試驗研究[J].解放軍理工大學學報(自然科學版),1988(4):72-82.YE Xu-shuang,DAI Zhen-hua,ZHUANG Zhao-ling.Study on the ideal rigid ground ball shape explosion shock wave overpressure experiments[J].Journal of PLA University of Science and Technology(Natural Science Edition),1988(4):72-82.
[6] 仲倩,王伯良,黃菊,等.TNT 空中爆炸超壓的相似律[J].火炸藥學報,2008,33(4):32-35.ZHONG Qian,WANG Bo-liang,HUANG Ju,et al.Study on the similarity law of TNT explosion overpressure in air[J].Chinese Journal of Explosives and Propellants,2008,33(4):32-35.
[7] 劉偉,鄭毅,秦飛.近地面TNT 爆炸的試驗研究和數值模擬[J].爆破,2012,29(1):5-9.LIU Wei,ZHENG Yi,QIN Fei.Experimental and numerical simulation of TNT explosion on the ground[J].Blasting,2012,29(1):5-9.
[8] 周保順,張立恒,王少龍,等.TNT 炸藥爆炸沖擊波的數值模擬與實驗研究[J].彈箭與制導學報,2010,30(3):88-90.ZHOU Bao-shun,ZHANG Li-heng,WANG Shaolong,et al.Numerical simulation and experimental research on TNT explosion shock wave[J].Journal of Projectiles,Rockets,Missiles and Guidance,2010,30(3):88-90.
[9] 李加貴,邊小華,張雷.爆炸沖擊波傳播的數值模擬與試驗數據對比[J].山西建筑,2006,32(8):106-107.LI Jia-gui,BIAN Xiao-hua,ZHANG Lei.Numerical simulation of blast wave propagation in tunnel compared with experiment data[J].Shanxi Architecture,2006,32(8):106-107.
[10]張廣福,劉玉存,王建華.爆炸沖擊波無限空氣領域傳播的數值模擬研究[J].山 西化工,2009,29(1):43-46.ZHANG Guang-fu,LIU Yu-cun,WANG Jian-hua.Numerical simulation study on propagation of shock wave in the air without obstacles[J].Shanxi Chemical Industry,2009,29(1):43-46.
[11]楊鑫,石少卿,程鵬飛.空氣中TNT 爆炸沖擊波超壓峰值的預測及數值模 擬[J].爆破,2008,25(1):15-19.YANG Xin,SHI Shao-qing,CHENG Peng-fei.Forecast and simulation of peak overpressure of TNT explosion shock wave in the air[J].Blasting,2008,25(1):15-19.
[12]李秀地.T 型坑道中爆炸沖擊波傳播規律的數值模擬[J].后勤工程學院學報,2011,27(4):8-12.LI Xiu-di.Numerical simulation for blast propagation and attenuation inside T-shaped tunnel from hecharges detonation[J].Journal of Logistical Engineering University,2007,30(4):55-57.
[13]盧紅琴,劉偉慶.空中爆炸沖擊波的數值模擬研究[J].武漢理工大學學報,2009,31(19):105-108.LU Hong-qin,LIU Wei-qing.Research on numerical simulation of blast wave in air[J].Journal of Wuhan University of Technology,2009,31(19):105-108.
[14]石磊,杜修力,樊鑫.爆炸沖擊波數值計算網格劃分方法研究[J].北京工業大學學報,2010,36(11):1465-1470.SHI Lei,DU Xiu-li,FAN Xin.A study on the mesh generation method for numerical simulation of blast wave[J].Journal of Beijing University of Technology,2010,36(11):1465-1470.
[15]GB6722-2003 爆破安全規程[S].北京:中國工程爆破協會,2003.
[16]張陶,惠君明,解立峰,等.FAE 爆炸場超壓與威力的實驗研究[J].爆炸與沖擊,2004,24(2):176-181.ZHANG Tao,HUI Jun-ming,XIE Li-feng,et al.Experimental research on the overpressure and power in the FAE blast field[J].Explosion Shock Waves,2004,24(2):176-181.
[17]董桂旭,杜茂華,黃雪峰.某型炸藥的沖擊波超壓峰值計算公式參數的修正[J].海軍航空工程學院學報,2010,25(5):542-544.DONG Gui-xu,DU Mao-hua,HUANG Xue-feng.Correction on parameters of calculation formula for shockwave overpressure peak value of one explosive[J].Journal of Naval Aeronautical and Astronautical University,2010,25(5):542-544.
[18]胡華權,裴明敬,許學忠,等.燃料空氣炸藥爆炸威力評價方法研究[C]∥第四屆全國爆炸力學實驗技術學術會議論文.合肥:安徽省力學學會,2006:313-318.HU Hua-quan,PEI Ming-jing,XU Xue-zhong,et al.Study on the evaluation method of fuel-air explosive power[C]∥The Fourth Session of the National Explosion Mechanics Experiment Technology of Academic Debate.Hefei:Anhui Society of Theoretical and Applied Mechanics,2006:313-318.
[19]郭煒,俞統昌,金朋剛.三波點的測量與實驗技術研究[J].火炸藥學報,2007,30(4):55-57.GUO Wei,YU Tong-chang,JIN Peng-gang.Test of triple point and study on its test technology[J].Chinese Journal of Explosives and Propellants,2007,30(4):55-57.
[20]王建靈,郭煒,馮曉軍.TNT、PBX 和Hexel空中爆炸沖擊波參數的實驗研究[J].火炸藥學報,2008,31(6):42-44.WANG Jian-ling,GUO Wei,FENG Xiao-jun.Experimental research on the air explosion shock wave parameters of TNT,PBX and Hexel[J].Chinese Journal of Explosives and Propellants,2008,31(6):42-44.
[21]高軒能,王書鵬.大空間柱殼結構爆炸動力響應的Ritz-POD數值模擬[J].土木建筑與環境工程,2010,32(2):64-70.GAO Xuan-neng,WANG Shu-peng.Numerical simulation for dynamic response of large-space cylindrical reticulated shell under internal explosion by Ritz-POD method[J].Journal of Civil,Architectural and Environmental Engineering,2010,32(2):64-70.
[22]姜宏,劉鵬.爆炸荷載作用下鋼框架動力響應有限元分析[J].湖南科技大學學報(自然科學版),2010,25(2):55-58.JIANG Hong,LIU Peng.Fea dynamic analysis of steel frame under explosive load[J].Journal of Hunan University of Science and Technology(Natural Science Edition),2010,25(2):55-58.
[23]田志敏,鄔玉斌,羅奇峰.隧道內爆炸沖擊波傳播特性及爆炸荷載分布規律研究[J].振動與沖擊,2011,30(1):21-26.TIAN Zhi-min,WU Yu-bin,LUO Qi-feng.Characteristics of in-tunnel explosion-induced air shockwave and distribution law of reflected shock wave load[J].Journal of Vibration and Shock,2011,30(1):21-26.
[24]Katayama M,Kibe S,Yamamoto T.Numerical and experimental study on the shaped charge for space debris assessment[J].Acta Astronauttca,2001,48(5-12):363-372.
[25]廖維張,杜修力,田志敏.爆炸荷載作用下部分埋置結構響應的數值模擬方法[J].北京工業大學學報,2007,33(2):155-159.LIAO Wei-Zhang,DU Xiu-li,TIAN Zhi-min.Numerical simulation methods on dynamic response of partially buried structure under blast loading[J].Journal of Beijing University of Technology,2007,33(2):155-159.
[26]尚曉江.ANSYS/LS-DYNA 動力分析方法與工程實例[M].北京:中國水利水電出版社,2005.