徐維錚，吳衛(wèi)國

(1. 武漢理工大學(xué) 高性能艦船技術(shù)教育部重點實驗室，湖北 武漢 430063；2. 武漢理工大學(xué) 交通學(xué)院船舶、海洋與結(jié)構(gòu)工程系，湖北 武漢 430063)

0 引 言

當(dāng)爆炸發(fā)生在艦船艙室內(nèi)部時，由于沖擊波傳播過程中受到艙室壁面限制，將產(chǎn)生沖擊的多次反射和匯聚，以及持續(xù)時間較長的準(zhǔn)靜態(tài)壓力，對內(nèi)部結(jié)構(gòu)和人員產(chǎn)生更加嚴(yán)重的毀傷效應(yīng)。近年來，國內(nèi)外學(xué)者針對內(nèi)爆炸特性和計算方法方面做了大量的研究工作。

WU等[1]開展了封閉艙室內(nèi)爆炸實驗，研究了炸藥形狀和起爆位置對爆炸沖擊波載荷的影響規(guī)律。侯海量等[2]實驗研究了艙室內(nèi)爆炸沖擊波載荷特性?？紫樯氐萚3]實驗研究了藥量、角隅結(jié)構(gòu)形式對艙內(nèi)爆炸載荷影響規(guī)律。侯海量等[4]采用MSC.DYTRAN軟件數(shù)值研究了艙室內(nèi)爆炸沖擊波載荷特性?？紫樯氐萚5]采用MSC.DYTRAN軟件數(shù)值模擬了艦船艙室內(nèi)部戰(zhàn)斗部爆炸及爆炸毀傷效應(yīng)。樊壯卿等[6]基于LS-DYNA軟件，數(shù)值模擬了大體積復(fù)雜艙室內(nèi)爆炸載荷傳播特性。陳攀等[7]采用AUTODYN軟件數(shù)值研究了艙室內(nèi)爆工況下不同區(qū)域爆炸沖擊波載荷特性及爆點位置對艙室內(nèi)爆炸沖擊波特性的影響規(guī)律。

炸藥爆炸屬于高壓比、高密度比問題，其對激波捕捉格式提出較為嚴(yán)格的要求。傳統(tǒng)的商用軟件（MSC.DYTRAN，LS-DYNA，AUTODYN）主要采用低階精度的黎曼求解器和FCT歐拉算法，針對爆炸載荷的計算比較粗糙，且當(dāng)前發(fā)展比較成熟的激波捕捉算法還沒有嵌入到商用程序中，因此自主開發(fā)高精度爆炸載荷數(shù)值計算程序顯得更有意義。

目前，激波捕捉格式眾多，如TVD格式[8]、ENO格式[9]、NND格式[10]等，本文主要關(guān)注被廣泛采用的WENO格式（weighted essentially non-oscillation scheme），LIU等[11]于1994年提出了該格式，SHU等[12–14]發(fā)展了該格式。WENO在格式的有效性、通量的光滑性和收斂解的穩(wěn)定性方面均優(yōu)于ENO格式。

本文基于FORTRAN平臺，采用3階、5階、7階WENO有限差分格式，開發(fā)了艙室內(nèi)爆炸載荷高精度三維數(shù)值計算程序。采用所開發(fā)的程序開展了封閉艙室、泄壓艙室內(nèi)爆炸載荷數(shù)值計算，研究了WENO格式精度對艙室內(nèi)爆炸載荷影響規(guī)律，為結(jié)構(gòu)抗爆設(shè)計提供一定的參考和指導(dǎo)。

1 控制方程

假定炸藥瞬時爆轟，將凝聚相炸藥等效為高壓、高密度氣體，流場的控制方程采用三維可壓縮歐拉方程進行描述，

其中：

2 數(shù)值計算方法

程序中采用3階、5階、7階WENO有限差分格式對歐拉方程進行數(shù)值離散和求解，下面給出不同精度WENO格式的離散過程。

2.1 三階WENO格式

3階WENO格式（WENO-JS3）的數(shù)值離散和推導(dǎo)過程如下：

式中：

上式即為本文程序中采用的3階WENO格式的最終形式，對于含激波間斷流場，式中的根據(jù)下式求得：

2.2 五階WENO格式

五階WENO格式（WENO-JS5）的數(shù)值離散和推導(dǎo)過程如下：

式中：

式(16)即為本文程序中采用的5階WENO格式的最終形式，對于含激波間斷流場，式中的根據(jù)下式求得：

2.3 七階WENO格式

7階WENO格式（WENO-JS7）的數(shù)值離散和推導(dǎo)過程如下，單元面中心處數(shù)值通量的4種重構(gòu)方式分別為[16]：

式中：

