프레임을 0.2로 밀었더니 점성이 5.97% 줄었다 — LBM 평형분포의 3차 모멘트 결손

LBM의 속도 상한은 안정성이 아니라 평형분포의 3차 모멘트가 먼저 정한다.

정지한 물과 흐르는 물에서 점성이 같아야 하는가#

같은 유체를 두 번 잰다고 하자. 한 번은 가만히 있는 상자 안에서, 한 번은 그 상자를 일정한 속도로 밀면서. 점성계수는 두 번 다 같아야 한다. 갈릴레이 불변성(등속으로 움직이는 두 좌표계에서 물리 법칙이 같다는 성질)이 그렇게 요구한다.

D2Q9 격자 볼츠만(LBM)으로 이 실험을 그대로 해보면 두 번째 값이 5.97% 작게 나온다. 격자를 촘촘히 해도, 완화시간 $\tau$ 를 바꿔도 이 비율은 거의 그대로다. 이 글은 그 5.97%를 평형분포의 모멘트 한 줄까지 따라간다. 결론부터 말하면 안정성 문제가 아니라 대수 문제다. 격자 속도 집합이 특정 항을 애초에 만들 수 없다.

격자 위의 한 줄을 Taylor로 펼치면#

LBM이 푸는 식은 한 줄이다.

f_i(\mathbf{x} + \mathbf{c}_i \Delta t,\; t + \Delta t) - f_i(\mathbf{x}, t) = -\frac{\Delta t}{\tau}\left[ f_i - f_i^{eq} \right]

$f_i$ 는 격자 속도 $\mathbf{c}_i$ 를 타고 가는 분포함수, $f_i^{eq}$ 는 그 자리의 평형분포, $\tau$ 는 완화시간이다.

왼쪽을 $\Delta t$ 로 Taylor 전개하면 1차항 말고 2차항이 남는다.

\Delta t \left( \partial_t + \mathbf{c}_i \cdot \nabla \right) f_i + \frac{\Delta t^2}{2} \left( \partial_t + \mathbf{c}_i \cdot \nabla \right)^2 f_i + O(\Delta t^3) = -\frac{\Delta t}{\tau} \left( f_i - f_i^{eq} \right)

이 2차항이 Chapman–Enskog 전개(Knudsen 수의 멱급수로 분포함수를 쪼개는 다중스케일 기법)에서 그 유명한 $\tau - 1/2$ 를 만든다. 그리고 그다음 단계에서, 얻어지는 연속체 방정식이 Navier–Stokes가 되려면 평형분포가 모멘트 네 개를 정확히 맞춰야 한다는 조건이 나온다.

\sum_i f_i^{eq} = \rho, \qquad \sum_i f_i^{eq} c_{i\alpha} = \rho u_\alpha, \qquad \sum_i f_i^{eq} c_{i\alpha} c_{i\beta} = \rho u_\alpha u_\beta + \rho c_s^2 \delta_{\alpha\beta}

\sum_i f_i^{eq} c_{i\alpha} c_{i\beta} c_{i\gamma} = \rho c_s^2 \left( u_\alpha \delta_{\beta\gamma} + u_\beta \delta_{\gamma\alpha} + u_\gamma \delta_{\alpha\beta} \right) + \rho u_\alpha u_\beta u_\gamma

앞의 세 개는 질량·운동량·압력을 정한다. 네 번째, 3차 모멘트는 점성 응력의 시간 미분을 대체할 때 쓰인다. 여기가 어긋나면 연속방정식도 Euler 방정식도 멀쩡한데 점성항만 오염된다.

아래 시뮬레이션에서 3차 모멘트가 맞는지 직접 확인해보자.

velocity u : 0.180 (Ma = 0.31)

Σ f c³ = 0.180000  |  Maxwell = 0.185832  |  gap = 0.005832  |  gap / u³ = 1.0000

속도 슬라이더를 밀면 노란 곡선(Maxwell 분포가 요구하는 값)과 파란 직선(D2Q9가 실제로 내놓는 값)이 벌어진다. sweep u로 자동 왕복시키면 벌어지는 폭이 $u$ 의 홀수 제곱을 따라간다는 것이 보인다. speeds ±2 버튼을 켜면 막대가 다섯 개로 늘고 간격이 0으로 닫힌다.

$c_x^3 = c_x$ — 격자가 만들 수 없는 항#

D2Q9의 $x$ 방향 성분은 $-1, 0, 1$ 세 값뿐이다. 이 세 수는 모두 세제곱해도 자기 자신이다. 그러므로 임의의 분포 $f_i$ 에 대해

