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

매 노드에서 완화율을 2로 두지 않고, 충돌 뒤 이산 H가 그대로인 두 번째 근만 남긴다.

엔트로피를 잃어버린 정보의 이름표로만 외우면, 이 계산이 왜 안 터지는지 설명이 안 된다. 격자 볼츠만(Lattice Boltzmann Method, LBM)에서 점성을 낮추려고 $\tau$ 를 $1/2$ 에 붙이면 BGK가 평형을 지나친다. 분포가 음수로 떨어진다. Karlin 계열 엔트로픽 LBM(Entropic LBM, ELBM)은 그 자리에서 이산 H를 충돌 전과 같은 값에 묶는다. 대가는 노드마다 비선형 근 하나다.

비평형 섭동 슬라이더를 오른쪽으로 밀면 노란 점이 점선 $\alpha=2$ 에서 떨어지고, $G(\alpha)=0$ 의 두 번째 근이 왼쪽으로 움직인다.

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.

섭동이 작으면 점이 2 위에 붙어 있다. 점선과 노란 점 사이가, BGK가 무시하는 과조화의 폭이다.

τ가 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차 절단은 Hermite를 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에서 벗어나 과조화를 깎는다.

점성 식을 숫자로 펼치면, 과조화가 벽에 얼마나 붙어 있는지 보인다. $c_s^2 = 1/3$ 이고 $\nu = c_s^2(\tau - 1/2)$ 다. $\tau = 0.51$ 이면 차이 $0.01$ 이라 $\nu = 1/300 \approx 0.00333$ 이다. 완화율은 $\omega = 1/0.51 \approx 1.961$ 이다. 평형에서 벗어난 길이 1을 1.961만큼 가므로, 평형을 0.961만큼 지나 반대편에 선다. $\tau = 0.501$ 로 더 붙이면 $\omega \approx 1.996$ 이다. 거울상 $\omega = 2$ 까지 남은 간격은 0.004뿐이다.

코드에 적힌 $\beta = 0.95$ 도 같은 식이다. $\tfrac{1}{2\beta}-\tfrac{1}{2} = \tfrac{1}{1.9}-\tfrac{1}{2} = \tfrac{1}{38}$ 이라 $\nu = 1/114 \approx 0.00877$ 이다. 같은 점성의 BGK라면 $\tau = 1/(2\beta) \approx 0.526$ 이다. 평형 근처에서 $\alpha \approx 2$ 이면 유효 완화율 $\alpha\beta = 1.9$ 라 이 BGK와 점성이 같다. 비평형이 커져 $\alpha$ 가 2 아래로 내려가면 $\alpha\beta$ 도 1.9보다 작아진다. 평형을 지나치는 길이 $\alpha\beta-1$ 이 0.9보다 짧아진다. 점성은 $\beta$ 가 정한 그대로고, 잘리는 것은 과조화뿐이다.

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를 보존하는 보정량이다.

섭동 두 줄은 모멘트로 읽힌다. $f_1$ 에 $0.03$ 을 더하고 $f_3$ 에서 $0.03$ 을 빼면 질량 합은 그대로다. 방향 1과 3의 격자 속도가 $+1$ , $-1$ 이라 $x$ 운동량은 $0.03-(-0.03)=0.06$ 만큼 는다. 처음 $\rho=1$ , $u_x=0.05$ 였으니 섭동 직후 $u_x=0.11$ 이다. 곱 형태 인자는 $(2u+\sqrt{1+3u^2})/(1-u)$ 다. $u=0.05$ 에서 약 1.162, $u=0.11$ 에서 약 1.391이다. 무차원 속도 0.3에서는 약 2.467까지 올라가고, 분모 $1-u$ 는 0.7이다. $|u| \lesssim 0.1$ 은 이 인자가 아직 완만한 쪽을 고른 것이다.

이분법의 끝도 부호에서 나온다. $d_i$ 가 음수인 성분이 0이 되는 지점은 $\alpha=-f_i/d_i$ 이고, 그중 최소가 $a_{\max}$ 다. 탐색은 $\min(0.999\,a_{\max},\,4)$ 에서 끊긴다. 4를 넘거나 분포의 부호를 바꾸는 근은 후보에서 빠진다.

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

앞의 곡선이 이 $G$ 다. 노란 점은 비자명 근 $\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를 넘는 순간이 온다. 그때 빨간 경고가 켜진다. 그 H 증가가 저점성에서 격자를 터뜨리는 시작이다.

구현하며 부딪힌 세 가지#

노드마다 비선형 풀이를 돌리면 비용이 점성보다 먼저 부담이 된다. 평형 근처, $\|f - f^{eq}\|$ 가 작을 때는 $\alpha \approx 2$ 로 두고, 비평형이 임계를 넘는 격자에서만 이분법을 켠다. 비용은 대략 1할까지 내려간다.

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

2차 절단 다항식을 평형으로 그대로 쓰면 H 최소라는 보증이 빠진다. 안정자의 이득을 다 가져가려면 평형 자체가 엔트로픽 곱이어야 한다. 다항식 위에 $\alpha$ 만 올리면 절반이다.

α를 2에 고정하면 β를 올려도 분포가 다시 음수가 된다#

평형 근처에서 $\alpha \approx 2$ 로 건너뛰는 것은, $\|f - f^{eq}\|$ 가 작을 때의 비용 절약이다. $\beta$ 를 1 쪽으로 밀면서 그 단축을 전 노드에 고정하면 유효 완화율 $\alpha\beta$ 가 다시 2에 붙는다. 거울상까지의 과조화가 돌아오고, 이산 H는 충돌 전보다 커질 수 있다. 점성만 낮추고 근은 풀지 않는 선택이다.

충돌이 H를 평평하게 잡아도, 이웃에서 이미 음수가 된 분포를 스트리밍이 실어 오면 그 근은 늦다. 부호가 이류 단계에서 깨질 때 양수 보존 플럭스 한정자가 아직 필요한지, 엔트로피 조건만으로 막히는지는 다음 갈림으로 남아 있다.