항력계수를 9% 깎았더니 액적이 35% 더 느려졌다 — 변형과 기화 보정이 붙는 자리

1983년, 불타는 액적 앞에서 항력계수가 맞지 않았다#

Renksizbulut와 Yuen은 고온 기류 속에서 기화하는 액적의 항력을 쟀다. 구형 상관식이 예측한 값보다 실측 항력이 꾸준히 낮았다. 액적이 구형을 유지하고 있었는데도 그랬다. 표면에서 뿜어져 나오는 증기가 경계층을 바깥으로 밀어내고 있었던 것이다.

같은 시기에 반대 방향의 어긋남도 보고되고 있었다. 고속 기류에 놓인 액적은 납작해지고, 그러면 항력은 구형보다 커진다. O'Rourke와 Amsden이 1987년에 TAB 모델로 정리한 것이 그쪽이다.

그래서 스프레이 코드의 액적 항력계수는 대개 세 겹이다. 구형 상관식 하나에 변형 보정 하나, 기화 보정 하나가 곱해진다. 이 글은 세 겹을 하나씩 벗겨서 각각이 궤적을 얼마나 움직이는지 같은 조건에서 잰다. 결론부터 말하면, 계수를 가장 많이 깎은 모델이 가장 먼저 멈춘다.

세 상관식은 같은 자리를 서로 다르게 고친다#

출발점은 구형 액적의 항력계수다. Schiller–Naumann 계열의 상관식을 쓴다.

$C_{D,\mathrm{sph}} = \frac{24}{Re}\left(1 + \tfrac{1}{6} Re^{2/3}\right), \qquad Re < 1000$

$Re = \rho_g u_{\mathrm{rel}} d / \mu_g$는 액적 지름과 상대속도로 만든 레이놀즈 수다. 앞의 $24/Re$가 Stokes 항력이고, 괄호 안이 관성 보정이다.

여기에 붙는 보정 두 개는 성격이 다르다.

$y$는 TAB 모델의 무차원 변형량이다. $y = 0$이면 구, $y = 1$이면 분열 판정이다. $B_M$은 Spalding 물질전달 수(표면과 원방의 증기 질량분율 차이를 무차원화한 값)로, 증발 소스항과 $d^2$ 법칙에서 다룬 그 수다.

변형이 실제로 어떻게 자라는지 아래에서 직접 조작해보자.

상대속도와 지름을 올리면 방울이 원반으로 눌리고 계수 막대가 회색 구형 값을 넘어선다. 관찰할 것은 점선으로 그은 정상상태 $y_{ss}$를 실제 $y$가 얼마나 넘어서는가다.

변형은 정면 면적이 아니라 계수에 얹힌다#

TAB은 액적을 감쇠 조화진동자로 본다. 기체 동압이 밀고, 표면장력이 되돌리고, 액체 점성이 감쇠한다.

$r$은 액적 반지름, $\sigma$는 표면장력, $\mu_l$은 액체 점성이다. 상수는 $C_F = 1/3$, $C_k = 8$, $C_d = 5$, $C_b = 1/2$를 쓴다. 이 방정식이 분열까지 밀고 가는 과정은 Weber 수와 TAB 분열에서 다뤘다.

미분항을 모두 0으로 두면 정상 변형량이 나온다.

$y_{ss} = \frac{C_F}{C_b C_k}We_r = \frac{We_r}{12}, \qquad We_r = \frac{\rho_g u_{\mathrm{rel}}^2 r}{\sigma}$

여기서 실무자가 자주 틀리는 지점이 하나 있다. $y_{ss}$를 그대로 항력 보정에 넣는 코드가 흔하다. 그런데 이 진동자는 거의 감쇠되지 않는다. n-헵탄 40 μm 액적의 감쇠비는 $10^{-2}$ 수준이다. 정지 상태에서 출발한 $y$는 $y_{ss}$를 지나쳐 거의 두 배까지 올라간다.

즉 정상상태 값으로 계산한 항력은 실제 첫 진동의 절반에 불과하다. 계수가 아니라 그 계수를 쓰는 시점이 문제다.

기화는 경계층을 밀어내서 계수를 깎는다#

Renksizbulut–Yuen 상관식은 기화를 곱셈 보정으로 처리한다.

$C_D = \frac{24}{Re}\left(1 + 0.2\,Re^{0.63}\right)(1+B_M)^{-0.2}$

