格子ペクレ数を2に上げたら解が-0.33まで落ちた — 最小化すべきエネルギーが消えた場所

安定化は対称性を取り戻さない。対称部分を膨らませて、代わりに離散最大値原理を買っているだけです。

一本の棒では四つの方法が同じ答えを出した

有限要素法の講義ノート第1章は、同じ棒の問題を四回解きます。直接剛性法、最小ポテンシャルエネルギー、重み付き残差法、そしてガラーキン法です。四回とも同じ5×5の剛性行列が出てきます。ノートにはこう書かれています。どの方法を使っても結果に大きな違いはなく、正確な値が得られる、と。

しかしこの一文には条件が付いています。その条件は棒の問題の中に隠れていて表に出てきません。流れの方程式はその条件を満たしません。この記事では条件が何であるか、そして条件が壊れたときに正確にどの保証を失うのかを測ります。さらに安定化手法が何を返してくれるのか — そして何は決して返してくれないのか — を数字で確認します。

最小化するには行列が対称でなければならない

最小ポテンシャルエネルギー法は、全ポテンシャルエネルギーを節点未知数で微分してゼロと置きます。離散化した形はこうなります。

\Pi(\boldsymbol{\phi}) = \frac{1}{2}\,\boldsymbol{\phi}^{T}\mathbf{K}\,\boldsymbol{\phi} - \mathbf{f}^{T}\boldsymbol{\phi}

$\boldsymbol{\phi}$ は節点未知数ベクトル、 $\mathbf{K}$ は剛性行列、 $\mathbf{f}$ は荷重ベクトルです。ここで $i$ 番目の成分で微分します。

\frac{\partial \Pi}{\partial \phi_i} = \frac{1}{2}\sum_{j}\left(K_{ij} + K_{ji}\right)\phi_j - f_i

$K_{ij} + K_{ji}$ が出てきたことに注目してください。この式が $\mathbf{K}\boldsymbol{\phi} = \mathbf{f}$ になるためには $K_{ij} = K_{ji}$ である必要があります。対称でなければ、最小化が実際に解いているのは対称化された行列 $\tfrac{1}{2}(\mathbf{K}+\mathbf{K}^{T})$ であり、それは元の方程式ではありません。

同じ場所からもっと根本的な話も出てきます。残差場 $\mathbf{r}(\boldsymbol{\phi}) = \mathbf{f} - \mathbf{K}\boldsymbol{\phi}$ が何らかのスカラー関数の勾配であるためには、そのヤコビアンが対称でなければなりません。対称でなければポテンシャル関数はそもそも存在しません。存在を確かめる一番速い方法は閉じた経路を一周することです。勾配場なら一周した仕事は必ずゼロになります。

下のシミュレーションで、自由度2個のおもちゃの行列を使って確かめてみましょう。

advection a1.00diffusion s1.00

Set advection a to 0: the pink dot goes all the way around and brings back 0.000 — the field is a gradient, and the amber ball slides straight down the ellipses. Push a up and the last lap returns 6.283, exactly 2πa. Now drag diffusion s across its whole range: the lap total does not budge. No amount of added diffusion buys back a potential.

advection aをゼロにすると、ピンクの点は一周して0.000を持ち帰ります。灰色の楕円は等エネルギー線で、黄色の球はその線を横切らずまっすぐ内側へ滑り込みます。aを上げると一周の仕事はちょうど $2\pi a$ になり、球は螺旋を描いて等エネルギー線を逆に登り始めます。この状態でdiffusion sを端まで動かしてみてください。一周の仕事はまったく動きません。

移流項が壊す量はちょうど $u$ です#

1次元移流拡散方程式 $u\,\phi' = \epsilon\,\phi''$ の弱形式を書くと、二つの項の性格の違いが見えてきます。

a(w, \phi) = \int_0^L \left( \epsilon\,\frac{dw}{dx}\frac{d\phi}{dx} + w\,u\,\frac{d\phi}{dx} \right) dx

$w$ は試験関数、 $\epsilon$ は拡散係数、 $u$ は移流速度です。第1項は $w$ と $\phi$ を入れ替えても変わりません — 対称です。第2項は部分積分すると符号が反転します。境界で消える試験関数に対して $\int w\,u\,\phi'\,dx = -\int \phi\,u\,w'\,dx$ ですから、移流項は純粋に反対称な成分です。

線形要素で組み立てると、この構造がそのまま係数に残ります。要素サイズを $h$ とすると、内部節点の行はこうなります。

K_{i,i-1} = -\frac{\epsilon}{h} - \frac{u}{2}, \qquad K_{i,i} = \frac{2\epsilon}{h}, \qquad K_{i,i+1} = -\frac{\epsilon}{h} + \frac{u}{2}

拡散は両隣に同じ値を入れ、移流は片側に $+u/2$ 、反対側に $-u/2$ を入れます。だから対称からのずれはちょうど一つの数になります。

K_{i,i+1} - K_{i+1,i} = u

格子をいくら細かくしてもこの値は $u$ のままです。 $h$ が入っていないからです。流れがある限り、最小ポテンシャルエネルギー原理は戻ってきません。曲線座標系で構成テンソルを変換した記事では、対称性のおかげで剛性を6×6のVoigt行列に畳めましたが、事情が違います。

