0.   引言

航空发动机与燃气轮机被誉为装备制造业“皇冠上的明珠”，是保证国防与能源安全，保持工业竞争力的核心战略产业。我国已启动的“两机专项”被列为国家“十三五”期间百个大项目之一，是工程研究的重中之重。目前我国航空发动机和燃气轮机的发展仍有较大不足，而燃烧室作为航空发动机和燃气轮机实现能量转化的核心部件之一，通常工作在高温、高压等极端环境之下，其设计水平直接影响到整机的性能^[1], 对其开展深入研究具有极大的现实意义。

近年来，燃烧室内旋流火焰结构和机理的研究越来越引起相关科研工作者的重视。在实验分析方面，Watanabe^[2]和Zhou^[3]等在常温常压下对旋流预混火焰进行了实验研究，分析了单旋流燃烧器中的火焰结构，阐释了组分扩散与热扩散间的作用关系，并对湍流流场进行了测量和分析；Li^[4]和Wang^[5]等对航空发动机贫燃预混预蒸发模型燃烧室进行了实验，研究了高温高压下航空发动机内的污染物生成机制；Worth^[6]和Singh^[7]等则着重分析了环形燃烧室内火焰和流场的相互作用关系。在数值模拟方面，苗文博^[8]等采用二维化学非平衡Navier-Stokes方程分析了空间发展的超声速燃烧反应剪切层的放热效应；Butz^[9]采用大涡模拟(Large Eddy Simulation, LES)方法对湍流预混旋流火焰进行了模拟，发现了内部回流区对燃尽产物的阻滞作用；Tanaka^[10]和Aoki^[11]等采用直接数值模拟(Direct Numerical Simulation, DNS)方法对常温常压下预混甲烷/空气旋流火焰进行了计算，从机理层面分析了旋流火焰和流场的结构以及湍动能的来源。

然而由于受高压强湍流条件下的实验诊断技术和DNS计算量的限制，大部分研究均未考虑高压强条件下湍流特征量与实际航空发动机燃烧室特征量的内在关系。与常压燃烧最大的区别在于，高压预混燃烧的火焰面极薄。在湍流预混燃烧理论中，常使用湍流均方根速度u′、湍流积分长度尺度l_t、湍流Kolmogorov尺度η、层流火焰速度S_L和层流火焰厚度δ_L等特征量来定义无量纲特征量如u′/S_L、l_t/δ_L、Karlovitz数(Ka)等。对于高压火焰，极薄的火焰厚度直接影响着上述参数的量级。对于航空发动机燃烧室内的湍流预混燃烧，Ka大致为10⁰~10²，在这一参数范围内，预混燃烧模式位于薄反应区和破碎火焰区中^[12]。

近期也有学者针对不同Ka下火焰开展了研究工作。Carlsson^[13]通过对各向同性湍流燃烧进行直接数值模拟，研究了高Ka下的破碎火焰结构特征；MacArt^[14]同样采用DNS方法模拟分析了不同Ka下平板预混燃烧射流中标量通量和梯度分布；Wang^[15]则对圆柱预混射流进行了直接数值模拟，分析了较高Ka下的火焰拉伸和火焰增厚现象。然而基于不同Ka的旋流火焰的研究还有所不足，对燃烧室内相应旋流火焰特征的理解仍需加强。

基于此，为了进一步理解航空发动机燃烧室内旋流预混燃烧的火焰结构特征和机理，本文设计构造了高温高压环形预混旋流模型燃烧室，并采用直接数值模拟方法研究了不同Ka的火焰和湍流结构特征，分析了旋流火焰中标量通量和梯度的对齐关系，讨论了对相关模型的影响，为发展适用于航空发动机燃烧的高精度模型和方法提供基础。

1.   直接数值模拟

1.1   计算模型

