混合分数稍有波动，平均温度就降 964 K —— presumed β-PDF 与 Z=0 处的无穷大

湍流火焰的平均温度不是平均混合分数处的火焰温度。一旦 β-PDF 的指数 a 小于 1，Z=0 附近的冷空气质量必须用不完全贝塔函数精确积分。

为什么不能用平均混合分数去读火焰温度

设想一个甲烷-空气湍流射流火焰。某一点上，时间平均的混合分数恰好是化学计量值 0.055。这一点的平均温度是多少 K？化学计量下的绝热火焰温度是 2230 K，于是很容易答 2230 K。

答案是 1266 K。这是混合分数波动处于中等水平时的值。本文说明这 964 K 去了哪里。我们亲手积分 presumed β-PDF（假定概率密度呈贝塔分布形状的方法）。最后追踪：为什么用 100 个均匀网格积分会偏差 115 K。

起点是 2018 年夏季燃烧研讨会上整理 FlameMaster 用法的一套幻灯片。在安装步骤、ScanMan 和 FlameMaster 的运行方法之间，夹着几页理论幻灯片：PDF 与 CDF、Beta 函数、混合分数与耗散率、S 曲线。幻灯片只留下了标题。这里补上其中"PDF 与 Beta 函数"这一空白。

Z 负责混合，T(Z) 负责化学#

混合分数 $Z$ （来自燃料流的质量比例）在空气中为 0，在燃料中为 1。若认为反应无限快，温度就只是 $Z$ 的函数。这就是 Burke–Schumann 解：在化学计量 $Z_{st}$ 处取峰值、两侧直线下降的屋顶形状。

在湍流中的某一点， $Z$ 随时间波动。空气团经过时它落到 0 附近，燃料团经过时它跳到较大的值。平均温度是 $T(Z)$ 在全部波动上的平均。

$\tilde T = \int_0^1 T(Z)\,P(Z)\,dZ \neq T(\tilde Z)$

$P(Z)$ 是该点 $Z$ 的概率密度， $\tilde Z$ 是平均混合分数。

如果 $T(Z)$ 是直线，两者相等。但屋顶在峰值处折弯。把平均值放在折点正上方，任何方向的波动都会拉低温度。在下面的模拟中亲自操作一下。

Z̄ = 0.055 s = Z″²/(Z̄(1−Z̄)) = 0.100 particles/frame 6

a = 0.000 · b = 0.00 · a<1: pdf → ∞ at Z=0 · particles 0 · running mean 0 K · exact T̃ 0 K · T(Z̄) 0 K · gap 0 K · P(Z<0.01) 0.0%

每个粒子都是从 β-PDF 中抽取的一个 $Z$ 。它落在屋顶曲线上的高度就是它的温度。粉色点是所有粒子的累计平均，很快贴到绿色线 $\tilde T$ 上。把 $s$ 设为 0.1，三分之一的粒子堆在 $Z<0.01$ 。再把 $\bar Z$ 移到 0.2：屋顶右侧坡度平缓，差距会小得多。

β-PDF：用均值和方差两个量确定形状#

LES 或 RANS 求解器只提供两个量： $\tilde Z$ 和方差 $\widetilde{Z''^2}$ 。没人去求解完整的 PDF。因此形状要事先确定。对于限制在 0 与 1 之间的变量，贝塔分布是自然的选择。

$P(Z) = \frac{Z^{a-1}(1-Z)^{b-1}}{B(a,b)}, \qquad a = \tilde Z\,\gamma,\quad b = (1-\tilde Z)\,\gamma,\quad \gamma = \frac{1}{s} - 1$

$B(a,b)$ 是贝塔函数（归一化常数）。 $s = \widetilde{Z''^2}/[\tilde Z(1-\tilde Z)]$ 是离析强度（segregation intensity）：0 为完全混合，1 为只剩 0 和 1 的完全分离。

在化学计量附近， $\tilde Z$ 很小。所以 $a$ 很快降到 1 以下。 $\tilde Z = 0.055$ 时，只要 $s > 0.052$ 就有 $a<1$ 。此时 $Z^{a-1}$ 在 $Z\to0$ 处趋于无穷大。积分仍然有限。但 PDF 本身的值在端点发散。正是这一性质在后面毁掉了均匀网格积分。

折线可以用不完全贝塔函数精确积分

若 $T(Z)$ 是分段线性的，就不需要数值积分。每一段只需要两个积分。

$\int_0^{x} P\,dZ = I_x(a,b), \qquad \int_0^{x} Z\,P\,dZ = \tilde Z\, I_x(a+1,b)$

$I_x(a,b)$ 是正则化不完全贝塔函数（贝塔分布的累积分布函数）。第二个式子成立，是因为 $Z\cdot Z^{a-1} = Z^{a}$ ，且 $B(a+1,b)/B(a,b) = a/(a+b) = \tilde Z$ 。

