基于动态分区概念的高超声速燃烧大涡模拟1)

姚卫 刘 杭 张 政 肖雅彬 岳连捷,†

* (中国科学院力学研究所高温气体动力学国家重点实验室,北京 100190)

† (中国科学院大学工程科学学院,北京 100049)

引言

2004 年,X-43 A 在35 km 高度的临近空间实现了接近马赫数10 的高超声速飞行[1],取得了大气层内吸气式飞行的最快速度.近年来,高马赫数飞行器引起了世界各国的高度重视,成为各空天科技强国竞相争夺的科技制高点.飞行马赫数上限达到12 甚至更高的超燃冲压发动机是实现可重复使用天地往返飞行器的关键[2].目前反射激波风洞和激波膨胀风洞是能够产生马赫数大于8 高焓气流的两种主要途径.由于高焓气流地面模拟的困难,目前世界上仅有为数不多的高马赫数风洞,例如美国的NASA HYPULSE 激波风洞[3]、日本的JAXA HIEST 自由活塞驱动激波风洞[4]和澳大利亚昆士兰大学T4 反射激波风洞[5].其中在T4 风洞中针对飞行马赫数7.5～ 12[6-12]工况开展了一系列的矩形转椭圆(REST)高超声速燃烧室缩比模型测试.由于高超声速实验测量手段和可获得数据均十分有限,计算流体力学(CFD)成为燃烧内部流场分析和发动机性能优化的必需手段.目前文献中已经开展了马赫数8,12,13 和15[13-15]等对应飞行条件下的高超声速超燃冲压发动机模拟,极大弥补了实验测量信息的不足.然而,随着飞行马赫数的提高,湍流强度增加,流动和化学反应特征时间相互重叠区域增加.同时,为了增加流动滞留时间以提高混合和燃烧效率,相比于马赫数4～ 8的常规超燃冲压发动机,马赫数大于8 的高马赫数超燃冲压发动机流道尺寸大幅增加,意味着巨大的全尺寸模拟计算代价.高马赫数条件下内流与外流的耦合效应也显著增强,为了提高模拟置信度一般计算区域需要至少包含部分飞行器前体.上述独特性意味着高马赫数超燃冲压发动机数值模拟在兼顾模型保真度(复杂多重湍流-化学反应交互作用模式的表征)与计算代价(耦合详细化学反应动力学的高解析度大涡模拟)两方面均面临着极大的挑战.

空间上复杂而时间上多变的湍流化学反应交互作用关系(turbulence-chemistry interaction,TCI)是超声速燃烧室内部湍流燃烧的主要特点.在超燃冲压发动机中,其内部流场不仅存在混合控制为主的快速反应区和化学反应控制为主的慢速反应区,甚至还存在预混、非预混和部分预混等复杂的湍流燃烧模式[16-20].跨越多个时间和空间尺度的湍流和化学反应不仅自身的模拟极为困难,而且其交互作用所引起的复杂耦合效应进一步加剧了其模拟挑战性.湍流燃烧模拟的核心在于对复杂的湍流-化学反应交互作用关系准确、高效建模.目前的湍流燃烧模型大致可分为两类:有限速率类和守恒标量类.有限速率类模型包括近似层流模型(QL)[21]、部分搅拌反应器模型(PaSR)[22]、涡耗散/破碎模型(EDC/EBU)[23]、增厚火焰面模型[24]和线性涡模型(LEM)[25]等.守恒标量类模型包括各类火焰面模型[26-28]、条件矩封闭模型(CMC)[29]、假定[30-31]/输运[32]概率密度函数模型等.有限速率类模型适用范围较广但一般需要在每个网格节点上求解所有的组分输运方程,因而计算量较大.守恒标量类模型通过对流动和化学反应的解耦操作实现了湍流燃烧的降维求解,一般可以显著降低计算量.然而,流动-化学反应可以解耦的前提是化学反应受湍流的影响模式单一,以便选择适合的守恒变量-热力学参数空间分布去表征当前化学反应状态.以火焰面模型为例,当湍流化学反应交互处于火焰面或薄反应区模式(戴姆克拉数Da＞10,卡尔洛维奇数Ka＜100)时,适宜采用计算简便的稳态火焰面模型;而当湍流化学反应模式处于慢反应区模式(Ka＞100,Da＜10)时,此时化学反应速率小于湍流混合速率,需要采用非稳态火焰面(UFM)或火焰面/进度变量(FPV)模型[33];在发动机燃烧室中还经常存在局部预混或部分预混模式,此时需要采用火焰面模型的其他变种如G 方程模型[34]或火焰面生成流形模型(FGM)[35]等.

与亚声速流动相比,超声速流动有如下特性:(1)滞留时间极短(毫秒量级),和碳氢燃料的化学反应特征时间相当;(2)雷诺数一般高于5 × 105,最小湍流尺度远小于火焰面厚度.这意味着超声速燃烧中湍流和化学反应因特征时间和特征尺度相当,存在强烈的耦合效应.超声速燃烧场在不同局部区域存在差异性极大的湍流-化学反应交互作用模式,反应状态参数对守恒标量的依赖存在极大的空间不均匀性,因此采用单一火焰面去描述整场湍流反应状态将会带来严重的误差,甚至一般认为火焰面模型不适用于超声速燃烧模拟.尽管目前文献中已经发展了基于多重火焰面的表征火焰面模型[36],但该模型仍然把整场作为一个火焰面分区,并不能很好地区分受当地湍流影响的化学反应状态.而如果能够根据局部流场状态自适应地采用单独的守恒标量反应状态参数空间去表征,即局部表征火焰面模型的概念,则可极大地提高整场计算的准确性.局部火焰面模型的理论依据是,尽管反应状态参数对守恒标量的依赖存在空间不均匀性,但针对局部流场可以假设该依赖关系具有统计均一性.分区模拟从物理本质上提高了火焰面模型在超声速燃烧模拟中的适用性,从而为高效高保真超声速燃烧大涡模拟提供了一个独特的解决思路.从针对各局部分区分别建立守恒标量与反应状态参数之间关联的角度看,分区湍流燃烧模型与条件矩封闭模型类似[29].但是分区湍流燃烧模型的分区是基于局部流场的湍流-化学反应交互作用模式并且随流场演化而动态调整,而条件矩封闭模型一般是单纯基于几何坐标位置(多为笛卡尔正交坐标)的固定分区.对于瞬态流场而言,局部湍流化学反应交互作用模式是动态变化的.在分区火焰面模型中,各分区的守恒标量-状态参数空间分布由一个表征火焰面(RIF)[37]描述.表征火焰面描述的并非实际物理上的火焰面结构,而是基于统计规律建立的近似映射关系.因此该分区表征火焰面模型并不局限于传统的薄火焰面假设,而是普遍适用于非平衡、瞬态和非均匀反应过程的模拟.分区湍流燃烧模拟的前提是分区内湍流化学反应交互作用模式局部近似均匀,不合适的分区会导致表征火焰面与实际的状态参数在守恒标量空间内分布差别较大.为了准确描述局部化学反应动力学状态,需要自适应地进行分区动态调整,以便使得各分区内化学反应受湍流的影响模式局部单一.

