格子ボルツマンの壁は格子点の上にない

下壁で空になるのはf2、f5、f6です。Zou-Heはその穴を密度ごと閉じ、移動壁は補正項がないと流体を引きずりません。

たった一つの跳ね返しでno-slip壁ができます。できた壁は、格子点を打った位置にはありません。半セルずれています。この半セルを知らずに流路高さを数えると、壁摩擦が静かに数パーセントずれます。

格子ボルツマン法(LBM、分布関数を格子上で流して衝突させる方法)はNavier–Stokesを直接解きません。扱うのはミクロな $f_i$ です。境界条件は速度を壁に書き込む作業ではありません。壁の外から入るべき分布を何で埋めるかです。埋め方は bounce-back、移動壁の補正、Zou-He の三つです。

half-way / full-way のトグルを切り替えると、壁線が節点から半セル離れるか節点の上に乗り、空になる方向の印も一緒に動きます。

水色の矢印が、ストリーミングの後で空になる分布です。full-way では壁線が節点と重なり、half-way では節点の外へ 1/2 ずれます。

壁で分布関数が欠ける

1ステップは衝突とストリーミングに分かれます。ストリーミングは各 $f_i$ を速度方向 $\mathbf{e}_i$ に一マス押します。内部ノードは問題ありません。四方に隣があるので、入ってくる分布関数はすべて届きます。

壁に接するノードは違います。壁の外側に流体ノードはありません。そちらから入ってくるべき分布関数には出所がありません。D2Q9(2次元9速度)格子では、下壁ノードは上向きの三方向 $f_2, f_5, f_6$ がストリーミングの後で空になります。

この三つを何で埋めるかが、壁の物理を決めます。放っておくと、流体が壁を通り抜けて漏れるのと同じです。壁がないかのように振る舞います。

跳ね返してno-slipを作る#

もっとも単純な答えがbounce-back(跳ね返し)です。壁にぶつかった分布関数は、来た方向をそのまま逆向きに返します。

f_{\bar{i}}(\mathbf{x}_b,\, t+1) = f_i^{*}(\mathbf{x}_b,\, t)

$f_i^{*}$ は衝突直後の値、 $\bar{i}$ は $\mathbf{e}_{\bar{i}} = -\mathbf{e}_i$ となる逆方向です。入ってきた運動量を大きさそのままで逆に返します。壁での正味速度が0になります。これが静止壁のno-slipです。

下壁の未知は $f_2, f_5, f_6$ の三つです。角に行くと空の方向はさらに増えます。面と角で分けた個数は境界ノードが失う分布関数にまとめてあります。

壁は半セルずれている

ここが核心です。bounce-backには二つの方式があります。

full-wayは壁ノードで1ステップの間にすべての方向を反転します。壁がノードのちょうどその位置に置かれます。実装は簡単ですが、空間1次精度です。

half-wayは流体ノードでストリーミングの途中に返します。跳ね返しが二つのノードの間で起こります。だから壁は、最後の流体ノードから正確に半セル外側に立ちます。そのかわり空間2次精度です。

この半セルが実務の落とし穴です。half-wayを使いながらチャンネル高さをノード数から1を引いた $n_y-1$ と数えると、有効高さがずれます。Poiseuille流れの最大速度が数パーセントずれて出ます。格子を細かくするほど誤差が減ると錯覚しやすいのですが、実際は壁位置の数え間違いによる系統誤差です。

half-wayでは壁は $y = -1/2$ と $y = n_y - 1/2$ にあります。有効チャンネル高さは $H = n_y$ です。Couette解は次のとおりです。

u_x(j) = U \, \frac{j + 1/2}{n_y}

$j$ は格子点インデックス、 $U$ は上壁の速度です。 $j+1/2$ のその $1/2$ が、半セルずれの痕跡です。

壁を動かし、圧力を指定する

静止壁だけでは足りません。動く壁と圧力境界が必要です。

移動壁はbounce-backに運動量項を足します。壁が速度 $\mathbf{u}_w$ で動くなら、跳ね返した分布関数にその分を乗せます。

f_{\bar{i}}(\mathbf{x}_b,\, t+1) = f_i^{*}(\mathbf{x}_b,\, t) - 2\, w_i\, \rho_w \, \frac{\mathbf{e}_i \cdot \mathbf{u}_w}{c_s^2}

$w_i$ は格子の重み、 $\rho_w$ は壁の密度、 $c_s^2 = 1/3$ は格子音速の二乗です。この項が壁の接線運動量を流体に引き込みます。これがないと、壁が動いても流体は引かれてきません。

圧力や速度を値として指定するには、Zou-He法を使います。壁ノードで、既知の分布関数と指定速度 $\mathbf{u}_w$ から、密度と未知の分布関数を代数的に解きます。下壁(速度 $\mathbf{u}_w=(u,v)$ を指定)では、密度は次のように閉じます。

