70
ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique PEREZ Florian 1 Résumé Type de document : Mémoire de Travail de Fin d’Etudes d’Ingénieur de l’Ecole Nationale du Génie de l’Eau et de l’Environnement de Strasbourg / Master Recherche 2 ème année « Ingénierie et Technologie » de l’Institut Professionnel des Sciences et Technologies, Option Mécanique, Spécialité Fluide et Environnement. Résumé : L’objet de cette étude est l’utilisation des méthodes des Eléments Finis Mixtes Hybrides et des Volumes Finis pour la modélisation de systèmes hydrothermaux situés au niveau des dorsales océaniques. L’objectif est de valider un jeu de données créé sous le logiciel de modélisation CAST3M développé depuis plus de vingt ans par le Commissariat à l’Energie Atomique. Deux schémas numériques temporels ont été testés, le schéma numérique implicite en temps et schéma de type Crank-Nicholson. La comparaison des résultats de calculs obtenus pendant cette étude avec les résultats publiés dans la littérature pour des systèmes identiques a permis de mettre en évidence des critères de précision telles que des gammes de pas de temps de calculs et de pas d’espace utilisables. Ces critères sont indispensables à la compréhension de systèmes plus complexes tels que les systèmes hydrothermaux naturels, pour lesquels il n’y a pas de description précise du comportement. Mots clés : hydrogéologie, hydrothermalisme, convection de Rayleigh-Bénard, modélisation, éléments finis mixtes hybrides, schéma implicite, schéma de Crank-Nicholson, CAST3M.

Rapport_de_Stage test

Embed Size (px)

Citation preview

Page 1: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 1

Résumé

Type de document : Mémoire de Travail de Fin d’Etudes d’Ingénieur de l’Ecole Nationale du Génie de l’Eau et de l’Environnement de Strasbourg / Master Recherche 2ème année « Ingénierie et Technologie » de l’Institut Professionnel des Sciences et Technologies, Option Mécanique, Spécialité Fluide et Environnement. Résumé : L’objet de cette étude est l’utilisation des méthodes des Eléments Finis Mixtes Hybrides et des Volumes Finis pour la modélisation de systèmes hydrothermaux situés au niveau des dorsales océaniques. L’objectif est de valider un jeu de données créé sous le logiciel de modélisation CAST3M développé depuis plus de vingt ans par le Commissariat à l’Energie Atomique. Deux schémas numériques temporels ont été testés, le schéma numérique implicite en temps et schéma de type Crank-Nicholson. La comparaison des résultats de calculs obtenus pendant cette étude avec les résultats publiés dans la littérature pour des systèmes identiques a permis de mettre en évidence des critères de précision telles que des gammes de pas de temps de calculs et de pas d’espace utilisables. Ces critères sont indispensables à la compréhension de systèmes plus complexes tels que les systèmes hydrothermaux naturels, pour lesquels il n’y a pas de description précise du comportement. Mots clés : hydrogéologie, hydrothermalisme, convection de Rayleigh-Bénard, modélisation, éléments finis mixtes hybrides, schéma implicite, schéma de Crank-Nicholson, CAST3M.

Page 2: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 2

Abstract Type of Report : Final Dissertation thesis at ENGEES/Master thesis in the Professional Institute of Sciences and Technologies, Mechanic Option, Fluid and Environment specialty. Abstract : The matter of this study is the use of the Hybrid Mixed Finite Element method and Finite Volume method for modelling of hydrothermal systems situated on oceanic ridge like the Rainbow site. The target is to validate a numerical code created with program CAST3M which has been developed by the “Commissariat à l’Energie Atomique” (France) for twenty years. We have tested two numerical schemes applied to this code, the implicit scheme, and the Crank-Nicholson scheme. The comparison of the results of our computations with other results from different authors for identical systems, permitted us to give rise to accuracy criterions like ranges of usable computation time steps and space steps. Those criterions are essential for the understanding of more complex systems such as hydrothermal systems whose behaviors are not precisely known. Key words : hydrogeology, hydrothermalism, Rayleigh-Bénard convection, modeling, hybrid mixed finite elements method, implicit scheme, Crank-Nicholson scheme, CAST3M.

Page 3: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 3

Remerciements

Mes remerciements vont tout particulièrement à Claude Mügler, chercheur au CEA-LSCE

qui m’a encadré lors de ce stage, à Philippe Jean-Baptiste, co-encadrant de ce stage, à Robert Mosé, directeur adjoint de l’IMFS et professeur à l’ENGEES, tuteur de mon stage. A toute l’équipe « Modélisation Hydrologique », un grand merci pour l’aide qu’ils m’ont apportée pendant ce stage. A mes proches qui m’ont soutenu. Aux équipes de l’ENGEES et de l’IPST pour la préparation de ce stage et la formation qu’ils m’ont donnée.

Page 4: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 4

Listes des notations Liste des abréviations µ : Viscosité dynamique du fluide [M][L-1][T -1] 2D : Bidimensionnel A : Rapport d’aspect, rapport largeur/hauteur des rouleaux de convection ou du

domaine étudié EFMH : Eléments Finis Mixtes Hybrides Fo : Nombre de Fourier g : Accélération gravitaire [L].[T-2]

h : Charge hydraulique [L], P

h zgρρρρ

= + dans le cas d’un écoulement souterrain

k : Perméabilité intrinsèque [L²] K : Perméabilité [L].[s-1] MAR : Mid-Atlantic Ridge, Dorsale Médio-Atlantique Nu : Nombre de Nusselt P : Pression [M].[L-1].[T -²] Pe : Nombre de Péclet Ra : Nombre de Rayleigh z : Altitude par rapport au niveau de référence [L] ρρρρ : Masse volumique du fluide [M].[L-3] VF : Volumes finis Liste des abréviations d’institutions ANR : Agence Nationale pour la Recherche CEA : Commissariat à l’Energie Atomique CNRS : Centre National de la Recherche Scientifique IFREMER : Institut Français de Recherche pour l’Exploitation de la MER IMFS : Institut de Mécanique des Fluides et des Solides IPST : Institut Professionnel des Sciences et Technologies LSCE : Laboratoire des Sciences du Climat et de l’Environnement UVSQ : Université Versailles St-Quentin

Page 5: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 5

Liste des tableaux

Tableau 1 : Tableau chronologique du déroulement du travail de fin d’études........................... 11 Tableau 2 : Concentrations des espèces chimiques mesurées dans les fluides prélevés lors de la campagne FLORES (Charlou et al. 2002)..................................................................................... 16 Tableau 3 : Compositions moyennes des gaz mesurées en sortie des fumeurs (Charlou et al. 2002)............................................................................................................................................... 16 Tableau 4 : Comparaison des gammes de pas de temps utilisables pour une fréquence d’oscillation de 500........................................................................................................................ 34 Tableau 5 : Paramètres du système............................................................................................... 36 Tableau 6 : Tableau récapitulatif des résultats de calculs obtenus avec CAST3M dans le cas d’un domaine fermé................................................................................................................................ 39 Tableau 7 : Comparaisons des résultats obtenus avec ceux publiés dans la littérature............... 42 Tableau 8 : Comparaisons des résultats obtenus sous CAST3M à l’aide de la méthode des EFMH et de la méthode des VF................................................................................................................. 45 Tableau 9 : Comparaison des résultats obtenus en boîte ouverte pour un nombre de Rayleigh croissant progressivement de 90, 108, 200, 300............................................................................ 47 Tableau 10 : Résultats obtenus en boîte ouverte pour des nombres de Rayleigh de 90, 108, 200, 300 et 400 avec une condition initiale de type CIO1.................................................................... 48 Tableau 11 : Comparaison des résultats obtenus en boîte ouverte pour un nombre de Rayleigh croissant progressivement de 400, 425, 430, 450 et 470............................................................... 49 Tableau 12 : Résultats obtenus en boîte ouverte pour des nombres de Rayleigh de 400, 425, 430, 450 et 470 avec une condition initiale de type CI1....................................................................... 49

Liste des figures et illustrations

Figure 1 : Mouvements de convection au sein du manteau........................................................... 12 Figure 2 : Modélisation des flux de chaleur émis au niveau de la croûte océanique. Les fortes valeurs de flux de chaleur trahissent la présence des dorsales océaniques (Planète terre, ENS Lyon (Stéphane Labrosse))............................................................................................................ 13 Figure 3 : Site de Rainbow (Douville et al., (2002))..................................................................... 14 Figure 4 : Principe général de la convection hydrothermale sur une coupe transversale d'une dorsale rapide (Figure d’Yves Fouquet, IFREMER)..................................................................... 15 Figure 5 : Evolution de la concentration en Chlorure en fonction de la concentration en Magnésium (Charlou et al., 2002)................................................................................................. 16 Figure 6 : Schéma de principe, modèle de convection avec séparation de phase inspiré du modèle de Kawada et al. (2004)................................................................................................................. 17 Figure 7 : Représentation schématique de la convection de Rayleigh-Bénard............................. 21 Figure 8 : Algorithme de résolution utilisé lors de cette étude..................................................... 35 Figure 9 : Schéma du domaine considéré et conditions aux limites............................................. 36 Figure 10 : Schéma du domaine considéré et conditions aux limites adimensionnées................. 36 Figure 11 : Champ de températures obtenu.................................................................................. 37 Figure 12 : Champ de températures publiés dans (Fontaine et al., 2007)................................... 37 Figure 13 : Champ de nombres de Péclet pour un maillage de 128×128, Fo = 0,5.................... 38 Figure 14 : Champ de températures permanent obtenu pour Ra = 950 et une condition initiale de type CI1.......................................................................................................................................... 40 Figure 15 : Evolution du nombre de Nusselt en fonction du temps, pour Ra = 950 et une condition initiale de type CI1......................................................................................................... 40

Page 6: Rapport_de_Stage test

Figure 16 : Champ de températures obtenu pour Ra = 950 et une condition initiale de type CIF800........................................................................................................................................... 41 Figure 17 : Evolution du nombre de Nusselt en fonction du temps, pour Ra = 950 et une condition initiale de type CIF800................................................................................................... 41 Figure 18 : Champ de températures obtenu pour Ra = 1200 et une condition initiale de type CIF950........................................................................................................................................... 41 Figure 19 : Evolution temporelle du nombre de Nusselt, pour Ra = 1200 et une condition initiale de type CIF950................................................................................................................... 41 Figure 20 : Champ de températures obtenu pour Ra = 800 et une condition initiale de type CI1......................................................................................................................................................... 43 Figure 21 : Evolution du nombre de Nusselt en fonction du temps, pour Ra = 800 et une condition initiale de type CI1......................................................................................................... 43 Figure 22 : Champ de températures permanent obtenu pour Ra = 950 et une condition initiale de type CI1.......................................................................................................................................... 44 Figure 23 : Evolution du nombre de Nusselt en fonction du temps, pour Ra = 950 et une condition initiale de type CI1......................................................................................................... 44 Figure 24 : Champ de températures permanent obtenu pour Ra = 950 et une condition initiale de type CIF800............................................................................................................................... 44 Figure 25 : Evolution du nombre de Nusselt en fonction du temps, pour Ra = 950 et une condition initiale de type CIF800................................................................................................... 44 Figure 26 : Schéma du domaine fermé et conditions aux limites adimensionnées....................... 45 Figure 27 : Evolution temporelle du nombre de Nusselt, obtenu à l’aide de la méthode des....... 46 Figure 28 : Champ de températures obtenu pour Ra = 470 et une condition initiale de type CIO425........................................................................................................................................... 47 Figure 29 : Evolution temporelle du nombre de Nusselt, obtenu à l’aide de la méthode des VF avec CAST3M, en faisant croître progressivement le nombre de Rayleigh.................................. 48

Page 7: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

7

SOMMAIRE

INTRODUCTION ............................................................................................................................................... 9

1. CONTEXTE GEOPHYSIQUE ET DESCRIPTION DES ECOULEMENTS HYDROTHERMAUX 12 1.1 Contexte géophysique ................................................................................................................ 12 1.2 Les écoulements hydrothermaux de la dorsale lente medio-atlantique, exemple du site de Rainbow ........................................................................................................................................... 13

1.2.1 Découverte du site de Rainbow....................................................................................... 14 1.2.2 Flux thermiques : Méthodes d’évaluation et résultats..................................................... 15 1.2.3 Géochimie : Résultats des prélèvements effectués sur le site de Rainbow..................... 16

Profondeur .......................................................................................................................16 1.3 Premières interprétations de ces résultats................................................................................... 17

1.3.1 Discussion sur une éventuelle séparation de phase......................................................... 17 1.3.2 Serpentinisation, production d’hydrogène et évolution de la perméabilité..................... 18

1.4 Morphologie du site ................................................................................................................... 19 1.5 Conclusion : Les inconnues du système, intérêt de la modélisation .......................................... 19

2. MISE EN PLACE DU MODELE................................................................................................................ 20 2.1 Hypothèses simplificatrices appliquées au modèle.................................................................... 20 2.2 Analogie des systèmes hydrothermaux avec la convection de Rayleigh-Bénard ...................... 20

2.2.1 La convection de Rayleigh-Bénard ................................................................................. 20 2.2.2 Le nombre de Rayleigh : explication physique de la convection.................................... 21 2.2.3 Convection de Rayleigh-Bénard en milieu poreux : description des résultats

courants 22 2.2.4 Pourquoi le nombre de Nusselt oscille-t-il ? ................................................................... 23 2.2.5 Conclusion....................................................................................................................... 24

2.3 Modèle mathématique général ................................................................................................... 24 2.3.1 Modèle hydraulique......................................................................................................... 24 2.3.2 Modèle de transfert de la chaleur .................................................................................... 25 2.3.3 Equations d’états ............................................................................................................. 26

2.4 Paramètres du modèle ................................................................................................................ 27 2.4.1 Paramètres de l’écoulement ............................................................................................ 27 2.4.2 Paramètres du transfert thermique................................................................................... 28 2.4.3 Hauteur caractéristique du domaine................................................................................ 28

Mémoire de Travail de Fin d’Etudes

réalisé au LSCE (CEA-CNRS-UVSQ Saclay, 91)

Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

Page 8: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

8

2.5 Simplification du modèle mathématique ................................................................................... 28 2.5.1 Formulation spécifique et simplification du modèle hydraulique................................... 28 2.5.2 Evolution de la masse volumique.................................................................................... 29 2.5.3 Formulation pour le calcul du transport de la chaleur..................................................... 29 2.5.4 Adimensionnalisation des équations ............................................................................... 29 2.5.5 Non-linéarité du système................................................................................................. 30

3. OUTILS NUMERIQUES ............................................................................................................................. 31 3.1 Outil de calcul en langage Gibiane : CAST3M.......................................................................... 31 3.2 Modélisation à l’aide des EFMH d’un champ hydrothermal dorsalien ..................................... 32

3.2.1 Choix du schéma spatial de discrétisation....................................................................... 32 3.2.2 Choix du schéma temporel de discrétisation................................................................... 32 3.2.3 Précision de calcul et stabilité : les critères..................................................................... 32 3.2.4 Discussion a priori sur la stabilité du schéma et la précision du calcul : pas

d’espace et pas en temps ............................................................................................................... 33 3.3 Modélisation à l’aide des VF d’un champ hydrothermal dorsalien ........................................... 34

3.3.1 Choix du schéma spatial de discrétisation....................................................................... 34 3.3.2 Choix du schéma temporel de discrétisation................................................................... 34 3.3.3 Les raisons d’une telle discrétisation .............................................................................. 34

3.4 Algorithme de résolution............................................................................................................ 35

4. VALIDATION DU JEU DE DONNEES – CAS-TEST ............................................................................. 36 4.1 Validation des modèles numériques en domaine fermé............................................................. 36

4.1.1 Description des conditions initiales et des conditions aux limites.................................. 36 4.1.2 Résultats obtenus à l’aide de la méthode des EFMH...................................................... 37 4.1.3 Résultats obtenus à l’aide de la méthode des VF............................................................ 43

4.2 Validation du jeu de données en domaine ouvert ...................................................................... 45 4.2.1 Conditions initiales et conditions aux limites ................................................................. 45 4.2.2 Résultats des calculs........................................................................................................ 46

4.3 Discussion .................................................................................................................................. 50 4.3.1 Approximation de Boussinesq, linéarité de l’équation de variation de la densité et

viscosité constante......................................................................................................................... 50 4.3.2 Linéarisation et précision du calcul................................................................................. 50 4.3.3 Importance de la condition initiale.................................................................................. 51

4.4 Conclusions et perspectives de modélisation d’un système hydrothermal réel ......................... 52 4.4.1 Conclusion sur l’utilisation des EFMH et des VF .......................................................... 52 4.4.2 Perspectives de modélisation : le champ hydrothermal réel ........................................... 53

BILAN DU STAGE ET CONCLUSION......................................................................................................... 54

BIBLIOGRAPHIE ............................................................................................................................................56

ANNEXES.......................................................................................................................................................... 60

Page 9: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

9

INTRODUCTION Contexte et objectifs de l’étude

A l’heure actuelle, un important enjeu scientifique et économique est de pouvoir faire face à la croissance de la demande énergétique mondiale tout en maîtrisant ses conséquences sur l’environnement. Les réserves de pétrole n’étant pas infinies, nombreuses sont les préoccupations concernant la recherche de nouvelles sources d’énergie. La combustion de l’hydrogène paraît sans conséquence immédiate sur l’environnement et la santé des populations puisqu’elle ne libère que de l’eau. La difficulté consiste en la production de cet hydrogène qui, elle, reste très coûteuse en énergie.

Cependant, des sources naturelles d’hydrogène existent, notamment au niveau des

dorsales océaniques, lieux de naissance de la croûte océanique. Ce stage de fin d’études s’intéresse plus particulièrement à la dorsale lente médio-atlantique. En effet, à son niveau, ont été découverts de nombreux champs hydrothermaux dont les cheminées, appelées fumeurs, laissent échapper des fluides chauds aux propriétés chimiques intéressantes puisqu’ils contiennent de nombreux métaux et pour ce qui nous concerne ici, dans certains cas, de l’hydrogène en forte concentration. C’est le cas pour un site particulièrement étudié par l’IFREMER Brest : le site de Rainbow.

Le Laboratoire des Sciences du Climat et de l’Environnement (LSCE), en la personne

de Claude Mügler, chercheuse du groupe modélisation hydrologique dirigé par Emmanuel Mouche, a pour but de répondre aux nombreuses interrogations en lien avec la nature même de l’écoulement hydrothermal de l’eau de mer au sein des roches constitutives de la dorsale médio-atlantique. En effet, la méconnaissance des propriétés des roches et des phases fluides qui les traversent, des réactions induisant les caractéristiques physico-chimiques des fluides hydrothermaux et le manque d’information sur la géométrie et les conditions aux limites en température imposent de faire des hypothèses fortes sur le fonctionnement du champ hydrothermal. La modélisation numérique, appuyée par des données de terrain et des résultats expérimentaux, devient alors indispensable afin d’infirmer ou confirmer certaines hypothèses.

