0 引　言

工业可燃性气云爆炸事故大多是由弱点火点燃可燃气体引起的，其传播形式大多是以亚音速传播的爆燃波。爆燃波的传播过程很复杂，受环境条件和物理因素的影响极大，当气体在传播过程中由于障碍物阻碍发生湍流加强效应或是在初始时刻由高能量直接起爆，其爆炸模式将由爆燃转为爆轰或直接以爆轰模式发生爆炸。气体发生爆轰以后将在空气中形成强度较大的带有负压区的空气冲击波，对工作人员和周边结构设施将造成较大的损害。针对空中可燃气云、云雾爆炸模拟，近年来徐胜利等学者进行过相应的数值研究^{[1 – 4]}，他们将可燃气云、云雾等效的简化为一团高压气体，而舱室内部可燃气云爆炸研究较少。

爆炸属于强间断问题，其数值模拟对激波捕捉格式提出较严格要求。目前，WENO格式作为一种典型的高精度激波捕捉格式，对流场内的激波间断具有较高的分辨率，适于求解包含激波、膨胀波以及接触间断等复杂结构的流场。LIU等^[5]于1994年提出了WENO（weighted essentially non-oscillation scheme）格式，之后，SHU等^{[6 – 7]}发展了该格式并扩展了其应用。

考虑到WENO格式高精度、稳定性较好的优势，基于FORTRAN平台，采用三阶WENO有限差分格式，自主开发了密闭多舱室内爆炸波二维数值计算程序。利用所开发的程序开展了多舱室内柱形气云爆炸波传播过程及不同舱室内部爆炸载荷特性研究。

1 控制方程

含激波流场采用可压缩欧拉方程进行描述，其具体形式如下：

${\overrightarrow Q _t} + {\overrightarrow E _x} + {\overrightarrow F _y} = 0{\text{。}}$

(1)

其中：

$\overrightarrow Q = \left( {\begin{array}{*{20}{c}} \rho \\ {\rho u} \\ {\rho v} \\ E \end{array}} \right){\rm{, }}\overrightarrow E = \left( {\begin{array}{*{20}{c}} {\rho u} \\ {\rho {u^2} + p} \\ {\rho uv} \\ {u\left( {E + p} \right)} \end{array}} \right){\rm{, }}\overrightarrow F = \left( {\begin{array}{*{20}{c}} {\rho v} \\ {\rho vu} \\ {\rho {v^2} + p} \\ {v\left( {E + p} \right)} \end{array}} \right),$

(2)

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

(3)

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

(4)

式中： $\rho $ 为密度； $u,v$ 为 $x,y$ 方向上的速度分量； $p$ 为流体压力； $E$ 为单位体积流体的总能量； $e$ 为比内能； $\gamma $ 为气体的绝热指数，本文取为1.4。

在方程（1）的每个方向上均可以看成是一个双曲守恒律方程：

$\frac{{\partial u}}{{\partial t}} + \frac{{\partial f\left( u \right)}}{{\partial x}} = 0{\text{，}}$

(5)

例如，针对x方向，方程（5）的数值离散形式为：

$ \begin{split} & \frac{{{\rm d}{u_i}}}{{{\rm d}t}} = - \frac{{\partial f}}{{\partial x}}\left| {_{x = {x_i}}} \right. = - \frac{{{h_{i + {1/2}}} - {h_{i - {1/2}}}}}{{\Delta x}} \simeq \\ &- \frac{1}{{\Delta x}}\left( {{{\overline f }_{i + {1/2}}} - {{\overline f }_{i - {1/2}}}} \right){\text{，}} \end{split}$

(6)

其中， ${\overline f _{i - {1/2}}}$ ， ${\overline f _{i + {1/2}}}$ 分别为单元（x_i-1/2，x_i+1/2）的左右界面处对流项数值通量， $\Delta x$ 为x方向的均匀网格间距，在本文的数值计算中，x，y两个方向的网格间距相同。

2 数值方法 2.1 三阶WENO-JS格式

三阶WENO格式的数值离散和推导过程如下^[6]：

右界面数值通量 ${\overline f _{i + {1/2}}}$ 的2种二阶精度通量重构方式分别为

$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}}{\text{。}}$

(7)

利用上述2种二阶通量的凸组合计算最终具有三阶精度的数值通量 ${\overline f _{i + {1/2}}}$ ，即