实际航空发动机燃烧室多为环形结构，且入口燃料温度高，室内压力和湍流强度都极高。考虑到直接数值模拟计算量的限制，对实际结构等比例建模较为困难，所以基于典型工况(p、T、U₀)及其对应的典型湍流预混火焰参数(u′/S_L、l_t/δ_L、Ka)设计了简化的计算模型。本研究中将航空发动机内的环形燃烧室简化为两个相邻的预混旋流。图 1所示为计算域示意图，在z方向上采用周期性边界条件，而在y方向上采用壁面边界条件以模拟实际环形燃烧器。单个环状旋流入口的内径D_in设置为0.1 mm，外径D_out设置为0.3 mm，两相邻旋流中心间距为2D_out，计算域整体尺寸为5.5D_out×2D_out×4D_out(x×y×z)。

图  1  计算域示意图(涡量等值面采用温度着色，左侧切片图为热释率)

Figure  1.  Schematic illustration of the computational domain (the isosurface colored by temperature presents vorticity and the slice on the left presents heat release rate)

下载: 全尺寸图片 幻灯片

为了研究不同Ka下的湍流火焰特征，设置了两个工况，其入口速度U₀分别为40 m/s和15 m/s，分别对应于航空发动机燃烧室内的高/低负荷工况，其入口雷诺数Re₀分别约为4000和1500，如表 1所示。其余参数将在后文详细介绍。

表  1  计算工况

Table  1.  Calculation conditions

	U0/(m·s-1)	Re0	Ka($\widetilde{C}$=0.8)	u′/SL($\widetilde{C}$=0.8)	Lt/δL($\widetilde{C}$=0.8)
1	40	4000	40	17.2	18.0
2	15	1500	8.5	5.7	12.1

下载: 导出CSV 

| 显示表格

为了获得具有实际意义的边界条件，首先单独计算了环形管道湍流，待其充分发展稳定后切片保存生成了入口文件，然后在正式计算时对入口文件进行读取从而得到流场入口信息。计算中，旋流数为1.0，燃料为航空煤油，当量比为0.7，入口预混气体温度为860 K，计算域压力为3 MPa。基于以上条件，通过计算一维层流火焰可以得到其层流火焰速度S_L=0.738 m/s, 层流火焰厚度δ_F=18.2 μm。

1.2   数值方法

本研究采用了变密度、低马赫数的燃烧求解器，求解了连续性方程、动量守恒方程、能量守恒方程和组分守恒方程^[16]。在气体预混燃烧直接数值模拟中，影响计算网格的因素主要有两个方面:从火焰角度，为了准确解析火焰面处组分扩散和燃烧放热效应，在火焰面处需要保持至少8~10个网格^[17]；从湍流角度，网格尺度需要能解析湍流最小特征尺度(η)的涡结构，通常而言，需要满足η>Δx_min/2^[18]。本文模拟中最小涡尺度在入口附近产生，约为1.5 μm。综合考虑火焰和湍流的精度要求，本研究中网格数量(N_x×N_y×N_z)设置为768×256×512≈9.3×10⁷。其中在x方向的前半段采用2 μm的均匀网格，在后半段采用拉伸率为1.05的拉伸网格，在y、z方向均采用大小为2.3 μm的均匀网格。为了降低计算量，采用高温高压条件下两步航空煤油/空气化学反应机理^[19]。

2.   旋流预混火焰分析

2.1   旋流预混燃烧状态

为了确定本文湍流预混燃烧火焰的参数特征，需在火焰面位置统计各湍流燃烧参数。本文中，火焰面根据过程变量C来定义：

式中Y_{O2, u}为未燃气体中氧气初始质量分数，Y_{O2, b}为燃尽物中氧气质量分数，Y_O2为当地氧气质量分数，C=0表示未燃物，C=1表示燃尽物。

图 2为工况1的y=0 mm截面上时间平均后的过程变量$\widetilde{C}$的云图。工况2与工况1基本一致，故不加赘述。可以发现，在旋流火焰上游，$\widetilde{C}$最小值接近0，而在火焰面下游，$\widetilde{C}$最小值均大于0.4，这表明随着旋流火焰的发展，主流中氧气质量分数逐渐降低，且主流逐渐沿着径向扩张。结合图 5得知，$\widetilde{C}$最小处为旋转主流所在位置，$\widetilde{C}$≈1处为旋流中心，$\widetilde{C}$处于中间值则为内外剪切层附近的位置。