Ce stage de fin d’études s’inscrit dans le cadre d’un projet scientifique soumis par le

CEA, l’ IFREMER et d’autres collaborateurs à l’ Agence Nationale pour la Recherche (ANR) : le projet MANTHY (MANTellic HYdrogen). Au sein de ce projet pluridisciplinaire, il m’a été demandé d’amorcer un modèle numérique permettant de décrire le fonctionnement d’un champ hydrothermal. Amorcer un modèle numérique de ce type demande de comprendre, dans un premier temps, la nature de l’écoulement et de traduire cette nature par des équations mathématiques décrivant avec le plus d’exactitude possible cet écoulement. La première partie de ce rapport est donc consacrée à la description d’un système hydrothermal à l’aide de l’exemple du site de Rainbow. Le modèle mathématique choisi est décrit en seconde partie. Par la suite, la troisième partie présente les méthodes numériques utilisées au cours de cette étude afin de résoudre les équations du modèle. Le type de schéma numérique choisi a une conséquence directe sur les résultats de calculs notamment sur leurs précisions. Il faut donc a priori faire un bon choix de schéma mais aussi, a posteriori, valider le modèle numérique en comparant nos résultats de calculs avec les résultats d’autres modèles utilisés pour des cas académiques, ce qui est fait en quatrième partie. Dans ce rapport, la méthode des éléments finis mixtes hybrides et la méthode des volumes finis sont testées.

Page 10: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

10

Présentation du LSCE

Le LSCE, Laboratoire des Sciences du Climat et de l’Environnement, est une des 54 unités mixtes de recherche du CEA, ce dernier étant, via cette unité, en partenariat avec l’Université de Versailles Saint-Quentin (UVSQ) et le Centre National de la Recherche Scientifique (CNRS). Le LSCE (UMR 1572) fait aussi partie des six laboratoires associés formant l’Institut Pierre Simon Laplace (IPSL) spécialisé

dans les sciences de l’environnement. Lors de mon stage, j’ai pu assister à l’évaluation du LSCE par un comité de l’Agence d’Evaluation de la Recherche et de l’Enseignement Scientifique (AERES) qui, dans son rapport, a mentionné sa satisfaction et a rappelé quelques points forts du LSCE, comme par exemple :

- « une position exceptionnellement forte sur la scène internationale » ; - « une maîtrise d’instrumentation et de méthodologie de classe mondiale pour l’analyse

des traceurs ». Les points à améliorer stipulés touchaient en particulier l’amélioration de la cohésion

entre les différentes équipes du LSCE, le développement des coopérations et l’intégration des laboratoires des unités mixtes tels que ceux de l’UVSQ.

La recherche au LSCE gravite autour de trois thèmes :

- Evolution du climat. - Cycles biogéochimiques. - Géosciences.

C’est au sein du thème Géosciences, et plus particulièrement de l’équipe « Modélisation

Hydrologique » dirigée par Emmanuel Mouche, que j’ai effectué mon stage. L’objectif de recherche de cette équipe consiste en la modélisation des hydrosystèmes continentaux appliquée aux différents thèmes que sont la gestion de la ressource en eau et le transfert de polluants, l’impact de la variabilité climatique sur l’hydrologie d’une région, la séquestration du CO2 dans des couches géologiques profondes etc…

Le CEA est un organisme public français de recherche scientifique créé en Octobre 1945. Il opère dans plusieurs domaines : l’énergie, les technologies pour l’information et la santé, la défense et la recherche fondamentale. Plus de 15 600 personnes collaborent au sein des neuf centres répartis dans toute la France, 5000 personnes pour le seul centre principal situé à Saclay, dans l’Essonne, sur lequel est

implanté l’un des deux centres du LSCE, le second centre du LSCE étant implanté non loin de là, à Gif-Sur-Yvette. En 2007, le CEA avait déposé 447 brevets prioritaires et était propriétaire de 1608 brevets prioritaires délivrés et en vigueur. Son budget annuel avoisinait les 3,4 milliards d’euros.

Le CNRS est un établissement public à caractère scientifique et technique (EPST) placé sous la tutelle du Ministère de l’enseignement supérieur et de la recherche. 11600 chercheurs, 14400 ingénieurs, techniciens et administratifs se répartissent parmi plus de 1200 unités de recherche et de service. Son budget pour l’année 2008 s’élève à 3,277 milliards d’euros dont 588 millions d’euros de ressources propres. Le LSCE fait partie du pôle de recherches « Environnement et Développement

Page 11: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

11

Durable » du CNRS.

L’ UVSQ, née en 1991 de la réunion de plusieurs annexes d’universités parisiennes, regroupe trente laboratoires et enseigne dans une multitude de disciplines : Sciences exactes, Sciences sociales, Sciences humaines, Sciences juridiques et politiques, Ingénierie et Technologie, Médecine.

Présentation de l’IFREMER

L’ IFREMER est l’un des collaborateurs du LSCE dans le projet MANTHY . Avec un budget annuel de près de 235 millions d’euros,

l’IFREMER possède 1500 salariés et collabore avec l’armateur Genavir pour sa flotte (8 navires, 1 submersible habité, 1 engin téléopéré et 2 AUV ou drones sous-marins). Cinq centres existent, Manche/Mer du Nord, Nantes, Méditerranée, Tahiti et l’IFREMER Brest.

Déroulement du stage

Tableau 1 : Tableau chronologique du déroulement du travail de fin d’études.

Dates Tâches Détails

Compréhension de la problématique du 26/01/09 au 6/02/09

Recherche Bibliographique

Initiation à CAST3M Décryptage d’un jeu de données CAST3M

-Problème d’ELDER – rapport explicatif

du 9/02/09 au 23/03/09

Recherche Bibliographique

Encodage sous CAST3M

Encodage pour un problème hydrothermal simple -Cas

d’une BOITE FERMEE – Schéma implicte en temps et Schéma de Crank-Nicholson.

Mardi 17 mars : Présentation de l’avancement du travail de stage à Philippe Jean-Baptiste (Equipe Gaz rares CEA-LSCE).

Vendredi 20 mars : Présentation de l’avancement du travail de stage à l’équipe « Modélisation hydrologique » (CEA-LSCE).

du 23/03/09 au 29/05/09

Recherche Bibliographique

Validation du Jeu de Données

Encodage sous CAST3M

Recherche des conditions de validité du Jeu de données. Encodage du cas d’une BOITE OUVERTE.

Jeudi 16 avril : Présentation de l’avancement du travail de stage à Jean-Luc CHARLOU (IFREMER Brest), Philippe

JEAN-BAPTISTE et Claude MUGLER (CEA-LSCE).

du 15/04/09 au 31/07/09

Recherche Bibliographique

Validation du Jeu de Données Ecriture d’un jeu de données (VF) pour les cas ouvert et fermé Validations des jeux de données.

Page 12: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

12

1. CONTEXTE GEOPHYSIQUE ET DESCRIPTION DES ECOULEMENTS HYDROTHERMAUX

1.1 Contexte géophysique

La Terre est une planète tellurique tout comme Mercure, Mars ou Vénus. Depuis sa

formation, il y a quelques 4,5 milliards d’années, la Terre a lentement acquis une organisation en succession d’enveloppes concentriques qui ont chacune des propriétés différentes.

Le manteau terrestre est le lieu de mouvements gigantesques de convection, phénomène permettant d’évacuer le plus efficacement possible la chaleur interne de la terre par transport de matière chaude. On estime à 42×1012 W la puissance dissipée par le noyau terrestre sous forme de chaleur, les trois quarts de cette énergie provenant de la décomposition isotopique des éléments du noyau.

Dans les profondeurs du manteau, les péridotites, roches constitutives de ce dernier, s’échauffent ; de gigantesques cellules de matière chaude, plus légère, se mettent alors en mouvement au travers des masses plus froides des zones supérieures.

Comme le montre la Figure 1, à l’opposé de ce mouvement ascendant, dont l’existence des dorsales océaniques est le témoin majeur, des cellules froides plongent vers les entrailles de la terre : au niveau des zones de subduction, la lithosphère océanique s’enfonce sous la lithosphère continentale plus légère et se délite dans l’asthénosphère. Parfois, les cellules descendantes atteignent la limite séparant le manteau inférieur du noyau externe. Il y a donc un échange de matière entre les différentes enveloppes externes et internes de la Terre.

A proximité de la surface, au niveau des dorsales, les cellules ascendantes rencontrent des conditions de pressions et de températures induisant la formation de magma qui se concentre et forme une chambre magmatique.

Figure 1 : Mouvements de convection au sein du manteau. Les flèches rouges représentent les roches chaudes entrant dans un mouvement ascendant, les flèches bleues les roches froides entrant dans un mouvement descendant.

Manteau inférieur

Lithosphère continentale

Asthénosphère

Lithosphère océanique

Noyau externe

Dorsale océanique

Cellule de convection

Cellule de convection

Chambre magmatique

Possible transfert de masse du manteau supérieur au manteau inférieur

Page 13: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

13

Au droit de ces dorsales, la croûte océanique, fine et peu dense, s’écarte et se fracture, tractée par son autre extrémité, plus froide et plus dense, subductée. Des failles normales divisent la croûte en d’énormes blocs rocheux. Le bombement observé provient de la poussée des matières chaudes (cf. Figures 1 et 2) qui remontent vers la surface terrestre pour former, en se refroidissant, une nouvelle croûte. Au centre de ce bombement, l’effondrement des blocs donne naissance à une vallée étendue appelée graben. Figure 2 : Modélisation des flux de chaleur émis au niveau de la croûte océanique. Les fortes valeurs de flux de chaleur trahissent la présence des dorsales océaniques (Planète terre, ENS Lyon (Stéphane Labrosse)).

Dans cet environnement complexe - une croûte non homogène, fracturée, à fort gradient thermique - de l’eau de mer s’infiltre, s’écoule, échange avec la roche des éléments chimiques : les sels précipitent ou se dissolvent, les minéraux s’appauvrissent ou s’enrichissent, s’hydratent puis réagissent et forment de nouvelles espèces minérales modifiant les caractéristiques des roches. L’eau de mer qui a pénétré la croûte en ressort, elle aussi, transformée physiquement et chimiquement et, à mesure que l’on s’approche du centre du bombement, donc des masses rocheuses chaudes, les écoulements s’intensifient. Ce travail de fin d’études porte plus particulièrement sur les écoulements d’eau de mer au travers de la croûte océanique formée au niveau de la dorsale lente medio-atlantique.

1.2 Les écoulements hydrothermaux de la dorsale lente medio-atlantique, exemple du site de Rainbow

En 1967, Palmason (Palmason 1967) avait suggéré que les variations locales des flux de chaleur mesurés au niveau du plancher océanique pouvaient avoir pour origine une circulation hydrothermale au sein de la croûte océanique. Depuis lors, de nombreuses expéditions ont confirmé cette hypothèse.

Page 14: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

14

1.2.1 Découverte du site de Rainbow

De nombreux sites hydrothermaux situés sur cette dorsale sont maintenant connus entre le 14°N et les îles des Açores. L’expression la plus visible et certainement la plus spectaculaire de ces écoulements, se traduit par l’existence de cheminées, appelées « fumeurs » qui émettent des fluides à hautes températures, de 70 à 370°C. Au niveau de ces cheminées, se sont développés des écosystèmes particuliers dont les organismes sont adaptés aux conditions de pression et de température régnant à de telles profondeurs, mais aussi aux caractéristiques chimiques spécifiques des fluides émis par ces cheminées.

Le site hydrothermal de Rainbow est particulièrement étudié par l’IFREMER de Brest et

illustre parfaitement la notion de champ hydrothermal. Il a pour la première fois été identifié lors de la « HEAT cruise – Charles Darwin 89 » réalisée du 19 Août au 13 Septembre 1994 (German et al., 1995). Long de 100 m et atteignant 200 m de large, il est localisé à 36°13,80 N et 33°54.12 W (cf. Figure 3), sur le flanc ouest de la dorsale à des profondeurs comprises entre 2270 et 2330 m (Charlou et al., 2002).

Figure 3 : Site de Rainbow, (a) et (b), localisation du site au droit de la dorsale médio-atlantique ; (c) principaux fumeurs du site témoins de la présence du champ hydrothermal (Douville et al., (2002)).

Dix fumeurs émettent très activement des fluides dont la composition chimique est

identique d’un fumeur à l’autre (Charlou et al., 2002). Comme le montre le tableau 2 (cf §1.2.3), cette composition est très différente de celle de l’eau de mer. La composition physico-chimique

Page 15: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

15

de ces fluides émis a permis de dégager un scénario général expliquant la formation de ces fluides : comme nous l’avons dit précédemment, au niveau de la dorsale, l’eau de mer s’infiltre dans des fissures et crevasses de la croûte et s’enfonce de plus en plus profondément ; puis cette eau se réchauffe en s’approchant des chambres magmatiques présentes au sein de la lithosphère (Figure 4). Les fluides chauds remontent en surface physiquement et chimiquement modifiés par leur passage au travers des roches ultrabasiques que sont les péridotites (Douville et al., 2002 ; Charlou et al., 2002).

Dans le cas du site de Rainbow, cette modification chimique entraîne une très forte

concentration en hydrogène, source potentielle d’énergie des années à venir.

Figure 4 : Principe général de la convection hydrothermale sur une coupe transversale d'une dorsale rapide (Figure d’Yves Fouquet, IFREMER).

L’intérêt d’un tel site semble donc évident et c’est dans ce cadre que ce travail de fin

d’études s’est inscrit.

1.2.2 Flux thermiques : Méthodes d’évaluation et résultats De nombreuses méthodes d’estimation des flux hydriques existent : les méthodes visuelles d’estimation des vitesses des particules sortant des fumeurs, des déterminations du champ de vitesses à l’aide de moulinets ou encore, la détermination des débits à l’aide d’un isotope de l’hélium, le 3He. Cette dernière méthode (Jean-Baptiste et al. 2004, 2008) permet dans certains cas non seulement de déterminer les flux hydriques focalisés, c'est-à-dire sortant au niveau des fumeurs, mais aussi le flux hydrique diffus s’échappant à l’interface croûte océanique-océan. Les flux thermiques en sont déduits par mesure de la température.

Page 16: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

16

1.2.3 Géochimie : Résultats des prélèvements effectués sur le site de Rainbow

Compositions ioniques

On constate que les cheminées sont plus riches en Cuivre, Zinc, Cobalt et Nickel qu’il

n’est habituellement observé notamment sur le site de Logatchev (14°45’N, Charlou et al., 1998 et Douville et al., 2002) et dans pratiquement chacun des sites à environnement basaltique tels que MARK-1/2, TAG, Broken spur, Lucky strike et Menez Gwen.

Tableau 2 : Concentrations des espèces chimiques mesurées dans les fluides prélevés lors de la campagne FLORES (Charlou et al. 2002) (R : Rainbow, SW : Eau de mer ambiante).

Figure 5 : Evolution de la concentration en Chlorure en fonction de la concentration en Magnésium (Charlou et al., 2002). Les concentrations sont en mM/kg.

Représente les valeurs mesurées pour l’eau de mer ambiante.

La figure 5 montre que les fluides hydrothermaux issus des fumeurs du site du Rainbow

sont enrichis en chlorure (cf. § 1.3.1). D’après Bischoff et Rosenbauer (1987), la présence d’une forte concentration en chlorure favorise la richesse en métaux (Cuivre, Fer, Nickel, Zinc…) via des réactions de complexation.

Composition des gaz

La composition en gaz des fluides hydrothermaux éjectés par les fumeurs du site de

Rainbow est caractérisée par une concentration relativement faible en H2S (1.2 mmol/kg) et une forte concentration en CH4 (2,5 mmol/kg), en H2 (16 mmol/kg), et CO (5 µmol/kg). Les concentrations des gaz sont ici fonctions linéaires de la concentration en magnésium. La similitude des concentrations d’un fumeur à un autre du site du Rainbow laisse penser que les fluides éjectés par ceux-ci possèdent tous la même origine.

Profondeur T°C H2S

(mM) CO2

(mM) CH4

(mM) CO

(mM) Ar

(mM) N2

(mM) H2

(mM) δδδδ13C(CO2) δδδδ13C(CH4)

R 2300 m 365 1.2 16 2.5 5000 - 1.8 16 -3.15 -15.8 SW 2°C 0 2,30 0,00003 0,3 16 0,59 0,0004 -5,1/-5,9 -

Tableau 3 : Compositions moyennes des gaz mesurées en sortie des fumeurs (Charlou et al. 2002).

Profondeur T°C pH Si(OH)4

(mM) Cl

(mM) Br

(mM) SO4

(mM) Na

(mM) Li

(mM) K

(mM) Rb

(mM) Cs

(mM) Ca

(mM) Sr

(µM) Ba

(µM) Fe

(µM) Mn

(µM) Cu

(µM) Zn

(µM)

R 2300 m 365 2,8 6,9 750 1178 0 553 340 20,4 36,9 333 66,6 200 >67 24050 2250 121/ 162

115/185

SW

2°C 7,8 <0,2 546 838 28,2 464 26 9,8 1,3 2,3 10,2 87 0,14 <0,001 <0,001 0,007 0,01

Page 17: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

17

Couche stagnante - saumure

Zone d’entraînement

Point chaud - intrusion magmatique

Couche de convection

Zone de rechargement, l’eau de mer à 2°C pénètre la croûte

Milieu fracturé

1.3 Premières interprétations de ces résultats

Ces données de terrain ont permis d’élaborer de nombreuses hypothèses sur le fonctionnement des sites hydrothermaux.

1.3.1 Discussion sur une éventuelle séparation de phase

La forte concentration en chlorure dans les fluides hydrothermaux laisse penser qu’il y a,

à proximité de la source de chaleur, une phase saumure, réservoir à chlorure qui alimente l’eau circulant au sein des fissures et remontant à la surface (cf. figure 6).

Figure 6 : Schéma de principe, modèle de convection avec séparation de phase inspiré du modèle de Kawada et al. (2004).

En effet, d’après Sourirajan et Kennedy (1962), l’enrichissement des fluides

