3차 모멘트를 정직하게 채웠더니 오차가 2.8배 늘었다 — 힘 텀이 지켜야 하는 줄

힘 텀은 몇 차 모멘트까지 맞춰야 하는가#

LBM 코드에 중력이나 표면장력을 넣을 때 우리가 맞추는 것은 세 줄이다. 0차 모멘트 0, 1차 모멘트 $\mathbf{F}$, 2차 모멘트 $\mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u}$. 그럼 3차는? 연속 볼츠만 방정식의 힘 항은 3차 모멘트도 분명한 값을 갖는다. 그 값을 격자 위에서도 맞춰 주면 정확해질 것 같다.

숫자를 찍어 보면 결과는 반대였다. 3차 모멘트를 운동론이 요구하는 값으로 채우자 오차가 0.053%에서 0.150%로 늘었다. 속도 슬라이더 $u_x$, $u_y$를 올리면 맞던 모멘트 줄이 갈라진다. 체크박스는 노란 줄만 닫는다.

속도를 0 쪽으로 당기면 줄이 다시 맞는다. 빨간 줄은 속도를 키워도 체크박스로 닫히지 않는다.

연속 볼츠만이 내놓는 힘의 3차 모멘트#

체적력 $\mathbf{F}$가 있는 볼츠만 방정식의 힘 항은 속도 공간의 미분이다.

$\xi$는 분자속도, $\rho$는 밀도다. 3차 모멘트는 부분적분 한 번으로 닫힌다.

$M_{\alpha\beta}=\int f^{eq}\xi_\alpha\xi_\beta\,d\boldsymbol{\xi} = \rho(u_\alpha u_\beta + c_s^2\delta_{\alpha\beta})$이다. 그래서 요구되는 3차 모멘트는 이렇게 적힌다.

$c_s^2$ 항과 $uuF$ 항이 나란히 들어 있다. 앞의 것은 속도와 무관하고, 뒤의 것은 속도의 2차다.

Guo forcing이 실제로 공급하는 값#

격자에 올리는 힘 텀은 보통 Guo 형태로 쓴다. 반 스텝 인자를 떼고 괄호 안만 $\hat{S}_i$라 부르자.

$w_i$는 격자 가중치, $\mathbf{c}_i$는 이산 속도다. 실제 $S_i$는 여기에 $(1-1/2\tau)$가 곱해진다. 그 인자가 왜 필요한지는 이산 LBE의 반 스텝에서 따로 다뤘다.

이 $\hat{S}_i$의 3차 모멘트를 D2Q9의 4차 등방성으로 계산하면 다음이 남는다.

$c_s^2$ 항은 정확히 맞는다. $uuF$ 항은 하나도 없다. 두 번째 항은 5차 모멘트를 만들어 D2Q9에서 0이 된다.

앞에서 속도를 올렸을 때 노란 두 줄만 열렸고, 빨간 줄은 잠긴 채로 남았다.

$c_x^3 = c_x$ — 격자가 잠가버린 두 성분#

빨간 줄의 정체는 격자 자체다. D2Q9의 속도 성분은 $-1, 0, 1$뿐이다. 그러면 모든 방향에서 $c_{ix}^3 = c_{ix}$가 성립한다. 두 벡터가 같은 벡터다.

$\Phi_{xxx}$는 고를 수 있는 값이 아니다. 1차 모멘트를 $F_x$로 정한 순간 함께 정해진다. 요구값은 $3F_x(u_x^2+c_s^2) = F_x + 3F_x u_x^2$이므로, $3F_x u_x^2$만큼이 영구 결손이다.

반면 $\Phi_{xxy}$와 $\Phi_{xyy}$는 독립이다. 9개 자유도 중 0~2차가 6개를 쓰고 3개가 남는다. 그래서 4개 성분 중 둘은 자유롭고 둘은 잠겨 있다. 텐서의 절반만 고칠 수 있다는 뜻이다.

질량도 운동량도 나르지 않는 항이 운동량을 옮긴다#

$\Phi_{xyy}$만 $d$만큼 바꾸는 보정 벡터 $v_i$를 만들 수 있다. 0차·1차·2차 모멘트가 모두 0이므로 밀도, 운동량, 응력을 그 자리에서는 건드리지 않는다. D2Q9에서 그 벡터는 이렇게 생겼다.

그런데 스트리밍을 한 번 하면 이야기가 달라진다. 각 성분이 서로 다른 이웃으로 날아간다. $x$방향 운동량만 모아보면 다음이 남는다.

