非構造格子で2次再構成が生む偽の極値 — Barth–Jespersen と Venkatakrishnan リミッター

無制限移流の最小は -0.0417 です。リミッターで床は 0 に戻りますが、その Φ が残差を止めることがあります。

問題なく回っていたソルバーです。2次を入れた瞬間、密度が負になりました。衝撃波の手前で、再構成した面の値がセル平均より上に出ていました。同じ格子の1次風上は壊れていません。精度を上げた側が計算を殺しています。非構造セルの中の直線が、隣接のどこにもない極値を作るからです。Barth–Jespersen はその極値を切ります。Venkatakrishnan は切った角を丸めます。

再構成がなかったはずの極値をつくる

有限体積法（FVM、セル平均を未知数とする離散化）で2次とは、セルの中の直線です。

u_f = u_i + \Phi_i\,\nabla u_i\cdot(\mathbf{x}_f-\mathbf{x}_i)

$u_i$ はセル $i$ の平均です。 $\nabla u_i$ は再構成した勾配、 $\mathbf{x}_f$ は面の中心、 $\Phi_i\in[0,1]$ はリミッターです。

$\Phi_i=1$ ではリミッターが切れています。勾配は最小二乗か Green–Gauss で、隣接を平均して作ります。不連続のそばでは、その平均が急すぎます。直線を面まで伸ばすと、どの隣接にもなかった最大・最小が生まれます。

勾配の gain を上げると、制限なし（Φ=1）の再構成が灰色の帯、つまり隣接の最大・最小を突き破って赤に変わります。

gradient gain: 1.40
Venkatakrishnan K: 0.30

Shaded band = allowed [min, max] from neighbors. Red segment = reconstruction escapes the band (new extremum). Φ=1 keeps full slope; Φ→0 flattens the cell to first order.

線が帯の中にあるセルでは、Φ は 1 のままです。ジャンプの隣だけ、線分が先に短くなります。同じ gain でリミッターを Barth–Jespersen に切り替えると、線は帯の外へ出ません。

偽の極値は次のステップのフラックスを汚します。振動の行き着く先は、負の密度と負のエネルギーです。再構成した面の値は、隣接がすでに持つ範囲を出てはいけません。 局所最大最小原理（LMP）です。 $\Phi_i$ は、その範囲を守るために勾配を削るスカラーです。

Barth–Jespersen：隣接の最大・最小で閉じ込める#

Barth と Jespersen（1989）は、許容範囲を先に固定します。セル $i$ と隣接が作る最大・最小です。

u_i^{\max}=\max\!\left(u_i,\max_{j\in N(i)}u_j\right),\quad u_i^{\min}=\min\!\left(u_i,\min_{j\in N(i)}u_j\right)

$N(i)$ はセル $i$ の面隣接の集合です。各面で、リミッターなしの増分 $u_f-u_i$ を見て、面ごとの係数 $\phi_f$ を求めます。

\phi_f= \begin{cases} \min\!\left(1,\dfrac{u_i^{\max}-u_i}{u_f-u_i}\right), & u_f-u_i>0\\[2mm] \min\!\left(1,\dfrac{u_i^{\min}-u_i}{u_f-u_i}\right), & u_f-u_i<0\\[1mm] 1, & u_f-u_i=0 \end{cases}

セルのリミッターは、全部の面のうち最も保守的な値です。

\Phi_i=\min_{f\in \text{faces}(i)}\phi_f

範囲の内側なら $\phi_f=1$ で、2次が残ります。はみ出そうとする面は、境界にちょうど触れる勾配だけを残します。極値のそばでは $\Phi_i\to 0$ になり、1次風上へ戻ります。

正の増分だけ読むと $\phi_f=\min(1,(u_i^{\max}-u_i)/(u_f-u_i))$ です。分子が分母より小さくなった瞬間、面の値は $u_i^{\max}$ に貼り付いて止まります。上の図で赤い線が灰色の境界で折れるのは、この比のせいです。

微分不可能という第二の罠

Barth–Jespersen は単調性を守ります。定常ソルバーに入れると、残差がある線から下がらず、そこで震えます。原因は $\min$ と割り算です。 $\phi_f$ は解に対して微分できません。どの面が最小を与えるかは、反復のたびに入れ替わります。リミッターが付いたり消えたりすると、残差が跳ね返ります。陰的（implicit）ソルバーのヤコビアンは、この不連続を嫌います。

滑らかな場所の Φ は 1 で足ります。本物の不連続でだけ起動し、その境界は滑らかであってほしいのです。鋭い折れを丸める式が次です。

Venkatakrishnan：滑らかに削る#

Venkatakrishnan（1993）は $\min(1,y)$ の折れを有理関数に置き換えました。ソースが触れている Michalak–Ollivier-Gooch リミッターも同じ系統です。面ごとの係数は次のように書きます。