지수가 $-0.2$로 작다. $B_M = 0.6$이면 인자는 $0.91$, $B_M = 2$면 $0.80$이다. 격렬하게 타는 액적이라도 계수는 20% 남짓 깎일 뿐이다.

원리는 막 이론(film theory)과 같다. 표면에서 나가는 증기 유속이 경계층을 두껍게 만든다. 속도 구배가 완만해지고, 벽면 전단이 줄고, 항력이 준다. 열전달의 Nusselt 수에도 같은 형태의 보정이 붙는다.

Python으로 네 액적을 같은 노즐에서 쏘았다#

정지 기체에 액적 하나를 쏘고 1 ms 동안 추적한다. 지배식은 1차원이다.

$\frac{du_d}{dt} = -\frac{3}{4}\frac{\rho_g}{\rho_l}\frac{C_D}{d}\,u_d|u_d|$

네 케이스는 항력계수만 다르다. A는 구형, B는 구형에 TAB 변형을 곱한 것, C는 기화 보정만, D는 기화 보정에 $d^2$ 법칙 수축까지 넣은 것이다.

import math
 
RHO_G, MU_G = 4.98, 3.3e-5      # 700 K, 10 bar 공기
RHO_L, SIG, MU_L = 620.0, 0.0145, 3.0e-4   # n-헵탄
D0, U0 = 40e-6, 25.0            # 분열 이후 액적, 분사 속도
K_EVAP = 1.0e-6                 # d^2 법칙 상수 [m^2/s]
CF, CK, CD_T, CB = 1.0 / 3.0, 8.0, 5.0, 0.5   # TAB 상수
 
def cd_sphere(re):
    if re < 1e-8:
        return 0.0
    return 24.0 / re * (1.0 + re ** (2.0 / 3.0) / 6.0) if re < 1000.0 else 0.424
 
def cd_distorted(re, y):
    return cd_sphere(re) * (1.0 + 2.632 * max(0.0, min(1.0, y)))
 
def cd_blowing(re, bm):
    if re < 1e-8:
        return 0.0
    return 24.0 / re * (1.0 + 0.2 * re ** 0.63) * (1.0 + bm) ** -0.2
 
def tab_oscillate(y, ydot, u, d, dt):
    r = 0.5 * d
    acc = (CF / CB) * RHO_G * u * u / (RHO_L * r * r) \
        - (CK * SIG / (RHO_L * r ** 3)) * y \
        - (CD_T * MU_L / (RHO_L * r * r)) * ydot
    ydot += acc * dt
    y += ydot * dt
    return y, ydot
 
def run_droplet(law, evaporating, bm=0.6, t_end=1.0e-3, dt=5e-8):
    u, x, d, t = U0, 0.0, D0, 0.0
    y, ydot, ymax = 0.0, 0.0, 0.0
    cd_sum, cd0, n = 0.0, None, 0
    while t < t_end and u > 1e-6:
        re = RHO_G * u * d / MU_G
        if law == 'sphere':
            cd = cd_sphere(re)
        elif law == 'tab':
            y, ydot = tab_oscillate(y, ydot, u, d, dt)
            ymax = max(ymax, y)
            cd = cd_distorted(re, y)
        else:
            cd = cd_blowing(re, bm)
        if cd0 is None:
            cd0 = cd
        cd_sum += cd
        n += 1
        u += -0.75 * RHO_G * cd * u * u / (RHO_L * d) * dt
        x += u * dt
        if evaporating:
            d = math.sqrt(max(1e-18, d * d - K_EVAP * dt))
        t += dt
    return dict(x_mm=x * 1e3, u=u, d_um=d * 1e6, cd0=cd0,
                cd_mean=cd_sum / n, ymax=ymax)
 
CASES = [
    ('A sphere            ', 'sphere', False),
    ('B sphere x TAB      ', 'tab', False),
    ('C blowing (B_M=0.6) ', 'blow', False),
    ('D blowing + d2-law  ', 'blow', True),
]
 
print('Re0 = %.0f   We_r = %.2f   t = 1.0 ms' %
      (RHO_G * U0 * D0 / MU_G, RHO_G * U0 * U0 * 0.5 * D0 / SIG))
print('case                  Cd(t=0)  <Cd>    u(1ms)   x(1ms)   d(1ms)')
base = None
for name, law, ev in CASES:
    r = run_droplet(law, ev)
    if base is None:
        base = r
    print('%s  %6.3f  %6.3f  %6.2f   %6.3f   %5.1f  (%+.1f%% x)' %
          (name, r['cd0'], r['cd_mean'], r['u'], r['x_mm'], r['d_um'],
           100.0 * (r['x_mm'] / base['x_mm'] - 1.0)))
 
