9개 화살표 안의 Maxwell–Boltzmann — LBM 평형 분포가 다항식인 이유

에르미트 다항식 2차 절단이 D2Q9를 만들고, 같은 절단이 LBM을 저속 흐름에 묶는다

1986년, Frisch·Hasslacher·Pomeau가 격자 위에 입자를 풀어 Navier–Stokes를 흉내 냈다. LGCA(Lattice Gas Cellular Automaton)다. 5년이 지나기 전에 죽었다 — 격자 한 칸에 입자가 0개 아니면 1개라는 정수 조건이 도저히 지울 수 없는 노이즈를 만들어냈기 때문이다. 1992년 McNamara와 Zanetti는 정수를 실수로 바꾼다. 입자 분포 함수 $f_i$ 로 격자를 채운다. 노이즈가 사라졌다 — 그 대신 새로운 질문이 따라왔다. 이 9개의 실수가 어떻게 Maxwell–Boltzmann이 들고 있던 무한차원 분포를 흉내 낼 수 있는가. 이 글은 그 답을 따라간다 — 에르미트 다항식 2차 절단이 LBM의 평형 분포를 만들고, 같은 절단이 LBM을 저속 흐름에 묶는다는 사실까지.

정수에서 실수로 — 1986년 LGCA의 죽음과 LBM의 탄생#

LGCA의 아이디어는 단순하고 우아했다. 정사각 격자, 6방향 입자, 충돌 규칙 하나. 거시적으로 평균하면 Navier–Stokes가 살아 돌아온다. 그러나 부울 입자(0 또는 1)는 statistical noise(통계적 잡음)를 키웠다. 노이즈를 평균으로 죽이려면 시간 평균을 길게 잡아야 했고, 그러면 본래의 transient 흐름이 사라졌다.

McNamara·Zanetti는 정수 입자 대신 실수 분포 $f_i(\mathbf{x}, t)$ 를 격자 노드마다 놓는다. $f_i$ 는 시각 $t$ , 위치 $\mathbf{x}$ 에서 격자 방향 $\mathbf{c}_i$ 를 향해 움직이는 입자의 기대 개수다. 노이즈가 본질적으로 사라진다. 한 스텝의 진행은 두 단계로 끝난다 — collision(충돌)과 streaming(전송).

f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) = f_i(\mathbf{x}, t) - \frac{f_i(\mathbf{x}, t) - f_i^{eq}(\rho, \mathbf{u})}{\tau}

$\tau$ 는 relaxation time(이완 시간), $f_i^{eq}$ 가 평형 분포 함수다. 충돌 항이 비평형 부분을 $\tau$ 시간 척도로 깎고, streaming은 그렇게 깎인 분포를 격자 한 칸 옆으로 넘긴다. 그런데 $f_i^{eq}$ 가 무엇이길래 9개 실수만으로 충분한가.

Maxwell–Boltzmann이 들고 있는 폭탄#

연속 운동론에서 평형 분포는 Maxwell–Boltzmann이다.

f^{MB}(\mathbf{c}) = \rho \, (2\pi R T)^{-D/2} \exp\!\left( -\frac{(\mathbf{c} - \mathbf{u})^2}{2 R T} \right)

$\mathbf{c}$ 는 입자 분자 속도, $\mathbf{u}$ 는 거시 유속, $D$ 는 공간 차원, $T$ 는 온도, $R$ 은 기체상수다. 이 식이 들고 있는 폭탄은 $\mathbf{c}$ 가 연속이라는 것이다. 격자 위에서 9개나 19개의 이산 방향만 남기려면 어떤 방식으로든 이 연속 분포를 추려야 한다.

직관적 시도 1 — 가우시안 적분을 그냥 직사각형 합으로 근사한다. 실패한다. Navier–Stokes에 필요한 모멘트 — $\rho$ , $\rho \mathbf{u}$ , 운동량 텐서 $\rho \mathbf{u}\mathbf{u}$ , 에너지 — 가 격자 합과 어긋난다. 보존이 깨지면 LBM은 그 자리에서 무용지물이 된다.

올바른 답은 함수 전개다. Maxwell–Boltzmann을 직교 다항식 급수로 펼쳐 놓고, 필요한 모멘트까지만 살리고 나머지를 자른다. 그 직교 다항식이 Hermite(에르미트)다.

에르미트 다항식이 분포를 잘게 썬다#

probabilist Hermite polynomial(확률론적 에르미트 다항식) $He_n(c)$ 는 표준 가우시안 가중 $w(c) = (2\pi)^{-1/2} e^{-c^2/2}$ 아래에서 직교한다.

