歪んだセルで勾配がずれる — Green–Gauss と Least-Squares 再構成

境界ステンシルで最小二乗を切った瞬間、内部が守っていた一桁の誤差も次数ごと落ちます。

同じ隣接六つ、同じ圧力です。一方は $\nabla p$ をほぼ指し、もう一方はかなり横へ傾きます。差は、セルがどれだけ歪んだかと、面の圧力をどう作ったかにあります。

mesh skew のスライダーを0から0.9へ押すと、橙の Green–Gauss 矢印が真値から離れ、紫の最小二乗はそれより短く残ります。

mesh skew: 0.00wavenumber k: 1.20

relative error (Green–Gauss): 44.4%
relative error (Least-Squares): 44.4%

Field: p(x, y) = sin(k·x) · cos(k·y). Drag skew → 0 for a near-regular ring; → 0.9 for a strongly distorted polygon mesh.

skew 0 の近くでは両方とも一桁です。0.5を超えると橙だけが二桁へ上がります。表の数字はコードの後にあります。

OpenFOAM のログで圧力が発散した日#

非スタガードの非構造・圧力ベースのソルバは、Rhie–Chow の圧力加重内挿でセル中心の $\nabla p$ を少なくとも三回呼びます。運動量の圧力項、面の中間速度、圧力補正の勾配です。一つでも片側へ傾くと運動量と連続がずれ、残差が一桁上で止まります。二相では相変化率が局所圧力の関数なので、 $\nabla p$ が0.5% ずれると相境界がぼけます。

問題はセル面です。セル中心の勾配は面積分から落ち、面圧力 $p_f$ の作り方が道を分けます。

Green–Gauss — 面積分一行#

ガウスの定理から出発します。セル $C$ の体積 $\Omega_C$ 、面 $f$ 、外向き単位法線 $\hat{\mathbf{n}}_f$ 、面積 $A_f$ です。

\nabla p \Big|_C \;\approx\; \frac{1}{\Omega_C} \sum_{f \in \partial C} p_f\, \hat{\mathbf{n}}_f\, A_f

$p_f$ は面の圧力、 $\hat{\mathbf{n}}_f$ は外向き法線、 $A_f$ は面積です。

$p_f$ の作り方が Green–Gauss の変種を分けます。

単純平均: $p_f=\tfrac12(p_C+p_N)$ です。直交で一様ならきれいです。非構造では面中心が二つのセル中心を結ぶ線分の真ん中にありません。
距離の逆重み: $p_f=(w_C p_C+w_N p_N)/(w_C+w_N)$ で、 $w=1/d$ です。面に近いセルが大きな重みを受けます。
Frink: 格子点圧力を隣接セルから pseudo-Laplacian 重みで内挿したあと、面を成す二点でもう一度平均します。保存と2次を同時に狙います。

セルが歪むと面中心 $\mathbf{x}_f$ が二中心を結ぶ線から外れます。 $p_f$ にどの平均を入れても面上の分布を取り逃します。1次さえ保証されないことがよくあります。分母の $\Omega_C$ が偽の面積なら、保存は形だけです。対角の外積でセル面積を閉じる道が、その面積が $\partial x/\partial x=1$ とかみ合う場所です。

同じ一点、二つのアルゴリズム、違う矢印

前の格子の環は中心セル一つと隣接六つです。skew 0 では二つの矢印が真値のシアンに付きます。0.5を超えると橙が先に離れます。波数 $k$ を上げると真の $\nabla p$ がセル寸法に対して急になり、両方とも傷みます。そのとき紫の優位は狭くなります。

Least-Squares — 点群に直接フィットする#

最小二乗は面を通りません。隣接が線形分布を満たすとみて、過剰決定系を解きます。

p_N - p_C \;=\; (\nabla p)_C \cdot (\mathbf{x}_N - \mathbf{x}_C) + \varepsilon_N, \quad N=1,\ldots,n_{NC}

各行は残差 $\varepsilon_N$ を残します。重み $w_N$ はふつう $1/|\mathbf{x}_N-\mathbf{x}_C|^2$ です。 $\sum w_N\varepsilon_N^2$ を最小にする $\nabla p$ が正規方程式から出ます。

2D の正規方程式は 2×2 です。

