Daedalus

Un integrador simpléctico híbrido para la dinámica de planetesimales.

Acerca del programa

Daedalus e Icarus por Jean Bouchet

Daedalus e Icarus volando fuera del laberinto, por Jean Bouchet (1476-1558)

DAEDALUS es un programa diseñado para resolver las ecuaciones de movimiento de un conjunto de cuerpos masivos gravitatoriamente interactuantes sujetos a la atracción de un cuerpo central de masa dominante. Típicamente, un disco protoplanetario orbitando una estrella central.

El programa implementa un algoritmo simpléctico híbrido capaz de resolver encuentros entre cuerpos y donde la evaluación de las fuerzas de interacción es efectuada implementando la estrategia de árbol jerárquico octal.

DAEDALUS está escrito en FORTRAN 77, con algunas extensiones que no forman parte del estandard (sentencias END DO, DO WHILE, y algunas subrutinas provistas por el compilador gfortran, soportadas también por el compilador ifort).

Descarga

Test

Con el fin de testear el desempeño del código recurrimos al escenario del crecimiento en fuga (runaway) de la acreción de planetesimales, situación modelada numéricamente por Kokubo & Ida (1996). La simulación consiste inicialmente de planetesimales de igual masa gr (con densidades de gr cm) distribuidos aleatoriamente en una región anular alrededor de UA con un ancho de 0.2 UA. Las excentricidades e inclinaciones de las órbitas osculatrices iniciales de los planetesimales están dadas por distribuciones gaussianas con dispersiones , siendo el radio de Hill de los cuerpos.

Situación no colisional

En el primer conjunto de experimentos numéricos, los cuerpos son tratados como masas puntuales, con lo cual las colisiones son ignoradas. Esto nos permite testear la precisión global del código (medida por la conservación de la energía total del sistema), la precisión del código de árbol en el cómputo de las interacciones entre los planetesimales, y la eficiencia en la detección de encuentros.

Histograma de errores relativos del potencial

Fig. 1. Histogramas de los errores relativos en los potenciales de los planetesimales (sin incluir el cuerpo central) en el paso inicial, para diversos valores del parámetro de expansión del árbol.

Para testear, en un dado paso, la precisión del árbol en el cómputo de las perturbaciones de los planetesimales (sin incluir al cuerpo central) podemos considerar el error cometido en la evaluación del potencial de un cuerpo por el algoritmo respecto de su valor exacto dado por la evaluación directa. Así, el error relativo en el potencial de un cuerpo está dado por, donde y son los potenciales para el cuerpo i-ésimo, evaluados con el método del árbol y el método directo partícula-partícula (PP), respectivamente. La Fig. 1 muestra, en el paso inicial, los histogramas de los errores relativos (esto es, el número de cuerpos para un dado rango del error) para diversos parámetros de expansión del árbol.

Una estimación globlal del error en el cálculo del potencial puede ser definida como sigue (Pfalzner & Gibbon, 1996): La Fig. 2 ilustra el comportamiento de tal error global con el parámetro de expansión del árbol para el primer paso de nuestras simulaciones. Por otra parte, la Fig. 3 da cuenta del tiempo de CPU invertido en los cómputos para un dado valor del parámetro de expansión (en dicha figura el tiempo es normalizado con respecto al tiempo de cómputo para el valor cero del parámetro de expansión).

Del conjunto de figuras, puede observarse que, para un parámetro de expansión tendiendo a cero, todos los cuerpos tienen un error relativo cercano a cero, en tanto que, para un parámetro de expansión tendiendo a la unidad, existe una proporción grande de cuerpos con errores relativamente grandes, y que, para parámetros de expansión mayores que la unidad, comenzamos a perder gran precisión para ganar poca velocidad de cómputo. El punto a destacar es que, para un valor del parámetro de expansión de 0.7 obtenemos una buena precisión conjuntamente con una buena eficiencia del método. Esto justifica que en las siguientes simulaciones de prueba escojamos dicho valor para el parámetro de expansión del árbol.

Fig. 2. Error global en el potencial en función del parámetro de expansión del árbol.