近年来,分区类燃烧模拟方法在发动机燃烧模拟中得到了较多关注[38-45].目前的分区燃烧模拟主要是依据化学反应热力学状态对燃烧场进行分区或聚类,以便采用简化化学反应计算模型(如冻结反应、均匀混合和平衡态假设等)或基于分区对化学反应进行整体求解.简单的分区一般基于几何空间分为未燃区、反应区/火焰锋面区和已燃区等[38-42,45],复杂的分区对具有相似反应状态的单元(可以几何不相邻)进行聚类操作[43-44,46-48].然而现有分区燃烧模拟并未考虑湍流效应,确切地说应该称之为“化学分区”.已经有部分研究[49-50]在多燃烧模型耦合计算方面开展了初步的探索,但是其所采用的方法需要提前初步评估各模型性能以指定分区,难以在后续计算中动态调整分区,因此不适用于瞬态性较强的超声速燃烧流场.对各模型结合当前算例进行预测评估也需要额外的计算代价,一定程度上降低了多模型耦合计算的效率优势.此外,有限速率类模型因计算代价过大,难以耦合详细的化学反应机理,与火焰面模型耦合使用时会产生“短板效应”,严重降低整体计算效率,因此具有不同理论基础的湍流燃烧模型分区耦合应用价值有限.

本文针对超声速燃烧模拟提出了一种基于当地湍流化学反应交互作用关系进行自适应流场分区和基于局部表征火焰面的分区湍流燃烧模拟方法.该方法不仅可以极大提高存在复杂多重湍流燃烧模式的超声速燃烧模拟的保真度,而且通过对局部湍流化学反应的解耦可取得极高的计算效率,从而突破长期以来制约超声速燃烧模拟的流动-化学反应多尺度作用耦合建模这一关键瓶颈.研究进一步结合该方法对不同构型的高超声速超燃冲压发动机进行了全尺寸内外流耦合一体化大涡模拟,并对其中的流动、混合和燃烧机理进行了深入分析,对比了两种构型高超声速燃烧室的性能,展示了动态分区湍流燃烧模拟方法在超声速燃烧模拟基础研究和超声速燃烧室优化设计方面的应用.

1 物理模型与数值方法

1.1 动态分区火焰面模型(DZFM)

定义当前局部分区内的表征火焰面结构函数为Qα=〈Yα|ξ(x,t)=η,x∈zone〉,其中 η 为混合分数空间的取样变量,x代表空间物理坐标,x∈zone 表示条件平均统计局限于当前区域.瞬态质量分数Yα与表征质量分数Qα的关联为

考虑相变效应和组分间差异扩散的描述质量守恒、瞬态混合分数 ξ 和组分质量分数Yα的完整控制方程如下

其中 ρ 为密度,U为速度向量,Dξ和Dα分别代表混合物、单独组分的扩散速率,Wα为化学反应速率,为相变/凝固等引起的质量增添,Yl,α为其他相(如液滴)中的组分 α 的质量分数,ξl为其他相的整体混合分数,对于纯燃料液滴可取 ξl=1.

将式(1)差分后代入式(4),并使用式(2)和式(3),可得

对式(5)取条件平均 ξ (x,t)=η,和位于当前分区内x∈zone,可得表征火焰面条件组分Qα的最终控制方程为

Evap表示相变效应贡献项,在不涉及两相相变的情况下为零.χ是标量耗散率定义为(∇ξ)2,ρη=〈ρ|η〉,分区条件扩散分区条件平均的引入意味着在每个分区内统计上遵循式(6)定义的守恒性.式(6)中左端第二项表示局部小火焰在相邻区域之间流动的对流输运,其物理意义为使得下游小火焰继承上游小火焰的反应状态,这是准确模拟点火过程和火焰抬升现象的关键.式(6)右端第1 项和第2 项分别模拟了因标量耗散和组分间差异扩散引起的混合分数空间内组分扩散,第3 项表示分区内组分因化学反应而发生的演化.

分区内的化学反应受分区条件平均温度控制,在本模型中采用基于当前分区内瞬态焓统计的方法获得[51].通过忽略焓脉动 (H′～O(0))[52-54],焓的概率密度分布 〈H|η〉 可以假设为以当前Favre 平均值为中心的狄拉克 δ 函数.因此条件焓的 〈H|η〉 分布可以基于历史统计方法获得

其中n为符合=η 的取样空间数据点.Thornber等[55]开展的高分辨率大涡模拟计算中采用了类似的统计方法以获取燃烧模型所需要的条件平均量.条件温度通过条件统计平均焓与混合分数空间组分信息获得,QT=f(〈H|η〉,Qα) .该统计模型极大节省了直接求解条件能量方程的计算成本,以及对超声速可压缩条件下条件能量源项建模的复杂度.该统计条件温度模型的适用前提是分区内包含足够多的时间或空间采样样本,以及数值模拟能够捕捉脉动瞬态值,因此统计模型方法仅适用于大涡模拟(LES)和直接数值模拟(DNS)等高解析度模拟方法[55].

上式表明式 (8) 右端第1 项eY的作用为在每个分区对应的采样空间内重新分布条件平均量,并且统计上重布效果相互抵消,分区内整体质量分数脉动守恒.通过动态更新火焰面分区使得每个火焰面分区对应一个较窄范围的混合分数子空间η ∈.其中为当前分区平均混合分数,Δ ξ 为瞬态混合分数的脉动范围.通过进一步细化分区,当前分区混合分数变化幅度趋近于零 Δ ξ →0,分区内标量概率密度分布坍缩为一个以为中心的狄拉克 δ 函数.据此有,这意味着式 (8) 中的条件脉动输运项通过动态自适应分区可以被忽略.

DZFM 模型在四维空间(三维物理空间加一维混合分数空间)内求解火焰面的演化.传统的火焰面模型可以看作是一个把整个计算区域当作一个单一分区的特殊分区火焰面模型,属于低维流型模型(LDMM) 的一种.全组分输运湍流燃烧模型(如PaSR[22])需在三维物理空间离散网格单元逐一求解组分输运方程和化学反应,计算代价一般远高于火焰面等流动-化学反应解耦模型.DZFM 模型实现了流动和化学反应在当前分区内的局部解耦,因而只需在火焰面分区单元上求解火焰面结构函数随空间的演化.由于火焰面分区的空间分辨率通常远小于CFD 网格数目,所需求解的化学反应体系大幅减少,因此整体上提高了计算效率.注意DZFM 模型的火焰面分区需随着流场而动态演化,以满足湍流化学反应交互模式局部单一的前提假设.

1.2 分区自适应化学(Z-DAC)

