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

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

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

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

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

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

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

爆炸力學;沖擊波;數值模擬;沖擊波創傷

0 引言

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

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

Greer[9]和Thom[10]通過約束計算模型外邊界上節點的運動方向和旋轉方向,采用理想氣體的狀態方程,將溫度曲線加載到沖擊波生成區域內,利用LS-DYNA計算獲得有效的沖擊波壓力曲線。但是計算時需要約束所有外邊界上節點的運動方向和旋轉方向,否則計算將無法進行。因為外邊界上的節點受到約束,沖擊波將無法透射出邊界,所以不能模擬無限大區域的沖擊波傳播特征。文獻[9-10]中沒有提供直接得到復雜沖擊波的方法,因此不夠完整。針對這些問題,本文給出了一種改進的方法,可以計算得到所需的沖擊波。

1 模型與方法

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

1.1 計算模型與邊界條件

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

圖1 計算模型Fig.1 Calculation model

圖1為模擬沖擊波的計算模型。設定空氣域的外邊界A、B、C為透射邊界(也可通過約束節點的運動和旋轉方向,設定為剛性邊界),邊界允許流體介質流出空氣域,以便模擬無限大空間的效應。控制沖擊波的傳播方向的方法:約束計算模型的上下表面(OXY平面)的所有節點(包括空氣域和沖擊波生成區域)在Z軸的運動方向與繞X、Y軸旋轉的方向;沖擊波生成區域與空氣域的接觸平面(OXZ平面)共節點,并約束了沖擊波生成區域余下的OXY、OXZ、OYZ平面上節點的所有運動和旋轉方向。

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

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

因為測量實驗中,通常所使用的壓電傳感器無法準確測量負壓,所以本文著重研究沖擊波的正壓相。Friedlander沖擊波的正壓相特征:壓力迅速升高的壓力峰值p+隨后以指數衰減到環境壓力。Friedlander沖擊波正壓相的數學方程式[11]為

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

計算模型中的空氣域和沖擊波生成區域均采用氣體的線性多項式的狀態方程[8]:

式中:p為沖擊波的超壓值;E為單位體積內能;μ= 1/V-1,V為相對體積。空氣域和沖擊波生成區域均采用空氣材料模型,具體參數如表1所示。

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

將表1中的參數代入到(2)式中,得到沖擊波生成區域和空氣域的狀態方程:

沖擊波生成區域內,設定了相對體積V為常量,所以爆轟壓力p與單位體積內能E呈線性關系。然而空氣域中狀態方程的相對體積V卻不為常量,所以空氣域中的沖擊波壓力p與單位體積內能E呈非線性關系。

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

復雜沖擊波的模擬方法:采用(3)式,將復雜沖擊波的實驗數據(p-t曲線)轉化為E-t曲線,最后將曲線加載到沖擊波生成區域的K文件中;利用LSDYNA有限元程序,即可在空氣域內模擬獲得到復雜沖擊波。

2 計算結果

2.1 計算實例

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

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

因為在LS-DYNA軟件中,使用的基本單位是厘米-克-微秒制,所以設定標準狀況下空氣域的單位體積的初始內能為E0為25 N·cm/(g·cm3).數值計算中沖擊波在空氣中傳播的規律,與邊界條件和計算區域的大小有關;如果需要在離沖擊波生成區域一定的目標位置產生沖擊波,則需要依據計算模型的大小、邊界條件和沖擊波的傳播規律,對加載的單位體積內能-時間曲線進行相應的校正。圖2是模擬沖擊波特征的計算結果,其中壓力監控單元位于空氣域中。

圖2 計算結果Fig.2 Calculated results

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

圖3 加載的內能-時間曲線Fig.3 Loaded internal energy vs.time

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

通過對計算結果的分析發現:空氣域內壓力監控單元計算獲得的沖擊波正壓作用時間與表2提供的沖擊波正壓作用時間保持一致,如圖3與圖4所示;但是空氣域內計算獲得的沖擊波超壓峰值p+與表2提供的沖擊波超壓峰值存在一定的差異,誤差范圍在5.3%～24%區間,如表3所示。

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

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

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

圖5 修正后加載的內能-時間曲線Fig.5 Modified curves of loaded internal energy-time

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

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

圖7 復雜沖擊波的實驗波形Fig.7 Experimental waveforms of complex shock waves

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

通過分析結果得出:模擬得到的沖擊波波形與沖擊波的實驗波形基本一致;由于復雜沖擊波的實驗數據只測量到3.1 ms,所以在模擬沖擊波的計算過程中,在3.1 ms后LS-DYNA軟件自動將超壓線性衰減到環境壓力。因此,本方法可以有效地模擬復雜沖擊波,模擬得到的波形與實驗波形基本一致。

2.2 應用實例

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

圖9 修正的Bowen創傷曲線Fig.9 Modified Bowen curve

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

圖10 a點的沖擊波波形Fig.8 Waveform of shock wave on Point a

以上所提供的沖擊波算例,只是針對沖擊波的正壓相。這是由于一般所使用的壓電傳感器無法準確測量負壓相,因為在實際操作和參數測量中,如何確定負壓對創傷的影響,存在很多問題。采用這里提供的數值方法,可以很容易形成所需要的負壓相,可研究其對人體目標的殺傷效應。這方面的內容在以后的研究中將給予關注。

3 結論

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

2)給出平面沖擊波方法,可以根據實驗數據或Friedlander沖擊波p-t曲線,通過空氣線性多項式的狀態方程,把p-t曲線轉變為EM-t曲線,并加載到有限元程序中,可計算獲得精度較高的模擬結果,誤差范圍在1%之內。

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

4)通過算例,給出了對應于Bowen創傷曲線的某一臨界創傷條件所需要的沖擊波p-t曲線,解決了現有其他數值模擬方法不能精確模擬實際沖擊波的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—),男,教授,博士生導師。E-mail:taogang@mail.njust.edu.cn