☰
超临界燃烧仿真核心难点拆解:物性突变与真实气体效应
2026/10/9 9:21:35 网站建设 项目流程

干燃烧仿真这行有几年了,刚做完一个超临界燃烧算例的复盘,正好借这个“超临界燃烧”的主题,把这条路从头捋一遍。所谓超临界燃烧,通俗说就是燃烧室压力高于燃料和氧化剂临界压力、温度也跨越临界温度之后的燃烧过程,在火箭发动机、超临界二氧化碳循环燃烧室、高压燃气轮机上几乎躲不掉。说句实在话,别把超临界燃烧理解成“普通高压燃烧”,高压只是入场券,真正难缠的是两件事:一是物性突变,二是真实气体效应。这篇文章就围绕这两个核心难点展开,把我踩过的网格、物性库、燃烧模型相关的坑也都交代一遍,给准备上手超临界燃烧仿真的朋友一条可以直接走的路径。文章偏实操,原理点到为止,重点放在“怎么设”“怎么调”“怎么排查”。

1. 先把“超临界燃烧”这几个字拆明白

1.1 临界状态到底是个什么状态

在纯物质的相图上,气液平衡线有一个终点,就是临界点。水的临界压力22.064 MPa、临界温度647.096 K;二氧化碳是7.377 MPa、304.13 K;甲烷是4.599 MPa、190.56 K;氧是5.043 MPa、154.59 K。压力一旦超过临界压力,加压不会让气态“液化”,体系进入超临界区。但超临界区内部并不是均匀的:沿着等压线升温,密度会连续下降,却在某个温度附近密度梯度极大、定压比热出现尖峰,这个温度叫伪临界温度。

伪临界温度不是固定值,它随压力升高而升高,每一条等压线都有自己对应的尖峰位置。拿25 MPa下的二氧化碳来说,伪临界温度大约在358 K附近,常温液态二氧化碳射入这个压力环境后,会被环境加热跨越伪临界线,密度从几百千克每立方米迅速掉到几十千克每立方米,这个转变过程在实验里叫伪沸腾,不是真的沸腾,没有气液界面,却有强烈的密度脉动。

仿真时最直观的感受是:伪临界温度附近,CFL数稍微大一点就震荡,残差曲线上上下下像心电图。这不是求解器抽风,而是物性在这个区间对温度极其敏感,密度、比热、黏度、导热系数几乎同时剧烈变化,任何数值误差都可能被放大成非物理振荡。

1.2 哪些工程场景绕不开超临界燃烧仿真

超临界燃烧不是实验室里猎奇的东西,工程场景非常多,而且都在关键装备上。

火箭发动机推力室是典型中的典型。液氧煤油发动机室压普遍在18到25 MPa,液氧喷注温度约90 K,虽然温度远低于氧的临界温度,但压力远高于氧的临界压力,所以液氧处于“冷态超临界”,不具备普通液体那样的表面张力,但密度仍然接近液相。液氧甲烷发动机更典型,甲烷临界温度仅190.6 K,喷注器附近甲烷温度比临界温度低很多,进入燃烧室后快速升温跨过伪临界线,经历从类液态到类气态的连续转变,喷射形态和混合过程完全不同于亚临界下的液柱破碎。

超临界二氧化碳循环燃烧室是另一个热点。整套系统压力常年在20到30 MPa,主流工质是超临界CO2,燃料通常为甲烷或合成气,氧化剂是高纯氧,燃烧产物直接和CO2工质掺混,温度从七八百K升到两千K以上。这个场景的难点在于主流工质CO2占了绝对多数,物性库和化学反应机理都要重新适配。

高压燃气轮机和超临界水氧化处理也在逐步进入这个范畴。前者燃烧室部分工况已经接近跨临界,后者在水的超临界点以上处理有机废水,涉及无焰燃烧和超临界水化学反应。虽然应用领域五花八门,仿真内核高度一致:真实气体物性、伪临界区捕捉、高压化学反应动力学。

2. 仿真的核心难点到底难在哪

