La estrella en reposo que empezó a moverse — FVM bien balanceado y conservación del momento angular en coordenadas polares

Una estrella polar en reposo sigue en reposo solo si la fuente de presión es el salto de área, no el 2P/r·V del manual.

Mismas ecuaciones, mismo limitador de pendiente, mismo solver de Riemann. Solo cambió la malla, de cartesiana a polar. El gas empezó a escapar del centro de una estrella en reposo. La presión era uniforme. Las velocidades, nulas.

Al subir $\Delta r$ , la fuente de manual $2P/r \cdot V$ se separa del salto de área $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.

El residuo baja como $O(\Delta r^{3})$ y no llega a cero. Una malla más gruesa abre esa brecha antes.

La culpa es de la malla. Las coordenadas curvilíneas añaden términos geométricos al momento, y una sola línea mal discretizada rompe el equilibrio estático. El cambio de momento-φ a momento angular queda al lado de esa línea.

Lo que las coordenadas polares añaden al momento#

En coordenadas polares $(r, \theta, \varphi)$ la ecuación de momento radial es

\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

donde $u, v, w$ son las velocidades en las direcciones $r, \theta, \varphi$ y $P$ la presión. Dos diferencias con la forma cartesiana: ① la divergencia queda envuelta por $1/r^{2}$ en forma conservativa; ② aparece una fuerza ficticia $-\rho(v^{2}+w^{2})/r$ en el lado derecho.

Estas pseudo-fuerzas no son aceleraciones reales. La base local cambia de dirección en cada punto, así que el mismo vector momento parece rotar cuando se expresa en esa base. Numéricamente, sin embargo, se suman a la actualización de la celda como si fueran fuerzas, y la manera de sumarlas decide el destino del cálculo.

El origen del 2P/r — diferencias de área disfrazadas de fuerza#

Si metemos $\partial P/\partial r$ dentro de la divergencia conservativa, la ecuación discreta de momento radial queda

\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

El flujo $F_{rr} = \rho u^{2} + P$ es exactamente lo que devuelve el solver de Riemann. El término $2P/r \cdot V$ del lado derecho es el protagonista de hoy.

Nace de reescribir $\partial P / \partial r$ como $\frac{1}{r^{2}} \partial (P r^{2})/\partial r - 2P/r$ . La primera parte se integra de forma natural en la divergencia de flujo; la segunda, $-2P/r$ , pasa al lado derecho como fuente. Discretizada celda a celda da

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

Parece inofensivo. Hasta que se introduce una solución estática — presión constante, velocidades cero — y se revisan las cuentas. El lado de los flujos arroja $\Delta_{r}(F_{rr} S_{r}) = P (S_{i+1/2} - S_{i-1/2})$ , que no es cero porque el área de la cara interna y la externa de una celda polar son distintas. La fuente, derivada por Taylor como $2P/r \cdot V$ , solo iguala esa diferencia de manera aproximada. La discrepancia se queda dentro de la celda y produce una aceleración espuria en cada paso.

De vuelta a los números impresos. $N=40$ , radio de 1 a 4, 200 pasos. El $\max|u|$ naive es $1.834\times 10^{-3}$ . La forma balanceada imprime $0.000\times 10^{0}$ . El estado inicial es $P=1$ , $\rho=1$ , $u=0$ , $\gamma=1.4$ . $\Delta r=(4-1)/40=0.075$ . El paso temporal del código es $\Delta t=0.4\,\Delta r/c_0$ con $c_0=\sqrt{\gamma P/\rho}$ . Doscientos de esos pasos acumulan una velocidad espuria de $1.834\times 10^{-3}$ . La forma balanceada no se mueve en el mismo bucle.

Al refinar a 200 celdas el residuo de velocidad cae como $(\Delta r)^{2}$ , no hasta cero. La figura cita un hueco $O(\Delta r^{3})$ . Un orden de diferencia. El error de la fuente en una celda es el resto de Taylor del desajuste de áreas, más pequeño frente al volumen. Dividido por $\rho V$ es una aceleración. Repetido con $\Delta t$ pasa a ser velocidad. Por eso el residuo de velocidad parece un orden más blando. Al sustituir $P(S_{i+1/2}-S_{i-1/2})$ se acaba la discusión de órdenes: la diferencia de flujos y la fuente son la misma expresión.

El volumen es $V=(r_R^{3}-r_L^{3})/3$ y el área $S=r^{2}$ . El factor $\sin\theta\,\Delta\theta\,\Delta\varphi$ se omite en ambos lados, como dice el comentario del código. Coeficientes iguales conservan la cancelación. Cambiar solo la definición de área hace que $P\cdot\Delta S$ deje de coincidir con la diferencia de flujos. El mismo agujero aparece cuando el área real de una celda curva sale de un producto cruzado de diagonales: $\partial x/\partial x$ se aparta de 1. La divergencia discreta y el área de cara tienen que ir en pareja, o la solución estática no se sostiene.

Una línea de código devuelve el equilibrio#

La solución es directa. Reemplazar la fuente por la misma expresión que la diferencia de flujos:

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

Ahora, en la solución estática, flujo y fuente se cancelan a precisión de máquina. La estrella deja de gotear.

El mismo truco aplica al término $\cos\theta P / (r \sin\theta)$ de la ecuación de momento $\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)

