3次モーメントを正直に埋めたら誤差が2.8倍になった — 力項が守るべき行

力項は何次モーメントまで合わせるべきか

LBMのコードに重力や表面張力を入れるとき、私たちが合わせるのは三行です。 0次モーメントは0、1次モーメントは $\mathbf{F}$、2次モーメントは $\mathbf{u}\mathbf{F}+\mathbf{F}\mathbf{u}$。 では3次はどうでしょうか。連続ボルツマン方程式の力項も、3次モーメントに明確な値を持ちます。 その値を格子上でも合わせれば、精度が上がりそうに見えます。

数値を打ってみると、結果は逆でした。 3次モーメントを運動論が要求する値で埋めると、誤差が0.053%から0.150%に増えました。 速度スライダー $u_x$、$u_y$ を上げると、合っていたモーメントの行が割れます。チェックボックスが閉じるのは黄色の行だけです。

速度を0側へ戻すと行は再び揃います。赤い行はスライダーを動かしても閉じません。

連続ボルツマンが要求する力の3次モーメント#

体積力 $\mathbf{F}$ があるボルツマン方程式の力項は、速度空間の微分です。

$\xi$ は分子速度、$\rho$ は密度です。3次モーメントは部分積分一回で閉じます。

ここで $M_{\alpha\beta}=\int f^{eq}\xi_\alpha\xi_\beta\,d\boldsymbol{\xi} = \rho(u_\alpha u_\beta + c_s^2\delta_{\alpha\beta})$ です。 したがって要求される3次モーメントはこう書けます。

$c_s^2$ の項と $uuF$ の項が並んでいます。前者は速度に依らず、後者は速度の2次です。

Guo forcingが実際に供給する値#

格子に載せる力項は通常Guoの形で書きます。半ステップ因子を外し、括弧の中だけを $\hat{S}_i$ と呼びます。

$w_i$ は格子重み、$\mathbf{c}_i$ は離散速度です。実際の $S_i$ にはさらに $(1-1/2\tau)$ が掛かります。 その因子が必要な理由は離散LBEが残す半ステップで別に扱いました。

この $\hat{S}_i$ の3次モーメントをD2Q9の4次等方性で計算すると、残るのは次だけです。

$c_s^2$ の項は正確に合います。$uuF$ の項は一つもありません。第二項は5次モーメントを作り、D2Q9では0になります。

先に速度を上げたとき開いたのは黄色の二行だけで、赤い行は閉じたまま残りました。

$c_x^3 = c_x$ — 格子が固定した二成分#

赤い行の正体は格子そのものです。D2Q9の速度成分は $-1, 0, 1$ しかありません。 すると全方向で $c_{ix}^3 = c_{ix}$ が成り立ちます。二つのベクトルが同じベクトルです。

$\Phi_{xxx}$ は選べる値ではありません。1次モーメントを $F_x$ に決めた瞬間に一緒に決まります。 要求値は $3F_x(u_x^2+c_s^2) = F_x + 3F_x u_x^2$ なので、$3F_x u_x^2$ が恒久的な欠損です。

一方 $\Phi_{xxy}$ と $\Phi_{xyy}$ は独立です。9個の自由度のうち0~2次が6個を使い、3個が残ります。 つまり4成分のうち二つは自由で、二つは固定されています。テンソルの半分だけ直せるということです。

質量も運動量も運ばない項が運動量を動かす

$\Phi_{xyy}$ だけを $d$ だけ変える補正ベクトル $v_i$ を作れます。 0次・1次・2次モーメントがすべて0なので、その場では密度も運動量も応力も触りません。 D2Q9でそのベクトルはこうなります。

ところがストリーミングを一回行うと話が変わります。各成分が別の隣接点へ飛びます。 $x$方向の運動量だけを集めると、次が残ります。

2次中心差分です。$d$ が一様なら厳密に0で、$d$ が振動すれば波数の2乗で大きくなります。 波形が一様なら格子ごとの変化は消え、正弦波で波長を縮めるとその変化が波数の2乗に沿って伸びます。格子の総和は0のままです。

