柱狀炸藥自由空氣爆炸的網格劃分方法研究

林文樂，王 碩，嚴 波

(國防科技大學 空天科學學院， 長沙 410072)

1 引言

網格劃分方法對空氣爆炸的沖擊波峰值數值模擬精度有較大影響，且不同炸藥當量下其網格尺寸取值亦有較大區別。爆炸沖擊波峰值的數值精度隨著網格尺寸的減小而增大，然而網格尺寸過小會導致單元數目呈指數倍增長，極大程度地降低了計算效率。因此，需要提出一種既綜合考慮計算精度和效率，又對不同當量TNT均有較強適用性的網格劃分方法，為大型結構采用流固耦合方法進行抗爆分析提供空氣域網格劃分方法的參考。

炸藥在空氣中爆炸后瞬間產生強烈的爆炸沖擊波，是爆炸對周圍目標產生破壞效應的重要因素之一，自20世紀40年代以來學者們便開展了廣泛的研究。柱狀炸藥作為主要的裝藥形式之一，其在空氣中的傳播規律吸引了眾多學者的研究目光。索強等利用試驗數據，比較了不同網格劃分比對爆炸沖擊波數值結果的影響，結果表明當網格劃分比為炸藥體邊長的1/32時，數值結果與試驗結果的誤差在5%以內；張社榮等采用AUTODYN軟件對不同網格劃分比的水下爆炸傳播規律進行了研究，結果表明采用炸藥半徑的1/3作為網格尺寸可滿足工程精度要求，并給出了不同網格劃分比下的誤差擬合公式；石磊等采用LS-DYNA軟件研究了網格劃分方法對爆炸沖擊波的影響，結果表明網格尺寸的漸進比例和漸進比例起始位置對計算結果有較大影響；姚成寶等利用LS-DYNA軟件研究了無反射邊界條件對沖擊波參量的影響，結果表明無反射邊界尺寸為考察區域2倍時數值結果與實驗結果相似；陳銘等采用LS-DYNA軟件對圓柱裝藥的徑高比和起爆方式進行了研究，結果表明徑高比對比沖量影響較大；高軒能等采用LS-DYNA軟件分析了炸藥材料參數、TNT藥量、網格尺寸、空氣域形狀和炸藥形狀等對沖擊波參量的影響；Luccioni等利用流體力學軟件研究了爆炸沖擊波的網格尺寸效應，結果表明10 cm的網格尺寸便可較為精確的模擬沖擊波的傳播規律。

由此可見，學者們對爆炸沖擊波的數值模擬研究較為全面，但是不同學者在研究中提出的網格劃分方法存在較大差別，而網格劃分方法對計算效率和計算精度的影響較大。本文采用LS-DYNA軟件，對比分析不同炸藥當量下網格劃分方法對計算精度及效率的影響，提出了一種對不同當量均有較強適用性的網格劃分方法。

2 基于經驗公式的超壓峰值

不同類型的炸藥爆炸產生的能量可用等效當量TNT炸藥來衡量。當爆炸發生在自由空氣時，其傳播規律由峰值超壓、沖量、持續時間等參量來描述。目前，最常用的實驗方法是采用Hopkinson比例定律，以小當量爆炸來研究大當量爆炸的特性，比例距離的表達式如下：

=13

(1)

式中，為測點與爆心之間的距離，為等效TNT藥量。

目前，具有較高權威性的實驗成果為UFC 3-340-02規范中給出的設計圖表；除此之外，其他學者以及我國的人防規范亦給出了相應的經驗公式，如圖1所示，不同學者提出的經驗公式結果存在較大差別，尤其是在靠近爆心差別越顯著。圖2所示，Henrych和baker提出的經驗公式與UFC 3-340-02規范中提出的圖表較為一致，當比例距離大于1時，各個公式預測的結果比較接近。文獻[12]表明：當沖擊波傳播距離>時(為測點距軸心的距離，為柱形炸藥的高度)，可近似看成球形裝藥的爆炸。為此，本研究采用Henrych公式對大當量柱狀炸藥的超壓數值結果進行驗證。

