周杰,陶鋼,張洪偉,潘保青

(1.南京理工大學(xué)能源與動(dòng)力工程學(xué)院,江蘇南京 210094;2.北京跟蹤與通信技術(shù)研究所,北京 100094)

模擬不同爆炸沖擊波的計(jì)算方法研究

周杰1,陶鋼1,張洪偉2,潘保青2

(1.南京理工大學(xué)能源與動(dòng)力工程學(xué)院,江蘇南京 210094;2.北京跟蹤與通信技術(shù)研究所,北京 100094)

給出了一種可獲得不同爆炸沖擊波的方法。利用SolidWorks軟件建立了沖擊波生成區(qū)域與沖擊波加載空氣域的幾何模型,采用Hypermesh軟件對(duì)模型進(jìn)行映射網(wǎng)格劃分,并賦予材料參數(shù)和狀態(tài)方程,通過設(shè)計(jì)Fortran程序,將所需沖擊波的壓力-時(shí)間曲線,轉(zhuǎn)換為單位體積內(nèi)能-時(shí)間的曲線,應(yīng)用修正公式對(duì)曲線進(jìn)行校正,然后加載到?jīng)_擊波生成區(qū)域內(nèi),最后通過LS-DYNA有限元程序計(jì)算獲得所需要的平面沖擊波場(chǎng)。采用該計(jì)算方法可以模擬獲得不同爆炸沖擊波,誤差范圍在1.0%之內(nèi)。該方法為深入研究爆炸沖擊波創(chuàng)傷效應(yīng)提供了一種有效的技術(shù)途徑。

爆炸力學(xué);沖擊波;數(shù)值模擬;沖擊波創(chuàng)傷

0 引言

不同爆炸源(猛炸藥爆炸、氣體源爆炸及不同規(guī)模的爆炸等)依照其具體情況會(huì)產(chǎn)生不同效果的爆炸沖擊波,其爆炸超壓、正壓作用時(shí)間和沖量等參數(shù)均不同。在研究爆炸沖擊波的殺傷和毀傷效應(yīng)時(shí),沖擊波的超壓、正壓作用時(shí)間和沖量等基本參數(shù)的組合效應(yīng)決定了創(chuàng)傷效果[1-3]。著名的Bowen創(chuàng)傷曲線[4]及沖量-壓力毀傷曲線[5],都是以爆炸沖擊波的超壓值和正壓作用時(shí)間或沖量值為創(chuàng)傷指示參數(shù)。因此,如何利用有限元計(jì)算程序模擬不同爆炸沖擊波,是爆炸創(chuàng)傷領(lǐng)域內(nèi)比較關(guān)注的問題。

爆炸實(shí)驗(yàn)中,炸藥尺寸、炸藥密度、起爆位置和環(huán)境等因素決定了爆炸沖擊波的波形,且隨著目標(biāo)的空間位置不同,其波形特征也要發(fā)生變化[6-7]。由于實(shí)驗(yàn)中爆炸沖擊波受多種因素制約,因此常規(guī)的數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)存在著一定的差異。常用的LS-DYNA顯式動(dòng)力學(xué)分析軟件,一般有兩種產(chǎn)生沖擊波的方法:一種是利用炸藥庫中的猛炸藥(TNT)模擬理想爆炸,選用MAT_HIGH_EXPLOSIVE_BURN材料模型和JWL狀態(tài)方程,通過炸藥起爆獲得沖擊波[8],然而此方法獲得的沖擊波的超壓峰值和正壓作用時(shí)間等參數(shù)是不可調(diào)控的;另一種方法是采用LOAD_BLAST,在實(shí)體表面(Lagrangian網(wǎng)格)加載沖擊波的壓力-時(shí)間曲線[8],這種方法雖然可以獲得計(jì)算需要的沖擊波波形參數(shù),卻不可以進(jìn)行流-固耦合計(jì)算分析,所以應(yīng)用范圍十分狹窄。

Greer[9]和Thom[10]通過約束計(jì)算模型外邊界上節(jié)點(diǎn)的運(yùn)動(dòng)方向和旋轉(zhuǎn)方向,采用理想氣體的狀態(tài)方程,將溫度曲線加載到?jīng)_擊波生成區(qū)域內(nèi),利用LS-DYNA計(jì)算獲得有效的沖擊波壓力曲線。但是計(jì)算時(shí)需要約束所有外邊界上節(jié)點(diǎn)的運(yùn)動(dòng)方向和旋轉(zhuǎn)方向,否則計(jì)算將無法進(jìn)行。因?yàn)橥膺吔缟系墓?jié)點(diǎn)受到約束,沖擊波將無法透射出邊界,所以不能模擬無限大區(qū)域的沖擊波傳播特征。文獻(xiàn)[9-10]中沒有提供直接得到復(fù)雜沖擊波的方法,因此不夠完整。針對(duì)這些問題,本文給出了一種改進(jìn)的方法,可以計(jì)算得到所需的沖擊波。

