Pasado y futuro del Sistema Solar: introducción a la Dinámica Molecular

1  Introducción.

Laplace.jpg Hace 200 años, Pierre Simon, Marqués de Laplace, comenzó a pensar sobre el origen del sistema solar. Su hipótesis fue presentada en 1796 en Exposition du Système du Monde. De acuerdo con la teoría de Laplace, el Sistema Solar evolucionó a partir de un gas incandescente que giraba alrededor de un eje. Según se iba enfriando, el gas se contraía, lo que originaba un rotación más rápida debido a la conservación del momento angular. En algún momento, la parte exterior del gas se separó formando anillos que posteriormente se condensaron formando planetas. La parte interior colapsó en un sólo cuerpo: el Sol. La idea detrás de esta hipótesis no se debe sólo a Laplace; ya fue propuesta antes en una forma menos elaborada por Emanuel Swedenborg e Immanuel Kant. En cualquier, caso Laplace desarrolló la idea cuyos detalles técnicos plasmó en su gran obra Traite de Mécanique Céleste (5 volúmenes publicados entre 1799 y 1825). Esta obra es la culminación del punto de vista newtoniano sobre la gravitación. De hecho, en ella Laplace mostró que el Sistema Solar era estable bajo la acción de pequeñas perturbaciones. Su trabajo demostraba la ausencia de necesidad de intervención divina en la formación y estabilidad del Sistema Solar. Se dice que Napoleón le comentó este hecho y Laplace le respondió: no necesito esa hipótesis.
Esta imágen genérica de Laplace es la que actualmente aceptamos como la teoría más viable para explicar la formación del Sistema Solar. En su versión más moderna se conoce como la hipótesis nebular. Esta hipótesis sostiene que el Sistema Solar se formó a partir del colapso gravitacional de un fragmento de nube gigante (de varios años luz de diámetro probablemente). Estudios recientes apuntan a que cierto número de explosiones de supernovas ocurrieron cerca de esta nube. Las ondas de choque procedentes de estas supernovas iniciaron la formación del Sistema Solar al crear fluctuaciones de densidad dentro de la nube, con el consiguiente colapso de las regiones de mayor densidad. Una de estas regiones de gas colapsado (llamadas nubes presolares) formaría lo que más tarde sería el Sistema Solar.
M42proplyds.jpg Debido a la conservación del momento angular, la nube presolar gira más rápido conforme se colapsa. El material del interior de la nube comienza a condensarse, convirtiendo energía potencial en cinética e incrementando la frecuencia de colisión de las partículas, lo que inicia reacciones de fusión. De esta manera, la zona central en la que se acumula la mayor parte de la masa se calienta mucho más que la región circundante. La competición entre la gravedad, la presión del gas, los campos magnéticos y la rotación da lugar a la contracción de la nube en un disco protoplanetario de unas 200 Unidades Astronómicas de diámetro y forma una protoestrella densa y caliente en el centro. En un período de 50 millones de años, la temperatura y la densidad en el núcleo de la protoestrella aumentan de manera sostenida hasta que el hidrógeno comienza a fusionarse, lo que marca la entrada del Sol en su primera fase de existencia, la secuencia principal.
Phobos-planetesimal.jpg Los diferentes planetas se forman entonces a partir del disco protoplanetario remanente que orbita alrededor del Sol recién formado, a través de un mecanismo conocido como acreción. Las partículas del gas colisionan entre sí de manera inelástica formando agregados cada vez mayores. En cuanto alcanzan tamaños suficientemente grandes, su atracción gravitatoria favorece su acumulación, dando lugar a objetos con un de tamaño que ronda los 5 km de diámetro que se denominan planetesimales. 'Estos incrementan gradualmente su tamaño con más colisiones, creciendo a un ritmo de varios centímetros al año a lo largo de varios millones de años. Los cometas y los asteroides se cree que son planetesimales sobrantes durante la formación del Sistema Solar.
La zona interior del Sistema Solar en formación está demasiado caliente para que moléculas volátiles tales como las de agua o metano se condensen, de manera que los planetesimales formados en esta región sólo pueden estar compuestos de elementos con puntos de fusión altos (como metales y silicatos). Estos cuerpos rocosos darán lugar a los planetas terrestres (Mercurio, Venus, la Tierra y Marte). Los compuestos que los forman son poco abundantes en la nube molecular inicial ( ∼ 0.6 %), por lo que los planetas terrestres no pueden crecer demasiado. Por otra parte, los planetas gaseosos gigantes (Júpiter, Saturno, Urano y Neptuno) se forman más lejos, más allá de la línea de congelación (el punto entre las órbitas de Marte y Júpiter donde el material está lo suficientemente frío como para que los compuestos volátiles congelados permanezcan en estado sólido). Estos materiales helados son mucho más abundantes que los metales y silicatos de los planetas terrestres, lo que permite que los planeta Jovianos crezcan lo suficiente como para capturar hidrógeno y helio, los elementos más ligeros y abundantes. Así, los planetesimales más allá de la línea de congelación pueden acumular hasta 4 veces la masa de la Tierra en 3 millones de años, y hoy en día los cuatro gigantes gaseosos acumulan el 99% de la masa total del Sistema Solar. La formación como tal del Sistema Solar termina cuando, después de varios millones de años (entre 3 y 10), los vientos solares de la joven estrella central han barrido todo el gas y el polvo sobrantes del disco protoplanetario, dispersándolos por el espacio interestelar.
Lhborbits.png
El futuro del sistema solar puede estudiarse a partir de simulaciones numéricas. Como ejemplo, se puede consultar el artículo de J. Sackmann, A.I. Boothroyd y K. Kraemer Our Sun. III. Present and Future en The Astrophysical Journal 418 457-468 (1993). Allí se nos realata como el Sol se convertirá en gigante roja, y despué de varios episodios de inestabilidad y pérdida de masa, acabará su vida como enana blanca.

