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

코드의 동쪽 치우침은 속도 0.0891, 마하 약 0.154에 선다. N=2 다항식이 버티라고 적힌 격자 속도는 약 0.173까지다.

1986년에 격자 가스가 죽었다. 칸마다 입자가 0개 아니면 1개라, 노이즈가 평균으로도 안 지워졌다. 1992년에 그 정수를 실수 분포 $f_i$ 로 바꿨다. 노이즈는 사라졌다. 질문은 남았다. 실수 9개가 어떻게 무한차원 Maxwell–Boltzmann을 대신하는가.

그림의 drift $u$ 와 절단 차수 $N$ 을 바꾸면, 실선 Maxwell–Boltzmann과 점선 다항식이 포개지거나 벌어진다.

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이 저속(낮은 마하수)에 묶이는 이유.

$N=2$ 에 $u=0.4$ 이면 밀도·운동량·에너지가 포개진다. $u$ 가 1을 넘으면 점선만 실선에서 떨어진다. 마하 천장 $|\mathbf{u}|/c_s \lesssim 0.3$ 의 자리가 이 벌어짐 근처가 아니다. 숫자는 뒤에서 나눈다.

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

LGCA(Lattice Gas Cellular Automaton)는 정사각 격자, 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은 그 분포를 한 칸 민다. 9개로 충분한 이유는 이 $f_i^{eq}$ 의 모양에 있다.

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개 방향만 남기려면 이 분포를 추려야 한다.

가우시안 적분을 직사각형 합으로 바꾸면 실패한다. 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차원 평형을 등온, 곧 $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은 이 합을 유한 차수에서 끊는다.

등온 전개의 $RT$ 는 1이라 열속도도 1이다. 그림의 $u$ 는 그 열속도로 나눈 값이다. 천장으로 적힌 0.3은 이 눈금의 $|u|$ 다. $u=0.4$ 는 천장보다 0.1 위인데도, $N=2$ 의 첫 세 모멘트는 맞는다. 점선이 갈라지는 $|u|>1$ 은 마하 1을 지난 자리다. 0.3과 1은 같은 실패가 아니다.

급수의 $n=3$ 항은 $(u^3/6)\,He_3(c)$ 이고 $He_3=c^3-3c$ 다. 계수만 보면 $u=0.3$ 에서 $0.3^3/6=0.027/6=0.0045$ 다. $u=0.4$ 에서는 $0.4^3/6=0.064/6\approx 0.0107$ 이라 천장의 약 2.37배다. $u=1$ 에서는 $1/6\approx 0.1667$ 이라 $0.0045$ 의 약 37배다. $1^3/0.3^3=1/0.027$ 이 그 비다. $N=2$ 는 이 항을 통째로 버린다. 그래서 $u=0.4$ 에서 밀도·운동량·에너지가 맞아도, 버린 3차의 크기는 이미 천장 값의 두 배를 넘는다. 점선이 갈라지는 $u=1$ 은 그 계수가 0.17 근처까지 올라온 자리다. 마하 천장 0.3은 곡선이 눈에 띄게 떨어지는 지점이 아니다. 3차 계수가 0.0045를 넘는 순간을 미리 그어 둔 선이다.

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

1차 점성 Navier–Stokes에 필요한 것은 세 층이다. $\rho$ , $\rho\mathbf{u}$ , 운동량 플럭스 $\Pi_{\alpha\beta}=\rho u_\alpha u_\beta+p\delta_{\alpha\beta}$ . $\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]

이것이 평형 분포의 연속 형태다. $N=3$ 이상은 열유속에 쓰인다. 표준 등온 LBM에는 남지 않는다. 버린 항이 등온 가정과 마하 0.3 천장을 같이 만든다. 그 3차를 도로 넣으면 등온이라는 가정이 먼저 깨진다. 열유속은 그 항이 나르는 조각이라, 2차에서 멈춘 코드는 온도가 흐르는 문제를 이 다항식으로 받지 못한다. 마하 천장과 등온 가정은 한 자름의 양면이다. 속도를 0.173 아래에 두더라도, 뜨거운 벽과 차가운 벽을 한 격자에 올리면 이 평형은 열을 나르지 않는다.

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

연속 적분 $\int f^{eq}(c)\,\phi(c)\,dc$ 를 격자 합 $\sum_i w_i f_i^{eq}\phi(\mathbf{c}_i)$ 로 바꾸려면 구적(quadrature)이 필요하다. $N$ 차까지 맞추려면, 그 차수의 다항식을 정확히 적분하면 된다.

2차원에서 $u^2$ 까지 맞추려면 가우시안 가중 아래 5차 다항식까지 정확해야 한다. 운동량 텐서가 $c^2\times u^2$ 를 만들기 때문이다. Gauss–Hermite 3점 구적이 그 적분을 맡는다.

