抗力係数を 9% 削ったら液滴が 35% 遅くなった — 変形と蒸発の補正が入る場所

蒸発補正は $C_D$ を下げますが、減速に効くのは $C_D/d$ です。直径が縮み始めた瞬間、削った係数は戻ってきます。

1983 年、燃える液滴の前で抗力係数が合わなかった#

Renksizbulut と Yuen は、高温気流の中で蒸発する液滴の抗力を測りました。実測された抗力は、球形相関式が予測する値より一貫して低く出ました。液滴が球形を保っていたにもかかわらず、です。表面から噴き出す蒸気が、境界層を外側へ押しやっていたのでした。

同じ時期に、逆向きのずれも報告されていました。高速気流に置かれた液滴は扁平になり、そうなると抗力は球形より大きくなります。O'Rourke と Amsden が 1987 年に TAB モデルとしてまとめたのがそちらです。

そのため、スプレーコードの液滴抗力係数はたいてい三重構造になっています。一つの球形相関式に、変形補正が一つ、蒸発補正が一つ掛かります。この記事では三重を一枚ずつ剥がし、それぞれが軌道をどれだけ動かすかを同じ条件で測ります。結論から言えば、係数を最も大きく削ったモデルが最初に止まります。

三つの相関式は同じ場所を別々の向きに直す

出発点は球形液滴の抗力係数です。Schiller–Naumann 系の相関式を使います。

$C_{D,\mathrm{sph}} = \frac{24}{Re}\left(1 + \tfrac{1}{6} Re^{2/3}\right), \qquad Re < 1000$

$Re = \rho_g u_{\mathrm{rel}} d / \mu_g$ は、液滴径と相対速度から作ったレイノルズ数です。前半の $24/Re$ が Stokes 抗力で、括弧の中が慣性補正にあたります。

ここに掛かる二つの補正は、性格が異なります。

補正	掛かる因子	向き	根拠
なし (球形)	$1$	—	剛体球、表面での物質移動なし
TAB 変形	$1 + 2.632\,y$	抗力増加	前面投影面積が広がる
蒸発 (R–Y)	$(1+B_M)^{-0.2}$	抗力減少	蒸気が境界層を押しのける

$y$ は TAB モデルの無次元変形量です。 $y = 0$ なら球、 $y = 1$ で分裂と判定します。 $B_M$ は Spalding 物質移動数(表面と無限遠の蒸気質量分率の差を無次元化した量)で、蒸発ソース項と $d^2$ 法則で扱ったあの数です。

変形が実際にどう育つのか、下で直接動かしてみましょう。

relative speed25 m/sdiameter40 umgas density5.0

y 0.000y_ss 0.000C_D x1.00We_r 0.00Re 0

Raise the speed and watch the drop flatten into a disc while the amber bar climbs past the grey sphere value. The oscillation overshoots the dashed steady line by almost 2x, so the peak drag is roughly twice what an equilibrium estimate would give. Turn liquid damping off and the ring never settles. Push far enough and y hits 1 — TAB stops and calls it breakup.

相対速度と直径を上げると、液滴が円盤状に潰れて係数のバーがグレーの球形値を超えていきます。観察すべきは、破線で引いた定常値 $y_{ss}$ を実際の $y$ がどれだけ行き過ぎるかです。

変形は前面投影面積ではなく係数に乗る

TAB は液滴を減衰調和振動子として捉えます。気体の動圧が押し、表面張力が戻し、液体の粘性が減衰させます。

\frac{d^2 y}{dt^2} = \frac{C_F}{C_b}\frac{\rho_g u_{\mathrm{rel}}^2}{\rho_l r^2} - \frac{C_k \sigma}{\rho_l r^3} y - \frac{C_d \mu_l}{\rho_l r^2}\frac{dy}{dt}

$r$ は液滴半径、 $\sigma$ は表面張力、 $\mu_l$ は液体粘性です。定数は $C_F = 1/3$ 、 $C_k = 8$ 、 $C_d = 5$ 、 $C_b = 1/2$ を使います。この方程式が分裂まで進む過程は、Weber 数と TAB 分裂で扱いました。

微分項をすべて 0 と置くと、定常変形量が出てきます。

$y_{ss} = \frac{C_F}{C_b C_k}We_r = \frac{We_r}{12}, \qquad We_r = \frac{\rho_g u_{\mathrm{rel}}^2 r}{\sigma}$

ここに実務者がよく踏む落とし穴が一つあります。 $y_{ss}$ をそのまま抗力補正に入れているコードは珍しくありません。ところがこの振動子はほとんど減衰しません。n-ヘプタンの 40 μm 液滴の減衰比は $10^{-2}$ 程度です。静止状態から出発した $y$ は $y_{ss}$ を通り過ぎ、ほぼ二倍まで跳ね上がります。

つまり定常値で計算した抗力は、実際の第一振動の半分にすぎません。問題は係数そのものではなく、その係数を使うタイミングです。

