油气爆炸是发生在化工领域和储油场所最严重的灾害之一^[1-2]。在实际的储油场所，比如油料洞库、覆土油罐巷道中不可避免地存在各种障碍物，如通风管道、各种阀门以及通道面积突然缩小的通道口(比如通道防护门)等。当火焰在传播过程中与此类障碍物相互作用时，火焰形态、火焰传播速度、流场结构及爆炸超压等与没有障碍物的空间相比都会发生显著的变化^[3-5]。因此，研究火焰与障碍物的相互作用机制及由此引发的流场结构变化对于有效预防油气爆炸事故的发生及降低灾害损失意义重大。

对置障管道内甲烷、氢气和丙烷等可燃气体的爆炸特性已有大量的研究^{[4, 6-13]}，研究重点主要为障碍物对火焰的湍流加速机理以及爆炸流场特性等。但受限于现有实验方法和测试手段，很难捕捉到详实的火焰传播过程和湍流流场。特别是在数值模拟方面，对于可燃气体爆炸过程的研究大都采用传统的雷诺平均模拟(RANS)方法^[14-16]，其主要缺点在于它只能提供湍流的平均信息，这很难满足研究置障管道内油气爆炸的复杂湍流过程的需求，无法有效预测油气爆炸火焰-湍流耦合的行为。

大涡模拟将耗散尺度的脉动进行过滤，只对大尺度脉动进行求解，其空间分辨率比雷诺平均更高，可以有效捕捉油气爆炸过程中的湍流特征。目前，关于置障管道内油气爆炸的大涡模拟研究偏少，油气爆炸火焰与湍流的耦合特性还亟待进一步研究。本文中针对连续圆孔障碍物条件下油气爆炸火焰与湍流耦合过程进行大涡模拟，并将数值模拟得到的火焰形态、火焰传播速度和爆炸超压等特征参数与实验结果进行比较，进而揭示油气爆炸过程中火焰-湍流耦合机理。

1.   数学模型

1.1   大涡模拟控制方程

由于在置障条件下的油气爆炸是一个强湍流爆燃过程，因此本文中采用较能捕捉湍流特征的壁面自适应局部涡粘模型(WALE)^[17]的大涡模拟(LES)对实验过程进行重现。大涡模拟对小尺度的脉动进行过滤，而只对大尺度的湍流脉动通过N-S方程进行计算求解，而对小尺度脉动产生的影响则通过亚网格模型进行模型假设来模拟，经过过滤后大涡模拟的控制方程为:

连续性方程:

动量方程:

能量方程:

反应进程变量c的输运方程:

式中: ρ为密度，p为压力，t为时间，u_i、u_j为速度分量，h_s为显焓，λ为热导率，D为扩散系数，${\dot \omega _{\rm{c}}}$为归一化的化学反应速率，根据${\dot \omega _{\rm{c}}} = \frac{{{S_{{\rm{h, chem}}}}}}{{{H_{{\rm{comb}}}}{Y_{f{\rm{uel}}}}}}$计算得到，其中S_{h, chem}为化学反应热，H_comb为1 kg燃料燃烧产生的热量，Y_fuel为混合物中燃料的质量分数，横线上标为物理空间过滤量，波浪线上标为质量权重过滤量。亚网格热焓通量通过梯度近似设为：$\rho \left( \widetilde{{{u}_{i}}{{h}_{\text{s}}}}-\text{ }{{{\tilde{u}}}_{i}}{{{\tilde{h}}}_{\text{s}}} \right)=-\frac{{{\mu }_{\text{sgs}}}{{c}_{p}}}{P{{r}_{\ \text{sgs}}}}\frac{\partial \tilde{T}}{\partial {{x}_{i}}}$，其中μ_sgs为亚网格黏度，Pr_sgs为亚网格普朗特数，T为温度，c_p为比定压热容。

建立亚网格模型使控制方程封闭, WALE模型^[17]通过大涡速度场动态地求出亚网格模型的系数，能够较为精确地捕捉层流到湍流的转变，且不需要显式过滤^[18]，计算量较小。因此本文采用WALE模型作为计算的亚网格模型。其中涡粘模型方程为:

模型中

式中：S_ij为其旋率张量，L_s为亚网格尺度混合长度。

1.2   亚网格燃烧模型

在湍流预混燃烧的模拟中，火焰的厚度一般较小，若直接计算，对网格的要求一般会比较高，因此本文中采用的燃烧模型为Zimont模型^[19]，在该模型中，通过加厚火焰前锋来计算湍流火焰速度：