1차원 노드는 $c\in\{-\sqrt{3},0,+\sqrt{3}\}$ , 가중은 $\{1/6,2/3,1/6\}$ 이다. 그 세 근이 왜 그 자리인지는 6차 모멘트를 9개 속도로 접은 에르미트 구적이 따로 센다. 2차원 텐서 곱은 $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]

$1/c_s^2=3$ , $1/(2c_s^4)=9/2$ , $1/(2c_s^2)=3/2$ 다. 식의 3, $9/2$ , $3/2$ 는 이 세 개다. $\mathbf{u}=0$ 이면 대괄호는 1이고 $f_i^{eq}=w_i\rho$ 다. 가중을 더하면 $4/9+4\cdot(1/9)+4\cdot(1/36)=1$ 이다. 정지한 질량의 $4/9$ 는 0번 방향에 남는다. 표의 세 가중은 1차원 가중 $\{1/6,\ 2/3,\ 1/6\}$ 의 곱이다. 정지 방향은 $(2/3)\times(2/3)=4/9$ 다. 축 방향은 $(2/3)\times(1/6)=1/9$ 다. 대각은 $(1/6)\times(1/6)=1/36$ 이다. 축 4개와 대각 4개를 더하면 $4/9+4\cdot(1/9)+4\cdot(1/36)=1$ 로 다시 닫힌다. 9라는 개수와 세 개의 분모는 같은 텐서 곱에서 나온다.

이 한 줄이 등온 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}")

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

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

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

초기 분포는 $f=w$ 에서 $f_1$ 에 0.05를 더하고 $f_3$ 에서 0.04를 뺀다. 질량은 $1+0.05-0.04=1.01$ 이다. $x$ 운동량은 $1\cdot 0.05+(-1)\cdot(-0.04)=0.09$ 다. 그래서 $u_x=0.09/1.01\approx 0.0891$ 이고 $u_y=0$ 이다. $c_s=1/\sqrt{3}$ 으로 나누면 마하는 약 0.154다. 천장 0.3의 격자 속도 $0.3/\sqrt{3}\approx 0.173$ 보다 아래다. 본문이 말하는 20스텝 뒤 $10^{-6}$ 아래의 잔차는, 이 보존된 $(1.01,\ 0.0891)$ 위의 평형으로 붙는 과정이다. 충돌이 남기는 것은 매 스텝의 $\rho$ 와 $\rho\mathbf{u}$ 뿐이다.

$\tau$ 를 0.55 쪽으로 내리면 잔차가 한두 스텝 만에 사라진다. 점성은 $\nu=c_s^2(\tau-1/2)$ 다. 코드의 $\tau=0.8$ 이면 $\nu=(1/3)\cdot 0.3=0.1$ 이다. $\tau=0.55$ 이면 $\nu=(1/3)\cdot 0.05=1/60\approx 0.0167$ 이라, 코드 점성의 6분의 1이다. $\tau$ 를 $1/2$ 로 더 밀면 점성은 0이 되고 격자는 불안정해진다. $\tau$ 를 키우면 점성은 늘고, 평형으로 끌리는 속도는 느려진다. 비평형 한 스텝의 배율은 $1-1/\tau$ 다. 코드의 $\tau=0.8$ 이면 $1/\tau=1.25$ 라 배율은 $-0.25$ 다. 부호가 뒤집히는 과이완이다. 절댓값 $0.25$ 를 20번 곱하면 $0.25^{20}\approx 9.1\times 10^{-13}$ 이다. 편차가 1 정도여도 20스텝 뒤에는 $10^{-6}$ 보다 한참 아래다. 본문이 20스텝을 $10^{-6}$ 의 눈금으로 잡은 것은, 이 배율만으로도 이미 넉넉한 횟수다. $\tau=1$ 이면 배율은 0이라 한 스텝에 평형으로 붙는다. 코드의 0.8은 그 한 스텝 평형에서 한 칸 비껴 선 값이다.

마하 0.3 너머에서는 이 다항식을 믿지 않는다#

$N=2$ 는 $\rho$ , $\rho\mathbf{u}$ , 운동량 플럭스까지만 맞춘다. 그 위는 버렸다. 격자 속도가 약 0.173을 넘으면 버린 항이 더 이상 작지 않다. $\tau$ 를 내려 점성을 $0.1$ 에서 $1/60$ 로 줄여도 평형 다항식은 그대로다. 속도를 천장 위로 올릴지, 점성만 내릴지는 다른 실패다.

정지한 틀에서 맞춘 점성이 틀을 옆으로 밀어도 그대로인지는 이 식이 답하지 않는다. 잘려 나간 3차가 그 답을 쥐고 있다. 프레임을 옮겼을 때 3차 모멘트 결손이 점성을 바꾸는지는 아직 이 다항식 밖에 있다.