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

정지 별의 압력 source는 2P/r·V가 아니라 면적 차 P·ΔS여야 중심이 비지 않는다.

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

Δ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})$ 로 줄어도 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 근사 값이다. 둘은 정확히 같지 않다. 잔차가 남고, 그 잔차가 매 스텝 가속이 된다.

출력 숫자로 돌아가자. $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}$ 로 줄지만 0은 아니다. 앞 그림의 잔차는 $O(\Delta r^{3})$ 다. 한 차수 차이. 한 셀의 source 오차는 면적 차의 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$ 과 같은 구멍이다. 이산 발산과 면적이 한 쌍이 아니면 정지 해는 버티지 못한다.

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

해법은 단순하다. 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

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}$ 이다.

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

기하 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. 그러나 별이 다시 정지한다.

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$ 방향은 두 면의 면적이 같아서 이 줄이 필요 없다.

회전이 세면 압력 source로 끝이 아니다. $\rho u w/r + \cos\theta\,\rho v w/r$ 가 $\varphi$ -운동량에 남는다. $l = w r \sin\theta$ 로 바꾸면 그 source가 식에서 사라진다. Kley 1998(A&A 338, L37)의 형태다. 회전이 반경마다 다르면 FARGO가 더 붙는다. 정지한 계면이 스스로 흐르는 병과는 기생 전류가 같다. 이산 발산이 만든 가짜 힘을 같은 이산식으로 빼지 않으면 평형은 자리로 돌아오지 않는다.