控制方程与物理建模

组分输运方程：原理与工程设置

组分方程的闭合缺口在多组分扩散、热扩散（Soret）与湍流—化学耦合三处。本文用 Stefan–Maxwell 方程说明 Fick 近似的误差来源，并算出 CH4/空气火焰厚度 57.9 μm、Kolmogorov 尺度 50 μm 与网格下限。

林工程师更新于 2026/9/22
CFDArrhenius工程设置与参数选择物理建模与适用边界控制方程与物理建模组分输运方程质量分数约束Stefan–Maxwell 扩散火焰厚度进阶

组分输运方程：原理与工程设置

组分输运方程的形式与标量输运方程一样，差别全在扩散通量的闭合上。把 ji\mathbf{j}_i 写成 −ρDi∇Yi-\rho D_i\nabla Y_i 是最省事的做法，但它丢掉了三件事：组分之间的交叉扩散、温度梯度驱动的热扩散、以及湍流脉动与化学反应的耦合。这三项的量级并不小——氢燃料火焰中忽略 Soret 效应会让火焰位置偏移数毫米，而层流火焰厚度只有几十微米。下面把三层闭合缺口分别量化。组分输运方程比动量方程多两个坑：一是扩散系数不是常数，选 Fick、恒 Lewis 还是多组分扩散会让结果差几十个百分点；二是反应源项的 Arrhenius 指数对温度极其敏感，活化温度填错一位数就直接把火焰算成惰性混合。设置时的顺序应当是先定扩散模型与物性，再填反应参数，最后用有界格式把质量分数锁在 [0,1][0,1] 内。

组分方程与质量分数约束

第 ii 个组分的质量守恒展开成随体形式

ρ∂Yi∂t+ρu⋅∇Yi=∇⋅(ρDi∇Yi)+ω˙i\rho\frac{\partial Y_i}{\partial t} + \rho\mathbf{u}\cdot\nabla Y_i = \nabla\cdot(\rho D_i\nabla Y_i) + \dot\omega_i

ω˙i\dot\omega_i 是单位体积的组分生成率，量纲 kg/(m3⋅s)\mathrm{kg/(m^3\cdot s)}。与动量方程不同，这组方程有一个硬约束

∑i=1NYi=1\sum_{i=1}^{N} Y_i = 1

逐组分求解时这个约束不会被自动满足，必须靠扩散通量的一致性修正（把 ∑iρDi∇Yi\sum_i \rho D_i\nabla Y_i 从对流通量中扣除）来维持。工程上最简单的检查是输出 ∑iYi\sum_i Y_i 的极值：收敛算例应在 10−610^{-6} 以内偏离 1，若偏离到 10−310^{-3} 量级，说明修正没做或源项没有守恒地对各组分分配。

组分质量守恒的微分形式与闭合缺口

第 ii 个组分的质量守恒写成

∂(ρYi)∂t+∇⋅(ρuYi)=−∇⋅ji+ω˙i\frac{\partial(\rho Y_i)}{\partial t} + \nabla\cdot(\rho\mathbf{u}Y_i) = -\nabla\cdot\mathbf{j}_i + \dot\omega_i

ji\mathbf{j}_i 是扩散质量通量，ω˙i\dot\omega_i 是化学反应源项。方程有 NN 个，但只有 N−1N-1 个独立，因为质量分数之和恒为 1。扩散通量必须满足

∑i=1Nji=0\sum_{i=1}^{N}\mathbf{j}_i = 0

这是 Fick 形式 ji=−ρDi∇Yi\mathbf{j}_i=-\rho D_i\nabla Y_i 不自动满足的：各组分独立扩散会给出非零的净质量通量。修正方法是引入修正速度

ji=−ρDi∇Yi+ρYiuc,uc=∑k=1NDk∇Yk\mathbf{j}_i = -\rho D_i\nabla Y_i + \rho Y_i\mathbf{u}_c,\qquad \mathbf{u}_c = \sum_{k=1}^{N} D_k\nabla Y_k

这一步不是可选的细节。缺少修正时，即使每个组分的方程都"守恒"，总和仍会漂移，表现为质量分数之和偏离 1 并随时间累积。

火焰厚度、Kolmogorov 尺度与网格下限

层流火焰厚度由热扩散率与火焰速度之比给出

δL=αSL\delta_L = \frac{\alpha}{S_L}

