衣服为什么在5℃也会晾干 — 扩散控制蒸发与热控制沸腾

晾在阳台上的衣服会干。气温是5℃。水的沸点是100℃。差了95度，水却照样消失。

同样的水在水壶里要等到100℃才消失。两种情况都是液体变成气体的相变。但决定速率的东西不同。CFD里也一样。给两种情况用同一个源项，其中一个必然算错。

本文把它们分开。分别为扩散控制的蒸发和热控制的蒸发建立方程，再把两者搬进Navier–Stokes方程作为源项。最后可以亲手挑选Lee模型的松弛系数 $r$ ，并看清文献中的 $r$ 为何散布在七个数量级上。

两种蒸发的瓶颈不同

蒸发需要两样东西：分子离开界面所需的能量，以及离开后被带走的通道。哪一样不够，就决定了属于哪个区制。

扩散控制。界面正上方的气体已经饱和。要继续蒸发，就必须把蒸汽输运到远处。瓶颈在气相侧的扩散。液体温度几乎不变。水洼、湿衣服、常温燃油喷雾都属于这一类。

热控制。界面附近的液体已经越过饱和温度。蒸汽的去路绰绰有余。瓶颈在于供应潜热（相变中隐性消耗的热量）的热流。沸腾、空化、壁面凝结属于这一类。

衡量两个区制的无量纲数也不同。扩散控制用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$ 表示界面。右端第二项就是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)$ ，直径的平方就线性减小。

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

这就是d²定律。画成直线的是 $d^2$ ，不是 $d$ 。寿命为 $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分钟延长到接近一小时。反过来提高温度， $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老实求解，就得追踪界面并在每一步满足上式。多数商业与开源求解器绕开了这条路。它们让界面模糊到一两个网格厚，再插入一个体积源项，使上述条件自行得到满足。

源项究竟落在哪几处

定义一个体积蒸发率 $\dot{m}'''$ [kg/m³/s]，它会出现在控制方程的四个位置。

它以符号相反的一对进入相连续方程。

\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

最容易漏掉的是第四处。相变会改变速度场的散度约束。

\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，温度剖面开始振荡，右侧历史曲线变成锯齿。绿色区间只有中间的两个数量级，而且一旦改变 $\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倍。