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

Si la difusión limita, la humedad lejana fija la vida de la gota. En el régimen térmico, r deja de servir en cuanto S cruza 2.

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.

Sin separarlas, el término fuente de una de las dos sale mal. Las ecuaciones cubren la evaporación limitada por difusión y la limitada por calor, y ambas pasan a Navier–Stokes como términos fuente. Al final queda cómo elegir el coeficiente de relajación $r$ de Lee, y por qué los valores publicados se reparten sobre siete órdenes de magnitud.

Al subir la humedad relativa del 30% al 95% la recta de vida se acuesta; al subir la temperatura, la misma gota de 1 mm se reduce a unos segundos.

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.

La pendiente de la derecha es $B_M$ , y el radio del disco de la izquierda usa el mismo $d(t)$ . Cuando la humedad se acerca a $Y_s$ , la vida se estira hacia una hora.

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í. El tramo de la curva de ebullición en el que una pared más caliente transporta menos calor cae en este régimen.

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$ .

$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.

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.

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

El $S$ de la tabla se vuelve a barrer con $r$ en la columna calentada.

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 ámbar — sobrecalentamiento que no pasó a vapor — se hincha y la interfaz se detiene. La pared sigue metiendo calor y nada hierve. Al pasar de $10^5$ , $S$ cruza 2, la temperatura oscila y la traza de la derecha se vuelve un diente de sierra. El verde cubre unas dos décadas, y esa ventana se desplaza entera si cambian $\Delta t$ o la malla.

La fila seca de la tabla, 1 mm a 25 °C y HR 0 %, imprime $B_M=0.02001$ , $K=4.69\times10^{-3}$ mm²/s y una vida de 3.55 min. Con $d_0^2=1$ mm² la vida es $1/K$ segundos: $1/0.00469\approx213$ s, y 213/60 da 3.55 min. HR 30 % son 5.04 min con $K=3.30\times10^{-3}$ . HR 97 % son 115.98 min con $K=1.44\times10^{-4}$ . El cociente de vidas 115.98/5.04 vale unos 23.0, y el cociente inverso de $K$ , $3.30\times10^{-3}/1.44\times10^{-4}$ , unos 22.9. La temperatura sigue en 25 °C. Entre el 30 % y el 97 % de humedad la vida se abre por un factor 23. Cuando un caso limitado por difusión no reproduce el experimento, suele fallar la humedad del campo lejano.

La columna $S$ es el otro régimen. Con $r=10^5$ , $S=1.435$ ; con $r=10^6$ , $S=14.352$ . El límite 2 queda entre ambos. Con $L/(T_{sat} c_p)\approx1.43$ y $\Delta t=10^{-5}$ , $2/(1.43\times10^{-5})\approx1.4\times10^{5}$ . El $r=10^5$ del código está dentro de ese tope y $10^6$ lo supera por diez. A 1 K de sobrecalentamiento, $\dot{m}'''$ escala con $r$ desde $2.67\times10^{2}$ hasta $2.67\times10^{6}$ . El $\rho_g=1.18$ de la ley d² es aire húmedo. El factor 1600 de la restricción de divergencia es el salto de volumen frente al vapor en ebullición. Mezclar esas dos densidades convierte la vida de la ropa y la expansión de la tetera en el mismo número.

Fallar la humedad de contorno, o dejar que S cruce 2#

Si una vida limitada por difusión no cuadra con el laboratorio, la humedad lejana va antes que el modelo. En la tabla, HR 30 % y 97 % comparten temperatura y aun así duran 5.04 min y 115.98 min. Copiar un $r$ de un artículo al régimen térmico hace que otro $\Delta t$ empuje $S$ más allá de 2, o que el sobrecalentamiento no se apague.

Meter una tasa de evaporación y dejar en cero el lado derecho de $\nabla\cdot\mathbf{u}$ quita masa y deja el volumen quieto. Para agua y vapor ese desajuste es de unas 1600 veces. Una gota que vuela mientras se seca hay que leerla junto con la corrección de arrastre que añade la vaporización. Poner el $r$ de Lee en el sitio de $B_M$ hace que la ropa espere a los 100 °C.