\int He_n(c)\, He_m(c)\, w(c)\, dc = n!\, \delta_{nm}

처음 몇 개는 $He_0 = 1$ , $He_1 = c$ , $He_2 = c^2 - 1$ , $He_3 = c^3 - 3c$ 다. 1차원 Maxwell–Boltzmann을 등온( $RT = 1$ 로 무차원화)으로 놓으면, 이항생성함수의 도움으로 다음 항등식이 곧장 따라온다.

f^{MB}(c) = w(c)\, \rho \exp\!\left( u c - \frac{u^2}{2} \right) = w(c)\, \rho \sum_{n=0}^{\infty} \frac{u^n}{n!} He_n(c)

각 항이 우주의 한 조각씩을 잡는다. $n=0$ 은 정지 가우시안. $n=1$ 은 작은 drift. $n=2$ 는 평균 운동에너지 보정. 더 높은 $n$ 은 점점 더 비대칭적인 꼬리 효과다. LBM이 노리는 것은 이 합을 유한 차수에서 끊는 일이다.

Maxwell–Boltzmann과 그 Hermite 절단의 비교. 표의 세 모멘트(밀도·운동량·에너지)가 같아지는 차수가 바로 LBM이 멈추는 자리다.

drift velocity u = 0.40 (= |u|/c_s)truncation order N = 2

N=2면 ρ·ρu·에너지 세 모멘트가 정확히 일치한다. u가 |u|/c_s ≈ 1을 넘어가면 절단오차가 폭주하며 격자가 깨진다 — LBM이 저속(낮은 마하수)에 묶이는 이유.

아래 시뮬레이션에서 직접 조작해보자. drift $u$ 를 0.4까지 올린 채 절단 차수 $N=2$ 를 고르면 첫 세 모멘트(밀도·운동량·에너지)가 Maxwell–Boltzmann과 정확히 일치한다. $u$ 를 1.0 너머로 밀면 점선이 곧장 떨어져 나간다. 이것이 LBM의 마하수 제약 $|\mathbf{u}|/c_s \lesssim 0.3$ 이 어디서 오는지 한 장에 보여준다.

2차 절단의 의미 — 보존되는 모멘트#

Navier–Stokes의 1차 점성 응력까지 살리려면 어떤 모멘트가 필요한가. 정확한 답이 있다 — $\rho$ , $\rho \mathbf{u}$ , 그리고 운동량 플럭스 $\Pi_{\alpha\beta} = \rho u_\alpha u_\beta + p \delta_{\alpha\beta}$ 까지 총 3개의 모멘트 차수다. $\rho \mathbf{u}\mathbf{u}$ 가 $u$ 의 2차이므로, 분포 함수에 $u^2$ 항이 들어가야 한다.

따라서 에르미트 급수를 $N=2$ 에서 자른다.

f^{eq}(c) = w(c)\, \rho \left[ 1 + uc + \frac{u^2}{2}(c^2 - 1) \right]

이것이 LBM 평형 분포의 연속 형태다. $N=3$ 이상은 비점성 흐름의 고차 모멘트(열전도)에는 도움이 되지만 표준 isothermal LBM에는 잉여다. 잘려나간 항들이 LBM의 두 가지 한계를 만든다 — 등온 가정과 마하수 0.3 천장.

D2Q9는 어떻게 만들어졌는가#

연속 적분 $\int f^{eq}(c) \phi(c)\, dc$ 를 격자 합 $\sum_i w_i f_i^{eq} \phi(\mathbf{c}_i)$ 로 옮기려면 quadrature(직교 적분 공식)가 필요하다. 핵심은 — $N$ 차 절단까지 정확한 quadrature는 그 차수의 다항식을 정확히 적분하면 충분하다는 점이다.

2차원에서 $u^2$ 항까지 정확하려면 가우시안 가중에 대해 5차 다항식까지 정확한 quadrature가 필요하다 (운동량 텐서가 $c^2 \times u^2$ 항을 만들기 때문). Gauss–Hermite 3점 quadrature가 정확히 그 일을 한다.

1D에서 점 3개 $c \in \{-\sqrt{3}, 0, +\sqrt{3}\}$ , 가중 $\{1/6, 2/3, 1/6\}$ . 2D에서 텐서 곱으로 $3 \times 3 = 9$ 점이 나온다 — 이것이 D2Q9다. 격자 단위로 다시 스케일하면 $c_s = 1/\sqrt{3}$ 이고 격자 속도 벡터는 정수 좌표 $(0,0), (\pm 1, 0), (0, \pm 1), (\pm 1, \pm 1)$ 이 된다. 격자 가중 $w_i$ 는 Gauss–Hermite 가중에 가우시안 인자를 흡수시킨 값이다.