hydrothermaux en chlorure serait dû à des processus incluant une séparation des phases vapeur et saumure à proximité de la source de chaleur. D’autres processus, l’hydratation des roches ignées (volcaniques) ou la précipitation/ dissolution des minéraux chlorés (Seyfried et al., 1986 ; Von Damm, 1988), peuvent être mis en jeu. Par contre, il n’existe pas de réaction entre la roche et l’eau de mer connue pour être prédominante dans le phénomène d’enrichissement de l’eau en chlorures (Bischoff et Rosenbauer, 1989). D’après Charlou et al., (2002), les concentrations en chlorures des fluides hydrothermaux de Rainbow, stables depuis 1997 (Charlou 2009), laissent supposer l’existence d’une séparation de phase, l’eau se trouvant à des conditions proches de son point critique à proximité de la source de chaleur. De plus, des inclusions de fluides à forte salinité ont été trouvées dans des gabbros originaires de la dorsale médio-atlantique dont les concentrations en NaCl sont les mêmes que celles issues d’une séparation de phase au point critique (Kelley et Gillis, 1993 ; Kelley, 1997).

Néanmoins, la stabilité de la concentration en chlorures des fluides hydrothermaux de

Rainbow, semble se rapprocher du modèle de Schoofs et Hansen (2000) qui présente une phase saumure lentement entraînée par la convection des fluides d’une couche supérieure (Charlou et al., 2002).

Il est difficile d’appliquer l’un ou l’autre des modèles ne serait-ce que parce que nous ne

connaissons ni l’impact de l’activité sismique ni celui de la précipitation ou de la serpentinisation sur l’évolution exacte de la perméabilité effective du manteau.

Page 18: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

18

La séparation de phase, bien qu’elle puisse influencer fortement l’écoulement, ne sera pas

prise en compte dans le modèle hydraulique développé lors du stage. D’autant plus que Lewis et Lowell (2004), ont montré, à l’aide d’un modèle à une dimension, que si l’on suppose que le transfert de la chaleur se fait majoritairement par diffusion au sein du champ hydrothermal l’épaisseur de la couche à forte densité représente 1% de la hauteur caractéristique de l’écoulement.

1.3.2 Serpentinisation, production d’hydrogène et évolution de la

perméabilité

Sous le terme serpentinisation l’on désigne les réactions d’hydratation des péridotites du manteau au contact de l’eau de mer. Il existe plusieurs réactions de serpentinisation en fonction des réactifs disponibles. D’après Moody (1976):

- 5Mg2SiO4+ FeSiO4 + 9H2O -> 3Mg3SiO5(OH)4 + Mg(OH)2 + 2Fe(OH)2

Olivine Serpentine Brucite - 3Fe(OH)2 -> Fe3O4 + H2 + 2 H2O.

Magnétite Le plancher du site est constitué d’une forte proportion de serpentine. Les calculs

s’appuyant sur les réactions théoriquement possibles faisant intervenir l’eau de mer et la péridotite prédisent un pH plus élevé, une teneur en SiO2 plus faible et une quantité en H2 bien supérieure à celles de l’eau de mer. De plus, une concentration en potassium [K] très sensible au ratio [K]eau/[K] roche est attendue (Charlou et al., 2002). Cela est vérifié pour le site de Rainbow sauf pour le pH qui est très faible pour ce site.

Cette dernière différence permet d’émettre l’hypothèse que l’intervention de la réaction de

serpentinisation dans les conditions du site conduit à l’abaissement du pH. En effet Allen et al. (1996), en supposant l’équilibre de la réaction de serpentinisation atteint à Trainbow et Prainbow, trouvent un pH de 5 qui s’abaisse à 3 lors du refroidissement. De plus, d’après Charlou et al., (2002), la faible disponibilité de Si et son incorporation lors d’altération de minéraux tels que la serpentine explique que celui-ci reste à des concentrations faibles.

De même, Charlou et al., (2002) précisent que d’après l’observation sur le terrain (Apps,

1985 ; Conveney et al. 1987 ; Abrajano et al 1988, 1990), la forte teneur en CH4 et H2 des fluides du site de Rainbow est issue d’une serpentinisation du type :

Fe2SiO4 (ou FeSiO3) + H2O + C (ou CO2) -> Fe3O4 + 3Mg3SiO5(OH)4 + Mg(OH)2 + H2 + CH4. Olivine ou Orthopyroxène Magnétite Serpentine Brucite

Pour cette dernière affirmation, Charlou et al. (2002) s’appuient de même sur des travaux

expérimentaux (Berndt et al., 1996) et des réactions théoriques (Wetzel et Shock, 2000). Enfin, une probable augmentation du pH et une perte de fer risquent d’être observées en

relation avec la baisse de disponibilité de l’Orthopyroxène étant donné que c’est la réaction de serpentinisation faisant intervenir ce dernier minéral qui est cinétiquement prépondérante dans les conditions (Trainbow, Prainbow) et qui est responsable de la forte libération de H+. La réaction de serpentinisation faisant intervenir l’Olivine est bien plus lente (Charlou et al., 2002). La production d’H2 risque donc de diminuer dans le temps. De plus, la serpentinisation s’accompagne d’une augmentation de volume des minéraux (Mével, 2003) et peut donc être

Page 19: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

19

responsable d’une forte chute de perméabilité des roches et par conséquent des débits d’éjection des fumeurs. Cette chute peut avoir d’autres causes telles que la précipitation/dissolution d’autres minéraux ou encore la dilatation des roches.

La réaction de serpentinisation, aboutissant à la production d’hydrogène, prend donc une place importante d’un point de vue hydraulique et chimique.

1.4 Morphologie du site A proximité du site de Rainbow, l’unique champ hydrothermal ayant fait l’objet d’études

par réflexion sismique est le champ nommé Lucky Strike (Singh 2008). Ces études ont permis de déduire la présence éventuelle d’une chambre magmatique à 3 kilomètres de profondeur sous la croûte océanique, mais il est difficile de transposer ceci au site de Rainbow. De plus, les données bathymétriques ne sont pas assez précises pour rendre compte de l’infiltration de petites quantités de matières chaudes dans les fractures rocheuses qui peuvent modifier le système hydrothermal.

Les nombreuses données de terrains récoltées notamment par l’IFREMER , comme les

flux de chaleur et de matières des champs hydrothermaux à l’interface océan - croûte océanique, permettent la délimitation de zones dites de décharge. Par contre, la superficie de l’aire d’entrée de l’eau de mer au sein de la croûte n’est absolument pas connue. En effet, l’eau de mer peut parfois pénétrer la croûte à plusieurs kilomètres de son point de sortie. La délimitation de la zone de recharge est aussi dépendante de la présence de failles qui jouent certainement un rôle important dans l’acheminement de l’eau au sein de la croûte puisqu’elles constituent de prime abord un chemin préférentiel.

1.5 Conclusion : Les inconnues du système, intérêt de la modélisation

Les inconnues du système sont très nombreuses et dues à la multitude de contraintes imposées par le contexte géomorphologique : en effet, aux difficultés habituelles de description des aquifères continentaux fracturés, s’ajoutent celles liées à la quasi-inaccessibilité du milieu étudié. Ces inconnues sont liées à la nature des roches, à leurs agencements et à celui des pénétrations de masses rocheuses chaudes, à la présence ou non de déformations géologiques telles que des failles ou à l’existence ou non d’une chambre magmatique. Elles peuvent être liées à la nature du fluide qui circule, à la présence d’une éventuelle séparation de phase en profondeur, à la précipitation de minéraux ou encore au passage possible de l’eau de mer à l’état subcritique. Il y a aussi les interactions entre la roche et l’eau de mer qui dépendent évidemment des champs de températures et de pressions qui sont eux-mêmes inconnus en profondeur. Enfin, pour être exploitable, un tel site doit posséder une longévité importante, une trop forte dépendance de son existence aux conditions sismiques le rendrait aléatoirement exploitable.

L’enjeu de la modélisation dans la description de tels systèmes paraît donc évident : la

modélisation est une aide à la compréhension. Couplée à l’exploitation des données de terrain, elle peut permettre de confirmer ou d’infirmer certaines hypothèses.

Page 20: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

20

2. MISE EN PLACE DU MODELE

2.1 Hypothèses simplificatrices appliquées au modèle

La complexité des écoulements hydrothermaux au niveau des dorsales impose de simplifier le système, tant du point de vue de sa géométrie que des phénomènes physiques pris en compte afin de répondre aux attentes précitées. Pour cette première approche, seul le couplage entre l’écoulement et le transfert de chaleur a été pris en compte.

Pour pouvoir valider le modèle et la résolution de ses équations, il est nécessaire de

comparer les résultats de calculs aux calculs publiés dans la littérature pour des situations identiques simples, puis pour des systèmes plus complexes. Pour ce stage, le domaine modélisé sera systématiquement carré, uniformément maillé et sera soit entièrement fermé, soit ouvert à son sommet. Les propriétés du milieu seront simplifiées : par exemple, la perméabilité du milieu sera considérée comme constante sur l’ensemble du domaine (milieu homogène).

Le modèle reste donc un modèle assez simple et ne fait intervenir que trois équations de

transfert d’énergie ou de masse, et des équations d’états de l’eau qui ferment le système :

L’équation de Darcy qui lie la vitesse et l’énergie potentielle du fluide, représentée en hydrogéologie par la charge, lorsque les vitesses sont inférieures au mètre par jour au sein du milieu poreux.

Une équation de conservation de la masse. Une équation de transfert de la chaleur. Des équations liant la masse volumique ou la viscosité de l’eau de mer aux

champs de températures ou de pressions.

Bien que le modèle soit présenté ici comme étant simple, les résultats présentés dans ce rapport montrent combien les phénomènes mis en jeu sont complexes et toujours sujets de recherches poussées, ce qui justifie de telles hypothèses.

2.2 Analogie des systèmes hydrothermaux avec la convection de Rayleigh-Bénard

2.2.1 La convection de Rayleigh-Bénard

Dans un fluide dilatable, l’instabilité thermique naît de tout gradient de température : cette

dernière tend à créer une différence de densité. Cependant même en présence d’un gradient thermique et soumis à la gravité, un fluide dilatable peut être en état de stabilité mécanique. Pour un liquide soumis à un gradient de température vertical, la condition de stabilité est généralement donnée sous la forme (Landau et Lifchitz, 1989) :

p

dT g Tdz C

αααα− < , (1)

où g représente la gravité, α le coefficient de dilatation thermique du fluide, T la température et cp la capacité thermique à pression constante.

Page 21: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

21

Si cette condition n’est pas satisfaite, la convection apparaît. Ainsi, si on place un fluide entre deux plaques planes horizontales et que l’on applique un gradient de température en chauffant légèrement plus la plaque inférieure, à partir d’un certain écart de température (∆T), des rouleaux de convection se forment : ce résultat est appelé convection de Rayleigh-Bénard (cf. Figure 7).

g

Thaut (°C)

Tbas(°C) = T haut +∆∆∆∆T Plaque plane inférieure (chaude) Fluide (ρ,ρ,ρ,ρ,µ,k therm )

d

Rouleaux de

convection

Plaque plane supérieure (froide)

g

Thaut (°C)

Tbas(°C) = T haut +∆∆∆∆T Plaque plane inférieure (chaude) Fluide (ρ,ρ,ρ,ρ,µ,k therm )

d

Rouleaux de

convection

Plaque plane supérieure (froide)

Figure 7 : Représentation schématique de la convection de Rayleigh-Bénard.

Dans le cas des dorsales océaniques, le moteur de l’écoulement des fluides

hydrothermaux est la forte dépendance de la masse volumique du liquide vis-à-vis de la température. Les principes théoriques et expérimentaux issus de la convection de Rayleigh-Bénard sont donc tout à fait applicables au cas d’un champ hydrothermal dorsalien. La plaque supérieure du dispositif est alors remplacée par une épaisse couche d’eau à température constante et le fluide convecté évolue dans un milieu poreux.

2.2.2 Le nombre de Rayleigh : explication physique de la

convection On se place dans le cas d’une convection de Rayleigh-Bénard. Lorsque la colonne de

fluide chaud remonte vers les zones hautes plus froides du fluide, elle perd une partie de sa chaleur par diffusion, atténuant ainsi le gradient local de température. On peut associer un temps caractéristique au phénomène de diffusion, noté ττττth, et défini comme suit :

²

ththerm

dk

ττττ = , (2)

avec d la distance entre les plaques horizontales du dispositif et ktherm la diffusivité thermique du fluide.

De même, la remontée de cette cellule chaude va dépendre de la viscosité dynamique (µ0)

du fluide qui ralentit le mouvement de convection, du coefficient de dilatation thermique (α), du gradient thermique (∆T), de la masse volumique caractéristique (ρ0) et de la gravité qui, quant à eux, favorisent la convection, ainsi que de la hauteur caractéristique du phénomène (d). On en déduit un temps caractéristique de convection, noté ττττm et défini comme suit :

0

0

=∆m

µcte

gd Tττττ

ρ αρ αρ αρ α . (3)

Page 22: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

22

Or l’apparition d’un mouvement durable n’est possible que lorsque la durée de vie de la cause du mouvement est plus élevée que la durée de la manifestation de l’effet. Autrement dit, pour qu’il y ait convection, on doit vérifier :

3

0

0

∆=h

m therm

g d Tµ k cste

τ ρ ατ ρ ατ ρ ατ ρ αττττ , (4a)

3

0

0

∆= > = ctherm

g d TRa cste Ra

µ kρ αρ αρ αρ α

, (4b)

où Ra représente une première écriture du nombre de Rayleigh. Rac correspond au nombre de Rayleigh critique au-delà duquel il y a convection.

En milieu poreux, le temps caractéristique de convection est défini de façon légèrement

différente :

0

0m

µ dcste

k gττττ

ρρρρ=

∆ . (5)

Le nombre de Rayleigh devient alors :

0

0

∆=therm

gdkRa

µ kρρρρ

. (6)

Ce nombre pilote l’écoulement (cf. § 4.).

2.2.3 Convection de Rayleigh-Bénard en milieu poreux : description des résultats courants

Dans les années 1980-1990, la modélisation de la convection de Rayleigh-Bénard au sein

de systèmes poreux a été l’objet de très nombreuses études. En effet, l’apparente simplicité de description d’un tel système instable a conduit de nombreux chercheurs à s’intéresser aux mécanismes de transition vers le chaos d’écoulements de fluides mis en mouvement par un différentiel de température. Ces études ont été menées le plus souvent sur un domaine carré fermé chauffé par le bas, possédant des parois latérales adiabatiques.

A partir d’un Rayleigh critique, égal à 4π² (Kimura et al. 1986, Caltagirone et Fabrie

1989, Graham et Steen 1992 et 1994, Cherkaoui et Wilcock 1999), on observe la mise en mouvement du fluide. La convection déforme alors le champ de températures.

Kimura et al. (1986), Steen et Aidun (1988) ou encore Caltagirone et Fabrie (1989)

ont étudié les mécanismes des transitions observées en fonction du nombre de Rayleigh. Après la première bifurcation du système vers un régime permanent par dépassement du Rayleigh critique, l’augmentation progressive du Rayleigh fait apparaître différents régimes repérables à l’aide de la courbe d’évolution temporelle du flux de température en limite inférieure du domaine. Ce flux est aussi appelé nombre de Nusselt et est défini par :

Page 23: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

23

**

** 0

0

|L

z

TNu dx

z =

∂=∂∫ , (7)

où L* est la longueur de la base du domaine adimensionnée et T* la température adimensionnée.

Selon Steen et Aidun (1988), lorsque le nombre de Rayleigh est supérieur à 390,7, le système quitte le régime permanent pour atteindre un premier régime d’oscillations périodiques (état P1). Puis pour une gamme de Rayleigh de 500 à 565, le régime devient quasi-périodique (état QP1) : deux fréquences indépendantes d’oscillations apparaissent puis se confondent vers Ra = 565. Le régime devient à nouveau périodique (état P2) jusqu’à atteindre un Rayleigh de 850 où une troisième fréquence d’oscillation indépendante des deux autres intervient (QP2). Par la suite, d’autres fréquences d’oscillation apparaissent jusqu’à l’apparition d’un régime chaotique.

De plus, il est important de noter que Caltagirone et Fabrie (1989) ont réussi à maintenir

un champ de températures correspondant à un seul rouleau de convection en augmentant le nombre de Rayleigh du régime périodique P1 jusqu’à un nombre de Rayleigh de 1500. Dans ce même article, les auteurs posent la question de l’existence ou non d’un régime chaotique.

2.2.4 Pourquoi le nombre de Nusselt oscille-t-il ?

En Annexe A est présentée l’évolution du champ de températures dans un milieu poreux

fermé, chauffé par le bas, pour un Rayleigh de 800, au cours d’une période de temps égale à l’inverse de la fréquence d’oscillation du nombre de Nusselt. On observe donc une convection de Rayleigh-Bénard dans le cas d’un régime périodique de type P2 (la fréquence principale f1 est de l’ordre de 280).

On remarque aisément à la base du domaine, un bombement de la couche limite qui

s’intensifie à mesure qu’il se déplace de la droite vers la gauche dans le sens général de l’écoulement, puis donne naissance à un bourgeonnement net de matière chaude à mi-parcours. Ce bourgeon, tout en se déplaçant vers la gauche possède une vitesse ascensionnelle liée à sa plus faible densité. (cf. Annexe B)

Par la suite, cette vague vient s’écraser contre la paroi latérale gauche, chasse une zone

plus froide pincée entre cette vague et la vague précédente qui remonte vers le sommet du domaine. Le système réagit symétriquement pour la zone froide aux perturbations engendrées par les vagues de chaleurs. Ces vagues correspondent à une tentative du système de mise en place d’un nouveau rouleau de convection afin d’évacuer le surplus d’énergie. Mais tant que le Rayleigh n’est pas supérieur à une valeur seuil, cette tentative est étouffée par le mouvement général de convection de l’unique rouleau.

Selon Caltagirone et al. (1989), les différentes fréquences apparaissant par exemple pour le régime QP1 sont dues aux différences entre la vitesse d’apparition du bourgeonnement et la vitesse d’étalement des cellules ascendantes (respectivement descendantes) le long de la couche limite froide (respectivement chaude). De plus, l’apparition d’une troisième fréquence serait due à la coalescence des vagues lors de leurs ascensions le long des parois latérales.

La validation du modèle développé lors de mon stage a donc dû tenir compte du nombre

de Nusselt, du nombre de Rayleigh et de la forme du champ de températures.

Page 24: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

24

2.2.5 Conclusion Les travaux présentés précédemment sont essentiels à l’étude des écoulements

hydrothermaux. D’un point de vue pratique, la méthode de validation des codes développés par les chercheurs qui modélisent les systèmes hydrothermaux au niveau des dorsales consiste à comparer leurs résultats de calculs avec des résultats déjà publiés. Puis le domaine modélisé, auparavant fermé, est ouvert pour simuler les échanges à l’interface océan-croûte océanique. Cependant, le comportement du système ouvert est assez différent de celui du système fermé et une telle validation des codes ne garantit pas l’exactitude des solutions trouvées pour un système ouvert.