Pythonで同じ格子に三つの手法を立てた#

要素10個、 $u = 1$ 、境界条件 $\phi(0)=0$ と $\phi(1)=1$ 。人工拡散係数 $\beta$ 一つで三つのスキームを作ります。 $\epsilon_{\text{eff}} = \epsilon + \beta\,u\,h/2$ として、 $\beta = 0$ なら純粋ガラーキン、 $\beta = 1$ なら完全風上差分、 $\beta = \coth(Pe_h) - 1/Pe_h$ ならSUPGです。

import numpy as np
 
def assemble_ad(n, vel, eps, beta):
    """線形要素上の1次元移流拡散。beta は人工拡散係数。"""
    h = 1.0 / n
    eps_eff = eps + beta * vel * h / 2.0
    kd = (eps_eff / h) * np.array([[1.0, -1.0], [-1.0, 1.0]])   # 拡散: 対称
    ka = (vel / 2.0) * np.array([[-1.0, 1.0], [-1.0, 1.0]])     # 移流: 反対称
    K = np.zeros((n + 1, n + 1))
    for e in range(n):
        K[e:e + 2, e:e + 2] += kd + ka
    return K
 
def solve_bvp(K):
    n = K.shape[0] - 1
    A, b = K.copy(), np.zeros(n + 1)
    A[0, :], A[0, 0], b[0] = 0.0, 1.0, 0.0
    A[n, :], A[n, n], b[n] = 0.0, 1.0, 1.0
    return np.linalg.solve(A, b)
 
def exact_ad(x, pe):
    return (np.exp(pe * (x - 1.0)) - np.exp(-pe)) / (1.0 - np.exp(-pe))
 
def skew_ratio(K):
    return np.linalg.norm(K - K.T) / np.linalg.norm(K + K.T)
 
def loop_work(K, m=20000):
    """自由度空間の単位円を一周したときの残差場 -K.phi の仕事。"""
    t = np.linspace(0.0, 2.0 * np.pi, m, endpoint=False)
    path = np.stack([np.cos(t), np.sin(t)])          # 位置
    tang = np.stack([-np.sin(t), np.cos(t)])         # dl / dt
    return float(np.sum(np.sum(-(K @ path) * tang, axis=0)) * (2.0 * np.pi / m))
 
n, vel = 10, 1.0
x = np.linspace(0.0, 1.0, n + 1)
print(f"{'Pe_h':>5} {'scheme':>9} {'max|err|%':>10} {'min phi':>9} {'skew/sym':>9} {'loop W':>8}")
for pe_h in [0.5, 1.0, 2.0, 5.0]:
    eps = vel / (2.0 * n * pe_h)
    ex = exact_ad(x, vel / eps)
    for name, beta in [("Galerkin", 0.0), ("upwind", 1.0),
                       ("SUPG", 1.0 / np.tanh(pe_h) - 1.0 / pe_h)]:
        K = assemble_ad(n, vel, eps, beta)
        phi = solve_bvp(K)
        print(f"{pe_h:5.1f} {name:>9} {100 * np.max(np.abs(phi - ex)):10.2f} "
              f"{phi.min():9.4f} {skew_ratio(K):9.4f} {loop_work(K[4:6, 4:6]):8.4f}")
 
K0 = assemble_ad(n, 0.0, 0.1, 0.0)
print(f"\nvel = 0 : skew/sym = {skew_ratio(K0):.2e},  loop W = {loop_work(K0[4:6, 4:6]):.2e}")
print(f"pi * u  = {np.pi * vel:.4f}")

 Pe_h    scheme  max|err|%   min phi  skew/sym   loop W
  0.5  Galerkin       3.45    0.0000    0.2924   3.1416
  0.5    upwind      13.17    0.0000    0.1954   3.1416
  0.5      SUPG       0.00   -0.0000    0.2704   3.1416
  1.0  Galerkin      13.53    0.0000    0.5774   3.1416
  1.0    upwind      19.80    0.0000    0.2924   3.1416
  1.0      SUPG       0.00   -0.0000    0.4428   3.1416
  2.0  Galerkin      35.17   -0.3334    1.1010   3.1416
  2.0    upwind      18.17   -0.0000    0.3885   3.1416
  2.0      SUPG       0.00    0.0000    0.5572   3.1416
  5.0  Galerkin      69.61   -0.6961    2.1517   3.1416
  5.0    upwind       9.09    0.0000    0.4836   3.1416
  5.0      SUPG       0.00    0.0000    0.5773   3.1416
 
vel = 0 : skew/sym = 0.00e+00,  loop W = -4.29e-16
pi * u  = 3.1416

境界値は0と1に固定されているのに、 $Pe_h = 2$ のガラーキン解は $-0.3334$ まで下がります。 $Pe_h = 5$ では $-0.6961$ です。物理的に出てはいけない値です。

$Pe_h = 1$ で係数が一つ符号を変える#

格子ペクレ数を $Pe_h = u h / (2\epsilon)$ と定義すると、先ほどの行の右側の係数はこう書き直せます。

