Indice | Previo: ProblemaNCuerpos.SolucionAnalitica | Siguiente: ProblemaNCuerpos.SolucionAnalitica.Energia

6.2.1. Aplicación de los teoremas de conservación

En la forma específica del problema de los N cuerpos presentada aquí, asumimos que el sistema de partículas esta aislado del Universo, es decir, sobre los cuerpos no actúa ninguna fuerzas externa.

Como vimos en la Sección Teoremas de conservación en sistemas de partículas, esta condición implica que tanto el momentum lineal total del sistema (ver Teo. ncuerpos_formulacion_ecuaciones), como su momentum angular total (ver Teo. ncuerpos_formulacion_ecuaciones) se conservan.

De acuerdo con esto, las funciones:

$$ \begin{array}{rcl} \sum m_i \dot{\vec r}_i & = & \vec P_\mathrm{CM}\\ \sum m_i {\vec r}_i\times\dot{\vec r}_i & = & {\vec L} \end{array} $$

son las primeras constantes de movimiento que reconocemos para el problema.

Dada la importancia que el método de cuadraturas tiene en el desarrollo de muchos de los resultados en este libro, en las siguientes secciones confirmaremos este resultado, lo interpretaremos en el contexto del problema de los N cuerpos y en particular estudiaremos su implicación para algunos sistemas astronómicos de interés, pero más importante, encontraremos dos constantes adicionales que no habíamos deducido hasta ahora en el caso general de sistemas de muchas partículas.

6.2.2. Momento lineal

Como vimos en la Sección Fuerzas y centro de masa, la manera de probar la conservación del momento angular en el problema de los N cuerpos, consiste en sumar las e.d.m. de cada una de las partículas (Ecs. ncuerpos_formulacion_ecuaciones):

\begin{equation} \label{eq:ncuerpos_suma_edm} \sum_i m_i\ddot{\vec r}_i = -\sum_i \sum_j \frac{m_i \mu_j}{r_{ij}^3} \vec{r}_{ij} \end{equation}

En el lado derecho de esta ecuación, el término en el denominador es simétrico, $r_{ji}=r_{ij}$, mientras el término en el numerador, $\vec{r}_{ji}=\vec{r}_i-\vec{r}_j$ es antisimétrico (el signo cambia al cambiar el orden de los índices). Como resultado, por cada término en la doble sumatoria (por ejemplo el término $i=1$, $j=4$) habra un término idéntico pero de signo contrario (el término con $i=4$, $j=1$) que lo anulará. Físicamente, esta propiedad matemática es equivalente a aplicar postulado de acción y reacción que usamos para deducir la Ec. (edm_sistema).

Así, el lado derecho de la Ec. (ncuerpos_suma_edm) siempre será nulo (sin importar el número de partículas) y la suma de las e.d.m. será:

\begin{equation} \sum m_i\ddot{\vec r}_i = 0. \end{equation}

Como era de esperarse esta ecuación es equivalente a la e.d.m. de un sistema de partículas sobre el que no actúan fuerzas externas (Ec. edm_sistema.)

En forma de cuadraturas la misma ecuación se escribe:

$$ \frac{d}{dt}\left(\sum_i m_i \dot{\vec{r}}_i\right) =0, $$

donde como siempre asumimos que las masas de todas las partículas del sistema $\{m_i\}$ son constantes.

Esta expresión no es otra cosa que una forma del teorema de conservación del momento lineal y provee la primera constante (vectorial) de movimiento del problema de N cuerpos:

\begin{equation} \label{eq:ncuerpos_momento} \sum_i m_i \dot{\vec{r}}_i={\vec P}_\mathrm{CM} \end{equation}

Por tratarse de una expresión vectorial, en esta cuadratura hay contenidas en realidad tres constantes de movimiento: $\sum_i m_i\dot{x}_i=P_{\mathrm{CM},x}$, $\sum_i m_i\dot{y}_i=P_{\mathrm{CM},y}$ y $\sum_i m_i\dot{z}_i=P_{\mathrm{CM},z}$.

Nota: El sistema de referencia del centro de masa. La constancia de ${\vec P}_{\mathrm{CM}}$ en el problema de los N cuerpos, nos permite introducir un sistema de referencia inercial ideal para la descripción simplificada del problema. Naturalemente nos referimos al sistema de referencia del centro de masa, que tiene, respecto al sistema de referencia general considerado hasta aquí, una velocidad constante e igual a ${\vec V}_{\mathrm{CM}}={\vec P}_{\mathrm{CM}}/M$ (donde $M=\sum m_i$).

