模拟舱室内部气云爆炸载荷的不同精度WENO格式比较

徐维铮; 吴卫国

doi:10.11858/gywlxb.20170689

模拟舱室内部气云爆炸载荷的不同精度WENO格式比较

doi: 10.11858/gywlxb.20170689

徐维铮^{1, 2,},
吴卫国^{1, 2}

1.
武汉理工大学高性能舰船技术教育部重点实验室, 湖北武汉 430063
2.
武汉理工大学交通学院船舶、海洋与结构工程系, 湖北武汉 430063

基金项目:

装备预研教育部联合基金（青年人才） 6141A020331

国家自然科学基金 51409202

中央高校基本科研业务费 2016-YB-016

详细信息

作者简介:
徐维铮(1991-), 男, 博士研究生, 主要从事束空间内炸药爆炸场高精度数值计算方法及三维程序开发研究.E-mail:xuweizheng@whut.edu.cn

中图分类号: O381
计量
- 文章访问数: 7129
- HTML全文浏览量: 3099
- PDF下载量: 200
出版历程
- 收稿日期: 2017-12-07
- 修回日期: 2017-12-19

Comparisons of Different Precision WENO Schemes for Simulating Blast Load of Gas Cloud Explosion inside a Cabin

XU Weizheng^{1, 2
,},
WU Weiguo^{1, 2}

1.
Key Laboratory of High Performance Ship Technology of Ministry of Education, Wuhan University of Technology, Wuhan 430063, China
2.
Departments of Naval Architecture, Ocean and Structural Engineering, School of Transportation, Wuhan University of Technology, Wuhan 430063, China

摘要

摘要: 为了研究加权本质无振荡（WENO）格式精度对舱室内部气云爆炸载荷的影响规律，基于FORTRAN平台，采用三、五、七、九阶WENO格式，开发了高精度舱室内部气云爆炸三维数值计算程序。选用Sod激波管、激波与熵波相互作用两个经典算例验证了程序编写的可靠性，并初步考察了不同精度WENO格式的计算性能。采用已验证的程序开展了封闭舱室和泄压舱室内部球体气云爆炸的数值模拟，并探讨了WENO格式精度对爆炸载荷的影响规律。研究表明：舱室内部气云爆炸载荷主要包含瞬态冲击波和持续时间较长的准静态超压；WENO格式精度对冲击波载荷影响较大，高阶格式给出更陡峭的峰值，而对形成的准静态超压影响较小。
- WENO格式 /
- 高精度 /
- 低耗散 /
- 舱室内部气云 /
- 爆炸载荷
Abstract: In the present work, we developed a high-resolution 3D code based on the FORTRAN platform to investigate the influence of accuracy of WENO schemes on the blast load generated by gas cloud explosion inside a cabin.The 3th, 5th, 7th and 9th WENO finite difference schemes were implemented in the code to capture the shock waves.Several one dimensional Riemann problems such as Sod shock tube and shock-entropy wave interaction were simulated to investigate the computing performance of the schemes preliminarily and validate the developed code.Then, the validated code was used to conduct the simulation of the gas cloud explosion inside a closed cabin and a venting cabin and the influence of the accuracy of WENO schemes on the blast load was discussed as well.The researches indicate that the blast load generated by gas cloud explosion inside the cabin mainly contains multiple peaks shock waves and long duration quasi-static overpressure.The accuracy of WENO schemes has considerable influence on the blast waves and higher order scheme gives sharper pressure peaks, whereas it has much less influence on the quasi-static overpressure.
- WENO scheme /
- high precision /
- low dissipation /
- gas cloud inside the cabin /
- blast load

HTML全文

可燃性气云爆炸是石化工业灾害预防等领域研究的焦点。开敞空间中气云爆炸后将形成强度较大的带有负压区的空气冲击波，对工作人员和周边结构设施造成较大的损害。近年来学者们针对可燃气云、云雾爆炸进行了相应的数值、实验研究^[1-4]。这些研究中可燃气云、云雾被等效地简化为一团均匀高压气体，且主要采用低阶精度数值方法进行模拟；而针对舱室内部可燃气云爆炸的研究较少。

高精度激波捕捉格式对含激波流场的数值模拟具有重要意义，不但可以降低网格的规模，而且能较好地分辨流场中复杂的波系结构。Liu等^[5]于1994年首次提出加权本质无振荡(Weighted Essentially Non-Oscillatory，WENO)格式。之后，Jiang和Shu^[6]通过引入新的光滑因子，提出三阶WENO格式和五阶WENO格式，并扩展了其应用^[7-8]。Balsara等^[9]将WENO格式推广到更高阶形式，给出七、九、十一阶WENO格式。本研究在参考上述文献的基础上，采用三、五、七、九阶WENO格式，基于FORTRAN平台自主开发了高精度舱室内部气云爆炸三维数值计算程序，研究WENO格式精度对舱室内部气云爆炸载荷的影响规律。