2.1 物性突变:伪临界温度附近“翻脸比翻书还快”

超临界燃烧仿真最坑的就是物性。以25 MPa的二氧化碳为例,在350到400 K这个区间内,密度从700 kg/m3量级一路降到100 kg/m3以下,定压比热在伪临界点附近出现一个非常窄的高峰,峰值可以是远离临界区的四五倍,导热系数和黏度也会在局部出现拐点。这种物性突变对数值格式是很大的考验。

我习惯把伪临界区的物性变化比作在悬崖边开车:温度是油门,密度是海拔,过了伪临界点,轻轻一脚油门,海拔就掉几百米。如果网格在崖边不够密,差个零点几个开尔文,密度计算值就完全不一样。

实际操作中要做的第一件事,就是在预处理阶段把工作压力附近的物性曲线全部拉出来看一遍。我常用NIST数据做基准,把密度、定压比热、黏度、导热系数随温度的变化曲线画在同一个坐标系里,标出伪临界温度,然后用真实气体状态方程去拟合。仿真里如果用的不是查表法,而是状态方程解析计算,需要确认拟合精度在关键区间是否够用,特别是cp的峰值能不能复现。峰值复现不了,温度场和火焰结构就会失真。

网格也必须在伪临界区加密。加密原则不是均匀加密,而是沿伪临界等值线附近做局部加密,让跨越伪临界带的网格层数至少达到10到20层。我早期做CO2射流算例时吃过亏,伪临界区只给了五六层网格,密度梯度根本拉不出来,射流核心区温度偏高两三百K,后面查问题查到网格上,后悔不已。

2.2 真实气体效应:理想气体假设彻底失效

常压燃烧仿真用理想气体状态方程,问题不大,误差两三个百分点,燃烧场趋势完全能抓住。但在25 MPa下,纯理想气体密度误差可以达到百分之二三十,比热、焓值的偏差也一样夸张,继续用理想气体模型就是在错误的地基上盖楼。

真实气体状态方程里,工程上最常用的是Peng-Robinson方程和Soave-Redlich-Kwong方程,两者区别主要在对偏心因子的处理上。PR方程形式如下:

p = RT/(v - b) - a / [v(v + b) + b(v - b)]

其中a和b由临界温度、临界压力和偏心因子确定,a还附加一个随温度变化的函数。混合物场景下,a和b通过混合规则从纯组分参数组合得到。组合时需要注意不同组分之间的交互作用系数,这个系数通常需要借助实验数据或经验关联式拟合,二氧化碳与烷烃、氮气与氧气等体系的交互参数在文献里都有推荐值。

计算物性时还有个容易忽略的点:真实气体效应不只是改变密度,还通过剩余焓、剩余熵进入能量方程,影响温度场和火焰传播。很多朋友只在密度上加了真实气体修正,焓值还是用理想气体表,结果能量不守恒,温度场出现莫名其妙的偏移。正确做法是热力学状态函数全部基于同一个状态方程计算,至少保证在主要组分范围内自洽。

2.3 湍流-燃烧-物性的强耦合

超临界燃烧仿真不是把各模块拼在一起就行,真正的挑战是湍流、燃烧、物性三者互相反馈。密度脉动直接影响浮力和剪切层发展,剪切层又决定燃料与氧化剂的混合,混合控制着化学反应放热,放热反过来再改变温度和物性。这个闭环在超临界压力下反应速度比常压快得多,因为高压下化学反应时间尺度更短,物性变化又随时在扰动流场。

湍流模型方面,RANS里SST k-omega用起来最稳,适合做初场和稳态预热;但超临界射流里的大尺度拟序结构,RANS根本解析不动,最终还是要上LES。LES的亚格子模型里,WALE和Sigma模型相对不容易在高剪切区域过度耗散,我用下来比Smagorinsky好。需要提醒的是,亚格子模型公式里通常隐含着“物性恒定或缓变”的假设,超临界场景里最好把密度和黏度的亚格子脉动也粗略考虑进去,否则射流穿透深度经常算得偏短。

