基于简化反应机理的甲烷-空气爆炸增厚火焰大涡模拟

2021-11-26郑凯蒋军成邢志祥吴凡吴洁余明高

化工学报 2021年11期

郑凯，蒋军成，邢志祥，吴凡，吴洁，余明高

（1 常州大学环境与安全工程学院，江苏常州213164； 2 重庆大学煤矿灾害动力学与控制国家重点实验室，重庆400044）

引言

甲烷作为天然气和煤矿瓦斯的主要成分，其在生产储运过程中一旦发生泄漏或积聚形成爆炸浓度范围内的可燃气体时，遇火源可发展为爆燃、爆轰以及从爆燃向爆轰的转变等爆炸事故[1-3]。一般而言，管道中甲烷-空气爆炸传播是涉及湍流流动的化学反应过程。随着高性能计算技术的发展，大涡模拟(large eddy simulation, LES)已经成为一种预测管道内甲烷-空气爆炸湍流反应的重要手段[4-6]。相对于大涡模拟，雷诺平均模拟(Reynolds average Navier-Stokes，RANS)由于对Navier-Stokes (N-S)方程进行了时间平均，计算准确性较差，难以精确预测爆炸火焰结构。然而由于计算能力的限制，当前大涡模拟中计算网格尺寸要远大于待求解火焰厚度。为了解决这一问题，学者们提出了三种亚网格模型：火焰表面密度模型(flame surface density，FSD)[7]，增厚火焰模型(thickening flame，TF)[8-9]和G 方程模型[10-11]。

在TF模型中，通过提高燃烧组分扩散系数并同比例降低反应速率，可以使层流火焰速度保持不变但火焰厚度增大，以保证火焰在大涡尺寸网格上求解。由于火焰厚度增加导致的亚网格火焰表面的丢失可用亚网格尺度褶皱因子来评估，其为火焰总表面与可求解的表面积的比值[12]。为了避免增厚过程改变组分的输运参数，采用火焰探测函数对火焰前锋进行动态增厚[13]。当前基于增厚火焰模型的大涡模拟技术已经取得长足进展，且在预测预混和非预混火焰的可行性方面开展了大量研究。Selle等[14]采用简化的甲烷-空气两步燃烧机理和增厚火焰模型开展了复杂几何结构中反应流的大涡模拟，表明增厚火焰模型能够预测复杂燃烧装置的反应特性。Wang 等[15]基于动态亚网格模型参数计算褶皱因子，结合动态增厚火焰开展大涡模拟分析了湍流本生灯火焰。张科等[16]采用基于释热率构造火焰探测函数的动态增厚火焰模型，对甲烷-空气同轴射流燃烧器的非预混燃烧进行可压缩大涡模拟。Guo 等[17]基于改进的动态火焰探测函数，开展了耦合多步反应机理的增厚火焰大涡模拟，分析了增厚因子对火焰结构的影响。

对于管道内甲烷-空气爆炸的大涡模拟，当前应用较多的为FSD 模型，研究结果表明其能够准确地捕捉爆炸预混火焰结构，精确预测超压变化和火焰传播流场特征[18-22]。但关于增厚火焰模型大涡模拟的当前报道相对较少，且相关研究主要侧重于分析障碍物引起的湍流对爆炸火焰加速过程的影响。Yu 等[23]将增厚火焰模型应用到可压缩燃烧预测并利用其捕捉了爆燃向爆轰转变过程。Emami 等[24]利用增厚火焰模型和27步反应机理，通过求解二维过滤后的Navier-Stokes 方程，研究了障碍物管道中预混氢气-空气的爆燃向爆轰转变过程。Volpiani等[25]利用动态火焰褶皱函数对爆燃火焰通过连续障碍物的传播过程进行了大涡模拟，结果表明大涡模拟能够预测爆炸超压等特征参数。

本文拟采用单步不可逆甲烷-空气简化燃烧机理开展三维大涡模拟，研究对象为出现郁金香结构[26]且基本以层流形式传播的管道内甲烷-空气爆炸火焰。模拟中采用Légier[13](TF1)与Durand[27](TF2)提出的两种火焰探测函数对火焰面进行动态增厚。基于计算得到的火焰三维等值面、二维切片分析火焰动态传播特征，将火焰传播速度、超压曲线与前期实验结果进行对比以验证模型可靠性。基于有效增厚因子分析火焰探测函数对大涡模拟结果的影响。

1 数值模型与计算步骤

1.1 增厚火焰模型

对于以预混火焰方式传播的管道内甲烷-空气爆炸，可以通过增加扩散速率与降低反应速率的方式实现火焰增厚。增厚后的组分输运方程如下：

