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

속도에 δt F/2ρ가 있는데 힘 항에 (1−1/(2τ))가 없으면, τ=0.6에서 주입 운동량이 1.83배가 된다.

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

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

모멘트 다이얼에서 스킴, $\tau$ , $|u|$ , prefactor를 바꾸면 M0·M1·M2 막대 높이가 같이 움직인다.

τ1.00F angle θ0°|u|0.080

다이얼을 돌려도 M0은 0이다. $|u|$ 를 올리면 plain의 M2만 목표에서 떨어지고, prefactor를 끄면 M1과 $\Delta j$ 가 한 번에 틀린다.

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

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

\partial_t f + \boldsymbol{\xi}\cdot\nabla_{\mathbf{x}} f + \frac{\mathbf{F}}{\rho}\cdot\nabla_{\boldsymbol{\xi}} f = -\frac{1}{\lambda}\left(f - f^{\rm eq}\right)

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

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

\frac{\mathbf{F}}{\rho}\cdot\nabla_{\boldsymbol{\xi}} f \simeq -\frac{(\boldsymbol{\xi}-\mathbf{u})\cdot\mathbf{F}}{\rho c_s^2}\,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차에서 자른다. 그러면 격자 위의 힘 항이 나온다.

F_i = w_i\left[\frac{\mathbf{c}_i-\mathbf{u}}{c_s^2} + \frac{(\mathbf{c}_i\cdot\mathbf{u})}{c_s^4}\mathbf{c}_i\right]\cdot\mathbf{F}

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

세 모멘트를 확인해 둔다.

\sum_i F_i = 0,\qquad \sum_i \mathbf{c}_i F_i = \mathbf{F},\qquad \sum_i \mathbf{c}_i\mathbf{c}_i F_i = \mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u}

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

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

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

f_i(\mathbf{x}+\mathbf{c}_i\delta t,\ t+\delta t) - f_i(\mathbf{x},t) = -\frac{\delta t}{2\lambda}\left[(f_i-f_i^{\rm eq})^{n} + (f_i-f_i^{\rm eq})^{n+1}\right] + \frac{\delta t}{2}\left[F_i^{\,n} + F_i^{\,n+1}\right]

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

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

\bar f_i = f_i + \frac{\delta t}{2\lambda}\left(f_i - f_i^{\rm eq}\right) - \frac{\delta t}{2}F_i

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

\bar f_i(\mathbf{x}+\mathbf{c}_i\delta t,\ t+\delta t) = \bar f_i - \frac{1}{\tau}\left(\bar f_i - f_i^{\rm eq}\right) + \delta t\left(1-\frac{1}{2\tau}\right)F_i

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

\rho\mathbf{u} = \sum_i \bar f_i \mathbf{c}_i + \frac{\delta t}{2}\mathbf{F}

이 둘은 한 몸이다. 같은 변환에서 같이 떨어져 나왔다. 사다리꼴 적분이 점성에 남긴 반 스텝도 이 반 칸과 한 뿌리에서 나온다. 하나만 챙기고 다른 하나를 빠뜨리면 회계가 맞지 않는다. 얼마나 안 맞는지는 아래에서 잰다.

네 스킴 비교표#

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

