0 引言

水倾倒入燃烧油中的危险性历来为人所知，英国皇家学院的Peter博士^[1]发布在国际互联网上的演示视频，则充分显示了这种危险性的严重程度。图 1给出了Peter所做实验的视频截图。水刚一倒入油中，燃烧火苗就急剧增长，形成远远超过油料稳定燃烧时的火焰高度，造成恐怖性后果。

相似问题，重质油品燃烧时遇水喷溅^[2-7]，是消防安全中的重要现象。1989年发生重大人员伤亡和财产损失的山东黄岛油库火灾事故，即由此造成。重质油品储油罐发生火灾事故时，轻质馏分会在油池表面进行稳定燃烧，高温重质馏分则会不断下沉，形成热波。当热波头与沉在油罐底部的水相遇时，水迅速蒸发汽化向上升腾，并夹带大量油品从油罐顶部冲出，形成非常剧烈的喷溅燃烧。

液态水遇到高温燃油，会瞬时在交界面形成沸腾型汽泡核，这些汽泡核快速增长，产生极高的泡内压力，并快速爆裂，这种现象叫做爆发沸腾，特别剧烈的情况也会被称为蒸汽爆炸。

由于这一现象非常重要，引起了广泛关注。Wohletz^[8]研究表明当液滴在高温液体中迅速膨胀时，液滴与周围液体分裂满足C-J爆轰条件。Manzello^[9]等通过实验研究了不同的液滴撞击速度、不同的液体表面温度，以及不同液池深度对撞击过程中形成的液柱高度、发生飞溅的临界韦伯数等的影响，发现水滴在220 ℃花生油内部发生了蒸汽爆炸现象。陈萍萍等^[10-11]也发现了类似现象，而且油品温度越高，蒸汽爆炸越剧烈。周源^[12]研究了水滴滴入熔融金属液体中，水滴破碎的情况，由于金属密度远大于水，金属液体内只出现了一个凹陷，并没有造成金属液滴喷溅。Koseki^[13]的实验表明喷溅现象发生后，热辐射是油池稳定燃烧的2.5~3倍，火焰高度是4~5倍，火焰温度上升约300 ℃。梁志桐^[14]进行了不同油水占比情况下的喷溅实验，测试了油层、油水混合层以及水层的温度变化。

计算机技术的飞速发展，使得数值模拟已经成为科学研究的重要手段，但是对于水倾倒入燃烧油中这样复杂现象的完全模拟，则还没有到达可行阶段。下落撞击和两相燃烧这两个过程的数值模拟相对成熟，具备数值模拟的可能，但是耦合汽化沸腾、喷溅蒸发这样剧烈的多时间尺度多空间尺度问题之后，完全数值模拟在现阶段就变得不可实现，因此为了能够推动这项研究，必须发展模型化方法，简化其中难以模拟的过程。

本文旨在建立数值模拟这一现象的简化模型，拟为进一步的研究奠定技术基础。首先，分析了倾倒入燃烧油中发生的物理化学过程；其次，对油滴喷溅建立了简化数学模型；然后采用数值计算的方法进行验证与分析。

1 现象分析 1.1 下落撞击

液态水倾倒后受重力作用，不断加速下落，与油池表面撞击，造成液态水的形态破碎。理想化的情况下，液态水会在撞击点部分穿透部分反弹。反弹的液态水高度很低，快速下落，然后再次与燃油掺混。考虑到液态水的密度略高于燃油密度，因此液态水会在重力作用下深入油面下方，缓慢下沉。而油池表面受到撞击则会凹陷，液面振荡，并带动整个油体的运动。当然在一定条件下，油面也会发生形态破碎，燃油涌溅。这个过程的主要影响因素是液态水的倾倒高度，这控制着液态水与油面撞击时的速度，决定了整个掺混过程的形态和规模。

1.2 汽化沸腾

油温高于液态水的沸腾温度，因此水油接触后会立即在接触面形成汽膜，汽膜受热膨胀，极速分裂破碎，形成爆炸式的沸腾现象。这一过程会将燃油炸离油池，在其上方形成弥漫的细碎油滴。这一过程极其复杂迅速，只要油温和油量足够，所有的液态水最终都会沸腾汽化，爆炸升腾，而油面激荡起伏，液滴飞腾。油温是这一过程的关键因素，它控制着汽化速率，决定了整个过程的持续时间。