甲烷—空气在当量比 1、298 K、1 atm 下的层流火焰速度 $S_L=0.38\ \mathrm{m/s}，代入得，代入得\delta_L=2.2\times10^{-5}/0.38=5.79\times10^{-5}\ \mathrm{m}，即57.9μm；绝热火焰温度约2223K。同样口径下氢气—空气的，即 57.9 μm；绝热火焰温度约 2223 K。同样口径下氢气—空气的S_L=2.1\ \mathrm{m/s}，，\delta_L=1.05\times10^{-5}\ \mathrm{m}，只有10.5μm。这两个数字直接决定网格：预混火焰至少需要10个单元落在火焰厚度内，甲烷要求，只有 10.5 μm。这两个数字直接决定网格：预混火焰至少需要 10 个单元落在火焰厚度内，甲烷要求\Delta x\le5.8\ \mu m，氢气要求，氢气要求\Delta x\le1.0\ \mu m$。

湍流一侧的下限由 Kolmogorov 尺度给出

η=L Re−3/4\eta = L\,Re^{-3/4}

取积分尺度 L=0.05 mL=0.05\ \mathrm{m}、Re=1.0×104Re=1.0\times10^4，得 η=0.05×(104)−0.75=0.05/1000=5.0×10−5 m=50 μm\eta=0.05\times(10^4)^{-0.75}=0.05/1000=5.0\times10^{-5}\ \mathrm{m}=50\ \mu m。于是卡洛维茨数

Ka=δLη=57.950=1.16Ka = \frac{\delta_L}{\eta} = \frac{57.9}{50} = 1.16

Ka<1Ka<1 属于皱褶火焰面区，可以用火焰面模型加湍流输运闭合；Ka>1Ka>1 进入薄反应区，火焰面内部结构被小尺度涡扰动，必须解析内层或使用加厚火焰模型。甲烷算例恰好落在 Ka≈1Ka\approx1 的分界上，说明用火焰面模型是可行的但余量很小，网格一旦放宽到 Δx>60 μm\Delta x>60\ \mu m 就会既解析不了火焰又解析不了 Kolmogorov 涡，两侧的物理同时丢失。

alpha, SL_ch4, SL_h2 = 2.2e-5, 0.38, 2.1        # m2/s, m/s
L, Re = 0.05, 1.0e4
eta = L * Re**-0.75
for name, SL in (("CH4", SL_ch4), ("H2", SL_h2)):
    dL = alpha / SL
    print(f"{name}: delta_L={dL*1e6:6.1f} um  dx_max={dL*1e5:5.2f} um  Ka={dL/eta:5.2f}")
# CH4: delta_L=  57.9 um  dx_max= 5.79 um  Ka= 1.16
# H2:  delta_L=  10.5 um  dx_max= 1.05 um  Ka= 0.21
print(f"eta={eta*1e6:.1f} um")                  # 50.0 um
print(f"Le(H2)={alpha/6.1e-5:.3f}  Le(CH4)={alpha/1.6e-5:.3f}")
# Le(H2)=0.361  Le(CH4)=1.375

注意氢气那一行的 Ka=0.21Ka=0.21 反而更小：火焰更薄但 $S_L更大，更大，\delta_L缩小的速度快于缩小的速度快于Ka$ 的下降，因此氢气火焰在同样的湍流场里反而更容易维持在火焰面区，但它对网格的要求（1.05 μm）比甲烷严格 5.5 倍。

Stefan–Maxwell 与 Fick 的差距

严格的扩散关系来自动量平衡，写成 Stefan–Maxwell 形式

∇Xi=∑j≠iXiXjDij(Vj−Vi)+∇TT∑j≠iXiXjDij(DjTρj−DiTρi)\nabla X_i = \sum_{j\neq i}\frac{X_iX_j}{D_{ij}}\left(\mathbf{V}_j-\mathbf{V}_i\right) + \frac{\nabla T}{T}\sum_{j\neq i}\frac{X_iX_j}{D_{ij}}\left(\frac{D_j^{T}}{\rho_j}-\frac{D_i^{T}}{\rho_i}\right)

第一项是浓度扩散，第二项是热扩散。DijD_{ij} 是二元扩散系数矩阵，DiTD_i^{T} 是热扩散系数。Fick 形式相当于把求和项压成一个有效 DiD_i，在两种情况下误差可以接受：混合物中某一组分浓度远高于其余（稀释近似），或所有二元系数彼此接近（如空气—燃烧产物体系）。误差不能接受的情况是含氢混合物——DH2-H2OD_{\mathrm{H_2\text{-}H_2O}} 与 DH2-N2D_{\mathrm{H_2\text{-}N_2}} 相差三倍以上，用单一 DiD_i 会让氢的优先扩散（preferential diffusion）完全消失，而正是这一效应决定了氢火焰的胞状不稳定性。

