La validation des codes HYDRO-MUSCL et HYDRO-COOL 

Le code HYDRO-MUSCL que j’ai utilisé pendant ma thèse, pour éaliser les simulations numériques de processus astrophysiques, est un code développé par notre quipeé. La philosophie d’avoir un code « maison » est

de pouvoir connaître les différentes méthodes numériquesmplémentéesi dans celui-ci et de pouvoir facilement modifier des parties de ce code afin de modéliser un système par ticulier ou de résoudre un problème numérique apparent. Afin d’étudier les jets de matière et les chocs radiatifs présents dans divers systèmes, notre équipe a commencé par utiliser le code public CLAW [Leveque, 2002a]qu’il est possible de trouver sur Internet1 .

Mais l’architecture de ce code et le solveur de Riemann implémenté ne permettent pas de simuler des phéno-mènes à grand nombre de Mach. En effet, nous avons montré queM ≈ 4 était la limite de convergence de ce code. C’est la raison pour laquelle nous avons développé nu tout nouveau code avec un schéma numérique plus correct et un solveur de Riemann plus robuste. L’exactitude est une propriété qui suppose que le schéma ne produit pas une trop grande quantité de diffusion numérique lors de la discrétisation des termes d’advec-tion (voir Sec. 2.1.5). La robustesse est une autre propriété qui suppose que le schéma produit tout de même une faible quantité de diffusion numérique pour compenser al dispersion numérique qui résulte de la propaga-tion d’erreurs numériques. Un schéma robuste implique aussi que la solution approchée obtenue reste stable. Afin d’acquérir ces propriétés, nous avons introduit un schéma qui s’appuie sur la méthode MUSCL-Hancock (MUSCL signifiant « Monotonic Upstream-Centered Scheme for Conservation Laws ») et qui est une améliora-tion de la méthode de Godunov. Quant à l’algorithme qui résout le problème de Riemann, il est de la même famille que le solveur HLL déterminé par Harten, Lax, et van Leer. Le code HYDRO-MUSCL permet à présent de traiter les problèmes d’hydrodynamique qui nous intéresent, à savoir à grand nombre de Mach. Nous avons aussi implémenté une extension de ce code afin de pouvoir étudier des systèmes hydrodynamiques où les pertes radiatives par refroidissement sont modélisées simplement (code HYDRO-COOL). Nous développons un autre code où le couplage de l’hydrodynamique (HYDRO-MUSCL) et du transfert radiatif est pris en compte (code HADES).

Le code HYDRO-MUSCL est donc le premier aboutissement de cette démarche et sa conception date de 2007.

Les développements réalisés par H. C. Nguyenportent essentiellement sur le schéma et le solveur. Cette ver-sion a été parallélisée avec MPI (« Message Passing Interface ») par F. Roy en 2008 afin de gagner du temps de calcul et elle fonctionne de manière efficace pour 4 coeurs de calcul (4 coeurs de calcul formant souvent un processeur), c’est-à-dire pour son utilisation sur le cl uster de calcul de l’Observatoire de Paris (Siolino). Le code HYDRO-COOL est la première extension du code et elle a été terminée en 2009. Dans cette version, nous avons introduit, en terme source, une fonction de refroidissement afin de traiter les pertes radiatives dans l’ap-proximation d’un milieu optiquement mince (approximation forte et restrictive du couplage hydrodynamique et rayonnement). De plus H. C. Nguyen a rendu cette version plus accessible à l’utilisateur et l’a parallélisée afin de pouvoir utiliser un grand nombre de coeurs de calcul. L e troisième code, HADES, est la version la plus importante car elle représente le travail de thèse de H. C. Nguyen. Il a réalisé pour cette version un couplage de l’hydrodynamique et du rayonnement selon une méthode basée sur la résolution des équations de transfert radiatif par la méthode des moments (méthode M1) et par groupe de fréquence (méthode multigroupe). Ce qui permet d’avoir une description correcte de l’action de l’hydrodynamique sur le rayonnement et de la rétroaction du rayonnement sur l’hydrodynamique.

