湍流状态下双燃料混合物自燃和火焰传播过程

周 磊，刘宗宽，袁成威，张 帆，卫海桥

(天津大学 内燃机燃烧学国家重点实验室，天津 300350)

随着化石能源日益枯竭，燃用石油质液体燃料带来的环境问题日渐严峻，内燃机节能减排技术成为研究热点.为控制船舶发动机对海洋环境的污染，各国政府和国际组织制定的排放法规更加严格.国际海事组织(IMO)出台的TierⅢ排放法规于2016 年开始生效，规定船用发动机NOx排放较上一阶段标准降低80%.

天然气最主要的成分是甲烷，具有高氢碳比、燃烧清洁和储量丰富等优点，有望成为发动机新型清洁替代能源的潜力.研究表明，柴油引燃天然气船用双燃料发动机可有效降低二氧化碳和颗粒物的排放，并可在采用稀薄燃烧技术时大幅降低NOx排放[1]，天然气船机只需简单地增加排放后处理设备就能满足TierⅢ排放要求.然而由于天然气具有比热容大、着火极限高和火焰传播速度较慢等理化特性，以及双燃料混合物复杂的化学动力学原理，共同导致双燃料低速机低负荷着火及燃烧稳定性差.同时，低速机缸内伴有湍流且空间尺度大，湍流作用下双燃料混合气的着火及火焰传播过程更加复杂；此外，缸内气流运动会造成燃料混合物浓度分层和温度分层.因此，探究湍流对着火及火焰发展的影响有助于理解双燃料发动机缸内自燃和燃烧模式转变的过程.

Mastorakos 等[2]通过直接数值模拟研究湍流混合层中燃烧过程，发现湍流会缩短混合物滞燃期，而标量耗散率决定着火位置，活性最强且标量耗散率低的混合分数位置最容易发生着火.此外，Pitsch 等[3]模拟研究了正庚烷喷雾的着火过程，发现标量耗散率降低可提高燃烧速度.Jin 等[4]研究了湍流对三重火焰结构的影响，发现由于湍流的存在，沿着化学当量比混合分数对应的位置会形成多个锋面火焰，并在传播过程中接触连结.Wang 等[5]通过直接模拟手段研究双燃料层流混合层中着火到火焰稳定传播过程的特征时间以及特征时间与压力、混合层厚度、预混气温度和当量比等参数的拟合关系，发现滞燃期与温度、压力密切相关，峰值温度与混合层厚度以及预混当量比相关，稳态传播特征时间和混合层厚度以及层流火焰速度有关.

综上所述，针对湍流状态下双燃料混合物着火及火焰传播过程的研究还鲜见报道，尤其在低速机缸内高温、高压热力学状态下，对湍流流场中双燃料着火及燃烧机制的分析不够全面.基于此，笔者采用直接数值模拟(DNS)方法分析湍流流场对双燃料混合层着火和火焰发展过程的影响，探究双燃料着火特性、混合层中火焰传播结构以及火焰发展的特征.模拟研究中分别以甲烷和正庚烷作为天然气和柴油的替代物，相关结论可为双燃料低速船机高效、清洁燃烧技术的开发提供一定理论依据.

1 数值模型设置

1.1 数值计算方法

双燃料发动机内的湍流燃烧由连续方程、N-S 方程以及组分和能量的输运方程组成，笔者使用一种高阶方法数值求解控制方程[6-7].求解过程中，时间积分项采用二阶对称分裂算子方法，即在两个半步的扩散和对流时间步之间积分一个完整的化学反应速率时间步；时间步的推进采用二阶龙格-库塔(Runge-Kutta)结合二阶Adam-Bashforth(AB2)格式；空间差分采用六阶中心差分格式；组分和能量方程中的对流项采用五阶 WENO(weighted essentially non-oscillatory scheme)格式求解.数值方法的具体描述参见文献[7].

