激波冲击火焰的涡量特性研究*

朱跃进,董 刚

(1.江苏大学能源与动力工程学院,江苏 镇江 212013; 2.南京理工大学瞬态物理重点实验室,江苏 南京 210094)



激波冲击火焰的涡量特性研究*

朱跃进1,2,董 刚2

(1.江苏大学能源与动力工程学院,江苏 镇江 212013; 2.南京理工大学瞬态物理重点实验室,江苏 南京 210094)

激波冲击火焰的现象涉及一系列复杂的物理化学过程，其中涡量的生成与演化对控制火焰发展起重要作用。为系统分析激波冲击火焰过程中的涡量特性，采用二维带化学反应的Navier-Stokes方程对平面入射激波及其反射激波与球形火焰作用的现象进行了数值研究，通过引入并行计算达到高网格分辨率的要求。计算结果表明，斜压项对火焰区内涡量生成起主导作用，压缩项和耗散项在火焰膨胀阶段抑制涡量生成，此外，火焰在激波压缩阶段主要受物理过程而非化学反应过程影响。

爆炸力学;涡量;Navier-Stokes方程;火焰;激波;斜压项

激波冲击火焰的现象常常出现在惯性约束聚变[1]、超声速燃烧推进[2]和工业爆炸灾害[3]等领域，涉及复杂的物理化学机制，相关研究有利于对激波-火焰作用现象的抑制或利用，因而具有重要的科学意义和应用价值。

G.H.Markstein[4]的经典实验首次清晰显示了火焰在弱激波及其反射激波作用下的变形失稳过程；英国Wales大学的G.O.Thomas等[5]实验研究了平面激波冲击球形火焰的过程，他们通过高速摄影发现，较强的入射激波及平面反射激波可使火焰形状发生严重变形。与此同时，随着计算机技术的发展，许多数值模拟工作得以开展。G.A.Bately等[6-7]数值研究了压力波及其反射波与柱形火焰的相互作用过程，计算结果表明，当反射波与火焰作用时，黏性效应对变形火焰的涡量演化有较大影响；Y.Ju等[8]研究了不同强度的入射激波对火焰变形的影响，结果发现火焰区的总燃烧速率随激波强度的增加而增加；美国海军研究实验室的E.S.Oran等[3]和A.M.Khokhlov等[9-11]对激波冲击火焰的问题展开了一系列研究，认为Richtmyer-Meshkov(RM)不稳定性是火焰变形失稳的主要机制, 且变形火焰的燃烧速率和放热率均有所提高，当条件适当时，反射激波与火焰的相互作用还会产生热点和爆轰。本文作者[12-15]采用Navier-Stokes(NS)方程，考察了化学反应机理、网格尺寸、激波强度、化学反应活性、火焰表面积等不同因素对激波-火焰作用过程的影响，并通过可视化图像、火焰积分量、火焰界面几何特性等方面的分析来探讨火焰变形、燃烧和混合等过程的变化规律及相互作用关系。

上述研究均表明，由于存在斜压效应，激波冲击火焰时可在火焰面处产生大量涡量，从而影响火焰界面变形、燃料与氧化剂混合、燃烧放热等一系列物理化学现象，因此有必要深入研究激波冲击火焰过程中的涡量特性，但目前尚缺乏对上述问题的系统研究。

本文中将采用高分辨率网格和基于MPI信息传递[16]的并行方法，开展入射及其反射激波与球形火焰多次作用的二维数值模拟，根据计算结果对变形火焰区域内的涡量特性进行深入分析研究，具体包括涡量演化过程、涡动力学分析、环量分析等，以期为理解激波冲击火焰这一经典问题提供新的视角和思路。

1 数理模型与计算

1.1 数理模型

采用的二维带化学反应的NS方程可表示为：

(1)

(2)

(3)

(4)

式中：ρ为密度；ui为i方向速度(i=1，2)；τij为黏性应力张量；E为单位体积总能量，E=p/(-1) ++ρqY；p为压力；为比热比；q为单位质量的总化学能；Y为反应物质量分数；qj为热通量，qj=-kT/xj；为化学反应速率；热传导系数k、扩散系数D、运动学黏性υ 的表达见文献[3]。