Etant donnée la chronologie de l’élaboration des différents codes, j’ai utilisé pendant ma thèse le code purement hydrodynamique HYDRO-MUSCL et le code incluant les pertes radiatives HYDRO-COOL. Mais c’est avec le premier code que j’ai obtenu le plus grand nombre de résultat que je présente dans les chapitres suivants. Avant de réaliser les simulations de systèmes particuliers, mon travail a consisté à tester ce nouveau code, afin de vérifier qu’il convergeait et qu’il présentait les différentes caractéristiques nécessaires à l’obtention de résultats physiques (robustesse, exactitude,. . .). C’est ce travail de test est présenté dans le chapitre sui-vant. Au cours de l’utilisation d’HYDRO-MUSCL, j’ai aussi pu évaluer les différents problèmes numériques qu’il pouvait engendrer et remonter à la source de ces erreur s. Toujours à l’aide de ce code, j’ai modélisé les coquilles de gaz choqué des RSN (Sec. 4.5), puis les mêmes coquilles perturbées par des hétérogénéités de densité comme celles présentes dans le MIS (Sec. 5.2.1). J’ai finalement étudié les coquilles perturbées par une fonction sinusoïdale (Sec. 5.2.2). C’est dans cette dernière modélisation que j’ai rencontré un important problème numérique, l’instabilité de grosseur (« carbuncle instability »), et c’est dans la partie correspondant à cette étude que je présenterai cet obstacle à l’étude de l’instabilité de Vishniac. En ce qui concerne les résultats obtenus avec le code HYDRO-COOL, j’ai réalisé des simulations du choc retour dans les colonnes d’accrétion des variables cataclysmiques magnétiques (dites polaires, chapitre III) et de propagation de jets d’étoile jeune (chapitre III ) , afin d’étudier l’effet du refroidissement sur la dynamique de ces systèmes. Mais pour l’instant, nous allons voir une description plus détaillée du code HYDRO-MUSCL.

La modélisation numérique d’un phénomène physique donne césac à des informations parfois inaccessibles à des modèles théoriques . Le code HYDRO-MUSCL, qui résout le système d’équations décrivant la 2 H. C. Nguyen est un doctorant de l’équipe, spécialisé en mathématiques appliquées.

