中心分级燃烧室冷态流场的大涡模拟

于 雷，王成军，赵冬凯，刘爱虢

(沈阳航空航天大学 航空航天工程学部(院)，沈阳 110136)

中心分级燃烧室冷态流场的大涡模拟

于 雷，王成军，赵冬凯，刘爱虢

(沈阳航空航天大学 航空航天工程学部(院)，沈阳 110136)

对中心分级燃烧室冷态流场进行大涡模拟，研究主燃级叶片安装角对流场中大尺度湍流结构运动特性的影响，亚格子尺度湍流及其对大尺度流动结构的影响使用动态Smagorinsky涡黏模型描述。结果显示：流场纵向截面上涡量云图揭示了流场上游出现内外两层涡结构，分别发源于文氏管尾缘和套筒末端，随旋流强度增加，涡层径向尺寸变宽，中心回流区长度缩短。主燃级叶片安装角从20°增加到40°时，流场上游角涡显著减少。下游低压核心与流线束环绕中心重合且偏离轴心，并在流场中下游产生进动涡核。中心轴线上的轴向速度功率谱密度表明，随着叶片安装角增加，涡核进动频率升高。

中心分级旋流器；大涡模拟；Q准则；进动涡核

航空发动机燃烧室中经常通过钝体或旋流的方式形成回流区来组织燃烧并稳定火焰，分层旋流燃烧可以提高燃料与空气的混合效率，增强火焰的稳定性和提高燃烧效率。旋流燃烧室中的流动处于强湍流状态，流动发展过程中常出现不稳定现象，如涡破碎和涡进动等[1]，是典型的非定常流动，这种运动行为对燃料混合效率、火焰稳定、燃烧效率、污染物排放等有重大影响。因此研究燃烧室中旋流不稳定现象具有重要的工程应用意义。

通过数值模拟方法获得燃烧室内的流动特性来研究旋流器旋流强度变化对流动的影响得到广泛应用。雷诺平均方法(RANS)对流动的控制方程进行统计平均来建立模型，忽略各尺度的湍流脉动情况，只能得到关于湍流的平均信息。在新出现的大涡模拟(LES)中，根据计算网格尺度，对流场中湍流结构进行过滤。对数值模拟中小于网格尺度的湍流脉动建立模型，进而能捕捉到流场中振荡、涡破碎、涡进动等非定常流动信息。

大涡模拟已经在湍流绝热旋流的研究中得到应用，并建立了多种关于冷态热态旋流场进动和不稳定性模型。Syred N 等人[2]使用Strouhal数和几何旋流数对绝热旋流的进动和不稳定性进行了深入研究。秦皓等人[3]的实验研究了旋流器出口的流动特点和流场对外界激励的响应。张济民等人[4-6]比较全面地研究了悉尼旋流器的流场，捕捉到进动涡核PVC(Precessing Vortex Core)出现在剪切层的下游。当旋流强度增加时，冷态流场回流区长度没有明显变化，涡破碎开始发生的位置向上游移动，进动涡核的运动强度沿流向逐渐减弱。郑韫哲[7]研究燃烧室内的流动及燃料空气混合时，发现涡团破碎程度和速度的增加推动燃料空气混合比较好的区域向上游移动，提高了全流场的燃料空气混合程度。Ranga Dinesh等人[8]研究了4个旋流数对悉尼旋流器流场的影响，捕捉到中心射流进动的低频振荡，叶片安装角增大，涡破碎的区域向下游移动，回流区剪切层的不稳定性增加。朱宇等人[9]的数值模拟结果发现叶片安装角增大回流区长度减小。Shanwu Wang等人[10-12]系统地研究了二级旋流器内部及其出口未受限的流场，比较了两个不同旋流数时流场的差异；两级叶片反向比起同向的情况，流场的中心回流区长度更短，剪切层强度更大，雾化效果更好。

本文对一种中心分级旋流器的冷态流场进行大涡模拟，主要讨论这种旋流器主燃级叶片安装角对燃烧室流场结构的影响。

1 物理模型和计算方法

1.1 中心分级旋流器简介