2차 중앙차분이다. $d$가 균일하면 정확히 0이고, $d$가 진동하면 파수의 제곱만큼 커진다. 파형이 균일이면 칸마다의 변화는 사라지고, 사인파에서 파장을 줄이면 그 변화가 파수 제곱을 따라 자란다. 칸의 합은 0으로 남는다.

힘 패턴을 고정하고 횡류를 태운 채널#

검증에는 정확해가 있는 문제가 필요하다. 주기 영역에 $x$방향 사인 체적력을 고정하고, $y$방향으로 균일한 횡류 $V$를 태운다. 다공질 벽으로 한쪽에서 불어넣고 반대쪽에서 빨아내는 상황이다.

$\nu$는 동점성, $k$는 파수, $a_0$는 힘 진폭이다. 단일 모드라 닫힌 해가 나온다.

횡류는 두 가지를 한다. 진폭을 깎고, 위상을 $\arctan\!\left(V/\nu k\right)$만큼 늦춘다. 동시에 $u_y = V$가 살아 있으니 $d = F_x u_y^2$가 공간에서 진동한다. 앞 절의 누출이 켜지는 조건이다.

네 조합을 같은 격자에 걸었다#

$16$노드 주기 격자, $\tau = 1$, $a_0 = 4\times10^{-3}$으로 네 가지를 돌렸다. 표준 Guo, 힘의 3차 채움, 평형분포의 3차 채움, 둘 다.

import numpy as np
 
CX  = np.array([0, 1, 0,-1, 0, 1,-1,-1, 1], float)
CY  = np.array([0, 0, 1, 0,-1, 1, 1,-1,-1], float)
WT  = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
CS2 = 1.0/3.0
 
def hermite3_vectors():
    """0~2차 모멘트는 0, 3차 xxy / xyy만 1로 만드는 보정 벡터."""
    rows = np.vstack([np.ones(9), CX, CY, CX*CX, CX*CY, CY*CY, CX*CX*CY, CX*CY*CY])
    pinv = np.linalg.pinv(rows)
    return pinv @ np.eye(8)[6], pinv @ np.eye(8)[7]
 
VXXY, VXYY = hermite3_vectors()
 
def feq_second_order(rho, ux, uy):
    cu = CX[:, None, None]*ux + CY[:, None, None]*uy
    return WT[:, None, None]*rho*(1 + cu/CS2 + cu*cu/(2*CS2**2)
                                  - (ux*ux + uy*uy)/(2*CS2))
 
def force_populations(ux, uy, fx, fy, tau):
    """Guo forcing: 1차 모멘트 F, 2차 모멘트 uF+Fu."""
    cF  = CX[:, None, None]*fx + CY[:, None, None]*fy
    gxx, gyy, gxy = ux*fx, uy*fy, 0.5*(ux*fy + uy*fx)
    quad = (gxx*(CX[:, None, None]**2 - CS2) + gyy*(CY[:, None, None]**2 - CS2)
            + 2*gxy*(CX*CY)[:, None, None])
    return WT[:, None, None]*(1 - 1/(2*tau))*(cF/CS2 + quad/CS2**2)
 
def analytic_wave(ny, nu, V, a0):
    """nu u'' - V u' + a0 sin(ky) = 0 의 주기 정상해."""
    k = 2*np.pi/ny
    y = np.arange(ny, dtype=float)
    c1 = a0/(nu*k*k + V*V/nu)
    return c1*np.sin(k*y) - (V*c1/(nu*k))*np.cos(k*y)
 
def run_cross_flow(ny, tau, V, a0, cubic_eq, cubic_force, nstep=6000, tol=1e-12):
    nx = 3
    k  = 2*np.pi/ny
    fx = np.tile(a0*np.sin(k*np.arange(ny)), (nx, 1))
    fy = np.zeros((nx, ny))
    rho = np.ones((nx, ny)); ux = np.zeros((nx, ny)); uy = np.full((nx, ny), V)
    f = feq_second_order(rho, ux, uy)
    half, prev = 1 - 1/(2*tau), ux.copy()
    for s in range(nstep):
        rho = f.sum(0)
        ux  = ((CX[:, None, None]*f).sum(0) + 0.5*fx)/rho
        uy  = ((CY[:, None, None]*f).sum(0) + 0.5*fy)/rho
        feq = feq_second_order(rho, ux, uy)
        if cubic_eq:                      # 평형분포의 3차 결손 rho*u*u*u 보충
            feq = feq + VXXY[:, None, None]*(rho*ux*ux*uy) \
                      + VXYY[:, None, None]*(rho*ux*uy*uy)
        src = force_populations(ux, uy, fx, fy, tau)
        if cubic_force:                   # 힘 텀의 3차 결손 F*u*u 보충
            src = src + half*(VXXY[:, None, None]*(2*fx*ux*uy + fy*ux*ux)
                              + VXYY[:, None, None]*(fx*uy*uy + 2*fy*ux*uy))
        f = f - (f - feq)/tau + src
        for i in range(9):
            f[i] = np.roll(f[i], (int(CX[i]), int(CY[i])), axis=(0, 1))
        if s % 1000 == 999:
            if np.max(np.abs(ux - prev)) < tol:
                break
            prev = ux.copy()
    return ux[0]
 
