τ=0.6で流速が6分の1になった — Guo forcingで半ステップ速度F/2を忘れると

Guoソース項の(1−1/2τ)係数は平衡速度のF/2補正と対です。片方だけだと、毎ステップFではなく(1−1/2τ)Fしか運動量が入りません。

同じ力をかけたのに、なぜ流速が6分の1なのか#

周期境界の箱に正弦波状の体積力をかけました。解析解は教科書の一行です。ところが格子ボルツマン法(LBM、分布関数を衝突・並進させて流体を解く方法)のコードが出した振幅は、解析解の0.16倍でした。粘性を変えると比率も変わりました。τ=1ではちょうど半分でした。

原因は一行です。巨視的速度を計算するときに $\mathbf{F}/2$ を足していませんでした。この記事では、その一行がなぜ正確に $(1-1/2\tau)$ 倍の誤差を生むのかを示します。Kolmogorov流(正弦体積力が駆動する周期せん断流)で実際に測り、自分のコードでこのバグを30秒で見つける方法で締めくくります。

1ステップの中で運動量は3回変わる#

LBMの1ステップは衝突と並進です。並進は運動量を運ぶだけです。運動量が生まれるのは衝突だけです。したがって、一様な流体の1点の帳簿を見れば十分です。

\mathbf{m}' = \mathbf{m} - \frac{1}{\tau}\left(\mathbf{m} - \rho\mathbf{u}\right) + \left(1-\frac{1}{2\tau}\right)\mathbf{F}

$\mathbf{m}=\sum_i f_i\mathbf{c}_i$ は衝突前の運動量、 $\mathbf{m}'$ は衝突後、 $\rho\mathbf{u}$ は平衡分布に入れた運動量、 $\tau$ は緩和時間、 $\mathbf{F}$ は格子単位の体積力です。

第2項は緩和項です。平衡速度が衝突前の速度と等しければゼロになります。第3項はGuoソース項の1次モーメントです。 $\rho\mathbf{u}$ に何を入れるかがすべてを決めます。

\rho\mathbf{u} = \mathbf{m} + \tfrac{1}{2}\mathbf{F} \;\Rightarrow\; \mathbf{m}' = \mathbf{m} + \frac{\mathbf{F}}{2\tau} + \mathbf{F} - \frac{\mathbf{F}}{2\tau} = \mathbf{m} + \mathbf{F}

\rho\mathbf{u} = \mathbf{m} \;\Rightarrow\; \mathbf{m}' = \mathbf{m} + \left(1-\frac{1}{2\tau}\right)\mathbf{F}

1行目では、緩和項が $\mathbf{F}/(2\tau)$ を補い、ちょうど $\mathbf{F}$ が入ります。2行目ではその補いが消えます。ソース項の $(1-1/2\tau)$ 係数だけが残ります。

下のシミュレーションで実際に操作してみましょう。

relaxation time τ = 0.60 ν = (τ−½)/3 = 0.0333

step 0 · ρu/(F·t): with F/2 0.0000 · without 0.0000 · theory 1−1/(2τ) = 0.1667

τスライダーを0.51側に動かしてください。赤い線の傾きがほぼゼロになり、右側の棒から濃い緩和の部分がまるごと抜け落ちるのが見えます。

半ステップ速度:衝突の前と後のちょうど中間

なぜ $\mathbf{F}/2$ なのでしょうか。衝突前の運動量は $\mathbf{m}$ 、力をすべて受けた後は $\mathbf{m}+\mathbf{F}$ です。物理的な速度はその時間区間の平均、つまり中間にあります。台形則(区間両端の値の平均で積分する2次精度の公式)と同じ理由です。

Guoら(2002)はソース項を次のように置きました。

S_i = w_i\left[\frac{\mathbf{c}_i-\mathbf{u}}{c_s^2} + \frac{(\mathbf{c}_i\cdot\mathbf{u})\,\mathbf{c}_i}{c_s^4}\right]\cdot\mathbf{F}

$w_i$ は格子重み、 $\mathbf{c}_i$ は離散速度、 $c_s^2=1/3$ は格子音速の2乗です。 $\sum_i S_i\mathbf{c}_i = \mathbf{F}$ 、 $\sum_i S_i = 0$ を満たします。

