1 μm 铝颗粒的升温比连续介质公式慢 4.6 倍 — Knudsen 过渡区的 Sherman 调和插值

过渡区颗粒的加热时间是连续介质值的 1 + Q_c/Q_fm 倍。在 1500 K 空气中,1 μm 铝颗粒升温慢 4.6 倍,10 μm 也慢 1.36 倍。

1 μm 颗粒为什么比计算熔化得晚#

把常温铝颗粒扔进 1500 K 的空气。直径 50 μm 和 1 μm 的颗粒,哪个先熔化?当然是小的。质量按 $d^3$ 减小,热传导只按 $d$ 减小。按连续介质公式,1 μm 颗粒 1.95 μs 就到达熔点。

可同一个颗粒用分子动理论来算,结果是 8.95 μs。相差 4.6 倍。原因是颗粒尺寸已经接近周围气体的平均自由程(分子两次碰撞之间飞行的平均距离)。本文讨论铝团聚颗粒燃烧模型的第一步——预热(pre-heating)阶段——这个差距从哪里来。我们分别建立连续介质端和自由分子端的公式,再用 Sherman 调和插值把两者连起来。最后会看到,整个结果可以压缩成一行: $t/t_c = 1 + Q_c/Q_{fm}$ 。

单个颗粒眼中的气体颗粒感

气体看起来是不是连续介质,取决于观察者的尺寸。这个比值就是 Knudsen 数。

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

$\lambda$ 是平均自由程, $d$ 是颗粒直径, $\sigma$ 是分子碰撞直径, $p$ 是压力。1500 K、1 atm 的空气取 $\sigma = 3.7\times10^{-10}$ m,得 $\lambda = 0.336$ μm。大约是常温的五倍。压力不变时温度升高,分子数密度下降。

通常把 $\mathrm{Kn} < 0.01$ 称为连续介质区, $\mathrm{Kn} > 10$ 称为自由分子区。两者之间是过渡区(transition regime)。50 μm 颗粒的 $\mathrm{Kn} = 0.013$ ,1 μm 颗粒是 0.67。相当一部分微米级铝粉落在过渡区。

在下面的模拟中亲手调整颗粒直径和适应系数。

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

把直径滑块往下拉,左侧分子的直线飞行距离会超过颗粒本身。被颗粒冷却的分子(蓝色)在连续介质一侧停留在颗粒附近,在自由分子一侧则径直飞远。右下角的竞赛中,两根条的到达时间之比等于红点纵坐标的倒数,也可以验证一下。

连续介质端:热量只靠扩散进来

$\mathrm{Kn} \to 0$ 时,颗粒周围的气体服从 Fourier 定律。静止气体中球体的稳态导热给出 Nusselt 数(Nu,对流换热与纯导热之比)等于 2。

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

$k$ 是气体导热系数, $T_g$ 和 $T_p$ 是气体和颗粒温度。热量只按直径正比进入。质量按 $d^3$ ,时间就按 $d^2$ 。直径缩到十分之一,升温快一百倍。

表里的 $t_c$ 把这个平方印了出来。50 μm 是 $4.883\times10^{-3}$ s,1 μm 是 $1.953\times10^{-6}$ s。直径比 50 的平方是 2500,时间比也是 2500。10 μm 对 1 μm、3 μm 对 0.3 μm,直径比 10,时间比 100。 $1.953\times10^{-4}/1.953\times10^{-6}=100$ 。 $1.758\times10^{-5}/1.758\times10^{-7}=100$ 。在连续介质端,只看直径就能把加热时间比算死。液滴变干的时间也是同一个 $d^2$ 。这个平衡在扩散控制的蒸发里又出现一次。

这个公式假设颗粒表面处的气体温度等于 $T_p$ 。实际上,表面以内约一个 $\lambda$ 的气体并没有经历足够多的碰撞。只有当这层厚度相对颗粒可以忽略时,这个假设才成立。

自由分子端:数一数撞上来的分子

$\mathrm{Kn} \to \infty$ 时,气体分子在到达颗粒之前彼此不碰撞。换热量等于单位时间撞击表面的分子数乘以每个分子交出的能量。

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$ 是热适应系数(thermal accommodation coefficient,碰撞分子离开时向壁面温度靠拢的程度), $\bar c$ 是平均分子速率, $\gamma$ 是比热比, $m$ 是分子质量。这一次 $Q_{fm} \propto d^2$ 。热量经表面积进入,加热时间只与 $d$ 成正比。

取两式之比,结构就清楚了。

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

代入 $k \sim \rho c_v \bar c \lambda$ ,这个比值大致与 $\alpha/\mathrm{Kn}$ 成正比。颗粒越小,自由分子公式相对连续介质公式越小。较小的那个就是瓶颈。

Sherman 调和插值:两个热阻串联#

在过渡区,两种机制先后起作用。从颗粒表面到约一个 $\lambda$ 处,分子靠飞行传热;再往外由扩散传热。把这两段看作串联热阻,就得到 Sherman(1963)的调和形式。

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

两端都是精确的。 $\mathrm{Kn} \to 0$ 时 $1/Q_{fm}$ 消失, $\mathrm{Kn} \to \infty$ 时 $1/Q_c$ 消失。原始报告除了这个形式,还整理了 Fuchs 的边界球(boundary sphere)方法,以及 Fuchs–Sutugin 和 Loyalka 的调和形式。边界球方法用半径约为 $\lambda$ 的虚拟球包住颗粒,内部按自由分子、外部按连续介质求解,在边界上匹配两个解。调和形式是把这一结果压缩成一个系数的近似。

代码:六个直径和四个适应系数