Cherkaoui et Wilcock (2001) ont tenté de réaliser des expériences de Rayleigh-Bénard

dans des cellules de Hele-Shaw de façon à modéliser expérimentalement des systèmes hydrothermaux. Ils expliquent la différence entre les résultats des simulations numériques et les résultats de leurs propres expériences notamment par la présence de parois non-adiabatiques. Le Rayleigh maximal exploré étant égal à 500, ils affirment donc que l’écoulement est resté laminaire. Or, par la suite, dans ce rapport, étant données les perméabilités et les conditions de température régnant au sein des systèmes hydrothermaux qui nous intéressent, le problème des hauts nombres de Rayleigh sera abordé. Cette expérience ne pourra donc pas servir à la validation des modèles développés lors de ce stage.

Le modèle mathématique correspond à un ensemble d’équations qui doivent traduire au mieux l’observation expérimentale ou de terrain. Les équations développées dans les prochains paragraphes sont des équations classiquement utilisées pour la modélisation des systèmes hydrogéologiques. Le nombre des équations développées est fonction du nombre d’inconnues prises en compte et donc de la précision de la description des systèmes étudiés.

2.3 Modèle mathématique général

2.3.1 Modèle hydraulique

La loi de Darcy Henri Darcy, ingénieur hydraulicien du XIXe siècle, proposa dans un rapport de 1856 (de

Marsily 2004), la loi suivante :

= =QU K i

A, (8)

où Q représente le débit hydraulique sortant (m 3.s-1), A, la section du cylindre (m²) dans lequel a lieu l’écoulement en milieu poreux, U, le débit moyen (m.s-1), K, une constante de

proportionnalité (m.s-1) et i la pente hydraulique : h

iL

∆= , avec h, la charge (m) et L, la longueur

du cylindre.

On écrit généralement cette loi sous une forme vectorielle :

( )U K gradh= −ur uuuuur

. (9)

Page 25: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

25

Cette loi n’est valide que sous certaines conditions. En effet, lorsque le gradient hydraulique est très faible, le coefficient de perméabilité ne peut plus être considéré linéaire et, lorsque le gradient hydraulique est trop important, un terme d’inertie intervient. La limite haute de validité de la loi de Darcy est souvent déterminée en fonction du nombre de Reynolds défini par :

ReΩ

=

urU ρρρρ

µµµµ, (10)

avec Ω, la taille caractéristique des pores, ρ, la masse volumique du fluide et µ, la viscosité du fluide. Lorsque le Reynolds est inférieur à 1, la loi de Darcy s’applique. On vérifiera par la suite que ces conditions sont remplies.

Enfin, étant données les grandes échelles de temps des écoulements étudiés, bien que la

loi de Darcy ne soit normalement valable qu’en régime permanent, on considérera que la variation de U par rapport au temps est négligeable par rapport au gradient de charge hydraulique. Cette simplification est utilisée dans la majeure partie des publications traitant de ce sujet.

Equation de continuité

L’équation de continuité, également appelée équation de conservation de la masse, s’obtient en écrivant le bilan entrée-sortie de masse de fluides sur un volume élémentaire :

( )t

div U Qρφρφρφρφ ρρρρ∂ + =∂

ur

, (11)

où Q est un terme puits/source, ρ représente la masse volumique de l’eau, U la vitesse de Darcy et φ la porosité du milieu.

2.3.2 Modèle de transfert de la chaleur

Si dans le manteau la chaleur est évacuée presque exclusivement par convection, au sein

du système hydrothermal, a priori, convection et diffusion peuvent toutes deux jouer des rôles importants. Au sein de la matrice rigide de la croûte, la conduction est largement majoritaire tandis qu’on ne peut négliger aucun des deux phénomènes de transport précités au sein du fluide hydrothermal. L’équation mathématique modélisant le transfert de la chaleur doit donc traduire l’existence des modes de conduction mais aussi des échanges de chaleur entre chacune des phases.

Pour la phase fluide, on peut donc écrire :

( ) ( (( ) ) ( ) )ll diffl displ l l l l

Tc div D D gradT c UT Q

tφ ρ φ ρφ ρ φ ρφ ρ φ ρφ ρ φ ρ∂ − + − =

uuuuur ur

, (12)

où Ql, représente le terme puits/source, (ρc)l, la capacité thermique volumique du fluide, T l, la température du fluide et U, la vitesse de Darcy. Ddiffl représente le coefficient de diffusion tandis que Ddispl représente le coefficient de dispersion liée aux hétérogénéités microscopiques du champ de vitesses.

Page 26: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

26

Pour la phase solide, cette équation se résume à :

(1 )( ) ((1 ) )ss diffs s s

Tc div D gradT Q

tφ ρ φφ ρ φφ ρ φφ ρ φ∂− − − =

uuuuur

, (13)

où Ts, est la température du solide, Ddiffs, est la diffusivité thermique du solide, (ρc)s, la capacité thermique volumique du solide et Qs, un terme puits/source.

Comme l’unique terme puits (respectivement source) du solide est le terme source

(respectivement puits) du fluide, s lQ Q= − . En injectant l’équation (12) dans l’équation (13) et en considérant que la roche et l’eau sont sensiblement à l’équilibre, l’équation du transfert de la chaleur peut alors s’écrire :

( ( ) (1 )( ) ) (( ( ) ))l s ldiv gradT UTt

φ ρ φ ρ λ ρφ ρ φ ρ λ ρφ ρ φ ρ λ ρφ ρ φ ρ λ ρ∂Τ+ − = −∂ eqc c c

uuuuur ur

, (14a)

( ) (( ( ) ))∂Τ = −∂ eqc c

uuuuur ur

eq ldiv gradT UTt

ρ λ ρρ λ ρρ λ ρρ λ ρ , (14b)

avec λeq la conductivité thermique du milieu saturé qui prend en compte un terme de diffusion et un terme de dispersion de la chaleur liée aux hétérogénéités du champ de vitesses du fluide convecté :

( ) (1 )eq diffl displ diffsD D Dλ φ φλ φ φλ φ φλ φ φ= + + − . (14c)

ρc correspond à la capacité thermique volumique équivalente du milieu saturé :

( ) ( ) (1 )( )= + −c c ceq l sρ φ ρ φ ρρ φ ρ φ ρρ φ ρ φ ρρ φ ρ φ ρ . (14d)

2.3.3 Equations d’états

Ces équations traduisent le comportement de l’eau placée dans des conditions

particulières, c'est-à-dire qu’elles lient les paramètres thermodynamiques intensifs du fluide (pression et température) à ses paramètres extensifs.

La diversité des systèmes hydrothermaux, la variabilité des températures aux limites, des

profondeurs ou encore des perméabilités qui existent même entre deux champs contigus, conduit l’eau de mer à passer par différents états. Dans certains systèmes, en profondeur, l’eau s’écoule à l’état subcritique, dans d’autres, une séparation de phase conduira à la formation d’une phase saumure, dense, et d’une phase aqueuse moins dense.

Page 27: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

27

2.4 Paramètres du modèle

2.4.1 Paramètres de l’écoulement

Paramètres propres au milieu poreux / fracturé

o La perméabilité intrinsèque

La facilité qu’a l’eau à s’écouler au sein d’un milieu poreux ou fracturé est représentée par le concept de perméabilité du milieu. Pour un milieu poreux, cette capacité de pénétration de l’eau va dépendre de la nature de la roche (sable, basalte …), du diamètre et des formes des grains solides, ou encore de la rugosité ainsi que des caractéristiques propres du fluide : sa viscosité et sa densité. Les fractures constituent des chemins préférentiels d’écoulement où la vitesse du fluide est plus élevée qu’au sein même de la roche. La perméabilité augmente donc avec le diamètre et la fréquence spatiale des fractures.

La nature du fluide variant beaucoup tout au long de l’écoulement hydrothermal, souvent,

la perméabilité mesurée et utilisée par les modèles numériques est la perméabilité intrinsèque (kint). Celle-ci ne dépend que de la nature du milieu fracturé et non du fluide qui le traverse : k int = K×µµµµ/ρρρρg (où µ est la viscosité du fluide, ρ sa masse volumique, g la gravité, K la perméabilité de Darcy). Cette formulation autorise la mesure de perméabilités intrinsèques très faibles via l’utilisation de gaz (voir Singh et al. (2008)). Elle permet donc de modéliser l’évolution de la perméabilité du milieu à l’aide d’un produit de variables indépendantes.

D’après Mével et al. (2003), la structure de la croûte océanique au niveau des dorsales

lentes n’obéit pas au modèle PENROSE de lithosphère litée, c'est-à-dire que la superposition des couches de Basaltes puis Gabbros et enfin les péridotites mantelliques – la discontinuité de MOHO étant située entre les Gabbros et les péridotites- n’est pas retrouvée pour ce type de dorsale. En effet, on trouve généralement un milieu hétérogène avec des variations de nature des péridotites plus ou moins serpentinisées ainsi que des inclusions gabbroïques. Les reliefs et les structures sont elles-mêmes variables d’un champ hydrothermal à l’autre.

Obtenir la caractérisation expérimentale d’un champ de perméabilités défini en tout point

du domaine modélisé paraît donc matériellement impossible. Cela imposerait la réalisation d’une multitude de forages pour chaque site étudié. Certaines méthodes, telle que l’approche inverse, permettent de reconstituer des champs de perméabilités à partir de la mesure expérimentale de charges ou de flux.

o La porosité

Nous l’avons vu précédemment, porosité et perméabilité sont étroitement liées. La

porosité représente le volume de vide sur le volume total de la roche étudiée. Les roches de la croûte réagissent au champ de températures en se dilatant, en se contractant ou, lorsque le choc thermique est trop grand en se fracturant. De même, l’activité sismique peut tenir un rôle majeur dans la modification des paramètres physiques des roches de la croûte océanique. La porosité est donc un paramètre très difficilement ajustable et lorsque l’on cherche à modéliser un système hydrothermal, afin de s’affranchir de cette incertitude, une porosité constante et l’approximation de Boussinesq sont souvent considérées (cf. § 2.5.1 et 4.3.1).

Page 28: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

28

Paramètres du fluide : masse volumique et viscosité La masse volumique et la viscosité de l’eau de mer dépendent de la concentration en

éléments solubles dans l’eau, de la température et de la pression. La dépendance de la viscosité en fonction de la pression est très souvent considérée négligeable vis-à-vis des deux autres dépendances.

2.4.2 Paramètres du transfert thermique

Capacité thermique volumique

La capacité thermique volumique ρρρρc correspond à la quantité d'énergie à apporter pour élever d'un degré la température de l'unité de volume d'une substance. Elle est égale au produit de la capacité thermique massique et de sa masse volumique. Souvent une capacité thermique équivalente (ρc)eq est définie pour l’ensemble du milieu poreux saturé. On considérera par la suite que le rapport (ρc)eq/(ρc)l est constant et proche de 1 (cf § 2.5.3, relation (18)).

Conductivité thermique La conductivité thermique (λ) représente la quantité de chaleur transférée par unité de surface pour un gradient de température imposé. Elle dépend évidemment de la nature du matériau, mais aussi de la température.

Diffusivité thermique La diffusivité thermique équivalente du milieu poreux saturé est déterminée par :

( )λλλλ

ρρρρ= eq

therm

eq

kc . Ce paramètre prend toute son importance en régime transitoire car il traduit la

vitesse de changement local de température du système.

2.4.3 Hauteur caractéristique du domaine La hauteur caractéristique du domaine ainsi que le rapport d’aspect (A = largeur/hauteur)

vont jouer sur le nombre de rouleaux de convection observés ainsi que sur la forme de ces rouleaux.

2.5 Simplification du modèle mathématique

Les équations présentées précédemment sont simplifiées pour permettre de mieux

identifier les causes des phénomènes observés. Ces simplifications sont utilisées par la plupart des auteurs. Celles-ci imposent la mise en place d’hypothèses fortes sur le comportement du système.

2.5.1 Formulation spécifique et simplification du modèle

hydraulique Dans l’équation de Darcy (9), il est possible de faire apparaître la masse volumique afin

de prendre en compte explicitement la dépendance de la vitesse vis-à-vis de la densité. La viscosité du fluide est supposée constante et on écrit :

Page 29: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

29

( )k

U gradP ggradzµ

= − + ρur uuuuur uuuuur

, (15)

où P représente la pression totale (Pa), k la perméabilité intrinsèque (m²), µ la viscosité du fluide (kg.m-1.s-1), ρ la masse volumique et g la gravité.

Comme cela est fait par de nombreux auteurs, on suppose que l’approximation de

Boussinesq est valable. La masse volumique du fluide ne dépend que de la température et n’intervient que dans le terme gravitaire de l’équation de Darcy. En supposant par ailleurs, que la variation temporelle de la porosité est négligeable, l’équation de continuité (11) se réduit alors à :

( ) 0div U =ur

. (16)

2.5.2 Evolution de la masse volumique

L’approximation de Boussinesq étant supposée valable, on considère une évolution linéaire de la masse volumique en fonction de la température :

( )( )0 0ρ = ρ 1 - T - Tαααα , (17)

où T0 est la température de référence pour laquelle ρ = ρ0.

2.5.3 Formulation pour le calcul du transport de la chaleur En considérant l’équilibre thermique vérifié, on peut écrire Tf = Ts = T. L’équation du transport thermique (14b) devient :

(( ))div gradT UTt

∂Τ = −∂ thermk

uuuuur ur

. (18)

2.5.4 Adimensionnalisation des équations

Il est en général intéressant d’adimensionner les équations afin de faire apparaître des

nombres caractéristiques tels que le nombre de Rayleigh et de restreindre la dépendance du résultat à un nombre limité de paramètres. Pour ce modèle simple, le milieu sera considéré homogène, l’adimensionnalisation est donc presque immédiate. On pose les constantes d’adimensionnalisation suivantes :

max minT T T∆ = − ,

min max( ) ( )T Tρ ρ ρρ ρ ρρ ρ ρρ ρ ρ∆ = − , *x d x= × , *y d y= × ,

*V K V= × , *T T T= ∆ × , *ρ ρ ρρ ρ ρρ ρ ρρ ρ ρ= ∆ × ,

*P gd Pρρρρ= ∆ × , . ²

*=therm

dt t

k

Page 30: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

30

Le nombre de Rayleigh est ici défini par la relation (6).

Les équations obtenues sont :

La vitesse de Darcy : * * *( )U grad P g

gρρρρ ∗

= − +ur uuuuur ur

. (19)

L’équation de conservation de la masse : ** ( ) 0div U =

ur

. (20)

L’équation de transfert de la chaleur : *

* * * * * **

. ( ) ( )T

RaU grad T div grad Tt

∂ + =∂

uuuuuur uuuuuurur

. (21)

L’équation de dépendance de la masse volumique : ( )( )0* * *

0ρ = ρ 1 - T T - Tαααα∗ ∆ . (22) Plusieurs intérêts existent à cette adimensionnalisation, elle permet :

la comparaison de nos résultats à d’autres publications pour un même Rayleigh ; d’éviter de répéter la simulation de scénarii similaires, par exemple cas d’une perméabilité

intrinsèque et d’une viscosité toutes deux multipliées par un même facteur ; de fixer clairement les limites de stabilité du système à partir d’un très petit nombre de

paramètres.

2.5.5 Non-linéarité du système

La formulation mathématique précédente met en évidence la forte non-linéarité du problème traité. En effet, le terme convectif de l’équation (21) du transfert de la chaleur est le produit de deux fonctions de la température (on rappelle que U est fonction de T via le terme de densité). La linéarisation du système par une méthode de Picard a été abordée lors de ce stage mais n’est pas développée dans ce rapport. Le traitement de cette non linéarité est présentée au paragraphe 4.1.2.

Page 31: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

31

3. OUTILS NUMERIQUES

3.1 Outil de calcul en langage Gibiane : CAST3M L’outil de calcul CAST3M a été développé à la Direction de l’Energie Nucléaire (DEN)

du CEA afin de répondre aux problèmes de mécanique des structures dans le domaine du nucléaire. Depuis, de nombreuses autres fonctionnalités ont été développées pour le calcul de problèmes de mécanique des fluides, d’hydrogéologie, d’électromagnétisme ou encore de thermique.

Ecrire sous CAST3M, c’est utiliser des directives via une syntaxe précise en langage

Gibiane. Ces directives font appel à des opérateurs codés en langage Esope, langage proche du Fortran. La simplicité d’utilisation cache néanmoins un problème, il est indispensable de connaître la nature exacte des opérations réalisées par les outils CAST3M. Lorsque l’on cherche à améliorer un schéma numérique, on est alors amené à décortiquer chacune des procédures codées en Esope, CAST3M perd alors son côté ludique.

Un calcul se déroule généralement en trois étapes :

La définition du modèle mathématique, c'est-à-dire le choix des équations gouvernant le système, de la géométrie et du maillage, des conditions initiales et aux limites, des propriétés matérielles etc.…

L’écriture du jeu de données : la construction du maillage, la discrétisation des équations, la

prise en compte des données caractérisant le système, le calcul. Le post-traitement qui correspond à l’exploitation numérique et graphique des résultats du

calcul.

CAST3M est un des logiciels de la plate-forme numérique ALLIANCES conjointement développée par le CEA, l’ANDRA et EDF. Au sein de cette plate-forme, CAST3M a été couplé au code de géochimie CHESS développé par le Centre de Géosciences de l’Ecole des Mines de Paris. C’est là tout l’intérêt de la modélisation hydrothermale sous CAST3M. Une fois les paramètres hydrauliques maîtrisés, il sera donc possible de prendre en compte la géochimie.

Dans de nombreuses études réalisées au sujet de la convection de Rayleigh-Bénard en milieu poreux, les codes développés sont des codes pseudo-spectraux fournissant une bonne précision de calcul. CAST3M ayant été développé dans le but de répondre à des problèmes de mécaniques, il n’existe pas de module de création de codes pseudo-spectraux.

La méthode numérique généralement utilisée sous CAST3M pour traiter des problèmes

hydrogéologiques est celle des éléments finis mixtes hybrides (EFMH). Cette méthode permet la résolution simultanée de l’équation de Darcy et de conservation de la masse. Par construction, elle permet d’assurer la continuité des flux aux interfaces entres mailles adjacentes et donc de respecter la conservation de la masse. Ce sont donc les EFMH qui ont été testés en premier lieu. La méthode des Volumes Finis (VF) a été plus récemment implémentée dans CASTEM, elle est d’ailleurs toujours activement développée.

Page 32: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

