TNT空中爆炸沖擊波的工程和數(shù)值計算

辛春亮，王俊林，余道建，李 通，宋師軍

（1. 北京航天長征飛行器研究所，北京，100076；2. 中國運(yùn)載火箭技術(shù)研究院，北京，100076）

0 引 言

炸藥在空氣中爆炸后瞬間形成高溫高壓的爆炸產(chǎn)物，產(chǎn)物強(qiáng)烈壓縮周圍靜止的空氣，在空氣中形成沖擊波向四周傳播，對結(jié)構(gòu)造成破壞。許多學(xué)者對于TNT空中爆炸沖擊波超壓和傳播規(guī)律進(jìn)行了研究[1～5]，總結(jié)出了 Brode、Henrych、Kingery-Bulmash、Kinney-Graham等超壓經(jīng)驗和半經(jīng)驗公式，并編寫了CONWEP、BEC、ATBLAST、BLAPAN、3D-BLAST、EBLAST、SHOCK、BEAM、BLASTX等工程計算程序，這些超壓計算公式和工程計算程序是在總結(jié)了大量試驗數(shù)據(jù)的基礎(chǔ)上發(fā)展起來的，大都將炸藥按爆熱折算成等效TNT當(dāng)量，不考慮負(fù)壓區(qū)和沖擊波的多次反射以及稀疏作用，并假設(shè)壓力以指數(shù)規(guī)律衰減，因此主要適用于自由場環(huán)境下比例距離較大時超壓的快速計算。數(shù)值模擬可以解決超壓計算公式和工程計算程序難以解決的更為復(fù)雜的問題，這方面的研究成果也很多[6～9]，大都采用基于有限差分、有限體積法（如 AUTODYN、Air3D、CTH、DYTRAN、DYSMAS、IFSAS、SHAMRC）以及有限元法軟件（如 LS-DYNA、EUROPLEXUS）中的顯式積分算法，數(shù)值計算超壓與試驗測試結(jié)果的吻合程度也各不相同，大體來說，基于有限差分和有限體積法的軟件計算準(zhǔn)確度較高，有限元軟件則差一些，其計算超壓曲線的顯著特點是：a）峰值之前有明顯的壓力爬升過程；b）超壓峰值容易被抹平，峰值過后壓力衰減過快；c）網(wǎng)格尺寸越小，計算結(jié)果越準(zhǔn)確。

本文將空氣沖擊波正壓區(qū)和負(fù)壓區(qū)工程計算公式及其參數(shù)結(jié)合起來，形成TNT空氣自由場爆炸沖擊波工程計算模型，并采用有限差分軟件AUTODYN和有限元軟件LS-DYNA對無限空間TNT空氣中爆炸沖擊波傳播過程進(jìn)行數(shù)值模擬，最后將三者得出的計算結(jié)果進(jìn)行對比分析。

1 空氣沖擊波工程計算模型

典型的TNT空氣中爆炸沖擊波曲線如圖1所示。

圖1 自由場空氣沖擊波壓力-時間曲線Fig.1 Schematic Profile of Airblast Shock Wave in Free Air

由圖1可知，在沖擊波到達(dá)之前，該處的壓力等于大氣壓力0Ρ，沖擊波在時間at到達(dá)該處后，壓力瞬間由大氣壓力突躍至最大值。壓力最大值與0Ρ的差值，通常稱為入射超壓峰值iΡ。波陣面通過后壓力迅速下降，經(jīng)過時間dt壓力衰減到大氣壓力并繼續(xù)下降，直至出現(xiàn)負(fù)超壓峰值，在一定時間內(nèi)又逐漸地回升到大氣壓力。負(fù)超壓與0Ρ的差值，通常稱為負(fù)超壓峰值nΡ。

空氣沖擊波壓力在正壓段大致按指數(shù)規(guī)律衰減，有許多經(jīng)驗公式可以描述此壓力衰減過程，其中修正的 Friedlander方程[2]較接近實際過程且又簡單易于計算：