控制方程(1)～(4)采用分裂算法求解：空间导数项的无黏部分采用高精度五点TVD格式[13-15]计算，黏性部分采用二阶空间中心差分计算，时间推进过程采用二阶Runge-Kutta方法求解。本文采用单步化学反应来描述燃烧过程，相关化学反应热力学和动力学参数可参见文献[14]。

1.2 计算条件

计算以文献[5]中的实验为依据，设计如图1所示的计算区域(灰色部分)与初始流场：计算区域流向(x)长0.17 m，法向(y)高0.038 m；球形火焰中心位于(x=0.035 m,y=0 m)处，初始时刻半径R0=0.019 m，内部密度ρ1=15.78 g/m3；火焰外部反应性预混气体组成为C2H4+3O2+4N2，密度ρ0=161.5 g/m3，初始温度T0=293 K；火焰内外压力均为p0=13.3 kPa；入射激波马赫数为1.7，初始时位于x=0.012 m处，激波沿x方向从左向右传播，当行至右端壁面后发生反射，左行的反射激波与火焰再次作用；激波后气体状态由Rankine-Hugoniot关系式给出。

数值计算采用的网格尺寸为47.5 μm×47.5 μm，因此本文中考虑的预混气体在入射激波和反射激波作用后的火焰反应区宽度可分别由大约8个、3个网格来刻画，网格尺寸足以精细反映火焰结构[12,15]。

图1 计算区域和初始流场示意图Fig.1 Computational domain and initial flow field

1.3 并行计算

图2 流场分区示意图Fig.2 Schematic of flow field block

由于采用较高的网格分辨率，计算规模较大，因此将采用MPI信息传递方法对已有程序进行并行化处理，即将计算区域划分为N个子分区并初始化，当程序开始运行后，系统同时启动N个进程进行计算。由于本文数值格式在计算网格界面通量时用到界面两侧各一层网格的值，因此在进行每一步的计算之前，需对每个分区的边界条件进行处理，而边界条件分为2种:物理边界和相邻子分区间的内边界。对于物理边界，可根据实际条件直接设置，其中，计算区域上边界(y=0.038 m)和x方向右端面为管壁，均采用无滑移的刚性绝热壁面边界；y=0 m处为对称面，采用对称边界；x方向左端面处采用零梯度边界；而分区内边界则需用MPI库函数对分区间的数据进行传递，如图2所示，边界0是分区A和B之间的内边界，边界1是分区A在分区B内的虚拟边界，边界2是分区B在分区A内的虚拟边界，具体交换方式如下：将分区A内落在边界2和0之间的这一层网格数据通过MPI库函数发送给区域B，作为区域B计算边界通量时需要的值，而将分区B内落在边界0和1之间的一层网格值发送给区域A，供区域A计算边界通量，这种方法保证了在分区内边界处的连续性。

2 结果与讨论

2.1 实验验证

数值方法和化学模型的可靠性通过计算结果与文献[5]中实验的对比来验证。图3显示了激波阵面距火焰点火位置(初始球形火焰中心)的距离、变形火焰宽度和高度三者随时间的变化，所有结果均采用火焰初始半径R0进行量纲一处理，通过图中对比可以看出，激波阵面距点火位置的距离、火焰尺寸的变化均与实验结果吻合较好，验证了本文所用计算方法和化学反应各参数的合理性与可靠性。除采用火焰面位置等宏观量验证，文献[13-14]中还对比了不同时刻的火焰结构形态，发现计算结果与实验结果一致，表明本文计算的可靠性。

图3 激波距点火位置距离、火焰界面宽度和高度计算结果(曲线)与实验结果(符号)[5]的对比Fig.3 Comparisons between calculated (curves) and experimental results (symbols) [5]about distance between shock wave and ignitor, height and width of deformed flame

2.2 涡量演化过程

为了考察激波冲击火焰过程中的涡量演化过程，图4首先给出了各时刻火焰区涡量ω的可视化结果，其中图4(a)～(d)为入射激波波后情形，图4(e)～(h)为反射激波波后情形。从图中可以看出，涡量首先出现在火焰界面上，但随着火焰的旋转变形，界面处的涡量逐渐向火焰中心集中，同时在黏性作用下发生衰减；当反射激波与火焰发生二次作用后，火焰内部涡量得到强化，形成红色与蓝色2个旋转方向相反的主涡量，从而诱导火焰结构向壁面发生翻转，同时火焰边界的小尺度涡量也得到极大增强，进而导致火焰界面皱褶增强。

