Teoría de perturbaciones y el achatamiento de la Tierra

1. Motivación

El plano orbital de la Estación Espacial Internacional gira hacia el oeste unos 5° cada día; los satélites Landsat y Sentinel-2 pasan por cada lugar siempre a la misma hora solar; los satélites Molniya mantienen su apogeo sobre Rusia durante años. Nada de esto ocurre en el problema de dos cuerpos: son consecuencias de que la Tierra no es esférica. Esta lección construye la herramienta general para tratar cualquier perturbación pequeña (las ecuaciones planetarias) y la aplica al término dominante del campo terrestre, el armónico zonal $J_2$.

2. Elementos osculadores y función perturbadora

Escribimos la ecuación de movimiento como el problema de dos cuerpos más una aceleración pequeña:

En cada instante definimos los elementos osculadores $(a, e, i, \Omega, \omega, M)$ como los de la cónica que tendría el cuerpo si $\mathbf f$ se anulara a partir de ese momento (lección 2). Si $\mathbf f = \mathbf 0$ son constantes; en general varían lentamente, y la idea de la teoría de perturbaciones es escribir ecuaciones diferenciales para ellos. Cuando $\mathbf f = \nabla\mathcal R$ deriva de un potencial, $\mathcal R$ se llama función perturbadora (con el signo tal que $U = -\mu/r - \mathcal R$).

3. Ecuaciones de Gauss

Descomponemos $\mathbf f = R\,\hat{\mathbf r} + T\,\hat{\boldsymbol\theta} + N\,\hat{\mathbf h}$ en las componentes radial, transversal (en el plano orbital, en el sentido del movimiento) y normal. Deduzcamos las dos primeras ecuaciones; las demás siguen el mismo método (Problema 1).

Semieje mayor. La energía osculadora $\varepsilon = -\mu/(2a)$ cambia según la potencia de $\mathbf f$: $\dot\varepsilon = \mathbf v\cdot\mathbf f = \dot r\,R + r\dot\theta\,T$. Con $\dot r = (\mu/h)\,e\sin\nu$ y $r\dot\theta = h/r = (\mu/h)(1+e\cos\nu)$ (lección 0):

Momento angular y excentricidad. $\dot{\mathbf h} = \mathbf r\times\mathbf f$, cuya componente a lo largo de $\hat{\mathbf h}$ es $rT$: $\dot h = rT$. Como $h^2 = \mu a(1-e^2)$, $2h\dot h = \mu\dot a(1-e^2) - 2\mu a e\dot e$; despejando $\dot e$ y sustituyendo $\dot a$ y $\dot h$:

Las componentes de $\dot{\mathbf h}$ perpendiculares a $\hat{\mathbf h}$ (solo $N$ las produce) inclinan el plano orbital y dan $\dot i$ y $\dot\Omega$; $\dot\omega$ y $\dot M$ se obtienen de la rotación de $\mathbf A$. El sistema completo (Murray & Dermott, ec. 2.145–2.150) es

Lectura física: solo $N$ cambia la orientación del plano ($i$, $\Omega$); $T$ es la componente eficaz para cambiar $a$ (base de la transferencia de Hohmann, lección 6); los denominadores $e$ y $\sin i$ son las mismas singularidades de la lección 2.

4. Ecuaciones de Lagrange

Cuando $\mathbf f = \nabla\mathcal R$, conviene expresar las variaciones mediante las derivadas parciales de $\mathcal R$ respecto a los propios elementos. Lagrange (1808) obtuvo, con $n = \sqrt{\mu/a^3}$:

Su estructura es hamiltoniana: los pares $(a,M)$, $(e,\omega)$, $(i,\Omega)$ aparecen cruzados con signos opuestos (más limpio aún en las variables de Delaunay $L = \sqrt{\mu a}$, $G = L\sqrt{1-e^2}$, $H = G\cos i$, conjugadas a $M$, $\omega$, $\Omega$). La ventaja decisiva: si $\mathcal R$ no depende de un ángulo, su conjugado se conserva. Y si promediamos $\mathcal R$ sobre $M$ (una órbita), $\partial\langle\mathcal R\rangle/\partial M = 0$ y por tanto $\langle\dot a\rangle = 0$: el semieje mayor no sufre variaciones seculares a primer orden para ninguna perturbación conservativa. Es el teorema de Laplace–Lagrange sobre la estabilidad de los semiejes.