1.3 喷溅蒸发

自油面飞腾而出的油滴，忽略气流曳力的情况下，主要受重力控制，因此油滴飞行减速并最终下落。在此过程中，油面燃烧产生的热气流使得油滴快速蒸发汽化，较小的油滴消耗殆尽，而较大的油滴会有液相残存，因此油面上方是气相油蒸汽和液相微小油滴共存的状态。这一过程的关键因素很多，油滴大小分布、飞腾速度分布、油的密度和饱和蒸气压等等都对油池上方的多相组成起重要作用。

1.4 两相燃烧

液态水倾倒而下自然将火焰从撞击区域排开，然后汽化沸腾，喷溅蒸发，而油面迅速回火，引燃上方气相油蒸汽和液相微小油滴，造成火焰形状和规模的急剧变化。待上方可燃物质消耗殆尽后，火焰才重新回归到油池表面稳定燃烧，整个爆发过程戛然而止。这一过程是前述过程的自然延续，其关键特征，如火焰爆发规模和火焰持续时间都是由前述过程实质决定的。

2 油滴喷溅模型

对于爆发沸腾这类现象，相关的实验研究和理论研究都有一些成果，基于分子运动论的数值模拟也获得了一些进展，但是详细可靠的微观机理还没有确立，因此还不可能实现耦合流体运动的直接数值模拟。同时，本文研究关注的是爆发火焰的规模和持续时间，并不把重点放在爆发沸腾过程的模拟上，因此可以将下落撞击、汽化沸腾和喷溅蒸发三个过程合并简化为

$ {M_{{\text{oil}}}} = C{M_w} $

(1)

式中：M_oil代表飞腾到油池之上的燃油质量，M_w代表倾倒入油池的液态水质量，系数C则直接代表单位质量液态水喷溅蒸发的燃油质量。

这样的一个模型非常简单，但抓住了物理问题数学描述的关键量，控制火焰爆发规模和持续时间的就是喷溅而起的油滴总质量，而这些油滴蒸发燃烧形成了猝发火焰。显然系数C目前无法通过理论推导或者数值方法得到，其也不可能是一个固定数值，一定是一个区间变化的随机量，具体数值会受到倾倒高度、油温、油的密度和黏性等因素的影响。本项研究计划通过真实试验确定系数C的取值范围，但本文数值模拟工作先期开展，以给试验研究提供方案设计的数据基础，因此本文中将系数C设定为若干值，考查了参数变化的影响程度。

油滴总质量之外，还需要假设油滴在空间中的分布。其实，燃油液滴从油池表面喷溅而出，穿越高温火焰或者火焰上部的热气流，会即时蒸发，形成弥漫的油蒸汽，并被热气流裹挟至高处，只有较大的油滴才可能在较长的时间内保持液态残留。但是假设油蒸汽的分布是不现实的，这需要在一个三维空间设定一个连续函数表征油蒸汽含量，而这个函数不可能是常值，必须指定空间分布，就目前而言，我们还不知道油蒸汽分布的任何信息，因此任何设定都是极不可靠的。

作为替代，假设油滴的空间分布相对可靠，毕竟爆发沸腾喷溅而出的是实实在在的油滴，将来真实实验时也可通过技术手段，测得油滴尺寸、油滴速度、喷溅角度、喷溅高度等关键信息。图 2给出了本文设定的油滴分布，这个分布是周向均匀的，底部和中央油滴密集，油滴体积也较大，顶部和侧部油滴分布逐渐稀疏，油滴体积也逐渐减小。这样的一个分布是在分析现象基础上作出的，与真实情况的吻合程度需要试验检验和修正。

3 数值计算

两相燃烧过程数值模拟是采用开源FDS(Fire Dynamics Simulator)软件完成的。这一软件是美国国家标准研究所建筑火灾研究实验室开发的模拟火灾中流体运动的计算流体动力学软件。该软件采用数值方法求解受火灾浮力驱动的低马赫数流动的Navier-Stokes方程，重点计算火灾中的烟气和热传递过程^[15]。

FDS所求解的基本控制方程如下：

连续性方程：