스킴	힘을 넣는 자리	평형에 쓰는 속도	2차 모멘트 $\sum\mathbf{c}_i\mathbf{c}_i F_i$
plain	$F_i = w_i(\mathbf{c}_i\cdot\mathbf{F})/c_s^2$	$\mathbf{j}/\rho$	0 — $\mathbf{u}\mathbf{F}$ 항이 통째로 없음
Shan–Chen (1993)	명시적 항 없음	$\mathbf{j}/\rho + \tau\mathbf{F}/\rho$	$\mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u} + \tau\mathbf{F}\mathbf{F}/\rho$
He (1998)	$(\mathbf{c}_i-\mathbf{u})\cdot\mathbf{F}\,f_i^{\rm eq}/(\rho c_s^2)$	$(\mathbf{j}+\delta t\mathbf{F}/2)/\rho$	$\mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u}$ (평형 절단 오차 내)
EDM (2004)	$f_i^{\rm eq}(\rho,\mathbf{u}+\delta\mathbf{u}) - f_i^{\rm eq}(\rho,\mathbf{u})$	$\mathbf{j}/\rho$	$\mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u} + \rho\,\delta\mathbf{u}\delta\mathbf{u}$
Guo (2002)	위 $F_i$ 에 $(1-\frac{1}{2\tau})$ 를 곱함	$(\mathbf{j}+\delta t\mathbf{F}/2)/\rho$	$(1-\frac{1}{2\tau})(\mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u})$

다섯 줄을 세로로 읽으면 공통점이 보인다. 계수를 곱하기 전 힘 항의 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}$ 에 필적할 때) 다상유동처럼 계면에서 힘이 급격히 변하는 자리다.

M1의 목표는 $\mathbf{F}$ 가 아니라 $(1-\frac{1}{2\tau})\mathbf{F}$ 다. 계수를 곱한 뒤의 힘 항을 재기 때문이다. 모자란 $\mathbf{F}/(2\tau)$ 는 반 칸 옮긴 평형이 되돌려 주고, $\Delta j$ 가 $\mathbf{F}$ 로 떨어진다. M2만 스킴마다 갈린다. plain은 모양이 틀려서 계수로 안 고쳐진다. $\tau$ 를 2 쪽으로 밀면 Shan–Chen의 편차가 자란다.

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

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

u(y) = \frac{F_x}{2\nu}\,y\,(H-y),\qquad \nu = c_s^2\left(\tau-\tfrac{1}{2}\right)

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

스킴을 갈아 끼워도 포물선의 진폭은 그대로다. 유체 노드 33개, $F_x = 10^{-5}$ 로 40,000스텝을 돌린 뒤 잰 최대 상대오차다.

$\tau$	plain	Shan–Chen	He	EDM	Guo
0.60	0.0875%	0.0875%	0.0875%	0.0875%	0.0875%
1.00	0.0306%	0.0306%	0.0306%	0.0306%	0.0306%
1.80	0.7358%	0.7358%	0.7358%	0.7358%	0.7358%

열끼리 비교할 것이 없다. 값이 같다. 남은 오차는 스킴이 아니라 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(\rho u) = \frac{1}{\tau}\cdot\frac{s\,\delta t F}{2} + g\,\delta t F

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

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

$\tau$	둘 다 켬	보정만 켬	예측 $1+\frac{1}{2\tau}$	계수만 켬	예측 $1-\frac{1}{2\tau}$
0.60	0.9991	1.8316	1.8333	0.1666	0.1667
0.80	0.9995	1.6240	1.6250	0.3751	0.3750
1.00	1.0003	1.5002	1.5000	0.5005	0.5000
1.40	1.0030	1.3609	1.3571	0.6452	0.6429
1.80	1.0074	1.2867	1.2778	0.7280	0.7222

측정과 예측이 소수점 셋째 자리까지 맞는다. $\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}$ 의 지문이다. 반대로 힘을 키워도 비율이 그대로면 회계 문제이지 비선형 문제가 아니다. 세 가지가 겹치면 스킴을 의심하지 말고 속도 정의부터 열어 보면 된다.

$\tau=0.6$ 이면 $2\tau=1.2$ 이라 $1/(2\tau)$ 는 $5/6$ , 곧 0.8333이다. 둘 다 켠 상태는 $s=1$ , $g=1-1/(2\tau)=1/6$ 이다. 한 스텝 주입량을 $\delta t F$ 로 나누면 $(1/\tau)\cdot(s/2)+g$ 가 된다. 대입하면 $(1/0.6)\cdot 0.5 + 1/6 = 5/6 + 1/6 = 1$ 이다. 목표와 같다. 표의 both 열 0.9991은 그 1에서 격자만큼만 빠진 값이다.

