水下爆炸數值仿真

吳國民，周心桃，肖漢林，段 宏

(中國艦船研究設計中心，湖北 武漢 430064）

水下爆炸數值仿真

吳國民，周心桃，肖漢林，段 宏

(中國艦船研究設計中心，湖北 武漢 430064）

基于MSC.Dytran軟件，通過一維球對稱數值仿真模型，模擬了水下爆炸沖擊波傳遞以及第一次氣泡脈動過程，其計算結果與經驗公式計算結果吻合較好，在此基礎上討論了網格劃分方式、網格密度以及計算區域大小對計算結果精度的影響。通過一維模型與三維模型數值計算結果的對比，論證了基于一維球對稱數值仿真模型，通過Remap映射技術進行三維模型數值仿真的計算方法。該方法在保證計算精度的同時較大地提高了計算效率，具有較強的工程應用價值。

MSC.Dytran;水下爆炸;數值仿真;映射

0 引言

艦艇在執行使命任務的過程中，會受到各種水下爆炸攻擊，有水下遠場爆炸、水下近場非接觸爆炸以及水下接觸爆炸等作用工況，因此，艦船抗水下爆炸性能是衡量其生命力的重要指標之一。在進行水下爆炸載荷作用下的艦艇抗爆抗沖擊設計時，水下爆炸數值仿真是重要的技術手段之一，通?？衫么笮偷纳逃糜邢拊浖M行建模計算，主要的軟件有AUTODYN、ABAQUS、LS-DYNA以及 MSC.Dytran，其中MSC.Dytran因為與功能強大的前后處理軟件MSC.Patran無縫連接，建模便捷，同時基于任意拉格朗日－歐拉求解方法(ALE）具有強大的流固耦合求解功能，因而被廣泛應用于水下爆炸數值仿真分析中［1－2］。國內外科研人員對于如何利用MSC.Dytran進行水下爆炸沖擊波傳遞以及氣泡脈動過程數值仿真也開展了深入的研究，如張振華等利用該軟件模擬了無限水域中球形藥包的水下爆炸沖擊波傳遞過程，并分析提出了通過修改水的狀態方程參數提升沖擊波壓力，以有效降低計算單元數量和計算時間的方法［3］;鄧文斌等也對水下爆炸沖擊波進行了數值仿真，并提出了優化計算文件節約計算時間提高解題效率的方法［4］;Matsumoto和方斌等則分別利用該軟件實現了對水下爆炸氣泡脈動的數值模擬，并重點分析了邊界條件對氣泡脈動壓力、脈動周期等特性的影響［5－6］。本文利用 MSC.Dytran軟件，通過一維球對稱模型對水下爆炸的沖擊波傳遞以及第一次氣泡脈動過程進行數值仿真，在此基礎上分析提出基于Remap映射技術提高三維模型數值仿真計算效率的計算方法，可為水下遠場非接觸爆炸作用下的艦船沖擊動響應分析提供借鑒。

1 水下爆炸沖擊波傳遞及氣泡脈動過程

在無限水域中爆炸沖擊波傳遞的主要過程可分為5個階段，如圖1所示，分別是指數衰減階段、倒數衰減階段、倒數衰減后端、第一次氣泡膨脹收縮段和脈動壓力段。在爆炸起始階段，沖擊波壓力在瞬間達到峰值，然后迅速衰減，同時爆炸產生的高壓氣體形成氣泡，開始第一次氣泡脈動，先是因內部壓力大而膨脹，達到氣泡半徑最大值后，在外部壓力的作用下，氣泡開始收縮，當氣泡半徑收縮到最小時，出現脈動壓力峰值，此壓力峰值不大于沖擊波壓力峰值的10%～20%，但脈動壓力持續的時間卻大大超過沖擊波壓力作用時間，因此，二者的沖量值大小是相當的［7］。在此之后氣泡脈動能量已所剩無幾，因此對于工程問題而言，主要考慮爆炸開始到第一次氣泡脈動結束這一過程［8］。

圖1 水下爆炸沖擊波壓力時歷曲線Fig.1 The time-pressure curve of underwater explosion

對于水下爆炸的前5個階段，庫爾在早期提出的經驗公式給出了很好的描述，并經試驗驗證具有良好的精度，因此一直沿用至今。其中對于沖擊波壓力峰值及指數衰減階段用下列公式表示［7］:

式中:Pm為沖擊波壓力峰值，Pa;W為TNT球狀藥包質量，kg;R為爆心到測點的距離，m;R0為藥包的初始半徑，m;θ為沖擊波時間衰減常數，s;而第一次氣泡的最大半徑rmax和脈動周期T分別用下列公式進行計算:

式中z為爆心處的流體靜壓力水柱值，m。

下節水下爆炸數值仿真將針對上述沖擊波傳遞以及第一次氣泡脈動壓力過程進行模擬，并與經驗公式的計算結果進行對比，以驗證數值仿真的精度。

2 數值仿真

對于水下爆炸數值仿真，MSC.Dytran采用歐拉方法進行求解，在時間域上采用顯示積分法，在空間上用歐拉單元離散，單元均為一階單元，可采用8節點任意六面體單元、6節點三棱柱單元和4節點四面體單元。對于歐拉網格中的炸藥和流體，分別利用狀態方程來定義其壓力與密度及比內能之間的函數關系。

2.1 水的狀態方程

在絕熱條件下，水的狀態方程為［9］

式中:p為壓力;ρ0為水在常溫狀態下的密度;ρ為整體密度;B和n為常數，B=3 045 kg/cm2，n=7.15。

在MSC.Dytran中，水的狀態方程用多項式狀態方程(EOSPOL）來模擬，其中壓力是相對體積及比內能的多項式函數［10］:

式中:e為單位質量比內能;μ=(ρ/ρ0）－1。計算中，取p=a1μ，體積模量 a1=2.2 GPa，并取 ρ0=1 000 kg/m3，e=83.950 kJ/kg。

2.2 TNT的狀態方程

在MSC.Dytran中，用JWL狀態方程(EOSJWL）來模擬炸藥爆炸過程:

式中:η=ρ/ρ0，ρ0為炸藥的初始密度，ρ為材料整體密度;e為單位質量比內能;A，B，ω，R1及R2為常數。對于TNT炸藥，密度為1 580 kg/m3，比內能為4.19MJ/kg，初始爆轟波速度為6 930m/s，A=371.2 GPa，B=3.231 GPa，ω =0.3，R1=4.15，R2=0.95。

2.3 計算工況

假定在靜止的水域中，水深40 m處，1 kg的TNT球形炸藥(藥包半徑為0.053 3 m）在水中爆炸，水面的大氣壓為1.01×105Pa，計算分析在離藥包中心0.4～1.0 m范圍內的沖擊波衰減特性以及水下爆炸第一次氣泡脈動的最大半徑和周期。

1）一維球對稱計算模型

在上述計算工況中，水深相對于藥包的直徑為大值，且在爆源附近的爆炸沖擊波作用過程為毫秒級的強間斷，可以近似地認為爆源附近的沖擊波傳遞為1個球對稱過程。因此，可用一維球對稱模型在MSC.Dytran中進行計算分析。

建立一維球對稱模型如圖2所示，坐標原點為爆心，沿x方向建模，楔形夾角為2°，水域長為30m，從爆心到水域邊界采用bias網格劃分方法，由密到疏共劃分6 000個歐拉單元，在x方向上最大單元與最小單元尺寸之比為3，在x=30 m的水域邊界處設置流入流出邊界條件。計算時間為0.10 s，對于一維球對稱模型在考慮計算時間步長時，只考慮徑向單元的尺寸。在x=0.4～1.0 m范圍內均勻設置7個考核點，以記錄沖擊波傳播過程和衰減規律。

圖2 一維球對稱模型局部放大圖Fig.2 The local view of1－D spherical symmetrymodel

在x=1.0 m處考核點記錄的沖擊波壓力時歷曲線如圖3所示。該曲線反映了沖擊波的衰減和第一次氣泡脈動的典型過程，圖4給出了x=1.0 m處沖擊波壓力衰減曲線與經驗公式計算結果的對比，圖5則給出了x=0.4～1.0 m共7個考核點的沖擊波壓力峰值與經驗公式結果的對比，從圖上可以看出二者很吻合。

圖6給出的是由第一次氣泡體積轉化得到的氣泡半徑的變化時歷曲線。從上述結果可知，當t=0.035 8 s時，第一次氣泡脈動的氣泡半徑達到最大值rmax=0.864 m，而第一次氣泡脈動周期為T=0.069 5 s。上述計算結果與經驗公式的比較見表1，誤差在10%以內。

圖6 氣泡半徑的時歷曲線Fig.6 The time-radius curve of the bubble's impulse

表1 第一次氣泡脈動特性對比Tab.1 The comparison of the first bubble's impulse

2）三維計算模型

一維球對稱模型由于單元少因此計算效率很高，同時上述分析表明其計算精度也在工程誤差范圍內，但應用到工程問題一般都需要建立三維模型進行計算分析，MSC.Dytran提供了從一維模型計算結果到三維模型的Remap映射技術。下面將首先建立三維模型進行計算，然后將一維模型與三維模型的計算結果進行對比，再利用Remap映射技術，將映射后的三維模型計算結果與直接用三維模型計算的結果進行對比，以論證一維球對稱模型在工程問題中應用的可行性。

為了使三維模型網格密度與一維球對稱模型一致，選取建立0.502 7 m×0.502 7 m×0.502 7 m的立方體水域三維模型，坐標原點為爆心，同樣是球形TNT裝藥1 kg，分別沿x，y及z向用bias網格劃分方法劃分181個六面體單元，共計5 929 741個單元，最大單元尺寸與最小單元尺寸之比為1.033 3，邊界條件及狀態方程與一維球對稱模型一致。

圖7給出了t=0.16 ms時刻一維模型與三維模型的計算結果對比，選取的是壓力、速度、比內能以及密度在x方向上的分布曲線對比，從圖中可以看出，在離爆心較近的區域尤其是1倍藥包半徑范圍內，2種模型中的計算結果有一定差異，這是由于爆轟波的初始模擬對網格較敏感，爆心越近，單元網格差異越大，尤其在爆心處，一維模型為4個節點共用的五面體錐形單元，而三維模型為規整的六面長方體單元，單元的差異導致離爆心較近區域的計算結果差異較大，而在1倍藥包半徑范圍外，2種模型的計算結果吻合得很好。

3）一維模型向三維模型的映射

所謂一維模型向三維模型的映射，是指先利用一維球對稱模型計算爆炸起始時刻t0到某一時刻t1，選取的t1必須是流體介質尚未完全流出計算區域的時刻，提取t1時刻一維球對稱模型計算結果中沿徑向的密度、比內能以及速度分布曲線(如圖7所示），并以數據對的形式分別存儲成*.xyd文件，然后在三維模型的*.dat文件中手工添加TABFILE語句調用上述*.xyd文件，計算時將上述密度、比內能以及速度分布映射到三維模型中，作為三維模型計算的第1個載荷步，即若三維模型需要求解的總的時間段為t0～t2時刻，前面的t0～t1時刻用一維模型計算，將t1時刻的計算結果導入三維模型，利用三維模型接著進行t1～t2時刻的計算分析，用上述接力方式達到大大節約計算時間、提高求解效率的目的［10］。

圖8給出的是t=0.16 ms時刻一維模型壓力分布云圖，其沖擊波陣面到達x=0.131 m處，將上述時刻的計算結果映射到三維模型中，如圖9所示，三維模型的第1個載荷步的壓力分布云圖與一維模型的完全一致，在此基礎上三維映射模型計算到t=0.02 ms時刻的壓力分布如圖10所示，此時相對于起爆點的絕對時刻為t=0.16+0.02=0.18 ms;直接用三維模型進行求解到t=0.18 ms時刻的壓力分布如圖11所示。由圖10與圖11進行對比，二者的壓力分布規律和最大值以及沖擊波陣面到達位置均非常一致。而在同樣的軟硬件條件下，上述直接用三維模型求解所耗的機時是利用映射技術求解所耗機時的10倍以上。因此，基于一維模型計算結果，利用映射技術，再建立三維映射模型進行求解，大大提高了計算效率，同時計算精度有所保證。

3 一維球對稱模型數值模擬影響因素

3.1 網格劃分方式對計算結果的影響

在建立一維球對稱模型時，網格劃分方式有非均勻劃分和均勻劃分2類，非均勻劃分即如上節中用bias網格劃分方法進行單元劃分，靠近爆心處的網格相對較密;均勻劃分即所有的單元在x方向的尺寸一致。如圖2所示的模型，在30 m區域內分別用非均勻和均勻2種方式劃分6 000個單元，在其他條件相同的情況下，在x=0.4～1.6 m的范圍內，二者的壓力計算結果對比如圖12所示。非均勻劃分網格方式的計算結果與經驗公式計算結果更接近，二者計算所得第一次氣泡脈動最大半徑和周期差別不大，說明用非均勻網格劃分方式的計算精度更高。

圖12 不同網格劃分方式壓力計算結果對比Fig.12 The comparison of the x－Pressure curves calculated with differentmethods ofmeshing

3.2 網格劃分密度對計算結果的影響

如圖2所示的一維球對稱模型，在30 m區域內采用同樣的網格劃分方式，即在x方向上最大單元與最小單元尺寸之比為3，分別劃分2 000，4 000，6 000，8 000以及10 000個單元，在其他條件相同的情況下，在x=0.4～1.6 m的范圍內，壓力計算結果對比如圖13所示，可以看出網格密度在6 000和8 000的壓力計算結果，比其他網格密度的結果更為接近經驗公式結果;氣泡脈動特性對比如表2所示，計算表明網格越密，氣泡脈動周期和半徑就越小。

圖13 不同網格密度壓力計算結果對比Fig.13 The comparison of the x－pressure curves calculated with different densities ofmesh

表2 不同網格密度的第一次氣泡脈動特性對比Tab.2 The comparison of the characters of the first bubble's impulse with different densities ofmesh

3.3 計算區域大小對計算結果的影響

分別建立10，20，30和40m區域的一維球對稱計算模型，每個模型在x=10 m以內的范圍按bias網格劃分方式劃分2 000個單元，即在x方向上最大單元與最小單元尺寸之比為3，x=10 m以外的范圍按同樣的單元尺寸進行均勻劃分。計算結果表明，在x=0.4～1.6 m的范圍內，壓力計算結果完全一致，不受計算區域大小的影響，但是氣泡脈動特性受計算區域影響較大，如表3所示，計算區域越大，氣泡半徑和脈動周期就越大，越接近經驗公式計算值，這說明雖然計算模型中在邊界處定義了流入流出的邊界條件，但是計算區域過小還是會有邊界反射效應，影響氣泡特性的計算精度。

表3 不同計算區域的第一次氣泡脈動特性對比Tab.3 The comparison of the characters of the first bubble's impulse with different ranges of the numericalmodel

4 結語

通過上面的計算分析，可以得到以下結論:

1）一維球對稱模型可以較好地模擬無限水域條件下水下爆炸沖擊波傳遞以及第一次氣泡脈動過程，數值仿真得出的沖擊波壓力、第一次氣泡脈動半徑和周期均與經驗公式能較好地吻合。

2）同等網格密度及其他邊界條件下，一維球對稱模型與三維模型在水下爆炸數值仿真中所得的沖擊波壓力、速度、比內能以及密度分布在藥包半徑范圍內有一定的差異，但在藥包半徑范圍外二者的計算結果基本一致。

3）基于一維模型計算結果，利用Remap映射技術進行三維模型數值仿真的計算方法，大大提高了求解效率，同時計算精度有所保證。上述計算方法在工程應用上很有意義，如對于圖14所示的水下遠場非接觸爆炸作用下艦船沖擊響應分析問題，在沖擊波陣面傳遞到離船體一定距離產生邊界反射效應前，其沖擊波傳遞過程可用無限水域下的水下爆炸沖擊波傳遞規律模擬，可利用上述計算方法，先用一維模型模擬上述沖擊波陣面的傳遞過程，然后利用映射技術進行后續的三維模型數值仿真計算，不但解決了因遠場水下爆炸仿真模型水單元數量巨大而無法計算的問題，而且計算效率大大提高。

圖14 沖擊波加載示意圖Fig.14 The illustration of loading shock wave

4）利用一維模型進行數值計算時，采用非均勻的bias網格劃分方法進行單元劃分，靠近爆心處的網格相對較密，計算精度更高;網格密度對沖擊波壓力和氣泡脈動特性計算值有一定影響，網格單元尺寸要適中，單元尺寸過小或過大，都會使在被考核的距離范圍內沖擊波壓力值偏離經驗公式的計算結果，而第一次氣泡脈動半徑和周期則是隨著網格密度的增大而減小，因此需要針對關注的具體問題選擇合適的網格密度;而在同等網格密度下，計算區域大小對沖擊波壓力計算值幾乎沒有影響，但對氣泡脈動特性有影響，計算區域越大，邊界反射效應影響越小，氣泡半徑和脈動周期計算值越準確。

［1］尹群，陳永念，等.水下爆炸研究的現狀和趨勢［J］.造船技術，2003，(6）:6－12.

［2］李磊，馮順山.水下爆炸對艦船結構毀傷效應的研究現狀及展望［J］.艦船科學技術，2008，30(3）:26－30.

LILei，FENG Shun-shan.Present state and perspectives of damage effect to ship construction subjected to UNDEX［J］.Ship Science and Technology，2008，30(3）:26 －30.

［3］張振華，朱錫，白雪飛.水下爆炸沖擊波的數值模擬研究［J］.爆炸與沖擊，2004，24(2）:182－188.

［4］鄧文彬，王國治.基于DYTRAN軟件的水下爆炸數值計算［J］.華東船舶工業學院學報(自然科學版），2003，17(6）:11－16.

［5］MATSUMOTO.Boundary curvature effects on gas bubble oscillations in underwater explosion［R］.ADA308087，1996.

［6］方斌，朱錫.不同邊界條件下水下爆炸氣泡的數值模擬［J］.海軍工程大學學報，2008，20(2）:85－90.

FANG Bin，ZHU Xi.Numerical simulation of underwater explosion bubble with different boundaries［J］.Journal of Naval University of Engineering，2008，20(2）:85 －90.

［7］庫爾.水下爆炸［M］.羅耀杰，等譯，北京:國防工業出版社，1960.6－9.

［8］姚熊亮.艦船結構振動沖擊與噪聲［M］.北京:國防工業出版社，2007.144－151.

［9］J.亨利奇.爆炸動力學及其應用［M］.熊建國，等譯，北京:科學出版社，1987.59 －62.

［10］MSC.Dytran User's Manual［M］.MSC Corporation，2008.

Numerical simulation of underwater explosion

WU Guo-min，ZHOU Xin-tao，XIAO Han-lin，DUAN Hong
(China Ship Development and Design Center，Wuhan 430064，China）

Using the software MSC Dytran，the shock wave's spread and the first bubble's impulse in underwater explosion were simulated by 1-D spherical symmetry numericalmodel.The results of numerical simulation were agreeable with the ones of experimental formulations.The factors influencing on the calculation resultswere compared and analyzed，such as the method ofmeshing，density ofmesh and the range of the numericalmodel.Through comparing the calculation results of 1-D model with the ones of 3-D model，the calculation method of translating 1-D spherical symmetry model into 3-D model by Remap was proved out.The above method improved the efficiency of numerical simulation，and also ensured the necessary precision，so itwas significant to the actual engineering problems.

MSC.Dytran;underwater explosion;numerical simulation;Remap

E925.2

A

1672－7649(2012）09－0020－07

10.3404/j.issn.1672－7649.2012.09.004

2012－03－19;

2012－04－10

吳國民(1980－），男，博士研究生，工程師，從事水面艦船結構抗爆抗沖擊設計工作。