下面的代码直接搬用上面三个公式。对颗粒能量平衡 $m c_p\, dT/dt = Q(T)$ 从 300 K 积分到铝的熔点 933 K。不按时间而按温度积分, $dt = m c_p\, dT / Q$ ,一个循环就够了。辐射和表面反应特意没有加入。

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 来自一行公式#

两个公式都与 $(T_g - T_p)$ 成正比,所以 $Q_{fm}/Q_c$ 与颗粒温度无关。于是 Sherman 加热率就是连续介质加热率乘以常数 $s = 1/(1 + Q_c/Q_{fm})$ 。加热时间之比是它的倒数。

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

1 μm 时 $Q_{fm}/Q_c = 0.28$ ,所以 $1 + 1/0.28 \approx 4.6$ ,与表中的 4.58 一致。0.1 μm 时是 $1 + 1/0.03$ ,约 37 倍。积分循环只是验证了这道算术。

调和那一侧是同一道除法。 $Q/Q_c=r/(1+r)$ , $r=Q_{fm}/Q_c$ 。1 μm 为 $0.28/1.28\approx0.219$ ,表上是 0.218。10 μm 为 $2.79/3.79=0.736$ ,与表相同。3 μm 为 $0.84/1.84\approx0.457$ ,表上是 0.456。50 μm 为 $13.96/14.96=0.933$ 。输出落在这个式子四舍五入的一格之内。

$Q_{fm}=Q_c$ 时时间比正好是 2。热量减半,时间加倍。正文里的 3.6 μm、 $\mathrm{Kn}\approx0.19$ 就是这一点。直径还早不到通常所说的连续介质,时间已经是两倍。

$\alpha$ 也直接乘在 $r$ 上。 $\alpha=0.5$ 时 $r=0.14$ , $1+1/0.14\approx8.14$ ,表上是 8.16。 $\alpha=0.1$ 时 $r\approx0.028$ , $1+1/0.028\approx36.7$ ,表上是 36.82。 $\alpha=0.02$ 时 $1+1/0.0056\approx180$ ,表上是 180.09。 $r/(1+r)$ 给出 0.123、0.027、0.006,贴着表上的 0.122、0.027、0.006。不是新的计算。只是把输出块再除一遍。

更值得注意的是 10 μm 那一行。 $\mathrm{Kn} = 0.067$ ,通常被归入连续介质范围,传入的热量却少了 26%。 $Q_{fm} = Q_c$ 的交叉直径约为 3.6 μm,即 $\mathrm{Kn} \approx 0.19$ 。交叉点位于比 $\mathrm{Kn} = 1$ 小五倍的地方。这来自 $(\gamma+1)/(\gamma-1) = 6$ 、 $1/16$ ,以及联系 $k$ 与 $\lambda$ 的系数。"Kn 小于 0.1 就是连续介质"这一经验,用在换热上太宽松了。

适应系数 α 的影响和直径一样大

表的下半部分是在 1 μm 下只改变 $\alpha$ 的结果。 $\alpha = 0.1$ 时加热时间变为 37 倍,0.02 时变为 180 倍。 $Q_{fm}$ 只以乘积 $\alpha d$ 的形式出现。把 $\alpha$ 减小十倍,等同于把直径减小十倍。

问题在于 $\alpha$ 并不清楚。金属表面与高温气体之间的报告值,随表面氧化膜、温度比和气体种类而分布很宽。Filippov & Rosner(2000)重新研究了大温度比条件下的能量传递本身。也有研究为了匹配纳米铝的点火延迟,把 $\alpha$ 取得远小于 1。在过渡区颗粒模型里, $\alpha$ 很容易沦为调参旋钮。它能让结果变化数倍,这一点应写进报告。

调和插值不稳的区间

Sherman 形式两端都对,中间的形状却是近似。在 $\mathrm{Kn} \sim 0.1$ –1,Fuchs 边界球解、Fuchs–Sutugin 公式和 Loyalka 修正彼此相差几十个百分点。上面 10 μm 的 26% 也应理解为落在这个模型差异范围内的数字。

$\lambda$ 的定义带来同样量级的不确定性。硬球公式还是基于黏性的定义,Kn 写成 $\lambda/r$ 还是 $2\lambda/d$ ,同一颗粒的 Kn 很容易相差几十个百分点。从别的论文借用关联式时,要连同那篇论文的 $\lambda$ 定义一起借来。否则边界球的半径会放错位置。

进入燃烧,瓶颈还在。颗粒旁边的氧化剂扩散同样要做传质一侧的 Knudsen 修正。团聚以后有效直径变大,Kn 又变小。那时能不能退回 Nu=2,要用这个有效直径把 Kn 重算一遍才知道。

给 Lagrangian 颗粒写下那一行换热公式之前#

Euler–Lagrange 代码中的颗粒换热,通常从一行 $\mathrm{Nu} = 2 + 0.6\,\mathrm{Re}^{1/2}\mathrm{Pr}^{1/3}$ 开始。这行公式站在连续介质端。尾部在 3 μm 时,表上的比值是 2.19。1 μm 则是 4.58。两三倍的感觉只贴在 3 μm 附近,到了 1 μm 已经超过四倍。小直径上的点火延迟短得反常,就先看这一行。

修正的运算很轻。单元的 $\lambda$ ,颗粒的 $Q_{fm}$ ,两者的调和平均。不记下 $\alpha$ 和 $\lambda$ 的定义,4.58 就复现不出来。这个倍数几乎就是这两个数。把热流再抬高,斜率也不会保持。沸腾曲线折下的临界热流一样,机制一换,传热量就弯了。这里的弯折是 $Q_{fm}$ 。