6次モーメントで15が9に落ちた — D2Q9を決めた三つのエルミート根

格子速度の集合は設計の好みではなくガウス・エルミート積分点です。ノードQ個は2Q−1次までしか正確になりません。

1864年、流体とは無関係に整理された多項式#

シャルル・エルミートは1864年、ガウス関数を重みとする直交多項式を整理しました。流体力学とはまったく関係のない純粋解析でした。その85年後、ハロルド・グラッドが同じ多項式でボルツマン方程式の分布関数を展開します。さらに半世紀後、その展開の積分点がLBMコードの速度配列になりました。

この記事では、どのLBMコードの冒頭にも書かれている cs2 = 1.0/3.0 という一行を最後まで辿ります。行き着く先は3次エルミート多項式の根 $\pm\sqrt{3}$ ただ一つです。重み 4/9、1/9、1/36 も同じ場所から出てきます。そして同じ導出が「D2Q9にはできないこと」の一覧も一緒に差し出します。

速度空間の積分には最初から重み関数が入っている

流体計算が必要とするのは分布関数 $f$ そのものではありません。モーメントです。

\rho = \int f \, d\mathbf{v}, \qquad \rho \mathbf{u} = \int f \, \mathbf{v} \, d\mathbf{v}, \qquad \mathbf{\Pi} = \int f \, \mathbf{v}\mathbf{v} \, d\mathbf{v}

$\rho$ は密度、 $\mathbf{u}$ は流速、 $\mathbf{\Pi}$ は運動量流束テンソルです。平衡の近くでは $f$ はマクスウェル分布に多項式を掛けた形になります。マクスウェル分布はガウス関数です。つまり上の積分はすべて「ガウス型の重み関数 × 多項式」の積分です。

この形なら展開の仕方は決まります。重み関数 $w(\mathbf{v}) = (2\pi)^{-D/2} e^{-v^2/2}$ に対して直交する多項式が、まさにエルミート多項式だからです。

f(\mathbf{v}) = w(\mathbf{v}) \sum_{n=0}^{N} \frac{1}{n!} \, \mathbf{a}^{(n)} : \mathbf{H}^{(n)}(\mathbf{v})

$\mathbf{H}^{(n)}$ は $n$ 階のエルミートテンソル、 $\mathbf{a}^{(n)} = \int f \mathbf{H}^{(n)} d\mathbf{v}$ はその係数です。積分は「ガウス重み付き多項式の積分」になり、そうした積分を有限個の点で厳密に終わらせる道具がガウス・エルミート求積法です。Qを1、3、5と動かすと黄色の破線 $2Q-1$ も一緒に動き、その右側だけが赤のまま残ります。

nodes Q = 3 degree n = 6 
n = 6: quadrature 0.0000 vs exact 15.0 → error 6.00e+0  ·  Q = 3 is exact through degree 5 (2Q-1 = 5)

緑の棒はその次数を正確に当てたもの、赤の棒は外したものです。Q=3なら5次まで緑で、6次で色が変わります。その境目が計算全体の主題です。

ノードQ個が2Q−1次で止まる理由#

ガウス求積法の規則は短いものです。 $Q$ 次直交多項式の根にノードを置けば、 $2Q-1$ 次以下の多項式を厳密に積分できます。根 $Q$ 個と重み $Q$ 個、合わせて $2Q$ 個の自由度を使い切った結果です。

$2Q$ ではなく $2Q-1$ なのは、次数を $0$ から数えるからです。 $0, 1, \dots, 2Q-1$ でちょうど $2Q$ 個になります。実際に測ってみるとはっきりします。

from math import factorial
 
def he(n, v):
    """確率論者のエルミート多項式 He_n(v): He_0=1, He_1=v, He_{k+1}=v*He_k - k*He_{k-1}"""
    a, b = 1.0, v
    if n == 0:
        return a
    for k in range(1, n):
        a, b = b, v * b - k * a
    return b
 
