洗濯物はなぜ5℃でも乾くのか — 拡散律速の蒸発と熱律速の沸騰

沸騰しない蒸発と沸騰する蒸発が別のソース項を必要とする理由

ベランダに干した洗濯物が乾きます。気温は5℃です。水の沸点は100℃です。95℃も足りないのに、水は消えていきます。

同じ水がやかんの中では100℃まで待ちます。どちらも液体が気体になる相変化です。しかし速度を決めているものが違います。CFDでも同じです。両方に同じソース項を使えば、どちらか一方は必ず外れます。

この記事では、その二つを切り分けます。拡散が律速となる蒸発と、熱が律速となる蒸発をそれぞれ別の式で立て、それをNavier–Stokes方程式のソース項へ移すところまで進みます。最後にLeeモデルの緩和係数 $r$ を自分で選んでみて、文献の $r$ が7桁にわたって散らばっている理由を確かめます。

二つの蒸発はボトルネックが違う

蒸発には二つが必要です。分子が界面を離れるだけのエネルギーと、離れた分子を界面から運び去る経路です。どちらが足りないかでレジームが分かれます。

拡散律速。界面直上の気体はすでに飽和しています。さらに蒸発させるには、その蒸気を遠くへ運ぶ必要があります。ボトルネックは気相側の拡散です。液体温度はほとんど変わりません。水たまり、濡れた洗濯物、常温の燃料噴霧がここに属します。

熱律速。界面近傍の液体がすでに飽和温度を超えています。蒸気の逃げ場は十分にあります。ボトルネックは潜熱(相変化に隠れて費やされる熱)を供給する熱です。沸騰、キャビテーション、壁面凝縮がここに属します。

二つのレジームを測る無次元数も異なります。拡散律速はSpalding物質移動数 $B_M$ (界面と遠方の蒸気濃度差)で、熱律速はJakob数 $\mathrm{Ja} = c_p \Delta T / L$ (顕熱と潜熱の比)で測ります。前者はSchmidt数 $\mathrm{Sc} = \nu/D$ と組み、後者はPrandtl数と組みます。

拡散律速 — 濃度勾配が速度を決める

球形液滴の表面から出る蒸気質量流束を立てます。Fickの法則(濃度勾配に比例する拡散)だけでは足りません。蒸気が出ていくとき、気体全体が外向きに押し出されるからです。この対流成分をStefan流と呼びます。

\dot{m}'' = -\rho_g D \frac{\partial Y_v}{\partial r}\bigg|_{s} + \dot{m}'' Y_{v,s}

$\dot{m}''$ は単位面積あたりの蒸発率、 $\rho_g$ は気体密度、 $D$ は蒸気の拡散係数、 $Y_v$ は蒸気質量分率、添字 $s$ は界面です。右辺第2項がStefan流が再び運んでくる分です。

準定常を仮定して半径方向に積分すると対数が現れます。

\dot{m}'' = \frac{\rho_g D}{r_s}\ln(1 + B_M), \qquad B_M = \frac{Y_s - Y_\infty}{1 - Y_s}

$r_s$ は液滴半径、 $Y_s$ は界面の飽和質量分率、 $Y_\infty$ は遠方の値です。 $B_M$ はSpalding物質移動数 — 蒸発を駆動する濃度差を一つの値に圧縮したものです。

ここに液滴の質量保存 $\dot{m} = -\frac{d}{dt}(\rho_l \pi d^3/6)$ を接続すると、直径の2乗が線形に減ります。

\frac{d(d^2)}{dt} = -K, \qquad K = \frac{8\rho_g D}{\rho_l}\ln(1 + B_M)

これがd²則です。直線になるのは $d$ ではなく $d^2$ です。寿命は $t_{life} = d_0^2 / K$ です。

要点は $Y_s$ が温度だけの関数だという点です。25℃の水の飽和蒸気圧は3.17 kPa、大気圧の3 %です。質量分率では $Y_s \approx 0.019$ 。ゼロではありません。だから乾きます。5℃では $Y_s \approx 0.0054$ まで下がりますが、それでもゼロではありません。沸点はこの話に登場しません。

