τ=0.51에서 점성을 51배 틀리게 읽었다 — 이산 LBE의 2차 전개가 남기는 반 스텝

연속식이 주는 점성과 격자가 내놓는 점성이 51배 달랐다#

연속 BGK 방정식만 보면 완화시간 $\tau$와 동점성계수 $\nu$의 관계는 $\nu = c_s^2 \tau$다. $\tau$를 0.51로 주면 $\nu$는 0.17이어야 한다. 그런데 같은 $\tau$를 D2Q9 코드에 넣고 Taylor–Green 와동의 감쇠율을 재면 0.00333이 나온다. 51배 차이다.

오차가 아니다. 부호도 아니다. 완화시간에서 정확히 $\Delta t/2$가 빠진 값이다. 반 스텝은 두 갈래로 나온다. 특성곡선 위의 2차 테일러 전개, 그리고 충돌적분의 중점 규칙. 두 길은 같은 자리에서 만난다. $\tau$ 슬라이더를 0.51 쪽으로 밀면 주황 점선 $c_s^2\tau$가 격자 진폭에서 떨어진다.

파란 곡선이 32×32 격자에서 살아 있는 와동 진폭이다. 초록 점선은 $c_s^2(\tau - 1/2)$다. $\tau = 0.51$에서 초록은 0.00333, 주황은 0.17로 갈라진다.

격자 위의 한 줄은 미분방정식이 아니다#

우리가 실제로 돌리는 것은 이 한 줄이다.

$f_i$는 속도 $\mathbf{c}_i$ 방향의 분포함수, $\Omega_i = f_i - f_i^{\rm eq}$는 비평형 부분(평형에서 벗어난 양), $\Delta t$는 한 스텝이다.

왼쪽은 미분이 아니다. 유한한 간격 $\Delta t$를 실제로 건너뛴 차이다. 이 차이를 미분으로 바꾸려면 전개를 해야 하고, 전개에는 2차항이 있다. CFD에서 흔히 하는 실수는 1차항만 남기고 나머지를 "절단오차"로 부르는 것이다. 여기서는 그 2차항이 오차가 아니라 물리계수를 바꾼다.

왜 무시할 수 없는지는 차수를 세보면 나온다. 2차항의 크기는 $\Delta t^2 D^2 f / 2$다. 오른쪽 충돌항은 $\Delta t / \tau$로 나뉜다. 격자단위에서 $\tau$는 0.5에서 2 사이의 $O(1)$ 수이고 $\Delta t = 1$이다. 즉 2차항과 1차항의 비는 작아지지 않는다. 무한소가 되어야 할 $\Delta t$가 1로 고정된 것이 격자기체 방법의 출발점이고, 그 대가가 이 항이다.

특성곡선을 따라 2차까지 펼치면 무엇이 남는가#

$D_i \equiv \partial_t + \mathbf{c}_i \cdot \nabla$로 두면 왼쪽 항은 $D_i$의 지수함수다. 2차까지 쓰고 $\Delta t$로 나눈다.

$D_i^2$ 항이 문제의 그 항이다. 이걸 버리면 $\nu = c_s^2\tau$가 나오고, 남기면 다른 값이 나온다.

Chapman–Enskog 전개를 걸어보자. $f_i = f_i^{(0)} + \varepsilon f_i^{(1)} + \varepsilon^2 f_i^{(2)}$, $\partial_t = \varepsilon \partial_{t_1} + \varepsilon^2 \partial_{t_2}$, $\nabla = \varepsilon \nabla_1$이고 $f_i^{(0)} = f_i^{\rm eq}$다. $\varepsilon$ 차수에서

가 나온다. 이 식을 $\varepsilon^2$ 차수에 그대로 재사용하는 것이 핵심이다. $D_{i1}^2 f_i^{(0)} = -\tfrac{1}{\tau} D_{i1} f_i^{(1)}$이므로

가 된다. $D_{i1} f_i^{(1)}$ 앞에 괄호 하나가 붙었다. 점성응력은 $f_i^{(1)}$의 2차 모멘트에 $\tau$를 곱해 나오므로, 그 $\tau$는 $\tau\left(1 - \frac{\Delta t}{2\tau}\right) = \tau - \frac{\Delta t}{2}$로 바뀐다.

