ノード1つで誤差が0.21%から76%に跳ねた — べき法則流体LBMの局所緩和時間

べき法則LBMにおけるτの上限は安定化オプションではなく、定常状態に到達するための条件です。

収束しないログに残った inf の一行#

血液をべき法則(power-law)流体として、格子ボルツマン法(LBM)でチャネル流れを回しました。速度場はもっともらしく見えました。残差だけが下がりません。2万ステップ回しても $10^{-3}$ 付近で止まったままでした。

緩和時間 $\tau$ をノードごとに出力すると、チャネル中心線のノード1つに inf がありました。タイプミスでもゼロ除算でもありません。そのノードは本当に無限大の粘度を要求していたのです。

発散はバグではありません。べき法則そのものの性質です。べき指数のスライダーを1.0から0.4へ下げると、中心ノードの $\tau$ が赤い帯へ押し出されます。

power-law index n = 0.60 tau cap = 5.00 nodes  cap tau

Uncapped, the centre node reads tau = infinity — the shear rate vanishes there, so a shear-thinning fluid asks for infinite viscosity. Green band 0.55 < tau < 2 is the usable range; red on the left is the stability edge, red on the right is where bounce-back walls start leaking.

左の速度分布は中心で平らになります。cap tau を外すと、そのノードは軸の外へ消えます。0.21%と76%を分けた場所がそこです。

LBMがひずみ速度をただでくれる場所#

有限体積法で非ニュートン粘度を扱うには、まず速度勾配が要ります。隣接セルを読み、最小二乗なりGreen–Gaussなりで再構成し、それから $\dot\gamma$ が出ます。並列計算ではここに通信が1回増えます。

LBMはその段階を飛ばします。Chapman–Enskog展開(クヌーセン数のべき級数で分布関数を分解する多重スケール手法)を1次追うだけで、非平衡分布の2次モーメントがひずみ速度テンソルと直接結びつきます。

S_{\alpha\beta} = -\frac{1}{2\rho c_s^2 \tau \Delta t} \sum_i c_{i\alpha} c_{i\beta}\, f_i^{neq}

$f_i^{neq} = f_i - f_i^{eq}$ はそのノードの平衡からのずれ、 $c_s^2 = 1/3$ は格子音速の2乗、 $\tau$ は緩和時間です。右辺はすべて1つのノード内の量です。隣を読みません。

同じ2次モーメントは、別の問題では欠損として顔を出すこともあります。ここでは逆に贈り物です。粘性応力はすでに分布関数の中にあり、取り出すだけで済みます。

べき法則はこの $\dot\gamma$ 1つで閉じます。

\dot\gamma = \sqrt{2 S_{\alpha\beta} S_{\alpha\beta}}, \qquad \nu = K \dot\gamma^{\,n-1}

$K$ は稠度係数、 $n$ はべき指数です。 $n<1$ ならせん断減粘(shear-thinning) — 速くかき混ぜるほど薄くなります。血液は $n \approx 0.6$ 付近です。 $n=1$ でニュートン流体に戻ります。

粘度が τ の中にあり τ が粘度の中にある

式を2つ並べると問題が見えます。BGK衝突では粘度を緩和時間が決めます。

\nu = c_s^2 \left( \tau - \frac{1}{2} \right) \Delta t

ところが上の $S_{\alpha\beta}$ の式にも $\tau$ が分母にあります。 $\tau$ を知らないと $\dot\gamma$ が測れず、 $\dot\gamma$ がないと $\nu$ が出ず、 $\nu$ がないと $\tau$ が決まりません。ノードごとに小さな非線形方程式が1つ置かれている形です。 $\nu = c_s^2(\tau - 1/2)\Delta t$ に半ステップがどこで付くかは離散LBEの2次展開に分けてあります。ここではその式をノードごとに解き直す側だけを扱います。

多くの実装はここで前ステップの $\tau$ をそのまま使います。この遅延(lagging)は定常状態では正しいのですが、界面や拍動のようにせん断が速く変わる場所では、1ステップ分ずれた粘度で衝突を回すことになります。

遅延せずに解くこともできます。鍵は $\tau$ に依存しない量を先に分離することです。非平衡2次モーメントの偏差(deviatoric)部分を $\Pi'_{\alpha\beta}$ と書き

Q \equiv \sqrt{2\, \Pi'_{\alpha\beta} \Pi'_{\alpha\beta}}

