τ만 바꿨는데 정지한 기체의 밀도가 바뀌었다 — Shan–Chen 속도 이동 forcing의 (τ−½)² 응력

점성만 바꿨는데 정지 상태가 바뀌었다#

정지한 기체는 점성을 모른다. 속도가 0이면 점성 응력도 0이다. 그래서 힘과 압력이 평형을 이룬 정지 상태는 점성과 무관해야 한다. 격자 볼츠만(LB, Lattice Boltzmann)에서 점성은 이완 시간 $\tau$가 정한다. 따라서 $\tau$를 바꿔도 정지 상태의 밀도 분포는 그대로여야 한다.

그런데 그렇지 않은 forcing 방식이 있다. Shan–Chen 계열 다상 모델이 쓰는 속도 이동(velocity shift) forcing이다. 같은 힘, 같은 격자에서 $\tau$만 0.8에서 3으로 바꾸면 밀도 분포가 달라진다. 이 글은 그 차이가 어디서 오는지 세 가지 forcing을 나란히 놓고 비교한다. 독자는 forcing마다 운동량 플럭스에 몰래 더해지는 응력의 계수를 한 줄로 얻어 간다. Guo는 0, Shan–Chen은 $(\tau-\tfrac12)^2$, EDM은 $\tfrac14$이다.

비교 기준 1: 평형 속도를 어디로 옮기는가#

세 방식 모두 한 스텝에 운동량 $F\,\delta t$를 더한다. 다른 점은 그 운동량을 어느 분포에 싣느냐다. 아래에서 $u'=\sum_i f_i c_i/\rho$는 분포에서 바로 읽은 속도다. 물리 속도는 세 방식 모두 $u = u' + F/(2\rho)$로 정의한다(반 스텝 보정).

정지 상태에서는 $u=0$이므로 $u'=-F/(2\rho)$다. 그러면 Shan–Chen이 평형에 넣는 속도는 다음과 같다.

$u_{eq}$는 평형 분포를 만들 때 쓰는 속도, $\tau$는 무차원 이완 시간, $F$는 단위 부피당 힘이다. 기체는 멈춰 있는데 평형은 $\tau$에 비례하는 속도로 움직이는 셈이다.

평형 분포의 2차 모멘트(운동량 플럭스)는 속도의 포물선이다.

$c_s^2=1/3$은 격자 음속의 제곱이다. $u_{eq}$가 0에서 벗어나면 플럭스가 $\rho u_{eq}^2$만큼 들린다. 아래 시뮬레이션에서 직접 조작해보자. 기둥을 따라 커서가 움직이며 그 위치의 힘 $F(y)$로 각 방식의 평형점을 포물선 위에 찍는다.

$\tau$ 슬라이더를 0.5 근처에서 3으로 올려보자. 주황 점(Shan–Chen)은 포물선을 타고 올라가지만 초록 점(Guo)은 바닥에 머문다. 보라 점(EDM)은 $\pm F/(2\rho)$ 두 곳에 앉아 $\tau$와 무관하게 같은 높이 $F^2/(4\rho)$를 유지한다.

비교 기준 2: 그 들림이 정지 상태까지 살아남는가#

평형이 들렸다고 거시 응력이 바로 들리는 것은 아니다. 실제로 운동량을 옮기는 것은 충돌 직후 분포의 2차 모멘트 $\Pi^*$다. 정상 상태에서 충돌 전후의 2차 모멘트가 같다고 놓으면 다음을 얻는다.