图  2  工况1中y=0 mm截面上的过程变量分布

Figure  2.  Distribution of progress variable at y=0 mm of case 1

下载: 全尺寸图片 幻灯片

如前文所述，在实际航空发动机燃烧器中，湍流预混燃烧特征参数处于一定的分布区间，其中具有代表性意义的无量纲特征量有u′/S_L、l_t/δ_L和Ka等。其中Ka从物理意义上表示层流火焰特征时间尺度与最小涡特征时间尺度的比值，定义为：

式中τ_L为层流火焰特征时间尺度，τ_η为最小涡特征时间尺度，ε为湍流耗散率，ν为运动黏度。

图 3为高/低负荷工况中不同流向位置上基于过程变量$\widetilde{C}$的Ka的统计。两个工况中随着$\widetilde{C}$的增大，Ka均逐渐减小，表明在剪切层和旋转主流中Ka较大，而在旋流中心Ka较小。可以发现，对于入口雷诺数较大的工况1，其Ka相比于工况2显著增大，即在湍流强度增加时火焰尺度与最小涡尺度的比值增大。这是因为较高入口雷诺数导致了较大的湍流脉动，从而增强了湍流耗散率。根据式(2)，Ka的变化主要源于湍流耗散率的变化，而湍流耗散率主要取决于当地密度、黏度和湍流脉动速度。在火焰面前后，温度的升高导致当地密度减小、运动黏度增大和脉动速度的减小，从而导致湍流耗散率的减小。因此，Ka的变化趋势出现转折处即为火焰面位置。在工况1中，转折点位于$\widetilde{C}$≈0.8处，而在工况2中，转折点位于$\widetilde{C}$≈0.9处。这表明在较大Ka时，小涡在火焰厚度尺寸内的扰动更加显著，火焰面位置更趋向于反应物侧。

图  3  不同x位置截面Ka基于过程变量的分布

Figure  3.  Variation of Ka based on progress value at different x location

下载: 全尺寸图片 幻灯片

在这两个工况中，Ka的值均在10⁰~10²范围内，与航空发动机燃烧室工况的参数一致。对于其他参数，图 4显示了湍流预混燃烧模式分区图^[20]，其中横纵坐标分别为湍流预混燃烧无量纲特征值l_t/δ_L和u′/S_L。可以看到，本文设计的两个工况参数均分布于“Aircraft engines”工况范围内。

图  4  湍流预混燃烧模式分区图(三角形：工况1；矩形：工况2)

Figure  4.  Regime diagram for premixed turbulent combustion(triangle: case 1; rectangle: case 2)

下载: 全尺寸图片 幻灯片

2.2   火焰和流场特征

图 5、图 6分别显示了某个时刻两个工况在y=0 mm处的瞬时速度和瞬时涡量云图。旋流预混气体从喷口进入燃烧室后，由于旋转离心导致产生张角，流场呈双“V”字形。在燃烧室内，由于旋转射流内外剪切层处存在较为强烈的剪切作用，从而产生速度小于0的回流区。回流区主要分布在旋流中心区域和两旋流间的角落处，如图 5所示。相比于单旋流火焰，双旋流火焰在火焰下游处由于存在相邻旋流的相互作用，导致下游湍流强度显著增强，另外由于环形燃烧室内相邻旋流之间有限的空间，导致双旋流火焰下游的扩张效应明显弱于单旋流火焰^[10]，即其内部回流区小于单旋流火焰，且外部回流区相比于单旋流而言也更小。相比于低Ka工况，高Ka工况下速度场褶皱更强烈，这与图 6中的涡量场相对应。在工况1中，湍流涡的分布更加细小，对火焰的影响更大。

图  5  y=0 mm截面处的瞬时速度云图(单位: m/s)

Figure  5.  Instantaneous velocity contours at y=0 mm (unit: m/s)

下载: 全尺寸图片 幻灯片