式中：A为模型常数，取值为0.5^[20]；l_t为湍流长度尺度，根据l_t=C_SΔ进行计算，其中C_S是Smagorinsky常数, Δ为单元特征长度；u′为亚网格脉动速度，根据u′=l_tτ_sgs^-1进行计算，其中τ_sgs^-1是亚网格混合率，由公式$\tau _{\text{sgs}}^{\text{-1}}=\sqrt{2{{{\tilde{s}}}_{ij}}{{{\tilde{s}}}_{ij}}}$计算得到，${{{\tilde s}_{ij}}}$是应变率张量，由${{\tilde s}_{ij}} = \frac{1}{2}\left( {\frac{{\partial {{\tilde u}_i}}}{{\partial {{\tilde x}_j}}} + \frac{{\partial {{\tilde u}_j}}}{{\partial {{\tilde x}_i}}}} \right)$计算得到；U_l为层流火焰速度，取0.43 m/s；α=k/ρc_p为未燃混合物的摩尔传热系数。

2.   数值模拟

2.1   物理模型

本文中大涡模拟以文献[11]中实验结果为基础进行，实验装置和实验系统如图 1所示。采用的物理模型和实验管道尺寸相同，如图 2所示，管道内部尺寸为100 mm×100 mm×1 000 mm，实验管道材质为有机玻璃，厚度为20 mm。在管道内安装了4个连续圆孔障碍物，第一个障碍物距离点火端100 mm，每个障碍物之间的间距为100 mm，圆孔半径为40 mm，障碍物阻塞率为49.8%。管道右侧端部采用盲板密封，盲板中心安装有点火头，距离点火头20 mm位置安装有一支压力传感器。左侧端部采用薄膜密封，薄膜的破口压力很小，对流场影响可以忽略^[6]。采用实验室自主研制的油气雾化装置进行配气，实验中初始压力和温度分别为实验室环境压力(约96.4~97.3 kPa)和温度(约298~303 K)。采用初始体积分数φ=1.7%的油气进行实验，并用电火花引爆油气。利用压力传感器记录管道内超压-时间曲线，高速摄影仪拍摄油气爆炸火焰发生和发展过程。为保证实验的可靠性，每组实验在相同工况下重复进行至少3次。

图  1  实验系统图

Figure  1.  Schematic diagram of experimental system

下载: 全尺寸图片 幻灯片

图  2  实验管道示意图

Figure  2.  Schematic diagram of the explosion pipe

下载: 全尺寸图片 幻灯片

在数值计算时，为了减小管道出口气体回流对管道内部压力的影响，物理模型在管道外沿x、y、z方向分别延伸500、500、1 000 mm，即设置了一个尺寸为500 mm×500 mm×1 000 mm的外部空间区域。这个外部空间区域和管道出口相连接，可以模拟管道内火焰冲出管道的反应进程，也可以模拟管道泄压过程。

2.2   网格划分

在大涡模拟中，亚网格尺度的物理扩散是随着网格的加密而减小的，因此，对于大涡模拟而言，并不存在通常数值计算中所谓的“网格无关性”概念，较好的做法是采用十分精细的网格和很小的时间步长，以确保其数值扩散远小于雷诺平均方法。但是，在实际计算中，由于受到计算机计算能力的限制，通常仍然会进行网格无关性检验，兼顾计算精度和计算成本。鉴于此，本文中在网格无关性检验的基础上，采用六面体网格对计算模型进行网格划分。整个计算区域网格总数为1 200 000个，其中管道内部大约有950 000个网格，网格尺寸大小为2~2.5 mm，外部空间网格总数为250 000左右，网格尺寸为2~10 mm。

2.3   边界条件和初始条件

管道右侧封闭端、管道内壁面和障碍物表面都设置成绝热无滑移壁面边界。另外，由于实验中管道左侧薄膜能够在很低的压力条件下破裂，对油气泄压爆炸过程影响较小，因此在大涡模拟计算时，忽略薄膜的影响^[7]。管道外部扩展区域的壁面设置为压力出口，压力出口表压设置为0 Pa。

未燃气体和已燃气体的比热容按照近似温度的五阶分段多项式函数来计算。气体分子黏度通过Sutherland法则来计算。初始反应物选用φ(CH)=1.7%的汽油蒸汽^[1-2]，层流火焰速度假设为随温度和压力变化保持常数不变，为0.43 m/s^[21]。

