열전도 τ만 따로 줬더니 점성 가열이 Pr=1에 갇혔다 — 2-population LB의 준평형 g*

Pr을 바꿨는데 왜 온도는 그대로일까#

Couette 채널을 하나 생각해보자. 아래 벽은 멈춰 있고 위 벽은 속도 $U$로 미끄러진다. 두 벽 온도는 같다. 그래도 가운데는 데워진다. 점성이 운동 에너지를 열로 바꾸기 때문이다(점성 가열, viscous heating). 가운데 온도 상승은 Prandtl 수(Pr, 운동량 확산과 열 확산의 비)에 정확히 비례한다.

이 글은 격자 볼츠만(LB, Lattice Boltzmann) 열 모델에서 Pr을 바꾸는 방법을 다룬다. 분포함수를 두 묶음으로 나누고 이완 시간을 따로 주면 Pr은 쉽게 바뀌는 것처럼 보인다. 실제로 돌려보면 열전도는 바뀌는데 점성 가열은 Pr=1 값에서 움직이지 않는다. 독자는 그 이유와 한 줄짜리 해법인 준평형(quasi-equilibrium) $g^*$를 얻어 간다.

Pr 슬라이더를 움직이면 응력 시계와 열전도 시계의 간격이 바뀌고, 점성 일 플럭스가 어느 곡선에 붙는지도 같이 바뀐다.

Pr 슬라이더를 0.71로 두면 주황(응력)과 초록 점선(g* 모델의 점성 일 플럭스)은 끝까지 겹친다. 빨강(단순 2-τ 모델)은 보라(열전도) 시계를 따라가서 응력과 갈라진다. Pr=1로 옮기면 네 선이 둘로 합쳐진다.

한 묶음의 분포함수는 Pr=1밖에 모른다#

BGK 충돌 연산자는 모든 비평형 모멘트를 같은 속도 $\omega$로 평형에 끌어당긴다. 응력도, 열 플럭스도 같은 시계로 이완된다. 그래서 점성 계수와 열전도 계수가 같은 $\omega$에 묶인다. 연속 BGK 방정식에서 Pr=1이 나오는 이유가 이것이다.

공기의 Pr은 0.71, 물은 약 7이다. Pr=1에 묶인 모델로는 경계층 두께 비조차 맞지 않는다. 처방은 두 갈래다. 하나는 충돌 연산자에서 열 플럭스 모멘트만 다른 속도로 이완하는 것이다(MRT, Shakhov형). 다른 하나는 에너지를 별도의 분포함수에 맡기는 것이다. 이 글은 두 번째 길을 따라간다.

두 묶음: f는 질량·운동량, g는 총에너지#

2-population 모델은 $f_i$가 밀도와 운동량을, $g_i$가 총에너지 $E=\rho(c_vT+u^2/2)$를 나른다($i=0,\dots,Q-1$). 운동 방정식은 다음과 같다.

$\omega_\mu$는 점성(동점성계수 $\nu$), $\omega_\kappa$는 열전도(열확산계수 $\alpha$)를 정하는 이완 계수다. $g_i^{eq}$는 국소 평형, $g_i^*$는 준평형 분포다.

두 계수는 확산 계수와 이렇게 이어진다.

$c_s^2=1/3$은 격자 음속의 제곱이다. $\tau_\mu=1/\omega_\mu$, $\tau_\kappa=1/\omega_\kappa$로 쓰면 $\mathrm{Pr}=(\tau_\mu-\tfrac12)/(\tau_\kappa-\tfrac12)$다.

코드는 $\tau_\mu=1$로 고정한다. 그러면 $\tau_\mu-\tfrac12=\tfrac12$라서 $\tau_\kappa=\tfrac12+(\tfrac12)/\mathrm{Pr}$이다. 표의 네 칸이면 이렇다. Pr=0.5는 $\tau_\kappa=1.5$. Pr=0.71은 $0.5/0.71\approx0.704$를 더해 $\tau_\kappa\approx1.204$. Pr=1은 $\tau_\kappa=1$. Pr=2는 $\tau_\kappa=0.75$. $\omega=1/\tau$이므로 Pr=2의 열 이완은 $\omega_\kappa\approx1.333$, 점성은 $\omega_\mu=1$이다. 점성 일은 느린 쪽에 있어야 한다. 단순 2-τ는 그 일을 빠른 $\omega_\kappa$에 싣는다. Pr=1 행만 두 $\omega$가 같아서 두 열 모두 1.000이다. 그 행만 보면 버그가 조용하다.

