|
Menu |
Rapport2Julien DURANCEAU et Alexandra PINARD Sciences de la nature Groupe 2203 Simulation d'une collision de galaxies Rapport - dépôt final présenté à M.Martin AUBÉ Projet de fin d'étude Département de physique Cégep de Sherbrooke 20 mai 2009 TABLE DES MATIÈRESINTRODUCTIONLes galaxies sont des superbes structures observables à l'aide de téléscopes. Le téléscope orbital Hubble nous a révélé la beauté de ces amas d'étoiles gigantesques. La forme des galaxies est l'une de leurs caractéristiques appréciables. Cette forme serait créée lors de la collision de deux galaxies. Il est intéressant de s'interroger à propos des conditions qui mènent à telle ou telle autre forme de galaxie. Il est à noter que d'autres auteurs ont précédemment étudié l'idée que la topologie des galaxies serait le fruit d'interactions entre les amas d'étoiles. C'est ce qui sera étudié dans le cadre de ce projet. Le but de ce projet est de tester différentes conditions de collisions entre deux galaxies pour observer la forme qui en résulte. Différents facteurs pourront être contrôlés. Premièrement, la forme initiale des galaxies déterminée arbitrairement en début de projet. Deuxièmement, la forme, la vitesse, ainsi que le paramètre d'impact de la deuxième galaxie pourront également être contrôlés. Finalement, la durée de la modélisation pourra être modifiée pour observer plusieurs possibilités de formes résultantes en fonction du temps écoulé. Il est évident qu'il serait impossible de calculer l'action de la gravité sur tous ces corps à la main. C'est pourquoi l'aide de l'ordinateur est requise. La simulation de la gravité est effectuée par un programme qui utilise la méthode de Verlet pour évaluer l'effet temporel de l'interaction entre les différents corps. Le programme a été programmé en langage fortran. Pour faciliter le projet, les galaxies initiales auront une forme sphérique et la répartition des étoiles à l'intérieur de celles-ci sera copiée sur celle d'une galaxie réelle. Les résultats obtenus pourront être observés à l'aide d'un logiciel de graphiques, Gnuplot, qui servira à afficher soit la position des étoiles ou leur vitesse relative en un instant donné. Il est prévu que le déroulement de la collision puisse être observé en format vidéo en faisant dérouler, une après l'autre, plusieurs images des galaxies lors de la collision. CADRE THÉORIQUELa simulation numérique de collision de galaxies fait actuellement l'objet d'une multitude de projets de recherche auprès des équipes en astrophysique du monde entier. Notamment, une équipe française d'astrophysiciens a su récemment démontrer par numérisation qu'à la suite d'une collision, les galaxies spirales fusionnent en galaxies elliptiques, entourées d'amas globulaires. Ces chercheurs ont eu recours à des supercalculateurs, afin de calculer (pendant environ une cinquantaine d'heures) l’évolution gravitationnelle des étoiles, du gaz, et de la matière noire qui entoure les galaxies, ayant été modélisés à l'aide de 36 millions de particules numériques. C'est donc l'utilisation des supercalculateurs du CCRT, le Centre de Calcul Recherche et Technologie du CEA (Commissariat à l'énergie atomique ) qui a permis d'atteindre une performance de résolution jamais atteinte jusqu'alors, pour une telle collision, c'est-à-dire utiliser un très grand nombre de particules dans la collision. La technologie moderne permet donc à ce jour une énorme puissance de calcul pour effectuer les collisions de galaxies. Évidemment, nous ne possédons ni les connaissances, ni l'équipement nécessaire, ni le temps pour réaliser une collision d'une telle envergure. Il est à noter que la simulation dont il a été question ci-haut aurait nécessité plus d'une trentaine d'ordinateurs conventionnels et les calculs auraient duré environ deux ans. Ainsi, le cas étudié pour ce projet sera plutôt celui d'une collision un peu plus modeste, soit la collision de deux galaxies sphériques. L'idée première est de se créer une galaxie de départ, dotée des paramètres les plus réalistes possibles. Il serait donc intéressant de clarifier de prime abord quels sont les caractéristiques des galaxies qui composent notre univers. Les galaxies se présentent sous trois formes principales: sphérique, spirale et elliptique. Une galaxie moyenne s'étend sur environ 100 000 années-lumières de diamètre et est composée d'environ 100 milliards d'étoiles. Évidemment,la vitesse des étoiles est très variable; la vitesse maximale moyenne se situe autour de 100 000 m/s. À titre de comparaison, le Soleil est une étoile moyenne dans l'univers, dont la masse est de 1,98x1030 kg. Les étoiles les plus massives ont généralement une masse d'environ 50 à 80 fois celle du Soleil. Il existe aussi des étoiles plus petites, appelées les naines rouges, pouvant avoir jusqu'au vingtième de la masse solaire. Les récentes observations ont démontré la présence de trous noirs au centre des galaxies. Ces trous noirs géants ont en moyenne une masse de 2,7 millions de fois la masse solaire. Cette observation pourrait s'avérer utile lors des calculs. La concentration des étoiles dans une galaxie n'est pas uniforme. Il existe plusieurs modèles qui la simulent. Parmi ceux-ci, se trouvent le modèle de la sphère uniforme, celui de la sphère isothermale,les modèles de Plummer’s et Hernquist’s, ainsi que le modèle de Jaffe. Le modèle qui a été utilisé ici pour la distribution aléatoire des corps est celui de Jaffe, puisque les calculs qu'il nécessite s'avèrent plus simples à programmer, d'autant plus qu'il s'agit d'un modèle assez réaliste. Le modèle de Jaffe s'exprime donc ainsi:  où p(r) est la densité des corps en fonction du rayon
M est la masse totale des étoiles dans la galaxie
a est une constante qui permet de faire varier la concentration des étoiles dans la région du centre de la galaxie
r est le rayon de la sphère
Afin de simuler la gravité entre les corps de la galaxie modélisée, il faudra être en mesure de calculer la force agissant sur chacune des étoiles de la galaxie. Cette force est le résultat de la somme des forces gravitationnelles entre chacune des autres étoiles et l'étoile étudiée (l'étoile i). La formule pour calculer la force exercée sur une étoile i s'exprime donc ainsi: {#\vec F_i = \sum_{j} \frac{G M_j m_i}{{r_{ji}}^3} \vec r_{ji} = \sum_{j} \frac{G M_j m_i}{{r_{ji}}^3} ((x_i-x_j) \vec x + (y_i-y_j) \vec y + (z_i-z_j) \vec z)#} avec {# r_{ji} = \sqrt{ (x_i-x_j)^2 + (y_i-y_j)^2 + (z_i-z_j)^2} #} où Fi est la force exercée sur l'étoile i
G est la constante de gravitation
mj est la masse d'une étoile j
mi est la masse de l'étoile i
rji est le rayon entre une étoile j et l'étoile i
xi est la position en x de l'étoile i
xj est la position en x de l'étoile j
yi est la position en y de l'étoile i
yj est la position en y de l'étoile j
zi est la position en z de l'étoile i
zj est la position en z de l'étoile j
Ensuite, selon la deuxième loi de Newton  où ai est l'accélération de l'étoile i causée par la force Fi
Comme  Il existe plusieurs méthodes pour calculer numériquement l'évolution de la position et de la vitesse d'un corps i dans le temps. Une méthode simple serait celle d'Euler, par contre, elle a le désavantage de générer une grande erreur. Il y a aussi la méthode de Runge-Kutta. Celle-ci est beaucoup plus précise, mais elle a le défaut d'être assez complexe à programmer. Finalement, c'est la méthode de Verlet qui sera utilisée ici. Elle offre une précision raisonnable tout en étant plutôt simple à programmer. Cette méthode peut s'exprimer comme suit:  où xi+1 est la position finale de l'étoile i après un temps Δt
xi est la position initiale de l'étoile i
x'i est la vitesse initiale de l'étoile i
Δt est le temps du pas de calcul
xi'' est l'accélération de l'étoile i pendant le temps Δt
 où x'i+1 est la vitesse finale de l'étoile après un temps Δt
xi+1'' est l'accélération finale de l'étoile après un temps Δt
On peut comprendre que la méthode de Verlet utilise les conditions de position, de vitesse et d'accélération à un temps t0 pour prévoir leurs nouvelles valeurs à t=t0+dt. Il est aussi possible de calculer l'énergie potentielle d'une galaxie. Pour se faire, il suffit de faire la sommation des énergies potentielles entre chacune des étoiles du système et de diviser cette valeur par 2.  où Epotgal est l'énergie potentielle totale de la galaxie
Epot est l'énergie potentielle entre deux étoiles de la galaxie
 m1 est la masse d'une certaine étoile
m2 est la masse d'une autre étoile
r1,2 est la distance entre les étoiles 1 et 2
Enfin, la détermination des conditions initiales de la collision nécessite certains éléments théoriques. Notamment, il sera utile de savoir que la distance moyenne réelle entre les galaxies dans l'univers équivaut à dix fois leur diamètre. Aussi, la loi de dispersion de Hubble sera utile pour déterminer les vitesses relatives entre les galaxies lors de notre simulation. Cette loi permet de prédire, de façon approximative, la vitesse entre deux galaxie éloignées d'une certaine distance. Le graphique suivant illustre la loi de Hubble: Figure 1: Graphique de la loi de Hubble servant à la détermination de la vitesse de collisions tiré de Edward L. Wright, 2007  La figure 1 permet de déduire la constante de Hubble qui est d'environ 64km/s/Mpc. 1Mpc représentant 3,086 X1022 mètres. CADRE MÉTHODOLOGIQUEPremièrement, un programme sera écrit pour créer une galaxie initiale aléatoire, dont la concentration en étoiles suivra le modèle de Jaffe. Pour se faire, la sphère représentant notre galaxie sera divisée en tranches, un peu à la manière d'un oignon. Il sera possible de calculer le nombre d'étoiles à l'intérieur de chacune des tranches de rayon delta r comme suit: et  où dm est la masse totale des étoiles contenues dans une couche de la galaxie
dv est le volume d'une couche de la galaxie.
Le volume de la couche peut être estimé avec une grande précision en multipliant l'aire de la sphère par la largeur de la couche  En manipulant les équations:  Pour obtenir le nombre d'étoiles contenues dans une couche de largeur delta r, il faut diviser la masse totale dans une couche par la masse d'une étoile:  où meest la masse d'une étoile
En divisant la masse totale des étoiles dans la galaxie par la masse d'une étoile (en supposant que toutes les étoiles possèdent la même masse), on obtient le nombre d'étoiles dans la galaxie.  où dN est le nombre d'étoiles contenues à l'intérieur d'une tranche de la galaxie
Ntot est le nombre d'étoiles contenues dans la galaxie
Pour raccourcir le temps de calcul, un nombre d'étoiles inférieur à celui que l'on retrouve dans les galaxies réelles sera utilisé. Toutefois, dans un souci de réalisme, la masse totale de la galaxie sera celle d'une galaxie réelle. Ainsi, en divisant la masse totale de la galaxie par le nombre d'étoiles utilisé, la masse de chacune des étoiles sera déterminée. Lorsque le nombre d'étoiles contenues à l'intérieur d'une tranche de la galaxie est connu, il faut ensuite que ces étoiles soient réparties de façon aléatoire à l'intérieur de cette tranche. Figure 2: graphique facilitant l'interprétation des angles  Un premier angle est tiré au hasard, c'est l'angle alpha. Cet angle se situe dans le plan xy. Ensuite, pour que la distribution soit bien aléatoire, il faut s'assurer que le nombre d'étoiles tirées par disque d'angle alpha soit proportionnel à la circonférence de ce disque. Si cette précaution n'est pas prise, il y aura un grand regroupement d'étoiles au niveau des pôles. Il est connu que la circonférence est proportionnelle au rayon d'un disque. Comme il est possible d'observer sur le schéma, le rayon du disque est donné par rcos(alpha). La méthode consiste donc à tirer plusieurs angles alpha par tranches de galaxie. Par la suite, un certain nombre d'angles phi est tiré, proportionnellement à la valeur du cosinus de chaque angle alpha. Finalement, dans toute cette banque d'étoiles possibles pour une tranche de la galaxie, un nombre d'étoiles déterminé précédemment à l'aide de la méthode de Jaffe est tiré au hasard. On peut comprendre que pour que les étoiles demeurent liées à l'intérieur de la galaxie par la gravité, il ne faut pas que leur énergie cinétique dépasse la valeur absolue de leur énergie potentielle gravitationnelle. Pour faciliter les calculs, lors du calcul de l'énergie potentielle d'une étoile, toutes les autres étoiles (ayant un rayon inférieur par rapport au centre de la galaxie que celui de l'étoile étudiée) sont assimilées à une géante masse placée au centre de la galaxie dont la masse est la somme des masses des étoiles ayant un rayon inférieur. C'est à partir de ce principe que la vitesse maximale des étoiles pourra être déterminée:  où Ek est l'énergie cinétique d'une étoile
Ep est l'énergie potentielle de l'étoile
 où v est la vitesse d'une étoile