Soret 效应与 Lewis 数不能取 1 的场合

热扩散的强弱用热扩散比衡量，量级上它与 Lewis 数直接相关

Le=αD=热扩散率质量扩散率Le = \frac{\alpha}{D} = \frac{\text{热扩散率}}{\text{质量扩散率}}

Le=1Le=1 时热扩散与质量扩散以相同速率进行，火焰结构最简单，这也是"恒 Lewis 数"模型默认取 1 的原因。298 K、1 atm 下甲烷在空气中的 D=1.6×10−5 m2/sD=1.6\times10^{-5}\ \mathrm{m^2/s}，空气热扩散率 α=2.2×10−5 m2/s\alpha=2.2\times10^{-5}\ \mathrm{m^2/s}，Le=1.375Le=1.375；氢气在空气中 D=6.1×10−5 m2/sD=6.1\times10^{-5}\ \mathrm{m^2/s}，Le=2.2×10−5/6.1×10−5=0.361Le=2.2\times10^{-5}/6.1\times10^{-5}=0.361。Le<1Le<1 表示质量扩散快于热扩散，轻组分优先向反应区输运，火焰温度高于绝热值、燃烧速度增大；Le>1Le>1 则相反。LeLe 偏离 1 的幅度超过 20% 时（即 Le<0.8Le<0.8 或 Le>1.2Le>1.2）就必须使用多组分或至少是逐组分 Fick 扩散，否则燃烧速度的预测误差会达到 20%～40%。

Soret 效应的独立判据是热扩散比 kTk_T。对氢气在 300 K 附近 kT≈−0.3k_T\approx-0.3，符号为负意味着氢向冷端富集；对氮气、甲烷这类较重分子 kTk_T 在 10−210^{-2} 量级，可以忽略。因此实用规则是：含 H2\mathrm{H_2} 或 He 的混合物、且存在 100 K 以上温差的算例必须打开热扩散；碳氢燃料—空气体系可以关闭。

边界条件与数值有界性

入口给定质量分数、出口用 inletOutlet（回流时取内部值）、不参与反应的壁面用 zeroGradient。若壁面有催化或吸附，需要换成对应的通量边界，此时壁面反应速率与体相速率用同一套 Arrhenius 参数会更自洽。

数值上，组分方程在反应区极薄、梯度极陡，普通二阶格式会产生负质量分数。做法是把对流项设成 Gauss limitedLinear 1，并在源项很强时启用 MULES 限制器做显式有界修正，同时保留隐式源项做时间推进。判据很直接：输出各 YiY_i 的最小值，出现负值就说明限制器没生效或时间步过大。

参数 取值 依据
DCH4-airD_{\mathrm{CH_4\text{-}air}} 1.6×10−5 m2/s1.6\times10^{-5}\ \mathrm{m^2/s} 298 K、1 atm
ScSc / LeLe 0.948 / 1.375 ν/D\nu/D 与 α/D\alpha/D
预混甲烷质量分数 0.055 化学当量：16.04 g CH₄ 配 275.9 g 空气
AA / β\beta / TaT_a 2.119×10112.119\times10^{11} / 0 / 24380 K 甲烷单步机理
化学时间步 1×10−7 s1\times10^{-7}\ \mathrm{s} τchem=0.108 ms\tau_{chem}=0.108\ \mathrm{ms} 的千分之一
达姆科勒数 139 τflow/τchem\tau_{flow}/\tau_{chem}

反应速率的 Arrhenius 参数怎么填

单步总包反应的速率常数写成

k=A Tβexp⁡(−EaRuT)k = A\,T^{\beta}\exp\left(-\frac{E_a}{R_u T}\right)

AA 的单位随反应级数变化，β\beta 是无量纲温度指数，EaE_a 是活化能，RuR_u 是通用气体常数。甲烷—空气单步机理的常用取值是 A=2.119×1011 m3/(kmol⋅s)A=2.119\times10^{11}\ \mathrm{m^3/(kmol\cdot s)}、β=0\beta=0、Ea=2.027×108 J/kmolE_a=2.027\times10^{8}\ \mathrm{J/kmol}。OpenFOAM 的 irreversibleArrheniusReaction 填的是活化温度 Ta=Ea/RuT_a=E_a/R_u，用 Ru=8314 J/(kmol⋅K)R_u=8314\ \mathrm{J/(kmol\cdot K)} 换算得