에너지 플럭스 안에 점성 일이 숨어 있다#

이 글의 장난감에서 $g^{eq}$의 모멘트는 셋만 정한다. 0차는 $E$, 1차(에너지 플럭스)는 $(E+p)\mathbf{u}$, 2차는 $c_s^2E\,\mathbf{I}+E\,\mathbf{u}\mathbf{u}$다. 1차 모멘트의 비평형 부분 $\mathbf{Q}$가 실제 열 플럭스다.

Chapman–Enskog 전개(작은 Knudsen 수로 비평형을 푸는 전개) 1차에서 $\mathbf{Q}\approx-\tau\,\nabla\!\cdot\!\mathbf{R}^{eq}$다. 여기서 $\nabla(c_s^2E)$가 두 조각으로 갈라진다.

첫 항은 Fourier 열전도다. 뒤쪽 항이 점성 일 $-\mathbf{u}\cdot\boldsymbol{\sigma}$의 재료다. 이 항은 $f$의 비평형 응력 $\boldsymbol{\Pi}^{neq}=\mathbf{P}-\mathbf{P}^{eq}\approx-\tau_\mu\rho c_s^2(\nabla\mathbf{u}+\nabla\mathbf{u}^T)$와 똑같은 모양이다. 문제는 앞에 붙는 $\tau$다.

$g$를 평형으로만 $\omega_\kappa$로 끌어당기면 두 조각 모두 $\tau_\kappa$를 달고 나온다. 점성 일의 계수가 $\nu$가 아니라 $\alpha$에 비례하게 된다. 열전도는 맞는데 점성 가열은 Pr=1일 때의 값으로 남는다. 앞 viz의 빨간 선이 이것이다.

준평형 g*: 플럭스만 응력의 시계에 묶는다#

해법은 $g^*$의 정의 한 줄이다. $g^*$는 $g^{eq}$와 모든 모멘트가 같고, 에너지 플럭스만 다르다.

$\mathbf{P}$는 $f$에서 잰 실제 운동량 플럭스, $\mathbf{P}^{eq}$는 그 평형값이다.

위 운동 방정식의 1차 모멘트를 취하면 두 비평형 모멘트의 한 스텝 갱신이 나온다.

$\mathbf{Q}$ 자체는 $\omega_\kappa$로 이완한다. 그래서 열전도는 $\alpha$를 따른다. 동시에 $\mathbf{Q}=\mathbf{u}\cdot\boldsymbol{\Pi}^{neq}$이면 다음 스텝에도 $\mathbf{Q}'=\mathbf{u}\cdot\boldsymbol{\Pi}'^{neq}$가 유지된다. 점성 일 조각은 응력의 시계 $\omega_\mu$를 탄다. Chapman–Enskog으로 계산해도 $\tau_\kappa$가 붙은 점성 일이 $(\omega_\kappa-\omega_\mu)$ 항과 상쇄되고 $\tau_\mu$만 남는다.

$\omega_\kappa=\omega_\mu$이면 $g^*$ 항이 사라진다. 이때는 단일 이완 모델과 같고 Pr=1이다. $g^*$는 Pr$\ne1$일 때만 일하는 보정이다. 한 분포만 쓰는 열 격자에서도 τ를 올려도 Pr이 1에 붙어 있었다. 그 두 상수가 잠긴 자리와 같은 벽이다. 분포를 둘로 쪼개도, 점성 일의 시계를 옮기지 않으면 그 벽에 남는다. 질문은 하나다. 어떤 모멘트를 어떤 속도로 이완시킬 것인가.

코드: 같은 Couette 채널, 두 가지 g#