보정만 켜면 $g$ 가 1로 고정된다. 같은 자리는 $5/6+1=1.8333$ 이다. 표는 1.8316을 준다. 차 0.0017은 예측 대비 약 0.09%다. 계수만 켜면 $s=0$ , $g=1/6$ 이라 주입량은 0.1667이다. 측정 0.1666과 넷째 자리에서 만난다.

$\tau=1$ 이면 $1/(2\tau)=1/2$ 다. 보정만 켜면 1.5, 계수만 켜면 0.5. 표는 1.5002와 0.5005다. 둘 다 켜면 1.0003이라 목표 1에서 0.03%가 남는다. $\tau=1.8$ 이면 $1/(2\tau)=1/3.6=0.2778$ 이다. 예측 1.2778과 0.7222에 대해 측정은 1.2867과 0.7280이다. $(1.2867-1.2778)/1.2778$ 은 약 0.007, 곧 0.7%다. 같은 $\tau$ 의 스킴 상대오차 0.7358%와 크기가 겹친다.

상대오차 표도 같이 읽힌다. $\tau=0.6$ 에서 0.0875%, $\tau=1.0$ 에서 0.0306%, $\tau=1.8$ 에서 0.7358%다. $\nu=c_s^2(\tau-1/2)$ 이고 $c_s^2=1/3$ 이니 $\tau=0.6$ 에서 $\nu=1/30$ , $\tau=1.0$ 에서 $\nu=1/6$ , $\tau=1.8$ 에서 $\nu=13/30$ 이다. $(\tau-1/2)$ 의 비는 0.1 대 0.5 대 1.3, 곧 1 대 5 대 13이다. 봉우리 $u(H/2)=F_x H^2/(8\nu)$ 라서 같은 $F_x=10^{-5}$ , $H=33$ 에서 높이는 $\nu$ 에 반비례한다. $\tau=1.8$ 의 봉우리는 $\tau=0.6$ 의 $1/13$ 이다. 퍼센트가 거기서 커져도 다섯 칸은 한 숫자로 붙어 있다. 스킴이 갈라진 자국이 아니다.

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$ , 반 칸 보정, prefactor를 위 출력의 both, shift_only, gain_only에 맞춰 바꿔 보라.

τ1.00fx8.00e-6

네 스킴을 눌러도 청록 곡선은 점선 포물선에 붙어 있다. 스위치 하나만 끄면 $\tau=0.6$ 에서 진폭이 1.83 또는 0.17 쪽으로 점프한다.

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

반 칸만 켜 두면 힘이 1.83배가 된다#

$\tau=0.6$ 의 shift_only는 1.8316이다. 예측 $1+\frac{1}{2\tau}=1.8333$ 과 셋째 자리까지 같다. Guo 항만 베끼고 속도는 $\mathbf{j}/\rho$ 로 두면, 보정은 켜지고 prefactor는 꺼진 이 칸에 떨어진다. 격자를 두 배로 해도 비율은 안 줄어든다. $\tau$ 를 1로 가져가면 1.5 근처, 0.5로 내리면 비율이 커진다. $1+\frac{1}{2\tau}$ 의 지문이다.

1차 모멘트는 다섯 스킴 모두 $\mathbf{F}$ 다. 갈리는 자리는 2차이고, 목표에 $(1-\frac{1}{2\tau})$ 가 붙어 있다. 정상 채널은 그 $xx$ 편차를 전단에 넣지 못해서 다섯 열이 0.0875%, 0.0306%, 0.7358%로 한 줄이다. 비정상이거나 계면에서 힘이 급할 때 M2를 본다. 2차를 넘어 3차까지 채우면 오차가 되레 커진다는 쪽은 힘 항의 3차 모멘트에 있다.