32

3.2 Modélisation à l’aide des EFMH d’un champ hydrothermal dorsalien

3.2.1 Choix du schéma spatial de discrétisation

La méthode de discrétisation des équations de transport est identique à celle appliquée pour la discrétisation de l’équation de Darcy et de la conservation de la masse (Dabbene 1993, 1995, Dabbene et al. 1994, Paillere et al. 1997). Le terme diffusif de l’équation de diffusion-convection pour le transport de la chaleur joue le rôle de la vitesse de Darcy :

* *U K gradh Q grad T= − ≡ = −

uuuuuuruuuuur

. (23) L’équation de diffusion-convection prend la place de l’équation de continuité, le terme

convectif étant traité comme un terme source.

3.2.2 Choix du schéma temporel de discrétisation Dans cette étude, sous CAST3M, la méthode de discrétisation temporelle est une théta-

méthode :

1 1

( ( ))

n n

n n

t t

t t

dTdt f T t dt

dt

+ +

=∫ ∫ , (24)

1

11

( ) (1 ) ( )n n

n nn n

T Tf T f T

t tθ θθ θθ θθ θ

++

+

− = + −−

. (25)

Pour passer d’un schéma implicite à un schéma explicite ou de Crank-Nicholson, il suffit

de changer la valeur du paramètre θ. D’une manière générale, dans cette partie, les deux schémas testés seront un schéma implicite en temps (θ = 1) et le schéma de Crank-Nicholson (θ = 1/2).

Le schéma implicite en temps est un schéma du premier ordre en temps :

1

1( )n n

nT Tf T

t

++− =

∆. (26)

Le schéma temporel de Crank-Nicholson est un schéma du second ordre en temps :

1 1( ) ( )

2

n n n nT T f T f Tt

+ +− +=∆

. (27)

3.2.3 Précision de calcul et stabilité : les critères

Certains critères ont été développés afin de pouvoir se placer a priori dans des conditions

de pas en temps et de pas en espace qui offrent la meilleure précision de calcul. Ces critères sont des nombres caractéristiques dépendant des données décrivant le modèle et le système. Les deux schémas numériques présentés précédemment doivent obéir à ces derniers.

Page 33: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

33

o Coefficient de Fourier

Pour le calcul du phénomène de diffusion en EFMH, la précision est déterminée par le coefficient de Fourier (Fo), par exemple pour le cas du transfert thermique :

2

²thermk t

Fox

∆=∆ , (28)

où ∆t est le pas de temps et ∆x, le pas d’espace du maillage. Ce coefficient est légèrement différent lorsque le problème est adimensionné :

2 *

* ²

tFo

x∆=

∆ . (29)

La précision du calcul est normalement optimale pour Fo = 1 (Maugis 2006).

o Coefficient de Courant-Friedrichs-Lewy

Le calcul du terme convectif à l’aide des EFMH pour un schéma implicite en temps implique une stabilité quelle que soit la valeur du critère de Courant-Friedrichs-Lewy (CFL). Le CFL est alors exclusivement un critère de précision. Celle-ci diminue à mesure que l’on s’éloigne d’un CFL de valeur 1. Le CFL est défini par :

* *

* therm

U t U t kCFL

x x k

∆ ∆= = ×

∆ ∆

uuur uuuur

. (30)

o Nombre de Péclet

Le nombre de Péclet traduit l’importance du phénomène de transfert par convection vis-à-

vis du transfert par diffusion. Il est défini par :

2CFL

PeFo

= . (31)

3.2.4 Discussion a priori sur la stabilité du schéma et la précision

du calcul : pas d’espace et pas en temps

Le système physiquement instable décrit est sujet à des oscillations périodiques du flux de chaleur aux bornes du domaine. La fréquence, l’amplitude et la valeur moyenne de ces oscillations sont dépendantes du Nombre de Rayleigh. Or Frédéric Dabbene, ingénieur numéricien au CEA, avait posé comme règle pratique que 80 pas de temps de calcul devaient s’écouler par demi-période d’oscillation afin de décrire précisément un phénomène de ce type. Les fréquences de ces oscillations varient de 30 à 500 pour la gamme de nombres de Rayleigh étudiée lors de la validation du code. Par exemple, une fréquence de 500, le pas de temps adimensionné maximal est d’environ 1,25×10-5. Hoteit et al. (2002) ont montré que, lorsqu’on utilise la méthode des EFMH pour résoudre une équation de transport diffusif, des oscillations peuvent apparaître en dessous d’une valeur seuil du nombre de Fourier, cette valeur seuil étant égale à 1/3. Ceci est d’autant plus contraignant que le nombre de Péclet est généralement très faible dans les cas étudiés. Le

Page 34: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

34

phénomène de diffusion peut donc occuper une place importante dans l’évolution du système. Aussi, on peut considérer que malgré la prise en compte de la convection dans nos calculs, le critère présenté par Hoteit et al. (2002) doit être retenu. Le tableau suivant présente une comparaison des gammes de pas de temps utilisables issus des critères précédemment présentés en fonction du maillage pour une fréquence d’oscillation de 500.

Maillage homogène carré : Fourier minimum suivant le critère d’Hoteit et al. (2002)

Fourier maximum suivant le critère du nombre de pas de

temps par période 64 × 64 1/3 0,10 96 × 96 1/3 0,23

128 × 128 1/3 0,41 256 × 256 1/3 1,64

Tableau 4 : Comparaison des gammes de pas de temps utilisables pour une fréquence d’oscillation de 500.

On constate qu’un maillage constitué de 64×64 mailles carrées est inutilisable pour l’obtention de résultats précis. Par contre, un maillage constitué de 128×128 mailles carrées paraît être un bon compromis mais peut conduire à des calculs très longs d’une à deux semaines. Dans cette étude, sauf dans le cas de la validation en boîte fermée pour un nombre de Rayleigh supérieur à 800, la fréquence d’oscillation du nombre de Nusselt est inférieure à 280 ce qui rend l’utilisation d’un maillage constitué de 96×96 mailles possible, le fourrier maximum étant alors de 0,41. Pour des fréquences d’oscillation plus élevées, ce même maillage permet l’obtention de résultats tout à fait satisfaisants (cf. § 3.4).

3.3 Modélisation à l’aide des VF d’un champ hydrothermal dorsalien

3.3.1 Choix du schéma spatial de discrétisation Seule l’équation du transfert de la chaleur est ici discrétisée à l’aide de la méthode des volumes finis. Les équations de conservation de la masse et de la vitesse du fluide restent discrétisées à l’aide de la méthode des EFMH.

3.3.2 Choix du schéma temporel de discrétisation Le schéma temporel de Crank-Nicholson a été le seul utilisé ici.

3.3.3 Les raisons d’une telle discrétisation D’après Bernard-Michel et Genty (2006), contrairement aux EFMH, dans le cas des VF, c’est la déformation de l’élément qui a une importance et non le nombre de Fourier de maille. Or l’approximation de Boussinesq revient à discrétiser un système d’équations hydrauliques définies pour un régime permanent. Donc l’utilisation des EFMH pour la partie hydraulique du système et des VF pour la partie transfert thermique permet de s’affranchir des contraintes sur le pas de temps.

Page 35: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

35

Calcul hydraulique

Un+1

Définition des conditions aux limites et initiales, et définition des paramètres

Calcul du transport thermique

Tn+1

Calcul de la masse volumique

ρρρρn+1

Création du maillage

3.4 Algorithme de résolution Au paragraphe 2.5.5 a été mis en évidence l’existence d’une forte non linéarité liée à la prise en compte du phénomène de convection au sein du modèle mathématique. De façon simple, cette non linéarité est traitée en résolvant de façon séquentielle les équations hydrauliques puis l’équation de transfert de la chaleur. Ceci est possible si l’on considère, dans le cadre de l’approximation de Boussinesq, que les transmissions de la pression, et donc de la vitesse sont instantanées par rapport au transfert de la chaleur.

Figure 8 : Algorithme de résolution utilisé lors de cette étude.

Le terme convectif de l’équation du transfert de la chaleur peut ainsi être considéré

comme terme source. Toutefois ce type de linéarisation induit des erreurs discutées au paragraphe 4.3.3.

ρρρρn+1 Tn+1

Mise en mémoire du champ de température

Variables entrantes à la n-ième itération en temps ρρρρn, Tn

Pas de temps suivant

Page 36: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

36

4. VALIDATION DU JEU DE DONNEES – CAS-TEST

Après l’écriture des jeux de données sous CAST3M, ceux-ci doivent être validés par comparaison des résultats de calculs avec des résultats issus de la littérature. Cette validation est primordiale, car les jeux de données développés vont ensuite servir à modéliser des systèmes dont le comportement n’est pas clairement connu.

4.1 Validation des modèles numériques en domaine fermé

4.1.1 Description des conditions initiales et des conditions aux limites

La configuration étudiée est dans un premier temps une boîte carrée bidimensionnelle

fermée de 1500 m de côté dont les conditions aux limites sont les suivantes :

Une condition de symétrie est imposée à chaque paroi latérale du domaine. Les paramètres du système sont regroupés dans le tableau suivant :

Valeurs Accélération gravitaire d’axe y : g 9,81 m.s-2

Masse volumique de l’eau de mer à 300 bars et 2°C : ρ0 130 kg/m3

Coefficient de dilatation volumique: α 1,17×10-2K-1 Viscosité de l’eau pure : µ 5×10-5 kg.m-1.s-1

Capacité thermique volumique de l’eau : (ρc)f 4,2×10-6 J.k-1.m-3

Capacité thermique volumique du milieu saturé : (ρc)eq ≈ 4,2×10-6 J.k-1.m-3

Conductivité thermique du milieu saturé λeq 2,5 W.m-1.K-1

Diffusivité thermique ktherm 6×10-7 m²/s Perméabilité intrinsèque du milieu poreux k : 4×10-13 à 4×10-15

m² Tableau 5 : Paramètres du système.

Figure 9 : Schéma du domaine considéré et conditions aux limites.

Figure 10 : Schéma du domaine considéré et conditions aux limites adimensionnées.

T = 1, u = v = 0, L = 1

T= 0, P= 2,25, u = v = 0

u = v = 0 d = 1

u = v = 0 Milieu poreux homogène isotrope

T = 600 °C, u = v = 0, L = 1500m

T= 2°C, P = 300 bars, u = v = 0

u = v = 0, d = 1500m

u = v = 0 Milieu poreux

homogène isotrope z

x

Page 37: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

37

Deux sortes de conditions initiales ont été utilisées selon les cas modélisés. La première est notée CI1 et, si Tdiff (z) est la température du système dans le cas d’un transport uniquement diffusif, la condition CI1 est définie par CI1 : T(t=0) = Tdiff (z) + 0,01 (cos(x) sin(z)) dans le cas du domaine adimensionné. Cette infime perturbation du champ de températures initiale est un artifice permettant de choisir le sens de rotation du rouleau de convection dans le cas du domaine fermé., ceci dans le but de pouvoir comparer aisément les champs de températures d’un calcul à l’autre. La deuxième condition initiale correspond au champ de températures obtenu lors de l’établissement du régime final pour une simulation réalisée avec un nombre de Rayleigh donné. Elle sera notée par la suite CIF suivi du nombre de Rayleigh (CIF800 si Ra = 800). Dans ce cas, par abus de langage, on parlera de croissance progressive du nombre de Rayleigh.

4.1.2 Résultats obtenus à l’aide de la méthode des EFMH

Conformément aux critères présentés précédemment en partie 3.2.4, des simulations ont été réalisées d’une part avec un maillage constitué de 96×96 mailles carrées, avec un nombre de Fourier variant de 0.34 à 0.42 et, d’autre part, avec un maillage constitué de 128×128 mailles carrées avec un nombre de Fourier variant de 0.34 à 1.

o Résultats obtenus pour un nombre de Rayleigh de 800

Résultats qualitatifs

Pour ce nombre de Rayleigh, une condition initiale de type CI1 a été utilisée. Pour chacun des schémas numériques testés, les résultats graphiques sont sensiblement les mêmes.

• Champ de températures

Les champs de températures obtenus avec CAST3M et ceux de Fontaine et al. (2007),

pour un nombre de Rayleigh de 800, sont donc sensiblement identiques. On remarque un étalement des masses chaudes (respectivement froides) à la fin de leur remontée (respectivement descente) le long de la paroi verticale. Au centre du domaine, la température est globalement homogène. Une couche limite chaude en partie inférieure et froide en partie supérieure du

Figure 11 : Champ de températures obtenu avec CAST3M (EFMH).

Figure 12 : Champ de températures publiés dans (Fontaine et al., 2007).

478 °C 359 239 120

478 °C 359 239 120

Page 38: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

38

domaine laissent entrevoir un léger bombement, puis un bourgeonnement, signe de l’amplification d’instabilité née au sein de cette couche.

• Champ de Pressions

Le champ de pressions est, en général, peu analysé par les auteurs car il est, comme on

peut le voir en Annexe C, peu déformé et ne permet pas une analyse fine des phénomènes de convection observés sur les champs de vitesses et de températures tels que le bourgeonnement.

• Champ de vitesses

A titre indicatif, les champs de vitesses et de normes de vitesses sont présentées en

Annexe D, mais sont peu exploitables par rapport au nombre de Nusselt. Ce champ de vitesses sera par contre très intéressant dans le cas de l’étude de la modélisation d’un champ hydrothermal réel.

• Champ de nombres de Péclet

Figure 13 : Champ de nombres de Péclet pour un maillage de 128×128, Fo = 0,5.

* *

therm

U xPe U x Ra

k

∆= = ∆

Le nombre de Péclet traduit la dominance du transport par convection sur la diffusion lorsqu’il est supérieur à 1. S’il est inférieur à 1, il y a prédominance de la diffusion.

• Conclusion sur les résultats qualitatifs Pour un nombre de Rayleigh de 800, les résultats qualitatifs obtenus avec CAST3M sont

très satisfaisants et ce quels que soient les maillages et nombres de Fourier utilisés. Cependant, une analyse quantitative des résultats consistant à étudier les variations temporelles du nombre de Nusselt s’avère nécessaire.

Résultats quantitatifs Le Tableau 6 répertorie les résultats des calculs réalisés pour un nombre de Rayleigh de

800 et des maillages de 96×96 et de 128×128 mailles. Les fréquences d’oscillation sont déduites par transformation de Fourier de l’évolution temporelle du nombre de Nusselt. Pour chacun des schémas numériques étudiés, les dispersions des valeurs moyennes et des fréquences d’oscillations sont inférieures à 2% dans la gamme de Fourier de 0,25 à 0,5.

1,615 1,09 0,480 2,61 10-02

Page 39: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

39

Tableau 6 : Tableau récapitulatif des résultats de calculs obtenus avec CAST3M dans le cas d’un domaine fermé

Domaine fermé

Méthode de discrétisation

utilisée

Schéma temporel

Nombre de

Rayleigh

Maillage Nombre de Fourier

∆∆∆∆t* ∆∆∆∆t (ans) Valeurs moyennes du nombre de

Nusselt

Fréquences d’oscillation

du nombre de Nusselt

Condition initiale

96x96 0,4 2,17x10-5 2,58 9,01 273 CI1

128x128 0,4 2,17x10-5 2,58 9,08 282 CI1 Implicite 800

128x128 0,7 2,14x10-5 2,54 9,06 275 CI1

96x96 0,4 2,17x10-5 2,58 9,16 281 CI1 96x96 0,5 2,71x10-5 3,23 9,16 277 CI1

128x128 0,25 7,63x10-6 0,91 9,1 277 CI1 128x128 0,4 1,22x10-5 1,45 9,16 286 CI1 128x128 0,5 1,53x10-5 1,81 9,07 277 CI1

800

128x128 1 3,05x10-5 3,63 9,04 260 CI1

96x96 0,4 2,17x10-5 2,58 10,90 399 & 212 CIF800 128x128 0,4 1,22x10-5 1,45 8,81 - CI1 950

128x128 0,4 1,22x10-5 1,45 11,00 400 & 217 CIF800 96x96 0,4 2,17x10-5 2,58 13,6 595 & 230 CIF950

EMFH

Crank-Nicholson

1200 128x128 0,4 1,22x10-5 1,45 13,6 599 & 232 CIF950

64x64 0,4 4,88x10-5 5,81 9,11 262 CI1

96x96 0,4 2,17x10-5 2,58 9,13 279 CI1 800

128x128 0,25 7,63x10-6 3,63 9,1 279 CI1

96x96 0,4 2,17x10-5 2,58 8,89 - CI1

VF Crank-Nicholson

950 128x128 0,4 1,22x10-5 1,45 9,04 260 CI1

Page 40: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 40

Les valeurs moyennes et les fréquences obtenues avec un maillage de 128×128 mailles convergent donc dans cette gamme de pas de temps. Par contre, au-delà d’un nombre de Fourier de 0,65, la fréquence d’oscillation du nombre de Nusselt chute.

Conclusion

Pour la suite des calculs, seul le schéma de Crank-Nicholson associé à un nombre de Fourier

de 0,4 sera utilisé, car il répond aux différents critères de stabilité et de précision de la majorité des cas traités dans ce rapport. Une réserve toutefois est faite pour des nombres de Rayleigh dépassant une valeur de 1000 car les fréquences du Nusselt associées dépassent alors une valeur de 500. Les résultats restent cependant satisfaisants (cf. § 4.1.4).

Remarque : On obtient de bons résultats même pour un nombre de Fourier inférieur à 1/3.

o Résultats obtenus pour un Rayleigh de 950

Dès les premiers calculs, il a semblé que pour un nombre de Rayleigh de 950, les résultats

obtenus par calcul numérique dépendaient de la condition initiale imposée. Deux configurations initiales ont été utilisées, la condition initiale CI1, puis la condition initiale CIF800 correspondant au régime oscillant établi pour un nombre de Rayleigh égal à 800.

Résultats obtenus pour une condition initiale de type CI1

Pour cette perturbation initiale faible, identique à celle du cas Ra = 800, les champs de

températures et de vitesses obtenus sont des champs permanents (cf. figure 15). La figure 14 présente le champ de températures obtenu ne présentant plus un unique rouleau de convection mais trois rouleaux.

Le nombre de Nusselt, d’abord oscillant, devient stable après une durée de calcul égale à

environ un cinquième du temps caractéristique. La valeur moyenne du Nusselt est de moitié inférieure à celle attendue.

Figure 14 : Evolution du nombre de Nusselt en fonction du temps, pour Ra = 950 et une condition initiale de type CI1.

Figure 15 : Champ de températures permanent obtenu pour Ra = 950 et une condition initiale de type CI1.

14.00

12.00

10.00

8.00

6.00

4.00

2.00