격자단위 $\Delta t = 1$에서는 $c_s^2(\tau - 1/2)$다. $c_s^2 = 1/3$이므로 $\tau = 0.51$이면 $\nu = 0.00333$이다. 위 시뮬레이션의 초록 점선이 바로 이 값이다.

남은 항을 분포함수 안으로 밀어 넣는다#

같은 결과를 전개 없이 얻는 길도 있다. 연속 BGK를 특성곡선 위에서 적분하되, 사다리꼴 규칙을 쓴다. 완화시간을 $\lambda$라 쓰자.

오른쪽에 미래 시각이 들어 있으니 음해식이다. 이대로는 못 돌린다. 그래서 변수를 바꾼다.

$\bar f_i$로 다시 쓰면 미래 시각 항이 상쇄되어 사라진다. $\bar\Omega_i = \bar f_i - f_i^{\rm eq} = \left(1 + \frac{\Delta t}{2\lambda}\right)\Omega_i$를 대입하면 남는 것은 이것이다.

완전히 양해식이고, 우리가 쓰는 그 한 줄과 글자까지 같다. 단 $\tau = \lambda + \Delta t/2$일 때만 그렇다. 즉 코드 안의 $\tau$는 물리 완화시간이 아니고, 물리 완화시간은 $\lambda = \tau - \Delta t/2$다. 반올림 문제가 아니다. 변수변환이 청구한 값이다. 사다리꼴 적분이 이 반 스텝을 남기는 자리는 완화시간에 붙은 $\Delta t/2$에서 이어진다.

같은 한 줄을 중점 규칙으로 다시 읽기#

사다리꼴 규칙은 $\theta$ 가중 규칙의 $\theta = 1/2$ 경우다. $\theta$를 움직여보면 왜 하필 절반인지가 보인다. 한 스텝 증폭률은 $A_\theta(z) = \frac{1 - (1-\theta)z}{1 + \theta z}$이고 정답은 $e^{-z}$다. 여기서 $z = \Delta t/\lambda$다.

$\theta$를 1/2에서 빼면 한 스텝 증폭률의 오차가 $(\theta - \tfrac12)z^2$로 기울기 2에 붙는다. $\theta = 1/2$에서만 기울기 3으로 떨어진다. 절반은 조율값이 아니라 $e^{-z}$의 2차 계수 $1/2$와 맞물리는 유일한 값이다.

Taylor–Green 와동을 τ 여섯 개로 감쇠시켰다#

주장은 여기까지다. 이제 격자에 직접 물어본다. 2D Taylor–Green 와동은 $\nu$를 정확히 알려주는 자 역할을 한다. 진폭이 $e^{-2\nu k^2 t}$로 줄기 때문이다.

import numpy as np
 
C = np.array([[0,0],[1,0],[0,1],[-1,0],[0,-1],[1,1],[-1,1],[-1,-1],[1,-1]])
W = np.array([4/9] + [1/9]*4 + [1/36]*4)
CS2 = 1/3
 
def equilibrium_d2q9(rho, ux, uy):
    cu = C[:,0,None,None]*ux + C[:,1,None,None]*uy
    u2 = ux*ux + uy*uy
    return W[:,None,None]*rho*(1 + cu/CS2 + cu*cu/(2*CS2**2) - u2/(2*CS2))
 
def taylor_green_init(n, u0):
    k = 2*np.pi/n
    x, y = np.meshgrid(np.arange(n), np.arange(n), indexing='ij')
    ux = -u0*np.cos(k*x)*np.sin(k*y)
    uy =  u0*np.sin(k*x)*np.cos(k*y)
    return np.ones((n,n)), ux, uy, k
 
def advance_lattice(f, tau):
    rho = f.sum(axis=0)
    ux = (C[:,0,None,None]*f).sum(axis=0)/rho
    uy = (C[:,1,None,None]*f).sum(axis=0)/rho
    f = f - (f - equilibrium_d2q9(rho, ux, uy))/tau      # 충돌
    for i in range(9):                                    # 이류
        f[i] = np.roll(f[i], (C[i,0], C[i,1]), axis=(0,1))
    return f, ux, uy
 