def gauss_hermite(Q):
    """ノード = He_Q の根、重み = Q! / (Q^2 He_{Q-1}(v)^2)"""
    nodes, lo, step = [], -8.0, 1e-3
    x0, f0 = lo, he(Q, lo)
    while x0 < 8.0:
        x1 = x0 + step
        f1 = he(Q, x1)
        if f0 == 0.0 or f0 * f1 < 0.0:
            a, b = x0, x1
            for _ in range(80):                 # 二分法
                m = 0.5 * (a + b)
                if he(Q, a) * he(Q, m) <= 0.0:
                    b = m
                else:
                    a = m
            nodes.append(0.5 * (a + b))
        x0, f0 = x1, f1
    wts = [factorial(Q) / (Q * Q * he(Q - 1, v) ** 2) for v in nodes]
    return nodes, wts
 
def exact_moment(n):
    """int v^n N(0,1) dv = 0 (奇数), (n-1)!! (偶数)"""
    if n % 2:
        return 0.0
    m = 1.0
    for k in range(1, n, 2):
        m *= k
    return m
 
def moment_errors(Q, nmax):
    v, w = gauss_hermite(Q)
    return [abs(sum(wi * vi ** n for vi, wi in zip(v, w)) - exact_moment(n))
            for n in range(nmax + 1)]
 
for Q in (3, 5):
    v, w = gauss_hermite(Q)
    print("Q=%d  nodes  %s" % (Q, "  ".join("%+9.6f" % x for x in v)))
    print("      wts    %s" % "  ".join("%9.6f" % x for x in w))
 
tab = {Q: moment_errors(Q, 10) for Q in (1, 3, 5)}
print()
print(" n |    Q=1    |    Q=3    |    Q=5    |  exact")
for n in range(11):
    print("%2d | %9.2e | %9.2e | %9.2e | %6.1f"
          % (n, tab[1][n], tab[3][n], tab[5][n], exact_moment(n)))
print()
for Q in (1, 3, 5):
    bad = [n for n in range(11) if tab[Q][n] > 1e-9]
    first = bad[0] if bad else 11
    print("Q=%d: exact through degree %d, first failure at %d   (2Q-1 = %d)"
          % (Q, first - 1, first, 2 * Q - 1))

Q=3  nodes  -1.732051  +0.000000  +1.732051
      wts     0.166667   0.666667   0.166667
Q=5  nodes  -2.856970  -1.355626  +0.000000  +1.355626  +2.856970
      wts     0.011257   0.222076   0.533333   0.222076   0.011257
 
 n |    Q=1    |    Q=3    |    Q=5    |  exact
 0 |  0.00e+00 |  0.00e+00 |  2.22e-16 |    1.0
 1 |  3.93e-28 |  0.00e+00 |  1.39e-17 |    0.0
 2 |  1.00e+00 |  2.22e-16 |  4.44e-16 |    1.0
 3 |  6.05e-83 |  0.00e+00 |  0.00e+00 |    0.0
 4 |  3.00e+00 |  8.88e-16 |  1.78e-15 |    3.0
 5 | 9.33e-138 |  0.00e+00 |  0.00e+00 |    0.0
 6 |  1.50e+01 |  6.00e+00 |  1.24e-14 |   15.0
 7 | 1.44e-192 |  0.00e+00 |  0.00e+00 |    0.0
 8 |  1.05e+02 |  7.80e+01 |  8.53e-14 |  105.0
 9 | 2.22e-247 |  0.00e+00 |  0.00e+00 |    0.0
10 |  9.45e+02 |  8.64e+02 |  1.20e+02 |  945.0
 
Q=1: exact through degree 1, first failure at 2   (2Q-1 = 1)
Q=3: exact through degree 5, first failure at 6   (2Q-1 = 5)
Q=5: exact through degree 9, first failure at 10   (2Q-1 = 9)

