耗散率提高202倍，温度只降了13% — 层流火焰面的熄火极限

熄火不是温度缓慢下降的结果。在扩散损失超过放热的那一点上，燃烧解整个消失。

扩散增加了，火焰为什么反而熄灭

构建火焰面库时必须手工定下一个网格：标量耗散率 $\chi_{st}$ 要拉到多远才截断表格。截得太短，湍流计算会查到表外；拉得太长，Newton 法会为不存在的解空转。为了亲眼看到那条边界，我直接求解了一维火焰面方程。

先写结果。把 $\chi_{st}$ 从 0.05 提到 10.12 1/s，即 202 倍，最高温度只从 2258.8 K 降到 1960.9 K，降幅 13.2%。再往前走一步，解就没有了。温度并不是缓慢趋于零，而是一直撑着的那支解整个消失。

这听起来像悖论。扩散是把反应物送到火焰面的唯一途径。在非预混火焰中，燃料和氧化剂从两侧而来，除了分子扩散没有别的相遇方式。可是把扩散加大，火却灭了。

把 strain rate a 滑块调高，混合层变薄，同一 $Z$ 上的 $\chi$ 变大。

strain rate a = 60 1/s Z_st = 0.055 

Layer thickness sqrt(2D/a) = 0.82 mm, chi_st = 1.48 1/s against the quench limit 10.12 1/s. Raising a squeezes the layer and the sheet dims; moving Z_st slides the sheet off the chi peak, so the same strain delivers far less dissipation to the reaction zone.

看火焰带变暗的那一侧。把 Z_st 从峰值旁推开，拉伸率不动，反应区接到的耗散却掉下去。

用一个标量写出火焰需要什么条件

假设所有组分的扩散系数相同。在这个假设下，流体微团的混合状态由一个守恒标量决定，即混合分数 $Z$ 。 $Z=0$ 为纯氧化剂， $Z=1$ 为纯燃料。

Peters 在 1984 年提出了以 $Z$ 为自变量的坐标变换。把方程写在 $Z$ 上而不是三维空间坐标上，就只剩下垂直于火焰面的那一个方向。

\rho \frac{\partial Y_k}{\partial t} = \rho \frac{\chi}{2} \frac{\partial^2 Y_k}{\partial Z^2} + \dot{\omega}_k

$Y_k$ 是组分 $k$ 的质量分数， $\dot{\omega}_k$ 是生成率， $\chi$ 是标量耗散率。能量方程具有相同的形式。

\rho \frac{\partial T}{\partial t} = \rho \frac{\chi}{2} \frac{\partial^2 T}{\partial Z^2} - \frac{1}{c_p}\sum_k h_k \dot{\omega}_k

$T$ 是温度， $c_p$ 是定压比热， $h_k$ 是组分焓。Ferreira 后来放宽了 $c_p$ 为常数的条件，提出更一般的形式；当 $c_p$ 为常数时它退回到上式。

这里消失的是流动。没有速度场，没有压力，没有三维网格。湍流对火焰所做的一切都被压进了 $\chi$ 这一个量。同一个 $\chi$ 上温度分成两支的地方，是火焰面表把 1759 K 和 818 K 挂在一起的位置。

标量耗散率究竟是什么的倒数

定义很简单。

\chi = 2D\,\left|\nabla Z\right|^2

$D$ 是扩散系数。量纲是 1/s，可读作扩散特征时间的倒数。

湍流拉伸火焰面时， $Z$ 的梯度变陡。 $\chi$ 正比于梯度的平方，因此涨得更快。进入反应面的分子扩散增加，同时离开反应面的热输运也按同样比例增加。可供反应的驻留时间则缩短。

这就是前面那个悖论的答案。扩散带来反应物，也带走热量。化学反应对温度呈指数依赖，所以失热的一方总是先赢。

对于对冲流层流火焰， $\chi$ 对 $Z$ 的依赖可以写成闭式，这正是 Pitsch 用 DNS 验证过的形式。

\chi(Z) = \chi_{st}\,\frac{\exp\!\left(-2\left[\mathrm{erfc}^{-1}(2Z)\right]^2\right)}{\exp\!\left(-2\left[\mathrm{erfc}^{-1}(2Z_{st})\right]^2\right)}

