密闭多舱室内柱形气云爆炸过程二维数值模拟

2018-11-20徐维铮吴卫国

舰船科学技术 2018年11期

徐维铮，吴卫国，况正

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

0 引言

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

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

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

1 控制方程

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

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

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

2 数值方法

2.1 三阶 WENO-JS 格式

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

利用上述2种二阶通量的凸组合计算最终具有三阶精度的数值通量，即

其中，本文参数ε取值为10-6。光滑因子的表达式如下：

2.2 时间项数值离散

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

3 试验验证

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

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

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

图1 计算区域尺寸及网格划分Fig.1 Size of the computational zone

图2 激波绕方腔流动初始条件Fig.2 Initial condition of the planar shock wave over a square cavity

图3 数值与文献[8]中试验结果对比Fig.3 Comparisons between numerical and experimental results[8]

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

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

4.1 舱室尺寸及测点布置

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

图4 多舱室布置及网格划分Fig.4 Size of multi-cabin and mesh distribution

4.2 初始条件及工况设置

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

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壁面测点NO1，NO2爆炸载荷时间历程，从图7可知舱室内爆炸载荷主要呈现峰值不断衰减的冲击波载荷和准静态压力。虽然测点NO1，NO2与爆源处于近似相同的距离，但是测点NO1超压载荷及冲量明显大于测点NO2，通过分析发现这主要是由于测点NO2靠近舱室2和舱室3，当爆炸波通过连接导管传播到舱室2和舱室3内部后将明显减弱测点NO2处的冲击波强度。

图5 初始条件及爆炸工况Fig.5 Initial condition and explosion cases

图6 工况 1 不同时刻爆炸波传播过程Fig.6 Blast wave propagation process at different time for explosion case one

图7 舱室 1 壁面测点 NO1，NO2 爆炸载荷时间历程（工况1）Fig.7 Blast load time histories of gauging point NO1, NO2 in cabin 1 for explosion case one

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

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

图8 舱室 2 壁面测点 NO3，NO4 爆炸载荷时间历程（工况1）Fig.8 Blast load time histories of gauging point NO3, NO4 in cabin 2 for explosion case one

图9 舱室 1 及舱室 2 壁面测点 NO1，NO2，NO3，NO4 爆炸载荷时间历程（工况1）Fig.9 Blast load time histories of gauging point NO1, NO2, NO3,NO4 in cabin 1, 2 for explosion case one

图10 工况 2 不同时刻爆炸波传播过程Fig.10 Blast wave propagation process at different time for explosion case two

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））。

图11 舱室 1 壁面测点 NO1，NO2 爆炸载荷时间历程（工况2）Fig.11 Blast load time histories of gauging point NO1, NO2 in cabin 1 for explosion case two

图12 舱室 2 壁面测点 NO3，NO4 爆炸载荷时间历程（工况2）Fig.12 Blast load time histories of gauging point NO3, NO4 in cabin 2 for explosion case two

4.4.2 爆炸载荷特性分析

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

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

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

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

图13 舱室 3 壁面测点 NO5，NO6 爆炸载荷时间历程（工况2）Fig.13 Blast load time histories of gauging point NO5, NO6 in cabin 3 for explosion case two

图14 舱室 1～舱室 3 壁面测点 NO1，NO2，NO3，NO4，NO5，NO6爆炸载荷时间历程（工况2）Fig.14 Blast load time histories of gauging point NO1, NO2, NO3,NO4, NO5, NO6 in cabin 1,2,3 for explosion case two