1984年，Kingery在對大量試驗數(shù)據(jù)總結(jié)分析的基礎(chǔ)上提出了 Kingery-Bulmash空氣沖擊波參數(shù)計算模型[4,10,11]，其研究成果已被廣泛采用，并被植入多種計算程序如CONWEP和LS-DYNA（即*LOAD_BLAST關(guān)鍵字）中。經(jīng)過一系列的試驗測試，1998年Kingery對Kingery-Bulmash超壓模型作了一些修改[12,13]，遠(yuǎn)場超壓預(yù)測值比以前高出許多，此修正模型已在BEC工程計算程序中實現(xiàn)。

Kingery沖擊波正壓區(qū)參數(shù)計算公式：

式中 FUNCTION代表沖擊波到達(dá)時間ta、入射超壓Ρi及正壓作用時間 td等參數(shù)；Z為比例距離， Z = d /；d為距離爆心的距離；W為裝藥質(zhì)量。當(dāng)比例距離Z＞40 m·kg-1/3時，對于大多數(shù)結(jié)構(gòu)已無破壞能力。A，B，C，D，E，F(xiàn)，G是系數(shù)，具體取值見表1。

表1 簡化的Kingery空氣沖擊波參數(shù)Tab.1 Shock Wave Parameters of Simplified Kingery Formula

對于衰減系數(shù)b，Martin Larcher[9]提出了下面的計算公式：

負(fù)壓區(qū)沖擊波參數(shù)雖然對剛性結(jié)構(gòu)（鋼筋混凝土）的設(shè)計不重要，但對柔性防護(hù)結(jié)構(gòu)（一般指鋼框架）尤其重要。Martin Larcher[9]采用雙折線方程來計算負(fù)壓區(qū)壓力：

式中0Ρ為大氣壓力；nΡ為負(fù)超壓；nt為負(fù)超壓持續(xù)時間。

對于負(fù)壓區(qū)參數(shù)（負(fù)超壓及其持續(xù)時間），Martin Larcher[9]對 Krauthammer[14]提供的圖表進(jìn)行擬合，得出下列計算公式。

負(fù)超壓計算公式：

負(fù)超壓持續(xù)時間計算公式：

將上述空氣沖擊波正壓區(qū)和負(fù)壓區(qū)工程計算公式結(jié)合起來，便可形成TNT空氣自由場爆炸沖擊波工程計算模型，該模型可以快速計算出不同TNT當(dāng)量在不同距離處的沖擊波壓力曲線。

2 數(shù)值計算模型

由于炸藥爆炸初期產(chǎn)生的沖擊波是高頻波，在數(shù)值計算模型中炸藥及其附近區(qū)域需要劃分細(xì)密網(wǎng)格才能反映出足夠頻寬的沖擊波特性，否則計算出的壓力峰值會被抹平。Martin Larcher[9]建議，對于1 kg TNT空氣中爆炸數(shù)值模擬，炸藥及其附近空氣的網(wǎng)格尺寸在1 mm左右，如此細(xì)密的網(wǎng)格劃分方式會導(dǎo)致網(wǎng)格數(shù)量極其龐大。無限空間TNT空中爆炸問題具有球?qū)ΨQ性質(zhì)，一維計算模型是最佳選擇，可顯著降低計算規(guī)模和計算時間。

AUTODYN軟件中的二維楔形網(wǎng)格軸對稱計算模型可用來等效一維計算模型，垂直于沖擊波傳播方向只有一個二維四邊形單元。而 LS-DYNA可采用一維梁單元球?qū)ΨQ計算模型，如圖2所示。在計算模型中采用多物質(zhì)Euler算法和缺省的人工粘性系數(shù)1×10-6，時間步長縮放因子均為0.9。

圖2 AUTODYN和LS-DYNA計算模型Fig.2 AUTODYN and LS-DYNA Computational Models

計算模型中的炸藥為 1 kg TNT，計算空氣域為12 m，網(wǎng)格劃分總數(shù)為6000，距離爆心1 m內(nèi)的網(wǎng)格尺寸為1 mm，然后由小到大逐漸向外圍漸變。