圖1 不同經驗公式下的沖擊波超壓峰值

圖2 Henrych、Baker與UFC 3-340-02超壓峰值

3 柱形炸藥的數值計算及對比分析

3.1 材料的狀態方程及參數

空氣采用*MAT_NULL模型，應用線性多項式狀態方式*EOS_LINEAR_POLUNOMIAL描述，表達式如下：

(2)

=1

(3)

式中：、、、、、、為常數，為單位初始體積內能，為初始相對體積。參數取值如表1所示。

表1 空氣的材料參數

炸藥采用*MAT_HIGH_EXPLOSIVE_BURN材料模型，定義壓力為相對體積和初始能量的函數，表達式如下：

(4)

式中:為炸藥的相對體積，為單位體積內能，、、、為表征炸藥材料特性的常數，上式的三項由左至右分別代表峰值超壓的高壓段、中壓段及低壓段。文獻[14-16]采用的TNT炸藥參數在CJ狀態下滿足約束守恒方程組，且能較好地模擬結構在爆炸工況下的動態響應過程，其參數取值如表2所示。

表2 TNT炸藥的材料參數

3.2 網格劃分方法對數值結果的影響

文獻[17]已驗證了建立1/8模型進行爆炸沖擊波研究的可行性。文獻[5]通過定義相對厚度=，其中為數值模型的最短邊尺寸，為柱狀裝藥的半徑，研究發現當>2時，無反射邊界對計算結果影響較小。文獻[3]研究表明網格尺寸小于炸藥半徑3倍以上時，其水中爆炸數值結果滿足工程精度要求。

本研究采用LS-DYNA軟件建立TNT爆炸的數值計算模型。為了提高計算效率，本研究采用1/8長方體模型，單元類型取8節點SOLID164單元，空氣域尺寸取為4 m×4 m×13m(高×寬×長)，通過關鍵字*CONSTRAINED_GLOBAL及*BOUNDARY_NON_REFLECTING定義模型的對稱邊界及無反射邊界，材料及狀態方程參數見表1和表2，以研究網格劃分方法對爆點在原點、長徑比為3的800kgTNT圓柱炸藥的峰值超壓的影響；通過關鍵字*DATABASE_ELOUT和*DATABASE_HISTORY_SOLID輸出以爆炸中心為起點，沿空氣域長度方向比例距離分別為0.1～1.1 m/kg、時間間隔為1的峰值超壓曲線，計算總時間為0.01 s。

3.2.1 網格劃分比對數值結果的影響

網格劃分比指“網格尺寸與炸藥半徑之比”?？紤]爆炸近區的沖擊波的衰減性較大，對比例距離小于0.3的爆炸區域進行加密(以爆炸中心為角點的3 m×3 m×3 m長方體區域)，加密區段網格尺寸分別為4 cm、5 cm、7.5 cm、10 cm,非加密區段網格尺寸均為10 cm。

圖3表明：當比例距離小于0.4時，數值結果與Henrych經驗公式相差較大；圖4表明：當比例距離大于0.4時，4 cm與5 cm網格尺寸計算結果與經驗公式吻合度較好。不同網格尺寸下的峰值超壓與Henrych經驗公式的誤差分析如表3所示，當比例距離大于0.4時，4 cm、5 cm、7.5 cm和10 cm網格尺寸的最大誤差分別為10.1%、-12.5%、-31.8%和-31%，單元數目分別為1 264 375、735 000、325 000和208 000。結果表明：當比例距離大于0.4時，峰值超壓的計算精度隨著網格尺寸的減小而增加，4 cm與5 cm網格尺寸的數值結果精度相差不大，誤差均在10%左右，但4 cm網格的單元數目5 cm網格的1.7倍，綜合考慮計算效率和計算精度的關系，加密區采用5 cm的網格尺寸。