蒸発は境界層を押し広げて係数を削る

Renksizbulut–Yuen の相関式は、蒸発を掛け算の補正として処理します。

$C_D = \frac{24}{Re}\left(1 + 0.2\,Re^{0.63}\right)(1+B_M)^{-0.2}$

指数が $-0.2$ と小さいのが効いています。 $B_M = 0.6$ なら因子は $0.91$ 、 $B_M = 2$ なら $0.80$ です。激しく燃える液滴でも、係数は 20% 程度しか削られません。

原理は膜理論(film theory)と同じです。表面から出ていく蒸気の流束が境界層を厚くします。速度勾配がなだらかになり、壁面せん断が減り、抗力が下がります。熱伝達の Nusselt 数にも同じ形の補正が掛かります。

Python で四つの液滴を同じノズルから撃った#

静止気体に液滴を一つ撃ち込み、1 ms のあいだ追跡します。支配方程式は 1 次元です。

$\frac{du_d}{dt} = -\frac{3}{4}\frac{\rho_g}{\rho_l}\frac{C_D}{d}\,u_d|u_d|$

四つのケースは抗力係数だけが違います。A は球形、B は球形に TAB 変形を掛けたもの、C は蒸発補正のみ、D は蒸発補正に $d^2$ 法則の収縮まで加えたものです。

import math
 
RHO_G, MU_G = 4.98, 3.3e-5      # 700 K, 10 bar の空気
RHO_L, SIG, MU_L = 620.0, 0.0145, 3.0e-4   # n-ヘプタン
D0, U0 = 40e-6, 25.0            # 分裂後の液滴、噴射速度
K_EVAP = 1.0e-6                 # d^2 法則の定数 [m^2/s]
CF, CK, CD_T, CB = 1.0 / 3.0, 8.0, 5.0, 0.5   # TAB 定数
 
def cd_sphere(re):
    if re < 1e-8:
        return 0.0
    return 24.0 / re * (1.0 + re ** (2.0 / 3.0) / 6.0) if re < 1000.0 else 0.424
 
def cd_distorted(re, y):
    return cd_sphere(re) * (1.0 + 2.632 * max(0.0, min(1.0, y)))
 
def cd_blowing(re, bm):
    if re < 1e-8:
        return 0.0
    return 24.0 / re * (1.0 + 0.2 * re ** 0.63) * (1.0 + bm) ** -0.2
 
def tab_oscillate(y, ydot, u, d, dt):
    r = 0.5 * d
    acc = (CF / CB) * RHO_G * u * u / (RHO_L * r * r) \
        - (CK * SIG / (RHO_L * r ** 3)) * y \
        - (CD_T * MU_L / (RHO_L * r * r)) * ydot
    ydot += acc * dt
    y += ydot * dt
    return y, ydot
 
def run_droplet(law, evaporating, bm=0.6, t_end=1.0e-3, dt=5e-8):
    u, x, d, t = U0, 0.0, D0, 0.0
    y, ydot, ymax = 0.0, 0.0, 0.0
    cd_sum, cd0, n = 0.0, None, 0
    while t < t_end and u > 1e-6:
        re = RHO_G * u * d / MU_G
        if law == 'sphere':
            cd = cd_sphere(re)
        elif law == 'tab':
            y, ydot = tab_oscillate(y, ydot, u, d, dt)
            ymax = max(ymax, y)
            cd = cd_distorted(re, y)
        else:
            cd = cd_blowing(re, bm)
        if cd0 is None:
            cd0 = cd
        cd_sum += cd
        n += 1
        u += -0.75 * RHO_G * cd * u * u / (RHO_L * d) * dt
        x += u * dt
        if evaporating:
            d = math.sqrt(max(1e-18, d * d - K_EVAP * dt))
        t += dt
    return dict(x_mm=x * 1e3, u=u, d_um=d * 1e6, cd0=cd0,
                cd_mean=cd_sum / n, ymax=ymax)
 
CASES = [
    ('A sphere            ', 'sphere', False),
    ('B sphere x TAB      ', 'tab', False),
    ('C blowing (B_M=0.6) ', 'blow', False),
    ('D blowing + d2-law  ', 'blow', True),
]
 
print('Re0 = %.0f   We_r = %.2f   t = 1.0 ms' %
      (RHO_G * U0 * D0 / MU_G, RHO_G * U0 * U0 * 0.5 * D0 / SIG))
print('case                  Cd(t=0)  <Cd>    u(1ms)   x(1ms)   d(1ms)')
base = None
for name, law, ev in CASES:
    r = run_droplet(law, ev)
    if base is None:
        base = r
    print('%s  %6.3f  %6.3f  %6.2f   %6.3f   %5.1f  (%+.1f%% x)' %
          (name, r['cd0'], r['cd_mean'], r['u'], r['x_mm'], r['d_um'],
           100.0 * (r['x_mm'] / base['x_mm'] - 1.0)))
 