力のパターンを固定して横流を通した流路

検証には厳密解のある問題が必要です。周期領域に $x$方向の正弦体積力を固定し、 $y$方向に一様な横流 $V$ を通します。多孔質壁で片側から吹き込み、反対側から吸い出す状況です。

$\nu$ は動粘性、$k$ は波数、$a_0$ は力の振幅です。単一モードなので閉じた解が出ます。

横流は二つのことをします。振幅を削り、位相を $\arctan\!\left(V/\nu k\right)$ だけ遅らせます。 同時に $u_y = V$ が生きているので $d = F_x u_y^2$ が空間で振動します。前節の漏れが入るスイッチです。

四つの組み合わせを同じ格子にかけた

16ノードの周期格子、$\tau = 1$、$a_0 = 4\times10^{-3}$ で四通り回しました。 標準Guo、力の3次を埋める、平衡分布の3次を埋める、両方です。

import numpy as np
 
CX  = np.array([0, 1, 0,-1, 0, 1,-1,-1, 1], float)
CY  = np.array([0, 0, 1, 0,-1, 1, 1,-1,-1], float)
WT  = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
CS2 = 1.0/3.0
 
def hermite3_vectors():
    """0~2次モーメントは0、3次の xxy / xyy のみ1にする補正ベクトル。"""
    rows = np.vstack([np.ones(9), CX, CY, CX*CX, CX*CY, CY*CY, CX*CX*CY, CX*CY*CY])
    pinv = np.linalg.pinv(rows)
    return pinv @ np.eye(8)[6], pinv @ np.eye(8)[7]
 
VXXY, VXYY = hermite3_vectors()
 
def feq_second_order(rho, ux, uy):
    cu = CX[:, None, None]*ux + CY[:, None, None]*uy
    return WT[:, None, None]*rho*(1 + cu/CS2 + cu*cu/(2*CS2**2)
                                  - (ux*ux + uy*uy)/(2*CS2))
 
def force_populations(ux, uy, fx, fy, tau):
    """Guo forcing: 1次モーメント F、2次モーメント uF+Fu。"""
    cF  = CX[:, None, None]*fx + CY[:, None, None]*fy
    gxx, gyy, gxy = ux*fx, uy*fy, 0.5*(ux*fy + uy*fx)
    quad = (gxx*(CX[:, None, None]**2 - CS2) + gyy*(CY[:, None, None]**2 - CS2)
            + 2*gxy*(CX*CY)[:, None, None])
    return WT[:, None, None]*(1 - 1/(2*tau))*(cF/CS2 + quad/CS2**2)
 
def analytic_wave(ny, nu, V, a0):
    """nu u'' - V u' + a0 sin(ky) = 0 の周期定常解。"""
    k = 2*np.pi/ny
    y = np.arange(ny, dtype=float)
    c1 = a0/(nu*k*k + V*V/nu)
    return c1*np.sin(k*y) - (V*c1/(nu*k))*np.cos(k*y)
 
def run_cross_flow(ny, tau, V, a0, cubic_eq, cubic_force, nstep=6000, tol=1e-12):
    nx = 3
    k  = 2*np.pi/ny
    fx = np.tile(a0*np.sin(k*np.arange(ny)), (nx, 1))
    fy = np.zeros((nx, ny))
    rho = np.ones((nx, ny)); ux = np.zeros((nx, ny)); uy = np.full((nx, ny), V)
    f = feq_second_order(rho, ux, uy)
    half, prev = 1 - 1/(2*tau), ux.copy()
    for s in range(nstep):
        rho = f.sum(0)
        ux  = ((CX[:, None, None]*f).sum(0) + 0.5*fx)/rho
        uy  = ((CY[:, None, None]*f).sum(0) + 0.5*fy)/rho
        feq = feq_second_order(rho, ux, uy)
        if cubic_eq:                      # 平衡分布の3次欠損 rho*u*u*u を補う
            feq = feq + VXXY[:, None, None]*(rho*ux*ux*uy) \
                      + VXYY[:, None, None]*(rho*ux*uy*uy)
        src = force_populations(ux, uy, fx, fy, tau)
        if cubic_force:                   # 力項の3次欠損 F*u*u を補う
            src = src + half*(VXXY[:, None, None]*(2*fx*ux*uy + fy*ux*ux)
                              + VXYY[:, None, None]*(fx*uy*uy + 2*fy*ux*uy))
        f = f - (f - feq)/tau + src
        for i in range(9):
            f[i] = np.roll(f[i], (int(CX[i]), int(CY[i])), axis=(0, 1))
        if s % 1000 == 999:
            if np.max(np.abs(ux - prev)) < tol:
                break
            prev = ux.copy()
    return ux[0]
 