\phi_f=\frac{1}{\Delta_-}\, \frac{(\Delta_+^2+\epsilon^2)\,\Delta_- + 2\,\Delta_-^2\,\Delta_+} {\Delta_+^2 + 2\,\Delta_-^2 + \Delta_-\Delta_+ + \epsilon^2}

$\Delta_-=u_f-u_i$ はリミッターなしの増分です。 $\Delta_+$ は増分の符号に応じて $u_i^{\max}-u_i$ か $u_i^{\min}-u_i$ です。手で選ぶ量は $\epsilon^2$ です。

\epsilon^2=(K\,\Delta x)^3

$K$ は調整パラメータ、 $\Delta x$ は格子サイズです。変動が $\epsilon$ より小さいと $\phi_f\to 1$ になり、リミッターは切れます。 $K$ を上げると閾値が上がり、リミッターが緩む領域が広がります。上げすぎると、衝撃波のそばの振動を見逃します。 $K$ を 0 に送ると Barth–Jespersen に戻ります。

コードはセルサイズを 1 として $\epsilon^2=K^3$ を計算します。既定の $K=0.3$ では $\epsilon^2=0.027$ です。 $\epsilon=0$ かつ $\Delta_+=\Delta_-=\Delta$ のとき、再構成が許容限界にちょうど触れる場所で、有理関数は $(3\Delta^2)/(4\Delta^2)=3/4$ になります。同じ場所で Barth–Jespersen は $\phi_f=1$ です。Venkatakrishnan は限界に触れる前に、勾配の 4 分の 1 を捨てます。角が一段早く鈍る理由です。

現場でよく試す範囲は $K\approx 0.1$ から $K\approx 5$ です。コードと同じ $\epsilon^2=K^3$ にこの両端を入れると、0.001 と 125 が出ます。閾値の比は 12 万 5 千です。 $K$ を少し動かすだけでも三乗なので、滑らかだと免除されるセルが急に増えます。逆に $K$ を 0.1 より下に置くと $\epsilon^2$ は実質 0 で、Barth–Jespersen の折れが戻ってきます。

上の図で Venkatakrishnan をオンにし、 $K$ を動かすと、Φ は段差なしで滑ります。

Python — スカラー移流で3つのリミッターを競わせる#

周期境界です。矩形波とガウス山を一緒に移流（advection）させます。無制限（Fromm）、Barth–Jespersen、Venkatakrishnan を並べ、最後の最小値を見ます。0 を下回れば、その深さだけ偽の谷が掘られています。

セル勾配は $s=(u_{i+1}-u_{i-1})/2$ です。面へ出す増分は $0.5s$ なので、左右の隣接の差が 1 なら $s=1/2$ 、面の増分は 0.25 です。その 0.25 が許容幅より大きいと、Φ は 1 を下回ります。1次元の中心差分は左右の隣接しか見ません。非構造セルで同じ勾配の重みを決めるのは、移動最小二乗がセルごとに合わせ直すやり方とは別です。リミッターは、勾配ができた後に付きます。

import numpy as np
 
NX, A, CFL = 200, 1.0, 0.4
dx = 1.0 / NX
dt = CFL * dx / A
 
def init_profile():
    x = (np.arange(NX) + 0.5) / NX
    u = np.where((x > 0.1) & (x < 0.3), 1.0, 0.0)   # 矩形波
    u += np.exp(-((x - 0.65) / 0.06) ** 2)          # ガウス山
    return u
 
def cell_slope(u):
    return (np.roll(u, -1) - np.roll(u, 1)) / 2.0    # 中心勾配（Fromm）
 
def barth_jespersen_phi(u, s):
    up, um = np.roll(u, -1), np.roll(u, 1)
    umax = np.maximum(u, np.maximum(up, um))
    umin = np.minimum(u, np.minimum(up, um))
    phi = np.ones_like(u)
    for du in (0.5 * s, -0.5 * s):                   # 2つの面
        f = np.ones_like(u)
        pos, neg = du > 1e-12, du < -1e-12
        f[pos] = np.minimum(1.0, (umax[pos] - u[pos]) / du[pos])
        f[neg] = np.minimum(1.0, (umin[neg] - u[neg]) / du[neg])
        phi = np.minimum(phi, f)
    return np.clip(phi, 0.0, 1.0)
 
def venkatakrishnan_phi(u, s, K=0.3):
    up, um = np.roll(u, -1), np.roll(u, 1)
    umax = np.maximum(u, np.maximum(up, um))
    umin = np.minimum(u, np.minimum(up, um))
    eps2 = K ** 3                                    # (K*h)^3、セルサイズ h=1 単位
    phi = np.ones_like(u)
    for du in (0.5 * s, -0.5 * s):
        d = np.where(du > 0, umax - u, umin - u)
        num = (d * d + eps2) * du + 2 * du * du * d
        den = d * d + 2 * du * du + d * du + eps2
        f = np.where(np.abs(du) < 1e-12, 1.0, num / (du * den))
        phi = np.minimum(phi, f)
    return np.clip(phi, 0.0, 1.0)
 