三つの行がすべて $2Q-1$ で正確に切れています。偶然ではなく定理です。

15が9に落ちる算数#

Q=3 が6次で外れる過程は暗算で追えます。ノードは $0, \pm\sqrt{3}$ 、重みは $2/3, 1/6, 1/6$ です。6次モーメントを求積法で計算するとこうなります。

\frac{1}{6}(\sqrt{3})^6 + \frac{2}{3}\cdot 0^6 + \frac{1}{6}(-\sqrt{3})^6 = 2 \times \frac{27}{6} = 9

真値は $5!! = 15$ です。6がまるごと欠けています。ノードが三つしかないので、 $v^6$ のように裾の重い関数を載せる場所がありません。中央のノードは $0^6 = 0$ で何も寄与せず、外側の二つは $\pm\sqrt{3}$ に固定されているため、それより遠い領域の寄与は永久に失われます。

4次までは完璧です。 $2 \times 9/6 = 3 = 3!!$ 。この「4次まで正確、6次から破綻」がD2Q9の能力の限界をそのまま規定します。

$\sqrt{3}$ を1にする換算が $c_s^2 = 1/3$ の正体#

ここまでは純粋な数値解析です。LBMが加える条件は一つだけです。移流ステップが厳密でなければならないこと。 $f_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) = f_i(\mathbf{x}, t)$ が補間なしで成り立つには、 $\mathbf{c}_i \Delta t$ がちょうど格子点である必要があります。

ノードはすでに $0, \pm\sqrt{3}$ に決まっています。これを $0, \pm 1$ に送るには単位を $\sqrt{3}$ で割ればよく、その換算で熱速度の単位も一緒に縮みます。

c_s = \frac{1}{\sqrt{3}} \frac{\Delta x}{\Delta t}, \qquad c_s^2 = \frac{1}{3}

$c_s$ は格子音速です。つまり $1/3$ は物理定数ではありません。3次エルミート多項式の根を格子点の上に座らせると決めた瞬間に付いてきた単位換算です。5点則では根が $\pm 1.3556, \pm 2.8570$ となり整数比になりません。D2Q25系が補間や非整数格子を受け入れるか、根を少しずらした変種を使うのはそのためです。

1次元の重みを掛けるとD2Q9の重みが出る#

2次元は新しい導出ではありません。1次元則のテンソル積です。 $(1/6, 2/3, 1/6)$ を自分自身と掛ければ終わりです。

\frac{2}{3}\cdot\frac{2}{3} = \frac{4}{9}, \qquad \frac{2}{3}\cdot\frac{1}{6} = \frac{1}{9}, \qquad \frac{1}{6}\cdot\frac{1}{6} = \frac{1}{36}

静止方向が 4/9、軸方向4本が 1/9、対角4本が 1/36。LBM教科書の最初のページの表が、掛け算三回で出てきます。3次元なら同じ積をもう一度とって、D3Q27の 8/27、2/27、1/54、1/216 になります。D3Q19とD3Q15はテンソル積ではなく、そこから方向を間引いて重みを合わせ直した縮約版です。

平衡分布を2次で切ると何が落ちるか#

求積法の次数とは別に、展開を $N$ 次で切るという選択がもう一つあります。標準的なLBMは $N = 2$ です。

f_i^{eq} = w_i \rho \left[ 1 + \frac{\mathbf{c}_i \cdot \mathbf{u}}{c_s^2} + \frac{(\mathbf{c}_i \cdot \mathbf{u})^2 - c_s^2 u^2}{2 c_s^4} \right]

$w_i$ は上で求めた重み、 $\mathbf{c}_i$ は格子速度です。この形で3次・4次モーメントを測ると、二つの誤差源がそれぞれいくら請求しているかを分離して見られます。

from itertools import product
from math import factorial
 
def he(n, v):
    a, b = 1.0, v
    if n == 0:
        return a
    for k in range(1, n):
        a, b = b, v * b - k * a
    return b
 
