엔트로피가 폭주를 막는다 — Karlin의 엔트로픽 격자 볼츠만

엔트로피는 보통 "잃어버린 정보"의 이름표로 배운다. 무질서가 늘어나는 방향, 되돌릴 수 없는 손실의 척도다. 그런데 격자 볼츠만(Lattice Boltzmann Method, LBM)에서는 바로 그 엔트로피가 폭주하는 계산을 붙잡아 세운다. 이 글은 Karlin 계열의 엔트로픽 LBM(Entropic LBM, ELBM)이 어떻게 H-정리를 이산화 단계에서 강제해 저점성 영역의 발산을 막는지 다룬다. 읽고 나면 충돌 연산자 한 줄에 왜 매 스텝 비선형 방정식을 푸는 비용을 감수하는지 납득하게 된다.

τ가 0.5에 다가가면 BGK가 무너진다#

표준 LBM은 BGK 충돌(단일 완화시간 근사)을 쓴다. 분포 $f_i$ 를 평형 $f_i^{eq}$ 쪽으로 $\omega$ 만큼 끌어당긴다.

f_i^* = f_i + \omega\,(f_i^{eq} - f_i)

여기서 $f_i$ 는 격자 방향 $i$ 의 분포, $f_i^{eq}$ 는 평형 분포, $\omega = 1/\tau$ 는 완화율이다.

점성은 $\nu = c_s^2(\tau - 1/2)$ 로 정해진다. $c_s^2 = 1/3$ 는 격자 음속의 제곱이다. 고Reynolds(저점성) 유동을 풀려면 $\tau \to 1/2$ , 즉 $\omega \to 2$ 로 밀어야 한다.

문제는 $\omega$ 가 2에 가까워질수록 충돌이 평형을 지나쳐 반대편으로 튄다는 점이다. 분포가 음수로 내려가고, 한 번 음수가 된 $f_i$ 는 다음 스텝에서 더 크게 진동한다. 격자가 터진다. BGK에는 이 과조화(over-relaxation)를 막을 장치가 없다.

볼츠만의 H-정리 — 시간의 화살#

볼츠만은 1872년에 다음 양이 충돌만으로는 절대 늘지 않음을 증명했다.

H = \int f \ln f \, d\mathbf{c}

$f$ 는 속도 분포, $\mathbf{c}$ 는 분자 속도다. $dH/dt \le 0$ 이며, 등호는 평형(Maxwell-Boltzmann 분포)에서만 성립한다.

H는 음의 엔트로피다. H가 단조 감소한다는 것은 엔트로피가 단조 증가한다는 말과 같다. 이것이 비가역성, 곧 시간의 화살이다. 핵심은 이 부등식이 안정성 보증서라는 사실이다. H가 절대 늘지 않는 수치 기법은 발산할 수 없다. 진폭이 무한정 커지려면 H가 늘어나야 하기 때문이다.

이산 H-함수와 엔트로픽 평형#

LBM은 연속 속도를 9개(D2Q9) 방향으로 잘라낸다. H도 그에 맞춰 이산화한다.

H(f) = \sum_i f_i \ln\frac{f_i}{w_i}

$w_i$ 는 격자 가중치(D2Q9에서 중심 $4/9$ , 면 $1/9$ , 대각 $1/36$ )다. 가중치로 나눈 비율의 로그를 취한 형태다.

엔트로픽 평형은 이 $H$ 를 보존량 제약 아래 최소화하는 분포로 정의한다. 질량 $\sum_i f_i = \rho$ 와 운동량 $\sum_i \mathbf{c}_i f_i = \rho\mathbf{u}$ 를 고정한 채 H를 최소화한다(라그랑주 승수 문제). D2Q9에서는 닫힌 곱(product) 형태가 존재한다.

f_i^{eq} = \rho\, w_i \prod_{\alpha=x,y}\left(2 - \sqrt{1+3u_\alpha^2}\right)\left(\frac{2u_\alpha + \sqrt{1+3u_\alpha^2}}{1 - u_\alpha}\right)^{c_{i\alpha}}