3 Si nous nous plaçons dans le cas de l’hydrodynamique, le code  numérique résout les mêmes équations que celles à la base d’un modèle analytique décrivant le problème. Mais à la différence de ce modèle théorique où une (ou plusieurs) solutions analytiques dynamique des écoulements fluides, permet la réalisation desimulations de divers processus hydrodynamiques comme ceux apparaissant dans les RSN, les jets d’étoile jeune ou encore les polaires. Nous synthétisons dans cette introduction l’architecture du code HYDRO-MUSCL afin de bien comprendre le rôle des différentes mé-thodes numériques qui le constituent. Nous présentons ces méthodes plus en détails dans les sous-sections suivantes. Les équations qui sont résolues dans le code sontles équations d’Euler écrites sous forme conserva-tive (forme explicite de la conservation de la densité, de l’impulsion et de la densité d’énergie totale). Un terme source géométrique est aussi introduit de manière optionnelle afin de pouvoir résoudre les équations d’Euler axisymétriques car, dans ce cadre particulier, des termes upplémentaires apparaissent. Ces termes sources per-mettent ainsi de modéliser une symétrie cylindrique bien utile dans le cas des jets d’étoile jeune et des RSN. Pour pouvoir prendre en compte d’autres processus physiques spécifiques, nous les plaçons aussi en terme source si leur nature le permet. La grille, sur laquelle sont discrétisées les équations de l’hydrodynamique par la méthode des volumes finis, est une grille cartésienne à deux dimensions. Pour pouvoir résoudre le système d’équations différentielles partielles (EDP) avec d’éventuels termes sources, nous devons étudier le problème par morceaux. La scission est faite selon le découplage de Strang (« Strang splitting ») en espace, en temps et en type d’équations. En effet, nous discernons la résolution des équations d’Euler dans les deux directions de la grille et nous différencions la résolution du système homogène (EDP où le terme source est annulé) et la résolution des équations différentielles ordinaires (EDO, seule la différentielle totale par rapport au temps est considérée). Ainsi la solution globale du système est la réunion des solutions particulières dans chacune des parties qui s’effectue selon une opération particulière (le calcul de la nouvelle valeur des grandeurs conservées est effectué sur un demi pas de temps). Pour déterminer la solution de l’EDO, nous utilisons un schéma de Runge-Kutta explicite du second ordre qui préserve la stabilité du système sur un nombre important de pas de temps. En ce qui concerne la détermination de la solution des équations d’Euler homogène, nous utilisons le schéma de MUSCL-Hancock qui est une amélioration du schéma de Godunov car il permet d’obtenir une solution approchée au second ordre en espace et en temps de lasolution exacte. Dans cette partie, le flux présent sur les interfaces des cellules est déterminé. Afin de permetre ce calcul, le schéma est construit en trois temps : d’abord la solution approchée aux bords est déterminéevia la reconstruction des données faisant appel à un limiteur de pente afin de minimiser l’erreur, ensuite les sol utions approchées aux bords évoluent temporelle-ment, et finalement le solveur résout le problème de Riemann ce qui permet de déterminer la vitesse des ondes du système. Les solveurs que nous avons implémentés dans leodec sont les solveurs HLLC et HLLE et ils sont adaptés pour les cas physiques que nous voulons étudier. Le modèle de parallélisme utilisé dans cette version est MPI et cette méthode permet de découper le domaine de simulation en sous-domaines traités chacun par un coeur de calcul différent. Nous allons maintenant détailler le rôle des différents modules implémentés dans le code.

peuvent être déterminées, dans les simulations numériques, une multitude de solutions peut être obtenue. Il reste toutefois au numéricien d’être critique sur ce panel de solutions obtenues afin de discerner les solutions non physiques.

La grille cartésienne à deux dimensions permet aux différentes grandeurs d’évoluer sur le plan (x,y). En effet, chacune de ces grandeurs caractéristiques correspond à une onde qui se propage d’une maille (i, j) à la maille suivante (i ±1, j) (ou encore (i, j ±1)) et qui résulte de la différence des solutions présentes dans chaque maille. Ainsi la vitesse est divisée en deux composantes quenous noterons dans ce chapitre u sur l’axe Ox et v sur l’axe Oy . Si nous considérons l’axeOx, les noeuds de maille (. . . , xi−1, xi, xi+1, . . . ) sont positionnés au centre des cellules comme nous le voyons sur la figure 2.1. L e pas spatial Δx est donné, pour la cellulexi, xi−1xixi+1 xi−1xi+ par l’écart présent entre deux bords de cette même celluleΔx = xi+1/2 − xi−1/2. Nous précisons que la grille à deux dimensions est constituée de mailles carrées ou rectangulaires et que les pas spatiaux Δx et Δy y sont réguliers. Cette grille nous permet de discrétiser les équations d’Euler que nous présentons maintenant.

Les équations d’Euler dépendantes du temps forment un système de lois de conservation hyperboliques et non-linéaires. Le code que nous avons implémenté permeta lrésolution de ces équations sous forme conser-vative. Le système impliqué s’écrit de manière générale soucette forme particulière : δU +∇ H=S, (2.1) où les quantités qui interviennent dans cette équation sontle vecteur des grandeurs conservées U , le vecteur des flux H de ces mêmes grandeurs et le terme sourceS. Le terme ∇ H est, comme nous l’avons déjà vu, un terme d’advection. Les lois de conservation énoncent que, quand nous intégrons ce système d’équations, nous obtenons naturellement la conservation de la densité, de laquantité de mouvement selon les deux directions et de la densité d’énergie. Il reste donc à calculer les flux de ces quantités et le traitement de ce calcul est un point important dans le code.

