A 1 μm Aluminum Particle Heated 4.6× Slower Than the Continuum Formula — Sherman Interpolation in the Knudsen Transition Regime

Why did the 1 μm particle melt later than predicted?#

Drop a room-temperature aluminum particle into 1500 K air. Which melts first, a 50 μm particle or a 1 μm one? The small one, obviously. Mass shrinks like $d^3$ while conduction shrinks only like $d$. The continuum formula says the 1 μm particle reaches its melting point in 1.95 μs.

Run the same particle through kinetic theory and you get 8.95 μs. That is 4.6 times longer. The particle has become comparable to the gas mean free path (the average distance a molecule travels between collisions). This post traces where that gap comes from in the pre-heating stage of an aluminum agglomerate combustion model. We build the continuum and free-molecular formulas separately, join them with Sherman's harmonic interpolation, and show that the whole result collapses to one line: $t/t_c = 1 + Q_c/Q_{fm}$.

How grainy the gas looks to a single particle#

Whether a gas looks like a continuum depends on the size of the observer. That ratio is the Knudsen number.

$\lambda$ is the mean free path, $d$ the particle diameter, $\sigma$ the molecular collision diameter, and $p$ the pressure. For air at 1500 K and 1 atm with $\sigma = 3.7\times10^{-10}$ m, $\lambda = 0.336$ μm. That is about five times longer than at room temperature, because molecular number density drops as temperature rises at fixed pressure.

The usual labels are continuum for $\mathrm{Kn} < 0.01$ and free molecular for $\mathrm{Kn} > 10$. Everything between is the transition regime. A 50 μm particle sits at $\mathrm{Kn} = 0.013$; a 1 μm particle at 0.67. A good share of micron-sized aluminum powder lives in the transition regime.

Try changing the particle diameter and the accommodation coefficient in the simulation below.

As you lower the diameter, molecules on the left fly straight paths longer than the particle itself. Molecules cooled by the particle (blue) linger near it on the continuum side and shoot away on the free-molecular side. In the race at bottom right, note that the ratio of arrival times equals the inverse of the red dot's height.

The continuum end: heat arrives only by diffusion#

As $\mathrm{Kn} \to 0$, the gas around the particle obeys Fourier's law. Steady conduction from a sphere in quiescent gas gives a Nusselt number (Nu, convective heat transfer divided by pure conduction) of 2.

$k$ is the gas thermal conductivity; $T_g$ and $T_p$ are the gas and particle temperatures. Heat comes in proportional to diameter only. Mass goes like $d^3$, so time goes like $d^2$. Cut the diameter by ten and the particle heats a hundred times faster.

The $t_c$ column prints that square. 50 μm is $4.883\times10^{-3}$ s and 1 μm is $1.953\times10^{-6}$ s. A diameter ratio of 50 squares to 2500, and the time ratio is 2500 too. 10 μm versus 1 μm, and 3 μm versus 0.3 μm, are a diameter ratio of 10 and a time ratio of 100. $1.953\times10^{-4}/1.953\times10^{-6}=100$. $1.758\times10^{-5}/1.758\times10^{-7}=100$. On the continuum end, diameter alone closes the heat-up ratio. A droplet drying under diffusion follows the same $d^2$. That balance shows up again in evaporation controlled by diffusion.

This formula assumes the gas at the particle surface sits at $T_p$. In reality, gas within one $\lambda$ of the surface has not had enough collisions to get there. The assumption holds only when that layer is negligible compared with the particle.

The free-molecular end: count the molecules that hit#

As $\mathrm{Kn} \to \infty$, gas molecules do not collide with each other before reaching the particle. Heat transfer is the number of molecules striking the surface per unit time times the energy each one hands over.

$\alpha$ is the thermal accommodation coefficient (how fully a reflected molecule adopts the wall temperature), $\bar c$ the mean molecular speed, $\gamma$ the heat capacity ratio, and $m$ the molecular mass. Now $Q_{fm} \propto d^2$. Heat enters through surface area, and heat-up time scales only with $d$.

Taking the ratio exposes the structure.

Substituting $k \sim \rho c_v \bar c \lambda$ shows the ratio scales roughly like $\alpha/\mathrm{Kn}$. The smaller the particle, the smaller the free-molecular rate relative to the continuum rate. Whichever is smaller is the bottleneck.

Sherman's harmonic form: two thermal resistances in series#

In the transition regime, both mechanisms act in sequence. Within about one $\lambda$ of the surface, molecules carry heat ballistically; beyond that, diffusion takes over. Treating the two as series thermal resistances gives Sherman's (1963) harmonic form.

It is exact at both ends. As $\mathrm{Kn} \to 0$, $1/Q_{fm}$ vanishes; as $\mathrm{Kn} \to \infty$, $1/Q_c$ vanishes. The source presentation also lists Fuchs' boundary-sphere approach and the harmonic forms of Fuchs–Sutugin and Loyalka. The boundary-sphere approach wraps the particle in a fictitious sphere of radius $\sim\lambda$, solves free-molecular transport inside and continuum transport outside, and matches the two at the boundary. The harmonic forms compress that result into a single coefficient.

Code: six diameters and four accommodation coefficients#