图  6  y=0 mm截面处的瞬时涡量云图

Figure  6.  Instantaneous vorticity contours at y=0 mm

下载: 全尺寸图片 幻灯片

图 7所示为瞬时热释率云图。火焰面在湍流作用下弯曲褶皱，表明在当前工况下湍流预混燃烧模式处于薄反应区中，如图 4的燃烧模式分区图所示。与常压旋流火焰相比，高压下的火焰面厚度极薄，因此存在于更小的湍流尺度范围内，火焰面受小尺度湍流的影响更显著，对应地，在预混湍流火焰燃烧模式分区图中，高压火焰的横坐标l_t/δ_L的值较常压火焰更大，火焰褶皱也就更剧烈。对比图 5可以发现，在靠近入口的上游燃烧区，火焰面沿着旋流内外剪切层双侧分布，而在入口的下游燃烧区，火焰面相互作用，变得更加褶皱甚至破碎，这在x=0.3 mm截面上更加突出。特别是工况1中，在入口速度更高的高负荷情形下，由于其较高的Ka，下游火焰的破碎更加明显。

图  7  瞬时热释率云图

Figure  7.  Instantaneous contours of heat release rate

下载: 全尺寸图片 幻灯片

2.3   速度分布统计

如图 5、图 7所示，预混旋流喷射入燃烧室后，在上游稳定燃烧，而在下游发生两旋流间的相互作用。为了理解旋流湍流预混燃烧中Ka的影响，分别对x=0.15 mm(上游)和x=0.30 mm(下游)两处的截面进行了速度分布的统计比较。

图 8显示了轴向、径向和周向的平均速度在过程变量空间上的分布。在上游x=0.15 mm且$\widetilde{C}$≈0处，两工况在三个方向上的初始平均速度均相等，排除了初始入口条件对流场发展的影响。整体而言，随着过程变量的增大，三个方向的平均速度均逐渐减小，甚至出现负值，这是因为$\widetilde{C}$=0位于旋转主流，各方向速度分量均较大，而在$\widetilde{C}$值较大的剪切层和中心回流区，剪切作用和回流作用导致速度逐渐减小。工况1中由于高负荷下Re较大，周围流体对旋转射流的剪切作用更强，导致速度降低得更为显著。

图  8  平均速度统计

Figure  8.  Statistics of mean velocity

下载: 全尺寸图片 幻灯片

从平均轴向速度来看，在x=0.15 mm、$\widetilde{C}$≈1处，无量纲轴向速度接近于0，说明在旋流上游中心处，不存在回流效应；而在下游x=0.30 mm处，无量纲轴向速度小于0，代表了回流区的存在。对比上下游平均速度可以发现，当$\widetilde{C}$较小时，下游轴向速度更大，而当$\widetilde{C}$较大时，下游轴向速度更小，表明燃烧导致的气体膨胀挤压作用会使旋转主流速度得到增强，而使周围流体的平均轴向速度降低。与平均轴向速度不同的是，随着旋流的发展，平均径向速度均得到增强，这是因为燃烧膨胀和内回流区的卷吸作用导致内回流区产生了较大的压强。而对平均周向速度而言，随着旋流的发展下游环状旋转主流的半径逐渐增大，根据角动量守恒原理，下游周向速度明显小于上游。

图 9所示为轴向、径向和周向的脉动速度统计。在上游处(x=0.15 mm)，由于工况1的雷诺数更大，湍流扰动效应也更强，所以脉动速度均大于工况2。而在下游情况则有所不同，对轴向和径向脉动而言，虽然在$\widetilde{C}$较小时同样是工况1更大，但是在$\widetilde{C}$较大时，工况2的轴向和径向速度更大，说明在下游的剪切层和回流区中，低Re时轴向和周向脉动会更强烈。

图  9  脉动速度统计

Figure  9.  Statistics of fluctuating velocity

下载: 全尺寸图片 幻灯片

3.   梯度输运检验

在雷诺平均(Reynolds Averaged Navier-Stokes, RANS)方法和大涡模拟(Large Eddy Simulation, LES)方法中，质量加权的标量通量F_{j, k}是守恒方程的重要项，其定义为：