$ \frac{{\partial \rho }}{{\partial t}} + \nabla \cdot \left( {\rho \mathit{\boldsymbol{u}}} \right) = {{\dot m'''}_b} $

(2)

动量守恒方程：

$ \frac{{\partial \left( {\rho \mathit{\boldsymbol{u}}} \right)}}{{\partial t}} + \nabla \cdot \rho \mathit{\boldsymbol{uu}} + \nabla \tilde p = \rho g + {f_b} + \nabla \cdot \tau $

(3)

组分输运方程：

$ \frac{{\partial \left( {\rho {Z_\beta }} \right)}}{{\partial t}} + \nabla \cdot \left( {\rho {Z_\beta }\mathit{\boldsymbol{u}}} \right) = \nabla \cdot \left( {\rho {D_\beta }\nabla {Z_\beta }} \right) + {{\dot m'''}_\beta } $

(4)

能量输运方程：

$ \begin{gathered} \frac{{\partial \rho {h_s}}}{{\partial t}} + \nabla \cdot \left( {\rho {h_s}\mathit{\boldsymbol{u}}} \right) = \frac{{{\text{D}}p}}{{{\text{D}}t}} - \hfill \\ \nabla \cdot \left( { - k\nabla T - \sum\limits_\beta {{h_{s,\beta }}\rho {D_\beta }\nabla {Z_\beta }} + {{\dot q'''}_r}} \right) + \dot q''' + {{\dot q'''}_b} \hfill \\ \end{gathered} $

(5)

理想气体状态方程：

$ \bar p = \frac{{\rho RT}}{{\bar W}} $

(6)

低速热浮力驱动压力表示方法：

$ p = \bar p + \tilde p $

(7)

式中：ρ为气体密度，u为气体速度，${{{\dot{m}}'''}_{b}}$为液滴蒸发的总质量，D_β为组分β扩散系数，k为导热系数，f_b为液滴阻力，p为气体压力，p为背压，$\tilde{p}$为扰动压力，Z_β为组分混合分数，R为理想气体常数，g为当地重力加速度，τ为黏性力张量，h_s为焓值，${{{\dot{q}}''}_{r}}$、${{{\dot{q}}'''}_{b}}$、${\dot{q}}'''$分别为辐射热通量、液滴传热量以及单位体积化学反应热释放率，T为气体温度，W为气体平均摩尔质量。

FDS采用燃烧模型为非预混混合组分反应模型，化学反应速率无限快，定义混合分数为独立于反应的守恒标量，即所有物质质量分数和温度仅是混合分数的函数，化学反应方程式如下：

$ \begin{gathered} {\nu _f}{{\text{C}}_a}{{\text{H}}_b}{{\text{O}}_m}{{\text{N}}_n} + {\nu _a}\left( {{v_{{{\text{O}}_2},1}}{{\text{O}}_2} + {v_{{{\text{N}}_2},1}}{{\text{N}}_2} + } \right. \hfill \\ \left. {{v_{{{\text{H}}_2}{\text{O}},1}}{{\text{H}}_2}{\text{O}} + {v_{{\text{C}}{{\text{O}}_2},1}}{\text{C}}{{\text{O}}_2}} \right) \to \hfill \\ {\nu _p}\left( {{v_{{\text{C}}{{\text{O}}_2},2}}{\text{C}}{{\text{O}}_2} + {v_{{{\text{H}}_2}{\text{O}},2}}{{\text{H}}_2}{\text{O}} + {v_{{{\text{N}}_2},2}}{{\text{N}}_2} + } \right. \hfill \\ \left. {{v_{{\text{CO}},2}}{\text{CO}} + {v_{S,2}}{\text{Soot}}} \right) \hfill \\ \end{gathered} $

(8)

式中：ν表示参加反应的集总数，ν的下标“f”表示燃料; υ表示反应物和生成物的计量数，υ的下标“1”表示反应前空气中各组分的计量数，下标“2”表示反应后气相中的各组分计量数; “a”表示空气; “p”表示产物。

$ \dot q''' \equiv \Delta {h_\beta }{{\dot m'''}_{\text{O}}} $

(9)

${{{\dot{m}}'''}_{\text{O}}}$为氧气的质量消耗率，可由组分输运方程(4)求得。

湍流采用Smagorinsky形式的大涡模拟方法，以过滤后的动量方程为例：

$ \frac{{\partial \bar \rho {{\tilde u}_i}}}{{\partial t}} + \frac{\partial }{{\partial {x_j}}}\left( {\bar \rho {{\tilde u}_i}{{\tilde u}_j}} \right) = - \frac{{\partial \bar p}}{{\partial {x_i}}} - \frac{{\partial \tau _{ij}^{{\text{dev}}}}}{{\partial {x_j}}} + \bar \rho {g_i} + {{\bar f}_b} $

(10)

$ \begin{gathered} \tau _{ij}^{{\text{dev}}} \equiv {{\bar \tau }_{ij}} + \tau _{ij}^{{\text{sgs}}} - \frac{1}{3}\tau _{kk}^{{\text{sgs}}}{\delta _{ij}} = \hfill \\ \;\;\;\;\;\; - 2\left( {\mu + {\mu _t}} \right)\left( {{{\tilde S}_{ij}} - \frac{1}{3}\left( {\nabla \cdot \tilde u} \right){\delta _{ij}}} \right) \hfill \\ \end{gathered} $

(11)

式中：τ_ij为过滤后的应力张量，τ_ij^sgs为亚格子应力项，τ_kk^sgs为亚格子动能项，${{\tilde{S}}_{ij}}$为应力偏量，δ_ij为Kolmogorov符号常数(当i=j时, δ_ij=1;当i≠j时, δ_ij=0)。μ_t为湍流黏度：

$ {\mu _t} = \rho {\left( {{C_s}\Delta } \right)^2}{\left[ {2{{\tilde S}_{ij}}{{\tilde S}_{ij}} - \frac{2}{3}{{\left( {\nabla \cdot \tilde u} \right)}^2}} \right]^{0.5}} $

(12)

式中：C_s为常数0.2，Δ=V_c^1/3, V_c=δxδyδz为过滤尺度。

FDS中标量和矢量分别在网格单元的中心和网格面上求解，第n个时间步长网格中心点的标量表示为ϕ_{i, j, k}ⁿ，时间离散上采用二阶预估校正法，以下用组分输运方程的离散为例说明：

$ \begin{gathered} \left( {\rho Z} \right)_\beta ^{{\text{es}}} = {\rho ^n}Z_\beta ^n - \delta t\left[ {\nabla \cdot {\rho ^n}Z_\beta ^n{\mathit{\boldsymbol{u}}^n} + } \right. \hfill \\ \;\;\;\left. {\nabla \cdot \left( {{\rho ^n}D_\beta ^n\nabla Z_\beta ^n} \right) + {{\left( {{{\dot m'''}_\beta } + {{\dot m'''}_{b,\beta }}} \right)}^n}} \right] \hfill \\ \end{gathered} $

(13)

$ \begin{gathered} \left( {\rho Z} \right)_\beta ^{(n + 1)} = \frac{1}{2}\left( {{\rho ^n}Z_\beta ^n + {\rho ^{{\text{es}}}}Z_\beta ^{{\text{es}}}} \right) - \frac{{\delta t}}{2}\left[ {\nabla \cdot {\rho ^{{\text{es}}}}Z_\beta ^{{\text{es}}}{\mathit{\boldsymbol{u}}^{{\text{es}}}}} \right. \hfill \\ \;\;\;\left. { + \nabla \cdot \left( {{\rho ^{{\text{es}}}}D_\beta ^{{\text{es}}}\nabla Z_\beta ^{{\text{es}}}} \right) + {{\left( {{{\dot m'''}_\beta } + {{\dot m'''}_{b,\beta }}} \right)}^n}} \right] \hfill \\ \end{gathered} $

(14)

空间上采用结构网格二阶有限差分方法:

$ \begin{gathered} {\left[ {\frac{{\delta \left( {\rho {Z_\beta }u} \right)}}{{\delta x}}} \right]_i} = \frac{{{{\overline {\left( {\rho {Z_\beta }} \right)} }_{i + \frac{1}{2}}}{u_{i + \frac{1}{2}}} - {{\left( {\rho {Z_\beta }} \right)}_{i - \frac{1}{2}}}{u_{i - \frac{1}{2}}}}}{{\delta x}} \hfill \\ \;\; = {\left( {\rho {Z_\beta }} \right)_i}\frac{{{u_{i + \frac{1}{2}}} - {u_{i - \frac{1}{2}}}}}{{\delta x}} + \frac{{{\Delta _{i + \frac{1}{2}}}{u_{i + \frac{1}{2}}} + {\Delta _{i - \frac{1}{2}}}{u_{i - \frac{1}{2}}}}}{{\delta x}} \hfill \\ \end{gathered} $

(15)

式中的${{\Delta }_{i+\frac{1}{2}}}$代表标量数据传输的有限斜率，由下式表示：

$ \begin{gathered} {\Delta _{i - \frac{1}{2}}} = \frac{1}{2}\left[ {{{\left( {\rho {Z_\beta }} \right)}_i} - {{\left( {\rho {Z_\beta }} \right)}_{i - 1}}} \right] \hfill \\ {\Delta _{i + \frac{1}{2}}} = \frac{1}{2}\left[ {{{\left( {\rho {Z_\beta }} \right)}_{i + 1}} - {{\left( {\rho {Z_\beta }} \right)}_i}} \right] \hfill \\ \end{gathered} $

(16)

可得：

$ \begin{gathered} {\left[ {\frac{{\delta \left( {\rho {Z_\beta }u} \right)}}{{\delta x}}} \right]_i} = {\left( {\rho {Z_\beta }} \right)_i}\frac{{{u_{i + \frac{1}{2}}} - {u_{i - \frac{1}{2}}}}}{{\delta x}} + \hfill \\ \;\;\;\frac{1}{2}{u_{i + \frac{1}{2}}}\left[ {\frac{{{{\left( {\rho {Z_\beta }} \right)}_{i + 1}} - {{\left( {\rho {Z_\beta }} \right)}_i}}}{{\delta x}}} \right] - \hfill \\ \;\;\;\frac{1}{2}{u_{i - \frac{1}{2}}}\left[ {\frac{{{{\left( {\rho {Z_\beta }} \right)}_i} - {{\left( {\rho {Z_\beta }} \right)}_{i - 1}}}}{{\delta x}}} \right] \hfill \\ \end{gathered} $

(17)

本文以室外池燃烧加水起火的危险演示图像资料做参考根据，油池尺寸选定为0.3 m×0.3 m×0.1 m，油池壁面为金属铁，厚度为0.003 m。整个计算空间选定为4 m×4 m×10 m，其四面及上方为开放环境，无风，环境温度25 ℃，气压1×10⁵ Pa。液态燃料选定亚麻酸，这是食用油中的一种主要成分，池内燃料高度设定为0.09 m。起始状态油池下方壁面温度设定为1800 ℃，四壁为800 ℃，待油池稳定燃烧后关闭油池壁面加热。倾倒水的质量假设为500 ml，即0.5 kg，单位质量液态水喷溅蒸发的燃油质量系数分别设定为1.2、0.9和0.6，每种情况下的燃油总质量并不是单次完全注入，而是短时间内分成三次注入，以刻画喷溅蒸发的持续时间尺度。图 3示意了计算空间以及油滴分布的相对尺度。

4 结果分析

图 4(a)给出了油池稳态燃烧的火焰形状图，清晰显示了稳态燃烧时的火焰规模和形状，这种图是根据热释放率绘出的，图中所示红区热释放率均超过200 kW/s。一般而言，稳态火焰可分为连续相火焰区，间断火焰区以及浮力羽流区三个区域，每段区域温度都不同，连续相区域火焰温度最高，间断区温度出现上下抖动，焰羽区温度较低。贴近油池上方有一个突出的火焰面，再往上火焰面变窄，这是由于空气从火焰表面进行卷吸，在富油区上方发生了燃烧，使火焰宽度从下往上逐渐变窄。由于是四面及上方敞开体系，热气流上升之后，压力降低，因此四周向体系内补充空气，使得火焰形状始终保持朝上。

图 4(b)~图 4(f)给出了油滴注入后五个时刻瞬态燃烧的火焰形状，其中油滴质量系数为0.9，油滴注入起始时刻为0时刻。很清楚，油池上方瞬时起火，从横向和纵向开始蔓延，火势很大，火焰高度以极快速度达到计算空间顶端，其后火焰即迅速消失，退回到稳定燃烧的状态，整个爆发火焰的持续时间不超过2 s。设定的油滴注入高度只有2.5 m，但是火焰高度在1 s内轻松达到10 m，这说明油滴蒸发产生的油蒸汽扩散速度非常快，迅即弥漫于广大空间。起初的火焰发展按着油滴分布展开，但部分油滴快速蒸发产生油蒸汽，随后在热浮力推动下在油池上方大面积扩散，与空气接触面积迅速增加，促进了火焰广度和高度的发展，形成爆发式燃烧。图 5给出了中心线上若干观测点的气流速度随时间变化历程。火焰稳定燃烧阶段，各观测点气流平均最高速度只能达到5.5 m/s，而爆发火焰阶段，气流速度超过20 m/s。

图 6给出了中心线上若干观测点的温度随时间变化历程。火焰生成之后，在稳定燃烧阶段，最高温度发生在距离油池表面0.5 m处，平均温度超过950 ℃，随着高度增加，平均温度逐渐降低，在5 m高度，平均温度只有80 ℃上下。油滴注入之后，火焰爆发，整个中心线上的温度都急剧升高，在短时内8m以下全部达到800 ℃以上的高温，8 m至10 m达到了400 ℃以上并快速回落，但维持在较高温度的时间都会持续几秒钟，例如5 m高度，超过800 ℃的高温都会持续接近2 s，这足以引燃一般常见的可燃物。火焰构图显示了火焰发展的形状变化，温度分布反映了不同高度燃烧反应的剧烈程度，如此一来，爆发火焰的发展过程显得更加具体。

图 7给出了距油池中心几个横向距离若干高度监测点的热辐射通量。火焰稳定燃烧时，距离中心0.5 m的轴线上0.5 m高度以下热辐射通量在2.5 kW/m²上下浮动，其后随着高度的增加而减小，在10 m高度热辐射通量仅有不到0.006 kW/m²，这在短时间内都不会对人员和设备造成损伤。爆发燃烧后，5 m以上热辐射通量的增量都超过了200 kW/m²，而5 m处热辐射通量则增加了230 kW/m²，增量达到最大，随着高度的增加，热辐射通量有所减小，但是在10 m处增量仍达到了220 kW/m²。随着离油池中心横向距离的增加，每个高度上的热辐射通量在量值上逐渐递减，但是爆发火焰之后的极速增长都是可观的，如此大范围高强度的热辐射是造成灾害损失的主要原因。

图 8和图 9给出了油滴质量系数为0.6和1.2的火焰形状图，显然参数影响是存在的，但影响程度有限。油滴质量系数的增加，可以产生周向规模更大的火焰，火焰持续时间也略长，但是火焰能够达到的最大高度则几乎没有变化。图 10和图 11给出了油滴质量系数为0.6和1.2时中心线上若干观测点的温度随时间变化历程，大的油滴质量系数产生的最高温度略高一些。

5 结论

本文对水倾倒入燃烧油中发生爆发燃烧的现象进行了深入分析，详细讨论了这一现象发生的物理化学过程，并对其中尚不能实施可信数值模拟的过程构造了数学模型，形成了模拟这一现象的数值模拟手段。数值算例测试表明，喷溅油滴蒸发产生的油蒸汽，在热气流的带动下，快速弥漫于广大空间，在极短的时间内就可以被携带至远大于喷溅油滴所能到达的位置，因而造成火焰规模的急剧扩大，而爆发火焰的辐射热通量则急剧增长，产生破坏性效果。通过本文设计的数值模拟手段，能够刻画真实情况发生时的爆发式火焰增长，但数学模型需要进一步完善，因为水不同的状态变化决定每个阶段喷溅油滴的质量，因此喷溅油滴质量系数不能完全具体反映每个阶段物理现象，接下来会细化喷溅油滴质量相关系数，将每个阶段水量损失和喷溅油滴质量联系起来，反映喷溅强度的变化，在理论上建立更具体的数学关系。此外，通过数值模拟结果所反应的细节来设计实验方案，通过实验手段，验证并改进模型，以得到与真实情况更相符的数值计算结果。