Una partícula de aluminio de 1 μm se calentó 4.6 veces más lento que con la fórmula del continuo — interpolación de Sherman en transición de Knudsen

En el régimen de transición, el tiempo de calentamiento es el del continuo multiplicado por 1 + Q_c/Q_fm. En aire a 1500 K, una partícula de 1 μm tarda 4.6 veces más y una de 10 μm, 1.36 veces más.

¿Por qué la partícula de 1 μm se fundió más tarde de lo previsto?#

Una partícula de aluminio a temperatura ambiente cae en aire a 1500 K. ¿Cuál se funde primero, una de 50 μm o una de 1 μm? La pequeña, claro. La masa disminuye como $d^3$ y la conducción solo como $d$ . Según la fórmula del continuo, la partícula de 1 μm alcanza su punto de fusión en 1.95 μs.

Con teoría cinética, la misma partícula tarda 8.95 μs. Son 4.6 veces más. La partícula ya tiene un tamaño comparable al recorrido libre medio del gas (la distancia promedio que recorre una molécula entre choques). Este artículo explica de dónde sale esa diferencia en la etapa de precalentamiento (pre-heating) de un modelo de combustión de aglomerados de aluminio. Se plantean por separado las fórmulas del extremo continuo y del extremo de moléculas libres, y se unen con la interpolación armónica de Sherman. Al final, todo el resultado se reduce a una línea: $t/t_c = 1 + Q_c/Q_{fm}$ .

Qué tan granular ve el gas una sola partícula#

Que un gas parezca continuo depende del tamaño del observador. Esa razón es el número de Knudsen.

\mathrm{Kn} = \frac{2\lambda}{d}, \qquad \lambda = \frac{k_B T}{\sqrt{2}\,\pi \sigma^2 p}

$\lambda$ es el recorrido libre medio, $d$ el diámetro de la partícula, $\sigma$ el diámetro de colisión molecular y $p$ la presión. Para aire a 1500 K y 1 atm con $\sigma = 3.7\times10^{-10}$ m, $\lambda = 0.336$ μm. Es unas cinco veces más largo que a temperatura ambiente. A presión fija, la densidad numérica de moléculas cae al subir la temperatura.

Se suele llamar continuo a $\mathrm{Kn} < 0.01$ y moléculas libres a $\mathrm{Kn} > 10$ . Lo que queda entre ambos es el régimen de transición (transition regime). Una partícula de 50 μm tiene $\mathrm{Kn} = 0.013$ ; una de 1 μm, 0.67. Buena parte del polvo de aluminio micrométrico cae en el régimen de transición.

En la simulación de abajo se pueden variar el diámetro de la partícula y el coeficiente de acomodación.

particle d = 1.00 um accommodation alpha = 1.00

transition | Kn = 0.672 | Q_sh/Q_c = 0.218 | t_sh/t_c = 1.00 | T_c = 300 K | T_sh = 300 K

Al bajar el diámetro, las moléculas de la izquierda vuelan en línea recta distancias mayores que la propia partícula. Las moléculas enfriadas por la partícula (azules) se quedan cerca de ella del lado continuo y salen disparadas del lado de moléculas libres. En la carrera de abajo a la derecha, la razón de tiempos de llegada coincide con el inverso de la altura del punto rojo.

El extremo continuo: el calor llega solo por difusión#

Cuando $\mathrm{Kn} \to 0$ , el gas alrededor de la partícula obedece la ley de Fourier. La conducción estacionaria desde una esfera en gas en reposo da un número de Nusselt (Nu, transferencia convectiva dividida entre conducción pura) igual a 2.

Q_c = 2\pi d\, k\, (T_g - T_p)

$k$ es la conductividad térmica del gas; $T_g$ y $T_p$ son las temperaturas del gas y de la partícula. El calor entra solo en proporción al diámetro. La masa va como $d^3$ , así que el tiempo va como $d^2$ . Reducir el diámetro a la décima parte calienta la partícula cien veces más rápido.

La columna $t_c$ imprime ese cuadrado. 50 μm da $4.883\times10^{-3}$ s y 1 μm da $1.953\times10^{-6}$ s. Un cociente de diámetros de 50, al cuadrado, es 2500, y el cociente de tiempos también es 2500. 10 μm frente a 1 μm, y 3 μm frente a 0.3 μm, son un cociente de diámetros de 10 y un cociente de tiempos de 100. $1.953\times10^{-4}/1.953\times10^{-6}=100$ . $1.758\times10^{-5}/1.758\times10^{-7}=100$ . En el extremo continuo, el diámetro solo cierra el cociente de tiempos. El secado de una gota sigue el mismo $d^2$ . Ese balance reaparece en la evaporación gobernada por difusión.