在流场求解过程中动态地依据当地热化学反应状态进行实时机理简化的方法称之为动态自适应化学(DAC)[56-57].自适应化学方法通过发展并即时应用适用于当前局部和瞬态流场状态的简化机理可以在保证化学反应模拟准确性的前提下提高化学反应求解效率.超声速燃烧中为了准确捕捉自点火和稳焰模态等流动-化学反应强耦合瞬态特征,一般需要采用尽可能详细的化学反应机理.对于需要直接积分(DI)的化学反应流求解,其计算代价近似与组分数目的平方成正比[58].限于当前的计算技术水平,目前超声速燃烧大涡模拟或直接数值模拟可接受的反应组分数目一般最多不超过50 个组分[59].过度简化的机理往往降低了机理的适用范围,而适用范围较广的机理尺寸又因尺寸过大无法应用于高解析度的三维反应流模拟.此外,流场中存在复杂而多变的热化学环境,难以提前针对所有环境发展准确的简化反应机理[56].自适应化学方法通过冻结对目标组分和特性关联较小的反应路径和组分实现了在不同时间和空间流场中应用不同的反应机理,其局部所应用的反应机理通常为更详细母反应机理的一个子集.由于局部子机理的应用目标热化学环境变窄,因此更容易在较少组分和反应步的情况下实现高化学反应保真度.这里局部子机理既可以采用当前局部分区状态进行实时简化的方式,也可以采用预估状态进行预先简化建库供后续查询的方式.

自适应化学准确应用的关键是选取合适的分区.以往的研究中一般采用根据空间坐标分区的方法甚至直接采用并行分区的方式.这种方式无法保证当前分区的热化学反应状态是均一的,所简化的局部子机理甚至无法适用于当前分区内所有CFD节点状态.对于采用并行分区的方法,虽然建模框架搭建较为简便,但是同样无法保证当前分区内热化学状态的均一性,此外分区方式(如剖分算法、剖分数目)的改变也会影响计算结果.

本文提出了一种基于当地自适应分区的动态自适应化学方法(zonal dynamic adaptive chemsitry,ZDAC).该方法依据流场的反应状态变量(如混合分数和温度等)将计算区域划分为若干分区,每个分区内的热化学反应状态相对均一.在本文中选取了混合分数来表征当前组分的变化,用温度来表征当前分区的化学反应活跃性.如果包含对压力敏感的化学反应,可进一步拓展加入压力指标体现压力不均匀性对反应路径的影响.每个分区均对应某一特定的混合分数子空间、较窄的温度子空间以及压力子空间,因此可认为分区内的热化学反应状态均匀,从而为自适应子机理的匹配性与适用性提供了保障.由于混合分数、温度和压力等反应状态参数是随时间动态变化的,因此Z-DAC 的分区同DZFM一样,也是随流场更新而不断动态更新的.通过判断当前CFD 网格单元热化学反应状态参数,自动匹配对应Z-DAC 分区.注意,Z-DAC 分区同样可以是不规则形状,甚至包含了离散的“孤岛”形CFD 网格单元.在应用中,反应状态参数空间的网格划分在等当量比、临界点/熄火温度等反应活跃性剧变的区间需要局部加密以提高关键燃烧反应区域流动-化学耦合求解的准确性.

Z-DAC 中的实时机理简化方法采用Lu 和Law[60]提出的直接关系图法(DRG),该方法基于组分之间的相互作用关系图谱,评估其他组分对于某几类关键目标组分生成或者消耗速率的影响大小,从而基于给定的误差阈值去除影响较弱的组分和相关反应步.DRG 法尤其对于规模较大的机理具有高效率和高可靠性.本文采用了基于误差传递的直接关系图法(DRGEP)[61],其考虑不同反应路径的误差传递过程,即考虑包含若干基元反应路径的整体重要性而非单个基元反应的重要性.根据反应过程选择合适的目标组分,如燃料、氧气、氮气和部分关键自由基.根据基元反应关系图确定目标组分与所有其他相关组分的反应路径,建立所有目标组分的依赖路径图谱.在每条路径的相邻组分之间计算直接相关系数(DIC)

其中Nr为当前反应机理中所有反应的步数,ωi为第i步反应的净反应速率,νAi为反应i中组分A 的计量系数,PA和CA分别为组分A 的生成和消耗速率.进而累计乘积路径上所有相邻组分之间的DIC得出路径相关系数(PIC).定义总相关系数(OIC)为所有反应路径PIC 的最大值.由于OIC 区分了受检组分与目标物组分的相对重要性,OIC 低于指定OIC 阈值的组分被视为不重要组分,涉及该组分的反应路径被移除.较之原始DRG 方法,DRGEP 方法中OIC 考虑了包含多步基元反应的反应路径的相对重要性,反应路径上距离目标组分越远的组分重要性越低.调节OIC 阈值可以确定原始机理中需要去除的不重要组分,直至使给定工况条件下简化机理与原始机理计算的点火延迟误差达到最大允许偏差.DRGEP 简化中用于起始路径搜索的目标组分选择对最终反应机理保留路径影响较大.目标组分一般选取为燃料、氧化剂和主要产物(一般为H2O 和CO2).如果有需要与实验对比的关键中间自由基组分如OH 和CH2O 等,也需要加入目标组分以在机理简化中保留这些中间组分及其关联反应路径.一般而言,当非最终产物(H2O 和CO2)中C/H-O 元素当量比高于某一特定值(如0.1)时则认为处于较活跃的燃烧反应区,此时需加入CO 和HO2等链式关键中间组分;当反应区温度高于1000 K 时则需加入NOx等氮氧化物组分.

1.3 分区并行自适应建表(Z-ISAT)

碳氢燃料的详细反应动力学模型往往包含成百上千个基元反应,其在实际的三维燃烧计算中的应用需要巨大的计算资源.即便采用骨架或简化机理,现有的计算机技术水平仍难以为继.为了进一步提高化学反应源项的计算速度,局部自适应建表(ISAT)方法[62-63]通过对局部热力学初始和最终状态建立映射存储从而大大节约了计算时间.对于固定的反应机理,特定时间步长之后反应热力学最终状态是初始状态的唯一解,因此初始和最终状态之间存在一一对应的映射关系.该方法首先按照常规方法通过积分计算某一化学热力学初始状态经过特定反应时间后的最终热力学状态,然后存储起始状态和最终状态.在后续的计算中,通过在存储链表中对比当前热力学状态参数与存储参数,将相似状态对应最终状态的插值近似作为当地特定反应时间之后的最终状态,以替代冗长费时的积分过程.理论上热力学状态的映射存储可以一次完整建立.然而映射表的尺寸随反应组分数成幂次增长,如果完整构建热力学状态参数空间,即便是较小尺寸的反应系统存储都需要巨大的存储空间,且后续的查找效率也较低.事实上,大多数反应系统仅涉及完整热力学状态参数空间的一个较小子空间.该子空间的范围取决于实际燃烧条件,如燃烧反应机理、热物性参数和湍流特性等.由于事先无法确定,必须在模拟过程中动态建立.随着映射表的逐步建立,化学反应的计算速率也逐步提高.对于准稳态反应流,在存储空间允许的范围内大部分化学反应都可以通过查表和插值近似得出,从而极大地加速了复杂化学反应计算.

ISAT 映射表的建立过程如下:

(1) 计算开始时,ISAT 链表为空.直接使用刚性反应求解器基于反应机理对初始组分Y0、温度T0和压力P0在特定时间步Δt上进行求解积分,得到最终组分Y1和温度T1(假设恒压反应).ISAT 树状链表中单个叶节点存储信息为:初始热力学参数ϕ0=(Y0,T0,P0),最终热力学参数 ϕ1=(Y1,T1,P1),映射梯度矩阵A=∂ϕ1/∂ϕ0和在热力学参数空间内原点位于 ϕ0的精度控制超椭球面(EOA).