El movimiento de las partículas en este sistema de referencia es el más simple que podemos percibir desde un sistema de referencia inercial. Por la misma razón, en lo sucesivo y siempre y cuando no se diga lo contrario, la dinámica de la mayoría de los sistemas de N cuerpos considerados en este libro se describirá en el sistema de referencia de su centro de masa.

6.2.3. Posición del centro de masa

Una segunda constante de movimiento puede obtenerse aplicando cuadraturas a la primera integral (Ec. ncuerpos_momentum):

$$ \frac{\mathrm{d}}{\mathrm{d}t}\left(\sum_i m_i \vec{r}_i\right)=\frac{\mathrm{d}}{\mathrm{d}t}\left({\vec P}_\mathrm{CM} t\right) $$

De donde podemos escribir:

\begin{equation} \label{eq:ncuerpos_centro_masa} \sum_i m_i \vec{r}_i-{\vec P}_\mathrm{CM} t=M\vec R_0 \end{equation}

Para expresar está última constante de movimiento, hemos elegido, arbitrariamente, llamar a su valor como $M\vec R_0$, donde $M$ es la masa total del sistema (que también es un valor constante). La elección de esta parametrización para el valor de la integral, no modifica en nada el hecho que el lado derecho de la Ec. (ncuerpos_centro_masa) es constante también.

¿Cómo podemos interpretar físicamente la constante de movimiento en la Ec. (ncuerpos_centro_masa)? Si dividimos la expresión de esta constante por la masa total del sistema $M$ y la escribimos como:

$$ \frac{\sum_i m_i \vec{r}_i}{M}=\vec R_0+\vec V t, $$

podemos identificar, del lado izquierdo, la posición del centro de masa del sistema en cualquier tiempo ${\vec R}_{\rm CM}=\sum m_i\vec{r}_i/M$ (ver Ec. centro_masa). Del lado derecho, por otro lado, encontramos a ${\vec V}_\mathrm{CM} t$, que no es otra cosa que el desplazamiento que sufre el centro de masa al moverse con velocidad constante ${\vec V}_\mathrm{CM}={\vec P}_\mathrm{CM}/M$. En consecuencia podemos entonces concluir que $\vec R_0$ no es otra cosa es la posición inicial del centro de masa. Así, la Ec. (ncuerpos_centro_masa) es la constante de movimiento asociada al centro de masa.

De nuevo, por tratarse de una expresión vectorial, la Ec. (ncuerpos_centro_masa) corresponde en realidad a tres constantes en lugar de una: $\sum_i m_i x_i-P_{\mathrm{CM},x} t=M x_{\mathrm{CM},0}$, $\sum_i m_i y_i-P_{\mathrm{CM},y} t=M y_{\mathrm{CM},0}$ y $\sum_i m_i z_i-P_{\mathrm{CM},z} t=M z_{\mathrm{CM},0}$.

En el sistema de referencia del centro de masa que habíamos mencionado en una nota anterior, esta constante de movimiento se puede escribir de forma simplificada como:

\begin{equation} \label{eq:ncuerpos_centro_masa} \sum_i m_i {{\vec r}_i}'=\vec 0, \end{equation}

donde, como hicimos en la Sección Dinámica en el centro de masa, las cantidades primadas son medidas respecto al centro de masa.

Esta expresión confirma la intuición expresada antes de que en este sistema de referenia la descripción del problema es la más simple posible.

6.2.4. Momentum angular

Como hicimos en la Sección *Teoremas de conservación en sistemas de partículas, si tomamos las e.d.m. del sistema de partículas (Ec. ncuerpos_formulacion_ecuaciones) y premultiplicamos vectorialmente por el vector posición y la masa de cada partícula (este será el factor integrante), el resultado es el conjunto de ecuaciones:

$$ \left\{m_i\vec{r}_i\times \ddot{\vec{r}}_i= -m_i\sum_j \frac{\mu_j}{r_{ij}^3} \vec{r}_i\times \vec{r}_{ij} \right\} $$