2  Modelo

Nuestro sistema está compuesto por el Sol y los nueve planetas, que suponemos sin estructura interna y caracterizados por su masa. La única interacción presente es la gravitatoria, y asumimos que el movimiento planetario se realiza en un plano. Usando datos reales para la excentricidad de las órbitas, las distancias y las masas de los diferentes planetas, podemos dar una condición inicial realista para el Sistema Solar y así estudiar su dinámica y estabilidad frente a perturbaciones.

3  Problemas

Obligatorio: simular la dinámica del sistema solar

Diseñar y escribir el programa que simula el comportamiento del Sistema Solar. Comprobar la validez de la simulación midiendo los períodos de rotación de los diferentes planetas y comparándolos con los reales. Mostrar la órbita relativa respecto a la Tierra del resto de planetas, en función del tiempo.

Voluntario A: futuro del sistemas solar

A partir de las condiciones iniciales actuales de los planetas (posición, velocidad y excentricidad de las órbitas), estudiar su comportamiento futuro sabiendo que el Sol irá sufriendo cambios estructurales. Realizar la simulación siguiendo el siguiente proceso:

Voluntario B: formación de los planetas a partir de planetesimales

Partimos de la condición inicial de un sol masivo y un gran número de pequeños cuerpos orbitando alrededor de éste con diferentes energías y excentricidades. Los pequeños planetesimales interaccionan entre sí gravitatoriamente. La formación de planetas se va realizando mediante choques totalmente inelásticos entre planetesimales. Dada una condición inicial de planetesimales, estudiar la distribución de tamaños de los planetas finales. Puesto que la colisión entre planetesimales es inelástica, parte de su energía cinética antes de su colisión se convierte en calor. Mientras que en ese mismo proceso las energías caloríficas que tenían se suman. Seguir la pista de esas energías y obtener las energías internas de la distribución de planetas finales. Comparar las propiedades obtenidas con las del sistema solar en la actualidad.

4  Movimiento de una partícula bajo la acción de fuerzas centrales

Las ecuaciones del movimiento para un planeta de masa m que sufre la atracción gravitatoria del Sol (de masa M) son, en 3 dimensiones,
m
d2
→
r
 
dt2
 =
→
F
 
 (
→
r
 
 )        con       
→
F
 
 (
→
r
 
 )=−G Mm
r3
→
r
 
 ,
