빨래는 왜 5도에서도 마르는가 — 확산 지배 증발과 열 지배 비등

끓지 않는 증발과 끓는 증발이 서로 다른 소스텀을 푸는 이유

베란다에 널어둔 빨래가 마른다. 기온은 5도다. 물의 끓는점은 100도다. 95도가 모자란데 물은 사라진다.

같은 물이 주전자 안에서는 100도가 되어야 사라진다. 두 경우 모두 액체가 기체로 바뀌는 상변화(phase change)다. 그런데 속도를 결정하는 것은 서로 다르다. CFD에서도 마찬가지다. 두 상황에 같은 소스텀을 쓰면 한쪽은 반드시 틀린다.

이 글은 그 둘을 갈라놓는다. 확산이 병목인 증발과 열이 병목인 증발을 각각 다른 수식으로 세우고, 그것을 Navier-Stokes 방정식의 소스텀으로 옮기는 방법까지 간다. 마지막에는 Lee 모델의 완화계수 $r$ 을 어떻게 골라야 하는지 — 그리고 왜 문헌의 $r$ 값이 일곱 자릿수에 걸쳐 흩어져 있는지 — 를 직접 조작해서 확인한다.

두 증발은 병목이 다르다#

증발이 일어나려면 두 가지가 필요하다. 분자가 계면을 떠날 만큼의 에너지, 그리고 떠난 분자를 계면에서 치워줄 통로. 어느 쪽이 부족하냐가 레짐을 가른다.

확산 지배(diffusion-limited). 계면 바로 위 기체는 이미 포화 상태다. 더 증발하려면 그 증기가 멀리 실려 나가야 한다. 병목은 기체 쪽 확산이다. 액체 온도는 거의 변하지 않는다. 물웅덩이, 젖은 빨래, 상온의 연료 분무가 여기 속한다.

열 지배(thermally-driven). 계면 근처 액체가 이미 포화 온도를 넘겼다. 증기가 실려 나갈 통로는 넘친다. 병목은 잠열(latent heat, 상변화에 드는 숨은 열)을 대줄 열 공급이다. 비등, 캐비테이션, 벽면 응축이 여기 속한다.

두 레짐을 구별하는 무차원수도 다르다. 확산 지배는 Spalding 물질전달수 $B_M$ (계면과 원방의 증기 농도 차)로, 열 지배는 Jakob 수 $\mathrm{Ja} = c_p \Delta T / L$ (현열/잠열 비)로 잰다. 앞의 것은 Schmidt 수 $\mathrm{Sc} = \nu/D$ 와 짝을 이루고, 뒤의 것은 Prandtl 수와 짝을 이룬다.

확산 지배 — 농도 구배가 속도를 정한다#

구형 액적 표면에서 나가는 증기 질량 플럭스를 세운다. Fick의 법칙(농도 구배에 비례하는 확산)만 쓰면 부족하다. 증기가 나가면서 기체 전체가 바깥으로 밀려나기 때문이다. 이 대류 성분을 Stefan 유동이라 부른다.

\dot{m}'' = -\rho_g D \frac{\partial Y_v}{\partial r}\bigg|_{s} + \dot{m}'' Y_{v,s}

$\dot{m}''$ 은 면적당 증발률, $\rho_g$ 는 기체 밀도, $D$ 는 증기의 확산계수, $Y_v$ 는 증기 질량분율, 아래첨자 $s$ 는 계면이다. 우변 둘째 항이 Stefan 유동이 다시 실어 오는 몫이다.

준정상 가정으로 반경 방향으로 적분하면 로그가 튀어나온다.

\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$ 는 액적 반경, $Y_s$ 는 계면 포화 질량분율, $Y_\infty$ 는 원방 질량분율이다. $B_M$ 은 Spalding 물질전달수 — 증발을 밀어내는 농도 낙차를 하나로 압축한 값이다.

여기에 액적 질량 보존 $\dot{m} = -\frac{d}{dt}(\rho_l \pi d^3/6)$ 을 붙이면 지름의 제곱이 선형으로 줄어든다.

\frac{d(d^2)}{dt} = -K, \qquad K = \frac{8\rho_g D}{\rho_l}\ln(1 + B_M)

이것이 d²-법칙이다. $d$ 가 아니라 $d^2$ 가 직선이다. 수명은 $t_{life} = d_0^2 / K$ .

핵심은 $Y_s$ 가 온도만의 함수라는 점이다. 25도 물의 포화 증기압은 3.17 kPa, 대기압의 3 %다. 질량분율로 $Y_s \approx 0.019$ . 0은 아니다. 그래서 마른다. 5도에서는 $Y_s \approx 0.0054$ 로 줄지만 여전히 0이 아니다. 끓는점은 이 이야기에 등장하지 않는다.

아래 시뮬레이션에서 직접 조작해보자.

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.