初始条件包括压力、流场初速度、初始温度和反应进程变量。初始温度设为300 K，其余初始参数设置为0。采用Patch功能实现点火引爆，具体操作方法为在管道右侧封闭端中心位置设置一个半径为5 mm的半球形区域，并且设置该区域的反应进程为1，以此来实现模拟点火功能^[7]。

2.4   计算方法

大涡模拟计算通过ANSYS公司的Fluent计算平台来实现。采用SIMPLE算法对压力和速度进行耦合计算，对流项和扩散项分别采用二阶迎风格式和二阶中心差分格式来进行离散。采用联想ThinkServer TD350服务器(64位、Xeon E5-2609 v4 CPU (8 processes)、8G RAM)实现模拟计算。迭代2计算步长设置为1×10^-5 s，每一个时间步长要求迭代20次，保证能量方程残差小于1×10^-6，反应进程变量方程残差小于1×10^-4，动量方程残差小于1×10^-5。完成一次模拟计算的时间大约为52 h。

3.   模型验证

以往大多数模拟研究，采用将模拟计算结果和实验结果进行对比来说明数值模拟模型的可靠性^[22-23]。本文中采用对比数值计算和实验所得火焰传播形态、压力-时序曲线、火焰锋面位置和火焰传播速度来判断数值模拟计算结果的精度。为了证明大涡模拟的精确性，在数值模拟计算时还采用了RNG k -ε湍流模型来进行计算，并将计算结果和大涡模拟计算结果进行对比，相关的湍流控制方程^[14]在此不做赘述。

图 3为选取的模拟计算和实验结果中8个不同时刻火焰传播形态对比图。从模拟结果可见，大涡模拟和RNG k-ε湍流模型都能比较直观地反映出火焰在管道内外的传播形态变化过程。但是，大涡模拟相比RNG k -ε湍流模型能更好地表现出火焰结构的褶皱、弯曲等精细化演变过程。图 4所示是数值模拟和实验所得火焰锋面位置随时间的变化关系曲线，从图 4中可见，大涡模拟计算结果和实验所得火焰锋面位置基本一致，体现出大涡模拟对火焰锋面位置较好的预测性。但是，采用RNG k -ε湍流模型计算所得的火焰锋面位置与实验结果相比较误差较大，尤其是预测的火焰锋面的移动速度明显偏高。

图  3  实验和模拟所得火焰结构对比

Figure  3.  Comparison between experimental and simulated flame structures at different times after ignition

下载: 全尺寸图片 幻灯片

图  4  实验和模拟火焰锋面位置对比

Figure  4.  Comparison between experi-mental and simulated flame front locations along pipe versus time

下载: 全尺寸图片 幻灯片

图 5为数值模拟和实验相同测点监测所得压力-时间曲线对比关系。从图中可见，实验、大涡模拟和RNG k -ε湍流模型计算结果中取得最大超压峰值的时间分别为27.4、25.5、16 ms，后两者与实验结果相比较误差分别为6.93%和41.6%。三者最大爆炸超压峰值分别为194.77、165.35、83 kPa，后两者与实验结果相比较误差分别误差为15.1%和57.4%。由此可见，大涡模拟对爆炸超压的预测较RNG k -ε湍流模型精度更高。值得注意的是，本文中实验和数值模拟得到的最大爆炸超压峰值相比其余大多数文献中获得的最大超压峰值更大，原因主要为：一是本文中采用的实验管道的长径比相较其余大多数文献中的管道长径比更大，因此油气爆炸火焰在管道内加速距离更长，导致爆炸发展更充分，引起更剧烈的爆炸；二是本文中采用了4个连续障碍物，显著增强了流场湍流度，提高了已燃气体和未燃气体的热交换速率。由此，增强了障碍物对火焰的加速效应，并能诱导更高的火焰传播速度，根据文献[11]中关于油气爆炸火焰传播速度和爆炸超压演变过程的耦合关系可知，最大爆炸超压峰值会变得更大。

图  5  实验和模拟超压-时间曲线对比

Figure  5.  Comparison between experimental and simulated overpressure time histories of monitor point at the closed end

下载: 全尺寸图片 幻灯片

图 6所示是实验和数值模拟所得火焰传播速度对比，火焰传播速度是通过计算火焰锋面在管道轴向移动速度得到，计算公式为：

图  6  实验和模拟所得火焰传播速度

Figure  6.  Comparison between experi-mental and simulated flame speeds