print("D2Q9 lock check:  max|cx^3 - cx| =", float(np.abs(CX**3 - CX).max()))
for tag, mom in (("xxy", CX*CX*CY), ("xyy", CX*CY*CY), ("xxx", CX**3), ("xy", CX*CY)):
    print(f"  correction vector v_xxy :  sum v_i c_{tag:3s} = {float(mom @ VXXY): .6f}")
 
NY, TAU = 16, 1.0
NU, A0  = CS2*(TAU - 0.5), 0.004
print(f"\nny={NY}  tau={TAU}  nu={NU:.5f}  k={2*np.pi/NY:.4f}  a0={A0}")
print("\n   V    |u|max   Guo       +force    +feq3     +both")
for V in (0.05, 0.10, 0.20):
    ref = analytic_wave(NY, NU, V, A0)
    errs = []
    for ce, cf in ((0, 0), (0, 1), (1, 0), (1, 1)):
        u = run_cross_flow(NY, TAU, V, A0, bool(ce), bool(cf))
        errs.append(100*np.linalg.norm(u - ref)/np.linalg.norm(ref))
    print(f"{V:5.2f}  {np.abs(ref).max():.5f}  " + "  ".join(f"{e:7.4f}%" for e in errs))

D2Q9 lock check:  max|cx^3 - cx| = 0.0
  correction vector v_xxy :  sum v_i c_xxy =  1.000000
  correction vector v_xxy :  sum v_i c_xyy = -0.000000
  correction vector v_xxy :  sum v_i c_xxx =  0.000000
  correction vector v_xxy :  sum v_i c_xy  = -0.000000
 
ny=16  tau=1.0  nu=0.16667  k=0.3927  a0=0.004
 
   V    |u|max   Guo       +force    +feq3     +both
 0.05  0.12258   0.5993%   0.5915%   0.0181%   0.0143%
 0.10  0.08374   1.6422%   1.6207%   0.0295%   0.0375%
 0.20  0.04826   3.7307%   3.6840%   0.0531%   0.1503%

表の +feq3 と +both の差は、波長を縮めるほど広がります。波形トグルと波長スライダーが、その漏れの棒です。

uniform では緑の棒が0になります。sine で波長を縮めると棒は波数の2乗に沿って伸び、格子の総和は0のままです。

0.053%を0.150%にした補正#

表を一行ずつ読むと順番が見えます。$V=0.20$ での値です。

力の3次モーメントを埋めた効果はせいぜい1.3%です。平衡分布の3次を埋めると誤差が70分の1になります。 そして両方入れると再び悪化します。三つの数字が同じことを語っています。

理由はChapman-Enskog展開の次数にあります。力の2次モーメントは $O(\varepsilon)$ の応力式に直接入ります。 力の3次モーメントは $O(\varepsilon)$ の3次モーメント式に入り、そこから応力へ一段、 応力から運動量へもう一段降りてきます。運動量方程式には $O(\varepsilon^3)$、つまりBurnett次数で到着します。 Navier-Stokesを目標とする等温LBMでは、合わせる対象ではありません。

一方、平衡分布の3次モーメントは $O(\varepsilon)$ の応力に直接入ります。それが 平衡分布の3次欠損が支配的になる理由です。

