熱伝導のτだけ分けたら粘性加熱が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)$が二つに分かれます。

第1項はFourierの熱伝導です。第2項が粘性仕事$-\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次モーメントを取ると、二つの非平衡モーメントの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分の1に過小予測します。

$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}$という欠陥があります。温度が変わったりマッハ数が大きくなったりすると、応力に誤差項が生じます。$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 に従わないなら、その前提は格子の中ですでに壊れています。