Por las propiedades del producto vectorial, el último término en el lado derecho será, $\vec{r}_i\times \vec{r}_{ij}=-\vec{r}_i\times \vec{r}_j$, que es también antisimétrico. En consecuencia, si sumamos todas las ecuaciones diferenciales resultantes, obtendremos la ecuación:

$$ \displaystyle \sum_i m_i\vec{r}_i \times \ddot{\vec{r}}_i = \vec{0}, $$

que podemos expresar en cuadraturas como:

$$ \frac{d}{dt}\left(\sum_i m_i \vec{r}_i \times \dot{\vec{r}}_i\right) = \vec 0 $$

Esta última ecuación permite identificar la constante de movimiento:

\begin{equation} \label{eq:ncuerpos_momento_angular} \sum m_i \vec{r}_i \times \dot{\vec{r}}_i= \vec L \end{equation}

que reconocemos, naturalmente, con el momento angular total del sistema.

Existe una interesante interpretación geométrica de esta constante de movimiento: si nos ubicamos en el sistema de referencia del centro de masa del sistema, la dirección (constante) del vector momento angular total define un plano invariable (en dirección perpendicular a él). No importa la posición de las partículas, ni el tiempo en el que se registren, dicho plano mantendrá siempre su orientación en el espacio. A este plano se lo conoce históricamente como el plano invariable de Laplace (Laplace, 1835).

Definición 11.10. Plano invariable (de Laplace) y sistema natural de referencia. Dado un sistema de N partículas puntuales $\{m_i\}$, llamamos plano invariable (de Laplace).. a aquel que tiene como vector normal $\hat n={\vec L}'/L'$, donde:

$${\vec L}'=\sum m_i {{\vec r}_i}' \times {\dot{\vec{r}}_i}'$$

es el momento angular total de las partículas referidas al sistema de referencia del centro de masa. Dado que ${\vec L}'$ es constante, es posible usar las posiciones y velocidades de las partículas en un tiempo arbitrario para calcular la dirección del plano invariable.

Llamamos sistema natural de coordenadas de un sistema de N cuerpos a aquel que tiene: 1) origen en el centro de masa y 2) plano $x-y$ sobre el plano invariable de Laplace (ver Figura Figura (plano_invariable)).)

Ilustración gráfica de la orientación del plano invariable de Laplace. El plano invariable esta definido en el sistema de referencia inercial del centro de masa y se mueve con él con velocidad $V_

Figura 6.65. Ilustración gráfica de la orientación del plano invariable de Laplace. El plano invariable esta definido en el sistema de referencia inercial del centro de masa y se mueve con él con velocidad $V_

Es importante insistir en que el vector que define la dirección del plano invariable de Laplace, no es, en general el vector en la Ec. (ncuerpos_momento_angular) que, como vimos en la Sección Dinámica referida al centro de masa se puede escribir como:

$$ \vec L=M{\vec R}_\mathrm{CM}\times{\vec V}_\mathrm{CM}+\sum m_i {{\vec r}_i}' \times {\dot{\vec{r}}_i}'. $$

El vector de la Definición ncuerpos_momento_angular, que es esencialmente el segundo término de la ecuación anterior, no tiene la misma dirección de $\vec L$ debido al término $M{\vec R}_\mathrm{CM}\times{\vec V}_\mathrm{CM}$ (momento angular del centro de masa.) Por lo tanto, es muy importante al construir el sistema de coordenadas natural, tener cuidado de asegurárse que las posiciones y velocidades de las partículas estén referidas al centro de masa.

Un poco de historia: El plano invariable del Sistema Solar. Durante más de 100 años, los astrónomos han intentado encontrar la orientación del plano invariable del Sistema Solar. El esfuerzo no ha sido sencillo, en tanto durante en el mismo lapso, la masa de los cuerpos no siempre se ha conocido con precisión e incluso, cuerpos enteramente nuevos se descubren de vez en cuando.

La determinación más precisa de la orientación del plano invariable se hizo recientemente (Souami & Souchay, 2012) y concluyó que este importante plano tiene una inclinación, respecto al plano del ecuador en el sistema ICRF (ver la Sección Sistema de referencia astronómicos) de 23$^\circ$ 0' 31".9 y una longitude del nodo ascendente de 3$^\circ$ 51' 9".4 (para una definición de estos ángulos ver la Sección Cónicas en el espacio.) Esto implica, que el plano invariable del Sistema Solar esta inclinado respecto al plano de la eclípica de J2000.0 (que es muy importante en la mecánica celeste al usarse como el plano de referencia fundamental para el cálculo de posiciones planetarias) en un ángulo de 1$^\circ$ 34' 43".3.

