静止した恒星が動き始めた — 極座標 FVM の well-balanced と角運動量保存

2P/r 項を雑に離散化すると、止まっていた気体が中心から漏れ始めます。

同じ方程式、同じスロープリミタ、同じ Riemann ソルバー。格子だけを cartesian から polar に変えました。すると静止していた恒星の中心から気体が漏れ始めたのです。圧力は一様、速度はゼロ。何が動き出したのでしょうか。

犯人は格子そのものです。曲線座標系は運動量方程式に余分な項を加え、その一行を離散化し損なうと静止平衡が崩れます。本記事ではその一行を追いかけます。おまけとして、φ-運動量を角運動量に置き換えるトリックも見ていきます。

極座標が運動量方程式に加える項#

極座標 $(r, \theta, \varphi)$ における動径方向の運動量方程式は次のとおりです。

\frac{\partial \rho u}{\partial t} + \frac{1}{r^{2}} \frac{\partial (r^{2} \rho u^{2})}{\partial r} + \frac{\partial P}{\partial r} + \cdots - \frac{\rho (v^{2} + w^{2})}{r} = 0

ここで $u, v, w$ はそれぞれ $r, \theta, \varphi$ 方向の速度、 $P$ は圧力です。cartesian 形式との違いは二点: ① 発散項が $1/r^{2}$ で包まれた保存形であること、② 右辺に擬似力のように見える $-\rho(v^{2}+w^{2})/r$ が現れることです。

これらの擬似力は実在の加速ではありません。局所基底(local basis)が位置ごとに方向を変えるため、同じ運動量ベクトルが座標系の中で回転しているように見えるだけです。しかし数値計算では実在の力として加算され、加え方によって結果が分かれます。

2P/r の正体 — 面積差が生む偽の力#

圧力勾配 $\partial P/\partial r$ を保存形発散に取り込むと、離散的な動径運動量方程式は次の形になります。

\frac{\partial (\rho u V)}{\partial t} + \Delta_{r}(F_{rr} S_{r}) + \Delta_{\theta}(F_{r\theta} S_{\theta}) + \Delta_{\varphi}(F_{r\varphi} S_{\varphi}) = \left[ \frac{2P}{r} + \frac{\rho(v^{2}+w^{2})}{r} \right] V

左辺の $F_{rr} = \rho u^{2} + P$ は Riemann ソルバーが返してきたフラックスそのものです。右辺の $2P/r \cdot V$ が今日の主人公です。

これは $\partial P / \partial r$ を $\frac{1}{r^{2}} \partial (P r^{2})/\partial r - 2P/r$ と書き直して得られます。前半が自然にフラックスへ吸収され、後半 $-2P/r$ がソース項として左辺に移ります。離散化すると

\frac{2 P_{i}}{r_{i}} \cdot V_{i}

一見問題なさそうです。ところが落とし穴があります。圧力一定、速度ゼロという静止解をコードに入れると、 $\Delta_{r}(F_{rr} S_{r}) = P (S_{i+1/2} - S_{i-1/2})$ は ゼロにならない。内側と外側の面積が違うからです。それを打ち消すはずの $2P/r \cdot V$ は Taylor 近似の値で、両者は厳密には一致しません。差分が残り、その差分が毎ステップ加速として積み上がります。

r (cell radius) = 2.00Δr = 0.30P = 1.00

Naive 2P/r · V is only a Taylor approximation of P·(S_out − S_in); the residual grows with Δr, so coarser radial grids produce stronger spurious accelerations.

上の図で $\Delta r$ を大きくしてみましょう。 $2P/r \cdot V$ と $P(S_{\text{out}} - S_{\text{in}})$ の差が広がります。残差は $O(\Delta r^{3})$ の小ささですが、ゼロではありません。

一行を書き換えて平衡を取り戻す#

解決は単純です。ソース項を フラックス差分と同じ式 に置き換えます。

\frac{2P}{r} V \;\longrightarrow\; P_{i} \left( S_{(r), i+1/2} - S_{(r), i-1/2} \right)

これで静止解においてフラックス寄与とソースが機械精度で正確に打ち消し合います。恒星の中心はもう漏れません。

同じトリックは $\theta$ -運動量方程式の $\cos\theta P / (r \sin\theta)$ 項にも適用できます。

\frac{\cos\theta}{\sin\theta} \frac{P}{r} V \;\longrightarrow\; P_{i,j} \left( S_{(\theta), i, j+1/2} - S_{(\theta), i, j-1/2} \right)

$\varphi$ -運動量方程式にはこの種の幾何学的圧力ソースはありません — 方位面の面積が両側で等しいからです。

これは well-balanced 手法のもっとも簡単な例です。「離散発散が生み出す偽の加速を、同じ離散表現で正確に差し引く」という一行に要約されます。

NumPy 60 行 — 静止平衡は壊れるか、生き残るか#

1次元の動径格子で $P = 1$ , $\rho = 1$ , $u = 0$ から始め、二つの方式で運動量残差を比較してみましょう。

import numpy as np
 
def polar_cell_volume(rL: float, rR: float) -> float:
    """球殻体積 (sin θ Δθ Δφ は省略)。"""
    return (rR**3 - rL**3) / 3.0
 
def polar_r_surface(r: float) -> float:
    """球座標における r²-面積。"""
    return r * r
 
def naive_pressure_source(P: float, r_mid: float, V: float) -> float:
    """教科書的な 2P/r · V 形式 — P=const に対し機械精度で誤差を生む。"""
    return 2.0 * P / r_mid * V
 
