Dveloppement dun code MHD 3D

1
Développement dun code MHD 3D pour le couplage
photosphère couronne G. Aulanier, P. Démoulin,
R. Grappin
1. De la physique au schéma numérique et aux
conditions aux limites
2. Etude paramétrique dun champ bipolaire
Stabilité des tubes de flux torsadés
Vecteur-B photosphérique
Modèle pour les sigmoïdes XUV
2
Objectifs
i Structuration 3D de la couronne Génération
de courants électriques alignés (taches-d,
protubérances, sigmoïdes) i Phénomènes
éruptifs Reconnexion magnétique topologies
complexes (flares confinées, CMEs) i Chauffage
par dissipation Ohmique Formation de courants
intermittents (mouvements BF braiding et HF
ondes)
Lobjectif est davoir un code MHD applicable à
des observations via des extrapolations
3
1. Sur la photosphère
4
Mouvement photosphérique Emergence depuis la
zone convective
i Ce quil se passe  en vrai  Les tubes de
flux émergent de façon ou torsadés depuis la
zone de convection Des mouvements de
cisaillement/torsion photosphériques cohérents
sont rarement observés i Contraintes
numériques dapplicabilité aux
observations Numérique Changts de régimes
brutaux du haut de la ZC à la couronne
Applicabilité Historique conduisant à un
magnétogramme observé ?
5
La photosphère parroi physique ?
i Contraste photosphère / couronne r 9
ordres de grandeur 1017 109 cm-3 b 3
ordres de grandeur 1 0.001 L
2 ordres de grandeur 150 km 10 Mm (sauf
pour ondes HF) gt Les ondes émises par la
couronne sont REFLECHIES par la photosphère ?
u(z0) 0 ? dbz/dt(z0) 0 gt Les mvts
photosphériques créent des VITESSES dans la
couronne ? u(z0) gt 0 ? u(zgt0) gt 0
gt Un point dancrage mouvant Ne déplace PAS
lautre point dancrage ? u(x1 z0) gt 0 ?
u(x2 z0) 0 i Approche forçage
photosphérique de la couronne On applique
uxuy(z0) gt0 là où on veut générer des
courants électriques uxuy(z0) 0 là où
la photosphère ne force pas On applique des
conditions infiniment conductrices (h 0) pour
ne pas faire jouer un rôle physique autre
quadvectif à la frontière
6
2. Présentation du code
7
Les équations de la MHD coronale
i Equation de continuité dt r (u .s) r
r (s.u) terme advectif terme compressible
i Equation dEuler dt u (u .s) u
(mr)1 (sx b) x b n ru terme advectif
force de Lorentz pseudo-force visqueuse i
Equation dinduction dt b (u .s) b b
(s.u) (b .s) u h rb terme idéal
terme résistif i gravité, pression, termes
Hall négligés (zéro-b collisionnel)


8
Discrétisation et dérivation spatiale
i Maillage fixe non uniforme Toutes les
quantités et les dérivées calculées aux bords des
mailles dxyz(n1) dxyz(n) x r r lt
1.03 200 km lt dxyz lt 2800 km maillage
nx ny nz 201 (8.12 x 106 points) i
Calcul des dérivées première et seconde pour les
s, r dzfi f ( fi-2 fi-1 fi fi1 fi2
dzi-2 dzi-1 dzi dzi1 ) Fonction f
calculée par combinaison de 4 devts de Taylor à
lordre 3 On voit tout de suite quil
faudra prescrire quelque chose aux frontières
(./..) i Méthode des différences finies au 3ème
ordre
fi-2 fi-1 fi
fi1 fi2
dzi-2
dzi-1 dzi
dzi1
z
9
Discrétisation et intégration temporelle
i Pas de temps variable contraint par la
propagation dinfo sur léchelle dune
maille Toutes les quantités et les dérivées
calculées aux bords des cellules dti C x min(
dxyz/(u cA) d2xyz/n d2xyz/h ) C
(facteur de Courant-Friedrich-Levy) 0.5
dti ( u0 cA1000 km/s) 0.1 s i Calcul des
champs rub à un temps tdt f Pi1 g (
dtfi-2 dtfi-1 dtfi
dti-2 dti-1 dti ) Adams-Bashforth
3-step f Ci1 h ( dtfi-1 dtfi
dtf Pi1 dti-1 dti )
Adams-Moulton 2-step Fonctions g et h
calculées par combinaison de 4 devts de Taylor à
lordre 3 Fonctions dtf Pi1 calculées par une
intégration des équations MHD i Méthode
prédicteur-correcteur
fi-2 fi-1 fi
fi1 dtfi-2 dtfi-1
dtfi dtfi1