式(25)即為本文程序中采用的7階WENO格式的最終形式，對于含激波間斷流場，式中的根據(jù)下式求得：

2.4 時間項數(shù)值離散

時間項采用3階TVD-RK法求解，格式如下[17]：

3 數(shù)值試驗

為了初步考察上述格式（WENO-JS3，WENO-JS5，WENO-JS7）的計算性能，選取3個經(jīng)典一維黎曼算例進行計算。

3.1 Sod激波管

該算例初始條件如式(30)所示[18]，網(wǎng)格數(shù)為400，計算結(jié)束時間為0.18。圖1給出計算結(jié)束時刻密度曲線圖及局部放大圖。

圖 1 計算結(jié)束后密度曲線及局部放大圖Fig. 1 Density curve and partial enlarged details

3.2 雙爆轟波碰撞

該算例初始條件如式(31)所示[19]，網(wǎng)格數(shù)為400，計算結(jié)束時間為0.038。圖2給出計算結(jié)束時刻密度曲線圖及局部放大圖。

3.3 激波與熵波相互作用

該算例初始條件如式(32)所示，網(wǎng)格數(shù)為400，計算結(jié)束時間為1.8。圖3給出計算結(jié)束時刻密度曲線圖及局部放大圖。

根據(jù)上述計算結(jié)果，評估各格式的計算性能發(fā)現(xiàn)：WENO-JS7格式計算性能最優(yōu)，WENO-JS5格式次之，WENO-JS3格式計算性能最低。即高精度WENO格式對含激波間斷流場具有更低的耗散特性和更高的分辨率。

圖 2 計算結(jié)束后密度曲線及局部放大圖Fig. 2 Density curve and partial enlarged details

圖 3 計算結(jié)束后密度曲線及局部放大圖Fig. 3 Density curve and partial enlarged details

4 艙室內(nèi)爆炸載荷計算

為了考察不同精度WENO格式對艙室內(nèi)爆載荷的影響規(guī)律，采用所開發(fā)的三維程序開展封閉艙室和泄壓艙室內(nèi)爆炸載荷數(shù)值計算。

4.1 密閉艙室內(nèi)爆炸

4.1.1 初始條件設(shè)置

密閉艙室尺寸為800×800×800 mm，壁面上設(shè)置2個測點，分別編號為NO1，NO2對爆炸超壓時間歷程進行輸出（見圖4）。

圖 4 艙室及測點分布Fig. 4 Cabin and gauging points

球形炸藥（200 g）位于艙室中間，等效后的高壓、高密度氣體參數(shù)：半徑為61.65 mm，密度為203.75 kg/m3，壓力為3.791×108Pa。計算初始條件見圖5（a）?？紤]計算時間及精度的要求，經(jīng)多次數(shù)值試驗，網(wǎng)格數(shù)為40×40×40，見圖5（b）。壁面邊界條件設(shè)置為反射邊界。

4.1.2 封閉艙室內(nèi)爆炸載荷特性

圖6給出封閉艙室內(nèi)爆炸壁面測點NO1和NO2超壓時間歷程曲線。從圖6可知，封閉艙室內(nèi)爆炸載荷主要包含多峰值、強度逐漸衰弱的沖擊波和持續(xù)時間較長、保持穩(wěn)定的準(zhǔn)靜態(tài)超壓。根據(jù)測點NO1和NO2的對比可知，測點的位置對爆炸前期爆炸波超壓峰值有較大的影響，而對最終形成的準(zhǔn)靜態(tài)超壓峰值幾乎沒有影響，說明封閉艙室內(nèi)部爆炸形成的準(zhǔn)靜態(tài)超壓峰值在空間上近似均勻。

4.1.3 不同精度爆炸載荷比較

為了探討不同精度WENO格式對封閉艙室內(nèi)爆炸載荷影響規(guī)律，選取壁面典型測點NO1對爆炸載荷進行輸出。

圖7（a）～圖7（c）給出不同精度WENO格式，測點NO1爆炸超壓時間歷程曲線，通過對比發(fā)現(xiàn)，高精度WENO格式能分辨出更精細(xì)的爆炸載荷特征，即高精度WENO格式具有更低的耗散特性和更高的分辨率。

圖 5 爆炸初場及網(wǎng)格分布Fig. 5 Initial condition and mesh distribution

圖 6 封閉艙室內(nèi)爆炸超壓時間歷程曲線Fig. 6 Overpressure time histories of all the gauging points in closed cabin

圖7（d）和圖7（e）分別給出3種精度WENO格式爆炸載荷對比圖及局部放大圖。從圖7（d）和圖7（e）可知，隨著WENO格式精度的提高，爆炸沖擊波載荷峰值得到較大的提高，即高精度WENO格式對于激波間斷具有更低的耗散特性。然而，不同精度WENO格式對最終形成的準(zhǔn)靜態(tài)超壓峰值影響較小，見圖7（d）中的粗實線。