$S_\Pi$는 추가 소스의 2차 모멘트다. Guo의 소스는 $u=0$에서 2차 모멘트가 0이다. EDM의 소스는 $\rho(u'+F/\rho)^2-\rho u'^2 = 0$이다. 그러니 세 방식의 차이는 평형에 넣은 속도에서만 나온다. 정지 기둥의 힘 평형은 이렇게 바뀐다.

$b$는 가짜 응력의 계수다. 연속체의 정수압 평형은 $b=0$이다.

이 항은 왜 교과서 유도에서 보이지 않을까. 표준 Chapman–Enskog 전개(작은 매개변수 $\epsilon$로 분포를 차수별로 나누는 다중 척도 전개)는 $F$를 1차 작은 양으로 둔다. 그러면 $F^2/\rho$는 $\epsilon^2$ 항이고 Navier–Stokes 수준에서 버려진다. 하지만 계수가 $\tau^2$로 커지면 "작다"는 보장이 사라진다. 소스 노트가 인용한 Revised Chapman–Enskog analysis for a class of forcing schemes in the lattice Boltzmann method가 바로 힘 항의 차수 배정을 다시 따지는 논문이다. 힘을 압력 기울기 형태로 넣을지 힘 형태로 넣을지는 비이상 LBM forcing의 두 형태에서 다뤘다. 이 글은 같은 힘을 어떤 방식으로 분포에 싣느냐를 본다.

코드: 주기 기둥, 세 가지 forcing, 다섯 개의 τ#

D1Q3 격자(속도 0, ±1) 32칸 주기 기둥에 $F = G\sin(ky)$를 건다. 처음엔 균일 밀도다. 음파가 점성으로 가라앉으면 정지 상태가 남는다. 수렴한 밀도에서 $c_s^2\rho = \lambda P(y) - bF^2/\rho + c$를 최소제곱으로 맞춰 $b$를 잰다. $P(y)=-\tfrac{G}{k}\cos(ky)$는 $P'=F$인 정수압 퍼텐셜이다. $\lambda$는 격자 이산화가 진폭에 주는 배율을 흡수한다.

import math
 
# D1Q3 격자: 속도 0, +1, -1 ; cs^2 = 1/3
C = (0.0, 1.0, -1.0)
W = (2/3, 1/6, 1/6)
CS2 = 1/3
 
def d1q3_eq(rho, u):
    return [W[i]*rho*(1 + C[i]*u/CS2 + 0.5*(C[i]*u)**2/CS2**2 - 0.5*u*u/CS2)
            for i in range(3)]
 
def collide_node(fn, F, tau, scheme):
    rho = sum(fn)
    mom = fn[1] - fn[2]
    if scheme == "guo":                 # Guo: u = (sum f c + F/2)/rho, 소스에 (1 - 1/2tau) 곱
        u = (mom + 0.5*F)/rho
        feq = d1q3_eq(rho, u)
        out = []
        for i in range(3):
            S = W[i]*((C[i] - u)/CS2 + C[i]*u*C[i]/CS2**2)*F
            out.append(fn[i] - (fn[i] - feq[i])/tau + (1 - 0.5/tau)*S)
        return out
    if scheme == "sc":                  # Shan-Chen: 평형 속도를 tau F / rho 만큼 이동
        feq = d1q3_eq(rho, mom/rho + tau*F/rho)
        return [fn[i] - (fn[i] - feq[i])/tau for i in range(3)]
    up = mom/rho                        # EDM: feq(u'+F/rho) - feq(u') 를 더함
    feq = d1q3_eq(rho, up)
    fsh = d1q3_eq(rho, up + F/rho)
    return [fn[i] - (fn[i] - feq[i])/tau + fsh[i] - feq[i] for i in range(3)]
 
def run_column(scheme, tau, G, N=32, tol=1e-12):
    k = 2*math.pi/N
    F = [G*math.sin(k*(j + 0.5)) for j in range(N)]
    f = [d1q3_eq(1.0, 0.0) for _ in range(N)]
    for it in range(1, 100001):
        post = [collide_node(f[j], F[j], tau, scheme) for j in range(N)]
        new = [[post[j][0], post[j-1][1], post[(j+1) % N][2]] for j in range(N)]
        if it % 100 == 0 and max(abs(sum(new[j]) - sum(f[j])) for j in range(N)) < tol:
            return F, [sum(fn) for fn in new], it
        f = new
    return F, [sum(fn) for fn in f], it
 
def stress_fit(F, rho, G):
    # 맞춤  cs^2 rho = lam*P(y) - b*F^2/rho + c,  P' = F  (정수압이면 b = 0)
    N = len(F); k = 2*math.pi/N
    X = [[F[j]**2/rho[j], -G/k*math.cos(k*(j + 0.5)), 1.0] for j in range(N)]
    r = [CS2*x for x in rho]
    A = [[sum(X[n][a]*X[n][b] for n in range(N)) for b in range(3)] for a in range(3)]
    v = [sum(X[n][a]*r[n] for n in range(N)) for a in range(3)]
    for i in range(3):
        for j in range(i + 1, 3):
            m = A[j][i]/A[i][i]
            A[j] = [A[j][c] - m*A[i][c] for c in range(3)]
            v[j] -= m*v[i]
    x = [0.0]*3
    for i in (2, 1, 0):
        x[i] = (v[i] - sum(A[i][c]*x[c] for c in range(i + 1, 3)))/A[i][i]
    return -x[0], x[1]
 
G = 0.02
print(f"periodic column N=32, F = {G} sin(ky), rho0 = 1")
print(" tau | b_Guo  b_SC   (tau-1/2)^2  b_EDM | lam   | max|rho_SC-rho_Guo|  iters")
for tau in [0.8, 1.0, 1.5, 2.0, 3.0]:
    out = {s: run_column(s, tau, G) for s in ("guo", "sc", "edm")}
    b = {s: stress_fit(out[s][0], out[s][1], G) for s in out}
    drho = max(abs(p - q) for p, q in zip(out["sc"][1], out["guo"][1]))
    bg = round(b["guo"][0], 3) + 0.0
    print(f"{tau:4.1f} | {bg:5.3f} {b['sc'][0]:6.3f}   {(tau-0.5)**2:6.3f}    {b['edm'][0]:5.3f} |"
          f" {b['guo'][1]:.4f} |  {drho:.2e}            {out['sc'][2]}")

periodic column N=32, F = 0.02 sin(ky), rho0 = 1
 tau | b_Guo  b_SC   (tau-1/2)^2  b_EDM | lam   | max|rho_SC-rho_Guo|  iters
 0.8 | 0.000  0.090    0.090    0.250 | 0.9968 |  5.49e-05            6100
 1.0 | 0.000  0.250    0.250    0.250 | 0.9968 |  1.53e-04            3600
 1.5 | 0.000  1.000    1.000    0.250 | 0.9968 |  6.11e-04            1900
 2.0 | 0.000  2.250    2.250    0.250 | 0.9968 |  1.38e-03            1300
 3.0 | 0.000  6.250    6.250    0.250 | 0.9968 |  3.84e-03            800

0, (τ−½)², 1/4가 말하는 것#

위 코드를 브라우저에서 직접 돌려보자. 세 기둥이 같은 힘 아래 동시에 진행하고, 매 프레임 $b$를 다시 맞춘다.

기본값 $\tau=1.5$, $G=0.02$에서 주황 곡선은 1.000, 보라는 0.250, 초록은 0.000에 붙는다. 출력 블록의 셋째 줄과 같다. $\tau$를 3으로 올리고 Reset하면 주황만 6.25로 올라간다. $G$를 바꾸면 점의 크기는 $G^2$로 변하지만 $b$는 그대로다.

출력에서 먼저 볼 것은 $b_{SC}$ 열과 $(\tau-\tfrac12)^2$ 열이다. 다섯 줄 모두 소수 셋째 자리까지 같다. 유도에서 버린 것이 없다는 뜻이다. Guo의 0.000도 중요하다. 정수압 평형이 $F^2$ 수준까지 맞는다.

$\tau=1$에서는 Shan–Chen과 EDM이 똑같이 0.250이다. 우연이 아니다. $\tau=1$이면 충돌 후 분포가 곧 평형이다. Shan–Chen은 $f^{eq}(u'+F/\rho)$가 되고, EDM도 $f^{eq}(u') + f^{eq}(u'+F/\rho) - f^{eq}(u')$로 같은 값이 된다. 두 방식은 $\tau=1$에서 한 점으로 만난다.

$\lambda=0.9968$은 forcing과 무관하다. 32칸 격자가 사인파 진폭에 주는 이산화 배율이고, 값은 $1-k^2/12$와 일치한다. 격자를 늘리면 1로 간다. $b$는 이 배율과 따로 잡히도록 맞췄다.

밀도 차이 자체는 작아 보인다. $\tau=3$에서 $3.84\times10^{-3}$이고 밀도 진폭 0.30의 1.3%다. 그런데 이 비율은 $G$에 비례해 커진다. 주 효과는 $F$에 비례하고 가짜 응력은 $F^2$에 비례하기 때문이다.

한 장으로 본 세 방식#

EDM은 $\tau$와 무관하지만 $\tfrac14$만큼 어긋난다. 어긋남이 상수라서 보정하기 쉽다. Shan–Chen은 $\tau=0.5$ 근처에서만 작다. 그 영역은 LB가 가장 불안정해지는 곳이다.

다상 계면에서는 F가 크다#

이 실험의 힘은 매끄럽고 약했다. 의사퍼텐셜(pseudopotential) 다상 모델의 힘은 다르다. 상호작용 힘 $F=-\mathcal{G}\psi\nabla\psi$($\mathcal{G}$는 상호작용 세기, $\psi$는 유효 밀도)는 계면 몇 칸 안에 몰려 있고 크다. 거기서 $bF^2/\rho$는 상태방정식이 만드는 압력과 같은 크기까지 자란다.

그래서 Shan–Chen forcing으로 평면 계면을 세우면 공존 밀도(coexistence density, 두 상이 평형일 때 각 상의 밀도)가 $\tau$에 따라 움직인다. 점성을 바꾸려고 $\tau$를 건드렸는데 액체 밀도가 바뀌는 셈이다. 이 기둥 실험은 그 현상을 계면 없이 분리해 보여준다. 가짜 응력은 계면의 기계적 평형 조건 안으로 그대로 들어간다.

반대로 이 항을 일부러 쓰는 방법도 있다. $b$를 조절 가능한 매개변수로 두면 공존 밀도를 열역학 일관성에 맞추는 손잡이가 된다. EDM의 $\tfrac14$나 Guo의 0 위에 원하는 만큼 더하는 식이다. 어떤 경우든 먼저 내 코드의 $b$가 얼마인지 알아야 한다. 같은 Chapman–Enskog 틀로 열 모델의 숨은 항을 찾은 예는 2-population LB의 준평형 g*에 있다.

새 forcing을 넣기 전에 돌려볼 32칸#

새 forcing 방식을 코드에 넣을 때 다상 문제부터 돌리지 말자. 계면에서는 원인이 섞인다. 대신 위의 주기 기둥을 먼저 돌린다. $\tau$를 두 값으로 바꿔 $b$를 잰다. $b$가 $\tau$에 따라 움직이면 속도 이동이 평형 안에 숨어 있다. $b$가 0이 아닌 상수라면 그 상수가 계면 압력에 더해질 양이다. 32칸, 몇 천 스텝이면 끝나는 확인이다.