在上述計算模型中，空氣密度為1.225 kg/m3，采用理想氣體狀態(tài)方程，γ=1.4。空氣中初始壓力為1個標(biāo)準(zhǔn)大氣壓（0.1 MPa）。

TNT炸藥爆炸產(chǎn)物壓力用JWL狀態(tài)方程來描述。

式中 E為單位質(zhì)量內(nèi)能；V為比容；A，B，1R，2R，ω為常數(shù)。其中，方程式右端第1項在高壓段起主要作用，第2項在中壓段起主要作用，第3項代表低壓段。在爆轟產(chǎn)物膨脹的后期，方程式前兩項的作用可以忽略，為了加快求解速度，將炸藥從JWL狀態(tài)方程轉(zhuǎn)換為更為簡單的理想氣體狀態(tài)方程（絕熱指數(shù)γ=ω+1）。TNT炸藥爆炸產(chǎn)物JWL狀態(tài)方程參數(shù)取值來自文獻(xiàn)[15]。TNT爆炸產(chǎn)物JWL狀態(tài)方程參數(shù)如表2所示。

表2 TNT爆炸產(chǎn)物JWL狀態(tài)方程參數(shù)Tab.2 JWL EOS Parameters of TNT Explosion Products

3 工程計算和數(shù)值計算結(jié)果的對比

圖3是距離爆心不同距離處沖擊波壓力計算曲線。圖3中AUTODYN和LS-DYNA數(shù)值計算曲線上均存多個峰值，第2個峰值是由TNT炸藥的爆心匯聚反射追趕造成的。由于空氣密度和壓強(qiáng)遠(yuǎn)小于爆轟產(chǎn)物的密度和壓強(qiáng)，在沖擊波形成的同時由界面向炸藥中心反射回一個稀疏波，使爆轟產(chǎn)物發(fā)生膨脹，降低內(nèi)部壓力，此稀疏波在炸藥中心匯聚后又向外傳播一個壓縮波，由于前導(dǎo)沖擊波已經(jīng)將空氣絕熱壓縮，此壓縮波的傳播速度將大于前導(dǎo)沖擊波，并逐漸向前追趕前導(dǎo)沖擊波。如果測點離炸藥不遠(yuǎn)，此二次壓力波峰值足夠高，就有可能將負(fù)壓區(qū)強(qiáng)行打斷，并再次衰減到負(fù)壓。

圖3 不同距離處沖擊波壓力-時間計算曲線Fig.3 Comparition of Airblast Shock Wave Profiles of Different Ranges

由圖3可知，AUTODYN和LS-DYNA數(shù)值計算曲線極為相似，走勢與工程計算曲線基本一致。距爆心越遠(yuǎn)，3條曲線符合程度越高，且數(shù)值計算曲線提供的信息更加豐富。

數(shù)值計算曲線與工程計算曲線的區(qū)別主要在于：

a）沖擊波到達(dá)時間：工程計算曲線沖擊波峰值最高，因此最早到達(dá)，其次是 AUTODYN，LS-DYNA計算曲線沖擊波到達(dá)時間最晚。由于數(shù)值計算模型中采用了細(xì)密網(wǎng)格和很小的人工粘性系數(shù)，計算壓力時間曲線爬升段接近強(qiáng)間斷。

b）正壓峰值：工程計算值最高，其次是AUTODYN計算值，LS-DYNA計算值最低，且距離爆心越近，也就是比例距離越小，差別越明顯。

c）正壓作用時間：工程計算值略高于數(shù)值計算值。

d）正壓區(qū)沖量：由 b）和 c）基本可以得出，工程計算值高于數(shù)值計算值，比例距離越小，差別越明顯。

e）負(fù)壓峰值：工程計算值與數(shù)值計算峰值基本相當(dāng)。負(fù)壓峰值到達(dá)時間提前。

f）負(fù)壓峰值到達(dá)時間：AUTODYN計算曲線負(fù)壓峰值到達(dá)時間最早，其次是LS-DYNA，工程計算曲線負(fù)壓峰值最晚到達(dá)。

