BGK가 발산하는 저점성에서 살아남는 법 — MRT 격자 볼츠만 충돌

같은 격자에서 같은 평형 분포를 향해 이완하는데, 한 충돌 모델은 매끄럽게 돌고 다른 하나는 체커보드 무늬로 터진다. 차이는 단 하나다. BGK는 아홉 개 분포함수를 하나의 완화시간으로 한꺼번에 이완시키고, MRT(Multiple Relaxation Time, 다중 완화시간)는 그것들을 모멘트로 바꿔 제각기 다른 속도로 이완시킨다. 이 글은 그 차이가 왜 안정성을 가르는지 D2Q9 격자 볼츠만으로 따라간다. 변환 행렬 $M$으로 분포를 모멘트로 옮기고, 점성과 무관한 자유도를 따로 조여 저점성에서 BGK가 터지는 영역을 MRT가 어떻게 버티는지 Taylor–Green 와류로 직접 확인한다.

한 박자로 모든 것을 이완시키는 BGK의 약점#

격자 볼츠만(LBM)은 분포함수 $f_i$가 격자 위를 이동(streaming)하고, 한 점에서 평형으로 이완(collision)하는 두 단계를 반복한다. 가장 흔한 BGK 충돌은 이렇게 쓴다.

여기서 $f_i$는 속도 $\mathbf{c}_i$ 방향 분포함수, $\tau$는 완화시간, $f_i^{\text{eq}}$는 국소 평형 분포다. 점성은 $\tau$ 하나로 정해진다.

문제는 여기서 시작한다. 점성을 낮추려면(Reynolds 수를 키우려면) $\tau$를 $1/2$에 바짝 붙여야 한다. 그런데 $\tau \to 1/2$가 되면 점성과 아무 관계 없는 고차 모드들까지 거의 이완되지 않은 채 격자를 돌아다닌다. 이 모드들이 격자 해상도 한계에서 진폭을 키우면 체커보드 진동이 생기고, 해석은 발산한다. BGK는 점성 조절과 안정성 조절이 하나의 손잡이에 묶여 있다.

충돌을 모멘트 공간에서 바라보면#

핵심 아이디어는 단순하다. 분포함수 아홉 개를 그대로 이완시키지 말고, 먼저 물리적 의미가 있는 모멘트로 바꾼 다음 이완시키자. 모멘트는 분포함수의 선형 결합이다.

$\mathbf{f} = (f_0,\dots,f_8)^\top$는 분포함수 벡터, $M$은 $9\times 9$ 변환 행렬, $\mathbf{m}$은 모멘트 벡터다. 충돌은 모멘트 공간에서 평형 모멘트 $\mathbf{m}^{\text{eq}}$로 이완시키고, 다시 속도 공간으로 되돌린다.

$S = \mathrm{diag}(s_0,\dots,s_8)$는 모멘트별 완화율을 담은 대각 행렬이다. 만약 모든 $s_i$를 같은 값 $1/\tau$로 두면 이 식은 정확히 BGK로 되돌아간다. MRT는 그 대각 성분을 모멘트마다 따로 푼다는 점이 전부다.

D2Q9의 아홉 모멘트와 변환 행렬 M#

D2Q9 격자(2차원, 9속도)에서 모멘트는 다음 아홉 개다. 밀도 $\rho$, 에너지 $e$, 에너지 제곱 $\varepsilon$, 운동량 $j_x, j_y$, 에너지 플럭스 $q_x, q_y$, 그리고 응력에 해당하는 $p_{xx}, p_{xy}$.

이 가운데 $\rho, j_x, j_y$는 충돌로 변하지 않는 보존량이다. 따라서 이들의 완화율은 $s=0$으로 고정한다. 나머지 여섯 개가 조절 가능한 손잡이다. Lallemand–Luo 표준 변환 행렬은 각 행이 직교하도록 설계되어, 역행렬을 norm으로 나눈 전치로 간단히 얻는다.

평형 모멘트는 밀도와 운동량만으로 닫힌다.