以 $Z_{st}$ 为界，分别代入贫燃侧和富燃侧的直线， $\tilde T$ 就有了闭式解。端点的无穷大在 $I_x$ 内部被解析处理。 $I_x$ 用连分式（Lentz 方法）计算，与 Numerical Recipes 的 betai 是同一算法。

火焰面库的情况也一样。标量耗散率与熄火极限中的 $T(Z;\chi)$ 曲线取代屋顶。如果表格的 $Z$ 网格是分段线性插值，逐段套用上面的公式即可。

代码：化学计量 Z̄ 下，方差从 0.01 到 0.8#

固定 $\tilde Z = Z_{st} = 0.055$ ，只改变 $s$ 。除精确积分外，再用"在单元中心取 PDF 值、再把总和重新归一化为 1"的均匀网格积分，分别取 $N=100$ 和 $1000$ 。最后两列是与精确值之差。

import math
 
# Burke-Schumann 甲烷-空气火焰温度 T(Z)
T_U, T_AD, Z_ST = 300.0, 2230.0, 0.055
 
def t_burke(z):
    if z <= Z_ST:
        return T_U + (T_AD - T_U) * z / Z_ST
    return T_U + (T_AD - T_U) * (1.0 - z) / (1.0 - Z_ST)
 
def betacf(a, b, x):
    # 不完全贝塔函数的连分式 (Lentz)
    qab, qap, qam = a + b, a + 1.0, a - 1.0
    c, d = 1.0, 1.0 - qab * x / qap
    d = 1.0 / (d if abs(d) > 1e-300 else 1e-300)
    h = d
    for m in range(1, 300):
        m2 = 2 * m
        aa = m * (b - m) * x / ((qam + m2) * (a + m2))
        d = 1.0 + aa * d; c = 1.0 + aa / c
        d = 1.0 / (d if abs(d) > 1e-300 else 1e-300); h *= d * c
        aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))
        d = 1.0 + aa * d; c = 1.0 + aa / c
        d = 1.0 / (d if abs(d) > 1e-300 else 1e-300)
        dl = d * c; h *= dl
        if abs(dl - 1.0) < 1e-14:
            break
    return h
 
def inc_beta(a, b, x):
    # 正则化不完全贝塔 I_x(a, b)
    if x <= 0.0: return 0.0
    if x >= 1.0: return 1.0
    lnf = (math.lgamma(a + b) - math.lgamma(a) - math.lgamma(b)
           + a * math.log(x) + b * math.log(1.0 - x))
    if x < (a + 1.0) / (a + b + 2.0):
        return math.exp(lnf) * betacf(a, b, x) / a
    return 1.0 - math.exp(lnf) * betacf(b, a, 1.0 - x) / b
 
def beta_ab(zm, s):
    # s = Z''^2 / (Zm (1 - Zm)), 离析强度
    g = 1.0 / s - 1.0
    return zm * g, (1.0 - zm) * g
 
def mean_t_exact(zm, s):
    # T 分段线性，每段用 I_x 积分
    a, b = beta_ab(zm, s)
    p_lo = inc_beta(a, b, Z_ST)                 # P(Z < Z_st)
    z_lo = zm * inc_beta(a + 1.0, b, Z_ST)      # int_0^Zst Z P dZ
    k1 = (T_AD - T_U) / Z_ST
    k2 = (T_AD - T_U) / (1.0 - Z_ST)
    lean = T_U * p_lo + k1 * z_lo
    rich = (T_U + k2) * (1.0 - p_lo) - k2 * (zm - z_lo)
    return lean + rich
 
def mean_t_grid(zm, s, n):
    # 朴素做法：在 n 个单元中心取 pdf 值并重新归一化
    a, b = beta_ab(zm, s)
    num = den = 0.0
    for i in range(n):
        z = (i + 0.5) / n
        p = z ** (a - 1.0) * (1.0 - z) ** (b - 1.0)
        num += p * t_burke(z); den += p
    return num / den
 
zm = Z_ST
print(f"Z_mean = {zm}, T(Z_mean) = {t_burke(zm):.1f} K")
print(" s      a       b      T_exact   gap    grid100  grid1000")
for s in [0.01, 0.03, 0.05, 0.1, 0.2, 0.4, 0.8]:
    a, b = beta_ab(zm, s)
    te = mean_t_exact(zm, s)
    g1, g2 = mean_t_grid(zm, s, 100), mean_t_grid(zm, s, 1000)
    print(f"{s:4.2f} {a:7.3f} {b:8.2f} {te:9.1f} {te - t_burke(zm):7.1f}"
          f" {g1 - te:+8.1f} {g2 - te:+8.1f}")