圖3 網格尺寸對比例距離小于0.4的超壓峰值曲線

圖4 網格尺寸對比例距離大于0.4的超壓峰值曲線

表3 網格尺寸對超壓峰值計算誤差的影響(%)

3.2.2 加密范圍對數值結果的影響

以3.2.1節的結果為依據，進一步探討加密范圍分別為4 m、5 m、6 m對峰值超壓結果的影響。網格劃分方法如下：以爆炸中心為角點，加密區域分別為4/5/6 m×4/5/6 m×4/5/6 m的長方體區域。若空氣域某一方向小于該加密長度，取空氣域該方向的最大長度作為加密長度，加密區網格尺寸取為5 cm，非加密區網格尺寸均為10 cm。

如圖5所示，圖例中如1-0.05-3表示加密范圍為3 m且加密范圍為5 cm的等比例網格。圖5表明：當比例距離小于0.4時，增大加密范圍對數值結果與Henrych經驗公式的誤差影響不大；圖6表明：當比例距離大于0.4時，增大加密范圍，能在一定程度上提高數值結果的精度。如表4所示，當比例距離大于0.4時，峰值超壓隨著加密范圍的增大而增大，當加密區長度由3 m增加至4 m時，最大誤差由-12.5%減少至-4.4%，誤差在5%以內，精度有了明顯提高；加密區長度進一步由4 m增大至5 m和6 m時，對計算精度的提高影響不明顯。由此可見，800 kg TNT當量下，加密區長度在4 m以內發生改變對精度有明顯的影響，據此推斷出其加密界限為0.4 m/kg。

圖5 加密范圍對比例距離小于0.4的超壓峰值曲線

表4 加密范圍對超壓峰值計算誤差的影響(%)

3.2.3 漸進比例對數值結果的影響

以3.2.2節研究結果為依據，為了提高計算效率、節省計算時間，進一步探討加密范圍內漸進比例網格對峰值超壓結果的影響。網格劃分方法如下：加密范圍取4 m，加密區網格尺寸取5 cm，加密范圍網格尺寸分別以1.01、1.02、1.03的漸進比例系數沿三個方向過渡到非加密范圍，非加密網格尺寸取10 cm。為了控制非加密區網格尺寸不超過10 cm和加密范圍為4 m，漸進系數為1.01的漸進網格由爆心開始，漸進系數為1.02和1.03的漸進網格分別從距離爆心三個方向的2.3 m和2.6 m開始，其單元數目分別為532 800、630 784、1 000 936。

如圖7、圖8所示，圖例中如“1.01-0.05”代表加密區漸進比例系數為1.01，網格尺寸為5 cm的數值模型。圖7表明當比例距離小于0.4時，漸進比例系數對近區沖擊波超壓影響不大；圖8表明當比例距離大于0.4時，3種漸進比例系數都與經驗公式有較好的吻合效果。表5展示了不同漸進比例系數下數值結果與Henruch的計算誤差，可以看出，3種漸進比例系數下，當比例距離大于0.4時，誤差都在10%左右，且當比例距離大于0.5時，誤差在5%以內，計算精度較高。此外，當漸進系數為1.01時，單元數目較其他2種漸進比例少，約為等比例網格劃分方式的0.5倍，在計算機處理器為Intel(R) Xeon(R) CPU E5-2620 v4 @ 2.10 GHz、內存為16 GB、LS-DYNA計算版本為ls971_s_R5.1.1等條件下，計算時間為23 min，較2.2.2節的加密區等比例網格劃分方法縮短了約1/2計算時間。

圖7 漸進比例對比例距離小于0.4的超壓峰值曲線

圖8 漸進比例對比例距離大于0.4的超壓峰值曲線

表5 漸進比例對超壓峰值計算誤差的影響(%)