La fórmula supone que el gas en la superficie está a $T_p$ . En realidad, el gas a menos de un $\lambda$ de la superficie no ha sufrido suficientes choques para llegar ahí. La suposición solo vale cuando esa capa es despreciable frente a la partícula.

El extremo de moléculas libres: contar las moléculas que golpean#

Cuando $\mathrm{Kn} \to \infty$ , las moléculas no chocan entre sí antes de llegar a la partícula. La transferencia de calor es el número de moléculas que golpean la superficie por unidad de tiempo por la energía que entrega cada una.

Q_{fm} = \alpha\, \pi d^2\, \frac{p\, \bar c}{8}\, \frac{\gamma + 1}{\gamma - 1}\, \frac{T_g - T_p}{T_g}, \qquad \bar c = \sqrt{\frac{8 k_B T_g}{\pi m}}

$\alpha$ es el coeficiente de acomodación térmica (thermal accommodation coefficient, qué tanto adopta la molécula reflejada la temperatura de la pared), $\bar c$ la rapidez molecular media, $\gamma$ la razón de calores específicos y $m$ la masa molecular. Ahora $Q_{fm} \propto d^2$ . El calor entra por el área superficial y el tiempo de calentamiento escala solo con $d$ .

El cociente de ambas expresiones muestra la estructura.

\frac{Q_{fm}}{Q_c} = \frac{\alpha\, d\, p\, \bar c\, (\gamma+1)}{16\, k\, T_g\, (\gamma-1)}

Sustituyendo $k \sim \rho c_v \bar c \lambda$ , el cociente escala aproximadamente como $\alpha/\mathrm{Kn}$ . Cuanto menor la partícula, menor la tasa de moléculas libres frente a la del continuo. La menor de las dos es el cuello de botella.

La forma armónica de Sherman: dos resistencias térmicas en serie#

En el régimen de transición ambos mecanismos actúan en secuencia. Hasta aproximadamente un $\lambda$ de la superficie, las moléculas transportan el calor en vuelo libre; más allá, lo hace la difusión. Tratar ambos tramos como resistencias térmicas en serie da la forma armónica de Sherman (1963).

\frac{1}{Q} = \frac{1}{Q_c} + \frac{1}{Q_{fm}}

Es exacta en ambos extremos. Si $\mathrm{Kn} \to 0$ , desaparece $1/Q_{fm}$ ; si $\mathrm{Kn} \to \infty$ , desaparece $1/Q_c$ . La presentación original también recoge el enfoque de esfera límite (boundary sphere) de Fuchs y las formas armónicas de Fuchs–Sutugin y Loyalka. El enfoque de esfera límite envuelve la partícula en una esfera ficticia de radio $\sim\lambda$ , resuelve moléculas libres dentro y continuo fuera, y empalma ambas soluciones en el borde. Las formas armónicas comprimen ese resultado en un solo coeficiente.

Código: seis diámetros y cuatro coeficientes de acomodación#

El código transcribe directamente las tres fórmulas. Integra el balance de energía $m c_p\, dT/dt = Q(T)$ desde 300 K hasta el punto de fusión del aluminio, 933 K. Integrar en temperatura en lugar de tiempo da $dt = m c_p\, dT / Q$ , así que basta un bucle. La radiación y las reacciones superficiales se omiten a propósito.

import math
 
KB = 1.380649e-23          # Boltzmann constant [J/K]
AMU = 1.66054e-27          # atomic mass unit [kg]
 
# gas: air at 1500 K, 1 atm
TG, P = 1500.0, 101325.0
K_GAS = 0.095              # thermal conductivity [W/m/K]
GAMMA = 1.4
M_GAS = 28.97 * AMU
SIGMA = 3.7e-10            # collision diameter [m]
 
# particle: aluminum
RHO_P, CP_P = 2700.0, 1100.0
T0, T_MELT = 300.0, 933.0
 
def mean_free_path(T, p):
    return KB * T / (math.sqrt(2.0) * math.pi * SIGMA**2 * p)
 
def q_continuum(d, Tp):
    return 2.0 * math.pi * d * K_GAS * (TG - Tp)
 
def q_free_molecular(d, Tp, alpha):
    cbar = math.sqrt(8.0 * KB * TG / (math.pi * M_GAS))
    return alpha * math.pi * d**2 * P * cbar / 8.0 \
        * (GAMMA + 1.0) / (GAMMA - 1.0) * (TG - Tp) / TG
 
