LBM에서 힘의 절반은 어디로 갔는가 — forcing 스킴 네 가지와 1−1/(2τ)

1993년 Shan과 Chen은 격자 위에서 물과 기름을 갈라놓았다. 처방은 소박했다. 충돌에 쓰는 평형 속도에 $\tau\mathbf{F}/\rho$를 얹은 것이 전부였다. 1998년 He는 힘을 평형분포에 직접 물렸다. 2002년 Guo는 앞의 처방들이 응력에 오차를 남긴다며 힘 항 앞에 $(1-\frac{1}{2\tau})$를 붙였다. 2004년 Kupershtokh는 평형분포 두 개의 차이만으로 같은 일을 해냈다.

같은 질문에 답이 네 개다. 그래서 무엇을 써야 하는가. 오늘 이 질문을 D2Q9 채널에 직접 물었다. 답은 예상과 달랐다. 정상 유동에서는 어느 스킴을 써도 소수점 넷째 자리까지 똑같은 오차가 나왔다. 진짜 함정은 스킴을 고르는 자리에 있지 않았다. 이 글은 유도로 그 함정의 위치를 짚고, 숫자로 함정의 크기를 잰다.

연속 방정식에서 힘 항을 잘라낸다#

외력이 있는 Boltzmann 방정식은 항이 하나 늘어난다.

$f$는 분포함수, $\boldsymbol{\xi}$는 입자 속도, $\mathbf{F}$는 단위 부피당 체적력, $\lambda$는 완화 시간이다.

말썽은 세 번째 항이다. 속도 공간의 미분 $\nabla_{\boldsymbol{\xi}} f$는 격자 위에 없다. 우리가 가진 것은 아홉 방향의 이산 속도뿐이다. 그래서 $f$를 $f^{\rm eq}$로 갈아 끼워 미분을 해석적으로 처리한다.

$c_s$는 격자 음속(D2Q9에서 $c_s^2=1/3$), $\mathbf{u}$는 유체 속도다.

$f$를 $f^{\rm eq}$로 바꿔치기해도 되는가. 저마하수에서는 된다. 비평형 부분 $f^{(1)}$은 Knudsen 수(평균자유행로/특성길이)의 1차 크기라 힘과 곱해지면 무시할 수 있는 차수로 내려간다. 대신 이 교체가 나중에 응력에 잔여 오차를 남긴다. 스킴들이 갈리는 지점이 바로 여기다.

이 표현을 이산 속도 $\mathbf{c}_i$에 사영하고 Hermite 급수를 2차에서 자른다. 그러면 격자 위의 힘 항이 나온다.

$w_i$는 격자 가중치다. 하필 2차에서 자르는 이유는 분명하다. Navier–Stokes를 되찾는 데 필요한 모멘트가 2차까지이기 때문이다.

세 모멘트를 확인해 둔다.

0차는 질량이다. 힘은 질량을 만들지 않으니 0이어야 한다. 1차는 운동량, 2차는 응력이다.

사다리꼴 적분이 남긴 절반#

이제 이산 속도 Boltzmann 방정식을 특성선을 따라 $\delta t$만큼 적분한다. 충돌항과 힘 항을 사다리꼴 규칙으로 적분하면 2차 정확도를 얻는다. 대신 우변에 $t+\delta t$가 끼어든다.

윗첨자 $n$과 $n+1$은 각각 $t$와 $t+\delta t$를 뜻한다. 이대로는 음함수다. 매 스텝 연립방정식을 풀어야 한다.

LBM은 이 음함수성을 변수 변환으로 지운다.

$\bar f_i$로 다시 쓰고 $\tau = \lambda/\delta t + 1/2$로 두면 익숙한 한 줄이 된다.

변환이 흔적을 두 개 남겼다. 하나는 힘 앞에 붙은 $(1-\frac{1}{2\tau})$다. 다른 하나는 눈에 덜 띈다. 우리가 배열에 저장하는 것은 $\bar f_i$이지 $f_i$가 아니다. 그래서 운동량도 절반만큼 어긋난다.