1. 控制方程

可燃气云简化为均匀高压气体，采用三维可压缩欧拉方程描述爆炸流场，其具体形式如下

$\frac{{\partial \mathit{\boldsymbol{U}}}}{{\partial t}} + \frac{{\partial \mathit{\boldsymbol{E}}}}{{\partial x}} + \frac{{\partial \mathit{\boldsymbol{F}}}}{{\partial y}} + \frac{{\partial \mathit{\boldsymbol{G}}}}{{\partial z}} = 0$

(1)

其中

$\mathit{\boldsymbol{U}} = \left[ \begin{array}{l} \rho \\ \rho u\\ \rho v\\ \rho w\\ E \end{array} \right],\;\;\;\;\mathit{\boldsymbol{E}} = \left[ \begin{array}{l} \rho u\\ \rho {u^2} + p\\ \rho uv\\ \rho uw\\ u\left( {E + p} \right) \end{array} \right],\;\;\;\;\mathit{\boldsymbol{F}} = \left[ \begin{array}{l} \rho v\\ \rho vu\\ \rho {v^2} + p\\ \rho vw\\ v\left( {E + p} \right) \end{array} \right], \\ \mathit{\boldsymbol{G}} = \left[ \begin{array}{l} \rho w\\ \rho wu\\ \rho wv\\ \rho {w^2} + p\\ w\left( {E + p} \right) \end{array} \right]$

(2)

$E = \rho e + \frac{1}{2}\rho {u^2} + \frac{1}{2}\rho {v^2} + \frac{1}{2}\rho {w^2}$

(3)

$p = \left( {\gamma - 1} \right)\rho e$

(4)

式中: ρ是密度；u、v、w是x、y、z方向上的速度分量；p为流体压力；E是单位体积流体的总能量；e是比内能；γ表示气体绝热指数，取为1.4。对于多维欧拉方程的计算，采用Strang维数分裂法，将欧拉方程分解为x、y、z 3个方向求解。

2. 数值方法

程序中采用三阶、五阶、七阶、九阶WENO有限差分格式对欧拉方程进行数值离散和求解，下面给出不同精度WENO格式的离散过程。

2.1 三阶WENO格式

三阶WENO格式(WENO-JS3)的数值离散和推导过程如下。单元面中心x_i+1/2处数值通量f_i+1/2的两种重构方式分别为^[6]