def q_sherman(d, Tp, alpha):
    qc, qf = q_continuum(d, Tp), q_free_molecular(d, Tp, alpha)
    return 1.0 / (1.0 / qc + 1.0 / qf)
 
def heat_up_time(d, law, n=4000):
    # integrate m*cp*dT/dt = Q(T) from T0 to T_MELT in T (dt = m cp dT / Q)
    m = RHO_P * math.pi * d**3 / 6.0
    dT = (T_MELT - T0) / n
    t = 0.0
    for i in range(n):
        T = T0 + (i + 0.5) * dT
        t += m * CP_P * dT / law(d, T)
    return t
 
lam = mean_free_path(TG, P)
print(f"lambda at {TG:.0f} K, 1 atm = {lam*1e6:.3f} um")
print()
print("alpha = 1.0")
print(f"{'d [um]':>8} {'Kn=2l/d':>8} {'Qfm/Qc':>8} {'Qsh/Qc':>8} {'t_c [s]':>10} {'t_sh [s]':>10} {'ratio':>6}")
for d_um in [50, 10, 3, 1, 0.3, 0.1]:
    d = d_um * 1e-6
    kn = 2 * lam / d
    r_fm = q_free_molecular(d, T0, 1.0) / q_continuum(d, T0)
    r_sh = q_sherman(d, T0, 1.0) / q_continuum(d, T0)
    tc = heat_up_time(d, q_continuum)
    ts = heat_up_time(d, lambda dd, T: q_sherman(dd, T, 1.0))
    print(f"{d_um:8.1f} {kn:8.3f} {r_fm:8.2f} {r_sh:8.3f} {tc:10.3e} {ts:10.3e} {ts/tc:6.2f}")
print()
print("d = 1 um, varying accommodation coefficient")
print(f"{'alpha':>6} {'Qsh/Qc':>8} {'t_sh/t_c':>9}")
d = 1e-6
tc = heat_up_time(d, q_continuum)
for a in [1.0, 0.5, 0.1, 0.02]:
    ts = heat_up_time(d, lambda dd, T: q_sherman(dd, T, a))
    print(f"{a:6.2f} {q_sherman(d, T0, a)/q_continuum(d, T0):8.3f} {ts/tc:9.2f}")

lambda at 1500 K, 1 atm = 0.336 um
 
alpha = 1.0
  d [um]  Kn=2l/d   Qfm/Qc   Qsh/Qc    t_c [s]   t_sh [s]  ratio
    50.0    0.013    13.96    0.933  4.883e-03  5.233e-03   1.07
    10.0    0.067     2.79    0.736  1.953e-04  2.653e-04   1.36
     3.0    0.224     0.84    0.456  1.758e-05  3.857e-05   2.19
     1.0    0.672     0.28    0.218  1.953e-06  8.949e-06   4.58
     0.3    2.240     0.08    0.077  1.758e-07  2.275e-06  12.94
     0.1    6.721     0.03    0.027  1.953e-08  7.191e-07  36.82
 
d = 1 um, varying accommodation coefficient
 alpha   Qsh/Qc  t_sh/t_c
  1.00    0.218      4.58
  0.50    0.122      8.16
  0.10    0.027     36.82
  0.02    0.006    180.09

El 4.58 sale de una fórmula de una línea#

Ambas tasas son proporcionales a $(T_g - T_p)$ , así que $Q_{fm}/Q_c$ no depende de la temperatura de la partícula. La tasa de Sherman es la del continuo multiplicada por una constante $s = 1/(1 + Q_c/Q_{fm})$ . La razón de tiempos de calentamiento es su inverso.

\frac{t_{sh}}{t_c} = 1 + \frac{Q_c}{Q_{fm}}

A 1 μm, $Q_{fm}/Q_c = 0.28$ , así que $1 + 1/0.28 \approx 4.6$ , igual que el 4.58 de la tabla. A 0.1 μm es $1 + 1/0.03$ , unas 37 veces. El bucle de integración solo confirma esta aritmética.

El lado armónico es la misma división. $Q/Q_c=r/(1+r)$ , con $r=Q_{fm}/Q_c$ . A 1 μm, $0.28/1.28\approx0.219$ frente al 0.218 impreso. A 10 μm, $2.79/3.79=0.736$ , igual que la tabla. A 3 μm, $0.84/1.84\approx0.457$ frente a 0.456. A 50 μm, $13.96/14.96=0.933$ . La salida cabe en un paso de redondeo de esta expresión.

Si $Q_{fm}=Q_c$ , el cociente de tiempos es exactamente 2. El calor baja a la mitad. El tiempo se duplica. Los 3.6 μm y $\mathrm{Kn}\approx0.19$ del texto son ese punto. Ya es un factor dos, en un diámetro que todavía se llamaría continuo.