と置くと、 $Q$ はその瞬間の $f_i$ だけで計算できます。 $\tau$ が入る余地がありません。そして $\dot\gamma = Q / (2\rho c_s^2 \tau)$ です。

不動点1つで閉じる局所方程式#

$m \equiv c_s^2 \tau$ と変数を替えると、3つの式が1行にまとまります。

m = \frac{c_s^2}{2} + B\, m^{1-n}, \qquad B \equiv K \left( \frac{Q}{2\rho} \right)^{n-1}

未知数は $m$ だけ、係数 $B$ は定数です。この形が良い理由はもう1つあります。右辺を $g(m)$ とすると

g'(m) = B(1-n)\, m^{-n}

であり、解では $m = c_s^2/2 + B m^{1-n} > B m^{1-n}$ なので $B m^{-n} < 1$ です。したがって $n<1$ のとき常に $|g'| < 1-n < 1$ が成り立ちます。せん断減粘流体なら単純代入反復が無条件に縮小写像になります。Newton法も安全装置も初期推定も要りません。

D2Q9で $\Pi_{xy}$ に寄与するのは対角リンク4本だけです。 $c_{ix}c_{iy}$ が0の軸方向リンクは、このモーメントに入りません。

同じチャネルに上限を5通り入れてみた#

体積力で駆動する完全発達チャネルです。壁は halfway bounce-back、幅は31ノード。べき法則チャネル流れには解析解があるので、誤差をそのまま測れます。

u(y) = \frac{n}{n+1} \left( \frac{g}{K} \right)^{1/n} \left( h^{\frac{n+1}{n}} - |y|^{\frac{n+1}{n}} \right)

$g$ は単位質量あたりの体積力、 $h$ は半幅です。 $n=1$ を入れると見慣れた放物線に戻ります。

import numpy as np
 
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, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])
OPP = np.array([0, 3, 4, 1, 2, 7, 8, 5, 6])
CS2, G, NY = 1/3, 2e-6, 31
 
def eq(rho, ux):
    cu = CX[:, None] * ux
    return W[:, None] * rho * (1 + cu/CS2 + cu*cu/(2*CS2**2) - ux*ux/(2*CS2))
 
def relax_time(Q, rho, n, K, sweeps=40):
    """m = cs2*tau。隣接ノードを一度も読まない不動点反復。"""
    with np.errstate(divide='ignore', over='ignore', invalid='ignore'):
        B = K * (Q / (2*rho))**(n - 1)
        m = np.ones_like(Q)
        for _ in range(sweeps):
            m = CS2/2 + B * m**(1 - n)
    return m / CS2
 
def exact(n, K):
    h = NY/2.0
    y = np.arange(NY) - (NY - 1)/2.0
    return n/(n+1) * (G/K)**(1/n) * (h**((n+1)/n) - np.abs(y)**((n+1)/n))
 
def consistency_K(n, umax=8e-3):
    h = NY/2.0
    return G / (umax*(n+1)/(n*h**((n+1)/n)))**n
 
def channel(n, K, cap, steps=20000, tol=1e-12):
    ue = exact(n, K)
    f = eq(np.ones(NY), ue)          # 解析解から出発し、定常状態への補正だけを緩和させる
    prev, tau = ue.copy(), np.full(NY, 0.8)
    for it in range(steps):
        rho = f.sum(0)
        ux = (CX[:, None]*f).sum(0)/rho + 0.5*G
        fe = eq(rho, ux)
        Pxy = (CX[:, None]*CY[:, None]*(f - fe)).sum(0)
        raw = relax_time(2*np.abs(Pxy), rho, n, K)
        tau = np.clip(np.nan_to_num(raw, nan=cap, posinf=cap), 0.5005, cap)
        cu = CX[:, None]*ux
        src = W[:, None]*(1 - 1/(2*tau))*rho*G*((CX[:, None]-ux)/CS2 + cu*CX[:, None]/CS2**2)
        post = f - (f - fe)/tau + src
        for i in range(9):
            f[i] = np.roll(post[i], CY[i])
        for i in range(9):                                  # halfway bounce-back
            if CY[i] == 1:  f[i, 0]  = post[OPP[i], 0]
            if CY[i] == -1: f[i, -1] = post[OPP[i], -1]
        if it % 200 == 199:
            res = np.abs(ux - prev).max()/ux.max()
            if res < tol: break
            prev = ux.copy()
    rho = f.sum(0)
    ux = (CX[:, None]*f).sum(0)/rho + 0.5*G
    err = np.linalg.norm(ux - ue)/np.linalg.norm(ue)
    return err, res, it + 1, tau, raw
 