下のシミュレーションで実際に操作してみてください。

RH30 %T25 °Cd₀1.00 mm

What to watch: the disk radius is d(t)/2 from the same d(t) the cyan trace draws, so the plot and the drop are one object. Push RH from 30 % to 95 % — B_M collapses, the amber line flattens, and the lifetime jumps from minutes to the better part of an hour. Raise T instead and Y_s runs away from Y_∞, so the same drop empties in seconds.

RHスライダーを30 %から95 %へ上げると、 $Y_\infty$ が $Y_s$ に近づいて $B_M$ が崩れ、右側の直線が寝てきます。1 mmの水滴の寿命が5分から1時間近くまで伸びます。逆に温度を上げると $Y_s$ が指数的に離れていき、同じ水滴が数秒で消えます。左の円の半径は、右のグラフと同じ $d(t)$ で描かれています。

熱律速 — 届いた熱が潜熱を支払う

次は界面が飽和温度に張り付いている場合です。蒸発率を決めるのは濃度ではなく、界面に届く熱の不均衡です。

\dot{m}'' L = k_l \frac{\partial T}{\partial n}\bigg|_l - k_v \frac{\partial T}{\partial n}\bigg|_v

$L$ は潜熱、 $k_l$ ・ $k_v$ は液相・気相の熱伝導率、 $n$ は界面法線です。液相側から入ってきた熱のうち気相側へ抜けなかった分が、すべて相変化に回ります。

これをStefan条件と呼びます。厄介な理由がはっきりしています。界面に条件が二つ同時に掛かるのです。温度は $T = T_{sat}$ に固定(Dirichlet)され、同時に流束の跳びが界面速度を決めます。条件が二つあって、境界の位置そのものが未知です。自由境界問題です。

Sharp-interfaceで正直に解くには界面を追跡し、毎ステップ上の式を満たさなければなりません。多くの商用・オープンソースソルバはその道を避けます。代わりに界面を1〜2セル厚に滲ませ、上の条件が自然に満たされるように設計した体積ソース項を入れます。

ソース項はどこに入るのか

体積蒸発率 $\dot{m}'''$ [kg/m³/s] を一つ定義すると、支配方程式の4か所に入ります。

相別連続方程式には符号を反転させた対で入ります。

\frac{\partial (\alpha_l \rho_l)}{\partial t} + \nabla\cdot(\alpha_l \rho_l \mathbf{u}) = -\dot{m}''' , \qquad \frac{\partial (\alpha_v \rho_v)}{\partial t} + \nabla\cdot(\alpha_v \rho_v \mathbf{u}) = +\dot{m}'''

$\alpha$ は体積分率です。二式を足せばソースが打ち消し合い、混合質量は保存されます。

エネルギー方程式には潜熱分のシンクとして入ります。蒸気質量分率を別に輸送する拡散律速ソルバなら、そこにも同じ値がソースとして入ります。

\frac{\partial (\rho e)}{\partial t} + \nabla\cdot(\rho \mathbf{u} h) = \nabla\cdot(k\nabla T) - \dot{m}''' L

最も落としやすいのは4番目です。相変化は速度場の発散拘束を変えます。

\nabla\cdot\mathbf{u} = \dot{m}''' \left( \frac{1}{\rho_v} - \frac{1}{\rho_l} \right)

非圧縮性ソルバの圧力Poisson方程式は通常、右辺がゼロです。相変化があればゼロではありません。大気圧で水が蒸気になると体積は約1600倍になります。蒸発率だけ入れてこの項を落とすと、質量は消えるのに何も押しのけられません。気泡は成長せず、圧力場だけがおかしくなります。

$r$ を選ぶ問題 — Leeモデルの罠#

最も広く使われるソース項がLeeモデルです。

\dot{m}''' = r\,\alpha_l \rho_l \frac{T - T_{sat}}{T_{sat}}

