τ=0.6에서 유속이 6분의 1로 나왔다 — Guo forcing에서 반 스텝 속도 F/2를 빼먹으면

Guo 소스항의 (1−1/2τ) 계수는 평형 속도의 F/2 보정과 짝이다. 하나만 넣으면 매 스텝 운동량이 F가 아니라 (1−1/2τ)F만 들어간다.

같은 힘을 걸었는데 왜 유속이 6분의 1일까#

주기 경계 상자에 사인 모양 체적력을 걸었다. 해석해는 교과서 한 줄이다. 그런데 격자 볼츠만(LBM, 분포함수를 충돌·이동시켜 유체를 푸는 방법) 코드가 낸 진폭은 해석해의 0.16배였다. 점성을 바꾸면 비율도 바뀌었다. τ=1에서는 정확히 절반이었다.

원인은 한 줄이다. 거시 속도를 계산할 때 $\mathbf{F}/2$ 를 더하지 않았다. 이 글은 그 한 줄이 왜 정확히 $(1-1/2\tau)$ 배의 오차를 만드는지 보인다. Kolmogorov 흐름(사인 체적력이 만드는 주기 전단 흐름)으로 직접 재고, 내 코드에서 이 버그를 30초 만에 찾는 법으로 끝낸다.

한 스텝 안에서 운동량은 세 번 바뀐다#

LBM의 한 스텝은 충돌과 이동이다. 이동은 운동량을 옮기기만 한다. 운동량이 생기는 곳은 충돌뿐이다. 그래서 균일한 유체 한 점의 장부만 보면 된다.

\mathbf{m}' = \mathbf{m} - \frac{1}{\tau}\left(\mathbf{m} - \rho\mathbf{u}\right) + \left(1-\frac{1}{2\tau}\right)\mathbf{F}

$\mathbf{m}=\sum_i f_i\mathbf{c}_i$ 는 충돌 전 운동량, $\mathbf{m}'$ 은 충돌 후, $\rho\mathbf{u}$ 는 평형 분포에 넣은 운동량, $\tau$ 는 이완 시간, $\mathbf{F}$ 는 격자 단위 체적력이다.

둘째 항은 이완 항이다. 평형 속도가 충돌 전 속도와 같으면 0이다. 셋째 항은 Guo 소스항의 1차 모멘트다. 이제 $\rho\mathbf{u}$ 에 무엇을 넣는지가 전부를 정한다.

\rho\mathbf{u} = \mathbf{m} + \tfrac{1}{2}\mathbf{F} \;\Rightarrow\; \mathbf{m}' = \mathbf{m} + \frac{\mathbf{F}}{2\tau} + \mathbf{F} - \frac{\mathbf{F}}{2\tau} = \mathbf{m} + \mathbf{F}

\rho\mathbf{u} = \mathbf{m} \;\Rightarrow\; \mathbf{m}' = \mathbf{m} + \left(1-\frac{1}{2\tau}\right)\mathbf{F}

첫 줄은 이완 항이 $\mathbf{F}/(2\tau)$ 를 메워 정확히 $\mathbf{F}$ 가 들어간다는 뜻이다. 둘째 줄은 그 메움이 사라진 경우다. 소스항의 $(1-1/2\tau)$ 계수만 남는다.

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

relaxation time τ = 0.60 ν = (τ−½)/3 = 0.0333

step 0 · ρu/(F·t): with F/2 0.0000 · without 0.0000 · theory 1−1/(2τ) = 0.1667

τ 슬라이더를 0.51 쪽으로 끌어 보자. 빨간 선의 기울기가 거의 0이 되고, 오른쪽 막대에서 진한 이완 조각이 통째로 빠진 것이 보인다.

반 스텝 속도: 충돌 전과 후의 한가운데#

왜 하필 $\mathbf{F}/2$ 일까. 충돌 전 운동량은 $\mathbf{m}$ , 힘을 다 받은 뒤는 $\mathbf{m}+\mathbf{F}$ 다. 물리적인 속도는 그 시간 구간의 평균, 즉 가운데에 있다. 사다리꼴 적분(구간 양 끝값의 평균으로 적분하는 2차 정확도 공식)과 같은 이유다.

Guo 등(2002)은 소스항을 다음처럼 두었다.

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

$w_i$ 는 격자 가중치, $\mathbf{c}_i$ 는 이산 속도, $c_s^2=1/3$ 은 격자 음속의 제곱이다. $\sum_i S_i\mathbf{c}_i = \mathbf{F}$ 이고 $\sum_i S_i = 0$ 이다.

