프레임을 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 방정식도 멀쩡한데 점성항만 오염된다.

속도 슬라이더를 0에서 밀어 올리면 노란 Maxwell 곡선이 파란 격자 직선에서 갈라진다.

velocity u : 0.180 (Ma = 0.31)

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

갈라진 폭이 $u$ 의 홀수 제곱을 따르는지는 sweep 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$ 다.

숫자로 한 번 밟아 두자. $c_s^2 = 1/3$ , $U = 0.2$ 이면 $U^2 = 0.04$ 이고 $U^2/c_s^2 = 0.12$ 다. $x$ 방향만 보면 점성이 12% 준다. Taylor–Green은 $k_x = k_y$ 라 감쇠율이 두 방향의 평균이다. 절반인 6%가 예측이다. 표의 $U_0 = 0.20$ 은 $\nu = 0.0940296$ 으로, 이론값 $0.100000$ 에서 $-5.9704\%$ 다. 6.00과 5.97 사이는 0.03퍼센트포인트다.

$\tau$ 를 흔들어도 같은 자리다. $\tau = 0.6$ 의 이론 점성은 $c_s^2 \times 0.1 = 0.0333333$ 이고, $\tau = 1.0$ 은 $c_s^2 \times 0.5 = 0.1666667$ 이다. 다섯 배다. 그런데 $U_0 = 0.15$ 의 상대 오차는 $-3.4076\%$ 에서 $-3.3099\%$ 로만 움직인다. $-U^2/(2c_s^2)$ 에 $U^2 = 0.0225$ 를 넣으면 $-0.03375$ , 곧 $-3.375\%$ 다. 세 줄이 이 숫자 옆에 붙어 있다. 점성 공식에서 $\tau - 1/2$ 가 어디서 나오는지는 이산 LBE를 2차까지 펼친 자리에 적어 두었다. 여기서는 그 인자가 분자와 분모에서 통째로 지워진다는 것만 남는다.

결손의 크기는 표의 감사 열에 이미 적혀 있다. $u = 0.20$ 에서 격자 합은 $0.20000000$ , Maxwell 쪽은 $0.20800000$ , 간격은 $0.00800000$ 이다. $0.2^3 = 0.008$ 이라 간격 대 $u^3$ 은 1.00000이다. $u = 0.05$ 로 내리면 간격은 $0.00012500$ 이 되지만 비는 그대로 1이다. $u = 0.30$ 이면 간격 $0.02700000$ 이 $0.3^3$ 과 같다. 속도를 낮춰 오차를 줄이는 일과, 평형이 $\rho u^3$ 을 못 담는 일은 같은 항이다. 난수를 넣어도 두 합의 차이는 $0.000\mathrm{e}{+}00$ 으로 남는다.

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%

frame velocity U를 올리면 오른쪽 패널의 소용돌이가 밀려 가며 감쇠가 느려진다.

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

같은 유체인데 오른쪽이 덜 사라진다. $\tau$ 슬라이더를 흔들면 아래 error 숫자가 거의 안 움직인다. 표의 $-3.41\%$ 와 $-3.31\%$ 가 그 자리다.

오차는 $\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\%$ 로 거의 움직이지 않는다. 점성을 키워 오차를 묻어버릴 수 없다는 뜻이다.

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

속도를 낮춘다. 오차가 $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 버튼이 그 계산이다. 대가는 스텐실이 넓어지는 것, 메모리가 늘어나는 것, 경계 처리가 복잡해지는 것이다. 일부 속도의 가중치가 음수가 될 수 있어 안정성도 다시 봐야 한다.

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

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

레벨마다 $\tau$ 를 다시 잡아도 이 결손은 사라지지 않는다. 상대 오차가 $\tau$ 와 무관하기 때문이다. 물성으로 넣지 않은 값이 격자 구조에서 자동으로 결정되는 자리와 같다.

진단은 싸다. 같은 문제를 평균류를 얹어서 한 번 더 돌리고 감쇠율이나 항력계수를 비교한다. 차이가 $U^2$ 에 비례해서 커지면 이 결손항이다. 경계 조건을 먼저 의심할 필요가 없다.

$u/c_s$ 를 낮출지, 속도를 $\pm 2$ 까지 늘릴지#

D2Q9는 3차에서 $\rho u^3$ 을 놓친다. $U_0 = 0.2$ 에서 청구서는 $-5.97\%$ 다. 격자를 촘촘히 해도, $\tau$ 를 0.6에서 1.0으로 올려도 이 비율은 거의 안 움직인다.

줄이는 길은 세 갈래다. $u/c_s$ 를 낮추면 오차는 $U^2$ 만큼 준다. 결손항 $\partial_\gamma(\rho u_\alpha u_\beta u_\gamma)$ 를 충돌에 되돌려 넣어도 된다. 아니면 속도 집합을 넓힌다. D2Q9의 세 속도가 왜 $c^3 = c$ 로 잠기는지는 에르미트 근이 격자 속도를 고르는 자리에서 이어진다. 그 선택을 바꾸지 않는 한, 프레임을 0.2로 밀면 점성은 다시 5.97% 줄어든다.