图4 火焰区涡量变化Fig.4 Vorticity variation of the flame zone

2.3 涡动力学分析

为澄清火焰区域内涡量随时间演变的内在机制，给出二维理想流体的涡量动力学方程：

(5)

图5 火焰区内各输运项绝对值之和随时间的变化Fig.5 Time histories of the sum of absolute values for per transport term within flame region

式中：D/Dt=∂/∂t+·u，等号右侧第1项为压缩项ωc，第2项为斜压项ωb，第3项为耗散项ωd。图5给出了火焰区域内涡量空间上的各输运项绝对值之和随时间的变化过程，分别记为|ωc|、|ωb|和|ωd|，值越大表示该输运项对涡量变化的影响越大，图中曲线的变化过程可分为4个阶段：Ⅰ为入射激波压缩阶段(0～94 μs)，Ⅱ为火焰膨胀阶段(94～328 μs)，Ⅲ为反射激波压缩阶段(328～374μs)，Ⅳ为火焰二次膨胀阶段。从图中可以看出，斜压效应一直是火焰区域内涡量变化的主要原因，特别是在Ⅲ和Ⅳ阶段，反射激波与变形火焰二次作用可引发更强的斜压效应，这表明图4(e)～(h)中涡量变化主要受斜压效应影响。需要注意的是，Ⅳ阶段中压缩项和耗散项变化趋势基本一致，但两者的变化原因并不相同：入射激波在流向右端固壁面反射后，激波强度增强，对变形火焰的压缩作用更强，因而引起火焰区内压缩项的大幅提升，此时，反射波后的温度压力进一步升高，气体黏性会随之增大，因此耗散项在反射波后也得到了较大增强。

图6 火焰区内各输运项正负值之和随时间的变化Fig.6 Time histories of the sum of positive and negative values for per transport term within flame zone

式(5)中各输运项在火焰区域内不同位置具有不同的正负值，除通过图5中各输运项绝对值之和来分析输运项对涡量生成的影响之外，还可利用各输运项正负值之和表征输运项对涡量生成所起的作用，分别记为ωc、ωb和ωd，即输运项正负值之和大于零时促进涡量生成，反之则抑制涡量生成，所得曲线见图6。结果表明，斜压项始终基本起促进涡量形成的作用，而压缩项、耗散项在火焰膨胀时(参见图5中Ⅱ和Ⅳ阶段)对涡量生成主要起抑制作用，因此，涡量强度在火焰膨胀期内降低主要受火焰快速膨胀的影响，尤其在反射波后(Ⅳ阶段)，火焰剧烈膨胀，并且气体黏性进一步增强，从而能够引发火焰区涡量强度快速下降，这与文献[14]中得到的涡量强度变化规律一致。

图5～6结果表明，由密度梯度和压力梯度方向不一致所引起的斜压项是促进火焰区内涡量生成的关键因素。为定量比较反射激波前后斜压项的变化，图7定量给出了入射及反射激波波后特定时刻火焰内斜压项最大值ωb,max和最小值ωb,min的变化。图中结果表明，反射激波波后(图5中Ⅳ阶段)斜压项的强度比入射激波波后(图5中Ⅱ阶段)斜压项强度强；入射激波波后斜压项的强度在缓慢降低；而反射激波波后，斜压项强度在火焰区剧烈的燃烧膨胀和增大的耗散作用下降低幅度较大，但从t=500.7s(图4(f)对应时刻)开始，变形火焰与右端固壁发生碰撞，使得斜压项强度得到一定程度的增强，随后斜压项强度受火焰膨胀和黏性的影响继续降低。

2.4 环量分析

涡量的强度可以通过环量来表达。为进一步研究火焰区涡量特性，考察火焰区内环量随时间的变化过程(见图8)，其中+和-分别代表火焰区正负环量的绝对值，而为正负环量之和，即火焰区的净环量。本文计算环量采用公式如下：

(6)

