乱流が火炎を速くする効果が折れ曲がる点 — Damköhler のしわ仮説と bending

乱流火炎速度が直線的に上昇し、折れ曲がり、消炎に至る三つの領域の物理

山火事は風が強いほど速く広がります。風速 5 m/s では、火の手は走る人が追いつけない速さで進みます。ところが風をさらに強めていくと、ある地点で火は自ら消えてしまいます。消防の現場で blow-off と呼ばれる現象です。予混合乱流火炎を実験室のバーナーで測っても、曲線の形はまったく同じになります。

乱流強度 $u'$ を上げると、乱流火炎速度 $S_T$ は最初ほぼ直線的に上昇します。そのあと折れ曲がります。この寝かせ方を bending 効果と呼びます。さらに押すと消炎します。この記事では、その曲線の三つの領域をそれぞれ別の物理で説明します。直線部分は Damköhler のしわ仮説で、折れ曲がりと消炎は Karlovitz 数で、そして火炎が乱流に返す力は baroclinic トルクで捉えます。

なぜ乱流は火炎を速くするのか

メタン-空気予混合火炎の層流火炎速度(静止した混合気の中を火炎面が自ら伝播する速度) $S_L$ は常温常圧で約 0.40 m/s です。火炎厚さ $\delta_L$ は 0.5 mm 程度です。この数字ではエンジンは回りません。ボア 40 mm を 2 ms で横切るには 20 m/s が必要です。 $S_L$ の 50 倍です。

その差を埋めるのが乱流です。だからこそ乱流火炎速度 $S_T$ という量を定義します。定義は反応率の積分から出てきます。平均燃料質量分率の輸送方程式を火炎ブラシの法線方向に積分し、拡散項を打ち消すと、こう残ります。

\rho_u S_T \left(Y_F^u - Y_F^b\right) = -\int_{-\infty}^{+\infty} \overline{\dot{\omega}_F}\ \mathrm{d}n

$\rho_u$ は未燃ガス密度、 $Y_F^u$ と $Y_F^b$ は未燃側・既燃側の燃料質量分率、 $\overline{\dot{\omega}_F}$ は平均燃料反応率、 $n$ は平均火炎ブラシの法線座標です。

押さえておきたいのは、これが層流火炎で使う関係式をそのまま拡張した形だという点です。層流版では $\overline{\dot{\omega}_F}$ の代わりに瞬時反応率が入ります。つまり $S_T$ は速度として定義された量ではありません。単位面積あたりどれだけ燃やせるかを積分した量です。この事実が以降のすべてを説明します。

誤った直感 — 乱流が燃料をよく混ぜてくれる?

いちばんよくある説明はこうです。乱流が燃料と空気をよく混ぜるから反応が速くなる。予混合火炎では誤りです。

予混合という言葉のとおり、燃料と空気は着火前にすでに混ざっています。火炎面の手前のガスはすでに完全な混合気です。乱流がこれ以上混ぜるものは残っていません。混合が速度を支配するのは拡散火炎(燃料と酸化剤が別々に供給される火炎)であって、予混合火炎ではありません。

第二の版もあります。乱流が有効拡散係数を大きくするから火炎が速くなる、という説明です。層流火炎速度は $S_L \sim \sqrt{D/\tau_c}$ でスケールします。ここで $D$ を乱流拡散係数 $u' l_t$ に差し替えると $S_T/S_L \sim \sqrt{u' l_t / D}$ になります。平方根です。ところが実験が示す初期勾配は直線です。指数が合いません。

両者にはもっと深いずれがあります。どちらも火炎が乱流をそのままにしておくと仮定しています。実際は逆です。着火直後、温度は 300 K から 2000 K 超へ跳ね上がります。動粘性係数 $\nu$ はおよそ $T^{1.7}$ で上がるので、20 倍以上に増えます。局所 Reynolds 数 $u' l_t / \nu$ はその分だけ下がります。よく発達した乱流が着火後に層流化しうるということです。乱流が火炎を揺さぶっている間、火炎のほうも乱流を殺しています。

Damköhler の答え — 速度ではなく面積である#

Damköhler は別の絵を出しました。乱流は局所の燃焼を速くしません。火炎面を折りたたみます。

仮定は一つです。しわの寄った火炎面のすべての点が、局所的には依然として $S_L$ で伝播する。化学には手を触れず、幾何だけを変えます。そうすると質量保存が答えを渡してくれます。しわの寄った実面積 $A_T$ を通る未燃ガスと、平均火炎ブラシの投影面積 $\bar{A}$ を通る未燃ガスが等しくなければなりません。