$u_\alpha$ 는 속도 성분, $c_{i\alpha}$ 는 격자 속도 성분(−1, 0, 1)이다. 표준 LBM이 쓰는 2차 절단 다항식 평형과 달리, 이 평형은 H의 진짜 최소점이다.

엔트로픽 안정자 — α를 매 스텝 푼다#

ELBM의 충돌은 BGK와 모양은 같지만 완화율을 둘로 쪼갠다.

f_i^* = f_i + \alpha\beta\,(f_i^{eq} - f_i)

$\beta \in (0,1)$ 는 점성을 정하고( $\nu = c_s^2(\tfrac{1}{2\beta} - \tfrac{1}{2})$ , $\beta\to1$ 이면 무점성), $\alpha$ 는 안정성을 정한다.

$\alpha$ 는 고정값이 아니다. 매 스텝, 매 격자점에서 다음 방정식의 비자명 근으로 푼다.

H\big(f + \alpha(f^{eq} - f)\big) = H(f)

이 식의 뜻은 명확하다. 충돌 후 H를 충돌 전과 정확히 같게 만드는 만큼만 평형을 향해 튄다. H가 늘어날 일이 원천 봉쇄된다. $\alpha=0$ 은 자명한 근이고, 우리가 찾는 것은 평형 너머의 두 번째 근이다. 평형 근처에서는 $\alpha \approx 2$ 로, BGK와 똑같아진다. 비평형이 클 때만 $\alpha$ 가 2에서 벗어나 과조화를 깎아낸다.

Python — 엔트로픽 α를 이분법으로#

곱 형태 평형과 H-함수, 그리고 이분법으로 $\alpha$ 를 푸는 한 스텝이다.

import numpy as np
 
# D2Q9 격자 가중치와 속도 (lattice weights & velocities)
W  = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
CX = np.array([0, 1, 0, -1,  0, 1, -1, -1,  1])
CY = np.array([0, 0, 1,  0, -1, 1,  1, -1, -1])
 
def h_function(f):
    # 이산 H-함수: H = Σ f_i ln(f_i / w_i)  (볼츠만 H-정리의 격자 버전)
    return np.sum(f * np.log(f / W))
 
def entropic_equilibrium(rho, ux, uy):
    # Ansumali–Karlin 곱(product) 형태의 정확한 엔트로픽 평형
    def factor(u):
        s = np.sqrt(1 + 3 * u * u)
        return (2 - s), (2 * u + s) / (1 - u)
    bx, rx = factor(ux)
    by, ry = factor(uy)
    return rho * W * bx * by * (rx ** CX) * (ry ** CY)
 
def entropic_alpha(f, feq):
    # H(f + α(feq − f)) = H(f) 의 비자명 근 α* 를 이분법으로 찾는다
    d = feq - f
    H0 = h_function(f)
    G = lambda a: h_function(f + a * d) - H0
    amax = np.min(np.where(d < 0, -f / d, np.inf))   # 분포 양수 보장 상한
    hi = min(0.999 * amax, 4.0) if np.isfinite(amax) else 4.0
    if G(hi) <= 0:
        return 2.0                                   # near-equilibrium 안전값
    lo = 1.0
    while G(lo) > 0 and lo > 1e-3:
        lo *= 0.5
    for _ in range(60):
        mid = 0.5 * (lo + hi)
        if G(mid) > 0: hi = mid
        else:          lo = mid
    return 0.5 * (lo + hi)
 
# --- 한 노드, 한 스텝 ---
feq = entropic_equilibrium(1.0, 0.05, 0.0)
f = feq.copy()
f[1] += 0.03; f[3] -= 0.03            # streaming 이 실어온 비평형 (non-equilibrium)
rho = f.sum(); ux = (CX * f).sum() / rho; uy = (CY * f).sum() / rho
feq = entropic_equilibrium(rho, ux, uy)
 
alpha = entropic_alpha(f, feq)
beta  = 0.95                          # ν = c_s²(1/(2β) − 1/2), β→1 이면 저점성
f_post = f + alpha * beta * (feq - f)
 