$j_x = \rho u_x$, $j_y = \rho u_y$는 운동량 성분이다. 응력 모멘트 $p_{xx}, p_{xy}$의 평형이 속도 구배와 연결되어, 이들의 완화율이 점성을 결정한다.

모멘트마다 다른 완화율 — 점성과 무관한 자유도#

전단 점성을 정하는 것은 응력 모멘트의 완화율 $s_\nu$다.

핵심은 나머지 완화율 $s_e$(체적), $s_\varepsilon$, $s_q$(에너지 플럭스)가 점성에 전혀 영향을 주지 않는다는 점이다. 이들은 자유 파라미터다. 저점성에서 $s_\nu$를 $2$에 붙여도, 고차 모드의 완화율은 안정 구간 $1.5\sim1.8$에 따로 둘 수 있다. BGK가 묶어 두었던 두 손잡이가 분리된다.

아래 막대에서 직접 완화율을 조작해보자.

$s_\nu$를 $1.99$까지 올리면 점성 $\nu$가 거의 $0$으로 떨어지는데, 이때도 $s_e, s_q$는 그대로 안정 구간에 둘 수 있음을 본다. 보존 모멘트($\rho, j_x, j_y$) 막대가 항상 0에 고정된 것도 확인하자.

충돌은 모멘트 공간, 스트리밍은 속도 공간#

알고리즘 흐름은 다음과 같다.

1. 분포 $f_i$에서 거시량 $\rho, u_x, u_y$ 계산
2. $\mathbf{m} = M\mathbf{f}$로 모멘트 변환
3. 평형 모멘트 $\mathbf{m}^{\text{eq}}$ 계산 (밀도·운동량으로)
4. 모멘트별 이완: $m_a \leftarrow m_a - s_a(m_a - m_a^{\text{eq}})$
5. $\mathbf{f}^{*} = M^{-1}\mathbf{m}$로 속도 공간 복원
6. 스트리밍: $f_i^{*}$를 $\mathbf{c}_i$ 방향 이웃으로 이동
7. 경계조건 적용 후 1로

충돌만 모멘트 공간에서 일어나고, 스트리밍은 그대로 속도 공간에서 진행된다. 추가 비용은 셀마다 $9\times 9$ 행렬 곱 두 번뿐이다.

Python — Taylor–Green 와류로 점성 감쇠를 잰다#

검증에는 Taylor–Green 와류를 쓴다. 이 유동은 점성 감쇠 해가 정확히 알려져 있어 구현이 맞는지 숫자로 확인할 수 있다. 운동에너지는 $\exp(-4\nu k^2 t)$로 줄어야 한다($k$는 와류 파수).

import numpy as np
 
# D2Q9 lattice
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])
 
# Lallemand-Luo transformation matrix (rows: rho,e,eps,jx,qx,jy,qy,pxx,pxy)
M = np.array([
    [ 1, 1, 1, 1, 1, 1, 1, 1, 1],
    [-4,-1,-1,-1,-1, 2, 2, 2, 2],
    [ 4,-2,-2,-2,-2, 1, 1, 1, 1],
    [ 0, 1, 0,-1, 0, 1,-1,-1, 1],
    [ 0,-2, 0, 2, 0, 1,-1,-1, 1],
    [ 0, 0, 1, 0,-1, 1, 1,-1,-1],
    [ 0, 0,-2, 0, 2, 1, 1,-1,-1],
    [ 0, 1,-1, 1,-1, 0, 0, 0, 0],
    [ 0, 0, 0, 0, 0, 1,-1, 1,-1],
], dtype=float)
Minv = np.linalg.inv(M)
 
def d2q9_equilibrium(rho, ux, uy):
    eq = np.zeros((9,) + rho.shape)
    for i in range(9):
        cu = cx[i] * ux + cy[i] * uy
        eq[i] = w[i] * rho * (1 + 3*cu + 4.5*cu**2 - 1.5*(ux**2 + uy**2))
    return eq
 
def relaxation_rates(nu):
    s_nu = 1.0 / (3*nu + 0.5)          # 점성을 정하는 완화율
    # 보존(rho,jx,jy)=0, 고차 모드는 안정 구간의 고정값
    return np.array([0, 1.64, 1.54, 0, 1.7, 0, 1.7, s_nu, s_nu])
 