\rho_u S_L A_T = \rho_u S_T \bar{A} \quad\Longrightarrow\quad \Sigma \equiv \frac{A_T}{\bar{A}} = \frac{S_T}{S_L}

$\Sigma$ はしわ係数(wrinkling factor)、すなわち単位投影面積あたりの火炎面積です。左辺は面積比、右辺は速度比です。二つは同じ数です。

しわの大きさは乱流が決めます。サイズ $l_t$ の渦が火炎面を速度 $u'$ で押しのける間、火炎面自身は $S_L$ で押し返します。この二つの変位の比がしわの傾きを決め、傾きが面積の増分を決めます。ここから Damköhler の線形則が出ます。

\frac{S_T}{S_L} = 1 + \frac{u'}{S_L}

$u'$ は未燃ガスの RMS 速度変動、 $S_L$ は層流火炎速度です。 $u' = 0$ なら $S_T = S_L$ に戻ります。

この絵が成り立つには、火炎面が渦より薄くなければなりません。二つの無次元数がその条件を測ります。

Da = \frac{\tau_t}{\tau_c} = \frac{l_t/u'}{\delta_L/S_L}, \qquad Ka = \frac{\tau_c}{\tau_\eta} = \left(\frac{u'}{S_L}\right)^{3/2}\left(\frac{l_t}{\delta_L}\right)^{-1/2}

$Da$ は Damköhler 数(乱流時間 / 化学時間)、 $Ka$ は Karlovitz 数(化学時間 / Kolmogorov 時間)です。 $\tau_t$ は渦の回転時間、 $\tau_c$ は化学反応時間、 $\tau_\eta$ は最小渦の時間スケールです。 $Da \gg 1$ なら火炎は渦が回るより速く反応します。 $Ka \lt 1$ なら最小の渦でさえ火炎厚さの中に入れません。この二つが同時に成り立つ領域が flamelet 領域で、Damköhler の線形則はその中でのみ有効です。

しわを自分で作ってみる

下のシミュレーションでしわを自分の手で作ってみましょう。

u'/S_L3.0l_t/δ_L20

What to watch: the drawn front's measured arc length is forced to equal S_T/S_L, so a faster flame is a literally longer line. Raise u'/S_L and the amber linear-law ghost runs away from the real front — that gap is the bending effect. Lower l_t/δ_L to drive Ka past 25 and the sheet thickens, desaturates, then breaks into holes.

$u'/S_L$ を 0 から 6 まで上げると火炎面が細かく折れますが、測定された $\Sigma$ は 4 付近でそれ以上伸びません。隣の直線的な Damköhler 予測は 7 まで上がっているので、二本の線がはっきり離れます。次に $l_t/\delta_L$ を 2 に下げたまま $u'/S_L$ を 12 まで押してみてください。 $Ka$ が 29 を超えてバッジが quenching に変わり、 $\Sigma$ は 3.5 から 2.0 まで崩れ落ちます。

実験が示した折れ曲がり — bending と消炎#

実験データは三つの領域を描きます。最初はほぼ直線。次に平坦化。最後に消炎。

平坦化する理由は三つ重なります。第一に、しわは生成されるだけでなく消滅もします。火炎面は自らの法線方向に伝播するので、凸のカスプ(cusp)を削り続けます。生成と消滅が釣り合うと $\Sigma$ は飽和します。第二に、 $Ka$ が 1 を超えると小さな渦が予熱帯の中に入り込み、火炎を厚くします。局所の $S_L$ 自体が下がります。第三に、強いひずみ速度が局所的に火炎片を消し、火炎シートに穴を開けます。

この三つを一つの式に圧縮した形を本記事では使います。

\Sigma = \frac{S_T}{S_L} = 1 + q\ \frac{u'/S_L}{1 + (u'/S_L)/\Sigma_\infty}, \qquad q = \exp\left[-\left(\frac{Ka}{Ka_q}\right)^{2}\right]

$\Sigma_\infty$ はしわの飽和上限(ここでは 6.0)、 $Ka_q$ は消炎 Karlovitz 数(ここでは 25.0)、 $q$ は 0 と 1 の間の消炎減衰係数です。分母の飽和項が bending を作り、 $q$ が消炎を作ります。

