정지된 별이 움직이기 시작했다 — 극좌표 FVM의 well-balanced와 각운동량 보존

같은 식, 같은 슬로프 리미터, 같은 Riemann solver. 격자만 cartesian에서 polar로 옮겼다. 그러자 정지된 별의 중심에서 가스가 새기 시작했다. 압력은 균일했고, 속도는 0이었다. 무엇이 살아 움직였는가.

대답은 격자 자체에 있다. 곡선 좌표계가 운동량 방정식에 덧붙이는 항들이 있고, 그것들을 이산화할 때 한 줄을 잘못 쓰면 정지 평형이 깨진다. 이 글은 그 한 줄을 추적한다. 부수적으로 φ-운동량을 각운동량으로 옮기는 트릭까지 본다.

극좌표가 운동량에 끼얹은 항들#

극좌표 $(r, \theta, \varphi)$ 에서 Euler 방정식의 반경 방향 운동량은 다음과 같다.

\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$ 을 보존형 발산 안으로 넣으면 식 (8.13)이 된다.

\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 solver가 그대로 돌려준 플럭스다. 우변 $2P/r \cdot V$ 항이 오늘의 주인공이다. 이 항이 어디서 왔는지 보자.

$\partial P / \partial r$ 을 $\frac{1}{r^{2}} \partial (P r^{2})/\partial r - 2P/r$ 로 다시 쓰면, 앞 부분이 자연스럽게 플럭스로 흡수되고 뒤의 $-2P/r$ 이 좌변으로 넘어가 source가 된다. 이산화하면

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

가장 자연스러운 표기다. 그런데 여기에 함정이 있다. 압력이 어디에서나 일정하고 속도가 0인 정지 해를 코드에 넣으면, $\Delta_{r}(F_{rr} S_{r}) = P (S_{i+1/2} - S_{i-1/2})$ 가 0이 아니다. 왜냐하면 안쪽 면과 바깥쪽 면의 표면적이 다르기 때문이다. 그것을 상쇄해야 할 우변의 $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})$ 로 작아지지만 0이 아니다.

한 줄을 바꿔 평형을 되돌리다#

해법은 단순하다. source 항을 다음으로 교체한다.

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

이제 정지 해에서 좌변 플럭스 차분과 우변 source가 비트 단위로 정확히 상쇄된다. 별의 중심에서 가스가 새지 않는다.

같은 트릭이 $\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$ -운동량에는 이런 기하 압력 source가 처음부터 없다 — 방위 방향 면적은 두 면이 같기 때문이다.

이 처방은 well-balanced 기법의 가장 단순한 예다. "디스크리트 발산이 만든 가짜 가속을 같은 디스크리트 방식으로 정확히 빼 준다"는 한 줄에 요약된다.

NumPy 60줄 — 정지 평형이 깨지나, 살아남나#

1D 반경 격자에서 $P = 1$ , $\rho = 1$ , $u = 0$ 로 시작해 두 방식의 운동량 잔차를 비교한다.

import numpy as np
 
def polar_cell_volume(rL: float, rR: float) -> float:
    """spherical-shell volume between rL and rR (sin θ Δθ Δφ omitted)."""
    return (rR**3 - rL**3) / 3.0
 
def polar_r_surface(r: float) -> float:
    """r²-weighted face area in spherical coords."""
    return r * r
 
def naive_pressure_source(P: float, r_mid: float, V: float) -> float:
    """The textbook 2P/r · V form — wrong at machine precision for P=const."""
    return 2.0 * P / r_mid * V
 
def balanced_pressure_source(P: float, S_in: float, S_out: float) -> float:
    """Same expression as the flux difference, so they cancel exactly."""
    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):
        # pressure flux through faces (P constant → exact value)
        flux_diff = P * (S[1:] - S[:-1])
        # source term (the line we want to scrutinise)
        if scheme == "naive":
            src = naive_pressure_source(P, r_mid, V)
        else:
            src = balanced_pressure_source(P, S[:-1], S[1:])
        # momentum update: d(ρ u V)/dt = -flux_diff + src
        # for static P, flux_diff is the only nontrivial term on LHS
        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}$ 로 줄지만 0은 되지 않는다. 균형 잡힌 형태는 처음부터 정확히 0이다.

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을 늘려 보자. 처음 100스텝 동안은 미세하지만 시간이 갈수록 누적된다. well-balanced로 토글하면 곡선이 0선에 붙어 떠나지 않는다. cells를 줄이면 격자가 거칠어지면서 naive의 잔차가 눈에 띄게 커진다.

φ-운동량을 각운동량으로#

기하 source는 압력에만 붙는 것이 아니다. $\varphi$ -운동량 방정식에는 $\rho u w / r + \cos\theta \rho v w / r$ 같은 항이 또 있다. 별이나 행성 disk처럼 회전이 강한 시스템에서 이 항들은 정확도와 안정성에 모두 해롭다.

해결은 변수를 바꾸는 것이다. $l = w r \sin\theta$ 를 정의하면 — 각운동량 그대로다 — $\varphi$ -운동량 식의 source가 전부 사라진다.

\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 solver가 돌려준 $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})$ 잔차	✅ 비트-단위	✅ (φ-방향만)
추가 코드 변경	없음	source 한 줄 교체	flux 정의 + 변수 추가
회전 시스템	⚠ 회전 누적 오차	△ 압력 source만 해결	✅ 본질적 보존
디버깅 비용	높음 (조용히 망가짐)	낮음	중간
어디서 만난다	"왜 별 중심이 비어 가는가?"	spherical hydro 표준	accretion disk

다시 별 앞에 서서#

격자가 직선이라면 보지 못했을 함정이다. 곡선 좌표계는 운동량 식에 항을 더 붙이고, 그것의 이산 표현을 플럭스 차분과 같은 형태로 적지 않으면 가만히 있던 가스가 흐른다. 이 글의 한 줄은 그래서

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

다. 한 줄, 한 commit. 그러나 별이 다시 정지한다.

노트북에 옮겨 적어 둘 것#

정지 해가 흔들리면 source 항 이산화부터 의심한다 — 보통 표면적 차이를 잘못 쓴 것이다.
"well-balanced"는 충격관 같은 비-정지 해보다 정지 해의 머신 정밀도 보존이 핵심이다.
회전이 강한 문제(별, accretion disk, 토네이도)는 $\varphi$ -운동량 대신 각운동량을 푼다.