0.00 0.00 0.50 1.00 1.50 2.00 2.50 3.00 Temps

Nusselt

Page 41: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 41

Résultats obtenus pour une condition initiale de type CIF800 Pour cette condition initiale, on parvient à maintenir une cellule de convection unique et le système reste périodique comme en témoignent les Figures 16 et 17.

Deux fréquences d’oscillation apparaissent, une fréquence principale f1 = 399 et une fréquence secondaire f2 = 212. La valeur moyenne du Nusselt est de 10,9.

o Résultats obtenus pour Ra = 1200. La condition initiale utilisée ici est une condition de type CIF950.

Figure 16 : Champ de températures obtenu pour Ra = 950 et une condition initiale de type CIF800.

Figure 17 : Evolution du nombre de Nusselt en fonction du temps, pour Ra = 950 et une condition initiale de type CIF800.

Figure 18 : Champ de températures obtenu pour Ra = 1200 et une condition initiale de type CIF950.

Figure 19 : Evolution temporelle du nombre de Nusselt, pour Ra = 1200 et une condition initiale de type CIF950.

0.00 0.50 1.00 1.50 2.00 2.50 3.00 Temps

12.00

10.00

8.00

6.00

4.00

2.00

0.00

Nusselt

Transition d’un nombre de Rayleigh de 800 à 950

478 °C 359 239 120

15.50

15.00

14.50

14.00

13.50

13.00

12.50

12.00

1.34 2.36 2.38 2.40 2.42 2.44 2.46 2.48 2.50 Temps

Nusselt

478 °C 359 239 120

Page 42: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 42

Deux fréquences d’oscillation existent pour ce cas, une fréquence principale élevée f1 = 595, et une fréquence secondaire f2 = 230. Ce jeu de données fournit donc des résultats très proches de ceux publiés par Fontaine et al. (2007) et paraît donc tout à fait valide pour le cas d’un domaine fermé.

o Conclusion sur la validité du code en boîte fermée pour les EFMH

En augmentant le nombre de Rayleigh au fur et à mesure, le système maintient un rouleau de convection unique. Le tableau 6 présente une comparaison entre les résultats obtenus avec CAST3M et les résultats de différents auteurs.

Domaine fermé Références

Ra 800 950 1200

Nu 9,1 10,9 13,6

f1 280 399 595 f2 212 230

Cette étude (CAST3M, EFMH)

Nu 9,2 10,9 13,06 f1 275 397 608 f2 200 228

Fontaine et al. (2007)

Nu 9,14 f1 299 f2

Cherkaoui et Wilcock (1999)

Nu 9,43 13,85 f1 296 695 f2 261 f3 69

Caltagirone et Fabrie (1989)

Tableau 7 : Comparaisons des résultats obtenus avec ceux publiés dans la littérature. Le modèle numérique utilisant la méthode des EFMH donne donc des résultats tout à fait

comparables aux résultats publiés dans la littérature, ce qui permet de valider le code. Remarque : Des calculs ont été réalisés avec un maillage de 64×64 mailles avec un schéma

de Crank-Nicholson. L’obtention de résultats exacts, stables en fonction du pas de temps, s’est faite avec de très faibles pas de temps correspondant à des nombres de Fourier allant de 0,02 à 0,1 (cf. Annexe E). Lorsque le pas de temps est abaissé en dessous de la valeur limite approximative du Fourier ≈ 0,02, la fréquence du nombre de Nusselt s’effondre en même temps que sa valeur moyenne et l’amplitude de ses oscillations. La condition présentée par Hoteit et al. (2002) semble donc fortement modifiée par la prise en compte de la convection dans le modèle numérique.

Page 43: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 43

4.1.3 Résultats obtenus à l’aide de la méthode des VF

Pour pouvoir comparer les deux méthodes, VF et EFMH, les maillages et les nombres de Fourier utilisés sont identiques à ceux utilisés précédemment.

o Résultats obtenus pour un nombre de Rayleigh de 800 Le calcul a été réalisé avec un maillage constitué de 96×96 puis de 128×128 mailles afin de vérifier la convergence spatiale des résultats. De la même façon, la convergence temporelle a été vérifiée en faisant varier le nombre de Fourier de 0,2 à 0,5 pour chaque type de maillage. On obtient un champ de températures identique à celui obtenu dans le cas de l’utilisation des EFMH comme en témoignent les Figures 20 et 21.

La fréquence d’oscillation du nombre de Nusselt obtenu est de 279 pour un maillage

constitué de 96×96 mailles et de 280 pour un maillage constitué de 128×128 mailles (cf. Tableau 6). Ces résultats sont donc très proches de ceux obtenus lors de l’utilisation des EFMH.

o Résultats obtenus pour un Rayleigh de 950

Résultats obtenus pour une condition initiale de type CI1

Pour cette perturbation initiale faible, identique à celle du cas Ra = 800, les champs de

températures et de vitesses obtenus sont des champs permanents, identiques à ceux obtenus dans le cas de l’utilisation de la méthode des EFMH. Le champ de températures ne présente plus un unique rouleau de convection mais trois rouleaux (cf. Figure 23).

Figure 20 : Champ de températures obtenu pour Ra = 800 et une condition initiale de type CI1.

Figure 21 : Evolution du nombre de Nusselt en fonction du temps, pour Ra = 800 et une condition initiale de type CI1.

11.00

10.00

9.00

8.00

7.00

6.00

5.00

4.00

3.00

2.00

1.00 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 Temps

Nusselt

478 °C 359 239 120

Page 44: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 44

Comme le montre la Figure 22, le nombre de Nusselt se stabilise à une valeur de 8,89 après

quelques oscillations. Cette valeur était de 8,81 dans le cas de l’utilisation de la méthode des EFMH. Le régime permanent est établi au bout d’un temps adimensionné de 0,30. On retrouve donc exactement le même comportement que lors du calcul utilisant la méthode des EFMH.

Résultats obtenus pour une condition initiale de type CIF800

Lorsque l’on prend pour condition initiale le champ de température obtenu pour un calcul utilisant un nombre de Rayleigh égal à 800, on parvient à maintenir un rouleau de convection unique (cf. Figure 24).

La Figure 25 laisse percevoir deux fréquences d’oscillation du nombre de Nusselt. La

transformation de Fourier permet d’en déduire ces deux valeurs de fréquence : f1 = 203 et f2 = 403.

Figure 22 : Champ de températures permanent obtenu pour Ra = 950 et une condition initiale de type CI1.

Figure 23 : Evolution du nombre de Nusselt en fonction du temps, pour Ra = 950 et une condition initiale de type CI1.

Figure 24 : Champ de températures permanent obtenu pour Ra = 950 et une condition initiale de type CIF800.

Figure 25 : Evolution du nombre de Nusselt en fonction du temps, pour Ra = 950 et une condition initiale de type CIF800.

12.00 11.60 11.20 10.80 10.40 10.00

0.95 1.00 1.05 1.10 1.15 1.20 Temps

Nusselt

16.00

14.00

12.00

10.00

8.00

6.00

4.00

2.00

0.00 0.00 0.20 0.40 0.60 0.80 Temps

Nusselt

478 °C 359 239 120

478 °C 359 239 120

Page 45: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 45

o Conclusion sur la validité du code en boîte fermée pour les VF

Les résultats obtenus à l’aide de la méthode des volumes finis sont très sensiblement identiques à ceux obtenus à l’aide de la méthode des éléments finis mixtes hybrides (cf. tableau 8).

Domaine fermé Références

Ra 800 950

Nu 9,1 10,9 f1 280 399 f2 212

CAST3M EFMH

Nu 9,2 11,1 f1 279 403 f2 203

CAST3M VF

Tableau 8 : Comparaisons des résultats obtenus sous CAST3M à l’aide de la méthode des EFMH et de la méthode des VF.

Une supposition couramment admise pour la validation des codes consiste à dire que si le comportement du système en boîte fermée est validé, alors étant donné qu’en domaine ouvert les équations à résoudre et les paramètres restent les mêmes, on peut valider de manière générale les résultats du code pour le système ouvert. Cependant, la disparité des résultats, même faible d’un auteur à l’autre dans le cas d’un domaine fermé, laisse penser qu’il est possible d’observer une amplification de la disparité des comportements des systèmes modélisés si la nature du domaine change. La validité du code dans le cas d’un domaine ouvert devrait donc être vérifiée.

4.2 Validation du jeu de données en domaine ouvert

Ce cas test utilise exactement le même jeu de données que le cas-test précédent. Seule la condition à la limite supérieure du domaine change. L’ensemble des résultats de calcul est présenté en Annexe F .

4.2.1 Conditions initiales et conditions aux limites

Figure 26: Schéma du domaine fermé et conditions aux limites adimensionnées.

T = 1, u = v = 0, L = 1

P= 2,25 Pa

u = v = 0, d = 1 u = v = 0 Milieu poreux homogène isotrope

∂T/∂z= 0 T= 0

Page 46: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 46

Dans ce cas, l’eau entre dans le domaine à une température de 0. La condition de sortie est une condition de gradient thermique nul, correspondant à un transfert thermique majoritairement convectif. Comme le montrent Cherkaoui et Wilcock (1999) (cf. Annexe G), la sensibilité aux conditions initiales dans le cas d’une boîte ouverte est plus importante qu’en boîte fermée. Comme dans le cas du domaine fermé, deux situations initiales pourront être utilisées (notées CIO1 pour un champ diffusif et CIO suivi du nombre de Rayleigh précédent pour une croissance progressive du nombre de Rayleigh).

4.2.2 Résultats des calculs

o Résultats obtenus à l’aide de la méthode des EFMH

Le maillage utilisé est un maillage de 96×96 mailles et un nombre de Fourier de 0,4. Dans un premier temps, pour chaque calcul à un nombre de Rayleigh donné, la condition initiale correspond au champ de températures obtenu à la fin du calcul précédent dont le nombre de Rayleigh est inférieur. Le résultat est une succession de champs de températures constitués de rouleaux de convection passant d’un rapport d’aspect de 1 à ½ autour d’un nombre de Rayleigh de 108. Ces champs de températures sont permanents jusqu’à une valeur du nombre de Rayleigh supérieure ou égale à 425 comme en témoigne la Figure 27. Certains champs de températures issus de ce calcul sont présentés en Annexes H, I, J, K.

Figure 27 : Evolution temporelle du nombre de Nusselt, obtenu à l’aide de la méthode des EFMH sous CAST3M, pour différents nombres de Rayleigh.

Le Tableau 9 compare le rapport d’aspect des rouleaux de convection et les caractéristiques

des nombres de Nusselt obtenus avec CAST3M à ceux issus de (Cherkaoui et Wilcock (1999)) et confirme l’absence de résultats aberrants.

12.00

10.00

8.00

6.00

4.00

2.00

0.00

0.00 1.00 2.00 3.00 4.00 5.00 6.00 7.00 Temps

Nusselt

Ra = 90 108 200 300 400 425 470

Page 47: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 47

Tableau 9 : Comparaison des résultats obtenus en boîte ouverte pour un nombre de Rayleigh croissant progressivement de 90, 108, 200, 300. *valeurs lues sur le graphique de l’article.

L’écart des valeurs moyennes et des fréquences entre les résultats des calculs de cette étude

et les résultats publiés par Cherkaoui et Wilcock (1999) est au maximum inférieur à 3%. A l’aide de cette méthode, un comportement apparemment similaire à celui décrit par Cherkaoui et Wilcock (1999) est obtenu. Il faut noter qu’ici, par manque de temps, la bifurcation présente autour d’un Rayleigh de 425 n’a pas été étudiée (voir le cas traité à l’aide des VF). Pour un nombre de Rayleigh d’une valeur d’environ 425, Cherkaoui et Wilcock (1999) observent la possibilité d’obtenir par le calcul un régime oscillant ou permanent. Malgré la forte croissance du nombre de Rayleigh, de 425 à 470, un régime oscillant identique à celui observé par Cherkaoui et Wilcock (1999) est obtenu pour un nombre de Rayleigh de 470.

Pour cette valeur, Ra = 470, au niveau de la limite basse du champ de température modélisé,

l’existence de faibles instabilités se traduit par l’apparition périodique de bourgeonnements de fluides chauds, disposés symétriquement par rapport à l’axe médian du domaine :

Figure 28: Champ de températures obtenu pour Ra = 470 et une condition initiale de type CIO425.

Une remarque importante toutefois : pour un nombre de Rayleigh inférieur à une valeur de 400, si la condition initiale utilisée correspond à un champ de température obtenu par la seule

Domaine ouvert Références Ra 90 108 200 300 400 425 470

Nu 3,21 3,6 6,49 8,16 9,16 9,38 9,85

f1 - - - - - - 121 Amplitude des

oscillations - - - - - - 1,82

rapport d'aspect A 1 1/2 1/2 1/2 1/2 1/2 1/2

Cette étude

Nu 3,20* 3,50* 6,40* 8,2* 9,10* 9,30* 9,80*

f1 - - - - - - 117,2 Amplitude des

oscillations - - - - - - 1,18*

rapport d'aspect A

1 1/2 1/2 1/2 1/2 1/2 1/2

Cherkaoui et Wilcock (1999)

478 °C 359 239 120

Page 48: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 48

diffusion de la chaleur en régime permanent (condition de type CIO1), on obtient exactement les mêmes résultats.

Domaine ouvert Références

Ra 90 108 200 300 400

Nu 3,21 3,6 6,49 8,16 9,16 f1 - - - - -

rapport d'aspect A

1 1/2 1/2 1/2 1/2 Cette étude

Tableau 10 : Résultats obtenus en boîte ouverte pour des nombres de Rayleigh de 90, 108, 200, 300 et 400 avec une condition initiale de type CIO1.

Or, pour un nombre de Rayleigh inférieur à 400, ces systèmes, solutions des calculs effectués

avec CAST3M, correspondent à ceux considérés comme répondant à l’hypothèse de Malkus (1954)* par Cherkaoui et Wilcock (1999). Autrement dit, les codes développés dans cette étude, offrant des résultats a priori moins sensibles à la condition initiale pour de faibles nombres de Rayleigh, convergent vers une solution physiquement stable. Ce résultat est discuté au paragraphe 3.6.2.

*D’après Malkus (1954), les systèmes convectifs adoptent la configuration qui maximise le transport de la chaleur.

o Résultats obtenus à l’aide de la méthode des VF De la même manière que lors de l’utilisation de la méthode des EFMH, deux types de

simulations ont été effectués. Le premier voit croître progressivement le nombre de Rayleigh :

Figure 29 : Evolution temporelle du nombre de Nusselt, obtenu à l’aide de la méthode des VF avec CAST3M, en faisant croître progressivement le nombre de Rayleigh.

Ra = 400 425 430 450 470

11.00

10.50

10.00

9.50

9.00

8.50

8.00 0.00 0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.00 Temps

Nusselt

Page 49: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 49

Domaine ouvert Références Ra 400 425 430 450 470

Nu 9 ,16 9,35 9,39 9,53 9,87 f1 - - - - 116,1

Amplitude des oscillations

- - - - 1,53*

rapport d'aspect A

1/2 1/2 1/2 1/2 1/2

Cette étude

Nu 9,10* 9,30* 9,42* 9,60* 9,80*

f1 - - - - 117,2 Amplitude des

oscillations - - - - 1,18*

rapport d'aspect A

1 1/2 1/2 1/2 1/2

Cherkaoui et Wilcock (1999)

Tableau 11 : Comparaison des résultats obtenus en boîte ouverte pour un nombre de Rayleigh croissant progressivement de 400, 425, 430, 450 et 470.

*valeurs lues sur le graphique de l’article.

Les valeurs moyennes du nombre de Nusselt sont donc sensiblement identiques à celles obtenues par Cherkaoui et Wilcock (1999) (cf. Annexe H). Cependant, malgré l’apparition d’un régime périodique pour un nombre de Rayleigh de 470, la bifurcation décelée par Cherkaoui et Wilcock (1999) autour d’un nombre de Rayleigh de 425 n’a pas été observée.

Enfin, lors du calcul pour un nombre de Rayleigh de 470, que ce soit avec la méthode des VF

ou celle des EFMH, l’amplitude des oscillations est supérieure de 30 à 50% à celle calculée par Cherkaoui et Wilcock (1999). Ces résultats sont discutés aux paragraphes 4.3.1 et 4.3.2.

Le deuxième type de simulation numérique est un calcul effectué avec une condition initiale

correspondant à un champ de température obtenu par diffusion seule en régime permanent, condition notée ici CIO1, les différents résultats sont présentés dans le tableau suivant :

Domaine ouvert Références

Ra 400 425 430 450 470

Nu 9 ,16 9,35 9,39 9,53 9,87

f1 - - - - 116,1 Amplitude des

oscillations - - - - 1,53*

rapport d'aspect A

1/2 1/2 1/2 1/2 1/2

Cette étude

Tableau 12 : Résultats obtenus en boîte ouverte pour des nombres de Rayleigh de 400, 425, 430, 450 et 470 avec une condition initiale de type CI1.

Avec une telle condition initiale, pour un nombre de Rayleigh inférieur ou égal à une valeur

de 470, les résultats sont exactement les mêmes que ceux obtenus lors de la croissance progressive du nombre de Rayleigh. L’étape de validation pour le cas ouvert est ici incomplète puisque seuls de faibles nombres de Rayleigh ont été testés et que nous n’obtenons pas exactement le comportement décrit par d’autres publications. De plus, l’unique expérimentation dont les résultats ont été publiés est celle de Cherkaoui et Wilcock (2001). Or les cellules de Hele-Shaw utilisées ne possèdent pas de parois adiabatiques et les résultats sont assez différents de ceux obtenus par les calculs numériques.

Page 50: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 50

4.3 Discussion

4.3.1 Approximation de Boussinesq, linéarité de l’équation de variation de la densité et viscosité constante

Comme le rappelle Zeytounian (2002), l'approximation de Boussinesq a été formulée par

Joseph Boussinesq (1903) dans son traité « Théorie Analytique de la chaleur » : « […] dans la plupart des mouvements provoqués par la chaleur sur nos fluides pesants, les

volumes ou les densités se conservent à très peu près, quoique la variation correspondante du poids de l’unité de volume soit justement la cause des phénomènes qu’il s’agit d’analyser.

De là résulte la possibilité de négliger les variations de la densité, là où elles ne sont pas multipliées par la gravité g, tout en conservant dans les calculs, leur produit par celle-ci. »

