채프먼-엔스코그 전개

1. 개요

채프먼-엔스코그 전개(Chapman-Enskog expansion)는 분자 하나하나의 통계를 다루는 볼츠만 수송방정식으로부터, 크누센수를 작은 매개변수로 삼는 다중척도 점근 전개를 통해 오일러 방정식·나비에-스토크스 방정식 같은 연속체 방정식과 그 수송계수를 유도해 내는 방법이다. 1916~17년 시드니 채프먼과 다비트 엔스코그가 각각 독립적으로 완성했으며¹, “왜 기체 분자의 아수라장이 매끄러운 유체 방정식으로 보이는가”에 대한 인류의 표준 답변이다.

CFD 하는 사람에게 이 전개가 중요한 이유는 따로 있다. 격자 볼츠만 방법이 “가상 입자 놀이”를 하는데도 진짜 유체 해를 내놓는 근거가 바로 이것이고, LBM 사용자가 매일 손으로 만지는 완화시간 $\tau$ 와 점성 $\nu$ 의 관계식 역시 이 전개의 1차 항에서 떨어져 나온 물건이기 때문이다. 즉 LBM 코드에서 tau = 3*nu + 0.5 한 줄을 아무 생각 없이 치는 순간, 우리는 이미 채프먼-엔스코그의 결과를 쓰고 있는 셈이다.

2. 크누센수와 척도 분리

전개의 소전제는 하나다. 분자 하나가 충돌 사이에 날아가는 거리(평균자유행로 $\lambda$)가, 우리가 보고 싶은 유동 구조의 크기 $L$ 보다 압도적으로 작을 것.

$$\mathrm{Kn} = \frac{\lambda}{L} \ll 1$$

이 무차원수가 크누센수다. $\mathrm{Kn} \lesssim 10^{-2}$ 면 연속체 영역이라 나비에-스토크스가 잘 맞고, $10^{-2} \sim 10^{-1}$ 은 미끄럼 유동, 그 위는 천이·자유분자 영역이라 연속체 방정식이 무너진다. MEMS 채널, 진공 장비, 재진입 기체가 후자에 해당하며, 이때는 아예 볼츠만 쪽에서 직접 푸는 DSMC 같은 도구로 갈아탄다.

전개의 두 번째 아이디어가 시간 척도의 분리다. 충돌은 빠르고 유동은 느리다. 그래서 시간 미분을 척도별로 쪼갠다.

$$f = f^{(0)} + \mathrm{Kn}\, f^{(1)} + \mathrm{Kn}^2 f^{(2)} + \cdots, \qquad \partial_t = \partial_{t_0} + \mathrm{Kn}\, \partial_{t_1} + \cdots$$

3. 차수별로 무엇이 나오는가

무차원화한 볼츠만 방정식 $\partial_t f + \boldsymbol{\xi}\cdot\nabla f = \mathrm{Kn}^{-1} Q(f,f)$ 에 위 전개를 대입하고 크누센수 차수별로 정리하면, 유체역학 교과서의 계보가 순서대로 튀어나온다.

• $\mathrm{Kn}^{-1}$: $Q(f^{(0)}, f^{(0)}) = 0$. 충돌항의 영공간이라 해는 국소 맥스웰 분포 하나뿐이다. 여기서 밀도·속도·온도라는 다섯 개의 거시 변수가 등장한다.
• $\mathrm{Kn}^{0}$: 압력이 등방적인 오일러 방정식. 점성도 열전도도 없다. 즉 비점성 유동은 볼츠만의 0차 근사다.
• $\mathrm{Kn}^{1}$: 선형화된 충돌 연산자를 가해성 조건(질량·운동량·에너지 모멘트가 0) 아래 뒤집어 $f^{(1)}$ 을 구하면, 변형률에 비례하는 뉴턴 점성응력과 온도구배에 비례하는 푸리에 열유속이 유도된다. 나비에-스토크스-푸리에가 여기서 완성된다.
• $\mathrm{Kn}^{2}$ 이상: 버넷·초버넷 방정식. 미끄럼 유동에서 정확도가 조금 올라가지만, 짧은 파장에서 선형 불안정(보빌레프 불안정)이 알려져 있어 그대로 쓰면 폭발한다. 그래서 실무는 정규화 버넷이나 그라드 13-모멘트, R13 같은 모멘트법 계열로 우회한다.²

