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

扩散控制时远场湿度决定寿命。热控制时一旦S超过2，Lee系数会把温度曲线打成锯齿。

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

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

不分开的话源项会对不上。扩散控制的蒸发和热控制的蒸发各写一套方程，再搬进 Navier–Stokes 作为源项。最后才是 Lee 模型松弛系数 $r$ 怎么取，以及文献里的 $r$ 为什么散在七个数量级上。

把 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$ 时，寿命会拉到接近一小时。

两种蒸发的瓶颈不同

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

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

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

衡量两个区制的无量纲数也不同。扩散控制用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$ ，依然不是零。沸点在这个故事里根本没有出场。

热控制 — 送达的热量偿付潜热

现在界面被钉在饱和温度上。决定蒸发率的不是浓度，而是抵达界面的热量的不平衡。

\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$ 不同。

基于物性的替代方案也有。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，温度开始振荡，右侧历史变成锯齿。绿色窗口只有中间两个数量级， $\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$ 是它的十倍。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}$ 的右端留成 0，质量在减少，体积却不动。水–蒸汽的这个错位大约是 1600 倍。边飞边干的液滴，要和阻力上附带气化修正的那一侧一起看寿命。把 Lee 的 $r$ 放进 $B_M$ 的位子，衣服会一直等到 100℃。