随机场TPDF在非结构网格平台的开发及检验

2023-05-13杨兴峰许笑颜高浩卜

航空发动机 2023年1期

王方，杨兴峰，方存，许笑颜，莫毅，高浩卜

（1.北京航空航天大学能源与动力工程学院，北京 100191；2.北京航空航天大学江西研究院，南昌 330096；3.北京航空航天大学成都航空动力创新研究院，成都 611930；4.中国航空发动机研究院，北京 101300）

0 引言

航空发动机燃烧室几何结构非常复杂，先进航空发动机设计对燃烧室提出更精确的燃烧组织要求，对燃烧室数值模拟提出了更高的精度要求[1]，对湍流燃烧模型也提出了更高的要求[2]。概率密度函数输运方程湍流燃烧模型（Transport Probability Density Function Turbulence Combustion Model，TPDF）可以耦合详细化学反应机理，高精度解析湍流燃烧过程[3]。

概率密度函数输运方程湍流燃烧模型[4]主要分为3种：速度-标量联合概率密度函数输运方程；标量联合概率密度函数输运方程；速度、标量和耗散率三者联合的概率密度函数输运方程。目前，3 种TPDF中应用最广泛的是标量联合的概率密度函数输运方程。在标量联合概率密度函数输运方程解法中，欧拉随机场法是一种概率密度函数的求解方法[5]。Jones等[5-7]采用射流、钝体、旋流、燃烧室等算例检验了欧拉随机场解法的TPDF 应用；王方等[8-9]在BOFFIN 程序的基础上，发展了适用于航空发动机两相湍流燃烧的AECSC 软件，并针对化学反应机理、两相蒸发模型等进行了研究；Gong 等[10]和Rihab 等[11]合作，在开源软件OpenFOAM 平台上开发了TPDF 燃烧模型代码，求解声速问题。TPDF可望成为燃烧室高精度模拟的基础模型之一[12]。

射流火焰和旋流火焰是燃烧室湍流燃烧情况的典型代表。Flame D 是非预混射流火焰，测量数据多且可靠，是目前湍流火焰预测方法检验的必经之路。比如Nik 等[13]检验速度标量联合滤波密度函数的方法；Renzo 等[14]验证FPV 湍流燃烧模型；Lysonko 等[15]检验涡耗散模型；Ge 等[16]检验稀疏拉格朗日MMC 模型。此外，还有学者用Flame D 研究辐射模型[17]影响、RCCE 应用[18]等；TECFLAM 是Dreizler 等[19]开发的甲烷/空气非受限强旋流贫油预混燃烧器；Freitag 等[20-22]利用直接数值模拟和大涡模拟（Large Eddy Simulation，LES）对TECFLAM 冷态和热态工况进行研究，并利用关联函数的方法对冷态流动工况中涡团尺度进行了分析；Wegner等[23]利用非稳态雷诺应力模拟（Unsteady Reynolds Averaged Navier-Stokes，URANS）和LES 对TECFLAM 冷态工况进行了模拟，分析了URANS 在模拟旋进涡核结构的适用性，并给出了URANS和LES 处理旋进频率的不同方法；隋春杰[24]用LES 和改进代数2 阶矩湍流燃烧模型，结合甲烷2 步化学反应机理，成功模拟了TECFLAM 冷态和热态工况，指出入口当量比扰动对预混燃烧模拟的影响；在Rolls-Royce 公司对PRECISE-UNS 软件发展测试工作中，对非结构网格下的火焰面类、TPDF 等模型进行了射流、旋流和燃烧室测试研究[12]。

Saturne 软件是开源软件，能够应用非结构网格。在Saturne 平台构建TPDF 并对射流、旋流、燃烧室进行模拟测试，可以检验TPDF 对航空发动机燃烧室的适用性，判断能否高精度求解复杂结构燃烧室内部燃烧过程。本文在非结构网格体系中研发随机场TPDF方法及配合的计算单元，测试模拟效果，以促进航空发动机燃烧室模拟技术难题的解决。

1 非结构网格平台、TPDF开发和算例介绍

1.1 非结构网格平台

Saturne 是由法国电力公司研究院（EDF R&D）设计开发的通用CFD 软件，可求解3 维Navier-Stokes 方程。Saturne软件基于有限体积离散，使用Fortran和C语言作为主要编程语言，支持非结构化网格并行计算。其涵盖包括湍流模型在内的多个物理模块，并且可以兼容目前大多数类型的网格文件，如cgns、gambit等格式，具有积极的工程应用价值。

1.2 TPDF简介

详细的欧拉随机场的TPDF 介绍可以参考文献[5]和[25]。考虑密度加权平均的标量联合概率密度函数输运方程为

式中：˜P为Favre 滤波后的概率密度函数；ρˉ为空间滤波后的密度；˜uj为Favre 滤波后的j方向速度；ψ为标量；α为组分；F为联合概率密度函数；ω̇α为组分α的化学反应源项；Ns为组分数量；ϕα为与α有关的标量；ψα为ϕα构成的样本空间；ψ为ψα构成的点的集合；μ为动力粘度；σ为施密特数。

