τ를 화학종마다 다르게 줬더니 총 질량이 16.8% 어긋났다 — 다화학종 LBM의 확산속도 보정

화학종을 셋으로 늘린 날 밀도 장부가 맞지 않았다#

연소 계산에 화학종을 하나 붙일 때마다 분포함수도 한 벌씩 늘어난다. 화학종마다 확산계수가 다르니 완화시간 $\tau_k$도 다르게 준다. 여기까지는 교과서대로다. 그런데 몇백 스텝을 돌리면 $\sum_k Y_k$가 1에서 벗어난다. 질량분율(mass fraction, 전체 질량 중 그 화학종이 차지하는 몫)의 합은 정의상 1이어야 한다.

이 글은 그 어긋남이 어디서 오는지 Chapman–Enskog 전개로 짚고, 한 줄짜리 보정으로 되돌린다. 3성분 혼합을 2000스텝 돌려 최대 결손이 0.168에서 0.0005로 줄어드는 과정을 숫자로 싣는다. 먼저 결론을 적자면 이건 버그가 아니라 모델이다. $\tau_k$를 서로 다르게 준 순간, 연속방정식이 깨지도록 되어 있다.

분포함수를 세 벌로 나눠 쓰는 구조#

압축성 다상 다화학종 LBM은 분포함수를 용도별로 쪼갠다. $f_i$가 유동장, $g_i$가 에너지, $h_{k,i}$가 $k$번째 화학종을 맡는다. 셋 다 같은 격자속도 $\mathbf{c}_i$ 위에서 같은 형태의 이산 Boltzmann 방정식을 푼다.

$h_{k,i}$는 방향 $i$에 실린 화학종 $k$의 분포함수, $\tau_k$는 그 화학종의 완화시간, $\delta t$는 시간스텝이다.

평형분포는 화학종 밀도 $\rho_k$와 유동속도 $\mathbf{u}$만 쓴다. 유동장 $f_i$와 달리 2차항이 필요 없다. 화학종 방정식은 이류-확산까지만 복원하면 되기 때문이다.

$w_i$는 격자 가중치, $c_s$는 격자 음속이다. 화학종 밀도는 $\rho_k=\sum_i h_{k,i}$, 질량분율은 $Y_k=\rho_k/\rho$로 읽는다.

$D_A,D_B,D_C$ 슬라이더를 벌리면 세 확산계수가 갈라지고, 아래 패널의 질량 합이 1에서 밀려난다.

세 슬라이더를 같은 값에 맞추면 아래 패널이 반올림까지 납작해진다. 벌린 채로 두면 매 스텝 오차가 쌓인다. 체크박스는 보정속도다. 켜면 그 합이 다시 1 쪽으로 접힌다.

Chapman–Enskog가 $h_k$에서 꺼내 놓는 식#

위 이산식을 $\delta t$로 Taylor 전개하고, 분포함수를 $h_{k,i}=h_{k,i}^{\mathrm{eq}}+\epsilon h_{k,i}^{(1)}+\cdots$로 펼친다. 0차 모멘트를 차수별로 모으면 화학종 하나짜리 이류-확산 방정식이 나온다.

$D_k$가 그 화학종의 확산계수다. $-1/2$는 이산화가 만들어내는 항이지 물리가 아니다. 같은 구조의 $-1/2$가 유동장에서는 점성, 에너지장에서는 열확산으로 나타난다. 그 두 장부를 하나로 묶으면 Prandtl 수가 상수로 잠긴다는 이야기는 열 LBM에 잠겨 있는 두 상수에 적어 두었다.

여기까지만 보면 아름답다. $\tau_k$ 하나를 돌려 확산계수를 원하는 값에 맞출 수 있다. 실제로 그렇게 맞는지는 뒤에서 잰다.

$\tau_k$가 서로 다르면 무엇이 남는가#

문제는 화학종 방정식 하나가 아니라 전부를 더할 때 생긴다. $k$에 대해 위 식을 합하고 $\sum_k\rho_k=\rho$를 쓰면 이렇게 된다.

좌변은 연속방정식이다. 우변은 0이어야 한다. 그런데 우변은 일반적으로 0이 아니다.

