做燃烧仿真这行也有几年了从最开始照着教程跑算例到后来自己改求解器、调离散格式回头再看整个知识体系发现最难啃、也最绕不开的一块硬骨头就是有限体积法。特别是做燃烧仿真不管是非预混火焰还是预混燃烧不管是RANS还是LES底层全是有限体积法在撑着。很多新手一上来就扑向湍流模型、化学反应机理结果算出来发散或者结果明显不对回头看往往都是栽在了离散格式、网格质量、时间推进这些有限体积法的基础环节上。这篇东西我把自己对有限体积法在燃烧仿真中的应用理解做个系统性整理从控制方程怎么变成代数方程到界面插值怎么选再到燃烧仿真特有的刚性问题怎么处理最后把实际算例中踩过的坑也一并列出来希望能给正在学燃烧仿真或者准备自己写求解器的朋友一点参考。1. 燃烧仿真到底在解什么先把方程组底牌摊开1.1 从连续介质到控制方程燃烧仿真的数学底盘燃烧现象在宏观尺度上本质上是流动、传热、传质和化学反应的耦合过程。我们做数值仿真并不是在模拟每一个分子而是在连续介质假设下把燃烧过程用一组偏微分方程描述出来然后用数值方法去求解。这组方程的核心包括连续性方程、动量方程N-S方程、能量方程以及组分输运方程。拿组分输运方程来说它的基本形式是∂(ρYk)/∂t ∇·(ρvYk) ∇·(ρDk∇Yk) ωk这个方程里左边第一项是组分质量分数的瞬态变化第二项是对流项右边第一项是扩散项第二项是化学反应源项。对整个燃烧仿真来说这串方程组合到一起描述的就是“流动带着组分跑扩散让组分混合化学反应把一种组分变成另一种组分”这个过程。我遇到不少朋友问为什么燃烧仿真一定要用有限体积法用有限差分或者有限元行不行这个问题其实背后有个很实际的原因。燃烧问题里经常出现强间断、大梯度尤其是火焰前锋附近温度和组分浓度变化非常剧烈。有限体积法是基于积分形式的方程它在控制体上做通量守恒这一点对保证“质量守恒、能量守恒”有先天的优势。燃烧数值结果如果守恒性不好算出来的火焰温度、组分分布都会失真。1.2 为什么偏偏是有限体积法三种离散方法横评CFD里有三大类离散思路有限差分法FDM、有限元法FEM、有限体积法FVM。很多人听说过名字但不清楚它们之间的实质差异这里我拿它们做燃烧仿真的效果和适用性来对比一下。有限差分法是最早发展起来的直接在网格点上用差分公式代替偏导数。它的数学推导简单直观计算效率高适合结构网格上的简单几何。但缺点也很明显要求网格足够规则对复杂几何边界难以处理而且在守恒性上需要额外注意。燃烧仿真里的火焰面、旋流燃烧室、多孔介质这类复杂几何用FDM会很吃力。有限元法在固体力学领域是霸主它的思路是把求解域分成有限个单元在每个单元上用试函数逼近真实解。它对复杂几何的适应性很好但当处理含对流占优的流动问题时经典的伽辽金格式会产生数值振荡需要额外的稳定化处理应用在燃烧这类强非线性问题上实现复杂度偏高。有限体积法的思想最贴近物理本质直接在控制体上做积分利用高斯散度定理把体积分转化为面积分把控制方程变成“通过控制体边界的净通量 内部源项 控制体内物理量变化率”这种形式。这种做法的最大优势就是局部守恒性天然满足不管网格多稀疏、多扭曲通量守恒在离散层面严格成立。对燃烧仿真来说在燃烧室内质量、动量、能量守恒是结果的底线有限体积法从根上保证了这一点。拿一个我实际用过的例子来说之前算一个带钝体稳焰的燃烧室网格质量不算好近壁面处扭曲比较大如果是有限差分格式早就开始振荡发散但有限体积法配合合适的限制器结果依然稳定温度场和实验数据对得上。正因为这样主流燃烧仿真工具OpenFOAM、Fluent、CFX包括很多自研求解器底层清一色是有限体积法。理解FVM就是理解这些工具求解机制的钥匙。2. 有限体积法的核心操作积分、散度定理与界面插值2.1 第一步在控制体上做积分有限体积法的起点是在一个控制体上对控制方程做积分。把计算域划分成一个个互不重叠的控制体可以理解为一个个小单元每个控制体中心有一个要计算的物理量值比如速度、压力、温度、组分质量分数。控制体可以是结构化网格的六面体也可以是非结构化的四面体、多面体。以不带源项的标量输运方程为例∂(ρφ)/∂t ∇·(ρvφ) ∇·(Γ∇φ)在一个控制体V上做体积分得到∫_V ∂(ρφ)/∂t dV ∫_V ∇·(ρvφ) dV ∫_V ∇·(Γ∇φ) dV这里逻辑很简单我们想知道控制体内总物理量怎么变化就需要把控制体内每一点的瞬态项、对流项、扩散项加起来。问题是控制体内各点的场是未知的连续函数这个积分没法解析算出来所以才需要数值离散。这里有个关键细节计算时我们并不是把体积分直接“求出来”而是通过高斯散度定理把它变换成面积分再用数值积分近似。这也是为什么我一直强调理解散度定理是理解FVM的关卡。2.2 第二步用高斯散度定理把体积分变成面积分高斯散度定理是有限体积法的理论支点。它说的是一个矢量场在控制体上的体积分等于这个矢量场在控制体表面上的面积分。∫_V ∇·F dV ∮_S F·n dS用这个定理上面的输运方程就变成∫_V ∂(ρφ)/∂t dV ∮_S (ρvφ)·n dS ∮_S (Γ∇φ)·n dS翻译成物理语言控制体内物理量的变化率等于通过控制体表面流入/流出的净对流通量加上净扩散通量。这也正是我前面说的“通量守恒”的由来。到了这一步每个控制体被拆成了它的各个面对流项和扩散项都变成面上通量的加和。把这个积分形式的方程写成离散形式就是(ρP φP - ρP_old φP_old) VP / Δt Σ_f (ρf vf φf)·Sf Σ_f Γf (∇φ)f·Sf这个方程里f表示控制体的面Sf是这个面的面积矢量P是控制体中心。可以看到要求解的就是控制体中心的值φP但方程里还出现了面上的物理量φf这就引出了下一步界面插值。2.3 第三步界面插值与离散格式选型界面插值是有限体积法里最考验功力的环节之一。控制体中心的值是已知待求的但面上的值如φf、ρf、vf需要从周围控制体中心的值插值得到。这个插值方式就决定了离散格式的精度和稳定性。常用的格式有三类一阶迎风Upwind取上游控制体中心的值作为面上的值。绝对稳定、不会振荡但数值耗散严重会把火焰前锋抹平导致火焰厚度被拉伸预测的火焰温度偏低。中心差分CD取两个相邻控制体中心值的平均。精度高但在强对流主导的燃烧问题里容易出现非物理振荡。二阶迎风 / QUICK在一阶迎风的基础上引入高阶修正兼顾稳定性和精度。QUICK格式对结构化网格上的层流扩散火焰效果不错但对非结构网格的适应性有限。我做燃烧仿真时最常用的组合是湍流场用二阶迎风组分和能量方程用MUSCL带限制器的格式压力用线性格式。这套组合在大多数工况下能做到精度和稳定性的平衡。选格式的时候有几个判断依据要先想清楚火焰附近的Peclet数对流强度/扩散强度是不是远大于1。如果对流占主导用中心差分大概率会振荡。Peclet数大于2时中心差分的精度优势就会变成不稳定的根源。网格够不够细。如果网格足够细高阶格式的误差优势才能体现网格粗的时候强行上QUICK容易算出带“波纹”的非物理解。是否多组分燃烧。组分方程数量多数值耗散会把火焰结构人为加宽多个组分方程的耗散效应叠加最终结果可能完全失真。我个人的一个习惯是先跑一阶迎风确认收敛性和宏观流场特征合理再切到高阶格式做正式计算。这样既能快速定位模型问题又能保证正式算例的精度。3. 时间推进、压力速度耦合与求解器框架3.1 显式还是隐式燃烧仿真的时间推进策略燃烧仿真的时间尺度跨度大。流动的时间尺度可能是毫秒级而化学反应的特征时间在某些机理里可以达到微秒级甚至更小这种刚度问题直接决定了时间推进格式的选型。显式格式如RK4实现简单、每步计算量小但稳定性受CFL条件限制。CFL数超过临界值就会发散。做低速燃烧流动声速远大于流速CFL条件由声速和网格尺寸控制显式格式的时间步需求往往小得离谱计算代价高到难以接受。隐式格式如欧拉隐式、BDF2求解时需要组装矩阵并迭代求解大型代数方程组每一步的代价高但时间步长可以大幅放宽。对燃烧仿真这种典型刚性问题工程计算几乎都走隐式路线再配合自适应时间步长控制。实际用CFD软件做燃烧仿真时还有一个选择会影响时间步长是否做双时间步推进。双时间步法的思路是在物理时间推进的每一层内部用一个“伪时间”做稳态迭代每一物理步内只求解稳态残差。好处是能兼顾物理时间精度和隐式稳定性OpenFOAM里做瞬态燃烧算例常用的PIMPLE算法就是这种思路。3.2 压力-速度耦合SIMPLE算法的核心思路燃烧仿真要同时解动量方程和连续性方程但压力和速度之间没有显式的联立方程需要通过连续性方程来约束。最经典的解法是SIMPLE算法Semi-Implicit Method for Pressure-Linked Equations。SIMPLE的思路分四步用当前压力场猜测速度场得到不满足连续性的“中间速度”。把动量方程代入连续性方程推导出压力修正方程。求解压力修正方程得到压力修正量。用压力修正量更新压力和速度重新进入下一轮迭代。整个过程中最考验数值经验的环节在第2步。推导压力修正方程时要把速度修正量用压力梯度的修正量来表达这一步涉及相邻网格点速度的耦合最后形成的方程是一个标准的拉普拉斯型方程它的系数矩阵特性对角线占优、对称性直接影响求解器的收敛速度。我做贫燃预混燃烧算例时发现SIMPLE算法在反应区附近收敛变慢是常态因为化学反应源项让温度和密度剧烈变化密度场反过来通过连续性方程强烈耦合到压力场上。这种情况下需要在压力修正方程求解时适当增加亚松弛否则压力场会在迭代过程中振荡导致残差曲线“颠簸”不下降。3.3 完整求解流程从方程到代数系统的落地步骤把上面这些串起来一个基于有限体积法的稳态燃烧求解流程大致是生成网格导入求解器确定边界条件和初始场。初始化工况参数如入口速度、温度、组分质量分数、湍流强度。开始外迭代循环求解动量方程得到新的速度场用当前压力场。求解压力修正方程更新压力场和速度场使其满足连续性。求解组分输运方程更新各组分质量分数。求解能量方程更新温度场。根据温度和组分计算新的密度、粘度、扩散系数、化学反应速率。检查残差是否满足收敛判据不满足则继续迭代。输出结果做后处理和验证。实际运行时每轮迭代内部的子方程是分别求解的这就是分离式求解器。每个代数方程组用什么方法来解也会直接影响整体效率我用GAMG比较多处理代数方程组的收敛速度很快。4. 燃烧仿真对有限体积法的特殊要求4.1 多组分输运方程组分数量带来的计算挑战燃烧仿真的组分输运不像单组分标量输运那么简单。以甲烷空气燃烧为例哪怕用简化机理也要涉及几十个组分反应。每个组分都有一个输运方程每个方程的离散形式都在对流项、扩散项、源项上和其它组分相互耦合。离散组分方程时有一个容易忽视的点组分质量分数之和必须恒等于1。数值误差累积会导致“质量不守恒假象”所有组分质量分数求和偏离1。处理办法是——选一个组分作为“剩余组分”不求解它的输运方程而用1减去其它组分之和来得到。我在OpenFOAM里做甲烷燃烧时就是这么处理N2的。如果你把所有组分都解一遍最后还要强制归一化反而容易引发局部振荡。4.2 化学反应源项的刚性问题燃烧方程里最有挑战性的就是化学反应源项的刚性。Arrhenius公式里的指数项对温度极其敏感温度差几十K反应速率就可能是几个数量级的差异。这个源项在离散方程里被放在显式位置时时间步长稍微大一点源项就会主导方程导致发散。一个实用的处理手段是对源项做雅可比线性化。把源项在当前位置做一阶泰勒展开线性部分隐含在方程左端系数里只有高阶残差部分留在显式端。这样源项的负斜率对温度的敏感度变成了左端矩阵的一部分增强了对角占优时间步长就能放宽很多。我在自研求解器里做这一步时踩过一个坑最初只线性化了源项对温度的导数没有线性化对组分浓度的导数。结果计算到点火阶段时组分振荡后来把组分导数也纳入隐式处理后问题就消失了。所以做燃烧仿真源项线性化一定要全面不能只做一半。4.3 湍流与化学反应的耦合燃烧仿真里还有一类难题湍流与化学反应的相互作用。湍流脉动影响当地混合速率化学反应又反过来改变密度场和湍流结构。如果直接用RANS平均方程化学反应源项的非线性会让平均源项不等于源项在平均温度、平均浓度下的值这就是湍流燃烧中的“未封闭项”问题。工程上常用EDM涡耗散模型、EDC涡耗散概念模型或者火焰面模型来处理这个未封闭问题。EDC模型在OpenFOAM里用得比较多它的思路是假设化学反应发生在湍流耗散尺度的微团内微团的体积分数由湍流参数决定化学反应在微团内按0维反应器推进再把微团内的反应结果映射回宏观平均场。从有限体积法的角度看EDC模型带来的是局部微团的反应源项计算它需要单独求解微团内的组分方程而且这个微团方程是常微分方程它的数值稳定性完全取决于前面的源项线性化和时间步长控制能力。所以很多朋友以为“燃烧仿真难在选对机理”其实机理选对了跑到有限体积法求解框架里一卡壳问题还是回到离散和刚性的处理能力上。4.4 火焰面模型和PDF方法的离散视角火焰面模型是另一种常用方法它的核心是假设火焰的结构由湍流混合时间尺度控制火焰内部化学反应时间足够快。它把火焰面做成一维的层流火焰结构库用混合分数作为控制变量。这个思路对非预混火焰尤其有效能够用廉价的输运方程解混合分数而不用解几十个组分方程。从离散实现角度火焰面模型的优势也很明显它把反应源项计算转移到了离线预处理的查表过程在线计算时只需要查插值表源项变成一个代数计算避免了刚性差分问题。代价是增加了表存储空间同时在有限体积法的离散过程中插值表的精度会影响最终结果的一致性。5. 实操记录搭一个有限体积法燃烧算例的关键节点5.1 网格生成与边界条件设置我以甲烷空气非预混射流火焰为例采用OpenFOAM框架。网格选用结构化分区网格燃烧器出口区域加密火焰高度估计在射流直径的十倍左右这个区域内的网格尺寸必须保证火焰前锋处至少有几个网格点否则一阶耗散就会把火焰结构吃掉了。边界条件设置要特别小心入口湍流参数。入口湍流强度设成5%还是2%火焰长度、抬升高度都会有显著差异。我做算例时习惯先用稳态非燃烧等温流动做一次初始化让湍流场先充分发展再开启反应项这比冷启动直接燃烧收敛快很多。5.2 离散格式、亚松弛与时间步长的配置参数离散格式选择上我推荐这套配置作为起步组合动量方程对流项二阶迎风组分方程对流项二阶迎风带Van Leer限制器能量方程对流项二阶迎风压力梯度插值线性格式梯度项高斯线性配合最小二乘法修正压力速度耦合稳态用SIMPLE瞬态用PISO或PIMPLE亚松弛因子的设置同样重要。压力亚松弛设0.3动量设0.7组分和能量设0.8。如果温度场和浓度场强烈耦合比如涉及快速点火温度亚松弛需要降到0.5左右避免每轮迭代中温度大起大落。时间步长的控制上瞬态计算先用固定时间步长试算观察残差和火焰前锋的运动状态再调整到自适应时间步长目标是将每个时间步内的最大库朗数控制在5以内同时保证化学反应时间尺度至少被时间步长覆盖5个点。5.3 收敛判断与结果验证的双重标准收敛判断不能只看残差数值。残差下降三到四个数量级虽然是常规判据但对燃烧问题来说局部火焰区域的残差可能顽固地停在一个平台上这种情况往往是因为化学反应源项和湍流耗散之间的耦合没有完全稳定。我的做法是双参数判断一是全局残差下降曲线进入平台期二是关键物理量出口平均温度、燃烧器出口处OH质量分数随迭代步数的变化率低于某个阈值。温度场每百步变化不超过0.5K才认为收敛。结果验证上也有一套流程不能只看温度云图顺眼就交差对比中轴线温度分布与实验数据对比火焰长度以OH或CH*分布判断与实验值检查组分质量分数和是否等于1检查出口质量流率与入口是否守恒偏差小于0.1%6. 常见问题与排查技巧实录6.1 算例发散的第一排查顺序燃烧算例发散是最常见的绝大多数人第一反应是调小时间步长但很多时候问题并不在时间步长。我建议按这个顺序排查先查边界条件入口压力或温度设错导致密度场突变。再查初始场如果初始温度场或组分场远离物理状态化学反应源项一开始就爆炸。然后查网格质量尤其关注是否有负体积、高度偏斜的网格。最后才看离散格式和亚松弛。其中初始场的问题最隐蔽尤其是点火模拟常常需要先给一个高温点火区域但点火区域大小和温度设定不当会在初始几步产生极高的反应速率直接把流场撕裂。我的习惯是把点火区域温度迭代过程中渐进提升而不是一步加到绝热火焰温度。6.2 组分浓度出现负值的处理有限体积法离散组分方程时如果对流项离散格式选择不当或时间步长过大即使物理上不可能的负组分浓度也会蹦出来。负组分一出现Arrhenius公式代入负浓度会导致反应速率为负进一步加剧震荡最终发散。处理方案通常有三板斧降低时间步长让每一步的通量变化限制在稳定区间内。更换限制器从无限制器切换到Van Leer或Minmod。如果程序中支持开启组分浓度的非负限制。但要注意限制器本身也是一种数值耗散加多了会让火焰位置偏移。能不加尽量不加靠网格和时间步长来保证稳定性是更健康的路径。6.3 温度振荡与燃烧不完全的调试经验还有一类问题是计算不发散但温度场始终在小范围波动或者火焰温度明显偏低。温度场波动往往与能量方程和组分方程的耦合时间滞后有关燃烧释放的热量在本时间步更新但组分输运已经提前消耗了燃料导致热量释放与组分消耗失配。这时需要检查组分方程和能量方程的源项是否同步更新如果求解器是逐个方程串行求解可以增加子迭代次数让它们收敛到一致状态。火焰温度偏低则要先区分是物理原因还是数值原因。物理原因包括机理不完善、热损失边界条件设置不当数值原因包括网格太粗导致的数值耗散、湍流燃烧模型参数不合适。我遇到过把EDC模型的C_eps系数从默认值调低后火焰温度显著上升的情况这个系数直接控制化学反应发生的时间尺度占比对热释放预测非常敏感。7. 几个值得长期留意的实操心得做了这些算例之后我一直觉得有限体积法对做燃烧仿真的人来说不只是“会用软件”层面的知识它更像是调试工具。当你看到残差不降、温度场失真、组分负值这些问题时能够回溯到离散格式、通量计算、源项线性化这些层面去理解根因而不是盲目调整参数这种能力比会用多少种燃烧模型更值钱。一个具体的体会是有限体积法里所有高深技巧的起点都是“通量守恒”这朴素的四个字。燃烧室入口多少质量出口就必须多少质量入口多少化学焓出口就必须带走多少焓除了壁面散热。只要守恒性出错火焰温度和位置一定错。因此我每一轮计算结束后都会先做进出口守恒校验再做结果分析。还有一个小技巧基于有限体积法编写自己的求解器代码时把一个简单标量输运方程先跑通再逐步加入组分、反应源项、湍流模型每一步都验证结果不要一下子把所有物理过程全部耦合进去。我发现把Stepwise验证贯彻到底做下来的代码稳定性远好于一次性堆叠所有模块的实现。最后多说一句做燃烧仿真不能只停留在“会操作”要把有限体积法的数学原理吃透遇到问题才有真正的底气去排查和优化。这篇内容的很多经验也是我从一次次发散、负浓度、温度振荡中积累出来的希望分享出来能让大家少走点弯路。
有限体积法在燃烧仿真中的核心作用与实操要点解析
做燃烧仿真这行也有几年了从最开始照着教程跑算例到后来自己改求解器、调离散格式回头再看整个知识体系发现最难啃、也最绕不开的一块硬骨头就是有限体积法。特别是做燃烧仿真不管是非预混火焰还是预混燃烧不管是RANS还是LES底层全是有限体积法在撑着。很多新手一上来就扑向湍流模型、化学反应机理结果算出来发散或者结果明显不对回头看往往都是栽在了离散格式、网格质量、时间推进这些有限体积法的基础环节上。这篇东西我把自己对有限体积法在燃烧仿真中的应用理解做个系统性整理从控制方程怎么变成代数方程到界面插值怎么选再到燃烧仿真特有的刚性问题怎么处理最后把实际算例中踩过的坑也一并列出来希望能给正在学燃烧仿真或者准备自己写求解器的朋友一点参考。1. 燃烧仿真到底在解什么先把方程组底牌摊开1.1 从连续介质到控制方程燃烧仿真的数学底盘燃烧现象在宏观尺度上本质上是流动、传热、传质和化学反应的耦合过程。我们做数值仿真并不是在模拟每一个分子而是在连续介质假设下把燃烧过程用一组偏微分方程描述出来然后用数值方法去求解。这组方程的核心包括连续性方程、动量方程N-S方程、能量方程以及组分输运方程。拿组分输运方程来说它的基本形式是∂(ρYk)/∂t ∇·(ρvYk) ∇·(ρDk∇Yk) ωk这个方程里左边第一项是组分质量分数的瞬态变化第二项是对流项右边第一项是扩散项第二项是化学反应源项。对整个燃烧仿真来说这串方程组合到一起描述的就是“流动带着组分跑扩散让组分混合化学反应把一种组分变成另一种组分”这个过程。我遇到不少朋友问为什么燃烧仿真一定要用有限体积法用有限差分或者有限元行不行这个问题其实背后有个很实际的原因。燃烧问题里经常出现强间断、大梯度尤其是火焰前锋附近温度和组分浓度变化非常剧烈。有限体积法是基于积分形式的方程它在控制体上做通量守恒这一点对保证“质量守恒、能量守恒”有先天的优势。燃烧数值结果如果守恒性不好算出来的火焰温度、组分分布都会失真。1.2 为什么偏偏是有限体积法三种离散方法横评CFD里有三大类离散思路有限差分法FDM、有限元法FEM、有限体积法FVM。很多人听说过名字但不清楚它们之间的实质差异这里我拿它们做燃烧仿真的效果和适用性来对比一下。有限差分法是最早发展起来的直接在网格点上用差分公式代替偏导数。它的数学推导简单直观计算效率高适合结构网格上的简单几何。但缺点也很明显要求网格足够规则对复杂几何边界难以处理而且在守恒性上需要额外注意。燃烧仿真里的火焰面、旋流燃烧室、多孔介质这类复杂几何用FDM会很吃力。有限元法在固体力学领域是霸主它的思路是把求解域分成有限个单元在每个单元上用试函数逼近真实解。它对复杂几何的适应性很好但当处理含对流占优的流动问题时经典的伽辽金格式会产生数值振荡需要额外的稳定化处理应用在燃烧这类强非线性问题上实现复杂度偏高。有限体积法的思想最贴近物理本质直接在控制体上做积分利用高斯散度定理把体积分转化为面积分把控制方程变成“通过控制体边界的净通量 内部源项 控制体内物理量变化率”这种形式。这种做法的最大优势就是局部守恒性天然满足不管网格多稀疏、多扭曲通量守恒在离散层面严格成立。对燃烧仿真来说在燃烧室内质量、动量、能量守恒是结果的底线有限体积法从根上保证了这一点。拿一个我实际用过的例子来说之前算一个带钝体稳焰的燃烧室网格质量不算好近壁面处扭曲比较大如果是有限差分格式早就开始振荡发散但有限体积法配合合适的限制器结果依然稳定温度场和实验数据对得上。正因为这样主流燃烧仿真工具OpenFOAM、Fluent、CFX包括很多自研求解器底层清一色是有限体积法。理解FVM就是理解这些工具求解机制的钥匙。2. 有限体积法的核心操作积分、散度定理与界面插值2.1 第一步在控制体上做积分有限体积法的起点是在一个控制体上对控制方程做积分。把计算域划分成一个个互不重叠的控制体可以理解为一个个小单元每个控制体中心有一个要计算的物理量值比如速度、压力、温度、组分质量分数。控制体可以是结构化网格的六面体也可以是非结构化的四面体、多面体。以不带源项的标量输运方程为例∂(ρφ)/∂t ∇·(ρvφ) ∇·(Γ∇φ)在一个控制体V上做体积分得到∫_V ∂(ρφ)/∂t dV ∫_V ∇·(ρvφ) dV ∫_V ∇·(Γ∇φ) dV这里逻辑很简单我们想知道控制体内总物理量怎么变化就需要把控制体内每一点的瞬态项、对流项、扩散项加起来。问题是控制体内各点的场是未知的连续函数这个积分没法解析算出来所以才需要数值离散。这里有个关键细节计算时我们并不是把体积分直接“求出来”而是通过高斯散度定理把它变换成面积分再用数值积分近似。这也是为什么我一直强调理解散度定理是理解FVM的关卡。2.2 第二步用高斯散度定理把体积分变成面积分高斯散度定理是有限体积法的理论支点。它说的是一个矢量场在控制体上的体积分等于这个矢量场在控制体表面上的面积分。∫_V ∇·F dV ∮_S F·n dS用这个定理上面的输运方程就变成∫_V ∂(ρφ)/∂t dV ∮_S (ρvφ)·n dS ∮_S (Γ∇φ)·n dS翻译成物理语言控制体内物理量的变化率等于通过控制体表面流入/流出的净对流通量加上净扩散通量。这也正是我前面说的“通量守恒”的由来。到了这一步每个控制体被拆成了它的各个面对流项和扩散项都变成面上通量的加和。把这个积分形式的方程写成离散形式就是(ρP φP - ρP_old φP_old) VP / Δt Σ_f (ρf vf φf)·Sf Σ_f Γf (∇φ)f·Sf这个方程里f表示控制体的面Sf是这个面的面积矢量P是控制体中心。可以看到要求解的就是控制体中心的值φP但方程里还出现了面上的物理量φf这就引出了下一步界面插值。2.3 第三步界面插值与离散格式选型界面插值是有限体积法里最考验功力的环节之一。控制体中心的值是已知待求的但面上的值如φf、ρf、vf需要从周围控制体中心的值插值得到。这个插值方式就决定了离散格式的精度和稳定性。常用的格式有三类一阶迎风Upwind取上游控制体中心的值作为面上的值。绝对稳定、不会振荡但数值耗散严重会把火焰前锋抹平导致火焰厚度被拉伸预测的火焰温度偏低。中心差分CD取两个相邻控制体中心值的平均。精度高但在强对流主导的燃烧问题里容易出现非物理振荡。二阶迎风 / QUICK在一阶迎风的基础上引入高阶修正兼顾稳定性和精度。QUICK格式对结构化网格上的层流扩散火焰效果不错但对非结构网格的适应性有限。我做燃烧仿真时最常用的组合是湍流场用二阶迎风组分和能量方程用MUSCL带限制器的格式压力用线性格式。这套组合在大多数工况下能做到精度和稳定性的平衡。选格式的时候有几个判断依据要先想清楚火焰附近的Peclet数对流强度/扩散强度是不是远大于1。如果对流占主导用中心差分大概率会振荡。Peclet数大于2时中心差分的精度优势就会变成不稳定的根源。网格够不够细。如果网格足够细高阶格式的误差优势才能体现网格粗的时候强行上QUICK容易算出带“波纹”的非物理解。是否多组分燃烧。组分方程数量多数值耗散会把火焰结构人为加宽多个组分方程的耗散效应叠加最终结果可能完全失真。我个人的一个习惯是先跑一阶迎风确认收敛性和宏观流场特征合理再切到高阶格式做正式计算。这样既能快速定位模型问题又能保证正式算例的精度。3. 时间推进、压力速度耦合与求解器框架3.1 显式还是隐式燃烧仿真的时间推进策略燃烧仿真的时间尺度跨度大。流动的时间尺度可能是毫秒级而化学反应的特征时间在某些机理里可以达到微秒级甚至更小这种刚度问题直接决定了时间推进格式的选型。显式格式如RK4实现简单、每步计算量小但稳定性受CFL条件限制。CFL数超过临界值就会发散。做低速燃烧流动声速远大于流速CFL条件由声速和网格尺寸控制显式格式的时间步需求往往小得离谱计算代价高到难以接受。隐式格式如欧拉隐式、BDF2求解时需要组装矩阵并迭代求解大型代数方程组每一步的代价高但时间步长可以大幅放宽。对燃烧仿真这种典型刚性问题工程计算几乎都走隐式路线再配合自适应时间步长控制。实际用CFD软件做燃烧仿真时还有一个选择会影响时间步长是否做双时间步推进。双时间步法的思路是在物理时间推进的每一层内部用一个“伪时间”做稳态迭代每一物理步内只求解稳态残差。好处是能兼顾物理时间精度和隐式稳定性OpenFOAM里做瞬态燃烧算例常用的PIMPLE算法就是这种思路。3.2 压力-速度耦合SIMPLE算法的核心思路燃烧仿真要同时解动量方程和连续性方程但压力和速度之间没有显式的联立方程需要通过连续性方程来约束。最经典的解法是SIMPLE算法Semi-Implicit Method for Pressure-Linked Equations。SIMPLE的思路分四步用当前压力场猜测速度场得到不满足连续性的“中间速度”。把动量方程代入连续性方程推导出压力修正方程。求解压力修正方程得到压力修正量。用压力修正量更新压力和速度重新进入下一轮迭代。整个过程中最考验数值经验的环节在第2步。推导压力修正方程时要把速度修正量用压力梯度的修正量来表达这一步涉及相邻网格点速度的耦合最后形成的方程是一个标准的拉普拉斯型方程它的系数矩阵特性对角线占优、对称性直接影响求解器的收敛速度。我做贫燃预混燃烧算例时发现SIMPLE算法在反应区附近收敛变慢是常态因为化学反应源项让温度和密度剧烈变化密度场反过来通过连续性方程强烈耦合到压力场上。这种情况下需要在压力修正方程求解时适当增加亚松弛否则压力场会在迭代过程中振荡导致残差曲线“颠簸”不下降。3.3 完整求解流程从方程到代数系统的落地步骤把上面这些串起来一个基于有限体积法的稳态燃烧求解流程大致是生成网格导入求解器确定边界条件和初始场。初始化工况参数如入口速度、温度、组分质量分数、湍流强度。开始外迭代循环求解动量方程得到新的速度场用当前压力场。求解压力修正方程更新压力场和速度场使其满足连续性。求解组分输运方程更新各组分质量分数。求解能量方程更新温度场。根据温度和组分计算新的密度、粘度、扩散系数、化学反应速率。检查残差是否满足收敛判据不满足则继续迭代。输出结果做后处理和验证。实际运行时每轮迭代内部的子方程是分别求解的这就是分离式求解器。每个代数方程组用什么方法来解也会直接影响整体效率我用GAMG比较多处理代数方程组的收敛速度很快。4. 燃烧仿真对有限体积法的特殊要求4.1 多组分输运方程组分数量带来的计算挑战燃烧仿真的组分输运不像单组分标量输运那么简单。以甲烷空气燃烧为例哪怕用简化机理也要涉及几十个组分反应。每个组分都有一个输运方程每个方程的离散形式都在对流项、扩散项、源项上和其它组分相互耦合。离散组分方程时有一个容易忽视的点组分质量分数之和必须恒等于1。数值误差累积会导致“质量不守恒假象”所有组分质量分数求和偏离1。处理办法是——选一个组分作为“剩余组分”不求解它的输运方程而用1减去其它组分之和来得到。我在OpenFOAM里做甲烷燃烧时就是这么处理N2的。如果你把所有组分都解一遍最后还要强制归一化反而容易引发局部振荡。4.2 化学反应源项的刚性问题燃烧方程里最有挑战性的就是化学反应源项的刚性。Arrhenius公式里的指数项对温度极其敏感温度差几十K反应速率就可能是几个数量级的差异。这个源项在离散方程里被放在显式位置时时间步长稍微大一点源项就会主导方程导致发散。一个实用的处理手段是对源项做雅可比线性化。把源项在当前位置做一阶泰勒展开线性部分隐含在方程左端系数里只有高阶残差部分留在显式端。这样源项的负斜率对温度的敏感度变成了左端矩阵的一部分增强了对角占优时间步长就能放宽很多。我在自研求解器里做这一步时踩过一个坑最初只线性化了源项对温度的导数没有线性化对组分浓度的导数。结果计算到点火阶段时组分振荡后来把组分导数也纳入隐式处理后问题就消失了。所以做燃烧仿真源项线性化一定要全面不能只做一半。4.3 湍流与化学反应的耦合燃烧仿真里还有一类难题湍流与化学反应的相互作用。湍流脉动影响当地混合速率化学反应又反过来改变密度场和湍流结构。如果直接用RANS平均方程化学反应源项的非线性会让平均源项不等于源项在平均温度、平均浓度下的值这就是湍流燃烧中的“未封闭项”问题。工程上常用EDM涡耗散模型、EDC涡耗散概念模型或者火焰面模型来处理这个未封闭问题。EDC模型在OpenFOAM里用得比较多它的思路是假设化学反应发生在湍流耗散尺度的微团内微团的体积分数由湍流参数决定化学反应在微团内按0维反应器推进再把微团内的反应结果映射回宏观平均场。从有限体积法的角度看EDC模型带来的是局部微团的反应源项计算它需要单独求解微团内的组分方程而且这个微团方程是常微分方程它的数值稳定性完全取决于前面的源项线性化和时间步长控制能力。所以很多朋友以为“燃烧仿真难在选对机理”其实机理选对了跑到有限体积法求解框架里一卡壳问题还是回到离散和刚性的处理能力上。4.4 火焰面模型和PDF方法的离散视角火焰面模型是另一种常用方法它的核心是假设火焰的结构由湍流混合时间尺度控制火焰内部化学反应时间足够快。它把火焰面做成一维的层流火焰结构库用混合分数作为控制变量。这个思路对非预混火焰尤其有效能够用廉价的输运方程解混合分数而不用解几十个组分方程。从离散实现角度火焰面模型的优势也很明显它把反应源项计算转移到了离线预处理的查表过程在线计算时只需要查插值表源项变成一个代数计算避免了刚性差分问题。代价是增加了表存储空间同时在有限体积法的离散过程中插值表的精度会影响最终结果的一致性。5. 实操记录搭一个有限体积法燃烧算例的关键节点5.1 网格生成与边界条件设置我以甲烷空气非预混射流火焰为例采用OpenFOAM框架。网格选用结构化分区网格燃烧器出口区域加密火焰高度估计在射流直径的十倍左右这个区域内的网格尺寸必须保证火焰前锋处至少有几个网格点否则一阶耗散就会把火焰结构吃掉了。边界条件设置要特别小心入口湍流参数。入口湍流强度设成5%还是2%火焰长度、抬升高度都会有显著差异。我做算例时习惯先用稳态非燃烧等温流动做一次初始化让湍流场先充分发展再开启反应项这比冷启动直接燃烧收敛快很多。5.2 离散格式、亚松弛与时间步长的配置参数离散格式选择上我推荐这套配置作为起步组合动量方程对流项二阶迎风组分方程对流项二阶迎风带Van Leer限制器能量方程对流项二阶迎风压力梯度插值线性格式梯度项高斯线性配合最小二乘法修正压力速度耦合稳态用SIMPLE瞬态用PISO或PIMPLE亚松弛因子的设置同样重要。压力亚松弛设0.3动量设0.7组分和能量设0.8。如果温度场和浓度场强烈耦合比如涉及快速点火温度亚松弛需要降到0.5左右避免每轮迭代中温度大起大落。时间步长的控制上瞬态计算先用固定时间步长试算观察残差和火焰前锋的运动状态再调整到自适应时间步长目标是将每个时间步内的最大库朗数控制在5以内同时保证化学反应时间尺度至少被时间步长覆盖5个点。5.3 收敛判断与结果验证的双重标准收敛判断不能只看残差数值。残差下降三到四个数量级虽然是常规判据但对燃烧问题来说局部火焰区域的残差可能顽固地停在一个平台上这种情况往往是因为化学反应源项和湍流耗散之间的耦合没有完全稳定。我的做法是双参数判断一是全局残差下降曲线进入平台期二是关键物理量出口平均温度、燃烧器出口处OH质量分数随迭代步数的变化率低于某个阈值。温度场每百步变化不超过0.5K才认为收敛。结果验证上也有一套流程不能只看温度云图顺眼就交差对比中轴线温度分布与实验数据对比火焰长度以OH或CH*分布判断与实验值检查组分质量分数和是否等于1检查出口质量流率与入口是否守恒偏差小于0.1%6. 常见问题与排查技巧实录6.1 算例发散的第一排查顺序燃烧算例发散是最常见的绝大多数人第一反应是调小时间步长但很多时候问题并不在时间步长。我建议按这个顺序排查先查边界条件入口压力或温度设错导致密度场突变。再查初始场如果初始温度场或组分场远离物理状态化学反应源项一开始就爆炸。然后查网格质量尤其关注是否有负体积、高度偏斜的网格。最后才看离散格式和亚松弛。其中初始场的问题最隐蔽尤其是点火模拟常常需要先给一个高温点火区域但点火区域大小和温度设定不当会在初始几步产生极高的反应速率直接把流场撕裂。我的习惯是把点火区域温度迭代过程中渐进提升而不是一步加到绝热火焰温度。6.2 组分浓度出现负值的处理有限体积法离散组分方程时如果对流项离散格式选择不当或时间步长过大即使物理上不可能的负组分浓度也会蹦出来。负组分一出现Arrhenius公式代入负浓度会导致反应速率为负进一步加剧震荡最终发散。处理方案通常有三板斧降低时间步长让每一步的通量变化限制在稳定区间内。更换限制器从无限制器切换到Van Leer或Minmod。如果程序中支持开启组分浓度的非负限制。但要注意限制器本身也是一种数值耗散加多了会让火焰位置偏移。能不加尽量不加靠网格和时间步长来保证稳定性是更健康的路径。6.3 温度振荡与燃烧不完全的调试经验还有一类问题是计算不发散但温度场始终在小范围波动或者火焰温度明显偏低。温度场波动往往与能量方程和组分方程的耦合时间滞后有关燃烧释放的热量在本时间步更新但组分输运已经提前消耗了燃料导致热量释放与组分消耗失配。这时需要检查组分方程和能量方程的源项是否同步更新如果求解器是逐个方程串行求解可以增加子迭代次数让它们收敛到一致状态。火焰温度偏低则要先区分是物理原因还是数值原因。物理原因包括机理不完善、热损失边界条件设置不当数值原因包括网格太粗导致的数值耗散、湍流燃烧模型参数不合适。我遇到过把EDC模型的C_eps系数从默认值调低后火焰温度显著上升的情况这个系数直接控制化学反应发生的时间尺度占比对热释放预测非常敏感。7. 几个值得长期留意的实操心得做了这些算例之后我一直觉得有限体积法对做燃烧仿真的人来说不只是“会用软件”层面的知识它更像是调试工具。当你看到残差不降、温度场失真、组分负值这些问题时能够回溯到离散格式、通量计算、源项线性化这些层面去理解根因而不是盲目调整参数这种能力比会用多少种燃烧模型更值钱。一个具体的体会是有限体积法里所有高深技巧的起点都是“通量守恒”这朴素的四个字。燃烧室入口多少质量出口就必须多少质量入口多少化学焓出口就必须带走多少焓除了壁面散热。只要守恒性出错火焰温度和位置一定错。因此我每一轮计算结束后都会先做进出口守恒校验再做结果分析。还有一个小技巧基于有限体积法编写自己的求解器代码时把一个简单标量输运方程先跑通再逐步加入组分、反应源项、湍流模型每一步都验证结果不要一下子把所有物理过程全部耦合进去。我发现把Stepwise验证贯彻到底做下来的代码稳定性远好于一次性堆叠所有模块的实现。最后多说一句做燃烧仿真不能只停留在“会操作”要把有限体积法的数学原理吃透遇到问题才有真正的底气去排查和优化。这篇内容的很多经验也是我从一次次发散、负浓度、温度振荡中积累出来的希望分享出来能让大家少走点弯路。