Fig. 3. Tiempo de cómputo de CPU en función del parámetro de expansión del árbol.

Con un paso de integración de 0.005 años (y un pámetro de expansión de 0.7) consideramos el comportamiento temporal de la energía del sistema a lo largo de 10 000 años. Si nuestro código es realmente simpléctico, el error relativo en la energía total del sistema no debería exhibir un comportamiento secular. Efectivamente, como se ilustra en la Fig. 4, si bien las variaciones de la energía sobre períodos cortos de tiempo son relativamente grandes, no resulta evidente un comportamiento secular de las mismas. De hecho, la dispersión de los valores respecto de la media es .

Fig. 4. Error relativo en la energía total del sistema en función del tiempo (expresado en años).

Situación colisional

Fig. 5. Evolución temporal del número de cuerpos en el disco. Abscisas: tiempo en años, ordenadas: número de cuerpos (escala logaritmica sobre ambos ejes).

Habiendo quedado satisfechos con los resultados de las simulaciones anteriores, consideramos ahora una nueva simulación, pero permitiendo esta vez la posibilidad de colisiones entre los planetesimales. En esta nueva simulación tal como lo hicieran Kokubo & Ida (1996), el radio inicial de los planetesimales es aumentado en un factor 5 respecto del valor realista que se obtiene asumiendo que la densidad de los mismos es 2 g cm. Como muestran Kokubo & Ida, ésto solo cambia el tiempo de escala de acreción, pero no el modo de crecimiento de los planetesimales. Con las mismas condiciones iniciales que antes, el sistema es integrado sobre 20000 años con un paso de integración de 0.01 años (tomando el parámetro de expansión del árbol igual a 0.7). Los resultados obtenidos son ilustrados en el conjunto de figuras 5,6 y 7, las cuales deberían compararse con las correspondientes al trabajo de Kokubo & Ida. La Fig. 5 muestra la evolución del número total de cuerpos en función del tiempo. Este decrece de 2000 a 833 en los 20000 años (en la simulación de Kokubo & Ida, el número final de cuerpos es 715). La Fig. 6 muestra la evolución temporal de la mayor masa en cada instante de tiempo y la masa media de los planetesimales (exceptuando la masa máxima). Al final de la integración, la masa final del cuerpo más masivo es unas 93 veces el valor de su masa inicial. Este valor difiere en un factor 2 del resultado de Kokubo & Ida, la cual se debe a un evento colisional aleatorio al final de la simulación. Aun así, la razón de la masa máxima a la masa media crece monótonamente con el tiempo a una tasa similar a la obtenida por Kokubo & Ida. La Fig 7 muestra la evolución temporal de las medias cuadráticas de la excentricidad e inclinación de las órbitas osculatrices de los planetesimales. Ambas se incrementan monótonamente con el tiempo y puede apreciarse que se reproduce la evolución esperada: la razón de una a la otra se mantiene aproximadamente igual a 2.

Fig. 6. Evolución temporal de la masa máxima (curva en rojo) y de la masa media (curva en verde). Abscisas: tiempo en años, ordenadas: masa en gr (escala logaritmica sobre ambos ejes).

Fig. 7. Evolución temporal de la media cuadrática de la excentricidad (curva en rojo) y de la inclinación (curva en verde) de la población del disco. Abscisas: tiempo en años (escala logaritmica sobre ambos ejes).

La Fig. 8 muestra cuatro instantáneas del sistema en diferentes tiempos de evolución. El tamaño de los símbolos indica tres rangos de masas: los puntos corresponden a cuerpos cuyas masas se mantienen iguales a la masa inicial, los círculos de tamaño intermedio a cuerpos con masas hasta diez veces la masa inicial, y los de mayor tamaño a cuerpos con masas superiores a diez veces la masa inicial.


a) t = 0 años, N = 2000.


b) t = 5000 años, N = 1285.


c) t = 10000 años, N = 1069.


d) t = 20000 años, N = 833.

Fig. 8. Evolución temporal del sistema en el plano x-y.