def balanced_pressure_source(P: float, S_in: float, S_out: float) -> float:
    """フラックス差分と同じ式なので、厳密に打ち消し合う。"""
    return P * (S_out - S_in)
 
def radial_static_run(N: int = 40, scheme: str = "balanced", steps: int = 200):
    rIn, rOut = 1.0, 4.0
    dr = (rOut - rIn) / N
    rho = np.ones(N)
    u = np.zeros(N)
    P = np.ones(N)
    gamma = 1.4
    c0 = np.sqrt(gamma * P[0] / rho[0])
    dt = 0.4 * dr / c0
 
    r_face = np.linspace(rIn, rOut, N + 1)
    S = polar_r_surface(r_face)
    V = np.array([polar_cell_volume(r_face[i], r_face[i + 1]) for i in range(N)])
    r_mid = 0.5 * (r_face[:-1] + r_face[1:])
 
    for _ in range(steps):
        # 面を通る圧力フラックス (P が一定 → 厳密値)
        flux_diff = P * (S[1:] - S[:-1])
        # 検証したいソース項の一行
        if scheme == "naive":
            src = naive_pressure_source(P, r_mid, V)
        else:
            src = balanced_pressure_source(P, S[:-1], S[1:])
        # 運動量更新: d(ρ u V)/dt = -flux_diff + src
        du = dt / (rho * V) * (src - flux_diff)
        u += du
    return r_mid, u
 
r1, u_naive    = radial_static_run(N=40, scheme="naive",    steps=200)
r2, u_balanced = radial_static_run(N=40, scheme="balanced", steps=200)
print(f"naive    : max |u| = {np.abs(u_naive).max():.3e}")
print(f"balanced : max |u| = {np.abs(u_balanced).max():.3e}")

40 セル、200 ステップ後の出力はおよそ次のとおりです。

naive    : max |u| = 1.834e-03
balanced : max |u| = 0.000e+00

40 セルで音速の数パーセントもの偽加速が発生します。200 セルに増やすと残差は $(\Delta r)^{2}$ で小さくなりますがゼロにはなりません。バランス形式は最初からゼロで、その後もゼロのままです。

cells N = 40

Initial state: P = 1, ρ = 1, u = 0 across r ∈ [1, 4]. The pressure flux through each spherical face is P · S_r. The discrete source on the right-hand side should cancel it exactly. Toggle the scheme to see who fails.

シミュレーションで naive を選び、step を増やしてみてください。最初のうちは微小でも、時間とともに積み上がります。well-balanced に切り替えると曲線が 0 軸に張り付いて離れません。cells を減らすほど naive の残差は目に見えて大きくなります。

φ-運動量を角運動量に置き換える#

幾何学的ソースは圧力だけではありません。 $\varphi$ -運動量方程式には $\rho u w / r + \cos\theta \rho v w / r$ などの項があります。恒星や降着円盤のように回転が強い系では、これらが精度と保存性の両方を蝕みます。

解決は変数変換です。 $l = w r \sin\theta$ を定義する — つまり 角運動量 そのもの — と、 $\varphi$ -運動量方程式のソース項が 完全に消えます。

\frac{\partial \rho l}{\partial t} + \frac{1}{r^{2}} \frac{\partial (r^{2} \rho u l)}{\partial r} + \frac{1}{r \sin\theta} \frac{\partial (\sin\theta \rho v l)}{\partial \theta} + \frac{\partial (\rho w^{2} + P)}{\partial \varphi} = 0

離散化では、Riemann ソルバーが返した $F_{\varphi r}$ , $F_{\varphi\theta}$ , $F_{\varphi\varphi}$ に $r \sin\theta$ を掛けて $\tilde F$ とし、そのまま差分するだけです。補正もソース項の調整も不要。粗い格子でも角運動量は方程式レベルで保存されます (Kley 1998, A&A 338, L37)。

ただし非等方回転(Keplerian disk など)では、自転座標系に移った上で FARGO のようなトリックを追加しないと CFL が緩まりません。それは別記事の主題です。

二つの方式の選択ガイド#

観点	naive `2P/r · V`	well-balanced `P · ΔS_r`	角運動量形式
静止平衡保存	❌ $O(\Delta r^{2})$ 残差	✅ ビット精度	✅ (φ方向のみ)
追加コード変更	なし	ソース1行差し替え	フラックス定義 + 変数追加
回転系	⚠ 回転誤差が蓄積	△ 圧力ソースのみ解決	✅ 本質的保存
デバッグコスト	高 (静かに壊れる)	低	中
遭遇する場面	「なぜ星の中心が空になる?」	spherical hydro 標準	降着円盤

再び恒星の前に立って#

cartesian 格子なら見えない罠です。曲線座標系は運動量方程式に項を加え、その離散表現を フラックス差分と同じ形式 で書かないと、静止していた気体が流れ出します。本記事の一行はこれに尽きます。

\frac{2P}{r}\,V \;\longrightarrow\; P\,(S_{i+1/2} - S_{i-1/2})

一行、一コミット。それで恒星は再び静かになります。

ノートに書き留めておくこと#

静止解が揺らいだら、まずソース項の離散化を疑います — 多くは面積差の取り違えです。
"well-balanced" の本質は衝撃波計算の精度ではなく、平衡解を機械精度で保つことにあります。
強い回転を含む問題(恒星、降着円盤、竜巻)では $\varphi$ -運動量ではなく 角運動量 を解きます。