print(f"엔트로픽 α  = {alpha:.4f}   (BGK 고정값 2.0)")
print(f"H 충돌 전   = {h_function(f):.8f}")
print(f"H 충돌 후   = {h_function(f_post):.8f}")
print(f"최소 분포   = {f_post.min():.3e}  (> 0)")

출력에서 $\alpha$ 가 2에 가깝지만 정확히 2는 아님을 확인할 수 있다. 그 미세한 차이가 H를 보존하는 보정량이다.

인터랙티브 — G(α)의 두 번째 뿌리#

아래 시뮬레이션에서 직접 조작해보자. 비평형 섭동을 키우면 곡선이 어떻게 휘는지 본다.

Non-equilibrium perturbation: 0.060

The yellow dot is the non-trivial root α* where the discrete H-function returns to its pre-collision value. Near equilibrium α* ≈ 2, which is exactly standard BGK. Push the perturbation up and α* drifts away from 2 — that gap is the entropic correction BGK ignores.

노란 점이 우리가 찾는 $\alpha^*$ 다. 섭동이 작을 때는 점선( $\alpha=2$ ) 바로 위에 붙어 있다. 슬라이더를 오른쪽으로 밀면 $\alpha^*$ 가 2에서 떨어진다. BGK는 이 이동을 무시하고 늘 2를 쓴다. 그 차이가 쌓이면 발산한다.

인터랙티브 — H가 늘어나는 순간을 막다#

같은 비평형 상태를 두 충돌 규칙으로 처리한 결과를 나란히 본다. 빨간 막대는 음수가 된 분포다.

Non-equilibrium shear: 0.070

Both panels show post-collision populations. The entropic collision lands exactly on the equal-entropy point, so H_ent = H_pre. Fixed-α BGK misses it; once the non-equilibrium part grows, BGK over-shoots and H actually increases — a violated discrete H-theorem, which is the seed of low-viscosity blow-up.

엔트로픽 쪽은 항상 $H_{ent} = H_{pre}$ 다. 설계가 그렇게 되어 있다. BGK 쪽은 비평형이 커지면 $H_{BGK}$ 가 충돌 전보다 커지는 순간이 온다. 그때 빨간 경고가 뜬다. 바로 그 H 증가가 저점성에서 격자를 터뜨리는 씨앗이다.

구현하며 부딪힌 세 가지#

첫째, 매 격자점 비선형 풀이는 비싸다. 그래서 실전에서는 평형 근처( $\|f - f^{eq}\|$ 가 작을 때)는 $\alpha \approx 2$ 로 두고, 비평형이 임계값을 넘는 격자에서만 이분법을 켠다. 비용이 1할 수준으로 떨어진다.

둘째, 곱 형태 평형은 $|u| < 1$ 에서만 정의된다. 분모에 $1 - u$ 가 있다. 무차원 속도가 0.3을 넘으면 압력 텐서가 비등방적으로 휘기 시작하므로, ELBM도 $|u| \lesssim 0.1$ 영역에서 쓰는 것이 안전하다.

셋째, 2차 절단 다항식 평형을 그대로 쓰면 H 최소화 보증이 깨진다. 엔트로픽 안정자의 이점을 온전히 누리려면 평형도 엔트로픽 곱 형태여야 한다. 다항식 평형에 $\alpha$ 만 얹는 것은 절반의 효과다.

핵심 3줄 요약#

BGK가 저점성에서 터지는 이유는 과조화가 이산 H를 늘리기 때문이다. ELBM은 충돌 후 H를 충돌 전과 같게 만드는 $\alpha$ 를 매 스텝 풀어 이를 막는다.
$\alpha$ 는 $H(f + \alpha(f^{eq}-f)) = H(f)$ 의 비자명 근이며, 평형 근처에서 $\alpha \approx 2$ 로 BGK를 자동 복원한다.
점성은 $\beta$ 가, 안정성은 $\alpha$ 가 분담한다. 비용은 매 스텝 비선형 풀이지만, 평형 근처 단축으로 실전에서는 감당할 만하다.