下载: 全尺寸图片 幻灯片

式中：S_f表示火焰传播速度, x_n+1-x_n表示相邻的两幅高速摄影照片中(在大涡模拟中指反应进程变量c=0.5的相邻两幅图像)火焰锋面的真实距离之差, Δt_n表示选取的两幅相邻图像的时间差(本文中Δt_n=0.001 s)。

从图 6中可见，实验和大涡模拟所得火焰传播速度的变化规律基本一致，在火焰锋面没有接触障碍物之前，火焰保持较低速度传播，大致为9 m/s；当火焰锋面经过障碍物之后，火焰速度开始急剧增大，并且在27 ms取得最大值。实验和大涡模拟所得最大火焰速度分别为410.71、390 m/s，大涡模拟和实验结果相对误差为5%。当火焰速度达到最大值之后，由于可燃气体量的减少和热量的大量损耗，火焰传播速度急剧下降。对于RNG k -ε湍流模型计算结果而言，其火焰传播速度在变化趋势上与实验结果和大涡模拟计算结果相比误差较大，不能体现出障碍物对火焰的加速过程，并且最大火焰速度仅为140 m/s，相比实验结果误差为65.9%。

由上述对比分析可见，本文中采用的大涡模拟计算模型能够较好地预测半受限空间汽油-空气混合物泄压爆炸过程火焰和爆炸超压特征。

4.   置障管道油气泄压爆炸火焰特性分析

4.1   障碍物对火焰传播形态影响

图 7所示为大涡模拟所得7个不同时刻火焰传播经过障碍物过程中的3维立体火焰图像，火焰结构采用反应进程变量为c=0.5的等值面渲染得到。从图 7可见，t=14 ms时，火焰还没受到管道侧壁和障碍物的影响，火焰锋面保持半球形向前传播；t=18 ms时，火焰锋面传播经过第1个障碍物，火焰形状由半球形转变成蘑菇形，并且火焰锋面面积显著增大；t=20~23 ms之间，火焰穿越剩余3个连续障碍物，在这一过程中，火焰锋面形态变得更加不规则，出现明显的褶皱和卷曲变形现象，并且当火焰锋面受到障碍物壁面阻碍后，还出现火焰向点火端回传现象(如图 7中20、21 ms火焰图像所示)。当火焰锋面穿越4个连续障碍物之后，火焰出现明显的破碎现象，如t=25，26 ms所示。此外，还可以观察到火焰锋面在传播经过障碍物时，在障碍物和管道内壁形成的腔体内还有一部分未燃气体滞留，文献[23-24]中的研究表明这一部分未燃气体与爆炸超压后期振荡有密切的关系。

图  7  火焰结构三维图像(反应进程等值面c=0.5)

Figure  7.  Three-dimensional flame structure images (iso-surface of progress variable c=0.5) at different times after ignition

下载: 全尺寸图片 幻灯片

4.2   火焰和流场耦合关系分析

4.2.1   内场火焰和流场耦合关系

图 8为管道内部置障条件下油气混合物预混火焰传播和流场结构图，图中绿色代表火焰锋面(c=0.5)，彩色箭头代表速度矢量。在t=14 ms之前，管道内气体流动是层流状态，流场强度较小，流线形态较规则，流场方向基本保持从点火端向管道开口端发展；t=18 ms时，火焰传播经过第1个障碍物，流场形态开始变得不规则，并且在相邻障碍物之间的流场区域出现涡旋结构，在涡旋结构作用下，火焰锋面形状从半球形向蘑菇状转变，增大了已燃气体和未燃气体的接触面积；t=20 ms时，火焰传播经过第2个障碍物，此时流场强度进一步增大，涡旋结构更加明显，并且在涡旋结构的卷曲作用下，火焰在第1~2个障碍物之间出现向点火端回传现象，同时可见此时火焰锋面前部的流场强度更加剧烈；t=21~23 ms之间，火焰传播经过第3~4个障碍物，在这个过程中，火焰结构变得更加褶皱、弯曲，增大了已燃气体与未燃气体的接触面积，提高了燃烧速率，增大了火焰传播速度，因此导致火焰锋面前方流场涡旋更加明显，流场强度急剧增大。同时，随着流场强度的增大，涡旋更加明显，对气体的卷曲作用明显增强，因此提高了已燃气体和未燃气体之间的交换速率，反过来提高燃烧速率和火焰传播速度。因此，在火焰传播和流场上述正激励作用的诱导下，火焰形态产生严重变形，火焰锋面前流场强度急剧增大。