中心分级旋流器结构如图1所示，结构包括值班级(第1级径向叶片、第2级径向叶片)、主燃级径向叶片、文氏管、套筒。其中值班级包括两级径向旋流器，1级旋流器主要作用是使喷嘴喷出打在其内表面的燃油颗粒进行二次雾化，其叶片安装角为45°，旋向沿进气方向为顺时针；2级旋流器作用是形成稳定回流区，以便稳定燃烧，其叶片安装角为45°，旋向与1级旋流器旋向相反；主燃级处在值班级的外层，空气通过主燃级径向叶片再形成一股旋转气流，满足大工况下燃油所需的空气量，旋向与2级旋流器相同。

图1 旋流器简图

1.2 计算模型

1.2.1 大涡模拟LES(Large Eddy Simulation)

大涡模拟[14]通过建立滤波函数，把尺度比滤波函数尺度小的涡从湍流瞬时运动方程中过滤掉，实际的湍流运动分解成可直接计算的大尺度湍流运动和小尺度脉动。通过建立亚格子应力模型来模化小尺度涡对大涡运动的影响。常用的滤波函数有3种：谱空间低通过滤、高斯过滤和盒式过滤，其中盒式过滤由于其过滤量处在单元中心而更符合有限体积的计算思想，被采用有限体积法的算法所采纳。

本文研究的流动情况进口是常温常压下的空气，进气速度较低，属于不可压缩流动。采用空间盒式过滤器过滤Navier-Stokes方程后，得到以下大尺度流动的方程：

(1)

(2)

1.2.2 涡黏模型

本文采用Smagorinsky涡黏模型模化亚格子尺度流动及其对大尺度湍流的影响，即

(3)

(4)

式(4)中：模型系数CS由局部动态模型[13]计算得到，较好地反映湍流的瞬时和局部动力学性质。

1.3 网格划分

模型计算时取环形燃烧室一个单元段，并以矩形代替扇形，这一近似是可靠的[15]。为消除壁面对旋流器流场的影响，燃烧室取较大的计算区域，截面边长90 mm，轴向200 mm。坐标原点在旋流器上游底圆中心。考虑到要在复杂的旋流器内部几何空间中流动数值模拟来获得接近真实的流动情况，对旋流器的内部空间采用Tet/Hybrid网格；燃烧室是单一的立方体结构，采用Tet网格。最终，网格总数为324万，整体网格模型如图 2所示。

1.4 边界条件定义及计算

在Fluent中求解瞬时不稳定状态时，粘性模型使用Large Eddy Simulation，流场介质为空气，旋流器进口速度给定10 m/s，燃烧室出口采用压力出口。压力-速度耦合求解采用SIMPLE算法，空间离散方法选取二阶迎风格式，并监控中心轴线上Z=90、120、150、180、210、240 mm 6个点的轴向速度。以标准初始化方式从速度进口计算。求解时间步长0.001 s，总共的时间步数2 000，计算时长为2 s的流动情况，保证流场在计算时间内获得充分发展，每个时间步长内最大迭代20次。

图2 划分网格结果

2 计算结果与讨论

2.1 涡结构

由于在顺流和回流分界处存在剪切层，引起回流区的边界产生一层漩涡，所以这个漩涡层指示了回流区边界和剪切层的位置。

图3显示了4个旋流强度下各流场纵切面上的涡量分布云图，旋流器出口扩张平面下游出现的内外两层大尺度涡结构中，大涡量值的涡结构分别发源于文氏管尾缘和套筒末端，这两处的涡是由这些位置的剪切层受到轴向和切向的Kelvin-Helmholtz不稳定性作用而引起，随后这些涡结构被对流到下游，并与周围流体相互作用能量耗散最终消失，而且涡量剧烈衰减的位置随旋流强度增加逐渐向上游收缩。