1 模型與方法

通過在沖擊波生成區(qū)域內(nèi)加載單位體積內(nèi)能-時(shí)間(E-t)曲線和設(shè)定空氣域的邊界條件(反射和透射邊界均可),計(jì)算獲得沖擊波。同時(shí)該計(jì)算模型也能應(yīng)用于流-固耦合計(jì)算,克服了LS-DYNA中一般方法在計(jì)算分析爆炸毀傷問題上存在的不足。

1.1 計(jì)算模型與邊界條件

采用Solidworks軟件建立的幾何模型由沖擊波生成區(qū)域和沖擊波加載空氣域(以下簡(jiǎn)稱:空氣域)構(gòu)成。為了節(jié)約計(jì)算時(shí)間,空氣區(qū)域的大小為3 m× 1.5 m×0.04 m,而沖擊波生成區(qū)域的尺寸為3 m× 0.01 m×0.04 m.通過Hypermesh軟件對(duì)幾何模型進(jìn)行映射網(wǎng)格(六面體網(wǎng)格)劃分,網(wǎng)格單元的尺寸為0.01 m,空氣域的網(wǎng)格數(shù)為180 000,網(wǎng)格總數(shù)為180 120,如圖1所示。

圖1 計(jì)算模型Fig.1 Calculation model

圖1為模擬沖擊波的計(jì)算模型。設(shè)定空氣域的外邊界A、B、C為透射邊界(也可通過約束節(jié)點(diǎn)的運(yùn)動(dòng)和旋轉(zhuǎn)方向,設(shè)定為剛性邊界),邊界允許流體介質(zhì)流出空氣域,以便模擬無限大空間的效應(yīng)?？刂茮_擊波的傳播方向的方法:約束計(jì)算模型的上下表面(OXY平面)的所有節(jié)點(diǎn)(包括空氣域和沖擊波生成區(qū)域)在Z軸的運(yùn)動(dòng)方向與繞X、Y軸旋轉(zhuǎn)的方向;沖擊波生成區(qū)域與空氣域的接觸平面(OXZ平面)共節(jié)點(diǎn),并約束了沖擊波生成區(qū)域余下的OXY、OXZ、OYZ平面上節(jié)點(diǎn)的所有運(yùn)動(dòng)和旋轉(zhuǎn)方向。

1.2 模擬沖擊波的生成原理和方法

Greer[9]和Thom[10]建立的計(jì)算模型包括:高溫區(qū)域和空氣域。計(jì)算模型采用理想氣體狀態(tài)方程pV=nRT,當(dāng)相對(duì)體積V為常量時(shí),壓力p與溫度T呈線性關(guān)系。將所需沖擊波的壓力-時(shí)間(p-t)曲線轉(zhuǎn)化為溫度-時(shí)間(T-t)曲線,將曲線加載到高溫區(qū)域,可計(jì)算獲得p-t曲線。但是由于此方法需要約束所有外邊界上節(jié)點(diǎn)的運(yùn)動(dòng)方向和旋轉(zhuǎn)方向,這樣處理會(huì)導(dǎo)致沖擊波無法透射出邊界,也就不能模擬無限大區(qū)域效應(yīng);而采用氣體線性多項(xiàng)式的狀態(tài)方程則不需要對(duì)所有的外邊界節(jié)點(diǎn)進(jìn)行約束。故采用氣體線性多項(xiàng)式的狀態(tài)方程代替理想氣體狀態(tài)方程,既能獲得所需沖擊波,又能夠有效地控制邊界條件。根據(jù)氣體線性多項(xiàng)式狀態(tài)方程計(jì)算沖擊波E-t曲線(設(shè)定相對(duì)體積為常量),可以將曲線直接加載到?jīng)_擊波生成區(qū)域,用來模擬不同爆炸類型的沖擊波。

因?yàn)闇y(cè)量實(shí)驗(yàn)中,通常所使用的壓電傳感器無法準(zhǔn)確測(cè)量負(fù)壓,所以本文著重研究沖擊波的正壓相。Friedlander沖擊波的正壓相特征:壓力迅速升高的壓力峰值p+隨后以指數(shù)衰減到環(huán)境壓力。Friedlander沖擊波正壓相的數(shù)學(xué)方程式[11]為