5. El achatamiento: potencial con $J_2$

El potencial exterior de un cuerpo con simetría axial se desarrolla en armónicos zonales:

con $R$ el radio ecuatorial, $\phi$ la latitud y $P_k$ los polinomios de Legendre. Para la Tierra $J_2 = 1.08262545\times10^{-3}$ (IERS 2010, vía JPL Horizons) y todos los demás son mil veces menores ($J_3 \approx -2.5\times10^{-6}$), así que nos quedamos con $P_2(x) = (3x^2-1)/2$:

donde la última igualdad es trigonometría esférica en el triángulo nodo–cuerpo–ecuador. $J_2 > 0$ significa que hay masa de más en el ecuador: un satélite siente una atracción extra hacia el plano ecuatorial.

5.1 Promedio sobre una órbita

Necesitamos $\langle r^{-3}\rangle$ y $\langle r^{-3}\sin^2(\omega+\nu)\rangle$ con $\langle\cdot\rangle = \frac1T\int_0^T(\cdot)\,dt$. Con $dt = r^2d\nu/h$ (ley de áreas) y $r = p/(1+e\cos\nu)$:

usando $hp = h^3/\mu$ y $h = na^2\sqrt{1-e^2}$. El mismo cálculo con $\sin^2(\omega+\nu)$ da la mitad (Problema 2), así que

que no depende de $M$, $\omega$ ni $\Omega$. Por las ecuaciones de Lagrange, $a$, $e$ e $i$ no tienen variación secular; solo precesan los tres ángulos.

5.2 Tasas seculares

Con $\partial\langle\mathcal R\rangle/\partial i = -\tfrac32 n^2J_2R^2(1-e^2)^{-3/2}\sin i\cos i$ y $\partial\langle\mathcal R\rangle/\partial e = \tfrac34 n^2J_2R^2\,e\,(1-e^2)^{-5/2}(2-3\sin^2i)$, las ecuaciones de Lagrange dan, con $p = a(1-e^2)$:

Son las fórmulas de Kozai (1959) y Brouwer (1959) a primer orden en $J_2$ (implementadas en src/lib/math/perturbaciones.ts). Tres consecuencias:

1. Regresión nodal. Para órbitas progradas ($\cos i > 0$) el nodo retrocede hacia el oeste; para las retrógradas avanza; para las polares no se mueve. El efecto es máximo en órbitas bajas ecuatoriales.
2. Inclinación crítica. $\dot\omega = 0$ cuando $5\cos^2i = 1$, es decir, $i = 63.435°$ o $116.565°$. Un satélite Molniya ($a \approx 26\,600$ km, $e \approx 0.74$, $i = 63.4°$, $T = 12$ h) congela así su apogeo sobre latitudes altas del hemisferio norte.
3. Órbita heliosíncrona. Si elegimos $\dot\Omega = +2\pi/\text{año tropical} = 0.98565°/\text{día}$ (positivo, por tanto $\cos i < 0$: retrógrada), el plano orbital gira al mismo ritmo que la dirección Tierra–Sol y el satélite cruza el ecuador siempre a la misma hora solar local. Para órbita circular:

6. Resultados con datos reales

Con $\mu = 398\,600.435436$ km³/s², $J_2 = 1.08262545\times10^{-3}$, $R = 6\,378.137$ km (JPL Horizons, IERS 2010) y año tropical de 365.2422 d, las fórmulas de la sección 5.2 dan:

• La ISS pierde $0.32°$ de nodo por órbita: en 73 días su plano da una vuelta completa, lo que fija las ventanas de lanzamiento desde cada cosmódromo.
• Landsat 8 (USGS: 705 km, $i = 98.2°$, 16 días de ciclo): la fórmula heliosíncrona da $i_{ss} = 98.21°$. Coincidencia a la centésima.
• Sentinel-2 (ESA: 786 km, $i = 98.62°$, 100.6 min, 10 días de ciclo): la fórmula da $98.54°$, $0.08°$ por debajo del valor de ESA. La sensibilidad es $di_{ss}/dh = 0.0042°/\text{km}$, así que la discrepancia equivale a 18 km de altitud; procede de que la órbita real se diseña con el modelo completo de campo (más armónicos, excentricidad «congelada» $e\approx10^{-3}$) y de que «786 km» es una altitud media sobre el geoide, no $a - R_{eq}$. La teoría $J_2$ de primer orden es exacta a $\sim0.1°$; el resto es ingeniería de detalle. Con 100.6 min, Sentinel-2 completa $143.14$ órbitas en 10 días, de modo que repite su traza cada 143 revoluciones.
• Inclinación heliosíncrona frente a altitud (órbita circular): 400 km → 97.03°; 600 km → 97.79°; 800 km → 98.60°; 1 000 km → 99.48°; 1 500 km → 101.96°. El límite es $\cos i = -1$: no existen órbitas heliosíncronas circulares por encima de 5 970 km de altitud.
• En Marte ($J_2 = 1.96045\times10^{-3}$, $R = 3\,396.2$ km, $\mu = 42\,828.376$ km³/s², año de 686.98 d) la misma fórmula da, a 400 km, $i_{ss} = 92.9°$: el $J_2$ marciano es mayor, pero su año es casi el doble y el resultado es una órbita casi polar.

El $J_2$ del Sol y Mercurio. Para un planeta, el ángulo relevante es la longitud del perihelio $\varpi = \Omega + \omega$ (con $i = 0$ respecto al ecuador solar, $\dot\varpi = \dot\Omega + \dot\omega = \tfrac32 nJ_2(R/p)^2$). Con $J_{2\odot} = 2.25\times10^{-7}$ (Park et al. 2017), $R_\odot = 695\,700$ km, $a = 0.38710$ UA y $e = 0.2056$, el avance es $3\pi J_2(R_\odot/a)^2/(1-e^2)^2 = 3.34\times10^{-10}$ rad por órbita: $0.029''$ por siglo, frente a los $42.98''$ relativistas. Por eso la medición de $J_{2\odot}$ y la prueba de la relatividad general con Mercurio van juntas.

7. Análisis numérico y de error

Teoría frente a integración. El test perturbaciones.test.ts integra con RK4 la ecuación completa $\ddot{\mathbf r} = -\mu\mathbf r/r^3 + \mathbf f_{J_2}$, con la aceleración cartesiana

durante 30 órbitas de un satélite a 800 km con $i = 98.6°$ (400 pasos por órbita) y compara el cambio de $\Omega$ osculador con $\dot\Omega_{sec}\,t$: coinciden al 1 %. La diferencia tiene tres orígenes, de distinto tamaño:

1. Términos periódicos omitidos al promediar: $\Omega$ osculador oscila con amplitud $\sim\tfrac32J_2(R/a)^2 \approx 1.3\times10^{-3}$ rad alrededor de la deriva secular. Es la cota de la comparación instantánea; promediando sobre una órbita desaparece.
2. Segundo orden en $J_2$: la teoría de Brouwer (1959) añade términos $\propto J_2^2 \approx 10^{-6}$, que para la ISS valen $\sim10^{-3}$°/día. Despreciable para diseñar, no para predecir a un año.
3. Error del integrador: RK4 con $h = T/400$ tiene error global $\propto h^4 \sim 10^{-10}$ relativo por órbita, muy por debajo de los anteriores; con $h = T/50$ el error del integrador ya es comparable a los términos periódicos. Compruébalo en el Problema 6.

Lo que $J_2$ no incluye. Para satélites bajos, el arrastre atmosférico sí cambia $a$ secularmente (la ISS pierde $\sim2$ km/mes y necesita reimpulsos); $J_3$ produce una oscilación de $e$ de largo periodo acoplada a $\omega$; el Sol y la Luna dominan la precesión nodal por encima de $\sim20\,000$ km (en GEO su efecto es comparable al de $J_2$: la inclinación de un satélite geoestacionario sin control deriva $\sim0.85°$/año). Toda cifra de la tabla es, por tanto, «solo $J_2$, primer orden».

8. Simulador

El simulador integra solo el problema de dos cuerpos, así que el periapsis no precesa. Úsalo para fijar las magnitudes de referencia: con $a = 1$ y $e = 0.74$ (Molniya a escala) mide la fracción del periodo que el cuerpo pasa en la mitad de la órbita más alejada del foco (debería ser $\approx 0.9$, Problema 5) y entiende por qué congelar el apogeo con $i = 63.4°$ vale tanto para las comunicaciones en latitudes altas.

Problemas