補正が誤差を作った経路

では力の3次を埋めたことが、なぜ害になるのでしょうか。前節の漏れの式が答えです。 $d = F_x u_y^2 = a_0 V^2 \sin(ky)$ なので、ストリーミングが次の大きさの偽の体積力を作ります。

$V=0.2$、$k=0.3927$、$\tau=1$ を入れると $a_0$ の $0.154\%$ です。 コードで測った二つの解の最大差は $|u|_{\max}$ の $0.156\%$ でした。二つの数字が一致します。 運動量方程式が求めてもいない項を、格子が忠実に配達したわけです。

$V^2 k^2$ 依存性もそのまま見えます。$V$ を $0.05 \to 0.10 \to 0.20$ と上げると、 二つの解の差は $0.010\% \to 0.039\% \to 0.156\%$ でした。毎回4倍です。

3次モーメントが本当に必要になる場所#

だからといって $\Phi_{\alpha\beta\gamma}$ の導出が無意味なわけではありません。必要な場所は確かにあります。

エネルギー方程式を一緒に解く熱・圧縮性LBMでは、3次モーメントが熱流束です。 そこでは力の3次モーメントが $O(\varepsilon^2)$ のエネルギー式に直接入ります。Burnettまで行く必要はありません。 非理想状態方程式で界面を作る多相LBMも似ています。力が大きく、空間的に急です。 $k$ が1に近づけば $V^2k^2$ はもう小さくありません。

$V=0.20$ で標準Guoは3.7307%です。力の3次だけ埋めると3.6840%です。差は0.0467%pで、3.7307で割ると約1.3%です。平衡分布の3次だけだと0.0531%です。3.7307を0.0531で割ると約70です。両方入れると0.1503%です。0.1503を0.0531で割ると約2.83で、表の2.8倍と合います。$V=0.05$ では両方入れた方が0.0143%で、平衡だけ0.0181%より小さいです。$V=0.10$ では0.0375%と0.0295%で順番が逆転します。$V=0.20$ の間隔は0.0972%pです。

漏れの式に $V=0.2$、$k=0.3927$、$\tau=1$ を入れると $1-1/(2\tau)$ は $1/2$ です。前の $1/2$ と掛けると $1/4$ で、$V^2=0.04$ です。$k^2$ は $0.3927$ の2乗で約 $0.1542$ です。$1/4$ 掛ける $0.04$ 掛ける $0.1542$ は $0.001542$、つまり $a_0$ の $0.154\%$ です。コードで測った差は $|u|_{\max}$ の $0.156\%$ でした。$0.002\%$p以内です。$V$ を $0.05$ から $0.10$、$0.20$ と2倍にすると $V^2$ は4倍です。間隔が $0.010\%$ から $0.039\%$、$0.156\%$ と4倍ずつ増えたことと合います。$|u|_{\max}$ は $0.12258$、$0.08374$、$0.04826$ と減ります。$C_1$ の分母は $\nu k^2+V^2/\nu$ で $\nu=0.16667$ なので、$V$ が大きいほど振幅は削れます。

力の半分をどこに置くかは、Guo、Shan–Chen、EDMを並べた比較が同じ $(1-1/2\tau)$ で分けています。目標の方程式が何か、その項が何次で届くかを先に数える必要があります。

格子にない項を半分だけ埋めてはいけない理由

今回の実験の本当の教訓は、次数の話より前にあります。$\Phi_{\alpha\beta\gamma}$ はテンソルです。 D2Q9はそのうち $xxy$、$xyy$ だけを許し、$xxx$、$yyy$ を固定します。 テンソルの半分だけ埋めれば、等方性の壊れた補正を入れたことになります。 打ち消し合っていた誤差が表に出ることがあり、今回がそれでした。

格子が運べない項を見つけたときの選択肢は二つです。 速度集合を増やしてそのモーメントを確保するか(D2Q17、D2Q37)、まったく触らないか。 中間はありません。力項に一行足す前に、その行の残りの成分が格子にあるかを数えるべきです。