图  8  管道内火焰传播和流场结构图

Figure  8.  Flame propagation and flow field in an obstructed pipe

下载: 全尺寸图片 幻灯片

4.2.2   外场火焰和流场耦合关系

图 9所示为外场火焰和流场结构图。t=26 ms以前，火焰还没传播至管道开口端，但是从t=21 ms起，管道外流场逐渐形成“双涡旋”结构，且随着时间的推移，涡旋结构越来越显著，并且以管道中轴线为分界线，中轴线上部涡旋方向为顺时针，中轴线下部涡旋方向为逆时针方向。t=26.5 ms时，火焰以射流火焰形式从管道内部传播到管道外部，并进入管道外部涡旋区域，在涡旋作用下，火焰锋面在t=27 ms时从柱状火焰向蘑菇状火焰转变，并于t=28 ms时形成形状较完整的蘑菇状火焰。通过图中流场矢量图颜色可以判断得到，在蘑菇状火焰区域，轴向流场速度明显大于径向流场速度，可以推断，蘑菇状火焰的形成和流场在轴向和径向的速度差存在密切的内在联系。此后，火焰锋面的蘑菇状火焰逐渐向哑铃状火焰转变(t=30 ms所示)，并且火焰锋面向前传播的同时，涡旋结构也向前移动，并且始终和火焰锋面重合。

图  9  外场火焰和流场耦合图

Figure  9.  Coupling results of the flame propagation and flow field outside the pipe

下载: 全尺寸图片 幻灯片

4.3   爆炸超压与火焰传播耦合关系

爆炸超压是油气爆炸过程中一个重要的爆炸参数，为了研究其演变规律与火焰传播之间的耦合关系，根据爆炸超压曲线的变化特点将其划分为图 10(a)中所示的4个变化阶段。在分析时，将大涡模拟计算得到的火焰传播速度、火焰面积与爆炸超压进行联合分析，对应的曲线如图 10所示。火焰面积根据下式计算：

图  10  爆炸超压与火焰传播速度、火焰面积耦合关系

Figure  10.  Coupling relationship between explosion overpressure and flame speed and flame surface area

下载: 全尺寸图片 幻灯片

式中：A_i代表等值面上微单元的面积；n代表等值面上微单元的数量(注：计算时等值面取c=0.5)。

由图 10可见，在15 ms以前(阶段Ⅰ)，爆炸超压平缓上升，对应的火焰传播速度和火焰面积变化幅度也很小，此时火焰还没受到障碍物的扰动，火焰处于层流传播阶段。从15~25 ms(阶段Ⅱ)，爆炸超压急剧上升，并在大约25 ms处取得最大值。在这一阶段，火焰穿越障碍物，火焰发生卷曲变形，层流火焰向湍流火焰转变，导致火焰传播速度和火焰面积急剧增大(如图 10(b)~(c)所示)，提高燃烧速率和热释放率，进而增大爆炸超压。在第Ⅲ阶段(大约25~30 ms)，爆炸超压和火焰传播速度都急剧下降，并且都在大约30 ms左右取得极小值，同时火焰面积也小幅下降。在前三个阶段，爆炸超压、火焰传播速度和火焰面积具有“同升同降”的相似变化趋势，耦合关系较好。在第Ⅳ阶段，火焰锋面传播到管道最左端，火焰逐渐熄灭，火焰面积在短时间内保持上升趋势，随后急剧下降。随着火焰面积的改变，爆炸超压呈现出振荡下降的趋势，并形成多个超压峰值，其中振荡期的第一个超压峰值对应火焰面积的最大峰值。由上述分析可见，油气在置障管道内泄压爆炸过程的爆炸超压、火焰传播速度和火焰面积之间耦合关系明显，三者的变化趋势存在较好的一致性。

5.   结论

(1) 通过比较大涡模拟和RNG k -ε湍流模型计算结果发现，大涡模拟在预测油气爆炸超压、火焰传播速度及火焰形态变化等方面比RNG k -ε湍流模型精确度更高，且能表现出更多流场的精细化结构；

(2) 障碍物诱导管道内形成湍流度较高的流场区域，导致火焰产生褶皱弯曲变形，增大火焰面积，加速火焰传播；

(3) 通过比较分析爆炸超压、火焰传播速度和火焰面积的演变规律，发现三者内在联系密切，具有显著的耦合性，且随时间的变化趋势存在高度的一致性。