$(1-1/2\tau)$ 계수는 Chapman–Enskog 전개(분포함수를 Knudsen 수로 전개해 거시 방정식을 끌어내는 방법)에서 나온다. 이산 시간의 Taylor 전개가 소스항에 $-\tfrac{1}{2}S_i$ 꼴의 찌꺼기를 남기기 때문이다. 이 찌꺼기는 반 스텝 점성 보정과 같은 뿌리다. 그리고 이 계수를 지우려면 평형 속도에 $\mathbf{F}/2$ 를 넣어야 1차 모멘트가 맞는다. 둘은 처음부터 한 쌍으로 유도된다.

코드: Kolmogorov 흐름, 세 가지 속도 정의#

D2Q9 격자를 $1\times N$ 주기 기둥으로 만든다. 힘은 $F_x = F_0\sin(ky)$ , $k=2\pi/N$ 이다. 정상 상태 해석해는 다음과 같다.

u_x(y) = \frac{F_0}{\nu k^2}\sin(ky),\qquad \nu = \frac{\tau - 1/2}{3}

$\nu$ 는 격자 단위 동점성계수다. 벽이 없으니 바운스백 경계 오차가 섞이지 않는다. 힘 처리만 따로 잴 수 있다.

세 경우를 돌린다. guo는 평형과 출력 모두 $\mathbf{F}/2$ 를 넣는다. nohalf는 Guo 소스항은 그대로 두고 속도에서 $\mathbf{F}/2$ 를 뺀다. rawout은 계산은 올바르게 하고 출력할 때만 $\sum f_i\mathbf{c}_i/\rho$ 를 쓴다. 진폭은 $\sin(ky)$ 푸리에 계수로 재고 해석해로 나눈다.

import math
 
# D2Q9 lattice: velocities, weights, sound speed squared
C = [(0,0),(1,0),(0,1),(-1,0),(0,-1),(1,1),(-1,1),(-1,-1),(1,-1)]
W = [4/9,1/9,1/9,1/9,1/9,1/36,1/36,1/36,1/36]
CS2 = 1/3
 
def f_equil(rho, ux):
    # uy = 0 everywhere: the flow is u_x(y) only
    out = []
    for (cx, cy), wi in zip(C, W):
        cu = cx*ux
        out.append(wi*rho*(1 + cu/CS2 + 0.5*(cu/CS2)**2 - 0.5*ux*ux/CS2))
    return out
 
def guo_term(ux, Fx, tau):
    out = []
    for (cx, cy), wi in zip(C, W):
        cu = cx*ux
        out.append((1 - 0.5/tau)*wi*((cx - ux)/CS2 + cu*cx/CS2**2)*Fx)
    return out
 
def kolmogorov_amp(tau, mode, N=16, F0=1e-6):
    """mode 'guo'    : u = (sum f c + F/2)/rho in equilibrium and output
       mode 'nohalf' : Guo source, but u = sum f c / rho everywhere
       mode 'rawout' : correct update, but report sum f c / rho"""
    k = 2*math.pi/N
    Fx = [F0*math.sin(k*(j + 0.5)) for j in range(N)]
    nu = (tau - 0.5)/3
    steps = int(12/(nu*k*k))                    # ~12 viscous decay times
    f = [f_equil(1.0, 0.0) for _ in range(N)]
    for _ in range(steps):
        post = []
        for j in range(N):
            fj = f[j]
            rho = sum(fj)
            mx = sum(fi*cx for fi, (cx, _) in zip(fj, C))
            ux = mx/rho if mode == 'nohalf' else (mx + 0.5*Fx[j])/rho
            feq = f_equil(rho, ux)
            src = guo_term(ux, Fx[j], tau)
            post.append([fi - (fi - fe)/tau + s for fi, fe, s in zip(fj, feq, src)])
        # streaming: only the y component moves anything (one cell in x, periodic)
        f = [[post[(j - C[i][1]) % N][i] for i in range(9)] for j in range(N)]
    amp = 0.0
    for j in range(N):
        rho = sum(f[j]); mx = sum(fi*cx for fi, (cx, _) in zip(f[j], C))
        u = mx/rho if mode in ('nohalf', 'rawout') else (mx + 0.5*Fx[j])/rho
        amp += 2*u*math.sin(k*(j + 0.5))/N
    return amp/(F0/(nu*k*k)), steps
 
print(" tau  steps   guo     nohalf  rawout   nohalf/guo  1-1/(2tau)")
for tau in [0.55, 0.6, 0.8, 1.0, 1.5]:
    r1, n = kolmogorov_amp(tau, 'guo')
    r2, _ = kolmogorov_amp(tau, 'nohalf')
    r3, _ = kolmogorov_amp(tau, 'rawout')
    print(f"{tau:4.2f} {n:6d}  {r1:.4f}  {r2:.4f}  {r3:.4f}   {r2/r1:.4f}      {1-0.5/tau:.4f}")

 tau  steps   guo     nohalf  rawout   nohalf/guo  1-1/(2tau)