式中，ρ为密度；u为速度矢量；Yk是组分k的质量分数；Dk为组分k的扩散速率；F为增厚因子；(Q)为组分k的反应速率，可以用Arrhenius 公式进行描述，其中Q表示用于计算反应速率的物理量。

考虑到前锋褶皱对火焰表面积的影响，将火焰前锋的扩散系数和反应速率乘以效率因子Ξ。本文中采用Zimont 提出的湍流火焰封闭模型(turbulent flame speed closure model，TFC)对效率因子进行评估[28]，该模型可改写成如下形式：

式中，Ut是湍流火焰速度；A是常数，对于甲烷-空气建议值为0.52；u'是亚网格脉动速度；Ul是层流火焰速度；δ是层流火焰厚度；lt为湍流长度尺度，可用下式计算得到：

式中，Cs是Smagorinsky常数，可用Smagorinsky-Lilly涡粘模型计算得到[29]；Δ是网格尺寸。

u'的计算公式为：

为应变率张量:

效率因子Ξ，为亚网格湍流速度在FΔ和Δ尺寸上的比值：

LES过滤后的组分输运方程为：

其中,-表示过滤值；～代表平均值。

1.2 火焰探测函数

为了避免非反应区域由于扩散系数的增加而导致的组分混合与传热问题，只在火焰前锋内进行动态增厚，火焰前锋的动态增厚系数由火焰探测函数Ω进行控制。理论上，在反应区域内Ω的值为1，在区域外则为0。本文中采用两种火焰探测函数(TF1和TF2)。

TF1为[13]:

TF2为[27]:

其中，c为反应进程变量：

基于火焰探测函数TF1 和TF2，可计算得到有效增厚因子：

1.3 反应机理

甲烷-空气爆炸的预混火焰反应过程用单步不可逆燃烧反应机理进行描述[15]：

反应速率可用Arrhenius公式进行描述：

式中，A是指前因子；Ea是活化能；ρ为密度；nk是组分k的反应指数；Mk是组分k的摩尔质量；R为理想气体常数；T为温度。反应机理的速率常数如表1中所示。根据该反应机理计算得到的层流火焰速度和绝热火焰温度分别为0.40 m/s和2328 K。

表1 反应机理的速率常数Table 1 Rate constants for the one step mechanism

1.4 参数设置

模拟在尺寸为1000 mm×100 mm×100 mm 的三维管道中进行。控制方程采用有限体积法离散，压力-速度场采用SIMPLE 算法耦合。采用二阶上迎风格式与中心差分格式离散对流项与扩散项。整个计算区域内初始网格尺寸为2.5 mm。为了确保最大限度求解湍动能，采用自适应网格对火焰前锋附近区域网格进行动态加密。加密过程是基于温度梯度来进行的，其中网格粗化与细化的阈值分别为30 K/m 与100 K/m。每个网格最多加密2 次，加密之后最小尺寸为0.625 mm。增厚因子Fmax的值为10。计算对象为当量比为1的甲烷-空气混合气体。各组分的比热和黏度用五次多项式和Sutherland 公式进行计算[30-31]。组分输运系数与热导率分别采用Chapman-Enskog 公式与动力学理论计算得到[17]。混合反应物的热物理特性由理想气体混合法则计算得到。初始温度和压力分别设定为295 K 和0.1 MPa。管道壁面采用无反射固体边界条件。压力监测点位于管道点火端中心处。在未燃气体中设置一个半径为5 mm 的半球形区域，假设该区域内甲烷/空气完全反应，区域内温度设置为2100 K。Wen等[4]的研究表明将点火半径设置为5 mm 时火焰锋面位置的大涡模拟与实验结果吻合度较高。管道中心轴线设为X轴(单位：m)，半球形点火区域球心坐标为(0，0，0)，即火焰沿X轴正向传播。时间步长设置为小于1×10-5s，并随着火焰速度增加逐渐减小。每个时间步长至少需要60 次迭代。动量方程的残差小于1×10-5，能量和物质输运方程的残差分别小于1×10-6。

2 结果与讨论

2.1 火焰传播结构与流场