$\chi_{st}$ 是化学当量混合分数 $Z_{st}$ 处的值，也是整个库的唯一索引。而 $Z_{st}$ 只由燃料与空气的组成决定。

$Z_{st} = \frac{1}{1 + s\,Y_{F,1}/Y_{O,0}}$

$s$ 是烧掉 1 kg 燃料所需的氧气质量， $Y_{F,1}$ 是燃料侧来流质量分数， $Y_{O,0}$ 是氧化剂侧的氧气质量分数。

熄火极限在代码里怎么抓

采用一步不可逆反应 $F + s\,O \to (1+s)P$ ，未知函数就减到一个。由于扩散系数相同， $Y_F + T/q$ 对 $Z$ 是线性的。

Y_F(Z) + \frac{T(Z)}{q} = \frac{T_0}{q} + Z\left(Y_{F,1} + \frac{T_1 - T_0}{q}\right)

$q$ 是单位质量燃料的放热量折算成温度单位。 $Y_O$ 也以同样方式与 $T$ 绑定。于是要解的只有 $T(Z)$ 一个场，一个三对角矩阵。

策略是延拓。在极小的 $\chi_{st}$ 处从 Burke–Schumann 解出发。每步把 $\chi_{st}$ 乘以 1.35 再跑 Newton，并把上一步的解交给下一步作初值。燃烧解消失的那一刻，Newton 要么不收敛，要么掉到常温。在这两个值之间做 40 次几何二分，极限值就出来了。

import math
 
N = 201
H = 1.0 / (N - 1)
ZG = [i * H for i in range(N)]
T_OX, T_FU = 300.0, 300.0      # 两侧来流温度 [K]
S_RATIO = 4.0                  # 每 kg 燃料所需氧气 [kg]
Y_OX = 0.232                   # 空气中的氧气质量分数
Q_HEAT = 36000.0               # 把放热量折算成温度 [K]
T_ACT = 24000.0                # 活化温度 [K]
A_PRE = 2.0e10                 # 指前因子 [1/s]
 
def inv_erfc(y):
    lo, hi = -6.0, 6.0
    for _ in range(80):
        mid = 0.5 * (lo + hi)
        if math.erfc(mid) > y: lo = mid
        else: hi = mid
    return 0.5 * (lo + hi)
 
def chi_profile(chi_st, z_st):
    # Pitsch 的对冲流形式。两端为零，所以加一个下限
    f = lambda z: math.exp(-2.0 * inv_erfc(2.0 * z) ** 2)
    fst = f(z_st)
    out = []
    for z in ZG:
        v = f(min(max(z, 1e-9), 1 - 1e-9)) / fst
        out.append(chi_st * max(v, 1e-4))
    return out
 
def couple_species(z, temp, y_fu):
    # 耦合函数：仅由 T 和 Z 还原 Y_F 与 Y_O
    yf = (T_OX / Q_HEAT + z * (y_fu + T_FU / Q_HEAT - T_OX / Q_HEAT)) - temp / Q_HEAT
    yo = (Y_OX + S_RATIO * T_OX / Q_HEAT
          + z * (S_RATIO * T_FU / Q_HEAT - Y_OX - S_RATIO * T_OX / Q_HEAT)) - S_RATIO * temp / Q_HEAT
    return max(yf, 0.0), max(yo, 0.0)
 
def tridiag_sweep(a, b, c, d):
    n = len(b)
    cp = [0.0] * n; dp = [0.0] * n
    cp[0] = c[0] / b[0]; dp[0] = d[0] / b[0]
    for i in range(1, n):
        m = b[i] - a[i] * cp[i - 1]
        cp[i] = c[i] / m
        dp[i] = (d[i] - a[i] * dp[i - 1]) / m
    x = [0.0] * n
    x[-1] = dp[-1]
    for i in range(n - 2, -1, -1):
        x[i] = dp[i] - cp[i] * x[i + 1]
    return x
 