式中，Y_k为组分k的质量分数。

然而，在RANS或LES中，根据式(3)所得到的标量通量项不能封闭，通常采用梯度扩散假设来建模，如：

式中：ν_t为涡黏度，Sc_t为湍流施密特数。

从式(4)中可以看出，梯度扩散假设是假设组分通量$\widetilde{u^{\prime \prime}_{j} Y^{\prime \prime}_{k}}$和组分梯度$\frac{\partial \widetilde{Y_{k}}}{\partial x_{j}}$之间存在完全反对齐关系，这种反对齐关系被称为梯度输运。而相反地，当组分通量和梯度为对齐关系时，为逆梯度输运。此前的研究表明，当Ka较小时，逆梯度输运发挥主导作用，而在Ka较大时，梯度输运效应更为显著^{[14, 21-22]}。

为研究航空发动机燃烧室内典型梯度输运假设的有效性，本节对组分扩散和组分梯度的对齐性进行了统计。

cosθ_k为组分k通量与组分梯度间夹角的余弦值。当cosθ_k→-1时，通量与梯度为反对齐关系，即梯度扩散占主导地位；当cosθ_k→1时，通量与梯度为对齐关系，即逆梯度扩散占主导地位。

图 10为基于主要反应物、生成物和中间组分的标量通量和梯度的对齐性统计。可以发现，反应物(O₂)和生成物(H₂O)的对齐性在各位置的过程变量空间中基本保持一致，夹角的余弦值均小于-0.4，且均在火焰面附近取得最大值，表明采用梯度扩散假设来对主要反应物和生成物通量进行建模是基本合理的。且由于工况1的Ka更大，其梯度扩散效应更明显。

图  10  基于反应物、生成物和中间组分的标量通量和梯度的对齐性统计

Figure  10.  Alignment statistics of scalar flux and gradient for fuel, product and radial mass fraction

下载: 全尺寸图片 幻灯片

然而对于中间产物，虽然同样当Ka较大时，夹角的余弦值更小，但是在火焰面附近，cosθ_CO的最大值约等于0.2，表明组分通量更偏向于与组分梯度呈对齐关系，即逆梯度扩散效应更为显著，传统梯度扩散模型失效，如图 10(c)所示。这表明传统上采用的基于梯度扩散的模型对于发动机燃烧室中间产物的模拟存在问题，在发展适用于航空发动机燃烧模拟的高精度模型和方法时，需要对相关通量项的封闭进行进一步的改进和优化。

4.   结论

基于航空发动机燃烧室湍流预混燃烧的特征参数，构造了环形旋流预混模型燃烧室，并研究了不同Ka条件下湍流火焰的区别及其内在机制，加强了对航空发动机燃烧室内预混湍流火焰和流场结构的理解。结果表明，发动机燃烧室内湍流预混燃烧模式位于薄反应区，且Ka越大时火焰的弯曲褶皱越强烈。火焰面沿着剪切层内外呈双侧分布，与单旋流燃烧不同的是，双旋流火焰在下游由于相邻旋流的碰撞而产生较为破碎的火焰结构，这在高Ka的工况中尤为明显。由于旋流剪切层的作用，导致在下游产生内部回流区，且Re和Ka越大，回流效应越显著。在剪切层附近，脉动速度最大，同样Re和Ka越大时，湍流扰动效应也越强。

在典型的RANS或LES方法中，常采用梯度扩散假设对组分通量进行建模，为了验证这一方法在发动机燃烧室内的适用性，对梯度输运效应进行了统计分析。分析表明，在航空发动机燃烧室工况下，反应物和生成物基本满足梯度输运效应，且Ka越大，梯度输运效应越显著，这与此前基于各向同性湍流和平板射流的研究结果相符。然而对于中间产物，如CO，梯度扩散假设在火焰面附近严重失效。因此在发展适用于航空发动机燃烧室的高精度数值模拟方法时，需要对标量通量项的封闭进行进一步的改进。