Usando SPICE podemos hacer una estimación de primer orden de la inclinación del plano invariable del Sistema Solar, respecto al sistema de la eclíptica de J2000.0 (ver la Sección Sistemas de referencia astronómicos.) Para ello podemos asumir, por simplicidad, que el momentum angular total del sistema esta concentrado en el Sol y el planeta más masivo, Júpiter (entre los dos cuerpos concentran el 99.96% de la masa del sistema.)

Para conseguirlo, primero debemos obtener de los kernels de SPICE, la masa, posición y velocidad de los cuerpos implicados en un tiempo de efemérides arbitrario (tomaremos $t_\mathrm{ef}=0$ por comodidad):

In [1]:
import spiceypy as spy

#Asumimos un tiempo cualquiera, en este caso t=J2000.0
tef=0

#Carga kernels con posiciones (bsp) y masas (tpc)
spy.furnsh('pymcel/data/de430.bsp')
spy.furnsh('pymcel/data/de430.tpc')

#Parámetro gravitacional, posiciones y velocidades de los cuerpos
musol=spy.bodvrd("SUN","GM",1)[1][0]
sol,tluz=spy.spkezr("SUN",tef,
                    "ECLIPJ2000","None","SSB")
rsol=sol[:3]
vsol=sol[3:]
mujupiter=spy.bodvrd("JUPITER_BARYCENTER","GM",1)[1][0]
jupiter,tluz=spy.spkezr("JUPITER_BARYCENTER",tef,
                        "ECLIPJ2000","None","SSB")
rjupiter=jupiter[:3]
vjupiter=jupiter[3:]

Sol:
	mu = 132712440041.9394 km^3/s^2
	Posición = [-1067598.50226456  -418234.39327422    30837.61810502] km
	Velocidad = [ 0.00931257 -0.01282475 -0.00016333] km/s
Jupiter:
	mu = 126712764.7999999 km^3/s^2
	Posición = [ 5.97499980e+08  4.39186501e+08 -1.51960586e+07] km
	Velocidad = [-7.90052522 11.14330827  0.13069904] km/s

Nótese que como señalamos antes nos hemos cuidado de calcular tanto la posición del Sol como de Júpiter, respecto al centro de masa (baricentro) del sistema solar (SSB por Solar System Barycenter en SPICE.)

Con esta información, el momento angular total del sistema será:

In [3]:
#Constante de gravitación universal:
G=6.67e-20 # km^3 / kg s^2

#Momentum angular total
from numpy import cross
L=(musol/G)*cross(rsol,vsol)+(mujupiter/G)*cross(rjupiter,vjupiter)

L = [4.31661898e+35 7.99455257e+34 1.92754426e+37] kg m^2/s^2

En el algoritmo anterior, para obtener la masa de cada cuerpo, hemos divido el valor de los parámetros gravitacionales $\mu$ por la constante $G$ (en las unidades apropiadas.)

Con este resultado, la inclinación estimada del plano invariable del Sistema Solar, se obtiene, finalmente, calculando el ángulo formado entre el momento angular total y el eje z (que en SPICE es perpendicular al plano de la eclíptica de J2000.0), es decir, $i=\cos^{-1}(L_z/L)$:

In [5]:
from numpy import arccos
from numpy.linalg import norm

i=arccos(L[2]/norm(L))

i = 1.3046988768626038 grados
  = 1:18:16

Aquí hemos usado la notación grados:minutos:segundos, mas conveniente en el contexto de los algoritmos.

La inclinación obtenida, $i=1^\circ\;18'\;16''$, difiere de la calculada con todos los cuerpos del Sistema Solar (ver el recuadro el Un poco de historia, el plano invariable del Sistema Solar) en casi $16$ minutos de arco (que corresponde a un error relativo de un poco más del 17\%.) El lector podría obtener una estimación mejor de esta inclinación si incluye otros cuerpos en el cálculo (ver sección de problemas al final de este capítulo.)

Indice | Previo: ProblemaNCuerpos.SolucionAnalitica | Siguiente: ProblemaNCuerpos.SolucionAnalitica.Energia