τ를 두 배로 바꿔도 Pr은 1.0000이었다 — 열 LBM에 잠겨 있는 두 상수

τ가 하나면 Pr = 1과 γ = 1 + 2/D는 유체의 물성이 아니라 격자가 정한 값이다. 둘 다 풀려면 에너지에 분포함수를 따로 세워야 한다.

물성으로 넣지 않은 값이 1.0000으로 나왔다#

D2Q9 격자에 사인 모양 전단파와 사인 모양 온도파를 하나씩 얹고, 각각의 진폭이 줄어드는 속도로 동점성계수 $\nu$ 와 열확산계수 $\alpha$ 를 쟀다. 이완 시간 $\tau$ 를 0.6, 0.8, 1.2로 두 배씩 바꿔가며 같은 측정을 반복했다. 세 번 모두 프란틀 수(Prandtl number, 운동량 확산과 열 확산의 비)가 1.0001, 1.0000, 1.0001로 나왔다.

물성으로 넣은 적이 없는 값이다. 격자가 정한 값이다.

이 글은 그 잠금장치가 코드의 어느 줄에 있는지 짚고, 에너지에 분포함수를 하나 더 세우면 무엇이 풀리는지를 같은 격자에서 다시 재본 기록이다. 잠긴 상수는 $\mathrm{Pr}$ 하나가 아니다. 비열비 $\gamma$ 도 같이 잠겨 있고, D2Q9에서 그 값은 2다. 공기의 1.4가 아니다. 그리고 풀어준 뒤에도 $\tau_g$ 를 무한정 밀 수는 없다 — 어디서 무너지는지도 숫자로 잰다.

아래에서 자물쇠를 직접 걸었다 풀었다 해보자.

tau_f (momentum)0.80tau_g (energy)0.80speed1.0x

nu 0.1000alpha 0.1000Pr 1.000step 0A_u 1.000A_T 1.000

Locked, the orange curve hides underneath the blue one at every tau_f you try — one relaxation time cannot hold two transport coefficients apart. Unlock it and drag tau_g: the thermal wave now outlives the shear wave (Pr > 1) or dies first (Pr < 1).

잠긴 상태에서는 $\tau_f$ 를 어떻게 움직여도 오른쪽 그래프의 주황색(온도) 곡선이 파란색(속도) 곡선 아래에 완전히 숨는다. 자물쇠를 풀고 $\tau_g$ 를 끌면 두 곡선이 갈라진다. 갈라지는 정도가 곧 $\mathrm{Pr}$ 이다.

τ 하나가 두 군데에서 쓰인다#

BGK 충돌 항을 쓴 격자 볼츠만 방정식은 이렇다.

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

$f_i$ 는 $i$ 번째 격자 방향의 분포함수, $\boldsymbol{e}_i$ 는 그 방향의 이산 속도, $\tau$ 는 이완 시간이다. Chapman–Enskog 전개(평형에서 멀어진 정도를 작은 매개변수로 두고 차수별로 나누는 전개)를 1차까지 하면 점성응력이 $f$ 의 1차 비평형 항에서 나온다. 그 결과가 잘 알려진 이 식이다.

\nu = c_s^{2}\left(\tau - \frac{\Delta t}{2}\right)

$c_s$ 는 격자 음속, $\Delta t/2$ 는 이산화가 남긴 보정이다. 이 반 칸이 어디서 오는지는 LBM 이산화가 남긴 Δt/2 세 자리에 따로 정리해뒀다.

문제는 그다음이다. 내부에너지를 같은 $f$ 의 2차 모멘트로 정의하면, 온도 방정식도 같은 전개에서 같은 1차 비평형 항을 통해 나온다. 열확산계수는 형태까지 똑같다.

\alpha = c_s^{2}\left(\tau - \frac{\Delta t}{2}\right)

두 식의 우변이 글자 하나까지 같다. 그러면 결론은 하나다.

\mathrm{Pr} = \frac{\nu}{\alpha} = 1

이유는 한 줄로 요약된다. 운동량 플럭스와 열 플럭스는 둘 다 같은 분포함수의 같은 비평형 항에서 나오고, 그 항 앞에 붙는 시간 상수는 $\tau$ 하나뿐이다. 상수가 하나면 비율은 선택할 수 없다.

기체에서 $\mathrm{Pr}$ 이 1 근처라는 것은 실험적 사실이고, 마찰과 열전달을 잇는 레이놀즈 상사도 거기서 나온다. 하지만 그것은 근사이고 선택이다. 여기서는 강제다. 물을 넣어도, 액체 금속을 넣어도 1이 나온다.