Cette approximation, couplée à l’approximation faite sur la variation de porosité en fonction du temps (∂φ/∂t =0), permet de s’affranchir de la définition de la porosité, de sa variation et des conséquences sur les perméabilités rencontrées. C’est une approximation forte, d’autant plus que les dimensions du système et l’existence sous certaines conditions d’un régime périodique ainsi que les pas de temps de calcul utilisés – à chaque pas de temps, une évolution du système est bien visible - ne permettent pas de dire que la variation temporelle de la masse volumique est négligeable. D’ailleurs, comme le relèvent Coumou et al. (2006), les simulations du système type « boîte ouverte », avec les différentes approximations que sont l’approximation de Boussinesq, la linéarisation de la densité et une viscosité constante donnent des résultats assez différents de ceux de modèles plus complexes.

Cependant, les approximations générales sur le comportement du fluide et de la roche sont

légion à la vue de la complexité du système. De plus, les résultats numériques des calculs en boîte fermée avec lesquels les résultats CAST3M ont été comparés utilisent souvent cette approximation (Graham et Steen 1992, Steen et Aidun 1988, Caltagirone et Fabrie 1989). Aussi, même si cette approximation n’est sûrement pas vérifiée dans le cas des systèmes hydrothermaux naturels, elle a permis lors de ce stage de pouvoir répondre aux interrogations posées.

4.3.2 Linéarisation et précision du calcul

L’algorithme de résolution, présenté au paragraphe 3.4, montre que la masse volumique, fonction de la température, est utilisée comme variable de couplage entre l’équation de Darcy et l’équation de transfert de la chaleur. Or, dans cette étude, lors du calcul de la vitesse de Darcy au temps n+1, on utilise une valeur de la masse volumique calculée au temps n. Par la suite, cette vitesse de Darcy est injectée dans le terme de convection de l’équation du transfert de la chaleur afin de calculer la masse volumique au temps n+1. Si rien n’est fait, il y a donc une discordance des valeurs au sein de la boucle en temps de l’algorithme. Si le pas de temps de calcul est trop grand, cela peut entraîner une imprécision de calcul. Le manque de précision décelé dans cette étude lors de l’utilisation du schéma implicite en temps pour l’équation de transfert de la chaleur en est probablement une des conséquences. L’augmentation de l’ordre en temps à l’aide de l’utilisation d’un schéma temporel de Crank-Nicholson semble pallier à ce défaut, en tout cas pour l’étude des domaines fermés, même à haut nombre de Rayleigh. En fait, le critère de 80 pas de temps par demi période d’oscillation du nombre de Nusselt conduit à l’utilisation de pas de temps petits par rapport à l’évolution de la masse volumique - la

Page 51: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 51

masse volumique étant une fonction linéaire de la température -. L’utilisation de tels pas de temps limite donc la discordance des valeurs de la masse volumique au temps tn et tn+1. La non linéarité est liée à la présence dans le terme de convection de l’équation du transfert de la chaleur d’une fonction linéaire de T - la vitesse de Darcy (U) -, multipliée par une fonction de T - le gradient spatial de T -. La discordance présentée auparavant entre l’utilisation de ρn et le calcul à partir de Tn+1 de ρ

n+1 est donc amplifiée par l’existence de cette non linéarité. Des erreurs s’ajoutent donc à un ordre supérieur ou égal à 2, d’un pas en temps de calcul à l’autre lors de l’utilisation notre méthode de linéarisation. Un problème de cinétique du système peut donc apparaître à partir d’une certaine valeur du pas de temps et conduire à des solutions aberrantes. Les méthodes de Picard ou de Newton-Raphson permettent normalement d’obtenir des résultats plus précis. Ces méthodes seront utilisées lors de la thèse qui fera la suite à ce stage.

4.3.3 Importance de la condition initiale

Le calcul, à l’aide des EFMH comme des VF, pour un nombre de Rayleigh de 950 montre

toute l’importance de la condition initiale dans les calculs réalisés à haute valeur du nombre de Rayleigh. Peut-on pour autant parler de régime chaotique ? Dans les cas des domaines clos étudiés, non. L’une des conditions nécessaires pour définir un régime chaotique est l’imprédictibilité du système. Or, pour un calcul avec un nombre de Rayleigh de 950, lorsque l’on impose une condition initiale correspondant au régime obtenu pour un Rayleigh de 800 à un temps t ou à t+T/2, T étant la période d’oscillation du régime périodique obtenu pour Ra = 800, les solutions obtenues sont identiques.

Lors du calcul pour un nombre de Rayleigh de 950 avec une condition initiale correspondant

à une faible perturbation du champ de températures, on remarque que dans un premier temps le système oscille autour d’une valeur du nombre de Nusselt plus faible que la valeur attendue, puis se stabilise et prend pour valeur 6,3 (cf. § 4.1.2 et 4.1.2).

Le chemin parcouru jusqu’à l’obtention d’un régime permanent lors du calcul pour un

nombre de Rayleigh de 950 peut s’expliquer de trois façons non-exclusives les unes par rapport aux autres. La première peut être l’existence d’un état d’équilibre thermodynamique intermédiaire, sorte « d’attracteur parasite » autour duquel, étant donné le chemin suivi naturellement par le système lors d’une telle évolution, ce même système vient s’agglomérer. La deuxième raison possible est numérique : il est possible que dans certains cas les erreurs de troncature liées au calcul numérique, ou la présence d’oscillations, conduisent à cet état permanent. Enfin la méthode utilisée pour traiter le caractère non linéaire du système d’équation peut conduire à un problème de cinétique du système modélisé (cf. § 4.1.2 et 4.3.2).

Cependant, Fontaine et al. (2007) ou encore Caltagirone et Fabrie (1989), constatent, eux

aussi, que pour maintenir un rouleau unique de convection en domaine fermé, pour de grands nombres de Rayleigh, il est nécessaire de faire croître progressivement ce nombre. Ceci n’accrédite pas systématiquement, mais du moins soulage, le choix de la méthode de linéarisation utilisée dans cette étude.

En domaine ouvert, Cherkaoui et Wilcock (1999) ont mis en évidence que le calcul pour un

nombre de Rayleigh de 550 avec une condition initiale correspondant au régime trouvé pour n’importe quel nombre de Rayleigh inférieur à 500, conduit à un régime permanent à quatre

Page 52: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 52

rouleaux de convection. Ce même régime est obtenu pour un Rayleigh de 300, 450 et 480 si la condition initiale imposée est le régime obtenu pour un Rayleigh de 700.

Dans le cas de la modélisation d’un système hydrothermal, une telle situation devrait donc

impliquer des comparaisons de résultats provenant de calculs effectués à partir de différentes conditions initiales réalistes.

Cependant, les calculs effectués dans notre étude, que ce soit à l’aide de la méthode des

EFMH ou des VF, n’exhibent pas de sensibilité aux conditions initiales en dessous d’un nombre de Rayleigh d’au moins 400 pour les EFMH et 470 pour les VF. De plus, les deux méthodes conduisent à des résultats tout à fait comparables et proches des résultats considérés comme physiquement stables par Cherkaoui et Wilcock (1999) (cf. § 4.2.2).

Donc, afin de lever l’ambiguïté liée à notre méthode de linéarisation, il est bien sûr

nécessaire de comparer les calculs obtenus avec d’autres méthodes de linéarisation. Des calculs utilisant des pas de temps très petits, de l’ordre de 5 à 40 fois plus petits que ceux utilisés dans cette étude, couplés à l’utilisation de notre méthode de linéarisation, peuvent apporter une aide précieuse. En effet, bien que dans la pratique de tels pas de temps ne soient pas utilisables de façon répétitive, cette levée simple du problème de discordance (cf. § 4.3.2) permettrait de valider la méthode de linéarisation finalement choisie.

La méthode de linéarisation choisie pour l’étude présentée dans ce mémoire, bien qu’a priori

moins précise, retrouverait alors toute son utilité.

4.4 Conclusions et perspectives de modélisation d’un système hydrothermal réel

4.4.1 Conclusion sur l’utilisation des EFMH et des VF Le temps de calculs, lorsqu’on utilise les EFMH, est assez conséquent : pour un maillage de

128×128 mailles et un nombre de Fourier de 0.4, un calcul avec CAST3M d’une durée d’un temps caractéristique prend environ 6 jours. Pour un maillage de 256×256 le temps de calcul devient supérieur à la vingtaine de jours, ce qui est en pratique très contraignant. On pourrait procéder à un relâchement du maillage dans la partie centrale du domaine où les instabilités sont étouffées par la diffusion dans le cas d’une boîte fermée, comme en témoigne le nombre de Péclet, mais la condition de stabilité présentée par Hoteit et al. (2002) limite immédiatement ce relâchement.

Après discussion avec d’autres chercheurs, il semble qu’à précision équivalente, la méthode

des volumes finis soit aussi coûteuse en temps. Cependant dans notre cas, le préconditionnement des matrices utilisées dans le cas des volumes finis, a quant à lui permis de réduire fortement le temps de calcul de près de 35%.

L’avantage qu’apporte la méthode des volumes finis est l’absence de valeur minimale pour

le pas de temps. Ce dernier avantage n’est pas négligeable : l’objectif étant de décrire le comportement du système pour des nombres de Rayleigh dépassant parfois une valeur de 2000, la multiplication des fréquences et l’élévation de leurs valeurs exigent des pas de temps très petits.

Finalement, les jeux de données développés lors de ce stage, malgré leur relative simplicité, donnent de bons résultats, mais la méthode de linéarisation doit absolument être comparée aux méthodes plus complexes, du type Newton-Raphson par exemple. L’objectif à présent est d’explorer le plus de configurations possibles en faisant varier le nombre de Rayleigh afin d’identifier

Page 53: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 53

d’éventuels problèmes dans le but d’enrichir ensuite le modèle et de pouvoir à terme modéliser un champ hydrothermal réel.

Les contraintes imposées sur les pas de temps et le fait que les VF soient une méthode

actuellement toujours développée par certains numériciens du CEA, ont amené à préférer la méthode des volumes finis utilisées pour la résolution de l’équation du transfert de la chaleur au détriment de la méthode des éléments finis mixtes hybrides. Celle-ci sera généralisée lors de la thèse à l’ensemble des équations du jeu de données dès lors que l’approximation de Boussinesq ne pourra plus être supposée valable. De même, la méthode de linéarisation doit très certainement être améliorée.

4.4.2 Perspectives de modélisation : le champ hydrothermal réel

Pour le projet MANTHY dans lequel s’inscrit ce stage, la modélisation d’un champ

hydrothermal doit être un support à la compréhension des mécanismes conduisant à la formation de l’hydrogène. Plusieurs étapes sont alors indispensables :

Le choix du site modélisé : pour ce mémoire, nous avons pris comme exemple le site de

Rainbow qui est certainement le site le plus documenté à l’heure actuelle. Mais justement, sa notoriété risque bien d’entraîner le projet vers une redite de publication.

La description la plus précise possible de l’hydraulique du système à partir des données de

terrain.

L’écriture d’un jeu de données pour la partie hydraulique et la validation de ce jeu, phase dans laquelle ce stage s’est inscrit. Cette phase n’est pas terminée. A présent, il est important d’étudier le comportement du modèle pour de grands nombres de Rayleigh dans le cas « ouvert » et de complexifier ensuite le modèle en affinant, par exemple, l’algorithme de résolution afin de linéariser au mieux le système et d’ainsi, peut-être augmenter la précision des résultats.

La modélisation du site et la détermination des différents scénarii hydrauliques qui peuvent

conduire aux données hydrauliques de terrain. Cela implique de bien choisir le type de conditions aux limites notamment sur la frontière inférieure du modèle. En effet, sur cette frontière, il pourrait être plus judicieux d’imposer un flux thermique plutôt qu’une température, car si l’on suppose la conservation de l’énergie vérifiée dans le plan modélisé, alors les flux sortants (mesurés sur le terrain) doivent être égaux aux flux entrants. Cette hypothèse impose toutefois de connaître la profondeur limite de l’écoulement hydrothermal. Ce choix sera fait en fonction des connaissances acquises du terrain. De même, pour pallier à la méconnaissance de l’agencement rocheux et donc des perméabilités, on sera peut-être amené à utiliser des méthodes inverses en combinaison avec un réseau de neurones. Il faudra de même prendre en compte le fait que l’on modélise en deux dimensions le site.

La modélisation géochimique et l’étude des rétroactions possibles de cette dernière sur la

partie hydraulique du système. Les expérimentations en laboratoire prévues par le projet MANTHY devraient apporter de nouvelles connaissances à ce sujet.

Ces travaux sont le sujet d’une thèse proposée par le CEA et l’IFREMER.

Page 54: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 54

BILAN DU STAGE ET CONCLUSION Les objectifs de ce stage étaient de comprendre le fonctionnement des sites hydrothermaux, d’amorcer un modèle numérique simple et de valider ce modèle. Cette étude montre que le modèle numérique développé lors de ce stage est une base solide pour le développement d’un modèle plus complexe. De plus, elle a permis de mettre en évidence une bonne précision des résultats malgré la méthode de linéarisation utilisée et les quelques points à améliorer (cf. § 4.3). La méthodologie de validation du code ayant été mise en place et appliquée, la validation d’un code plus complexe n’en sera que plus rapide.

Bien sûr, dans le cadre de la modélisation hydrothermale, le travail mené lors de ce stage n’a pas apporté de nouvelles connaissances scientifiques. Le modèle développé et les schémas numériques sont restés simples mais la somme des problèmes déjà rencontrés et l’étendue des domaines concernés ont été une grande source d’enrichissement, d’autant plus qu’une grande liberté d’action m’a été laissée pendant ce stage et m’a permis d’apprendre à avoir une démarche scientifique pour répondre aux problèmes rencontrés. Les recherches bibliographiques se sont portées tant sur les problèmes numériques que sur les problèmes de transition vers le chaos en milieu poreux, sur les types de modèles développés pour la circulation hydrothermale, les données de terrain et leurs analyses par de multiples chercheurs.

Les apports aux laboratoires, issus d’un tel stage, sont multiples : d’une part l’application de

méthodes en court de développement permet aux numériciens d’avoir un « retour sur investissement ». Par exemple, l’application à un cas concret de la méthode des volumes finis, telle qu’elle est développée au CEA, a permis de montrer son efficacité et sa précision malgré la simplicité des schémas utilisés. Cette application a d’ailleurs mis en évidence la nécessité de développer de nouveaux outils pour le traitement de l’équation de continuité en régime transitoire. D’autre part, la somme des connaissances recueillies aux sujets des systèmes hydrothermaux, tant sur un plan numérique qu’expérimental, les références bibliographiques associées, la méthodologie développée et les problèmes soulevés dans ce stage, présentés dans ce rapport, permettent à quiconque voudrait poursuivre ce travail, d’acquérir une avance considérable. Les problèmes rencontrés ont surtout été d’ordre numérique. Je pensais d’ailleurs, au début de ce stage, que les problèmes de stabilité et de précision des schémas étaient mieux maîtrisés, qu’aux problèmes rencontrés correspondaient un schéma « parfait », bien connu. Mais les nombreuses discussions avec des numériciens ou modélisateurs et la multitude de méthodes présentées dans les publications traitant de la modélisation des sites hydrothermaux (volumes finis couplés aux éléments finis, différences finies « améliorées», codes pseudo-spectraux) m’ont fait comprendre que bien souvent le choix d’une méthode ou d’un schéma découle d’un historique d’utilisation, les schémas à disposition étant très souvent modifiés et enrichis en fonction des applications.

Autour de ce sujet, des réunions très enrichissantes ont été organisées notamment avec Jean-Luc Charlou (IFREMER), Philippe Jean-Baptiste (CEA-LSCE) et Claude Mügler (CEA-LSCE), pour rendre compte de l’avancée du travail et du sujet de thèse qui normalement suit ce stage. Les réunions d’équipe, quant à elles, ont touché beaucoup de sujets de modélisation assez différents, tels que la modélisation du ruissellement de surface ou du transfert thermique dans le sol, l’équipe ayant un large panel de sujets auxquels la modélisation peut apporter une réponse. Le gain personnel issu d’un tel stage est donc indéniable.

Page 55: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 55

La source première de ma motivation lors de ce stage a incontestablement été le sujet de la thèse qui débutera en octobre. Tout le problème d’un tel sujet peut être résumé dans l’expression « exploration numérique de domaines physiques peu connus ». En effet, bien que les méthodes numériques voient leurs performances augmenter significativement, l’apport expérimental et les données de terrain sont indispensables à la compréhension plus ou moins précise de tels systèmes. C’est la force de cette thèse qui attachée au projet MANTHY pourra bénéficier de l’apport des connaissances de différents laboratoires de pointe dans le domaine expérimental.

Page 56: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 56

BIBLIOGRAPHIE

Ouvrages, publications et rapports internes : Apps, J.A., (1985). Methane formation during hydrolysis by mafic rocks. Univ. of Calif., Lawrence Berkeley Lab., Berkeley, Calif., annual report, pp. 13–17. Abrajano , T.A., Sturchio, N.C., Bohlke, J.K., Lyon, G.L., Poreda, R.J., Stevens, C.M., (1988). Methane– hydrogen gas seeps, Zambales Ophiolite, Philippines: deep or shallow origin? Chem. Geol. 71, 211 –222. Abrajano , T.A., Sturchio, N.C., Kennedy, B.M., Lyon, G.L., Muehlenbachs, K., Bohlke, J.K., (1990). Geochemistry of reduced gas related to serpentinization of the Zambales ophiolite, Philippines. Appl. Geochem. 5, 625– 630. Allen, D.E., Berndt, M.E., Seyfried, W.E., (1996). Reduction of CO2 during serpentinization of olivine at 300 jC and 500 bar. Geology 24 (4), 351–354. Bergé, P., Y. Pomeau, Ch. Vidal , L’ORDRE dans le CHAOS, Vers une approche déterministe de la turbulence.-Paris : Hermann, 1984. Bernard-Michel, G., Genty, A., (2006). Modules d’écoulement-modèle de Richards – et de transport en milieux poreux non saturés dans CAST3M. Rapport interne CEA-DM2S 11. Berndt, M.E., Allen, D.E., Seyfried Jr., W.E., (1996). Reduction of CO2 during serpentinisation of olivine at 300°C and 500 bar. Geology 24, 351-354. Bischoff, J.L., Rosenbauer, R.J., (1987). Phase separation in seafloor geothermal systems: an experimental study of the effects on metal transport. Am. J. Sci. 287, 953–978. Bischoff, J.L., Rosenbauer, R.J., (1989). Salinity variations in submarine hydrothermal systems by layered double-diffusive convection. J. Geol. 97, 613–623. Boussinesq J., Théorie Analytique de la chaleur, Vol. II Gauthier-Villars, Paris, 1903. Caltagirone J.P., P. Fabrie, (1989), Natural in a porous medium at high rayleigh numbers, Part I – Darcy’s model, Eur. J. Mech., B/Fluids, Vol 8, n°3, 207-227. Charlou, J.L., Donval, J.P., Douville, D., Knoery, J., Fouquet, Y., Jean-Baptiste, P., Bourg, C., Prieur, D., German, C., (1998). Fluides hydrothermaux en contexte de roches ultrabasiques : exemples de Logatchev (14j45VN) et Rainbow (36j14VN) sur la dorsale médio-Atlantique. Réunion des Sciences de la Terre RST 98, Brest, 31 Mars-3 Avril 1998, p. 89. Charlou, J.L., Donval J.P., Fouquet Y., Jean-Baptiste P., Holm N., (2002), Geochemistry of high H2 and CH4 vent fluids issuing from ultramafic rocks at the Rainbow hydrothermal field (36°14’N, MAR), Chem. Geol. 191, 345-359. Charlou, J.L, (2009). Communication personnelle.