式中左侧3 项分别为时间项、对流项和化学反应项，均是可以精确计算的。右侧2 项则需要通过梯度扩散模型和小尺度混合模型进行模化，得到最终封闭形式的TPDF方程为

式中：μT为湍流粘度；σT为湍流施密特数；τ为湍流混合时间尺度；ϕ͂α为Favre滤波后的ϕα。

设有N个随机场的集合ξm(x,t)=[ξm1,…,ξmN]，该随机场对于所有标量都是欧拉形式的。采用随机场解法的Ito形式求解TPDF方程，得到最终可以应用到程序中的TPDF方程形式为

1.3 基于Saturne软件的TPDF研发

TPDF 构建框架如图1 所示。从图中可见，仅用红框标出的文件为修改的文件，主要目的是实现TPDF的独立；全部标红的文件为新增的文件，这是湍流燃烧模型的主要求解文件，用于增加组分输运方程和能量方程以及数据的传输和其他功能构建。本文在Saturne 开发的TPDF 构建了质量分数输运场和总焓输运场，另外创建了一些不参与输运方程求解的包括物种生成率、温度等在内的属性场。主要求解文件是pdfini、pdftss、pdfphy、pdftcl，分别用来实现模型初始化、源项求解、物理属性更新、边条处理。

图1 TPDF构建框架

1.4 算例介绍和边界条件

1.4.1 Flame D算例

Flame D 的几何结构见文献[23、26]。FUEL 区为主喷射口，直径d=7.2 mm，喷射工质为体积分数25%的甲烷和75%的空气组成的混气。PILOT 区为高温值班火焰，主要成分是已燃气，具体的各物质质量分数以及不同轴向位置截面的试验数据可以在Sandia实验室官方网站[27]查到。Flame D 算例边界条件见表1。

表1 Flame D算例边界条件

1.4.2 TECFLAM算例

TECFLAM 燃烧器几何结构如图2 所示。甲烷和空气先在燃烧器进口位置充分掺混，然后经旋流叶片进入主要的燃烧区域。本文选定的热态试验工况为30 kW的热负荷下的试验工况，旋流数为0.75，甲烷空气当量比为0.83，预混气总体积流量为37.92 m3/s[28]。冷态和热态工况的旋流和伴流进口边界条件见表2、3，其余所有壁面均为无滑移绝热壁面条件，出口为压力出口条件。

图2 TECFLAM 燃烧器几何结构[6]

表2 TECFLAM 冷态算例进口边界条件

表3 TECFLAM 热态算例进口边界条件

1.5 计算网格和化学反应机理

用ICEM 软件绘制网格划分O-BLOCK，Flame D网格轴向长度0.54 m，径向长度0.3 m，总网格单元数目约为30 万，对进口和中心射流区进行加密，Flame D 算例网格如图3 所示。TECFLAM 算例采用文献[24]的网格划分策略，TECFLAM 算例网格如图4 所示。文献[22]采用长度2D的环形通道使流动发展，计算区域总轴向长度为10D。D代表TECFLAM 中心钝体的直径30 mm。最终的网格单元总数目约为80万。本文选用Zhou 等[29]的甲烷单步化学反应机理，在mol-cm-g-s 单位制标准Arrhenius 表达式中，A=3.2×1014cm3/（mol·s），b=0.0，E=36560 kcal/mol，甲烷和氧气的反应级数均为1.0。

图3 Flame D算例网格

图4 TECFLAM 算例网格

2 模拟结果与讨论

2.1 Flame D算例

RANS-TPDF 模拟的Flame D 时均温度如图5 所示。Flame D 火焰核心射流进口流速较高，以喷嘴出口速度为特征速度的雷诺数超过20000，属于典型高雷诺数部分预混射流火焰。其PILOT 区提供了1880 K 的高温，甲烷被高温点燃发生化学反应释放出大量热，瞬时高温可达2200～2300 K。随着射流发展甲烷不断被消耗，到中部截面处高温区汇合，形成“蜡烛形”火焰。之后，高温气体继续扩散。

图5 Flame D时均温度

分别将模拟平均温度和甲烷浓度与试验数据进行对比，Flame D 不同轴向位置时均温度径向分布如图6 所示，Flame D 不同轴向位置时均甲烷质量分数径向分布如图7 所示。从图6、7 中可见，模拟结果总体与试验值有较好一致性。前3 个截面是燃烧发展区，流动主导而反应刚开始，吻合较好；中间3 个截面化学反应起主要作用，放热较为集中且有大量热量累积。在测试中不同化学反应机理有不同的温度分布和最高温度，Saturne 单步机理与多步机理时均温度结果对比如图8 所示，多步化学反应机理预测结果明显优于单步化学反应机理的。RANS 模拟对流动预测存在较大误差，单步甲烷机理放热快速，二者在中间截面有相当大的误差。下游3 个截面主要受对流扩散影响，吻合较好。

图6 Flame D不同轴向位置时均温度径向分布

图7 Flame D不同轴向位置时均甲烷质量分数径向分布

