Shom 3 Final

VALIDATION/DE
VELOPPEMENT DE PARAME
TRISATIONS DU ME
LANGE POUR DES MODE LES DE CIRCULATION GE
NE
RALE OCE
ANIQUE, A PARTIR DE SIMULATIONS DES GRANDES E
CHELLES (LES) NON-HYDROSTATIQUES. Marche EPSHOM-UJF No 00.87.070.00.470.29.25

RAPPORT SCIENTIFIQUE DE FIN D'E
TUDE (Premiere tranche) Periode: Juin 2001 - Fevrier 2002

Responsable Scienti que: Olivier METAIS Correspondant EPSHOM: Yves MOREL

1

Table des matieres

1 Fiche technique

1.1 Cadre des Travaux . . . . . . . . . . . . . . . . . 1.2 Renseignements Administratifs . . . . . . . . . . 1.2.1 Periode d'activite concernee par le rapport 1.2.2 Responsable scienti que . . . . . . . . . . 1.2.3 Correspondant EPSHOM . . . . . . . . . 1.2.4 Personnes ayant contribue aux travaux . . 1.3 Objet du rapport . . . . . . . . . . . . . . . . . . 1.4 Perspectives . . . . . . . . . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

2 Rapport d'activites scienti ques

2.1 Prise en main des outils numeriques . . . . . . . . . . . 2.1.1 Con guration de l'ecoulement . . . . . . . . . . 2.1.2 Equations regissant l'ecoulement . . . . . . . . . 2.1.3 Les termes de forcage . . . . . . . . . . . . . . . 2.1.4 Adimensionalisation des equations . . . . . . . . 2.1.5 Modelisation sous-maille . . . . . . . . . . . . . 2.1.6 Calculs d'ilustration . . . . . . . . . . . . . . . 2.2 Implementation des eets de la Salinite . . . . . . . . . 2.2.1 Nouveau systeme d'equations adimensionalisees 2.2.2 Ilustration des eects de salinite . . . . . . . . . 2.3 Implementation des conditions limites de ux impose . 2.3.1 Nouvelles conditions limites adimensionalisees . 2.3.2 Ilustration de la condition de ux impose . . . . 2.4 Implementation d'un ux de chaleur/salinite variable . 2.4.1 Nouvelles equations adimensionalisees . . . . . . 2.4.2 Ilustration des eects de ux variable . . . . . . 2.5 Implementation des eets du vent en surface . . . . . . 2.5.1 Nouvelles equations adimensionalisees . . . . . . 2.5.2 Ilustration des eets du vent de surface . . . . .

3 Manuel d'utilisation du code 3.1 3.2 3.3 3.4

Le Systeme nal d'equations adimensionalisees . Les conditions aux frontieres . . . . . . . . . . . L'etat de base . . . . . . . . . . . . . . . . . . . Discretization des equations . . . . . . . . . . . 3.4.1 La discretisation temporelle . . . . . . . 3.4.2 La discretisation spatiale . . . . . . . . . 3.4.3 Le systeme d'equations nal discretise . 3.4.4 Les conditions limites . . . . . . . . . . . 3.5 Description des routines . . . . . . . . . . . . . 2

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

3 3 3 3 3 3 3 5

6

6 6 7 8 9 11 12 12 17 17 19 19 20 21 21 24 24 26 27

29

29 30 32 32 32 34 39 41 41

1 Fiche technique

1.1 Cadre des Travaux

Marche EPSHOM-UJF No 00.87.070.00.470.29.25 sur la validation et le developpement de parametrisations du melange pour des modeles de circulation generale oceanique.

1.2 Renseignements Administratifs

1.2.1 P
eriode d'activit
e concern
ee par le rapport La premiere tranche de l'etude (9mois), de Juin 2001 a Fevrier 2002. La duree totale de l'etude est de 24 mois (9mois + 15mois).

1.2.2 Responsable scienti que Olivier Metais, Directeur de Recherche, CNRS Laboratoire des Ecoulements Geophysiques et Industriels BP53, 38041 Grenoble Cedex 9 Tel: 04-76-82-51-22 Telecopie: 04-76-82-52-71 Email: [email protected]

1.2.3 Correspondant EPSHOM Yves Morel

1.2.4 Personnes ayant contribu
e aux travaux Carlos B. da Silva, Bernard Barnier, Olivier Metais.

1.3 Objet du rapport

Les neuf premiers mois du present contrat ont ete consacres a l'implementation dans le code initial (Navier-Stokes dans l'aproximation de Boussinesq) des eets de la salinite, d'une condition de ux imposee (de chaleur/salinite), des ux de chaleur/salinite variables, et des eets du vent en surface. L'annexe scienti que (voir ci-dessous) detaille l'implementation, les nouvelles equations resolues par le code et quelques nouveaux calculs dans des con gurations interessantes. Les t^aches prevues pour la premiere tranche de ce contrat etaient les suivantes:

T^ache 1.1 termin
ee] Prise en main des outils num
eriques: Cette

t^ache correspond a la familiarisation avec l'outil numerique de depart (le code developpe par Padilla-Barbosa 6]) et a la revision bibliographique du sujet de la convection turbulente en milieu tournant en particulier dans le phenomene de formation des eaux profondes dans l'ocean. 3

T^ache 1.2 termin
ee] Impl
ementation des eets de la Salinit
e: Dans

cette t^ache les eets de la salinite ont ete inclus et valides dans le code de depart par l'ajout d'une equation de transport de cette variable. La densite nale de l'eau de la mer est calculee en prenant en compte la temperature et la salinite. L'analyse detaillee des eets de la salinite dans la physique de la formation des eaux profondes sera etudiee dans la deuxieme annee du contract. En particulier la relation lineaire entre densite, temperature et salinite implementee dans le code sera remplace par une relation nonlineaire, plus proche de la physique de l'ocean.

T^ache 1.3 termin
ee] Impl
ementation des conditions limites de ux impos
e: Cette tranche concerne le remplacement des termes de source dans

les equations de temperature et salinite par de nouvelles conditions limites dans ces equations, visant a imposer des ux donnes de chaleur/salinite. Ces nouvelles conditions limites ont ete validees et elles sont un grand pas en avant en ce qui concerne le realisme de la description physique de la formation des eaux profondes, notament dans la phase de preconditionement. Cependant, ces conditions limites peuvent ^etre encore ameliorees, car la viscosite turbulente a la surface, utilisee dans ces conditions, est tres dicile a calculer (la parametrisation LES n'est pas valable a la surface). Une autre facon d'imposer les ux de chaleur/salinite a la surface sera analysee dans le deuxieme annee du contrat.

T^ache 1.4 termin
ee] Impl
ementation d'un ux de chaleur/salinit
e variable: Un ux de chaleur/salinite variable dans le temps a ete implemente et valide dans cette t^ache. Ceci permetra d'etudier l'eet de plusieurs cycles de refroidissement dans la formation des eaux profondes.

T^ache 1.5 termin
ee] Impl
ementation des eets du vent en surface:

Dans cette t^ache les eets du vent en surface ont ete implementes et valides. Ceci peut ^etre acompli de deux facons dierentes: Soit par un terme de source dans les equations de quantite de mouvement, soit comme condition de ux impose a la surface.

T^ache 1.6 en cours] Impl
ementation d'une topographie de fond simple (avec conditions limites de inow/outow): Pour etudier les eets de la topographie du fond, il faut developper un nouvel outil de calcul. Plusieurs strategies possibles pour etudier ce probleme sont en cours d'analyse. Cette tache est la principale de la 2eme phase de l'etude.

4

1.4 Perspectives

La proposition initiale du contrat de cette etude prevoyait: Premiere annee: Formation des eaux profondes: Implantation dans le code SPECOMPACT d'un forcage realiste par le vent en surface implantation des eets de salinite validation du modele sur des cas idealises de convection profonde, et comparaisons des resultats avec les modeles turbulents simplies classiques. Deuxieme annee: Etude d'un courant gravitaire sur un talus continental: Transformation du code SPECOMPACT pour la prise en compte de geometries complexes d'ecoulements, et validation sur le cas idealise d'un ecoulement d'eau dense le long d'une pente topographique.