燃烧模型的选型更敏感。稳态火焰面模型在常压湍流扩散火焰里非常成熟,但超临界压力下,火焰面内部温度极高、局部密度极低、标量耗散率范围又宽,准稳态假设可能失效。如果用稳态火焰面,建议至少做一步验证:把涡耗散时间与化学时间对比,看Damköhler数是否远大于1,不够大的话就别硬撑着。

更稳妥的方案是输运PDF模型或EDC模型。输运PDF对CPU的消耗确实大,但它把化学反应和湍流混合解耦,物性变化不影响化学反应子步的求解,我们组在超临界CO2循环燃烧室算例里就是用LES加输运PDF跑下来的,虽然一个算例烧了一个月,但温度场和实验对标误差控制在5%以内。EDC计算量小很多,但高压下E系数需要重新标定,不标定会高估反应速率。

2.4 数值稳定性:算得动比算得准更重要

压力高、密度变化剧烈、化学反应快,三个条件叠加,数值稳定性就变得非常棘手。超临界区里声速变化很大,密度基求解器在低压区表现好,到伪临界区常常面临刚性极大的问题;压力基求解器在低压可压缩流里精度又不够。这里需要全速度算法或者预处理方法,把两种求解器的优势结合起来。

我用过的组合方案有这么几个,供参考:

  • OpenFOAM系:reactingFoam做常压没问题,超临界工况要改物性接口,工作量大;rhoReactingCentralFoam对激波和接触间断的分辨率更好,但计算开销高。国外有课题组直接基于rhoCentralFoam扩展了超临界物性和化学反应源项,我照着类似思路改造过一版,烧得慢但稳定。
  • Fluent:用density-based solver加PR真实气体模型,能跑通,但如果化学反应机理带了几十个组分,计算速度感人;pressure-based solver加real-gas模型在超临界射流中容易振荡,需要开低松弛。
  • 专业燃烧软件如CFD++和气动热化学软件也有超临界模块,火箭领域用得多,通用性差一些。

网格质量方面,第一层网格高度按y+小于1去控制,超临界下密度变化大,壁面附近的速度梯度比常压更陡,近壁网格宁密勿疏。剪切层和混合区要做局部加密,格间膨胀比控制在1.2以内。伪临界区可以先用各向同性加密,跑通后再用自适应加密进一步优化。

时间步长控制是稳定性的核心。显式方案下,全局Courant数在准稳态阶段建议不超过0.5,初始阶段0.1到0.3更安全。超临界射流发展初期,伪临界区里局部速度可能异常高,需要把局部Courant数监控起来,设置钳制上限。隐式方案虽然理论上稳定,但伪临界区物性剧烈变化时也会出现非物理振荡,建议打开双时间步推进,内部子迭代步数提高到10到20步。

3. 实操流程与关键设置,照这个走不迷路

3.1 先定工况与状态空间

仿真开始前,先把工况定死。我在做超临界CO2循环燃烧室算例时,用的是这样的设置:燃烧室压力25 MPa,主流工质CO2,氧化剂入口是O2/CO2混合气,温度800 K,速度60 m/s;燃料入口是CH4/CO2混合气,温度600 K,速度20 m/s;出口环境压力25 MPa,壁面采用等温条件。

这个工况的好处在于,入口温度已经高于CO2临界温度,但压力远高于临界压力,主流流体处于超临界区,而燃料和氧化剂混合后在火焰附近会跨越伪临界温度,燃烧场内有丰富的物性变化特征,非常能考验仿真设置。做二维轴对称算例时,计算域取喷射器下游一定长度,一般取喷口直径的40到60倍,出流边界才能基本不受数值反射影响。

定好工况后,把所有入口组分的物性曲线在目标压力下算一遍,存成参考数据。这一步别省,后面排查温度尖峰、密度负值都靠它。

3.2 网格:伪临界区和剪切层就是命根子

网格画得好不好,决定这个算例能不能收敛。我的经验是先粗网格跑一个稳态场,观察伪临界等值线大致位置和火焰面位置,再围绕这两个区域做加密。不要一开始就在全计算域铺细网格,成本太高,伪临界区位置都不知道,加密了也白加。