Ta=2.027×1088314=2.438×104 KT_a = \frac{2.027\times10^{8}}{8314} = 2.438\times10^{4}\ \mathrm{K}

这一步换算是输入正确性的关键：把 TaT_a 误填成 EaE_a 的数值（2.027e8）会让 exp⁡(−Ea/RuT)\exp(-E_a/R_uT) 直接下溢为零，反应完全不发生；反过来把 EaE_a 当成 TaT_a 填入又会让反应瞬时完成。两个温度的速率对比很能说明指数敏感性：T=1500 KT=1500\ \mathrm{K} 时 Ea/(RuT)=2.027×108/(8314×1500)=16.25E_a/(R_uT)=2.027\times10^8/(8314\times1500)=16.25，k=2.119×1011×e−16.25=2.119×1011×8.716×10−8=1.85×104 m3/(kmol⋅s)k=2.119\times10^{11}\times e^{-16.25}=2.119\times10^{11}\times8.716\times10^{-8}=1.85\times10^{4}\ \mathrm{m^3/(kmol\cdot s)}；T=2000 KT=2000\ \mathrm{K} 时 Ea/(RuT)=12.19E_a/(R_uT)=12.19，k=2.119×1011×5.109×10−6=1.08×106k=2.119\times10^{11}\times5.109\times10^{-6}=1.08\times10^{6}。温度只升 33%，速率放大 58.6 倍。

这个倍数决定了数值策略。设燃料摩尔浓度 CCH4=0.5 kmol/m3C_{\mathrm{CH_4}}=0.5\ \mathrm{kmol/m^3}，化学时间尺度 τchem=1/(kC)=1/(1.85×104×0.5)=1.08×10−4 s\tau_{chem}=1/(kC)=1/(1.85\times10^4\times0.5)=1.08\times10^{-4}\ \mathrm{s}，即 0.108 ms；流动时间尺度取 τflow=L/U=0.3/20=0.015 s\tau_{flow}=L/U=0.3/20=0.015\ \mathrm{s}，达姆科勒数

Da=τflowτchem=0.0151.08×10−4=139Da = \frac{\tau_{flow}}{\tau_{chem}} = \frac{0.015}{1.08\times10^{-4}} = 139

Da≫1Da\gg1 表示反应远快于流动，火焰薄且强耦合，必须用小时间步或隐式化学积分（seulex、EulerImplicit）配合 10−7 s10^{-7}\ \mathrm{s} 量级的初始化学时间步。

// constant/combustionProperties
reactions
{
    methaneAir
    {
        type        irreversibleArrheniusReaction;
        reaction    "CH4 + 2O2 = CO2 + 2H2O";
        A           2.119e11;      // m^3/(kmol s)
        beta        0;
        Ta          24380;         // = Ea/Ru，K；不是 Ea
    }
}

// constant/chemistryProperties
chemistry       on;
initialChemicalTimeStep 1e-7;
odeCoeffs { solver seulex; absTol 1e-12; relTol 1e-4; }

扩散系数模型：Fick、恒 Lewis 与多组分

层流扩散的严格形式是 Stefan–Maxwell 方程组，工程实现中通常退化为三种模型。Fick 模型逐组分给出二元扩散系数 DiD_i，最接近物理但需要 N(N−1)/2N(N-1)/2 个系数；恒 Lewis 模型只给一个 Lewis 数，反算 Di=α/LeD_i=\alpha/Le；恒 Schmidt 模型直接给 ScSc。三者可以通过两个无量纲数互相换算

Sc=νD,Le=αD=ScPrSc = \frac{\nu}{D},\qquad Le = \frac{\alpha}{D} = \frac{Sc}{Pr}