En accord avec ce contrat, la premiere annee a ete consacree a l'implementation/teste de plusieurs nouvelles equations qui rendent le modele numerique plus proche de la realite physique de l'ocean. Ce travail ouvre des portes a de nouvelles recherches sur la dynamique de la formation des eaux profondes. Le nouveau code de calcul permet a la fois d'etudier en detail la physique (a grande et petite echelle) de la formation des eaux profondes et de construire une base de donnees de simulations qui peut ^etre utilisee pour ameliorer les parametrisations utilisees dans les modeles de circulation oceanique. Ces travaux seront realises dans la deuxieme tranche du contrat. En particulier une comparaison entre les resultats LES et ceux obtenus avec la parametrisation KPP sera l'objet d'une etude speci que. Egalement dans cette 2eme tranche, un nouvel outil numerique capable de prendre en compte les eets d'une topographie de fond sera developpe et valide. Dans un premier temps, les con gurations topographiques realisees seront simples (talus de pente constante, variation lineaire avec la profondeur). Nous validerons le code dans des cas idealisees de l'ecoulement d'eau dense. A plus longue terme (dans le cadre d'etudes a venir) nous chercherons a etablir une nouvelle parametrisation pour decrire le melange le long d'une pente topographique.

5

2 Rapport d'activites scientiques Cette partie decrit en detail les equations, implementations et resultats des t^aches eectuees dans cette etude. Les t^aches sont les suivantes:

2.1 2.2 2.3 2.4 2.5

Prise en main des outils numeriques Implementation des eets de la Salinite Implementation des conditions limites de ux impose Implementation d'un ux de chaleur/salinite variable Implementation des eets du vent en surface

2.1 Prise en main des outils numeriques

Cette etude a comence par la prise en main des outils numeriques, et par la repetition de certains calculs de Padilla-Barbosa et Metais 7].

2.1.1 Con guration de l'
ecoulement

La con guration de l'ecoulement de l'etude est schematisee dans la gure 1. Le systeme est constitue d'une partie de l'ocean de largeur Lx = Ly = 32 km et de profondeur egale a H = 2 km. Comme condition initiale, l'ocean est considere comme etant au repos et possedant une densite initiale constante egale a 0 = 1000 kg:m;3. Un refroidissement localise en surface, sur un disque de rayon R = 8 km, fait declencher des mouvements de convection verticale. Le systeme est soumis a une rotation avec un taux xe par le parametre de Coriolis f . Bo g

Ω D = 2R

H ρ

Z

Lx = Ly Y

X

1 { Con guration des simulations numeriques visant a reproduire le phenomene de la convection des eaux profondes. Fig.

6

Dans les simulations de Padilla-Barbosa et Metais 7], le refroidissement local en surface correspond a un ux d'energie egal a Q = 800 Wm;2 . Ce refroidissement est applique pendant un periode de tforc = 48 h et le temps total de la simulation est de tsim = 96 h.

2.1.2 E
quations r
egissant l'
ecoulement

Dans ce systeme l'evolution des champs de vitesse est bien decrite par les equations de Navier-Stokes dans l'approximation de Boussinesq:

@u + @ (uu) + @ (uv) + @ (uw) = fv ; 1 @p0 +  @ 2 u + @ 2 u + @ 2 u @t @x @y @z @x2 @y2 @z2 0 @x

!

@v + @ (vu) + @ (vv) + @ (vw) = ;fu ; 1 @p0 +  @ 2 v + @ 2 v + @ 2 v @t @x @y @z @x2 @y2 @z2 0 @y

!

(1) (2)

! @w + @ (wu) + @ (wv) + @ (ww) = ; 0 g ; 1 @p0 +  @ 2 w + @ 2 w + @ 2 w (3) @t @x @y @z @x2 @y2 @z2 0 0 @z @u + @v + @w = 0 (4) @x @y @z Ou u, v et w sont les vitesses selon les directions x, y et z (voir gure 1).  est la viscosite cinematique f = 2 sin ' est le parametre de Coriolis. Le champ de densite est donne par, (x y z t) = 0 + 0 (x y z t) (5) ou 0 represente une densite de reference et 0 une petite uctuation autour de 0, 0 << 0 . Le champ de pression peut ^etre decompose en deux parties: une partie \hydrostatique" et une partie uctuante.

p(x y z t) = p0 (z) + p0 (x y z t) ou la pression hydrostatique ne depend que de la profondeur,

(6)

p0(z) = P0 ; 0 gz: (7) P0 est la pression a la surface. Dans le code de depart les eets de salinite ne sont pas pris en compte et la temperature est liee a la densite par la relation lineaire, 7

= 0 1 ; (T ; T0)] (8) ou  est le coecient d'expancion thermique de l'eau, et T0 est une temperature de reference. De ce fait l'equation classique d'energie devient une equation d'evolution de la densite,

! @ 0 + @ 0 + @ 0 + @ 0 =  @2 0 + @2 0 + @2 0 (9) @t @x @y @z @x2 @y2 @z2 qu'il faut ajouter au syteme d'equations (1), (2), (3) et (4), pour fermer le systeme. Une quantite importante caracterisant ce systeme est le ux de ottaison B0 , de nit comme, _ B0 = gCQ : (10) 0 w g est l'acceleration de la gravite, Q_ le ux de chaleur a la surface et Cw est la chaleur speci que de l'eau.

2.1.3 Les termes de forcage Dans la version initiale du code, les eets de refroidissement en surface sont inclus dans le probl^eme par un terme de forcage dans l'equation d'evolution de densite. Si on considere une couche de hauteur h << H , pres de la surface, la variation temporele de la temperature dans cette couche de uide est donnee par, @T = Q_  (11) @t 0 Cw h ou Q_ est le ux de chaleur en surface. Une fois que la densite et la temperature sont liees par une relation lineaire (equation (8)), l'evolution de la densite est donnee par, @ = ; Q_ : (12) @t Cw h Alors, l'equation pour la densite (9) devient,

! @ 0 + @ 0 + @ 0 + @ 0 =  @ 2 0 + @ 2 0 + @ 2 0 + Q_ (13) @t @x @y @z @x2 @y2 @z2 Cw h Le refroidissment en surface appara^$t ainsi comme une source qui augmente la densite des elements de uide dans une couche de uide d'epaisseur h pres de 8

la surface. A noter que ce forcage est applique seulement dans un disque de rayon R de ni par une fonction a, qui s'ecrit:

 h 2 2 i2 8 > 1 ; exp ; R;(xR=8+y ) < a=> :0 0

0

pour z = 0 ;0:05H

(14)

ailleurs La gure 2 montre a en fonction de la coordonnee horizontale y, pour une coupe verticale a la moitie du domaine (x = Lx=2).

Fig.

a.

2 { L'adoucissement de la distribution du forcage represente par la fonction

2.1.4 Adimensionalisation des
equations

Prenant comme parametres de reference la profondeur H et le ux d'energie Bo, les dimensions caracteristiques du systeme peuvent ^etre ecrites en fonction uniquement de H et Bo. La longueur, la vitesse et la gravite reduites, caracterisant le systeme, sont donnees par, 2= 3 B o lnorot = H unorot = (BoH ) (15) norot = H 1=3 : Le systeme d'equations adimensioneles a resoudre est donc donne par: 1=3

g0

@ u + @ (u u) + @ (u v) + @ (u w) = ; @  + 1 v @t @x @y @z @ x Ronorot ! 1 @2u + @2u + @2u + Re @ x2 @ y 2 @ z 2 9

(16)

@ v + @ (v u) + @ (v v ) + @ (v w) = ; @  ; 1 u @t @x @y @z @ y Ronorot ! 1 @2v + @2v + @2v + Re @ x2 @ y 2 @ z 2 @ w + @ (w u) + @ (w v) + @ (w w) = ; @  @t @x @y @z ! @z 2 2 2 1 @ w+@ w+@ w ; + Re @ x2 @ y 2 @ z 2

(17)

(18)

@u + @v + @w = 0 @x @y @z

(19)

! @ + @ + @ + @ = 1 @2 + @2 + @2 + ; @ t @ x @ y @ z RePr @ x2 @ y2 @ z 2

(20)

ou les caracteres en gras representent les variables (ou operateurs) adimensionnelles, de nies a partir des expressions suivantes:

t xi u = ui t = xi = lnorot i unorot lnorot =unorot

 = uP2 0 norot

=

0g 0 0 gnorot

Apres adimensionalisation, le terme de forcage devient, ; = lnorot (21) h a ou la fonction a est de ni par l'equation (14). L'adimensionalisation des equations fait appara^$tre des nombres adimensionnels, tels que le nombre de Reynolds Re, le nombre de Rossby Ro et le nombre de Prandtl Pr,

Re = unorotlnorot

unorot Ronorot = fl

norot

10

Pr =  :

2.1.5 Mod
elisation sous-maille La modelisation sous-maille utilise le modele de la Fonction de Structure (Lesieur & Metais 5]), tenant compte de l'hypothese de Boussinesq. La viscosite turbulente t (~x t) est calculee par,