喷口附近剪切层的网格尺寸,按喷射直径的百分之一到千分之二控制。举个例子,喷口直径10 mm,剪切层最小网格取0.1到0.2 mm,这个量级在RANS下够用,LES下最好再细一点。近壁网格按y+小于1去估算第一层高度,超临界下密度高、黏度低,边界层更薄,我吃过亏的第一版网格第一层高度就是按常压工况放的,结果y+跑到5以上,壁面温度算得偏高。

伪临界区网格要顺着等值线走。25 MPa的CO2伪临界温度在358 K附近,火焰外围和入口低温区之间必然有一条伪临界带,顺着这条带加密,厚度方向上至少放15层网格。我在后处理里常用密度梯度幅值找这条带的位置,确认它是不是和火焰面相交、是不是被数值扩散抹平了。

3.3 物性模型:状态方程和输运性质要配套

物性模型是超临界仿真的灵魂,这一步没做对,后面全白搭。国内做燃烧仿真的同行常用Fluent,设置时在流体材料里切换成real-gas,选择Peng-Robinson状态方程,输入各组分的临界温度、临界压力、偏心因子和二元交互作用参数。OpenFOAM里通常是自定义一个thermo类,把PR方程写进去,再通过NASA多项式拟合理想气体部分的热容。

输运性质也不能拉下。黏度和导热系数在伪临界区有异常行为,最简单的做法是让它们随温度压力按对应态原理计算,虽然精度不如实验关联式,但数值连续性比查表法好。扩散系数方面,超临界压力下组分扩散系数不是简单的低压二元扩散公式,需要引入高压修正项,否则组分输运会失真。

我建议每一步都做物性对标:用状态方程计算出的密度、cp和NIST数据对比,偏差控制在3%以内。偏差一大,优先检查混合规则的交互参数是不是给错了,其次检查临界参数的来源是不是有误。曾经有一个算例,计算密度比实验低8%,查来查去发现CO2和CH4的交互作用系数取的是默认值0,换成文献值0.091之后误差立刻降到1%以内。

3.4 湍流与燃烧模型:稳妥的组合方案

模型选择没有银弹,但有一个组合区间相对稳妥:RANS预热加LES正式计算,燃烧模型用稳态火焰面或EDC,物性统一走PR方程。

先跑RANS-SST,拿一个初步的冷态混合场和热态温度场,这个阶段不需要太高的网格精度,收敛标准也放低,出口流量平衡误差小于1%就行。接着用RANS结果作为LES的初场,切换WALE亚格子模型,时间推进用双时间步,物理时间步按喷射直径和喷射速度估算。

燃烧模型方面,如果用稳态火焰面,火焰面库生成时一定要用真实气体状态方程去算组分和温度的关系,不能用理想气体火焰面库。另一个关键是把标量耗散率的上限设够。超临界扩散火焰的耗散率通常比亚临界高,上限设太小,火焰面库会把高温熄火区截断,算出来火焰过分旺盛。我习惯先把工况下的标量耗散率范围做一维对冲火焰扫一遍,找到熄火极限,再乘以一个1.5到2倍的安全系数作为火焰面库的上边界。

用EDC的话,需要把反应区内的化学反应理解为发生在精细结构里,高压下精细结构和湍流能量耗散的关系与常压不同,E系数从默认值往下调。曾经有论文推荐在超临界压力下把E系数打五折,我个人试下来不敢说普适,但确实比默认值稳定。

3.5 边界条件与工况设定:细节决定成败

下面给出我那个25 MPa超临界CO2循环燃烧室算例的边界条件表格,参考时可以按自己的工况线性缩放:

边界类型温度速度组分
燃料入口质量入口600 K20 m/sCH4 20%, CO2 80%
氧化剂入口质量入口800 K60 m/sO2 30%, CO2 70%
燃烧室出口压力出口——压力25 MPa
燃烧室壁面等温壁1200 K—无滑移