$(1-1/2\tau)$ 係数はChapman–Enskog展開(分布関数をKnudsen数で展開して巨視方程式を導く方法)から出てきます。離散時間のTaylor展開が、ソース項に $-\tfrac{1}{2}S_i$ 型の残りかすを残すためです。この残りかすは半ステップ粘性補正と同じ根を持ちます。そしてこの係数を打ち消すには、平衡速度に $\mathbf{F}/2$ を入れて1次モーメントを合わせる必要があります。両者は最初から一対として導かれます。

コード:Kolmogorov流と3通りの速度定義#

D2Q9格子を $1\times N$ の周期柱にします。力は $F_x = F_0\sin(ky)$ 、 $k=2\pi/N$ です。定常解析解は次のとおりです。

u_x(y) = \frac{F_0}{\nu k^2}\sin(ky),\qquad \nu = \frac{\tau - 1/2}{3}

$\nu$ は格子単位の動粘性係数です。壁がないので、バウンスバック境界の誤差が混ざりません。力の扱いだけを測れます。

3つのケースを回します。guo は平衡と出力の両方に $\mathbf{F}/2$ を入れます。nohalf はGuoソース項はそのままで、速度から $\mathbf{F}/2$ を外します。rawout は計算は正しく行い、出力時だけ $\sum f_i\mathbf{c}_i/\rho$ を使います。振幅は $\sin(ky)$ のフーリエ係数で測り、解析解で割ります。

import math
 
# D2Q9 lattice: velocities, weights, sound speed squared
C = [(0,0),(1,0),(0,1),(-1,0),(0,-1),(1,1),(-1,1),(-1,-1),(1,-1)]
W = [4/9,1/9,1/9,1/9,1/9,1/36,1/36,1/36,1/36]
CS2 = 1/3
 
def f_equil(rho, ux):
    # uy = 0 everywhere: the flow is u_x(y) only
    out = []
    for (cx, cy), wi in zip(C, W):
        cu = cx*ux
        out.append(wi*rho*(1 + cu/CS2 + 0.5*(cu/CS2)**2 - 0.5*ux*ux/CS2))
    return out
 
def guo_term(ux, Fx, tau):
    out = []
    for (cx, cy), wi in zip(C, W):
        cu = cx*ux
        out.append((1 - 0.5/tau)*wi*((cx - ux)/CS2 + cu*cx/CS2**2)*Fx)
    return out
 
def kolmogorov_amp(tau, mode, N=16, F0=1e-6):
    """mode 'guo'    : u = (sum f c + F/2)/rho in equilibrium and output
       mode 'nohalf' : Guo source, but u = sum f c / rho everywhere
       mode 'rawout' : correct update, but report sum f c / rho"""
    k = 2*math.pi/N
    Fx = [F0*math.sin(k*(j + 0.5)) for j in range(N)]
    nu = (tau - 0.5)/3
    steps = int(12/(nu*k*k))                    # ~12 viscous decay times
    f = [f_equil(1.0, 0.0) for _ in range(N)]
    for _ in range(steps):
        post = []
        for j in range(N):
            fj = f[j]
            rho = sum(fj)
            mx = sum(fi*cx for fi, (cx, _) in zip(fj, C))
            ux = mx/rho if mode == 'nohalf' else (mx + 0.5*Fx[j])/rho
            feq = f_equil(rho, ux)
            src = guo_term(ux, Fx[j], tau)
            post.append([fi - (fi - fe)/tau + s for fi, fe, s in zip(fj, feq, src)])
        # streaming: only the y component moves anything (one cell in x, periodic)
        f = [[post[(j - C[i][1]) % N][i] for i in range(9)] for j in range(N)]
    amp = 0.0
    for j in range(N):
        rho = sum(f[j]); mx = sum(fi*cx for fi, (cx, _) in zip(f[j], C))
        u = mx/rho if mode in ('nohalf', 'rawout') else (mx + 0.5*Fx[j])/rho
        amp += 2*u*math.sin(k*(j + 0.5))/N
    return amp/(F0/(nu*k*k)), steps
 
print(" tau  steps   guo     nohalf  rawout   nohalf/guo  1-1/(2tau)")
for tau in [0.55, 0.6, 0.8, 1.0, 1.5]:
    r1, n = kolmogorov_amp(tau, 'guo')
    r2, _ = kolmogorov_amp(tau, 'nohalf')
    r3, _ = kolmogorov_amp(tau, 'rawout')
    print(f"{tau:4.2f} {n:6d}  {r1:.4f}  {r2:.4f}  {r3:.4f}   {r2/r1:.4f}      {1-0.5/tau:.4f}")

 tau  steps   guo     nohalf  rawout   nohalf/guo  1-1/(2tau)
