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

静止した星を保つには、圧力ソースを教科書の 2P/r·V ではなく面積差 P·ΔS で書く必要があります。

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

$\Delta r$ のスライダを上げると、教科書の $2P/r \cdot V$ と面積差 $P(S_{\text{out}} - S_{\text{in}})$ が離れます。

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.

残差は $O(\Delta r^{3})$ まで小さくなってもゼロにはなりません。格子が粗いほど、その隙間が先に見えます。

犯人は格子そのものです。曲線座標系は運動量方程式に余分な項を加えます。その一行を離散化し損なうと静止平衡が崩れます。φ-運動量を角運動量に置き換える操作も、その隣にあります。

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

極座標 $(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 近似の値で、両者は厳密には一致しません。差分が残り、その差分が毎ステップ加速として積み上がります。

出力の数字に戻ります。 $N=40$ 、半径は 1 から 4、ステップは 200 です。naive の $\max|u|$ は $1.834\times 10^{-3}$ です。balanced は $0.000\times 10^{0}$ です。初期状態は $P=1$ 、 $\rho=1$ 、 $u=0$ 、 $\gamma=1.4$ です。 $\Delta r=(4-1)/40=0.075$ です。時間刻みはコードどおり $\Delta t=0.4\,\Delta r/c_0$ で、 $c_0=\sqrt{\gamma P/\rho}$ です。その刻みを 200 回積むと、偽の速度が $1.834\times 10^{-3}$ まで届きます。バランス形式は同じループで動きません。

セルを 200 に増やすと速度残差は $(\Delta r)^{2}$ で減りますが、ゼロにはなりません。先の図の残差は $O(\Delta r^{3})$ です。一桁違います。一セルのソース誤差は面積差の Taylor 剰余なので、体積に対してはより小さくなります。それを $\rho V$ で割ると加速度になり、 $\Delta t$ を繰り返すと速度になります。速度側の次数が一段低く見える理由はここです。 $P(S_{i+1/2}-S_{i-1/2})$ に置き換えると次数の話は終わります。左辺のフラックス差分と右辺が同じ式だからです。

体積は $V=(r_R^{3}-r_L^{3})/3$ 、面積は $S=r^{2}$ です。 $\sin\theta\,\Delta\theta\,\Delta\varphi$ はコードの注釈どおり両側から外しています。係数が同じなら相殺は残ります。面積の定義だけを変えると $P\cdot\Delta S$ がフラックス差分とずれます。曲線セルの真の面積を外積で入れると、 $\partial x/\partial x$ が 1 から外れます。極座標の $2P/r$ と同じ穴です。離散発散と面積が対でなければ、静止解は持ちこたえません。

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

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

\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

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 を増やすと、曲線は 0 から離れます。well-balanced に切り替えると 0 軸に張り付きます。cells を減らすと naive の残差だけが大きくなります。40 セル、200 ステップの $\max|u|=1.834\times 10^{-3}$ がそのずれで、balanced は $0.000\times 10^{0}$ です。

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

幾何学的ソースは圧力だけではありません。 $\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})

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

2P/r を残すか、ΔS に換えるか#

圧力だけが揺れるなら、 $2P/r\,V$ を $P(S_{i+1/2}-S_{i-1/2})$ に換えます。40 セルで $1.834\times 10^{-3}$ だった $\max|u|$ が 0 になる箇所です。 $\theta$ 方向の $\cos\theta\,P/(r\sin\theta)$ も $P\cdot\Delta S_{\theta}$ への同じ交換です。 $\varphi$ 方向は両面の面積が等しいので、この一行は要りません。

回転が強いときは圧力ソースだけでは終わりません。 $\rho u w/r + \cos\theta\,\rho v w/r$ が $\varphi$ -運動量に残ります。 $l = w r \sin\theta$ に変えると、そのソースは式から消えます。Kley 1998 (A&A 338, L37) の形です。回転が半径ごとに違うなら FARGO がさらに要ります。止まった界面が勝手に流れ出す症状は、寄生電流と同じ病気です。離散発散が作った偽の力を、同じ離散式で引かなければ、平衡は元の位置に戻りません。