dti-2
dti-1 dti
t
10
Définition du filtre diffusif pour u
i Pourquoi ne pas utiliser le  vrai  terme
visqueux Numériquement on doit dissiper à
léchelle de la maille les vrais n sont trop
faibles ! Adapter n sur chaque maille détendue
diminue le pas de temps en d2 avec un r i
Utilisation dun pseudo-laplacien r ?
filtre On calcule les dérivées seconde sur
i1,nx pas sur idx(1),dx(nx) rA
d2iA d2jA d2kA Ai1 2Ai Ai1 Aj1
2Aj Aj1 Ak1 2Ak Ak1 Puis on choisit
chaque terme diffusif en fonction dune vitesse
caractéristique un tel que n un /
minboîte(dxdydz) Ce qui revient à léchelle
de chaque maille à Re u / un



11
3. Intégration de conditions aux
limites réalistes et numériquement stables
12
Quelques mots sur les conditions aux limites

• i Conditions périodiques
• Très avantageuses, car bien définies et
  utilisation possible de la FFT
• - Problème pour traiter la photosphère et le
  haut de la boîte
• i Conditions de Dirichlet-Neumann
• Neumann Pair Dirichlet Impair (implique
  que bnorm 0 ou btrans 0)
• - Bien adaptées pour des parrois
  réfléchissantes glissantes
• - Difficultés car en général, les 6 faces dune
  boite ont bxbybz 0
• la photosphère naura pas tjs uxuy
  0
• Conditions ouvertes
• Méthode des caractéristiques
• - possible mais bon pour ondes seulement

13
Conditions aux limites avec cellules fantômes 1.
Photosphère  paroi rigide conductrice 
z y x
i Conditions sur r Neumann symétrique i
Conditions sur b ux uy Pseudo-Dirichlet anti
symétrique autour de la frontière Stable.
Erreurs sur les dn minimisées i Conditions sur
uz Dirichlet antisymétrique (réfléchissant) i
Conditions sur les termes MHD diffusifs n h
0 advection idéale Limite lamplitude de ces
termes dans le domaine
r z

bx,by,bz,ux,uy z

uz z

cellules fantôme
domaine
14
Conditions aux limites avec cellules fantômes 2.
Couronne  ouverte  aux 5 faces
z y x
i  Zero gradient extrapolation  Méthode
 classique  en (M)HD si près des frontières
l échelle des gradients A (dnA) 1 gt
quelques mailles Transmet le fluide et les
ondes (codes VAC, ZEUS, NRL) Mais réflexion
dondes parasites HF de  faible  amplitude i
une pseudo-frontière Couronne, près des
frontières A (dnA) 1 taille dune
maille Forces de Lorentz jxb  parasites 
conduisant à u cA Ajout dune
 pseudo-frontière  pour u - calculer jxb
seulement là ou les cellules fantômes
ninterviennent pas - 2 cellules dans le
domaine pour les dérivées à 5 points
r,bx,by,bz z

ux,uy,uz z

cellules fantôme
domaine
15
4. Calcul du comportement des tubes de flux
torsadés avec photosphère infiniment conductrice
inertielle Problème déjà étudié mais
avec dautres types de schémas numériques et
conditions aux limites

• Amari et al. (1996) Implicite (filtrage
  des ondes), pas déq pour r, grand domaine fermé
