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

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

水下爆炸數(shù)值仿真

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

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

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

MSC.Dytran;水下爆炸;數(shù)值仿真;映射

0 引言

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

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

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

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

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

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

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

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

2 數(shù)值仿真

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

2.1 水的狀態(tài)方程

在絕熱條件下，水的狀態(tài)方程為［9］

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

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

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

2.2 TNT的狀態(tài)方程

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

式中:η=ρ/ρ0，ρ0為炸藥的初始密度，ρ為材料整體密度;e為單位質(zhì)量比內(nèi)能;A，B，ω，R1及R2為常數(shù)。對于TNT炸藥，密度為1 580 kg/m3，比內(nèi)能為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范圍內(nèi)的沖擊波衰減特性以及水下爆炸第一次氣泡脈動的最大半徑和周期。

1）一維球?qū)ΨQ計算模型

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

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

圖2 一維球?qū)ΨQ模型局部放大圖Fig.2 The local view of1－D spherical symmetrymodel

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

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

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

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

2）三維計算模型

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

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

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

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

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

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

3 一維球?qū)ΨQ模型數(shù)值模擬影響因素

3.1 網(wǎng)格劃分方式對計算結(jié)果的影響

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

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

3.2 網(wǎng)格劃分密度對計算結(jié)果的影響

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

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

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

3.3 計算區(qū)域大小對計算結(jié)果的影響

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

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

4 結(jié)語

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

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

2）同等網(wǎng)格密度及其他邊界條件下，一維球?qū)ΨQ模型與三維模型在水下爆炸數(shù)值仿真中所得的沖擊波壓力、速度、比內(nèi)能以及密度分布在藥包半徑范圍內(nèi)有一定的差異，但在藥包半徑范圍外二者的計算結(jié)果基本一致。

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

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

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

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

［2］李磊，馮順山.水下爆炸對艦船結(jié)構(gòu)毀傷效應(yīng)的研究現(xiàn)狀及展望［J］.艦船科學(xué)技術(shù)，2008，30(3）:26－30.

LILei，F(xiàn)ENG 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］張振華，朱錫，白雪飛.水下爆炸沖擊波的數(shù)值模擬研究［J］.爆炸與沖擊，2004，24(2）:182－188.

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

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

［6］方斌，朱錫.不同邊界條件下水下爆炸氣泡的數(shù)值模擬［J］.海軍工程大學(xué)學(xué)報，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］.羅耀杰，等譯，北京:國防工業(yè)出版社，1960.6－9.

［8］姚熊亮.艦船結(jié)構(gòu)振動沖擊與噪聲［M］.北京:國防工業(yè)出版社，2007.144－151.

［9］J.亨利奇.爆炸動力學(xué)及其應(yīng)用［M］.熊建國，等譯，北京:科學(xué)出版社，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－），男，博士研究生，工程師，從事水面艦船結(jié)構(gòu)抗爆抗沖擊設(shè)計工作。