空氣自由場裝藥爆炸荷載一維建模數(shù)值算法

楊允鍇， 武雙章， 高振儒， 毛益明

（ 陸軍工程大學(xué)野戰(zhàn)工程學(xué)院， 江蘇 南京 210007）

0 引言

當下結(jié)構(gòu)與構(gòu)件的抗爆性能及防護研究任務(wù)仍舊繁巨，然而結(jié)構(gòu)與構(gòu)件爆炸試驗成本昂貴，且安全風(fēng)險大，因此數(shù)值模擬已成為爆炸研究領(lǐng)域的有效手段。準確可靠的爆炸荷載計算，是分析結(jié)構(gòu)爆炸響應(yīng)和毀傷評估的基礎(chǔ)與關(guān)鍵。如何利用數(shù)值模擬軟件準確計算爆炸荷載并模擬空氣沖擊波傳播過程，是進行結(jié)構(gòu)抗爆研究的首要問題。

關(guān)于自由空氣中TNT 炸藥爆炸荷載的計算，已有不同學(xué)者得出一些預(yù)測方法及經(jīng)驗公式， 雖然已經(jīng)廣泛應(yīng)用于結(jié)構(gòu)爆炸響應(yīng)分析和構(gòu)件損傷評估的研究中， 但采用不同方法計算結(jié)果存在較大差異。 在有限元分析軟件LS-DYNA 中，爆炸荷載的施加主要有三種方式：第一種是三角形沖擊荷載，該方法計算效率最高，建模簡單，但此法無法模擬爆炸荷載的非均布加載， 也無法模擬空氣沖擊波與構(gòu)件相互作用的問題， 因此得到的誤差較大； 第二種是CONWEP 法，即通過*LOAD_BLAST、*LOAD_SEGMENT_SET 和*DEFINE_SEGMENT 等關(guān)鍵字來施加一些特定的爆炸荷載，此法可以略過炸藥與空氣建模步驟，簡化建模流程、提高計算效率， 但無法考慮空氣沖擊波反射及繞射等現(xiàn)象，僅適于遠距離爆炸下結(jié)構(gòu)的動態(tài)響應(yīng)分析；第三種是流固耦合法，建立空氣域和炸藥模型，直接模擬炸藥爆炸過程、空氣沖擊波傳播過程以及沖擊波荷載加載過程，因此考慮了空氣沖擊波與目標結(jié)構(gòu)的相互作用， 效果最為準確直觀，缺點在于網(wǎng)格尺寸、初始條件對計算精度影響大，建模復(fù)雜、計算效率低[1]。 姚成寶[2]通過研究得出保證三維流固耦合法計算精度的兩個結(jié)論： ①為保證爆炸荷載計算精度，空氣網(wǎng)格劃分應(yīng)小于2mm；②空氣域模型的尺寸至少應(yīng)為構(gòu)件尺寸的2 倍以上， 以消除無反射邊界條件的影響。因計算機性能限制，以上兩點結(jié)論很難在結(jié)構(gòu)與構(gòu)件毀傷數(shù)值模擬研究中實現(xiàn)。

為解決采用三維流固耦合法開展構(gòu)件毀傷數(shù)值模擬時計算效率低下、計算誤差大的問題，首先利用有限元分析軟件LS-DYNA 建立裝藥爆炸與空氣沖擊波一維仿真模型， 計算了6.82g 和1000g 球形TNT 裝藥在自由空氣中爆炸的空氣沖擊波超壓值， 并分別與經(jīng)驗公式計算結(jié)果、三維模型計算結(jié)果[3]進行了對比，結(jié)果表明利用一維仿真模型計算爆炸荷載計算效率與計算精度更高。 再按照“兩步法”將一維仿真模型的計算結(jié)果映射到三維，進一步模擬空氣沖擊波傳播并毀傷構(gòu)件的過程。 實際算例表明，“兩步法” 能夠在保證計算效率與計算精度的前提下，模擬空氣沖擊波與構(gòu)件的相互作用過程，在結(jié)構(gòu)與構(gòu)件抗爆領(lǐng)域具有廣泛的應(yīng)用前景。

1 爆炸荷載工程算法