tab = run_droplet('tab', False)
print('TAB peak distortion y_max = %.3f  ->  Cd multiplier %.2fx' %
      (tab['ymax'], 1.0 + 2.632 * tab['ymax']))

Re0 = 151   We_r = 4.29   t = 1.0 ms
case                  Cd(t=0)  <Cd>    u(1ms)   x(1ms)   d(1ms)
A sphere               0.910   1.708    3.36    9.261    40.0  (+0.0% x)
B sphere x TAB         0.910   2.233    2.66    7.490    40.0  (-19.1% x)
C blowing (B_M=0.6)    0.828   1.537    3.68    9.732    40.0  (+5.1% x)
D blowing + d2-law     0.828   2.096    2.18    8.868    24.5  (-4.3% x)
TAB peak distortion y_max = 0.632  ->  Cd multiplier 2.66x

変形補正は貫通距離を 19.1% 縮めます。第一振動で係数が 2.66 倍まで上がるためです。蒸発補正だけを掛けると、逆に 5.1% 伸びます。ここまでは各相関式の符号どおりです。

係数が最も低い方が最初に止まった理由

問題は D です。C とまったく同じ抗力係数の式を使っています。初期係数も 0.828 で同じです。ところが 1 ms 後の速度は 2.18 m/s と四つの中で最も低くなっています。A の 3.36 m/s より 35% 遅い値です。

減速の式をもう一度見れば答えが出ます。右辺には $C_D/d$ があります。抗力係数を 9% 削るあいだに、直径は 40 μm から 24.5 μm へ 39% 縮みました。割り算の側がはるかに大きく動いたわけです。

物理的にはこうなります。抗力は面積( $d^2$ )に比例し、慣性は質量( $d^3$ )に比例します。比を取れば $1/d$ です。液滴が痩せるほど、同じ抗力がより大きな減速を生みます。応答時間 $\tau_d \simeq \rho_l d^2 / 18\mu_g$ が $d^2$ で縮む、と言い換えても同じことです。

$\frac{a}{m} \propto \frac{d^2}{d^3} = \frac{1}{d}$

つまり蒸発は、抗力係数を下げる効果と減速を強める効果を同時に持ちます。前者は $-0.2$ の指数に縛られ、後者は直径に直結しています。少し長く計算すれば、後者が勝ちます。

下で四つの液滴を同じノズルから同時に撃ってみましょう。

injection u025 m/sd040 umB_M0.60K1.00 mm2/s

A 0.00 mm (--)B 0.00 mm (-100.0%)C 0.00 mm (-100.0%)D 0.00 mm (-100.0%)

Push B_M up: lane C's drag coefficient drops and it pulls ahead of the grey sphere. Now push K up and lane D uses that same reduced coefficient yet falls behind — the shrinking diameter multiplies the deceleration faster than blowing removes it. Setting K to 0 collapses D onto C, which is the cleanest way to see which of the two effects is doing the work.

$B_M$ を上げると C レーンがグレーの球形を追い抜きます。その状態で $K$ を上げると、同じ係数を使っている D レーンだけが遅れていきます。 $K = 0$ に戻せば D は C にぴったり重なるので、二つの効果のどちらが効いているのかが一度に切り分けられます。

どの相関式をいつ有効にするか

状況	変形補正	蒸発補正	理由
$We_r < 0.5$ 、低温	オフ	オフ	$y_{ss} < 0.04$ 、補正は 4% 未満
ディーゼル噴霧の近傍場	オン	オン	相対速度が最大、変形が支配的
噴霧の遠方場	オフ	オン	減速後に $We$ が下がり、蒸発だけが残る
固体粒子	オフ	オフ	変形も物質移動もない
超臨界近傍	オフ	注意	$\sigma \to 0$ で TAB が無意味・ $B_M$ が発散

最後の行が実際によく問題になります。表面張力が 0 に向かうと TAB の復元項が消え、 $y$ が発散します。その領域では変形補正を切って、拡散界面系のモデルへ移る必要があります。

$We$ が十分低いのに抗力が合わないなら、見るべきは補正ではなく別の場所です。壁面近くの抗力危機と境界層剥離のように、 $Re$ そのものが相関式の有効範囲を外れている場合です。

相関式を一つ変えると検証ケースが丸ごと動く

この計算で貫通距離は -19.1% から +5.1% のあいだを行き来しました。格子も、時間刻みも、乱流モデルも触っていません。液滴一つに掛かる無次元因子を二つ、オンオフしただけです。

ですからスプレー検証で実験とずれたとき、まず格子を詰めるのは順序が違います。先に確認するのは、抗力係数の式が何か、その中に $y$ と $B_M$ がどんな値で入っているか、 $y$ を定常値で入れたのか振動を解いたのか、です。直径が時間とともに縮んでいるなら、 $C_D$ ではなく $C_D/d$ を出力して見ます。係数だけを眺めていては、なぜ液滴が先に止まるのかは最後まで見えてきません。