混合分率が少し揺らぐだけで平均温度は964 K下がる — presumed β-PDFと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)$ が直線なら両者は一致します。しかし屋根は頂点で折れています。折れ目の真上に平均を置くと、揺らぎはどちら向きでも温度を下げます。下のシミュレーションで実際に操作してみましょう。

粒子の一つひとつが、β-PDFから引いた $Z$ です。屋根の曲線に落ちた高さがその粒子の温度です。ピンクの点は粒子の累積平均で、すぐに緑の線 $\tilde T$ に張り付きます。$s$ を0.1にすると、粒子の3分の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$ で回します。最後の2列は厳密値との差です。

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

上のコードをブラウザで直接動かしてみましょう。

Playを押すと、$s$ を0.01から0.8まで対数間隔で掃引しながら3本の曲線を描きます。既定値では、右の表が上の出力と同じ値を出します。赤い領域が $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より大きいか、1に近いからです。$s=0.1$ で $a=0.495$ になると、誤差は+115 Kに跳ねます。1000セルに増やしても+32 Kが残ります。$s=0.8$ では1000セルでも+507 Kです。格子を10倍細かくしても、誤差は3.6分の1にしか減りません。

$a<1$のとき、一様格子は冷たい空気を取りこぼす#

符号がいつもプラスであることが、原因を教えてくれます。格子積分は平均を熱くしています。冷たい側の質量を少なく数えているということです。

$s=0.1$ で、最初のセル $[0, 0.01]$ に実際に入っている確率は $I_{0.01}(a,b) = 0.320$ です。空気の塊の3分の1がこの1セルにあります。セル中心 $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 ファイルを書き出してライブラリを作る、という1行です。そのライブラリをSLF・FPVテーブルとして使うとき、β-PDF積分がすぐ次の段階です。フレームレット解は普通、$Z_{st}$ 付近に点を集めます。ところがPDFが発散する場所は $Z=0$ です。

テーブルを作ったら、一度だけ確かめてみてください。$\tilde Z = Z_{st}$、$s=0.1$ で、$\tilde T$ を格子積分と $I_x$ 積分でそれぞれ計算します。二つの値が数十K以上違えば、その差はそのままLESの平均温度場に乗ります。そして反応率やNO生成のように温度に指数的に敏感な量では、数倍に膨らみます。