tab = run_droplet('tab', False)
print('TAB peak distortion y_max = %.3f  ->  Cd multiplier %.2fx' %
      (tab['ymax'], 1.0 + 2.632 * tab['ymax']))

Re0 = 151   We_r = 4.29   t = 1.0 ms
case                  Cd(t=0)  <Cd>    u(1ms)   x(1ms)   d(1ms)
A sphere               0.910   1.708    3.36    9.261    40.0  (+0.0% x)
B sphere x TAB         0.910   2.233    2.66    7.490    40.0  (-19.1% x)
C blowing (B_M=0.6)    0.828   1.537    3.68    9.732    40.0  (+5.1% x)
D blowing + d2-law     0.828   2.096    2.18    8.868    24.5  (-4.3% x)
TAB peak distortion y_max = 0.632  ->  Cd multiplier 2.66x

변형 보정은 관통거리를 19.1% 줄인다. 첫 진동에서 계수가 2.66배까지 올라가기 때문이다. 기화 보정만 걸면 반대로 5.1% 늘어난다. 여기까지는 각 상관식의 부호 그대로다.

계수가 가장 낮은 쪽이 가장 먼저 멈춘 이유#

문제는 D다. C와 똑같은 항력계수 식을 쓴다. 시작 계수도 0.828로 같다. 그런데 1 ms 뒤 속도가 2.18 m/s로 넷 중 가장 낮다. A의 3.36 m/s보다 35% 느리다.

감속식을 다시 보면 답이 나온다. 우변에 $C_D/d$가 있다. 항력계수를 9% 깎는 동안 지름은 40 μm에서 24.5 μm로 39% 줄었다. 나눗셈 쪽이 훨씬 크게 움직인 것이다.

물리적으로는 이렇다. 항력은 면적($d^2$)에 비례하고 관성은 질량($d^3$)에 비례한다. 비율은 $1/d$다. 액적이 마를수록 같은 항력이 더 큰 감속을 만든다. 응답시간 $\tau_d \simeq \rho_l d^2 / 18\mu_g$가 $d^2$로 줄어든다는 말과 같다.

$\frac{a}{m} \propto \frac{d^2}{d^3} = \frac{1}{d}$

즉 기화는 항력계수를 낮추는 효과와 감속을 키우는 효과를 동시에 가진다. 앞의 것은 $-0.2$ 지수에 묶여 있고, 뒤의 것은 지름에 직접 걸려 있다. 조금만 오래 계산하면 뒤가 이긴다.

아래에서 네 액적을 같은 노즐에서 동시에 쏘아보자.

$B_M$을 올리면 C 레인이 회색 구형을 앞지른다. 그 상태에서 $K$를 올려보면 같은 계수를 쓰는 D 레인만 뒤처진다. $K = 0$으로 되돌리면 D가 C에 정확히 겹치므로, 두 효과 중 어느 쪽이 일하고 있는지 한 번에 분리된다.

어느 상관식을 언제 켜는가#

마지막 줄이 실제로 자주 문제가 된다. 표면장력이 0으로 가면 TAB의 복원항이 사라지고 $y$가 발산한다. 그 영역에서는 변형 보정을 끄고 확산 계면 쪽 모델로 넘어가야 한다.

$We$가 충분히 낮은데도 항력이 안 맞으면 보정이 아니라 다른 데를 봐야 한다. 벽면 근처의 항력 위기와 경계층 박리처럼 $Re$ 자체가 상관식의 유효 범위를 벗어난 경우다.

상관식 하나를 바꾸면 검증 케이스가 통째로 움직인다#

이 계산에서 관통거리는 -19.1%에서 +5.1% 사이를 오갔다. 격자도, 시간간격도, 난류 모델도 건드리지 않았다. 액적 하나에 곱해지는 무차원 인자 두 개를 켰다 껐다 한 것이 전부다.

그래서 스프레이 검증에서 실험과 어긋날 때 격자부터 조이는 것은 순서가 아니다. 먼저 항력계수 식이 무엇인지, 그 안에 $y$와 $B_M$이 어떤 값으로 들어가는지, $y$를 정상상태로 넣었는지 진동을 풀었는지를 확인한다. 지름이 시간에 따라 줄고 있다면 $C_D$보다 $C_D/d$를 출력해서 본다. 계수만 봐서는 왜 액적이 먼저 멈추는지 끝내 보이지 않는다.