以 298 K、1 atm 下甲烷在空气中的扩散为例，DCH4-air=1.6×10−5 m2/sD_{\mathrm{CH_4\text{-}air}}=1.6\times10^{-5}\ \mathrm{m^2/s}，空气的运动黏度 ν=1.516×10−5 m2/s\nu=1.516\times10^{-5}\ \mathrm{m^2/s}，热扩散率 α=2.2×10−5 m2/s\alpha=2.2\times10^{-5}\ \mathrm{m^2/s}，于是 Sc=1.516×10−5/1.6×10−5=0.948Sc=1.516\times10^{-5}/1.6\times10^{-5}=0.948，Pr=1.516×10−5/2.2×10−5=0.689Pr=1.516\times10^{-5}/2.2\times10^{-5}=0.689，Le=0.948/0.689=1.375Le=0.948/0.689=1.375。若按 OpenFOAM 的 constantLewis 默认取 Le=1Le=1，扩散系数会被设成 α=2.2×10−5 m2/s\alpha=2.2\times10^{-5}\ \mathrm{m^2/s}，比真实值高 37.5%；对甲烷这类 LeLe 接近 1 的燃料误差可接受，但对氢气（Le≈0.3Le\approx0.3）会显著改变火焰位置与燃烧速度，必须改用 Fick 或多组分扩散。

故障模式与判定试验

现象 根因 判定试验
质量分数之和随时间单调偏离 1 扩散通量缺少修正速度，∑ji≠0\sum\mathbf{j}_i\neq0 输出 ∑iYi\sum_i Y_i 的时间序列，修正后应在 10−610^{-6} 内波动
氢火焰燃烧速度比实验低 30% 用 Le=1Le=1 抹掉了优先扩散 输出 LeLe 场并改逐组分 Fick 扩散，比较 $S_L$ 变化
富氢混合物中火焰前锋出现非对称偏移 未打开 Soret 热扩散 打开热扩散复算，前锋位置应发生可测的毫米级移动
网格加密到 10 μm 后火焰速度仍在变 火焰厚度 57.9 μm，10 μm 只给 0.17 个单元/μm 的分辨率不足 按 Δx≤δL/10\Delta x\le\delta_L/10 加密到 5.8 μm 以下
湍流燃烧算例的火焰面模型在 Ka>1Ka>1 区失效 火焰面假设被小尺度涡破坏 由 LL、ReRe 算 η\eta 与 KaKa，超过 1 时换加厚火焰或有限速率模型
燃烧室温度峰值比绝热火焰温度高 200 K 用 Le<1Le<1 燃料时超绝热效应被误判为数值误差 核对 LeLe：Le<1Le<1 时局部超绝热是物理结果，不应做限制器截断
温度场升到 2000 K 但组分完全不变 TaT_a 填成了 EaE_a，指数项下溢 打印 TaT_a 与 Ea/RuE_a/R_u 的比值，应为 1
质量分数之和偏离 1 达 1% 扩散通量未做一致性修正，各组分独立扩散 输出 ∑iYi\sum_i Y_i 极值，应在 10−610^{-6} 内
火焰面上出现负质量分数 对流格式无界或时间步过大 检查 YiY_i 最小值；把格式换成 limitedLinear 1 并把 Δt\Delta t 减半
氢燃料火焰位置比实验靠下游 用 Le=1Le=1 代替真实 Le≈0.3Le\approx0.3，扩散被低估 改成 Fick 扩散并输出 LeLe 场，观察火焰前锋是否上移
反应区温度出现网格相关的尖峰 化学时间步未自适应，刚性方程显式积分失稳 输出化学积分器接受的步长，应小于 10−6 s10^{-6}\ \mathrm{s}
封闭腔内燃料缓慢单调减少 出口边界泄漏，或源项符号写反 积分总体燃料质量随时间，质量守恒应闭合到 10−610^{-6}

参考资料

  1. Poinsot T., Veynante D., Theoretical and Numerical Combustion, 2nd ed., Edwards, 2005.
  2. Williams F.A., Combustion Theory, 2nd ed., Benjamin/Cummings, 1985.
  3. Law C.K., Combustion Physics, Cambridge University Press, 2006.
  4. Peters N., Turbulent Combustion, Cambridge University Press, 2000.
  5. Westbrook C.K., Dryer F.L., “Simplified reaction mechanisms for the oxidation of hydrocarbon fuels in flames”, Combustion Science and Technology, 27:31–43, 1981.
  6. Bird R.B., Stewart W.E., Lightfoot E.N., Transport Phenomena, 2nd ed., Wiley, 2002.
  7. Kee R.J., Rupley F.M., Miller J.A., CHEMKIN-II: A Fortran Chemical Kinetics Package for the Analysis of Gas-Phase Chemical Kinetics, Sandia National Laboratories, SAND89-8009, 1989.
  8. Turns S.R., An Introduction to Combustion, 3rd ed., McGraw-Hill, 2012.