Indice | Previo: Problema3Cuerpos.RegionesExclusion | Siguiente: Problema3Cuerpos.PuntosEquilibrioLagrange
En la sección anterior comprobamos teorica y numéricamente que las superficies de cero velocidad actúan como espejos contra los que se reflejan las trayectorias de partículas de prueba en el CRTBP. La razón de ello es que, por definición, cuando la partícula alcanza estas superficies imaginarias se detiene. Esto es consistente con el hecho además de que más allá de ellas se encuentran regiones que no pueden ser atravesadas por la partícula. Pero surge la pregunta de ¿por qué una vez la partícula llega a la superficie de cero velocidad y se detiene, vuelve después a ganar velocidad para desprenderse de ella? ¿en qué dirección se acelera después de alcanzar la superficie de cero velocidad?.
Para resolver estas preguntas podemos apelar a la ecuación de movimiento de la partícula que escribimos en la Sección El problema circular restringido de los tres cuerpos:
\begin{equation} \label{eq:crtbp_edm} \ddot{\vec{r}} = -(1-\alpha)\frac{\vec{r}_{1}}{r^3_{1}} - \alpha\frac{\vec{r}_{2}}{r^3_{2}} - \hat{e}_z\times(\hat{e}_z\times \vec{r}) -2\hat{e}_z\times \dot{\vec{r}} \end{equation}Para el asunto particular que nos ocupa, nos interesa la forma que adopta esta ecuación de movimiento unicamente cuando la partícula esta en reposo, es decir se encuentra sobre una superficie de cero velocidad. En estas condiciones el término de la aceleración de Coriolís desaparece y la ecuación se convierte en:
\begin{equation} \label{eq:edm_crtbp_cero_velocidad} \ddot{\vec{r}} = -(1-\alpha)\frac{\vec{r}_{1}}{r^3_{1}} - \alpha\frac{\vec{r}_{2}}{r^3_{2}} - \hat{e}_z\times(\hat{e}_z\times \vec{r}) \end{equation}donde el subíndice 0 en $\ddot{\vec{r}}_0$ lo hemos usado para indicar que esta ecuación solo aplica para partículas que se encuentran en reposo en el sistema rotante.
Nota: Reposo no implica equilibrio. Aunque para este punto el lector debe seguramente haber reflexionado ampliamente sobre el malentendido común de que el reposo es lo mismo a la ausencia total de movimiento, no sobra, dada la importancia en este punto del libro, recordar esta distinción conceptual fundamental.
Por reposo entendemos el estado de movimiento caracterizado exclusivamente por la condición $\vec v=\vec o$. Ahora bien, el estado completo de movimiento de un cuerpo no solo esta determinado por su velocidad. Otras propiedades cinemáticas también son importantes, siendo la más relevante para esta discusión la aceleración $\vec a$.
Si una partícula esta en reposo, $\vec v=\vec o$ pero también tiene una aceleración nula $\vec a=\vec o$, entonces, en virtud del teorema de inercia la partícula se mantendrá en reposo hasta que una fuerza la saque de esa situación. A esta condición la llamamos equilibrio.
Sin embargo hay estados de movimiento en los que $\vec v=\vec o$ pero $\vec a\neq\vec o$ (reposo pero no equilibrio). Este es el caso por ejemplo de una partícula lanzada hacia arriba en un campo gravitacional casi uniforme (el campo gravitacional muy cerca a la superficie de la Tierra). En la parte alta de su trayectoria la partícula alcanza el reposo (se detiene) pero su aceleración (la aceleración de la gravedad) es distinta de cero (no hay equilibrio); un instante de tiempo después, la partícula estará nuevamente en movimiento hacia abajo.
Esto es justamente lo que pasa cuando en el CRTBP la partícula de prueba alcanza una superficie de cero velocidad: al hacerlo, se detiene, pero si la aceleración en aquel lugar es distinta de cero (dependiendo de lo que dicte justamente la Ec. edm_crtbp_cero_velocidad) entonces un momento después estará otra vez en movimiento.
Es fácil mostrar (ver Problemas al final del capítulo) que la ecuación de movimiento en (edm_crtbp_cero_velocidad) puede escribirse en la forma:
\begin{equation} \label{eq:edm_crtbp_reposo} \ddot{\vec{r_0}} = - \vec\nabla \left[-\frac{(1-\alpha)}{r_{1}}-\frac{\alpha}{r_{2}}-\frac{1}{2}(x^2 + y^2)\right] \end{equation}es decir, en el sistema de referencia rotante, una partícula en reposo experimentará una aceleración producto de un potencial escalar diferente del simple potencial gravitacional newtoniano. A este potencial lo llamaremos aquí el potencial modificado del CRTBP y tiene la forma:
\begin{equation} \label{eq:Vmod_crtbp} V_\mathrm{mod} = -\frac{(1-\alpha)}{r_{1}} - \frac{\alpha}{r_{2}} - \frac{1}{2}(x^2 + y^2) \end{equation}En términos del potencial modificado las ecuación de movimiento del sistema (Ec. crtbp_edm) incluso cuando la velocidad de la partícula de prueba no es 0, se pueden escribir como:
\begin{eqnarray} \label{eq:crtbp_edm_Vmod} \ddot{x}-2\dot y & = & -\frac{\partial V_\mathrm{mod}}{\partial x} \\ \nonumber \ddot{y}+2\dot x & = & -\frac{\partial V_\mathrm{mod}}{\partial y} \\ \nonumber \ddot{z} & = & -\frac{\partial V_\mathrm{mod}}{\partial z} \\ \end{eqnarray}¿Cómo luce este potencial modificado?. Naturalmente el potencial $V_\mathrm{mod}$ depende de todas $(x,y,z)$ y visualizarlo sería muy complejo. Sin embargo podemos, como hicimos en secciones anteriores, concentrarnos en el valor del potencial sobre el plano $x-y$ (el plano de la órbita de las partículas masivas). Visualizarlo allí puede darnos alguna intuición de su estructura 3d.
En el algoritmo a continuación se calcula el valor del potencial modificado en una malla coordenada como la que hemos usado a lo largo de este capítulo y se visualiza en el espacio de tres dimensiones:
#Parámetros del sistema
alfa=0.3
#Malla coordenada
from numpy import meshgrid,zeros_like,linspace
rango=1.5
NG=100
X,Y=meshgrid(linspace(-rango,rango,NG),
linspace(-rango,rango,NG))
Z=zeros_like(X)
#Factor de suavizado (ver nota)
eps=0.8*alfa
#Distancia relativa
from numpy import sqrt
r1=sqrt((X+alfa)**2+Y**2+Z**2+eps**2)
r2=sqrt((X-1+alfa)**2+Y**2+Z**2+eps**2)
#Calcula el potencial
Vmod=-(1-alfa)/r1-alfa/r2-0.5*(X**2+Y**2)
Nótese que el valor absoluto del potencial es del orden de la mitad de la constante de Jacobi. Esto era de esperarse del hecho de que las expresiones que definen ambas cantidades son casi idénticas (Ecs. Vmod_crtbp y constante_Jacobi).
Nota: el factor de suavizado. En el cálculo del potencial hemos usado un artefacto matemático, que no procede de la teoría del CRTBP. Puesto que el valor del potencial es teóricamente infinito cuando $r_1=0$ o $r_2=0$ hemos introducido en las expresiones para el potencial una cantidad ficticia eps
proporcional al parámetro $\alpha$ del sistema (relacionado con la masa de las partículas), y que hace que el potencial mantenga un valor finito cerca al centro de las partículas. La introducción de esta cantidad artificial implica también que los valores calculados del potencial no son los correctos a una cierta distancia de las partículas (típicamente similar al valor de eps
).
Este artefacto numérico facilita la representación gráfica del potencial, pero también es importante en otras aplicaciones de la mecánica celeste. Así por ejemplo cuando se simula numéricamente la interacción gravitacional de partículas en sistemas de N cuerpos tales como galaxias o nubes de gas, el facto
eps
recibe el nombre de factorde suavizado y se usa para simular de alguna manera el hecho de que las *partículas de materia** tienen en realidad una dimensión finita (no son realmente puntuales).
Una gráfica del potencial modificado sobre el plano $x-y$ en el CRTBP se puede crear con el algoritmo a continuación:
%matplotlib inline
#Grafico
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig=plt.figure()
ax=fig.add_subplot(111, projection='3d')
ax.plot_surface(X,Y,Vmod,cmap="jet",rstride=2,cstride=2)
ax.set_zlim((-3,0))
ax.view_init(30,-15)
ax.set_xlabel("$x$")
ax.set_ylabel("$y$")
ax.set_zlabel("$V_\mathrm{mod}$")
ax.set_title(f"Potencial modificado en el CRTBP con $\\alpha={alfa}$")
fig.tight_layout()
Otra manera de visualizar el mismo potencial es usando un gráfico de contornos:
fig=plt.figure()
ax=fig.gca()
#Contornos
ax.contour(X,Y,Vmod,levels=15,colors='k',linestyles='solid',alpha=0.3)
c=ax.contourf(X,Y,Vmod,levels=100,cmap="jet")
fig.colorbar(c)
#Decoración
ax.set_xlabel("$x$")
ax.set_ylabel("$y$")
ax.set_title(f"Potencial modificado en el CRTBP con $\\alpha={alfa}$")
fig.tight_layout()
#Puntos de interés
ax.plot([0],[0],'ko')
ax.text(0,0,'A',fontsize=18)
ax.plot([1],[1],'ko')
ax.text(1,1,'B',fontsize=18)
ax.plot([0.8*alfa],[-1+0.7*alfa],'ko')
ax.text(0.8*alfa,-1+0.7*alfa,'C',fontsize=18)
ax.plot([1+0.5*alfa],[0],'ko')
ax.text(1+0.5*alfa,0,'D',fontsize=18);
¿Cómo podemos interpretar la estructura del potencial modificado?. Como es corriente en mecánica, el potencial asociado a una fuerza, dicta la manera como la partícula se acelera. En particular, la dirección de la aceleración es perpendicular a las curvas de equipotencial (dirección del gradiente) y va en el sentido en el que el potencial disminuye.
En la gráfica de equipotenciales de la Figura (code:Vmod_contornos) hemos resaltado algunos puntos que merecen atención y que nos permiten interpretar mejor el significado de $V_\mathrm{mod}$:
Indice | Previo: Problema3Cuerpos.RegionesExclusion | Siguiente: Problema3Cuerpos.PuntosEquilibrioLagrange