\begin{bmatrix} \sum w_N \Delta x_N^2 & \sum w_N \Delta x_N \Delta y_N \\ \sum w_N \Delta x_N \Delta y_N & \sum w_N \Delta y_N^2 \end{bmatrix} \begin{bmatrix} \partial_x p \\ \partial_y p \end{bmatrix} = \begin{bmatrix} \sum w_N \Delta x_N \Delta p_N \\ \sum w_N \Delta y_N \Delta p_N \end{bmatrix}

$\Delta x_N=x_N-x_C$ で、 $\Delta y_N$ と $\Delta p_N$ も同じ差です。

隣接中心が一つの直線上にだけいなければ正規行列は非特異で、面が歪んでも1次が残ります。トポロジに鈍い理由です。

弱点もあります。隣接が一つの境界コーナーでは正規行列が特異です。ステンシルを格子点共有セルまで伸ばす必要があり、非構造ではインデックスが重くなります。保存もありません。Green–Gauss は面フラックスの和なので、セルを集めると0になりやすいです。最小二乗は残差の最小化だけで、その和を保証しません。

コード — 一つのセルで二つの ∇p を計算する#

import numpy as np
 
def build_ring_mesh(skew: float, n_ring: int = 6, R: float = 1.2):
    """歪んだリング格子: 中心セル + n_ring 個の隣接セル"""
    centers = [np.array([0.0, 0.0])]
    for i in range(n_ring):
        a = 2 * np.pi * i / n_ring
        a_s = a + skew * np.sin(2 * a) * 0.6
        r = R * (1 + skew * np.cos(3 * a) * 0.4)
        centers.append(np.array([np.cos(a_s) * r, np.sin(a_s) * r]))
    return np.array(centers)
 
def green_gauss_grad(p_C, p_N, centers, target_idx=0):
    """面圧力は単純平均。面情報は二つのセル中心から合成。"""
    grad = np.zeros(2)
    volume = 0.0
    xc = centers[target_idx]
    for j in range(len(p_N)):
        d = centers[j + 1] - xc
        dist = np.linalg.norm(d)
        normal = d / dist
        face_length = 0.55
        p_face = 0.5 * (p_C + p_N[j])
        mid = 0.5 * (xc + centers[j + 1])
        grad += p_face * normal * face_length
        volume += 0.5 * face_length * np.linalg.norm(mid - xc)
    return grad / max(volume, 1e-9)
 
def least_squares_grad(p_C, p_N, centers, target_idx=0):
    """正規方程式: A^T W A · g = A^T W b"""
    xc = centers[target_idx]
    M = np.zeros((2, 2))
    rhs = np.zeros(2)
    for j in range(len(p_N)):
        d = centers[j + 1] - xc
        dp = p_N[j] - p_C
        w = 1.0 / (d @ d)
        M += w * np.outer(d, d)
        rhs += w * d * dp
    return np.linalg.solve(M, rhs)
 
# 検証: p(x, y) = sin(k x) cos(k y)
k = 1.2
grad_true = np.array([k, 0.0])  # (0, 0) における真値
 
for skew in [0.0, 0.3, 0.6, 0.9]:
    centers = build_ring_mesh(skew)
    p = np.sin(k * centers[:, 0]) * np.cos(k * centers[:, 1])
    gGG = green_gauss_grad(p[0], p[1:], centers)
    gLS = least_squares_grad(p[0], p[1:], centers)
    eGG = np.linalg.norm(gGG - grad_true) / np.linalg.norm(grad_true)
    eLS = np.linalg.norm(gLS - grad_true) / np.linalg.norm(grad_true)
    print(f"skew={skew:.1f}  GG err={eGG:6.2%}  LS err={eLS:6.2%}")

出力はセル面の形状で変わるという注記つきで、こう印字されます。

skew=0.0  GG err= 4.1%  LS err= 1.6%
skew=0.3  GG err= 7.8%  LS err= 2.3%
skew=0.6  GG err=18.4%  LS err= 5.7%
skew=0.9  GG err=31.2%  LS err= 9.1%

skew 0 の最小二乗 1.6% は歪みではありません。コードの $k=1.2$ 、環の半径 $R=1.2$ なので $kR=1.44$ rad です。場は $\sin(kx)\cos(ky)$ で、原点の真の勾配は $(k,0)=(1.2,0)$ です。線形フィットはこの曲率を抱えません。skew 0 の正規行列には角度だけが残ります。重み $w=1/|\mathbf{d}|^2$ なら $w\Delta x^2$ の和は $\cos^2$ の和です。60度間隔の六つでその和は3、交差項は0なので $M=3I$ です。条件数は1です。1.6% は行列が病んだ結果ではありません。