이 둘은 한 몸이다. 같은 변환에서 같이 떨어져 나왔다. 하나만 챙기고 다른 하나를 빠뜨리면 회계가 맞지 않는다. 얼마나 안 맞는지는 아래에서 잰다.

네 스킴 비교표#

$\mathbf{j}=\sum_i \bar f_i\mathbf{c}_i$를 저장된 분포함수의 생 운동량, $\delta\mathbf{u}=\delta t\,\mathbf{F}/\rho$를 힘이 한 스텝에 주는 속도 증분이라 하자.

다섯 줄을 세로로 읽으면 공통점이 보인다. 계수를 곱하기 전 힘 항의 1차 모멘트는 전부 정확히 $\mathbf{F}$다. 그렇게 설계했으니 당연하다. 갈리는 것은 2차 모멘트, 곧 응력이다.

그런데 변수 변환이 요구하는 목표값은 $\mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u}$가 아니라 $(1-\frac{1}{2\tau})(\mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u})$다. 이 계수를 명시적으로 짊어진 것은 맨 아랫줄 하나뿐이다. 나머지는 $\frac{1}{2\tau}$만큼의 응력 편차를 안고 간다.

편차의 크기는 힘과 속도의 곱, 즉 $u F$ 차수다. 저마하수에서 $u$가 작으면 이 항 자체가 작다. Shan–Chen과 EDM이 남기는 $\mathbf{F}\mathbf{F}$ 항은 더 작아서 $F^2$ 차수다. 그래서 약한 중력으로 미는 문제에서는 어느 것을 써도 별 차이가 없다. 차이가 드러나는 곳은 힘이 세거나($\tau\mathbf{F}/\rho$가 $\mathbf{u}$에 필적할 때) 다상유동처럼 계면에서 힘이 급격히 변하는 자리다.

아래 도형에서 직접 조작해보자.

힘 방향 다이얼을 아무리 돌려도 M0(질량) 막대는 0에 붙어 있다. 힘은 질량을 만들지 않으니 당연하다. M1(운동량)의 목표는 $\mathbf{F}$가 아니라 $(1-\frac{1}{2\tau})\mathbf{F}$다. 계수가 곱해진 뒤의 힘 항을 재고 있기 때문이다. 모자란 $\mathbf{F}/(2\tau)$는 반 칸 옮긴 평형이 되돌려주고, 맨 아랫줄 $\Delta j$가 정확히 $\mathbf{F}$로 떨어지는 것이 그 확인이다.

스킴이 갈리는 곳은 M2(응력) 한 줄뿐이다. $|u|$ 슬라이더를 올리면 plain의 M2는 0에 머문 채 목표를 놓친다. 크기가 아니라 모양이 틀렸으니 계수로는 고쳐지지 않는다. $\tau$를 2 쪽으로 밀면 Shan–Chen의 편차가 눈에 띄게 자란다. prefactor 스위치를 끄면 M1과 M2와 $\Delta j$가 한꺼번에 어긋난다.

정상 Poiseuille은 이들을 구별하지 못한다#

가장 흔한 검증 문제로 가 보자. 위아래 벽 사이를 체적력으로 미는 채널이다. 정상 상태의 해는 포물선으로 알려져 있다.

$H$는 채널 높이, $\nu$는 동점성계수다. half-way bounce-back을 쓰면 벽이 격자점에서 반 칸 밖에 놓이므로 $H$는 유체 노드 수와 같고 $j$번 노드의 좌표는 $y=j+0.5$다.

아래 시뮬레이션에서 직접 조작해보자.

스킴 버튼을 네 개 다 눌러 보라. 점선 포물선 위에서 청록색 곡선이 움직이지 않는다. $\tau$를 0.6에서 2.0까지 밀어도 마찬가지다. 그다음 아래 두 스위치 중 하나만 꺼 보라. 그 순간 포물선의 진폭이 통째로 갈아탄다.

숫자로도 같다. 유체 노드 33개, $F_x = 10^{-5}$로 40,000스텝을 돌린 뒤 잰 최대 상대오차다.

열끼리 비교할 것이 없다. 값이 같다. 남은 오차는 스킴이 아니라 bounce-back의 이산화에서 온다.