1차 결과의 위력은 수송계수가 공짜로 따라 나온다는 데 있다. 지름 $d$ 인 강체구 기체의 1차 근사 점성은

$$\mu = \frac{5}{16 d^2}\sqrt{\frac{m k_B T}{\pi}}$$

로, 밀도에 무관하고 $\sqrt{T}$ 로 증가한다. 압력을 두 배 올려도 점성이 그대로라는 맥스웰의 반직관적 예측이 여기서 나오고, 실험이 그걸 지지했다. 단원자 기체의 프란틀수가 $\mathrm{Pr} = 2/3$ 로 떨어지는 것도 같은 계산의 부산물이다.

4. LBM에서 ν = c_s²(τ − 1/2)Δt 가 나오는 자리

격자 볼츠만 방법의 BGK 갱신식은 완전히 이산적이다.

$$f_i(\mathbf{x} + \mathbf{e}_i \Delta t,\, t + \Delta t) - f_i(\mathbf{x}, t) = -\frac{1}{\hat\tau}\left(f_i - f_i^{\mathrm{eq}}\right)$$

좌변을 $\Delta t$ 로 테일러 전개하면 1차 항 $\Delta t(\partial_t + \mathbf{e}_i\cdot\nabla)f_i$ 뒤에 2차 항 $\tfrac{\Delta t^2}{2}(\partial_t + \mathbf{e}_i\cdot\nabla)^2 f_i$ 가 남는다. 여기에 위와 똑같이 $f_i = f_i^{\mathrm{eq}} + \varepsilon f_i^{(1)} + \varepsilon^2 f_i^{(2)}$, $\partial_t = \varepsilon\partial_{t_1} + \varepsilon^2\partial_{t_2}$ 를 넣고 차수를 맞추면, $\varepsilon^2$ 방정식의 유속 항 앞에 인자 $\left(1 - \frac{1}{2\hat\tau}\right)$ 가 붙는다. 이 인자를 2차 모멘트로 옮기면 그대로

$$\nu = c_s^2\left(\hat\tau - \frac{1}{2}\right)\Delta t, \qquad c_s^2 = \frac{1}{3}\frac{\Delta x^2}{\Delta t^2}\ (\text{D2Q9})$$

가 된다. 그 악명 높은 $-1/2$ 는 물리가 아니라 이산 격자의 2차 테일러 잔여항이다. 이걸 흡수해 준 덕에 LBM은 1차가 아니라 2차 정확도를 갖고, 동시에 $\hat\tau > 1/2$ 라는 하한이 강제된다. 레이놀즈수를 올릴수록 $\hat\tau \to 1/2^{+}$ 로 붙으면서 유효 점성이 0에 수렴하고, 그 지점에서 SRT-LBM이 불안정해지는 것도 같은 식이 말해 주는 사실이다.³ SRT에서는 체적점성까지 같은 $\hat\tau$ 하나에 묶여 따로 못 만지는데, MRT가 완화율을 모멘트별로 쪼개 이 족쇄를 푸는 것도 채프먼-엔스코그로 각 모멘트의 역할을 확인했기에 가능한 설계다.

5. 격자 기체 오토마타와 갈릴레이 불변성

이 전개가 진짜 무기로 쓰인 첫 무대는 LBM이 아니라 격자 기체 오토마타였다. FHP류 LGCA는 점유수 $n_i \in \{0,1\}$ 만 다루는 셀룰러 오토마타인데, 앙상블 평균 후 채프먼-엔스코그를 돌려 보면 이류항 앞에 밉살스러운 인자가 하나 붙는다.

$$\partial_t \mathbf{u} + g(\rho)\,(\mathbf{u}\cdot\nabla)\mathbf{u} = -\nabla P + \nu\nabla^2\mathbf{u}, \qquad g(\rho) = \frac{D}{D+2}\cdot\frac{1-2d}{1-d}$$