La grille étant bi-dimensionnelle, les équations traitéesdans le code décrivent aussi un système à deux dimen-sions. Ainsi le vecteur de flux se divise en deux composantes F et G et l’Eq. 2.1 s’écrit [Toro, 1999] : Ut + F (U )x + G(U )y = S(U ) , (2.2) où les indices notent la dérivée partielle par rapport au temps t et à l’espace selon  Ox  et Oy . Les fonctions s’expriment :ρ U = ρu , ρv E √ ρu F = , ρuv ρu2 + p u(E + p) ρv G = , (2.3) ρv2 + p ρuv v(E + p) où E = p/(γ − 1) + (1/2)ρ u2 + v2 est toujours la densité d’énergie totale prenant en compte esl deux composantes de la vitesse, et F et G sont les deux fonctions de flux dans chacune des directions. D ans ce système, les fonctions décrivant les grandeurs physiques onts ρ = ρ(x, y, t), u = u(x, y, t), v = v(x, y, t), et p = p(x, y, t). Quant au vecteur des grandeurs conservées, sa fonction estnotéeU (x, y, t). Nous précisons que dans ce système, comme nous traitons encore un cas général où nous n’avons pas spécifié de propriétés géométriques, le terme sourceS est inexistant ou seulement un terme source physique. Par contre si le système est axisymétrique, c’est-à-dire s’il existe dans un système à trois dimensions une symétrie par rapport à un des axes ce qui permet de ramener le système à deux dimensions, alors un terme source géométrique apparaît en supplément d’un terme source physique le cas échéant.

Dans la configuration d’une symétrie cylindrique, nous déterminons les propriétés géométriques qu’il faut inclure à ce système. L’axe le long duquel s’effectue la symétrie est l’axe Oz qui se nomme maintenant la direction axiale. Nous introduisons aussi la coordonnée r qui mesure la distance à l’axe de symétrie Oz et l’axe qui résulte de cette variable est l’axe radialOr . Ainsi les équations d’Euler décrivant un espace à trois dimensions font normalement appel à trois coordonnées mais dans le cas axisymétrique, le système se réduit à deux dimensions. Ce système axisymétrique s’écrit [Toro, 999]1 : Ut + F (U )r + G(U )z = S(U ) , (2.4)

où U , F et G sont les mêmes vecteurs que précédemment Set = Sphys + Sgeo est le terme source global qui peut contenir une partie physique Sphys à laquelle vient s’ajouter la partie géométriqueSgeo qui s’écrit : ρuv 1ρu2 u(E + p)ρu Sgeo = −r . (2.5)

Si le terme source est non nul, la résolution du système d’équation n’est pas directe. En effet il faut diviser le problème en deux parties et résoudre chacune des parties éparéments. De la même manière, le système d’équations doit être résolu sur chacune des directions de’espacel. Ces opérations sont rendues possibles par le découplage de Strang que nous présentons maintenant.

Pour résoudre les équations d’Euler sous forme conservativ et avoir un schéma du second ordre en espace et en temps, nous utilisons le découplage de Strang du même ordre [Strang, 1968]. Mais avant de présenter cette méthode, nous discrétisons le vecteur des randeursg conservées selon l’espace et le temps : U (x, y, t) devient Ui,jn que nous notons U n lorsque seule une évolution temporelle du vecteur est effectuée. La nouvelle fonction Ui,jn est une solution approchée du vecteur de solutions physique. Dans cette formulation, le n N . Quant à l’espace, il est discrétisé par les temps est discrétisé par les éléments où l’entier n varie de 0 à points de maille xi,j où les entiers (i, j) varient respectivement de 1 − nf à mx + nf et de 1 − nf à my + nf où nf est le nombre de cellules fantômes et mx, my sont les maximums du nombre de cellules présentes dans la boîte de simulation pour les deux directions. Nous précisons que la valeur n = 0 détermine les conditions initiales et que les valeurs i = 1 −nf , 1 −nf + 1, . . . , 0 et i = mx + 1, . . . , mx + nf −1, mx + nf déterminent les conditions aux bords (même chose pour l’entierj). Nous formulerons dans ce qui suit les conditions aux bords qu’il est nécessaire de déterminer afin de pouvoir résoudre le système d’équations.