选用由Liu 等[8]提出的正庚烷简化机理，包含45组分、112 步反应.机理涵盖正庚烷自燃、高温化学以及甲烷氧化动力学，可准确预测高压下正庚烷滞燃期[8-9]、火焰传播速度和绝热火焰温度[10]等关键参数.此外，该机理广泛用于预混甲烷混合气下正庚烷喷雾、着火及燃烧的相关研究，与试验结果一致性良好[11].

1.2 计算区域及参数定义

图1 所示二维矩形计算区域尺寸为3.2 mm×3.2 mm，计算区域由2 0482个均匀网格划分，大约420 万个网格，网格解析度可达到1.56 μm[12].计算区域中正庚烷与甲烷/空气均匀混合气之间设置层流混合层，混合层中热分层和组分分层根据双曲正切函数进行设定，即

图1 二维湍流计算区域示意Fig.1 2D schematic of turbulence-like computational domain

式中：f (x)代表组分质量分数或温度随空间位置x变化的通用参数；fR和 fL分别为右侧值和左侧值；cx 代表中间位置；δ为初始混合层厚度.

根据双曲正切函数的特性，正庚烷质量分数为1%～90%的实际分布范围约为初始混合层厚度的3倍.计算区域初始条件设定见表1.依据合成湍流方法[13]添加湍流场，初始湍流场使用Passot 等[14]提出的各向同性能谱方程设置.平均流速设置为0，均方根流速(RMS) u＇ ＝0.5 m/s[15].积分长度尺度lλ＝0.8 mm，为区域长度的1/4，柯尔莫哥洛夫长度尺度lk＝22µm，积分尺度雷诺数为120.

表1 二维湍流混合层参数Tab.1 Main parameters for 2D turbulence-like mixing layer

需要指出，混合层厚度δ和混合区域厚度不是同一概念，混合区域大约为3δ厚度.计算区域各组分的空间分布见图2.计算区域左、右(垂直方向)为非反射的出流条件，上、下(水平方向)为周期性边界条件.

图2 初始状态下计算区域主要组分的空间分布Fig.2 Initial distribution of major species for computational domain

研究中所用的参数混合分数ξ根据氮气质量分数进行定义，公式为

2 模拟结果与数据分析

2.1 湍流场双燃料着火和火焰发展

运用CHEMKIN 软件中零维均匀反应器计算所得混合气滞燃期在混合分数空间的分布如图3 所示.第一阶段滞燃期极小值τ1,min=0.278 ms，对应混合分数为ξMR,1=0.240.第二阶段滞燃期极小值τ2,min=0.374 ms，对应的混合分数为ξMR,2=0.275.由初始当量比的分布可知，第一、二阶段滞燃期极小值均发生在浓混合物中(混合物总体当量比φ＞1).第二阶段着火位置对应温度低于负温度系数(NTC)温度区间的下限为850 K，这是因为初始温度继续升高，NTC 现象将导致计算区域混合物滞燃期变长，计算区域反应活性最强的混合分数对应ξMR,2位置，此处混合物浓度合适且不受NTC 效应影响，故自燃倾向最大、滞燃期最短.

图3 滞燃期、初始温度和当量比在混合分数空间上分布Fig.3 Distribution of ignition delay times(IDs)，initial temperature and equivalence ratio in mixture fraction ξ space

温度梯度的变化可表征化学反应空间位置和放热强度的变化，图4 展示了不同时刻下计算区域温度梯度的二维空间分布.计算区域局部温度和关键组分的变化如图5 所示.