여기서 $d$ 는 방향당 점유확률이다. $g \neq 1$ 이라는 것은 갈릴레이 불변성이 깨졌다는 뜻이다. 등속으로 움직이는 좌표계에서 물리가 달라진다는, 유체 코드가 절대 하면 안 되는 짓을 격자가 하고 있는 것이다. LBM이 페르미-디랙형 점유수를 실수 분포함수로 바꾸고 2차 다항식 평형분포를 도입한 결정적 이유가 이 $g(\rho)$ 를 1로 되돌리기 위해서였다.

그런데 완전히 고쳐진 것은 아니다. D2Q9 격자는 4차 모멘트 텐서까지는 등방적이지만 3차 모멘트 $\sum_i w_i e_{i\alpha}e_{i\beta}e_{i\gamma}$ 를 연속체 값과 일치시키지 못한다. 그 결과 운동량 방정식에 $O(u^3)$ 짜리 오차항(입방 속도 결함)이 남고, 이것이 표준 LBM에 저마하수 제약이 붙는 근본 원인이다. 마하수를 0.1~0.2 아래로 유지하라는 국룰은 경험칙이 아니라 이 전개가 계산해 준 오차 항의 크기다.⁴

6. 한계와 오해

• 점근 급수이지 수렴 급수가 아니다. 차수를 올린다고 정확해지는 게 아니라, 어느 지점부터 나빠진다. 버넷이 불안정한 것이 그 대표적 증상이다.
• 가해성 조건이 본질이다. 각 차수에서 $f^{(n)}$ 의 보존 모멘트를 0으로 두는 선택이 없으면 해가 유일하지 않다. 이 선택이 “밀도·속도·온도는 전부 $f^{(0)}$ 이 정의한다”는 물리적 규약이다.
• 강한 비평형에는 안 통한다. 충격파 내부, 재진입 유동, 진공 장비처럼 국소 $\mathrm{Kn}$ 이 커지는 곳에서는 1차 근사가 무너진다. 압축성 유동 코드가 충격파 두께를 물리적으로 재현하지 못하는 이유이기도 하다.
• LBM은 볼츠만 방정식을 푸는 게 아니다. 정확히는 “채프먼-엔스코그를 거쳐 나비에-스토크스와 점근적으로 일치하도록 설계된 이산 동역학”을 푸는 것이다. 미시 물리가 진짜인 척하면 곤란하다.

7. 관련 문서

• 격자 볼츠만 방법 · 격자 기체 오토마타
• 나비에-스토크스 방정식 · 오일러 방정식
• 점성 · 무차원수 · 크누센수
• 기체동역학 · 압축성 유동
• 모멘트법 · 셀룰러 오토마타
• 전산유체역학 · 수치해석

8. Footnotes

1. 두 사람은 서로의 작업을 모른 채 거의 동시에 도달했고, 방법론도 달랐다. 채프먼은 맥스웰식 수송방정식을, 엔스코그는 볼츠만 방정식의 체계적 섭동 해법을 썼다. 지금 교과서에 실린 정돈된 형태는 엔스코그 쪽에 가깝다. ↩

2. 버넷 방정식은 “차수를 올렸으니 더 정확하겠지”라는 순진한 기대를 배신하는 대표 사례다. 짧은 파장 교란에 대해 선형 불안정이라, 격자를 촘촘히 할수록 더 빨리 터진다. 격자 수렴 연구를 돌리다 발산하는 걸 보고 코드를 의심하면 안 되고 방정식을 의심해야 하는, 흔치 않은 상황. ↩

3. 그래서 LBM 튜닝의 국룰은 “레이놀즈수를 올리고 싶으면 $\tau$ 를 1/2로 밀지 말고 격자를 늘려라”이다. 물론 격자를 늘리면 시간 스텝도 같이 줄어들어 계산량이 순식간에 갈려 나간다. 세상에 공짜는 없다. ↩

4. 3차 모멘트 결함을 없애려면 속도 방향을 더 준 격자(D2Q17, D2Q37 등)를 써야 한다. 열유동·고마하수용 고차 격자가 존재하는 이유가 이것인데, 방향 수가 늘면 메모리와 대역폭이 그대로 늘어난다. 결국 “어느 모멘트까지 살릴 것인가”가 격자 설계의 전부다. ↩