Le découplage de Strang s’effectue sur le temps, sur l’espace et sur le type d’équations d’Euler, comme ces différents aspects sont difficiles à traiter simultanément (mais ce n’est pas impossible). Cette méthode permet donc de simplifier grandement la résolution des équations. Pour le découplage spatial, le vecteur de solutions U n prend la valeur suivante au temps tn+1 : U n+1 = Y( 1 Δt)X(Δt) Y( 1 Δt)(U n) , (2.6)

où Δt est le pas de temps, et X et Y sont respectivement les opérateurs de solution de l’axeOx et de l’axe Oy .

Pour le découplage du type d’équation, le vecteur de solutions U n prend la valeur suivante au temps tn+1 : U n+1 = S( 1 Δt)Z(Δt)S( 1 Δt)(U n) , (2.7)

où S et Z sont respectivement les opérateurs de solution de l’EDO et de l’équation d’Euler homogène (EDP). Ainsi, pour chacun des découplages, la méthode de Strang permet de découper la résolution du système global en la résolution de deux systèmes séparément. Cette méthodefacilite donc la résolution du système d’équations. Nous détaillons la résolution des deux types d’équations dans le paragraphe suivant.

L’équation différentielle ordinaire s’écrit : d U ( , t) = S (U ( , t)) , (2.8)

où le point ( ) désigne n’importe quelle position de l’espace. Cette équation est résolue par un Runge-Kutta d’ordre 2 à condition que le terme source soit non raide (présence de gradients pas trop forts). Pour déterminer la méthode optimale, nous avons testé aussi la méthode de Heun, la méthode du point milieu, et des Runge-Kutta d’ordre 3 et 4 dans le cadre de simulations de jets d’étoile jeune (voir Nguyen et al. [2010]). Les résultats montrent que ces différentes méthodes ne changent pas la valeur des grandeurs caractéristiques (la variation est négligeable) étant donné le second ordre présent dans usto les schémas énoncés. Techniquement, nous calculons ce terme dans un module différent en faisant appelaux valeurs des grandeurs fluides calculées lors de la résolution des équations d’Euler homogènes présentéeci-après. Mais d’abord, nous rappelons les deux types de termes sources qui vont nous être utiles.

Le terme source géométrique

Nous avons vu que la symétrie pouvant être introduite dans laboîte de calcul est une symétrie cylindrique et qu’elle est permise par le terme source géométriqueSgeo. En effet ce terme permet d’introduire un axe de révolution sur un des bords de la boîte de simulation et donc d’être dans un système axisymétrique.

Le terme source physique : la fonction de refroidissement

Pour les différentes études effectuées dans cette thèse, usno n’avons utilisé qu’un seul type de terme source physique qui est lié à l’émission de rayonnement. Mais il existe d’autres termes sources de cette nature qui permettent d’étudier un système où un effet physique suplémentaire est pris en compte (par exemple, les équations d’Euler avec de la gravité). Le terme source physique Sphys que nous avons introduit est un vecteur contenant la fonction de refroidissement Λ présentée dans la Sec. 1.3. En effet, nous avons vu qu’une manière de traiter, dans un milieu optiquement mince, les processus radiatifs qui ont une rétroaction sur la dynamique du système est d’introduire une perte ou un gain d’énergie dans l’équation de conservation de l’énergie. Dans le cas du processus de refroidissement, la perte d’énergie se modélise par une fonction de refroidissement. Le terme source physique s’écrit dans ce cas : 0 0Sphys = − Λ(ρ, p) (2.9)