τだけ変えたら静止した気体の密度が変わった — Shan–Chen速度シフトforcingの(τ−½)²応力

平衡速度を τF/ρ だけずらす forcing は運動量フラックスに (τ−½)²F²/ρ を加えるため、静止状態の密度分布が緩和時間によって変わります。

粘性を変えただけで静止状態が変わった

静止した気体は粘性を知りません。速度がゼロなら粘性応力もゼロです。したがって、力と圧力がつり合った静止状態は粘性と無関係なはずです。格子ボルツマン法（LB, Lattice Boltzmann）では粘性を緩和時間 $\tau$ が決めます。だから $\tau$ を変えても静止状態の密度分布は変わらないはずです。

ところが、そうならない forcing 方式があります。Shan–Chen 系の多相モデルが使う速度シフト（velocity shift）forcing です。同じ力、同じ格子で $\tau$ だけを 0.8 から 3 に変えると、密度分布が動きます。この記事では三つの forcing を並べ、その差がどこから来るかを比較します。読み終えると、各方式が運動量フラックスにこっそり加える応力の係数が一行で手に入ります。Guo は 0、Shan–Chen は $(\tau-\tfrac12)^2$ 、EDM は $\tfrac14$ です。

比較基準1：平衡速度をどこへずらすか#

三方式とも、1ステップで運動量 $F\,\delta t$ を加えます。違うのは、その運動量をどの分布に載せるかです。以下で $u'=\sum_i f_i c_i/\rho$ は分布から直接読んだ速度です。物理速度は三方式とも $u = u' + F/(2\rho)$ と定義します（半ステップ補正）。

方式	平衡に入れる速度	追加ソース
Guo	$u$	$(1-\tfrac{1}{2\tau})S_i(u,F)$
Shan–Chen	$u' + \tau F/\rho$	なし
EDM (Kupershtokh)	$u'$	$f_i^{eq}(u'+F/\rho)-f_i^{eq}(u')$

静止状態では $u=0$ なので $u'=-F/(2\rho)$ です。すると Shan–Chen が平衡に入れる速度は次のようになります。

u_{eq} = u' + \frac{\tau F}{\rho} = \left(\tau - \tfrac12\right)\frac{F}{\rho}

$u_{eq}$ は平衡分布を作るときの速度、 $\tau$ は無次元緩和時間、 $F$ は単位体積あたりの力です。気体は止まっているのに、平衡は $\tau$ に比例する速度で動いていることになります。

平衡分布の2次モーメント（運動量フラックス）は速度の放物線です。

\Pi_{eq}(u_{eq}) = \sum_i f_i^{eq} c_i c_i = c_s^2\rho + \rho\, u_{eq}^2

$c_s^2=1/3$ は格子音速の2乗です。 $u_{eq}$ がゼロからずれると、フラックスが $\rho u_{eq}^2$ だけ持ち上がります。下のシミュレーションで実際に操作してみましょう。カーソルが柱に沿って動き、その位置の力 $F(y)$ で各方式の平衡点を放物線上に打ちます。

τ = 1.50 force amplitude G = 0.120

y = 0.0 · F = 0.0000 · SC u_eq = (τ−½)F/ρ = 0.0000 · SC extra 0.00e+0 · EDM extra 0.00e+0 · Guo 0

$\tau$ スライダーを 0.5 付近から 3 まで上げてみてください。オレンジの点（Shan–Chen）は放物線を登りますが、緑の点（Guo）は底にとどまります。紫の点（EDM）は $\pm F/(2\rho)$ の二か所に座り、 $\tau$ と無関係に同じ高さ $F^2/(4\rho)$ を保ちます。

比較基準2：その持ち上がりは静止状態まで残るか#

平衡が持ち上がっても、巨視的な応力がすぐ持ち上がるとは限りません。実際に運動量を運ぶのは衝突直後の分布の2次モーメント $\Pi^*$ です。定常状態で衝突前後の2次モーメントが等しいとおくと、次が得られます。

\Pi^* = \Pi_{eq} + \tau\, S_\Pi

$S_\Pi$ は追加ソースの2次モーメントです。Guo のソースは $u=0$ で2次モーメントがゼロです。EDM のソースは $\rho(u'+F/\rho)^2-\rho u'^2 = 0$ です。つまり三方式の差は、平衡に入れた速度だけから生じます。静止した柱の力のつり合いはこう変わります。

\frac{d}{dy}\left(c_s^2\rho + b\,\frac{F^2}{\rho}\right) = F, \qquad b_{Guo}=0,\quad b_{SC}=\left(\tau-\tfrac12\right)^2,\quad b_{EDM}=\tfrac14

$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$ をフィットし直します。

τ = 1.50 G = 0.020 steps/frame 10

step 0 · b: Guo 0.000 · SC 0.000 (theory 1.000) · EDM 0.000 · max|ρ_SC−ρ_Guo| 0.00e+0 · Δρ/step 1.0e+0

既定値 $\tau=1.5$ 、 $G=0.02$ では、オレンジの曲線が 1.000、紫が 0.250、緑が 0.000 に張り付きます。出力ブロックの3行目と同じです。 $\tau$ を 3 に上げて Reset すると、オレンジだけが 6.25 まで上がります。 $G$ を変えると点の大きさは $G^2$ で変わりますが、 $b$ は変わりません。

出力でまず見るべきは $b_{SC}$ の列と $(\tau-\tfrac12)^2$ の列です。5行すべて小数第3位まで一致します。導出で捨てたものがないということです。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$ に比例するからです。

三方式を一枚で

	Guo	Shan–Chen	EDM
平衡速度（静止時）	$0$	$(\tau-\tfrac12)F/\rho$	$-F/(2\rho)$
偽の応力の係数 $b$	$0$	$(\tau-\tfrac12)^2$	$\tfrac14$
$\tau$ 依存	なし	2次	なし
静止状態が粘性と無関係か	はい	いいえ	はい
静水圧平衡と一致	はい	$\tau=\tfrac12$ 付近のみ	いいえ（ $\tfrac14$ 固定）

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$ がゼロでない定数なら、その定数が界面圧力に上乗せされる量です。32セル、数千ステップで終わる確認です。