Instabilité elliptique

Approche locale, approche globale

L’approche locale consiste à considérer de façon lagrangienne un écoulement nonconfiné perturbé par une onde plane. Cela revient à ne considérer que des perturbations de courtes longueurs d’onde. Cette approche, initiée par Bayly (1986) et Craik & Criminale (1986) pour l’instabilité elliptique, s’est progressivement raffinée jusqu’aux travaux sur  la méthode WKB dans le cadre de l’hydrodynamique de Friedlander & Vishik (1991) et Lifschitz & Hameiri (1991). L’approche locale permet d’obtenir les taux de croissance, souvent sous forme analytique explicite, et permet de prouver qu’une simple ligne de courant elliptique est intrinsèquement instable.
L’approche globale considère un écoulement de base, éventuellement confiné, et le décompose en ses modes normaux. Cette approche permet donc de tenir compte des effets de confinement et de suivre l’évolution de l’instabilité mode par mode. Il est possible d’étendre cette approche à une analyse faiblement non-linéaire : l’évolution temporelle de l’instabilité peut alors être calculée, mode par mode, depuis sa croissance initiale exponentielle jusqu’à sa saturation. Dans le cadre de l’instabilité elliptique, cette approche a été développée par Eloy et al. (2000, 2003) en géométrie confinée cylindrique, en se basant sur les modes de Kelvin (section 1.3.2.3). En géométrie confinée sphérique, Lacaze et al. (2004, 2005b) ont confirmé la pertinence d’une telle approche. Les instabilités inertielles résultant d’une résonance triadique de modes normaux n’apparaissent alors que si les conditions suivantes, appelées conditions de résonance, sont vérifiées :
m2 − m1 = m (2.11)
ω2 − ω1 = ω (2.12)
k2 − k1 = k (2.13) où mi , ωi et ki sont respectivement les nombre d’onde azimutaux, les pulsations et les nombres d’onde radiaux des deux modes normaux en résonance avec le forçage de nombre d’onde azimutal m, de pulsation m et de nombre d’onde radial k.

Approche numérique 

l’instabilité multipolaire par exemple, les conditions de résonance pour un élongationnel fixe dans le référentiel inertiel sont données par (e.g. Eloy & Le Dizès, 2001) : ω = 0 et k = 0, m étant l’ordre de la symétrie de la perturbation élongationnelle (m = 2 pour l’instabilité elliptique). Dans le cas de l’instabilité inertielle de précession en géométrie cylindrique (e.g. Lagrange et al., 2008), le forçage de précession donne m = 1 tandis que ω et k sont la pulsation et le nombre d’onde radial du mode de Kelvin forcé qui constitue l’écoulement de base (e.g. Lagrange et al., 2008, 2011).D’un point de vue numérique, la plupart des études de l’instabilité elliptique considèrent la dynamique 3D de deux vortex déformés initialement 2D (e.g. Lundgren & Mansour, 1996; Sipp & Jacquin, 1998; Lacaze et al., 2007; Roy et al., 2008). En effet, le travail pionnier de Pierrehumbert (1986) montre numériquement la présence de cette instabilité pour un écoulement plan associé à un vortex elliptique au sein d’une boîte, avec des conditions de glissement et une périodicité supposée dans la direction axiale : le problème linéaire aux valeurs propres est alors résolu par une méthode spectrale, et le taux de croissance de l’instabilité obtenu. Les études numériques en milieu confiné sont moins nombreuses. Mason & Kerswell (1999) utilisent un système de coordonnées elliptico-polaire non-orthogonal afin de résoudre l’écoulement dans un cylindre déformé par des méthodes spectrales ; l’évolution temporelle non-linéaire de deux modes différents de l’instabilité elliptique est alors calculée avec des conditions d’adhérence sur les parois latérales et des conditions de glissement sur les parois inférieure et supérieure du cylindre. Seyed-Mahmoud et al. (2000) obtiennent numériquement les fréquences et les taux de croissance de l’instabilité dans des ellipsoïdes et des coquilles ellipsoïdales en utilisant une méthode de Galerkin linéaire, en projetant l’écoulement sur une sélection d’ondes inertielles. Finalement, Ou et al. (2007) étudient la stabilité de domaines fluides ellipsoïdaux compressibles auto-gravitants et soulignent l’apparition de l’instabilité elliptique. Cependant, jusqu’ici, l’évolution non-linéaire visqueuse de l’instabilité elliptique au sein d’un conteneur n’a pas encore été simulée numériquement. Cette section présente une méthode numérique permettant de telles simulations. Le modèle numérique est validé par comparaison avec les résultats de la littérature sur un fluide en rotation au sein d’un ellipsoïde triaxial faiblement déformé.

Méthode numérique utilisée