關(guān)于球形TNT 裝藥自由空中的爆炸荷載計算，已有學(xué)者進行了大量的相關(guān)研究， 并獲得了可信度較高的成果。爆炸荷載主要用沖擊波壓力、峰值超壓、比沖量、持續(xù)時間、到達時間等沖擊波參數(shù)來描述，其中沖擊波超壓是研究結(jié)構(gòu)與構(gòu)件毀傷的重要參數(shù)。 沖擊波超壓通常用比例距離來表達。 比例距離可以定義為：

貝克（W.E.Baker）[6]通過實驗得到的經(jīng)驗公式是：

式中：F—沖擊波到達時間ta， 自由空中爆炸時沖擊波的峰值超壓ΔPm及正壓作用時間t+等參數(shù)；A、B、C、D、E、F、G—對應(yīng)的實參數(shù)系數(shù)，具體取值見表1；其余參數(shù)與式（2）相同。

表1 簡化的Kingery 空氣沖擊波參數(shù)Tab.1 Shock wave parameters of simplified kingery formula

為了將不同經(jīng)驗公式的預(yù)測結(jié)果進行直觀的對比，在對數(shù)坐標系中給出各經(jīng)驗公式預(yù)測的峰值超壓值隨比例距離的變化關(guān)系，見圖1。

圖1 不同經(jīng)驗公式?jīng)_擊波峰值超壓-比例距離對比圖Fig.1 Comparison of peak pressure attenuation against scaled distance

從圖中可以看出，在比例距離較小時，式（3）～（5）的計算結(jié)果接近，而式（2）因超出公式適用范圍，計算結(jié)果偏差偏大；隨著比例距離增大，式（2）～（4）結(jié)果趨于一致，但式（5）主要用于為結(jié)構(gòu)抗爆設(shè)計提供參考，計算結(jié)果偏大。 通過比較，式（3）整體范圍內(nèi)計算結(jié)果均處于可信范圍，因此本文將亨利奇公式作為推薦公式，將其計算結(jié)果與仿真計算結(jié)果進行對比。

2 一維仿真模型

選用大型通用顯式有限元分析軟件LS-DYNA 建立球形TNT 裝藥在自由空中爆炸的一維仿真模型，模擬空氣沖擊波傳播并計算爆炸荷載。

2.1 物理模型

在建立一維模型時，使用BEAM161 梁單元建立TNT炸藥和空氣物理模型， 模型采用cm-g-μs 單位制。 為便于將計算結(jié)果與三維模型[3]進行對比，TNT 藥量確定為6.82g 和1000g 兩種工況。

針對2 種不同工況，分別建立相應(yīng)的一維仿真模型。工況一的TNT 裝藥模型長度1cm，空氣域模型長度200cm；工況二中TNT 裝藥模型長度5.3cm， 空氣域模型長度2000cm。TNT 與空氣兩種材料之間采用共節(jié)點接觸，并劃分為1mm均勻網(wǎng)格， 分別在比例距離為0.1、0.5、1、2、3、4、5、6、7、8m/kg1/3位置處設(shè)置追蹤點，物理模型見圖2。

圖2 一維模型示意圖Fig.2 The one-dimensional model

2.2 材料模型與狀態(tài)方程

TNT 炸藥選用008 號材料模型*MAT_HIGH_EXPLOSIVE_BURN，其具體材料參數(shù)在表2[9]中給出。爆炸力學(xué)中用狀態(tài)方程來描述爆轟產(chǎn)物系統(tǒng)中各物理量（壓力、體積、溫度等）之間的關(guān)系，在進行一維數(shù)值模擬時選用*EOS_JWL 狀態(tài)方程， 該方程定義壓力為相對比容V 和爆轟產(chǎn)物比內(nèi)能e 的函數(shù)：

表2 反演驗證計算結(jié)果

表2 TNT 炸藥材料模型參數(shù)Tab.2 Parameters of TNT model parameters

表3 JWL 狀態(tài)方程參數(shù)Tab.3 Parameters of the*EOS_JWL

采用009 號材料模型*MAT_NULL 以及線性多項式狀態(tài)方程*EOS_LINEAR_POLYNOMIAL[9]描述空氣的材料特性，線性多項式狀態(tài)方程表述為：