Z_mean = 0.055, T(Z_mean) = 2230.0 K
 s      a       b      T_exact   gap    grid100  grid1000
0.01   5.445    93.55    1895.7  -334.3     +5.4     -0.0
0.03   1.778    30.55    1663.1  -566.9     -5.9     -0.2
0.05   1.045    17.95    1513.3  -716.7     +3.0     -0.2
0.10   0.495     8.50    1266.2  -963.8   +114.6    +32.3
0.20   0.220     3.78     980.9 -1249.1   +335.4   +167.4
0.40   0.083     1.42     682.6 -1547.4   +585.9   +382.0
0.80   0.014     0.24     391.7 -1838.3   +664.2   +506.9

在浏览器中直接运行上面的代码。

Z̄ = 0.055 N₁ = 100 N₂ = 1000 s-points/step 1

point 0/80 · s = 0.000 · a = 0.000 (singular at Z=0) · T̃_exact 0.0 K · gap to T(Z̄) -2230.0 K · N₁ error +0.0 K · N₂ error +0.0 K

按下 Play， $s$ 以对数间隔从 0.01 扫到 0.8，同时绘出三条曲线。默认参数下，右侧表格给出与上面输出相同的数值。红色区域是 $a<1$ 区间，蓝色误差曲线一越过边界就猛然上升。把 $\bar Z$ 调到 0.2，红色区域向右退， $N_1=100$ 的误差也随之减小。

964 K 和 +115 K 说明了什么#

第一个数字是物理。 $s=0.1$ 时，rms 波动为 0.072，比均值 0.055 还大，在混合层边缘很常见。此处 $\tilde T$ 为 1266 K，比 $T(\tilde Z)$ 低 964 K。即使是 $s=0.01$ 的微小波动，也会损失 334 K。位于峰值之上的平均值对波动最敏感。

把 $s$ 提高到 0.8， $\tilde T$ 降到 392 K。在 $s\to1$ 的极限下，PDF 分裂为 $Z=0$ 和 $Z=1$ 两个 δ 函数。两端都是 300 K。平均温度回到室温，这正是应有的结果。

第二个数字是数值问题。 $s\le0.05$ 时，100 单元网格的误差不超过 6 K，因为 $a$ 大于或接近 1。 $s=0.1$ 时 $a=0.495$ ，误差跳到 +115 K。增加到 1000 单元仍残留 +32 K。 $s=0.8$ 时，即使 1000 单元也偏高 +507 K。网格加密 10 倍，误差只减为 1/3.6。

$a<1$ 时，均匀网格漏掉了冷空气#

符号始终为正，这就是线索。网格积分把平均值算热了，说明冷的一侧质量被少算了。

$s=0.1$ 时，第一个单元 $[0, 0.01]$ 中的真实概率是 $I_{0.01}(a,b) = 0.320$ 。三分之一的空气团都在这一个单元里。在单元中心 $Z=0.005$ 处取 PDF 值，这个单元只分到 0.247。因为 $Z^{a-1}$ 在单元内远高于中心值。缺失的 0.073 经过重新归一化流向了更热的单元。

中点公式的误差本应按 $h^2$ 收敛。但被积函数像 $Z^{a-1}$ 这样奇异时，第一个单元的相对误差与网格无关、保持不变。第一个单元的质量正比于 $h^a$ ，所以总误差也按 $h^a$ 减小。 $a=0.495$ 时，网格加密 10 倍，误差变为 $10^{-0.495}\approx 0.32$ 倍，与输出中的 115 → 32（0.28 倍）基本吻合。

修正方法有两种。一是像上面那样，用 $I_x$ 的差值把每个单元的质量精确放进去。二是把网格按几何级数向 $Z=0$ 加密。前者更便宜、更可靠。另外，当 $s$ 极小、 $\gamma$ 很大时，直接计算 $B(a,b)$ 会溢出，应使用 lgamma 在对数空间计算。

构建火焰面表格时该重新检查 Z 网格的地方#

FlameMaster 幻灯片的最后只有一行：用 ListTool 导出 .kg 文件来建库。当这个库作为 SLF/FPV 表格使用时，β-PDF 积分就是紧接着的一步。火焰面解通常把网格点集中在 $Z_{st}$ 附近。然而 PDF 发散的地方却是 $Z=0$ 。

建好表格后检查一次即可。在 $\tilde Z = Z_{st}$ 、 $s=0.1$ 处，分别用网格积分和 $I_x$ 积分计算 $\tilde T$ 。若两者相差几十 K 以上，这个差值就会原样进入 LES 的平均温度场。而对于反应速率、NO 生成这类对温度呈指数敏感的量，它会放大好几倍。