def gauss_hermite(Q):
    nodes, x0, step = [], -8.0, 1e-3
    f0 = he(Q, x0)
    while x0 < 8.0:
        x1 = x0 + step
        f1 = he(Q, x1)
        if f0 == 0.0 or f0 * f1 < 0.0:
            a, b = x0, x1
            for _ in range(80):
                m = 0.5 * (a + b)
                if he(Q, a) * he(Q, m) <= 0.0:
                    b = m
                else:
                    a = m
            nodes.append(0.5 * (a + b))
        x0, f0 = x1, f1
    wts = [factorial(Q) / (Q * Q * he(Q - 1, v) ** 2) for v in nodes]
    return nodes, wts
 
def tensor_lattice(Q):
    """1D ガウス・エルミート則 -> 2D 直積則。熱単位では c_s = 1。"""
    n, w = gauss_hermite(Q)
    cs, ws = [], []
    for (vx, wx), (vy, wy) in product(zip(n, w), repeat=2):
        cs.append((vx, vy))
        ws.append(wx * wy)
    return cs, ws
 
def feq_hermite(rho, ux, c, w, order):
    """打ち切ったエルミート平衡分布, u = (ux, 0), c_s = 1"""
    out = []
    for (cx, _cy), wi in zip(c, w):
        s = 1.0 + cx * ux
        if order >= 2:
            s += 0.5 * ux * ux * (cx * cx - 1.0)
        if order >= 3:
            s += ux ** 3 * (cx ** 3 - 3.0 * cx) / 6.0
        out.append(wi * rho * s)
    return out
 
def raw_moment(f, c, p):
    return sum(fi * cx ** p for fi, (cx, _cy) in zip(f, c))
 
print("D2Q9 from the 3-node rule (rescaled by c_s = 1/sqrt(3)):")
c9, w9 = tensor_lattice(3)
for (cx, cy), wi in zip(c9, w9):
    ex, ey = round(cx / 3 ** 0.5), round(cy / 3 ** 0.5)
    print("   e = (%+2.0f,%+2.0f)   w = %9.6f  ~ %.0f/36" % (ex, ey, wi, round(wi * 36)))
 
c25, w25 = tensor_lattice(5)
print("\n   u    |  Pi_xxx exact | D2Q9 (N=2) |   gap    | D2Q25 (N=3) |  gap")
for u in (0.02, 0.05, 0.10, 0.20):
    ex = 3.0 * u + u ** 3
    m9 = raw_moment(feq_hermite(1.0, u, c9, w9, 2), c9, 3)
    m25 = raw_moment(feq_hermite(1.0, u, c25, w25, 3), c25, 3)
    print("  %.2f  |   %10.6f  | %10.6f | %+8.1e |  %10.6f | %+8.1e"
          % (u, ex, m9, m9 - ex, m25, m25 - ex))
 
print("\n   u    | Pi_xxxx exact | D2Q9 (N=2) |   gap    | D2Q25 (N=3) |  gap")
for u in (0.02, 0.05, 0.10, 0.20):
    ex = 3.0 + 6.0 * u ** 2 + u ** 4
    m9 = raw_moment(feq_hermite(1.0, u, c9, w9, 2), c9, 4)
    m25 = raw_moment(feq_hermite(1.0, u, c25, w25, 3), c25, 4)
    print("  %.2f  |   %10.6f  | %10.6f | %+8.1e |  %10.6f | %+8.1e"
          % (u, ex, m9, m9 - ex, m25, m25 - ex))

