Indice | Previo: Problema3Cuerpos.CRTBP.Numerico | Siguiente: Problema3Cuerpos.RegionesExclusion
8.7. La constante de Jacobi¶
¿Qué podemos "entender" del movimiento de la partícula de prueba en el CRTBP sin resolver numéricamente la ecuación de movimiento?. Como lo hemos hecho a lo largo de este libro, y lo hicimos en particular en los casos del problema de los N cuerpos y el problema de los dos cuerpos relativo, es posible, como mínimo, reconocer algunas propiedades generales del sistema buscando las cuadraturas de su ecuación de movimiento.
Empecemos por considerar nuevamente la ecuación en su forma más general (Eq. crtbp):
\begin{equation} \ddot{\vec{r}} = -\mu_1\frac{\vec{r}_{1}}{r^3_{1}} - \mu_2\frac{\vec{r}_{2}}{r^3_{2}} - \vec{\omega}\times(\vec{\omega}\times \vec{r}) -2\vec{\omega}\times \dot{\vec{r}} \end{equation}Con $\vec{r}_{1}=[x(t) - x_1]\hat{e}_x + y(t)\hat{e}_y + z(t) \hat{e}_z$ y $\vec{r}_{2}=[x(t) - x_2]\hat{e}_x + y(t)\hat{e}_y + z(t) \hat{e}_z$
Si usamos el factor integrande $\dot{\vec{r}}\cdot$ obtenemos en el lado izquierdo de la ecuación de movimiento:
\begin{equation} \dot{\vec{r}} \cdot \ddot{\vec{r}} = \frac{d}{dt}\left(\frac{1}{2}\dot{\vec{r}}^2 \right) \end{equation}De otra parte, el primer término del lado derecho: queda como
\begin{equation} -\frac{\vec{r}_{1} \cdot \dot{\vec{r}}}{r_{1}^3} = -\left[\frac{(x - x_1)\dot{x} + y\dot{y} + z \dot{z}}{r_{1}^3}\right] \equiv \frac{d}{dt}\left(\frac{1}{r_{1}}\right) \end{equation}siendo $r_{1} = \sqrt{(x-x_1)^2 + y^2 + z^2}$. El segundo término será análogo.
Analicemos ahora el término $\vec{\omega}\times(\vec{\omega}\times \vec{r})\cdot \dot{\vec{r}}$. Notemos primero que:
\begin{eqnarray} \vec{\omega} \times (\vec{\omega} \times \vec{r}) & = & \omega \hat{k}\times (x \omega \hat{j} - \omega y \hat{i}) \nonumber \\ & = & - x \omega^2 \hat{e}_x - y \omega^2 \hat{e}_y \nonumber \end{eqnarray}y por lo tanto:
$$ \vec{\omega}\times(\vec{\omega}\times \vec{r})\cdot \dot{\vec{r}} = - \omega^2 (x\dot x+y\dot y) = -\frac{1}{2}\omega^2\frac{d}{dt}(x^2+y^2) $$Por propiedades del producto punto es claro que multiplicando escalarmente el último término de la derecha (aceleración de Coriolís) por $\dot{\vec{r}}$ obtenemos:
\begin{equation} \displaystyle 2(\vec{\omega} \times \dot{\vec{r}}) \cdot \dot{\vec{r}} = 0 \end{equation}Reuniendo todos estos resultados obtenemos finalmente que:
\begin{equation} \frac{d}{dt}\left(\frac{1}{2}\dot{\vec{r}}^2\right) = \frac{d}{dt}\left[\frac{\mu_1}{r_{1}} + \frac{\mu_2}{r_{2}} + \frac{1}{2}\omega^2 (x^2 + y^2) \right] \end{equation}Agrupando términos y multiplicando ambos lados de la ecuación por dos se halla finalmente:
\begin{equation} \frac{d}{dt}\left[-\dot{\vec{r}}^2 + 2\frac{Gm_1}{r_{1}} + 2\frac{Gm_2}{r_{2}} + \omega^2 (x^2 + y^2) \right] = 0 \end{equation}De aquí encontramos que una cuadratura al CRTBP es:
\begin{equation} \label{eq:cuadratura_Jacobi} \frac{2\mu_1}{r_{1}} + \frac{2\mu_2}{r_{2}} + \omega^2 (x^2 + y^2)-v^2\equiv C_J \end{equation}A esta expresión general la llamamos la cuadratura de Jacobi. Al valor $C_J$ (que puede obtenerse por ejemplo con las condiciones iniciales del problema) y que es estrictamente constante en el CRTBP la llamamos la constante de Jacobi. Usaremos de forma intercambiable los términos cuadratura y constante para referirnos a la Ec. (cuadratura_Jacobi) aunque el lector debe recordar que son conceptos sutilmente diferentes.
¿Existen otras cuadraturas del CRTBP?. Lamentablemente no. La constante de Jacobi resulta ser la única cantidad conservada en el problema. Naturalmente esto es cierto cuando se describe la dinámica en el sistema rotante y bajo las particulares condiciones en el que se lo define (cuerpos masivos siguiendo una trayectoria totalmente independiente de la partícula de prueba).
¿Esta relacionada $C_J$ con la energía total del sistema?. No realmente. En primer lugar debe notarse que no es posible definir una energía cinética para los dos objetos masivos del sistema puesto que están, por construcción, en reposo. Aún así si es posible hablar de una energía potencial entre ellos que se mantiene constante, a pesar del movimiento de la partícula de prueba. Es claro que si se considera el problema en un sistema inercial, la energía total del sistema debería ser constante (como lo es en cualquier problema de N cuerpos). Pero al pasarnos al sistema rotante y fijar el movimiento de las partículas masivas hemos arruinado las condiciones que conducían a dicha conservación.
¿Qué podemos decir de la energía mecánica total de la partícula de prueba?. En el sistema rotante la energía mecánica de la partícula no es conservada debido a que sobre ella actúa la fuerza de coriolis que depende de la velocidad y que por la misma razón no pueden hacerse derivar de una función escalar (una de las condiciones para que una fuerza sea conservativa). En este sentido podemos decir que la fuerza de Coriolís es disipativa y que no hay conservación de la energía mecánica.
¿Qué podemos decir sobre el momento angular de la partícula de prueba? ¿podría conservarse como en otros sistemas estudiados antes?. La fuerza sobre la partícula de prueba en el CRTBP tampoco es una fuerza central, de modo que la conservación del momento angular no esta garantizada como lo esta tanto en el problema general de los N cuerpos como en el problema relativo de los dos cuerpos.
En resumen, de por sí podríamos considerar afortunado encontrar al menos una cuadratura en un problema que por su preparación ha roto con todas las convenientes propiedades que la forma específica de la fuerza gravitacional le había otorgado a los sistemas que habíamos estudiado hasta ahora.
8.8. El valor de la constante de Jacobi¶
Si la constante de Jacobi no es la energía del sistema ¿qué unidades tiene? ¿cuál es su orden de magnitud? Es interesante analizar estas cuestiones antes de seguir profundizando en el CRTBP. Entre más intución desarrollemos sobre el significado de esta importante cantidad, mejor preparados estaremos para entender las propiedades globales de un sistema descrito con el formalismo del CRTBP una vez conocido el valor de $C_J$.
Examinando la definición de la constante en la Ec. (constante_Jacobi) reconocemos $C_J$ tiene unidades de velocidad al cuadrado, o lo que es lo mismo unidades de energía específica (energía por unidad de masa). La cuadratura, en las unidades que introdujimos en la Sección Unidades canónicas del CRTBP se puede escribir como:
\begin{equation} \label{eq:constante_Jacobi} \frac{2(1-\alpha)}{r_{1}} + \frac{2\alpha}{r_{2}} + (x^2 + y^2)-v^2=C_J \end{equation}Como vemos la constante puede tener un valor real arbitrario, positivo o negativo (esto debido al signo del cuadrado de la velocidad que siempre es una cantidad positiva). Si asumimos una situación en la que $r_1\sim r_2\sim (x^2+y^2)\sim 1$, vemos que independiente del valor de $\alpha$, la constante tiene un valor:
$$ C_J\sim 5-v^2 $$o lo que es lo mismo, su valor estará, en situaciones realistas y en las unidades canónicas, entre 0 y unos pocos, y muy seguramente alrededor de $2-4$ como lo atestiguan algunos de los experimentos numéricos que realizaremos a continuación.
Para hacernos a una idea del valor de $C_J$ podemos calcular su valor para muchos puntos alrededor de un sistema con un valor de $\alpha$ específico y asumiendo, por sencillez que la partícula se encuentra en reposo en esos puntos. Con ese propósito construyamos primero una malla de puntos (matriz de pares) en un plano paralelo al plano $x-y$:
#Valor de z en el que se encuentra el plano
zplano=0
#Tamaño de la malla y número de filas y columnas
rango=1.5
NG=80
#Malla de puntos
from numpy import meshgrid,linspace,ones_like
X,Y=meshgrid(linspace(-rango,rango,NG),linspace(-rango,rango,NG))
Z=zplano*ones_like(X)
Las matrices X, Y y Z contienen los valores de las respectivas coordenadas de los puntos sobre la malla. Así, por ejemplo, la componente X[3,10] tiene el valor de la coordenada $x$ del punto en la cuarta fila (índice 3) y en la décima columna (índice 10).
Con las matrices de las coordenadas de la malla podemos proceder a calcular los valores de la constante de Jacobi correspondiente a un determinado valor (fijo) de la rapidez:
#Parámetro del sistema
alfa=0.3
#Rapidez de la partícula de prueba
v=0.0
#Valor de CJ
from numpy import sqrt
CJ=2*(1-alfa)/sqrt((X+alfa)**2+Y**2+Z**2)+\
2*alfa/sqrt((X-1+alfa)**2+Y**2+Z**2)+\
(X**2+Y**2)-\
v**2
Como vemos CJ es una matriz que tiene el valor de la constante en cada uno de los puntos de la malla coordenada. El resultado esbozado, coindice con nuestras expectativas: la constante tiene un valor de unos pecos y en el borde de la malla tiene un valor cercano a $5-v^2$.
Podemos representar gráficamente este resultado, construyendo a partir de estas matrices, contornos numéricos (líneas de igual valor de $C_J$) que nos permiten, nuevamente, ganar un poco de intuición sobre esta importante cantidad. El algoritmo a continuación lleva a cabo ese cometido:
%matplotlib inline
import matplotlib.pyplot as plt
fig=plt.figure(figsize=(6,5))
ax=fig.gca()
#Curva de cero velocidad
from numpy import linspace
c=ax.contourf(X,Y,CJ,levels=linspace(0,5,10))
plt.colorbar(c)
ax.plot([-alfa],[0],'o',color='blue',markersize=10*(1-alfa))
ax.plot([1-alfa],[0],'o',color='orange',markersize=max(2,10*alfa))
ax.set_xlabel("x")
ax.set_xlabel("y")
ax.set_title(f"z = {zplano}, v = {v}")
ax.set_xlim((-rango,rango))
ax.set_ylim((-rango,rango))
fig.tight_layout()
Debe notarse que hemos excluído de la malla los puntos que están muy cercanos a las partículas masivas y para los cuales el valor de $C_J$ es muy grande.
Busque las figuras interactivas y las animaciones incluídas en el sitio en línea del libro.
8.9. Cuadratura de Jacobi de un sistema simulado¶
Contrastermos ahora el resultado teórico obtenido en las secciones anteriores con lo que podemos obtener de la solución numérica a las ecuaciones de movimiento. Primero consideremos la solución a la Eq. (crtbp) que sabemos tiene tiene la cuadratura expresada por la constante de Jacobi:
#Propiedades del sistema y condiciones iniciales
alfa=0.3
ro=[1.0,0.0,0.0]
vo=[0.0,0.45,0.0]
#Tiempos de integración
from numpy import linspace
Nt=1000
ts=linspace(0,10,Nt)
#Resuelve numéricamente la ecuación de movimiento
from pymcel.export import crtbp_solucion
solucion=crtbp_solucion(alfa,ro,vo,ts)
#Extrae las posiciones y velocidades en el sistema rotante
rs=solucion[0]
vs=solucion[1]
Podemos ahora calcular el valor de la constante de Jacobi para cada uno de los puntos de la trayectoria de la partícula. Para ello, escribamos una rutina general que permita calcular la constante para el conjunto de valores de la posición y la velocidad de la partícula de prueba:
def constante_jacobi(alfa,r,vel):
from numpy import array
r=array(r)
vel=array(vel)
#Valor de x, y, z
x=r[:,0]
y=r[:,1]
z=r[:,2]
#Rapidez
from numpy.linalg import norm
v=norm(vel,axis=1)
#Posiciones relativas
from numpy import sqrt
r1=sqrt((x+alfa)**2+y**2+z**2)
r2=sqrt((x-1+alfa)**2+y**2+z**2)
#Valor de la constante
CJ=2*(1-alfa)/r1+2*alfa/r2+(x**2+y**2)-v**2
return CJ
Ahora aplicamos la rutina para calculas los valores obtenidos con la solución numérica de las ecuaciones de movimiento:
CJ=constante_jacobi(alfa,rs,vs)
Como era de esperarse, dentro del margen de las incertidumbres numéricas, el valor de la cuadratura es el mismo, independiente del tiempo en el que se calcule.
8.1. Cuadratura de Jacobi de un sistema real¶
El resultado al final de la última secciín, más allá de ofrecer una comprobación numérica de la teoría, es bastante predecible: el sistema simulado tiene exactamente todas las propiedades necesarias para que la cuadratura de Jacobi (Ec. constante_Jacobi) resulte constante.
Aún menos trivial sería comprobar si esa misma cuadratura se satisface incluso en el caso de un sistema de tres cuerpos real, es decir, uno en el que todas las partículas interactúan.
Consideremos el siguiente sistema, descrito usando la notació que introdujimos en la Sección *Solución numérica al problema de los N cuerpos:
sistema=[
dict(
m=1e-4,
r=[0.5,0.0,0.5],
v=[0,0.5,0],
),
dict(
m=3.0,
r=[0,0,0],
v=[0,0,0]),
dict(
m=1.0,
r=[1.0,0,0],
v=[0,2.0,0],
),
]
Resolvamos numéricamente las ecuaciones de movimiento correspondientes:
from numpy import linspace
Nt=1000
ts=linspace(0,10,Nt)
from pymcel.export import ncuerpos_solucion
rs,vs,rps,vps,constantes=ncuerpos_solucion(sistema,ts)
Ahora veamos una gráfica de las trayectorias en el sistema de referencia del centro de masa (sistema de referencia inercial):
from pymcel.export import plot_ncuerpos_3d
fig=plot_ncuerpos_3d(rps,vps);
Como era de esperarse el sistema se comporta aproximadamente como el CRTBP: hay dos partículas masivas que orbitan una respecto a otra en trayectorias aproximadamente circulares, mientras que una tercera más liviana se mueve siguiendo una trayectoria muy compleja (ciertamente no kepleriana) afectada por las otras dos. El efecto de esta última sobre la trayectoria de las dos primeras parece despreciable.
Para examinar si el comportamiento de este sistema es similar al que hemos visto en el CRTBP, pasemos las posiciones y velocidades a un sistema rotante que tenga la misma frecuencia angular que el movimiento orbital medio de las partículas más masivas. Para ello usemos como aproximación la ley armónica $n^2 a^3=\mu$ que aplicada al sistema formado por las dos partículas masivas sería:
$$ n=\sqrt{\frac{G(m_1+m_2)}{a^3}} $$donde $a$ es la distancia promedio entre las partículas. Numéricamente:
#Distancia entre las partículas 1 y 2
from numpy.linalg import norm
r12=norm(rps[1]-rps[2],axis=1)
#Promedio
a=r12.mean()
Con ello la velocidad angular promedio será:
from numpy import sqrt
m1=sistema[1]["m"]
m2=sistema[2]["m"]
n=sqrt((m1+m2)/a**3)
Rotemos ahora todas las posiciones y velocidades de las partículas usando los mismos métodos de la Sección Un ejemplo numérico de sistemas rotantes:
#Velocidad angular
w=n
from numpy import array
omega=array([0,0,w])
#Vectores en el sistema rotado
from numpy import zeros_like
rps_rot=zeros_like(rps)
vps_rot=zeros_like(vps)
for i in range(Nt):
#Matriz de rotación
from spiceypy import rotate,mxv
R=rotate(w*ts[i],3)
#Rota las posiciones y velocidades de cada partícula
for n in range(3):
rps_rot[n,i]=mxv(R,rps[n,i])
#v' = v + w x r
from spiceypy import vcrss
vps_rot[n,i]=mxv(R,vps[n,i]-vcrss(omega,rps[n,i]))
Veamos como se ve la dinámica en el sistema rotante después de esta transformación:
from pymcel.export import plot_ncuerpos_3d
fig=plot_ncuerpos_3d(rps_rot,vps_rot,marker='o',ms=3);
Como era de esperarse, las partículas más masivas parecen quietas, mientras que la partícula de prueba describe una trayectoria compleja en el sistema.
La prueba de fuego de nuestra teoría consiste en averiguar si incluso en un sistema que no satisface estrictamente las condiciones del CRTBP (la partícula de prueba tiene una masa distinta de cero, la órbita de las partículas masivas no es exactamente circular debido a la perturbación de la partícula), la cuadratura de Jacobi sigue siendo aproximadamente constante. Para ello calculemos el valor de $C_J$ en este sistema usando la fórmula original en la Ec. (constante_Jacobi):
#Posiciones de la partícula de prueba
r=rps_rot[0]
#Posiciones y velocidades relativas a las partículas masivas
from numpy.linalg import norm
r1=norm(rps_rot[0]-rps_rot[1],axis=1)
r2=norm(rps_rot[0]-rps_rot[2],axis=1)
#Rapideces
v=norm(vps_rot[0],axis=1)
#Parametros gravitacionales
mu1=m1
mu2=m2
#Cuadratura de Jacobi
CJ=2*mu1/r1+2*mu2/r2+w**2*(r[:,0]**2+r[:,1]**2)-v**2
Veamos un gráfico de $C_J$ como función del tiempo:
fig=plt.figure()
ax=fig.gca()
ax.plot(ts,CJ)
ax.set_xlabel("$t$ [u.c.]");
ax.set_ylabel("$C_J$ [u.c.]");
La primera observación de interés que podemos hacer sobre la Figura (code:constante_jacobi_sistema_real) es que efectivamente se verifica que incluso en un sistema real, el valor de $C_J$ se mantiene aproximadamente constante, lo que nos da confianza frente a la aplicación del formalismo aproximado desarrollado para el CRTBP en casos de sistemas reales.
El valor de la constante en este caso ($C_J\approx 11.3$) no coincide con las expectativas de las secciones anteriores, en las que vimos que en condiciones normales esta cantidad está entre 2 y 5. Sin embargo debemos tener en cuenta que todas las cantidades en la simulación del sistema real que vimos en esta sección están expresadas en un sistema de unidades canónicas que no coincide necesariamente con el del CRTBP.
Indice | Previo: Problema3Cuerpos.CRTBP.Numerico | Siguiente: Problema3Cuerpos.RegionesExclusion