(1)
donde G es la constante de gravitación universal. El momento angular se define como →L=→r×→p, o equivalentemente Li=ϵijkrjpk, donde ϵijk es el símbolo de Levi-Civita (o de las permutaciones), tal que ϵ123=1, ϵ213=−1, etc. La conservación del momento angular implica que
d
→
L
 
dt
 =
→
v
 
 ×
→
p
 
 +
→
r
 
 ×
→
F
 
 =0.
(2)
Las condiciones iniciales fijan entonces →L, que es perpendicular al plano de movimiento del planeta. Por tanto, tenemos un problema puramente bidimensional.
El lagrangiano del planeta, L=[1/2]m(· x2+· y2 )+G[Mm/r], se puede escribir en coordenadas polares como L=[1/2]m(· r2+r2· θ2)+G[Mm/r], con lo que las ecuaciones del movimiento quedan
d
dt
 ⎛
⎜
⎝ 		∂L
∂
⋅
r
 
⎞
⎟
⎠ 		− ∂L
∂r
 =0    
⇒
    m
⋅⋅
r
 
 =−G Mm
r2
 + l2
mr3
  ,
d
dt
 ⎛
⎜
⎝ 		∂L
∂
⋅
θ
 
⎞
⎟
⎠ 		− ∂L
∂θ
 =0    
⇒
    d
dt
 ⎛
⎝ 		mr2
⋅
θ
 
⎞
⎠ 		=0  .
La solución de la segunda ecuación es trivial, mr2· θ=l=cte. Es importante notar que parametrizar nuestro problema en función del tiempo, (r(t),θ(t)), es equivalente a parametrizarlo en función del ángulo θ, (r(θ),t(θ)), pues ambas descripciones se relacionan trivialmente usando la segunda ecuación. Así, podemos reescribir la primera ecuación en función de θ,
⎛
⎝ 		d
dt
 = dθ
dt
 d
dθ
 = l
mr2
 d
dθ
⎞
⎠ 		      ⇒        l
r2
 d
dθ
 ⎛
⎝ 		l
mr2
 dr
dθ
⎞
⎠ 		− l2
mr3
 =−G Mm
r2
 ,
(3)
Si ahora hacemos primero un cambio de variable u=1/r, que lleva a,
d2u
dθ2
 +u=G Mm2
l2
 ,
(4)
y definimos una nueva variable y tal que u=y+G[(Mm2)/(l2)], obtenemos la siguiente ecuación diferencial
d2y
dθ2
 +y=0,
(5)
cuya solución es sencilla,
1
r
 =G Mm2
l2
 [1+εcos(θ−θ′)],
(6)
donde ε ≥ 0 y θ′ son constantes que dependen de las condiciones iniciales. En particular, ε se puede relacionar con la energía total del sistema, E, que también es una magnitud conservada. Resulta sencillo demostrar que
ε =   ⎛
√
1+ 2El2
G2M2m3
 
 .
(7)
La solución obtenida es una curva cónica genérica: Finalmente, dejamos como ejercicio dibujar para valores diferentes de ε la curva (x,y): r=1/(1+εcosθ), x=rcosθ, y=rsinθ.

5  Método numérico

Nuestro sistema se compone en general de N discos de diferentes masas mi, i ∈ [1,N]. Los cuerpos interaccionan entre si mediante fuerzas gravitatorias newtonianas. La ecuación de evolución del cuerpo i es:
mi
d2
→
r
 
i 
dt2
 =−Gmi
∑
j ≠ i 
mj (
→
r
 
i 
 −
→
r
 
j 
 )
|
→
r
 
i 
 −
→
r
 
j 
 |3
(8)
Existen métodos muy precisos para resolver numéricamente este conjunto de N ecuaciones diferenciales. Sin embargo, resultan muy costosos computacionalmete, en especial cuando N es muy grande. En estos casos, uno puede renunciar en parte a la precisión del cálculo de cada una de las trayectorias de las partículas por la rápidez del cálculo. Este hecho, per se, justifica el uso de algoritmos simplificados. El uso de estos algoritmos simplificados sólo tendrá sentido cuando nos interesen propiedades globales o macroscópicas del sistema, puesto que en ese caso lo que importa no es el detalle del camino seguido sino sus propiedades promedio. En cualquier caso, el algoritmo utilizado debe cumplir ciertas propiedades de estabilidad. Por ejemplo, si la dinámica conserva la energía y el momento angular, es importante y necesario que el algoritmo sea diseñado para que los conserve.
Desde el punto de vista computacional los algoritmos más eficaces son los llamados algoritmos de un paso, esto es, dadas las posiciones y momentos de la partículas en un instante dado, podemos a partir de ellas calcular sus respectivas posiciones y momentos en un paso posterior. Veamos como construir uno de estos algoritmos, el llamado algoritmo de Verlet en velocidad. Este algoritmo se basa en el desarrollo de Taylor de las posición y velocidad de la partícula,
→
r
 
 (t+h)=