La ecuación de momento $\varphi$ no tiene fuente geométrica de presión: las dos caras azimutales comparten área.

Este es el ejemplo más simple de un esquema well-balanced: la aceleración ficticia que produce la divergencia discreta se resta con la misma expresión discreta al otro lado de la ecuación.

Sesenta líneas de NumPy — ¿sobrevive el equilibrio?#

En una malla radial 1D arrancamos con $P = 1$ , $\rho = 1$ , $u = 0$ y comparamos el residuo de momento con cada formulación.

import numpy as np
 
def polar_cell_volume(rL: float, rR: float) -> float:
    """volumen de cáscara esférica (sin θ Δθ Δφ omitidos)."""
    return (rR**3 - rL**3) / 3.0
 
def polar_r_surface(r: float) -> float:
    """área de cara ponderada por r² en coordenadas esféricas."""
    return r * r
 
def naive_pressure_source(P: float, r_mid: float, V: float) -> float:
    """Forma de manual 2P/r · V — errónea a precisión de máquina con P=const."""
    return 2.0 * P / r_mid * V
 
def balanced_pressure_source(P: float, S_in: float, S_out: float) -> float:
    """Misma expresión que la diferencia de flujos, así que se cancela exacto."""
    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):
        # flujo de presión a través de las caras (P constante → valor exacto)
        flux_diff = P * (S[1:] - S[:-1])
        # la línea de fuente que queremos auditar
        if scheme == "naive":
            src = naive_pressure_source(P, r_mid, V)
        else:
            src = balanced_pressure_source(P, S[:-1], S[1:])
        # actualización de momento: 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}")

Tras 200 pasos en una malla de 40 celdas, la salida es aproximadamente

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.

Con el interruptor en naive, al sumar pasos la curva se aparta de cero. Al pasar a well-balanced se queda en el eje. Al bajar cells solo crece el residuo naive. En 40 celdas y 200 pasos ese apartamiento es $\max|u|=1.834\times 10^{-3}$ ; la forma balanceada imprime $0.000\times 10^{0}$ .

Sustituir el momento-φ por momento angular#

La presión no es la única fuente geométrica. La ecuación de momento- $\varphi$ contiene $\rho u w / r + \cos\theta \rho v w / r$ . En sistemas con rotación fuerte — estrellas, discos de acreción, atmósferas — esos términos erosionan precisión y conservación.

El remedio es un cambio de variable. Definir $l = w r \sin\theta$ — el momento angular sin más — hace que las fuentes en la ecuación de momento- $\varphi$ desaparezcan por completo:

\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

En forma discreta basta con tomar los flujos del solver de Riemann $F_{\varphi r}$ , $F_{\varphi\theta}$ , $F_{\varphi\varphi}$ , multiplicarlos por $r \sin\theta$ para obtener $\tilde F$ y diferenciarlos. Sin correcciones, sin malabares con fuentes. Incluso en mallas gruesas el momento angular se conserva por construcción (Kley 1998, A&A 338, L37).

Los sistemas que rotan diferencialmente (discos de Kepler) necesitan además el truco FARGO del sistema rotante para aflojar la CFL. Eso queda para otro post.

Guía de decisión entre formas#

Cuestión	naive `2P/r · V`	well-balanced `P · ΔS_r`	forma angular
equilibrio estático	❌ residuo $O(\Delta r^{2})$	✅ exacto a nivel de bit	✅ (solo dirección φ)
código extra	ninguno	reemplazar una línea de fuente	nuevo flujo + variable adicional
sistemas rotantes	⚠ deriva acumulativa	△ solo arregla la fuente de presión	✅ conservación intrínseca
coste de depuración	alto (falla en silencio)	bajo	medio
dónde aparece	"¿por qué se vacía el núcleo de la estrella?"	estándar en hidro esférica	discos de acreción

De vuelta frente a la estrella#

Una malla cartesiana esconde esta trampa. Las coordenadas curvilíneas añaden términos a la ecuación de momento, y a menos que su forma discreta coincida exactamente con la diferencia discreta de flujos, un gas en reposo no permanece en reposo. La sola línea editada en este artículo es

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

Una línea, un commit. La estrella vuelve a estar quieta.

¿Se deja 2P/r o se cambia por ΔS?#

Si solo tiembla el balance de presión, $2P/r\,V$ pasa a $P(S_{i+1/2}-S_{i-1/2})$ . Ese es el cambio que convierte $\max|u|=1.834\times 10^{-3}$ en 40 celdas en 0. La fuente en $\theta$ , $\cos\theta\,P/(r\sin\theta)$ , admite el mismo reemplazo, $P\cdot\Delta S_{\theta}$ . La dirección $\varphi$ no lo necesita: las dos caras azimutales comparten área.

Con rotación fuerte la fuente de presión no basta. $\rho u w/r + \cos\theta\,\rho v w/r$ sigue en el momento- $\varphi$ . Definir $l = w r \sin\theta$ saca esa fuente de la ecuación. Es la forma de Kley 1998 (A&A 338, L37). Si la rotación cambia con el radio, todavía hace falta FARGO. Una interfaz en reposo que empieza a fluir sola es la misma enfermedad que las corrientes parásitas. Si la fuerza falsa de la divergencia discreta no se resta con la misma expresión discreta, el equilibrio no vuelve a su sitio.