• ? gonflement plus rapide quune exponentielle
• Törok Kliem (2003) Lax-Wendroff, lissage
  des quantités, grand domaine fermé
• ? perte déquilibre à Ntour gt 1.4 1.5

16
Génération dun champ bipolaire théorique
i Magnétogramme synthétique S 2
gaussiennes bz (xyz0t0) bo exp
(xxc) y 2 / a2 bo exp (xxc) y2 /
a2 bo 650 G xc 8 Mm a 15 Mm i
Calcul du champ potentiel sx b 0 avec BLFFF
(Démoulin et al., 1996) Calcul avec conditions
périodiques en x200,200, y200,200 Domaine
dintégration choisi x99,99, y99,99,
z0,200
b (t0) z 3.8
17
Champ de vitesse de torsion appliqué

• Champ de vitesse photosphérique imposé
  dtbz(z0) uz 0
• on souhaite imposer une champ de vitesse
  transverse ut incompressible st . ut 0
• dbz/dt st (bz ut) ut . st (bz) 0
  ltgt ut st y(bz) x k
• i Choix du potentiel y(bz) pour la torsion
• y(xy) exp ( bz2(xyz0) max
  bz2(xyz0) ) / db2 db 10.
  bzmax (xyz0)
• 0.9 bzmax (xyz0)
• 0.5 bzmax (xyz0)
• ux (xyz0t) uo a(t) dy y(xy)
• uy (xyz0t) uo a(t) dx y(xy)
• i Profil daccélération graduelle avant une
  phase ut cst
• a(t) ½ tanh( 2 t ta / tw ) ½
  ta 300 s tw 100 s