图4a 中，t＝0.313 ms 时刻白色虚线代表混合物分数ξMR,1=0.240 等值线.图4b 所示t＝0.383 ms 时刻白色实线代表ξMR,2=0.275 等值线.第一阶段着火发生后的t＝0.313 ms 时刻、沿ξMR,1等值线左侧更浓混合物位置逐渐发生反应放热，形成多个连接的放热区域.结合图5a 中t＝0.313 ms 时刻线条A 上温度和组分分布特征，发现温度梯度峰值位置和低温化学表征组分过氧化氢酮KET(OC7OOH)峰值重合，且H2O2和HO2质量分数逐渐升高.文献[16]指出第一阶段着火后，H2O2逐渐积累增多，HO2在第一阶段和第二阶段着火后快速上升形成峰值.综合可知，此处火焰为正庚烷低温化学反应放热主导形成的冷焰.文献[17]研究混合层中低温化学反应时指出，冷焰逐渐形成于比ξMR,1更浓的混合分数位置，与笔者分析结论一致.

第二阶段着火(τig=0.372 ms)发生后，根据图4b中t＝0.383 ms 时刻的空间放热情况分布可知，沿着ξMR,2等值线右侧稀混合物中生成多个膨胀核心，正庚烷侧的浓混合物中形成薄条状的火焰结构.由图5b 中线条B 上的温度、温度梯度和组分分布可知，膨胀核心为着火后形成的高温膨胀核心.膨胀核心边缘放热形成两个温度梯度峰值，核心内部低温反应组分KET 快速消耗，高温反应区H2O2和HO2也在第二阶段着火后开始反应消耗.薄条状的放热区域KET 快速生成，定义KET 组分质量分数峰值超过20%区域为低温化学反应区域(绿色区域)，与窄条状区域基本重合，可知其本质是低温化学反应形成的冷焰区域.

随着反应进行至图4c 所示t＝0.403 ms 时刻，多个高温膨胀核心逐渐连结，形成高温的火焰传播前沿.正庚烷侧和甲烷侧高温火焰前沿区域为高温已燃区.在正庚烷侧(浓混合物侧)火焰下游仍存在低温化学形成的冷焰区，甲烷侧火焰前沿反应更强，温度也快速升高.该时刻的突出特点是出现同时存在浓预混冷焰和高温火焰的双火焰结构，如图5c 所示.

图5 标志线条上温度、温度梯度和关键组分质量分数分布Fig.5 Spatial distribution of temperature，temperature gradient and key species along lines

反应继续进行如图4d 和图4e 所示，正庚烷侧传播方向高温火焰前沿和低温冷焰之间的距离逐渐缩短，更浓混合物中的火焰前沿处当量比高、氧气含量少；而浓引燃油侧反应逐渐变缓，发生中温化学反应.相较于层流混合层，湍流使反应前沿变得更加褶皱，反应持续进行，甲烷侧火焰锋面褶皱程度低于正庚烷侧.

图4 不同时刻下温度梯度在计算区域的空间分布Fig.4 Distribution of temperature gradient for a sequence of time instants in the computational domain

第二阶段着火后在ξMR,2右侧生成多个高温膨胀核心，为探究湍流对自燃点和膨胀核心出现位置的影响，检测了着火时刻、膨胀核心在混合分数梯度场及ξMR,2等值线上混合分数梯度的分布特征.其中，混合分数梯度表示混合分数在物理空间上的梯度变化，即式中：ξ表示混合物分数.混合分数梯度是标量耗散率的组成部分，可用于解释湍流等因素导致的混合不均对热量和组分的耗散作用[2,18].图6 为湍流混合层中高温膨胀核心分布.其中，标记区域Z1～Z4用于现象解释，Z2和Z4为高温膨胀核心区域.

由图6a 可知，在ξMR,2等值线左侧，混合分数梯度变动很大，而右侧相对平缓，只在Z1和Z3两个区域存在大的梯度波动.而该区域并不存在高温膨胀核心，即该区域不适合着火和化学反应发展.Z2和Z4两个高温膨胀核心生成在当量比梯度平缓区域.分析可知，混合分数梯度大的区域，基团和热量快速耗散，不利于着火核心的形成，即使着火核心形成，热量也很快被耗散掉，无法继续燃烧膨胀，故Z1和Z3区域无法形成高温膨胀核心.初始火核更容易生成在ξMR等值线上标量耗散率小的区域.

