INSTABILITES ET BIFURCATIONS ASSOCIEES A
LA MODELISATIONDE LA GEODYNAMO
Instabilité convective
Modélisation
Système d’équations
Dans nos simulations, nous adoptons une géométrie simplifiée pour le noyau fluide terrestre. Celui-ci est en effet approché par une coquille sphérique (volume compris entre deux sphères concentriques). Ce choix nous amène donc à négliger les trois points suivants : – Les frontières noyau-manteau (CMB) et noyau interne-noyau externe (ICB) sont légèrement elliptiques, l’ellipticité de la CMB est d’environ 2.67 · 10−3 [71] et est plus faible pour l’ICB. – Ces deux frontières comportent des irrégularités de petites échelles [31]. – l’ICB se déplace en rayon à cause de la cristallisation de la graine solide. Par souci de comparaison avec la théorie, nos modèles numériques ne prendront en compte par la suite que la convection thermique comme source d’énergie (celle-ci est formellement équivalente à la convection solutale aux valeurs de diffusivité près). Dans nos calculs, le fluide est considéré comme Newtonien et Boussinesq. L’approximation Boussinesq signifie que le fluide est homogène et incompressible sauf pour la force d’Archimède o`u sa compressibilité est prise en compte au premier ordre. D’après le modèle PREM la compressibilité du mélange dans le noyau externe est d’environ 20%. Cette hypothèse d’incompressibilité est une hypothèse forte puisqu’elle nous amène à ne pas prendre en compte les variations de la densité avec la température et exclue le gradient thermique adiabatique présent dans le noyau fluide terrestre. Pour un fluide compressible si une particule fluide est moins dense que son environnement, elle va monter (direction radiale) vers des pressions plus faibles donc se dilater et se refroidir. Au contraire si une particule fluide est plus lourde, elle va descendre donc se contracter et s’échauffer. Il existe donc un gradient thermique tel que lorsqu’on déplace une particule fluide sans échange de chaleur elle reste en équilibre avec son environnement.
Géostrophie et modèle quasi-géostrophique
Considérons pour l’instant la première équation du système (2.22) sans le terme de for¸cage thermique, on a ∂u ∂t = −(u · ∇)u − ∇Π + ∆u − E −1 (ez ∧ u). (2.24) On sait que pour le noyau terrestre Ω ≃ 7.3 · 10−5 rad/s, avec re ≃ 3.5 · 106 m et ν ≈ 10−6m2 .s −1 (Stacey, 1992) on obtient un nombre d’Ekman E ≈ 10−15 . Pour un nombre d’Ekman aussi petit, on peut, en première approximation, négliger tous les termes devant le terme de Coriolis, sauf bien sˆur le gradient de pression qui lui est asservi afin de garantir ∇ · u = 0. On obtient l’équation dite de l’équilibre géostrophique E −1 (ez ∧ u) = −∇Π . (2.25) En prenant le rotationnel de cette équation, on obtient (ez · ∇)u = 0 ⇒ ∂u ∂z = 0 . (2.26) Cette contrainte constitue le théorème de Proudman-Taylor : dans un écoulement géostrophique la vitesse est invariante selon la direction de l’axe de rotation z. Pla¸cons nous en coordonnées cylindriques (s, ϕ, z), un écoulement géostrophique dans une sphère aura nécessairement us et uz nulles pour satisfaire ∇ · u = 0 et la non-pénétration aux bords. Prenons l’exemple d’un us non nul, étant invariant selon z, us sera aussi non nul à la frontière inférieure et supérieure. La non-pénétration à la paroi externe impliquera alors une vitesse verticale non nulle et opposée pour ces deux frontières. L’invariance selon z ne sera ainsi plus vérifiée pour uz. Les seuls mouvements respectant le théorème sont donc ceux pour lesquels seule uϕ est non nulle (mouvements angulaires ou zonaux), uϕ ne dépend évidemment pas de z et ne peut pas non plus dépendre de ϕ pour satisfaire la contrainte de divergence nulle, on a donc u = uϕ(s)eϕ. Le mouvement s’organise en cylindres concentriques, inclus dans la sphère et coaxiaux avec l’axe de rotation de la sphère, appelés contours géostrophiques (cf fig.2.1). Si maintenant on réintroduit le for¸cage thermique, on note d’emblée que les mouvements convectifs ne vont pas pouvoir respecter cette contrainte. En effet des mouvements purement zonaux ne permettent pas au fluide le transport radial de la température. Pour autant la géométrie de l’écoulement restera fortement influencée par la rotation, nous resterons donc pour la suite de notre étude en coordonnées cylindriques (s, ϕ, z). Bien que le mouvement obtenu ne soit plus exactement géostrophique, les expériences (Carrigan et Busse, 1983 ; Cardin et Olson, 1994) ont montré que la quasi-invariance selon z des mouvements convectifs est une caractéristique robuste de l’écoulement à bas nombre d’Ekman. Cependant les mouvements verticaux ne peuvent pas ˆetre z-invariant. Nous avons donc utilisé un modèle à deux dimensions dit modèle quasigéostrophique (QG) o`u seuls les mouvements horizontaux sont invariant selon z mais pas les mouvements verticaux, ce modèle, dit quasi-géostrophique sphérique fut introduit pour ce problème par Busse (1970). Cette hypothèse implique, pour satisfaire à un divergence nulle, une vitesse verticale linéaire en z. Nous supposons que les mouvements verticaux ont une variation à l’échelle de la sphère, et nous allons donc par la suite négliger les gradients verticaux devant les gradients horizontaux ( ∂ ∂z ≪ ∂ ∂s , ∂ ∂ϕ). L’équation de conservation de la matière devient donc ∇H · u ≃ 0 . (2.27) o`u l’indice H désigne les composantes horizontales. En coordonnées cylindriques le système (2.22) s’écrit ∂u ∂t + (u · ∇)u = −∇Π − E −1 (ez ∧ u) + ∆u + Raθ(ses + zez), ∂θ ∂t + (u · ∇)θ = Pr−1 (sus + zuz + ∆θ), (2.28) ∇ · u = 0. Prenons maintenant le rotationnel de la première équation de ce système, afin d’éliminer le gradient de pression. Il en résulte l’équation de la vorticité. Nous ne nous intéressons qu’à ω, la composante selon z de cette équation (ω = (∇ ∧ u) · ez), puisqu’elle ne met en jeu que les composantes horizontales de la vitesse, cette équation s’écrit ∂ω ∂t + (u · ∇)ω = ∆ω + (E −1 − ω) ∂uz ∂z − Ra ∂θ ∂ϕ. (2.29) Le deuxième terme du membre de droite de cette équation peut ˆetre simplifié. Il comporte en effet un terme en E−1 , qui est très grand puisque nous allons nous intéresser à la limite des petits nombres d’Ekman, et un terme en ω, qui à priori est d’ordre un. Ce dernier terme est donc négligeable devant le terme précédent. En moyennant ensuite en z l’équation (2.29) à l’aide de < >z= 1 2H R H −H dz, o`u H est la demi hauteur d’une colonne de fluide, il vient ∂ω¯ ∂t + (u · ∇)¯ω = ∆¯ω + E−1 2H [uz] H −H − Ra ∂ ¯θ ∂ϕ . (2.30) De plus, nous avons ¯ω = ω car ω, la composante selon z de la vorticité, est indépendante de z d’après l’hypothèse quasigéostrophique. Le terme en E −1 , qui correspond à l’injection ou au pompage de fluide aux extrémités inférieure et supérieure d’une colonne fluide, agit comme une source ou un puit de vorticité. La vitesse verticale uz a une double origine. La première est due à la conversion des mouvements radiaux en mouvements verticaux, du fait de la condition de non pénétration à la sphère externe. Cette condition de non pénétration aux parois implique en effet u · nsup = uscosλ + uzsinλ = 0 ⇒ uzsup = −ussup s H = −ussup s √ 1 − s 2 , u · ninf = uscosλ − uzsinλ = 0 ⇒ uzinf = usinf s H = usinf s √ 1 − s 2 , o`u λ correspond à la latitude à laquelle une colonne de fluide, de demi hauteur H et située au rayon s, touche la sphère externe. La normale supérieure (inférieure) à la sphère externe est désignée par nsup (respectivement ninf). La deuxième est due au pompage induit par la couche d’Ekman (voir la section suivante pour une définition). La vitesse verticale vp induite loin du bord par ce pompage a pour expression (Greenspan, 1968) vp = − 1 2 E 1 2 ez · ∇ ∧ n ∧ u + u p |n · ez| ! ez . (2.31) On extrait de cette formule l’expression de uzsup induite par le pompage uzsup = (vp · ez)(H) = − 1 2 E1/2 p |n · ez| ω. (2.32) 2.1. MODELISATION ´ 29 Le pompage induit à l’extrémité inférieure produit une vitesse uzinf opposée. Nous allons par la suite négliger ce terme de pompage devant le terme de viscosité dans l’équation de ψ. En effet le terme de pompage est en E − 1 2 et pour que ce terme devienne du mˆeme ordre que le terme viscosité, E− 1 2 ω ∼ ∆ω, il faut que l’échelle caractéristique l de variation de ω vérifie l ∼ E 1 4 . Or comme nous le verrons dans la section 2.2, la convection a pour échelle caractéristique E 1 3 , à cette échelle plus petite la viscosité est dominante. L’équation (2.30) s’écrit donc ∂ω ∂t + (u · ∇)ω = ∆ω − E −1 s 1 − s 2 us − Ra ∂ ¯θ ∂ϕ . (2.33) La deuxième équation du système (2.22) (équation de la chaleur) est elle aussi moyennée en z avec <>z ∂ ¯θ ∂t = −(u · ∇)θ + Pr−1 (sus + ∆¯θ). (2.34) Le terme zuz est quant à lui négligé puisqu’il implique une vitesse verticale. Notons que la géométrie de notre domaine a été modifiée lors de cette intégration. Le rˆole de la graine est ainsi étendu sur l’ensemble du cylindre tangent. Il en résulte un domaine d’intégration à deux dimension qui n’est plus simplement connexe. Nous négligeons donc tous les phénomènes convectifs à l’intérieur du cylindre tangent. Ceci est raisonnable pour la convection près du seuil. En effet, Roberts (1965) a calculé le seuil pour ce mode convectif, ce seuil est bien plus élevé que celui que nous allons étudier (en terme de nombre de Rayleigh). D’après l’hypothèse quasi-géostrophique on a un champ de vitesse u qui dans le plan (s, ϕ) est à divergence horizontale nulle. Celui-ci peut donc ˆetre exprimé en terme de fonction de courant ψ comme suit us = 1 s ∂ψ ∂ϕ et uϕ = − ∂ψ ∂s , d’o`u ω = −∆ψ + (∇ ∧ u0eϕ) · ez , (2.35) o`u u0 est le mouvement zonal moyenné en ϕ. La condition de non pénétration implique que la vitesse normale soit nulle à la frontière interne et externe, ψ est donc constante à la paroi. La fonction de courant étant définie à une constante près on peut, dans le cas d’un domaine simplement connexe, prendre celle-ci comme nulle. Or, notre domaine n’est pas simplement connexe et ψ ne peut donc pas ˆetre prise comme nulle aux deux frontières sans introduire un mouvement azimutal axisymétrique (ou vent zonal) arbitraire qui serait sans origine physique (voir [47]). Pour contourner cette difficulté, nous avons choisi de résoudre séparément les équations des mouvements axisymétriques et non-axisymétriques. Les équations sont résolues en modes de Fourier selon ϕ. Pour tous les modes m 6= 0, ψ peut effectivement ˆetre prise comme nulle aux deux parois (voir [2] pour une démarche similaire). 2.1.3 Couches d’Ekman et pompage Dans le cas de conditions aux limites de non-glissement la vitesse du fluide doit se raccorder à zéro à la paroi. Ceci a pour conséquence la présence d’une couche limite vis-queuse, le fluide étant en rotation cette couche limite est ici une couche d’Ekman. Alors que par définition dans notre modèle la force de Coriolis est dominante dans le coeur de la coquille sphérique, dans la couche d’Ekman c’est un équilibre entre la force de Coriolis et les forces visqueuses qui est réalisé (et bien sˆur le gradient de pression). Cette couche limite est d’épaisseur E1/2 sauf dans la région proche de l’équateur de la sphère d’extension E1/5 o`u elle est d’épaisseur E2/5 . La couche d’Ekman est importante pour l’écoulement dans son ensemble. En effet, le profil de vitesse horizontale créé dans la couche (la spirale d’Ekman) induit, pour que la conservation de la masse soit assurée, un écoulement normale à la couche appelé pompage d’Ekman. Cet écoulement secondaire créé par le pompage est non nul loin de la couche d’Ekman et influe sur la circulation générale du fluide (voir en annexe pour plus de détails sur la couche d’Ekman). Contrairement aux modes non-axisymétrique (m 6= 0), le mode axisymétrique (m=0) est de grande échelle en ϕ, nous allons donc établir l’équation des mouvements axisymétriques en tenant compte de ce pompage d’Ekman. En prenant la composante ϕ de l’équation de Navier-Stokes du système (2.22), on obtient ∂uϕ ∂t + (u · ∇)uϕ = − 1 s ∂Π ∂ϕ − E −1us + ∆uϕ. (2.36) En moyennant ensuite cette équation en ϕ et en z (us signifiant par exemple la moyenne en ϕ de us) il vient ∂uϕ ∂t + (u · ∇)uϕ = − E−1 2H Z H −H usdz + ∆uϕ − uϕ s 2 . (2.37) Le gradient de pression a disparu identiquement. Le terme de Coriolis est nul à l’ordre dominant d’après l’incompressibilité et en utilisant la formule de Green (D.33). Pour trouver l’expression du pompage il nous faut donc aller aux ordres supérieurs. On part de l’incompressibilité intégrée sur le volume, ce qui d’après le théorème de la divergence nous donne ZZ S u · dS = 0. (2.38) Considérons S la surface fermée formée par le cylindre de rayon si (Σ1), le cylindre concentrique de rayon s (Σ2) et la calotte sphérique trouée fermant au-dessus et en-dessous l’espace compris entre ces cylindres (Σ3). On applique le mˆeme moyennage que précédemment ZZ Σ1 us dS + ZZ Σ2 us dS + ZZ Σ3 [n · u] H −H dS = 0 . (2.39) Le premier terme est nul car us est nulle à la graine et donc nulle sur tout le cylindre de rayon si (us est indépendante de z). On obtient donc 2πs Z H −H usdz + ZZ Σ3 [n · u] H −H dS = 0 . (2.40) Le terme [n · u] H −H traduit l’existence d’une vitesse normale à la paroi engendrée par le pompage d’Ekman. Cette vitesse, pour le couvercle supérieur, s’exprime par (Greenspan, 1968) u(H) · n = − 1 2 E 1 2∇e ∧ n ∧ u + u p |n · ez| ! · n, (2.41) o`u ∇e n’implique que des dérivés horizontales. En injectant cette expression dans l’équation (2.40), il vient 2πs Z H −H usdz − E 1 2 ZZ Σ3 ∇e ∧ n ∧ u + u p |n · ez| ! dS = 0 . (2.42) Il nous faut maintenant appliquer le théorème du rotationnel (Stokes) sur le deuxième terme de cette égalité. Pour utiliser ce théorème il faut que le contour de la circulation soit connexe ce qui n’est pas le cas pour la calotte sphérique trouée que nous considérons. On peut cependant écrire ZZ Σ3 = ZZ Σtotale dS − ZZ Σtrou dS , (2.43) o`u Σtotale, la surface englobant le trou, et Σtrou sont deux surfaces connexes. En posant comme nulle la vitesse dans le cylindre tangent ce qui annule la contribution du trou, on obtient ZZ Σ3 ∇e ∧ n ∧ u + u p |n · ez| ! = I Γ (n ∧ u) · eϕ p |n · ez| sdϕ + I Γ u · eϕ p |n · ez| sdϕ , (2.44) o`u Γ est le contour de rayon s. Ce résultat n’est valable que pour des conditions aux limites de type non-glissement et pour des bords au repos dans le référentiel considéré. Il devra ˆetre adapté pour d’autres configurations (par exemple pour une super rotation de la graine). Le premier terme du membre de droite de l’équation précédente contient un produit mixte avec n la normale contenue dans un plan (s, z) et eϕ, il est donc nul, on obtient l’égalité suivante Z H −H usdz = E 1 2 uϕ p |n · ez| . (2.45) Le système d’équations à résoudre est le système quasi-géostrophique ∂ ∂t − ∆ ∆ψ = (u · ∇)ω + E −1 1 − s 2 ∂ψ ∂ϕ + Ra ∂ ¯θ ∂ϕ ∂ ∂t − ∆ − 1 r 2 u0 = − (u · ∇)uϕ 0 − E−1/2 2(1 − s 2) 3/4 u0 (2.46) ∂ ∂t − Pr−1∆ ¯θ = −(u · ∇)θ + Pr−1 ∂ψ ∂ϕ . Dans le cadre de notre modèle quasi-géostrophique et pour des conditions aux limites de contraintes libres le terme de pompage sur le vent zonal n’existe pas. En effet, de part l’hypothèse quasi-géostrophique le vent zonal est invariant selon z, la contrainte verticale est donc nulle et son raccordement à zéro aux parois ne crée pas de pompage. Dans le cas de conditions aux limites de non-glissement le pompage sur le vent zonal sera tout d’abord négligé puis réintroduit pour étudier son effet sur les caractéristiques de l’écoulement.
1 Introduction |