4.2 泄壓艙室內(nèi)爆炸

4.2.1 初始條件設(shè)置

計算初始條件同4.1.1小節(jié)，唯一的不同是，泄壓艙室壁面上設(shè)置一個邊長為160 mm的方形泄壓口（見圖8）。壁面邊界條件設(shè)置為反射邊界，泄壓口處邊界條件設(shè)置為透射邊界。

圖 7 封閉艙室內(nèi)爆炸不同精度WENO格式測點NO1爆炸超壓時間歷程曲線Fig. 7 Overpressure time histories of gauging point NO1 for different schemes in closed cabin

圖 8 艙室及測點分布Fig. 8 Cabin and gauging points

圖 9 泄壓艙室內(nèi)爆炸超壓時間歷程曲線Fig. 9 Overpressure time histories of all the gauging points in venting cabin

4.2.2 泄壓艙室內(nèi)爆炸載荷特性

圖9給出泄壓艙室內(nèi)爆炸壁面上測點NO1和NO2超壓時間歷程曲線。從圖9可知，泄壓艙室內(nèi)爆炸載荷主要包含多峰值、強度逐漸衰弱的沖擊波和持續(xù)時間較長、呈現(xiàn)指數(shù)衰減規(guī)律[20]的準(zhǔn)靜態(tài)超壓。根據(jù)測點NO1和NO2的對比可知，測點的位置對爆炸前期爆炸波超壓峰值有一定的影響，而對準(zhǔn)靜態(tài)超壓時間歷程幾乎沒有影響，說明泄壓艙室內(nèi)爆炸形成的準(zhǔn)靜態(tài)超壓時間歷程在空間上近似均勻。

4.2.3 不同精度爆炸載荷比較

為了探討不同精度WENO格式對泄壓艙室內(nèi)爆炸載荷影響規(guī)律，選取壁面典型測點NO1對爆炸載荷進行輸出。

圖10（a）～圖10（c）給出不同精度WENO格式，測點NO1爆炸超壓時間歷程曲線，通過對比發(fā)現(xiàn)，高精度WENO格式能分辨出更精細(xì)的爆炸載荷特征，即高精度WENO格式具有更低的耗散特性和更高的分辨率。

圖 10 泄壓艙室內(nèi)爆炸不同精度WENO格式測點NO1爆炸超壓時間歷程曲線Fig. 10 Overpressure time histories of gauging point NO1 for different schemes in venting cabin

圖10（d）和圖10（e）分別給出3種精度WENO格式爆炸載荷對比圖及局部放大圖。從圖10（d）和圖10（e）可知，隨著WENO格式精度的提高，爆炸沖擊波載荷峰值得到較大的提高，即高精度WENO格式對于激波間斷具有更低的耗散特性。然而，不同精度WENO格式對逐漸形成的、呈現(xiàn)指數(shù)衰減規(guī)律的準(zhǔn)靜態(tài)超壓影響較小，見圖10（d）中的粗實線。

5 結(jié) 語

基于自主開發(fā)的高精度艙室內(nèi)爆炸三維數(shù)值計算程序研究了WENO格式精度對艙室內(nèi)爆炸載荷影響規(guī)律，通過本文的研究主要得到以下4點結(jié)論：

1）高精度WENO格式對含激波間斷流場具有更低的耗散特性和更高的分辨率。

2）WENO格式精度對艙室內(nèi)爆炸沖擊波載荷影響較大，高精度WENO格式給出更陡峭的激波峰值。

3）WENO格式精度對艙室內(nèi)爆炸準(zhǔn)靜態(tài)超壓（封閉艙室內(nèi)爆炸形成的保持不變的準(zhǔn)靜態(tài)超壓、泄壓艙室內(nèi)爆炸形成的呈指數(shù)衰減的準(zhǔn)靜態(tài)超壓）載荷影響較小。

4）本文所開發(fā)的三維程序可用于艦船艙室內(nèi)爆炸載荷數(shù)值計算，為結(jié)構(gòu)抗爆設(shè)計提供載荷依據(jù)。

[1]WU C, LUKASZEWICZ M, SCHEBELLA K, et al.Experimental and numerical investigation of confined explosion in a blast chamber [J]. Journal of Loss Prevention in the Process Industries, 2013, 26(4): 737–750.

[2]侯海量, 朱錫, 李偉, 等. 艙內(nèi)爆炸沖擊載荷特性實驗研究[J]. 船舶力學(xué), 2010, 14(8): 901–907.