에너지에 분포함수를 따로 세운다#

풀어주는 방법은 구조적으로 단순하다. 상수가 하나뿐이어서 생긴 문제니 상수를 하나 더 만든다. $f$ 는 질량과 운동량만 담당하고, 에너지는 두 번째 분포함수 $g_i$ 에 맡긴다. 이중 분포함수 (double-distribution-function, DDF) 방법이다.

$g$ 가 만족해야 할 모멘트 조건은 두 줄이다.

\sum_i g_i = \rho E, \qquad \sum_i \boldsymbol{e}_i\, g_i = \rho E \boldsymbol{u} + p\,\boldsymbol{u}

$E = c_v T + |\boldsymbol{u}|^{2}/2$ 는 단위질량당 총에너지, $p\boldsymbol{u}$ 는 압력이 하는 일이다. 두 번째 줄에 $p\boldsymbol{u}$ 가 들어가는 것이 핵심이다. 에너지의 대류 플럭스는 $\rho E \boldsymbol{u}$ 하나가 아니라 압력 일까지 합친 $\rho H \boldsymbol{u}$ 이고, 그래서 $g$ 의 평형 분포는 $f$ 의 평형 분포를 그대로 베낄 수 없다.

$g$ 를 $\tau_g$ 로 이완시키고 같은 Chapman–Enskog 전개를 돌리면 열확산계수는 이제 $\tau_g$ 를 본다.

\mathrm{Pr} = \frac{\nu}{\alpha} = \frac{\tau_f - \Delta t/2}{\tau_g - \Delta t/2}

$\tau_f$ 는 레이놀즈 수가 정하고, $\tau_g$ 는 프란틀 수가 정한다. 두 요구가 서로 다른 손잡이를 잡게 된 것이다.

Python으로 잠금과 해제를 같은 격자에서 쟀다#

말로 끝내지 않고 재봤다. $f$ 는 D2Q9, $g$ 는 D2Q5(온도 전용 5속도 격자)로 세우고, $y$ 방향으로만 변하는 사인파를 초기조건으로 준다. 전단파의 진폭은 $\exp(-\nu k^{2} t)$ 로, 온도파의 진폭은 $\exp(-\alpha k^{2} t)$ 로 줄어든다. 진폭 로그를 최소제곱으로 맞추면 $\nu$ 와 $\alpha$ 가 나온다.

import math
 
NY, CS2 = 64, 1.0 / 3.0
K = 2.0 * math.pi / NY
 
EX9 = [0, 1, 0, -1, 0, 1, -1, -1, 1]
EY9 = [0, 0, 1, 0, -1, 1, 1, -1, -1]
W9 = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]
 
EY5 = [0, 0, 1, 0, -1]
W5 = [1/3, 1/6, 1/6, 1/6, 1/6]
 
 
def feq_d2q9(rho, ux, uy):
    u2 = ux * ux + uy * uy
    return [w * rho * (1 + 3 * (ex * ux + ey * uy)
                       + 4.5 * (ex * ux + ey * uy) ** 2 - 1.5 * u2)
            for w, ex, ey in zip(W9, EX9, EY9)]
 
 
def geq_d2q5(temp):
    return [w * temp for w in W5]
 
 
def fit_decay_rate(samples):
    """samples = [(스텝, 진폭)] -> ln(진폭)의 기울기를 최소제곱으로."""
    n = len(samples)
    xs = [s for s, _ in samples]
    ys = [math.log(a) for _, a in samples]
    mx, my = sum(xs) / n, sum(ys) / n
    num = sum((x - mx) * (y - my) for x, y in zip(xs, ys))
    den = sum((x - mx) ** 2 for x in xs)
    return -num / den
 
 
def run_shear_wave(tau_f, steps=3000, amp=1e-3):
    f = [[0.0] * NY for _ in range(9)]
    for y in range(NY):
        for i, v in enumerate(feq_d2q9(1.0, amp * math.sin(K * y), 0.0)):
            f[i][y] = v
    log = []
    for step in range(steps + 1):
        rho = [sum(f[i][y] for i in range(9)) for y in range(NY)]
        ux = [sum(f[i][y] * EX9[i] for i in range(9)) / rho[y] for y in range(NY)]
        if step % 200 == 0:
            a = 2.0 / NY * sum(ux[y] * math.sin(K * y) for y in range(NY))
            log.append((step, a))
        post = [[0.0] * NY for _ in range(9)]
        for y in range(NY):
            eq = feq_d2q9(rho[y], ux[y], 0.0)
            for i in range(9):
                post[i][y] = f[i][y] - (f[i][y] - eq[i]) / tau_f
        for i in range(9):
            for y in range(NY):
                f[i][(y + EY9[i]) % NY] = post[i][y]
    return fit_decay_rate(log) / (K * K)
 
 