Dans notre étude de l’instabilité elliptique, le problème à résoudre numériquement peut se décrire génériquement sous la forme suivante : au sein d’une géométrie nonaxisymétrique à parois fixes, quel est l’écoulement d’un fluide en rotation ? La configuration de référence dans ce travail est issue des études expérimentales dans un ellipsoïde faiblement déformé (Lacaze et al., 2004, 2005b; Herreman et al., 2009; Le Bars et al., 2010). Dans ces expériences, une cavité sphérique de rayon R, moulée dans un cylindre de silicone, est remplie avec du liquide et mise en rotation à vitesse angulaire constante Ω autour de l’axe (Oz). Le cylindre est alors légèrement compressé d’une quantité s selon (Ox), perpendiculairement à l’axe de rotation. La cavité devient alors un ellipsoïde triaxial constante le long des parois dans chaque plan perpendiculaire à l’axe de rotation Ω.d’axes (a, b, c) = (R − s, R + s, R), avec une ellipticité équatoriale β = (b2 − a2)/(a2 + b2) et une vitesse tangentielle constante le long de la paroi déformée, égale à ΩR à l’équateur. Une telle configuration est le modèle le plus simple d’un noyau planétaire liquide sans graine solide interne, au sein d’un manteau déformé par les forces de marées, avec une vitesse tangentielle constante. De façon similaire, notre modèle numérique de référence résout l’écoulement d’un fluide en rotation au sein d’un ellipsoïde d’axes (a, b, c) liés au référentiel (Ox, Oy, Oz), avec une vitesse tangentielle constante le long de la paroi dans chaque plan perpendiculaire à l’axe de rotation Ω (fig. 2.2). La longueur c de l’axe polaire peut être choisie indépendamment des autres longueurs a et b (avec b > a), ce que ne permettait pas le dispositif expérimental. De plus, l’axe de rotation de l’ellipsoïde peut être incliné par rapport à l’axe polaire, ce qui nous permettra d’étudier l’influence de l’obliquité en section 2.3.3 ou l’interaction avec la précession en section 2.5. Notons cependant que, par défaut, c sera égal à la moyenne des axes équatoriaux Req = (a+b)/2, et l’axe de rotation sera selon (Oz), de même que dans le dispositif expérimental. Dans toutes nos simulations, le fluide est initalement au repos, puis soudainement, une vitesse angulaire constante est imposée de sorte que la vitesse tangentielle le long de la paroi déformée dans chaque plan perpendiculaire à l’axe de rotation soit égale à Ω (a′ + b′)/2, où a′ et b′ sont les axes de la paroi elliptique dans ce plan. Nos résultats sont adimensionnés par l’échelle de longueur Req et l’échelle de temps Ω−1. Ainsi, cinq nombres adimensionnels permettent de caractériser le système : le nombre d’Ekman E = ν/(Ω R2eq), où ν est la viscosité cinématique du fluide, l’ellipticité β = (b2 − a2)/(a2 + b2) de la déformation, le rapport d’aspect c/b qui quantifie l’aplatissement de l’ellipsoïde, et finalement l’inclinaison θ et la déclinaison φ de l’axe de rotation. Le problème résolu numériquement est donc régi par le système d’équations ∂u ∂t + u · ∇u = −∇p + E △ u − 2 Ω∗ c × u , (2.14) ∇ · u = 0 , (2.15) où des conditions d’adhérence à la paroi sont utilisées pour le fluide. Notons que nous
travaillons toujours dans un référentiel où les parois de l’ellipsoïde sont fixes, ce qui correspond au référentiel inertiel de référence dans la plupart des cas. La force de Coriolis −2 Ω∗c × u est uniquement utilisée en section 2.3.2 où l’ellipsoïde est entièrement soumis à une rotation supplémentaire Ω∗ c ez.
Dans la plupart des travaux, les études numériques des noyaux liquides planétaires
approximent les noyaux par une sphère ou un sphéroïde, ce qui permet d’exploiter l’axisymétrie par décomposition spectrale. Les méthodes spectrales, à la fois rapides et précises,ne peuvent cependant pas être appliquées aisément à notre cas où l’axisymétrie est brisée.Nos simulations sont donc menées avec un code basé sur la méthode des éléments finis.Ce type de code est largement utilisé en ingénierie car des géométries complexes peuvent être considérées, ce qui constitue d’ailleurs un des points forts de ce type de méthode numérique. Ainsi, l’implémentation de notre ellipsoïde triaxial et des conditions aux limites nécessaires se fait assez simplement, le prix à payer étant un coût de calcul accrpar rapport à des méthodes spectrales. Avec le logiciel commercial COMSOL Multiphysicsr, un maillage non-structuréd’éléments tétrahédriques est créé. Les éléments utilisés sont des é ments standards de Lagrange P1 − P2, linéaires pour la pression et quadratiques pour
le champ de vitesses. Notons qu’aucune technique de stabilisation n’est utilisée dans ce travail. L’avance en temps est régie par un solveur IDA (Implicit Differential-Algebraic), basé sur un schéma BDF (Backward Differencing Formulas), décrit par Hindmarsh et al.(2005). A chaque pas de temps, le système est résolu avec le solveur direct linéaire pour matrices creuses PARDISO 1.L’instabilité elliptique induit une déstabilisation 3D des lignes de courant initalement 2D et elliptiques. Pour étudier ses propriétés globales, il est donc naturel d’introduire la valeur moyenne de la vitesse verticale en valeur absolue W =1VRRRV |w| dτ , avec w la
vitesse adimensionnelle verticale et V le volume de l’ellipsoïde. L’évolution typique de Wen fonction du temps est représenté en figure 2.3a pour E = 1/500 et β = 0.317. À t = 0,le fluide est au repos au sein de l’ellipsoïde, et la condition d’adhérence à la paroi met progressivement le fluide en rotation. Le premier pic sur W, juste après t = 0, provient du pompage d’Ekman (section 1.3.3) qui apparaît durant le spin-up (section 1.3.5) i.e. sur une durée de l’ordre du temps d’Ekman tE = E−1/2 Ω−1 , beaucoup plus court que le temps typique de diffusion visqueuse tv = R2
eq/ν = E−1 Ω −1 (Benton & Clark, 1974).En figure 2.3a par exemple, le temps adimensionnel d’Ekman donne ΩtE ≈ 22, en accord avec l’évolution de W. Une fois le fluide en rotation, l’écoulement est essentiellement bidimensionnel avec des lignes de courant elliptiques, ce qui correspond à l’écoulement de base de l’instabilité elliptique