Por qué la ropa se seca a 5 °C — evaporación limitada por difusión frente a ebullición térmica

Por qué la evaporación sin ebullición y con ebullición piden términos fuente distintos

La ropa tendida en el balcón se seca. El aire está a 5 °C. El agua hierve a 100 °C. Faltan noventa y cinco grados y el agua desaparece igual.

Esa misma agua dentro de una tetera espera hasta los 100 °C. Ambos casos son cambio de fase: líquido que pasa a vapor. Pero lo que fija la velocidad es distinto en cada uno. En CFD ocurre lo mismo. Usar un único término fuente para las dos situaciones garantiza que una de ellas salga mal.

Este artículo las separa. Plantea las ecuaciones de la evaporación limitada por difusión y de la limitada por calor, y lleva ambas a las ecuaciones de Navier–Stokes como términos fuente. Al final queda elegir a mano el coeficiente de relajación $r$ del modelo de Lee, y ver por qué los valores publicados de $r$ se reparten sobre siete órdenes de magnitud.

Las dos evaporaciones tienen cuellos de botella distintos#

La evaporación necesita dos cosas: energía suficiente para que una molécula abandone la interfaz, y un camino que se lleve esa molécula una vez fuera. Cuál de las dos escasea define el régimen.

Limitada por difusión. El gas justo encima de la interfaz ya está saturado. Para seguir evaporando hay que transportar ese vapor lejos. El cuello de botella es la difusión del lado del gas. La temperatura del líquido apenas cambia. Los charcos, la ropa mojada y los aerosoles de combustible a temperatura ambiente viven aquí.

Impulsada térmicamente. El líquido cercano a la interfaz ya superó la temperatura de saturación. Al vapor le sobra sitio para irse. El cuello de botella es el calor que debe pagar el calor latente (la energía oculta que consume el cambio de fase). La ebullición, la cavitación y la condensación en pared viven aquí.

Incluso los números adimensionales que miden cada régimen son distintos. La evaporación por difusión se escala con el número de transferencia de masa de Spalding $B_M$ (la caída de concentración de vapor entre interfaz y campo lejano); la térmica, con el número de Jakob $\mathrm{Ja} = c_p \Delta T / L$ (calor sensible frente a latente). El primero se empareja con el número de Schmidt $\mathrm{Sc} = \nu/D$ ; el segundo, con el de Prandtl.

Limitada por difusión — el gradiente de concentración fija la velocidad#

Se plantea el flujo másico de vapor que sale de una gota esférica. La ley de Fick (difusión proporcional al gradiente de concentración) no basta por sí sola, porque el vapor que escapa empuja hacia afuera a toda la mezcla gaseosa. Esa componente convectiva se llama flujo de Stefan.

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

$\dot{m}''$ es la tasa de evaporación por unidad de área, $\rho_g$ la densidad del gas, $D$ la difusividad del vapor, $Y_v$ la fracción másica de vapor y el subíndice $s$ indica la interfaz. El segundo término del lado derecho es la parte que el flujo de Stefan vuelve a arrastrar.

Al integrar radialmente bajo hipótesis cuasiestacionaria aparece un logaritmo.

\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$ es el radio de la gota, $Y_s$ la fracción másica saturada en la interfaz y $Y_\infty$ el valor lejano. $B_M$ es el número de transferencia de masa de Spalding: toda la caída de concentración que impulsa la evaporación comprimida en un solo valor.

Si se añade la conservación de masa de la gota, $\dot{m} = -\frac{d}{dt}(\rho_l \pi d^3/6)$ , resulta que el cuadrado del diámetro decrece linealmente.

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

Esta es la ley d². Lo que se dibuja como recta es $d^2$ , no $d$ . La vida útil vale $t_{life} = d_0^2 / K$ .

La clave está en que $Y_s$ depende solo de la temperatura. El agua a 25 °C tiene una presión de saturación de 3,17 kPa, un 3 % de la atmosférica, lo que da $Y_s \approx 0{,}019$ . No es cero. Por eso las cosas se secan. A 5 °C baja hasta $Y_s \approx 0{,}0054$ , y sigue sin ser cero. El punto de ebullición nunca entra en esta historia.

Conviene manipular los parámetros directamente en la simulación siguiente.

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.

Al llevar el control de humedad relativa del 30 % al 95 %, $Y_\infty$ sube hacia $Y_s$ , $B_M$ se derrumba y la recta de la derecha se acuesta. Una gota de 1 mm pasa de durar 5 minutos a casi una hora. Si en cambio se sube la temperatura, $Y_s$ se dispara exponencialmente y esa misma gota desaparece en segundos. El radio del disco de la izquierda se dibuja con el mismo $d(t)$ que traza la gráfica.