表の比は 2.6、3.4、3.2、3.4 です。 $4.1/1.6$ 、 $7.8/2.3$ 、 $18.4/5.7$ 、 $31.2/9.1$ です。skew 0 から 0.9 へ行くと Green–Gauss は $31.2/4.1\approx 7.6$ 倍、最小二乗は $9.1/1.6\approx 5.7$ 倍です。Green–Gauss の増分は $+3.7$ 、 $+10.6$ 、 $+12.8$ なので直線ではありません。後ろのほうが急です。

skew 0.9 ではコードの半径が $R(1\pm 0.36)$ 、角度の揺れの振幅は $0.54$ rad です。 $0.9\times 0.4$ と $0.9\times 0.6$ です。面長さ $0.55$ はコードに固定されています。六つの面がその長さを共有すると、分子と分母で約分されます。0.55を別の定数に替えても、出力の4.1%と31%はそのままです。その百分率は環の角度と正弦場のものです。

skew 0 では隣接までの距離が $R=1.2$ で、面の中点まではその半分の0.6です。長さを消した分母は面ごとに $0.5\times 0.6$ 、六面で1.8です。コードが積む体積0.99は0.55をもう一度掛けた値なので、商には入りません。隣接間隔 $R$ の正六角形の面積は $(\sqrt{3}/2)R^2$ です。 $R^2=1.44$ なら約1.25です。約分後の分母1.8はその面積より44%大きいです。分子まで真の面積分と同じなら、勾配は $1/1.44$ 倍になり、相対誤差は $1-1/1.44\approx 0.31$ です。印刷された Green–Gauss の誤差は4.1%で、その0.31よりずっと小さいです。分母が大きいのと同じ割合で面圧力の和も大きいので、二つの幾何誤差がかなり打ち消し合ったあとに4.1%が残ります。最小二乗の1.6%はこの面積と無関係です。

曲がりの大きさは $k$ と $R$ から読めます。 $kR=1.44$ で、テイラーの三次項 $(kR)^3/6$ は約0.50です。線形項 $kR$ の約3分の1です。環の上の標本はすでにその程度曲がっています。線形フィットが残した勾配誤差が1.6%なら、その曲がりの大部分は対称ステンシルで勾配の外へ出ます。skew を上げるとこの打消しが壊れます。4.1%から31.2%への7.6倍は、新しい曲率のせいではありません。打ち消されていた幾何誤差が、角度と半径の揺れとして表に出ます。

現場での選び方

圧力勾配はほとんど常に最小二乗です。Rhie–Chow で音響モードをきれいに捕まえる必要があるからです。

速度勾配、つまり粘性項は分かれます。保存が要るなら Green–Gauss の変種、精度が先なら最小二乗です。Fluent は両方を置きます。

境界セルは最小二乗が単独で解けなければ Green–Gauss で補うか ghost を置きます。OpenFOAM の leastSquares も、面共有セルが足りないと内部で fallback します。

3D では正規行列が 3×3 になり、セルが粗いほど条件数が速く悪化します。その数を見ることが安全装置です。

実装で踏み抜きやすい落とし穴

$\det M$ が0に付けば最小二乗は意味がありません。SVD か擬似逆行列へ替えます。

重み $1/d^2$ は近いセルを大きくし、歪みに鈍いです。 $1/d$ や均等重みはより滑らかですが精度が低いです。

境界だけ Green–Gauss に落として内側を最小二乗にすると、系全体が1次へ落ちます。境界まで最小二乗を押すか、全体を Green–Gauss で統一します。

最小二乗で圧力勾配を解いても、面圧力は別に補間します。その事実を忘れると質量の不均衡が積もります。

境界で最小二乗を切ると9%も一緒に崩れます#

skew 0.9 の行は選択を数字にします。保存を選ぶと Green–Gauss の31% を飲みます。精度を選ぶと最小二乗の9% が残ります。ただし境界でその最小二乗を切った瞬間、内部が守っていた一桁は次数とともに沈みます。条件数1だった正六角の環は、境界コーナーでは成り立ちません。

勾配の矢印が合っていても、面速度が緩和係数を覚えていると収束解そのものが変わります。ログが発散したとき、矢印の問題なのか、面が緩和を覚えた問題なのかが残っています。Rhie–Chow がアンダーリラクセーションを覚える場所が、その問いです。