图7 火焰区内斜压项最大值和最小值随时间的变化Fig.7 Time histories of max and min values of baroclinic term within flame zone

图8 火焰区内环量随时间的变化Fig.8 Time histories of circulationwithin flame zone

式中：下标F代表火焰区域(定义为反应物质量分数Y≤0.99)，为涡量。从图8中可以看出，正负环量在入射及反射激波与火焰作用时(分别对应图5中Ⅰ和Ⅲ阶段)均快速上升，在入射激波及反射激波波后变化幅度放缓。由式(6)可知，环量变化同时受火焰面积和涡量两者变化的影响，而文献[15]表明火焰面积在Ⅰ和Ⅲ阶段内减小，在Ⅱ和Ⅳ阶段内增加，特别是Ⅳ阶段内增加速度较快，因此结合图8可知，火焰区内涡量强度在Ⅰ、Ⅲ阶段内上升，在Ⅱ和Ⅳ阶段内下降，尤其在Ⅳ阶段内下降幅度较大，这与2.3节及文献[14]中的结论一致。

此外，本文中还将上述环量计算结果与激波与惰性气泡作用下的环量模型(PB模型[17]和YKZ模型[18])的计算结果进行了对比。PB和YKZ模型选取激波刚扫过气泡下游极点的时刻，在简化的假设前提下，对半圆型气泡区域内的涡量进行积分。图8中垂直实线分别对应t=36.94、361.2 μs，即入射和反射激波刚好扫过火焰区域的时刻，与PB和YKZ模型计算结果的对比见表1。结果表明，YKZ模型预测结果与本文计算结果吻合较好，而PB模型计算的环量值大致为数值计算结果的2倍，这与G.Layes等[19]的结论一致。传统的PB和YKZ模型均用于入射激波与惰性气泡相互作用的过程，而本文中研究发现YKZ模型还适用于入射及反射激波扫过火焰后的环量预测，这实际上表明在激波压缩火焰阶段(对应图5中I和III阶段)，火焰区的变化主要受物理过程而非化学反应过程控制。

表1 火焰区域内的环量计算结果

3 结 论

基于二维带化学反应的NS方程，对平面入射激波及其反射激波冲击球形火焰的现象进行了数值研究，采用并行计算克服高精度时空分辨率带来的计算量过大问题，主要从涡量可视化、涡动力学及环量的角度分析火焰区内的涡量特性。计算结果表明，火焰区内涡量在激波压缩阶段快速增长，在火焰膨胀阶段下降；斜压项在火焰区内涡量变化过程中起主导作用，主要促进涡量生成，而压缩项和耗散项在火焰燃烧膨胀阶段均抑制涡量生成；此外，环量计算还表明在激波压缩阶段，火焰区变化主要受物理过程控制。

[1] Lindl J D, McCrory R L, Campbell E M. Progress toward ignition and burn propagation in inertial confinement fusion[J]. Physics Today, 1992,45(9):32-40.

[2] Marble F E, Hendrick G J, Zukoski E E. Progress toward shock enhancement of supersonic combustion process[R]. AIAA, 1987:87-1880.

[3] Oran E S, Gamezo V N. Origins of the deflagration-to-detonation transition in gas-phrase combustion[J]. Combustion and Flame, 2007,148(1/2):4-47.

[4] Markstein G H. A shock-tube study of flame front-pressure wave interaction[C]∥6th Symposium (International) on Combustion. Pittsburgh, USA: The Combustion Institute, 1957:387-398.

[5] Thomas G O, Bambrey R, Brown C. Experimental observations of flame acceleration and transition to detonation following shock-flame interaction[J]. Combustion Theory and Modeling, 2001,5(4):573-594.

[6] Batley G A, Mcintosh A C, Brindley J, et al. A numerical study of the vorticity field generated by the baroclinic effect due to the propagation of a planar pressure wave through a cylindrical premixed laminar flame[J]. Journal of Fluid Mechanics, 1994,279:217-237.

[7] Batley G A, Mcintosh A C, Brindley J. The baroclinic effect in combustion[J]. Mathematical and Computer Modelling, 1996,24(8):165-176.