1. (Demostración) Deduce $di/dt$ y $d\Omega/dt$ de Gauss. Sugerencia: $\dot{\mathbf h} = \mathbf r\times\mathbf f = rT\,\hat{\mathbf h} - rN\,\hat{\boldsymbol\theta}$; descompón $\hat{\boldsymbol\theta}$ en la base de la línea de nodos y proyecta el cambio de $\hat{\mathbf h}$ sobre las direcciones que cambian $i$ y $\Omega$.

2. (Demostración) Calcula $\langle r^{-3}\sin^2(\omega+\nu)\rangle$ y completa la deducción de $\langle\mathcal R_{J_2}\rangle$.

3. (Inclinación crítica) Demuestra que $\dot\omega = 0$ exactamente en $i = \arccos(1/\sqrt5) = 63.435°$ y que en ese caso $\dot\Omega \ne 0$. ¿Cuánto tarda el nodo de un Molniya en dar una vuelta?

4. (Heliosíncrona) Deduce la altitud máxima de una órbita heliosíncrona circular y calcula $i_{ss}$ para una órbita circular a 600 km. ¿Qué ocurre con $i_{ss}$ si se mantiene $a$ pero $e = 0.1$?

5. (Molniya) Calcula, con la anomalía excéntrica, la fracción del periodo que un cuerpo pasa con $\nu\in[\pi/2, 3\pi/2]$ (la mitad «lejana» de la órbita, más allá del semilado recto) y evalúala para $e = 0.74$.

6. (Computacional) Integra con RK4 el movimiento de un satélite a 800 km, $i = 98.6°$, $e = 0$, con la aceleración $\mathbf f_{J_2}$ de la sección 7 (o usa campoJ2 de src/lib/math/perturbaciones.ts), convierte cada estado a elementos osculadores y representa $\Omega(t)$, $\omega(t)$, $a(t)$ durante 30 órbitas. Ajusta una recta a $\Omega(t)$ y compara su pendiente con $\dot\Omega_{sec}$. Repite con 400, 100 y 50 pasos por órbita y explica qué cambia.

import numpy as np
mu, J2, R = 3.98600435436e14, 1.08262545e-3, 6378137.0
def f(s):
    x, y, z, vx, vy, vz = s; r2 = x*x+y*y+z*z; r = np.sqrt(r2)
    k = 1.5*J2*mu*R*R/(r2*r2*r); z2 = 5*z*z/r2
    ax = -mu*x/r**3 + k*x*(z2-1); ay = -mu*y/r**3 + k*y*(z2-1); az = -mu*z/r**3 + k*z*(z2-3)
    return np.array([vx, vy, vz, ax, ay, az])
def rk4(s, h):
    k1 = f(s); k2 = f(s+h/2*k1); k3 = f(s+h/2*k2); k4 = f(s+h*k3)
    return s + h/6*(k1+2*k2+2*k3+k4)
# estado inicial desde elementos (lección 2), bucle, estado_a_elementos en cada paso, np.polyfit(t, Omega, 1)

7. (Mercurio) Reproduce el cálculo de $0.029''$/siglo para la precesión del perihelio de Mercurio por $J_{2\odot}$ y determina con qué precisión relativa habría que conocer $J_{2\odot}$ para que su incertidumbre fuera inferior al $10^{-4}$ del efecto relativista ($42.98''$/siglo).

Autoevaluación

Lecturas

• Murray & Dermott (1999), §2.9 y cap. 6: ecuaciones de Gauss y de Lagrange con deducción completa; §6.11 trata el achatamiento con la misma función perturbadora promediada.
• Vallado (2013), cap. 9: la referencia de ingeniería para $J_2$, órbitas heliosíncronas, Molniya y teorías de Brouwer/Kozai; §8.7 da la aceleración cartesiana.
• Curtis (2020), cap. 10: deducción elemental de la regresión nodal y ejemplos numéricos con satélites reales.
• Kozai (1959), DOI 10.1086/107957, y Brouwer (1959), DOI 10.1086/107958: las dos teorías originales del satélite artificial, a primer y segundo orden en $J_2$.
• Park et al. (2017), DOI 10.3847/1538-3881/aa5be2: $J_{2\odot}$ y la precesión del perihelio de Mercurio medidos con MESSENGER.
• IERS Conventions (2010), Technical Note 36: constantes de referencia del campo terrestre.