m1 est la masse de l'étoile
m2 est la somme des masses des étoiles ayant un rayon inférieur
r1,2 est le rayon entre l'étoile et le centre de la galaxie.
En isolant la vitesse:  Cette vitesse correspond à la vitesse maximale que peuvent avoir les étoiles sans théoriquement être éjectées. Ainsi, notre stratégie pour déterminer la vitesse des étoiles consiste à donner à celles-ci une vitesse comprise entre 0 et la dite vitesse maximale. Par contre, même en respectant cette stratégie de détermination de la vitesse, il se peut que certaines étoiles soient tout de même éjectées de la galaxie. En effet, lorsque deux étoiles passent trop près l'une de l'autre, il peut se produire un effet de fronde, amplifié par les erreurs de calculs possibles, ce qui peut provoquer l'éjection des étoiles concernées. Pour remédier à la situation, une distance critique a été déterminée entre deux étoiles, en-dessous de laquelle les étoiles seront agglomérées. Pour déterminer cette distance critique, le critère de roche sera utilisé: {#R_{roche}= R_{roche1}+R_{roche2}#} où Rroche est la distance critique entre deux étoiles
Rroche1 est le rayon de roche de la première étoile
Rroche2 est le rayon de roche de la deuxième étoile
{#R_{roche1}={m_1}^{1/3*0,14}#} où m1 est la masse de l'étoile 1
{#R_{roche2}={m_2}^{1/3*0,14}#} où m2 est la masse de l'étoile 2
Stratégie expérimentale:Lorsqu'une galaxie numérique aura été créée correctement, il faudra premièrement faire des tests pour tenter de déterminer l'influence du pas de calcul sur les résultats pour que l'erreur ne soit pas trop importante, tout en gardant un temps de calcul raisonnable. Ensuite, il faudra tester les paramètres de la galaxie pour s'assurer qu'elle demeure bien en équilibre quand elle est soumise, seule, à l'action de la gravité. Si ce n'est pas le cas, quelques paramètres pourront être modifiés, par exemple, le paramètre a, influant sur la concentration des étoiles. Aussi, le nombre d'étoiles a peut-être une influence sur l'état d'équilibre, il peut être modifié. Ensuite, une modification sur le diamètre de la galaxie pourrait influer sur sa capacité à conserver un certain état d'équilibre. Par ailleurs, l'ajout d'un trou noir au centre de la galaxie comme dans la réalité pourrait aider à la stabilisation des étoiles. Pour les premiers calculs, une galaxie d'environ 500 étoiles sera utilisée. Si les vitesses démontrent un équilibre satisfaisant, un objet de grande masse (s'apparentant à un trou noir) sera lancé dans la galaxie, selon plusieurs angles et vitesses différents. Ensuite, deux galaxies différentes pourront être lancées l'une vers l'autre. Pour ce qui est des conditions initiales de la collisions, quatre paramètres seront à déterminer, soit la distance verticale initiale entre les deux galaxies, le paramètre d'impact, la vitesse des galaxies ainsi que la durée des collisions. Premièrement, pour être plus réaliste, il faudrait idéalement que la distance verticale initiale entre les galaxies soit de dix diamètres, mais par contrainte de temps de calcul, cette distance a été fixée à 6x1020 m, ce qui correspond approximativement à 6 rayons de galaxie. Deuxièmement, le paramètre d'impact, correspondant à la distance perpendiculaire à la vitesse entre les deux galaxies, constitue la principale variable de l'expérimentation. Il a été fixé à approximativement deux rayons de galaxie (2x1020m),un rayon de galaxie (1x1020m) et finalement un troisième paramètre d'impact nul. La vitesse a quant à elle été fixée à l'aide de la loi de Hubble comme suit: {#v = h r#} où v est la vitesse entre les deux galaxies
h est une constante dont la valeur est 2,07 x 10-21 km/(s*m)
r est est la distance entre les deux galaxies
{#v = 2,07 \times 10^{-21} * 6 \times 10^{20}#} v = 1242 m/s Encore une fois, la contrainte du temps de calcul nous force à prendre des vitesses plus rapides qu'en réalité. Avec une vitesse de 1244 m/s, le temps réel de la collision aurait été extrêmement long, donc le temps nécessaire au calcul d'une telle collision aurait aussi été très long. C'est pourquoi il a été décidé que les deux vitesses testées allaient être 200 000 et 300 000 m/s. Finalement, la durée de collision est déterminée en fonction de la vitesse et de la distance verticale initiale des galaxies, de manière à ce qu'un certain temps s'écoule à la suite de la collision pour que l'on puisse observer les changements engendrés. De plus, une animation de la collision entre les deux galaxies sphériques sera produite. Il s'agit d'utiliser une grande quantité de photos de la collision et ce, à différents moments. Ensuite, les photos, se succédant rapidement une à la suite de l'autre, donneront une idée visuelle de la collision. Finalement, tout dépendant de la qualité des résultats obtenus et du temps qu'il nous reste, le projet pourra être perfectionné en modifiant certains paramètres. Par exemple:
MATÉRIEL, INSTRUMENTATION ET EXPÉRIMENTATIONLe matériel utilisé est:
Expérimentation:
Pour créer des animations en format gif des collisions:
Pour analyser les valeurs énergétiques de la galaxie:
OBSERVATIONS, INTERPRÉTATION ET RÉSULTATSFigure 3: Simulation de l'orbite terrestre. Cette figure démontre le bon fonctionnement de l'algorithme dynamique puisque la terre suit une orbite cyclique sur une année.  Figure 4: graphique du rayon d'une étoile en fonction du pas de calcul d'une simulation de 10 000 ans.  Figure 5: Résultat pour une galaxie d'environ 600 étoiles en gravité durant 100 millions d'années démontrant la stabilité de cette galaxie   Figure 5: résultat de l'essai 1  Figure 6: résultat de l'essai 2: (vue de devant la collision)  Figure 7: résultat de l'essai 3  Figure 8: résultat de l'essai 3: (vue rapprochée)  Figure 9: résultat de l'essai 4  Figure 10: résultat de l'essai 5  Figure 11: résultat de l'essai 6:   DISCUSSIONComparaison et/ou validation des résultatsValidation de la méthodologie:Premièrement, afin de valider l'algorithme numérique de simulation de la gravité, l'orbite de la Terre autour du Soleil a été simulé. Pour ce faire, les paramètres spécifiques du système Terre-Soleil ont été insérés dans le programme "galaxie" et le temps de simulation a été fixé à un an. Les résultats démontrent (voir animation 1) que la Terre, un an après le début de sa révolution autour du Soleil, est revenue sensiblement à sa position initiale. Ainsi, comme le programme est parvenu à simuler la révolution de la Terre autour du Soleil d'une manière réaliste, l'algorithme numérique de la simulation de la gravité a été validé. Deuxièmement, il a fallu déterminer un pas de calcul raisonnable pour la suite des calculs. Le graphique 3 représentant le rayon d'une étoile à la fin d'une courte simulation (10 000 ans) démontre l'erreur qui est commise au fur et à mesure que la valeur du pas de calcul augmente. Il est possible de constater que l'erreur augmente rapidement lorsque le pas de calcul augmente. Il a été décidé que le pas de calcul allait être de 100 ans. cette décision peut sembler bizarre, puisque l'erreur est relativement élevée à ce pas de calcul. Par contre, il faut tenir compte de temps de calcul. En utilisant un pas de calcul inférieur à 100 ans, le temps de calcul ne convient pas avec le temps qui nous est alloué pour l'expérimentation. De plus, avant de procéder à la simulation d'une collision, il fallait valider la stabilité d'une seule galaxie. La simulation de l'effet de gravité a donc été effectuée sur une galaxie comportant 608 étoiles, et ce pour une durée de l'ordre de celle d'une collision, soit 100 million d'années. En observant l'animation résultante (animation 2), il est possible de remarquer que la grande majorité des étoiles demeurent à l'intérieur du rayon initial de la galaxie. Toutefois, une dizaine d'étoiles sont carrément éjectées de la galaxie, phénomène qui pourrait fortement être dû à un pas de calcul trop élevé pour la densité d'étoiles à ces endroits. Néanmoins, comme seules une dizaine d'étoiles ont été éjectées sur les 608 étoiles, le phénomène a été considéré négligeable, et la galaxie a été qualifiée comme suffisamment stable. Analyse qualitative:Essai 1: Il est premièrement possible d'observer un changement considérable au niveau de la trajectoire de la deuxième galaxie. La vitesse de cette deuxième galaxie, initialement orientée uniquement vers le haut, gagne une composante horizontale à mesure que la collision se déroule. L'effet d'attraction entre les deux galaxies est donc flagrant, d'autant plus que la première galaxie acquiert une certaine vitesse, alors qu'elle était initialement au repos. Un certain changement dans la forme des galaxies peut également être observé. Ce changement se traduit par une dispersion des étoiles, particulièrement dans l'axe horizontal ou perpendiculaire à la collision. Les galaxies tendent donc apparemment à adopter une forme elliptique plus ou moins aplatie, plutôt que sphérique. Essai 2: Dans ce deuxième essai, en comparaison avec le premier essai, les mêmes changements au niveau de la trajectoire des galaxies sont notables. Toutefois, la dispersion s'effectue d'une manière particulière, en raison d'un paramètre d'impact inférieur. Il est en fait possible de remarquer que la dispersion s'accentue davantage du côté opposé à la collision, et ce pour les deux galaxies. Cela pourrait s'expliquer par le fait que les étoiles près du site de la collision subissent davantage les interactions avec les étoiles de la galaxie voisine. Autrement dit, plus les étoiles sont éloignées du site de collision, moins elles suivent le mouvement de la galaxie engendré par la collision. Comme le mouvement des étoiles n'est pas uniforme, la dispersion qui en résulte ne peut l'être; les étoiles plus éloignées sont plus dispersées. Essai 3: Le changement observable dans la forme des galaxies pour ce troisième essai concerne encore une fois la dispersion des étoiles. Avec cette fois-ci un paramètre d'impact nul, les galaxies se rencontrent de plein fouet. Ainsi, un maximum d'étoiles entrent en interaction avec celles de la galaxie voisine, ce qui entraîne une dispersion encore plus forte que dans les deux premiers essais. Aussi, le troisième essai présente une caractéristique particulière que nous n'avons pas repérée dans les cinq autres tests. En fait, il est possible d'apercevoir qu'il y a un certain échange d'étoiles d'une galaxie à l'autre lors de la collision; il semblerait que l'effet de gravité ait été assez puissant pour créer les dites échanges. Essai 4,5 et 6: Comme le tableau 1 l'indique, les quatrième, cinquième et sixième essais comprennent exactement les mêmes paramètres que les essais 1, 2 et 3 respectivement, à l'exception près d'une plus grande vitesse dans les trois derniers essais. Ainsi, comme le paramètre d'impact et la distance entre les galaxies de l'essai 4 et de l'essai 1 sont les mêmes, les changements dans la trajectoire et la dispersion des galaxies sont très similaires. La nuance ici est que les effets sont atténués avec une plus grande vitesse. Ceci s'explique par le fait que les étoiles disposent de moins de temps pour interagir avec les autres. Analyse quantitative:Peu d'analyse peut être faite quantitativement à partir des données recueillies suite aux collisions de galaxies. Cependant, des résultats étonnants ont pu être observés au niveau de l'énergie potentielle des galaxies avant et après leur collision. Comme le montre le tableau 2, dans chacun des 6 essais, l'énergie potentielle finale des galaxies est supérieure à leur énergie potentielle initiale. Ce gain d'énergie potentielle est observable visuellement par la dispersion des étoiles à l'extérieur du rayon initial de la galaxie suite à la collision. Un autre fait notable est que non seulement les énergies finales sont supérieures aux énergies initiales, mais de plus, la valeur de l'énergie potentielle de chacune des galaxies en interaction tend vers une valeur semblable. Initialement, les galaxies ont une énergie potentielle respective assez différente, en raison du tir aléatoire, mais suite à la collision, ces valeurs sont beaucoup plus rapprochées. Il est possible d'en conclure qu'il y a eu échange d'énergie entre les galaxies durant la collision. ConclusionCe projet avait comme objectif principal de réussir, à l'aide d'une suite de programmes informatiques, à modéliser une collision entre 2 galaxies. Le principal défi était de faire en sorte que l'action de la gravité soit calculée étoile par étoile et non à l'aide d'une quelconque approximation. Cette modélisation allait permettre de tirer des conclusions au sujet des changements géométriques et dynamiques des galaxies initiales suite à leur collision. Elle allait aussi permettre de faire varier différents paramètres de la collision et d'observer le résultat de ces variations sur les galaxies finales. Somme toutes, les programmes mis au point sont assez efficaces et permettent une belle visualisation des collisions. Au niveau quantitatif, les résultats sont moins probants, mais quelques conclusions ont tout de même pu être tirées. Il a été démontré que moins la vitesse de collision était grande, plus les interactions entre les galaxies étaient importantes. Aussi, plus le paramètre d'impact est petit, plus la dispersion des étoiles est importante durant la collision. Au niveau énergétique, l'énergie potentielle des galaxie est plus grande suite à leur collision et il a été conclu qu'un certain échange d'énergie était possible entre les galaxies durant les collision. Malheureusement, aucune formation géométrique significative n'a pu être observée dans une galaxie suite à sa collision avec une autre galaxie. Finalement, une des pistes de solutions envisageables pour l'obtention de forme géométrique remarquable suite à une collision de galaxies serait de faire en sorte que cette galaxie ait un mouvement initial rotatif très important. Cette condition pourrait mener à l'obtention de bras spiraux suite à une collision. Aussi, il serait intéressant de tester le résultat de collision entre 2 galaxies de forme particulière, par exemple des galaxies de forme elliptique. MÉDIAGRAPHIE1.Auteur inconnu. cea (page consultée le 26 février2009), http://www.cea.fr/le_cea/actualites/simulation_formation_galaxie_elliptique-8947 2.Auteur inconnu. Annual Reviews (page consultée le 4 février 2009), http://nedwww.ipac.caltech.edu/level5/March04/Binney/Binney3.html 3.COMBES, Françoise. interstices (page consultée le 25 mars 2009), http://interstices.info/jcms/c_19147/levolution-des-galaxies 4.HAWKING, Stephen. 5 Galaxy models (page consultée le 4 février 2009), http://www.kof.zcu.cz/st/dis/schwarzmeier/galaxy_models.html AnnexeProgramme aleatoire: c programme pour creer une galaxie spherique aleatoire c declaration des variables double precision alpha(1000000),phi(1000000),alph,x,y,deltar, + a,f,g double precision alphaf(1000000),phif(1000000),posx(1000000), + posy(1000000) double precision posz(1000000),vx(1000000),vy(1000000),x0,z0, + vz(1000000) double precision xtmp,ytmp,ztmp,v0 c martin < double precision vmax,v,facteur,masse,rayonact,rayon,massei c > martin integer k,i,j,n,jmax,nbre,nn,nmax,m,pos,ii,nbmasse,nbgal,ng,ne c a est un parametre fixe arbitrairement qui influe sur la concentration des etoiles au centre de la galaxie. c nn est le nombre total d'etoiles dans la galaxie c deltar est la largeur d'une couche de la galaxie c k permet de faire varier le nombre de couches de la galaxie c posx, posy et posz sont les positions en coordonnees cartesiennes des etoiles c m est un compteur qui compte le nombre d'etoiles total c vmax est la vitesse maximale d'une etoile masse=4.e+37 a=19.46e+19 deltar=4.73e+17 facteur=1000000000000000000. m=0 ne=0 open(unit=1,file='3d',status='unknown') open(unit=2,file='vector',status='unknown') print*,'Nombre de galaxies?' read*,nbgal do ng=1,nbgal print*,'Position initiale de la galaxie ',ng,' en x et z' read*,x0,z0 print*,'Vitesse initiale de la galaxie ',ng,' en z' read*,v0 nn=1000 do k=1,1000 n=0 nbre=dnint((deltar*dble(nn)*a)/((a+(dble(k)*deltar))**2.)) nmax=dnint(sqrt(dble(nbre)/0.5914)) do i=1,100 c avec une boucle de 100 et 100 cos alpha on obtiens 6227 valeurs c On tire ici la banque des angles pour la position call RANDOM_NUMBER(x) alph=(x-0.5)*3.1415926 jmax=dnint(100.*cos(alph)) do j=1,jmax n=n+1 alpha(n)=alph call RANDOM_NUMBER(y) phi(n)=(y*6.2831853) c print*,alpha(n),phi(n) enddo enddo c print*,'n=',n do i=1,nbre call RANDOM_NUMBER(f) pos=dnint(f*dble(n)) m=m+1 alphaf(m)=alpha(pos) phif(m)=phi(pos) c print*,k,m, alphaf(m),phif(m) posx(m)=dble(k)*deltar*cos(alphaf(m))*cos(phif(m)) posy(m)=dble(k)*deltar*cos(alphaf(m))*sin(phif(m)) posz(m)=dble(k)*deltar*sin(alphaf(m)) c write(1,*) posx(m),posy(m),posz(m) c martin < enddo c do i=1,m c > martin c martin < pos=dnint(g*dble(n)) alphaf(m)=alpha(pos) phif(m)=phi(pos) c > martin c alphaf(m)=alpha(v) c phif(m)=phi(v) c vx(m)=dble(k)*deltar*cos(alphaf(m))*cos(phif(m)) c vy(m)=dble(k)*deltar*cos(alphaf(m))*sin(phif(m)) c vz(m)=dble(k)*deltar*sin(alphaf(m)) c write(1,*) posx(m),posy(m),posz(m),vx(m),vy(m),vz(m) c print*,k,m, alphaf(m),phif(m),vx(m),vy(m),vz(m) c martin < print*,m,'/',nn c > martin enddo enddo nn=m m=0 n=0 do i=1,100 c avec une boucle de 100 et 100 cos alpha on obtiens 6227 valeurs call RANDOM_NUMBER(x) alph=(x-0.5)*3.1415926 jmax=dnint(100.*cos(alph)) do j=1,jmax n=n+1 alpha(n)=alph call RANDOM_NUMBER(y) phi(n)=(y*6.2831853) c print*,alpha(n),phi(n) enddo enddo do i=1,nn call RANDOM_NUMBER(f) pos=dnint(f*dble(n)) m=m+1 alphaf(m)=alpha(pos) phif(m)=phi(pos) call RANDOM_NUMBER(g) rayonact=sqrt(posx(m)**2.+posy(m)**2.+posz(m)**2.) nbmasse=0 do ii=1,nn rayon=sqrt(posx(ii)**2.+posy(ii)**2.+posz(ii)**2.) if (rayon.lt.rayonact) then nbmasse=nbmasse+1 endif enddo massei=dble(nbmasse)*masse vmax=sqrt((2.*6.67428e-11*massei)/rayonact) v=g*vmax vx(m)=v*cos(alphaf(m))*cos(phif(m)) vy(m)=v*cos(alphaf(m))*sin(phif(m)) vz(m)=v*sin(alphaf(m))+v0 xtmp=posx(m)+x0 ytmp=posy(m) ztmp=posz(m)+z0 write(2,*) xtmp,ytmp,ztmp,vx(m)*facteur, + vy(m)*facteur,vz(m)*facteur write(1,*) xtmp,ytmp,ztmp,vx(m),vy(m),vz(m),masse ne=ne+1 print*,'Etoile no.',ne print*,i,'/',nn enddo m=0 enddo close(unit=1) close(unit=2) stop end Programme galaxie: c c programme pour simuler le probleme du mouvement a N corps sous c l action de la force de gravite c programme par Martin Aube c utilisation de la methode de verlet c c ===== c Declaration des variables c integer size parameter (size=2048) double precision xa(size),xn(size),ya(size),yn(size),za(size), + zin(size) double precision dt,dti,G,rroche1,rroche2,rroche double precision m(size),pi,ray double precision vxa(size),vya(size),vxn(size),vyn(size), + vzia(size),vzin(size) double precision axa(size),aya(size),axn(size), + ayn(size),aza(size),azin(size) integer step,nloop,i,j,mout,n,ne,nm,nn,ii character*30 nomvect c c ===== c Identification des variables c c m(n)=masse de la masse n c nm=nombre de masses c xa(n) ya(n) za(n) les position de la particile n a t_i c xn(n) et yn(n) zn(n) sont les position de la particile n a t_n c vxa(n) et vya(n) vzia(n) sont les vitesses de la particile n a t_i c vxn(n) et vyn(n) vzn(n) sont les vitesses de la particile n a t_n c axa(n) et aya(n) aza(n) sont les accelerations de la particile n a t_i c axn(n) et ayn(n) azn(n) sont les accelerations de la particile n a t_n c ray=rayon entre deux masses c dt=pas de calcul c dti=pas de calcul initial c Tm,To,Ts periodes orbitalesza c nloop=nombre de boucles de calcul c mout=saut pour l impression de sorties c c c Commentaires: c Le programme attend en entree un fichier nomme ci.txt (conditions initiales) c et retourne les conditions finales cf.txt dans le meme format c le format est le nombre de masses sur la premiere ligne suivie sur les lignes c suivantes (une ligne par masse) des x y z vx vy vz m c c Le programme fait 2E9 boucles de calculs temporels. La variable dti (a ajuster c avant de compiler) est le pas de calcul. Pour connaitre le temps apres les c 2E9 boucles il faut faire tf = 2E9 x dti - 0.99dti c Dans sa version originale dti=1jour ce qui implique un tf>5E6 annees c c ===== c Valeurs initiales c c dti=86400*7 => une semaine c dti=86400.*365.25*50. nloop=4000000 mout=1000000 n=0 pi=3.14159 G=6.6725985E-11 dt=dti nn=0 c c ===== c Lecture des positions et vitesses initiales c open(unit=1,file='ci.txt',status='unknown') print*,'Reading initial conditions: file ci.txt' read(1,*) nm if (nm.gt.size) then print*,'Stars number exceed max value! Abort.' stop endif do i=1,nm read(1,*) xa(i),ya(i),za(i),vxa(i),vya(i),vzia(i),m(i) enddo close(unit=1) c c c c c ===== c c c ouverture du fichier de sortie open(unit=2,file="cf.txt",status="unknown") do n=1,nloop do i=1,nm axn(i)=0. ayn(i)=0. azin(i)=0. do j=1,nm if (i.ne.j) then ray=((xa(j)-xa(i))**2.+(ya(j)-ya(i))**2.+ + (za(j)-za(i))**2.)**0.5 c calcul du rayon de roche c (16^.333*rayon soleil/Masse du soleil ^0.333)*Masse ^0.333=0.14 rroche1=m(i)**0.33333333333*0.14 rroche2=m(j)**0.33333333333*0.14 rroche=rroche1+rroche2 if (ray.lt.rroche) then print*,'limite de roche atteinte' xa(i)=(m(i)*xa(i)+m(j)*xa(j))/(m(i)+m(j)) ya(i)=(m(i)*ya(i)+m(j)*ya(j))/(m(i)+m(j)) za(i)=(m(i)*za(i)+m(j)*za(j))/(m(i)+m(j)) vxa(i)=(m(i)*vxa(i)+m(j)*vxa(j))/(m(i)+m(j)) vya(i)=(m(i)*vya(i)+m(j)*vya(j))/(m(i)+m(j)) vzia(i)=(m(i)*vzia(i)+m(j)*vzia(j))/(m(i)+m(j)) m(i)=m(i)+m(j) m(j)=0. xa(j)=100. ya(j)=0. za(j)=0. vxa(j)=0. vya(j)=0. vzia(j)=0. endif axn(i)=axn(i)+G*m(j)*(xa(j)-xa(i))/ray**3. ayn(i)=ayn(i)+G*m(j)*(ya(j)-ya(i))/ray**3. azin(i)=azin(i)+G*m(j)*(za(j)-za(i))/ray**3. endif enddo if (n.eq.1) then c c ===== c Ajuste la premiere valeur prededente a la valeur actuelle c dt=dti/100. axa(i)=axn(i) aya(i)=ayn(i) aza(i)=azin(i) else dt=dti endif vxn(i)=vxa(i)+0.5*(axa(i)+axn(i))*dt vyn(i)=vya(i)+0.5*(aya(i)+ayn(i))*dt vzin(i)=vzia(i)+0.5*(aza(i)+azin(i))*dt xn(i)=xa(i)+vxa(i)*dt+0.5*axa(i)*dt**2. yn(i)=ya(i)+vya(i)*dt+0.5*aya(i)*dt**2. zin(i)=za(i)+vzia(i)*dt+0.5*aza(i)*dt**2. nn=nn+1 c c ===== c Impression des resultats c if (nn.eq.mout) then nn=0 print*," temps(jours)=",(dble(n-1)*dti+dti/100.)/3600./24. + ,'/',dble(nloop)*dti/3600./24., + nint(100.*(dble(n-1)+1/100.)/dble(nloop)),'%' endif nomvect='toto' if (nloop/50.eq.n) nomvect='vecteur01.txt' if (2*nloop/50.eq.n) nomvect='vecteur02.txt' if (3*nloop/50.eq.n) nomvect='vecteur03.txt' if (4*nloop/50.eq.n) nomvect='vecteur04.txt' if (5*nloop/50.eq.n) nomvect='vecteur05.txt' if (6*nloop/50.eq.n) nomvect='vecteur06.txt' if (7*nloop/50.eq.n) nomvect='vecteur07.txt' if (8*nloop/50.eq.n) nomvect='vecteur08.txt' if (9*nloop/50.eq.n) nomvect='vecteur09.txt' if (10*nloop/50.eq.n) nomvect='vecteur10.txt' if (11*nloop/50.eq.n) nomvect='vecteur11.txt' if (12*nloop/50.eq.n) nomvect='vecteur12.txt' if (13*nloop/50.eq.n) nomvect='vecteur13.txt' if (14*nloop/50.eq.n) nomvect='vecteur14.txt' if (15*nloop/50.eq.n) nomvect='vecteur15.txt' if (16*nloop/50.eq.n) nomvect='vecteur16.txt' if (17*nloop/50.eq.n) nomvect='vecteur17.txt' if (18*nloop/50.eq.n) nomvect='vecteur18.txt' if (19*nloop/50.eq.n) nomvect='vecteur19.txt' if (20*nloop/50.eq.n) nomvect='vecteur20.txt' if (21*nloop/50.eq.n) nomvect='vecteur21.txt' if (22*nloop/50.eq.n) nomvect='vecteur22.txt' if (23*nloop/50.eq.n) nomvect='vecteur23.txt' if (24*nloop/50.eq.n) nomvect='vecteur24.txt' if (25*nloop/50.eq.n) nomvect='vecteur25.txt' if (26*nloop/50.eq.n) nomvect='vecteur26.txt' if (27*nloop/50.eq.n) nomvect='vecteur27.txt' if (28*nloop/50.eq.n) nomvect='vecteur28.txt' if (29*nloop/50.eq.n) nomvect='vecteur29.txt' if (30*nloop/50.eq.n) nomvect='vecteur30.txt' if (31*nloop/50.eq.n) nomvect='vecteur31.txt' if (32*nloop/50.eq.n) nomvect='vecteur32.txt' if (33*nloop/50.eq.n) nomvect='vecteur33.txt' if (34*nloop/50.eq.n) nomvect='vecteur34.txt' if (35*nloop/50.eq.n) nomvect='vecteur35.txt' if (36*nloop/50.eq.n) nomvect='vecteur36.txt' if (37*nloop/50.eq.n) nomvect='vecteur37.txt' if (38*nloop/50.eq.n) nomvect='vecteur38.txt' if (39*nloop/50.eq.n) nomvect='vecteur39.txt' if (40*nloop/50.eq.n) nomvect='vecteur40.txt' if (41*nloop/50.eq.n) nomvect='vecteur41.txt' if (42*nloop/50.eq.n) nomvect='vecteur42.txt' if (43*nloop/50.eq.n) nomvect='vecteur43.txt' if (44*nloop/50.eq.n) nomvect='vecteur44.txt' if (45*nloop/50.eq.n) nomvect='vecteur45.txt' if (46*nloop/50.eq.n) nomvect='vecteur46.txt' if (47*nloop/50.eq.n) nomvect='vecteur47.txt' if (48*nloop/50.eq.n) nomvect='vecteur48.txt' if (49*nloop/50.eq.n) nomvect='vecteur49.txt' if (50*nloop/50.eq.n) nomvect='vecteur50.txt' c c ===== c Ecrire les donnees en format vector pour gnuplot c x,y,z,deltax,deltay,deltaz if (nomvect.ne.'toto') then open(unit=3,file=nomvect,status="unknown") write(3,*) nm do ii=1,nm write(3,*) xa(ii),ya(ii),za(ii),vxa(ii)*1.e15,v + ya(ii)*1.e15,vzia(ii)*1.e15 enddo close(unit=3) endif c c ===== c store la valeur precedentes c axa(i)=axn(i) aya(i)=ayn(i) aza(i)=azin(i) vxa(i)=vxn(i) vya(i)=vyn(i) vzia(i)=vzin(i) xa(i)=xn(i) ya(i)=yn(i) za(i)=zin(i) enddo enddo c c ===== c Ecrire les donnees de sortie c write(2,*) nm do i=1,nm write(2,*) xa(i),ya(i),za(i),vxa(i),vya(i),vzia(i),m(i) enddo close(unit=2) stop end Programme sepgal: #!/bin/bash # usage: sepgal n_stars #file="vecteur10.txt" let ngal=$1/2 let ndemi=ngal+1 #head -$ndemi $file | tail -$ngal > "gal1"$file #exit 0 list=`ls -1 vecteur*.txt` for i in $list do echo $i echo $ngal > "gal1"$i head -$ndemi $i | tail -$ngal >> "gal1"$i echo $ngal > "gal2"$i tail -$ngal $i >> "gal2"$i done Programme mkanim2gal: #!/bin/bash # usage: mkanim axis_limits list=`ls -1 vecteur*.txt` for i in $list do echo "set term gif ; set multiplot ; set pointsize 0.4; set data style points; set view 0,0; set xrange [-"$1":"$1"]; set yrange [-"$1":"$1"]; set zrange [-"$1":"$1"]; splot 'gal1"$i"' ; replot 'gal2"$i"'" >toto.1 gnuplot -persist toto.1 > $i.gif done convert -delay 10 -loop 100 vecteur*gif animation.gif animate animation.gif Programme vcmsansm: c programme pour calculer le vitesse et la position du centre de masse c d'une galaxie a partir d'un fichier ci.txt ou cf.txt c double precision vx1(1000),vy1(1000),vz1(1000),x1(1000), +y1(1000),z1(1000),vxcm1,vycm1,vzcm1,xcm1,ycm1,zcm1,m,mcm,v1, +bidon,mcin1,vmax double precision ecm1,epot1,r1,G,ecin1,l1,lx1,ly1,lz1,normx1, +normy1,normz1 double precision norme1,erot1,miner1,rcm1,vcm1,nomoy1 double precision vx2(2000),vy2(2000),vz2(2000),x2(2000), +y2(2000),z2(2000),vxcm2,vycm2,vzcm2,xcm2,ycm2,zcm2,v2, +mcin2 double precision ecm2,epot2,r2,ecin2,l2,lx2,ly2,lz2,normx2, +normy2,normz2 double precision norme2,erot2,miner2,rcm2,vcm2,nomoy2,ecin double precision r,epot,etot,ecint,ek integer i,netoile,j,nmax character*2 nom character*40 nom1,nom2 G=6.6725985e-11 print*,'Numero du fichier (vecteur##.txt)' read*,nom nom1='gal1vecteur'//nom//'.txt' nom2='gal2vecteur'//nom//'.txt' xcm1=0. ycm1=0. zcm1=0. vxcm1=0. vycm1=0. vzcm1=0. xcm2=0. ycm2=0. zcm2=0. vxcm2=0. vycm2=0. vzcm2=0. c print*,'Entrez la masse d une etoile:' c read*,m open(unit=2,file='ci.txt',status='old') read(2,*) bidon read(2,*) bidon,bidon,bidon,bidon,bidon,bidon,m close(unit=2) print*,'Masse d une etoile=',m mcm=0 open(unit=1,file=nom1,status='old') read(1,*) netoile ecin1=0. miner1=0. lx1=0. ly1=0. lz1=0. nomoy1=0. nmax=0 do i=1,netoile read(1,*) x1(i),y1(i),z1(i),vx1(i),vy1(i),vz1(i) xcm1=xcm1+x1(i)*m ycm1=ycm1+y1(i)*m zcm1=zcm1+z1(i)*m vx1(i)=vx1(i)/1.e15 vy1(i)=vy1(i)/1.e15 vz1(i)=vz1(i)/1.e15 vxcm1=vxcm1+vx1(i)*m vycm1=vycm1+vy1(i)*m vzcm1=vzcm1+vz1(i)*m ecin1=ecin1+0.5*m*(vx1(i)**2.+vy1(i)**2.+vz1(i)**2.) c calcul du moment cinetique normx1=(y1(i)-ycm1)*(vz1(i)-vzcm1)-(z1(i)-zcm1)*(vy1(i)-vycm1) normy1=(z1(i)-zcm1)*(vx1(i)-vxcm1)-(x1(i)-xcm1)*(vz1(i)-vzcm1) normz1=(x1(i)-xcm1)*(vy1(i)-vycm1)-(y1(i)-ycm1)*(vx1(i)-vxcm1) norme1=sqrt(normx1**2.+normy1**2.+normz1**2.) nomoy1=nomoy1+norme1 c normx=normx/norme c normy=normy/norme c normz=normz/norme rcm1=sqrt((x1(i)-xcm1)**2.+(y1(i)-ycm1)**2.+(z1(i)-zcm1) +**2.) vcm1=sqrt(vx1(i)**2.+vy1(i)**2.+vz1(i)**2.) lx1=lx1+m*normx1 ly1=ly1+m*normy1 lz1=lz1+m*normz1 miner1=miner1+m*rcm1**2. mcm=mcm+m enddo xcm1=xcm1/mcm ycm1=ycm1/mcm zcm1=zcm1/mcm vxcm1=vxcm1/mcm vycm1=vycm1/mcm vzcm1=vzcm1/mcm v1=sqrt(vxcm1**2.+vycm1**2.+vzcm1**2.) ecm1=0.5*dble(netoile)*m*v1**2. l1=sqrt(lx1**2.+ly1**2.+lz1**2.) erot1=l1**2./(2.*miner1) nomoy1=m*nomoy1/dble(netoile) print*,'Position CM gal1=',xcm1,ycm1,zcm1 print*,'Vitesse CM gal1=',vxcm1,vycm1,vzcm1 print*,'grandeur vitesse CM gal1=',v1 print*,'energie cinetique du CM gal1=',ecm1 close(unit=1) ek=0. do i=1,netoile vmax=sqrt((vx1(i)-vxcm1)**2.+(vy1(i)-vycm1)**2.+ +(vz1(i)-vzcm1)**2.) if (vmax.gt.2.e5) then ek=ek+0.5*m*vmax**2. print*,'On depasse c!',ek, vmax nmax=nmax+1 print*,nmax endif enddo ecin1=ecin1-ek c calcul de l'energie potentielle c epot1=0. do i=1,netoile do j=1,netoile if (i.ne.j) then r1=sqrt((x1(i)-x1(j))**2.+(y1(i)-y1(j))**2.+(z1(i)- +z1(j))**2.) epot1=epot1-G*m*m/r1 endif enddo enddo epot1=epot1/2. print*,'Energie potentielle propre galaxie 1',epot1 mcm=0 open(unit=1,file=nom2,status='old') read(1,*) netoile ecin2=0. miner2=0. lx2=0. ly2=0. lz2=0. nomoy2=0. do i=1,netoile read(1,*) x2(i),y2(i),z2(i),vx2(i),vy2(i),vz2(i) xcm2=xcm2+x2(i)*m ycm2=ycm2+y2(i)*m zcm2=zcm2+z2(i)*m vx2(i)=vx2(i)/1.e15 vy2(i)=vy2(i)/1.e15 vz2(i)=vz2(i)/1.e15 vxcm2=vxcm2+vx2(i)*m vycm2=vycm2+vy2(i)*m vzcm2=vzcm2+vz2(i)*m ecin2=ecin2+0.5*m*(vx2(i)**2.+vy2(i)**2.+vz2(i)**2.) c calcul du moment cinetique normx2=(y2(i)-ycm2)*(vz2(i)-vzcm2)-(z2(i)-zcm2)*(vy2(i)-vycm2) normy2=(z2(i)-zcm2)*(vx2(i)-vxcm2)-(x2(i)-xcm2)*(vz2(i)-vzcm2) normz2=(x2(i)-xcm2)*(vy2(i)-vycm2)-(y2(i)-ycm2)*(vx2(i)-vxcm2) norme2=sqrt(normx2**2.+normy2**2.+normz2**2.) nomoy2=nomoy2+norme2 c normx=normx/norme c normy=normy/norme c normz=normz/norme rcm2=sqrt((x2(i)-xcm2)**2.+(y2(i)-ycm2)**2.+(z2(i)-zcm2)**2.) vcm2=sqrt(vx2(i)**2.+vy2(i)**2.+vz2(i)**2.) lx2=lx2+m*normx2 ly2=ly2+m*normy2 lz2=lz2+m*normz2 miner2=miner2+m*rcm2**2. mcm=mcm+m enddo xcm2=xcm2/mcm ycm2=ycm2/mcm zcm2=zcm2/mcm vxcm2=vxcm2/mcm vycm2=vycm2/mcm vzcm2=vzcm2/mcm v2=sqrt(vxcm2**2.+vycm2**2.+vzcm2**2.) ecm2=0.5*dble(netoile)*m*v2**2. l2=sqrt(lx2**2.+ly2**2.+lz2**2.) erot2=l2**2./(2.*miner2) nomoy2=m*nomoy2/dble(netoile) print*,'Position CM gal2=',xcm2,ycm2,zcm2 print*,'Vitesse CM gal2=',vxcm2,vycm2,vzcm2 print*,'grandeur vitesse CM gal2=',v2 print*,'energie cinetique du CM gal2=',ecm2 close(unit=1) ek=0. do i=1,netoile vmax=sqrt((vx2(i)-vxcm2)**2.+(vy2(i)-vycm2)**2.+ + (vz2(i)-vzcm2)**2.) if (vmax.gt.5.e5) then ek=ek+0.5*m*vmax**2. print*,'On depasse c!',ek, vmax nmax=nmax+1 print*,nmax endif enddo ecin2=ecin2-ek c calcul de l'energie potentielle c epot2=0. do i=1,netoile do j=1,netoile if (i.ne.j) then r2=sqrt((x2(i)-x2(j))**2.+(y2(i)-y2(j))**2.+(z2(i)- +z2(j))**2.) epot2=epot2-G*m*m/r2 endif enddo enddo epot2=epot2/2 print*,'Energie potentielle propre galaxie 2',epot2 c calcul de l energie potentielle totale do i=1,netoile do j=1,netoile if (i.ne.j) then r=sqrt((x1(i)-x1(j))**2.+(y1(i)-y1(j))**2.+(z1(i)- +z1(j))**2.) epot=epot-G*m*m/r endif enddo do j=1,netoile r=sqrt((x1(i)-x2(j))**2.+(y1(i)-y2(j))**2.+(z1(i)- +z2(j))**2.) epot=epot-G*m*m/r enddo enddo do i=1,netoile do j=1,netoile if (i.ne.j) then r=sqrt((x2(i)-x2(j))**2.+(y2(i)-y2(j))**2.+(z2(i)- +z2(j))**2.) epot=epot-G*m*m/r endif enddo do j=1,netoile r=sqrt((x2(i)-x1(j))**2.+(y2(i)-y1(j))**2.+(z2(i)- +z1(j))**2.) epot=epot-G*m*m/r enddo enddo epot=epot/2. ecin=ecin1+ecin2 etot=epot+ecin1+ecin2 ecint=ecin1-ecm1+ecin2-ecm2 print*,'Bilan total' print*,'Energie potentielle totale (les deux galaxies) =', +epot print*,'Energie cinetique totale (les deux galaxies) =',ecin +,ecin1,ecin2 print*,'Energie totale (les deux galaxies) =',etot print*,'Energie cinetique interne (rot et disp.)=', +ecint c print*,'Energie de rotation interne',erot print*,vxcm2,vycm2,vzcm2 stop end RésuméDans ce projet, il est question de la simulation numérique de collisions entre deux galaxies sphériques, dans le but d'analyser les caractéristiques géométriques et énergétiques des galaxies suite à leur collision. Pour y parvenir, une suite de cinq programmes a été créée en langage fortran. Ces programmes ont été mis au point dans le but de créer les galaxies initiales, de simuler sur ces dernières l'effet de la gravité, de mettre en collisions deux galaxies et de visualiser les collisions sous forme d'animations. Lors de ces collisions, les deux galaxies utilisées sont identiques et de formes sphériques. Les étoiles qui les constituent sont tirées de façon aléatoire, suivant le modèle de distribution aléatoire de Jaffe. De plus, dans un souci de précision et de fidélité à la réalité, l'effet de la gravité dans les galaxies est simulé en faisant la sommation de toutes les interactions gravitationnelles étoile par étoile. À partir de cette force, il est possible de déduire l'accélération d'une étoile au début de la simulation. À partir de cette accélération, la méthode de Verlet est utilisée. Cette méthode permet de prédire les valeurs de position et de vitesse d'une particule après un temps delta t à partir de ses conditions initiales de position, de vitesse et d'accélération. Les collisions étudiées impliquent des galaxies de 329 étoiles et sont simulées sur une durée variant de 130 à 200 millions d'années. Les vitesses relatives entre chacune des galaxies sont soit de 200 000 ou de 300 000 m/s. Les collisions ont été simulées avec succès, les programmes permettant de visualiser les collisions fonctionnent bien et certaines conclusions ont pu être tirées au niveau énergétique. Les résultats montrent que les collisions font augmenter la dispersion des étoiles dans les galaxies et que plus la vitesse est lente, plus cette dispersion est importante. Au niveau énergétique, les galaxies finales ont plus d'énergie potentielle qu'avant la collision et il a été démontré qu'il y a échange d'énergie entre les galaxies durant la collision. AbstractThis project consists in simulating a digital collision between two spherical galaxies. The goal of this simulation is to analyze geometric and energetic characteristics of the two galaxies after their collision. To be able to do it, five programs were created in fortran language. These programs were finalized to create initial galaxies, simulate gravity effect on them, realize a collision between them and visualize the collision in a video format. During these collisions, both galaxies are spherical and identical. The stars in these galaxies are placed randomly following Jaffe’s random model. Moreover, in a worry of being precise and accurate, the gravity effect in the galaxies is simulated by adding up every gravitational interaction between each star. Using this force, it is possible to calculate the acceleration of a star at the beginning of the simulation. Using this acceleration, Verlet’s method is used. This method makes it possible to predict the position and the speed of a particle after a period delta t of time from its initial conditions of position, speed and acceleration. The created collisions imply 329 stars per galaxy and are simulated on a duration which varies between 130 and 200 billion of years. The speed between each galaxy is either 200 000 or 300 000 m/s. The collisions were simulated successfully, the programs used to visualize the collision work well and we were able to draw some facts about energy. The results show that the dispersion of the stars in the galaxies is higher after the collision then before. Moreover, when the speed is slower, this dispersion is more important. Concerning energy, the final galaxies have more potential energy than before the collision and it was demonstrated that there is energy exchange between the galaxies during the collision. |