质量入口的好处是流量固定,稳态更容易收敛。压力出口必须给定目标压力,同时设置回流温度和回流组分,防止迭代初期回流把非物理组分带进去。壁面温度我取了1200 K,比火焰温度低,能模拟燃烧室壁面被冷却的效果,后处理时也能观察近壁热分层。

燃料入口组分里掺了大量CO2,这是模拟废气再循环的效果,真实系统里为了控制燃烧温度也会这么做。别小看稀释比例,它直接影响火焰温度和火焰面位置,抄别人的边界条件前一定要先确认稀释比是否一致。

3.6 求解与监控:盯住该盯的物理量

求解器参数设置上,我习惯先全隐式低库朗数跑2000步,等火焰锚定住之后再慢慢放开。RANS阶段残差降到1e-4以下还不够,必须看出口温度是否平稳、燃烧室最高温度是否还在缓慢爬升。燃烧问题最怕看起来收敛了,实际上火焰还在缓慢移动,这时候看残差没用,得看积分量。

LES阶段监控的就更多了。除了残差,我还会放几个监测点:射流轴线上的一列温度探头、火焰面附近的组分浓度探头,以及整个计算域的瞬时最高温度。超临界工况下瞬时最高温度可能有周期性脉动,这正常,但要确保脉动不冲破物理上限。如果最高温度持续高于绝热火焰温度两三倍,说明物性或者燃烧模型出了问题,要停下来查,而不是把松弛因子调小硬压。

出口设置为压力出口后,要定期检查质量守恒:入口总流量和出口总流量的偏差控制在0.5%以内。超临界流体会在某个局部区域短暂出现回流,只要回流幅度不大、不影响入口条件,是可以接受的。

3.7 后处理:用伪临界等值线解读火焰结构

后处理很多时候比计算还花时间,因为超临界燃烧场的特征不像常压燃烧那样直观。我必看的三张图:密度场、温度场、伪临界等值线。密度场能直接反映射流的类液态核心长度和扩散混合程度;温度场看火焰形状;伪临界等值线则告诉我流体在哪个位置发生了类液态到类气态的转变。

把伪临界等值线和温度场叠加,能判断火焰面是位于伪临界区上游还是下游。很多超临界燃烧工况里,燃料射流先跨过伪临界线再着火,如果仿真里火焰位置出现在伪临界线之前,说明混合或者反应速率设置有问题。从密度场提取核心长度也是个好办法:沿射流轴线找密度从接近液相值掉到半液相值的位置,这个长度和实验阴影图里的核心长度直接对应。

组分后处理上,我会看OH质量分数和CO2质量分数分布,CO2这个组分的分布能说明主火焰区和后燃区各自烧掉了多少燃料,尤其在CO2循环燃烧室里,CO2既是产物又是工质,辨析起来更有意义。

4. 常见问题与排查技巧实录

4.1 残差降不下去,火焰飘忽不定

残差不降是最常见的问题,原因通常不在求解器而在物性。伪临界区物性变化太陡,网格又不够密,局部密度估算跳变,导致压力修正方程不稳定。排查思路很直接:把残差最高的区域云图画出来,看是不是集中在伪临界带附近,如果是,加密那个区域,同时把物性表插值改成高精度分段插值。

另一个典型原因是燃烧模型和物性的时间尺度不匹配。稳态火焰面表在物性剧烈变化时,查表出来的温度、密度可能与当前流场状态不一致,表现为火焰位置左右摇摆。解决方法是把火焰面库的网格在标量耗散率方向加密,或者改用非稳态火焰面。

4.2 温度场出现局部尖峰

温度尖峰,尤其是超过绝热火焰温度两三倍的尖峰,多半是组分混合和反应速率匹配出了问题。常压燃烧里化学反应速率足够快,湍流混合是控制因素;超临界燃烧里高压让反应速率更快,混合控制依然占主导,但物性变化可能让局部组分的扩散系数变得特别小,组分堆积导致局部反应放热过度。