(2) 设有某初始热力学参数空间向量 ϕq,0,在ISAT 树形链表中查询与之最接近的节点设为 ϕ0,根据映射梯度矩阵线性插值计算近似最终结果ϕq,1=.如果 ϕq,1位于EOA 范围内,则上述线性插值的结果在制定误差 εtot允许范围之内,所计算结果将被取用.否则,则直接积分计算得到 ϕDI,1,确定映射误差,其中B为缩放矩阵.若 ε ＜εtot,则选用 ϕq,1,同时扩张EOA 范围以包括状态向量 ϕq,0;否则 ϕq,0生成新叶节点.

(3) 计算中按上述步骤依次查询或建立新的叶节点,最终形成一个二叉树形链表.当新的叶节点插入时,最相似原始叶节点的位置将变成一个二叉结点,重新连接原始叶节点和新的叶节点.在二叉结点上将计算并存储其所属两片叶节点的超级分割平面.计算中将以这个超级分割平面作为判据从树形链表的根部二叉结点依次往下层分支二叉结点上判断与 ϕq,0最相似的叶节点.

初始计算时,大部分树形链表的操作都是增加叶节点(增操作)或扩张原有节点EOA(扩操作).当大部分状态参数子空间都被树形链表的EOA 区域覆盖时,读取操作将更加频繁.尽管初期的增、扩操作会耗费额外的计算时间,但随着链表覆盖范围的完善,计算效率将逐步提高.注意,EOA 扩操作包含了一定程度的近似,即并不能保证EOA 范围内的所有点的误差均在指定误差范围之内,这是由于映射关系的非线性特性.实际计算表明EOA 范围内绝大部分状态参数子空间的误差都在指定的控制误差范围之内[64-65],并且局部极个别的误差越界对整体计算精度影响可以忽略,因此ISAT 的EOA 误差控制方法是符合计算需求的.

随着计算中增操作的进行,ISAT 链表的尺寸会越来越大,尤其是对于非稳态流场.过大的二叉树链表尺寸一方面耗费存储空间另一方面增加了二叉树搜索深度.由于先前历史中添加的叶节点被提取的概率逐渐降低,可以定期对链表进行清理维护以保证每个叶节点都有较高的提取率.本研究中实时记录每个节点被提取的次数并动态排名,根据排名确定最常用叶节点双向链表.当整个链表的尺寸超过设定上限后则清空整个二叉树形链表,再以最常用双向链表中的节点信息重新建立二叉树形链表.计算中发现,一般最常用链表的尺寸上限设为整个树形链表尺寸上限20%～ 50%可保持较好的加速比.

在并行计算中,一般采取针对各并行分区建立独立的ISAT 链表,相互之间没有数据交换,即纯粹本地建表(purely local processing,PLP) 策略[66].PLP 策略的缺陷在于在强瞬态、非均匀燃烧流场中无法保证反应计算载荷均衡性.特别是在超声速燃烧流场中,不同分区之间热化学反应状态分布差异极大,并且随时间演化而剧烈变化.此时不同分区ISAT 表的尺寸和成功取回率均差异较大.一般主要反应区所需维护的ISAT 链表尺寸较大,而无反应区的ISAT 链表尺寸甚至最小仅包含一个叶节点.此外,当并行分区方式或数目改变时,原有的ISAT 表需重新建立,无法重复利用.

本文提出了基于热化学反应状态进行动态分区的ISAT 建表策略,即Z-ISAT (zonal ISAT).该方法首先依据混合分数、温度、压力、反应进度变量等当地反应状态参数对流场进行动态分区,并在每个分区建立单独的ISAT 表.由于每个分区的热化学反应状态相对单一,因此ISAT 表一旦建成即可保持尺寸长期稳定,热化学状态点覆盖率和成功取回率效率也可大幅提升,降低了ISAT 链表树清理和重整频率.不同于从初始起便固定的并行分区,该基于热化学状态的动态分区建表策略尤其适用于强瞬态和空间非均匀燃烧场.计算中为保持各建表分区与热化学状态子范围空间的严格对应关系,需定期动态更新分区包含的CFD 计算节点.同样地,ISAT 分区为不规则区域并且可以包含孤立的CFD 节点.对应每种热力学状态分区的ISAT 表可以存储以便在计算重启或相似工况计算时再利用.

无通讯PLP-ISAT 策略适用于反应区较集中而稳定的燃烧情形.超燃冲压发动机中,反应区一般集中于射流尾迹或凹腔,但瞬态反应区存在剧烈的时空变化,甚至在特殊的配置条件下还存在振荡燃烧模态[67].主反应区在各热力学状态分区之间频繁移动,意味着相应的ISAT 表存在频繁的节点插入和增长操作,一方法增加了ISAT 表的尺寸即内存需求,另一方面降低了查表获取率,加剧了直接积分需求.为此本研究中发展了基于云计算的ISAT 并行策略.该云计算策略将本地热化学反应状态分为频繁取用和临时两类状态,只有频繁取用的节点才存储ISAT 链表,而临时状态点则通过打包分配到其余空闲节点进行直接积分计算求解.这里区分频繁取用的判据为ISAT 中被调用5 次以上的状态点.临时状态点的打包分发同样会占据内存带宽与通讯时间,为此研究中发展了主节点统计各从节点待分配状态点数目,绘制状态点重布向量地图,各从节点再依据此重布向量地图直接采用点对点通讯交互数据的方式,从而避免了集中分发所带来的额外通讯代价.该云计算策略既保证了ISAT 表较小的尺寸和极高的查询获取率,又通过并行分发策略提高了强瞬态流场中化学反应计算求解的载荷均衡性.

1.4 控制方程

其中 ν 为依赖于温度的运动黏性系数,可解尺度的应变率张量为

其中Sct=1 为湍流施密特数.

动态分区火焰面模型[70]中求解了四维空间(三维物理空间加一维混合分数空间)中的条件组分输运方程 (17) 而不是传统的Favre 平均组分输运方程.平均组分质量分数通过对混合分数空间内火焰面变量Qα进行概率密度函数加权平均得到

其中模型参数Cvar在亚声速燃烧中校准为0.1[73],但超声速燃烧中该值需进一步校准.

1.5 物理模型与求解器细节

大部分发动机燃烧室模拟都必须要考虑壁面的摩擦和冷却效应,相应地近壁区域的湍流模拟对发动机性能预测尤其重要.目前近壁边界层模拟主要有3 种方法:(1)直接数值模拟求解所有湍流尺度、(2)近壁湍流模型、(3)壁面函数方法.在存在惯性子区的湍流中,直接求解最小湍流尺度的单维度方向网格量需求为雷诺数的幂次方Re3/4,求解最小耗散时间步需求为Re1/2[74].然而壁面附近不存在明显的惯性子区,因此近壁边界层的模拟对网格有更加严苛的需求.在保持相同无量纲近壁距离(=Δynuτ/ν)的情况下,近壁边界层直接数值模拟求解所需要的网格数量为