이유는 2차 모멘트가 들어가는 자리에 있다. 이 유동은 정상이고 단방향이다. 속도는 $u_x(y)$ 하나뿐이고 힘도 $x$ 성분만 있다. 스킴이 갈리는 항은 $\mathbf{u}\mathbf{F}$, 곧 $xx$ 성분이다. 그런데 이 유동의 운동량 균형을 실제로 지배하는 것은 $xy$ 전단응력이다. 힘이 만드는 $xx$ 편차는 균형식에 들어갈 통로가 없다.

$\partial_t\mathbf{u}=0$인 것도 같이 작용한다. 힘 항의 오차가 시간에 누적되어 나타날 여지가 없다. 스킴 차이를 보려면 비정상 유동이나 공간적으로 변하는 힘으로 가야 한다. 적어도 이 시험으로는 다섯을 가려낼 수 없다.

이건 나쁜 소식이 아니다. 검증 문제를 고를 때 알아야 할 정보다. Poiseuille가 통과했다고 forcing 구현이 옳다고 말할 수 없다는 뜻이니까.

절반을 두 번 세면 힘이 커진다#

그러면 무엇이 실제로 틀어지는가. 앞에서 한 몸이라 한 두 가지, 즉 속도의 반 칸 보정과 힘 앞의 계수다.

한 스텝에 실제로 주입되는 운동량을 세어 보자. 평형 속도에 $\delta t\mathbf{F}/2\rho$를 얹으면 충돌이 그중 $1/\tau$만큼을 매 스텝 실어 나른다. 여기에 힘 항이 직접 주는 몫이 더해진다. 반 칸 보정을 켰는지를 $s\in\{0,1\}$, 힘 앞의 계수를 $g$라 하면

이것이 $\delta t F$와 같아야 한다. 조건은 $g = 1 - s/(2\tau)$ 하나다. 두 스위치는 독립이 아니다.

어긋내면 진폭이 예측 가능한 비율로 어긋난다. Guo 힘 항을 고정하고 두 스위치만 바꿔 측정한 진폭비다.

측정과 예측이 소수점 셋째 자리까지 맞는다. $\tau=0.6$에서 반 칸 보정만 켜면 힘이 83% 커진다. 계수만 켜면 힘이 6분의 1로 줄어든다.

이건 미묘한 정확도 문제가 아니다. 논문에서 Guo의 $F_i$를 그대로 베껴 오고 속도는 예전 코드의 $\mathbf{j}/\rho$를 그냥 둘 때 정확히 이 일이 벌어진다. 점성이 이상하다며 $\tau$를 만지기 시작하면 더 깊이 들어간다. $\tau$를 바꿀 때마다 오차 비율이 같이 움직이기 때문이다.

증상으로 구별하는 법이 있다. 격자를 두 배로 늘려도 오차 비율이 줄지 않으면 이산화 오차가 아니다. $\tau$를 1에 가깝게 가져갈수록 비율이 1.5 근처로 수렴하고 $\tau$를 0.5로 내릴수록 발산하듯 커지면 $1+\frac{1}{2\tau}$의 지문이다. 반대로 힘을 키워도 비율이 그대로면 회계 문제이지 비선형 문제가 아니다. 세 가지가 겹치면 스킴을 의심하지 말고 속도 정의부터 열어 보면 된다.

Python — 두 스위치를 어긋내 본다#

두 스위치를 인자로 뺀 D2Q9 채널 솔버다. $x$ 방향으로 균질하므로 열 하나만 남긴다. 위 표의 다섯 스킴은 source_terms와 평형에 넣는 속도만 갈아 끼우면 되므로 여기서는 Guo 힘 항 하나로 고정한다.

import numpy as np
 
EX = np.array([0, 1, 0, -1, 0, 1, -1, -1, 1])
EY = 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])
OPP = np.array([0, 3, 4, 1, 2, 7, 8, 5, 6])
CS2 = 1.0 / 3.0
 
 
def lattice_equilibrium(rho, ux, uy):
    feq = np.empty((9, rho.size))
    usq = ux**2 + uy**2
    for i in range(9):
        eu = EX[i] * ux + EY[i] * uy
        feq[i] = W[i] * rho * (1 + eu/CS2 + eu**2/(2*CS2**2) - usq/(2*CS2))
    return feq
 
 