gt
18
Conditions initiales champ potentiel
i Représentation de b à z gt 0
partie centrale du domaine
domaine complet
bz contours /- 50 100 150 200 250 300
350 400 450 G max(bz) 480 G max(uy)
20 km/s
19
Conditions initiales vortex étendus (Run
A) db 10. bzmax (xyz0)
partie centrale du domaine
bz (y0,z0)
uy (y0,z0)
bz contours /- 50 100 150 200 250 300
350 400 450 G max(bz) 480 G max(uy)
20 km/s
20
Conditions initiales vortex moyens (Run
B) db 0.9 bzmax (xyz0)
partie centrale du domaine
bz (y0,z0)
uy (y0,z0)
bz contours /- 50 100 150 200 250 300
350 400 450 G max(bz) 480 G max(uy)
20 km/s
21
Conditions initiales vortex concentrés (Run
C) db 0.5 bzmax (xyz0)
partie centrale du domaine
bz (y0,z0)
uy (y0,z0)
bz contours /- 50 100 150 200 250 300
350 400 450 G max(bz) 480 G max(uy)
20 km/s
22
Conditions initiales densité termes diffusifs
i Paramètres de calcul Densité et vitesse
dAlfvèn initiales cA(x,y,z,t0) 1000
km/s r (x,y,z,t0) m-1 b2(x,y,z,t0) /
cA(x,y,z,t0) Termes diffusifs n un d
3 1010 m2/s un 150 km/s
min(dx,dy,dz) 200 km h uh d 3 109 m2/s
uh 15 km/s
Re (l 30 Mm) 100 Rm (l 30 Mm) 1,000
évalués pour u 100 km/s 0.1 cA(x,y,z,t0)
Lu (l 30 Mm) 10,000
23
Evolution du tube de flux db 10. bzmax
(xyz0)
ligne de champ axiale lignes ancrées à bz 475
G lignes ancrées à bz 400 G lignes ancrées à bz
300 G lignes ancrées à bz 250 G lignes
ancrées à bz 150 G lignes ancrées à bz 50 G
Champ potentiel
50 100 150
t 0630 s
24
Evolution du tube de flux db 10. bzmax
(xyz0)
ligne de champ axiale lignes ancrées à bz 475
G lignes ancrées à bz 400 G lignes ancrées à bz
300 G lignes ancrées à bz 250 G lignes
ancrées à bz 150 G lignes ancrées à bz 50 G
50 100 150
t 0630 s
25
Evolution du tube de flux db 10. bzmax
(xyz0)
ligne de champ axiale lignes ancrées à bz 475
G lignes ancrées à bz 400 G lignes ancrées à bz
300 G lignes ancrées à bz 250 G lignes
ancrées à bz 150 G lignes ancrées à bz 50 G
50 100 150
t 0840 s
26
Evolution du tube de flux db 10. bzmax
(xyz0)
ligne de champ axiale lignes ancrées à bz 475
G lignes ancrées à bz 400 G lignes ancrées à bz
300 G lignes ancrées à bz 250 G lignes
ancrées à bz 150 G lignes ancrées à bz 50 G
Torsion du tube rose 1 tour
50 100 150
t 1020 s
27
Evolution du tube de flux db 10. bzmax
(xyz0)
ligne de champ axiale lignes ancrées à bz 475
G lignes ancrées à bz 400 G lignes ancrées à bz
300 G lignes ancrées à bz 250 G lignes
ancrées à bz 150 G lignes ancrées à bz 50 G
50 100 150
t 1110 s
28
Evolution du tube de flux db 10. bzmax
(xyz0)
ligne de champ axiale lignes ancrées à bz 475
G lignes ancrées à bz 400 G lignes ancrées à bz
300 G lignes ancrées à bz 250 G lignes
ancrées à bz 150 G lignes ancrées à bz 50 G
50 100 150
t 1170 s
29
5. Analyse de la dynamique du système (lentement)
forcé
30
Energétique
i Ce quil se passe Saturation de EB
- expulsion du domaine  ouvert  - champ
partiellement ouvert Chute de dt -
vidage de la photosphère - vitesses
croissantes
31
Dynamique
i Comportement en deux phases Expansion
lente( uz uxy lt 0.1 cA) , puis
quasi-Alfvénique ( uz gt 0.5 cA gtgt
uxy) Sensibilité du système à la taille des
vortex i Instabilité (e.g. kink) ? Perte
déqulibre (e.g. ligne de courant dans b
potentiel) ?
32
6. Comportement du système en relaxation MHD
33
Méthode(s) de relaxation
i Décélération graduelle du forçage Rampe de
diminution de uxy(z0) g 0 sur quelques temps
dAlfvén Utilisé par Amari et al. (1996) et
Törok Kliem (2003)
i Redémarrage en coupant les vitesses
Recommencer un calcul en prenant pour
conditions b (initial) b (au
temps t du calcul à torsion continue) u
(initial) 0 uxy(x, y, z0 , t)
0 Laisser les forces j x b initialement
présentes agir diminuer ? g relaxation vers
un équilibre force-free non linéaire ?
ou samplifier ? g perte déquilibre ?
34
Energétique dynamique de quelques relaxations
Ntour 1.05
Ntour 0.92
Ntour 1.48
Ntour 1.41
35
Energétique dynamique de quelques relaxations
i Oscillations autour dun équilibre
Expansion ralentissement redescente rapide
Jamais dexpansion monotone, seulement des
OSCILLATIONS (viscosité pas trop
forte) Altitude d équilibre  z (dtuz
dzuz 0 )
36
7. Construction des courbes déquilibre
37
Positions des  équilibres 
i Courbe générique déquilibre z z0 exp A
N2tour A (gros vortex) 2.16
A (moyens vortex) 1.12 A (petits
vortex) 0.50 Cest sévère ! Cette courbe
est prédite analytiquement en cisaillement sphéri
que axisymétrique !
38
Existe til une bifurcation ?
Törok Kliem (2003)
A 2.16
A 1.12
A 0.75
A 0.50
i Bifurcation ou pas ? Amari et al. (1996)
gonflement de en rapide et relaxations
de en longues ? Törok Kliem (2003)
perte déquilibre à Ntour gt 1.4 1.5 au
moment où la courbe déquil nest pas très raide
? Nous courbe
déquilibre Z générique z z0 exp A N2tour
jamais de perte déquilibre même là où
Z est raide jxb z (pendant les
oscillations) faibles 0.1 bidzbk
39
Arguments contre une instabilité de kink
i Instabilité de kink en géométrie cylindrique
(rqz) Pour un tube aussi large que haut
instabilité ? N gt 1.25 tour (Hood, 1992) Pour
des rapports daspect variables instabilité ?
bq / bz gt 1 (Baty, 2001) gt tubes long
étroits instabilité ? N par unité de longueur
i En géométrie cartésienne (xyz) Pour un
nombre de tour N donné gt 1, si la taille des
vortex k ? Pourcentage du flux torsadé
k ? Expansion du tube selon z
k ? longueur du tube torsadé
k ? Torsion par unité de longueur
m Pour une taille de vortex pas trop petite, si
N k ? Altitude déquilibre
k en exp N2 !!! ? Expansion du
tube selon z k ? Torsion par
unité de longueur m
40
8. Propriétés observationnelles Vecteurs b et j
photosphériques Simulation et topologie de
sigmoïdes XUV
41
Formation dun sigmoïde en torsion continue