高马赫数发动机中存在剧烈变化的温度和压力,理想气体模型仅在一定范围内适用.为此物性计算采用在热力学状态空间内分区计算,具体为低于300 K,低于1 kPa 或高于临界压力时密度计算采用真实气体模型,如广义对应状态法则[79];高温和中等压力条件下则采用计算较为简便的理想气体模型.热扩散率和组分扩散率计算分别采用JANAF 格式热物性数据库[80]和CHEMKIN 格式输运属性库[81].混合物平均黏性和热传导率的计算分别采用修正版Wilke 定律和摩尔加权平均[82].

计算平台采用基于OpenFOAM 开源框架[83]开发的可压缩密度基超声速燃烧求解器Amber (曾用版本名AstroFoam).该求解器采用低耗散修正[84-85]二阶Kurganov-Tadmor 中心迎风求解格式.求解器已经在超声速燃烧工况模拟中得到广泛验证[16-18,70,86-89].

图1 显示了计算中流动求解和湍流燃烧求解的耦合关系.其中流动模块Amber 求解式(12)～ 式(16) 基础N-S 方程以及基于Spalart-Allmaras 的IDDES 湍流模型[90].与PaSR[22]等直接求解三维物理空间内的平均组分输运方程不同,DZFM 湍流燃烧模型在四维空间内求解分区火焰面组分输运方程,平均组分则通过对条件组分进行概率密度函数加权平均得到.分区化学反应速率由条件温度控制,其通过统计方法估算得到,从而避免了复杂的可压缩非绝热条件温度建模与求解.为体现可压缩效应,平均温度由平均焓和平均组分反推得出.化学反应采用Jachimowski[91]针对超声速燃烧开发的13 组分33 步反应详细氢气机理.本研究中整个计算域被分割成18 000 个火焰面分区,其中在主流动方向分割为200 个薄片状区域,进而对每个薄片区域依据当地混合分数分割为90 个“年轮”状环带形子区域.火焰面分区可以是不规则、非连续区域,其随混合分数流场而动态更新以保证每个分区对应一个窄域的混合分数子空间.

图1 流动模型与湍流燃烧模型耦合示意图Fig.1 Coupling diagram of flow model and turbulent combustion model

1.6 测试工况

本文中选取了高超声速飞行中两种典型的发动机燃烧室构型进行对比测试.该发动机目标飞行马赫数10,飞行高度40 km.如图2(a) 所示,模型发动机包括进气道、隔离段、燃烧室和尾喷管4 部分.两种构型发动机的主要区别在于燃烧室中采用不同的被动增混方式,分别为中心支板方式(如图2(b) 所示)和壁面撑挡方式(如图2(c) 所示).发动机的整体长度为5.61 m,计算中包括了一部分长度为0.67 m的外流部分以模拟自由流条件下的进气道特性.进气道长度2.49 m,采用基于Busemann 基准流场的流线追踪法[92]设计,其中边界层黏性修正采用参考温度法.进气道的几何收缩比为10,静压压缩比约为62.隔离段平滑过渡进气道出口到圆形超声速燃烧室,燃烧室包括一段0.985 m 长的微扩张圆筒以过渡到单边扩张的尾喷管.中心支板燃烧室共计有24 个燃料喷孔,其中6 个燃料喷孔位于连接中心支板的支架侧面上,其余18 个燃料喷孔以3 个为一组位于中心支板侧壁.喷孔直径从0.8～ 1.5 mm 直接变化以适应于当地横向流动速度增强射流穿透深度.壁面撑挡超声速燃烧室同样采用了24 个喷孔,2 个为一组均匀环绕周向分别从撑挡部件顶部(3 组)和燃烧室壁面(9 组)直接喷注.

图2 燃烧室中心支板和壁面撑挡结构几何尺寸Fig.2 Geometric dimensions of the strut and pylon structures

本文中数值测试了上述两种典型构型高马赫数超燃冲压发动机在飞行马赫数10 和海拔40 km高度的性能.所模拟飞行工况参数如表1 所示,其中来流温度和压力均依据对应高度处大气层参数.燃料流均为声速喷注,总温为室温,流量保持总包当量比为1.0.

表1 测试飞行工况参数Table 1 Flight test conditions

计算域包括外流和发动机内部区域,其中中心支板发动机总网格数量为1.25 亿,壁面撑挡发动机总网格数量为1.4 亿.考虑到复杂的壁面突起结构和中心支板部件,大部分区域(96.6%体积)采用无结构四面体网格剖分,仅在拐角处采用少量楔形(3.39%体积分数)、棱柱形和金字塔形(共计0.01%体积分数)填充.平均网格尺度0.8 mm,在燃料喷孔附近局部加密至0.1 mm,近壁边界层采用15 层向内逐渐收缩的棱柱膨胀层,无量纲近壁距离y+＜1.网格质量分析表明99.4% 的网格正交性超过0.5,99%的网格偏斜度小于0.5.研究针对中心支板燃烧室开展了5413 万、7176 万、1.047 7 亿和1.251 亿共4 套网格的独立性验证,其余3 套较粗网格与最精细网格的平均壁面压力误差分别为1.7%,1.2%和0.8%,随网格尺度呈指数递减,因此可认为最精细网格的计算结果具有较好的网格独立性,如图3 所示.

图3 中心支板燃烧室壁面平均压力的网格独立性分析Fig.3 Grid independence analysis on average wall pressure

2 结果与讨论

2.1 REST 算例验证

基于本研究的数值模型和方法框架对国际上已公开发表的澳大利亚昆士兰大学(简称UQ)矩形变椭圆高超声速燃烧室(REST)开展了计算验证.昆士兰大学基于T4 反射激波风洞 (RST)开展了一系列的自由和半自由高马赫数超燃冲压发动机测试,在马赫数 7.5[6],8[7],8.7[8],10[9],直到12[10-12]均观察到了正净推力.与本研究类似,REST 超燃冲压发动机的自由射流模拟[15]均包含了部分飞行器前体以模拟内外流一体化耦合过程.验证REST 算例对应飞行马赫数12 和海拔36 km,燃料为当量比1.0 常温声速喷注氢气.模拟中采用了1.15 亿单元的无结构网格.如图4(a)所示,基于当前计算平台和数值框架所预测的平均压力与实验测量值吻合较好,同时捕捉到了与文献中预测类似的双峰结构.压力峰值位于燃烧室喷注点下游的微扩张段,同时进气道喷注点附近也有初始燃烧但产生的压升弱于弓形激波反射产生的第二峰值压升.如图4(b)～ 图4(e) 所示,计算中包含了0.5 m 长的飞行器前体部分以模拟边界层等外流效应对发动机内流的影响.上游喷注燃料的燃烧较为微弱,从温度和产物分布看剧烈燃烧反应主要发生于下游喷注点之后的微扩张燃烧段.马赫数分布表明燃烧室处于超燃冲压模式,近壁面处特别是燃烧室上侧壁面处的燃烧造成了连续而狭长的局部亚声速区域.数值纹影表明下游喷注点之前的进气道和隔离段存在多重反射激波并与剪切层相互作用,而在剧烈燃烧下游因体积膨胀效应无明显激波存在但可见丰富的湍流漩涡结构.