式中:p(t)為沖擊波超壓值的時(shí)間函數(shù);p0為環(huán)境壓力;p+為沖擊波正壓相的超壓峰值;t+為正壓作用時(shí)間;b為衰減常數(shù)。超壓峰值和作用時(shí)間是由實(shí)驗(yàn)得來的,衰減常數(shù)根據(jù)爆炸的類型而定。

計(jì)算模型中的空氣域和沖擊波生成區(qū)域均采用氣體的線性多項(xiàng)式的狀態(tài)方程[8]:

式中:p為沖擊波的超壓值;E為單位體積內(nèi)能;μ= 1/V-1,V為相對(duì)體積?？諝庥蚝蜎_擊波生成區(qū)域均采用空氣材料模型,具體參數(shù)如表1所示。

表1 沖擊波生成區(qū)域和空氣域的參數(shù)Tab.1 Loading area of shock wave and air basin parameters

將表1中的參數(shù)代入到(2)式中,得到?jīng)_擊波生成區(qū)域和空氣域的狀態(tài)方程:

沖擊波生成區(qū)域內(nèi),設(shè)定了相對(duì)體積V為常量,所以爆轟壓力p與單位體積內(nèi)能E呈線性關(guān)系。然而空氣域中狀態(tài)方程的相對(duì)體積V卻不為常量,所以空氣域中的沖擊波壓力p與單位體積內(nèi)能E呈非線性關(guān)系。

Friedlander理想沖擊波的模擬方法:利用Fortran程序編制Friedlander沖擊波的p-t曲線。采用(3)式將沖擊波的p-t曲線,轉(zhuǎn)化為E-t曲線,最后將曲線加載到?jīng)_擊波生成區(qū)域的K文件中;利用LSDYNA有限元程序,即可在空氣域內(nèi)計(jì)算獲得所需的Friedlander理想沖擊波。

復(fù)雜沖擊波的模擬方法:采用(3)式,將復(fù)雜沖擊波的實(shí)驗(yàn)數(shù)據(jù)(p-t曲線)轉(zhuǎn)化為E-t曲線,最后將曲線加載到?jīng)_擊波生成區(qū)域的K文件中;利用LSDYNA有限元程序,即可在空氣域內(nèi)模擬獲得到復(fù)雜沖擊波。

2 計(jì)算結(jié)果

2.1 計(jì)算實(shí)例

針對(duì)Friedlander理想沖擊波,表2提供了所需要模擬的沖擊波參數(shù),其中超壓峰值和正壓作用時(shí)間參數(shù)已知。

表2 加載的沖擊波參數(shù)Tab.2 Parameters of loaded shock waves

因?yàn)樵贚S-DYNA軟件中,使用的基本單位是厘米-克-微秒制,所以設(shè)定標(biāo)準(zhǔn)狀況下空氣域的單位體積的初始內(nèi)能為E0為25 N·cm/(g·cm3).數(shù)值計(jì)算中沖擊波在空氣中傳播的規(guī)律,與邊界條件和計(jì)算區(qū)域的大小有關(guān);如果需要在離沖擊波生成區(qū)域一定的目標(biāo)位置產(chǎn)生沖擊波,則需要依據(jù)計(jì)算模型的大小、邊界條件和沖擊波的傳播規(guī)律,對(duì)加載的單位體積內(nèi)能-時(shí)間曲線進(jìn)行相應(yīng)的校正。圖2是模擬沖擊波特征的計(jì)算結(jié)果,其中壓力監(jiān)控單元位于空氣域中。

圖2 計(jì)算結(jié)果Fig.2 Calculated results

由表2參數(shù),采用Friedlander理想沖擊波(1)式計(jì)算獲得沖擊波的p-t曲線,通過(3)式將沖擊波的p-t曲線轉(zhuǎn)換為E-t曲線,如圖3所示;將曲線加載到計(jì)算模型中的沖擊波生成區(qū)域中,利用LS-DYNA計(jì)算獲得空氣域內(nèi)沖擊波p-t曲線,如圖4所示。

圖3 加載的內(nèi)能-時(shí)間曲線Fig.3 Loaded internal energy vs.time

圖4 沖擊波的波形Fig.4 Waveforms of shock waves

通過對(duì)計(jì)算結(jié)果的分析發(fā)現(xiàn):空氣域內(nèi)壓力監(jiān)控單元計(jì)算獲得的沖擊波正壓作用時(shí)間與表2提供的沖擊波正壓作用時(shí)間保持一致,如圖3與圖4所示;但是空氣域內(nèi)計(jì)算獲得的沖擊波超壓峰值p+與表2提供的沖擊波超壓峰值存在一定的差異,誤差范圍在5.3%～24%區(qū)間,如表3所示。