g）負(fù)壓作用時間：數(shù)值計算峰值與工程計算值基本相當(dāng)。

對于不同的比例距離，上述規(guī)律均相同。此外，由于工程計算模型簡單，考慮因素少，耗時短，瞬間即可完成計算。LS-DYNA計算耗時 4 min 52 s。AUTODYN軟件本身計算速度比LS-DYNA慢，再加上二維四邊形單元比一維梁單元計算耗費(fèi)大，因此AUTODYN計算耗時最長，為34 min，是LS-DYNA的7倍。

4 結(jié) 論

a）本文將把 Kingery、Bulmash、Martin Larcher以及Krauthammer等人在空氣沖擊波正壓區(qū)和負(fù)壓區(qū)參數(shù)計算研究成果結(jié)合起來，形成TNT空氣自由場爆炸沖擊波工程計算模型，通過該模型可以快速計算出TNT空氣自由場爆炸沖擊波壓力衰減演化曲線。該模型的正壓區(qū)參數(shù)計算公式由大量實驗數(shù)據(jù)擬合而成，準(zhǔn)確可靠。由于過于簡單地采用雙折線模型將負(fù)壓區(qū)持續(xù)時間均分，與實際會存在較大偏差。

b）AUTODYN數(shù)值計算曲線與工程計算曲線吻合程度最高，負(fù)壓區(qū)包含的信息比工程計算曲線更為豐富，也更接近實際。

c）LS-DYNA計算曲線與工程計算曲線吻合程度稍差，但LS-DYNA計算時間遠(yuǎn)低于AUTODYN。

[1] Brode H L. Numerical solutions of spherical blast wave[J]. JAP, 1955,26(6)： 766-775.

[2] Baker W E. Explosions in air[M]. Austin： University of Texas Press, 1973.

[3] Henrych J. The dynamics of explosion and its use[M]. Amsterdam：Elsevier, 1979.

[4] Kingery, Charles N, Bulmash, Gerald. Airblast parameters from TNT spherical air burst and hemispherical surface burst[R]. Maryland： Defence Technical Information Center, Ballistic Research Laboratory, Aberdeen Proving Ground, 1984.

[5] Kinney G F, Graham K J. Explosive shocks in air (2ndEdition)[M]. New York： Springer-Verlag, 1985.

[6] Borgers J B. Vantomme J. Towards a parametric model of a planar blastwave created with detonating cord[C]. Ralston： 19th military aspects of blast and shock, 2006.

[7] Clutter J K. Stahl M. Hydrocode simulations of air and water shocks for facility vulnerability assessments[J]. Journal of Hazardous Materials,2004(106)： 9-24.

[8] Alia A, Souli M. High explosive simulation using multi-material formulations[J]. Applied Thermal Engineering, 2006, 26(10)：1032-1042.

[9] Martin L. Simulation of the effects of an air blast wave[R]. Ispra：JRC41337, 2007.

[10] Army Engineer Waterways Experiment Station. Fundamentals of protective design for conventional weapons[R]. Vicksburg： US Department of the Army, WES/IR/SL-1, 1987.

[11] Conwep Software, Hyde D W. Conventional weapon effects[R].Albuquerque： US Army Engineering Waterways Experimental Station,919167, 1992.

[12] TB 700-2, NAVSEAINST 8020.8 B. DoD ammunition and explosives[R].Washington： Hazards Classification Procedures, DOD-4145.26-M-REV,1998.

[13] Swisdak M M. Simplified kingery airblast calculations[C]. Florida：Proceedings of the 26th DoD Explosives Safety Seminar, 1994.

[14] Krauthammer T, Altenberg A. Negative phase blast effects on glass panels[J]. International Journal of Impact Engineering, 2000, 24(1)： 1-18.

[15] Dobratz B M, Crawford P C. LLNL explosives handbook - properties of chemical explosives and explosive stimulants[M]. California： Lawrence Livermore National Laboratory, UCRL-52997, 1985.