해석해는 간단하다. 정상 상태 에너지식은 $kT''+\mu(u')^2=0$이다. 속도가 선형이면 온도는 포물선이 된다.

가운데 최댓값은 $\theta_{max}=\mathrm{Pr}$이다. 아래 코드는 D2Q9를 $x$ 방향 1칸 주기로 두고 $y$로만 푼다. 벽은 비평형 외삽(Guo) 경계다. 단순화를 위해 $f$는 등온 평형을 쓰고 에너지는 한 방향으로만 결합한다. 저마하($U=0.1$) Couette에서는 압력이 균일해서 이 단순화가 결과를 바꾸지 않는다.

import numpy as np
 
# D2Q9, x 방향 1칸 주기 → y 방향 1D 배열
cx = np.array([0, 1, 0, -1, 0, 1, -1, -1, 1])
cy = np.array([0, 0, 1, 0, -1, 1, 1, -1, -1])
w = np.array([4/9] + [1/9]*4 + [1/36]*4)
CV, T0 = 1.0, 1.0
 
def feq_lattice(rho, ux):
    cu = cx[:, None]*ux
    return w[:, None]*rho*(1 + 3*cu + 4.5*cu**2 - 1.5*ux**2)
 
def geq_energy(E, rho, ux):
    p = rho*(E/rho - 0.5*ux**2)/CV          # p = rho*T
    qx = (E + p)*ux
    cu = cx[:, None]*ux
    return w[:, None]*(E + 3*cx[:, None]*qx + 4.5*E*cu**2 - 1.5*E*ux**2)
 
def couette_heating(Pr, quasi, N=21, U=0.1, tau_mu=1.0, steps=8000):
    om_mu = 1/tau_mu
    om_ka = 1/(0.5 + (tau_mu - 0.5)/Pr)
    rho = np.ones(N); ux = np.zeros(N)
    E = rho*CV*T0*np.ones(N)
    f = feq_lattice(rho, ux); g = geq_energy(E, rho, ux)
    uw = np.zeros(N); uw[-1] = U
    for _ in range(steps):
        rho = f.sum(0); ux = (cx[:, None]*f).sum(0)/rho; E = g.sum(0)
        fe = feq_lattice(rho, ux); ge = geq_energy(E, rho, ux)
        Pxy = ((cx*cy)[:, None]*(f - fe)).sum(0)
        f = f + om_mu*(fe - f)
        if quasi:   # g*: 에너지 플럭스에 u·(P - P_eq)를 더한 준평형
            gs = ge + w[:, None]*3*cy[:, None]*(ux*Pxy)
            g = g + om_mu*(ge - g) + (om_ka - om_mu)*(gs - g)
        else:       # 단순 2-τ: g는 평형으로만 ω_κ
            g = g + om_ka*(ge - g)
        for k in range(9):
            f[k] = np.roll(f[k], cy[k]); g[k] = np.roll(g[k], cy[k])
        for wall, nb in ((0, 1), (N-1, N-2)):   # 비평형 외삽 벽
            r = f[:, nb].sum(); un = (cx*f[:, nb]).sum()/r; En = g[:, nb].sum()
            fn = feq_lattice(np.array([r]), np.array([un]))[:, 0]
            gn = geq_energy(np.array([En]), np.array([r]), np.array([un]))[:, 0]
            Ew = r*(CV*T0 + 0.5*uw[wall]**2)
            f[:, wall] = feq_lattice(np.array([r]), np.array([uw[wall]]))[:, 0] + f[:, nb] - fn
            g[:, wall] = geq_energy(np.array([Ew]), np.array([r]), np.array([uw[wall]]))[:, 0] + g[:, nb] - gn
    rho = f.sum(0); ux = (cx[:, None]*f).sum(0)/rho
    T = (g.sum(0)/rho - 0.5*ux**2)/CV
    theta = (T - T0)/(U**2/(8*CV))
    return theta.max(), ux
 
print(" Pr   exact   2-tau(g_eq)   g* quasi-eq")
for Pr in (0.5, 0.71, 1.0, 2.0):
    a, _ = couette_heating(Pr, False)
    b, u = couette_heating(Pr, True)
    print(f"{Pr:4.2f}  {Pr:6.3f}   {a:9.3f}    {b:9.3f}")
print(f"u(mid)/U = {u[10]/0.1:.4f}")

 Pr   exact   2-tau(g_eq)   g* quasi-eq
0.50   0.500       1.000        0.498
0.71   0.710       1.000        0.710
1.00   1.000       1.000        1.000
2.00   2.000       1.000        1.997
u(mid)/U = 0.5000

위 코드를 브라우저에서 직접 돌려보자. 두 모델이 같은 채널에서 나란히 시간 전진한다.

Pr 슬라이더를 움직이면 초록 곡선의 정점은 점선 해석해 Pr을 따라가고, 빨강은 1.000에 붙어 있다. $U$나 $\tau_\mu$를 바꿔도 빨강의 $\theta_{max}$는 1에서 떠나지 않는다. 출력 블록에 없는 Pr=3에서도 확인해보자.

1.000이라는 숫자가 말하는 것#

단순 2-τ 열은 Pr과 무관하게 1.000이다. 오차가 아니라 구조다. 점성 일 계수가 $\tau_\kappa-\tfrac12$에, 열전도 계수도 $\tau_\kappa-\tfrac12$에 비례한다. 둘의 비인 $\theta_{max}$에서 $\tau_\kappa$가 정확히 약분된다. 공기(Pr=0.71)라면 벽 사이 온도 상승을 41% 과대 예측한다. 물 같은 Pr=7 유체라면 7배 과소 예측이다.

$U=0.1$, $c_v=1$이면 분모는 $U^2/(8c_v)=0.00125$다. 공기의 참 $\Delta T$는 $0.71\times0.00125\approx8.9\times10^{-4}$, 단순 2-τ는 $1.25\times10^{-3}$이다. 격자 온도는 1.0009와 1.0013. 소수점 두 자리면 둘 다 1.00이다. Pr=7이면 참값 $0.00875$, 모델 $0.00125$. 1.009 대 1.001. $\theta$로 봐야 1과 7이 갈린다.

$g^*$ 열은 0.498, 0.710, 1.000, 1.997이다. 0.498은 0.500보다 0.4% 작고, 1.997은 2보다 0.15% 작다. 0.710은 exact와 같다. 격자 21칸에서 남은 차이는 벽과 해상도다. 가운데 속도 $u/U=0.5000$은 두 모델이 같은 $f$를 쓰므로 똑같다. 온도 오류는 속도장에 안 보인다.

압축성으로 가면 f에도 준평형이 붙는다#

이 글의 $f$는 등온 평형을 썼다. 압축성 2-population 모델에서는 $f$ 쪽에도 준평형 $f^*$가 들어간다. D2Q9의 3차 모멘트에는 $c_{ix}^3=c_{ix}$라는 결함이 있다. 온도가 변하거나 Mach 수가 커지면 응력에 오차 항이 생긴다. $f^*$는 그 결함을 보정하는 압력 텐서를 담는다. 운동량 플럭스를 바로잡는 쪽이 $f^*$, 에너지 플럭스에 점성 일을 붙이는 쪽이 $g^*$다.

원리는 같다. 보존량(밀도·운동량·총에너지)은 평형과 준평형이 공유한다. 다른 것은 이완 속도를 따로 줘야 하는 모멘트 하나뿐이다. 비이상 유체의 압력 항을 어디에 넣을지도, 결국 그 모멘트의 자리다.

열 LB 코드에서 τκ를 따로 줄 때 확인할 한 줄#

에너지 분포의 충돌식을 열어 $\omega_\kappa$가 어디에 곱해지는지 본다. $g$ 전체가 $\omega_\kappa$로 평형에 끌려가면 점성 가열은 Pr=1 값이다. 확인은 Couette 한 번이면 된다. $\theta_{max}$가 Pr을 따라가지 않으면 준평형 항이 빠진 것이다.

반대로 $\omega_\kappa=\omega_\mu$로만 검증했다면 이 버그는 조용히 남는다. 열 LB 검증에는 Pr$\ne1$ 점성 가열을 하나 넣는다. 정지 유체의 온도 확산만으로는 이 항이 돌지 않는다. 마찰과 열을 한 비로 묶는 레이놀즈 상사도, 그 비가 Pr에 매여 있다는 전제에서 시작한다. $\theta_{max}$가 Pr을 따라가지 않으면 그 전제가 격자 안에서 이미 깨져 있다.