def measure_decay_rate(n, tau, u0=0.01, nstep=400, skip=100):
    rho, ux, uy, k = taylor_green_init(n, u0)
    f = equilibrium_d2q9(rho, ux, uy)
    amp = []
    for s in range(nstep):
        f, ux, uy = advance_lattice(f, tau)
        amp.append(np.sqrt((ux**2 + uy**2).mean()))
    t = np.arange(skip, nstep, dtype=float)
    a = np.log(np.array(amp[skip:]))                      # 초기 천이 구간은 버린다
    slope = np.polyfit(t, a, 1)[0]
    return -slope/(2*k*k), k
 
def viscosity_scan(n=32):
    print(f"{'tau':>6} {'nu_measured':>12} {'cs2*(tau-1/2)':>14} {'cs2*tau':>10} "
          f"{'err_half[%]':>11} {'err_naive[%]':>12}")
    for tau in (0.51, 0.55, 0.6, 0.8, 1.0, 1.5):
        nu, k = measure_decay_rate(n, tau)
        half = CS2*(tau - 0.5)
        naive = CS2*tau
        print(f"{tau:6.2f} {nu:12.6f} {half:14.6f} {naive:10.6f} "
              f"{100*abs(nu-half)/half:11.2f} {100*abs(nu-naive)/naive:12.2f}")
 
viscosity_scan()

   tau  nu_measured  cs2*(tau-1/2)    cs2*tau err_half[%] err_naive[%]
  0.51     0.003333       0.003333   0.170000        0.01        98.04
  0.55     0.016667       0.016667   0.183333        0.00        90.91
  0.60     0.033337       0.033333   0.200000        0.01        83.33
  0.80     0.100115       0.100000   0.266667        0.12        62.46
  1.00     0.167205       0.166667   0.333333        0.32        49.84
  1.50     0.337688       0.333333   0.500000        1.31        32.46

theta를 1/2에서 빼면 오른쪽 기울기가 2에 붙고, 1/2에서만 3으로 떨어진다.

왼쪽 파란 영역이 $e^{-s/\lambda}$의 면적이고, 분홍 점선이 $\theta$ 규칙이 대신 쓰는 높이다. $\theta = 1/2$만 $e^{-z}$의 2차 계수 $1/2$와 맞물린다. 표의 반 스텝 오차 0.01%는 그 중점이 맞을 때의 잔차다.

$\tau - 1/2$ 쪽 오차는 0.011.31%다. 격자 이산화가 남기는 고차 오차 수준이다. $\tau$ 쪽은 3298%다. 오차라고 부를 수 없는 크기다.

초기 진폭을 0.01로 낮게 잡은 것은 압축성 효과를 지우기 위해서다. 마하수가 0.02 수준이면 $f^{\rm eq}$의 잘린 3차항이 감쇠율에 섞이지 않는다. 앞쪽 100스텝을 버리는 것도 같은 이유다. 초기조건을 평형분포로 주면 $f^{(1)}$이 아직 0이므로, 비평형 성분이 자리를 잡기까지 몇십 스텝이 필요하다. 그 구간을 넣고 기울기를 재면 점성이 실제보다 작게 나온다.

$\tau$가 커질수록 $\tau - 1/2$ 쪽 오차가 조금씩 자라는 것도 읽을 거리다. $\tau = 1.5$에서 1.31%가 된다. 전개에서 버린 $O(\Delta t^2)$ 항이 $\tau$와 함께 커지기 때문이다. 반 스텝 보정은 완결이 아니라 2차까지의 이야기다.

표를 $\tau/(\tau - 1/2)$로 다시 읽으면 51배가 한 줄이다. $\tau = 0.51$이면 $0.51/0.01 = 51$이다. $0.170000 / 0.003333$도 51.0이다. $\tau = 0.55$는 $0.55/0.05 = 11$이고, $0.183333/0.016667$도 11이다. $\tau = 0.60$은 $0.60/0.10 = 6$, $\tau = 1.0$은 $1.0/0.5 = 2$, $\tau = 1.5$는 $1.5/1.0 = 1.5$다. 나이브 오차 98.04%, 90.91%, 83.33%, 62.46%, 49.84%, 32.46%는 이 비가 1로 다가가며 줄어든 자국이다.