${\overline f _{i + {1/2}}} = \sum\limits_{k = 0}^1 {{\omega _k}f_{i + {1/2}}^k} = {\omega _0}f_{i + {1/2}}^0 + {\omega _1}f_{i + {1/2}}^1{\text{。}}$

(8)

对于光滑情形，有 ${\omega _0} = {d_0} = \displaystyle\frac{1}{3}$ ， ${\omega _1} = {d_1} = \displaystyle\frac{2}{3}$ 。式（8）所给出的形式既适合光滑流场也适合含间断流场，对于含激波间断流场，式中的非线性权函数 ${\omega _k}$ 需要根据下式求得：

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

(9)

其中，本文参数ε取值为10^–6。光滑因子 ${\beta _k}\left( {k = 0,1} \right)$ 的表达式如下：

${\beta _0} = {\left( {{f_{i - 1}} - {f_i}} \right)^2},{\beta _1} = {\left( {{f_i} - {f_{i + 1}}} \right)^2}{\text{。}}$

(10)

2.2 时间项数值离散

采用三阶TVD-RK法对欧拉方程（1）中的时间项进行数值离散，具体的离散格式如下^[6]：

$\left\{ \begin{aligned}& {u^{\left( 1 \right)}} = {u^n} + \Delta tL\left( {{u^n}} \right) {\text{，}}\\& {u^{\left( 2 \right)}} = \frac{3}{4}{u^n} + \frac{1}{4}{u^{\left( 1 \right)}} + \frac{1}{4}\Delta tL\left( {{u^{\left( 1 \right)}}} \right) {\text{，}}\\& {u^{n + 1}} = \frac{1}{3}{u^n} + \frac{1}{3}{u^{\left( 2 \right)}} + \frac{2}{3}\Delta tL\left( {{u^{\left( 2 \right)}}} \right) {\text{。}}\end{aligned} \right.$

(11)

3 试验验证

为了检验所开发数值程序的可靠性和正确性，选取文献[8]开展的平面激波绕方腔流动试验进行对比分析。图1给出计算域的尺寸和网格划分示意图（本文计算中，网格尺寸为1 mm）。

在标准状况下，平面激波马赫数为1.3，根据冲击波物理可以确定波后物理参数。图2给出该算例初始条件。

图3给出不同时刻数值模拟纹影图与试验的对比，图中不同时刻的波系结构和台阶处绕射形成的漩涡结构均证明了文中程序开发的可靠性和正确性，为进行多舱室内爆炸模拟奠定了基础。

4 多舱室内爆炸过程数值计算

本节主要采用验证的程序对密闭多舱室内平面柱形高压气云爆炸过程进行数值模拟研究，主要探讨爆炸波传播路径及不同舱室内部爆炸载荷特性，为后期的结构抗爆设计及毁伤评估提供一定的参考价值。

4.1 舱室尺寸及测点布置

多舱室有3个舱室组成，分别标记为舱室1、舱室2、舱室3，舱室1与舱室2及舱室3之间通过方形连接导管相连通，每个舱室的具体尺寸参见 图4。每个舱室的壁面中点处均设置2个压力测点（分别编号为NO₁，NO₂，NO₃，NO₄，NO₅，NO₆）对爆炸超压载荷时间历程进行输出（见图4（a））。壁面边界条件设置为刚性边界（不考虑舱室结构的变形），均匀网格尺寸取为1 mm（见图4（b））。

4.2 初始条件及工况设置

柱形均匀高压气云的具体参数：半径为20 mm，密度为1.99 kg/m³，压力为13.8×10⁵ Pa。为了探讨爆炸发生位置对爆炸过程的影响规律，本文选取2种典型爆炸工况：工况1表示爆炸发生在舱室1内部，如图5（a）所示，工况2表示爆炸发生在舱室3内部，如图5（b）所示。均匀高压气云均设置在舱室中间（见图5中的浅色区域），周围深色区域表示空气域，空气密度为1.225 kg/m³，压力为1.0×10⁵ Pa。

4.3 工况一爆炸过程分析

本节主要根据数值模拟结果分析工况1爆炸后爆炸波传播过程及不同舱室爆炸载荷分布特性。

4.3.1 爆炸波传播过程分析

为了探讨工况1气云爆炸后爆炸波在多舱室内部的传播过程，这里给出不同时刻爆炸压力云图和数值纹影图（见图6）。根据图6可知工况1多舱室内爆炸过程爆炸波的传播路径：爆炸波首先在舱室1中进行自由膨胀运动（图6（a））爆炸波初次到达舱室壁面时发生正规则反射，壁面处形成局部高压区域（图6（b）），之后爆炸波沿着连接导管进入舱室2和舱室3内部（图6（c）），并紧接着在舱室2和舱室3上下壁面处进行系列正规则斜反射和马赫反射（图6（d）和图6（e）），当爆炸波传播到舱室2和舱室3的最右端壁面处时再次形成复杂的壁面反射波（图6（f））。由于多舱室处于封闭空间，因此可以预知爆炸波经过上述在不同舱室之间的传播过程，爆炸波强度将不断衰弱，爆炸气体的动能逐渐转化为气体的内能，最终舱室内部将趋近于一个准静态过程。

4.3.2 爆炸载荷特性分析

为了探讨工况1气云爆炸后不同舱室内部爆炸载荷分布特性，这里给出不同舱室内部壁面测点超压时间历程曲线和冲量历程时间曲线。针对工况1，由于舱室2和舱室3上下对称，因此，这里只需要讨论舱室1和舱室2内部壁面测点的爆炸载荷特性。图7给出舱室1壁面测点NO₁，NO₂爆炸载荷时间历程，从图7可知舱室内爆炸载荷主要呈现峰值不断衰减的冲击波载荷和准静态压力。虽然测点NO₁，NO₂与爆源处于近似相同的距离，但是测点NO₁超压载荷及冲量明显大于测点NO₂，通过分析发现这主要是由于测点NO₂靠近舱室2和舱室3，当爆炸波通过连接导管传播到舱室2和舱室3内部后将明显减弱测点NO₂处的冲击波强度。

图8给出舱室2壁面测点NO₃，NO₄爆炸载荷时间历程，从图8可知在爆炸初期，测点NO₃超压载荷及冲量稍微大于测点NO₄，而在爆炸后期，测点NO₃超压载荷及冲量演变为小于测点NO₄，通过分析发现这主要是由于爆炸初期测点NO₃相较测点NO₄更靠近爆源，因此其接触到的初始冲击波强度较大。而在爆炸后期由于测点NO₄与连接导管正对，在后续的爆炸波传播过程中，其将会受到爆炸波的连续直接冲击，反而造成其超压载荷、冲量时间历程稍微高于测点NO₃。

图9给出舱室1，舱室2壁面测点NO₁，NO₂，NO₃，NO₄爆炸载荷时间历程，从图9可知舱室1内部测点NO₁，NO₂爆炸载荷强度明显强于舱室2内部测点NO₃，NO₄，通过分析发现这主要是由于爆炸发生在舱室1内部，因此其主要能量集中在舱室1内部。同时根据图9（b）中的冲量历程曲线还得知，不同舱室内部最终的准静态压力相同，如果将冲量积分曲线开始变为直线段的转折点作为准静态压力开始形成的时间，则在该爆炸工况下不同舱室内部测点到达准静态压力的时间基本相同（近似为0.5 ms）。

4.4 工况2爆炸过程分析

本节主要根据数值模拟结果分析工况2爆炸后爆炸波传播过程及不同舱室爆炸载荷分布特性。

4.4.1 爆炸波传播过程分析

为了探讨工况2气云爆炸后爆炸波在多舱室内部的传播过程，这里给出不同时刻爆炸压力云图和数值纹影图（见图10）。根据图10可知工况2多舱室内爆炸过程爆炸波的传播路径：爆炸波首先在长方形舱室3中进行自由膨胀运动，初次到达舱室上下壁面时发生正规则反射，壁面处形成局部高压区域（图10（a）），之后爆炸波沿着连接导管进入舱室1内部并以球形方式进行扩张（图10（b）和图10（c）），在扩张过程中由于舱室1上壁面的限制，形成局部马赫反射波（图10（d）），紧接着，爆炸波传播到舱室2内部并形成复杂的壁面反射波（图10（e）和图10（f））。

4.4.2 爆炸载荷特性分析

为了探讨工况2气云爆炸后不同舱室内部爆炸载荷分布特性，这里给出不同舱室内部壁面测点超压时间历程曲线和冲量历程时间曲线。图11给出舱室1壁面测点NO₁，NO₂爆炸载荷时间历程，从图11可知测点NO₁超压载荷及冲量稍微大于测点NO₂，通过分析发现这主要是由于测点NO₁靠近舱室3的连接导管，当爆炸波通过连接导管传播到舱室1内部时，测点NO₁将遭受到更强的爆炸波冲击。

图12给出舱室2壁面测点NO₃，NO₄爆炸载荷时间历程，从图12可知测点NO₃，NO₄虽然超压载荷时间历程差异较大，但是冲量载荷历程基本相同，通过分析发现这主要是由于测点NO₃，NO₄位于舱室2内部，相较于舱室1中的爆源处于远场区域，当爆炸波通过连接导管扩张到舱室1再传播到舱室2内部后，由于舱室2容积也较舱室1偏小，传入的爆炸波冲击很快使得舱室2内部达到准静态过程。

图13给出舱室3壁面测点NO₅，NO₆爆炸载荷时间历程，从图13可知测点NO₆初始超压载荷明显强于测点NO₅，后期强度衰减速度快于测点NO₅，这主要是由于舱室3为长方形舱室，测点6距离爆源较近，因此初期冲击波载荷强度大，由于测点6距离连接导管近，因此其后期超压衰减速率更快。舱室3容积较小，虽然测点NO₆距离连接导管较近，对该处的冲量历程有一定的影响，但是测点6距离爆源较近，最终使得NO₅，NO₆冲量历程曲线基本重合。

图14给出舱室1~舱室3壁面测点NO₁~NO₆爆炸载荷时间历程，从图14可知舱室3内部测点NO₅，NO₆爆炸载荷强度明显强于舱室1内部测点NO₁，NO₂，通过分析发现这主要是由于爆炸发生在舱室3内部，因此其主要能量集中在舱室3内部。舱室1内部测点NO₁，NO₂爆炸载荷强度稍微强于舱室2内部测点NO₃，NO₄，通过分析发现这主要是由于NO₁，NO₂相较测点NO₃，NO₄距离爆炸源更近一些，冲击波强度随着传播距离的增大，强度逐渐衰弱。同时根据图14（b）中的冲量历程曲线还得知，尽管爆炸波达到各测点的时间存在差异性，但是该爆炸工况下不同舱室内部测点到达准静态压力的时间基本相同（近似为0.5 ms），这一结论与4.3.2小节中的分析结果相同。因此，通过这2个工况的爆炸分析，可以初步得出如下结论：对于特定的封闭舱室，准静态压力形成时间主要与舱室容积和爆源能量相关。其实，从准静态压力形成的物理过程来分析也很容易阐述这种关系，由于准静态压力全场均匀，因此其在全场空间内必然应同时到达。

4.5 工况1与工况2对比分析

为了进一步分析爆源位置对多舱室内部不同测点处爆炸载荷特性的影响规律，图15～图20分别给出2种工况下测点NO₁~NO₆的超压载荷和冲量时间历程曲线图。

根据图15和图16可知，针对舱室1内部测点NO₁，NO₂，工况1相较工况2具有更强的超压载荷和更大的冲量，这主要是由于工况1中爆炸发生在舱室1中。

根据图17和图18可知，针对舱室2内部测点NO₃，NO₄，2种工况超压载荷历程差异较大，但是冲量载荷历程差别较小，这主要是由于舱室2中的测点相较2种工况中的爆源距离较大，属于远场的范围，冲量的历程主要由准静态压力决定。

根据图19和图20可知，针对舱室3内部测点NO₅，NO₆，工况2相较工况1具有更强的超压载荷和更大的冲量，这主要是由于工况2中爆炸发生在舱室3中。

5 结　语

采用三阶WENO有限差分格式，自主开发了密闭多舱室内爆炸波二维数值计算程序。利用所开发的程序开展了多舱室内柱形气云爆炸波传播过程及不同舱室内部爆炸载荷特性研究。通过本文的研究主要得到以下几点结论：

1）多舱室内爆炸过程中，爆炸所在舱室内部爆炸载荷强度明显强于相邻舱室，尤其需要重点关注爆炸所在舱室通过连接导管对直接相邻舱室爆炸载荷的影响；爆源位置对不同舱室内部超压载荷时间历程影响较大，而对处于远场的舱室内冲量时间历程影响较小。

2）对于特定的封闭舱室，准静态压力峰值及准静态压力形成时间均全场均匀，准静态压力形成时间主要与舱室容积和爆源能量相关。

3）本文研究可为工程抗爆结构设计及爆炸毁伤评估提供一定的参考和指导。