对于本文大涡模拟所采用管道中传播的甲烷-空气爆炸预混火焰，前期实验结果表明其传播会经历经典“郁金香”火焰形成的四个阶段，即：半球形、指形、火焰接触管道壁面到形成平板形状、郁金香形状[32]。图1 分别给出了TF1 和TF2 时郁金香火焰形成过程中温度二维切片，代表了火焰传播前锋结构演变过程。从图1 中可看出，无论是对于TF1 还是TF2，大涡模拟均能够完整再现郁金香火焰形成过程。点火完成后，在11 ms 时可观测到火焰以半球形传播。此时，TF1 和TF2 火焰传播结构和距离相差不大。此后，由于管道壁面的限制，半球形火焰会转变成指形。在28 ms时，TF1 和TF2 中均可观测到火焰以指形结构传播。但是相对于TF1 (0.11 m)，TF2 (0.13 m)中火焰前锋传播的要稍远些，代表了更大的火焰传播速度。随着火焰的进一步传播，在接触到管道壁面后，火焰前锋面积开始减少，在图1 中39 ms 可看到TF1 和TF2 中靠近管道壁面处火焰前锋逐渐消失，表明进入火焰接触管道壁面到形成平板形状阶段。此时，TF2 仍然出现了相对较大的传播距离。火焰前锋面积的进一步减少最终导致平板火焰前锋的出现。在实验中，平板火焰出现时间为44 ms，而模拟结果分别为47 ms (TF1)和45 ms(TF2)，实验和模拟结果基本吻合。此后，火焰前锋出现逆转，火焰靠近管道壁面处速度大于管道中心位置，火焰前锋出现凹陷，最终形成经典的“郁金香”火焰结构。图1 中86 ms 时火焰切片代表了TF1 和TF2 预测的郁金香火焰结构。“郁金香”火焰出现后会保持这一结构传播直到整个爆炸过程结束。

图1 郁金香火焰形成过程中不同时刻的温度二维切片Fig.1 2D-slices of the temperature field at different instants during the tulip shaped flame formation

图2 为TF1 中不同时刻大涡模拟预测的火焰前锋三维1200 K 温度等值面结构(左侧)与流场演化特征(右侧)。对于TF1 与TF2，除传播速度外，大涡模拟预测的火焰传播结构与流场基本一致。在流场结构中，两条曲线为火焰探测函数Ω=0.1 时的等值线，以划分已燃区、爆炸反应区和未燃区。在图2中11 ms 时刻，可以看到火焰以光滑的半球形传播，此时火焰处于自由膨胀阶段。因此，在流场中靠近火焰前锋流线由反应区呈扇形指向未燃区。在靠近管道壁面处流线方向发生改变，最终在远离火焰前锋的未燃区所有流线均平行于管道轴线。在28 ms火焰由半球形转变成指形后，受到管道壁面的限制，三维结构中火焰封面在靠近管道壁面时被压平。此时，已燃产物膨胀作用推动火焰前锋高速传播。对应流场中流线起源于靠近管道壁面处的火焰前锋，指向管道轴线方向并流向未燃气体区域。火焰前锋接触到管道壁面后，在39 ms 可以看出火焰在接触到管道壁面处的锋面消失，进而前锋总面积大大减小。火焰前锋面积的减小意味着反应速率的降低，导致该阶段火焰尖端传播速度出现下降。此时，火焰前锋附近流线主要起源于火焰前锋与管道壁面的接触点。随着火焰前锋面积的逐步减少，平板火焰在47 ms 时刻出现。流场中在未燃气体区域流线保持正向流动，在已燃产物区出现两个涡旋，其旋转方向均为由管道壁面指向中心轴线。涡旋的旋转推动靠近壁面处气体向前但拖拽中心轴线处气体向产物流动，最终导致了“郁金香”火焰的出现。图中86 ms 为完全形成的经典“郁金香”火焰结构。在“郁金香”火焰形成后，在三维图中可看到火焰呈现四个锥形结构前锋。在流场中未燃气体保持正向流动，而已燃气体区域基本上保持逆向流动。

图2 TF1中不同时刻的温度三维等值面和二维切片流场结构Fig.2 3D iso-surface by temperature of 1200 K and flow field with TF1 at different instants

2.2 火焰传播速度与超压

对于管道中可燃气体爆炸事故，火焰传播速度和超压是衡量爆炸危险程度的重要参数。本文中大涡模拟可靠性可通过作者前期实验与本文模拟得到的火焰传播速度与超压曲线定量对比来进行验证[32]。