图4 基于相同计算框架的REST 马赫数12 高超声速燃烧室模拟Fig.4 Simulation of REST Mach 12 hypersonic combustor based on the same computational framework

2.2 中心支板与壁面撑挡燃烧室性能对比

图5 展示了两种构型高超声速燃烧室的内部流场结构.从马赫数分布来看,两种构型燃烧室均在支板或撑挡上下游出现了较大的亚声速区域,并在支板或撑挡对侧壁面诱发了明显的边界层分离.由于中心支板结构缩减,加大了流通面积,因此其对来流的减速效应更加明显,表现为支板处的超声速射流核心缩窄至近乎壅塞.壁面撑挡燃烧室由于凹腔的动量交换和稳焰助燃均进一步降低了局部马赫数,因此其下游的亚声速区域延伸更长距离.从温度分布来看,两种构型均出现了较长区域的喷注点前部燃烧,主要集中在燃料喷注撑挡一侧,表明部分燃料在回流裹挟下被带到上游并在较高的边界层滞止温度(＞1500 K)下被点燃.中心支板燃烧室主要燃烧区域位于较靠下游区域,表明燃料混合距离较长.而壁面撑挡燃烧室中凹腔的稳焰作用比较明显,上半部分在凹腔剪切层中出现了火焰,下半部分在凹腔内部回流区实现了稳焰.可见即使对于当地马赫数整体大于2.5 的高超声速燃烧室,凹腔的稳焰作用依然较为显著.壁面撑挡燃烧室的H2O 含量略高于中心支板回流区,同时整体2400 K 以上的温度也略高于整体2000 K 左右的中心支板回流区,表明壁面撑挡回流区反应强度更高.从两个构型流场对比来看,需要在燃烧室设计中一方面缩短两个喷注点的横向空间间距,另一方面尽可能让射流喷注点远离外扩形壁面,因其流动转向会增加下游尾迹的横向距离.CEMA 指数λe是衡量燃烧化学反应剧烈程度的重要指标[15,93-94],为化学反应速率雅可比矩阵的最大特征值,通常正值表示爆炸模式,负值表示趋于平衡的缓燃模式.正值发生于预混点火区域,而负值发生于扩散控制和后点火区域.图中所示为带符号对数值signλe× lg|λe|,本研究中燃烧室中均为负值,表明燃烧主要由扩散控制,高速流动中极短的滞留时间通常难以维持高温预混气团至点火状态.其中对应温度高于2000 K 的上游回流区、支板后主反应区和凹腔等燃烧剧烈处的λe均为108s-1量级的特征反应速率.尾喷管内因温度普遍低于1500 K,即便仍有未完全反应的燃烧,反应特征速率降为105～ 106s-1.由于反应机理中包含了部分O2,N2和H2O 的离解反应,因此空气来流区域也存在特征时间104s-1量级的反应.

图5 中心支板与壁面撑挡高超声速发动机内部流场Fig.5 The internal flow fields of strut and pylon high-Ma scramjets

图6 展示了沿程流向各截面提取参数对比.壁面撑挡燃烧室静压峰值85 kPa,为中心支板燃烧室压力峰值60 kPa 的1.5 倍.继进气道压缩引起的压升之后,回流区燃烧进一步提高了压力;其中壁面撑挡燃烧室对应回流区反应更加剧烈因此压升效应更加明显并产生了一个较小的峰值.壁面撑挡燃烧室的压力峰值位于凹腔处,由剧烈燃烧诱导;并在其后缓慢下降至尾喷管入口.而中心支板燃烧室的压力峰值由支板或撑挡诱导的斜激波及其在壁面的反射激波引起,支板下游未见明显的燃烧诱导峰值;其后燃烧室压力出现了一个40 kPa 的平台区,燃烧近乎等压燃烧.与冷态燃烧室不同,燃烧条件下混合效率的计算应采用元素法[15],即统计氢元素与氧元素的比例,以考虑已反应部分的燃料以及组分间差异扩散效应.回流区两种构型燃烧室的混合效率接近.说明支板或撑挡之后,壁面撑挡燃烧室的混合效率迅速升高远超中心支板燃烧室,这是由于壁面撑挡燃烧室的射流穿透深度较高,并且从图6(c) 凹腔前后燃料浓度分布可知凹腔顶部剪切层和内部回流区也起到了一定程度的增混作用.混合效率在燃烧室下游靠近尾喷管入口处(x=4.2 m)趋于稳定值,表明在尾喷管内平均马赫数5.0 以上的高速流动发生的混合极少.壁面撑挡燃烧室的最终混合效率为85.3%,大幅高于中心支板燃烧室的61.2%.燃烧效率的曲线与混合效率高度相似,最终燃烧效率分别为83.9%和60.7%,仅略低于混合效率,表明高超声速燃烧室因滞止温度较高且氢气的特征反应时间极短(微秒级)存在混合即燃烧的特点.这也与CEMA 分析中指出的燃烧多为扩散混合控制的观察结果相符.从这个角度上说,高超声速燃烧室提高燃烧效率的关键在于增强混合.总压计算中考虑了真实气体比热随温度变化的特性,而理想情况下的非等熵关系式.根据等熵关系式积分可得

图6 沿流向的准一维性能指标Fig.6 Quasi-one-dimensional performance indices along the flow path

其中p和T为当前静压、静温,Tt为总温,R为气体常数.高超声速燃烧室的总压损失极高,中心支板燃烧室的总压损失为99.3%,壁面撑挡燃烧室的总压损失为99%,意味着产生了接近2 倍的熵增.两种构型燃烧室的总压损失曲线基本一致,其中壁面撑挡附近因局部动量交换更为剧烈导致总压损失更加迅速.在设计中,进气道型面微调为平滑过渡到燃烧室入口,因此中心支板燃烧室的进气道略短导致总压损失起始延后,但在保持相同几何压缩的条件下进气道末端的总压损失一致.定义轴向冲量为

其中A为流向面积,ρ为密度,U为流向速度.根据动量守恒,轴向推力F=ΔI,图6(e)中下降曲线(负梯度)表示阻力,上升曲线段表示正推力.一个很明显的结果是进气道因压缩来流是产生阻力的主要部件;对应x=3.48 m 处,中心支板和壁面撑挡均产生了较大的阻力;壁面撑挡燃烧室x=3.64 m 处的轴向冲量突跃位于凹腔,同时产生的推力和阻力,总体上表现为阻力;近等直燃烧段产生了微弱阻力;尾喷管是获得推力的主要部件.从以入口冲量无量纲化的冲量曲线对比可见,两种构型燃烧室进气道阶段产生的阻力和尾喷管产生的推力均比较接近,主要差别在于燃烧段的阻力差别较大:特别是中心支板产生的阻力显著高于壁面撑挡,凹腔虽然产生了一定程度的阻力但也有效增加了混合和燃烧,其利弊需要进一步探索.