$\left\{ \begin{array}{l} f_{i + 1/2}^0 = - \frac{1}{2}{f_{i - 1}} + \frac{3}{2}{f_i}\\ f_{i + 1/2}^1 = \frac{1}{2}{f_i} + \frac{1}{2}{f_{i + 1}} \end{array} \right.$

(5)

利用上述两种模板的凸组合计算数值通量f_i+1/2，即

${f_{i + 1/2}} = {\omega _0}f_{i + 1/2}^0 + {\omega _1}f_{i + 1/2}^1$

(6)

对于含激波间断流场，(6)式中的ω_k根据下式求得

${\omega _k} = {\alpha _k}/\sum\limits_{s = 0}^1 {{\alpha _s}} ,{\alpha _k} = \frac{{{d_k}}}{{{{\left( {\varepsilon + {\beta _k}} \right)}^2}}},\;\;\;\;k = 0,1$

(7)

式中：d_k表示WENO格式的线性权值，对于三阶WENO格式，其值为d₀=1/3，d₁=2/3。为避免分母为零，取ε=1.0×10^-6。光滑因子β_k(k=0, 1)的表达式如下^[6]

$\left\{ \begin{array}{l} {\beta _0} = {\left( {{f_i} - {f_{i - 1}}} \right)^2}\\ {\beta _1} = {\left( {{f_{i + 1}} - {f_i}} \right)^2} \end{array} \right.$

(8)

2.2 五阶WENO格式

五阶WENO格式(WENO-JS5)的数值离散和推导过程如下。单元面中心x_i+1/2处数值通量f_i+1/2的3种重构方式分别为^[6]

$\left\{ \begin{array}{l} f_{i + 1/2}^0 = \frac{1}{3}{f_{i - 2}} - \frac{7}{6}{f_{i - 1}} + \frac{{11}}{6}{f_i}\\ f_{i + 1/2}^1 = - \frac{1}{6}{f_{i - 1}} + \frac{5}{6}{f_i} + \frac{1}{3}{f_{i + 1}}\\ f_{i + 1/2}^2 = \frac{1}{3}{f_i} + \frac{5}{6}{f_{i + 1}} - \frac{1}{6}{f_{i + 2}} \end{array} \right.$

(9)

利用上述3种模板的凸组合计算数值通量f_i+1/2，即

${f_{i + 1/2}} = {\omega _0}f_{i + 1/2}^0 + {\omega _1}f_{i + 1/2}^1 + {\omega _2}f_{i + 1/2}^1$

(10)

对于含激波间断流场，(10)式中的ω_k根据下式求得

${\omega _k} = {\alpha _k}/\sum\limits_{s = 0}^2 {{\alpha _s}} ,{\alpha _k} = \frac{{{d_k}}}{{{{\left( {\varepsilon + {\beta _k}} \right)}^2}}},\;\;\;\;k = 0,1,2$

(11)

对于五阶WENO格式，d₀= $\frac{1}{{10}}$ ，d₁= $\frac{3}{{5}}$ ，d₂= $\frac{3}{{10}}$ 。光滑因子β_k(k=0, 1, 2)的表达式如下

$\left\{ \begin{array}{l} {\beta _0} = \frac{{13}}{{12}}{\left( {{f_{i - 2}} - 2{f_{i - 1}} + {f_i}} \right)^2} + \frac{1}{4}{\left( {{f_{i - 2}} - 4{f_{i - 1}} + 3{f_i}} \right)^2}\\ {\beta _1} = \frac{{13}}{{12}}{\left( {{f_{i - 1}} - 2{f_i} + {f_{i + 1}}} \right)^2} + \frac{1}{4}{\left( {{f_{i - 1}} - {f_{i + 1}}} \right)^2}\\ {\beta _2} = \frac{{13}}{{12}}{\left( {{f_i} - 2{f_{i + 1}} + {f_{i + 2}}} \right)^2} + \frac{1}{4}{\left( {3{f_i} - 4{f_{i + 1}} + {f_{i + 2}}} \right)^2} \end{array} \right.$

(12)

2.3 七阶WENO格式

七阶WENO格式(WENO-JS7)的数值离散和推导过程如下。单元面中心x_i+1/2处数值通量f_i+1/2的4种重构方式分别为^[9]

$\left\{ \begin{array}{l} f_{i + 1/2}^0 = - \frac{1}{4}{f_{i - 3}} + \frac{{13}}{{12}}{f_{i - 2}} - \frac{{23}}{{12}}{f_{i - 1}} + \frac{{25}}{{12}}{f_i}\\ f_{i + 1/2}^1 = \frac{1}{{12}}{f_{i - 2}} - \frac{5}{{12}}{f_{i - 1}} + \frac{{13}}{{12}}{f_i} + \frac{1}{4}{f_{i + 1}}\\ f_{i + 1/2}^2 = - \frac{1}{{12}}{f_{i - 1}} + \frac{7}{{12}}{f_i} + \frac{7}{{12}}{f_{i + 1}} - \frac{1}{{12}}{f_{i + 2}}\\ f_{i + 1/2}^3 = \frac{1}{4}{f_i} + \frac{{13}}{{12}}{f_{i + 1}} - \frac{5}{{12}}{f_{i + 2}} + \frac{1}{{12}}{f_{i + 3}} \end{array} \right.$

(13)

利用上述4种模板的凸组合计算数值通量f_i+1/2，即

${f_{i + 1/2}} = {\omega _0}f_{i + 1/2}^0 + {\omega _1}f_{i + 1/2}^1 + {\omega _2}f_{i + 1/2}^2 + {\omega _3}f_{i + 1/2}^3$

(14)

对于含激波间断流场，(14)式中的 ω_k 根据下式求得

${\omega _k} = {\alpha _k}/\sum\limits_{s = 0}^3 {{\alpha _s}} ,{\alpha _k} = \frac{{{d_k}}}{{{{\left( {\varepsilon + {\beta _k}} \right)}^2}}},\;\;\;\;k = 0,1,2,3$

(15)

对于七阶WENO格式，d₀= $\frac{1}{{35}}$ ，d₁= $\frac{12}{{35}}$ ，d₂= $\frac{18}{{35}}$ ，d₃= $\frac{4}{{35}}$ 。光滑因子β_k(k=0, 1, 2, 3)的表达式如下^[9]

$\left\{ \begin{array}{l} {\beta _0} = {f_{i - 3}}\left( {547{f_{i - 3}} - 3882{f_{i - 2}} + 4642{f_{i - 1}} - 1854{f_i}} \right) + \\ \;\;\;\;\;\;\;{f_{i - 2}}\left( {7043{f_{i - 2}} - 17246{f_{i - 1}} + 7042{f_i}} \right) + \\ {f_{i - 1}}\left( {11003{f_{i - 1}} - 9402{f_i}} \right) + 2107f_i^2\\ {\beta _1} = {f_{i - 2}}\left( {267{f_{i - 2}} - 1642{f_{i - 1}} + 1602{f_i} - 494{f_{i + 1}}} \right) + \\ \;\;\;\;\;\;\;{f_{i - 1}}\left( {2843{f_{i - 1}} - 5966{f_i} + 1922{f_{i + 1}}} \right) + \\ {f_i}\left( {3443{f_i} - 2522{f_{i + 1}}} \right) + 547f_{i + 1}^2\\ {\beta _2} = {f_{i - 1}}\left( {547{f_{i - 1}} - 2522{f_i} + 1922{f_{i + 1}} - 494{f_{i + 2}}} \right) + \\ \;\;\;\;\;\;\;{f_i}\left( {3443{f_i} - 5966{f_{i + 1}} + 1602{f_{i + 2}}} \right) + \\ {f_{i + 1}}\left( {2843{f_{i + 1}} - 1642{f_{i + 2}}} \right) + 267f_{i + 2}^2\\ {\beta _3} = {f_i}\left( {2107{f_i} - 9402{f_{i + 1}} + 7042{f_{i + 2}} - 1854{f_{i + 3}}} \right) + \\ \;\;\;\;\;\;{f_{i + 1}}\left( {11003{f_{i + 1}} - 17246{f_{i + 2}} + 4642{f_{i + 3}}} \right) + \\ {f_{i + 2}}\left( {7043{f_{i + 2}} - 3882{f_{i + 3}}} \right) + 547f_{i + 3}^2 \end{array} \right.$

(16)

2.4 九阶WENO格式

九阶WENO格式(WENO-JS9)的数值离散和推导过程如下。单元面中心x_i+1/2处数值通量f_i+1/2的5种重构方式分别为^[9]

$\left\{ \begin{array}{l} f_{i + 1/2}^0 = \frac{1}{5}{f_{i - 4}} - \frac{{21}}{{20}}{f_{i - 3}} + \frac{{137}}{{60}}{f_{i - 2}} - \frac{{163}}{{60}}{f_{i - 1}} + \frac{{137}}{{60}}{f_i}\\ f_{i + 1/2}^1 = - \frac{1}{{20}}{f_{i - 3}} + \frac{{17}}{{60}}{f_{i - 2}} - \frac{{43}}{{60}}{f_{i - 1}} + \frac{{77}}{{60}}{f_i} + \frac{1}{5}{f_{i + 1}}\\ f_{i + 1/2}^2 = \frac{1}{{30}}{f_{i - 2}} - \frac{{13}}{{60}}{f_{i - 1}} + \frac{{47}}{{60}}{f_i} + \frac{9}{{20}}{f_{i + 1}} - \frac{1}{{20}}{f_{i + 2}}\\ f_{i + 1/2}^3 = - \frac{1}{{20}}{f_{i - 1}} + \frac{9}{{20}}{f_i} + \frac{{47}}{{60}}{f_{i + 1}} - \frac{{13}}{{60}}{f_{i + 2}} + \frac{1}{{30}}{f_{i + 3}}\\ f_{i + 1/2}^4 = \frac{1}{5}{f_i} + \frac{{77}}{{60}}{f_{i + 1}} - \frac{{43}}{{60}}{f_{i + 2}} + \frac{{17}}{{60}}{f_{i + 3}} - \frac{1}{{20}}{f_{i + 4}} \end{array} \right.$

(17)

利用上述5种模板的凸组合计算数值通量f_i+1/2，即

${f_{i + 1/2}} = {\omega _0}f_{i + 1/2}^0 + {\omega _1}f_{i + 1/2}^1 + {\omega _2}f_{i + 1/2}^2 + {\omega _3}f_{i + 1/2}^3 + {\omega _4}f_{i + 1/2}^4$

(18)

对于含激波间断流场，(18)式中的ω_k根据下式求得

${\omega _k} = {\alpha _k}/\sum\limits_{s = 0}^3 {{\alpha _s}} ,\;\;\;{\alpha _k} = \frac{{{d_k}}}{{{{\left( {\varepsilon + {\beta _k}} \right)}^2}}},\;\;\;k = 0,1,2,3,4$

(19)

对于九阶WENO格式，d₀= $\frac{1}{{126}}$ ，d₁= $\frac{10}{{63}}$ ，d₂= $\frac{10}{{21}}$ ，d₃= $\frac{20}{{63}}$ ，d₄= $\frac{5}{{126}}$ 。光滑因子β_k(k=0, 1, 2, 3, 4)的表达式如下^[9]

$\left\{ \begin{array}{l} {\beta _0} = {f_{i - 4}}\left( {22658{f_{i - 4}} - 208501{f_{i - 3}} + 364863{f_{i - 2}} - 288007{f_{i - 1}} + 86329{f_i}} \right) + \\ \;\;\;\;\;\;\;{f_{i - 3}}\left( {482963{f_{i - 3}} - 1704396{f_{i - 2}} + 1358458{f_{i - 1}} - 411487{f_i}} \right) + \\ \;\;\;\;\;\;\;{f_{i - 2}}\left( {1521393{f_{i - 2}} - 2462076{f_{i - 1}} + 758823{f_i}} \right) + \\ \;\;\;\;\;\;\;{f_{i - 1}}\left( {1020563{f_{i - 1}} - 649501{f_i}} \right) + 107918f_i^2\\ {\beta _1} = {f_{i - 3}}\left( {6908{f_{i - 3}} - 60871{f_{i - 2}} + 99213{f_{i - 1}} - 70237{f_i} + 18079{f_{i + 1}}} \right) + \\ \;\;\;\;\;\;\;{f_{i - 2}}\left( {138563{f_{i - 2}} - 464976{f_{i - 1}} + 337018{f_i} - 88297{f_{i + 1}}} \right) + \\ \;\;\;\;\;\;\;{f_{i - 1}}\left( {406293{f_{i - 1}} - 611976{f_i} + 165153{f_{i + 1}}} \right) + \\ \;\;\;\;\;\;\;{f_i}\left( {242723{f_i} - 140251{f_{i + 1}}} \right) + 22658f_{i + 1}^2\\ {\beta _2} = {f_{i - 2}}\left( {6908{f_{i - 2}} - 51001{f_{i - 1}} + 67923{f_i} - 38947{f_{i + 1}} + 8209{f_{i + 2}}} \right) + \\ \;\;\;\;\;\;\;{f_{i - 1}}\left( {104963{f_{i - 1}} - 299076{f_i} + 179098{f_{i + 1}} - 38947{f_{i + 2}}} \right) + \\ \;\;\;\;\;\;\;{f_i}\left( {231153{f_i} - 299076{f_{i + 1}} + 67923{f_{i + 2}}} \right) + \\ \;\;\;\;\;\;\;{f_{i + 1}}\left( {104963{f_{i + 1}} - 51001{f_{i + 2}}} \right) + 6908f_{i + 2}^2\\ {\beta _3} = {f_{i - 1}}\left( {22658{f_{i - 1}} - 140251{f_i} + 165153{f_{i + 1}} - 88297{f_{i + 2}} + 18079{f_{i + 3}}} \right) + \\ \;\;\;\;\;\;\;{f_i}\left( {242723{f_i} - 611976{f_{i + 1}} + 337018{f_{i + 2}} - 70237{f_{i + 3}}} \right) + \\ \;\;\;\;\;\;\;{f_{i + 1}}\left( {406293{f_{i + 1}} - 464976{f_{i + 2}} + 99213{f_{i + 3}}} \right) + \\ \;\;\;\;\;\;\;{f_{i + 2}}\left( {138563{f_{i + 2}} - 60871{f_{i + 3}}} \right) + 6908f_{i + 3}^2\\ {\beta _4} = {f_i}\left( {107918{f_i} - 649501{f_{i + 1}} + 758823{f_{i + 2}} - 411487{f_{i + 3}} + 86329{f_{i + 4}}} \right) + \\ \;\;\;\;\;\;\;{f_{i + 1}}\left( {1020563{f_{i + 1}} - 2462076{f_{i + 2}} + 1358458{f_{i + 3}} - 288007{f_{i + 4}}} \right) + \\ \;\;\;\;\;\;\;{f_{i + 2}}\left( {1521393{f_{i + 2}} - 1704396{f_{i + 3}} + 364863{f_{i + 4}}} \right) + \\ \;\;\;\;\;\;\;{f_{i + 3}}\left( {482963{f_{i + 3}} - 208501{f_{i + 4}}} \right) + 22658f_{i + 4}^2 \end{array} \right.$

(20)

2.5 时间项数值离散

欧拉方程时间项采用三阶TVD-RK(Total Variation Diminishing Runge-Kutta)法进行数值离散，具体离散格式如下^[10]

$\left\{ \begin{array}{l} {\varphi ^{\left( 1 \right)}} = {\varphi ^n} + \Delta tL\left( {{\varphi ^n}} \right)\\ {\varphi ^{\left( 2 \right)}} = \frac{3}{4}{\varphi ^n} + \frac{1}{4}{\varphi ^{\left( 1 \right)}} + \frac{1}{4}\Delta tL\left( {{\varphi ^{\left( 1 \right)}}} \right)\\ {\varphi ^{\left( {n + 1} \right)}} = \frac{1}{3}{\varphi ^n} + \frac{1}{3}{\varphi ^{\left( 2 \right)}} + \frac{2}{3}\Delta tL\left( {{\varphi ^{\left( 2 \right)}}} \right) \end{array} \right.$

(21)

式中：φⁿ表示n时刻的守恒通量，φ⁽¹⁾、φ⁽²⁾为中间变量，Δt为时间步长，L为运算算子。

3. 经典算例

为了初步考察上述格式(WENO-JS3、WENO-JS5、WENO-JS7、WENO-JS9)的计算性能并验证所开发程序的可靠性，选取两个经典一维黎曼算例进行计算。

3.1 Sod激波管

该算例初始条件如(22)式所示^[11]，网格数为400，计算结束时间为0.18。图 1为计算结束时刻的密度曲线及其局部放大图。

${\left( {\rho ,u,p} \right)^{\rm{T}}} = \left\{ \begin{array}{l} {\left( {1,0,1} \right)^{\rm{T}}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;0 \le x < 0.5\\ {\left( {0.125,0,0.100} \right)^{\rm{T}}}\;\;\;\;\;\;\;\;\;0.5 \le x \le 1.0 \end{array} \right.$

(22)

图 1 Sod激波管算例密度曲线及其局部放大图

Figure 1. Density curve and its partially enlarged details for Sod shock tube

下载: 全尺寸图片幻灯片

3.2 激波与熵波相互作用

该算例初始条件如(23)式所示，网格数为400，计算结束时间为1.8。图 2为计算结束时刻的密度曲线及其局部放大图。

${\left( {\rho ,u,p} \right)^{\rm{T}}} = \left\{ \begin{array}{l} {\left( {3.857143,2.629369,10.333330} \right)^{\rm{T}}}\;\;\;\;\;\;\;\;\; - 5 \le x < - 4\\ {\left( {1 + 0.2\sin 5x,0.1} \right)^{\rm{T}}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; - 4 \le x \le 5 \end{array} \right.$

(23)

图 2 激波与熵波相互作用算例密度曲线及其局部放大图

Figure 2. Density curve and its partially enlarged details for shock-entropy wave interaction

下载: 全尺寸图片幻灯片

3.3 小结

根据图 1和图 2的计算结果，评估各格式的计算性能发现：WENO-JS9格式的计算性能最优，WENO-JS7、WENO-JS5格式次之，WENO-JS3格式的计算性能最低。即高精度WENO格式对含激波间断流场具有更低的耗散特性和更高的分辨率。

4. 舱室内气云爆炸载荷计算

为了考察不同精度WENO格式对舱室内部气云爆炸载荷的影响规律，采用所开发的三维程序开展封闭舱室和泄压舱室内爆炸载荷数值计算。

4.1 初始条件设置

舱室尺寸为800 mm×800 mm×800 mm，泄压舱室壁面上开有一个边长为160 mm的正方形泄压口。壁面上设置两个测点，分别编号为No.1、No.2，对爆炸超压时间历程进行输出，如图 3所示。

图 3 封闭舱室和泄压舱室及其测点分布(单位：mm)

Figure 3. Closed cabin and venting cabin and their gauging points (Unit:mm)

下载: 全尺寸图片幻灯片

半径为200 mm的球体气云位于舱室中间，气云参数为^[12-13]：密度1.99 kg/m³，压力1.38 MPa。周围空气参数为：密度1.225 kg/m³，压力0.1 MPa。计算初始条件见图 4(a)。考虑计算时间及精度的要求，经多次数值试验，网格数为40×40×40(见图 4(b))。壁面边界条件设置为反射边界，泄压口处边界条件设置为透射边界条件^[14]。

图 4 爆炸初场及网格分布

Figure 4. Initial condition and mesh distribution

下载: 全尺寸图片幻灯片

4.2 舱室内爆炸载荷特性

图 5给出了舱室内部气云爆炸壁面测点No.1和No.2处的超压时间历程曲线。由图 5可知：封闭舱室内爆炸载荷主要包含多峰值、强度逐渐衰弱的冲击波和持续时间较长、峰值保持稳定的准静态超压，泄压舱室内爆炸载荷主要包含多峰值、强度逐渐衰弱的冲击波和持续时间较长、呈指数衰减规律的准静态超压。

图 5 舱室内爆炸超压时间历程曲线

Figure 5. Overpressure histories of all the gauging points inside the cabin

下载: 全尺寸图片幻灯片

根据测点No.1和No.2的对比可知，测点的位置对爆炸前期的爆炸波超压峰值有较大的影响，而对形成的准静态超压时间历程几乎没有影响，说明舱室内部爆炸形成的准静态超压时间历程在空间上是近似均匀的。

4.3 封闭舱室内爆炸载荷比较

选取壁面典型测点No.1对爆炸载荷进行输出，探讨不同精度WENO格式对封闭舱室内气云爆炸载荷的影响规律。

通过对比图 6(a)~图 6(d)发现，高精度WENO格式能分辨出更精细的爆炸载荷特征，即高精度WENO格式具有更低的耗散特性和更高的分辨率。根据图 7可知，高精度WENO格式对于激波间断具有更低的耗散特性，给出更陡峭的激波峰值。然而，不同精度WENO格式对最终形成的准静态超压峰值影响较小，如图 7(a)中的绿色粗实线所示。

图 6 不同精度WENO格式下封闭舱室内气云爆炸测点No.1处的超压时间历程曲线

Figure 6. Overpressure histories at gauging point No.1 for different schemes in closed cabin

下载: 全尺寸图片幻灯片

图 7 不同精度WENO格式下封闭舱室内气云爆炸测点No.1处的超压比较及其局部放大图

Figure 7. Comparison of overpressure at gauging point No.1 for different schemes in closed cabin and its partially enlarged details

下载: 全尺寸图片幻灯片

4.4 泄压舱室内爆炸载荷比较

选取壁面典型测点No.1对爆炸载荷进行输出，探讨不同精度WENO格式对泄压舱室内爆炸载荷的影响规律。

通过对比图 8(a)~图 8(d)发现，高精度WENO格式能分辨出更精细的爆炸载荷特征，即高精度WENO格式具有更低的耗散特性和更高的分辨率。根据图 9可知，高精度WENO格式对于激波间断具有更低的耗散特性，给出更陡峭的激波峰值。然而，不同精度WENO格式对逐渐形成的、呈现指数衰减规律的准静态超压影响较小，如图 9(a)中的绿色粗实线所示。

图 8 不同精度WENO格式下泄压舱室内气云爆炸测点No.1处的超压时间历程曲线

Figure 8. Overpressure histories at gauging point No.1 for different schemes in venting cabin

下载: 全尺寸图片幻灯片

图 9 不同精度WENO格式下泄压舱室内气云爆炸测点No.1处的超压比较及其局部放大图

Figure 9. Comparison of overpressure at gauging point No.1 for different schemes in venting cabin and its partially enlarged details

下载: 全尺寸图片幻灯片

5. 结论

基于自主开发的高精度舱室内部气云爆炸三维数值计算程序研究了WENO格式精度对舱室内部气云爆炸载荷的影响规律，主要得到以下结论。

(1) 高精度WENO格式对含激波间断流场具有更低的耗散特性和更高的分辨率。

(2) 舱室内部气云爆炸载荷主要包含多峰值的瞬态冲击波和持续时间较长的准静态超压；封闭舱室内爆炸形成峰值保持不变的准静态超压，泄压舱室内爆炸形成呈指数衰减的准静态超压。

(3) WENO格式精度对舱室内爆炸冲击波载荷影响较大，高精度WENO格式给出更陡峭的激波峰值，而对舱室内爆炸准静态超压载荷的影响较小。

图 1 Sod激波管算例密度曲线及其局部放大图

Figure 1. Density curve and its partially enlarged details for Sod shock tube

下载: 全尺寸图片幻灯片

图 2 激波与熵波相互作用算例密度曲线及其局部放大图

Figure 2. Density curve and its partially enlarged details for shock-entropy wave interaction

下载: 全尺寸图片幻灯片

图 3 封闭舱室和泄压舱室及其测点分布(单位：mm)

Figure 3. Closed cabin and venting cabin and their gauging points (Unit:mm)

下载: 全尺寸图片幻灯片

图 4 爆炸初场及网格分布

Figure 4. Initial condition and mesh distribution

下载: 全尺寸图片幻灯片

图 5 舱室内爆炸超压时间历程曲线

Figure 5. Overpressure histories of all the gauging points inside the cabin

下载: 全尺寸图片幻灯片

图 6 不同精度WENO格式下封闭舱室内气云爆炸测点No.1处的超压时间历程曲线

Figure 6. Overpressure histories at gauging point No.1 for different schemes in closed cabin

下载: 全尺寸图片幻灯片

图 7 不同精度WENO格式下封闭舱室内气云爆炸测点No.1处的超压比较及其局部放大图

Figure 7. Comparison of overpressure at gauging point No.1 for different schemes in closed cabin and its partially enlarged details

下载: 全尺寸图片幻灯片

图 8 不同精度WENO格式下泄压舱室内气云爆炸测点No.1处的超压时间历程曲线

Figure 8. Overpressure histories at gauging point No.1 for different schemes in venting cabin

下载: 全尺寸图片幻灯片

图 9 不同精度WENO格式下泄压舱室内气云爆炸测点No.1处的超压比较及其局部放大图

Figure 9. Comparison of overpressure at gauging point No.1 for different schemes in venting cabin and its partially enlarged details

下载: 全尺寸图片幻灯片

参考文献(14)

[1]	徐胜利, 汤明钧, 糜仲春.近地空中气云爆炸波遇地面反射的研究[J].爆炸与冲击, 1996, 16(4):298-304. http://www.cqvip.com/QK/94778X/1996004/2291573.html XU S L, TANG M J, MI Z C.Studies on reflection of blast waves for symmetric cloud explosion close to the ground[J].Explosion and Shock Waves, 1996, 16(4):298-304. http://www.cqvip.com/QK/94778X/1996004/2291573.html
[2]	徐胜利, 彭金华.多爆源云雾爆炸波相互作用的三维数值研究[J].爆炸与冲击, 2000, 20(1):1-6. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=bzycj200001001 XU S L, PENG J H.Three dimensional computaion on the interaction of blast waves generated by multi-sources of FAE[J].Explosion and Shock Waves, 2000, 20(1):1-6. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=bzycj200001001
[3]	岳鹏涛, 彭金华.FAE爆炸波对地面目标作用的三维数值研究[J].爆炸与冲击, 2000, 20(2):97-102. http://www.cnki.com.cn/Article/CJFDTOTAL-BZCJ200002000.htm YUE P T, PENG J H.3D numerical simulations on the interaction between FAE blast waves and ground targets[J].Explosion and Shock Waves, 2000, 20(2):97-102. http://www.cnki.com.cn/Article/CJFDTOTAL-BZCJ200002000.htm
[4]	陈明生, 李建平, 白春华.非圆截面云雾爆炸超压场数值模拟[J].含能材料, 2015, 23(5):484-489. doi: 10.11943/j.issn.1006-9941.2015.05.015 CHEN M S, LI J P, BAI C H.Simulation of explosion overpressure distribution for non-circular cross-section cloud[J].Chinese Journal of Energetic Materials, 2015, 23(5):484-489. doi: 10.11943/j.issn.1006-9941.2015.05.015
[5]	LIU X D, OSHER S, CHAN T.Weighted essentially non-oscillatory schemes[J].Journal of Computational Physics, 1994, 115(1):200-212. doi: 10.1006/jcph.1994.1187
[6]	JIANG G S, SHU C W.Efficient implementation of weighted ENO schemes[J].Journal of Computational Physics, 1996, 126(1):202-202. doi: 10.1006/jcph.1996.0130
[7]	SHU C W.Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[M].Heidelberg:Springer Berlin Heidelberg, 1998:325-432.
[8]	SHU C W.High order weighted essentially nonoscillatory schemes for convection dominated problems[J].Siam Review, 2009, 51(1):82-126. doi: 10.1137/070679065
[9]	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. doi: 10.1006/jcph.2000.6443
[10]	SHU C W, OSHER S.Efficient implementation of essentially non-oscillatory shock-capturing schemes[J].Journal of Computational Physics, 1988, 77(2):439-471. doi: 10.1016/0021-9991(88)90177-5
[11]	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. http://www.sciencedirect.com/science/article/pii/0021999178900232
[12]	赵衡阳.气体和粉尘爆炸原理[M].北京:北京理工大学出版社, 1996. ZHAO H Y.Principle of gas and dust explosion[M].Beijing:Beijing Institute of Technology Press, 1996.
[13]	王建, 段吉员, 黄文斌, 等.乙炔-氧气混合气体强爆轰参数的理论估算与实验研究[J].高压物理学报, 2011, 25(4):365-369. http://www.gywlxb.cn/CN/abstract/abstract1386.shtml WANG J, DUAN J Y, HUANG W B, et al.Calculation and experiment of overdriven detonation parameters of C₂H₂-O₂ mixture[J].Chinese Journal of High Pressure Physics, 2011, 25(4):365-369. http://www.gywlxb.cn/CN/abstract/abstract1386.shtml
[14]	徐维铮, 吴卫国.泄压口大小对约束空间爆炸准静态超压载荷的影响规律[J].高压物理学报, 2017, 31(5):619-628. doi: 10.11858/gywlxb.2017.05.016 XU W Z, WU W G.Effects of size of venting holes on the characteristics of quasi-static overpressure in confined space[J].Chinese Journal of High Pressure Physics, 2017, 31(5):619-628. doi: 10.11858/gywlxb.2017.05.016