1. 超临界燃烧仿真:为什么说这是燃烧计算里的硬骨头
做燃烧仿真的人,大概都听过一句话:层流燃烧好算,湍流燃烧难算,超临界燃烧最难算。这几年我陆陆续续接触过几个超临界燃烧相关的项目,从液氧煤油发动机的喷注燃烧到燃气轮机预混燃烧室的高压工况,说实话,每一次都被超临界状态下的物理行为折磨得不轻。但反过来讲,也正是这些折磨,让我把很多教科书上写得模糊、甚至干脆没写的东西彻底弄明白了。
先说清楚一个概念:超临界燃烧,简单说就是在压力超过流体临界压力、温度接近或超过临界温度的条件下发生的燃烧过程。这时候氧化剂和燃料往往以“稠密气体”的状态喷入燃烧室,不再有清晰的液滴表面,也没有传统意义上的蒸发过程。很多人第一次接触这个领域时会觉得奇怪:液体没有了蒸发现象,那燃料怎么和氧化剂混合?答案其实藏在“伪沸腾”这个现象里。
伪沸腾是超临界流体最反直觉的行为:当压力超过临界压力后,流体在跨越伪临界温度时,定压比热会出现一个很尖锐的峰值,密度急剧下降,膨胀率猛增,看起来非常像液体的沸腾,但实际上是单相流体的剧烈热膨胀。这种密度的大幅变化会直接改变喷注射流的破碎和混合机制,进而影响火焰的稳定性和燃烧效率。
超临界燃烧仿真之所以难,难在四个层面:第一,物性描述不再能用理想气体定律,必须依赖真实流体状态方程;第二,流场中密度变化剧烈,湍流和压缩性效应耦合非常复杂;第三,化学反应机理在高压下表现出和常压完全不同的反应路径;第四,喷注器附近的微小尺度现象(射流破碎、局部混合、火焰厚度极小)对整体燃烧影响巨大,而常规工程网格根本捕捉不到。
这篇文章适合谁?如果你正在做高压燃烧室设计、火箭推力室仿真、超临界喷注混合研究,或者准备从常规燃烧仿真转向超临界工况,这些内容就是给你写的。我自己用的工具包括开源CFD求解器和商业软件,后处理常用Python做数据处理,文里涉及的参数和步骤都是实际项目里验证过的,不是纸上谈兵。
2. 物性模型:超临界燃烧仿真的第一道生死关
2.1 理想气体为什么在超临界区彻底失效
刚开始接触超临界燃烧时,很多人习惯性地沿用常压燃烧仿真的思路,把气相燃烧室里的工质当成理想气体来处理。压力一上来,理想气体假设的误差会迅速放大。
理想气体模型的核心假设是分子间没有相互作用力、分子本身不占体积。在常压高温燃烧环境下,这个假设误差很小。但在超临界条件下,流体密度接近甚至超过液相密度(比如液氧在超临界态下密度可达800kg/m³以上),分子间距被压缩到很小,分子间引力和分子本身体积再也不能忽略。按照理想气体状态方程计算出来的密度、焓和比热,可能和真实值相差30%到一倍以上。
有一个非常典型的例子:在超临界态下,流体跨越伪临界温度时比热会猛增数十倍,形成我们前面说的“伪沸腾”。理想气体模型给出的比热是一条平缓的曲线,完全捕捉不到这个尖峰。燃烧室里的温度场、火焰面位置对热流和比热极其敏感,比热算不准,温度场就全歪了,整个仿真就失去了意义。
所以超临界燃烧仿真的第一步,必须换用真实流体状态方程。这也是我在第一轮超临界仿真中踩过最大的一次坑:用理想气体跑完整个流场,温度场看着挺漂亮,但跟实验对照时热流密度偏差超过40%,返工几乎重来了大半轮计算。
2.2 状态方程怎么选:PR、RS还是其他
真实流体状态方程是超临界仿真的基石。目前工程上用得最多的是两三类:Peng-Robinson(简称PR)、Redlich-Kwong及其改进型(比如RK-Soave,也就是RS)、以及针对复杂碳氢燃料的多参数状态方程。
我个人的习惯是:对于甲烷、乙烷、氧气、二氧化碳这类简单分子体系,PR方程足够用,计算效率和精度平衡得比较好。对于煤油、航空煤油这类多组分复杂燃料,纯PR方程对近临界区的密度预报偏差会偏大,这时候通常用RK-Soave或者能自定义二元交互参数的状态方程做修正。
这里有一个很关键、但经常被忽略的点:状态方程里的纯组分参数(临界温度、临界压力、偏心因子)取错,整个物性计算就完蛋了。这些参数看起来只是三个数,但不同数据库给的值并不一致。我之前收到过一份来自某合作方的网格文件,燃料选的是煤油替代物,临界温度填的是文献里某个早期数据,结果燃烧室里流体的压缩性完全不对,火焰几乎贴在喷注器出口。后来把临界参数改成该替代物官方推荐值,流场才正常。
混合规则也值得多说两句。多组分流体在超临界状态下的混合行为不是简单的摩尔加权平均,而是需要用混合规则把纯组分参数组合起来。最常用的是van der Waals混合规则,需要指定二元交互参数kij。kij看起来只是一个0到0.1之间的小数,但它对混合物临界点的预测影响极大。kij取0和取0.08,计算出来的混合物临界密度可能相差20%以上。如果文献里查不到目标混合物的kij,至少要做一次敏感性分析,看看它对关键结果(比如火焰面位置)的影响有多大。
2.3 热物性和输运系数的耦合处理
有了状态方程,密度和相平衡问题解决了,但燃烧仿真还需要比热、焓、粘度、导热系数。这几个量在超临界区的行为也不循常理。
比热和焓的处理方式是:用理想气体比热加上偏离函数修正。理想气体比热一般用NASA多项式拟合,偏离函数则由状态方程积分得出。这样处理的好处是:在高温燃烧区(几千K),非理想效应减弱,计算结果自动回归到理想气体值;而在低温稠密区,修正项起主导作用,精确捕捉伪沸腾的比热尖峰。
粘度、导热系数的超临界修正比较麻烦。工程上常用的做法是基于对应态原理,在低压输运系数基础上叠加上密度的修正项。但务必注意,这类修正只适用于远离临界滞止点的区域。在临界点附近,输运系数会出现奇异性,几乎所有工程修正模型都会失真。好在燃烧室工况虽然压力超过临界压力,但温度大部分时候远离临界温度,输运系数的误差对整体流场影响相对可控。
操作层面,如果你用的是商业CFD软件,内置的NIST物性数据库通常有比较完整的纯工质物性,但计算成本很高。一个折中方案是:在准备阶段把目标压力下的物性表算好、存成表格,仿真时用查表插值代替实时物性计算。这套做法能节省大概一到两倍的计算时间,而且物性精度几乎没有损失。我自己在长时长燃烧仿真中基本都走这条路线。
提示:物性模型的选择和验证,不应该等整个仿真算完再做。开工前单独跑一个零维或一维算例,验证伪临界温度附近的比热尖峰是否出现、密度变化是否平滑、临界点参数是否和NIST数据对得上。这一步只需要几分钟,能省下后面几周的返工时间。
3. 网格与计算域:超临界火球薄到让你怀疑人生
3.1 尺度估算:为什么常规网格会被秒杀
燃烧仿真中,网格要捕捉的关键特征尺度是火焰面厚度和射流混合层的宽度。常压甲烷空气火焰的层流火焰面厚度大约在0.5到1毫米量级,工程网格做到零点几毫米,配合较好的湍流模型,基本能看个大概。
超临界工况完全不同。高压下密度升高,化学反应速率加快,火焰面大幅变薄。在典型的超临界燃烧室压力(比如15到25MPa)下,层流火焰面厚度可以缩到几十微米甚至更低。也就是说,超临界火焰面比常压稀薄一到两个数量级。
而超临界喷注射流还有一个特点:射流以稠密流体状态喷入高温高压燃烧室,由于密度差巨大,射流核心区会保持很长的“液态针”结构,破碎和蒸发过程被抑制,混合主要由剪切层中的湍流扩散决定。剪切层厚度在喷注出口附近可能只有喷孔直径的百分之几,喷孔直径通常0.5mm到1mm,剪切层初始厚度可能只有几十微米。
于是问题来了:喷注器近场你需要几十微米的网格来分辨混合层,火焰面需要更细的网格来分辨反应区,而整个燃烧室可能有半米长、几十厘米直径。所有尺度塞进一套网格,直接算DNS(直接数值模拟)的计算量是个天文数字。据我粗算,一个全尺寸工程燃烧室若做全分辨DNS,网格量奔着上千亿去了,以现有算力根本不现实。
所以在工程仿真中,我们必须接受网格分辨不足的现实,然后用合适的模型去弥补。
3.2 工程网格策略:让粗网格发挥最大价值
我常用的策略分三步。
第一步,喷注器近场局部加密。喷孔出口、剪切层起始区、火焰稳定区是全场最敏感的区域,网格尺寸控制在射流直径的1/20到1/50,并且用渐进过渡网格逐步放大,避免突变网格导致的伪反射。以1mm喷孔为例,近场最小网格做到20到30微米,在很多工程算例中已经能獲得比較合理的射流穿透和混合趋势。
第二步,中远场适度放粗,但仍需保证火焰面附近网格密度足够支撑温度梯度的捕获。一个实用的判断标准是:网格能否平滑地分辨出温度场从低温射流区到高温燃气区的过渡,如果温度出现非物理的台阶状跳跃,说明网格还不够细。
第三步,利用网格自适应加密。新型求解器基本都支持基于温度梯度或密度梯度的自适应加密,能把精细网格“长”到火焰面和剪切层上,其他区域保持粗网格。这是解决超临界燃烧计算量矛盾的最有效工程手段。不过自适应加密也不是万能的,在三维非稳态计算中频繁加密和粗化会带来额外开销,对于算力紧张的项目,我会优先做近场静态加密,把自适应留到后期精细化阶段。
3.3 近壁区与壁面热流的陷阱
超临界燃烧室的壁面热流极高,推力室壁面热流动辄几MW/m²起步,燃气侧壁面温度的控制是设计成败的关键之一。但近壁区的网格处理,恰恰是很多人出错的地方。
壁面的网格除了要满足y+要求外,还必须注意超临界流体的密度分层效应。近壁区温度梯度极大,燃气在壁面附近温度骤降,密度剧烈上升。这种密度壁面层对换热系数有决定性的影响。如果壁面网格太粗或者壁面函数不适用于强变密度流动,壁面热流会低估百分之二三十以上。
我的做法是:壁面第一层网格厚度按目标y+小于1来设置,并保证近壁区有至少10层以上网格覆盖边界层内的密度变化区域。必须说明,经典的壁面函数建立在常密度、对数律速度分布基础上,在强变密度超临界流动中可能失效。如果软件支持,应选择可压缩、变密度的壁面处理,或者干脆用低雷诺数壁面模型。
注意:在超临界燃烧仿真中,我不是很建议一上来就开标准壁面函数图省事。壁面热流算不准,直接导致冷却结构设计和寿命评估全盘漂移。宁可多花一倍网格量,把边界层分辨出来,也不要在这里省钱。
4. 化学机理与湍流-化学反应耦合:从骨架机理到高保真模型
4.1 高压化学反应机理的选取思路
常压甲烷燃烧的GRI-Mech 3.0,共53个组分、325个基元反应,是很多人默认的选择。但这个机理是在1个大气压附近标定的,直接外推到20MPa以上,不少反应速率常数会出错。超临界高压环境下,第三体效应、压力依赖反应、碰撞效率都会变化,高保真机理的表现并不一定好。
工程上更推荐使用专为高压工况标定或验证过的简化机理。哪怕是“骨架机理”级别,比如包含十几到二十几个组分的简化甲烷氧化机理,只要压力适用范围覆盖目标工况,在工程精度要求下往往比照搬常压全机理更可靠、计算更快。
这里我遇到过一个真实的坑:有个项目用GRI机理跑超临界甲烷燃烧,火焰结构和温度场看着都正常,但点火延迟时间明显偏短。换成高压标定的简化机理后,点火位置向后移动,与实验的高速摄影吻合度大幅提升。原因就在于GRI机理中几个关键链分支反应在高压下的速率常数外推误差累积起来,改变了火焰传播速度。
选机理时还有一点值得注意:燃烧室里的实际燃料往往不是纯物质,而是多组分混合物(比如煤油替代物、裂解气混合物)。这种情况下建议先做机理简化或替代物标定,在动力学层面把混合物的氧化过程等效映射到可计算的反应网络上。需要提醒的是,替代物的组分比例变了,机理可能就不再适用,所以替代物选取的依据必须充分,不能拍脑袋凑。
4.2 湍流-化学反应耦合:湍流火焰模型的选择逻辑
湍流和化学反应在超临界状态下的耦合是燃烧仿真最核心也最深层的问题之一。这里没法回避一个事实:多数稳态工程简化模型(比如涡耗散模型EDM),以“化学反应速率无限快、由湍流混合率控制”为前提,在超临界条件下往往不够准确。
超临界火焰的一个特点是局部熄火和再着火现象频繁。射流以稠密冷流体喷入高温环境,混合层内局部当量比和温度的梯度极其陡峭,流动拉伸率也很高,火焰很容易在局部被吹熄。EDM这类混合控制模型天生不会预测局部熄火,它要么给出全燃,要么给出不燃,中间态是失真的。
我自己做超临界燃烧工程仿真的主力方案有两类。
第一类是火焰面模型(FGM或FPV),核心思路是:把湍流火焰看成一组一维“层流火焰片”的统计集合,火焰内部结构由混合分数和一个进展变量(或热释放率变量)参数化。计算时不需要直接求解每个反应组分,而是把组分浓度和温度做成火焰面数据库,在流场中通过混合分数和进展变量的输运方程来插值。这个方法计算成本低、对稳态喷射火焰效果好,很适合喷注燃烧室内的扩散火焰。
第二类是输运概率密度函数方法,用蒙特长枪或随机粒子方法求解组分和温度的联合概率密度输运方程,对湍流-化学相互作用的表现力强,能捕捉局部熄火和再着火,但计算成本比火焰面模型高十倍甚至更多。这通常用于研究型小算例、单射流喇叭燃烧室构型的精细分析,不适合大型三维全程仿真。
梯度陡峭条件下,还有一个模型失效的共同原因:梯度假设。无论是涡扩散还是火焰面模型,都依赖于梯度输运的正常预测。但在超临界射流近场,密度梯度极大,标量耗散率强烈各向异性,梯度输运会出偏差。我一般会在喷注近场的后处理中单独查看局部标量耗散率分布,确认模型是否在合理范围内工作。
4.3 稳态与非稳态仿真的取舍
超临界燃烧中火焰天然具有非稳态特征。射流内部的剪切涡、大尺度相干结构、火焰面抖振,都对燃烧和传热有影响。但如果一上来就做非稳态大涡模拟(LES),你很快就会面临“时间步长微小、模拟物理时间极短、算一天只走了几毫秒”的窘境。
合理的推进路线是分步走:
- RANS稳态起步:先拿RANS模型加火焰面或EDM快速把流场结构、喷注穿透深度、整体燃烧效率、壁面热流摸一遍。这个阶段重点是物性和边界条件的合理性。一般十天到两周能完成一轮。
- 在RANS结果基础上做LES:把RANS解作为初场,换用大涡模拟,近场加密网格,非稳态推进。只关注火焰稳定区、局部熄火概率、压力脉动这些RANS给不了的东西。
- 数据后处理中重点同步实验:高速摄影的火焰结构、动态压力传感器的频谱、热电偶温度曲线,这些都是验证非稳态仿真的金标准。
需要留意的是,LES从RANS解出发有一个初始化适应期。RANS流场的湍动能和耗散率分布是“平均”意义上的,LES需要一段时间让大尺度涡结构成长起来。这个时期一般要若干倍穿越时间,具体多长没有定论,但如果你只推进了1ms就拿来出结论,那基本是拿初场和噪声在骗自己。
5. 超临界燃烧仿真实战调试:压力震荡、温度尖峰和收敛陷阱
5.1 压力场震荡:一个不解决就全盘皆输的麻烦
超临界燃烧仿真在工程应用中出现的第一类麻烦,通常是压力场在局部区域出现非物理震荡或锯齿状分布。具体表现是:明明边界条件和物性设置都合理,残差却怎么都降不下去,或者压力场出现严重的棋盘式振荡。
这种震荡的根本原因,在于密度-压力-温度三者的耦合关系在超临界区高度非线性。密度对温度极其敏感,而温度又受燃烧放热强烈调制,求解器中压力修正方程在这种强耦合下很容易失稳。我自己排查这类问题时的顺序如下:
第一步,检查物性表的光滑性。如果你用查表法,物性表本身在伪临界温度附近的数据点如果取得太疏,插值会引入锯齿状的不光滑性,哪怕物理上合理,数值上也会激发出震荡。把伪临界温度附近的物性表加密一倍以上,经常立竿见影。
第二步,检查压力-速度耦合算法。超临界工况下传统的SIMPLE类算法在密度修正环节可能比较吃力,遇到强烈密度变化时优先尝试PISO类算法,并缩小时间步长。在LES中时间步长控制在库朗数0.3到0.5以内,基本是硬要求。
第三步,检查边界条件的初始化质量。超临界燃烧室的初始场不能拿常温常压场做线性外推,而是要先做“纯混合冷流”或“无反应等温高压”预备计算,等压力场和密度场在燃烧室构型下先达到一个合理的稳定态,再开启反应源项。直接把冷流初始加上反应项点火,几乎必然产生压力冲击波,有时候数值直接爆掉。
5.2 火焰面温度的“尖峰”问题和伪沸腾陷阱
火焰面温度尖峰是另一个高频问题。在正确的物理图景下,超临界扩散火焰的峰值温度应该在某个合理的绝热燃烧温度附近。但仿真结果里常出现火焰面温度高出理论绝热温度好几百度、甚至上千度的局部尖峰。这不是真实的燃烧现象,几乎都是数值伪影。
出现尖峰的常见原因是混合分数插值表的分辨率不足。火焰面数据库在富燃侧往往是温度非线性最剧烈的区域(因为有些组分在这个区域发生大量吸热裂解反应),数据库在这一侧点取得稀疏,插值就容易过冲,在火焰面局部造出假的温度尖峰。
另一个常见原因是热物性耦合断裂:温度剧烈升高时,物性算出的比热和焓没有平滑过渡,造成局部焓值失衡,温度被虚假“顶高”。这两个原因的排查方法,都是在后处理中取火焰面上一个典型位置,把组分浓度、温度、焓值沿法向截面画出来,对照物性表和火焰面数据库逐点检查。如果焓值曲线在某个节点出现跳跃,问题多半出在物性表或数据库的插值连续性上。
5.3 判断收敛的正确姿势:别只盯残差
超临界燃烧仿真的收敛判断,比常规燃烧仿真要多留几个心眼。残差降下来不代表算对了,残差降不下来也不一定代表算法坏了。
我的经验是,至少同时检查三组指标:
第一组是全局守恒量,包括燃烧室出口质量流量、总焓进出口平衡、燃烧室平均压力。如果这些量在数百步迭代内没有明显漂移,说明稳态解有了基本骨架。
第二组是局部关键物理量,比如喷注近场的射流穿透长度、某监测点的温度时间序列、火焰面位置。这些量最能暴露模型和网格的缺陷,比全局量敏感得多。
第三组是数值噪声水平。观察压力场和密度场云图有没有非物理的高频波纹。哪怕残差已经很低,如果云图上有肉眼可见的棋盘格纹,那就是数值噪声,这个结果是不可信的。
在做LES时,判断“统计收敛”还需要做时间平均。取一段时间窗口,统计平均温度场和脉动速度场,如果窗口长度加大一倍而平均结果变化很小,统计才算平稳。很多论文里给出的漂亮的LES平均场,仔细一查可能只平均了很短一段模拟时间,其统计置信度是可疑的。
5.4 从实验到仿真再回到实验的闭环验证
最后聊聊验证。超临界燃烧实验本身难度极高,高压、高温、透明窗口又少,能拿到的实验数据往往只有壁面热流、燃烧室压力、推力、高速摄影的火焰形态。但恰恰是这些“少得可怜”的数据,用好了就能把仿真质量“钉死”。
我做仿真验证时有个固定流程:先做冷流(无反应)对比,验证喷注射流的穿透深度和混合场结构。再做热态验证,对照燃烧室压力、出口温度分布和壁面热流。最后做火焰形态的定性对比:仿真得出的火焰长度、抬举高度是否和高速摄影轮廓匹配。三层对照下来,模型哪里有问题基本就暴露了。
这里有个我反复提醒项目组同事的细节:不要拿调整过边界条件去“凑”实验数据。很多边界条件(如入口温度、壁温、燃烧室背压)本身有测量不确定度,调整它们使仿真和实验对上,是合理的参数标定过程。但如果为了“好结果”反复压网格、换机理、调模型常数,最后得到的是过拟合数据,不是可信预测。超临界燃烧本身计算量大、调试周期长,坚持这个底线不容易,但这是仿真工作真正产生价值的前提。
我个人的切身体会:超临界燃烧仿真项目几乎没有一次成型的时候,每次都会在物性表、机理选择或者网格方案上栽跟头。但恰恰是这些返工,逼着人把真实流体的物理机制想透了。如果你刚开始做这个方向,别怕慢,先从一小段射流、一个简化机理、一平方米的精细网格做起,把物性和火焰结构跑通,再逐步放大计算域。这样积累的每一步经验,都能在后续全尺寸仿真中替你省下成倍的时间和代价。