格子が歪むとラプラシアンが漏れる — 非直交拡散フラックス補正

直交補正の壁は70°ではなく、cosθが1/2になる60°にあります。

立方体で安定していた拡散ソルバーを、歪んだ三角形へ移すと、答えが静かにずれます。方程式は同じです。コードも同じです。温度場だけがぼけ、収束は遅くなります。境界条件でも物性でもありません。面の上で勾配を測るやり方です。

非直交角のスライダーを20°から65°へ押すと、過緩和の青緑ベクトルは面積より長くなり、最小補正のピンクのベクトルはさらに長くなります。

non-orthogonality θ = 35°

over-relaxed   |E_f|/|S_f| = 1.221   |T_f|/|S_f| = 0.700
explicit/implicit ratio ρ = |T_f|/|E_f| :   min 0.70  |  orth 0.60  |  over 0.57

青緑がピンクより短くなると、画面下の ρ が 1 を超えます。45°で最小補正だけが先に赤へ向かうかを見ます。

面の勾配が厄介者

有限体積法(FVM、保存量をセルごとに積分する方法)では、拡散項は面を通るフラックスの和です。面 $f$ 一つの拡散フラックスは次のように書きます。

F_f = \Gamma_f \, \nabla\phi_f \cdot \mathbf{S}_f

$\Gamma_f$ は面の拡散係数です。 $\nabla\phi_f$ は面の勾配です。 $\mathbf{S}_f$ は外向きの面積ベクトル、面積かける法線です。

$\nabla\phi_f \cdot \mathbf{S}_f$ をセル値へ落とすのが問題です。構造直交格子では、二つの中心を結ぶ $\mathbf{d}$ が法線と平行です。 $(\phi_N-\phi_P)/|\mathbf{d}|$ がそのまま法線勾配になります。

非構造格子ではその仮定が崩れます。 $\mathbf{d}$ と $\mathbf{S}_f$ がずれます。その角が非直交角 $\theta$ です。 $(\phi_N-\phi_P)/|\mathbf{d}|$ は、もう法線勾配ではなく $\mathbf{d}$ 方向の勾配です。差を無視すると拡散が漏れ、解はぼけます。

S_f を二つに割る#

面積ベクトルを二つに分けます。一つは $\mathbf{d}$ 方向、一つは残りです。

\mathbf{S}_f = \mathbf{E}_f + \mathbf{T}_f

$\mathbf{E}_f$ は $\mathbf{d}$ に沿う成分です。 $\mathbf{T}_f$ は接線側の残りです。フラックスへ入れると項が二つに割れます。

F_f = \Gamma_f |\mathbf{E}_f| \frac{\phi_N - \phi_P}{|\mathbf{d}|} + \Gamma_f \, \overline{\nabla\phi}_f \cdot \mathbf{T}_f

前項は $\phi_N-\phi_P$ だけを使います。陰的に行列へ入れます。セル中心の二点だけです。後項は面平均勾配 $\overline{\nabla\phi}_f$ が要ります。隣接勾配の平均です。係数に入れにくいので、陽的補正として右辺へ投げます。その勾配を最小二乗で組み直す側は非構造格子で勾配を残す移動最小二乗に別にあります。遅延補正が揺れるとき、犯人は ρ ではなくその勾配であることがあります。

$\theta=0$ なら $\mathbf{T}_f=0$ です。昔の中心差分へ戻ります。

最小・直交・過緩和：三つの分解

$\mathbf{E}_f$ を $\mathbf{d}$ 方向に置く長さは一つではありません。 $\hat{\mathbf{d}}=\mathbf{d}/|\mathbf{d}|$ 、 $\cos\theta=\hat{\mathbf{d}}\cdot\mathbf{S}_f/|\mathbf{S}_f|$ とすると三つあります。

最小補正: $|\mathbf{E}_f|=|\mathbf{S}_f|\cos\theta$ 。 $\mathbf{T}_f$ が最も短いです。
直交補正: $|\mathbf{E}_f|=|\mathbf{S}_f|$ 。陰的な大きさを面積のままにします。
過緩和: $|\mathbf{E}_f|=|\mathbf{S}_f|/\cos\theta$ 。非直交が大きくなるほど陰的成分を伸ばします。

上の画面で θ を 60°まで上げると、過緩和の $\mathbf{E}_f$ は面積ベクトルより長くなります。対角係数が大きくなった合図です。最小補正は逆で、 $\mathbf{E}_f$ が縮み $\mathbf{T}_f$ が伸びます。

コード：遅延補正ラプラスソルバー

陽的補正は直前の反復の勾配を使います。一回解いて終わりではありません。補正を更新しながら何度も掃く、遅延補正のループです。

import numpy as np
 