RH 슬라이더를 30 %에서 95 %로 올려보면 $Y_\infty$ 가 $Y_s$ 를 향해 올라오면서 $B_M$ 이 무너지고, 오른쪽 직선의 기울기가 눕는다. 1 mm 물방울의 수명이 5분에서 한 시간 가까이로 늘어난다. 반대로 온도를 올리면 $Y_s$ 가 지수적으로 달아나 같은 방울이 수 초 만에 사라진다. 왼쪽 원의 반지름은 오른쪽 그래프와 같은 $d(t)$ 로 그려져 있다.

열 지배 — 도착한 열이 잠열을 갚는다#

이제 계면이 포화 온도에 붙어 있는 경우다. 증발률을 정하는 것은 농도가 아니라 계면에 도착하는 열의 불균형이다.

\dot{m}'' L = k_l \frac{\partial T}{\partial n}\bigg|_l - k_v \frac{\partial T}{\partial n}\bigg|_v

$L$ 은 잠열, $k_l$ · $k_v$ 는 액체·증기 열전도도, $n$ 은 계면 법선이다. 액체 쪽에서 들어온 열에서 증기 쪽으로 빠져나간 열을 뺀 나머지가 전부 상변화로 간다.

이것을 Stefan 조건이라 부른다. 까다로운 이유가 있다. 계면에 조건이 둘 걸리기 때문이다. 온도는 $T = T_{sat}$ 로 고정(Dirichlet)이고, 동시에 플럭스 점프가 계면 속도를 정한다. 조건이 두 개인데 경계 위치는 미지수다. 자유경계 문제다.

Sharp-interface로 정직하게 풀려면 계면을 추적하고 매 스텝 위 식을 만족시켜야 한다. 대부분의 상용·오픈소스 solver는 그 길을 피한다. 대신 계면을 한두 셀 두께로 번지게 두고, 위 조건이 저절로 만족되도록 유도하는 체적 소스텀을 넣는다.

소스텀으로 옮기면 어디에 무엇이 들어가는가#

체적 기화율 $\dot{m}'''$ [kg/m³/s]을 하나 정의하면, 지배방정식 네 곳에 각각 들어간다.

상별 연속방정식에 부호를 뒤집어 한 쌍으로 들어간다.

\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$ 는 체적분율이다. 두 식을 더하면 소스가 상쇄되어 혼합 질량은 보존된다.

에너지 방정식에는 잠열만큼 싱크로 들어간다. 증기 질량분율 방정식을 따로 푸는 확산 지배 solver라면 거기에도 같은 값이 소스로 들어간다.

\frac{\partial (\rho e)}{\partial t} + \nabla\cdot(\rho \mathbf{u} h) = \nabla\cdot(k\nabla T) - \dot{m}''' L

가장 자주 빠뜨리는 곳은 네 번째다. 상변화는 속도장의 발산 제약을 바꾼다.

\nabla\cdot\mathbf{u} = \dot{m}''' \left( \frac{1}{\rho_v} - \frac{1}{\rho_l} \right)

비압축성 solver의 압력 Poisson 방정식은 보통 우변이 0이다. 상변화가 있으면 0이 아니다. 대기압에서 물이 증기로 바뀌면 부피가 약 1600배가 된다. 이 항을 빠뜨린 채 기화율만 넣으면 질량은 사라지는데 아무것도 밀려나지 않는다. 기포는 자라지 않고 압력만 이상해진다.

$r$ 을 고르는 문제 — Lee 모델의 함정#

가장 널리 쓰이는 소스텀이 Lee 모델이다.

\dot{m}''' = r\,\alpha_l \rho_l \frac{T - T_{sat}}{T_{sat}}

$r$ 은 완화계수(relaxation coefficient) [1/s]다. 여기서 오해가 시작된다. $r$ 은 물성이 아니다. 실험으로 잴 수 있는 값도 아니다. $r$ 의 역할은 계면 셀의 온도를 $T_{sat}$ 에 붙여놓는 것뿐이다. 온도가 $T_{sat}$ 에 붙으면 그 셀의 에너지 수지가 저절로 Stefan 조건이 된다. 즉 $r$ 은 앞 절의 자유경계 조건을 페널티 방식으로 강제하는 수치 장치다.

그래서 $r$ 은 크면 클수록 이론적으로는 좋다. 문제는 명시적으로 적분할 때 생긴다. 소스 항만 떼어 보면 $T - T_{sat}$ 에 대한 선형 감쇠고, 그 감쇠율이 $r L/(T_{sat} c_p)$ 다. 명시적 Euler의 안정 조건은 다음과 같다.

S = r\,\frac{L}{T_{sat}\,c_p}\,\Delta t < 2

물–증기에서 $L/(T_{sat} c_p) \approx 1.43$ 이다. $\Delta t = 10^{-5}$ s면 $r$ 은 $1.4\times10^{5}$ 를 넘을 수 없다. 격자를 줄여 $\Delta t$ 가 작아지면 허용 $r$ 이 커진다. 문헌의 $r$ 값이 0.1부터 $10^7$ 까지 흩어져 있는 이유가 이것이다. 저자마다 $\Delta t$ 가 다르다.

아래 열 지배 컬럼에서 $r$ 을 직접 쓸어보자.

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.