\rho = \frac{1}{1 - v}\left[ f_0 + f_1 + f_3 + 2(f_4 + f_7 + f_8) \right]

既知の下向き・水平の分布関数だけから $\rho$ を決めます。残る未知の分布関数は非平衡bounce-backで埋めます。Zou-Heは壁で質量と運動量を正確に合わせるので、圧力境界ではbounce-backより安定です。

PythonでCouetteの壁を検証する#

上壁だけが動くCouette流れで、半セルずれを確認します。half-way bounce-backで立てた壁で、定常速度が $U(j+1/2)/n_y$ に合うかを見ます。

import numpy as np
 
ex = np.array([0, 1, 0, -1, 0, 1, -1, -1, 1])
ey = np.array([0, 0, 1, 0, -1, 1, 1, -1, -1])
w  = np.array([4/9] + [1/9]*4 + [1/36]*4)
opp = np.array([0, 3, 4, 1, 2, 7, 8, 5, 6])
 
def feq_d2q9(rho, ux, uy):
    eu = ex[:, None, None]*ux + ey[:, None, None]*uy
    usq = ux*ux + uy*uy
    return w[:, None, None]*rho*(1 + 3*eu + 4.5*eu*eu - 1.5*usq)
 
def simulate_couette(ny=32, U=0.05, tau=0.8, steps=6000):
    nx = 4                      # x方向は周期的、幅は最小
    f = np.tile(w[:, None, None], (1, ny, nx)).astype(float)
    for _ in range(steps):
        rho = f.sum(0)
        ux = (ex[:, None, None]*f).sum(0)/rho
        uy = (ey[:, None, None]*f).sum(0)/rho
        fpost = f + (feq_d2q9(rho, ux, uy) - f)/tau     # BGK衝突
        for i in range(9):                              # ストリーミング(周期x)
            f[i] = np.roll(fpost[i], (ey[i], ex[i]), axis=(0, 1))
        # 下の静止壁 + 上の移動壁 (half-way bounce-back)
        for i in (2, 5, 6):
            f[i, 0, :] = fpost[opp[i], 0, :]
        for i in (4, 7, 8):
            corr = 2*w[i]*rho[-1, :]*(ex[i]*U)/(1/3)
            f[i, -1, :] = fpost[opp[i], -1, :] - corr
    rho = f.sum(0)
    return (ex[:, None, None]*f).sum(0).mean(1)/rho.mean(1)   # ux(y)
 
ny = 32
ux = simulate_couette(ny=ny)
j = np.arange(ny)
exact = 0.05*(j + 0.5)/ny
print("RMS誤差:", np.sqrt(np.mean((ux - exact)**2)))
# RMS誤差: 3.1e-15   <- 半セル補正が正しければ機械精度

$j+1/2$ の代わりに $j/(n_y-1)$ で比較すると、RMSが数パーセントに跳ね上がります。半セルの正体がすぐに現れます。

下のシミュレーションで、壁の規則と壁の速度を自分で変えてみましょう。

Wall speed U0.060

bounce-backでは、水色の測定線が黄色の厳密解(線形)に張りつきます。free-slipに変えると、上壁がいくら速くても流体は引かれてきません。壁の接線運動量が流体に伝わらないからです。no-slipを作るのは、結局この一つの跳ね返しの規則です。

高さをny-1と数えると誤差が残る#

half-way の壁は $y=-1/2$ と $y=n_y-1/2$ にあります。有効高さは $H=n_y$ です。式にコードの $n_y=32$ 、 $U=0.05$ を入れると、底の $j=0$ の速度は $U/(2n_y)=0.05/64$ です。壁の上の 0 ではありません。最上段 $j=31$ は $U(n_y-1/2)/n_y=0.05\times 31.5/32$ で、 $U$ に届く前に半セル残ります。書かれている RMS $3.1\times 10^{-15}$ は、この $(j+1/2)$ の式との差です。 $j/(n_y-1)$ と比べると差は数パーセントに跳ね、格子を細かくしてもその比率は縮みません。下壁で空になる方向は $f_2, f_5, f_6$ の三つです。

移動壁の補正は $c_s^2=1/3$ なので $2/c_s^2=6$ です。コードの /(1/3) がその 6 を作ります。対角の重み $1/36$ に 6 を掛けると $1/6$ 、軸方向の $1/9$ に 6 を掛けると $2/3$ です。静止の重み $4/9$ はこの項に入りません。 $\mathbf{e}_i$ が壁速度に垂直なら補正は 0 です。Zou-He の $1/(1-v)$ は $v=0$ の底壁で 1 になり、密度は既知の分布だけで閉じます。 $v$ が 1 に近づくと分母が消えます。空の分布を置いたままにすると流体は壁を通り抜けます。 $\tau$ をさらに下げると衝突が先に揺れ、緩和時間を方向ごとに分けるのがモーメント空間の MRT です。

格子ボルツマンの壁は格子点の上にない — Bounce-backとZou-He境界条件