K1 = consistency_K(1.0)
e, r, it, tau, _ = channel(1.0, K1, 1e6)
print(f"Newtonian n=1.0  L2={e:.4%}  tau={tau.min():.4f}..{tau.max():.4f}  steps={it}")
 
n, K = 0.6, consistency_K(0.6)
_, _, _, _, raw = channel(n, K, 4.0, steps=200)
print(f"n=0.6  tau at wall {raw[0]:.4f}   tau at centre {raw[NY//2]:.4e}")
print(f"{'cap':>9} {'L2 err':>9} {'residual':>10} {'steps':>7} {'capped':>8}")
for cap in [0.8, 1.5, 4.0, 40.0, 1e6]:
    e, r, it, tau, _ = channel(n, K, cap)
    print(f"{cap:9.1f} {e:9.4%} {r:10.2e} {it:7d} {int((tau >= cap - 1e-9).sum()):5d}/{NY}")

Newtonian n=1.0  L2=0.1362%  tau=0.5901..0.5901  steps=20000
n=0.6  tau at wall 0.5690   tau at centre inf
      cap    L2 err   residual   steps   capped
      0.8   0.2112%   1.08e-06   20000     3/31
      1.5   0.2121%   1.14e-06   20000     1/31
      4.0   0.2134%   1.07e-06   20000     1/31
     40.0   0.4850%   1.02e-08   20000     1/31
1000000.0  75.9752%   7.69e-03   20000     1/31

Pi_xy を上げると緑の曲線が下がり、不動点の階段が対角線に届く段数が減ります。

Pi_xy = 4.00e-4 n = 0.60 K = 2.5e-3  run

The staircase is the whole solver: 0 steps of m → cs2/2 + B m^(1-n) land on tau = 0.5894 (gamma = 2.04e-3, nu = 2.98e-2). Freezing tau at the previous step instead would answer 0.6011. No neighbour node is read at any point.

$n$ を1より上にすると傾きが負になり、階段は螺旋を描きます。 $n<2$ ならそれでも縮小します。ロゼットの軸方向リンクが灰色なのは、 $c_{ix}c_{iy} = 0$ で $\Pi_{xy}$ に乗らないからです。

0.21%と76%を分けたのはノード1つだった#

ニュートン流体の基準線が0.1362%です。ソルバ自体は正しく動いています。同じ格子に $n=0.6$ を入れ、上限だけを変えたところ、誤差が0.2112%から75.9752%まで動きました。

capped の列を見てください。上限が1.5以上のとき、引っかかるノードは31個中ちょうど1個です。チャネル中心線のノードです。その1ノードの $\tau$ を $10^6$ にするか4.0にするかの間に、誤差75.76ポイントが挟まっています。

理由は衝突の式にあります。 $f_i \leftarrow f_i - (f_i - f_i^{eq})/\tau$ で $\tau \to \infty$ にすると、括弧の前の係数が0になります。そのノードは平衡に戻りません。初期分布が持っていた誤差がそのまま残り、周囲のノードがストリーミングで持ち去ります。残差が $7.69 \times 10^{-3}$ で止まったのはその結果です。

なぜ中心線で $\dot\gamma$ がちょうど0なのか。対称性です。チャネル中心で $\partial u / \partial y$ が符号を変えるので、その位置で0を通ります。離散的にも $\Pi_{xy}$ が相殺してちょうど0になります。すると $\nu = K \dot\gamma^{\,n-1}$ で $n-1 = -0.4$ ですから $\nu \to \infty$ です。連続体でも同じです。べき法則はせん断速度0で定義されません。

壁の $\tau = 0.5690$ を粘性に戻すと $\nu = c_s^2(\tau - 1/2) = (1/3)(0.0690) = 0.0230$ です。ニュートン側の $\tau = 0.5901$ は $(1/3)(0.0901) = 0.030033$ です。同じ $G = 2\times 10^{-6}$ 、幅31ノードでも、壁の粘性はすでに違います。コードの下限クリップは0.5005です。そこでは $\nu = (1/3)(0.0005) = 1.67\times 10^{-4}$ まで下がり、壁の0.0230より二桁以上小さくなります。 $n = 0.6$ なら $1-n = 0.4$ なので、不動点の傾きは $|g'| < 0.4$ の中に収まります。縮小はします。ただし中心で $Q = 0$ になると、指数 $n-1 = -0.4$ が $B$ を発散させます。