图8 Saturne单步机理与多步机理时均温度结果对比

进一步与相同条件下Fluent 中涡破碎（Eddy Break-Up，EBU）模型、E-A（EBU-Arrehnius）模型和有限速率（Finite-Rate，F-R）模型计算结果对比，不同燃烧模型在不同轴向位置的时均温度径向分布如图9 所示，可见新软件精度与E-A 和F-R 湍流燃烧模型相近，优于EBU湍流燃烧模型的。

图9 不同燃烧模型在不同轴向位置的时均温度径向分布

2.2 TECFLAM 冷态工况研究

URANS 瞬时轴向速度和时均轴向速度如图10所示。瞬时轴向速度在不停地波动，时均轴向速度存在因强旋产生的中心回流区以及在喷嘴出口壁面附近因台阶效应产生的小回流区。

图10 URANS瞬时轴向速度和时均轴向速度

TECFLAM 冷态p=1.5 Pa 压力等值面如图11 所示，显示为旋进涡核结构（Precessing Vortex Core，PVC）。对于PVC 旋进频率，可在试验和LES 中由湍流能谱[20]通过时间自相关得到，也可在URANS 中通过对固定点的速度时间序列进行傅里叶变换得到[19]。根据文献和测试结果，选定距喷嘴出口中心轴向30 mm 和径向20 mm 为固定点。

图11 TECFLAM 冷态p =1.5 Pa压力等值面

点（30 mm，20 mm）瞬态轴向速度时域与频域如图12 所示。从图中可见，PVC 结构稳定，旋进周期约为0.0225 s，旋进频率约为45 Hz。Freitag 等[20]用LES预测的旋进频率约为42 Hz；Wegne 等[23]用LES 预测的旋进频率为50 Hz，URANS 为35 Hz；隋春杰等[24]用LES预测约为40 Hz。本文结果与前人研究一致。

图12 点（30 mm，20 mm）瞬态轴向速度时域与频域

TECFLAM 冷态4 个轴向位置的时均速度径向分布如图13所示。从图中可见，对比距喷嘴出口1、30、60、90 mm 截面的时均轴向速度，模拟与试验结果在峰值大小和整体趋势上接近。说明新软件对复杂旋流的适用性好，为热态模拟奠定了基础。

图13 TECFLAM 冷态4个轴向位置的时均速度径向分布

2.3 TECFLAM 热态工况研究

TECFLAM 热态时均温度如图14 所示。从图中可见，高温区主要在强旋流产生的中心回流区内，稳定在中心钝体下游。中心回流区将高温燃气从下游卷回，结合中心钝体，给甲烷和空气提供了1 个稳定、低速、高温的稳定燃烧区域。从上游向下游，高温区变宽后温度不断降低，宽度渐渐变窄。从中心高温区到常温空气伴流区，有非常明显的温度梯度下降。

图14 TECFLAM 热态时均温度

TECFLAM 热态4 个轴向位置的时均速度径向分布如图15 所示。从图中可见，模拟的最大轴向速度值和整体趋势与试验值较为相符，特别是在上游截面，这与射流模拟结果一致。热态工况和冷态相比，时均轴向速度正向最大值和负向最大值都有了提高，说明燃烧的存在使局部流场温度提高密度减小，对局部气流产生了加速效果。下游2 个截面的速度误差一部分来源于与射流类似的情况，即RANS 模拟方法和简单反应机理；另一部分来源于热态旋流各向异性的输运特征。热态工况远离中心区域，在径向大于30 mm 的位置，模拟值与试验值有比较好的贴合。在化学反应剧烈和旋流强的区域，湍流模型、湍流燃烧模型和化学反应机理都受到挑战。

图15 TECFLAM 热态4个轴向位置的时均速度径向分布

TECFLAM 热态4 个轴向位置的时均速度径向分布与文献[24]中改进代数2 阶矩模型（Modified Algebraic Second Order Moment，MASOM）和部分搅拌反应器模型（Partially Stirred Reactor，PaSR），文献[30]火焰面生成流型模型（Flamelet Generated Manifolds，FGM）的结果对比如图16 所示。本文模拟的时均温度与试验值在最高温和整体趋势上都保持了较好的一致性。Saturne 模拟结果从中心火焰区到低温伴流区温度降低滞后且降低速度过快，说明模拟的火焰面结构过薄，与试验存在差异。本文基于RANS 的模拟结果优于文献结果，特别是一些基于LES 的结果，说明TPDF的优越性。

图16 TECFLAM 热态4个轴向位置的时均速度径向分布与文献[24]、[30]对比

3 真实燃烧室算例测试

基于一致性等研究基础[31]，应用新软件模拟某型两相燃烧室，采用377万非结构化网格，如图17所示。空气速度进口为121.07 m/s、温度为696.01 K，出口为充分发展条件，壁面绝热无滑移。初始压力为790216.58 Pa，温度为696.01 K。中心截面速度、温度和出口温度分布如图18～20 所示。从图中可见，出口温度分布的预测数值和试验测量得到的数值比较接近，相对误差不超过3.8%。