0.55   4668  0.9746  0.0885  0.9734   0.0908      0.0909
0.60   2334  0.9754  0.1621  0.9729   0.1662      0.1667
0.80    778  0.9836  0.3660  0.9759   0.3721      0.3750
1.00    466  1.0001  0.4936  0.9872   0.4936      0.5000
1.50    233  1.0772  0.7010  1.0515   0.6508      0.6667

上のコードをブラウザで直接動かしてみましょう。

τ = 0.60 N = 16 steps/frame = 20

step 0/0 · amplitude/exact: guo 0.0000 · nohalf 0.0000 · rawout 0.0000 · nohalf/guo 0.0000 (1−1/2τ = 0.1667)

τを0.6にすると、赤の nohalf の振幅は0.16付近で止まり、赤い破線 $1-1/2\tau$ に張り付きます。Nを32以上に上げると緑の guo は1に近づきますが、nohalf/guo の比はそのままです。

0.166と0.973が語ること#

nohalf/guo 列は $1-1/2\tau$ と重なります。τ=0.55で0.0908対0.0909、τ=0.6で0.1662対0.1667です。運動量帳簿の結論が、実際の流れでもそのまま現れます。定常状態では、粘性拡散が毎ステップ入ってくる運動量とつり合います。入る量が $(1-1/2\tau)$ 倍なら、振幅もその分になります。

τが大きくなるほど、比は理論値から少しずれます。τ=1.5で0.6508対0.6667です。 $N=16$ だと $k\approx0.39$ で、1波長あたりの格子点が少なくなります。この領域では $k^2$ 次の格子打ち切り誤差が混ざります。guo 自身が1.0772と1から離れているのも同じ理由です。

rawout は guo とほぼ同じです。出力速度が $F/(2\rho)$ だけ遅れるだけです。相対誤差は $\nu k^2/2$ 程度で、τ=0.55では約0.1%です。デバッグ時に混同しやすい点です。出力側の $\mathbf{F}/2$ 欠落は小さく、平衡側の欠落は大きいのです。

最も危険なのは高Reynolds数の計算です。粘性を下げようとτを0.5に近づけると、実効的な力もゼロに向かいます。流速が出ないからと力を大きくすると、隠れた係数がすべての結果に残ります。

疑似コードのどこに F/2 が入るのか#

よくあるLBMループの順序は、初期化 → 巨視量の計算 → 平衡分布 → 力 → 衝突 → 並進です。バグはこの順序から生まれます。巨視量の段階で $\mathbf{u}$ を求めるとき、力はまだ計算されていません。そのため $\mathbf{u}=\sum f_i\mathbf{c}_i/\rho$ で済ませてしまいがちです。

解決策は順序を少し入れ替えることです。

まず $\rho = \sum_i f_i$ を求めます。
力 $\mathbf{F}$ を計算します。重力のように $\rho$ だけで済むなら、ここで終わりです。
$\rho\mathbf{u} = \sum_i f_i\mathbf{c}_i + \mathbf{F}/2$ で速度を求めます。
同じ $\mathbf{u}$ で $f^{eq}$ と $S_i$ を作ります。出力にもこの $\mathbf{u}$ を使います。

力が速度に依存する場合(例:抗力)、3番目は陰的な式になります。線形抗力なら代数的に解けます。それ以外は固定点反復を1〜2回行えば十分です。

$(1-1/2\tau)$ を消して「直す」のも答えではありません。1次モーメントは合いますが、2次モーメントにτ依存の応力が残ります。この問題は力の項の3次モーメントの記事とShan–Chen速度シフトforcingの記事で扱いました。

自分のLBMループで F/2 を探す1ステップテスト#

周期箱を作り、 $\rho=1$ 、 $\mathbf{u}=0$ で満たします。一様な力 $\mathbf{F}$ をかけて、ちょうど1ステップだけ回します。そして $\sum f_i\mathbf{c}_i$ の変化量を $\mathbf{F}$ で割ります。

結果が1なら対がそろっています。 $1-1/2\tau$ なら、平衡速度から $\mathbf{F}/2$ が抜けています。τを2つ選んで回せば確実です。比がτに合わせて動くならバグです。検証ケースを回す前に、まずこの1ステップを確かめましょう。