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

拡散律速では遠方の湿度が寿命を決める。熱律速ではSが2を超えたrで温度履歴がのこぎり波になる。

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

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

二つを分けないとソース項が合いません。拡散が律速の蒸発と熱が律速の蒸発を別の式で立て、それをNavier–Stokes方程式のソース項へ移します。最後はLeeモデルの緩和係数 $r$ の選び方と、文献の $r$ が7桁に散らばる理由です。

RHを30%から95%へ上げると寿命の直線が寝て、温度を上げると同じ1 mmの水滴が数秒まで縮みます。

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.

右の傾きが $B_M$ で、左の円の半径はそのグラフと同じ $d(t)$ です。湿度が $Y_s$ に近づくと寿命は1時間側へ伸びます。

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

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

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

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

二つのレジームを測る無次元数も異なります。拡散律速は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$ まで下がりますが、それでもゼロではありません。沸点はこの話に登場しません。

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

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

\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$ が違います。

物性ベースの代替もあります。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  発散

表の $S$ を、加熱カラムで $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$ や格子が変わると窓ごと動きます。

表の 1 mm、25℃、RH 0 % は $B_M=0.02001$ 、 $K=4.69\times10^{-3}$ mm²/s、寿命 3.55 分です。 $d_0^2=1$ mm² なので寿命は $1/K$ 秒です。 $1/0.00469\approx213$ 秒、213/60=3.55 分と一致します。RH 30 % は 5.04 分、 $K=3.30\times10^{-3}$ です。RH 97 % は 115.98 分、 $K=1.44\times10^{-4}$ です。寿命の比は 115.98/5.04≈23.0 で、 $K$ の逆比 $3.30\times10^{-3}/1.44\times10^{-4}\approx22.9$ と同じです。温度は 25℃ のままです。湿度 30 % と 97 % の間で寿命は 23 倍開きます。拡散律速で実験と合わないときは、遠方の湿度境界が原因であることが多いです。

$S$ の列は別のレジームです。 $r=10^5$ で $S=1.435$ 、 $r=10^6$ で $S=14.352$ です。限界の 2 はその間にあります。 $L/(T_{sat} c_p)\approx1.43$ 、 $\Delta t=10^{-5}$ なら $2/(1.43\times10^{-5})\approx1.4\times10^{5}$ です。コードの $r=10^5$ は限界の内側で、 $10^6$ はその 10 倍です。1 K 過熱の $\dot{m}'''$ は $r$ に比例して $2.67\times10^{2}$ から $2.67\times10^{6}$ まで進みます。d²則の $\rho_g=1.18$ は湿り空気です。発散拘束の 1600 倍は沸騰蒸気との体積比なので、記号を混ぜると洗濯物の寿命とやかんで膨張が同じ数字になります。

湿度境界を外すか、 $S$ が2を超えるか#

拡散律速で寿命が実験とずれるなら、モデルより先に遠方の湿度を疑います。表の RH 30 % と 97 % は温度が同じなのに、寿命は 5.04 分と 115.98 分です。熱律速で $r$ を論文から写すと、 $\Delta t$ が変わった瞬間に $S$ が2を超えるか、過熱が消えずに残ります。

蒸発率を入れて $\nabla\cdot\mathbf{u}$ の右辺をゼロのままにすると、質量は減るのに体積はそのままです。水と蒸気ではそのずれは約1600倍です。飛びながら乾く液滴は、抗力に気化補正が付く側と寿命を一緒に見る必要があります。 $B_M$ の席に Lee の $r$ を置くと、洗濯物は100℃まで待ちます。