def taylor_green_init(N, U0):
    x = np.arange(N)
    X, Y = np.meshgrid(x, x, indexing='ij')
    k = 2*np.pi / N
    ux = -U0 * np.cos(k*X) * np.sin(k*Y)
    uy =  U0 * np.sin(k*X) * np.cos(k*Y)
    rho = 1 - 0.75 * U0**2 * (np.cos(2*k*X) + np.cos(2*k*Y))
    return d2q9_equilibrium(rho, ux, uy), k
 
def macros(f):
    rho = f.sum(axis=0)
    ux = (cx[:, None, None] * f).sum(axis=0) / rho
    uy = (cy[:, None, None] * f).sum(axis=0) / rho
    return rho, ux, uy
 
def mrt_collide(f, s):
    rho, ux, uy = macros(f)
    jx, jy = rho*ux, rho*uy
    m = np.tensordot(M, f, axes=([1], [0]))        # 모멘트 공간으로
    meq = np.zeros_like(m)
    meq[0] = rho
    meq[1] = -2*rho + 3*(jx**2 + jy**2)
    meq[2] = rho - 3*(jx**2 + jy**2)
    meq[3], meq[4] = jx, -jx
    meq[5], meq[6] = jy, -jy
    meq[7] = jx**2 - jy**2
    meq[8] = jx*jy
    m -= s[:, None, None] * (m - meq)               # 모멘트별 이완
    return np.tensordot(Minv, m, axes=([1], [0]))    # 속도 공간 복원
 
def stream(f):
    for i in range(9):
        f[i] = np.roll(f[i], (cx[i], cy[i]), axis=(0, 1))
    return f
 
def run_mrt_lbm(N=64, nu=0.004, U0=0.04, steps=2000):
    f, k = taylor_green_init(N, U0)
    s = relaxation_rates(nu)
    e0 = None
    for t in range(steps):
        f = mrt_collide(f, s)
        f = stream(f)
        if t % 400 == 0:
            _, ux, uy = macros(f)
            E = np.mean(ux**2 + uy**2)
            e0 = E if e0 is None else e0
            print(f"t={t:5d}  KE/KE0={E/e0:.4f}  analytic={np.exp(-4*nu*k**2*t):.4f}")
 
run_mrt_lbm()

출력은 수치 감쇠가 해석 곡선 $\exp(-4\nu k^2 t)$와 소수점 둘째 자리까지 맞는 것을 보여준다. 같은 격자·같은 $\nu$로 BGK를 돌리면 더 낮은 점성에서 먼저 발산한다.

BGK와 MRT를 같은 Re에서 나란히#

아래 시뮬레이션에서 점성 슬라이더를 내리고 BGK/MRT 버튼으로 충돌 모델을 바꿔보자.

점성을 $0.002$ 부근까지 낮춘 뒤 BGK로 두면 격자에 거친 무늬가 끼며 "diverged"가 뜬다. 같은 점성에서 MRT로 바꾸고 reset하면 와류가 매끄럽게 감쇠한다. 두 모델이 같은 Reynolds 수를 목표로 하는데도 안정 한계가 다른 것이 한눈에 보인다.

다음에 격자 볼츠만을 만질 때#

MRT는 마법이 아니다. 점성을 정하는 손잡이($s_\nu$)와 안정성을 정하는 손잡이($s_e, s_q$)를 물리적으로 분리했을 뿐이다. 그 분리가 저점성·고Reynolds에서 격자 볼츠만의 작동 범위를 넓힌다.

세 가지만 남긴다. 첫째, BGK가 저점성에서 터지는 이유는 점성과 무관한 고차 모드가 함께 느려지기 때문이다. 둘째, MRT는 충돌만 모멘트 공간으로 옮겨 그 모드를 따로 조이고, 추가 비용은 셀당 행렬 곱 두 번이다. 셋째, 구현이 맞는지는 Taylor–Green 와류의 $\exp(-4\nu k^2 t)$ 감쇠로 항상 검증할 수 있다.