图3中给出了火焰传播速度随时间变化曲线的实验和大涡模拟对比结果。爆炸发生后，实验与大涡模拟中都可观测到火焰传播速度的指数形式增长。这一现象主要是由于火焰在半球形和指形阶段前锋面积的指数增长所引起[33]。火焰的加速传播过程因前锋接触到管道壁面而结束。图3 中，实验观测到的火焰传播速度峰值和相应时间分别为17.5 m/s 和35 ms，TF1 和TF2 计算得到的对应值分别为15.4 m/s、35 ms 和15.7 m/s、32 ms。显然，TF1和TF2 预测的火焰传播速度峰值基本一致，稍小于实验结果，但实验与大涡计算结果最大相差仅12%。而TF1 预测的到达火焰传播速度峰值时间与实验结果相同，TF2 相应值则稍大。在火焰接触到管道壁面后，锋面面积减少，火焰传播速度也逐渐降低，并在形成平板火焰时达到最小值。伴随着郁金香火焰的形成，火焰再次出现加速。最终伴随着火焰传播速度的逐渐减小，火焰前锋以郁金香形状传播直到结束。火焰传播过程中，实验和大涡模拟得到的火焰传播速度曲线基本吻合，但TF1 预测结果准确度相对较高。

图3 火焰传播速度随时间的变化Fig.3 The time curves of flame speed

图4对比了实验与大涡模拟得到的超压随时间变化曲线。相对于火焰传播速度，超压的实验与大涡模拟曲线在“郁金香”火焰阶段存在较大差异。在图4中，与火焰传播速度变化趋势相一致的是，爆炸超压在火焰传播初期也呈指数增长。这是因为火焰传播速度的指数增长表明燃烧反应速率呈指数增长，导致产物生成速率和放热速率也相应提高，最终造成了超压的指数上升。因此，实验中观测到的较大火焰速度导致超压相对于大涡模拟结果要上升的更快，且两者之间差值在火焰传播速度达到最大值时达到最大。在火焰锋面接触管道壁面后，超压曲线斜率出现下降，该过程持续到平板火焰的形成时刻。此后，在“郁金香”火焰形成后，超压基本上呈线性增长。TF1 和TF2 预测结果均与实验基本吻合，但后者增长要稍微快于前者。

图4 超压随时间的变化Fig.4 The time curves of pressure

2.3 有效增厚因子

在增厚火焰模型中，有效增厚因子是影响大涡计算结果的重要参数。图5 给出了39 ms 时TF1 和TF2的Ω、F、Ξ以及Ξ×F分布二维切片，代表了火焰传播过程中两种火焰探测函数对计算结果的影响效果。根据Guo 等[17]的观点，一个理想的Ω值应该在火焰前锋内为1，而在火焰前锋外则为0。从图5中可看出，TF1 和TF2 预测结果均能达到该标准。在靠近火焰前锋中心位置，Ω的值基本接近1。随着远离火焰前锋中心，Ω的值逐渐减小，并在已燃区和未燃区中减小到0。显然，TF1捕捉到的反应区域要明显大于TF2。TF1 预测的Ω值大多为1，只有在靠近已燃与未燃区域边界处少量值小于0.5，但总体上是从1突然下降到0的。而在TF2中，只有在预测区域中心位置处少量值为1，且从1 下降到0 的过程是渐进的。

在大涡模拟中，效率因子Ξ为在过滤网格上火焰总表面与已求解火焰表面积的比率，其值与火焰增厚导致的亚网格尺度火焰面积丢失有关[15]。根据流场分析结果可知，火焰在以指形结构传播过程中管道内基本上保持层流流动。所以在图5中整个管道大部分区域内Ξ均为1。只有在火焰前锋及靠近火焰前锋附近的已燃区域内Ξ的值大于1。在已燃区域内，Ξ主要是受到火焰前锋与管道壁面接触后导致该区域流线方向改变的影响，在该区域内Ξ的值较小(小于1.18)。而在火焰前锋尖端，Ξ的分布则主要受到火焰拉伸效应的影响，故其在该区域内Ξ的值较大(接近于1.5)。总体上，TF1 和TF2 两种工况中计算得到的Ξ分布基本保持一致，表明火焰捕捉函数对Ξ基本上没有影响。

由于TF1和TF2中Ξ的分布基本一致，可知Ξ×F的分布只受到F的影响。而根据有效增厚因子计算公式可知，当Fmax一定时，F是Ω的函数。因此，图5中F与Ξ×F的分布规律均与Ω相似，且TF1和TF2预测的F最大值均小于Fmax= 10。对于TF1和TF2，受到Ξ分布的影响，Ξ×F的最大值均分布在火焰前锋上。在反应区域内，TF1 中Ξ×F大部分数值要大于10，而在TF2 中只有靠近火焰前锋中心位置少部分区域的Ξ×F大于10。因此，总体上TF1对爆炸火焰的增厚效应要强于TF2。而增厚效应的增强，可能会导致火焰前锋燃烧反应求解的不充分，进而导致TF2 预测的火焰传播速度与超压上升均稍大于TF1。