$\alpha$ multiplica $r$ de forma directa. Con $\alpha=0.5$ , $r=0.14$ y $1+1/0.14\approx8.14$ , junto al 8.16 de la tabla. Con $\alpha=0.1$ , $r\approx0.028$ y $1+1/0.028\approx36.7$ , junto a 36.82. Con $\alpha=0.02$ , $1+1/0.0056\approx180$ , junto a 180.09. El mismo $r/(1+r)$ da 0.123, 0.027 y 0.006 frente a 0.122, 0.027 y 0.006. No es un cálculo nuevo. Es el bloque de salida dividido otra vez.

La fila que merece atención es la de 10 μm. Con $\mathrm{Kn} = 0.067$ , un rango que casi todos llamarían continuo, el calor entregado cae un 26%. El diámetro de cruce donde $Q_{fm} = Q_c$ es de unos 3.6 μm, es decir $\mathrm{Kn} \approx 0.19$ . El cruce queda cinco veces por debajo de $\mathrm{Kn} = 1$ . Se debe a $(\gamma+1)/(\gamma-1) = 6$ , al factor $1/16$ y al coeficiente que liga $k$ con $\lambda$ . La regla práctica "Kn menor que 0.1 significa continuo" es demasiado laxa para transferencia de calor.

El coeficiente de acomodación pesa tanto como el diámetro#

La parte inferior de la tabla varía solo $\alpha$ a 1 μm. Con $\alpha = 0.1$ , el calentamiento tarda 37 veces más; con 0.02, 180 veces. $Q_{fm}$ entra solo a través del producto $\alpha d$ . Dividir $\alpha$ entre diez equivale a dividir el diámetro entre diez.

El problema es que $\alpha$ se conoce mal. Los valores reportados entre superficies metálicas y gases calientes varían mucho según la capa de óxido, la razón de temperaturas y la especie gaseosa. Filippov & Rosner (2000) revisaron la propia transferencia de energía a razones de temperatura grandes. Algunos estudios fijan $\alpha$ muy por debajo de 1 para ajustar retardos de ignición de nanoaluminio. En modelos de partículas en transición, $\alpha$ se convierte fácilmente en una perilla de ajuste. Que cambie los resultados varias veces debe constar en el informe.

Dónde se tambalea la forma armónica#

La forma de Sherman acierta en ambos extremos, pero la curva intermedia es una aproximación. Para $\mathrm{Kn} \sim 0.1$ –1, la solución de esfera límite de Fuchs, la fórmula de Fuchs–Sutugin y la corrección de Loyalka difieren entre sí en decenas de por ciento. El 26% a 10 μm debe leerse como un número dentro de esa dispersión de modelos.

La definición de $\lambda$ añade una incertidumbre del mismo orden. Fórmula de esfera rígida o definición basada en viscosidad, $\lambda/r$ o $2\lambda/d$ : el Kn de la misma partícula cambia fácilmente decenas de por ciento. Al tomar una correlación de otro artículo, conviene traer también su definición de $\lambda$ . De lo contrario, la esfera límite queda en el lugar equivocado.

El cuello de botella sigue en la combustión. La difusión del oxidante hacia la partícula también pide una corrección de Knudsen del lado de la masa. Al aglomerarse, el diámetro efectivo crece y Kn baja otra vez. Volver a Nu = 2 en ese punto exige recalcular Kn con ese diámetro efectivo.

Antes de escribir una línea de transferencia de calor para partículas lagrangianas#

En los códigos Euler–Lagrange, el calor de las partículas suele empezar con una sola línea: $\mathrm{Nu} = 2 + 0.6\,\mathrm{Re}^{1/2}\mathrm{Pr}^{1/3}$ . Esa línea es el extremo continuo. Una cola en 3 μm lleva el cociente 2.19 de la tabla. En 1 μm el cociente es 4.58. La idea de dos o tres veces solo cabe cerca de 3 μm. En 1 μm el desvío ya pasa de cuatro. Si el retardo de ignición sale corto en los diámetros pequeños, esta línea es el primer sitio donde mirar.

La aritmética es barata. $\lambda$ por celda, $Q_{fm}$ por partícula, media armónica. Sin anotar de dónde salen $\alpha$ y la definición de $\lambda$ , el 4.58 no vuelve. Ese factor es casi esos dos números. Subir el flujo de calor tampoco mantiene la pendiente. Igual que en la curva de ebullición donde el flujo de calor se quiebra, el mecanismo cambia y el calor entregado se dobla. Aquí el quiebre es $Q_{fm}$ .