3.3 網格劃分方法對不同當量TNT的適用性

網格劃分方法的結果為：加密范圍不小于0.4 m/kg，加密區網格劃分比不大于2/15、非加密區網格劃分比不大于1/4或10 cm(取兩者的較小值)，加密區漸進比例為1.01。根據以上研究成果，分別驗證該網格劃分方法對500 kg和1 500 kg當量，長徑比為3的柱狀TNT炸藥的適用性。500 kg當量TNT和1 500 kg當量TNT的網格劃分方法如下：

1) 500 kg當量?？諝庥虺叽鐬?.5 m×3.5 m×11 m(高×寬×長)，加密范圍為3.5 m，加密區網格取4 cm，由爆心開始以1.01漸進比例沿三個方向過渡到非加密區，非加密區網格尺寸均為8 cm。

2) 1 500 kg當量。空氣域尺寸為4.5 m×4.5 m×15 m(高×寬×長)，加密范圍為4.5 m，加密區網格取6 cm，以1.01漸進比例系數沿三個方向過渡到非加密區，非加密區網格尺寸均為10 cm。

如圖9所示，當比例距離小于0.4時，該網格劃分方法所得的峰值超壓與Henrych經驗公式的誤差較大；當比例距離大于0.4時，500 kg當量和1 500 kg當量TNT均與Henrych經驗公式吻合度較好。由此驗證了當比例距離大于0.4時，該網格劃分方法對不同當量TNT均有較強的適用性。

圖9 500 kg/1 500 kg當量的超壓峰值曲線

4 與文獻提出的網格劃分方法的對比

文獻[2-4]提出的網格劃分方法表現為：在全局范圍內劃分等比例網格，網格劃分比分別為1/32、3/80和1/3。如表6所示，索強和石磊提出的網格劃分方法誤差在5%以內，精度較高，但在空氣域尺寸為4 m×4 m×13 m(高×寬×長)的幾何模型下，單元數目在千萬級以上，計算資源要求高，容易導致在大當量(800 kg)爆炸工況下由于計算時間過長、內存較大等原因而難以進行結構的爆炸動態響應分析。張社榮提出的網格劃分方法單元數目少，但是計算精度較低。對比文獻[2-4]提出的網格劃分方法，本文所提出的網格劃分方法在比例距離大于0.4 m/kg時，誤差在10%以內，雖然較文獻提出的方法精度低，適用的比例距離范圍較小，但該方法在同等的空氣域尺寸下，劃分的單元數目低于1百萬、計算內存占比低、計算效率較高，且在比例距離大于0.4 m/kg時精度在10%以內，符合工程誤差的標準，可以更好地應用于幾何模型較大的結構在大當量TNT柱狀炸藥工況下的數值計算。

表6 本文的網格劃分方法與文獻[2-4]方法的精度及效率對比

5 結論

本文使用LS-DYNA軟件，對大當量柱狀TNT炸藥在自由空氣中的沖擊波傳播過程進行了數值模擬，通過分析網格劃分比、網格加密范圍和爆炸近區漸進系數對不同當量下長徑比為3的TNT柱狀炸藥的計算精度及效率的影響，提出了一種“加密區網格劃分比不大于2/15，非加密區網格劃分比不大于1/4或10 cm(取兩者的較小值)，加密范圍不小于0.4 m/kg，加密區漸進比例系數為1.01”的網格劃分方法，研究結果表明：

1)院該網格劃分方法對不同當量長徑比為3的TNT柱狀炸藥均有較強的適用性，當比例距離大于0.4時，與Henrych經驗公式相比，超壓峰值誤差均在10%以內。

2) 該網格劃分方法單元數目較少、計算內存占比低、計算效率較高，當比例距離大于0.4 m/kg時，精度在10%以內，適用于幾何模型較大的結構在大當量TNT柱狀炸藥工況下的數值計算。