ここは正直であるべきところです。この式は導出された法則ではありません。観測された曲線の形を再現するように合わせた現象論的な形です。 $\Sigma_\infty = 6.0$ と $Ka_q = 25.0$ は物理定数ではなくチューニング値で、燃料・圧力・バーナー形状が変われば一緒に変わります。出典も同じことを言っています。曲線の折れ曲がりを予測するのは難しく、消炎限界を予測するのは「ほぼ不可能」だと。物理から出てきたのは直線部分だけで、残りの二領域はいまも相関式です。

火炎も乱流を変える — baroclinic トルク#

火炎面を通るとガスは膨張します。質量流束は保存されなければならないので

\rho_u S_L = \rho_b u_b \quad\Longrightarrow\quad \frac{u_b}{S_L} = \frac{\rho_u}{\rho_b} \approx 7

$\rho_b$ は既燃ガス密度、 $u_b$ は火炎面を抜けるガス速度です。火炎面は流れを 7 倍に加速する装置です。この加速がそのまま乱流に働きます。

渦度輸送方程式の四つの項がその経路を見せてくれます。

\frac{\mathrm{D}\boldsymbol{\omega}}{\mathrm{D}t} = \underbrace{(\boldsymbol{\omega}\cdot\nabla)\mathbf{u}}_{\text{I}} - \underbrace{\boldsymbol{\omega}(\nabla\cdot\mathbf{u})}_{\text{II}} + \underbrace{\nu\nabla^{2}\boldsymbol{\omega}}_{\text{III}} + \underbrace{\frac{1}{\rho^{2}}\nabla\rho\times\nabla p}_{\text{IV}}

$\boldsymbol{\omega}$ は渦度ベクトル、 $\mathbf{u}$ は速度、 $\rho$ は密度、 $p$ は圧力です。I は渦伸長、II は膨張(dilatation)、III は粘性散逸、IV は baroclinic トルク(密度勾配と圧力勾配がずれたときに生じる回転力)です。

冷たい乱流では II と IV はゼロです。火炎があると両方とも生き返ります。火炎面では $\nabla \cdot \mathbf{u} \gt 0$ なので、II は渦度を減らします。一方 IV は新しい渦度を作ります。火炎面を横切る $\nabla \rho$ は面に垂直で大きく、加速が作る $\nabla p$ は向きが違います。二つの外積はゼロではありません。これが「flame-generated turbulence」の出どころです。

下のシミュレーションで渦を一つ火炎面に投げ込んでみましょう。

u_θ/S_L2.0r/δ_L10ρ_u/ρ_b6.0

What to watch: at the defaults (u_θ/S_L = 2, r/δ_L = 10) the vortex dies on the front and a violet counter-rotating pair is all that is left. Push u_θ/S_L past ~9 and it punches through, dragging a cusp that stays. The arrows get ρ_u/ρ_b times longer on the burnt side — that jump is the dilatation term; set ρ_u/ρ_b = 1 and the baroclinic pair disappears entirely.

$u_\theta/S_L$ を 2 付近にして Release を押すと、渦は火炎面で壊され、その場所に逆回転の渦度が新しく生まれます。同じ値を 10 付近まで上げると、渦はほぼそのまま通り抜け、深いカスプだけを残します。 $\rho_u/\rho_b$ を 1 に下げると $\nabla\rho$ が消え、baroclinic 生成がまるごと切れることも確認できます。