表2 对比了两种构型燃烧室的总体性能.其中壁面撑挡燃烧室的未装机净推力为344 N,比冲比较理想为1234 s;而中心支板燃烧室净推力较弱仅为99 N,比冲437 s 与常规火箭发动机性能接近,体现不出吸气式发动机的优势.壁面撑挡燃烧室的燃烧效率也较为理想,高于通常认为的可产生推力的80% 准则[95].壁面撑挡燃烧室的压比16 也显著高于中心支板燃烧室的11,由此相比中心支板产生了增幅36% 的无黏(压力)推力.同时,中心支板燃烧室的阻力比壁面撑挡燃烧室高23%.中心支板燃烧室进气道外形在保持几何压缩比的条件下做了一定的调整以适应微扩的隔离段,因此其进气道捕获流量略低,但计算中保持总当量比恒定.

表2 壁面撑挡和中心支板燃烧室总体性能对比Table 2 Overall performance comparison between pylon and strut combustors

2.3 动态分区火焰面特性

如图7 所示,动态分区火焰面模型中火焰面结构函数随时间和空间而不断演化以准确表征当地燃烧化学反应状态.分区条件平均温度取决于当前分区的焓增/减程度和反应进度,用于控制当前分区反应速率.随着流向距离增加,分区条件平均温度随着反应趋向平衡递增,并在燃烧室末端x=4.478 m 处达到峰值,后续在尾喷管内随着部分总焓转换为动能又逐渐降低.上游的火焰面状态会被下游分区继承并继续反应,因此随着反应趋于平衡产物H2O 的浓度持续增加.自由基O 由O2在高于2500 K 温度的离解反应生成,并为在H2-O2的链式反应中起到重要作用的中间组分 d [O]/dt=k1[ H][O2]-k2[ H2][O],(k1和k2为反应速率常数).一般认为燃烧反应过程中O 的浓度在宽温度范围内相对恒定[96],离解产生的过量O 在燃烧段被逐渐消耗,在温度低于离解温度(2500 K)仅有燃烧反应的尾喷管中保持相对恒定.火焰面的演化除了在当前分区控制温度作用下演化以外,还不断与相邻分区通过对流交互信息.以O2浓度为例,其在燃烧段被不断消耗的同时也不断从相邻未燃来流分区卷吸高浓度的O2,因此在燃烧段其浓度甚至会短暂上升.与H2O 的不断生成对应,O2的浓度在自上游至下游的流动过程中被逐渐消耗,而在尾喷管中因卷吸作用减弱,O2在反应中被逐渐消耗.

图7 不同流向距离处火焰面函数变化Fig.7 Representative flamelet variations changes at different flow distance

动态分区火焰面的动态过程一方面是指火焰面的时空演化,另一方面火焰面分区也不断自适应于流场而变化,如图8(a)所示.首先根据流向进行空间分区以区别不同流动阶段的火焰面.精细的流向分区是准确捕捉点火距离的关键,因其准确区分了从纯混合、点火到充分燃烧反应甚至熄火等不同阶段的反应状态.分区独立性分析[70]表明,当分区空间间隔小于一定程度时,流场计算结果不再依赖于分区数目.一般在超声速燃烧室模拟中,一个合适的分区准则为分区间隔应不大于其点火距离,并在燃烧段保证至少50 个空间分区.即空间上对流场进行分区后,进一步在各“树墩”状剖分段依据当地混合分数进行分区.图8(b) 展示了支板前、中、后典型位置处的混合分数分区情形.其中在中段距离燃料喷注点附近,当地混合分数覆盖从0 到1 的分布,因此分区较为密集,而前后段由于扩散效应仅覆盖低混合分数段,因此分区数目减少.再往上游或下游的更远端,当地混合物状态为完全未预混的纯来流或已经混合较均匀的状态,此时混合分数分区数目缩减为单个,即单一火焰面即可充分表征当地反应状态.在依据混合分数分区时,为提高主反应区的表征保真度,一般需要在当量比混合分数附近加密分区,如图8(b) 所示.从图中还可以看到,混合分数分区形状不一定是规则的,甚至可以是分散的孤岛集合,在计算空间对流时依据格林公式沿交界面进行积分.当流场不均匀性增大时,可进一步增加分区指标.例如根据当地温度、马赫数等进一步进行分区剖分以提高分区条件均一性.

图8 火焰面分区示意图Fig.8 Diagram of flamelet division

2.4 分区自适应化学特性统计

图9 展示了分区自适应化学(Z-DAC)中各分区反应机理的统计.从图9(a) 中可见前13 步反应为在Z-DAC 简化中几乎被始终保留的反应步,而14～ 19步反应被使用的概率略低,约为90% 左右.20～ 33步反应则在63.5% 的分区机理简化中被忽略.前13 步反应为H2的氧化反应,在燃烧反应区被完整保留,仅在纯空气来流区域被忽略.第14～ 19 步涉及过氧化氢(H2O2)的反应仅在OH 浓度较高的活跃反应区才被激活.而后14 步涉及N 元素的反应则在1800 K 以上的高温区域才会被激活.从各分区的反应机理尺寸统计来看,64%的区域采用了31～ 33步反应,27%的区域采用了19 步反应,8%的区域采用了11～ 13 步反应.另有极少的区域采用了7 步反应,为涉及N2和O2离解反应高温且无燃料区域.可见分区自适应化学方法在将近一半的计算域上降低了反应求解计算代价,特别是在无燃料区反应机理的简化幅度更加明显.

图9 分区自适应化学(Z-DAC)统计数据Fig.9 Zonal dynamic adaptive chemistry (Z-DAC) statistics

2.5 与传统燃烧模型效率对比

表3 对比了DZFM 模型与以PaSR 为代表的传统有限速率模型的计算效率.PaSR 模型除了需要对流场中所有CFD 网格单元逐一积分求解当地化学反应外,还需要求解全组分输运方程,因此计算代价一般远高于流动-化学反应解耦模型.超燃冲压发动机中的特征流通时间一般以毫秒为量级,这里选取1 ms 流动物理时间考察不同燃烧模型及加速方法的计算效率.在相同并行规模504 核(CPU 型号Intel Xeon CPU E5-2640,单核主频2.40 GHz)的计算测试中,纯PaSR 模型的单毫秒的计算周期需要24 d,采用Z-DAC 加速后计算周期缩短为约20 d,采用ZISAT 加速缩短为14 d.由于Z-DAC 机理简化与ISAT 建表查询方法近乎相互独立,因此二者耦合使用的加速效率近似于二者单独加速比的乘积:332.92 h × (468.16 h/582.61 h)=267.52 h=11 d,仅略小于其耦合使用的12 d.ISAT 与DAC 的耦合使用的组合在部分文献[97]中被称为建表自适应化学(TDAC),本研究中的方法为基于分区的Z-TDAC.纯DZFM的单毫秒计算周期缩短为52 h,对比纯PaSR 的加速比高达11 倍.采用DZFM 模型,每个火焰面分区上仅有一维混合分数空间的若干状态点需要更新化学反应状态.Z-DAC 尽管对当前分区反应机理进行了最大程度的简化,然而因反应状态点数量大幅降低,因此其加速效率不如与传统有限速率模型耦合时明显,仅实现了5%的加速.同样地,Z-ISAT 的加速率也为5%.二者耦合使用的Z-TDAC 加速率略高,也仅为8%.表3 列出了不同模型计算1 ms 物理时间对应的CPU 时间.可见分区解耦湍流燃烧模拟将三维物理空间千万量级的化学反应状态点求解近似为较少的火焰面分区数目(本文中为1.8 万个分区)乘以一维混合分数空间内的离散状态点上的反应求解,从而极大地提升了加速比,同时使得传统的反应计算加速方法加速效率相对不明显.由于机理简化和建表查询均会产生一定的近似误差,因此在加速比较低的情况下可以考虑优先提高精度.