모든 $D_k$가 같을 때만 $\sum_k D_k\nabla Y_k=D\nabla\sum_k Y_k=D\nabla(1)=0$이 된다. 화학종마다 $\tau_k$를 다르게 준다는 것은 이 취소를 스스로 깨는 일이다. 같은 말을 플럭스로 적으면 더 분명하다. Fick 플럭스 $\mathbf{J}_k=-\rho D_k\nabla Y_k$의 합이 0이 아니다.

물리적으로 이건 알려진 사실이다. Fick 법칙은 2성분 혼합의 근사이고, 3성분 이상에서는 Stefan–Maxwell 관계가 맞다. LBM은 이 근사를 화학종마다 독립된 $\tau_k$로 옮겨 적으면서 결손을 그대로 물려받는다.

3성분 혼합을 2000스텝 돌려 장부를 맞춰봤다#

D1Q3 격자에 화학종 세 벌을 올렸다. $\rho=1$, $\mathbf{u}=0$인 정지 혼합이다. 확산 말고는 아무 일도 일어나지 않으므로 $\sum_k\rho_k$는 영원히 1이어야 한다. 먼저 감쇠하는 사인 모드로 $D_k$가 정말 $c_s^2(\tau_k-1/2)$인지 재고, 그다음 $D_k$ 조합을 바꿔가며 결손을 쟀다.

import numpy as np
 
N, CS2 = 256, 1.0/3.0
W = np.array([2.0/3.0, 1.0/6.0, 1.0/6.0])   # D1Q3 가중치: c = 0, +1, -1
C = np.array([0, 1, -1])
 
def feq_species(rho_k, u):
    return rho_k[None, :] * W[:, None] * (1.0 + C[:, None] * u[None, :] / CS2)
 
def correction_velocity(rho_k, D):
    dY = np.array([(np.roll(r, -1) - np.roll(r, 1)) * 0.5 for r in rho_k])
    return np.sum(D[:, None] * dY, axis=0)           # V_c = sum_k D_k grad Y_k
 
def lbm_species_step(h, D, use_correction):
    rho_k = h.sum(axis=1)                            # (K, N)
    u = correction_velocity(rho_k, D) if use_correction else np.zeros(N)
    tau = D / CS2 + 0.5
    for k in range(len(D)):
        h[k] += (feq_species(rho_k[k], u) - h[k]) / tau[k]
        for i in range(3):
            h[k, i] = np.roll(h[k, i], C[i])
    return h
 
def diffusion_leak_run(D, steps, use_correction):
    x = np.arange(N)
    r1 = 0.30 + 0.20 * np.sin(2 * np.pi * x / N)
    r2 = 0.45 + 0.15 * np.cos(4 * np.pi * x / N)
    r3 = 1.0 - r1 - r2                               # 시작할 때 합은 정확히 1
    h = np.array([feq_species(r, np.zeros(N)) for r in [r1, r2, r3]])
    worst = 0.0
    for _ in range(steps):
        h = lbm_species_step(h, D, use_correction)
        worst = max(worst, np.abs(h.sum(axis=(0, 1)) - 1.0).max())
    return worst
 
def measure_D_eff(tau, steps=400):
    x = np.arange(N)
    amp0 = 0.01
    rho = 0.5 + amp0 * np.sin(2 * np.pi * x / N)     # 사인 모드 하나를 넣고 감쇠율을 잰다
    h = feq_species(rho, np.zeros(N))
    for _ in range(steps):
        r = h.sum(axis=0)
        h += (feq_species(r, np.zeros(N)) - h) / tau
        for i in range(3):
            h[i] = np.roll(h[i], C[i])
    amp = np.abs(np.fft.rfft(h.sum(axis=0))[1]) * 2.0 / N
    k2 = (2 * np.pi / N) ** 2
    return -np.log(amp / amp0) / (k2 * steps)
 
print("tau    D_theory=cs^2(tau-1/2)   D_measured   error[%]")
for tau in (0.55, 0.80, 1.20, 2.00):
    dth, dm = CS2 * (tau - 0.5), measure_D_eff(tau)
    print(f"{tau:.2f}   {dth:.6f}                {dm:.6f}     {100*(dm-dth)/dth:+.3f}")
 
print()
print("D1,D2,D3            correction   max|sum(rho_k)-1|")
for D in (np.array([0.05, 0.05, 0.05]),
          np.array([0.05, 0.10, 0.02]),
          np.array([0.01, 0.10, 0.30]),
          np.array([0.005, 0.10, 0.50])):
    for corr in (False, True):
        print(f"{D[0]:.3f},{D[1]:.3f},{D[2]:.3f}     {str(corr):5s}        "
              f"{diffusion_leak_run(D, 2000, corr):.3e}")