print("D2Q9 lock check:  max|cx^3 - cx| =", float(np.abs(CX**3 - CX).max()))
for tag, mom in (("xxy", CX*CX*CY), ("xyy", CX*CY*CY), ("xxx", CX**3), ("xy", CX*CY)):
    print(f"  correction vector v_xxy :  sum v_i c_{tag:3s} = {float(mom @ VXXY): .6f}")
 
NY, TAU = 16, 1.0
NU, A0  = CS2*(TAU - 0.5), 0.004
print(f"\nny={NY}  tau={TAU}  nu={NU:.5f}  k={2*np.pi/NY:.4f}  a0={A0}")
print("\n   V    |u|max   Guo       +force    +feq3     +both")
for V in (0.05, 0.10, 0.20):
    ref = analytic_wave(NY, NU, V, A0)
    errs = []
    for ce, cf in ((0, 0), (0, 1), (1, 0), (1, 1)):
        u = run_cross_flow(NY, TAU, V, A0, bool(ce), bool(cf))
        errs.append(100*np.linalg.norm(u - ref)/np.linalg.norm(ref))
    print(f"{V:5.2f}  {np.abs(ref).max():.5f}  " + "  ".join(f"{e:7.4f}%" for e in errs))

D2Q9 lock check:  max|cx^3 - cx| = 0.0
  correction vector v_xxy :  sum v_i c_xxy =  1.000000
  correction vector v_xxy :  sum v_i c_xyy = -0.000000
  correction vector v_xxy :  sum v_i c_xxx =  0.000000
  correction vector v_xxy :  sum v_i c_xy  = -0.000000
 
ny=16  tau=1.0  nu=0.16667  k=0.3927  a0=0.004
 
   V    |u|max   Guo       +force    +feq3     +both
 0.05  0.12258   0.5993%   0.5915%   0.0181%   0.0143%
 0.10  0.08374   1.6422%   1.6207%   0.0295%   0.0375%
 0.20  0.04826   3.7307%   3.6840%   0.0531%   0.1503%

표의 +feq3와 +both 간격은 파장에 따라 커진다. 파형 토글과 파장 슬라이더가 그 누출의 크기 막대다.

uniform이면 초록 막대가 0이다. sine에서 파장을 줄이면 막대가 파수 제곱을 따라 자라고, 칸의 합은 0으로 남는다.

0.053%를 0.150%로 만든 보정#

표를 한 줄씩 읽으면 순서가 보인다. $V=0.20$ 기준이다.

힘의 3차 모멘트를 채운 효과는 기껏 1.3%다. 반면 평형분포의 3차를 채우면 오차가 70배 줄어든다. 그리고 둘을 함께 넣으면 다시 나빠진다. 세 숫자가 같은 이야기를 한다.

이유는 Chapman–Enskog 전개의 차수에 있다. 힘의 2차 모멘트는 $O(\varepsilon)$ 응력식에 바로 들어간다. 힘의 3차 모멘트는 $O(\varepsilon)$ 3차 모멘트식에 들어가고, 거기서 응력으로 한 번, 응력에서 운동량으로 또 한 번 내려온다. 운동량 방정식에는 $O(\varepsilon^3)$, 즉 Burnett 차수에서 도착한다. Navier–Stokes를 목표로 하는 등온 LBM에서 그 항은 맞출 대상이 아니다.

반면 평형분포의 3차 모멘트는 $O(\varepsilon)$ 응력에 직접 들어간다. 그것이 평형분포의 3차 결손이 지배적인 까닭이다.

보정이 오차를 만든 경로#