式中：c0～c6—與空氣性質(zhì)有關(guān)的常數(shù)；pa—空氣的壓力；μ—相對密度，μ=ρ/ρ0-1，ρ、ρ0分別為空氣的密度、 初始密度；E0，a—空 氣 的 初 始 體 積 內(nèi) 能。 空 氣 的 初 始 密 度 取0.00125g·cm-3，狀態(tài)方程中具體參數(shù)取值見表4[10]，其中V0為初始的相對體積。

表4 LINEAR_POLYNOMIAL 狀態(tài)方程參數(shù)Tab.4 Parameters of the*EOS_LINEAR_POLYNOMIAL

2.3 仿真計算與數(shù)據(jù)分析

通過*SECTION_ALE_1D 定義一維模型的歐拉算法，并使用Van Leer+HIS（二階精度）方法進行計算，以應(yīng)對炸藥與空氣網(wǎng)格的嚴重畸變。 2 種不同工況的求解時間 分 別 為30000μs 和90000μs。 利 用 后 處 理 軟 件Ls-Prepost 對仿真結(jié)果進行處理和分析， 得到各追蹤點位置處空氣沖擊波的計算結(jié)果， 并繪制了工況二中比例距離和比例距離處的沖擊波壓力時程曲線，分別如圖3（a）和圖3（b）。

圖3 工況二不同位置處壓力時程曲線Fig.3 Overpressure response time history of condition two

將2 種工況條件下數(shù)值模擬的計算結(jié)果與亨利奇公式的計算結(jié)果進行對比， 并得到其相對誤差， 具體見表5?？梢钥闯? 種工況下數(shù)值模擬計算結(jié)果與經(jīng)驗公式預(yù)測值均存在一定誤差，且比例距離時誤差較大。主要原因在于推導(dǎo)經(jīng)驗公式的實驗數(shù)據(jù)誤差。 空氣中TNT 爆炸沖擊波近場參數(shù)測量難， 早先的工程算法研究的試驗裝置靈敏度和精度不夠理想，造成試驗數(shù)據(jù)誤差，使得經(jīng)驗公式預(yù)測值存在一定偏差。

2.3.1 經(jīng)驗公式計算值對比

通過對表5 空氣沖擊波壓力值計算結(jié)果的對比分析，可以發(fā)現(xiàn)：數(shù)值仿真結(jié)果符合爆炸相似原理；2 種工況的最大相對誤差值均出現(xiàn)在r=0.1m/kg1/3處， 可見不宜采用數(shù)值仿真模擬構(gòu)件在裝藥近距離爆炸時的毀傷效果；當比例距離r>0.1m/kg1/3時，各追蹤點的數(shù)值仿真計算值與亨利奇公式計算值的誤差全部均小于20%之內(nèi)；裝藥量增加能夠減小數(shù)值模擬計算值與經(jīng)驗公式預(yù)測值的相對誤差。

表5 不同比例距離處空氣沖擊波超壓一維仿真模型計算值與經(jīng)驗公式計算值對比Tab.5 Commanded overpressure comparison with simulated overpressure

通過對比可以看出，當比例距離r>0.1m/kg1/3時，一維仿真模型計算的空氣沖擊波峰值超壓誤差較小， 符合工程使用要求。

2.3.2 三維模型仿真計算值對比

楊鑫[3]運用LS-DYNA 軟件建立三維模型計算球形裝藥在自由空氣中爆炸產(chǎn)生的沖擊波超壓值， 其三維仿真模型使用的裝藥藥量與本文完全一致， 分別為6.82g和1000g 兩種工況，并計算了比例距離0.1、0.5、1、2、3、4、5、6、7、8m/kg1/3處的超壓峰值。 兩種工況下三維仿真模型計算結(jié)果與經(jīng)驗公式計算值的相對誤差范圍分別為32.2%～94.8%和25.7%～92.9%，而一維仿真模型的兩種工況的誤差范圍分別為-61.3%～12.3%和-43.0%～19.2%。

通過對比可以看出， 一維仿真模型對空氣沖擊波峰值壓力的計算結(jié)果與經(jīng)驗公式的預(yù)測值相對誤差更小，一維模型求解空氣沖擊波荷載參數(shù)的計算精度更高。