def run_thermal_wave(tau_g, amp=1e-3):
    """진폭이 항상 e^-2.5만큼 줄어들도록 스텝 수를 tau_g에 맞춰 잡는다."""
    alpha_th = CS2 * (tau_g - 0.5)
    steps = max(240, min(12000, int(2.5 / (alpha_th * K * K))))
    every = max(1, steps // 12)
    g = [[0.0] * NY for _ in range(5)]
    for y in range(NY):
        for i, v in enumerate(geq_d2q5(amp * math.sin(K * y))):
            g[i][y] = v
    log = []
    for step in range(steps + 1):
        temp = [sum(g[i][y] for i in range(5)) for y in range(NY)]
        if step % every == 0:
            a = 2.0 / NY * sum(temp[y] * math.sin(K * y) for y in range(NY))
            log.append((step, a))
        post = [[0.0] * NY for _ in range(5)]
        for y in range(NY):
            eq = geq_d2q5(temp[y])
            for i in range(5):
                post[i][y] = g[i][y] - (g[i][y] - eq[i]) / tau_g
        for i in range(5):
            for y in range(NY):
                g[i][(y + EY5[i]) % NY] = post[i][y]
    return fit_decay_rate(log) / (K * K)
 
 
def tau_for(nu, target_pr):
    return 0.5 + nu / (target_pr * CS2)
 
 
TAU_F = 0.8
nu = run_shear_wave(TAU_F)
print(f"tau_f = {TAU_F}   nu(theory) = {CS2*(TAU_F-0.5):.6f}   nu(measured) = {nu:.6f}")
print()
print("[A] single distribution: one tau relaxes momentum AND energy")
print(f"{'tau':>6} {'nu':>10} {'alpha':>10} {'Pr':>8}")
for t in (0.6, 0.8, 1.2):
    n_, a_ = run_shear_wave(t), run_thermal_wave(t)
    print(f"{t:6.2f} {n_:10.6f} {a_:10.6f} {n_/a_:8.4f}")
print()
print("[B] double distribution: tau_g chosen for a target Pr (tau_f = 0.8)")
print(f"{'gas':>8} {'Pr(target)':>11} {'tau_g':>8} {'alpha':>10} {'Pr(meas)':>9} {'err%':>7}")
for name, pr in (("mercury", 0.025), ("air", 0.71), ("Pr=1", 1.0), ("water", 7.0)):
    tg = tau_for(nu, pr)
    a_ = run_thermal_wave(tg)
    prm = nu / a_
    print(f"{name:>8} {pr:11.3f} {tg:8.4f} {a_:10.6f} {prm:9.4f} {100*(prm-pr)/pr:7.2f}")
 
print()
print("[C] how far can tau_g be pushed?  (theory: alpha = cs2*(tau_g-0.5))")
print(f"{'tau_g':>7} {'alpha(th)':>10} {'alpha(meas)':>12} {'err%':>7}")
for tg in (0.51, 0.55, 0.7, 1.0, 2.0, 4.0, 8.0, 12.5):
    th = CS2 * (tg - 0.5)
    ms = run_thermal_wave(tg)
    print(f"{tg:7.2f} {th:10.5f} {ms:12.5f} {100*(ms-th)/th:7.2f}")

tau_f = 0.8   nu(theory) = 0.100000   nu(measured) = 0.100057
 
[A] single distribution: one tau relaxes momentum AND energy
   tau         nu      alpha       Pr
  0.60   0.033368   0.033363   1.0001
  0.80   0.100057   0.100060   1.0000
  1.20   0.233144   0.233124   1.0001
 
[B] double distribution: tau_g chosen for a target Pr (tau_f = 0.8)
     gas  Pr(target)    tau_g      alpha  Pr(meas)    err%
 mercury       0.025  12.5069   2.047239    0.0489   95.50
     air       0.710   0.9228   0.140963    0.7098   -0.03
    Pr=1       1.000   0.8002   0.100117    0.9994   -0.06
   water       7.000   0.5429   0.014308    6.9931   -0.10
 
[C] how far can tau_g be pushed?  (theory: alpha = cs2*(tau_g-0.5))
  tau_g  alpha(th)  alpha(meas)    err%
   0.51    0.00333      0.00334    0.16
   0.55    0.01667      0.01668    0.10
   0.70    0.06667      0.06672    0.08
   1.00    0.16667      0.16667   -0.00
   2.00    0.50000      0.49624   -0.75
   4.00    1.16667      1.11274   -4.62
   8.00    2.50000      1.94812  -22.08
  12.50    4.00000      2.04752  -48.81

[A] 표가 잠금장치다. $\tau$ 를 0.6에서 1.2로 두 배 키우면 $\nu$ 는 일곱 배가 되는데 $\mathrm{Pr}$ 은 소수점 넷째 자리까지 1이다. [B] 표는 해제다. 공기는 0.7098, 물은 6.9931로 목표값에서 0.1% 안쪽에 들어왔다. 단, 맨 윗줄 수은은 95% 틀렸다. 그 줄은 뒤에서 따로 다룬다.

잠긴 상수는 하나 더 있다 — 비열비#

$\mathrm{Pr}$ 만 보고 넘어가면 두 번째 자물쇠를 놓친다. 격자 위의 입자는 $D$ 개의 병진 방향으로만 움직인다. 회전도, 진동도 없다. 그러면 등적비열은 병진 자유도만 세어 $c_v = DR/2$ 가 되고, 기체상수는 $R = c_v / (D/2)$ 로 되받아진다. 비열비는 자동으로 결정된다.

c_v = \frac{(D + n_0)\,R}{2}, \qquad \gamma = 1 + \frac{2}{D + n_0}

$n_0$ 은 병진 외에 추가로 실어주는 내부 자유도의 개수다. 아무것도 하지 않으면 $n_0 = 0$ 이다.

격자	$n_0$	$c_v$	$\gamma$	대응하는 기체
D2Q9	0	$1.0R$	2.000	없음
D3Q19	0	$1.5R$	1.667	단원자(Ar, He)
D2Q9	3	$2.5R$	1.400	공기
D3Q19	2	$2.5R$	1.400	공기
D2Q9	4	$3.0R$	1.333	수증기

2차원에서 손대지 않은 격자의 $\gamma$ 는 2다. 음속은 $\sqrt{\gamma R T}$ 이므로 공기 대비 $\sqrt{2/1.4} = 1.195$ 배, 즉 19.5% 빠르다. 마하수, 충격파 각도, 노즐의 초킹 조건이 전부 그만큼 어긋난다. 압축성 계산에서는 이쪽이 $\mathrm{Pr}$ 보다 먼저 문제가 된다.

n0 extra DOF0.0speed1.0x

gamma 2.000c_v 1.00 Rsound speed 0.0%

Start at n0 = 0 with D2Q9: gamma = 2, and the blue pulse runs 19% fast against air. Drag n0 until the blue marker lands on the dashed line — the number you land on is how many internal degrees of freedom the energy distribution has to carry.

D2Q9와 $n_0 = 0$ 에서 출발해 두 음향 펄스가 얼마나 벌어지는지 본 다음, $n_0$ 을 끌어 파란 표식을 점선 위에 올려보자. 올라간 자리의 숫자가 에너지 분포함수에 실어야 하는 내부 자유도 개수다. 3차원에서 공기를 맞추려면 2개, 2차원에서는 3개가 필요하다.

τ_g를 어디까지 밀 수 있는가#

[C] 표는 $\alpha = c_s^{2}(\tau_g - \Delta t/2)$ 가 언제까지 믿을 만한지를 잰 것이다. $\tau_g$ 가 2까지는 오차가 0.8% 안쪽이다. 4에서 4.6%, 8에서 22%, 12.5에서는 49%가 벌어진다. 이론값의 절반밖에 안 나온다.

이유는 그 식의 출처에 있다. $\alpha = c_s^{2}(\tau_g - \Delta t/2)$ 는 Chapman–Enskog 전개의 1차 결과다. 전개가 성립하려면 이완 시간이 유동의 시간 규모보다 짧아야 한다. $\tau_g$ 가 커지면 버린 2차 항 — 파수의 네제곱에 비례하는 초확산 항 — 이 1차 항과 같은 크기로 자란다. 실질 상한은 격자 단위로 $\alpha \approx 0.5$ 근처다.

그래서 [B]의 수은 줄이 깨졌다. $\mathrm{Pr} = 0.025$ 를 $\tau_f = 0.8$ 에서 얻으려면 $\alpha = 4.0$ 이 필요한데, 그 값은 위 상한의 여덟 배다. 처방은 $\tau_g$ 를 더 올리는 것이 아니라 $\nu$ 를 내리는 것이다. $\alpha \le 0.5$ 를 지키면서 $\mathrm{Pr} = 0.025$ 를 얻으려면 $\nu \le 0.0125$ , 즉 $\tau_f \le 0.5375$ 여야 한다. 이번에는 반대쪽 벽이다. $\tau_f$ 가 0.5에 붙으면 BGK가 불안정해진다.

정리하면 이렇다. 높은 $\mathrm{Pr}$ 은 $\tau_g \to 0.5$ 쪽에서 안정성 벽을 만나고, 낮은 $\mathrm{Pr}$ 은 $\tau_g$ 가 커지는 쪽에서 정확도 벽을 만난다. DDF가 준 것은 무제한의 자유가 아니라 창이다.

점성 가열은 어느 장부에 오르는가#

$f$ 와 $g$ 를 갈라놓으면 새 문제가 하나 생긴다. 에너지 방정식에는 점성 소산 항 $\boldsymbol{u} \cdot (\nabla \cdot \boldsymbol{\Pi})$ 가 들어 있다. 그런데 $\boldsymbol{\Pi}$ 는 $f$ 의 1차 비평형에서 나오는 양이고, $\tau_f$ 로 이완한다.

\boldsymbol{\Pi}^{(1)} = -2\rho\, c_s^{2}\left(\tau_f - \frac{\Delta t}{2}\right) \boldsymbol{S}

$\boldsymbol{S}$ 는 변형률 텐서다. $g$ 는 $\tau_g$ 로 이완하므로, $g$ 의 전개가 되돌려주는 소산 항 앞에는 $\tau_g$ 가 붙는다. 두 이완 시간이 다른 순간 계수가 어긋난다. $\mathrm{Pr} = 1$ 일 때만 저절로 맞는다. 결합형 DDF가 $g$ 방정식에 보정항을 하나 더 붙이는 이유가 이것이다.

보정항을 격자점마다 국소적으로 계산하려면 $\boldsymbol{\Pi}^{(1)}$ 을 모멘트로 닫아야 한다. Grad의 13-모멘트 근사가 그 자리를 메운다.

f_i \simeq w_i\left[\rho + \frac{\rho\, \boldsymbol{e}_i \cdot \boldsymbol{u}}{c_s^{2}} + \frac{(\boldsymbol{e}_i \boldsymbol{e}_i - c_s^{2}\boldsymbol{I}) : \left(\rho \boldsymbol{u}\boldsymbol{u} + \boldsymbol{\Pi}^{(1)}\right)}{2 c_s^{4}}\right]

$\boldsymbol{\Pi}^{(1)} = 0$ 을 넣으면 익숙한 평형 분포가 그대로 나온다. 그 다항식이 어디서 왔는지는 9개 화살표 안의 Maxwell–Boltzmann에 적어뒀다. Grad 근사는 그 전개를 한 차수 더 끌고 가서, 비평형 응력까지 분포함수 안에 되돌려 넣는 장치다. 2007년의 결합형 DDF 논문(Li 등의 압축성 Navier–Stokes용 모델)과 저마하 디커플링 모델이 갈라지는 지점도 여기다. 전자는 보정항을 명시적으로 붙이고, 후자는 점성 가열을 아예 버린다.

분포함수 하나가 감당할 수 있는 물리의 크기#

$f$ 하나로 실을 수 있는 것은 $\rho$ , $\boldsymbol{u}$ , 그리고 이완 시간 하나까지다. 온도를 같은 $f$ 의 2차 모멘트에 얹는 순간 $\mathrm{Pr}$ 과 $\gamma$ 는 유체의 물성이 아니라 격자 상수가 된다. 저마하 부시네스크 계산에서는 온도가 어차피 수동 스칼라라 이 문제가 드러나지 않는다. 압축성 열유동으로 넘어가면 두 상수가 동시에 청구된다.

넘겨받은 열 LBM 코드에서 세 줄만 먼저 찾으면 된다. 온도를 $f$ 의 모멘트에서 꺼내는가 아니면 별도 배열에서 꺼내는가. $\gamma$ 가 상수로 하드코딩돼 있는가 아니면 $D + n_0$ 에서 계산되는가. 그리고 $g$ 의 충돌 항 옆에 $\boldsymbol{\Pi}^{(1)}$ 을 쓰는 보정항이 있는가. 세 번째가 없는데 $\tau_f \ne \tau_g$ 라면, 그 코드의 점성 가열은 계산된 적이 없는 값이다.