• SXR ? heat
• I intégration verticale du terme de chauffage
  Joule (à la louche !)
• Images ci-dessous pour Ntours 1.25 tours ltgt
  f 2.5 p

S
j 2 dz
étendus
moyens
concentrés
Sigmoïdes S-inversés
42
Jz photospherique
 Bras spiraux  pour Jz a lt 0 amax 0.25
Mm-1
Extrémités des sigmoïdes z0
correspondance a lt 0 / S-inversé
43
Lignes de champs torsadées ?
B lines around vortex centers Flux rope B lines
on sigmoïd ends B lines on brightest parts
Le sigmoïde Nest PAS le tube de flux torsadé
44
Lignes de champs du sigmoïde
B lines around vortex centers Flux rope B lines
on sigmoïd ends B lines on brightest parts
Le sigmoïde est présent dans des lignes CISAILLEES
45
Jz et bxy photospherique
Le magnétogramme force-free non-linéaire présente
des asymétries
46
3 propriétés de bxy photospherique
Ligne dinversion cisaillement
bxy non-radial existe dans les régions de forts
courants dans chaque polarité


centre des vortex en bxy décalé par rapport
aux centre des vortex sur bzmax
47
Conclusions
48
Conclusions 1 / 2 (pour les théoriciens)

• i Ca marche !
• Comportement similaire à Amari et al. (1996) et
  Törok Kliem (2003)
• Traitement de la frontière inférieure robuste
  (pas ou peu de couche limite)
• Frontières ouvertes utilisables ? diminution
  sensible de la taille du domaine
• Expansion en deux phases sans transition
  abrupte et SANS perte déquilibre
• Torsion à grande échelle quasi-statique continue
  umax(z0) 0.02 cA
• Réponse des grandes lignes de champ et expansion
  u(zzmax) 0.3 - 0.7 cA
• Expansion du flux torsadé TRES influencée par le
  flux externe faiblement torsadé
• Calcul de relaxation suggérant courbe
  déquilibre GENERIQUE à la Sturrock et al. (1995)
• Extension des vortex et du flux torsadé va à
  lencontre du kink
• Courbe déquilibre difficile à reconstruire
  perte déquilibre  numérique  chez Törok ?

49
Conclusions 2 / 2 (pour les observateurs)

• i Simulation dobservations X (à tester avec
  SOLAR-B)
• I terme intégré (optiquement mince) du terme
  Ohmique (température)
• Formes sigmoïdales sens du S compatible avec
  le signe des courants
• PAS le tube de flux torsadé
• Enveloppe de lignes cisaillées en 3D en forme de
  J sous-jacentes
• Magnétogrammes vectoriels synthétiques (à
  tester avec THEMIS SOLAR-B)
• Des profils de bu (z0) très simples résultent
  de structures complexes por bj (z0)
• Structures spirales pour jz
• Co-existence de courants direct de retour
• Régions de champ vectoriel non-radial où bxy //
  contours de bz
• Séparation entre bz(max) et régions de champ
  purement vertical