散逸率を202倍にしても温度は13%しか下がらなかった

消炎は温度がゆっくり下がって起きるのではありません。拡散損失が発熱を上回る一点で、解そのものが消えます。

拡散が増えたのに、なぜ火炎が消えるのか

フレームレットライブラリを作るとき、格子を一つ手で決めなければなりません。スカラー散逸率 $\chi_{st}$ をどこまで伸ばして表を切るか、です。短く切れば乱流計算が表の外を参照し、長く取れば存在しない解を探してNewton法が空回りします。その境界を目で見るために、1次元フレームレット方程式を直接解いてみました。

結果を先に書きます。 $\chi_{st}$ を0.05から10.12 1/sまで、202倍に増やす間に最高温度は2258.8 Kから1960.9 Kへ、13.2%しか下がりませんでした。そこからもう一歩進めると解が消えました。温度がゆっくり0に向かうのではなく、支えていた解が丸ごと無くなります。

逆説のように聞こえます。拡散は反応物を火炎面へ運ぶ唯一の手段です。非予混合火炎では燃料と酸化剤が反対側から来るので、分子拡散以外に出会う道がありません。その拡散を強めたら火が消えたのです。

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$ を独立変数とする座標変換を提案しました。3次元の空間座標の代わりに $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$ が定数なら上の式に戻ります。

ここで消えたのは流れです。速度場も圧力も3次元格子もありません。乱流が火炎に対して行うことの全部が $\chi$ 一つに圧縮されています。同じ $\chi$ で温度が二股に折れる場所は、フレームレット表が1759 Kと818 Kを同時に掛けるところです。

スカラー散逸率は何の逆数なのか

定義は単純です。

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

$D$ は拡散係数です。次元は1/sで、拡散にかかる特性時間の逆数と読めます。

乱流がフレームレットを引き伸ばすと $Z$ の勾配が大きくなります。 $\chi$ は勾配の2乗に乗るので、さらに速く大きくなります。反応面へ入る分子拡散が増え、同時に反応面から出ていく熱伝達も同じ割合で増えます。反応できる滞留時間は短くなります。

これが先ほどの逆説の答えです。拡散は反応物を連れてきますが、熱も持ち去ります。化学反応は温度に指数的に依存するので、熱を失う側が必ず先に勝ちます。

対向流層流火炎では $\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                  # 燃料1 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の対向流形。両端で0になるので下限を置く
    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):
    # 結合関数: Y_F, Y_OをTとZだけから復元する
    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}$ を左へ戻しても上がってきません。凍結解も定常解だからです。ヒステリシスこそがこの曲線の本質です。

表の最初の6行が平らなことに注目してください。 $\chi_{st}$ が90倍になる間に温度は94.8 Kしか下がっていません。最後の一マスで残りの203 Kが一気に抜けます。これがS字曲線の膝です。

$Z_{st}=1/(1+s Y_{F,1}/Y_{O,0})$ に $s=4$ 、 $Y_{O,0}=0.232$ を入れます。 $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 から 2164.0 へ、94.8 K しか下がりません。10.1232 の一マスが 1960.9 K まで残り 203 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$ を詰める理由です。

燃料を希釈すると、なぜ限界がより大きく崩れるのか

出力の2番目の表がそれを示しています。燃料質量分率を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}$ 自体が乱流モデルの出力です。索引が 2 倍ずれれば膝の左右が反転します。常温へ落ちた計算は、反応機構より先に $\chi_{st}$ の入力を見ます。その表を作る三重対角一回の費用は、陰的拡散の Thomas 消去と同じ場所にあります。

散逸率を202倍にしても温度は13%しか下がらなかった — 層流フレームレットの消炎限界