t (~x 'c t) = 0:105Ck;3=2'cF2 (~x 'c t)]1=2 (22) ce qui permet de calculer le nombre de Reynolds turbulent, lnorot Ret (~x t) = unorot (23) t(~x t) Dans l'expression (22), F2(~x 'c t) est la fonction de structure d'ordre 2 des vitesses resolues, de nie comme, F 2(~x 'c t) = hk~u(~x t) ; ~u(~x + ~r t)k2ik~rk= (24) Dans la pratique, on n'a pas acces aux moyennes d'ensemble. Prenant en compte l'hypothese d'ergodicite, ou les moyennes d'ensemble et spatiales se confondent, et en supposant que l'ecoulement est localement isotrope, on remplace la moyenne d'ensemble par une moyenne spaciale locale. Ainsi <> est determinee sur les points voisins du point ~x. Par exemple, pour des mailles cubiques, la moyenne peut se faire sur les six points voisins, de la facon suivante: c

3 X F 2(~x 'c t) = 16 k~u(~x + 'xi~ei ) ; ~u(~x)k2 + k~u(~x ; 'xi~ei) ; ~u(~x)k2] (25)

i=1

Si le maillage est non-uniforme, ce modele nous permet de faire des interpolations. Soit 'c = ('x1 'x2 'x3 )1=3 , on aura alors:

 2=3 3 X F 2(~x 'c t) = 61 k~u(~x + 'xi~ei ) ; ~u(~x)k2 + k~u(~x ; 'xi~ei) ; ~u(~x)k2] ''xc i (26) i

i=1

ou 'xi~ei, designe le pas d'espace oriente dans la direction i par le vecteur unitaire ~ei. Le nombre de Prandtl turbulent est pris egale a: (27) Prt = t = 0:6 t Le systeme d'equations a resoudre devient alors,

@ ui + @ (u u ) = ; @  ; 1 " u i3k k @ t "@ xj i j @ xi Ronorot !#  + @@x Re +1 Re @@xui + @@uxj ; i3 j

t

j

11

i

(28)

@ ui = 0 (29) @ xi " # @ + @ ( u ) = @  1 + 1 @ + ; (30) j @ t @ xj @ xj RePr Ret Prt @ xj ou la barre sur les variables represente le ltrage spatial implicite dans la simulation des grandes echelles.

2.1.6 Calculs d'ilustration Certains des calculs menes par Padilla-Barbosa et Metais 7] ont ete refaits. Les gures 3 a) a f) montrent des isosurfaces de basse temperature d'un calcul correspondant a f = 1:0
10;4s;1 (H4) a plusieurs instants. Les gures 3 montrent la formation d'un reservoir thermique (3 a)) qui croit jusqu'a devenir instable ce qui provoque la formation des plumes d'eau froide (3 b)) qui descendent dans l'ocean jusqu'a attendre le fond (3 c)). Pendant ce temps, le uide froid reste dans un cylindre ou des mouvements quasi-verticaux dominent l'ecoulement. Une fois que le uide froid est arrive au fond de l'ocean, on observe des mouvements de transport horizontal, domines par la rotation, qui nissent par melanger le uide dans les directions horizontales ( gures 3 d) a f). Les gures 4, 5 et 6 montrent des details des vitesses dans l'ecoulement dans la phase de melange violent cause par la formation des plumes d'eau froide.

2.2 Implementation des eets de la Salinite

L'implementation des eets de la salinite se fait en ajoutant une equation de transport de cette variable (S ). Il faut aussi remplacer l'ancienne equation de transport de la densite , (equation (30)) par une equation de transport de la temperature (T ). La densite nale de l'eau de mer est calculee en considerant la relation lineaire, = 0 1 ; (T ; T0) +  (S ; S0)] (31) au lieu de l'equation (8). S0 est la valeur de reference pour la salinite et  est la valeur du coecient de contraction saline de l'eau. La prise en compte des eects de salinite intervient aussi dans la de nition du ux de ottaison B0 . Au lieu de l'equation (10) on a maintenant, _ E ; P )} B0 = gCQ + gS s({z (32) | 0 w | {z } B2

B1

ou E ; P est la dierence entre l'evaporation et la precipitation, et Ss est la salinite en surface. Le terme B1 correspond a la contribution d^ue aux pertes 12

a)

d)

b)

e)

c)

f)

3 { Isosurfaces de basse temperature dans un calcul de convection des eaux profondes avec f = 1
10;4s;1 . Le seuil de temperature correspond a T 0 = ;0:006K . a) Apres 3 heures( b) Apres 9 heures( c) Apres 24 heures( d) Apres 48 heures( e) Apres 72 heures( f) Apres 96 heures. Fig.

13

a)

b)

4 { Isosurfaces de vitesse verticale dans un calcul de convection des eaux profondes avec f = 1
10;4s;1. Le seuil des isosurfaces de vitesse correspond a uz = ;0:02ms;1 (vert) et uz = +0:02ms;1 (gris). a) Apres 9 heures( b) Apres 24 heures. 14

Fig.

5 { Detail de l'ecoulement des eaux profondes avec f = 1
10;4s;1 apres 24 heures. La gure montre des vecteurs de vitesse dans un plan (x y) traversee par des plumes a basse temperature indiquees par des isosurfaces grises. Le seuil de temperature est T 0 = ;0:02K . Les vecteurs vitesse sont dans un plan a z = ;1000m de profondeur. Les courbes d'iso-vorticite sont indiques dans le m^eme plan. Fig.

15

a)

b)

6 { Composante horizontale des vecteurs de vitesse a 1000m de profondeur dans un calcul de convection des eaux profondes avec f = 1
10;4s;1. a) Apres 24 heures( b) Apres 96 heures. Fig.

16

d'energie en surface et B2 est d^ue aux eets du ux de salinite. On de nit aussi la constante f telle que, 2 f = B (33) B0 :

2.2.1 Nouveau systeme d'
equations adimensionalis
ees

Dans le nouveau systeme d'equations a resoudre on remplace l'equation (30) par,

" # @ T + @ (T u ) = @  1 + 1 @ T ; (1 ;  ) ; (34) j f @ t @ xj @ xj RePr RetPrt @ xj et on ajoute l'equation de transport de salinite, " # @ S + @ (S u ) = @  1 + 1 @ S +  ;: (35) f j @ t @ xj @ xj RePr Ret Prt @ xj La temperature T et la salinite S adimensionelles sont de nies par, 0 gS 0 T = ggT S = (36) 0 0 gnorot norot ou T 0 et S 0 representent respectivement, l'ecart de la temperature et salinite par rapport a ses valeurs de reference: T = T0 + T 0 S = S0 + S 0 (37) Alors, en utilisant l'equation (31) la densite adimensionelle est obtenue simplement avec, = ;T + S : (38) A noter qu'on suppose aussi que la difusivite thermique et de salinite sont egales: T = S = .

2.2.2 Ilustration des eects de salinit
e