인덱스	$\mathbf{c}_i$	$w_i$
0	$(0, 0)$	$4/9$
1–4	$(\pm 1, 0), (0, \pm 1)$	$1/9$
5–8	$(\pm 1, \pm 1)$	$1/36$

격자 위 평형 분포는 연속 형태에 $c_s^2 = 1/3$ 을 대입한 것이다.

f_i^{eq} = w_i\, \rho \left[ 1 + 3 (\mathbf{c}_i \cdot \mathbf{u}) + \frac{9}{2}(\mathbf{c}_i \cdot \mathbf{u})^2 - \frac{3}{2}\, \mathbf{u} \cdot \mathbf{u} \right]

이 한 줄이 모든 isothermal LBM 코드의 심장이다. 단순 다항식이지만 그 안에 Maxwell–Boltzmann의 첫 세 모멘트가 봉인되어 있다.

Python — BGK 충돌 한 스텝#

격자 한 칸에서 충돌 한 스텝을 돌려본다.

import numpy as np
 
cx = np.array([0, 1, 0, -1, 0, 1, -1, -1, 1])
cy = np.array([0, 0, 1, 0, -1, 1, 1, -1, -1])
w  = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
 
def equilibrium(rho, ux, uy):
    """D2Q9 평형 분포 — 9개 다항식"""
    ciu = cx * ux + cy * uy            # 점곱
    usq = ux * ux + uy * uy
    return w * rho * (1 + 3 * ciu + 4.5 * ciu**2 - 1.5 * usq)
 
def bgk_collision(f, tau):
    rho = f.sum()
    ux  = (cx * f).sum() / rho
    uy  = (cy * f).sum() / rho
    feq = equilibrium(rho, ux, uy)
    return f - (f - feq) / tau, rho, ux, uy
 
# 비평형 초기 분포
f = w * 1.0
f[1] += 0.05; f[3] -= 0.04          # 동쪽으로 살짝 치우침
 
for step in range(20):
    f, rho, ux, uy = bgk_collision(f, tau=0.8)
    err = np.linalg.norm(f - equilibrium(rho, ux, uy))
    print(f"step {step:2d}: ρ={rho:.4f}, u=({ux:+.4f},{uy:+.4f}), ||Δ||={err:.2e}")

20 스텝쯤 돌리면 잔차가 $10^{-6}$ 아래로 떨어지면서 9개 $f_i$ 가 평형 다항식에 맞춰 자리를 잡는다. 보존되는 것은 무엇인가 — 매 스텝마다 $\rho$ 와 $\rho \mathbf{u}$ 다. 충돌은 비평형 모멘트만 잘라낸다.

D2Q9 격자 위에서 9개 분포 함수가 BGK 충돌을 통해 평형으로 끌려가는 과정. 노란 화살표는 거시 속도 u, 가는 화살표는 각 격자 방향 c_i의 분포 크기.

u_x = 0.12u_y = 0.06τ (이완 시간) = 0.80

τ→0.5에 가까울수록 한 스텝에 평형으로 끌려가지만, 점성 ν = c_s²(τ−1/2)도 작아져 안정성이 흔들린다.

위 시뮬레이션에서 $\tau$ 슬라이더를 0.55 가까이 끌어내려 보면, 한두 스텝 만에 잔차가 0에 닿는다. 점성은 $\nu = c_s^2 (\tau - 1/2)$ 인데, $\tau \to 1/2$ 로 가면 점성이 0이 되며 격자가 불안정해진다. 반대로 $\tau$ 를 키우면 점성은 커지지만 평형으로 끌려가는 속도가 느려진다. 이 trade-off가 LBM 솔버 튜닝의 핵심이다.

핵심 3줄 요약#

LBM 평형 분포 $f_i^{eq}$ 는 Maxwell–Boltzmann을 에르미트 다항식 급수로 펼친 뒤 2차에서 자른 형태다 — 그 절단이 $\rho$ · $\rho \mathbf{u}$ ·운동량 플럭스 세 모멘트를 정확히 보존한다.
D2Q9 격자는 2차 절단 평형 분포를 정확히 적분하는 Gauss–Hermite 3점 quadrature의 2D 텐서 곱으로 유도된다. 9는 자연수가 아니라 직교 적분의 산물이다.
절단 차수가 평형 분포의 한계를 정한다 — isothermal 가정, 마하수 $\lesssim 0.3$ 천장은 모두 " $N=2$ 에서 잘랐다"는 한 줄에서 자동으로 따라온다.