def decompose_face(S_f, d, scheme="over"):
    """面積ベクトル S_f = E_f（陰的、d 方向） + T_f（陽的補正）。"""
    d_hat = d / np.linalg.norm(d)
    Smag = np.linalg.norm(S_f)
    dS = d_hat @ S_f                      # = |S_f| cos(非直交角)
    if scheme == "min":                   # 最小補正
        E = dS * d_hat
    elif scheme == "orth":                # 直交補正
        E = Smag * d_hat
    else:                                 # 過緩和 (over-relaxed)
        E = (Smag**2 / dS) * d_hat
    return E, S_f - E                     # E_f, T_f
 
def sweep_residual(S_f, d, scheme, sweeps=25):
    """E_f は陰的、T_f は遅延。残差は rho = |T|/|E| の比で縮む。"""
    E, T = decompose_face(S_f, d, scheme)
    rho = np.linalg.norm(T) / np.linalg.norm(E)
    r, hist = 1.0, [1.0]
    for _ in range(sweeps):
        r *= rho                          # 遅延補正 1 回 = 残差 x rho
        hist.append(r)
    return rho, hist
 
S_f = np.array([1.0, 0.0])                # 面法線 (面積ベクトル, |S|=1)
for ang in (20, 45, 65):
    d = np.array([np.cos(np.radians(ang)), np.sin(np.radians(ang))])
    for s in ("min", "orth", "over"):
        rho, hist = sweep_residual(S_f, d, s)
        tag = "収束" if rho < 1 else "発散"
        print(f"theta={ang:2d}  {s:5s}  rho={rho:.3f}  {tag}  (25回後 {hist[-1]:.1e})")

出力は、角が大きくなるほど三方式が割れるという意味です。20°では三つとも穏当です。65°で最小補正は $\rho=\tan 65°\approx 2.14$ で発散し、過緩和は $\rho=\sin 65°\approx 0.91$ でまだ収束します。同じ格子、同じ方程式なのに、分解の一つが収束と発散を分けます。

三つの ρ は角度だけで閉じます。面積の長さを 1 とすると、最小補正は E の長さが cosθ、T の長さが sinθ なので ρ は tanθ です。tan 45° は 1 です。45°が収束の壁です。その角では遅延ループが残差に 1 を掛けます。25回回しても数字はそのままです。20°の tan は約 0.364 です。65°の tan は約 2.14 です。発散側の数字はこの正接から出ます。

直交補正は E の長さが 1 で、T の長さは 2 から 2cosθ を引いた値の平方根です。ρ が 1 になるには cosθ が 1/2、つまり 60°です。65°ではすでに約 1.07 です。70°まで待って直交が危ないと言うのは遅いです。壁は 60°にあります。

過緩和の ρ は sinθ です。65°で約 0.906 です。スニペットは残差にこの比を 25 回掛けます。0.906 の 25 乗は約 0.085 です。25回の後も陽的誤差の約 8 パーセントが残ります。同じ式を 20°に使うと sin は約 0.342 で、25 乗は $10^{-12}$ の水準まで消えます。補正回数を 2 か 3 にする習慣は、65°の面では 0.906 の 3 乗、約 0.74 です。誤差の約 4 分の 3 がその時間ステップに残ります。

コードの角は d と S の間の角です。S を単位ベクトルに固定すると、20°、45°、65°が比較点です。最小補正の 45°は tan がちょうど 1 なので、経験で選んだ閾値ではありません。

なぜ ρ が T の長さ割り E の長さなのか。遅延補正は陰的係数、E の長さで解き、掃くたびに T の長さだけの誤差を陽的に戻します。一回掃くごとに残差がその比で減るか増えます。比が 1 より小さくないとループは閉じません。

非直交が増すと何が起こるか

三方式の残差を一つの画面で一緒に進めます。角のスライダーが ρ を変え、掃くたびに曲線がその比で下がります。

mesh non-orthogonality θ = 55°

● minimum ρ=1.43 diverges● orthogonal ρ=0.92 converges● over-relaxed ρ=0.82 converges

30°より下では三つが一緒に落ちます。45°では最小補正だけが水平になります。60°を越えると直交も ρ が 1 を越え、過緩和だけが sinθ のまま残ります。

代償はあります。過緩和は陰的成分を伸ばす分、陽的補正も大きいです。同じ残差へ届くには掃きを何度も回します。nNonOrthogonalCorrectors がその回数です。格子が良ければ 0 か 1 で足ります。ひどく潰れていると 2 か 3 まで上げます。65°ではその二、三回が誤差の 4 分の 3 を残します。むやみに増やすと時間ステップが高くなります。角そのものを減らす方が、ほとんどいつも安いです。

45°で最小補正を捨てる理由#

45°で最小補正は残差を減らせません。ρ がちょうど 1 だからです。直交は 60°で同じ壁に当たり、65°ではすでに 1 を越えています。過緩和は sinθ なので、90°の直前まで ρ は 1 を越えません。代わりに一歩の中に残る誤差が大きいです。

回数を足すのか、角が 60°を越える前に格子を敷き直すのか。三角形で覆う順序がその角を作ります。点群を三角形で覆う二つの生成法は、その順序がどこで非直交を残すかを問います。補正ループを伸ばす方が安いのか、フロントを再び押す方が安いのか。