在流场中上游，涡破碎引起了中心回流区的建立，内涡层内侧低涡量值的涡结构显示的正是中心回流区边缘，这里的漩涡层在流动中翻滚拉伸最终破裂成小的涡泡。强烈的旋流在径向方向上建立了很大的压力梯度，并绕着中心轴线诱导出一个低压区，切向速度随着旋流强度增加，需要更大的低压区来平衡离心力，这就是流场旋流强度增加时中心回流区径向变宽的原因[16]，表现在图中就是内涡层径向尺寸变宽。回流区阻塞效应的存在，减小了燃烧室内有效流动通道的截面积，导致顺流区的流动速度增加，剪切层强度增加。

图3 X-Z平面上瞬时涡量分布

朱宇等人[7]的文章认为中心回流区长度随叶片安装角的增加向上游缩短。结合表 1可知该旋流器的中心回流区长度同样随着角度的增加逐渐缩短。在实际应用中为了使燃烧室出口温度分布满足涡轮叶片对温度要求，避免涡轮前温度分布不合理，造成涡轮叶片根部和叶稍热应力急剧增加，损坏叶片或降低叶片寿命，要求燃烧火焰要最短，燃烧室流场中，在保证燃烧稳定完全的情况下，使回流区长度尽量短。通过以上讨论，增大主燃级叶片安装角是实现要求的方法之一，但旋流强度太强会导致回流区过短使热态流场发生回火。为了达到回流区的风速要求，冷态流场回流区主体要有合适的长度，保证火焰的稳定而且不发生回火。

表1 中心回流区轴向长度

从中心回流区到文氏管内的涡旋在管内也形成强烈的剪切层。当喷嘴满负荷工况运行时，部分燃料会喷射到文氏管内表面形成一层薄膜，在强烈的剪切层作用下，这层薄膜上使其雾化破裂，促进燃料完全雾化燃烧。旋流器出口截面突扩，四周的流体由于卷吸效应的负压形成角落回流区。角落回流区受主燃级叶片角度影响显著，旋流在离心力的作用下快速膨胀向外流动，旋流越强流体膨胀流动越快，再加上外涡层径向尺寸扩大，结果使角落回流区急剧减小。

2.2 拟序结构

2.2.1 涡结构的Q准则计算方法

Q准则是速度梯度张量的二次不变量，定义式是：

Q=(ΩijΩij-SijSij)/2

式中：Sij和Ωij分别是速度梯度张量中的对称张量和反对称张量。二者定义式如下：

Sij=(uij+uji)/2

Ωij=(uij-uji)/2

Q表示转动速率的ΩijΩij大过应变率SijSij的程度。在Q>0的位置，转动速率ΩijΩij占主导地位，即在该区域涡旋结构占主导地位。

分别使用涡核与Q准则来捕捉流场中的涡结构。涡核是与拟序结构密切相关的物理量，其计算方法有两种：涡向量法和速度梯度特征值法，后一种计算方法复杂，但能减少伪涡核得到更真实的涡核分布，本文就选用这种方法来提取。图4上部分是涡核的分布，下部分是涡核与Q准则叠加的情况，取值Q=100 000。图4中显示Q准则等值面很好地包围了涡核，说明这两种识别涡分布的方法有效。主燃级叶片角度越小时，流场中上游部分的涡旋结构越丰富，尤其在角落回流区附近的涡旋明显增多。随着叶片角度增加，角落回流区的涡旋结构显著减少，说明主燃级角度改变对旋流器出口扩张平面处的流场影响显著，不但角落回流减小，涡结构也减少。叶片角度大于等于40°时，比较图4中虚线椭圆所指示的涡结构，随旋流强度增加涡结构逐渐发育规则并向下游生长伸长，而且叶片角度在从20°增加到40°过程中，流场涡结构有重大变化，下游涡结构显著减少。

图4 涡核分布与Q准则等值面

中下游轴线附近的涡不同于上游涡团，它有明显的周期性振荡特性，是进动涡核(PVC)。进动涡旋的典型特征是周围的流体绕着涡核旋转，涡核又另外绕自身的涡轴。Syred文中[2]指出PVC不仅有以上两个特征，而且涡轴处切向速度为零。本文以这些特征来识别中下游的PVC。图5上面显示的是流场中Z=240 mm横截面上的瞬时压强与流线叠加图，下面是切向速度云图。该平面上流线束环绕中心与低压核心均偏离燃烧室轴心且基本重合，实际流动中二者不能完全重合是流体不稳定变形和粘性引起的。而且当叶片安装角大于等于40°时低压核心处的速度云图均显示是零，综合以上特征表明中下游的涡是PVC。