我排查温度尖峰时会先看那个区域的组分分布:如果局部氧和燃料同时大量共存,说明混合已经完全发生,反应速率也没有被湍流限制住,问题在燃烧模型。如果共存区内标量耗散率特别高,火焰面库查出来的温度和密度可能已经到了熄火边界,插值会给出不合理的高温,这时候需要扩展火焰面库的耗散率范围。

另一个冷门但真实存在的原因是,差分格式在伪临界区产生了数值振荡,压力出现锯齿状波动,温度跟着被拉高。遇到这种情况,降低局部Courant数,或者把对流项格式从二阶迎风切换成带限制器的混合格式,一般能压住。

4.3 密度出现负值

密度负值是燃烧仿真里最让人头疼的错误之一。超临界压力下,PR方程在极端组分组合或极端温度区间可能求根失败,返回非物理的密度值。这种情况经常出现在火焰面附近的贫氧侧和富燃侧,局部组分与空气差别很大,二元交互参数又覆盖不到,状态方程的根就飞了。

处理办法有几个:一是对组分浓度做钳制,确保任何单元格的组分质量分数不低于零、不高于一;二是对温度和压力做钳制,防止进入非物理区间;三是采用混合状态方程的保守求根策略,在Newton迭代发散时自动切换到二分法。OpenFOAM里可以在物性类里加保护逻辑,Fluent那边通常靠限制器处理。我把这个保护逻辑称为“最后的防呆机制”,没有它,任何算例早晚会在某个角落翻车。

4.4 火焰面库生成失败或熄火区被错误截断

稳态火焰面库生成时,最常遇到的问题是熄火极限标量耗散率算不准。超临界压力下,化学动力学和输运性质同时改变,一维对冲火焰的可燃上限和常压完全不同。我自己遇到过火焰面库在某个耗散率之后出现温度跳跃,查下来是对冲火焰求解器没有用真实气体物性,火焰面温度场不连续。

解决方法是重新生成火焰面库时把物性模型和主求解器统一,输入组分和压力也严格一致。生成后要扫几条标量耗散率下的温度曲线,确认火焰结构和熄火位置合理,再装入主仿真。

4.5 常用问题速查表

把上面几个问题整理成速查表,方便现场排查:

现象可能原因处理建议
残差周期性振荡伪临界区网格不足、物性插值不连续加密伪临界带,改用分段三次插值,降低松弛
温度尖峰燃烧模型耗散率上限不足、组分扩散系数错误扩展火焰面库范围,检查高压扩散修正
密度负值状态方程求根失败、组分越界添加组分钳制,TRUNC限值,状态方程求根加上限保护
火焰面库温度跳变对冲火焰未用真实气体物性统一物性模型,重新生成库
射流穿透深度偏短LES亚格子耗散过大换WALE或Sigma模型,加密剪切层
稳态“假收敛”火焰面缓慢漂移,残差却达标监控最高温度和出口温度,延长计算时长

5. 给新手的路径建议和一点个人体会

提到这里,多说几句个人实操体会。每次开新算例,我都要求自己和组里同事先做一维和零维验证,再进三维。超临界燃烧仿真变量太多,直接上三维满配,出了问题连排查入口都找不到。零维验证点火延迟时间,一维验证火焰面和熄火极限,二维轴对称验证射流穿透和火焰锚定,全部通过了才敢展开三维。

另外,物性模型一定不要迷信默认值。不管是Fluent还是OpenFOAM,默认的物性参数都是按常见工况标定的,超临界压力下一定要逐个检查临界参数、偏心因子、交互作用系数。我见过太多算例卡在物性库上,最后发现只是某个组分的偏心因子填错了。

最后分享一个小技巧。超临界燃烧的后处理里,随手画一张“当地密度 vs 当地温度”的散点图,所有网格点都叠上去。正常情况下,这个散点图会形成一条清晰的“伪沸腾走廊”:低温高密度区、伪临界密度骤降区、高温低密度区三段分明。如果散点一团乱麻,密度和温度关系完全不成体系,说明数值振荡已经污染了流场,这时候别纠结后处理,回头查物性和数值格式才是正解。这条散点图我用了很多年,屡试不爽。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询