1 μm 알루미늄 입자가 연속체 공식보다 4.6배 늦게 데워졌다 — Knudsen 전이영역의 Sherman 조화 보간

1 μm 입자는 왜 계산보다 늦게 녹았을까#

1500 K 공기에 상온의 알루미늄 입자를 던진다. 지름 50 μm 입자와 1 μm 입자 중 어느 쪽이 먼저 녹을까. 답은 당연히 작은 쪽이다. 질량은 $d^3$로 줄고 열전도는 $d$로만 줄기 때문이다. 연속체 공식으로는 1 μm 입자가 1.95 μs 만에 녹는점에 닿는다.

그런데 같은 입자를 운동론으로 계산하면 8.95 μs가 나온다. 4.6배 차이다. 입자가 주변 기체의 평균자유행로(mean free path, 분자가 다음 충돌까지 날아가는 평균 거리)와 비슷한 크기가 됐기 때문이다. 이 글은 알루미늄 응집입자 연소 모델의 첫 단계인 예열(pre-heating)에서 이 차이가 어디서 오는지 다룬다. 연속체 끝과 자유분자 끝의 공식을 각각 세우고, Sherman의 조화 보간으로 둘을 잇는다. 마지막에 그 결과가 한 줄짜리 식 $t/t_c = 1 + Q_c/Q_{fm}$로 줄어드는 것을 확인한다.

입자 하나가 보는 기체의 입자성#

기체가 연속체로 보이는지는 관찰자 크기에 달려 있다. 그 비율이 Knudsen 수다.

$\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이다. 마이크로 알루미늄 분말의 상당 부분이 전이영역에 걸친다.

아래 시뮬레이션에서 입자 지름과 적응계수를 직접 바꿔보자.

지름 슬라이더를 내리면 왼쪽 분자의 직선 비행 거리가 입자보다 길어진다. 입자에 닿아 식은 분자(파란색)가 연속체 쪽에서는 입자 곁에 머물고, 자유분자 쪽에서는 곧장 멀리 날아간다. 오른쪽 아래 경주에서 두 막대의 도착 시간 비가 빨간 점의 세로값의 역수와 같다는 점도 확인해보자.

연속체 끝: 열은 확산으로만 들어온다#

$\mathrm{Kn} \to 0$이면 입자 주변 기체는 Fourier 법칙을 따른다. 정지한 기체 속 구의 정상 열전도는 Nusselt 수(Nu, 대류 열전달을 전도로 나눈 무차원 수) 2를 준다.

$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$면 기체 분자는 입자에 닿기 전까지 서로 부딪히지 않는다. 열전달은 단위 시간에 표면에 닿는 분자 수와 분자 하나가 넘겨주는 에너지의 곱이다.

$\alpha$는 열 적응계수(thermal accommodation coefficient, 부딪힌 분자가 벽 온도로 얼마나 맞춰져 나가는지의 비율), $\bar c$는 평균 분자 속력, $\gamma$는 비열비, $m$은 분자 질량이다. 이번에는 $Q_{fm} \propto d^2$다. 열이 표면적을 통해 들어온다. 가열 시간은 $d$에 비례할 뿐이다.

두 식의 비를 구하면 구조가 드러난다.

$k \sim \rho c_v \bar c \lambda$를 넣으면 이 비는 대략 $\alpha/\mathrm{Kn}$에 비례한다. 작은 입자일수록 자유분자 식이 연속체 식보다 작아진다. 둘 중 작은 쪽이 병목이다.

Sherman 조화 보간: 두 열저항을 직렬로 잇는다#

전이영역에서는 두 기구가 차례로 작동한다. 입자 표면에서 한 $\lambda$까지는 분자가 날아서 열을 나르고, 그 바깥은 확산이 나른다. 두 단계를 직렬 열저항으로 보면 Sherman(1963)의 조화 형태가 나온다.

양 끝에서 정확하다. $\mathrm{Kn} \to 0$이면 $1/Q_{fm}$이 사라지고, $\mathrm{Kn} \to \infty$면 $1/Q_c$가 사라진다. 원문 발표는 이 형태 외에 Fuchs의 경계구(boundary sphere) 접근과 Fuchs–Sutugin, Loyalka의 조화형을 함께 정리한다. 경계구 접근은 입자를 반지름 $\sim\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})$를 곱한 것이다. 가열 시간의 비는 그 역수다.

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. $Q/Q_c=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$ 정의도 같은 크기의 불확실성을 만든다. 강체구 공식, 점성 기반 정의, $\mathrm{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}$이다.