그럼 힘 3차를 채운 것이 왜 해가 되는가. 앞 절의 누출식이 답이다. $d = F_x u_y^2 = a_0 V^2 \sin(ky)$이므로 스트리밍이 다음 크기의 가짜 체적력을 만든다.

$V=0.2$, $k=0.3927$, $\tau=1$을 넣으면 $a_0$ 대비 $0.154\%$다. 코드에서 재본 두 해의 최대 차이는 $|u|_{\max}$ 대비 $0.156\%$였다. 두 숫자가 맞는다. 운동량 방정식이 원하지도 않은 항을 격자가 성실하게 배달한 것이다.

$V^2 k^2$ 의존성도 그대로 보인다. $V$를 $0.05 \to 0.10 \to 0.20$으로 올렸을 때 두 해의 차이는 $0.010\% \to 0.039\% \to 0.156\%$였다. 4배씩이다.

3차 모멘트가 정말 필요해지는 자리#

그렇다고 $\Phi_{\alpha\beta\gamma}$ 유도가 무의미한 것은 아니다. 필요한 곳이 분명히 있다.

에너지 방정식을 함께 푸는 열·압축성 LBM에서는 3차 모멘트가 열유속이다. 그때 힘의 3차 모멘트는 $O(\varepsilon^2)$ 에너지식에 직접 들어간다. Burnett까지 갈 필요가 없다. 비이상 상태방정식으로 계면을 만드는 다상 LBM도 비슷하다. 힘이 크고 공간적으로 급하다. $k$가 1에 가까워지면 $V^2k^2$가 더는 작지 않다.

$V=0.20$에서 표준 Guo는 3.7307%다. 힘의 3차만 채우면 3.6840%다. 빼면 0.0467%p이고, 3.7307로 나누면 약 1.3%다. 평형분포 3차만 채우면 0.0531%다. 3.7307을 0.0531로 나누면 약 70이다. 둘을 같이 넣으면 0.1503%다. 0.1503을 0.0531로 나누면 약 2.83이라, 표에 적힌 2.8배와 맞다. $V=0.05$에서는 둘 다 넣은 쪽이 0.0143%로, 평형만 채운 0.0181%보다 작다. $V=0.10$에서 0.0375%와 0.0295%로 순서가 뒤집힌다. $V=0.20$의 간격은 0.0972%p다.

누출식에 $V=0.2$, $k=0.3927$, $\tau=1$을 넣으면 $1-1/(2\tau)$는 $1/2$다. 앞의 $1/2$와 곱하면 $1/4$이고, $V^2=0.04$다. $k^2$는 $0.3927$의 제곱으로 약 $0.1542$다. $1/4$ 곱하기 $0.04$ 곱하기 $0.1542$는 $0.001542$, 곧 $a_0$의 $0.154\%$다. 코드가 재본 차이는 $|u|_{\max}$ 대비 $0.156\%$였다. $0.002\%$p 안쪽이다. $V$를 $0.05$에서 $0.10$, $0.20$으로 두 배씩 올리면 $V^2$는 네 배다. 간격 $0.010\%$에서 $0.039\%$, $0.156\%$로 네 배씩 커진 것과 맞다. $|u|_{\max}$는 $0.12258$, $0.08374$, $0.04826$으로 줄어든다. $C_1$의 분모는 $\nu k^2+V^2/\nu$이고 $\nu=0.16667$이라, $V$가 커질수록 진폭이 깎인다.

힘의 절반을 어디에 둘지는 Guo, Shan–Chen, EDM을 나란히 둔 비교가 같은 $(1-1/2\tau)$를 두고 가른다. 목표 방정식이 무엇인지, 그 항이 몇 차에서 도착하는지를 먼저 세어야 한다.

격자에 없는 항을 반만 채우면 안 되는 이유#

이번 실험의 진짜 교훈은 차수보다 앞에 있다. $\Phi_{\alpha\beta\gamma}$는 텐서다. D2Q9는 그중 $xxy$, $xyy$만 허락하고 $xxx$, $yyy$는 잠근다. 텐서의 절반만 채우면 등방성이 깨진 보정이 들어간다. 상쇄되고 있던 오차가 드러날 수 있고, 이번 경우가 그랬다.

격자가 못 나르는 항을 발견하면 선택지는 둘이다. 속도 집합을 늘려 그 모멘트를 확보하거나(D2Q17, D2Q37), 아예 손대지 않는 것이다. 가운데는 없다. 힘 텀에 한 줄을 더하기 전에, 그 줄의 나머지 성분이 격자에 있는지 먼저 세어볼 일이다.