Ejemplos

Fig 9. Evolución temporal del semieje mayor de un planeta de 10 masas terrestres localizado inicialmente a 5 UA en un disco de unas 150 masas terrestres. Abscisas: tiempo en años.

Cionco & Brunini (2002) analizan la interacción planeta-disco desde un punto de vista resonante, efectuando un conjunto de simulaciones numéricas de la migración orbital de un planeta masivo (en el rango de 10 masas terrestres a una masa de Júpiter) sumergido en un disco de planetesimales. A partir de dichas simulaciones, los autores determinan la variación del semieje mayor del planeta y de la misma derivan la velocidad de migración y el torque ejercido sobre el planeta. La Fig. 9 muestra la evolución temporal del semieje para un planeta de 10 masas terrestres obtenida con nuestro código para las mismas condiciones del trabajo mencionado. La tasa de migración derivada es compatible con los resultados de dicho trabajo.

Gomes et al. (2004) muestran que, en el caso de nuestro Sistema Solar, dependiendo de la masa del disco de planetesimales, Neptuno podría haber tenido una migracción orbital amortiguada, limitada a unas pocas unidades astronómicas, o bien una migracción forzada que lo lleva rápidamente hacia el borde del disco. La Fig. 10 ilustra esta última situación en una simulación obtenida con nuestro código bajo condiciones iniciales semejantes a la del trabajo mencionado.

Fig 10. Evolución temporal del semieje mayor de los cuatro planetas gigantes en un disco de unas 200 masas terrestres que se extienden inicialmente entre las 20 y 50 UA. Abscisas: tiempo en años, ordenadas: semieje mayor en UA.

Gomes et al. (2005) analizan el origen del Bombardeo Lunar Tardio en el marco del modelo de Nice, específicamente proponen que dicho pico de material entrante al Sistema Solar Interior es debido a la desestabilicación del disco de planetesimales durante el cruce de la resonancia 1:2 de Jupiter y Saturno. La Fig 11. muestra la evolución temporal de los semiejes de los planetas gigantes obtenida con nuestro código, con condiciones iniciales semejantes a las del trabajo mencionado. En nuestra simulación el cruce de la resonancia se produce a años.

Fig 11. Migración planetaria de los cuatro planetas gigantes en el marco del modelo de Nice para una simulación del origen del bombardeo lunar tardío. El disco se extiende inicialmente desde 15.5 UA hasta 34 UA con un perfil que decae inversamente con la distancia y una masa total de 35 masas terrestres. Ordenadas: distancia al perihelio q y al afelio Q, abscisas: tiempo en años.

Bibliografía/Referencias

  • Brunini, A. & Viturro, H., A tree code for planetesimal dynamics: comparison with a hybrid direct code, 2003, MNRAS, 346, 924-932.
  • Chambers, J. E., A hybrid simplectic integrator that permits close encounters between massive bodies, 1999, MNRAS, 304, 793.
  • Duncan, M. J., Levison, H. F., & Lee, M. H., A multiple time step symplectic algorithm for integrating close encounters, 1998, AJ, 116, 2067.
  • Pfalzner, S. & Gibbon, P., 1996, Many-Body Tree Methods in Physics, Cambridge University Press.
  • Santamaría, P., Un código de árbol para la dinámica de planetesimales, 2004, Tesis de Licenciatura.
  • Cionco, R. & Brunini, A., Orbital migrations in planetesimal discs: N-body simulations and the resonant dynamical friction, 2002, MNRAS, 334, 77-86.
  • Gomes, R., Morbidelli, A., Levison, H., Planetary migration in a planetesimal disk: why did Neptuno stop at 30 UA?, 2004, Icarus 170 492-507
  • Gomes, R., Levison, H., Tsiganis, K. & Morbidelli, A, Origin of the cataclysmic Late Heavy Bombardment period of the terrestrial planets, 2005, Nature, Vol 435|26.
  • Kokubo, E. & Ida, S., On runaway growth of planetesimals, 1996, Icarus 123, 180.
Última modificación: 2012/03/05 13:52