图6 湍流混合层中高温膨胀核心分布Fig.6 Distribution of high-temperature expanding kernels in the turbulence-like mixing layer

为进一步分析湍流混合层中燃料着火和火焰发展过程，图7 展示着火后温度梯度、OH 质量分数、温度峰值和CH4质量分数在混合分数空间的分布特征.由温度梯度变化和高温化学表征组分OH 质量分布可知，高温化学发生在ξMR,2=0.275 位置(浓混合物且当量比大于1)，然后向两侧传播，分别进入正庚烷侧和甲烷侧，甲烷侧的反应强度远大于正庚烷侧火焰前沿的反应强度.稀引燃油侧传播的火焰前沿逐渐向化学当量比混合分数(ξst=0.024 9)位置传播.2D 湍流混合层算例由于混合层厚度较大以及着火后膨胀作用，ξst位置超出计算区域外，未发现边缘火焰火焰结构和预混甲烷的引燃过程[4].图7c 显示着火后计算区域温度峰值在混合物空间的变化，发现温度峰值逐渐向化学当量比混合分数位置传播.图7d 展示了甲烷在浓混合物逐渐生成、在稀混合物中快速消耗的过程.第二阶段着火后，浓混合物侧火焰前沿由于氧含量和其他活性基团不足，甲烷和正庚烷无法完全燃烧，正庚烷裂解反应生成甲烷.稀混合物传播分支，氧气含量增多，反应温度逐渐升高，甲烷成分逐渐燃烧消耗.

图7 不同时刻下温度梯度、温度峰值、羟基和甲烷质量分数随混合分数的变化Fig.7 Conditional averaged temperature gradient，peak temperature，mass fraction of CH4 and OH against mixture fraction at different time instants

2.2 湍流场化学爆炸模式分析

化学爆炸模式分析(chemical explosion model analysis，CEMA)是由Lu 等[19]提出的一种定量检测方法，可以用来识别临界火焰特征(熄火极限和火焰前沿)、火焰和着火结构等，已在相关研究中广泛应用[19-22].CEMA 基于化学源项雅可比矩阵 Jw的特征值确定混合物自燃倾向.

对反应系统离散守恒由源项和混合项构成，即

式中：y 为因变量向量，包括组分浓度和其他稳态变量，在CEMA 中，组分摩尔浓度和温度包括在其中；为物质导数算子；w 为源项；S 为混合项.

对控制方程的两边取雅可比矩阵，主要包括化学源项的雅可比矩阵以及混合项的雅可比矩阵，即

式中：Ja为因变量向量物质导数的雅克比矩阵；Js代表混合项的雅克比矩阵.CEMA 分析中，通过化学雅可比矩阵确定临界的火焰特征.通过求解雅可比矩阵特征值λe，并利用特征值分析化学反应特征[19,22-23]，其表达式为

进而可得CEMA 特征值Λe，且当Λe＞0 对应化学爆炸模式(CEM)，Λe＜0 对应非化学爆炸模式(Non-CEM)，Λe计算公式为

式中：be和αe分别是雅可比矩阵的左、右特征向量；eλ为雅可比矩阵的特征值；sgn()为符号函数；Re()表示取特征值的实数部分.

化学爆炸模式(CEM)下的反应物活性较高，自燃倾向大，且Λe越大表示反应物活性越高，自燃倾向越大，着火后燃烧剧烈，通常对应未燃反应区；非化学爆炸模式(Non-CEM)对应的区域反应活性较低，无自燃倾向.研究表明，CEM 向Non-CEM 的转变过程与临界火焰结构、着火以及预混火焰前沿位置密切相关，可用来确定着火和火焰前沿的位置.同时，Λe的零临界点可用于区分已燃区域和未燃区域.笔者采用CEMA 分析混合层中自燃和火焰传播特征，关于CEMA 方法的详细介绍可参见文献[21].CEMA 分析主要针对混合物中单点的组分浓度以及温度，未考虑混合物中的扩散混合作用，但并不意味混合层中组分和热量扩散作用不重要.