tau    D_theory=cs^2(tau-1/2)   D_measured   error[%]
0.55   0.016667                0.016874     +1.243
0.80   0.100000                0.100137     +0.137
1.20   0.233333                0.233122     -0.091
2.00   0.500000                0.498135     -0.373
 
D1,D2,D3            correction   max|sum(rho_k)-1|
0.050,0.050,0.050     False        1.865e-13
0.050,0.050,0.050     True         2.091e-13
0.050,0.100,0.020     False        5.050e-02
0.050,0.100,0.020     True         5.201e-05
0.010,0.100,0.300     False        1.155e-01
0.010,0.100,0.300     True         1.354e-04
0.005,0.100,0.500     False        1.683e-01
0.005,0.100,0.500     True         5.215e-04

$D_k$ 슬라이더를 벌리면 빨간 합 화살표의 길이가 달라진다.

$D_k$와 기울기를 아무렇게나 놓아도 빨간 합 화살표는 거의 0에 서지 않는다. 체크박스를 켜면 초록 구간이 붙으면서 합이 정확히 0으로 접힌다.

윗 표부터 읽자. $D_k=c_s^2(\tau_k-1/2)$는 0.1% 수준으로 맞는다. $\tau=0.55$에서 1.24%로 튀는 건 $\tau$가 $1/2$에 붙을수록 커지는 격자 오차다. 이 관계 자체는 믿어도 된다는 뜻이다.

아랫 표가 본론이다. 세 확산계수를 모두 0.05로 맞추면 결손이 $1.9\times10^{-13}$, 배정도 반올림이다. 2.5배만 벌려도 $5\times10^{-2}$가 된다. 100배로 벌리면 0.168까지 간다. $\rho=1$이니 이건 총 질량의 16.8%다.

$c_s^2=1/3$을 $D_k=c_s^2(\tau_k-1/2)$에 넣으면 윗 표의 이론값이 그대로 나온다. $\tau=0.80$이면 $(0.30)/3=0.100000$이다. 측정은 0.100137, 오차 $+0.137\%$다. $\tau=2.00$이면 $(1.50)/3=0.500000$, 측정 0.498135, 오차 $-0.373\%$다. 부호가 바뀐다. $\tau=0.55$는 $(0.05)/3=0.016667$인데 측정 0.016874, $+1.243\%$다. $\tau=1.20$은 이론 0.233333, 측정 0.233122, $-0.091\%$다. $1/2$에서 멀어질수록 퍼센트 오차는 작아진다.

결손 비도 아랫 표에서 바로 나뉜다. $D=(0.050,\,0.100,\,0.020)$에서 $5.050\times10^{-2}$를 $5.201\times10^{-5}$로 나누면 971이다. $D=(0.005,\,0.100,\,0.500)$에서 $1.683\times10^{-1}$을 $5.215\times10^{-4}$로 나누면 323이다. 확산계수 비가 2.5배일 때보다 100배일 때 보정 배율이 오히려 작다. 세 값이 모두 0.05이면 결손은 $1.865\times10^{-13}$, 보정을 켜도 $2.091\times10^{-13}$이다. 취소가 이미 성립한 자리라 보정속도가 할 일이 없다. 2000스텝, $\rho=1$에서 0.168은 총질량의 16.8%이고, 보정 뒤 $5.215\times10^{-4}$는 그 자리의 잔차다.

0.168을 0.0005로 바꾼 한 줄#

고치는 방법은 오래 알려져 있다. 모든 화학종에 같은 보정속도 $\mathbf{V}_c$를 얹어 플럭스 합을 강제로 0으로 만든다.

이 $\mathbf{V}_c$를 넣으면 $\sum_k\mathbf{J}_k=-\rho\sum_k D_k\nabla Y_k+\rho\mathbf{V}_c\sum_k Y_k=0$이 된다. $\sum_k Y_k=1$이라는 조건을 그대로 쓴 것이다.

LBM에서는 이걸 소스항으로 붙이지 않는다. 평형분포의 속도를 $\mathbf{u}$에서 $\mathbf{u}+\mathbf{V}_c$로 바꾸기만 하면 된다. 위 코드의 correction_velocity 한 줄과 use_correction 분기가 전부다. 속도를 올렸을 때 확산이 줄어드는 여분 플럭스는 대류-확산 LBM이 남기는 항과 같은 종류다.