$r$ は緩和係数 [1/s] で、ここから誤解が始まります。 $r$ は物性値ではありません。実験で測れる量でもありません。 $r$ の役目は界面セルの温度を $T_{sat}$ に張り付かせることだけです。温度が $T_{sat}$ に張り付けば、そのセルのエネルギー収支が自然にStefan条件になります。つまり $r$ は、前節の自由境界条件をペナルティで強制する数値的な仕掛けです。

理屈の上では大きいほど良いことになります。問題は陽的に積分するときに出ます。ソース項だけを取り出すと $T - T_{sat}$ に対する線形減衰で、その減衰率が $r L/(T_{sat} c_p)$ です。陽的Eulerの安定条件は次のようになります。

S = r\,\frac{L}{T_{sat}\,c_p}\,\Delta t < 2

水と蒸気では $L/(T_{sat} c_p) \approx 1.43$ です。 $\Delta t = 10^{-5}$ s なら $r$ は $1.4\times10^{5}$ を超えられません。格子を細かくして $\Delta t$ が小さくなれば、許容される $r$ は大きくなります。文献の $r$ が0.1から $10^7$ まで散らばっているのはこのためです。著者ごとに $\Delta t$ が違います。

下の加熱カラムで $r$ を自分で掃引してみてください。

log₁₀ r2.0e+3ΔT_wall10 K

What to watch: the amber shaded area is superheat the model failed to convert into vapor. Drag r down to 10² and the area swells while the interface stalls — the wall keeps pouring in heat and nothing boils. Push past 10⁵ and S crosses 2: the profile starts ringing and the history trace turns into a sawtooth. Green sits in between — roughly 10³ to 10⁵ here — and that window moves whenever you change Δt or the mesh.

$r$ を $10^2$ まで下げると、黄色の網掛け — 蒸気に変換できなかった過熱 — が膨らみ、界面が止まります。壁は熱を注ぎ続けているのに、何も沸きません。 $10^5$ を超えると $S$ が2を跨ぎ、温度プロファイルが振動して右の履歴グラフがのこぎり波になります。緑の区間はその間の2桁分しかなく、 $\Delta t$ や格子を変えればその窓ごと移動します。

物性ベースの代替もあります。Hertz–Knudsen–Schrage関係は気体分子運動論から出発し、界面蒸発率を適応係数(accommodation coefficient) $\sigma$ 一つで表します。Tanasawaモデルはそれを線形化して $\dot{m}'' \propto \sigma (T - T_{sat})$ の形にします。 $\sigma$ は物性なので実測値があります。ただし水の $\sigma$ の報告値も0.01から1まで散らばっており、不確かさは消えるというより場所を移すだけです。

Pythonで両レジームを並べる#

一方は湿度が、もう一方は時間刻みが答えを決めます。同じスクリプトで確認します。

import math
 
M_V, M_A, P_ATM = 18.015, 28.96, 101325.0
RHO_L, RHO_G, D_AB = 997.0, 1.18, 2.5e-5     # 水、湿り空気、水蒸気の拡散係数
L_VAP, CP_L, T_SAT = 2.26e6, 4220.0, 373.15
 
def p_sat(t_c):
    """Antoine式(水、mmHg) -> Pa"""
    return 10 ** (8.07131 - 1730.63 / (233.426 + t_c)) * 133.322
 
def mass_fraction(x):
    """蒸気のモル分率 -> 質量分率"""
    return x * M_V / (x * M_V + (1 - x) * M_A)
 
def spalding_number(t_c, rh):
    """Spalding物質移動数 B_M = (Y_s - Y_inf) / (1 - Y_s)"""
    x_s = p_sat(t_c) / P_ATM
    y_s = mass_fraction(x_s)
    y_inf = mass_fraction(x_s * rh)
    return (y_s - y_inf) / (1 - y_s)
 