图8～图10 分别展示不同时刻下计算区域内Λe、O H 质量分数和温度的分布情况.图8a 所示t＝0.313 ms 时刻计算区域中央呈现窄条状的青色区域(Λe＜0)，存在非化学爆炸模式(Non-CEM)，其形成主要与冷焰反应有关，主要集中在第一阶段着火后浓混合物中缓慢发生低温化学反应的区域.分析该时刻下OH 组分质量分数和温度分布，如图9a 和图10a 所示，发现Non-CEM 区域右侧区域温度逐渐上升，形成多个高活性核心，但处于化学反应初期，化学反应活性仍很强，故对应的位置为化学爆炸模式(CEM).

图8 不同时刻下计算区域Λe 空间分布Fig.8 Spatial distribution ofΛe in the computational domain at different times

图9 不同时刻下计算区域OH空间分布Fig.9 Spatial distribution of OH concentration in the computational domain at different times

图10 不同时刻下计算区域温度空间分布Fig.10 Spatial distribution of temperature in the computational domain at different times

图8b 所示t＝0.383 ms 时刻第二阶段着火发生后，着火核心燃烧后发热膨胀，形成高温膨胀核心，因而离散的高温膨胀核心着火后形成蓝色(Λe＜0)的Non-CEM 区域.实际上，难以找到合适的OH 组分质量分数限值来确定反应前沿，在图9b 中该时刻下OH 质量分数较高处和已燃区并非完全匹配，而使用CEMA 分析方法可以量化方式准确地确定火焰前沿和反应区域.

反应继续进行至t＝0.403 ms 时刻，多个位置相继进入第二阶段着火时期，已存在的高温膨胀核心膨胀连结，开始在空间上形成两个独立的区域.同时，正庚烷侧下游存在冷焰形成的条状Non-CEM 区域.高温已燃区和冷焰区中间的中温反应区仍然保持较高的自燃倾向.t＝0.430 ms 时刻低温冷焰和高温已燃Non-CEM 区域已经连结，继续向两侧传播.

3 结论

(1) 第一阶段着火后，湍流作用下在ξMR,1等值线边缘偏浓区域逐渐生成冷焰；第二阶段着火后，在ξMR,2等值线偏甲烷侧生成多个高温的膨胀核心；高温核心在混合分数平缓的区域更易生成；混合分数梯度过大会增大标量耗散率，加强热量和活性基团的耗散，不利于燃烧反应的稳定进行.

(2) 不同于层流混合层中火焰前沿，湍流作用下火焰前沿形成褶皱向两侧传播；位于热膨胀核心位置的火焰前沿传播较快，混合分数梯度大的区域反应进行缓慢，反应前沿传播滞后；甲烷侧氧含量高，前沿传播加快，火焰前沿褶皱程度逐渐降低；正庚烷侧火焰前沿在传播下游存在冷焰反应区域，呈现“双火焰”结构，冷焰位置低温化学组分KET 大量生成，火焰前沿处温度高，高温化学反应组分OH 大量生成；随着反应进行，火焰前沿传播至更浓混合物中，双火焰面之间的距离逐渐缩短.

(3) 由混合分数空间上温度、组分和温度梯度时间变化特征可知，甲烷侧的化学反应强度更大，反应前沿逐渐向ξst位置传播；另外，甲烷组分在混合分数空间上同时存在生成和消耗反应区，甲烷侧的温度高且氧气充分，甲烷分子氧化消耗，而正庚烷侧氧含量低，正庚烷无法完全氧化，部分裂解反应生成甲烷分子；CEMA 分析发现特征值可以准确地区分冷焰和高温膨胀核心区，相较于表征组分等值线的方式，可以量化确定反应区、火焰前沿和双燃料NTC 效应影响的中温反应区.