3 一維模型計算結(jié)果的映射應(yīng)用

一維仿真模型無法直接模擬空氣沖擊波與結(jié)構(gòu)構(gòu)件的相互作用過程，但可以使用“兩步法”將一維數(shù)值模擬結(jié)果映射到三維流固耦合模型中，此時一維模型計算得到的空氣沖擊波荷載將以三維形式完成對構(gòu)件的毀傷模擬。

第一步：在進行一維仿真計算時，首先應(yīng)設(shè)置合適求解時間， 使得空氣沖擊波在求解時間內(nèi)傳播距離大致為裝藥到達構(gòu)件迎爆面的距離， 然后在模型求解時生成相應(yīng)的映射（map）文件。

第二步：建立空氣域及構(gòu)件的三維模型，構(gòu)件迎爆面方向的空氣域尺寸可以減小， 且空氣域網(wǎng)格可進行稀疏劃分，能夠提高計算效率且對荷載計算精度的影響很小。通過*INITIAL_ALE_MAPPING 以及*DEFINE_VICTOR 關(guān)鍵字將map 文件映射至三維模型中， 并定義映射起始點及方向，此時空氣沖擊波荷載將在三維模型中繼續(xù)傳播。

以裝藥爆炸毀傷鋼板為例，驗證兩步法的可行性。將第2 節(jié)中1000g 球形TNT 裝藥爆炸的一維仿真結(jié)果，映射到建立的空氣和鋼板所在三維模型， 其中空氣域尺寸1200mm×1200mm×500mm，鋼板尺寸1000mm×1000mm×2mm，兩者采用流固耦合法模擬相互作用關(guān)系。 映射方向為鋼板迎爆面法線的負方向， 映射點為鋼板中心點在空氣域表面的正投影點，三維模型及映射關(guān)系見圖4。

圖4 三維模型及映射關(guān)系圖Fig.4 3D model andmapping relationship

由圖5 可以看出， 沖擊波實現(xiàn)了在三維模型中的映射與傳播， 并對構(gòu)件造成了毀傷。 可見“兩步法”在保證沖擊波荷載計算精度和計算效率的前提下，準確模擬沖擊波與構(gòu)件的相互作用過程。

圖5 鋼板毀傷圖Fig.5 Damage effect of steel plate

4 結(jié)論

通過建立一維模型進行球形TNT 炸藥自由空氣中爆炸工況的數(shù)值模擬研究， 對比經(jīng)驗公式計算值和前人所做三維數(shù)值模擬結(jié)果，得出以下結(jié)論：

現(xiàn)有的空中TNT 爆炸沖擊波荷載的經(jīng)驗公式有很多， 但不同學(xué)者的經(jīng)驗公式預(yù)測結(jié)果存在一定的適用范圍及差別。研究近爆作用下構(gòu)件毀傷效應(yīng)時，亨利奇公式的預(yù)測結(jié)果更加準確。

通過對比一維仿真模型、 三維仿真模型獲得的空氣沖擊波超壓值與經(jīng)驗公式預(yù)測值的相對誤差， 可知一維仿真模型的計算結(jié)果誤差更小， 即采用一維模型求解爆炸空氣沖擊波荷載時計算精度更高。 其原因在于數(shù)值模擬計算裝藥爆炸荷載時， 網(wǎng)格尺寸能夠顯著影響求解精度，網(wǎng)格尺寸越大，誤差就越大，網(wǎng)格尺寸小于2mm 時計算精度能夠得到保證； 而三維模型的網(wǎng)格尺寸減小會導(dǎo)致計算效率的大幅下降，受計算機性能的限制，實際應(yīng)用中難以實現(xiàn)?！皟刹椒ā?是將一維模型計算所得高精度空氣沖擊波荷載映射到三維模型中，并模擬其與構(gòu)件的相互作用的一種方法。 利用“兩步法”能夠把一維模型建模簡單、計算精度高、計算效率高的優(yōu)點，與三維流固耦合模型真實模擬空氣沖擊波與構(gòu)件相互作用過程的優(yōu)點相結(jié)合，為結(jié)構(gòu)與構(gòu)件爆炸毀傷研究提供了一種高效的數(shù)值模擬方法。