\sum_i f_i c_{ix}^3 = \sum_i f_i c_{ix} = \rho u_x

가 항상 성립한다. 평형분포를 어떻게 설계하든, 심지어 난수를 채워 넣어도 마찬가지다. 3차 모멘트는 이미 1차 모멘트에 묶여 있다.

한편 Maxwell 분포가 요구하는 값은 $\rho(u_x^3 + 3 c_s^2 u_x)$ 이고, D2Q9는 $c_s^2 = 1/3$ 이므로 $3c_s^2 = 1$ 이다. 결국 요구값은 $\rho(u_x + u_x^3)$ 이 된다. 두 값의 차이는 정확히 $\rho u_x^3$ 하나다. 오차 항이 근사적으로 그 정도라는 말이 아니라, 나머지 없이 그 항 하나다.

결손항이 운동량 방정식에 남기는 얼굴#

빠진 3차 모멘트는 점성 응력에 그대로 더해진다. Chapman–Enskog 전개를 끝까지 밀면 응력에 항이 하나 붙는다.

\sigma_{\alpha\beta}^{\text{err}} = -\left( \tau - \tfrac{1}{2} \right) \Delta t \, \partial_\gamma \left( \rho\, u_\alpha u_\beta u_\gamma \right)

이제 평균류 $U$ 가 $x$ 방향으로 흐르고 그 위에 작은 교란 $u'$ 이 얹혀 있다고 하자. 위 항에서 $u'$ 에 1차인 부분만 남기면 $-(\tau - 1/2)\Delta t\, U^2 \partial_x^2 (\rho u'_\alpha)$ 가 되고, 이건 점성항과 형태가 같다. 즉 $x$ 방향 점성계수가 통째로 바뀐다.

\nu_{xx}^{\text{eff}} = \left( \tau - \tfrac{1}{2} \right) \Delta t \left( c_s^2 - U^2 \right)

상대 오차는 $-U^2 / c_s^2$ 다. 여기서 $\tau$ 가 분자와 분모에서 동시에 빠진다. 격자를 아무리 정교하게 완화시켜도 상대 오차는 그대로라는 뜻이다. 그리고 오차는 속도의 제곱에 붙는다. 마하수로 쓰면 $-\mathrm{Ma}^2$ 다.

Python으로 같은 소용돌이를 두 프레임에서 쟀다#

Taylor–Green 소용돌이를 64×64 주기 격자에 올리고, 여기에 등속 $U_0$ 를 더한다. 진폭의 감쇠율을 로그 기울기로 재면 점성계수가 역산된다. $U_0$ 만 바꾸고 나머지는 전부 고정한다.

import numpy as np
 
CX  = np.array([0, 1, 0, -1,  0, 1, -1, -1,  1], dtype=float)
CY  = np.array([0, 0, 1,  0, -1, 1,  1, -1, -1], dtype=float)
W   = 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 f_equilibrium(rho, ux, uy):
    # 2차까지 전개한 표준 평형분포
    cu = CX[:, None, None] * ux + CY[:, None, None] * uy
    usq = ux * ux + uy * uy
    return W[:, None, None] * rho * (1 + cu / CS2 + cu * cu / (2 * CS2**2) - usq / (2 * CS2))
 
print("[1] third moment audit:  sum f_i c_ix^3   vs   Maxwell rho*(u^3 + 3*cs2*u)")
print("   u        lattice      Maxwell        gap        gap/u^3")
for u in (0.05, 0.10, 0.20, 0.30):
    feq = f_equilibrium(np.ones((1, 1)), np.full((1, 1), u), np.zeros((1, 1)))
    lat = float((feq * (CX**3)[:, None, None]).sum())
    exact = u**3 + 3 * CS2 * u
    print(f"  {u:4.2f}  {lat:11.8f}  {exact:11.8f}  {exact-lat:11.8f}   {(exact-lat)/u**3:8.5f}")
 
rng = np.random.default_rng(7)
frand = rng.random((9, 1, 1))
d = float((frand * (CX**3)[:, None, None]).sum() - (frand * CX[:, None, None]).sum())
print(f"  any random f:  sum f c_x^3 - sum f c_x = {d:.3e}   (c_x^3 = c_x, so always 0)")
 
def taylor_green_run(U0, tau=0.8, N=64, steps=900, amp=0.04, warmup=300, sample=25):
    # 같은 소용돌이에 등속 U0를 얹고 감쇠율로 점성을 되잰다
    k = 2 * np.pi / N
    x = np.arange(N)[:, None] * np.ones(N)[None, :]
    y = np.ones(N)[:, None] * np.arange(N)[None, :]
    ux = U0 - amp * np.cos(k * x) * np.sin(k * y)
    uy = amp * np.sin(k * x) * np.cos(k * y)
    rho = np.ones((N, N))
    f = f_equilibrium(rho, ux, uy)
    ts, logs = [], []
    for n in range(steps + 1):
        rho = f.sum(axis=0)
        ux = (f * CX[:, None, None]).sum(axis=0) / rho
        uy = (f * CY[:, None, None]).sum(axis=0) / rho
        if n % sample == 0:
            up, vp = ux - ux.mean(), uy - uy.mean()   # 평균류를 뺀 변동 성분
            ts.append(n)
            logs.append(0.5 * np.log(2 * np.mean(up * up + vp * vp)))
        f += (f_equilibrium(rho, ux, uy) - f) / tau   # BGK 충돌
        for i in range(9):                            # 스트리밍
            f[i] = np.roll(np.roll(f[i], int(CX[i]), axis=0), int(CY[i]), axis=1)
    ts, logs = np.array(ts, dtype=float), np.array(logs)
    m = ts >= warmup
    slope = np.polyfit(ts[m], logs[m], 1)[0]
    return -slope / (2 * k * k)
 
print("\n[2] same vortex measured in a uniformly moving frame  (tau=0.8, nu_theory=0.100000)")
nu_th = CS2 * (0.8 - 0.5)
print("   U0      nu_eff      rel.err     rel.err/U0^2")
for U0 in (0.00, 0.05, 0.10, 0.15, 0.20):
    nu = taylor_green_run(U0)
    e = nu / nu_th - 1
    tail = f"{e/U0**2:10.4f}" if U0 > 0 else "         -"
    print(f"  {U0:4.2f}  {nu:10.7f}  {e:+9.4%}  {tail}")
 
print("\n[3] does the error depend on tau?  (U0=0.15 fixed)")
print("   tau     nu_theory    nu_eff      rel.err")
for tau in (0.6, 0.8, 1.0):
    nu = taylor_green_run(0.15, tau=tau)
    th = CS2 * (tau - 0.5)
    print(f"  {tau:4.2f}  {th:10.7f}  {nu:10.7f}  {nu/th-1:+9.4%}")

[1] third moment audit:  sum f_i c_ix^3   vs   Maxwell rho*(u^3 + 3*cs2*u)
   u        lattice      Maxwell        gap        gap/u^3
  0.05   0.05000000   0.05012500   0.00012500    1.00000
  0.10   0.10000000   0.10100000   0.00100000    1.00000
  0.20   0.20000000   0.20800000   0.00800000    1.00000
  0.30   0.30000000   0.32700000   0.02700000    1.00000
  any random f:  sum f c_x^3 - sum f c_x = 0.000e+00   (c_x^3 = c_x, so always 0)
 
[2] same vortex measured in a uniformly moving frame  (tau=0.8, nu_theory=0.100000)
   U0      nu_eff      rel.err     rel.err/U0^2
  0.00   0.1000175   +0.0175%           -
  0.05   0.0996433   -0.3567%     -1.4268
  0.10   0.0985207   -1.4793%     -1.4793
  0.15   0.0966495   -3.3505%     -1.4891
  0.20   0.0940296   -5.9704%     -1.4926
 
[3] does the error depend on tau?  (U0=0.15 fixed)
   tau     nu_theory    nu_eff      rel.err
  0.60   0.0333333   0.0321975   -3.4076%
  0.80   0.1000000   0.0966495   -3.3505%
  1.00   0.1666667   0.1611502   -3.3099%

오차는 $\tau$ 가 아니라 $U^2$ 에 붙어 있다#

[1]의 마지막 열이 전부 1.00000이다. 결손이 $\rho u^3$ 하나라는 주장이 소수점 다섯 자리까지 맞는다는 뜻이다.

[2]에서 정지 프레임의 오차는 +0.0175%다. 측정 방법 자체의 잡음이 그 정도라는 기준선이다. $U_0$ 를 올리면 오차가 $-0.36\%, -1.48\%, -3.35\%, -5.97\%$ 로 커진다. 마지막 열은 오차를 $U_0^2$ 로 나눈 값이고, $U_0$ 가 작아질수록 $-1.5$ 에 붙는다.

이 $-1.5$ 가 앞 절의 예측이다. 상대 오차 $-U^2/c_s^2$ 는 $x$ 성분에만 걸리는데, Taylor–Green 모드는 $k_x = k_y$ 인 대각 모드라 감쇠율이 두 방향의 평균으로 나온다. 그래서 절반, 즉 $-U^2/(2c_s^2) = -1.5 U^2$ 다. $U_0 = 0.2$ 에서 예측은 $-6.00\%$ , 측정은 $-5.97\%$ 다.

[3]이 더 중요하다. $\tau$ 를 0.6에서 1.0까지 올려 점성을 다섯 배로 키워도 상대 오차는 $-3.41\%$ 에서 $-3.31\%$ 로 거의 움직이지 않는다. 점성을 키워 오차를 묻어버릴 수 없다는 뜻이다.

frame velocity U : 0.18 (Ma = 0.31)relaxation τ : 0.80

ν = 0.10000  |  νeff = 0.09514  |  error -4.86 %  |  step 0  |  amplitude ratio 1.000

frame velocity U를 올리면 오른쪽 소용돌이가 흘러가면서 왼쪽보다 천천히 흐려진다. 같은 유체인데 오른쪽이 덜 사라진다. $\tau$ 슬라이더를 흔들어보면 아래쪽 error 값이 거의 고정이라는 것도 같이 보인다.

고치는 방법 세 가지와 각각의 청구서#

속도를 낮춘다. 오차가 $U^2$ 이므로 격자 속도를 절반으로 줄이면 오차는 4분의 1이 된다. 대신 같은 물리 시간을 재현하는 데 필요한 스텝 수가 늘어난다. LBM에서 "마하수를 0.1 아래로"라는 경험칙은 대개 안정성 이야기로 소개되지만, 실제로 먼저 걸리는 것은 이 결손항이다.

보정항을 더한다. $\partial_\gamma(\rho u_\alpha u_\beta u_\gamma)$ 를 유한차분으로 계산해서 충돌 단계에 반대 부호로 넣는다. 추가 비용은 기울기 한 번이고, 격자 볼츠만의 장점인 완전 국소 충돌이 조금 깨진다. 그래도 회전 좌표계처럼 평균류가 큰 문제에서는 가장 싼 선택이다.

속도 집합을 넓힌다. $\pm 2$ 를 포함하는 다속도 격자는 $c^3 \neq c$ 이므로 3차 모멘트를 정확히 맞출 수 있다. 위 viz의 speeds ±2 버튼이 그 계산이다. 대가는 스텐실이 넓어지는 것, 메모리가 늘어나는 것, 경계 처리가 복잡해지는 것이다. 일부 속도의 가중치가 음수가 될 수 있어 안정성도 다시 봐야 한다.

이 결손이 실제 계산에서 얼굴을 내미는 자리#

평균류가 큰 곳이 전부 후보다. 회전 기계의 회전 좌표계, 슬라이딩 메시, 빠른 균일류 위에 얹힌 난류, 이동하는 물체를 따라가는 프레임. 정지 유동 벤치마크에서는 멀쩡하던 코드가 이런 문제에서만 유효 레이놀즈 수가 어긋나는 경우가 여기에 해당한다.

LBM 격자 세밀화의 비평형 재스케일에서 레벨마다 $\tau$ 를 다시 잡는 이야기를 했는데, 이 결손은 그 조정으로 사라지지 않는다. 상대 오차가 $\tau$ 와 무관하기 때문이다. 열 LBM에 잠겨 있는 두 상수에서 본 것과 같은 구조이기도 하다. 물성으로 넣지 않은 값이 격자 구조에서 자동으로 결정된다.

진단은 싸다. 같은 문제를 평균류를 얹어서 한 번 더 돌리고 감쇠율이나 항력계수를 비교한다. 차이가 $U^2$ 에 비례해서 커지면 이 결손항이다. STL 복셀화의 홀짝 판정처럼 경계 쪽을 의심할 필요가 없다.

다음에 프레임마다 점성이 다르게 나오면#

격자 볼츠만에서 속도 상한을 정하는 것은 안정성만이 아니다. 평형분포가 맞출 수 있는 모멘트의 차수가 먼저 상한을 정한다. D2Q9는 3차에서 $\rho u^3$ 을 놓치고, 그 대가는 $-\mathrm{Ma}^2$ 의 점성 오차로 청구된다.

수치가 이상할 때 격자를 촘촘히 하거나 $\tau$ 를 만지는 것이 먼저 떠오른다. 이 오차는 둘 다로 줄지 않는다. 줄어드는 것은 $u/c_s$ 를 낮추거나, 결손항을 직접 되돌려 넣거나, 속도 집합을 넓힐 때뿐이다.