このフィードバックには厄介な副作用が一つあります。 $S_T(u')$ 曲線でいう「 $u'$ 」がどこの $u'$ なのか曖昧になる点です。未燃ガスで測った値、既燃ガスで測った値、火炎ブラシ全体で測った値は互いに別の数字です。密度重み付けをしないブラシ平均は、未燃・既燃という二状態を行き来する間欠性まで変動として数えてしまい、値が水増しされます。この区別はモデルにあまり反映されません。RANS コードが火炎近傍で非物理的な挙動を見せる理由の一つです。

Python で描く bending 曲線#

上のモデルをそのままコードに移します。二つの乱流スケールで線形予測と bending 予測を並べて出力します。

import math
 
SIGMA_INF = 6.0     # しわの飽和上限
KA_QUENCH = 25.0    # 消炎 Karlovitz 数
 
 
def damkohler_number(u_ratio, l_ratio):
    """Da = 乱流時間 / 化学時間。"""
    return l_ratio / u_ratio
 
 
def karlovitz_number(u_ratio, l_ratio):
    """Ka = 化学時間 / Kolmogorov 時間。"""
    return u_ratio**1.5 * l_ratio**-0.5
 
 
def quench_factor(ka):
    """Ka が大きくなると 0 へ落ちる消炎減衰係数。"""
    return math.exp(-(ka / KA_QUENCH) ** 2)
 
 
def wrinkling_factor(u_ratio, l_ratio):
    """Sigma = S_T / S_L(飽和と消炎を含む)。"""
    ka = karlovitz_number(u_ratio, l_ratio)
    return 1.0 + quench_factor(ka) * u_ratio / (1.0 + u_ratio / SIGMA_INF)
 
 
def turbulent_flame_speed(s_laminar, u_ratio, l_ratio):
    """乱流火炎速度 [m/s]。"""
    return s_laminar * wrinkling_factor(u_ratio, l_ratio)
 
 
U_RATIOS = [0.5, 1, 2, 3, 4, 6, 8, 10, 12]
 
for l_ratio in (10.0, 2.0):
    print(f"l_t/delta_L = {l_ratio:.0f}")
    print(f"{'u/SL':>6}{'Da':>9}{'Ka':>9}{'linear':>9}{'bending':>9}")
    for u_ratio in U_RATIOS:
        da = damkohler_number(u_ratio, l_ratio)
        ka = karlovitz_number(u_ratio, l_ratio)
        lin = 1.0 + u_ratio                      # Damkohler の線形則
        ben = wrinkling_factor(u_ratio, l_ratio)  # 飽和 + 消炎
        print(f"{u_ratio:6.1f}{da:9.2f}{ka:9.2f}{lin:9.2f}{ben:9.2f}")
    print()
 
print("S_T at S_L = 0.40 m/s, u'/S_L = 6, l_t/delta_L = 10 :",
      round(turbulent_flame_speed(0.40, 6.0, 10.0), 3), "m/s")

出力はこうなります。

l_t/delta_L = 10
  u/SL       Da       Ka   linear  bending
   0.5    20.00     0.11     1.50     1.46
   1.0    10.00     0.32     2.00     1.86
   2.0     5.00     0.89     3.00     2.50
   3.0     3.33     1.64     4.00     2.99
   4.0     2.50     2.53     5.00     3.38
   6.0     1.67     4.65     7.00     3.90
   8.0     1.25     7.16     9.00     4.16
  10.0     1.00    10.00    11.00     4.20
  12.0     0.83    13.15    13.00     4.03
 
l_t/delta_L = 2
  u/SL       Da       Ka   linear  bending
   0.5     4.00     0.25     1.50     1.46
   1.0     2.00     0.71     2.00     1.86
   2.0     1.00     2.00     3.00     2.49
   3.0     0.67     3.67     4.00     2.96
   4.0     0.50     5.66     5.00     3.28
   6.0     0.33    10.39     7.00     3.52
   8.0     0.25    16.00     9.00     3.28
  10.0     0.20    22.36    11.00     2.68
  12.0     0.17    29.39    13.00     2.00
 
S_T at S_L = 0.40 m/s, u'/S_L = 6, l_t/delta_L = 10 : 1.559 m/s

二つの列を比べると分かれ目が見えます。 $u'/S_L = 1$ までは線形と bending が 15% 以内で重なります。 $u'/S_L = 6$ ではすでに 7.00 対 3.90 に開きます。そして大きなスケール( $l_t/\delta_L = 10$ )では $u'/S_L = 12$ まで押しても $\Sigma$ が 4 付近を保ちます。小さなスケール( $l_t/\delta_L = 2$ )では同じ $u'$ で $Ka$ が 29 を超え、2.0 まで落ちます。乱流強度が同じでも、スケールが小さいほうが先に消えます。

次に乱流火炎に出会ったら

$S_T$ を速度ではなく面積として読む習慣が半分を解決します。乱流が化学を速くしたのではありません。同じ $S_L$ で燃える面をより広く広げただけです。実験で $S_T$ が伸びないときは、反応速度を疑う前に火炎面の面積がなぜ飽和したのかを先に見ます。

$u'$ 一つで判断しないことが残りの半分です。 $u'/S_L$ が同じでも $l_t/\delta_L$ が違えば $Ka$ は 10 倍変わります。上の表では $u'/S_L = 12$ のとき $\Sigma$ が 4.03 だったり 2.00 だったりしました。分けたのは乱流強度ではなく乱流スケールです。 $Ka$ を併記しない $S_T$ 相関式は半分しか語っていません。

そして計算機が吐き出した bending 曲線を法則として扱わないことです。直線部分は幾何から出ましたが、折れ曲がりと消炎はデータに合わせたものです。RANS の結果が火炎近傍でおかしいときは、乱流モデル定数に手を伸ばすより $u'$ をどこで測っているかを確認するほうが速いです。