表3 沖擊波超壓峰值的誤差分析Tab.3 Error analysis on overpressures of shock waves

計(jì)算產(chǎn)生誤差的原因:網(wǎng)格的質(zhì)量(數(shù)量、單元尺寸、形狀)還不夠高;有部分能量轉(zhuǎn)化為氣體的動(dòng)能;沖擊波由生成區(qū)域向空氣域中傳播時(shí)會(huì)導(dǎo)致能量的損失。因此在實(shí)際計(jì)算過程中需要對(duì)加載的E-t曲線進(jìn)行必要的修正,采用的修正公式為

式中:EM(t)為修正后的單位體積內(nèi)能;E(t)為理論計(jì)算單位體積的內(nèi)能;E0為單位體積的初始內(nèi)能; A為修正系數(shù),在0.3～0.03區(qū)間。圖5、圖6分別為修正后的單位體積內(nèi)能-時(shí)間(EM-t)曲線和修正后計(jì)算獲得的沖擊波波形。修正后的沖擊波超壓峰值與表2提供的沖擊波超壓峰值的誤差控制在1%范圍之內(nèi)。

圖5 修正后加載的內(nèi)能-時(shí)間曲線Fig.5 Modified curves of loaded internal energy-time

圖6 修正后沖擊波的波形Fig.6 Modified waveforms of shock waves

通常在復(fù)雜環(huán)境(如墻壁或有限空間)中由于沖擊波反射和折射而產(chǎn)生混合疊加的鋸齒形波形,形成復(fù)雜沖擊波。本文提供的方法也適用于復(fù)雜沖擊波的模擬。圖7為復(fù)雜沖擊波的實(shí)驗(yàn)曲線,將實(shí)驗(yàn)測(cè)得的沖擊波壓力-時(shí)間曲線,轉(zhuǎn)換為E-t曲線,并根據(jù)修正公式校正曲線,最后將EM-t曲線加載到?jīng)_擊波生成區(qū)域內(nèi),采用LS-DYNA可計(jì)算獲得復(fù)雜沖擊波波形。數(shù)值模擬獲得的復(fù)雜沖擊波波形如圖8所示,壓力監(jiān)控單元如圖2所示。

圖7 復(fù)雜沖擊波的實(shí)驗(yàn)波形Fig.7 Experimental waveforms of complex shock waves

圖8 數(shù)值模擬獲得的復(fù)雜沖擊波波形Fig.8 Simulated waveforms of complex shock waves

通過分析結(jié)果得出:模擬得到的沖擊波波形與沖擊波的實(shí)驗(yàn)波形基本一致;由于復(fù)雜沖擊波的實(shí)驗(yàn)數(shù)據(jù)只測(cè)量到3.1 ms,所以在模擬沖擊波的計(jì)算過程中,在3.1 ms后LS-DYNA軟件自動(dòng)將超壓線性衰減到環(huán)境壓力。因此,本方法可以有效地模擬復(fù)雜沖擊波,模擬得到的波形與實(shí)驗(yàn)波形基本一致。

2.2 應(yīng)用實(shí)例

利用本方法,調(diào)控超壓峰值和正壓作用時(shí)間,可以加載不同爆炸類型的沖擊波,這對(duì)爆炸沖擊波毀傷的研究有著較重要意義。圖9為修正的Bowen創(chuàng)傷曲線[9]。

圖9 修正的Bowen創(chuàng)傷曲線Fig.9 Modified Bowen curve

參照修正的Bowen創(chuàng)傷曲線,根據(jù)人體胸前入射沖擊波超壓峰值p+和正壓作用時(shí)間t+,可獲得創(chuàng)傷程度。圖9中a點(diǎn)的沖擊波作用于人體時(shí),人處于臨界創(chuàng)傷程度;a點(diǎn)的超壓峰值p+為150 kPa,正壓作用時(shí)間t+為3 ms,通過本文提供的方法,可模擬a點(diǎn)的沖擊波波形,如圖10所示。如在沖擊波加載區(qū)域中加入人體模型,就可進(jìn)一步研究不同類型的沖擊波與人體作用的創(chuàng)傷機(jī)理。

圖10 a點(diǎn)的沖擊波波形Fig.8 Waveform of shock wave on Point a