K_{i,i+1} = \frac{\epsilon}{h}\left(Pe_h - 1\right)

$Pe_h > 1$ になると非対角成分が正になります。この瞬間に行列はM行列であることをやめ、離散最大値原理 — 内部の解が境界値の作る区間を出ないという保証 — も一緒に消えます。表でmin phiが $Pe_h = 1.0$ まではゼロで、2.0で負になる理由がこれです。

下のシミュレーションは同じ系を時間前進で解きます。定常解ができあがる過程を見ると、振動がどこから育つのかが見えます。

Pe_h2.0beta0.00

Leave beta at 0 and drag Pe_h past 1: a_E flips sign, and the marching profile starts ringing — at Pe_h = 2 the node next to the outlet dives to 0.000 with the boundary values still pinned at 0 and 1. Hit SUPG and the error goes to 0.00 %. The pink line at bottom right is the skew part of the matrix: none of the three buttons moves it.

betaをゼロにしたままPe_hを1の向こうへ押してみてください。右のa_Eバーが左に越えて赤くなり、次の時間ステップで節点値が0の下に沈みます。SUPGボタンを押すと誤差が0 %になります。ところが右下のピンクの行 — 行列の反対称成分 — は三つのボタンのどれでも動きません。

安定化は対称性を取り戻さない

ここでよくある誤解を整理しておきます。風上差分やSUPGが「安定性を回復させる」という言い方は、対称性と最小化原理を取り戻すという意味ではありません。表のloop W列を見てください。三つのスキーム、四つの $Pe_h$ — 全12行がすべて3.1416です。この値は $\pi u$ であり、 $\epsilon$ も $\beta$ も入っていません。

理由は単純です。 $\beta$ が大きくするのは $\epsilon_{\text{eff}}$ だけで、それは行列の対称部分です。反対称部分の $u$ には手を触れません。skew/sym比が $Pe_h = 5$ で2.1517から0.4836に減ったのも、反対称が小さくなったのではなく、分母の対称部分が大きくなったからです。

つまり安定化が実際に買っているのは、より弱い保証です。最小化原理(エネルギーノルムでの最良近似)を失った代わりに、M行列の性質と離散最大値原理を手に入れます。代金は精度で払います。 $Pe_h = 2$ で風上差分は振動を消しましたが、誤差は18.17 %です。一方SUPGは $\beta$ を必要な分だけ入れるので節点で正確です。

\beta_{\text{opt}} = \coth(Pe_h) - \frac{1}{Pe_h}

この値は $Pe_h \to 0$ でゼロへ、 $Pe_h \to \infty$ で1へ向かいます。拡散が支配するなら安定化を切り、移流が支配するなら完全風上差分へ行くという意味です。ただし節点での正確さは1次元定数係数問題の特権です。2次元では流線方向にだけ人工拡散を入れる本来のSUPGの形が必要で、そのときも「正確」は消えます。

有限体積法が同じ場所で使う別の名前

この計算は有限要素法の言葉で行いましたが、結論は離散化手法に依りません。有限体積法でも中心差分の対流項はまったく同じステンシルを作り、同じ $Pe_h = 1$ で係数が符号を変えます。そこで登場するのが1次風上差分で、それが持ち込む数値拡散は $\beta = 1$ のときの $u h / 2$ と同じ値です。

二つの世界が分かれるのは保存性です。有限体積法の風上差分は面フラックスを変える方式なので、総量保存が保たれます。保存形と原始形が衝撃波速度で分かれた話がここでも繰り返されます。有限要素法の人工拡散は剛性行列に項を足すので、何を保存しているかを別に確認する必要があります。

重み付き残差の系列の中で別の選択肢を選ぶと、また別の代価を払います。最小二乗法は残差ノルムを最小化するので、常に対称正定値行列を作ります — 最小化原理が戻ってきます。その代わり条件数が二乗され、線形要素では要素内部の2階微分がゼロなので拡散項がまるごと消えます。不連続ガラーキンの積分点の下限を検討した記事と同じく、基底と積分則が定式化の性質を静かに変える場所です。

非対称な行列を引き継いだとき

講義ノートの「どの方法でも同じ」という文の有効範囲は自己随伴作用素です。拡散、弾性、ポテンシャル流はその中にあります。移流が入った瞬間に外へ出ます。

実務では三つを順に確認すれば十分です。第一に、組み立てた行列は対称か。対称ならCG系が使え、エネルギーノルムでの最良近似という保証が付いてきます。非対称ならGMRES系で、その保証はありません。第二に、セルペクレ数は1を超えているか。超えていれば振動はバグではなく、スキームの定義された挙動です。第三に、安定化を入れたなら、それが買ったのは精度か有界性か。たいていは有界性で、精度は支払った側です。

解が境界値の区間の外へ出たとき、格子を細かくするのは応急処置ではなく正攻法です。 $h$ を小さくすれば $Pe_h$ も一緒に小さくなるからです。ただし3次元で $Pe_h$ を5から1まで下げるには、セル数が125倍になります。その計算をしてから安定化項を選ぶほうがよいでしょう。