반 스텝 열은 훨씬 조용하다. $\tau = 0.51$에서 측정과 예측이 둘 다 0.003333이라 오차는 0.01%다. $\tau = 0.80$은 측정 0.100115, 예측 0.100000, 차 0.000115라 0.12%다. $\tau = 1.50$은 0.337688 대 0.333333, 차 0.004355라 1.31%다. 버린 $O(\Delta t^2)$가 $\tau$와 같이 큰다는 말이 이 세 칸에 들어 있다. 초기 진폭 0.01, 앞쪽 100스텝, 격자 32는 그 잔차를 압축성에서 떼려고 고른 값이다. 마하 0.02면 잘린 3차항이 감쇠율에 섞이지 않는다. $\tau = 0.55$의 나이브 오차 90.91%는 $1 - 0.016667/0.183333$과 같다. 그 비는 1/11이라, 연속식 $\nu = c_s^2\tau$로 읽으면 열한 배 큰 점성을 코드에 적는 셈이다. $\tau = 1.5$에서도 비는 1.5라 나이브 오차 32.46%를 작은 오차라고 부르기 어렵다.

경계와 힘항이 같은 반 스텝을 다시 청구한다#

이 반 스텝은 한 번만 나타나지 않는다.

bounce-back 경계가 벽을 노드 위가 아니라 두 노드의 중간에 놓는 이유가 같다. 분포함수가 반사되는 시점이 스텝의 중간이기 때문이다. $\tau$가 0.5에서 멀어지면 이 위치가 $\tau$에 의존해 흔들리고, 그것이 유명한 2차 정확도 붕괴다.

Guo 형태의 힘항에 $\left(1 - \frac{1}{2\tau}\right)$가 붙는 이유도 같다. 힘도 특성곡선 위에서 적분되고, 같은 사다리꼴 규칙을 받는다. $\bar f$ 변환이 $\Omega_i$에만 작용하고 힘항에는 작용하지 않으므로 계수가 남는다. 이 괄호를 빼먹으면 정지해 있어야 할 유체가 천천히 흐르거나, 중력 하의 정수압 분포가 미세하게 기운다.

격자를 반으로 줄일 때 $\tau$를 그냥 두 배 하면 안 되는 이유도 여기 있다. $\nu$가 보존되어야 하니 $\tau_{\rm fine} - 1/2 = 2(\tau_{\rm coarse} - 1/2)$를 풀어야 한다. 세밀화 루프 안에서 이 한 줄을 어디에 두느냐가 오차를 가른다.

$\tau$를 노드마다 다르게 주는 경우에는 하한이 곧 이 반 스텝이다. $\tau \to 0.5$에서 $\nu \to 0$이므로, 점성이 국소적으로 0에 가까워지는 문제는 격자 자체의 벽에 부딪힌다.

0.5는 하한이 아니라 좌표원점이다#

$\tau > 0.5$를 안정조건으로 외우면 원인을 놓친다. 물리 완화시간은 $\tau - 1/2$다. $\tau = 0.51$에서 $c_s^2\tau = 0.170000$이고 $c_s^2(\tau - 1/2) = 0.003333$이라 비는 51이다. 음의 완화시간은 물리가 아니다.

코드에 적을 것은 $\tau$가 아니라 $\nu$다. $\tau = 1/2 + \nu/c_s^2$를 입력단 한 곳에 둔다. 힘항에 붙는 $\left(1 - \frac{1}{2\tau}\right)$도 같은 반 스텝의 청구서다. 그 괄호를 힘의 3차 모멘트와 같이 두는 한계는 힘 텀이 지켜야 하는 줄에서 이어진다. 감쇠가 51배로 느려지면, 어딘가에서 $\tau$를 $\nu$로 읽고 있는 것이다.