图5 瞬时压力云图与流线叠加图以及切向速度云图

涡核进动会引起周围流场物理量的周期性变化。监测变化的物理量并作做傅里叶分析即可确定附近流场运动是否有特定的周期性。本文对燃烧室中心轴线上Z=90、120、150、180、210、240 mm共6个点的轴向速度做傅里叶变换得到各点的功率谱密度PSD(Power Spectral Density)如图 6所示。主燃级叶片20°时的PSD图表明特征峰存在的轴向空间很小，结合图 5中低压核切向速度不为零的情况得知在20°流场中下游没有进动的涡结构。主燃级叶片角度大于等于40°时，流场从Z=150mm的位置开始出现各自的特征峰，各流场特征峰所对应的特征频率依次是f=31、35、38 Hz。在Shanwu Wang等[11]的文中，旋流强度加大时，对压力波动监测的频谱图上显示低频成分显著增加。而在本文的旋流器几何结构条件下，主燃级叶片角度增加对下游涡结构的进动产生的影响是使进动频率上升。而图中的低频振荡可能是由于顺流与回流之间的涡层破碎成小涡泡引起的。

图6 燃烧室中心轴线上6个点的轴向速度功率谱密度(PSD)

特征峰强度在流场中出现与衰减消失的位置指示了PVC的轴向长度，安装角40°时PVC终结于燃烧室内部，60°、80°时则一直影响到燃烧室出口。 进动会影响流场的特性，一定强度的涡结构进动能促进燃料空气以及热态流场中高活化能组分混合，加强燃烧器中大尺度湍流可以提高燃烧效率，热态流场中一定频率的进动可能与燃烧室内的声波产生耦合共振，改变火焰结构影响燃烧室的释热行为，造成燃烧不稳定。流场内的有规则行为也可能会与机械结构发生共振减损发动机的寿命。在实际设计燃烧室时，为保证燃烧的稳定性，应尽量避免流场中PVC的产生。

3 结论

本文采用大涡模拟的涡黏模型对中心分级旋流器的冷态流场进行了数值模拟，研究了主燃级叶片安装角对流场结构的影响。

(1)流场纵向截面上涡量分布揭示了流场中内外两层高能量涡分别发源于文氏管尾缘和套筒后缘，内外两涡层随旋流强度增加而径向变宽，中心回流区轴向长度减小。外涡层向燃烧室上游扩张平面靠近，角落回流区受到越来越明显的抑制。

(2)涡核与Q=100 000的等值面捕捉到流场中的大尺度涡结构。显示主燃级叶片安装角从20°增加到40°的过程中流场发生了显著变化，中下游大尺度湍流结构减少且开始组织出规则的流动结构，扩张平面附近的涡旋结构剧烈减少。

(3)叶片安装角度大于等于40°时，Z=240 mm平面上的流线图、压力云图、切向速度云图中流线环绕中心、低压核心、零切向速度位置重合且偏离几何中心，说明中下游出现了进动涡核。

(4)轴向速度的功率谱密度图表明进动涡核从Z=150 mm的位置处开始出现，且进动频率随主燃级旋流强度增加而升高。

[1]LUCCA-NEGRO O,O′DOHERTY T.Vortex breakdown:a review[J].Progress in Energy and Combustion Science,2001,27(4):431-481.

[2]SYRED N,BEER JM.Combustion in swirling flows:a review[J].Combust Flame,2009，23(2)：143-201.

[3]秦皓，丁志磊，林宇震，等.同心分级旋流结构的动态响应特征[J].航空动力学报,2015,30(4):793-799.

[4]张济民，张宏达，韩超，等.分层旋流燃烧器冷态流场的大涡模拟.[J]航空动力学报,2015,29(10):2369-2376.

[5]张济民，韩超，张宏达，等.钝体绕流中回流区与进动涡核的大涡模拟.[J]推进技术,2014,35(8):1070-1079.