def source_terms(rho, ux, uy, fx):
    """Guo forcing term, before the (1 - 1/2tau) gain."""
    src = np.empty((9, rho.size))
    for i in range(9):
        eu = EX[i]*ux + EY[i]*uy
        src[i] = W[i] * ((EX[i] - ux)/CS2 + eu*EX[i]/CS2**2) * fx
    return src
 
 
def run_forced_channel(tau, half_shift, prefactor, ny=33, fx=1.0e-5, steps=40000):
    rho = np.ones(ny)
    f = lattice_equilibrium(rho, np.zeros(ny), np.zeros(ny))
    gain = (1.0 - 1.0/(2*tau)) if prefactor else 1.0
    for _ in range(steps):
        rho = f.sum(axis=0)
        jx = (f*EX[:, None]).sum(axis=0)
        jy = (f*EY[:, None]).sum(axis=0)
        ux = (jx + (0.5*fx if half_shift else 0.0)) / rho   # switch 1
        uy = jy / rho
        fpost = f - (f - lattice_equilibrium(rho, ux, uy))/tau \
                + gain*source_terms(rho, ux, uy, fx)        # switch 2
        for i in range(9):                                   # streaming in y
            f[i] = fpost[i] if EY[i] == 0 else np.roll(fpost[i], EY[i])
        for i in range(9):                                   # half-way bounce-back
            if EY[i] > 0:
                f[i, 0] = fpost[OPP[i], 0]
            elif EY[i] < 0:
                f[i, -1] = fpost[OPP[i], -1]
    jx = (f*EX[:, None]).sum(axis=0)
    return (jx + 0.5*fx) / f.sum(axis=0)      # physical velocity, always shifted
 
 
NY, FX = 33, 1.0e-5
y = np.arange(NY) + 0.5
for tau in (0.6, 1.0, 1.8):
    exact = FX / (2*CS2*(tau - 0.5)) * y * (NY - y)
    ok = run_forced_channel(tau, True, True).max() / exact.max()
    m1 = run_forced_channel(tau, True, False).max() / exact.max()
    m2 = run_forced_channel(tau, False, True).max() / exact.max()
    print(f"tau={tau:4.2f}  both={ok:.4f}  shift_only={m1:.4f} (pred {1+1/(2*tau):.4f})"
          f"  gain_only={m2:.4f} (pred {1-1/(2*tau):.4f})")

출력은 이렇다.

tau=0.60  both=0.9991  shift_only=1.8316 (pred 1.8333)  gain_only=0.1666 (pred 0.1667)
tau=1.00  both=1.0003  shift_only=1.5002 (pred 1.5000)  gain_only=0.5005 (pred 0.5000)
tau=1.80  both=1.0074  shift_only=1.2867 (pred 1.2778)  gain_only=0.7280 (pred 0.7222)

$\tau=1.8$에서 예측과 0.7% 어긋나는데, 같은 조건에서 스킴 자체의 이산화 오차가 0.74%다. 크기가 같다. 어긋남의 출처가 회계가 아니라 격자라는 뜻이다.

기억할 점#

• 힘 항 앞의 $(1-\frac{1}{2\tau})$와 속도의 $+\delta t\mathbf{F}/2\rho$는 사다리꼴 적분의 변수 변환에서 함께 나온 한 쌍이다. 둘 중 하나만 쓰면 힘이 $1\pm\frac{1}{2\tau}$배로 어긋난다. $\tau=0.6$이면 83% 초과다.
• 정상 단방향 Poiseuille은 forcing 스킴을 가려내지 못한다. 다섯 스킴이 넷째 자리까지 같은 오차를 낸다. 스킴 비교는 비정상 유동이나 비균질 힘에서 해야 한다.
• 어떤 스킴이든 힘 항 자체의 1차 모멘트는 $\mathbf{F}$로 맞춰져 있다. 차이는 2차 모멘트, 곧 응력에 있고 목표값에 $(1-\frac{1}{2\tau})$가 들어 있다. 새 스킴을 검토할 때 M2부터 확인하면 빠르다.