Impulsada térmicamente — el calor que llega paga el latente#

Ahora la interfaz está clavada en la temperatura de saturación. Lo que fija la tasa no es la concentración, sino el desequilibrio de calor que llega a la interfaz.

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

$L$ es el calor latente, $k_l$ y $k_v$ las conductividades del líquido y del vapor, y $n$ la normal a la interfaz. Todo el calor que entra por el lado líquido y no sale por el lado vapor se va íntegro al cambio de fase.

Esto se conoce como condición de Stefan, y resulta incómoda por un motivo concreto: sobre la interfaz cuelgan dos condiciones a la vez. La temperatura queda fijada en $T = T_{sat}$ (Dirichlet) y, simultáneamente, el salto de flujo determina la velocidad de la interfaz. Dos condiciones, y la posición del contorno es ella misma una incógnita. Es un problema de frontera libre.

Resolverlo con honestidad mediante interfaz aguda exige rastrear esa interfaz y satisfacer la ecuación anterior en cada paso. La mayoría de los solvers comerciales y de código abierto evita ese camino. En su lugar dejan que la interfaz se difumine sobre una o dos celdas e introducen un término fuente volumétrico diseñado para que la condición anterior se cumpla por sí sola.

Dónde aparece realmente el término fuente#

Basta definir una tasa volumétrica de evaporación $\dot{m}'''$ en kg/m³/s para que aparezca en cuatro sitios.

Entra en las ecuaciones de continuidad por fase como un par de signo opuesto.

\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$ es la fracción volumétrica. Al sumar ambas los fuentes se cancelan, de modo que la masa de la mezcla se conserva.

En la ecuación de energía aparece como sumidero por el valor del calor latente. Si el solver además transporta una fracción másica de vapor (la vía limitada por difusión), el mismo valor aparece allí como fuente.

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

El cuarto sitio es el que suele olvidarse. El cambio de fase modifica la restricción de divergencia del campo de velocidad.

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

La ecuación de Poisson de presión de un solver incompresible suele tener cero a la derecha. Con cambio de fase deja de tenerlo. A presión atmosférica, el agua que se convierte en vapor se expande unas 1600 veces. Si se mete la tasa de evaporación pero se omite este término, la masa desaparece sin que nada sea desplazado: la burbuja no crece y el campo de presión se vuelve extraño.

Elegir $r$ — la trampa del modelo de Lee#

El término fuente más extendido es el modelo de Lee.

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

$r$ es el coeficiente de relajación, en 1/s, y ahí empieza el malentendido. $r$ no es una propiedad del material. Tampoco es algo medible. Su única función es sujetar la temperatura de la celda interfacial en $T_{sat}$ . Una vez clavada ahí, el balance de energía de esa celda se convierte en la condición de Stefan. Dicho de otro modo, $r$ es un artificio numérico que impone por penalización la condición de frontera libre de la sección anterior.

En teoría, entonces, cuanto mayor mejor. El problema surge al integrar de forma explícita. Aislando el término fuente queda un decaimiento lineal en $T - T_{sat}$ con tasa $r L/(T_{sat} c_p)$ . La condición de estabilidad del Euler explícito es por tanto

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

Para agua y vapor, $L/(T_{sat} c_p) \approx 1{,}43$ . Con $\Delta t = 10^{-5}$ s, $r$ no puede pasar de $1{,}4\times10^{5}$ . Si se refina la malla y $\Delta t$ se reduce, el $r$ admisible crece. Ese es el motivo de que los valores publicados de $r$ vayan de 0,1 hasta $10^7$ : cada autor trabajó con un $\Delta t$ distinto.

Conviene barrer $r$ personalmente en la columna calentada siguiente.

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.

Al bajar $r$ hasta $10^2$ , la zona sombreada en ámbar — el sobrecalentamiento que el modelo no logró convertir en vapor — se hincha mientras la interfaz se detiene. La pared sigue inyectando calor y nada hierve. Al pasar de $10^5$ , $S$ cruza 2: el perfil de temperatura oscila y la traza histórica de la derecha se vuelve un diente de sierra. El verde ocupa apenas dos décadas intermedias, y esa ventana se desplaza entera en cuanto se cambia $\Delta t$ o la malla.

Existen alternativas basadas en propiedades. La relación de Hertz–Knudsen–Schrage parte de la teoría cinética y expresa la tasa interfacial mediante un único coeficiente de acomodación $\sigma$ . El modelo de Tanasawa la lineariza hasta la forma $\dot{m}'' \propto \sigma (T - T_{sat})$ . Como $\sigma$ sí es una propiedad, tiene valores medidos, aunque los reportados para el agua se dispersan entre 0,01 y 1: la incertidumbre más que desaparecer cambia de sitio.

Los dos regímenes en un mismo script de Python#

En un lado la humedad decide la respuesta; en el otro, el paso temporal. Se comprueban en el mismo script.

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     # agua, aire húmedo, difusividad del vapor
L_VAP, CP_L, T_SAT = 2.26e6, 4220.0, 373.15
 
def p_sat(t_c):
    """Ecuación de Antoine (agua, mmHg) -> Pa"""
    return 10 ** (8.07131 - 1730.63 / (233.426 + t_c)) * 133.322
 
def mass_fraction(x):
    """fracción molar de vapor -> fracción másica"""
    return x * M_V / (x * M_V + (1 - x) * M_A)
 
def spalding_number(t_c, rh):
    """número de 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):
    """vida útil de la ley 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):
    """tasa volumétrica de evaporación del modelo de Lee [kg/m^3/s]"""
    if temp > T_SAT:
        return r * alpha_l * RHO_L * (temp - T_SAT) / T_SAT
    return 0.0
 
print("limitada por difusion - gota de agua de 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"  HR {rh*100:4.0f} %   B_M = {spalding_number(25.0, rh):.5f}"
          f"   K = {k:.2e} mm^2/s   vida = {life/60:7.2f} min")
 
print("\nimpulsada termicamente - coeficiente de Lee y limite explicito")
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 de sobrecalentamiento) = {rate:9.2e} kg/m^3/s"
          f"   S = {s:8.3f}  {'diverge' if s > 2 else 'estable'}")

La salida es esta.

limitada por difusion - gota de agua de 1 mm, 25 degC
  HR    0 %   B_M = 0.02001   K = 4.69e-03 mm^2/s   vida =    3.55 min
  HR   30 %   B_M = 0.01406   K = 3.30e-03 mm^2/s   vida =    5.04 min
  HR   60 %   B_M = 0.00806   K = 1.90e-03 mm^2/s   vida =    8.77 min
  HR   90 %   B_M = 0.00202   K = 4.78e-04 mm^2/s   vida =   34.85 min
  HR   97 %   B_M = 0.00061   K = 1.44e-04 mm^2/s   vida =  115.98 min
 
impulsada termicamente - coeficiente de Lee y limite explicito
  r =    1e+02   m''' (1 K de sobrecalentamiento) =  2.67e+02 kg/m^3/s   S =    0.001  estable
  r =    1e+03   m''' (1 K de sobrecalentamiento) =  2.67e+03 kg/m^3/s   S =    0.014  estable
  r =    1e+04   m''' (1 K de sobrecalentamiento) =  2.67e+04 kg/m^3/s   S =    0.144  estable
  r =    1e+05   m''' (1 K de sobrecalentamiento) =  2.67e+05 kg/m^3/s   S =    1.435  estable
  r =    1e+06   m''' (1 K de sobrecalentamiento) =  2.67e+06 kg/m^3/s   S =   14.352  diverge

Entre el 30 % y el 97 % de humedad la vida útil se abre por un factor 23, con la temperatura sin tocar. Cuando un caso limitado por difusión se niega a reproducir el experimento, el culpable suele ser la condición de contorno de humedad en el campo lejano.

Para la libreta de laboratorio#

Primero conviene determinar qué limita la evaporación y después elegir el modelo. Si es el gradiente de concentración del lado gas, tocan $B_M$ y la ley d²; si es el desequilibrio térmico interfacial, tocan la condición de Stefan y Lee o Tanasawa. Invertir ese orden hace que un resultado correcto sea pura suerte.
El $r$ del modelo de Lee es un coeficiente de penalización, no una propiedad. No conviene copiar un valor de un artículo: se toma tan grande como permita $S = r L \Delta t/(T_{sat} c_p) < 2$ y después se verifica que el sobrecalentamiento de la celda interfacial cae realmente a cero.
Si se añadió una tasa de evaporación, hay que añadir también $\nabla\cdot\mathbf{u} = \dot{m}'''(1/\rho_v - 1/\rho_l)$ . Omitirlo hace que la masa desaparezca mientras el volumen se queda donde estaba. Para agua y vapor ese error es un factor 1600.