def d2_lifetime(d0_mm, t_c, rh):
    """d^2則の寿命 t = d0^2 / K,  K = 8 rho_g D / rho_l * ln(1 + B_M)"""
    b_m = spalding_number(t_c, rh)
    k = 8 * RHO_G * D_AB / RHO_L * math.log(1 + b_m)     # m^2/s
    return (d0_mm * 1e-3) ** 2 / k, k * 1e6              # s, mm^2/s
 
def lee_source(r, alpha_l, temp):
    """Leeモデルの体積蒸発率 [kg/m^3/s]"""
    if temp > T_SAT:
        return r * alpha_l * RHO_L * (temp - T_SAT) / T_SAT
    return 0.0
 
print("拡散律速 - 1 mm 水滴、25 degC")
for rh in (0.0, 0.3, 0.6, 0.9, 0.97):
    life, k = d2_lifetime(1.0, 25.0, rh)
    print(f"  RH {rh*100:4.0f} %   B_M = {spalding_number(25.0, rh):.5f}"
          f"   K = {k:.2e} mm^2/s   寿命 = {life/60:7.2f} min")
 
print("\n熱律速 - Lee緩和係数と陽的安定限界")
dt = 1e-5
for r in (1e2, 1e3, 1e4, 1e5, 1e6):
    s = r * L_VAP / (T_SAT * CP_L) * dt
    rate = lee_source(r, 1.0, T_SAT + 1.0)
    print(f"  r = {r:8.0e}   m''' (1 K 過熱) = {rate:9.2e} kg/m^3/s"
          f"   S = {s:8.3f}  {'発散' if s > 2 else '安定'}")

出力はこうなります。

拡散律速 - 1 mm 水滴、25 degC
  RH    0 %   B_M = 0.02001   K = 4.69e-03 mm^2/s   寿命 =    3.55 min
  RH   30 %   B_M = 0.01406   K = 3.30e-03 mm^2/s   寿命 =    5.04 min
  RH   60 %   B_M = 0.00806   K = 1.90e-03 mm^2/s   寿命 =    8.77 min
  RH   90 %   B_M = 0.00202   K = 4.78e-04 mm^2/s   寿命 =   34.85 min
  RH   97 %   B_M = 0.00061   K = 1.44e-04 mm^2/s   寿命 =  115.98 min
 
熱律速 - Lee緩和係数と陽的安定限界
  r =    1e+02   m''' (1 K 過熱) =  2.67e+02 kg/m^3/s   S =    0.001  安定
  r =    1e+03   m''' (1 K 過熱) =  2.67e+03 kg/m^3/s   S =    0.014  安定
  r =    1e+04   m''' (1 K 過熱) =  2.67e+04 kg/m^3/s   S =    0.144  安定
  r =    1e+05   m''' (1 K 過熱) =  2.67e+05 kg/m^3/s   S =    1.435  安定
  r =    1e+06   m''' (1 K 過熱) =  2.67e+06 kg/m^3/s   S =   14.352  発散

湿度30 %と97 %の間で寿命が23倍開きます。温度はそのままです。拡散律速のケースで計算が実験と合わないときは、たいてい遠方の湿度境界条件が原因です。

実験ノートに写しておくこと

蒸発速度を決めているものを先に見極めてからモデルを選びます。気相側の濃度勾配なら $B_M$ とd²則、界面の熱不均衡ならStefan条件とLee/Tanasawaです。順序を逆にすると、正しい答えが出るのは偶然です。
Leeモデルの $r$ は物性値ではなくペナルティ係数です。論文から値を写すのではなく、 $S = r L \Delta t/(T_{sat} c_p) < 2$ の範囲でできるだけ大きく取り、界面セルの過熱が実際にゼロへ張り付くかを確認します。
蒸発率を入れたなら $\nabla\cdot\mathbf{u} = \dot{m}'''(1/\rho_v - 1/\rho_l)$ も一緒に入れます。これを落とすと質量は消えるのに体積はそのままです。水と蒸気ではその誤差は1600倍です。