[8] Ju Y, Shimano A, Inoue O. Vorticity generation and flame distortion induced by shock flame interaction[C]∥27th Symposium (International) on Combustion. Pittsburgh, USA: The Combustion Institute, 1998:735-741.

[9] Khokhlov A M, Oran E S, Chtchelkanova A Y, et al. Interaction of a shock with a sinusoidally perturbed flame[J]. Combustion and Flame, 1999,117(1/2): 99-116.

[10] Khokhlov A M, Oran E S, Thomas G O. Numerical simulation of deflagration-to-detonation transition: The role of shock-flame interactions in turbulent flame[J]. Combustion and Flame, 1999,117(3):323-339.

[11] Khokhlov A M, Oran E S. Numerical simulation of detonation initiation in a flame brush: The role of hot spots[J]. Combustion and Flame, 1999,119(4):400-416.

[12] Dong G, Fan B C, Ye J F. Numerical investigation of ethylene flame bubble instability induced by shock waves[J]. Shock Waves, 2008,17(6):409-419.

[13] 朱跃进，董刚，范宝春.受限空间内激波与火焰作用的三维计算[J].推进技术,2012,33(3):405-411. Zhu Yue-jin, Dong Gang, Fan Bao-chun. Three-dimensional computation of the interactions between shock waves and flame in a confined space[J]. Journal of Propulsion Technology, 2012,33(3):405-411.

[14] Zhu Y J, Dong G, Liu Y X. Three-dimensional numerical simulations of spherical flame evolutions in shock and reshock accelerated flows[J]. Combustion Science and Technology, 2013,185(10):1415-1440.

[15] 朱跃进，董刚，刘怡昕，等.激波诱导火焰变形、混合和燃烧的数值研究[J].爆炸与冲击,2013,33(4):430-437. Zhu Yue-jin, Dong Gang, Liu Yi-xin, et al. A numerical study on shock induced distortion, mixing and combustion of flame[J]. Explosion and Shock Waves, 2013,33(4):430-437.

[16] 都志辉.高性能计算之并行编程技术----MPI并行程序设计[M].北京:清华大学出版社,2001.

[17] Picone J M, Boris J P. Vorticity generation by shock propagation through bubbles in a gas[J]. Journal of Fluid Mechanics, 1988,189:23-51.

[18] Yang J, Kubota T, Zukoski E E. A model for characterization of a vortex pair formed by shock passage over a light-gas inhomogeneity[J]. Journal of Fluid Mechanics, 1994,258:217-244.

[19] Layes G, Jourdan G, Houas L. Experimental study on a plane shock wave accelerating a gas bubble[J]. Physics of Fluids, 2009,21(7):074102.

(责任编辑 曾月蓉)

A study of vorticity characteristics of shock-flame interaction

Zhu Yue-jin1,2, Dong Gang2

(1.SchoolofEnergyandPowerEngineering,JiangsuUniversity,Zhenjiang212013,Jiangsu,China; 2.KeyLaboratoryofTransientPhysics,NanjingUniversityofScienceandTechnology,Nanjing210094,Jiangsu,China)

The phenomenon of shock wave interacting with a flame involves a series of complicated physical and chemical processes, in which the generation and evolution of vorticity play an important role in controlling flame development. To systematically analyze the vorticity characteristics in the course of shock-flame interaction, a numerical study of a planar incident shock wave and its reflected wave interaction with a spherical flame was carried out by using the two-dimensional Navier-Stokes equations coupled with chemical reaction, and the requirement of high-resolution grid was met via the parallel computation. It is found that the baroclinic term plays a dominant role in the generation of vorticity within the flame zone, and the compression and dissipation terms restrain the generation of vorticity in the flame expanding stages. Besides, in the compression stages, the evolution of flame is mainly affected by the physical-rather than chemical-process.

mechanics of explosion; vorticity; Navier-Stokes equations; flame; shock wave; baroclinic term

10.11883/1001-1455(2015)06-0839-07

2014-04-29；

2014-08-21

国家自然科学基金项目(11372140,11402102);江苏省自然科学基金青年项目(BK20140524); 江苏大学高级专业人才科研启动基金项目(14JDG031);江苏省博士后基金项目(1402013B)

朱跃进(1986— )，男，博士，讲师,zyjwind@163.com。

O382 国标学科代码： 13035

A