→
r
 
 (t) + h
→
v
 
 (t) + h2
2
→
a
 
 (t) + O(h3)  ,
→
v
 
 (t+h)=
→
v
 
 (t) + h
→
a
 
 (t) + h2
2
→
b
 
 (t) + O(h3)  ,
donde →b(t)=[(d →a)/d t](t), y h es el paso temporal de integración del algoritmo. La primera derivada de la aceleración se puede escribir como
→
b
 
 (t)=
→
a
 
 (t+h) −
→
a
 
 (t)
h
 + O(h).
(9)
Usando esta expresión para →b(t) en la ecuación para la velocidad, y truncando ambos desarrollos a orden h3, obtenemos
→
r
 
 (t+h) ≈
→
r
 
 (t) + h
→
v
 
 (t) + h2
2
→
a
 
 (t)  ,
(10)
→
v
 
 (t+h) ≈
→
v
 
 (t) + h
2
 ⎡
⎣
→
a
 
 (t) +
→
a
 
 (t+h)⎤
⎦ 		.
(11)
El algoritmo de Verlet de un solo paso queda totalmente definido por las ecuaciones anteriores. Este algoritmo conserva la energía y el momento angular de nuestro sistema. Sin embargo, es importante notar que esto no significa que cuando introduzcamos el algoritmo en el ordenador se conserven dichas cantidades, puesto que los errores de redondeo se irán acumulando y afectarán a dichos valores. Sin embargo, el conocimiento en todo momento de su evolución nos dará idea de cuan lejos evolutivamente podemos llegar manteniendo el sentido físico de los resultados. En resumen, el algoritmo de Verlet consiste en : Como ejercicio, podemos aplicar el algoritmo de Verlet a la resolución del movimiento de un péndulo real, m[(d2θ)/(dt2)]=−mgsinθ.
Para aplicar el algoritmo de Verlet en nuestro caso, hemos de tener en cuenta que vamos a trabajar con números muy grandes y muy pequeños simultáneamente (por ejemplo, G=6.67×10−11  N  m2/kg2 y la Masa solar es Ms=1.99×1030  kg). Por tanto resulta conveniente reescalar las ecuaciones del movimiento. En particular, podemos elegir el siguiente reescalamiento:
r′= r
c
        ;        t′=⎡
⎣ 		G Ms
c3
⎤
⎦ 		1/2

 
 t     ;       m′= m
Ms
 ,
(12)
donde elegimos c=1.496×1011 m, que no es más que la distancia entre la tierra y el Sol. En estas unidades, Neptuno se encuentra aproximadamente a 30 unidades c de distancia al Sol, y la unidad temporal t′ corresponde a unos 58.1 dias. Las ecuaciones del movimiento se convierten ahora en
d2
→
ri′
 
dt′2
 =−
∑
j ≠ i 
mj′(
→
ri′
 
 −
→
rj′
 
 )
|
→
ri′
 
 −
→
rj′
 
 |3
(13)
Por último, para el caso del problema de la formación de los planetas a partir de los planetesimales suponemos que la agregación de materia se produce mediante colisiones totalmente inelásticas. Así, el algoritmo de colisión entre las partículas i y j lo definimos de la siguiente manera: Cabe destacar que al introducir el choque inelástico, la energía cinética no se conserva, generándose calor en el proceso. Por tanto, este calor disipado hay que tenerlo en cuenta si queremos utilizar la energía total del sistema para controlar los errores de redondeo.


File translated from TEX by TTH, version 4.03.
On 3 Feb 2014, 13:22.