def newton_flamelet(chi, guess, y_fu, itmax=60):
    T = list(guess)
    n = N - 2
    for _ in range(itmax):
        a = [0.0] * n; b = [0.0] * n; c = [0.0] * n; r = [0.0] * n
        for k in range(n):
            i = k + 1
            yf, yo = couple_species(ZG[i], T[i], y_fu)
            ex = math.exp(-T_ACT / T[i])
            w = A_PRE * yf * yo * ex
            dw = A_PRE * ex * (-yo / Q_HEAT - yf * S_RATIO / Q_HEAT + yf * yo * T_ACT / T[i] ** 2)
            if yf <= 0.0 or yo <= 0.0: w = 0.0; dw = 0.0
            cf = chi[i] / (2.0 * H * H)
            r[k] = -(cf * (T[i + 1] - 2 * T[i] + T[i - 1]) + Q_HEAT * w)
            a[k] = cf if k > 0 else 0.0
            b[k] = -2.0 * cf + Q_HEAT * dw
            c[k] = cf if k < n - 1 else 0.0
        d = tridiag_sweep(a, b, c, r)
        lim = max(abs(v) for v in d)
        scale = min(1.0, 200.0 / lim) if lim > 200.0 else 1.0   # 限制温度更新量
        for k in range(n):
            T[k + 1] = max(T[k + 1] + scale * d[k], 280.0)
        if lim < 1e-6:
            return T, True
    return T, False
 
def burke_schumann(y_fu):
    # 无限反应速率极限 —— 延拓的起点
    z_st = 1.0 / (1.0 + S_RATIO * y_fu / Y_OX)
    T = []
    for z in ZG:
        if z <= z_st: T.append(T_OX + (T_FU - T_OX) * z + Q_HEAT * y_fu * z)
        else: T.append(T_OX + (T_FU - T_OX) * z + Q_HEAT * y_fu * z_st * (1 - z) / (1 - z_st))
    return z_st, T
 
def extinction_scan(y_fu, chi0=0.05, growth=1.35, steps=60):
    z_st, T = burke_schumann(y_fu)
    chi_st = chi0
    hist = []
    last_ok = None
    while steps > 0:
        chi = chi_profile(chi_st, z_st)
        Tn, ok = newton_flamelet(chi, T, y_fu)
        tmax = max(Tn)
        if not ok or tmax < 1000.0:
            break
        T = Tn; last_ok = (chi_st, tmax, list(Tn))
        hist.append((chi_st, tmax))
        chi_st *= growth; steps -= 1
    lo, hi = last_ok[0], chi_st
    base = last_ok[2]
    for _ in range(40):                      # 用几何二分法逼近最后一个燃烧解
        mid = math.sqrt(lo * hi)
        Tn, ok = newton_flamelet(chi_profile(mid, z_st), base, y_fu)
        if ok and max(Tn) > 1000.0: lo, base = mid, Tn
        else: hi = mid
    return z_st, hist, lo, base
 
if __name__ == "__main__":
    z_st, hist, chi_q, Tq = extinction_scan(1.0)
    print("Z_st = %.4f   T_ad(Burke-Schumann) = %.1f K" % (z_st, max(burke_schumann(1.0)[1])))
    print("  chi_st [1/s]    T_max [K]   drop [K]")
    tad = hist[0][1]
    for c, t in hist[::3]:
        print("    %9.4f    %8.1f    %7.1f" % (c, t, tad - t))
    print("    %9.4f    %8.1f    %7.1f   <- quench" % (chi_q, max(Tq), tad - max(Tq)))
    print("chi range %.4f -> %.4f = x%.1f,  T_max drop = %.1f %%"
          % (hist[0][0], chi_q, chi_q / hist[0][0], 100.0 * (tad - max(Tq)) / tad))
    print()
    print("  Y_F1     Z_st    T_ad [K]   chi_q [1/s]")
    for yf in (1.0, 0.5, 0.25, 0.125):
        zs, h, cq, _ = extinction_scan(yf)
        print("   %5.3f   %6.4f    %7.1f    %9.3f" % (yf, zs, max(burke_schumann(yf)[1]), cq))

Z_st = 0.0548   T_ad(Burke-Schumann) = 2273.2 K
  chi_st [1/s]    T_max [K]   drop [K]
       0.0500      2258.8        0.0
       0.1230      2251.1        7.7
       0.3027      2245.4       13.4
       0.7447      2233.0       25.8
       1.8322      2210.9       47.9
       4.5079      2164.0       94.8
      10.1232      1960.9      297.9   <- quench
chi range 0.0500 -> 10.1232 = x202.5,  T_max drop = 13.2 %
 
  Y_F1     Z_st    T_ad [K]   chi_q [1/s]
   1.000   0.0548     2273.2       10.123
   0.500   0.1039     2168.8        8.662
   0.250   0.1883     1991.3        3.942
   0.125   0.3169     1719.8        0.561