D2Q9 from the 3-node rule (rescaled by c_s = 1/sqrt(3)):
   e = (-1,-1)   w =  0.027778  ~ 1/36
   e = (-1,+0)   w =  0.111111  ~ 4/36
   e = (-1,+1)   w =  0.027778  ~ 1/36
   e = (+0,-1)   w =  0.111111  ~ 4/36
   e = (+0,+0)   w =  0.444444  ~ 16/36
   e = (+0,+1)   w =  0.111111  ~ 4/36
   e = (+1,-1)   w =  0.027778  ~ 1/36
   e = (+1,+0)   w =  0.111111  ~ 4/36
   e = (+1,+1)   w =  0.027778  ~ 1/36
 
   u    |  Pi_xxx exact | D2Q9 (N=2) |   gap    | D2Q25 (N=3) |  gap
  0.02  |     0.060008  |   0.060000 | -8.0e-06 |    0.060008 | +1.8e-16
  0.05  |     0.150125  |   0.150000 | -1.3e-04 |    0.150125 | +2.8e-16
  0.10  |     0.301000  |   0.300000 | -1.0e-03 |    0.301000 | +2.2e-16
  0.20  |     0.608000  |   0.600000 | -8.0e-03 |    0.608000 | +4.4e-16
 
   u    | Pi_xxxx exact | D2Q9 (N=2) |   gap    | D2Q25 (N=3) |  gap
  0.02  |     3.002400  |   3.001200 | -1.2e-03 |    3.002400 | -1.6e-07
  0.05  |     3.015006  |   3.007500 | -7.5e-03 |    3.015000 | -6.2e-06
  0.10  |     3.060100  |   3.030000 | -3.0e-02 |    3.060000 | -1.0e-04
  0.20  |     3.241600  |   3.120000 | -1.2e-01 |    3.240000 | -1.6e-03

QとNを一緒に上げると、表の $-8.0\times 10^{-3}$ がどの列で $10^{-16}$ まで落ちるかが分かれます。

Q per axis = 3 Hermite order N = 2 u = 0.20 
Pi_xxx gap 0.00e+0  ·  Pi_xxxx gap 0.00e+0

Qだけを5にすると格子はD2Q25になりますが、 $N=2$ の平衡は $u^3$ をまだ落とします。Nだけを3にしてもノードが三つなら $c_x^3=c_x$ がその項を消します。両方上げて初めて、D2Q25の列のように隙間が機械精度まで下がります。

重みの表がテンソル積からそのまま再現されています。その下の二つの表が本題です。

$u = 0.2$ で3.2416が3.1200になる#

$\Pi_{xxx}$ の欠損はちょうど $-\rho u_x^3$ です。 $u = 0.2$ で $-8.0 \times 10^{-3}$ 、 $u = 0.1$ で $-1.0 \times 10^{-3}$ 。速度を半分にすると誤差は8分の1になります。3次です。

理由は格子が $c_x^3 = c_x$ を強制するからです。 $c_x$ が $-1, 0, 1$ しかないので、3乗しても自分自身に戻ります。そのため $\sum f_i c_x^3$ は $\sum f_i c_x = \rho u_x$ を返しますが、連続系が要求する値は $\rho(3c_s^2 u_x + u_x^3) = \rho(u_x + u_x^3)$ です。 $u_x^3$ が丸ごと抜けています。同じ制約は力の項にも等しくかかります。そちらは力の項が合わせるべき3次モーメントで別に扱いました。

4次はもっと大きくなります。 $u = 0.2$ で真値3.2416に対して3.1200、相対誤差3.75%です。欠損は $-\rho(u_x^2 + u_x^4)$ で、2次から始まります。

$u$ を 0.02、0.05、0.10、0.20 と四つ取った欄がその冪を見せます。 $\Pi_{xxx}$ の隙間は $-8.0\times 10^{-6}$ 、 $-1.3\times 10^{-4}$ 、 $-1.0\times 10^{-3}$ 、 $-8.0\times 10^{-3}$ です。0.05から0.10へ速度を2倍にすると、隙間は $-1.3\times 10^{-4}$ から $-1.0\times 10^{-3}$ へ、およそ8倍になります。0.10から0.20も同じです。 $-1.0\times 10^{-3}$ が $-8.0\times 10^{-3}$ になります。連続系が要求する 0.608000 に対して格子は 0.600000 を返します。差は $u^3=0.008$ です。