표의 숫자로도 같다. 결손 $5.05\times10^{-2}$가 $5.20\times10^{-5}$로, $1.68\times10^{-1}$이 $5.22\times10^{-4}$로 내려간다. 각각 971배, 323배다.

보정을 넣어도 0이 되지 않는 잔차#

눈에 걸리는 게 하나 있다. 확산계수가 모두 같을 때는 $1.9\times10^{-13}$인데, 보정을 켠 경우는 $5\times10^{-4}$에서 멈춘다. 기계 정밀도가 아니다.

이유는 이산 연산자가 두 개이기 때문이다. correction_velocity는 중심차분으로 $\nabla Y_k$를 계산한다. 반면 LBM이 실제로 만들어내는 확산 플럭스는 비평형 부분 $h_{k,i}^{(1)}$의 1차 모멘트, 즉 스트리밍이 정의하는 또 다른 이산 기울기다. 두 기울기가 $O(\delta x^2)$만큼 어긋난다. 그 차이가 취소되지 않고 남는다.

없애려면 보정속도도 같은 비평형 모멘트에서 뽑으면 된다. $\mathbf{V}_c=\frac{1}{\rho}\sum_k\big(1-\frac{1}{2\tau_k}\big)\sum_i\mathbf{c}_i h_{k,i}^{(1)}$처럼 쓰는 방식이다. 중심차분이 한 스텝 뒤진 값을 쓰는 것도 같은 이유로 잔차를 만든다. 프레임을 밀면 점성 자체가 달라지는 자리는 평형분포의 3차 모멘트 결손이다.

다만 실무에서는 $5\times10^{-4}$와 $1.7\times10^{-1}$의 차이가 크다. 전자는 축적되지 않고 격자를 조이면 줄어든다. 후자는 시간에 따라 자란다.

연소 코드에서 이 항이 얼굴을 내미는 자리#

메탄 공기 연소의 Lewis 수는 화학종마다 0.3에서 1.8까지 흩어진다. H와 H$_2$는 극단적으로 빠르고, 큰 탄화수소는 느리다. 확산계수로 보면 20배 이상 벌어진다. 위 표의 세 번째 행이 딱 그 범위다.

결손은 밀도 오차로만 끝나지 않는다. 화학 반응률은 $Y_k$의 지수에 들어간다. 합이 1을 넘으면 반응률이 과대평가되고, 그 열이 다시 밀도를 밀어 발산으로 간다. 많은 코드가 매 스텝 $Y_k$를 합으로 나눠 정규화하는데, 그건 결손을 지우는 게 아니라 화학종 사이에 임의로 나눠 담는 것이다. 어느 화학종이 손해를 볼지는 격자 번호에 달린다.

가장 흔한 실무 처방은 가장 무거운 화학종 하나를 잉여종(bath species)으로 지정해 $Y_{N_2}=1-\sum_{k\neq N_2}Y_k$로 직접 계산하는 것이다. 이러면 합은 정확히 1이지만 결손이 전부 N$_2$에 쌓인다. 보정속도는 그 부담을 $Y_k$ 비율대로 나눈다.

결손을 N$_2$에 몰면 장부는 맞아도 종이 틀린다#

$\tau_k$를 화학종마다 다르게 주는 순간 길이 셋이다. 결손을 그대로 두는 길, 합으로 나눠 한 종에 몰아넣는 길, 보정속도로 $Y_k$ 비율로 나누는 길. 그대로 두는 길은 없다. 몰아넣는 길은 잉여종이 많을 때만 버틴다. 희석 연소의 $Y_{N_2}\approx 0.7$이 그 조건이다. 기본은 보정속도다. 비용은 스텝마다 기울기 한 벌과 평형분포 속도 한 줄이다.

$\tau$의 범위도 같이 잠긴다. 표에서 $\tau=0.55$의 확산계수 오차는 $+1.243\%$다. $\tau=2.00$에서는 $-0.373\%$로 부호가 바뀐다. 그 바닥은 $\tau=0.51$에서 점성을 51배 틀린 반 스텝과 같은 자리다. 화학종 축으로 그 바닥이 한 번 더 깔린다. 로그에 $\sum_k Y_k=1.0003$이 찍히면 $D_k$부터 나란히 놓는다. 값이 다르면 원인은 그 줄에 있다.