$r$ 을 $10^2$ 까지 내리면 노란 음영 — 기화로 바뀌지 못한 과열 — 이 부풀고 계면이 멈춘다. 벽은 계속 열을 붓는데 아무것도 끓지 않는다. $10^5$ 을 넘기면 $S$ 가 2를 지나 온도 프로파일이 울리고, 오른쪽 이력 그래프가 톱니가 된다. 초록 구간은 그 사이 두 자릿수뿐이고, $\Delta t$ 나 격자를 바꾸면 그 구간도 통째로 움직인다.

물성 기반 대안이 없지는 않다. Hertz–Knudsen–Schrage 관계는 기체 운동론에서 출발해 계면 증발률을 수용계수(accommodation coefficient) $\sigma$ 하나로 표현한다. Tanasawa 모델은 그것을 선형화해 $\dot{m}'' \propto \sigma (T - T_{sat})$ 형태로 만든다. $\sigma$ 는 물성이라 실험값이 있다. 다만 물의 $\sigma$ 실측치도 0.01에서 1 사이로 흩어져 있어서, 불확실성이 사라진다기보다 자리를 옮긴다.

Python으로 겹쳐 보는 두 레짐#

한쪽은 습도가, 다른 쪽은 시간 스텝이 답을 정한다. 같은 스크립트에서 확인한다.

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     # 물, 습공기, 수증기 확산계수
L_VAP, CP_L, T_SAT = 2.26e6, 4220.0, 373.15
 
def p_sat(t_c):
    """Antoine 식 (물, mmHg) -> Pa"""
    return 10 ** (8.07131 - 1730.63 / (233.426 + t_c)) * 133.322
 
def mass_fraction(x):
    """증기 몰분율 -> 질량분율"""
    return x * M_V / (x * M_V + (1 - x) * M_A)
 
def spalding_number(t_c, rh):
    """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):
    """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):
    """Lee 모델 체적 기화율 [kg/m^3/s]"""
    if temp > T_SAT:
        return r * alpha_l * RHO_L * (temp - T_SAT) / T_SAT
    return 0.0
 
print("확산 지배 — 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"  RH {rh*100:4.0f} %   B_M = {spalding_number(25.0, rh):.5f}"
          f"   K = {k:.2e} mm^2/s   수명 = {life/60:7.2f} min")
 
print("\n열 지배 — Lee 완화계수와 명시적 안정 한계")
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 과열) = {rate:9.2e} kg/m^3/s"
          f"   S = {s:8.3f}  {'발산' if s > 2 else '안정'}")

출력은 이렇다.

확산 지배 — 1 mm 물방울, 25 degC
  RH    0 %   B_M = 0.02001   K = 4.69e-03 mm^2/s   수명 =    3.55 min
  RH   30 %   B_M = 0.01406   K = 3.30e-03 mm^2/s   수명 =    5.04 min
  RH   60 %   B_M = 0.00806   K = 1.90e-03 mm^2/s   수명 =    8.77 min
  RH   90 %   B_M = 0.00202   K = 4.78e-04 mm^2/s   수명 =   34.85 min
  RH   97 %   B_M = 0.00061   K = 1.44e-04 mm^2/s   수명 =  115.98 min
 
열 지배 — Lee 완화계수와 명시적 안정 한계
  r =    1e+02   m''' (1 K 과열) =  2.67e+02 kg/m^3/s   S =    0.001  안정
  r =    1e+03   m''' (1 K 과열) =  2.67e+03 kg/m^3/s   S =    0.014  안정
  r =    1e+04   m''' (1 K 과열) =  2.67e+04 kg/m^3/s   S =    0.144  안정
  r =    1e+05   m''' (1 K 과열) =  2.67e+05 kg/m^3/s   S =    1.435  안정
  r =    1e+06   m''' (1 K 과열) =  2.67e+06 kg/m^3/s   S =   14.352  발산

습도 30 %와 97 % 사이에서 수명이 23배 벌어진다. 온도는 그대로다. 확산 지배 레짐에서 계산이 실험과 안 맞는다면 대개 원방 습도 경계조건이 범인이다.

실험실 노트로 옮길 것#

증발 속도를 정하는 것이 무엇인지 먼저 정하고 모델을 고른다. 기체 쪽 농도 구배면 $B_M$ 과 d²-법칙, 계면 열 불균형이면 Stefan 조건과 Lee/Tanasawa다. 순서를 바꾸면 맞는 답이 우연히 나올 뿐이다.
Lee 모델의 $r$ 은 물성이 아니라 페널티 계수다. 값을 문헌에서 베끼지 말고 $S = r L \Delta t/(T_{sat} c_p) < 2$ 안에서 가능한 크게 잡은 뒤, 계면 셀의 과열이 실제로 0에 붙는지 확인한다.
기화율을 넣었으면 $\nabla\cdot\mathbf{u} = \dot{m}'''(1/\rho_v - 1/\rho_l)$ 도 같이 넣는다. 이걸 빠뜨리면 질량은 사라지는데 부피는 그대로다. 물–증기에서 그 오차는 1600배짜리다.