[6]张宏达，张济民，韩超，等.大涡模拟钝体有旋流流场的拟序结构.[J]航空学报,2014,35(7):1854-1864.

[7]郑韫哲,朱民,姜羲.旋流预混燃烧室流动混合的大涡模拟[J].航空动力学报,2012,27(1):33-40.

[8]K K J RANGA DINESH,M P KIRKPATRICK.Study of jet precession,recirculation and vortex breakdown in turbulent swirling jets using LES[J].Computers & Fluids，2009，38(5)：1232-1242.

[9]朱宇，张群，徐华胜，等.轴向旋流器几何对回流区的影响.[J]航空动力学报,2014,29(11):2684-2693.

[10]SHANWU WANG,VIGOR YANG,GEORGE HSIAO,et al.Large-eddy simulations of gas-turbine swirl injector flow dynamics[J].Fluid Mech,2007,583(9):99-122.

[11]SHANWU WANG,SHIH-YANG HSIEH,VIGOR YANG.Numerical simulation of gas turbine swirl-stabilized injector dynamics[R].AIAA-2001-0334,2001.

[12]SHANWU WANG,VIGOR YANG.Modeling of gas turbine swirl cup dynamics,part 5:Large eddy simulation of cold flow[R].AIAA-2003-485,2003.

[13]YING HUANG,VIGOR YANG.Swirling flow structures and flame characteristics in a lean-premixed combustor[R].AIAA-2004-978,2004.

[14]张兆顺，崔桂香，许春晓.湍流大涡数值模拟的理论和应用[M].北京：清华大学出版社，2008.

[15]JONG-CHAN KIM,HONG-GYE SUNG,DAI-KI MIN,et al.Large eddy simulation of the turbulent flow field in a swirl stabilized annular combustor[R].AIAA-2009-645,2009.

[16]KIRAN MANOHARAN,SAMUEL HANSFORD,JACQUELINE O′CONNOR,et al.Instability Mechanism in a swirl flow combustor precession of vortex core and influence of density gradient[C].Proceedings of ASME Turbo Expo 2015:Turbine Technical Conference and Exposition,GT2015-42985.

(责任编辑：吴萍 英文审校：赵欢)

Large-eddy simulation of non-reacting flow in a central staged combustor

YU Lei,WANG Cheng-jun,ZHAO Dong-kai,LIU Ai-guo

(Faculty of Aerospace Engineering,Shenyang Aerospace University,Shenyang 110136,China)

Large-eddy simulation (LES) was used to explore non-reacting flow fields of the central staged combustor and investigate effects of the vane angle of swirler primary stage on the large-scale flow structures.Smagorinsky eddy viscosity model was used to describe sub-grid turbulence flow and its influences on the large-scale flow structures.The contours of transient vorticity magnitude on x-z plane show that there are two vortex layers in upstream of the flow field.The inside layer derives from the end of the venturi and the outside one from the outer tube tail.With the increase of swirl strength,radial distance of the two vortex layers becomes larger,and the axial length of center recirculation zone is reduced.The upstream corner vortex is significantly decreased when the vane angle increases from 20° to 40°.In downstream,the position of low pressure core and the streamline surrounded center overlaps.The overlap region deviates from the axis of chamber,and precessing vortex core generates in downstream of flow field.The axial velocity power spectrum density on the chamber axis shows that precessing frequency of the vortex rises with the increase of angle.

central staged swirler;large-eddy simulation;Q criterion;precessing vortex core

2017-04-26

国家自然科学基金(项目编号：51476106)；辽宁省自然科学基金(项目编号：2015020639)

于雷(1993 - )，男，辽宁抚顺人，硕士研究生，主要研究方向：航空发动机燃烧技术，E-mail：1243892925@qq.com；王成军(1967- )，男，辽宁沈阳人，副教授，博士，主要研究方向：航空发动机燃烧技术，E-mail:wangchengjun22@sina.com。

2095-1248(2017)03-0043-07

V231.2

A

10.3969/j.issn.2095-1248.2017.03.006