The code below transcribes the three formulas directly. It integrates the particle energy balance $m c_p\, dT/dt = Q(T)$ from 300 K to aluminum's melting point of 933 K. Integrating in temperature instead of time gives $dt = m c_p\, dT / Q$, so one loop does it. Radiation and surface reactions are deliberately left out.

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

4.58 comes from a one-line formula#

Both rates are proportional to $(T_g - T_p)$, so $Q_{fm}/Q_c$ does not depend on particle temperature. The Sherman heating rate is just the continuum rate times a constant $s = 1/(1 + Q_c/Q_{fm})$. The heat-up time ratio is its inverse.

At 1 μm, $Q_{fm}/Q_c = 0.28$, so $1 + 1/0.28 \approx 4.6$, matching the 4.58 in the table. At 0.1 μm it is $1 + 1/0.03$, about 37 times. The integration loop only confirms this arithmetic.

The harmonic side is the same division. $Q/Q_c=r/(1+r)$ with $r=Q_{fm}/Q_c$. At 1 μm, $0.28/1.28\approx0.219$ against the printed 0.218. At 10 μm, $2.79/3.79=0.736$, the same as the table. At 3 μm, $0.84/1.84\approx0.457$ against 0.456. At 50 μm, $13.96/14.96=0.933$. The print sits inside one rounding step of this expression.

When $Q_{fm}=Q_c$ the time ratio is exactly 2. Heat is cut in half. Time doubles. The 3.6 μm and $\mathrm{Kn}\approx0.19$ in the text are that point. That is already a factor of two, at a diameter most people would still call continuum.

$\alpha$ multiplies $r$ directly. At $\alpha=0.5$, $r=0.14$ and $1+1/0.14\approx8.14$, next to the table's 8.16. At $\alpha=0.1$, $r\approx0.028$ and $1+1/0.028\approx36.7$, next to 36.82. At $\alpha=0.02$, $1+1/0.0056\approx180$, next to 180.09. The same $r/(1+r)$ gives 0.123, 0.027, 0.006 against the printed 0.122, 0.027, 0.006. Not a new run. Just the output block divided again.

The row worth staring at is 10 μm. At $\mathrm{Kn} = 0.067$, a range most people would call continuum, delivered heat drops by 26%. The crossover diameter where $Q_{fm} = Q_c$ is about 3.6 μm, or $\mathrm{Kn} \approx 0.19$. The crossover sits five times below $\mathrm{Kn} = 1$. That comes from $(\gamma+1)/(\gamma-1) = 6$, the factor $1/16$, and the coefficient linking $k$ to $\lambda$. The rule of thumb "Kn below 0.1 means continuum" is too loose for heat transfer.

The accommodation coefficient moves things as much as diameter#

The bottom of the table varies only $\alpha$ at 1 μm. With $\alpha = 0.1$, heat-up takes 37 times longer; with 0.02, 180 times. $Q_{fm}$ enters only through the product $\alpha d$. Cutting $\alpha$ tenfold is the same as cutting the diameter tenfold.

The catch is that $\alpha$ is poorly known. Reported values between metal surfaces and hot gases spread widely depending on oxide layers, temperature ratio, and gas species. Filippov & Rosner (2000) revisited energy transfer itself at large temperature ratios. Some studies set $\alpha$ well below 1 to match nano-aluminum ignition delays. In transition-regime particle models, $\alpha$ easily becomes a tuning knob. The fact that it shifts results by multiples belongs in the report.

Where the harmonic form wobbles#

The Sherman form nails both ends, but the shape in between is an approximation. For $\mathrm{Kn} \sim 0.1$–1, the Fuchs boundary-sphere solution, the Fuchs–Sutugin formula, and Loyalka's correction disagree by tens of percent. Read the 26% at 10 μm as a number inside that model spread.

The definition of $\lambda$ adds uncertainty of the same size. Hard-sphere formula versus viscosity-based definition, $\lambda/r$ versus $2\lambda/d$: the Kn of the same particle can easily shift by tens of percent. When borrowing a correlation from another paper, bring that paper's $\lambda$ definition along. Otherwise the boundary sphere lands in the wrong place.

The bottleneck stays into combustion. Oxidizer diffusion toward the particle also needs a Knudsen correction on the mass side. Agglomerates grow the effective diameter and Kn drops again. Whether Nu = 2 is safe then depends on recomputing Kn with that effective diameter.

Before adding one heat-transfer line to Lagrangian particles#

Particle heat transfer in Euler–Lagrange codes usually starts with one line: $\mathrm{Nu} = 2 + 0.6\,\mathrm{Re}^{1/2}\mathrm{Pr}^{1/3}$. That line is the continuum end. A tail at 3 μm carries the table ratio 2.19. At 1 μm the ratio is 4.58. "Two or three times" only fits near 3 μm. At 1 μm the miss is already past four. If ignition delay comes out short on the small diameters, this line is the first place to look.

The arithmetic is cheap. $\lambda$ per cell, $Q_{fm}$ per particle, harmonic mean. Leave out where $\alpha$ and the $\lambda$ definition came from and 4.58 will not come back. That factor is almost those two numbers. Cranking the heat flux does not keep the slope either. As on the boiling curve where the heat flux turns over, the mechanism changes and the delivered heat bends. Here the bend is $Q_{fm}$.