把 $\chi_{st}$ 旋钮向右推，橙色分支上的点就跟着耗散率走。

chi_st = 0.500 1/s fuel Y_F1 = 1.000 

Z_st = 0.055, T_ad = 2274 K, quench limit chi_q = 10.20 1/s, current T_max = 2273 K (burning). Push chi_st past the red line and the dot falls off the orange branch. Pull it back down — it stays on the floor, because the frozen branch is also a steady solution. Only re-ignite lifts it.

把点沿橙色分支往右推，越过红线，它就掉到底部。此时再把 $\chi_{st}$ 拉回左边，它也不会回升。因为冻结解同样是一个定常解 —— 迟滞正是这条曲线的本质。

注意表格前六行有多平。 $\chi_{st}$ 涨了 90 倍，温度只让出 94.8 K。最后一格里剩下的 203 K 一次性全掉光。这就是 S 形曲线的膝点。

把 $s=4$ 、 $Y_{O,0}=0.232$ 代进 $Z_{st}=1/(1+s Y_{F,1}/Y_{O,0})$ 。 $Y_{F,1}=1$ 时分母是 $1+4/0.232$ ，倒数就是输出里的 $Z_{st}=0.0548$ 。 $Y_{F,1}=0.125$ 时同一式子给出第二张表的 0.3169。Burke–Schumann 的 2273.2 K 和 $\chi_{st}=0.05$ 的 2258.8 K 已经相差 14.4 K。那是无限反应速率极限和第一支燃烧解之间的间隔。

接着 $\chi_{st}$ 从 0.0500 走到 4.5079，温度只从 2258.8 K 降到 2164.0 K，就是那 94.8 K。再到 10.1232 的一格把剩下的 203 K 带到 1960.9 K。表上的 drop 297.9 K 就是这两段之和。经过 0.1230、0.3027、0.7447、1.8322 时，降幅分别是 7.7 K、13.4 K、25.8 K、47.9 K。膝点之前几乎是直线。每次乘 1.35 的延拓会粗粗地跨过最后这一格。40 次几何二分把 $\chi_q$ 收到 10.123，原因就在这里。

稀释燃料为什么让极限垮得更厉害

输出中的第二张表说明了这一点。把燃料质量分数从 1.0 稀释到 0.125，相差 8 倍。绝热火焰温度从 2273.2 K 降到 1719.8 K，降了 24%。而熄火极限 $\chi_q$ 从 10.123 跌到 0.561，跌了 18 倍。

温度降 24% 换来极限降 1800%。原因在 Arrhenius 项。反应速率挂在 $\exp(-T_a/T)$ 上，这里 $T_a = 24000$ K。 $T$ 从 2273 走到 1720，指数从 $-10.6$ 走到 $-14.0$ 。指数差 3.4 在乘积上就是 30 倍。

$Z_{st}$ 从 0.055 移到 0.317 也朝同一方向起作用。在第一个模拟里把 Z_st 向右推就能看到。 $\chi(Z)$ 的峰值在 $Z=0.5$ 附近， $Z_{st}$ 越大，反应区离峰值越近。同样的拉伸率会送来更多耗散。

工程上出现这种组合的场合是废气再循环(EGR)。用惰性气体稀释燃料或氧化剂， $T_{ad}$ 缓缓下降，熄火极限却急剧下降。先出现局部熄火，然后是火焰抬举(lift-off)，最后是吹熄(blow-off)。

查到膝点右边，温度就掉到常温

把 $\chi_{st}$ 网格铺到 $\chi_q=10.123$ ，表是密了，可 4.5079 处的 2164.0 K 到 10.1232 处的 1960.9 K，一格就掉 203 K。这个膝用几个表点来写，才是真正在换的东西。 $\chi_{st} > \chi_q$ 的单元返回什么，由表之外的模型决定。送去冻结解，用熄火概率加权，或者转到非稳态火焰面。

$\chi_{st}$ 本身是湍流模型的输出。索引偏两倍，膝的左右就对调。掉到常温的计算，先看 $\chi_{st}$ 的输入，再看反应机理。造这张表时一次三对角消元的代价，和隐式扩散的 Thomas 消元是同一类活。