Page 57: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 57

Cherkaoui, A.S.M.,. Wilcock, W.S.D., (1999). Characteristics of high Rayleigh number two-dimensional convection in an open-top porous layer heated from below, J. Fluid Mech., 394, 241–260. Cherkaoui A.S.M., Wilcock W.S.D., (2001), Laboratory studies of high Rayleigh number circulation in an open-top Hele-Shaw cell : An analog to mid-ocean ridge hydrothermal systems, J. Geophys. Res., Vol. 106, n°B6, 10,983-11,000. Conveney, R.M., Goebel, E.D.,Zeller, E.J., Dreschoff, G.A.M., Angino, E.E., (1987) Serpentinization and the origin of hydrogen gas in Kansas. AAPG Bull. 71, 39-48. Coumou D., Driesner T., Geiger S., Heinrich C.A., Matthäi S., (2006), The dynamics of mid-ocean ridge hydrothermal systems: splitting plumes and fluctuating vent temperatures, Earth and Planet. Sci. Let. 245, 218-231. Dabbene F., (1993), Résolution des équations de Darcy par une méthode d’éléments finis mixtes hybrides. Rapport CEA-DMT 637. Dabbene F., X. Nouvellon, (1994), Extension de la méthode des éléments finis mixtes hybrides à la résolution des équations de Darcy en régime transitoire. Rapport CEA-DMT 566. Dabbene F., (1995), Schéma de diffusion-convection en éléments finis mixtes hybrides. Rapport CEA-DMT 613. De Marsily, G., (2004), Cours d’Hydrogéologie, Université Paris VI. Douville E., Charlou J.L., Oelkers E.H., Bienvenu P., Colon C.F.J., Donval J.P., Fouquet Y., Prieur D., Appriou P., (2002), The rainbow vent fluid : the influence of ultramafic rocks and phase seperation on trace metal content in Mid-Atlantic Ridge hydrothermal fluids, Geol. 184, 37-48.

Fontaine F.J., Wilcock W.S.D., Butterfield D.A., (2007), Physical controls on the salinity of mid-ocean ridge hydrothermal vent fluids, Earth and Planet. Sci. Let. 257, 132-145. German, C.R. et al., (1995). RRS Charles Darwin Cruise 89, 19 Aug-13 Sep 1994. Hydrothermal exploration at the Azores Tiple junction : HEAT, Institue of Oceanographic Sciences Deacon Laboratory, Cruise Report, 246, 58pp. Graham, M.D., P.H. Steen, (1992), Strongly interacting travelling waves and quasiperiodic dynamics in porous medium convection, Physica D 54, 331-350. Graham M.D., P.H. Steen, (1994), Plume formation and resonant bifurcations in porous-media convection, J. Fluid Mech., Vol. 272, 67-89. Hoteit. H., R. Mosé, B. Philippe, Ph. Ackerer, J. Erhel, (2002), The maximum principle violations of the mixed-hybrid finite-element method applied to diffusion equations, Int. J. Numer. Meth. Engng, 55, 1373-1390. Jean-Baptiste, P., Fourre, E., Charlou, J.L., German, C.R., Radford-Knoery, J., (2004). Helium isotopes at the Rainbow hydrothermal site (Mid-AtlanticRidge, 6°14′N). Earth Planet. Sci. Lett. 221, 325–335.

Page 58: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 58

Jean-Baptiste, P., E. Fourre´, A. Dapoigny, J. L. Charlou, and J.-P. Donval (2008), Deepwater mantle 3He plumes over the northern Mid-Atlantic Ridge (36_N–40_N) and the Azores Platform, Geochem. Geophys. Geosyst., 9, Q03010, doi:10.1029/2007GC001765. Kawada Y., Yoshida S., Watanabe S., (2004). Numerical simulations of mid-ocean ridge hydrothermal circulation including the phase separation of seawater, earth planets space, 56, 193-215. Kelley, D.S., Gillis, K.M., (1993). Fluid evolution in submarine magma-hydrothermal systems at the Mid-Atlantic Ridge. J. Geophys. Res. 98, 19579–19596. Kelley, D.S., (1997). Fluid evolution in slow-spreading environments. Proc. Ocean Drill. Program Sci. Results 153, 309–415. Kimura , S., G. Schubert, J.M. Straus, (1986), Route to chaos in porous-medium thermal convection, J. Fluid. Space, 56, 193-324. Landau, L., E. Lifchitz , Physique Théorique, Mécanique des fluides, Vol.6. - Deuxième édition. Moscou : MIR, 1989. Lewis, K.C., Lowell, R.P., (2004). Mathematical modelling of phase separation of seawater near an igneous dike, Geofluids, 4, 198 – 207. Malkus, W. V. R. (1954) The heat transport and spectrum of thermal turbulence. Proc. R. Soc. Lond. A 225, 196-212. Maugis P., (2006), Transferts complexes en milieu poreux : Quelques approaches physiques et numériques, Thèse de Doctorat, Université Paris 6, Paris. Mével C., (2003). Serpentinization of abyssal peridotites at mid-ocean ridges, Comptes-Rendus Geoscience 335 825-852. Moody, J.B., (1976). Serpentinization: a review. Lithos 9, 125–138. Neal, C., Stanger, G., 1983. Hydrogen generation from mantle source rocks in Oman. Earth Planet. Sci. Lett. 66, 315– 320. Paillere H., F.Dabbene, (1997), Initiation à la simulation numérique en mécanique des fluides à l’aide de CAST3M2000, recueil d’exemples commentés. Cours ENSTA MF307, Rapport DMT 308. Palmason, G., (1967), On heat flow in Iceland in relation to the mid-Atlantic ridge, Iceland and Mid-ocean Ridges, Publ. 38, Soc. Sci. Island., Reykjavik, 111-127. Ribando, R.J., K.E. Torrance, D.L.Turcotte, (1976), Numerical Models for Hydrothermal Circulation in the Oceanic crust, J. Geophys. Res. Vol.81, n°17, 3007-3012. Schoofs, S., Hansen, U., (2000). Depletion of a brine layer at the base of ridge-crest hydrothermal systems. Earth Planet. Sci. Lett. 180, 341– 353. Seyfried, W.E., Berndt, M.E., Janecky, D.R., (1986). Chloride depletions and enrichments in seafloor hydrothermal fluids: constraints from experimental basalt alteration studies. geochim. Cosmochim. Acta. 50, 469– 475.

Page 59: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 59

Singh S.C., Wayne C.C., H. Carton, T. Seher, V. Combier, M. Cannat, J.P. Canales, D. Düsünür, J. Escartin, J.M. Miranda, (2008), Discovery of a magma chamber and faults beneath a Mid-Atlantic Ridge hydrothermal field, Nature, Vol. 442, 1029-1032. Sourirajan , S., Kennedy, G.C., (1962). The system H2O–NaCl at elevated temperatures and pressures. Am. J. Sci. 260, 115– 141. Steen P.H., C.K. Aidun, (1988), Time-periodic convection in porous media : transition mechanism, J. Fluid Mech. Vol. 196, 263-290. Von Damm, K.L., (1988). Systematics of and postulated controls on submarine hydrothermal solution chemistry. J. Geophys. Res. 93, 4551– 4561. Von Damm, K.L., Bray, A.M., Buttermore, L.G., Oosting, S.E., (1998). The geochemical controls on vent fluids from the Lucky Strike vent field, Mid-Atlantic Ridge. Earth Planet. Sci. Lett. 160, 521– 536. Wetzel, L.R., Shock, E.L., (2000). Distinguishing ultramafic-from basalt-hosted submarine hydrothermal systems by comparing calculated vent fluid compositions. J. Geophys. Res. 105, 8319– 8340. Zeytounian R. Kh., (2002), Joseph Boussinesq and his approximation : a contemporary view, C.R. Mecanique 331, 575 – 586. Site internet : Stéphane Labrosse, Planète terre, ENS Lyon : http://planet-terre.ens-lyon.fr/planetterre/XML/db/planetterre/metadata/LOM-geophysique-fonds-oceaniques.xml consulté le 21 juin 2009.

Page 60: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 60

ANNEXES

Tables des annexes

Annexe A : Naissance et transport d’une vague de fluide chaud.................................................................. 61

Annexe B : Compléments sur l’analyse de la convection............................................................................... 62

Annexe C : Exemple de champ de pressions obtenu pour le cas Ra = 800 .................................................. 62

Annexe D : Champ de vitesses et normes de vitesses ..................................................................................... 63

Annexe E : Résultats des simulations numériques effectuées pour un nombre de Rayleigh de 800, avec un maillage de 64x64 mailles et des nombres de Fourier inférieurs à 1/3. .......................................... 64

Annexe F : Tableau récapitulatif des résultats de calculs obtenus avec CAST3M dans le cas d’un domaine ouvert .................................................................................................................................................. 65

Annexe G : Graphique issue de (Cherkaoui et Wilcock, 1999)..................................................................... 66

Annexe H : Résultats obtenus pour une nombre de Rayleigh de 90, une condition initiale de type CI1, en domaine ouvert, avec un maillage de 96×96 mailles.................................................................................. 67

Annexe I : Résultats obtenus pour une nombre de Rayleigh de 108, une condition initiale de type CIO90, en domaine ouvert, avec un maillage de 96×96 mailles .................................................................... 68

Annexe J : Résultats obtenus pour un nombre de Rayleigh d’une valeurs de 200, une condition initiale de type CIO108, en domaine ouvert, avec un maillage de 96×96 mailles ........................................ 69

Annexe K : Résultats obtenus pour une nombre de Rayleigh de 300, une condition initiale de type CIO200, en domaine ouvert, avec un maillage de 96×96 mailles .................................................................. 70

Page 61: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 61

Annexe A : Naissance et transport d’une vague de fluide chaud

t = 26212 ans t = 26270 ans t = 26328 ans

t = 26386 ans

t = 26619 ans t = 26561 ans

t = 26444 ans t = 26503 ans

478 °C 359 239 120

Page 62: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 62

Annexe B : Compléments sur l’analyse de la convection La superposition du champ de vitesses et du champ de températures permet de compléter

l’analyse des instabilités.

Champ de températures et de vitesses, zoom de la couche

limite au niveau de la formation d’un bourgeon, Ra = 800.

Cette superposition montre que la croissance des bourgeons chauds est effectivement couplée à une modification du champ de vitesses puisque ce dernier voit sa composante verticale fortement augmenter au sein du bourgeon. De plus, la conservation du débit massique est respectée en aval du bourgeon comme en témoigne la chute des valeurs des normes de vitesse. Le déplacement de ce bourgeon est similaire à l’avancée d’une vague de surface sur l’océan. Annexe C : Exemple de champ de pressions obtenu pour le cas Ra = 800

386 Bars 362 339 311 288

478 °C 359 239 120

Page 63: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 63

Annexe D : Champ de vitesses et normes de vitesses

Champ de vitesses au régime P2 (Ra = 800),

temps = 29945 ans (Maillage 64×64).

Champ de normes de vitesses, régime P2 (Ra = 800), temps = 29945 ans (Maillage 64×64).

La zone centrale et les angles du domaine sont à vitesses faibles de l’ordre de 6,33.10-2 m.j-1. Les masses se déplaçant verticalement le long des parois ont une vitesse de départ importante puis ralentissent à proximité des angles du domaine. Cela est cohérent avec les observations faites sur le champ de températures. En effet, l’étalement des masses chaudes ou froides remarqué au paragraphe 4.1.3 correspond au ralentissement observé couplé à une relative homogénéisation des températures ; tandis que les zones à forte vitesse correspondent à l’écrasement des masses contre les parois latérales et donc à une concentration de l’écoulement couplé à un fort différentiel de densité.

3,92 m.j -1 2,63 1,16 6,33.10-02

Page 64: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 64

Annexe E : Résultats des simulations numériques effectuées pour un nombre de Rayleigh de 800, avec un maillage de 64x64 mailles et des nombres de Fourier inférieurs à 1/3.

Tableau rassemblant les résultats obtenus pour un maillage de 64×64 mailles,

un nombre de Rayleigh de 800 et différents nombres de Fourier inférieurs à 1/3.

Evolution du nombre de Nusselt en fonction du temps pour un maillage de 64×64 mailles – Rayleigh 800 – Fourier 0,01.

Fourier utilisé ∆∆∆∆t* ∆∆∆∆t

(ans)

Valeurs moyennes du

Nusselt

Fréquences du Nusselt

Maillage 64×64 Crank-Nicholson Fo = 0,25 3,05x10 -5 3,63 9,11 273 Fo = 0,10 1,22x10 -5 1,45 9,1 280 Fo = 0,05 6,10x10 -6 0,73 9,09 283 Fo = 0,02 2,44x10 -6 0,29 9,08 286 Fo = 0,01 1,22x10 -6 0,15 6,31 149

9.0

8.0

7.0

6 .0

5.0

4 .0

3.0

2.0

1.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 temps

Nusselt

Page 65: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

65

Annexe F : Tableau récapitulatif des résultats de calculs obtenus avec CAST3M dans le cas d’un domaine ouvert

Domaine

ouvert

Méthode de

discrétisation Schéma temporel

Nombre de

Rayleigh

Maillage Nombre de

Fourier

∆∆∆∆t* ∆∆∆∆t (ans) Valeurs moyennes du nombre de Nusselt

Fréquences d’oscillation du nombre de Nusselt

Condition initiale

90 96x96 0,4 2,17x10-5 2,58 3,22 - CIO1 108 96x96 0,4 2,17x10-5 2,58 3,6 - CIO90 ou CIO1 200 96x96 0,4 2,17x10-5 2,58 6,49 - CIO108 ou CIO1 300 96x96 0,4 2,17x10-5 2,58 8,16 - CIO200 ou CIO1 400 96x96 0,4 2,17x10-5 2,58 9,16 - CIO300 ou CIO1 425 96x96 0,4 2,17x10-5 2,58 9,38 - CIO400 ou CIO1

EFMH Crank-Nicholson

470 96x96 0,4 2,17x10-5 2,58 9,85 121 CIO425 400 96x96 0,4 2,17x10-5 2,58 9,16 - CI1 425 96x96 0,4 2,17x10-5 2,58 9,35 - CIO400 ou CIO1 430 96x96 0,4 2,17x10-5 2,58 9,39 - CIO425 ou CIO1 450 96x96 0,4 2,17x10-5 2,58 9,53 - CIO430 ou CIO1 470 96x96 0,4 2,17x10-5 2,58 9,87 116 CIO450 ou CIO1 400 96x96 0,2 1,09 x10-5 1,29 9,16 - CIO1

VF Crank-Nicholson

470 96x96 0,2 1,09 x10-5 1,29 9,87 - CIO1

Page 66: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian

66

Annexe G : Graphique issue de (Cherkaoui et Wilcock, 1999) Ce graphique présente les résultats obtenus par Cherkaoui et Wilcock (1999) lors de la simulation d’un système identique à celui de notre étude dans le cas d’un domaine ouvert. Le nombre de Rayleigh croît progressivement à partir d’un nombre de Rayleigh initial variable. Les comportements observés diffèrent suivant la condition initiale. Les résultats de nos calculs pour un système identique sont représentés en rouge.

Page 67: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 67

Annexe H : Résultats obtenus pour une nombre de Rayleigh de 90, une condition initiale de type CI1, en domaine ouvert, avec un maillage de 96×96 mailles

Champ de températures en régime permanent,

un seul rouleau de convection est visible.

Evolution du nombre de Nusselt en fonction du temps (adimensionné), pour Ra = 90 et une condition initiale de type CI1.

3.50

3.00

2.50

2.00

1.50

1.00

0.00 0.20 0.40 0.60 0.80 1.00 1.20 Temps

Nusselt

Régime permanent établi

478 °C 359 239 120

Page 68: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 68

Annexe I : Résultats obtenus pour une nombre de Rayleigh de 108, une condition initiale de type CIO90, en domaine ouvert, avec un maillage de 96×96 mailles

Champ de températures en régime permanent, deux panaches de fluides chauds montent le long des parois latérales.

Evolution du nombre de Nusselt en fonction du temps (adimensionné), pour Ra =108 et une condition initiale de type CIO90.

4.20

4.00

3.80

3.60

3.40

3.20

3.00

0.80 1.00 1.20 1.40 1.60 1.80 2.00 Temps

Nusselt

Régime permanent établi

478 °C 359 239 120

Page 69: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 69

Annexe J : Résultats obtenus pour un nombre de Rayleigh d’une valeurs de 200, une condition initiale de type CIO108, en domaine ouvert, avec un maillage de 96×96 mailles

Champ de températures en régime permanent, deux panaches de

fluides chauds montent le long des parois latérales.

Evolution du nombre de Nusselt en fonction du temps (adimensionné), pour Ra =200 et une condition initiale de type CIO108.

7.00

6.50

6.00

5.50

5.00

4.50

4.00

3.50

3.00

2.50

1.00 1.20 1.40 1.60 1.80 2.00 2.20 2.40 2.60 2.80 3.00 Temps

Nusselt

Régime permanent établi

478 °C 359 239 120

Page 70: Rapport_de_Stage test

ENGEES – IPST : Mémoire de Travail de Fin d’Etudes Modélisation bidimensionnelle de la circulation hydrothermale au niveau de la dorsale lente médio-atlantique

PEREZ Florian 70

Annexe K : Résultats obtenus pour une nombre de Rayleigh de 300, une condition initiale de type CIO200, en domaine ouvert, avec un maillage de 96×96 mailles

Champ de températures en régime permanent, deux panaches de

fluides chauds montent le long des parois latérales.

Evolution du nombre de Nusselt en fonction du temps (adimensionné), pour Ra =300 et une condition initiale de type CIO200.

8.50

8.00

7.50

7.00

6.50

6.00

5.50

2.80 3.00 3.20 3.40 3.60 3.80 4.00 Temps

Régime permanent établi

478 °C 359 239 120