以上所提供的沖擊波算例,只是針對(duì)沖擊波的正壓相。這是由于一般所使用的壓電傳感器無法準(zhǔn)確測(cè)量負(fù)壓相,因?yàn)樵趯?shí)際操作和參數(shù)測(cè)量中,如何確定負(fù)壓對(duì)創(chuàng)傷的影響,存在很多問題。采用這里提供的數(shù)值方法,可以很容易形成所需要的負(fù)壓相,可研究其對(duì)人體目標(biāo)的殺傷效應(yīng)。這方面的內(nèi)容在以后的研究中將給予關(guān)注。

3 結(jié)論

1)本文使用氣體線性多項(xiàng)式的狀態(tài)方程代替理想氣體狀態(tài)方程,既可以獲得所需的沖擊波,又能正確的處理邊界條件,改進(jìn)了文獻(xiàn)[9-10]方法的不足。

2)給出平面沖擊波方法,可以根據(jù)實(shí)驗(yàn)數(shù)據(jù)或Friedlander沖擊波p-t曲線,通過空氣線性多項(xiàng)式的狀態(tài)方程,把p-t曲線轉(zhuǎn)變?yōu)镋M-t曲線,并加載到有限元程序中,可計(jì)算獲得精度較高的模擬結(jié)果,誤差范圍在1%之內(nèi)。

3)利用該方法可以有效地模擬爆炸沖擊波(簡(jiǎn)單波或復(fù)雜波),為研究不同爆炸類型的沖擊波創(chuàng)傷機(jī)理提供了有效的技術(shù)途徑。

4)通過算例,給出了對(duì)應(yīng)于Bowen創(chuàng)傷曲線的某一臨界創(chuàng)傷條件所需要的沖擊波p-t曲線,解決了現(xiàn)有其他數(shù)值模擬方法不能精確模擬實(shí)際沖擊波的p-t曲線問題。

References)

[1] Stuhmiller J H.Biological response to blast overpressure:a summary of modeling[J].Toxicology,1997,121(33):91-103.

[2] Stuhmiller J H,Ho K,Vorst M.A model of blast overpressure injury to the lung[J].Journal of Biomechanics,1996,29(2): 227-234.

[3] Bowen I G,Fletcher E R,Richmond D R.Estimate of man's tolerance to the direct effects of air blast,DASA-2113[R]. Washington,DC:Defense Atomic Support Agency,1968.

[4] Argyros G J.Management of primary blast injury[J].Toxicology, 1997,121(1):105-115.

[5] Fernando D A,Enrique G F,Marta D M,et al.Consequence analysis by means of characteristic curves to determine the damage to humans from the detonation of explosive substances as a function of TNT equivalence[J].Journal of Loss Prevention in the Process Industries,2008,20(1):74-81.

[6] Johansson C H,Person P A.Demonic of high explosives[M]. New York:Academic Press,1970.

[7] Kinney G F.Explosive shocks in air[M].New York:The Macmillan Company,1962.

[8] LS-DYNA key word user's manual[EB].Livermore,CA:Livermore Software Technology Corporation,2003.

[9] Greer A D.Numerical modeling for the prediction of primary blast injury to the lung[D].Canada:University of Waterloo,2005.

[10] Thom C G.Soft materials under air blast loading and their effect on primary blast injury[D].Canada:,University of Waterloo, 2009.

[11] Wu C,Hao H.Modeling of simultaneous ground shock and air blast pressure on nearby structures from surface explosions[J]. International Journal of Impact Engineering,2005,31(6): 699-717.

Research on Simulating Method of Different Blast Waves

ZHOU Jie1,TAO Gang1,ZHANG Hong-wei2,PAN Bao-qing2
(1.School of Energy and Power Engineering,Nanjing University of Science and Technology,Nanjing 210094,Jiangsu,China; 2.Beijing Institute of Tracking and Telecommunications Technology,Beijing 100094,China)

A method to acquire different blast waves is proposed.The geometrical models of the loading area of shock wave and air basin are established using SolidWorks,and then are divided to the structured grids by using Hypermesh.Material parameters and equation of state are assigned to the model.The curve of internal energy per unit volume and time is calculated with FORTRAN program and revised with the modified formula,then imported into the loading area of shock wave.The required shock wave is acquired by use of finite element program LS-DYNA.The proposed method can used to simulate the different blast waves with error of 1.0%.The method provides an effective route for research on blast injury.

explosion mechanics;shock wave;numerical simulation;blast injury

TJ011.1

A

1000-1093(2014)11-1846-05

10.3969/j.issn.1000-1093.2014.11.016

2013-09-22

周杰(1986—),男,博士研究生。E-mail:Beijihu1986@163.com;

陶鋼(1962—),男,教授,博士生導(dǎo)師。E-mail:taogang@mail.njust.edu.cn