Certaines images du calcul correspondant au cas f = 4
10;4s;1, sont montrees sur la gure 7. Ceci montre la plongee des eaux froides dans le processus de convection et leur expansion au fond de l'ocean. Dans le cas de gauche (a),b),c)) il n'y a pas d'eet de salinite, tandis que dans le cas de droite (d),e),f)) l'eet de la salinite est ajoute en considerant le taux de evaporation moins precipitation (E ; P ) egale `a 20 mm:jour;1. Les images montrent que l'ajout de la salinite cause une forcage additionnel qui accelere la dynamique des plongees d'eau. Ceci 17

a)

d)

b)

e)

c)

f)

7 { Isosurfaces de basse temperature dans un calcul de convection des eaux profondes. Dans le cas de gauche (a),b),c)) il n'y a pas d'eets de salinite tandis que dans le cas de droite (d),e),f)) l'equation de transport de salinite etait prise en compte. Le seuil de temperature correspond a T 0 = ;0:007K . a),b) Apres 9 heures( c),d) Apres 24 heures( e),f) Apres 72 heures. Fig.

18

8 { Isosurfaces de salinite dans un calcul de convection des eaux profondes apres 9 heures. Le seuil de salinite correspond a S 0 = +0:00026 = . Fig.

peut ^etre observe par la vitesse a laquelle les deux ecoulements se developpent. Celui de droite evolue beaucoup plus rapidement. En outre la salinite, qui est un scalaire passif, exhibe des isosurfaces similaires a celles de la temperature. Ceci peut ^etre compris en comparant la temperature et la salinite au m^eme instant (voir gures 7 d) et 8). Dans cette implementation des eets de la salinite, on a considere que les coecients d'expansion thermique  et de contraction de salinite  sont constants dans l'equation (31). Dans la realite, ces coecients varient de facon complexe avec, entre autres, la profondeur, ce qui rend le probleme fortement non-lineaire. Il serait donc interessant d'etudier ses eets dans la dynamique des eaux profondes. Ce probleme sera analyse dans une prochaine etude.

2.3 Implementation des conditions limites de ux impose 2.3.1 Nouvelles conditions limites adimensionalis
ees

Dans la version initiale du code les eets de la perte d'energie en surface sont inclus dans les termes de source de l'equation (30). Dans une premiere phase, on a utilise la m^eme procedure pour les equations (34) et (35). Une condition limite plus realiste consiste a eliminer ces termes de source et a imposer un ux d'energie et de salinite a la surface (z = ;H ): 0 _ ( + t ) @@z z=;H = ; CQh ; (39) w 0  Q_ ;  = ; (40) ( + t ) @T @z z=;H 0 Cw 19

0

 ( + t ) @S @z z=;H = Ss(E ; P ) ; Sous la forme adimensionelle ces equations s'ecrivent: @  @ z z=;H = + (RePr + Ret Prt) ; @ T  @ z z=;H = ;(1 ; f ) (RePr + Ret Prt) ; @ S  @ z z=;H = +f (RePr + Ret Prt) ;

(41) (42) (43) (44)

2.3.2 Ilustration de la condition de ux impos
e Pour illustrer les dierences entre les deux cas, les gures 9 montrent des contours de temperature dans deux calculs de convection des eaux profondes, utilisant des termes de source dans les equations de temperature et salinite (gauche) et en imposant un ux donne de temperature et salinite (droite). Dans ces calculs, la viscosite cinematique de l'eau de mer est pris a  = 1:0m2:s;1 ce qui rend l'ecoulement laminaire et permet d'analyser l'evolution de la temperature en detail. Ceci permet en plus d'eliminer d'autres problemes. En fait si rien n'est fait pour le prevenir, en debut de calcul les derivees selon la direction verticale peuvent ^etre beaucoup trop elevees, ce qui peut causer des problemes numeriques (Gibbs). Ceci resulte du fait que la viscosite turbulente soit presque nulle en debut de calcul, et par consequent, les gradients de vitesse en surface selon la direction verticale doivent ^etre anormalement eleves, si on veut imposer des ux de chaleur/salinite importants. Pour eliminer ce probleme, en debut de calcul les uxs de chaleur/salinite doivent ^etre imposes avec une fonction croissante dans le temps. La m^eme procedure est utilisee couramment dans la litterature des eaux profondes. Les premieres images montrent des diferences dans la forme des isolignes de temperature pres de la surface. Dans le cas d'une source (negative) de chaleur la temperature reste presque constante dans une couche de uide pres de la surface (voir gure 9 a)). Un gradient de temperature s'etablit en dessous de cette premiere couche. Utilisant une condition limite de ux impose, on observe au m^eme instant qu'un gradient vertical de temperature constant s'etablit pres de la surface (voir gure 9 e)). En consequence, la suite des evenements est aectee dieremment dans les deux cas (voir les gures 9 b),c),d) et f),g),h)). Dans le cas d'une source de chaleur (gauche), on voit que le champ de temperature pres de la surface reste presque constant, provocant la plongee des eaux denses en agissant dans la partie de uide qui l'est en dessus. Dans le cas de ux impose (droite) la dynamique des 20

plongees d'eaux semble se declencher plus t^ot a cause des instabilites au bord de la cheminee. Les gures 10 montrent des isosurfaces de vitesse verticale negative/positive a plusieurs instants pour ces deux calculs. On observe la formation de cellules de Rayleigh-Benard de forme cylindrique, en dessous du disque de refroidissement. On voit que, en ce qui concerne la vitesse verticale la topologie de l'ecoulement dans les deux cas est similaire, mais le cas de ux impose se developpe plus rapidement. Sur l'implementation des conditions limites de ux impose, on a obtenu des resultats plus realistes qu'auparavant, en ce qui concerne le gradient de temperature proche de la surface. Cependant, un probleme resulte du fait que la viscosite turbulente a la surface non seulement n'est pas connue au debut du calcul, mais il n'est pas possible de determiner a partir du modele de LES, car la physique de la turbulence a la surface est completement diferentte et malleureusement tres mal comprise. De facon a decrire encore mieux la realite physique de ce probleme, des solutions alternatives, que dispensent le calcul de la valeur exacte du coecient de difusion thermique turbulent t a z = ;H , seront analysees.

2.4 Implementation d'un ux de chaleur/salinite variable

Le fait d'imposer un ux ou une source de chaleur/salinite variable dans le temps, cause seulement des petites modi cations dans le systeme d'equations a resoudre. On de nit les quantites, Z Tq Z Ts 1 1 _ Q(t)dt Is = T E (t) ; P (t)]dt (45) Iq = T q 0 s 0 ou Tq et Ts sont les temps caracteristiques des forcages appliques respectivement en chaleur et salinite1. Le ux de ottaison devient alors, q B0 = gI + gS sIs} {z | C (46) 0 w | {z } B2

B1

2.4.1 Nouvelles
equations adimensionalis
ees Les equations de temperature et salinite doivent etre changees pour prendre en compte les eets du ux de chaleur/salinite variable. On de nit d'abord, les quantites suivantes: _ ;qt (t) = QI(t) ;st(t) = E (t) I; P (t) (47) q s 1

Typiquement, Tq et Ts sont les durations des forcages de chaleur et de salinit
e

21

a)

e)

b)

f)

c)

g)

d)

h)

9 { Contours de temperature dans deux calculs de convection des eaux profondes. Le parametre de Coriolis est f = 4
10;4s;1 . Utilisant des termes de sources dans les equations de temperature et salinite (a),b),c),d))( Utilisant des conditions limites de ux impose en surface (e),f),g),h)). a),e) Apres 6 heures( b),f) Apres 12 heures( c),g) Apres 25 heures. d),h) Apres 31 heures. Fig.

22

a)

b)

c)

d)

e)

f)

10 { Isosurfaces de vitesse verticale dans deux calculs de convection des eaux profondes utilisant des termes de sources dans la temperature (a),b)) ou des conditions limites de ux impose (c),d),e),f)). Le seuil de vitesse correspond a w = ;0:002m:s;1 (gris) et w = +0:002m:s;1 (argent). Apres 8 heures (c))( Apres 16 heures (d))( Apres 25 heures (a),e))( Apres 31 heures (b),f)). Fig.

23

Utilisant les termes de forcage dans les equations de temperature (34) et salinite (35), ces termes deviennent alors,

Equation de temp
erature : ; (1 ; f ) ; ;qt (t) Equation de salinit
e : + f ; ;st(t):

(48)

Dans le cas ou on utilise des conditions limites de ux imposes, les anciennes conditions de ux exprimees dans les equations (44), deviennent,

@ T  q = ; (1 ;  f ) (RePr + Ret Prt ) ; ;t (t) z = ; H @z @ S  s = +  f (RePr + Ret Prt ) ; ;t (t) z = ; H @z

(49) (50)

2.4.2 Ilustration des eects de ux variable Pour montrer les eets d'un forcage variable, deux calculs dierents ont ete realises, avec un forcage constant ( gures 11 a) a d)) et un forcage variable ( gures 11 e) a f)). L'evolution temporel du forcage utilisee dans les deux calculs est montree dans la gure 12. Dans les deux calculs, la quantite de chaleur retiree a la surface est la m^eme. Comme prevu, pour les premieres heures du calcul, l'ecoulement se developpe plus rapidement dans le cas d'un forcage constant, car le ux d'energie au depart est plus important. Cependant, apres 45 heures, les deux ecoulements semblent ^etre arrives au m^eme etat, d^u a un plus fort ux de chaleur applique dans le cas de ux variable apres  10h.

2.5 Implementation des eets du vent en surface

Les eets du vent se traduisent par des tensions de cisaillement (x y ) a la surface de l'eau de la mer. Ceci peut ^etre implemente soit par un terme de source dans les equations de transport de quantite de mouvement, soit en imposant les tensions de cisaillement dans la condition limite de ces equations a la surface. Dans le premier cas, les equations pour u et v deviennent, au lieu de (51) et (52),

@u + @ (uu) + @ (uv) + @ (uw) = fv ; 1 @p0 @t @x @y @z ! 0 @xs 2 2 2 + @@xu2 + @@yu2 + @@zu2 + xh 0

@v + @ (vu) + @ (vv) + @ (vw) = ;fu ; 1 @p0 @t @x @y @z 0 @y 24

(51)

a)

b)

c)

d)

e)

f)

11 { Isosurfaces de basse temperature dans un calcul de convection des eaux profondes utilisant une source de chaleur constante (gauche) ou une source de chaleur variable (droite). Le seuil de temperature correspond a T 0 = ;0:006K . a),b) Apres 5 heures( c),d) Apres 15 heures( e),f) Apres 45 heures. Fig.

25

12 { Variation temporel de la quantite de chaleur retiree a la surface dans les deux cas de l'etude. Fig.

! s 2 2 2 @ v @ v @ v + @x2 + @y2 + @z2 + yh 0

(52)

xs (x y) at ys(x y) sont les tensions de cisaillement qui resultent de l'interaction ocean/atmosphere. Dans le deuxieme cas, les conditions limites pour les vitesses u, et v sont,  xs(x y) ( + t ) @u  = (53) @z z=;H 0   s(x y) (54) ( + t ) @v z=;H = y @z 0

2.5.1 Nouvelles
equations adimensionalis
ees Dans la pratique, les tensions de cisaillement dues au vent en surface sont appliquees dans le disque de refroidissement. On de nit alors les quantites adimensionelles, s (x y ) ;1 = xhg a (55) 0 0 norot  s(x y) ;2 = yhg0 a (56) 0 norot et les equations de quantite de mouvement, ltrees et adimensionalisees sont, (au lieu de (28)) 26

@ ui + @ (u u ) = ; @  ; 1 " u i j norot i3k k @ t @ x @ x Ro j i " @ ui @ uj !# 1 @ + @ x Re + Re @ x + @ x ; i3 + ;i( i1 + i2 ) (57) j t j i Si les eets du vent sont mis en tant que ux de quantite de mouvement les conditions limites sous la forme adimensionelle, s'ecrivent,  1 + 1 @ u z=;H = ;1 (58) Re Ret @ z 1 @ v  1 2 (59) Re + Ret @ z z=;H = ;

2.5.2 Ilustration des eets du vent de surface

Pour ilustrer les eets du vent en surface trois calculs ont ete menes (voir gure 13). Dans le premier on impose une tension de cisaillement d^u au r vent egale a (x y ) = (0:042 0)N:m;2. Ceci correspond a une vitesse de u = jx0 j = 0:0065m:s;1. Le deuxieme cas correspond a (x  y ) = (0:42 0)N:m;2. Comparant les isosurfaces de ces deux cas a 40 heures ( gures 13 a) et b)) on observe presque le m^eme comportement au niveau des grandes echelles. Il semble que le cas du vent plus intense ( gure 13 b) soit legerement plus developpe, ce qui indique que le vent en surface agis surtout au niveau des mouvements a petite echelle, en augmentant le melange dans la phase de preconditionement. Par contre, lorsque le vent en surface est tres intense (voir gure 13 c) on observe que la cheminee est convectee dans la direction du vent. La gure 13 d) montre les vecteurs vitesse a 1000 m de profondeur dans le cas de vent tres intense. On observe une direction preferentielle des vitesses dans la direction du vent et davantage de mouvements a petite echelle, dans la region de la cheminee, que pour le cas de reference. La physique a grande et petite echelle causee par les eets du vent en surface sera analysee en detail dans les travaux futurs.

27

a)

b)

c)

d)

13 { Isosurfaces de basse temperature dans un calcul de convection des eaux profondes utilisant une forcage du vent en surface apres 40 heures. Le seuil de temperature correspond a T 0 = ;0:006K . a) Avec (x y ) = (0:042 0) N:m;2( b) Avec (x y ) = (0:42 0) N:m;2 ( c) Avec (x y ) = (4:42 4:42) N:m;2. d) Composante horizontale des vecteurs de vitesse a 1000m de profondeur correspondant au cas c). Le parametre de Coriolis est f = 4
10;4s;1. Fig.

28

3 Manuel d'utilisation du code

3.1 Le Systeme nal d'equations adimensionalisees

Le systeme etudie correspond a une situation de convection profonde generee par une source de ottaison isolee Bo, a travers une surface circulaire de rayon R, sous l'inuence d'un taux de rotation constant f . La profondeur du domaine est H , et les dimensions horizontales, Lx = Ly (voir gure 1). Ce systeme est bien decrit par les equations de Navier-Stokes, dans l'approximation de Boussinesq. Sous la forme adimensionnelle, ces equations s'ecrivent:  Conservation de la quantite de mouvement

@ ui ; " u  = ; @  ; 1 " u @ t ijk j k !# @ xi Ronorot i3k k " + @@x Re +1 Re @@xui + @@uxj ; i3 + ;i( i1 + i2 ) (60) j t j i  Conservation de la masse @ ui = 0 (61) @ xi  Conservation de l'energie Si les eets de salinite ne sont pas pris en compte on resoud l'equation de conservation de l'energie dans la forme, " # @ + @ ( u ) = @  1 + 1 @ + ; ;q (t) (62) j t @ t @ xj @ xj RePr Ret Prt @ xj Alors que si on prend en compte les eets de salinite l'equation de conservation de l'energie est resolue dans sa forme originale,

" # @ T + @ (T u ) = @  1 + 1 @ T ; (1 ;  ) ; ;q (t) (63) f j t @ t @ xj @ xj RePr Ret Prt @ xj et on ajoute l'equation de transport de salinite, " # @ S + @ (S u ) = @  1 + 1 @ S +  ; ;s(t): j f t @ t @ xj @ xj RePr Ret Prt @ xj et la densite adimensionel est obtenu par, =T +S 29

(64) (65)

Les caracteres en gras representent les variables (ou operateurs) adimensionnelles de nies a partir des expressions suivantes:

~x = H~x

 2=3 t = T t = BH t o

~u = U~u = (BoH )1=3u~ Bo2 g H

0= o

T0 =

1 Bo2 g H

!1=3

!1=3

T

P = (BoH )2=3 S0 =

1 Bo2 g H

!1=3

S

Le ux de ottaison est de nit par, q + gSsIs B0 = gI 0 Cw avec, Z Tq Z Ts Iq = T1 Q_ (t)dt Is = T1 E (t) ; P (t)]dt q 0 s 0 Finalement, on a de nit aussi les variables suivantes, sIs f = gS B0 ; = Ha h

_ ;qt (t) = QI(t) ;1 =

q s x (x y)a 0 0 hgnorot

;st(t) = E (t) I; P (t) s s  (x y)a ;2 = y hg0 : 0 norot

(66) (67)

(68)

Les nombres adimensionnels sont respectivement, le Reynolds, Rossby, et Prandtl,

Re = UH 

U Ret (~x t) =  UH Ro = fH t (~x t)

3.2 Les conditions aux fronti eres

Pr =  Prt(~x t) =  (~x t)

Les conditions aux frontieres du systeme sont les suivantes: 30

Direction horizontale x: Condition de glissement sans frottement, et parois

verticales adiabatiques: A x = 0 1:

ux = 0

@ uz = @ uy = 0 @x @x

@ =0 @x

@T = 0 @x

@S = 0 @x

(69)

Direction horizontale y: Les frontieres y = 0 Ly presentent des conditions de periodicite. Cette periodicite se traduit par une condition a veri er par toutes les variables dynamiques a chaque instant, a savoir:

f (x 0 z) = f (x Ly  z)

(70)

Direction verticale z: Pour les conditions aux frontieres horizontales, on a

considere deux situations. La premiere correspond au cas ou on applique un forcage par un terme de source dans les equations de transport de densite, ou de temperature et salinite. Dans ce cas on a impose des condition limites de glissement sans frottement et parois horizontales adiabatiques, soit: A z = 0 ;1: uz = 0 @@uzz = @@uzy = 0 @@z = 0 @@Tz = 0 @@Sz = 0 (71) Dans le cas de ux imposee, les termes de source dans les equations de densite, temperature et salinite sont mis a zero, et les conditions limites deviennent, A z = 0:

@ uz = @ uy = 0 @z @z @ = + (RePr + Re Pr ) ; ;q (t) t t t @z @ T = ;(1 ;  ) (RePr + Re Pr ) ; ;q (t) f t t t @z @ S = + (RePr + Re Pr ) ; ;s(t) f t t t @z

uz = 0

et A z = ;1: uz = 0 @@uzz = @@uzy = 0

@ =0 @z 31

@T = 0 @z

@S = 0 @z

3.3 L'etat de base

L'etat de base du systeme consiste en un champ de vitesse, et des deviations de densite, temperature et salinite nulles partout:

u~ (~x t = 0) = ~0

(~x t = 0) = 0

T (~x t = 0) = 0 S (~x t = 0) = 0

3.4 Discretization des equations 3.4.1 La discr
etisation temporelle

La discretisation des derivees temporelles de nos equations est realisee par un schema explicite de Runge-Kutta (Williamson 8]) d'ordre trois, aussi bien pour les termes d'advection que les termes de diusion. L'obtention du champ de pression est decouplee du champ de vitesse a l'aide d'un schema a pas fractionnaire. Ainsi, le calcul est fractionne en deux etapes. Dans un premier temps, on obtient le champ de vitesse independamment du champ de pression. Posterieurement, on corrige le champ de vitesse par un champ de vitesse incompressible a l'aide d'un gradient de pression, solution de l'equation de Poisson. A chaque pas de temps p, les equations discretisees s'expriment par:

u~ p ; u~ p;1 = pA(u~ p;1 p;1) + pA(u~ p;2 p;2) ; r~ p (72) 't
p ;
p;1 = pB (u~ p;1
p;1) + pB(u~ p;2
p;2) (73) 't r~  u~ p = 0 (74) ou le sous-pas de temps p peut ^etre p = 1 2 3 et
represente , T ou S . Notons que p ; 2 est ignore lors d'une valeur de p = 1. A deux pas de temps consecutifs n et n + 1 on a donc u~ pi =0 = u~ ni et u~ pi =3 = u~ ni +1 avec et

 1  1 1 ~ A(u~  ) = u~
 + Ro u~
~ez + r: Re + Ret ru~ ; ~ez

B(u~ 
) = ;u~  r~
+ r:

 1  1 RePr + Ret Prt r
+ 

(75) (76)

La precision en temps de troisieme ordre des termes non-lineaires est obtenue en utilisant les coecients p et p, tels que (Williamson 8]):

32

1 = 0

1 = 158

2 = 125 2 = ; 1760 3 = 43 3 = ; 125 L'operation visant a incompressibiliser le champ de vitesse est executee a chaque sous-pas de temps. Cette procedure a ete proposee par Kim & Moin 2] et posterieurement generalisee a une methode du type Runge-Kutta par Le & Moin 3]. Elle est realisee en deux etapes a chaque sous-pas de temps. Tout d'abord, elle resout l'equation de conservation de quantite de mouvement sans considerer le champ de pression, a savoir:

u~p ; u~ p;1 = pA(u~ p;1 't

p;1 ) +  p

A(u~ p;2

p;2 )

(77)

De m^eme elle resout l'equation de conservation de l'energie:


p ;
p;1 = pB (u~ p;1
p;1) + pB(u~ p;2
p;2) 't

(78)

r~  u~ p = 0

(79)

Generalement, a cette etape, la condition d'incompressibilite:

n'est pas veri ee. La deuxieme etape consiste a reprendre le champ u~ p des vitesses, a n d'obtenir le champ corrige de vitesse u~ p, a partir de l'expression suivante: u~ p ; u~ p = ;r~ p (80) 't A partir de l'equation (80) et de la condition de conservation de masse (eq. 74), on peut deduire une equation de Poisson discrete, a partir de laquelle on obtient le champ  qui est posterieurement utilise pour eectuer la correction du champ de vitesse de l'expression (80). Cette expression s'ecrit: 1r ~  u~ p = r2  (81) 't Dans la pratique, la discretisation temporelle du code de calcul utilise resoud les equations de Navier-Stokes dans l'ordre suivant: 1. Dans un premier temps, on obtient une vitesse intermediaire u~ p sans prendre en compte, le champ de pression, a savoir:

u~p ; u~ p;1 = pA(u~ p;1 't 33

p;1 ) + 

p

A(u~ p;2

p;2)

(82)

2. On determine le scalaire
, de l'equation de conservation de l'energie:


p ;
p;1 = pB(u~ p;1
p;1) + pB (u~ p;2
p;2) 't

(83)

3. On obtient une vitesse u~ p, a partir de:

u~ p ; u~ p = ;r~ 0p 't

(84)

4. Une derniere etape consiste a rendre incompressible le champ u~ p a n d'obtenir le champ nal u~ p, a partir de l'expression:

u~ p ; u~ p = ;r~ p 't

(85)

Les scalaires  et 0 sont obtenus a partir de: r2 0p = '1t r~  u~ p r2 p = '1t r~  u~ p

(86) (87)

3.4.2 La discr
etisation spatiale

Pour les derivees premieres et secondes dans la direction verticale z et une des directions horizontales x, on utilise des schemas implicites a dierences nies, dits schemas compacts (voir Lele 4]). La troisieme direction (y) utilise des methodes pseudo-spectrales (Canuto et al 1]) pour les derivees premieres et secondes.

Les sch
emas spectraux

Si une variable donne de l'ecoulement (x y z t) est periodique, on peut le developper, en utilisant la transformee inverse de Fourier, 2 X;1 ^ (x y z t) = (x ky  z t)ek y ny

y

j =;

ny

2

ou ky est le nombre d'onde de ni par,

ky = L2 j y 34

(88)

p

i = ;1 est l'unite imaginaire, Ly est la longeur spatiale selon la direction y, et ny le nombre de points dans cette direction. Chaque coecient de Fourier ^(x y z t), peut ^etre obtenu en utilisant la transformee de Fourier (directe) discrete, nX ;1 (x y z t)e;k y ^(x ky  z t) = n1 y

y j =0

y

(89)

La derivation dans l'espace physique revient a une multiplication dans l'espace spectral:

d @ = k ^ y @y

(90)

Les sch
emas compacts

L'idee de depart des schemas compacts est d'exprimer une derivee en un point comme une fonction des derivees aux points voisins. Le methode est donc implicite, ce qui dans la pratique revient a resoudre un systeme d'equations pour calculer la derivee dans chaque point. Considerant une maillage uniforme avec coordonnees spatial xi = (i ; 1)'x, avec 'x = Const: (i 2 1 nx], xi 2 0 Lx]), et utilisant le developpement en serie de Taylor de la function fi = f (xi ) = f (x), avec premiere derivee fi0 = f 0(xi) = df (xi ) et derivee seconde fi00 = f 00(xi ) = dx d2f (x ), on peut arriver a la formule exacte, dx2 i p X

j =;p

j fi0+j =

q X

k=;q

ak fi+k + O('xn)

(91)

avec p q 2 N . L'ordre de l'approximation n, depend des restrictions imposees sur j and ak . Des cas particuliers sont:

   

8 j 6= 0  j = 0 ; > schema explicite. 9 j 6= 0  with j = 0 ; > schema Hermitian (Compact). 8 (j k)  ;j = j and a;k = ;ak ; > schema centre. 9 (j k) with ;j =6 j and/or a;k =6 ;ak ; > schema decentre.

Si on met p = 2 et q = 3 dans l'equation 91 on obtient,

fi0;2 + fi0;1 + fi0 + fi0+1 + fi0+2 = ; fi;1 + b fi+2 ; fi;2 + c fi+3 ; fi;3 + O('xn) a fi+12' x 4'x 6'x

Dans ce cas les restrictions dans (  a b c) sont:  a + b + c = 1 + 2 + 2 if n 2. 35

(92)

 a + 22b + 32c = 2 3!2! ( + 22 ) if n 4.  a + 24b + 34c = 2 5!4! ( + 24 ) if n 6.  ... Dans ce code toutes les derivees aux points interieurs sont calculees par un schema Compact centre, (93) fi0;1 + fi0 + fi0+1 = a fi+1';xfi;1 + b fi+2';xfi;2 : Pour les points interieurs entre i = 3 ::: nx ; 2, on utilise un schema de textrmme 6 ordre. Ceci est obtenu avec  = 1=3 et en utilisant les restrictions suivantes sur a et b,

a = ( + 2)=3 = 7=9  b = (4 ; 1)=12 = 1=36: (94) Pres des deux frontieres (i = 2 and i = nx;1 ) un schema classique de 4eme ordre Pade a ete utilisee. Ceci est obtenu avec  = 1=4 et, a = ( + 2)=3 = 3=4  b = (4 ; 1)=12 = 0: (95) Finallement, un schema implicite, decentre de 3eme odre est utilise pour les points aux frontieres (i = 1 and i = nx), (96) f10 + 2f20 = 2'1 x (;5f1 + 4f2 + f3) fn0 + 2fn0 ;1 = 2'1 x (5fn ; 4fn ;1 ; fn ;2 ) (97) Dans la pratique, le calcul d'une derivee dans un seul point x~p = (xp yp zp) consiste a resoudre une matrice avec tous les points dans la direction x. Le systeme s'ecrit, x

x

x

x

x

0 1 0 f0 1 1 1 BB 2 1 2 CC BB f120 CC BB CC BB f30 CC  1  BB CC BB .. CC ... ... ... BB CC BB . CC BB CC BB fi0 CC =  1  BB CC BB .. CC ... ... ... BB CC BB 0 . CC  1  BB CC BB fn ;2 CC n ;1 1 n ;1 A B @ @ fn0 ;1 CA n 1 fn0 x

x

x

x

36

x

x

0a 1 B ; a B 2 B B ; b B B B 1 B B 'x B B B B B B B @

b1 0 ;a ...

c1 a2 0 a ... ... ;b ;a ...

b ... ... 0 a ... ... ;b ;a

b ... ... 0 a b ;anx;1 0 anx;1 cnx bnx anx

1 0 f1 1 CC BB f2 CC CCC BBB f3 CCC CC BB ... CC CC BB C CC BB f.i CCC (98) CC BB .. CC CC BB fn ;2 CC CA B@ f CA n ;1 fn x x

x

avec les coecients,

1 = n = 2  a1 = ;an = ; 25  b1 = ;bn = 2  c1 = ;cn = 12

2 = n;1 = 14  a2 = an;1 = 34 1  = 13  a = 97  b = 36 Pour la deuxieme derivee, l'equation analogue a (99) avec fi0; > fi00 est, fi00;2 + fi00;1 + fi00 + fi00+1 + fi00+2 = 2fi + fi;1 + b fi+2 ; 2fi + fi;2 + c fi+3 ; 2fi + fi;3 + O('xn) (99) a fi+1 ; ' x 4'x 9'x Les restrictions en (  a b c) sont:

   

a + b + c = 1 + 2 + 2 if n 2. a + 22b + 32c = 2 2!4! ( + 22 ) if n 4. a + 24b + 34c = 2 6!4! ( + 24 ) if n 6. ...

Une fois de plus, tous les points interieurs sont calcules avec un schema Compact centre, ce qui pour la seconde derive s'est,

fi00;1 + fi00 + fi00+1 = a fi+1 ;'2fxi2+ fi;1 + b fi+2 ;'2fxi2+ fi;2

(100)

Pour les points interieurs (i = 3 ::: nx ; 2), on utilise un schema de 6eme ordre. Ceci correspond a  = 2=11 et, 37

a = 4(1 ; )=3 = 12=11  b = (4 ; 1)=12 = 3=44: (101) Comme auparavant, les points i = 2 et i = nx;1 sont traites avec un schema Pade de 4eme ordre. Pour la derivee seconde, ceci est obtenu avec  = 1=10 et, a = 4(1 ; )=3 = 12=10  b = (4 ; 1)=12 = 0: (102) Un schema decentre, implicite d'ordre 3 est utilise dans les points des frontieres i = 1 and i = nx, f100 + 11f200 = 2'1x2 (13f1 ; 27f2 + 15f3 ; f4) (103) fn00 + 11fn00 ;1 = 2'1x2 (13fn ; 27fn ;1 + 15fn ;2 ; fn ;3) (104) Le systeme nal d'equations est maintenant, x

x

x

x

x

x

0 1 0 f 00 1 1 1 BB f100 CC B C 2 1 2 B CCC BB f2300 CC B B  1  CC BB .. CC B . . . B . . . . . . CC BB . CC B B CC BB fi00 CC = B  1  B CC BB .. CC B ... ... ... B CC BB 00. CC B B  1  CC BB fn ;2 CC B B n ;1 1 n ;1 A B @ @ fn00 ;1 CA n 1 fn00 x

x

x

x

0 a1 b1 c1 B a2 ;2a2 a2 B B B b a ; 2( a + b) B B . . .. .. B 1 B B B b 'x2 B B B B B B B @

x

x

d1 a b ... ... ... a ;2(a + b) a b ... ... ... ... ... b a ;2(a + b) a anx;1 ;2anx ;1 dnx cnx bnx

avec les coecients,

10 1 CC B ff12 C CC CC BB B f CC B .3 CC CC BB .. CC CC BB fi CC CC BB .. CC CC BB . CC B C b C CC BB fn ;2 CC an ;1 A @ fn ;1 A fn an x x

x

x

(105)

1 = n = 11  a1 = an = 13  b1 = bn = ;27  c1 = cn = 15  d1 = dn = ;1 38

x

1  a = a = 12 2 = n;1 = 10 2 n;1 10 2  a = 12  b = 3  = 11 11 44 On peut ecrire les equations (98) et (105) sous la forme, A1f = '1x B1f (106) (107) A2f = '1x2 B2f ou A1, A2, B1 et B2 sont des matrices (nx
nx) et f , f et f sont des vecteurs (nx
1). La procedure pour obtenir les vecteurs f et f a partir de f , est la suivante: 1. Inversion des matrices A1 et A2 , 0

00

0

0

00

00

A1;1  A2;1 2. Multiplication des matrices B1 et B2 par le vecteur 1x f ,  1 1 'x B1f  'x2 B2f

(108) (109)

3. Multiplication des matrices A1;1 et A2;1 par les vecteurs 1x B1f et 1x2 B2 f ,

A1

;1

1  1 ; 1 B f  A2 'x2 B2f : 'x 1

(110)

Ceci reduit considerablement le temps de calcul de l'obtention des derivees.

3.4.3 Le syst
eme d'
equations nal discr
etis
e Les equations (82) a (87) s'ecrivent dans l'espace spectral de la facon suivante: Pour le champ de vitesse sans considerer le champ de pression:

ubp ; ubp;1 =  A p;1 p;1 ) +  A p;2 p;2 ) p c(~u  p c(~u  't Pour le champ du scalaire  ( , T ou S ): bp ; bp;1 =  B p;1 p;1 c(~up;2 p;2) p c (~u   ) + pB 't L'obtention de ubp, dans une etape intermediaire: 39

(111) (112)

d b 0p + ubp ubpx = 't dx x ubpy = iky 'tb 0p + ubpy ubpz = 't dzd b 0p + ubpz L'obtention du champ de vitesse corrige se fait a partir de: d b p + ubp ubpx = 't dx x p p uby = iky 'tb + ubpy d b p + ubp ubpz = 't dz z Les scalaires 0p et p sont obtenus a partir de: ! ! d2 ; k2 + d2 b 0p = 1 d ubp + ik ubp + d ubp y y dx2 y dz2 't dx x dz z ! ! d2 ; k2 + d2 b p = 1 d ubp + ik ubp + d ubp y y dx2 y dz2 't dx x dz z avec:

ou

( 2 ) 2 p 1 1 d d p p p p p 2 c Ax(~u ) = Ffuy!z ; uz !y g ; Ro ubx + Re dx2 ; ky + dz2 ubpx ( 2 ) 2 p 1 d d 1 p 2 p p p p c Ay (~u ) = Ffuz !x ; ux!z g + Ro uby + Re dx2 ; ky + dz2 ubpy ( 2 ) 2 p 1 d d p p p p 2 c Az (~u ) = Ffux!y ; uy!xg + Re dx2 ; ky + dz2 ubpz ; b ( 2 ) 2 p d d 1 2 Bc (~u ) = Ff~u  r~ g + RePr dx2 ; ky + dz2 b + b (

(113)

(114) (115) (116)

(117) (118)

! ) d + ik + d b (119) y dx dz La condition d'incompressibilite a veri er s'exprime dans l'espace de Fourier: d ubp + ik ubp + d ubp = 0 (120) dx x y y dz z

r~ = F ;1

40

3.4.4 Les conditions limites Les conditions aux limites pour la pression sont les conditions de Neumann, elles s'ecrivent:

! @P 't @x = ; (uxj1n ; uxj1n ) (121) 1nx ! @P = ; (uz j1n ; uzj1n ) (122) 't @x 1nz avec uxj1n = uz j1n = 0. On impose ces conditions aux frontieres lors de la correction des vitesses. Cette correction a lieu une fois que l'on a calcule la pression a partir de l'equation de Poisson, puis le gradient de pression. Les conditions aux limites pour les composantes de vitesses tangentielles aux parois non periodiques (x = 0 Lx et z = 0 ;H ), ainsi que le scalaire passif  (conditions 69 et 71) sont imposes a l'aide des schemas compacts d'ordre 2, soit: (123) f1 = f2 ; '2x f 02 Pour la condition a la frontiere horizontale de la surface (z = 0), on considere que la paroi est adiabatique (eq. (69)). x

x

x

z

z

z

3.5 Description des routines Routines principales (descrition): main2.f : C'est la routine de demarrage du code. Elle lit les donnees d'entree et les parametres de calcul (dans param)( ouvre les chiers d'ecriture de donnees( initialise les routines des ts (dans tfax)( calcule les nombres d'onde (dans wave)( et appelle la routine shear. shear2.f : C'est la routine que prepare les tableaux du code pour le calcul. Initialise les matrices du schema Compact pour les derivees premieres et secondes (dans derive1 et derive2)( de nie la fonction a (equation 14) pour le disque de refroidissement (dans idisq)( Initialise les champs de vitesses u~ , et les champs de uctuation de densite , temperature T , et salinite S , ou, s'il y a une reprise d'un calcul precedent, ces champs sont initialises par leurs valeurs precedentes( et ensuite apelle la routine rkutta.

C'est la routine qui demarre le boucle de Runge-Kutta pour l'avancement temporel. Dans cette boucle la routine calcule les termes convectifs des equations de transport de quantite de mouvement et de scalaire passif (densite, temperature ou salinite) (dans nonlear)( calcule les termes visqueux et de rkuttap2.f :

41

forcage des equations (dans right)( resoud l'equation de poisson pour la pression (dans pres). Routines auxiliaires (descrition): matrix2.f : Calcul des matrices A1 , B1 , A2 et B2 pour les derivees premieres et secondes des schemas Compact. param3.f : Lecture et adimensionalisation des parametres du calcul. wave2.f : Calcul des nombres d'onde pour les ts. sepx4.f : Resoulution de l'equation de poisson pour la pression. boundy2.f : Imposition des conditions limites pour les equations de transport de quantite de mouvement et la pression. Les parametres du calcul sont choisis dans le chier dat_in2.dat. La signi cation des plusieurs parametres est donnee dans le chier. Un exemple est donne en desous: 'les' --> 2 ------> 1 ------> 0 ------> 1 ------> 1.e-6 --> 1. -----> 0.6 ----> 1. -----> 1. -----> 0.0625 -> 1. -----> 50. ----> 48. ----> 0 ------> 1. -----> 5.e-2 --> 32. ----> 2. -----> 4.e-4 --> 800. ---> 16. ----> 1.e-6 --> 0. -----> 0.042 ---> 0. ------>

CURDAT: nom generique des fichiers resultat ITEMP: =0 no transp. eq. =1 density eq. =2 temperature and salinity eq. FTEMP: =1 body forcing in temp/salt eq. =2 imposed flux in temp/salt eq. IGAMAT: =0 const heat/salt flux =1 variable heat/salt flux. IGAMAUV: =0 no wind forcing at surface =1 wind forcing at surface ALPHA: ecart caracteristique a la densite moyenne PRANDTL: nombre de Prandtl PRAT: nombre de Prandtl turbulent CHI : =1 si calcul par LES - SF model (=0 sinon) AX: rapport d'aspect Lx/Ly AZ: rapport d'aspect Lz/Ly DT: pas de temps (en minutes) DELTAT: duree de la simulation (en heures-reelles) TCUT (en heures-reelles) NRESUVW: pas a partir duquel on repart (0 sinon) TSAVE: periode d'ecriture des champs (en heures) AMPBR: amplitude du bruit sur la temperature LY: largeur de boite (en kilometres) L echelle caracteristique (en kilometres) F: parametre de coriolis (en sec-1) FORC: forcage (en watt.metre-2) DDISQ: diametre de la zone de forcage (en kilometres) NU: viscosite (en metre2.sec-1) EMP: Evaporation minus precipatation (mm/day) TAUX: Shear stress caused by wind in the x direction (N.m-2) TAUY: Shear stress caused by wind in the y direction (N.m-2)

42

~ ~ ~

References 1] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang. Spectral Methods in Fluid Dynamics. Springer-Verlag, 1987. 2] J. Kim and P. Moin. Application of a fractional-step method to incompressible navier-stokes equations. J. Comp. Phys., 59:308{323, 1985. 3] H. Le and P. Moin. An improovement of fractional-step methods for the incompressible navier-stokes equations. J. Comp. Phys., 92:369{379, 1991. 4] S. K. Lele. Compact nite dierence schemes with spectral-like resolution. J. Comp. Phys., 103:15{42, 1992. 5] M. Lesieur and O. Metais. New trends in large-eddy simulations of turbulence. Annu. Rev. Fluid Mech., 28:45{98, 1999. 6] J. Padilla. Simulation des grandes echelles de convection turbulente en milieux tournant. PhD thesis, INPG, Grenoble., 1999. 7] J. Padilla-Barbosa and O. Metais. Large-eddy simulations of deep-ocean convection: analysis of the vorticity dynamics. Journal of Turbulence, 009, 2000. 8] J. H. Williamson. Low-storage runge-kutta schemes. J. Comp. Phys., 35:48{ 56, 1980.

43