HOU Hai-liang, ZHU Xi, LI Wei, et al. Experimental studies on characteristics of blast loading when exploded inside ship cabin [J]. Journal of Ship Mechanics, 2010, 14(8): 901–907.

[3]孔祥韶, 吳衛(wèi)國, 李俊, 等. 角隅結(jié)構(gòu)對艙內(nèi)爆炸載荷影響的實驗研究 [J]. 中國造船, 2012(3): 40–50.

KONG Xiang-shao, WU Wei-guo, LI Jun, et al. Experimental research of influence of corner structure on blast loading under inner explosion[J]. Chinese Journal of Ship Research, 2012(3):40–50.

[4]侯海量, 朱錫, 梅志遠(yuǎn). 艙內(nèi)爆炸載荷及艙室板架結(jié)構(gòu)的失效模式分析 [J]. 爆炸與沖擊, 2007, 27(2): 151–158.

HOU Hai-liang, ZHU Xi, MEI Zhi-yuan. Study on the blast load and failure mode of ship structure subject to internal explosion [J]. Explosion and Shock Waves, 2007, 27(2):151–158.

[5]孔祥韶, 吳衛(wèi)國, 李曉彬, 等. 艦船艙室內(nèi)部爆炸的數(shù)值模擬研究 [J]. 中國艦船研究, 2009, 4(4): 7–11.

KONG Xiang-shao, WU Wei-guo, LI Xiao-bin, et al.Numerical simualtion of cabin structure under inner explosion[J]. Chinese Journal of Ship Research, 2009, 4(4):7–11.

[6]樊壯卿, 王偉力, 黃雪峰, 等. 典型艙室內(nèi)爆炸仿真分析 [J].工程爆破, 2015, 21(3): 13–17.

FAN Zhuang-qing, WANG Wei-li, HUANG Xue-feng, et al.Simulation analysis on typical cabin internal explosion[J].Engineering Blasting, 2015, 21(3): 13–17.

[7]陳攀, 劉志忠. 艙室內(nèi)爆沖擊波載荷特性及影響因素分析[J]. 艦船科學(xué)技術(shù), 2016, 38(3): 43–48.

CHEN Pan, LIU Zhi-zhong. Research on loading of explosive and influencing factors inside closed cabin[J]. Ship Science and Technology, 2016, 38(3): 43–48.

[8]HARTEN A. High resolution schemes for hyperbolic conservation laws[J]. Journal of Computational Physics, 1983,135(3): 357–393.

[9]HARTEN A, ENGQUIST B, OSHER S, et al. Uniformly high order accurate essentially non-oscillatory schemes, III [J].Journal of Computational Physics, 1987, 71(2): 231–303.

[10]ZHANG H, ZHUANG F. NND schemes and their applications to numerical simulation of two-and three-dimensional flows*[M]//JOHN W H, THEODORE Y W. Advances in Applied Mechanics. Elsevier. 1991: 193–256.

[11]LIU X D, OSHER S, CHAN T. Weighted essentially nonoscillatory schemes [J]. Journal of Computational Physics,1994, 115(1): 200–212.

[12]JIANG G S, SHU C W. Efficient implementation of weighted ENO schemes [J]. Journal of Computational Physics, 1996,126(1): 202–202.

[13]SHU C W. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws [M].Springer Berlin Heidelberg. 1998: 325–432.

[14]SHU C W. High order weighted essentially nonoscillatory schemes for convection dominated problems [J]. Siam Review,2009, 51(1): 82–126.

[15]謝昌坦. 氣—水可壓縮流物質(zhì)界面的R-M不穩(wěn)定性研究[D]. 哈爾濱: 哈爾濱工程大學(xué), 2011.

XIE Chang-tan. The study on R-M instability of the material interface in gas-water compressible flow[D]. Harbin: Harbin Engineering University Master of Engineering, 2011.

[16]BALSARA D S, SHU C W. Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy [J]. Journal of Computational Physics, 2000,160(2): 405–452.

[17]SHU C W, OSHER S. Efficient implementation of essentially non-oscillatory shock-capturing schemes [J]. Journal of Computational Physics, 1988, 77(2): 439–471.

[18]SOD G A. A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws [J]. Journal of Computational Physics, 1978, 27(1): 1–31.

[19]WOODWARD P, COLELLA P. Colella, P. The numerical simulation of two-dimensional fluid flow with strong shocks[J]. Journal of Computational Physics, 1984, 54(1):115–173.

[20]ANDERSON C E, BAKER W E, WAUTERS D K, et al.Quasi-static pressure, duration, and impulse for explosions (e.g. HE) in structures [J]. International Journal of Mechanical Sciences, 1983, 25(6): 455–464.