上限0.8と4.0のL2は0.2112%と0.2134%です。差は0.0022パーセントポイントです。1.5は0.2121%で、0.8より0.0009パーセントポイント大きいだけです。引っかかったノードは3個から1個に減ります。上限を $10^6$ まで開くと誤差75.9752%、残差 $7.69\times 10^{-3}$ です。上限4.0の残差 $1.07\times 10^{-6}$ より三桁以上大きい値です。ステップ数はどちらも20000です。上限を閉じれば $1.07\times 10^{-6}$ まで下がり、開くと $7.69\times 10^{-3}$ で止まります。冒頭の $10^{-3}$ はその開いた側です。引っかかったノードは相変わらず31個中1個です。

上限値をどこに置いてもほぼ同じ答えになる理由

ここに実務上重要な事実があります。0.8、1.5、4.0の誤差が0.2112%、0.2121%、0.2134%です。3桁まで一致しています。

そのノードはせん断応力をほとんど運んでいないからです。完全発達チャネルで運動量収支が決める応力は $\sigma(y) = \rho g y$ です。中心では $y=0$ なので $\sigma=0$ です。粘度をいくら間違えても、掛かるせん断速度が0なので応力への寄与がありません。

つまり発散するのは $\nu$ であって $\nu\dot\gamma$ ではありません。流れが実際に感じるのは後者です。だから上限は精度を損なわずに入れられます。上限を40まで上げると誤差が0.4850%と倍になりますが、これは $\tau$ が大きくそのノードの緩和が遅いためです。精度ではなく収束速度の問題です。

逆に上限を下げすぎると別の代償があります。0.8では3ノードが引っかかります。中心近くの2ノードが実際より薄くなり、その誤差は応力に乗ります。表ではまだ0.2112%と小さいのですが、 $n$ がさらに小さくなると引っかかるノード数が急に増えます。最初のシミュレーションで nodes を81に上げ、 $n$ を0.4まで下げると、赤い帯に入るノードが目に見えて増えます。

血液を解くときに話が変わるところ

べき法則がせん断速度0で発散することは、レオロジーでも知られた欠点です。実際の血液はせん断速度が0に向かっても粘度が有限値(ゼロせん断粘度)に留まります。Carreau–Yasudaモデルがその部分を直します。

\nu = \nu_\infty + (\nu_0 - \nu_\infty)\left[ 1 + (\lambda \dot\gamma)^a \right]^{\frac{n-1}{a}}

$\nu_0$ と $\nu_\infty$ が両端を塞ぎ、 $\lambda$ はせん断減粘が始まるせん断速度の逆数です。このモデルなら $\dot\gamma \to 0$ で $\nu \to \nu_0$ と有限です。上限を手で入れる必要がありません。 $\nu_0$ がそのまま上限です。

ではなぜべき法則を使うのか。管中央を除く大半の領域で両モデルの差が小さく、べき法則には検証できる解析解があるからです。さらに上の不動点が1行で閉じるのはべき法則のときだけです。Carreau–Yasudaでも同じ反復になりますが、縮小性の証明はそこまできれいになりません。

LBM側で $\tau$ が大きくなるときにもう1つ注意すべきことがあります。bounce-back壁の実効位置が $\tau$ とともに動く点です。せん断減粘流体は壁近くで $\tau$ が小さいので、この誤差はむしろ小さくなります。代わりに $\tau$ が $0.5$ に近づくと安定性の余裕が減ります。上の計算では壁ノードが0.5690でした。余裕はありますが、十分とまでは言えません。 $\tau - 0.5 = 0.069$ が、下限までの間隔です。

発散することと使えないことは別である

inf はゼロ除算のバグではありませんでした。せん断速度0で $n-1 = -0.4$ なら $\nu \to \infty$ が正しい答えです。計算が止まったのは、BGKが $\tau^{-1}$ を係数に使うからです。粘度が無限大だという文と、そのノードを永久に緩和しないという文が、同じ式になります。

上限は後者を止める離散化です。表では1.5と4.0のL2が0.2121%と0.2134%で、ほとんど同じです。 $n=1$ なら同じチャネルの式が放物線に戻り、その基準はポアズイユの解析解の側にあります。ログに引っかかったノードが1を大きく超えたら、直すのは上限ではなく $K$ か格子です。