表3 DZFM 模型与PaSR 模型计算1 ms 物理时间对应的CPU 时间Table 3 CPU time of DZFM and PASR model calculation of 1 millisecond physical time

图10 基于Borghi 图统计分析了以氢气为燃料的高超声速燃烧室中各分区的湍流化学反应交互作用关系(TCI).其中戴姆克拉数Da定义为湍流特征时间 τt与化学反应特征时间 τc之比,Da=τt/τc.基于大涡模拟数据计算湍流特征时间时需要包括可解部分的湍流脉动kres[98],即 τt=(kres+kt)/ε .化学反应特征时间 τc为CEMA 指数的导数.卡尔洛维奇数Ka定义为化学反应特征时间与最小湍流尺度之比Ka=τc/τk,其中Kolmogorov 特征时间 τk=(ν/ε)1/2.据此有Re～(τt/τk)2=Da2·Ka2.根据Ka和Da可将湍流化学反应交互作用关系分为3 种模式:(1) 火焰面模式Ka＜1 且Da＞10,(2) 薄反应区模式:1 ＜Ka＜100且Da＞10,(3) 慢反应模式:Ka＞100 且Da＜10.从两种燃烧室的统计分析中可以看出,超过90% 的绝大部分分区处于快速反应假设成立的火焰面模式,5%左右的少部分分区处于薄反应区模式,小于3%的分区处于慢速反应区.氢气超声速燃烧与碳氢燃料超声速燃烧最大的区别在于,碳氢超声速燃烧绝大部分处于薄反应区(大于50%)和慢速反应区(小于30%)[99].可见即便对于来流马赫数极高的高超声速燃烧室,氢气反应速率仍然远快于滞留时间,并且燃烧室中极高的滞止温度也进一步提高了氢气的反应速率,在混合即燃烧的情况下决定燃烧效率的瓶颈在于高效增混.高超声速燃烧室中雷诺数高达108,最小湍流特征尺度为微米量级,这意味直接求解高超声速燃烧的直接数值模拟计算代价仍然是目前计算机技术无法承受的,同时也预示高超声速流场中跨越极广尺度的湍流与速率较慢的碳氢燃烧化学反应存在极其复杂的交互作用.

图10 分区湍流化学反应交互作用关系Borghi图Fig.10 Borghi diagram of zonal turbulent chemical reaction interaction relationship

图10 分区湍流化学反应交互作用关系Borghi 图(续)Fig.10 Borghi diagram of zonal turbulent chemical reaction interaction relationship (continued)

3 结论

本文展示了基于动态分区概念的高马赫数全尺寸超燃冲压发动机的内外流耦合一体化改进延迟分离涡(IDDES)模拟研究.研究建立了完整的动态分区燃烧模拟框架,包括动态分区火焰面湍流燃烧模型(DZFM),以及分区自适应化学(Z-DAC)和分区并行自适应建表(Z-ISAT)等化学反应加速模型.前者通过分区解耦的思想既准确表征了当地湍流化学交互作用关系,又有效提升了整场湍流燃烧的计算效率.后两者通过在分区框架内对化学反应机理进行动态实时简化和建表查询,可进一步提升当前分区内化学反应的求解效率.沿不同空间位置处的火焰面结构函数的变化表明动态分区火焰面模型可以捕捉反应混合物微团内温度和组分等参数随反应进程逐渐演化的物理过程.精细的流向分区准确区分了从纯混合、点火到充分燃烧反应甚至熄火等不同阶段的反应状态,因而是准确捕捉点火距离的关键.混合分数分区形状不必是规则的,甚至可以是分散的孤岛集合,其在不同流向距离处的分布密集程度也自适应于流场而变化.

首先通过对比有完备实验数据的马赫12 REST标准高超声速燃烧室模型,验证了分区模拟框架的保真性.分别基于1.25 和1.4 亿网格动态分区框架对马赫数10 和海拔40 km 条件下全尺寸的中心支板(strut)和壁面撑挡型(pylon)两类构型氢气高超声速燃烧室进行了深入的数值分析,揭示了高超声速燃烧中一系列独特的规律和机理.

壁面撑挡燃烧室的未装机净推力为344 N,比冲比较理想为1234 s;而中心支板燃烧室净推力较弱仅为99 N,比冲437 s.壁面撑挡燃烧室的燃烧效率也较为理想,高于通常认为的可产生推力的80% 准则.壁面撑挡燃烧室的压比为16,也显著高于中心支板燃烧室的压比11,由此相比中心支板燃烧室产生了增幅为36%的无黏(压力)推力.同时,中心支板燃烧室的阻力也比壁面撑挡燃烧室高了23%.两种构型燃烧室的总压损失均高达99%,其中壁面撑挡附近因局部动量交换更为剧烈导致总压损失更加迅速.

基于Borghi 图的数值分析表明,上述两类燃烧室内燃烧的模式均为扩散主要控制的火焰面模式.超过90% 的绝大部分分区处于快速反应假设成立的火焰面模式,远大于碳氢超声速燃烧室中的约10% 的比例.在混合即燃烧的情况下,提高燃烧效率的关键在于增强混合.壁面撑挡燃烧室中撑挡结构的燃料喷注由于处于亚声速区因而其穿透深度较中心支板钝体表面喷注更深,从而取得了更好的混合效率,进而提高了燃烧效率.

分区自适应化学Z-DAC 在准确地区分了各分区的化学反应状态,其中H2的氧化反应在燃烧反应区被完整保留,涉及过氧化氢(H2O2) 的反应仅在OH 浓度较高的活跃反应区才被激活,而涉及N 元素的反应则在1800 K 以上的高温区域才会被激活.64% 的区域采用了近乎完整的33 步反应机理,27%的区域采用了19 步反应,8%的区域采用了11～13 步反应,另有极少的高温纯空气离解区域采用了7 步反应机理.可见分区自适应化学方法在将近一半的计算域上降低了反应求解计算代价,特别是在无燃料区反应机理的简化幅度更加明显.

相比于传统的有限速率PaSR 模型,DZFM 模型将三维物理空间千万量级的化学反应状态点求解近似为较少的火焰面分区数目(本研究中为1.8 万个分区)乘以一维混合分数空间内的离散状态点上的反应求解,实现了高达11 倍的加速比,同时使得传统的反应计算加速方法加速效率相对不明显.由于机理简化和建表查询均会产生一定的近似误差,因此在加速比较低的情况下可以考虑优先提高精度.