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 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. PEREZ Florian 1 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 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 CrankNicholson 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. PEREZ Florian 2 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 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. PEREZ Florian 3 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 Listes des notations Liste des abréviations µ 2D A EFMH Fo g h k K MAR Nu P Pe Ra z ρ VF : Viscosité dynamique du fluide [M][L-1][T-1] : Bidimensionnel : Rapport d’aspect, rapport largeur/hauteur des rouleaux de convection ou du domaine étudié : Eléments Finis Mixtes Hybrides : Nombre de Fourier : Accélération gravitaire [L].[T-2] P : Charge hydraulique [L], h = ρ g + z dans le cas d’un écoulement souterrain : Perméabilité intrinsèque [L²] : Perméabilité [L].[s-1] : Mid-Atlantic Ridge, Dorsale Médio-Atlantique : Nombre de Nusselt : Pression [M].[L-1].[T-²] : Nombre de Péclet : Nombre de Rayleigh : Altitude par rapport au niveau de référence [L] : Masse volumique du fluide [M].[L-3] : Volumes finis Liste des abréviations d’institutions ANR CEA CNRS IFREMER IMFS IPST LSCE UVSQ : Agence Nationale pour la Recherche : Commissariat à l’Energie Atomique : Centre National de la Recherche Scientifique : Institut Français de Recherche pour l’Exploitation de la MER : Institut de Mécanique des Fluides et des Solides : Institut Professionnel des Sciences et Technologies : Laboratoire des Sciences du Climat et de l’Environnement : Université Versailles St-Quentin PEREZ Florian 4 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 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 PEREZ Florian 5 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 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 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 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 PEREZ Florian 7 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 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 PEREZ Florian 8 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 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édioatlantique. 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. PEREZ Florian 9 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 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 PEREZ Florian 10 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 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 Dates du 26/01/09 au 6/02/09 Tâches Recherche Bibliographique Initiation à CAST3M Détails Compréhension de la problématique Décryptage d’un jeu de données CAST3M -Problème d’ELDER – rapport explicatif Encodage pour un problème hydrothermal simple -Cas d’une BOITE FERMEE – Schéma implicte en temps et Schéma de Crank-Nicholson. du 9/02/09 au 23/03/09 Recherche Bibliographique Encodage sous CAST3M 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. Tableau 1 : Tableau chronologique du déroulement du travail de fin d’études. PEREZ Florian 11 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 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. Dorsale océanique Chambre magmatique Lithosphère océanique Lithosphère continentale Cellule de convection Asthénosphère Possible transfert de masse du manteau supérieur au manteau inférieur Manteau inférieur Cellule de convection Noyau externe 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. 12 PEREZ Florian 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 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. PEREZ Florian 13 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 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édioatlantique ; (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 PEREZ Florian 14 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 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. PEREZ Florian 15 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 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. Si(OH)4 (mM) 6,9 67 0,14 Fe (µM) 24050 Mn (µM) 2250 Cu Zn (µM) (µM) 121/ 115/ 162 185 0,01 Profondeur R S W 2300 m T°C pH 365 2°C 2,8 7,8 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 PEREZ Florian 18 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 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. PEREZ Florian 19 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 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 RayleighBé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) : − dT gα T < dz Cp , (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. PEREZ Florian 20 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 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). Thaut(° C) Plaque plane supérieure (froide) d Rouleaux de convection g Tbas(° = T haut +∆T C) ∆ Plaque plane inférieure (chaude) Fluide (ρ, therm) ρ,µ,k ρ, 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 RayleighBé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 : τ th = d² ktherm , (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 : τ m = cte µ0 ρ 0 gdα∆T . (3) PEREZ Florian 21 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 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 : τ h ρ 0 gα d 3∆T = , τ m µ0 ktherm cste Ra = (4a) ρ 0 gα d 3∆T µ0 k therm > cste = Rac , (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 : µd τ m = 0 cste . (5) k0 g ∆ρ Le nombre de Rayleigh devient alors : Ra = ∆ρ gdk0 µ0 k therm . (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 : PEREZ Florian 22 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 L* Nu = ∫ 0 ∂T * dx * , * |z =0 ∂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. PEREZ Florian 23 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 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 : U= Q =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 : i = L , avec h, la charge (m) et L, la longueur du cylindre. On écrit généralement cette loi sous une forme vectorielle : ur uuuuu r U = − K ( gradh) . ∆h (9) PEREZ Florian 24 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 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 = ur U Ωρ µ , (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 : ur ∂ρφ + div ( ρ U ) = Q , ∂t (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 : φ (ρc)l uuuuu r ur ∂Tl − div(φ((Ddiffl + Ddispl ) gradT )l − (ρc)l UTl ) = Ql , ∂t (12) où Ql, représente le terme puits/source, (ρc)l, la capacité thermique volumique du fluide, Tl, 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. 25 PEREZ Florian 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 Pour la phase solide, cette équation se résume à : (1 − φ )( ρ c ) s uuuuu r ∂Ts − div ((1 − φ ) Ddiffs gradTs ) = Qs , ∂t (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, Qs = −Ql . 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 : (φ ( ρ c) l + (1 − φ )( ρ c) s ) ( ρ c) eq uuuuu r ur ∂Τ = div ((λeq gradT − ( ρ c) l UT )) , ∂t (14a) uuuuu r ur ∂Τ = div ((λeq gradT − ( ρ c) l UT )) , ∂t (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é : λeq = φ ( Ddiffl + Ddispl ) + (1 − φ ) Ddiffs . ρc correspond à la capacité thermique volumique équivalente du milieu saturé : ( ρ c)eq = φ ( ρ c) l + (1 − φ )( ρ c) s . (14c) (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. PEREZ Florian 26 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 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 : kint = 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). PEREZ Florian 27 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 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 = ( ρ c )eq . Ce paramètre prend toute son importance en régime transitoire car il traduit la k therm 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 : PEREZ Florian 28 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 ur r uuuuu r k uuuuu U = − ( grad P + ρg gradz ) , µ (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 à : ur div (U ) = 0 . 2.5.2 (16) 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 ( 1 - α (T - T0 ) ) , où T0 est la température de référence pour laquelle ρ = ρ0. 2.5.3 (17) 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 : uuuuu r ur ∂Τ = div ((k therm gradT − UT )) . ∂t 2.5.4 (18) 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 : ∆ρ = ρ (Tmin ) − ρ (Tmax ) , x = d × x *, y = d × y *, V = K ×V * , T = ∆T × T * , ρ = ∆ρ × ρ * , P = ∆ρ gd × P * , .t = d ² t * k therm ∆T = Tmax − Tmin , PEREZ Florian 29 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 Le nombre de Rayleigh est ici défini par la relation (6). Les équations obtenues sont : ur* uuuuu * r r ρ∗ u U = −( grad P * + g) . La vitesse de Darcy : (19) g ur * div * (U ) = 0 . L’équation de conservation de la masse : (20) * uuuuuu r uuuuuu r ur * ∂T + RaU . grad * (T * ) = div * ( grad *T * ) . (21) L’équation de transfert de la chaleur : * ∂t * ∗ * * L’équation de dépendance de la masse volumique : ρ = ρ0 1 - α∆T (T - T0 ) . (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. PEREZ Florian 30 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 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. PEREZ Florian 31 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 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 diffusionconvection pour le transport de la chaleur joue le rôle de la vitesse de Darcy : uuuuuu r uuuuu r U = − K gradh ≡ Q = − grad *T * . (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étaméthode : t n +1 ∫ t n dT dt = dt t n +1 ∫ t n f (T ( t )) dt , (24) T n+1 − T n = θ f (T n+1 ) + (1 − θ ) f (T n ) . n +1 n 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 : T n +1 −T ∆t n = f (T n +1 ). (26) Le schéma temporel de Crank-Nicholson est un schéma du second ordre en temps : T n +1 −T ∆t n = f (T n +1 ) + f (T n ) . 2 (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. PEREZ Florian 32 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 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 : 2ktherm ∆t , (28) ∆x ² 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é : Fo = Fo = 2 ∆t * . ∆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 : CFL= uuu r U ∆t ∆x = uuuur U * ∆t * ∆x * × k k therm . (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 : Pe = 2CFL . Fo (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 PEREZ Florian 33 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 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. 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. Maillage homogène carré : Fourier minimum suivant le critère d’Hoteit et al. (2002) 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. PEREZ Florian 34 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 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. Création du maillage Définition des conditions aux limites et initiales, et définition des paramètres Variables entrantes à la nième itération en temps ρn, Tn Calcul hydraulique Un+1 ρn+1 Pas de temps suivant Tn+1 Calcul du transport thermique Tn+1 Calcul de la masse volumique ρn+1 Mise en mémoire du champ de température 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. PEREZ Florian 35 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 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 : T= 2°C, P = 300 bars, u = v = 0 T= 0, P= 2,25, u = v = 0 u = v = 0, d = 1500m u=v=0 Milieu poreux homogène isotrope z u=v=0 d=1 Milieu poreux homogène isotrope u=v=0 x T = 600 °C, u = v = 0, L = 1500m T = 1, u = v = 0, L = 1 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. 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 : Accélération gravitaire d’axe y : g Masse volumique de l’eau de mer à 300 bars et 2°C : ρ0 Coefficient de dilatation volumique: α Viscosité de l’eau pure : µ Capacité thermique volumique de l’eau : (ρc)f Capacité thermique volumique du milieu saturé : (ρc)eq Conductivité thermique du milieu saturé λeq Diffusivité thermique ktherm Perméabilité intrinsèque du milieu poreux k : Valeurs 9,81 m.s-2 130 kg/m3 1,17×10-2K-1 5×10-5 kg.m-1.s-1 4,2×10-6 J.k-1.m-3 ≈ 4,2×10-6 J.k-1.m-3 2,5 W.m-1.K-1 6×10-7 m²/s 4×10-13 à 4×10-15 m² Tableau 5 : Paramètres du système. PEREZ Florian 36 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 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 478 ° C 478 ° C 359 359 239 239 120 120 Figure 11 : Champ de températures obtenu avec CAST3M (EFMH). Figure 12 : Champ de températures publiés dans (Fontaine et al., 2007). 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 PEREZ Florian 37 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 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 1,615 Figure 13 : Champ de nombres de Péclet pour un maillage de 128×128, Fo = 0,5. Pe = U ∆x ktherm = U * ∆x * Ra 1,09 0,480 2,61 10-02 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. PEREZ Florian 38 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 Domaine fermé Méthode de Schéma Nombre Maillage Nombre de Fourier discrétisation temporel de utilisée Rayleigh 96x96 128x128 128x128 96x96 96x96 128x128 128x128 128x128 128x128 96x96 128x128 128x128 96x96 128x128 64x64 96x96 128x128 96x96 128x128 0,4 0,4 0,7 0,4 0,5 0,25 0,4 0,5 1 0,4 0,4 0,4 0,4 0,4 0,4 0,4 0,25 0,4 0,4 ∆t* ∆t (ans) Valeurs moyennes du nombre de Nusselt 2,58 2,58 2,54 2,58 3,23 0,91 1,45 1,81 3,63 2,58 1,45 1,45 2,58 1,45 5,81 2,58 9,01 9,08 9,06 9,16 9,16 9,1 9,16 9,07 9,04 10,90 8,81 11,00 13,6 13,6 9,11 9,13 9,1 8,89 9,04 Implicite 800 800 EMFH CrankNicholson 950 1200 VF CrankNicholson 800 2,17x10-5 2,17x10-5 2,14x10-5 2,17x10-5 2,71x10-5 7,63x10-6 1,22x10-5 1,53x10-5 3,05x10-5 2,17x10-5 1,22x10-5 1,22x10-5 2,17x10-5 1,22x10-5 4,88x10-5 2,17x10-5 Fréquences d’oscillation du nombre de Nusselt 273 282 275 281 277 277 286 277 260 399 & 212 400 & 217 595 & 230 599 & 232 262 279 279 260 Condition initiale CI1 CI1 CI1 CI1 CI1 CI1 CI1 CI1 CI1 CIF800 CI1 CIF800 CIF950 CIF950 CI1 CI1 CI1 CI1 CI1 950 7,63x10-6 3,63 2,17x10-5 2,58 1,22x10-5 1,45 Tableau 6 : Tableau récapitulatif des résultats de calculs obtenus avec CAST3M dans le cas d’un domaine fermé 39 PEREZ Florian 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 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. Nusselt 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 Figure 15 : Champ de températures permanent obtenu pour Ra = 950 et une condition initiale de type CI1. Figure 14 : Evolution du nombre de Nusselt en fonction du temps, pour Ra = 950 et une condition initiale de type CI1. 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. PEREZ Florian 40 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 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. Nusselt 12.00 478 ° C 10.00 359 8.00 6.00 Transition d’un nombre de Rayleigh de 800 à 950 239 4.00 2.00 0.00 0.00 0.50 1.00 1.50 2.00 2.50 3.00 Temps 120 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. 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. 478 ° C Nusselt 15.50 15.00 359 14.50 14.00 239 13.50 13.00 120 12.50 12.00 1.34 2.36 2.38 2.40 2.42 2.44 2.46 2.48 2.50 Temps 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. PEREZ Florian 41 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 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 Nu f1 f2 Nu f1 f2 Nu f1 f2 Nu f1 f2 f3 800 9,1 280 9,2 275 9,14 299 9,43 296 950 10,9 399 212 10,9 397 200 1200 13,6 595 230 13,06 608 228 Cherkaoui et Wilcock (1999) 13,85 695 261 69 Caltagirone et Fabrie (1989) Fontaine et al. (2007) Cette étude (CAST3M, EFMH) 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. PEREZ Florian 42 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 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. Nusselt 11.00 478 ° C 10.00 9.00 359 8.00 7.00 6.00 239 5.00 4.00 3.00 2.00 1.00 0.60 0.70 Temps Figure 21 : Evolution du nombre de Nusselt en fonction du temps, pour Ra = 800 et une condition initiale de type CI1. 0.00 0.10 0.20 0.30 0.40 0.50 120 Figure 20 : Champ de températures obtenu pour Ra = 800 et une condition initiale de type CI1. 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). PEREZ Florian 43 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 478 ° C Nusselt 16.00 14.00 12.00 359 10.00 8.00 6.00 4.00 2.00 0.00 0.00 0.80 Temps Figure 23 : Evolution du nombre de Nusselt en fonction du temps, pour Ra = 950 et une condition initiale de type CI1. 0.20 0.40 0.60 239 120 Figure 22 : Champ de températures permanent obtenu pour Ra = 950 et une condition initiale de type CI1. 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). Nusselt 12.00 478 ° C 11.60 359 11.20 10.80 10.40 239 120 10.00 1.20 Temps Figure 25 : Evolution du nombre de Nusselt en fonction du temps, pour Ra = 950 et une condition initiale de type CIF800. 0.95 1.00 1.05 1.10 1.15 Figure 24 : Champ de températures permanent obtenu pour Ra = 950 et une condition initiale de type CIF800. 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. PEREZ Florian 44 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 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 Nu f1 f2 Nu f1 f2 800 9,1 280 9,2 279 950 10,9 399 212 11,1 403 203 CAST3M EFMH 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 P= 2,25 Pa T= 0 ∂T/∂z= 0 u = v = 0, d = 1 Milieu poreux homogène isotrope u=v=0 T = 1, u = v = 0, L = 1 Figure 26: Schéma du domaine fermé et conditions aux limites adimensionnées. PEREZ Florian 45 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 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. Nusselt 12.00 10.00 Ra = 90 108 200 300 400 425 470 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 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. PEREZ Florian 46 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 Ra Nu f1 Amplitude des oscillations rapport d'aspect A Nu f1 Amplitude des oscillations rapport d'aspect A 90 3,21 1 3,20* - 108 3,6 1/2 3,50* - Domaine ouvert 200 300 400 8,16 9,16 6,49 1/2 6,40* - Références 425 9,38 1/2 9,30* - 470 9,85 121 1,82 1/2 9,80* 117,2 1,18* 1/2 Cherkaoui et Wilcock (1999) Cette étude 1/2 8,2* - 1/2 9,10* - 1 1/2 1/2 1/2 1/2 1/2 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 : 478 ° C 359 239 120 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 PEREZ Florian 47 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 diffusion de la chaleur en régime permanent (condition de type CIO1), on obtient exactement les mêmes résultats. Domaine ouvert 200 6,49 1/2 Références 300 8,16 1/2 400 9,16 1/2 Ra Nu f1 rapport d'aspect A 90 3,21 1 108 3,6 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 : Nusselt 11.00 10.50 Ra = 10.00 400 425 430 450 470 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 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. PEREZ Florian 48 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 Ra Nu f1 Amplitude des oscillations rapport d'aspect A Nu f1 Amplitude des oscillations rapport d'aspect A 400 9 ,16 1/2 9,10* - Domaine ouvert 425 430 9,35 1/2 9,30* - Références 450 9,53 1/2 9,60* - 470 9,87 116,1 1,53* 1/2 9,80* 117,2 Cherkaoui et Wilcock 1,18* (1999) 1/2 Cette étude 9,39 1/2 9,42* - 1 1/2 1/2 1/2 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 Ra Nu f1 Amplitude des oscillations rapport d'aspect A 400 9 ,16 1/2 425 9,35 1/2 430 9,39 1/2 450 9,53 1/2 470 9,87 116,1 1,53* 1/2 Cette étude Références 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. PEREZ Florian 49 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 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 CrankNicholson 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 PEREZ Florian 50 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 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 PEREZ Florian 51 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 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 PEREZ Florian 52 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 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. PEREZ Florian 53 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 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 JeanLuc 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. PEREZ Florian 54 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 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. PEREZ Florian 55 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 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. PEREZ Florian 56 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 Cherkaoui, A.S.M.,. Wilcock, W.S.D., (1999). Characteristics of high Rayleigh number twodimensional 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 midocean 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. PEREZ Florian 57 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 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, 193215. 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. PEREZ Florian 58 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 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-fondsoceaniques.xml consulté le 21 juin 2009. PEREZ Florian 59 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 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 PEREZ Florian 60 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 Annexe A : Naissance et transport d’une vague de fluide chaud t = 26212 ans t = 26270 ans t = 26328 ans t = 26386 ans t = 26444 ans 478 ° C t = 26503 ans 359 239 120 t = 26561 ans t = 26619 ans PEREZ Florian 61 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 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. 478 ° C 359 239 120 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 PEREZ Florian 62 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 Annexe D : Champ de vitesses et normes de vitesses Champ de vitesses au régime P2 (Ra = 800), temps = 29945 ans (Maillage 64×64). 3,92 m.j-1 2,63 1,16 6,33.10-02 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é. PEREZ Florian 63 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 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. ∆t (ans) -5 Fourier utilisé Maillage 64×64 Crank-Nicholson ∆t* Valeurs moyennes du Nusselt 3,63 1,45 0,73 0,29 0,15 9,11 9,1 9,09 9,08 6,31 Fréquences du Nusselt 273 280 283 286 149 Fo = 0,25 Fo = 0,10 Fo = 0,05 Fo = 0,02 Fo = 0,01 3,05x10 -5 1,22x10 -6 6,10x10 -6 2,44x10 -6 1,22x10 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. Nusselt 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 Evolution du nombre de Nusselt en fonction du temps pour un maillage de 64×64 mailles – Rayleigh 800 – Fourier 0,01. PEREZ Florian 64 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 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 EFMH CrankNicholson VF CrankNicholson Nombre Maillage de Rayleig h 90 96x96 108 96x96 200 96x96 300 96x96 400 96x96 425 96x96 470 96x96 400 96x96 425 96x96 430 96x96 450 96x96 470 96x96 400 96x96 470 96x96 Nombre de Fourier 0,4 0,4 0,4 0,4 0,4 0,4 0,4 0,4 0,4 0,4 0,4 0,4 0,2 0,2 ∆t* ∆t (ans) 2,17x10-5 2,17x10-5 2,17x10-5 2,17x10-5 2,17x10-5 2,17x10-5 2,17x10-5 2,17x10-5 2,17x10-5 2,17x10-5 2,17x10-5 2,17x10-5 1,09 x10-5 1,09 x10-5 2,58 2,58 2,58 2,58 2,58 2,58 2,58 2,58 2,58 2,58 2,58 2,58 1,29 1,29 Valeurs moyennes du nombre de Nusselt 3,22 3,6 6,49 8,16 9,16 9,38 9,85 9,16 9,35 9,39 9,53 9,87 9,16 9,87 Fréquences d’oscillation du nombre de Nusselt 121 116 - Condition initiale CIO1 CIO90 ou CIO1 CIO108 ou CIO1 CIO200 ou CIO1 CIO300 ou CIO1 CIO400 ou CIO1 CIO425 CI1 CIO400 ou CIO1 CIO425 ou CIO1 CIO430 ou CIO1 CIO450 ou CIO1 CIO1 CIO1 65 PEREZ Florian 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 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. 66 PEREZ Florian 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 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 478 ° C 359 239 120 Champ de températures en régime permanent, un seul rouleau de convection est visible. Nusselt 3.50 Régime permanent établi 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 Evolution du nombre de Nusselt en fonction du temps (adimensionné), pour Ra = 90 et une condition initiale de type CI1. PEREZ Florian 67 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 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 478 ° C 359 239 120 Champ de températures en régime permanent, deux panaches de fluides chauds montent le long des parois latérales. Nusselt 4.20 4.00 3.80 Régime permanent établi 3.60 3.40 3.20 3.00 0.80 1.00 1.20 1.40 1.60 1.80 2.00 Temps Evolution du nombre de Nusselt en fonction du temps (adimensionné), pour Ra =108 et une condition initiale de type CIO90. PEREZ Florian 68 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 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 478 ° C 359 239 120 Champ de températures en régime permanent, deux panaches de fluides chauds montent le long des parois latérales. Nusselt 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 Régime permanent établi Evolution du nombre de Nusselt en fonction du temps (adimensionné), pour Ra =200 et une condition initiale de type CIO108. PEREZ Florian 69 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 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 478 ° C 359 239 120 Champ de températures en régime permanent, deux panaches de fluides chauds montent le long des parois latérales. 8.50 8.00 Régime permanent établi 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 Evolution du nombre de Nusselt en fonction du temps (adimensionné), pour Ra =300 et une condition initiale de type CIO200. PEREZ Florian 70
Please download to view
All materials on our website are shared by users. If you have any questions about copyright issues, please report us to resolve them. We are always happy to assist you.
...

Rapport_de_Stage test

by chwarzi-nari

on

Report

Category:

Documents

Download: 0

Comment: 0

112

views

Comments

Description

Download Rapport_de_Stage test

Transcript

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 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. PEREZ Florian 1 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 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 CrankNicholson 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. PEREZ Florian 2 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 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. PEREZ Florian 3 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 Listes des notations Liste des abréviations µ 2D A EFMH Fo g h k K MAR Nu P Pe Ra z ρ VF : Viscosité dynamique du fluide [M][L-1][T-1] : Bidimensionnel : Rapport d’aspect, rapport largeur/hauteur des rouleaux de convection ou du domaine étudié : Eléments Finis Mixtes Hybrides : Nombre de Fourier : Accélération gravitaire [L].[T-2] P : Charge hydraulique [L], h = ρ g + z dans le cas d’un écoulement souterrain : Perméabilité intrinsèque [L²] : Perméabilité [L].[s-1] : Mid-Atlantic Ridge, Dorsale Médio-Atlantique : Nombre de Nusselt : Pression [M].[L-1].[T-²] : Nombre de Péclet : Nombre de Rayleigh : Altitude par rapport au niveau de référence [L] : Masse volumique du fluide [M].[L-3] : Volumes finis Liste des abréviations d’institutions ANR CEA CNRS IFREMER IMFS IPST LSCE UVSQ : Agence Nationale pour la Recherche : Commissariat à l’Energie Atomique : Centre National de la Recherche Scientifique : Institut Français de Recherche pour l’Exploitation de la MER : Institut de Mécanique des Fluides et des Solides : Institut Professionnel des Sciences et Technologies : Laboratoire des Sciences du Climat et de l’Environnement : Université Versailles St-Quentin PEREZ Florian 4 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 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 PEREZ Florian 5 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 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 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 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 PEREZ Florian 7 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 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 PEREZ Florian 8 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 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édioatlantique. 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. PEREZ Florian 9 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 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 PEREZ Florian 10 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 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 Dates du 26/01/09 au 6/02/09 Tâches Recherche Bibliographique Initiation à CAST3M Détails Compréhension de la problématique Décryptage d’un jeu de données CAST3M -Problème d’ELDER – rapport explicatif Encodage pour un problème hydrothermal simple -Cas d’une BOITE FERMEE – Schéma implicte en temps et Schéma de Crank-Nicholson. du 9/02/09 au 23/03/09 Recherche Bibliographique Encodage sous CAST3M 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. Tableau 1 : Tableau chronologique du déroulement du travail de fin d’études. PEREZ Florian 11 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 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. Dorsale océanique Chambre magmatique Lithosphère océanique Lithosphère continentale Cellule de convection Asthénosphère Possible transfert de masse du manteau supérieur au manteau inférieur Manteau inférieur Cellule de convection Noyau externe 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. 12 PEREZ Florian 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 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. PEREZ Florian 13 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 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édioatlantique ; (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 PEREZ Florian 14 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 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. PEREZ Florian 15 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 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. Si(OH)4 (mM) 6,9 67 0,14 Fe (µM) 24050 Mn (µM) 2250 Cu Zn (µM) (µM) 121/ 115/ 162 185 0,01 Profondeur R S W 2300 m T°C pH 365 2°C 2,8 7,8 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 PEREZ Florian 18 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 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. PEREZ Florian 19 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 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 RayleighBé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) : − dT gα T < dz Cp , (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. PEREZ Florian 20 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 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). Thaut(° C) Plaque plane supérieure (froide) d Rouleaux de convection g Tbas(° = T haut +∆T C) ∆ Plaque plane inférieure (chaude) Fluide (ρ, therm) ρ,µ,k ρ, 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 RayleighBé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 : τ th = d² ktherm , (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 : τ m = cte µ0 ρ 0 gdα∆T . (3) PEREZ Florian 21 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 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 : τ h ρ 0 gα d 3∆T = , τ m µ0 ktherm cste Ra = (4a) ρ 0 gα d 3∆T µ0 k therm > cste = Rac , (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 : µd τ m = 0 cste . (5) k0 g ∆ρ Le nombre de Rayleigh devient alors : Ra = ∆ρ gdk0 µ0 k therm . (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 : PEREZ Florian 22 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 L* Nu = ∫ 0 ∂T * dx * , * |z =0 ∂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. PEREZ Florian 23 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 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 : U= Q =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 : i = L , avec h, la charge (m) et L, la longueur du cylindre. On écrit généralement cette loi sous une forme vectorielle : ur uuuuu r U = − K ( gradh) . ∆h (9) PEREZ Florian 24 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 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 = ur U Ωρ µ , (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 : ur ∂ρφ + div ( ρ U ) = Q , ∂t (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 : φ (ρc)l uuuuu r ur ∂Tl − div(φ((Ddiffl + Ddispl ) gradT )l − (ρc)l UTl ) = Ql , ∂t (12) où Ql, représente le terme puits/source, (ρc)l, la capacité thermique volumique du fluide, Tl, 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. 25 PEREZ Florian 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 Pour la phase solide, cette équation se résume à : (1 − φ )( ρ c ) s uuuuu r ∂Ts − div ((1 − φ ) Ddiffs gradTs ) = Qs , ∂t (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, Qs = −Ql . 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 : (φ ( ρ c) l + (1 − φ )( ρ c) s ) ( ρ c) eq uuuuu r ur ∂Τ = div ((λeq gradT − ( ρ c) l UT )) , ∂t (14a) uuuuu r ur ∂Τ = div ((λeq gradT − ( ρ c) l UT )) , ∂t (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é : λeq = φ ( Ddiffl + Ddispl ) + (1 − φ ) Ddiffs . ρc correspond à la capacité thermique volumique équivalente du milieu saturé : ( ρ c)eq = φ ( ρ c) l + (1 − φ )( ρ c) s . (14c) (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. PEREZ Florian 26 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 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 : kint = 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). PEREZ Florian 27 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 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 = ( ρ c )eq . Ce paramètre prend toute son importance en régime transitoire car il traduit la k therm 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 : PEREZ Florian 28 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 ur r uuuuu r k uuuuu U = − ( grad P + ρg gradz ) , µ (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 à : ur div (U ) = 0 . 2.5.2 (16) 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 ( 1 - α (T - T0 ) ) , où T0 est la température de référence pour laquelle ρ = ρ0. 2.5.3 (17) 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 : uuuuu r ur ∂Τ = div ((k therm gradT − UT )) . ∂t 2.5.4 (18) 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 : ∆ρ = ρ (Tmin ) − ρ (Tmax ) , x = d × x *, y = d × y *, V = K ×V * , T = ∆T × T * , ρ = ∆ρ × ρ * , P = ∆ρ gd × P * , .t = d ² t * k therm ∆T = Tmax − Tmin , PEREZ Florian 29 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 Le nombre de Rayleigh est ici défini par la relation (6). Les équations obtenues sont : ur* uuuuu * r r ρ∗ u U = −( grad P * + g) . La vitesse de Darcy : (19) g ur * div * (U ) = 0 . L’équation de conservation de la masse : (20) * uuuuuu r uuuuuu r ur * ∂T + RaU . grad * (T * ) = div * ( grad *T * ) . (21) L’équation de transfert de la chaleur : * ∂t * ∗ * * L’équation de dépendance de la masse volumique : ρ = ρ0 1 - α∆T (T - T0 ) . (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. PEREZ Florian 30 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 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. PEREZ Florian 31 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 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 diffusionconvection pour le transport de la chaleur joue le rôle de la vitesse de Darcy : uuuuuu r uuuuu r U = − K gradh ≡ Q = − grad *T * . (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étaméthode : t n +1 ∫ t n dT dt = dt t n +1 ∫ t n f (T ( t )) dt , (24) T n+1 − T n = θ f (T n+1 ) + (1 − θ ) f (T n ) . n +1 n 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 : T n +1 −T ∆t n = f (T n +1 ). (26) Le schéma temporel de Crank-Nicholson est un schéma du second ordre en temps : T n +1 −T ∆t n = f (T n +1 ) + f (T n ) . 2 (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. PEREZ Florian 32 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 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 : 2ktherm ∆t , (28) ∆x ² 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é : Fo = Fo = 2 ∆t * . ∆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 : CFL= uuu r U ∆t ∆x = uuuur U * ∆t * ∆x * × k k therm . (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 : Pe = 2CFL . Fo (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 PEREZ Florian 33 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 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. 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. Maillage homogène carré : Fourier minimum suivant le critère d’Hoteit et al. (2002) 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. PEREZ Florian 34 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 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. Création du maillage Définition des conditions aux limites et initiales, et définition des paramètres Variables entrantes à la nième itération en temps ρn, Tn Calcul hydraulique Un+1 ρn+1 Pas de temps suivant Tn+1 Calcul du transport thermique Tn+1 Calcul de la masse volumique ρn+1 Mise en mémoire du champ de température 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. PEREZ Florian 35 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 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 : T= 2°C, P = 300 bars, u = v = 0 T= 0, P= 2,25, u = v = 0 u = v = 0, d = 1500m u=v=0 Milieu poreux homogène isotrope z u=v=0 d=1 Milieu poreux homogène isotrope u=v=0 x T = 600 °C, u = v = 0, L = 1500m T = 1, u = v = 0, L = 1 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. 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 : Accélération gravitaire d’axe y : g Masse volumique de l’eau de mer à 300 bars et 2°C : ρ0 Coefficient de dilatation volumique: α Viscosité de l’eau pure : µ Capacité thermique volumique de l’eau : (ρc)f Capacité thermique volumique du milieu saturé : (ρc)eq Conductivité thermique du milieu saturé λeq Diffusivité thermique ktherm Perméabilité intrinsèque du milieu poreux k : Valeurs 9,81 m.s-2 130 kg/m3 1,17×10-2K-1 5×10-5 kg.m-1.s-1 4,2×10-6 J.k-1.m-3 ≈ 4,2×10-6 J.k-1.m-3 2,5 W.m-1.K-1 6×10-7 m²/s 4×10-13 à 4×10-15 m² Tableau 5 : Paramètres du système. PEREZ Florian 36 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 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 478 ° C 478 ° C 359 359 239 239 120 120 Figure 11 : Champ de températures obtenu avec CAST3M (EFMH). Figure 12 : Champ de températures publiés dans (Fontaine et al., 2007). 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 PEREZ Florian 37 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 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 1,615 Figure 13 : Champ de nombres de Péclet pour un maillage de 128×128, Fo = 0,5. Pe = U ∆x ktherm = U * ∆x * Ra 1,09 0,480 2,61 10-02 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. PEREZ Florian 38 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 Domaine fermé Méthode de Schéma Nombre Maillage Nombre de Fourier discrétisation temporel de utilisée Rayleigh 96x96 128x128 128x128 96x96 96x96 128x128 128x128 128x128 128x128 96x96 128x128 128x128 96x96 128x128 64x64 96x96 128x128 96x96 128x128 0,4 0,4 0,7 0,4 0,5 0,25 0,4 0,5 1 0,4 0,4 0,4 0,4 0,4 0,4 0,4 0,25 0,4 0,4 ∆t* ∆t (ans) Valeurs moyennes du nombre de Nusselt 2,58 2,58 2,54 2,58 3,23 0,91 1,45 1,81 3,63 2,58 1,45 1,45 2,58 1,45 5,81 2,58 9,01 9,08 9,06 9,16 9,16 9,1 9,16 9,07 9,04 10,90 8,81 11,00 13,6 13,6 9,11 9,13 9,1 8,89 9,04 Implicite 800 800 EMFH CrankNicholson 950 1200 VF CrankNicholson 800 2,17x10-5 2,17x10-5 2,14x10-5 2,17x10-5 2,71x10-5 7,63x10-6 1,22x10-5 1,53x10-5 3,05x10-5 2,17x10-5 1,22x10-5 1,22x10-5 2,17x10-5 1,22x10-5 4,88x10-5 2,17x10-5 Fréquences d’oscillation du nombre de Nusselt 273 282 275 281 277 277 286 277 260 399 & 212 400 & 217 595 & 230 599 & 232 262 279 279 260 Condition initiale CI1 CI1 CI1 CI1 CI1 CI1 CI1 CI1 CI1 CIF800 CI1 CIF800 CIF950 CIF950 CI1 CI1 CI1 CI1 CI1 950 7,63x10-6 3,63 2,17x10-5 2,58 1,22x10-5 1,45 Tableau 6 : Tableau récapitulatif des résultats de calculs obtenus avec CAST3M dans le cas d’un domaine fermé 39 PEREZ Florian 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 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. Nusselt 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 Figure 15 : Champ de températures permanent obtenu pour Ra = 950 et une condition initiale de type CI1. Figure 14 : Evolution du nombre de Nusselt en fonction du temps, pour Ra = 950 et une condition initiale de type CI1. 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. PEREZ Florian 40 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 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. Nusselt 12.00 478 ° C 10.00 359 8.00 6.00 Transition d’un nombre de Rayleigh de 800 à 950 239 4.00 2.00 0.00 0.00 0.50 1.00 1.50 2.00 2.50 3.00 Temps 120 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. 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. 478 ° C Nusselt 15.50 15.00 359 14.50 14.00 239 13.50 13.00 120 12.50 12.00 1.34 2.36 2.38 2.40 2.42 2.44 2.46 2.48 2.50 Temps 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. PEREZ Florian 41 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 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 Nu f1 f2 Nu f1 f2 Nu f1 f2 Nu f1 f2 f3 800 9,1 280 9,2 275 9,14 299 9,43 296 950 10,9 399 212 10,9 397 200 1200 13,6 595 230 13,06 608 228 Cherkaoui et Wilcock (1999) 13,85 695 261 69 Caltagirone et Fabrie (1989) Fontaine et al. (2007) Cette étude (CAST3M, EFMH) 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. PEREZ Florian 42 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 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. Nusselt 11.00 478 ° C 10.00 9.00 359 8.00 7.00 6.00 239 5.00 4.00 3.00 2.00 1.00 0.60 0.70 Temps Figure 21 : Evolution du nombre de Nusselt en fonction du temps, pour Ra = 800 et une condition initiale de type CI1. 0.00 0.10 0.20 0.30 0.40 0.50 120 Figure 20 : Champ de températures obtenu pour Ra = 800 et une condition initiale de type CI1. 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). PEREZ Florian 43 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 478 ° C Nusselt 16.00 14.00 12.00 359 10.00 8.00 6.00 4.00 2.00 0.00 0.00 0.80 Temps Figure 23 : Evolution du nombre de Nusselt en fonction du temps, pour Ra = 950 et une condition initiale de type CI1. 0.20 0.40 0.60 239 120 Figure 22 : Champ de températures permanent obtenu pour Ra = 950 et une condition initiale de type CI1. 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). Nusselt 12.00 478 ° C 11.60 359 11.20 10.80 10.40 239 120 10.00 1.20 Temps Figure 25 : Evolution du nombre de Nusselt en fonction du temps, pour Ra = 950 et une condition initiale de type CIF800. 0.95 1.00 1.05 1.10 1.15 Figure 24 : Champ de températures permanent obtenu pour Ra = 950 et une condition initiale de type CIF800. 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. PEREZ Florian 44 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 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 Nu f1 f2 Nu f1 f2 800 9,1 280 9,2 279 950 10,9 399 212 11,1 403 203 CAST3M EFMH 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 P= 2,25 Pa T= 0 ∂T/∂z= 0 u = v = 0, d = 1 Milieu poreux homogène isotrope u=v=0 T = 1, u = v = 0, L = 1 Figure 26: Schéma du domaine fermé et conditions aux limites adimensionnées. PEREZ Florian 45 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 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. Nusselt 12.00 10.00 Ra = 90 108 200 300 400 425 470 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 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. PEREZ Florian 46 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 Ra Nu f1 Amplitude des oscillations rapport d'aspect A Nu f1 Amplitude des oscillations rapport d'aspect A 90 3,21 1 3,20* - 108 3,6 1/2 3,50* - Domaine ouvert 200 300 400 8,16 9,16 6,49 1/2 6,40* - Références 425 9,38 1/2 9,30* - 470 9,85 121 1,82 1/2 9,80* 117,2 1,18* 1/2 Cherkaoui et Wilcock (1999) Cette étude 1/2 8,2* - 1/2 9,10* - 1 1/2 1/2 1/2 1/2 1/2 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 : 478 ° C 359 239 120 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 PEREZ Florian 47 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 diffusion de la chaleur en régime permanent (condition de type CIO1), on obtient exactement les mêmes résultats. Domaine ouvert 200 6,49 1/2 Références 300 8,16 1/2 400 9,16 1/2 Ra Nu f1 rapport d'aspect A 90 3,21 1 108 3,6 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 : Nusselt 11.00 10.50 Ra = 10.00 400 425 430 450 470 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 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. PEREZ Florian 48 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 Ra Nu f1 Amplitude des oscillations rapport d'aspect A Nu f1 Amplitude des oscillations rapport d'aspect A 400 9 ,16 1/2 9,10* - Domaine ouvert 425 430 9,35 1/2 9,30* - Références 450 9,53 1/2 9,60* - 470 9,87 116,1 1,53* 1/2 9,80* 117,2 Cherkaoui et Wilcock 1,18* (1999) 1/2 Cette étude 9,39 1/2 9,42* - 1 1/2 1/2 1/2 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 Ra Nu f1 Amplitude des oscillations rapport d'aspect A 400 9 ,16 1/2 425 9,35 1/2 430 9,39 1/2 450 9,53 1/2 470 9,87 116,1 1,53* 1/2 Cette étude Références 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. PEREZ Florian 49 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 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 CrankNicholson 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 PEREZ Florian 50 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 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 PEREZ Florian 51 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 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 PEREZ Florian 52 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 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. PEREZ Florian 53 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 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 JeanLuc 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. PEREZ Florian 54 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 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. PEREZ Florian 55 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 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. PEREZ Florian 56 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 Cherkaoui, A.S.M.,. Wilcock, W.S.D., (1999). Characteristics of high Rayleigh number twodimensional 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 midocean 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. PEREZ Florian 57 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 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, 193215. 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. PEREZ Florian 58 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 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-fondsoceaniques.xml consulté le 21 juin 2009. PEREZ Florian 59 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 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 PEREZ Florian 60 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 Annexe A : Naissance et transport d’une vague de fluide chaud t = 26212 ans t = 26270 ans t = 26328 ans t = 26386 ans t = 26444 ans 478 ° C t = 26503 ans 359 239 120 t = 26561 ans t = 26619 ans PEREZ Florian 61 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 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. 478 ° C 359 239 120 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 PEREZ Florian 62 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 Annexe D : Champ de vitesses et normes de vitesses Champ de vitesses au régime P2 (Ra = 800), temps = 29945 ans (Maillage 64×64). 3,92 m.j-1 2,63 1,16 6,33.10-02 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é. PEREZ Florian 63 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 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. ∆t (ans) -5 Fourier utilisé Maillage 64×64 Crank-Nicholson ∆t* Valeurs moyennes du Nusselt 3,63 1,45 0,73 0,29 0,15 9,11 9,1 9,09 9,08 6,31 Fréquences du Nusselt 273 280 283 286 149 Fo = 0,25 Fo = 0,10 Fo = 0,05 Fo = 0,02 Fo = 0,01 3,05x10 -5 1,22x10 -6 6,10x10 -6 2,44x10 -6 1,22x10 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. Nusselt 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 Evolution du nombre de Nusselt en fonction du temps pour un maillage de 64×64 mailles – Rayleigh 800 – Fourier 0,01. PEREZ Florian 64 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 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 EFMH CrankNicholson VF CrankNicholson Nombre Maillage de Rayleig h 90 96x96 108 96x96 200 96x96 300 96x96 400 96x96 425 96x96 470 96x96 400 96x96 425 96x96 430 96x96 450 96x96 470 96x96 400 96x96 470 96x96 Nombre de Fourier 0,4 0,4 0,4 0,4 0,4 0,4 0,4 0,4 0,4 0,4 0,4 0,4 0,2 0,2 ∆t* ∆t (ans) 2,17x10-5 2,17x10-5 2,17x10-5 2,17x10-5 2,17x10-5 2,17x10-5 2,17x10-5 2,17x10-5 2,17x10-5 2,17x10-5 2,17x10-5 2,17x10-5 1,09 x10-5 1,09 x10-5 2,58 2,58 2,58 2,58 2,58 2,58 2,58 2,58 2,58 2,58 2,58 2,58 1,29 1,29 Valeurs moyennes du nombre de Nusselt 3,22 3,6 6,49 8,16 9,16 9,38 9,85 9,16 9,35 9,39 9,53 9,87 9,16 9,87 Fréquences d’oscillation du nombre de Nusselt 121 116 - Condition initiale CIO1 CIO90 ou CIO1 CIO108 ou CIO1 CIO200 ou CIO1 CIO300 ou CIO1 CIO400 ou CIO1 CIO425 CI1 CIO400 ou CIO1 CIO425 ou CIO1 CIO430 ou CIO1 CIO450 ou CIO1 CIO1 CIO1 65 PEREZ Florian 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 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. 66 PEREZ Florian 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 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 478 ° C 359 239 120 Champ de températures en régime permanent, un seul rouleau de convection est visible. Nusselt 3.50 Régime permanent établi 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 Evolution du nombre de Nusselt en fonction du temps (adimensionné), pour Ra = 90 et une condition initiale de type CI1. PEREZ Florian 67 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 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 478 ° C 359 239 120 Champ de températures en régime permanent, deux panaches de fluides chauds montent le long des parois latérales. Nusselt 4.20 4.00 3.80 Régime permanent établi 3.60 3.40 3.20 3.00 0.80 1.00 1.20 1.40 1.60 1.80 2.00 Temps Evolution du nombre de Nusselt en fonction du temps (adimensionné), pour Ra =108 et une condition initiale de type CIO90. PEREZ Florian 68 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 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 478 ° C 359 239 120 Champ de températures en régime permanent, deux panaches de fluides chauds montent le long des parois latérales. Nusselt 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 Régime permanent établi Evolution du nombre de Nusselt en fonction du temps (adimensionné), pour Ra =200 et une condition initiale de type CIO108. PEREZ Florian 69 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 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 478 ° C 359 239 120 Champ de températures en régime permanent, deux panaches de fluides chauds montent le long des parois latérales. 8.50 8.00 Régime permanent établi 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 Evolution du nombre de Nusselt en fonction du temps (adimensionné), pour Ra =300 et une condition initiale de type CIO200. PEREZ Florian 70
Fly UP