$\Pi_{xxxx}$ はもっと早く開きます。 $u=0.02$ では真値 3.002400、D2Q9 は 3.001200、隙間は $-1.2\times 10^{-3}$ です。相対では 0.04% です。 $u=0.20$ では 3.241600 が 3.120000 になり、隙間は $-1.2\times 10^{-1}$ です。D2Q25 に $N=3$ を与えても $-1.6\times 10^{-3}$ が残ります。 $0.2^4=0.0016$ と同じ大きさで、符号だけが逆です。ノードを増やした代償はこの項に届きません。

二つの誤差源は互いに独立

表のD2Q25の列がこの点を分離して見せてくれます。 $\Pi_{xxx}$ の隙間は $10^{-16}$ 、機械精度まで落ちます。求積次数9と展開次数3のどちらも十分だからです。

ところが同じD2Q25でも $\Pi_{xxxx}$ は $u = 0.2$ で $-1.6 \times 10^{-3}$ が残ります。これはちょうど $-u^4$ です。求積法のせいではなく、展開を3次で切ったために生じた値です。ノードをいくら増やしてもこの項は減りません。 $N$ を4に上げて初めて消えます。

誤差の出どころは二つに分かれます。展開次数 $N$ は平衡分布が持てる $u$ のべきを決め、求積次数 $2Q-1$ はその多項式を厳密に積分できる上限を決めます。 $M$ 次モーメントを正しく得るには $2Q-1 \ge N + M$ が必要です。D2Q9は $2Q-1 = 5$ 、 $N = 2$ なので、3次モーメントがすでに境界線上にあります。

この不等式が実際に塞いでいること

最初に塞がるのはガリレイ不変性です。 $\Pi_{xxx}$ の $u^3$ の欠損はChapman–Enskog展開で $\partial_x(\rho u_x^3)$ の形の偽の項として蘇ります。平均流を乗せた状態でせん断波の減衰率を測ると、粘性が流速によって変わって見えます。離散化が請求する半ステップと重なって、測定粘性と入力 $\tau$ の関係が二重に曇ります。

熱流動は一次数ぶん高くつきます。エネルギー流束は3次モーメントであり、その平衡値を合わせるには $N = 3$ 、したがって $2Q-1 \ge 6$ が必要です。D2Q9は5で止まります。標準のD2Q9が等温モデルとしてしか使われないのは、この一行が理由です。平衡分布を多項式で切る側は9本の速度の中のマクスウェル分布に続きます。

高マッハでは欠損項が $u^3$ 、 $u^4$ で育ちます。 $\Pi_{xxxx}$ の誤差は $u = 0.02$ で0.04%ですが、 $u = 0.2$ では3.75%です。「格子マッハ数を0.1以下に」という経験則は安定性だけの問題ではなく、この表の数字です。D2Q9のまま $u^3$ の欠損だけを補正項で戻す道もあります。 $c_x^3 = c_x$ が消した成分を空間微分で近似して足し直します。微分が一つ増えるので境界付近で次数が落ち、非平衡部分との干渉も改めて検討する必要があります。

格子表を選ぶときに実際に署名しているもの

e[9][2] と w[9] をコードに書く作業は、速度集合を選ぶことではありません。積分次数に署名する作業です。ノードを三つ選んだ瞬間に「5次までしか正確でない」という契約が結ばれ、その契約が等温・低マッハ・ガリレイ近似という三条件になって戻ってきます。

次にD3Q19とD3Q27で迷ったら、最初に問うべきはメモリではありません。このコードが必要とする最も高いモーメントは何次か、そして平衡分布を何次で切るか。 $2Q-1 \ge N + M$ の一行があれば、残りは付いてきます。