0.55   4668  0.9746  0.0885  0.9734   0.0908      0.0909
0.60   2334  0.9754  0.1621  0.9729   0.1662      0.1667
0.80    778  0.9836  0.3660  0.9759   0.3721      0.3750
1.00    466  1.0001  0.4936  0.9872   0.4936      0.5000
1.50    233  1.0772  0.7010  1.0515   0.6508      0.6667

위 코드를 브라우저에서 직접 돌려보자.

τ = 0.60 N = 16 steps/frame = 20

step 0/0 · amplitude/exact: guo 0.0000 · nohalf 0.0000 · rawout 0.0000 · nohalf/guo 0.0000 (1−1/2τ = 0.1667)

τ를 0.6에 두면 빨간 nohalf 진폭이 0.16 근처에서 멈추고 빨간 점선 $1-1/2\tau$ 에 붙는다. N을 32 이상으로 올리면 초록 guo가 1에 더 가까워지지만 nohalf/guo 비율은 그대로다.

0.166과 0.973이 말하는 것#

nohalf/guo 열이 $1-1/2\tau$ 와 겹친다. τ=0.55에서 0.0908 대 0.0909, τ=0.6에서 0.1662 대 0.1667이다. 운동량 장부의 결론이 실제 유동에서도 그대로 나온다. 정상 상태에서는 점성 확산이 매 스텝 들어오는 운동량과 맞선다. 들어오는 양이 $(1-1/2\tau)$ 배면 진폭도 그만큼이다.

τ가 커질수록 비율이 이론값에서 조금 벗어난다. τ=1.5에서 0.6508 대 0.6667이다. $N=16$ 이면 $k\approx0.39$ 로 파장당 격자가 적다. 이 영역에서는 $k^2$ 차수의 격자 절단 오차가 섞인다. guo 자체도 1.0772로 1에서 떨어져 있다. 같은 이유다.

rawout은 guo와 거의 같다. 출력 속도가 $F/(2\rho)$ 만큼 늦을 뿐이다. 상대 오차는 $\nu k^2/2$ 정도라 τ=0.55에서 0.1%다. 진단할 때 혼동하기 쉬운 지점이다. 출력 쪽의 $\mathbf{F}/2$ 누락은 작다. 평형 쪽의 누락은 크다.

가장 위험한 곳은 고 Reynolds 수 계산이다. 점성을 낮추려고 τ를 0.5에 붙이면 유효 힘도 0으로 간다. 유속이 안 나와서 힘을 키우면 숨은 계수가 모든 결과에 남는다.

수도코드 어디에 F/2가 들어가는가#

흔한 LBM 루프 순서는 초기화 → 거시량 계산 → 평형 분포 → 힘 → 충돌 → 이동이다. 버그는 이 순서 때문에 생긴다. 거시량 단계에서 $\mathbf{u}$ 를 구할 때는 아직 힘이 계산되지 않았다. 그래서 $\mathbf{u}=\sum f_i\mathbf{c}_i/\rho$ 로 끝내기 쉽다.

해법은 순서를 조금 바꾸는 것이다.

$\rho = \sum_i f_i$ 를 먼저 구한다.
힘 $\mathbf{F}$ 를 계산한다. 중력처럼 $\rho$ 만 필요하면 여기서 끝난다.
$\rho\mathbf{u} = \sum_i f_i\mathbf{c}_i + \mathbf{F}/2$ 로 속도를 구한다.
같은 $\mathbf{u}$ 로 $f^{eq}$ 와 $S_i$ 를 만든다. 출력도 이 $\mathbf{u}$ 를 쓴다.

힘이 속도에 의존하면(예: 항력) 3단계가 암묵식이 된다. 선형 항력이면 대수적으로 풀린다. 그 밖의 경우는 고정점 반복 한두 번이면 충분하다.

$(1-1/2\tau)$ 를 지워 "고치는" 것도 답이 아니다. 1차 모멘트는 맞지만 2차 모멘트에 $\tau$ 의존 응력이 남는다. 이 문제는 힘 텀의 3차 모멘트 글과 Shan–Chen 속도 이동 forcing 글에서 다뤘다.

내 LBM 루프에서 F/2를 찾는 한 스텝 테스트#

주기 상자를 만들고 $\rho=1$ , $\mathbf{u}=0$ 으로 채운다. 균일한 힘 $\mathbf{F}$ 를 걸고 딱 한 스텝을 돌린다. 그리고 $\sum f_i\mathbf{c}_i$ 의 변화량을 $\mathbf{F}$ 로 나눈다.

결과가 1이면 짝이 맞다. $1-1/2\tau$ 면 평형 속도에서 $\mathbf{F}/2$ 가 빠졌다. τ를 두 개 골라 돌리면 확실하다. 비율이 τ를 따라 움직이면 버그다. 검증 케이스를 돌리기 전에 이 한 스텝부터 확인하자.