def muscl_rhs(u, limiter):
    s = cell_slope(u)
    phi = limiter(u, s) if limiter else np.ones_like(u)
    uL = u + 0.5 * phi * s              # 面 i+1/2 左側の状態
    flux = A * uL                       # 風上フラックス (a > 0)
    return -(flux - np.roll(flux, 1)) / dx
 
def advance_muscl(u, limiter):         # SSP-RK2
    k1 = muscl_rhs(u, limiter)
    u1 = u + dt * k1
    k2 = muscl_rhs(u1, limiter)
    return 0.5 * (u + u1 + dt * k2)
 
def run_limiter_race(steps=160):
    fields = {"none": init_profile(), "bj": init_profile(), "venk": init_profile()}
    lims = {"none": None, "bj": barth_jespersen_phi, "venk": venkatakrishnan_phi}
    for _ in range(steps):
        for k in fields:
            fields[k] = advance_muscl(fields[k], lims[k])
    for k, u in fields.items():
        print(f"{k:5s}  min={u.min():+.4f}  max={u.max():.4f}")
 
run_limiter_race()

出力は次のとおりです。

none   min=-0.0417  max=1.0231
bj     min=+0.0000  max=1.0000
venk   min=-0.0004  max=1.0009

Unlimited (Fromm)Barth–JespersenVenkatakrishnan

steps/frame: 3t = 0.00

Watch the trailing edge of the square wave: the unlimited scheme grows ripples below zero, while both limiters stay monotone.

無制限の最小 -0.0417 は、高さ 1 の矩形波が 0 の下に掘れた深さです。最大 1.0231 は、同じ高さから上にはみ出した量です。Barth–Jespersen は最小・最大を $[0,1]$ に固定しました。Venkatakrishnan の最小 -0.0004 は、その谷の約 104 分の 1 です。1 を超える幅は 0.0231 から 0.0009 に縮みます。

160 ステップの移動距離は $A \cdot dt \cdot 160$ です。 $dx=1/200$ 、CFL が 0.4 なら $dt=0.002$ なので、距離は 0.32 です。200 セルのうち 64 セル分です。矩形波の区間 $(0.1, 0.3)$ は $(0.42, 0.62)$ へ動き、ガウスの中心 0.65 は 0.97 に着きます。両方とも同じ速度なので、相対位置はそのままです。出力の最小・最大は、二つの波がぶつかった結果ではありません。それぞれの角で別々に生まれた値です。矩形波の長さ 0.2 は 40 セル、ガウスの式のスケール 0.06 は 12 セルです。最小・最大を分けたのは、その角の近くの数セルです。

矩形波の後ろ側の角だけ見れば足ります。無制限（赤）は 0 の下にさざ波を育てます。二つのリミッターは単調なままです。Venkatakrishnan（黄）は Barth–Jespersen（シアン）より、角が少し鈍くなります。 $\epsilon=0$ で勾配の 3/4 だけが残った計算と同じ向きです。コードの $K=0.3$ はその折れを消さず、超過だけを 0.0009 の水準まで下げます。

現場でリミッターをオンにするとき

隣接をどう束ねるかで、許容範囲が分かれます。Barth–Jespersen の $\max/\min$ を面隣接で取るか、頂点（vertex）隣接で取るかで結果が変わります。非構造格子で頂点隣接、つまり頂点を共有するすべてのセルを使うと、方向の偏りが減ってバランスが良くなります。

次は $K$ です。小さすぎると Barth–Jespersen のように収束が停滞します。大きすぎると衝撃波で振動を見逃します。 $\epsilon^2$ に $\Delta x^3$ が入るので、格子を細かくすると同じ $K$ でもリミッターはより頻繁に起動します。

ベクトル系（Euler 方程式）では、保存変数ごとに別々にリミッターをかけると、成分のあいだのずれから新しい振動が生まれます。特性変数（characteristic variable）へ射影してからかける方が安全です。費用は増えます。

残差が止まったら、勾配を先に疑わない

2次ソルバーが負の密度で死んだら、面の値が隣接範囲の外に出たかを先に見ます。外なら、悪いのは勾配ルーチンではなく、Φ を作った隣接集合です。面隣接だけでは、対角方向の極値を落とします。再構成はその隙間からまた漏れます。

定常で残差だけが足踏みするなら、 $K$ を開きます。 $K=0$ は Barth–Jespersen の折れなので、ヤコビアンが反復ごとに飛びます。ジャンプで Φ が 0 まで下がると、そのセルの次数は 1 次に張り付きます。滑らかな領域は 2 次の近くに残るでしょうか。残るなら、 $\epsilon^2$ がその領域を放しています。この分岐を一つの計算で同時に測ると何が出るかは、衝撃波と滑らかな渦を並べて次数を測った記録の側にまだ開いています。