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

같은 격자, 같은 평형인데 한쪽은 매끄럽게 돌고 다른 쪽은 체커보드로 터진다. BGK는 아홉 개 분포함수를 하나의 완화시간으로 민다. MRT(Multiple Relaxation Time, 다중 완화시간)는 모멘트로 바꾼 뒤 제각기 다른 속도로 민다.

$s_\nu$ 막대를 1.99까지 밀면 점성 숫자가 0 쪽으로 꺼진다. $s_e$와 $s_q$ 막대는 그 자리에 남는다.

보존량 막대 $\rho$, $j_x$, $j_y$는 0에 붙어 있다. BGK는 이 막대들을 한 값으로 묶는다. 저점성에서 체커보드가 나오는 자리가 그 묶음이다.

한 박자로 모든 것을 이완시키는 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가 한 손잡이에 묶던 것이 여기서 갈라진다.

코드의 기본 점성은 $\nu=0.004$다. $s_\nu = 1/(3\nu + 1/2)$라서 $3\times 0.004 + 0.5 = 0.512$, $s_\nu = 1/0.512 = 1.953125$다. 앞의 막대를 1.99까지 밀면 $\nu = (1/3)(1/1.99 - 1/2)$이다. $1/1.99 \approx 0.5025$에서 0.5를 빼면 약 0.0025, 3으로 나누면 $\nu \approx 8.4\times 10^{-4}$다. 0에 가깝지만 0은 아니다. 그 점성에서 BGK는 $s_e$, $s_\varepsilon$, $s_q$까지 1.99로 같이 민다. 코드의 MRT는 그 세 자리를 $1.64$, $1.54$, $1.7$에 고정한다. 돌아가는 배열은 [0, 1.64, 1.54, 0, 1.7, 0, 1.7, s_nu, s_nu]라서 0은 보존량 세 칸뿐이다.

평형 모멘트가 밀도·운동량만으로 닫히는 것은 분포를 낮은 차수로 자른 뒤의 일이다. 그 자르기는 9개 속도 안의 Maxwell–Boltzmann 쪽에 있다.

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

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

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()

print가 나란히 놓는 해석값은 $\exp(-4\nu k^2 t)$다. 기본값은 $N=64$, $\nu=0.004$, $k=2\pi/N$이라 $4\nu k^2 = 0.016\,(\pi/32)^2 \approx 1.54\times 10^{-4}$다. $t=400$이면 지수가 약 $-0.0617$, 해석 인자는 약 0.940이다. $t=1600$이면 약 $-0.247$, 인자는 약 0.781이다. range(2000)이라 찍히는 마지막 $t$는 1600이다. 소수점 둘째 자리까지 맞는다는 말은 이 인자와 수치 운동에너지의 비교다. 같은 격자·같은 $\nu$로 BGK를 돌리면 더 낮은 점성에서 먼저 발산한다.

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

점성 슬라이더를 내리면 감쇠가 느려지고, BGK/MRT 버튼은 같은 $\nu$에서 격자가 버티는지를 가른다.

점성을 $0.002$ 부근까지 낮춘 뒤 BGK로 두면 격자에 거친 무늬가 끼며 "diverged"가 뜬다. 같은 점성에서 MRT로 바꾸고 reset하면 와류가 매끄럽게 감쇠한다. $\nu=0.002$면 $s_\nu=1/(3\times 0.002+0.5)=1/0.506\approx 1.976$이다. BGK는 고차 모드까지 이 값으로 민다. 코드의 $1.64$, $1.54$, $1.7$은 그 옆에 남아 있다.

s_ν만 2에 붙이면 고차 모드가 남는다#

점성 손잡이 $s_\nu$와 안정성 손잡이 $s_e$, $s_q$를 물리적으로 분리한 것이 MRT의 전부다. $\nu=0.004$의 $s_\nu=1.953$을 모든 모드에 복사하면 BGK다. 배열의 $1.64$, $1.54$, $1.7$을 그대로 두고 $s_\nu$만 올리는 쪽이 MRT다. 추가 비용은 셀당 $9\times 9$ 행렬 곱 두 번이다.

그 분리가 안 되면 저점성 Taylor–Green은 체커보드로 죽는다. 엔트로피 제한으로 같은 폭주를 막는 다른 길은 H 정리를 충돌에 넣은 격자 볼츠만이다. 감쇠 확인은 $\exp(-4\nu k^2 t)$ 한 줄로 충분하다.