Je tiens à exprimer mes sincères remerciements à toutes les personnes qui ont contribué à la réalisation de cette thèse de doctorat.

Il faudrait sûrement plus que quelques mots pour remercier mes encadrants, qui, durant trois ans, ont eu à la fois la motivation, le courage et la patience, de m'indiquer là ou il fallait que j'aille et là ou il valait mieux que j'arrête de m'entêter. Merci pour l'intérêt que vous avez porté à mon travail, votre expertise et vos commentaires constructifs qui m'ont permis d'améliorer la qualité de cette thèse. Tout au long de cette thèse, j'ai pu m'inspirer et apprendre d'eux, chacun m'apportant à sa manière une part de son expérience pour ma prochaine étape dans ce monde, fascinant mais parfois cruel, qu'est la recherche scientifique.

Je suis reconnaissante envers Aurélien Larcher pour son encadrement, ses conseils avisés et sa disponibilité tout au long de ces années (je pense à nos réunions dont la durée syndicale minimale était équivalente à trois pauses déjeuner françaises). Sa rigueur scientifique, sa patience et sa bienveillance ont été des atouts majeurs pour mener à bien ce projet de recherche. Il m'a transmis son envie de savoir, de comprendre et de recherche à expliquer que; derrière tout résultat; il y a une explication justifiable et plausible; qu'il est parfois possible de se tromper dans notre interprétation des choses, mais que s'arrêter à des simples observations par peur d'avoir tord ne sera jamais suffisant. C'est le genre de directeur de thèse qui ne se contente pas d'une relation nécessaire et codifiée avec ses doctorants: il cherche également à les connaître, les mettre en confiance, et à comprendre ce qu'ils ont derrière la tête. Enfin, s'il y a une dernier aspect du personnage dont je vais me souvenir, c'est bien son humanité et sensibilité. En effet, il est doté d'une grande capacité à sentir quand la mise en valeur de la part de travail bien faite vaut mieux que la critique de celle à améliorer; sans oublier sa facilité à étroitement parsemer quelques graines d'humour dans un océan de sérieux scientifique. S'il y a une qualité que j'aimerais apprendre de lui, c'est bien sa capacité à encourager ses étudiants sans cesse pour qu'ils tiennent jusqu'au bout du chemin.

Elie Hachem fut une directeur de thèse à l'emploi du temps souvent chargé et parfois quelque peu contraignant. Cependant, cela ne l'a pas empêché de se rendre disponible par vidéo ou audio conférences, ou pour m'accueillir pendant quelques petites minutes dans son bureau. Si son humour m'a laissé parfois quelque peu iii perplexe, il m'a cependant énormément marqué par sa détermination à faire fructifier la recherche au sein de l'équipe, assumant totalement et fièrement son originalité. C'est aussi quelqu'un qui m'a bien montré comment organiser mes idées afin de rendre mon manuscrit le plus claire possible tout en restant cohérente et rigoureuse.

Je remercie également les membres de mon jury de thèse; Alvaro COUTINHO, Joan BAIGES, Nissrine AKKARI ; pour leur intérêt pour mon travail, leur expertise et leurs commentaires constructifs qui ont permis d'améliorer la qualité de cette thèse. Leurs remarques ont été précieuses pour affiner mes analyses et conclusions.

Je suis reconnaissant envers mes collègues de laboratoire pour leur collaboration, leur aide et leur amitié tout au long de cette aventure. Leurs discussions scientifiques et leurs encouragements ont été des moments forts de cette expérience de recherche.

Je remercie également mes proches, ma famille et mes amis, pour leur soutiens indéfectible, leur encouragement et leur compréhension pendant ces années de travail intense. Leur présence et leur soutien ont été essentiels pour me permettre de mener à bien ce projet.

Je tiens à remercier toutes les personnes qui ont contribué des près de de loin à la réalisation de cette thèse, ainsi que les institutions qui ont soutenu financièrement ce projet. Je pense notamment aux partenaires industriels. 5.7 Set-up for the 3D industrial quenching(left) and the boundary layer mesh used (right). The solid is represented in red . . . . . . . . . . . 5.8 Sensors' position inside the cylindrical workpiece. . . . . . . . . . . . 5.9 Industrial quenching. Results at t=1.55s, 6.

... À mes parents, Il a suffit que je pense à toi, Papa, pour que mes peurs s'émiettent comme du biscuit .. Chaque fois que j'emprunte un chapitre, comme on emprunte un passage secret, je perçois une présence penchée par dessus mon épaule, je me retourne, et c'est toi, ma maman adorée .. obtained by: (a,b) by Li et al. [START_REF] Li | Lattice Boltzmann Method Simulation of 3-D Natural Convection with Double MRT Model[END_REF]; (c,d) by [START_REF] Nee | Hybrid meso-macroscopic simulation of threedimensional natural convection combined with conjugate heat transfer[END_REF], (e,f) This study . . . 4.17

List of Figures

Industrial motivation

Due to their complexity, the simulation of industrial processes allows to predict how they actually behave in a given set of operating conditions. The design of industrial processes using simulation software allows, inter alia, to increase factories efficiency and especially to reduce their environmental impact. In our work, we are especially interested in quenching process defined as the soaking of a metal at a high temperature, above the recrystallization phase, followed by a rapid cooling process to obtain a material with some desirable properties. The fields of application are multiple: the cooling of automotive, nuclear and aerospace parts (e.g landing gear or turbines for airplanes and rocket's parts). This process has a direct impact on the evolution of mechanical properties, microstructure control and release of residual stresses. In order to build an optimal process, it is important to correctly control the phase transformations that take place in the alloy, and thus obtain a microstructure with the desired thermomechanical properties.

Quenching is a highly non-linear process because of strong coupling between fluid mechanics, heat transfer at different interfaces, phase transformations of the solid and boiling of the quenching bath. Despite advances in numerical simulation, this process remain very expensive to be simulated because it is a multi-physics process involving a complicated pattern of couplings among heat transfer, phase transformations and stress evolution.

Problematic

The driving physical event is the heat transfer that takes place when a solid is soak in a fluid because it triggers other processes and involve very complex features with steep gradients. The finite element approach is adopted for the numerical treatment of these physical problems with a finer discretization of the computational domain allowing to get a better representation of the approximated solution. However, this approach tends to drastically increase the number of elements in the mesh to well capture anisotropic features and consequently leads to a considerable increase of the computational cost. Nowadays, the aim is to simulate these relatively complex problems while obtaining sufficiently accurate results and by limiting computational cost. The idea is to discretize the computational domain in a finer way in the regions where it is required. The domain is refined according to the evolution of the different interfaces: water/vapour/solid. The notion of mesh refinement was introduced by [22] and improved later by [23]. However, this refinement of the mesh has a price and the obtained meshes can become very expensive in terms of number of nodes.

The more the number of nodes is important, the more the computational time is longer.

Ambition of this work and the state of the art

Anisotropic meshes adaptation seems to be an interesting candidate to deal with the solution accuracy, they are formed of stretched and oriented elements in order to minimize the interpolation error of the domain while respecting a number of nodes. These particular meshes allow to get solutions whose precision is equivalent to those obtained with conventional isotropic meshes and often contain much less nodes. These savings in terms of nodes and elements are not negligible especially with the solution computational time and the required memory space. These types of meshes were first introduced in [5,6] and developed by several other authors [1], [2], [3], [4].These methods rely on a posteriori error estimates, to identify the regions of interest to be adapted and provide the best suitable local orientation and stretching of the elements. The new mesh is obtained applying local modifications to the previous one using a metric field resulting from an error analysis. The idea of using an error estimator to drive the mesh adaptation procedure is proposed for the first time by [24], [25]. Despite the level of accuracy obtained with anisotropic mesh adaptation, it turns out that when we deal with large computational domain and long time simulations, the mesh resolution would be more complex inducing a considerable computational cost. The comprehension of real world phenomena and industrial processes is today largely based on partial differential equations (PDEs) allowing to model the different physical and mechanical processes involved. The numerical resolution of any PDE most often requires an iterative algorithm. In general, the discretization of PDEs gives rise to large-scale systems. As the resolution of large systems is very expensive in terms of computational cost, an important question arises. In order to obtain a good quality approximated solution, when to stop linear solver resolution to avoid unnecessary iterations ? This key point has been addressed in several publications, we refer, inter alia, to Chaillou A. L., Suri M. [19] et L. El Aloui , Ern A., M. Vohralík, [21].

It turns out that for a given simulation, up to 80 % of the time is spent to solve iteratively the linear system that stems from the discretization of the governing PDE. Therefore, to increase efficiency and performance we need to provide an optimal control of the convergence of the iterative algorithm used to reduce its number of iterations while keeping a good precision of the numerical method. The question of computation time saving and numerical method accuracy in the framework of finite elements has been addressed in [7,8] where an a priori error estimator is used as a stopping criterion to drive the convergence of the iterative algorithm. More recent contributions based on a posteriori error estimators has been considered in [9,10]. In fact, this approach provide stopping criterion allowing a global control of the error. This error is not controlled according to the regularity of the exact solution as in the case of a priori estimators but using of the discrete solution that we can compute.

The a posteriori error estimators are not only basic tool to build an adaptive mesh, because they also allow to optimize the convergence of iterative algorithm by avoiding unnecessary iterations. The a posteiori approach meets several objectives. First of all, it globally control the discretization error of the considered problem. Furthermore, it allows to use computer resources in a suitable and efficient way. Indeed, the advantage of the a posteriori approach is to provide explicit bounds on the error between the numerical solution and the exact one as soon as the approached solution is know. This analysis was initiated by Babuška [11] for parabolic problems and developed by Verfürth [12,13]. We also quote Ainsworth and Oden [14], Zienkiewicz and Zhu [15,16], Strouboulis and Hague [17] and Repin [18]. The a posteriori analysis can provide stopping criterion that ensure a global control of the error. A posteriori error estimators have been designed by Chaillou and Suri [19,20] allowing to distinguish linearization and discretization errors for non linear problems and have been developed in the framework of iterative algorithms by L. El Alaoui, A. Ern and M. Vohralík [21].

Objective of the thesis

In this work, we are interested in improving the adaptive stopping criterion based on the a posteriori error estimator available from the mesh adaptation procedure developed in the works of G. Manzinalli [26]. In order to fulfil this challenge, a mathematical analysis has been addressed to estimate the scaling coefficient involved within the adaptive criterion previously determined through a parametric study. The obtained stopping criterion will be compared with the previous one through validation tests, it has to be performant and efficient for the different PDEs allowing to describe the complex physical phenomena involved in quenching (heat transfers and fluid dynamics). In addition, we implemented different gradient recovery approaches [28], [29], [30] allowing the computation of the edge-based a posteriori error estimator and we compared them to the already implemented technique developed in [27], it is based on a Least Square minimization technique. The most optimal interpolation operator will finally be used to define a performant stopping criterion allowing to obtain an accurate representation of the solution variations with the least computational cost. We finally evaluate the validity, accuracy and efficiency of the proposed framework through 2D and 3D analytical validation tests and for industrial numerical simulations including workpieces' quenching processes.

In the following chapter, we will explore several mathematical models and numerical methods presented in the literature allowing to simulate the underlying process. To model the fluid dynamics, the incompressible Navier-Stokes equation will be considered and the convection-diffusion-reaction problem for heat transfer. The stabilized finite element method implemented in CIMLIB -CFD will be used in order to yield accurate and oscillation free numerical solution. The CIMLIB-CFD is a finite element library developed at CEMEF by the research team of T. Coupez [31,32].

Résumé du chapitre en français

Ce premier chapitre propose une introduction générale des enjeux industriels et numériques du procédé de trempe industrielle. La trempe est une technique de traitement thermique qui consiste à plonger un matériau porté à haute température dans un fluide de quelques dizaines de degrés Celsius afin de le doter de propriétés nouvelles. Les domaines d'application de cette opération métallurgique sont multiples: le refroidissement des pièces automobiles, nucléaires et aérospatiales (e.x trains d'atterrissage ou turbines d'avions et pièces de fusée). La trempe a un impact direct sur l'évolution des propriétés mécaniques, le contrôle de la microstructure et la libération des contraintes résiduelles. Afin de réaliser un processus optimal, il est important de contrôler correctement les transformations de phase qui ont lieu dans l'alliage, et ainsi obtenir une mictrostructure avec les propriétés thermomécaniques souhaitées. La trempe est un processus hautement non linéaire en raison du couplage fort entre la mécanique des fluides, le transfert de chaleur à différentes interfaces, les transformations de phase du solide et l'ébullition du bain de trempe. Malgré l'amélioration continue des méthodes numériques d'une part et des puissances de calcul informatique d'autre part, des difficultés persistent quand au coût de simulation du procédé de trempe en raison de la grandeur des installations (dizaines de mètres), la complexité des phénomènes physiques traités et la longue durée du processus (dizaines d'heures). Nous avons adopté une discrétisation du domaine de calcul par éléments finis combinée avec une méthode d'adaptation de maillage anisotrope qui visent à améliorer les précisions de calcul d'une part et réduire le temps de calcul d'autre part. L'idée consiste à discrétiser le domaine de calcul d'une manière plus fine dans les régions où c'est nécessaire en particulier dans les voisinages des forts gradients. En effet, le domaine est raffiné en fonction de l'évolution des différentes interfaces: eau/vapeur/solide. L'objectif scientifique de ce projet de thèse est d'améliorer le critère d'arrêt adaptatif basé sur l'estimateur d'erreur a posteriori, initialement utilisé pour conduire le processus d'adaptation de maillage anisotrope, développé dans les travaux de G. Manzinalli [26]. Afin de relever ce défi, une analyse mathématiques a été effectué pour estimer le coefficient de scaling préalablement déterminé grâce à une étude paramétrique. Nous avons implémenté ce scaling dans notre bibliothèque Cimlib-CFD à l'aide de la bibliothèque externe PETSc. Le critère d'arrêt ainsi obtenu a été comparé au précédent à travers différents benchmarks tirés de la littérature. Nous avons également implémenté différentes techniques de reconstruction du gradient [28], [29], [30] permettant de calculer l'estimateur d'erreur a posteriori le long des arêtes du maillage. L'implémentation a été effectué dans notre bibliothèque Cimlib-CFD. L'objectif étant de simuler des applications 3D réelles et complexes dans de grands espaces et pour de longues durées, les calculs doivent alors être lancés sur plusieurs processeurs. Ainsi, nous avons décrit l'implémentation parallèle de ces méthodes de reconstructions. Le cadre ainsi développé a été validé à travers plusieurs exemples.

Introduction

The present chapter focuses on the model equations and numerical methods allowing to achieve an accurate and a realistic simulation of the quenching process. This requires the knowledge and the understanding of the involved multi-physical phenomena. In addition, the interactions between these phenomena are difficult to be determined, because they happen depending on each other.

During heat treatment, temperature and cooling speed are essential factors to make metallic solids respecting the desired mechanical and microstructure properties. These factors determine most of material final characteristics. Numerical simulation is an excellent tool to monitor the evolution of the temperature of the metal part, the heat exchanges between the quenching liquid and the metal parts and the agitation speed of the quenching liquid. Therefore, it is essential to model and resolve the heat transfer and the fluid flow. The computations consist in simultaneously solving the Navier-Stokes and the Convection Diffusion Reaction equation. Stabilized finite element method (SUPG) implemented in CIMLIB-CFD library is used to discretize the governing equation. The library is developed by CFL research group in CEMEF Mines Paris. The proposed finite element solvers have been validated in several previous works. For that purpose, we refer to [1,4] for more details.

After a brief review of equations relating to thermal (heat equation) and hydrodynamics (Navier-Stokes equation) aspects, we then describe numerical and stabilized techniques, in case of convection dominated regimes, used to solve these equations.

Convection-diffusion reaction equation

The Convection Diffusion Reaction (CDR) equation describes the concentration u of a certain substance and its variations in space and time resulting from diffusion, convection and reaction processes. In fact, the heat equation requires only convection and diffusion, but for sake of generality we present the full CDR equation

The solution of the CDR equation covers a large spectrum of fluid mechanics applications in particular the heat transfer through a permeable medium or the transportation of a pollutant through the interface. Historically, the Galerkin Finite Element method is the first method applied to solve this type of equation.

In this section, we will start by recalling the convection-diffusion reaction equation. After that, we will introduce standard Galerkin finite element method. And we will finish with pointing out the need of stabilized finite element method in case of convection dominated regimes. We will also provide the formulation of the so-called Streamline Upwind Petrov-Galerkin (SUPG) stabilization technique.

Governing equation

Let Ω ∈ R d , d = 2, 3 a polygonal domain (polyhedral for d = 3) open and bounded, ∂Ω its boundary and [0, T ] a time interval. We consider the unsteady convectiondiffusion-reaction equation modelling the transport of a quantity u:

   ∂ t u + v • ∇u -κ∆u + σu = f in Ω × [0, T ], u(., 0) = u 0 in Ω, u = g on ∂Ω, (2.1)
where

v ∈ [W 1,∞ (Ω)] d is the divergence-free velocity field, κ ∈ L ∞ (Ω) is the diffusion coefficient, σ ∈ L 2 (Ω) is a reaction coefficient, u 0 is the initial condition,

Standard Galerkin finite element formulation

In order to give the weak formulation of the Equation (2.1), let us define some Hilbert-Sobolev spaces L 2 and H 1 , where L 2 is the Lebesgue space of square integrable functions given by:

L 2 (Ω) = u ∈ Ω : Ω |u(x)| 2 dx < ∞
and H 1 the Sobolev space of square integrable first order derivatives functions defined by:

H 1 (Ω) = {v ∈ L 2 (Ω) : ∥∇v∥ ∈ L 2 (Ω)} . (2.2)
We also define the subspace H 1 0 (Ω) ∈ H 1 (Ω) the set of functions vanishing on the boundary. The Hilbert-Sobolev space H 1 is a natural choice to "measure" functions involved in the weak formulation of PDEs as the existence of the integrals relies on the fact that integral of power 2 of u and weak derivative ∇u exists. The Galerkin variational formulation of the convection-diffusion-reaction equation is obtained by multiplying Equation (2.1) by a test function w ∈ H 1 0 (Ω) and integrating over the computational domain Ω. Thus the weak formulation associated to the convection diffusion reaction is written:

   Find u ∈ H 1 0 (Ω) such that: (∂ t u, w) + a(u, w) -l(w) = 0 ∀w ∈ H 1 0 (Ω) (2.3) where    a(u, w) = (v • ∇u, w) + (κ∇u, ∇w) + (σu, w), l(w) = (f, w) (2.4)
We define by (u, w) = Ω u(x)v(x) the L 2 scalar product. Now, let us consider a finite element subdivision T h of the continuous domain Ω into simplex elements K. Using this partition, the functional spaces are now approximated by discrete spaces called H 1 0,h (Ω) and H 1 h (Ω) Therefore, the discrete weak formulation or the Galerkin finite element formulation of convection diffusion reaction problem reads:

   Find u h ∈ H 1 0,h (Ω) such that : (∂ t u h , w h ) + a(u h , w h ) -l(w h ) = 0 ∀w h ∈ H 1 0,h (Ω) (2.5)
The Standard Galerkin formulation is unstable and gives birth to oscillations when the convective term of the equation is dominant. Therefore, we stabilize the formulation using the Streamline Upwind Petrov-Galerkin (SUPG) numerical scheme [35]. This approach has proved its efficiency in eliminating spurious oscillations related to the Galerkin formulation.

Streamline Upwind Petrov Galerkin method (SUPG)

The SUPG method is formulated by adding an additional weight τ K u • ∇w h to the standard Galerkin test function w h in the upwind direction for all terms in the equation. This modification in interpreted as adding artificial diffusion which acts only in the direction of the flow. In other words, we add more weight to the nodes in the upstream direction and reduces the weight of the nodes in the downstream direction.

The new test functions are defined by

wh = τ K u • ∇w h (2.6)
τ K is a piecewise-defined stabilization coefficient. For the details on the choice of stabilization parameter τ K we refer to [14] and references therein. The stabilized variational formulation of Equation 2.5 reads:

   Find u h ∈ H 1 0,h (Ω) such that: (∂ t u h , w h ) + a τ (u h , w h ) -l τ (w h ) = 0 ∀w h ∈ H 1 0,h (Ω) (2.7)
where

       a τ (u h , w h ) = a(u h , w h ) + K∈τ h τ K (v • ∇u h -κ∆u h + σu h , u • ∇w h ) K l τ (w h ) = l(w h ) + K∈τ h τ K (f, u • ∇w h ) K (2.8)
The SUPG method used in this work is globally stable, and is a good compromise between stability and accuracy of the solution.

The incompressible Navier-Stokes equations

In fluid mechanics, the Navier-Stokes equations are non-linear partial differential equations that describe fluids motion as continuous media. For example, they govern the air motion in the atmosphere, ocean currents, flow water through a pipe, and many others fluid flow phenomena. These equations were introduced in the 19th century and were named by two physicists, Claude Navier and George Stokes. Their resolution allows to obtain a velocity/pressure pair (v, p) the solution of the considered problem.

In this section, we will give an overview of the mathematical expression of the governing equations of incompressible fluid flows. The incompressibility condition is introduced when the fluid density changes are ignored. We set the problem in a stabilized finite element framework to deal with the instabilities. To get a compatibility between velocity and pressure approximation spaces, the inf-sup condition must be verified, its also known as the Brezzi-Babuska condition, for more information refer to [5,6]. Several methods to deal with instabilities problems of Stokes and Navier-Stokes problems can be found in the literature. Among the various existing methods, we have: 1 The Unusual Stabilized finite element method introduced by Franca and Farhat in [7], 2 The residual based stabilized methods are well described by several authors such as Brezzi and al in [8] and Codina et al [9][10][11]. In this work, The variational Multiscale (VMS) method will be used to deal with instabilities. It was implemented and validated in CIMLIB by E. Hachem [1]. Its purpose is to deal with convection dominated regimes problems by satisfying the inf-sup Brezzi-Babuška condition.

Governing equations

We are interested in a fluid's evolution over time through the study of its velocity field in every point of space and at every moment: the chosen formulation is then an Eulerian formulation, contrary to a Lagrangian formulation where one would rather seek to describe the position of fluid's particles at every moment.

We assume that the fluid is of positive constant density ρ, incompressible (the space occupied by a certain amount of fluid at each moment may change shape, but not of volume) and we will assume that it is viscous.

Let us now describe the considered equations: let Ω ∈ R d a bounded domain of dimension d ∈ {2, 3} and [0, T ] the time interval. The fluid's incompressibility is translated into Eulerian language by a free-divergence velocity field at every time step. The incompressible Navier-Stokes equations in strong form are written:

   ρ(∂ t v + (v • ∇)v) -∇ • σ = f in Ω × (0, T ), ∇ • v = 0, in Ω × (0, T ), (2.9)
where f represents the external forces acting on the flow and σ is the Cauchy stress tensor for a Newtonian fluid, it is defined by:

σ = 2µϵ(v) -pI d (2.10)
where µ and ρ are the viscosity and the density, p the pressure and ϵ(v) is the strain rate tensor, defining the symmetric part of the velocity gradient:

ϵ(v) = 1 2 [∇v + (∇v) T ] (2.11)
Using the definitions of stress tensor σ and the strain rate ϵ, Equation 2.9 can be rewritten as follows

   ρ(∂ t v + v • ∇v) -2µ∇ • ϵ(v) + ∇p = f in Ω × (0, T ), ∇ • v = 0 in Ω × (0, T ), (2.12)
To close the system, we use appropriate initial and boundary conditions

   v(x, 0) = v 0 (x) in Ω v(., t) = v D in ∂Ω D × (0, T ) ∇v • n = ϕ N in ∂Ω H × (0, T ) (2.13)
where ∂Ω = ∂Ω D ∪ ∂Ω H and ∂Ω D ∩ ∂Ω H = ∅.

The Galerkin formulation

To get the variational formulation corresponding to the Navier-Stokes Equation 2.9, we introduce some functional spaces

V = v ∈ (H 1 (Ω)) d | v = v d on ∂Ω D V 0 = v ∈ (H 1 (Ω)) d | v = 0 on ∂Ω D Q = q ∈ L 2 (Ω)
They define respectively the the function spaces for the velocity, the scalar function fo the pressure and the weighted function space for the velocity: We obtain the Galerkin variational formulation of the problem (2.9) by multiplying it by a test function (w, q) ∈ (V 0 , Q). Thus the weak formulation associated to the incompressible Navier-Stokes equation reads:

                   Find (v, p) ∈ (V, Q) such that: ρ(∂ t v + v • ∇v, w) Ω + (2µϵ(v), ϵ(w)) Ω -(p, ∇ • w) Ω = (f , w) Ω + (ϕ N , w) ∂Ω N ∀w ∈ V 0 (∇ • v, q) Ω = 0 ∀q ∈ Q (2.14)
For the spatial discretization we consider the finite element admissible partition T h of the domain Ω into a family simplex elements {K}. The previously defined functional spaces V and V 0 are approached by finite dimensional discrete spaces generated by piecewise linear functions defined bellow

V h = w h ∈ V h | v h | K ∈ P 1 (K) d , ∀K ∈ T h V h,0 = v h ∈ (C 0 (Ω)) d | w h | ∂Ω = 0 Q h = q h ∈ C 0 (Ω) | q h | K ∈ P 1 (K), ∀K ∈ T h
Thus, Galerkin formulation associated to (2.14) reads:

                   Find (v h , p h ) ∈ (V h , Q h ) such that : ρ(∂ t v h + (v h • ∇)v h , w h ) Ω + (2µϵ(v h ), ϵ(w h )) Ω -(p h , ∇ • w h ) Ω = (f , w h ) + (ϕ N , w h ) ∂Ω N ∀w h ∈ V h,0 (∇ • v h , q h ) Ω = 0 ∀q h ∈ Q h (2.15)
The use of the standrad Galerkin finite element method to simulate the incompressible Navier-Stokes equations at higher Reynolds number is challenging and often failed due to the domination of convective term, it leads to a loss of coercivity of the discrete operator associated to the problem (2.15). This issue can be overcomed by using a stabilization technique. As said previously, we use the Variational Multiscale method (VMS), it was originally introduced and developed by Hughes et al [15,16] . The main idea is to rewrite the governing equations as a coupled system with two types of scales: coarse and fine. First, we solve the fine scales in an approximate manner and then we replace their effects into the large-scale equation.

Variational Multiscale method(VMS)

The VMS approach provides a natural stabilization by an orthogonal decomposition of the solution spaces (velocity and pressure). We present here only an outline of the method, and the reader is referred to [32] for extensive details about the formulation. The partition of functional spaces is performed as follows:

Ṽ = V h ⊕ V ′ and Q = Q h ⊕ Q ′
The unknowns (velocity and pressure) can be split into two components, a resolvable coarse-scale and unresolved fine-scale as follows:

ṽ = v h + v ′ p = p h + p ′
The fine-scales provide additional stabilization in the regions where the velocity gradients are strong.

In the same way, we obtain the following test functions:

w = w h + w ′ q = q h + q ′
Subscript h is used to denote the finite element coarse component, whereas the prime is used for the so called fine component of the unknown. Consequently, the mixed-finite element weak formulation associated to problem (2.15) reads

                          
Find (ṽ, p) ∈ ( Ṽ , Q) such that:

(ρ∂ t (v h + v ′ ) + ρ((v h + v ′ ) • ∇)(v h + v ′ ), (w h + w ′ )) + (2µϵ((v h + v ′ )), ϵ((w h + w ′ )) -((p h + p ′ ), ∇ • (w h + w ′ )) = (f , (w h + w ′ )) + (ϕ N , (w h + w ′ )) ∂Ω N (∇ • (v h + v ′ ), (q h + q ′ )) = 0
for any pair (w h , q h ) ∈ ( Ṽ , Q)

(2.16) Thanks to the linearity of test function, we split equation 2.16 into sub-problems, namely the coarse-scale problem and the fine-scale problem

                   Coarse-Scale problem (ρ∂ t (v h + v ′ ) + ρ(v h + v ′ ) • ∇)(v h + v ′ ), w h ) + (2µϵ((v h + v ′ )), ϵ(w h ) -((p h + p ′ ), ∇ • w h ) = (f, w h )) + (h N , w h ) ∂Ω N ∀w h ∈ V h,0 (∇ • (v h + v ′ ), q h ) = 0 ∀q h ∈ Q h,0 (2.17)
                  
Fine-scale problem

(ρ∂ t (v h + v ′ ) + ρ(v h + v ′ ) • ∇(v h + v ′ ), w ′ ) + (2µϵ((v h + v ′ )), ϵ(w ′ ) -((p h + p ′ ), ∇ • w ′ ) = (f, w ′ ) + (h N , w ′ ) ∂Ω N ∀w ′ ∈ V ′ 0 (∇ • (v h + v ′ ), q ′ ) = 0 ∀q ′ ∈ Q ′ 0 (2.18)
The fine scales acts as a local stabilization on coarse-scale problem. Thus, we begin by approximating the fine scale problem and then, we substitute the fine scale solution back into the coarse problem. To do so, let start by rearrange the terms of the fine scale formulation:

           ρ(∂ t v ′ , w ′ ) + (ρ(v h + v ′ ) • ∇v ′ , w ′ ) + (2µϵ(v ′ ), ϵ(w ′ )) + (∇p ′ , w ′ ) = (R m , w ′ ) ∀w ′ ∈ V ′ 0 (∇ • v ′ , q ′ ) = (R c , q ′ ) ∀q ′ ∈ Q 0 (2.19)
where R m , and R c are respectively the residuals of the momentum and continuity large-scale problem projected onto fine-scale spaces

R m = f -ρ∂ t v h -(ρv) h • ∇v h -∇p h (2.20) R c = -∇ • v h (2.21)
In order to obtain an approximate solution of the small-scale problem, we will introduce some simplifying assumptions to deal with the time-dependency and non-linear terms appearing in Equation (2.19).

• We will neglect the contribution of the time dependent term. In fact, tracking time evolution of fine-scales causes a huge computational cost especially when we deal with complex 3D problems, check [17] for further information about the legitimacy of such a choice.

Considering these assumptions true, the fine-scale problem (2.19) becomes

   (ρv h • ∇v ′ , w ′ ) K + (2µϵ(v ′ ), ϵ(w ′ )) K + (∇p ′ , w ′ ) K = (R m , w ′ ) K ∀w ′ ∈ V ′ 0 (∇ • v ′ , q ′ ) K = (R c , q ′ ) K ∀q ′ ∈ Q 0 (2.22
) However, it is hard to get both velocity and pressure associated to the smallscale problem. The strategy is to account their contribution in the coarse-problem by approximating them within each element K ∈ T h using stabilization parameters and residuals terms R m and R c as follows

v ′ = τ v R m p ′ = τ p R c (2.23)
Several definitions are given in the literature for the stabilization coefficient, for instance Codina in [9] proposed the following stabilization parameters:

τ v | K = c 1 µ ρh 2 m 2 + c 2 ρ∥v h ∥ K h K 2 -1 2
(2.24)

τ p | K = µ ρ 2 + c 2 ρ∥v h ∥ K c 1 h K 2 -1 2
(2.25)

where h K the characteristic length of the element K, ∥v h ∥ K the norm of v h in K and c 1 , c 2 two constants independent from h K . We take c 1 = 4 and c 2 = 2 for linear elements (see [18]). We can introduce the time step size of the temporal disretization in the definition of τ v to improve convergence of the algorithm and to deal with the non-linearity of the problem:

τ v|K = c 1 µ ρh 2 m 2 + c 2 ρ||v h || K h K 2 + 2ρ ∆t 2 -1 2
(2.26)

Now, we have all the necessary ingredients to solve the coarse problem 2.17. Taking into account the assumptions on the fine-scale problem, and integrating by parts certain terms of the fine-scale equation and replacing the expressions of both v ′ and p ′ , the coarse scale equation becomes:

         (ρ∂ t v h , w h ) + (ρ(v h • ∇)v h , w h ) - K∈τ h (τ v R m , ρv h • ∇w h ) K + (2µϵ(v h , ϵ(w h ))) -(p h , ∇ • w h ) + K∈τ h (τ p R c , ∇ • w h ) = (f, w h ), ∀w h ∈ V h,0 (∇ • v h , q h ) - K∈τ h (τ v R m , ∇q h ) = 0 ∀q h ∈ Q 0 (2.27)
Finally, substituting R m and R c by their explicit expressions,the coarse-scale problem can be expressed using coarse scale functions

                                                               (ρ∂ t v h + ρv h • ∇v h , w h ) + (2µϵ(v h , ϵ(w h ))) -(p h , ∇ • w h ) + (∇ • v h , q h ) -(f, w h ) Standard Galerkin terms + K∈τ h τ v (ρ(∂ t v h + v h • ∇v h ) + ∇p h -f, ρv h • ∇w h ) K fine-scale velocity upwind stabilization terms + K∈τ h τ v (ρ(∂ t v h + v h • ∇v h ) + ∇p h -f, ∇q h ) K fine-scale pressure stabilization terms + K∈τ h τ p (∇ • v h , ∇ • w h ) K = 0 grad-div stabilization term ∀w h ∈ V h,0 , q ∈ Q 0 (2.28)
This formulation can be subdivided into 4 parts: the standard Galerkin formulation, the second and third are fine-scale velocity contributions and the last part models the effects of small-scale pressure. In convection dominated regimes, the added terms provide additional stabilization and allow to circumvent the inf-sup condition known as the Babuška-Brezzi condition [5,6].

It is important to mention that the element length h K influence significantly the behaviour of the stabilization parameter. For each element mesh K, the expression of the stabilization coefficient depends on this parameter. Since we are working on convection dominated problems, we used the definition proposed in [31] to compute h K . The idea is to compute h K as the diameter of the element K in the direction of the velocity field u h to deal with any mesh type, especially strongly anisotropic meshes. The expression of h K reads

h K = 2∥u h ∥ K ϕ i |(u h • ∇)ϕ i (x K )| (2.29)
with P 1 (K) = span{φ i } and x K the cell centroid.

Overview on Turbulence modelling

The majority of daily observed "natural" flows are turbulent (cigarette smoke, milk cream poured into the coffee, clouds). They are furthermore very diversified: biological flows (blood), Geofluids motion (wind, sea current), motion of stellar fluids (gas circulation around planets). Despite this, these flows have common properties which we will set out later. The equations that govern the fluids motion, whether turbulent or not, are The Navier-Stokes equations. These equations had been previously written by Euler but Navier has had the ingenuity to add a friction term between the diverse fluid layers. The turbulence became an experimental science toward the end of the 19th century when the Englishman Osborne Reynolds was able to observe the transition from the laminar regime to the turbulent one. He thus highlighted some rather simple laws and introduced a dimensionless number bearing its name to characterize this transition. Nevertheless, before the 1950s, turbulence was an obscure subject. The only way for the engineer was to experiment on physical models in order to improve his know-how. Fortunately, after the 1960s, the situation was about to unravel with some accomplished progresses in modelling, at the same time as the capacity of numerical treatments was increasing strongly. However, the predetermination of local statistical properties was still impossible. Currently, numerical computations allow a good estimation of the average flow in the presence of a developed turbulence. Large numbers of models have been developed and studied in the last three decades. We can classify these modelling approaches into three categories: DNS, LES and RANS. The Direct Numerical Simulation (DNS) attempts to simulate and resolve all the scales of motion without approximation or need of any additional modelling. Since this approach aims to solve all spatial and temporal gradients, its application requires huge computational cost. Unlike DNS, Large Eddy simulations (LES) are based on filtering and decomposing the scales into large-scale and small-scale. Its aim is to solve and simulate more precisely large-scales since they are more energetic than small-sales and are the main contributors to the transport of conserved quantities. Consequently, they are much less costly than DNS simulations. The Reynolds-averaged Navier-Stokes equation (RANS) determines flow parameters, such as water velocity and temperature using some turbulence modelling. The idea here is that all turbulent motions will be modelled. As turbulent fluctuations are removed, coarser grids can be used , consequently two and three-dimensional simulations can be performed with significantly less computer resources compared to DNS and LES. In Figure (2.1), we compare the various turbulence approaches mentioned above. However, while small-scale eddies are note solved, some additional models are required to approximate the turbulence effects. One of the most widely used model is two-equations k -ϵ model. More details can be found in [33], [34], [33].

Water-cooling of a rectangular ingot by natural convection

To assess the validity, accuracy and efficiency of the proposed framework of stabilization techniques applied to the Navier-Stokes equations for velocity coupled with convection diffusion reaction equation for temperature, we consider the watercooling by natural convection of a heated workpiece. The numerical simulations were carried out using C++ CIMLIB finite element library.

The problem studies the heat transfer by conjugate natural convection involved in 18 seconds of cooling process of an Aluminium rectangular shape ingot initially heated to 110 • and embedded into a water filled cubical cavity. This test case was carried out in the framework of the ANR HECO [36] project, the tests were carried out at the LTeN (Laboratoire de Thermique et Energie de Nantes).

The interest from an industrial perspective is to understand the characteristics of the flow and ingot's thermal history in the view of benefiting to the maximum from this naturally occurring quenching process without resorting to cooling devices. The numerical analysis could allow for optimizing the heat transfer effects in terms of enclosure's volume and to determine the needed spacing between ingots when several workpieces are to be cooled in the same times. The computational domain, the initial mesh and the immersed solid geometry are give in Figure 2.3. The water inside the enclosure is initially at rest as a temperature of 25 • . We note that the thermophysical properties of the water and the workpiece are updated dynamically in accordance with their temperatures. We compare in Figure 2.2 the temperature variations at sensors 1, 2 and 3 respectively located at 2.5mm, 4.5mm and 6.5mm of the ingot surface. Despite no consideration of boiling, resulting in a slower cooling during the first few seconds, the obtained temperature evolution are in very good agreement with experimental results because the cooling rate is almost identical. We point out that the temperature profiles are not polluted with spurious oscillations (undershoots or overshoots) that are frequently observed in the presence of high temperature gradients at the interface or in a convection dominated problems. This can be attributed to the stabilized finite element discretization applied to the coupled system of fluid flow and heat transfer. Summing up, the combination of local mesh adaptation and stabilized iterative solvers with the smoothed distribution of the thermo-physical properties across the interface enables the simulation to overcome the numerical instabilities and lead to a good numerical behaviour. Natural convection in an enclosed container has received a great deal of attention by research community due to its importance in many engineering devices.

Conclusion

Despite developments in computational fluid dynamics, the resolution of equations described above poses scientific complexities, their resolution requires considerable computational resources and long lasting calculations. To tackle this problem, we will see in the next chapter how developments in mesh adaptation techniques allowed to devise realistic configurations while limiting the shortcomings of lack of solution's accuracy and computational efficiency. In fact, the resolution of these unsteady equations on a uniform mesh with a limited number of degrees of freedom often fails to capture fine scales, has drastic computational cost and might give incorrect solutions. We aim in this chapter to present the anisotropic mesh adaptation technique developed within our CFL research team based on an error estimator computed on the edges of the mesh, we will also present the development that we implemented to improve it.

Résumé du chapitre en français

Introduction

A posteriori error estimators and mesh adaptation procedures became over the last decades, one of the most important development directions in the field of numerical analysis and an indispensable tool for the discretization of partial differential equations. Their purpose is, on the one hand, to measure the quality of the approximate solution obtained and, on the other hand, to reduce the computational time by building a fine mesh at the desired location and coarser elsewhere.

The aim of a posteriori error indicators is to measure the quality of the approximated solution by giving upper bounds of the approximation error defined as the difference between the exact solution which is unknown and the numerical approximation in a certain norm. These bounds are computable once the approximated solution is known. A local computation of these indicators, for example on every mesh element or over edges, will allow to predict the error localization. Thanks to the latter, the resolution, and in an equal way, adaptive mesh refinement, can be concentrated in this region of computation.

Error indicators studied through this thesis are based on the reconstruction of derivatives, gradient and Hessian, at the mesh vertices. Indeed, as we use Lagrange finite elements of the first order P 1 , the numerical solution u h is only piecewise linear, the gradient ∇u h is constant on each cell K ∈ T h , and ∆ h is zero on each cell K ∈ T h . In order to reconstruct derivatives, we will recall, analyze and compare some non-smooth interpolation operators techniques based on different strategies (L 2 projection, Least-Square minimization, ...). These methods aim at recovering the solution's derivatives from the piecewise linear approximate solution u h defined at the vertices of the mesh.

Gradient recovery techniques

The main objective of interpolation operators is to give a "good" estimation of a piecewise continuous and derivable function u h defined on a mesh. A lot of researchers have published articles that allow to give a solution to the irregularity problem of ∇u h . Their methods for derivatives reconstruction are also called derivatives recovery methods.

Our main objective in this work will be to describe the different methods of derivatives reconstruction present in the literature. Those methods represent an important aspect of adaptation methods. In fact, in general, we do not know the solution of a finite element problem but its necessary to know the Hessian matrix to construct a metric or just first derivatives which intervene in the adaptation procedure. The methods that we present bellow allow to obtain the same convergence order for u h and ∇u h . We note by V (T h ) the set of mesh vertices. We associate to every node V i a "patch" defined as the set of simplices containing V i .

Hereafter some examples of patches for linear finite elements (left) and quadratic ones (right) Elements patches 3 node elements (Linear)

6 node elements (quadratic) △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △ △

Clément interpolation operator: L 2 projection method

Let us assume in this section that the polynomial approximation degree is equal to

q = 1. The function Π h u is in the approximation space V 1 h . We denote by (ϕ i ) 1≤i≤N the basis functions of V 1 h , so that Π h u is written Π h u(X) = n i=1 u(V i )ϕ i (X) (3.3)
Let us denote by K a triangle and [V 1 , V 2 , V 3 ] its three vertices. For X ∈ K, we have

∇Π h u(X) = Π h u(V 1 )∇ϕ 1 (X) + Π h u(V 2 )∇ϕ 2 (X) + Π h u(V 3 )∇ϕ 3 (X) (3.4)
In 2D, if (x i , y i ) are the coordinates of the node V i , the basis functions gradients are given by ∂ϕ

1 ∂x = y 2 -y 3 2|K| ∂ϕ 1 ∂y = x 3 -x 2 2|K| ∂ϕ 2 ∂x = y 3 -y 1 2|K| ∂ϕ 2 ∂y = x 1 -x 3 2|K| ∂ϕ 3 ∂x = y 1 -y 2 2|K| ∂ϕ 3 ∂y = x 2 -x 1 2|K|
We deduce that ∇Π h u is piecewise constant and discontinuous, it means that ∇Π h u ∈ V 0 h . We denote by (∇Π h u) |K the ∇Π h u gradient which is constant on the element K. However, ∇Π h u is not defined at the mesh vertices. We propose to reconstruct the nodal gradients from the piecewise constant representation of ∇Π h u by using an L 2 projection. The L 2 projection operator is based on the Clément operator [1].

The idea is to project v ∈ L 2 into the set of piecewise constant functions V 0 h , independently on each patch S i ⊂ T h . More precisely, for v ∈ L 2 (Ω), we define Π 0 v ∈ V 0 h by:

∀S i ∈ Ω h ,    (Π 0 v)| S i is constant S i (Π 0 v -v)w = 0, ∀w constant (3.5)
The Clément interpolation operator is defined by Π c :

V 0 h -→ V 1 h by: Π c v := n i=1 Π 0 v(V i )ϕ i (3.6)
Since u h is piecewise linear, we apply the Clément operator to v = ∇u h which is piecewise constant (i.e ∈ V 0 h ). For every support S i ∈ Ω h , we have for the particular choice w = 1:

S i (Π 0 (∇u h ) -∇u h )dx = 0 equivalently S i Π 0 (∇u h )dx = S i ∇u h dx Then |S i |Π 0 (∇u h )| S i = K j ∈S i K j ∇u h dx Finally: Π 0 (∇u h )| S i = K j ∈S i |K j |(∇u h )| K j |S i |
where |K| and |S i | are respectively the area (volume) of the element K and the support S i , i.e.

|S i | = j,K j ⊂S i |K j |
We thus obtain, for each vertex V i of the mesh, the following nodal gradient recovery:

∇ R u h (V i ) = K j ∈S i |K j |(∇u h )| K j |S i |
This reconstruction is nothing other than the average of the gradients of the patch elements S i weighted by the volume of the patch elements. This corresponds to the notion of "P 1 -Galerkin nodal gradient in finite elements". Hence, we obtain a piecewise P 1 representation of the gradient. The same procedure can be applied, once again, to each component of the gradient to recover the Hessian.

Clément interpolation operator verifies some regularity properties allowing to justify its use to exhibits an a posteriori error control. These properties are stated in the following, we will also analyze their validity in the case of anisotropic meshes.

Theorem 1. The operator Π c satisfies the following properties:

1. ∃c 1 > 0 ∀v ∈ V 0 h , ∥∇Π c (v)∥ L 2 ≤ c 1 ∥∇v∥ L 2 (Ω) (stability) 2. ∃c 2 > 0 ∀v ∈ V 0 h , ∥v -Π c (v)∥ L 2 (Ω) ≤ c 2 h∥∇v∥ L 2 (Ω) (Interpolation inequality)
where Ω is a bounded Lipschitz domain in R n and h is the mesh-size.

Remark 1. The stability property mentioned above is even better when the constant c 1 admits a sufficiently small majorant. Indeed, the smaller is the constant c 1 , the better the stability and continuity of the Clément interpolation operator. According the stability proof addressed in [30], we have c 1 := M 1/p 3 M 4 M 5 such that:

M 3 := card {V j ∈ V (T h ) | ∃K ∈ S i V ij ∈ K}. M 4 := max z∈V (T h ) (1 + e z
) where e z is the constant of the Poincaré-Wirtinger inequality, its explicit bounds depends on diam(S i ). The bound for ∥.∥ L 2 was first established in [31].

M 5 := sup x∈Ω z∈V (T h ) h(x) q |Dϕ z (x)| q 1/q
, (ϕ z : z ∈ V (T h )) denote the nodal basis functions of all nodes with respect to the first-order finite element.

Notice that M 3 and M 4 are robust for small aspect ratios and M 5 may be very large for small angles in the triangulation.

Moreover, in [30] Theorem 2.6. an interesting result of convergence has been demonstrated, it is briefly reminded bellow

Theorem 2. For f ∈ W 1,2 (Ω), we have | Ω f • (v -Π c (v))dx| ≤ O(∥h∥ 2 ∞ )∥∇v∥ (Ω)
This theorem illustrates that the error of the recovered solution is of order 2 and thus of higher order compared to the first-order finite element method whose convergence is of order 1.

1. One of limitations of the Clément operator is that it does not preserve homogeneous boundary conditions, i.e., if u vanishes on the border, this is not generally the case for the Clément operator ( Π c u).

2. The technique presented above, can be generalized for others finite elements as well as for areas with curved boundaries [32], [33].

Scott-Zhang operator

As mentioned previously, Clément operator does not preserve boundary conditions. In [2], Scott and Zhang have addressed this issue and have defined an alternative interpolation operator. For instance, boundary condition of the form:

{v ∈ V h | ∇ • v = 0}
are preserved by Scott-Zhang interpolation operator. The construction of this operator is carried out as follows:

For every node V i , if V i ∈ ∂Ω, we restrict the interpolation to external facets S i connected to V i , (cf Figure 3.3 left). Otherwise, if V i is internal, we then interpolate on the set of cells S i connected to V i , (c.f Figure 3.3 right.)

We denote by n i the number of nodes belonging to S i and {ϕ i,q } 1≤q≤n i the restriction to S i of local basis functions associated to S i vertices.

△ △

▲ ▲ ▲ ▲ ▲ ▲ One can simplify the notation defining ϕ i,1 = ϕ i ∀V i ∈ V (T h ). Now we build a family {γ i,q } 1≤q≤n i as: for an integer q, 1 ≤ q ≤ n i , we define {γ i,q } 1≤q≤n i ∈ Vect {ϕ i,1 , ..., ϕ i,n i } as the unique function such that:

S i γ i,q ϕ i,r = δ qr 1 ≤ q, r ≤ n i (3.7)
where δ qr is the Kronecker delta. To simplify, we denote

γ i = γ i,1 ∀V i ∈ V (T h ) (3.8)
We then have

S i γ i (x)ϕ i (x)dx = δ ij i, j = 1, 2, ..., N (3.9)
where ϕ i is a nodal basis of V h . Thus, The Scott-Zhang interpolation operator is defined as:

SZ h : V 0 h (Ω) -→ V 1 h (Ω) v -→ SZ h v(x) = N i=1 ϕ i (x) S i γ i (ξ)v(ξ)dξ By taking v = ∇u h and x = V i , we obtain SZ h ∇u h (V i ) = S j γ i (ξ)∇u h (ξ)dξ Equivalently SZ h ∇u h (V i ) = n j=1 K j ∩S j ̸ =∅ γ i (ξ)∇ϕ j (ξ)u h|K j ∩S j dξ
Thanks to the change of variable, we obtain:

SZ h ∇u h (V i ) = n j=1 Kj | Kj | |K j | γ i ( ξ)J -T K j ∇ φj ( ξ)u h|K j ∩S j |detJ T K j |d ξ ( * * )
Here, J T K j is the Jacobian of the affine mapping T K mapping K to K, its determinant is equal to d!|K|. Thus

SZ h ∇u u (V i ) = n j=1 ( * * ) Kj γ i ( ξ)J -T K j ∇ φj ( ξ)u h|K j ∩S j |detJ K j |d ξ d! n j=1 |K j | By using a quadrature rule, we obtain ( * * ) ≈ N q=1 | Kj |w q γ i (x q )J -T K j ∇ φj (x q )u h|K j ∩S j d!|K j | (3.10) By simplifying | Kj | with d! we obtain ( * * ) ≈ N q=1 w q γ i (x q )J -T K j ∇ φj (x q )u h|K j ∩S j |K j | (3.11)
By definition of Gauss weights on the reference cell, we have

N q=1 w q = | K|
However, in our implementation, the weights are normalized, i.e

N q=1 w q = 1
Therefore, the factor | K| was added in front of w q in Equality 3.10. Just like Clément, Scott-Zhang interpolation operator admits stability properties, we recall them bellow through the following theorem:

Theorem 3. The following properties holds for Scott-Zhang interpolation operator:

1. ∃c 1 > 0 ∀v ∈ V 0 h , ∥SZ h (v)∥ L 2 ≤ c 1 ∥v∥ L 2 (Ω) (stability) 2. ∃c 2 > 0 ∀v ∈ V 0 h , ∥v-SZ h (v)∥ L 2 (Ω) ≤ c 2 h∥∇v∥ L 2 (Ω) (Interpolation inequality)
This theorem states that the Scott-Zhang operator is stable and consistent. In other words, the operator is well-posed and allows to estimate the error by passing to the limit in the interpolation inequality.

Zienkiewicz-Zhu interpolation operator (Least Square minimisation)

The Zienkiewicz-Zhu (ZZ) recovery method is based on a Least-Square minimisation. Like the two methods presented previously, ZZ recovery uses"patchs" (see Figure 3.4) around each node to reconstruct the gradient. We consider u a function defined on Ω ⊂ R n , we denote by u h the u interpolation of type P 1 on T h . We would like to find σ 1 an approximation of ∇u which is also of type P 1 . Clearly we can not use σ := ∇u h because this approximation is piecewise constant P 0 .

In order to reconstruct a P 1 gradient in a given node V i , we use x i , i ∈ 1, n where n is the number of sampling points, the set of elements K midpoints (those elements set up the patch S i associated to the node V i ) as sampling points.

With this new recovery procedure, we use a linear polynomial u * P = Pa in the node V i in a Least-Square sens as follows

F (a) = n i=1 (u h (x i ) -u * p (x i )) 2 = n i=1 (u h (x i ) -P(x i )a) 2 (3.12)
Typically, in the case of P 1 -finite elements and in 2D, the polynomial basis used reads as:

P = [1, x, y] and a = [a 1 , a 2 , a 3 ]
Respectively, in 3D, the polynomial basis is given by:

P = [1, x, y, z] and a = [a 1 , a 2 , a 3 , a 4 ]
The minimization condition of F (a) implies that the unknown a verifies:

n i=1 P T (x i )P(x i )a = n i=1 P T (x i )u h (x i )
This can be solved as a linear system as follows

P T (x i )u h (x i )
Thus, the recovered gradient ∇u h in the node V i is defined as follows in 2D

∇u h = ∇P V i = a 2 a 3
and in 3D

∇u h = ∇P V i =   a 2 a 3 a 4  
According to [36] the recovered solution u * converges faster toward u than u h . We can thus write:

∥u -u * ∥= O(h 1+η ) and ∥u -u h ∥= O(h) with η ∈ R > 0 (3.13)
∥.∥ is a Sobolev norm. Similarly, according to [36] and applying the triangular inequality to the error estimator ∥ē∥= ∥u

* -u h ∥≡ ∥(u -u h ) -(u -u * )∥, we obtain ∥u -u h ∥ + ∥u -u * ∥ ≤ ∥u * -u h ∥ ≤ ∥u -u h ∥ + ∥u -u * ∥ (3.14)
Using the notations above

∥e∥ + ∥e * ∥ ≤ ∥ē∥ ≤ ∥e∥ + ∥e * ∥ (3.15)
which is equivalent to

1 - ∥e * ∥ ∥e∥ ≤ θ ≤ 1 + ∥e * ∥ ∥e∥ (3.16)
θ is the effectivity index allowing to measure the quality of the error estimator ∥ē∥.

In fact, the error estimator is said to be asymptotically exact if θ approaches the unity as ∥e∥ tends to 0. In other words, the error estimator will convergences to the exact error ∥e∥ while this decreases. Obviously, the reliability of the error estimator ∥ē∥ depend on the quality of the recovered solution u * . Further details about error estimator are provided in Section 3.5.

Numerical results: Comparison of the implemented interpolation operators

In this section, we will present a series of numerical examples to assess the validity, accuracy and the implementation efficiency of the previously presented interpolation operators. The aim is to illustrate the convergence of the proposed gradient recovery techniques. To this end, we consider, on the one hand, a unit square domain Ω = [0, 1] 2 with f (x, y) = x + y, and on the other hand, the same domain with

f (x, y) = x 2 + 2y 2 .
As the exact gradient of functions above is known, we can determine how accurate are the presented gradient recovery techniques. To evaluate the accuracy of the proposed reconstructions, we compute the L 2 norm of the error between the exact and the recovered gradient on increasingly fine meshes. Figures 3.5 and 3.6 present the variation of the interpolation errors' L 2 norm induced by the reconstruction techniques on isotropic meshes with increasing number of elements. In Figure 3.5, as f is linear and we represent it using piecewise linear functions we obtain an interpolation error close to machine precision: around 10 -28 for Clément and 10 -16 for Zhu-Zienkiewicz. We also notice a slight increase of the error proportionally to the increase of the number of elements used for the computation. This is due to an accumulation of the truncation error. In Figure 3.6, we can clearly see how the interpolation error decreases as the number of cells is doubled.

Truncation errors linked to system inversion for Zhu-Zienkiewicz interpolate causes a loss of accuracy of interpolation errors compared to those obtained with Clément. N umber of elements N umber of elements

||(∇u -∇ R u h )||
∥(∇u -∇ R u h )∥

Error indicators for anisotropic mesh refinement

In this section, we are interested in building an error estimator used to build a metric tensor that modifies the mesh sizes in all directions in a way to provide a good capture of considered physical entities. The aim is to ensure that the discrepancy between the exact and the approximated solution is minimal by increasing locally the mesh density at the regions of high approximation error. In the following, we present details of the derivation of the a posteriori error estimator allowing to predict the regions of high solution variations.

Overview on error estimators techniques

When we solve a partial differential equation, the aim is to obtain a numerical solution within a prescribed tolerance using minimal amount of work. The main tool are error estimators in an a priori and a posteriori sense. The idea was first proposed in [26]. Recent developments based on anisotropic interpolation estimates have proved to be successful for generating meshes with high aspect ration [3,[7][8][9][10].

Traditionally, the quality of approximated numerical solutions is expressed through a priori error estimators. These estimators can be evaluated before the computation, they give upper-bounds on the difference between exact and approximated solution. This bound typically depends on the mesh size (which tends toward 0 in the mesh refinement) and an unknown constant that depends on the exact solution.

Unfortunately, in practice the upper bound can almost never be evaluated because it depends on the exact solution that is unknown. Consequently, they do not allow to quantify the finite element solution. However, they proved to be essential for convergence studies of numerical methods, see [19] and [12] for further details.

On the other hand, we havea posteriori error estimators allowing to give global and local information on the error of the numerical solution using only the computed numerical solution and known data of the problem. In fact, the aim of a posteriori error estimators is to provide upper bounds on the error between the numerical approximation and the unknown exact solution that can be computed in practice, once the approximation solution is known, cf, Verfürth [24], Ainsworth and Oden [25], Babuška and Strouboulis [26], Neittaanmäki and Repin [27], Han [28], or Repin [29]. Consequently, a posteriori error estimations can be used to : 1. estimate the error magnitude between the exact and the approximated solution. 2. know where the error is located. Thus, we can also hope to use them to design algorithms that meet the following requirements: 3. The accuracy given before the start of the simulation is obtained at the end of the simulation (error control ). 4. The amount of work required is reduced as possible (efficient computation).

A posteriori error estimators are characterized by the following properties:

Provide a fully computable error's upper-bound (without any unknown constant) using the approximated solution (guaranteed upper bound ).

Give the expression of this estimation locally, for example in each mesh element, and ensure that this estimation is locally an error's lower-bound up to a multiplicative constant(local efficiency).

Ensure that the effectivity index, given as the ratio between the estimated and the current error tends toward 1.

Ensure the three previous properties independently of the problem parameters and their variation.

Give estimations that can be evaluated locally (negligible cost).

We can distinguish four main frequently used a posteriori error estimators: the residual based, the hierarchical, the recovery based and the goal oriented methods. The residual based error estimators [13,14], as suggested by their name, are derived from the residual of the discretized problem, thus they are problem dependent and are easy to implement. However, they lack of precision due to successive approximation in their derivation. The hierarchical error estimators [15], [16], based on considering high order approximation of the exact solution as a good representation of the exact solution and evaluating the error (of the actual computed solution) with respect to it. Since these techniques require two solutions of different order they imply high computational costs. Moreover, the Goal-oriented error estimators [17][18][19][20] consist in exploiting the dual problem whose source term is the quantity of interest. In fact, the solution of the dual problem is used to optimize a given functional of interest. Consequently, the obtained mesh would be adapted to capture the solution features that are relevant to the considered functional, thus different choices of functional generate different adapted meshes. The last well known method of error estimations is based on gradient recovery techniques [1][2][3][4]21].

A posteriori error estimations

Given a Sobolev norm ∥.∥, the aim of a posteriori error estimators is to approximate the (unknown) error ∥u -u h ∥ by a computable quantity η caller error indicator depending on data such as the discrete solution u h , Ω and ∂Ω. To be reliable and and efficient, the a posteriori error estimator must provide lower and upper bounds for the approximation error [26] in the following sense:

C 1 η ≤ ∥u -u h ∥≤ ∥C 2 η (3.17)
where multiplicative constants C 1 and C 2 are strictly positive, they may depend on the mesh T h and are supposed to be bounded. The two bounds above, allow to proof the equivalence between the estimator η and the approximation error ∥.∥.

We will try to estimate C 1 and C 2 in the case of Lagrange P 1 .

1. Estimation of C 1 .

According to Céa's Lemma [44,45] we have

∥u -u h ∥ a ≤ M α ∥u -Π h u∥ a (3.18)
where M > 0 is the continuity constant and α > 0 is the coercivity constant associated to the bilinear form a(•, •), ∥.∥ a = a(•, •) 1/2 is the energy norm. On the other hand, according to Interpolation Inequality for Lagrange P 1 , we have

∥u -Π h u∥ L 2 (Ω) ≤ h 2 Ω ∥∆u∥ L 2 (Ω) (3.19)
where

h Ω = max K∈T h (h K )
Consequently, C 1 depends on h Ω , the continuity constant M and the coercivity constant α.

Estimation of C 2 .

Thanks to inverse inequality for piecewise-polynomials, there exists a positive constant C such that:

∥∇(u -u h )∥ L 2 (K) ≤ C∥β -1 K ∥ 2 ∥u -u h ∥ L 2 (K) (3.20)
where ∥β K ∥≤ h K ρ K with h K = diam(K) and ρ K := sup diam(S)|S is the inscribed circle contain When we deal with a quasi-uniform family of triangulation,

∃ 0 < C < 1 | h K ρ K ≤ C.
In this work, special emphasis in placed on error estimates prescribed to mesh edges also called edge-based error estimates. They comes from an a posteriori error analysis of the discrete solution. The main advantage of edge-based error estimates over cell-based error estimates is that they provide local directional information about the error. This approach will be discussed further in the section bellow.

Edge-based error estimation

We consider V h a P 1 finite element approximation space defined as follows:

V h = v h ∈ C 0 (Ω) | v h|K ∈ P 1 , K ∈ T h (3.21)
where the triangulation T h of Ω is a family of disjoint elements K (segment, triangle, tetrahedron, ...) such that Ω = K∈T h

K.

Let u h ∈ V h the P 1 the interpolated solution obtained by applying the Lagrange interpolation operator to a regular function u ∈ C 2 (Ω). Let π h : C 2 (Ω) → V h a Lagrange interpolation operator. At each vertex of the mesh, we denote by u i the nodal value of π h at V i , for all i ∈ 1, ..., N where N is the number of nodes, we have:

u i := π h (V i ) = u h (V i )
for the sake of simplicity, we introduce the following notations for a generic node V i :

V ij = V j -V i and u ij = u j -u i
In the case of P 1 finite element, the gradient of u h is a piecewise constant vector field discontinuous from element to element. However, its projection onto nodes is continuous and depends only on the nodal values of u, therefore , by considering a Taylor expansion of the variable on interest on the node V i , one can write

u j = u i + ∇ ij u h • V ij (3.22)
Using notations introduced previously and rearranging the terms, we obtain

∇ ij u • V ij = u ij (3.23)
Following the works of [4], an a posteriori error estimate based on the length disctribution tensor approach and the associated edge-based error analysis is defined:

∥∇u h • V ij -∇u(V i ) • V ij ∥ L 2 ≤ max y∈[V i ,V j ] |H u (y)V ij • V ij | (3.24)
where H u (y) = ∂u ∂x i ∂x j ij is the associated Hessian of u evaluated at a point y of the edge V i , V j . We denote by G i = ∇ R u(V i ) and G j = ∇ R u(V j ) the recovered gradient of u at the node V i and V j respectively, thus we have

∇u • V ij = G j -G i = G ij (3.25)
Consequently, the projection of the Hessian based on the recovered gradient at the extremities of the edge is obtained as follows:

G ij • V ij = H u (V i )V ij • V ij (3.26)
We denote the error along the edges using the following expression:

e ij = |G ij • V ij |
For P 1 -finite elements, we use this projection as an evaluation of the L 2 norm of the interpolation error along the edge. However, this error cannot be evaluated exactly as it requires that the gradient of u to be known and continuous at the nodes of the mesh. For that reason, we resort to a gradient recovery procedure such as those presented in Section 3.2. We will give here bellow, another gradient recovery technique based on a least square minimization, which it is developed in [35].

Edge-based gradient recovery

The gradient recovery technique that we will present here comes from the work in [4], it is based on a least squares approximation of ∇u h along the edges of the mesh. Using the length distribution tensor, a continuous gradient can be defined directly at the nodes of the mesh and depending only on the values of the interpolation solution.

For each node V i , we seek G i ∈ R d , the proposed gradient reconstruction satisfying:

G i = arg min G∈R d j∈Γ(i) |(G -∇u h ) • V ij c.f 3.25 | 2 = arg min G∈R d j∈Γ(i) |(G • V ij -U ij )| 2 (3.27)
The minimum can be obtained by setting the derivative of the equality above to zero i.e.

j∈Γ(i)

(G i V ij -U ij )V ij = 0 ⇐⇒ G i j∈Γ(i) V ij ⊗ V ij = j∈Γ(i) U ij V ij (3.28)
where ⊗ identifies the tensor product, and

X i = j∈Γ(i) V ij ⊗ V ij is the length distri-
bution tensor. To simplify the notation, we define

U i = j∈Γ(i) U ij V ij (3.29)
We get:

G i = (X i ) -1 U i (3.30)
where X i is the length distribution tensor. The reconstructed gradient GG i is defined under the assumption that X i is a positive definite tensor. Consequently, the estimated error is thus written as:

e ij = |G ij • V ij | (3.31)
This quantity can be used to account for several sources of error by applying error estimator, cf Equation 3.31, to each component separately. This error estimator will be related to a metric field in order to drive the mesh adaptation procedure.

Estimators quality

Once an estimator is defined, we must derive criteria and means allowing to judge the performances of this estimator. In general, we must seek to build an estimator whose asymptotic behavior, this means when the mesh size tends toward 0, follows that of the error. In other words, this behaviour is translated into the existence of two constants C 1 and C 2 depending on the problem's data and the discretization and not the elements' size verifying the Inequality 3.17 as h approaches 0. Thus, to measure estimators' performance, there are different intrinsic criteria.

Efficiency index

The efficiency index is defined as the ratio between the computed error via an estimator e estimated and the true error e true :

θ = e est e true (3.32)
Unless we have an analytical solution, the true error is computed as the difference between a solution obtained using a very fine mesh and that computed on a given mesh. An efficiency index close to unity characterize a good estimator. If this property is satisfied when the mesh size h approaches 0, the estimator is said to be asymptotically exact. The estimation must be sufficiently accurate (θ close to 1) for coarse meshes in order to be used for meshes' refinement. However, it is advisable to overestimate the error (θ > 1 ) in order to be usable as a stopping criterion in an iterative process.

Studies have been conducted by Strouboulis and Haque [37] and by Babuška [38] based on different test-cases with an analytical solution or for which a reference solution can be determined on a very fine mesh. Generally, it appears that the quality of estimators depends on the mesh's topology, solution's regularity and elements' regularity (anisotropy level, flattening). The global index θ only does not fully allow to determine the way the local error behaves: it may be possible, in some cases, that the estimator is satisfactory according to this criteria without managing to localize the regions where the error in the energy norm is high [39].

Robustness index

To free one self from the local behaviour of the error estimator, another criteria has been defined by Babuška and Rodriguez [46]. For a given triangulation T h , the local efficiency index θ k related to the cell k is defined by:

θ k = e k est e k true (3.33)
Babuška and Strouboulis have shown that we can have access to the asymptotic variation's range of θ k [38]:

0 < C k inf ≤ θ ω ≤ C k sup < +∞ (3.34)
The deviation to 1 of θ k is measured by the index R k defined by:

R k = max |1 -C k sup | + |1 -C k inf |, |1 - 1 C k sup | + |1 - 1 C k inf | (3.35)
The robustness index R is the greatest R k obtained by browsing through all the mesh's cells belonging to a class of triangulation. This approach allows to highlight the discretization error on a subdomain but also the solution's regularity and the distortion of mesh's elements. However, this technique seems complicated to be implemented for complex meshes used to deal with engineering problems.

Numerical validation: error estimator and Effectivity index

In the following, we will compare the performance of the error estimator, defined as the difference between the recovered gradient and that given directly by the finite element resolution, using the different recovery techniques presented previously. We consider a two dimensional problem as in [36] -∆u = f in Ω (3.36)

with boundary conditions u = ū on ∂Ω (3.37)

where Ω = [0, 1] × [0, 1] and f is computed so that the exact solution is given by:

u(x, y) = (x 3 + y 3 )sin(xy) (3.38)
The error in the energy norm is written as

∥e∥= Ω (∇u -∇u h ) T (∇u -∇u h )dx 1/2 (3.39)
The error estimator is defined by

∥ē∥= K∈T h (∥ē∥ K ) 2 1/2 (3.40)
where the element error estimator is given by

∥ē∥ K = K (∇ R u -∇u h ) T (∇ R u -∇u h )dx 1/2 (3.41)
where the subscript R represents the recovery procedure and will be replaced by CL, AM, ZZ and ER for respectively Clément, Arithmetic-mean, Zhu-Zienkiewicz and Edge-Recovery. The numerical results of the performance of the error estimators for this example are plotted in Figures 3.7 and 3.8. We can observe that θ CL and θ AM are close to one for all meshes used even for very coarse ones. In fact, θ CL and θ AM has nearly constant values close to 1 for all mesh requirements. On the other hand, θ ZZ converges to one rapidly for all the meshes tested while the mesh is refined, the farthest value from 0 is 1.1 obtained with a very coarse mesh containing 2988 elements. consequently, recovery techniques used will overall provide reliable a posteriori error estimators.

In order to study the mesh convergence rates of the previously presented recovery techniques, we have conducted convergence studies of the L 2 norm of interpolation errors obtained respectively Lagrange, Clément (CL), Arithmetic-mean (AM), Zhu-Zienkiewicz(ZZ) and Edge-Recovery interpolation techniques. Figure 3.7 summarizes the convergence order estimates. We can observe that, as expected from theoretical analysis (see Section 3.2.3), the convergence rate of ∥∇u -∇u h ∥ L 2 is of first order whereas that of ∥∇u -∇ R u h ∥ L 2 is slightly better than the first order.

The performance of gradient recovery techniques for an arbitrarily generated anisotropic mesh shown in Figure 3.9 are now considered. The numerical results of θ CL , θ AM and θ ZZ are reported in Table 3.1. It is observed that effectivity indices are again very close to unity. This shows that the error estimator ∥ē∥ will provide reliable error estimation not only for uniform meshes but also for arbitrary irregular meshes, an important feature in adaptive analysis procedure.

MeshElements

θ CL θ AM θ ZZ ∥ē∥ CL ∥ē∥ AM ∥ē∥ ZZ 1336 elements 9.66e -1 9.68e -1 9.65e -1 1.04e -1 1.03e -1 4.22e -1
Table 3.1: The effectivity index on an irregular mesh for various gradient recovery techniques.

10 -2.5 10 -2 10 -1.5

-5 The performance of gradient recovery techniques for arbitrarily generated irregular meshes shown in Figure 3.10 are now considered. As for regular meshes, Figure 3.11 compares the evolution of ∥∇u -∇ R u h ∥ L 2 with that of ∥∇u -∇u h ∥ L 2 using the different gradient recovery techniques. Indeed, as for regular meshes, we can observe that the convergence rate of ∥∇u -∇u h ∥ L 2 is faster than that of ∥∇u -∇ R u h ∥ L 2 no matter the non-smooth interpolation technique used. The numerical results of θ CL , θ AM , θ ZZ and θ ER are reported in Figure 3.12. It is observed that the effectivity indices are again very close to unity.

-4 -3 -2 h ∥∇u -∇u h ∥ L 2 ∥∇u -∇ R u h ∥ L 2 (ZZ) ∥∇u -∇ R u h ∥ L 2 (CL) ∥∇u -∇ R u h ∥ L 2 (AM ) ∥∇u -∇ R u h ∥ L 2 (ER) h h 2 Figure 3.7: ∥∇u-∇ R u h ∥ L 2 vs ∥∇u-∇u h ∥ L 2 using Zhu-Zienkiewicz,
The performances of the implemented recovery techniques has also been studied in parallel. Indeed, as one can see through Figures 3.13 and 3.14, where we compute the interpolation error and the efficiency index obtained using different numbers of processes, values are the same regardless the number or processors used. The edge error estimator will be followed by the derivation of an optimal metric construction used to drive anisotropic mesh adaptation. Further details are provide in the upcoming section. 3 Error indicators and anisotropic mesh refinement

10 -2.4 10 -2.2 10 -2 10 -1.8 10 -1.6 10 -1.4 10 -1.2 -5 -4 -3 -2 -1 h ∥∇u -∇u h ∥ L 2 ∥∇u -∇ R u h ∥ L 2 (ZZ) ∥∇u -∇ R u h ∥ L 2 (CL) ∥∇u -∇ R u h ∥ L 2 (AM ) ∥∇u -∇ R u h ∥ L 2 (ER) h h 2 Figure 3.11: ∥∇u -∇ R u h ∥ L 2 vs ∥∇u -∇u h ∥ L 2 using
∥(∇ R u h -∇u h )∥ L 2 ∥(∇ R u h -∇u h )∥ L 2 [Cl] ∥(∇ R u h -∇u h )∥ L 2 [ZZ] ∥(∇ R u h -∇u h )∥ L 2 [AM ] ∥(∇ R u h -∇u h )∥ L 2 [ER]
Figure 3.13: The evolution of error estimator ∥∇ R u h -∇u h ∥ L 2 obtained using different numbers of procs. Within the context of finite elements method numerical simulations, non structured mesh adaptation has amply proven its usefulness. Indeed, this technique allows in theory, (i) significantly reduce the computational cost (by reducing the number of degrees of freedom) and to achieve the numerical solution with the desired precision. This has the effect of allowing the simulation of complex phenomena in 3D without resorting to heavy computing means (architectures and parallel algorithms, for example). In principle, a mesh adaptation method is based on the estimation of the error made and its translation into a numerical information named Error estimator used as an indicator for the adaptation procedure to determine the mesh sizes accordingly in the different regions [46,47,49,50]. Thus, we converge gradually toward a "limit" mesh corresponding to a numerical solution with the desired precision. Anisotropic mesh adaptation allows to control the size, shape and the orientation of elements to better represent the numerical solution with a reduced computational cost. The first studies on anisotropic mesh adaptation were proposed in the late 1980s [42], [43]. A new method based on the computation of the edge lengths in a discrete metric defined on the mesh nodes was adopted later. This anisotropic metric is build using an a posteriori error estimator based on a discrete approximation of the Hessian of the approximated solution. The aim is to equidistribute the interpolation error on the mesh edges. In the following, we will detail the construction of the metric following the studies developed by the CFL research group at CEMEF.

Theoretical framework

Let us provide some theoretical framework and basic definitions that are useful for the development of anisotropic mesh adaptation procedure. In the following, the physical domain Ω is an open subset of R d , its discretized using a finite element mesh Ω h = K∈T h K, with K being a d-simplex (segment, triangle, tetrahedron, ...)

and T h an admissible mesh.

For any v ∈ V h , we call a metric M i ⊂ R d×d , a d × d symmetric positive definite matrix such that m i is the M i -conjugate scalar product.

∀x ∈ R d , < x, y > m i =< M i x, y > Or, in algebraic form, ∀x ∈ R d , < x, y > m i = y T M i x
A metric M i can be diagonalized into

M i = R    λ 11 • • • 0 . . . . . . . . . 0 • • • λ n    R T
where R is a rotation matrix whose columns are the eigenvectors of M v and (λ i ) i=1,...,d ≥ 0 are the corresponding eigenvalues.

The main feature feature of anisotropic mesh adaptation techniques is the control of elements orientation and their sizes. This is exactly the information that is brought by metric tensors. Indeed, the metric eigenvalues λ i provide the mesh sizes while eigenvectors give directions.

Metric construction at the mesh vertices

Let us consider the star-topology on the mesh T h , that is defined as the vertexvertex incidences. We note by Γ(i) (see Figure 3.17) the set of nodes V j sharing an edge with V i :

Γ(i) = {V j ∈ V (T h ), ∃K ∈ T h , V ij ∈ K}
The aim is to define a discrete metric in V i such that the length of each edge connected to V i is equal to 1 with respect to the metric field M i . By definition of the M i -conjugate scalar product, we have:

(M i V ij , V ij ) = 1, ∀j ∈ Γ(i) (3.42)
As stated in [4], the problem of finding a unitary metric is equivalent to the following least squares problem:

M i = arg min M ∈R d×d sym j∈Γ i M V ij • V ij -card(Γ(i)) 2 (3.43)
Then summing up the Equality 3.42 over j, we obtain

j∈Γ i (M i V ij , V ij ) = j∈Γ i 1 = card(Γ(i)) (3.44)
In addition, by using the Einstein notation for the tensor scalar product, the previous equation can be written equivalently as follows: M i :

j∈Γ i V ij ⊗ V ij = card(Γ(i)) (3.45)
Thus

M i : 1 card(Γ(i)) j∈Γ i V ij ⊗ V ij X i = 1 (3.46)
For a finite element admissible mesh, Γ(i) has to form at least d non co-linear edges with vertex V i (see Figure 3.17). X i is the length distribution tensor, it gives an Reynolds number). The obtained mesh will be used all over the simulation. In this case the mesh adaptation criteria are geometric and do not evolve with the solution profile. The implicit strategy, meanwhile, consists in dynamically adapting the mesh and minimizing as much as possible the global equidistributed error. In this case, the anisotropic adaptation is driven through the behaviour of an a posteriori error estimator whose theory will be developed in the next section, taking into account the geometry as well as the evolving solutions on the desired indicator fields. Therefore, the following vector depending on different sources of error of CFD problems can be used as a criterion for mesh adaptation:

e(x i ) =    U i ∥U ∥ L 2 , U i max j |U j | , T i T max , ϕ i ϕ max   
Because all adaptation fields are normalized (the velocity components by the velocity norm, the temperature and the level set function by their respective maximum), a field that is much larger in absolute value does not dominate the error estimator, consequently, the variation of all fields are fairly taken into account.

As an illustration, Figure 3.19 presents the zero isovalues of the level set function for two immersed solid bodies. From an isotropic mesh and after several iterations, one can clearly see through Figure 3.19 that the interface between these objects and the surroundings is well adapted. The node distribution is highly dense in the vicinity of the interface allowing its accurate capture. This validates how the developed algorithm optimizes the distribution of the nodes to produce a sharp anisotropic mesh that is well adapted based on a given variable (here, the level set of the immersed body). Figure 3.20 shows that after few iterations, a well respected shape in terms of curvatures, sharp angles and thin layers is obtained for both ingots.

This step in commonly used to initialize an industrial setup and was extended to deal with boundary layers capturing in complex flows. We will propose another numerical example to illustrate the anisotropic mesh adaptation based on the normalized velocity field U i max j |U j | (see Figure 3.22). The geometry of the test case if given through Figure 3.21. The air is pumped from the left inlet at 1m/s and vented out the enclosure through the two outlets positioned at the vertical wall. The shape of the fluid flow is well captured by applying the anisotropic mesh adaptation. The fluid flow is described thanks to the discretization in space of incompressible Navier-Stokes equations and the turbulence equations by an unstructured grid stabilized finite element. Consequently, the numerical oscillations are well captured and smooth solutions are obtained. As the adaptation is performed on the velocity field, one can detect steep anisotropic characteristics as a response to the sharp gradients of the velocity field, see Figure 3.22 right. We can also observe the high density of the nodes in the wake and its gradual decrease as we move away. More test problems on this subject are treated in [34].

An important point to mention here is that this mesh adaptation technique works under the constraint of a fixed number of nodes. With such an advantage, we avoid a drastic increase of mesh complexity and hence computational cost when dealing with complex CFD problems.

However, one serious drawback of many boundary values CFD problems is that they may contain singularities and/or have layers where the solution changes rapidly. For such non-smooth functions, the application of interpolation operators for nonsmooth functions over polytopal meshes will be introduced and analysed later. We depict in Figure 3.24 the general diagram summarizing our work. We mention that a posteriori error estimators also allow to refine the mesh where the error is high even if it is not the aim of this thesis.

Conclusion

In this chapter, we have devised a posteriori error estimates between the exact and the numerical solution. These estimators are based on gradients recovery which makes them convenient because they does not request information on the nature of the starting equation. They are based on comparing the numerical solution's gradient and the recovered gradient, more regular then the first one. This technique is known as "Superconvergent Patch Recovery", it has been initially developed by Zhu-Zienkiewicz [40] for a linear elasticity problem. It has been studied in electromagnetism recently by Nicaise [41]. Several gradient recovery techniques has been introduced and analysed: those based on L 2 projection techniques and those based on Least-squares minimization. Their efficiency and robustness in a posteriori estimates has been compared through different numerical tests. These indicators allow to devise adaptive stopping criteria for linear solvers, they guarantee the global error control. Consequently, they are used to reduces the number of iterations and the computational time. We will analyse this point in more detail in the next chapter. We will show in the following chapter how an optimal control of iterative solvers through a stopping criterion depending on the a posteriori error estimator defined above allows to make the resolution procedure faster and more efficient.

Résumé du chapitre en français

Malgré le développement régulier des puissances de calcul, l'adaptation de maillage anisotrope semble être un moyen incontournable pour réduire la complexité des simulations numériques, en particulier lorsqu'il s'agit d'applications industrielles 3D complexes. En effet, il est hautement souhaitable d'obtenir la solution la plus précise en réduisant le coût de calcul au maximum. L'objectif de l'adaptation de maillage anisotrope est de s'assurer que l'écart estimé entre la solution exacte et la solution approchée est minime. A cet effet, la densité du maillage doit être augmentée localement dans les régions où l'erreur d'approximation est élevée de manière à prévoir avec précision la solution. Dans ce chapitre, l'estimateur d'erreur a posteriori calculé le long des arêtes du maillage est utilisé pour mesurer localement l'écart entre la solution exacte et la solution approchée. De plus, l'estimation de l'erreur nécessite que le gradient soit connu et continu aux noeuds du maillages. En revanche, dans le cadre des éléments finis P 1 le gradient est constant par maille. Pour cette raison, nous nous sommes appuyer sur des méthodes de reconstruction du gradient au niveau des noeuds du maillage. Les techniques d'interpolation étudiées et implémentées sont basées sur deux stratégies principales: projection L 2 pour Clément [1] et Scott-Zhang [2] et minimisation au sens des moindres carrées pour Zienkiewicz-Zhu [3] et Edge-based [4].

Nous avons ensuite présenté divers exemples numériques et benchmarks tirés de littérature afin d'évaluer la qualité de l'estimateur d'erreur calculé avec les différentes techniques d'interpolation implémentées. Nous avons également détaillé l'implémentation parallèle des différents opérateurs d'interpolation étudiés. En utilisant le principe d'équidistribution de l'erreur et sous la contrainte d'un nombre fixe de noeuds dans le maillage, un ensemble de facteurs d'étirement associés aux arêtes du maillage sont déterminées. Une simple modification du tenseur de distribution de longueur pour tenir compte de ces facteurs d'étirement aboutira à une nouvelle métrique conduisant à un maillage anisotrope et bien adapté.

Introduction

Approximated solutions obtained from the finite element discretization of the Convectiondiffusion-reaction and Navier-Stokes equations are obtained thanks to an iterative procedure. Their accuracy is controlled by a stopping criteria used to drive the convergence of the iterative algorithm. In this Chapter we will give a brief introduction of iterative solvers. For more details, we refer the reader to Appendix A. Furthermore, we will see how a posteriori error estimates developed in the previous Chapter 3 allow to formulate the so-called adaptive stopping criterion for our solvers. This leads to an important time saving while ensuring accurate results. We will show through different numerical examples that we can reduce significantly the number of iterations and the computational time required.

Iterative solvers for sparse linear systems

The resolution of Navier-Stokes or Convection-diffusion-reaction equations implies the resolution of a linear system. This system presents several particularities. First of all, the linear system is of high dimension. A standard 3D Navier-Stokes mesh consists of tens of millions of nodes. As we have five unknowns per node (the resolution of CDR and Navier-Stokes coupled), the linear system contains more than thirty million unknowns. Furthermore, this system is obtain through a discretization of equations using the finite element method. Considering P 1 Lagrange finite elements, the fill graph of the matrix [1] is identical to the mesh graph. The number of neighbour vertices of another defined the number of non-zero blocks of every line and column of the matrix, its of the order of several tens. This is to be compared to the tens of millions of the mesh nodes. The matrix is then extremely sparse.

There are two main categories of methods to solve large size linear systems: iterative methods or direct ones. Iterative methods, as their name suggests, create a sequence of approximation of the solution. The user define a resolution tolerance, defined as the maximal norm of the residual r = b -Ax from which the user is satisfied with the solution. The direct method compute the exact solution (except for numerical errors) x = A -1 b. Coutinho and Jouglard [START_REF] Jouglard | A Comparison of Iterative Multi-level Finite Element Solvers[END_REF], proposed a revisited termination condition of the Conjugate Gradient algorithm based on the energy norm of the increment instead of the residual norm, which appears to be more convenient. Indeed, it has been deduced through different numerical examples that acceptable solutions can be obtained with less iterations than those needed to reach the numerical convergence using the residual norm. The goal would be to extend this strategy to the General Minimal residual Method (GMRES).

Advantages and disadvantages of these two classes of methods are easy to understand. Iterative methods are based on the creation of approximation sequence either by iterations of the matrix ( as for Krylov spaces based methods), or by multiplication by parts of the matrix (relaxation methods as Jacobi, Gauss-Seidel and SOR). These methods are based on matrix-vector products, for which storage formats can be optimized [1], Chapter 3. The memory cost of these algorithms is only that of the sparse matrix and of a number of vectors. For problems addressed in this thesis, the number of vectors composing the Krylov space is around a hundred.

Direct solvers compute the application of A -1 to a vector. In general, this requires a decomposition of the matrix A. For example, the LU decomposition [2,3]. One this decomposition is computed, we can solve almost instantly several different source terms for the same matrix. The resolution time does not depend on the matrix conditioning κ(A) = ||A||||A -1 ||. Conversely, some of iterative methods, like the conjugate gradient for symmetric positive definite matrices, have a convergent rate directly dependent on κ(A). Direct methods have an important disadvantage of the required memory. For big systems, the LU decomposition costs as much in memory as the matrix. A direct method is then very suitable. However, the LU decomposition of a sparse matrix gives factors L and U much fuller then the initial matrix A. There are direct solvers as MUMPS [6,7] adapted to sparse matrices. The storage cost of full factors is then identical of that of the matrix. Unlike iterative methods, the parallelization of direct methods is difficult, as shown by the parallelization of a ILU preconditioning method [1].

Preconditioning

The preconditioning consists on transforming the linear system to an equivalent form that is numerically easier to solver with the aim of accelerating the linear resolution. It consists on transforming the matrix A of the linear system Ax = b to bring it as closer as possible from the identity. To do so, we need to look through a matrix M , whose inverse is simple to calculate and approach A -1 . Thus, AM -1 ≈ I and the linear system will be easier to solve. Left preconditioning consists on multiplying the left hand side of the system by M -1 . We then obtain

M -1 Ax = M -1 b (4.1)
The right preconditioning is simply deduced from Ax = AM -1 M x. A variable change allows to obtain

AM -1 y = b x = M -1 y (4.2)
As one can see through Appendix [? ] , the GMRES method is very sensitive to small eigenvalues of the matrix. A good preconditioning allows to to avoid having too small eigenvalues harmful to the convergence. A preconditioning matrix must be quick to assemble and the product of this matrix by a vector is inexpensive. In fact, in the GMRES method, the preconditioning matrix is applied after every matrix-vector product in the case of a left preconditioning ( or before in the case of a right preconditioning). The resolution of the GMRES method requires hundreds or even thousands of iterations. The application of preconditioning must then be fast.

Similarly, preconditioning creation time must be limited. The limiting case M -1 = A -1 leads to use a direct method to inverse the matrix A. The GMRES algorithm will then need one iteration to converge. The inverse of a sparse matrix is generally full. A good preconditioning will then require a lot of memory.

It is therefore necessary to find a good compromise between the preconditioning quality, and creation and application costs, both in memory that in time. This makes Saad say in [1] that "Finding a good preconditioner to solve a given sparse linear system is often viewed as a combination of art and science. Theoretical results are rare and some methods work surprisingly well, often despite explications" There are many methods to approach the inverse of a matrix. They are been addressed in [1,8]. We can mention the so-called SOR (Successive Over Relaxation) methods, ILU methods addressed in [4], are an adaptation of the LU factorization to sparse matrices. Finally, every linear system resolution method can be used to approach the effect of the inverse of a matrix on a vector. Consequently, the use of the algebraic Multi Grid methods (AMG for Algebraic Multi Grid) as preconditioner for others iterative solvers is very successful [5]. The main disadvantage is that an approximate resolution method does not always guarantee to obtain the same preconditioner. We then must use a solver that accept the variation of the preconditioner during iterations, such as GMRES [9]. Others preconditioning methods, not based on linear system resolution and factorization, exists. For example, the SPAI method (SParse Approximate Inverse) [10,11] is based on a least-square minimization of the difference between AT and the identity for every column of the matrix A, T is the preconditionner.

Finally, the preconditioning can be based on a simpler matrix than A. Finite volume methods of higher order lead to fuller matrices those with linear finite elements. The factorization process is then much more expensive. Using a matrix of lower order discretization allows to have a lighter and more efficient preconditioning [12,13].

Stopping criterion for iterative algorithms

Iterative algorithms provide a sequence to the solution u k k = 0, 1, 2, ... such that r k = (b -Au k ) → ∞ 0 as k → ∞. In practice, the number of iterations of an iterative algorithm is finite. Thus, we do not know how many mathematical operations are required to obtain a satisfactory solution. At some point, we merely declare that the solution is close enough. What is a rational criterion for deciding when to stop? It is necessary to find a compromise between quickly resolving linear systems and not perform many iterations. Naturally, the numerical resolution of a PDE generates errors which are typically of three natures: 1) error related to the numerical scheme used and to the mesh chosen, 2) error linked to the use of the non linear solver and the error linked to the linear algebraic solver. When these errors are very small, it means that the system has been solved almost exactly. We will try to estimate the three sources of error caused by the numerical simulation. This estimation of error components will allow in particular to reduce the number of iterations of the iterative algorithm without affecting the quality of the solution.

When we use an iterative algorithm, the idea is to find a solution u n h obtained at the n th step of the iterative method, that will be a good approximation of the exact solution u.

∥u -

u h ∥≈ ∥u -u n h ∥ (4.3)
where u h is the approximated solution obtained with a direct solver, ∥.∥ is a norm (or seminorm).

The a priori approach

In the case of an a priori error estimation the exact solution of a PDE and the approximated one u h are linked through the following inequality

∥u -u h ∥ L 2 ≤ Ch s ∥u∥ L 2 s ∈ N * (4.4)
where h is the mesh step and C is a constant depending on the exact solution u. This type of estimate provides important information about the convergence of the numerical method used. For further example of a priori error estimates, the interested readers can refer to Strang [16], Johnson [17], Ern et Guermond [18] or Maday, Bernardi, Rapetti [19]. In the context of conforming finite elements, a priori error estimates rely heavily of the famous Céa lemma [20] ∥u -

u h ∥ L 2 ≤ C∥u -Π h u∥ L 2 (4.5)
The main disadvantage of these a priori estimates is that the upper bound of the error is not computable in practice as it depends on the exact solution u, which is unknown.

The a posteriori approach

Since the 1970's, with the work of [21], inspired by Prager and Synge [22], estimates named a posteriori error estimates develop. By "a posteriori " we mean that the estimation takes place after the computation of the approximated solution u h . Among the long list of works on the subject, we can cite the works of Babuška and Rheinboldt [23], Ainsworth and Oden [24], Repin [32], Verfürth [28], and the contributions of Ern and Vohralík [27] for elliptic problems. For parabolic problems, we can cite the contributions of Verfürth [28], Nicaise and Nicaise and Soualem [29], Picasso and Prachittham [31] and Ern, Swears and Vohralík [30]. These estimates are defined through following inequality:

∥u -u h ∥ L 2 ≤ η(u h ) (4.6)
where η(u h ) is a computable quantity depending only on the approximated solution u h , ∥u-u h ∥ L 2 is the approximation error. According to Vohralik [START_REF] Vohralík | A posteriori error estimates, stopping criteria, and adaptivity for two-phase flows[END_REF], an a posteriori error estimator is optimal if it satisfies the following properties:

1. The upper bound of the error is guaranteed : The condition [4.6] must be verified and the estimator η(u h ) must be completely computable in function of u h .

2. Efficiency: For steady problems, for every element K of the mesh T h , the local estimator η k (u h ) must represents a lower bound of the error in the neighbourhood of K (represented by ϵ K )

η K (u h ) ≤ C ef f ∥u -u h ∥ ϵ K (4.7)
For unsteady problem, at every time step 1 ≤ n ≤ N t , we have

η n K (u h ) ≤ C ef f ∥u -u h ∥ ϵ K ×(t n-1 ,t n ) (4.8)
3. Asymptotic accuracy: the efficiency index I ef f is defined by

I ef f = η(u h ) ∥u -u h ∥
For steady problems (4.9)

or

I ef f = Nt n=1 K∈T h (η n K ) 2 1 2
∥u -u hτ ∥ For unsteady problems (4.10)

tends to 1 when we refine the mesh T h and the time step.

Robustness:

The constant C ef f must be independent of the problem parameters.

5. Low computational cost: a lower computational cost of local estimators η K .

Distinction of error components:

We must be able to estimate all error components of the numerical simulation.

In our study, the idea is to find a solution u n h obtained at the n th step of the iterative solver allowing to minimize the real error defined by ∥u -u n h ∥ L 2 . However, this error verify the following control:

∥u -u n h ∥ L 2 ≤ ∥u -u h ∥ L 2 +∥u h -u n h ∥ L 2 (4.11)
Thanks to anisotropic mesh adaptation, the approximation error is controlled as follows:

||u -u h || ≤ C K∈T h η 2 K 1/2 . (4.12)
Indeed, the anisotropic adaptive mesh is obtained by measuring the approximation error through the a posteriori error estimator and by refining the mesh where the error is the greatest. The a posteriori error estimator is locally defined as reads

η p K = K H(u h (x K ))(x -x K ) • (x -x K ) p dx . (4.13)
Consequently, the algebraic error must be smaller than a fraction of the approximation error to stop solving

||u h -u n h || L 2 ≤ C 2 (0.01) K∈T h η 2 K 1/2 . (4.14)
The choice of the scaling factor 1% is purely heuristic. In the literature, it has been suggested that the algebraic error can be approximated by the Euclidean norm of the residual (see [START_REF] Manzinali | Adaptive control for iterative solvers in an FE framework with mesh adaptation, for CFD simulations of industrial processes[END_REF] for more details) ∥r∥.

According to the following result, C 2 depends on

∥A -1 ∥ ∥e h ∥ L 2 (Ω) = ∥u h -u∥ L 2 (Ω) = ∥A -1 (Au h -b)∥≤ ∥A -1 ∥∥r n ∥
Using the definition of condition number κ(A)

κ(A) = ∥A∥∥A -1 ∥= |λ max | |λ min | = ρ(A)ρ(A -1 ) (4.15)
where ρ(A) = max {|λ i |, λ i eigenvalues ofA} is the spectrum of the operator A. On the other hand, we have ρ(A) ≤ ||A||.

Then ||A -1 || ≤ 1 |λ min | (4.16)
Using this result, the adaptive stopping criterion is defined as

∥r n ∥ ⩽ 1 |λ min | (0.01) K∈T h η 2 K 1/2 . (4.17)
This approach has been implemented in our library Cimlib-CFD. The Portable Extensible Toolkit for Scientific computation (PETSc) library has been used to obtain an estimation of λ min . PETSc, [START_REF] Balay | Efficient Management of Parallelism in Object-Oriented Numerical Software Libraries[END_REF] is an open-source suite od data structures and routines to solve scientific applications modelled by PDEs. The General Minimal RESidual (GMRES) method introduced by Saad and al. in [14] implemented in PETSc has been used to solve our non-symmetric problems. In the upcoming chapter, we will evaluate the efficiency of the developed numerical tools in resolving large scale and long time industrial simulations.

Numerical results

In this section, we will present several numerical examples to demonstrate the efficiency of the proposed adaptive stopping criterion in terms of computational time and error control. The aim is to validate the method on known problems by studying convergence, accuracy and efficiency of the considered test cases. We compare the new adaptive

∥r n ∥ ⩽ 1 |λ min | (0.01) K∈T h η 2 K 1/2 . (4.18)
with a fixed stopping criterion, where we impose a given precision, e.g, ϵ = 10 -7 ,

2D Taylor-Green problem

The first benchmark used for the validation and the convergence characterization of the code is the test case of Taylor-Green vortex in 2D proposed in [43] The non-linear system that stems from the finite element discretization of this (parabolic/hyperbolic) problem is solved using a GMRES algorithm with an Incomplete LU factorization (ILU) as a preconditioner.

The TG vortex was mainly used to test the sensitivity of the adaptive stopping criterion to the inherent instabilities present in vortices.

3D Taylor-Green problem

We perform the same test on a 3D Taylor-Green test case with a periodic domain, x, y ∈ Ω = [0, 2π] 3 .

The initial condition for the simulation is We depict in 4.11 the numerical velocity obtained and the corresponding anisotropic adapted mesh of 800k elements. 4.6 depict the evolution of the resolution time (left) and the number of iterations needed to solve the linear system (right). As we can see, the obtained results show a significant drop in the computational time and number of iterations required: for a mesh with 10 6 elements, the resolution time 94 is reduced by 47% for Taylor-Green 2D (TG 2D) and 40% for Taylor-Green 3D (TG 3D), the number of iterations is reduced by 75% for TG 2D and 65% for TG 3D. The L 2 norm of the approximation error and the Euclidean norm of the final residual norm are reported in Figures 4.4 and 4.7, the use of the adaptive stopping criterion clearly does not effect the accuracy of the resulting solution. In fact, the accuracy level is similar in both cases.

u(x, y) =   sin(x) cos(y) cos(z) 2 √ 3 -cos(x) sin(y) cos(z) 2 √ 3 0   With the following analytical solution u(x, y) =    sin(x) cos(y) cos(z) 2 √ 3 exp(-π 2 t Re ) -cos(x) sin(y) cos(z) 2 √ 3 exp(-π 2 t Re ) 0    ( 4

2D natural convection benchmark

We aim to evaluate the accuracy of the adaptive stopping criterion that we use for solving coupled problems. We compare the performance of the adaptive criterion with a classical one on results found in the literature on the 2D Natural convection benchmark. It consists in solving a classical flow in a cavity with differentially heated vertical walls. Natural convection in closed enclosures is due to the presence of a temperature gradient. This gradient is manifested by a non-uniform distribution of the temperature in the domain. Under gravity effects, the fluid moves. Natural convection in rectangular cavities has been the subject of numerous experimental and numerical studies, ( [44], [45]), [46], its a classical problem to test and validate many techniques and numerical models. In this type of enclosure, the two vertical walls are heated at a constant and uniform temperature, while both other walls are adiabatic. The left wall is maintained at a temperature T h and the right wall at a temperature T c where T h > T c . The heated walls are parallel to the Y axis, the adiabatic walls are in the direction of the X axis. The Figure 4.8 illustrates the geometric configuration studies.

The Natural convection flow in assumed to be laminar. It is governed by conservation equations mass, amount of motion and energy. The considered fluid is assumed to be newtonian, incompressible and viscous and obeys the Boussinesq approximation [47] ρ

= ρ 0 [1 -β(T -T 0 )] (4.23)
where: ρ 0 : density at T 0 , β: thermal volumetric coefficient and β > 0 in general. The Boussinesq approximation allows to take the effect of thermal dilatation only as a source term in the momentum equation. Note that here the fluid motion and the temperature field are fully coupled: the incompressible Navier-Stokes equations for fluid flow are solved coupled with the thermal convection-diffusion-reaction Equations.

To assess the accuracy of the obtained numerical solution, we compare the obtained Nusselt number with very accurate reference solutions taken from the literature. Furthermore, we will check and highlight the efficiency of the error estimator used in the stopping criterion for coupled problems. The Nusselt number is given by:

L = 1 ∂Ω ∂Ω bottom ∂T ∂n = 0 ∂Ω top ∂T ∂n = 0 ∂Ω lef t ∂Ω right g T h T c
N u = L 0 - ∂T ∂x dy at x = 0 (4.24)
Through Table 4.1, we can clearly identify that the results are in good agreement with the literature references though few differences can be detected. Figure 4.9 highlights how the temperature and the flow responds to expectations. In fact, an upward flow is observed along the left wall, and a downward flow along the the right wall(see Figure4.9), the presence of horizontal walls which are waterproof forcing the fluid to turn right at the top of the hot wall and symmetrically to turn left at the bottom of the cold wall. In other words, the fluid motion inside the cavity is a consequence of different densities between the hot and the cold walls. A sensitivity study of the results to the stopping criterion used (fixed vs adaptive) was carried out. The Table 4.1 reports the data corresponding to the averaged Nusselt number computed on the cold wall shows that the obtained results, with both stopping criterion, are very close to the literature results. Furthermore, the gap between the values of the averaged Nusselt number obtained with the fixed and the adaptive stopping criterion is about 0.1%. In the other hand, Figure 4.13 shows that both stopping criterion produce similar results in terms of temperature convergence in function of time. This highlight one again the potential of the adaptive stopping criterion in not deteriorating the accuracy aspect of the computation. In addition, we perform some investigations to determine some other advantages of the adaptive 96 Ra 10 stopping criterion. As reported in Figures 4.12, the use of the adaptive stopping criterion allows to drastically reduce the computation time and the number of iterations needed to solve the coupled 2D heat transfer and fluid flow problem. For a mesh made up from 10 5 elements the resolution time of all solvers involved in the computation is reduced by 47% and the number of iterations by about 70%.

3D natural convection benchmark

We continue our investigation by considering the three dimensional case of this classical benchmark on the cubic domain [0, 1] 3 . No slip boundary conditions are imposed at the walls, Dirichlet boundary conditions are imposed on the velocity field. The left wall is kept at a high temperature T h whereas the right wall in maintained at a cold temperature T c and horizontal walls are adiabatic (see Figure 4.14). Numerical experiments are done by discretizing the cavity to 293406 nodes for Ra = 10 6 (Ra 99 is the Reynolds number). The mesh is adapted on the normalized velocity field

U i max j |U j | .
The solution is evolved using a time step equal to 0.5s.

In Figure 4.15 we presents different close-ups on the mesh, corresponding to Ra = 10 6 , to show how the elements can highly stretched along the directions of the layers. This reflects well the accuracy and details of the resolution caused by the discontinuity of the boundary conditions and the nature of the flow. For high Rayleigh number Ra = 10 8 the complexity of the flow is clearly shown by the streamlines displayed in Figure 4.17. When the Rayleigh number increases, the flow loses its symmetry indicating that it becomes unsteady. Indeed, the observed flow trajectory manifest a qualitatively new feature by passing from one side to the other of the enclosure and, consequently, by losing its mirror symmetry. As we can see through the front view if Figure 4.17 the spiralling flow structure explores the full width of the cavity. Although, now quantitative measurements were reported in the given references, similar behaviour of the velocity was reported in [START_REF] Zhou | A 3D thermal LB model on non-orthogonal grid and its application for natural convection in irregular domains[END_REF][START_REF] Ugurlubilek | 3D laminar natural convection in a cubical enclosure with gradually changing partitions[END_REF].

Again, in 3D where the flow characteristics are much more complicated, the solution quality is not destroyed and the computational cost is significantly reduces. In Table 4.2 we summarize results of the computational cost and solvers iterations obtained with the fixed SC and the adaptive one for Ra = 10 6 and an anisotropic mesh of 1668464 elements. We can clearly identify a significant decrease of both thanks to the adaptive stopping criterion. In Table 4.3 we reported results of the convergence analysis four different grid resolutions. We compute the relative error of the velocity using the L Error norm 5754 50143 109385 1668464

T h T c x z y
∥u F -u A ∥ L 2 ∥u F ∥ L 2 3.13 • 10 -3 4.54 • 10 -4 3.52 • 10 -4 1.28 • 10 -4 ∥∇(u F -u A )∥ L 2 ∥∇(u F )∥ L 2
2.51 • 10 -3 2.91 • 10 -4 2.70 • 10 -4 9.12 • 10 -5

∥u F -u A ∥ L 2 ∥u F ∥ L 2 = Ω (u F -u A ) 2 dΩ Ω u 2 F dΩ 1 2 and ∥∇(u F -u A )∥ L 2 ∥∇(u F )∥ L 2 = Ω ∇(u F -u A ) 2 dΩ Ω ∇u 2 F dΩ 1 2
where u F is the numerical velocity obtained with the fixed criterion and u A is the numerical velocity obtained with the adaptive criterion. As expected, the finer the mesh the lower the error. On the other hand, the new adaptive criterion yield good results when comparing to the "analytical" solution obtained with the fixed criterion. It does not effect the convergence of the method with respect to the number of elements used, resulting is the same accuracy between the compared results.

Flow past a circular cylinder

In this example, a turbulent unsteady flow past a circular cylinder is analysed to verify the proposed adaptive stopping criterion. The main focus of this numerical example is from one hand to offer large diversity of numerical benchmarks and from the other hand it can be seen as an illustration of a flow past a heated ingot inside an industrial furnace. Recall that the objective is to demonstrate the validity of the proposed adaptive framework on a turbulent unsteady flow with a high Reynolds number. Even for relatively simple geometries, in this case a cylinder, simulating such flows and the loading imposed on the bodies is a difficult task and is computationally considerably more expensive and time consuming.

Physically, when the Reynolds number exceeds some critical value, vortex shedding past a circular cylinder appears naturally and immediately. Numerically, vortex shedding requires a long run-time and is proportional to the Reynolds numbers. It was shown that the quality of vortex-shedding prediction depend strongly on the turbulence model used [START_REF] Zang | A Dynamic Mixed Subgrid-Scale Model and Its Application to Turbulent Recirculating Flows[END_REF], [START_REF] Rodi | On The Simulation of Turbulent Flow past Bluff Bodies[END_REF]. Full details concerning this test using different turbulence models are given by [START_REF] Shao-Shi & Jia | Three-dimensional numerical simulation of the flow past a circular cylinder based on LES method[END_REF], [START_REF] Patel | Numerical Investigation of Flow Past a Circular Cylinder and in a Staggered Tube Bundle Using Various Turbulence Models[END_REF]. Here, we would like to investigate the contribution of the adaptive stopping criterion in term of computational cost. The computational domain for this example as well as the boundary conditions are illustrated in Figure 4.18. The computational domain is discretized by 658286 elements and 112724 nodes. The inflow velocity is constant and given equal to one and the Reynolds number is set equal to different values, 10 3 , 10 4 , 10 5 , 10 6 , so that the sensitivity of the adaptive stopping criterion to an excessive production of turbulent kinetic energy can be tested. We choose the time step of 0.001, which is found to be sufficient to track unsteady characteristics of the flow. Our aim through this example is to investigate the sensitivity of the adaptive stopping criterion to increasing Reynolds numbers. In fact, the advantage of using the adaptive stopping criterion is that it could decrease considerably the computational cost as well as solver iterations, see 4.20, without affecting the solution accuracy as depicted in Figure 4.21. Indeed, as you can see through Figure 4.21, the evolution of the mean value of drag coefficient is almost identical in both cases. This confirms again the validity of the implemented adaptive stopping criterion. This plot shows how the drag coefficient decreases as the Reynolds number increases. This behaviour is expected because the drag is due to skin friction which dominates at low Reynolds numbers. What is also important to note is that, even if the adaptive stopping criterion does not allow going too far in the resolution (see Figure 4.21 right), the quality of the solution is not destroyed as explained previously.

Conclusion

In this chapter, we proposed an adaptive stopping criterion for iterative solvers in the framework of anisotropic adaptive finite elements, applied to convection diffusion reaction and Navier-Stokes equations. The proposed framework is designed from a posteriori error indicator allowing to refine the mesh where the error is the greatest thanks to the estimations of the approximated solution. The idea is to stop solving when the algebraic error is bellow the level of the estimated approximation error.

All numerical tests performed demonstrate that our adaptive strategy based on an early stopping of solver's iterations effects neither the precision nor the quality of the obtained numerical solution. The results proved that, following the provided approach, the total number of iterations needed can be reduced by more than 50%, without no significant effect on the accuracy of the computed solution. More realistic numerical tests based on industrial data will be presented in the next Chapter 5 to illustrate the efficiency of the developed numerical tools in resolving large scale and long time realistic industrial applications.

Résumé du chapitre en français

Introduction

In this chapter, we focused our attention on the simulation on the quenching process on an industrial point of view. In fact, we are interested in evaluating the potential of the the developed numerical tools in resolving long time computation 3D enclosures. The study will combine thermal cooling inside quenching pools by either natural or forced convection. We recall that this work was developed in the context of the ANR industrial chair INFINITY gathering twelve industrial partners whose objective is the establishment of a software allowing to get an accurate and efficient simulation of quenching for arbitrary geometries and conditions. The validity of the framework thus established is demonstrated by comparing the obtained results of industrial quenching simulation through high performance computing with experimental data.

In this work, a multiphase framework allowing the consideration of all scales and all the physics of the problem is proposed. The Level Set method is adopted to track solid/liquid/gaz interfaces, we describe it in the following subsection. This method is combined to an a posteriori error estimator allowing an anisotropic mesh adaptation and the control of linear solvers used by deriving an adaptive stopping criterion. A mixed formulation of finite elements stabilized with a Variational Multiscale method (VMS) (cf section 2.3.3) is considered allowing an accurate resolution of turbulent flows with a high Reynolds number without explicitly resolving small-scales. The surface tension, allowing to describe the separation of the vapour film and thus the hydrodynamic behaviour of the system will be developed later. A simplified model of the phase transformation will also be presented. Finally, the high performance computing environment, resolution algorithms and the computation resources are also given.

The Level Set method

The interface tracking methods allow to take into account complex topological changes that may appear at the interface of multiphase flows. In this work, the Level Set method implemented in CIMLIB CFD finite element library developed at CEMEF is used. The Level Set function, denoted by α, it is a signed distance function whose zero isovalue corresponds to the interface Γ between two domains. Its definition is as follows:

α(x) =    -d(x, Γ) if x ∈ Ω 1 0 if x ∈ Γ d(x, Γ) if x ∈ Ω 2 (5.1)
Ω is the entire domain, x is a point of this domain, Γ is an interface between two sub-domains Ω 1 and Ω 2 .

The resolution of a convection equation allows to predict the time evolution of the interface in a velocity field,

∂α ∂t + (u • ∇)α = 0 (5.2)
u is the velocity field obtained by the resolution of Navier-Stokes equations.

The application of this method to the simulation of multiphase flows is not immediate because of some defects. For example, dissipative numerical errors in the resolution of the transport equation lead to mass losses all the more as the velocity field is under-resolved. It is necessary to reinitialize the Level Set by solving the following equation:

∂α ∂t + ∇(α(1 -α)n -ϵ(∇α • n)n) = 0 (5.3)
An important property of this function is that the computation of geometrical features (the normal ⃗ n and the curvature κ) of the problem is particularly easy.

     n = ⃗ ∇α ∥ ⃗ ∇α∥ κ = ∇ • ( ∇α ∥∇α∥ ) (5.4) 0 0 0 0 0 0 -1 -1 -1 -1 -1 -2 -2 -2

Immersed Volume Method for quenching

In this section, we will introduce the Immersed Volume Method (IVM) used to embed solid geometries inside fluid domains. Indeed, this method combined with anisotropic mesh adaptation, is intensively used for fluid modelling in the context of multiphase flows and for fluid-structure interactions in the context of heat and mass transfer [1][2][3][4][5]. The proposed approach demonstrates its efficiency both for 2D 2D and 3D unsteady heat transfer and turbulent flows inside complex geometries in the presence of conducting solids and can be applied in a wide range of multi-material applications. In the context of quenching process, problems and simulations can be viewed as a fluid structure interaction whose resolution necessitates a technique for modelling the exchanges between the solid and the fluid parts. The obtained results will be used by the industrial partners in the view of improving microstructure and thermo-mechanical properties of the treated workpieces.

The key feature of the IVM method is to combine the level Set method, presented in the previous section, to delimit the different subdomains and the anisotropic mesh adaptation described in Chapter 3 to provide a high resolution at the interface and the mixing laws. The considered subdomains are treated as a single fluid with different material properties, and thus the solution is evolved with a single set of equations on the global domain. The governing equations are solved using the methods presented in Chapter 2 with minor modifications for the non-constant aspect of the thermal material-properties. The anisotropically adapted mesh have to be very fine at the interface level in order to permit a high accuracy in the distribution of the physical properties, without forgetting a good capture of the flow changes and heat transfer with an affordable computational cost. We will introduce in the section bellow the third main component of the IVM, the mixing laws, the two first ones were detailed previously.

The mixing laws

Once the mesh is well adapted around the zero isovalues of the level Set function to have an accurate capture of the solid/fluid interface, we resort to use the mixing laws in order to well distribute the physical properties of the different subdomains.

The assignment of the physical properties is based on the level Set function's sign through a Heaviside function defined by:

H(α) = 1 if α ≥ 0 0 if α < 0 (5.5)
However, defining the Heaviside function as above might be devastating for the numerical solver because of sharp variations across the interface. To that purpose, 120 the idea is to smooth the Heaviside function defined above and, consequently, the physical properties, by artificially thickening the interface layer [6]

H ϵ (α) =        1 if α > ϵ 1 2 α ϵ + 1 π sin πα ϵ if |α| ≥ ϵ 0 if α < ϵ (5.6)
where ϵ is a thickness parameter. In numerical experiments, we take ϵ = 2h d where h d is the mesh size in the normal direction to the interface. Now, we can, take advantage of the regularity of the smoothed Heaviside function to define physical properties of the heat transfer and fluid flow problems such as the density ρ, the temperature T , the dynamic viscosity µ, the heat capacity C p and the conductivity κ as follows:

ρ = ρ f H ϵ (α) + ρ s (1 -H ϵ (α)) µ = µ f H ϵ (α) + µ s (1 -H ϵ (α)) κ = κ f H ϵ (α) + κ s (1 -H ϵ (α))
(5.7) solid properties are referred using the subscript 's' and the fluid ones using the subscript 'f'. Figure 5.2 illustrates the distribution of the conductivity over the domain, obtained using the above described mixing laws, on a gradually refined mesh obtained using layer-based approach (left) and on an anisotropic adaptive mesh (right) using the approach described in 3. Both meshes are made up of 200K nodes. We can clearly see how the transition from the solid subdomain to the fluid one is precisely captured using the anisotropic meshing tools. We note that if the sharp discontinuity of the material properties across the interface is not properly represented, the quality of the numerical solution will be effected. The Immersed Volume Method can be applied to any heat transfer and fluid flow problem without resorting to experimental investigations in order to determine the heat transfer coefficients, the only needed parameters are the physical properties of the different material involved. Numerical simulation of quenching is still challenging because of the wide range of scales to consider and the physics involved. From the nucleation of a vapour bubble to the film boiling, one needs to consider scales for µm to the m. In the case of industrial quenching process, three phases interact with different material properties: the immersed solid, the vapour and the liquid since a vapour film surrounds the part because of a high thermal gradient at the interface. Indeed, this variety of configurations and the complexity of the surrounding flow must be taken into account because they have direct impact of changing mechanical properties and controlling the workpieces' microstructure. Therefore, the study of two-phase flows with phase change, thermal transfer and solid-fluid interactions is a first step to a better understanding of quenching process. The classical well know method of thermal treatment of a solid is based on the use of the heat transfer coefficient experimentally deduced. This technique allows to simply deduce the cooling of a solid, without any knowledge of the fluid behaviour without giving information on the underlying physics. Furthermore, this approach is only useful for a particular configuration in term of geometry, range of temperature and fluid motion. Only a direct numeri-cal simulation can take into account all the physics and all scales involved in boiling.

The majority of approaches devoted to the simulation of boiling depend of the thickness of the interface. A signed distance function (Level Set) and a Volume fraction (Volume-Of-Fluid) is conducted using a transport equation [7][8][9][10][11][12]. These methods are popular due the simplicity of implementation and the fact that complex topology changes are handled naturally. Other approaches in the literature can be mentioned like Phase-field that seems the most promising one [13][14][15]. This is based on an order parameter which is used to distinguish the different phases, and its evolution gives access to the interface location. The interface is assumed to have a thickness of few nanometers.

To be able to correctly model quenching process, several challenging questions have to be tackled. On the one hand, the complexity comes from the modelling of the liquid/vapour phase change at the interface. On the other hand, the surface tension force at the liquid/vapour interface playing an important role in the mechanical behaviour of this interface. Thus, we use phase change model developed in CEMEF CFL research group to predict precisely the liquid to vapour phase change during boiling and the semi-implicit time integration of the surface tension force designed to circumvent the time step restriction. This combined with the Variational Multiscale Method for the Navier-Stokes equation, developed in Chapter 2.3, accounting for the mass transfer between the phases at the interface, the a posteriori error estimator for dynamic anisotropic mesh adaptation addressed in Chapter 3.6 and all numerical tools previously presented allow to deal with immersed solids as in quenching process.

Phase transformation during quenching process

During the quenching process, the workpiece materials will undergo several complicated phenomena as thermal expansion, elastic and plastic deformations as well as phase transformation depending on the cooling rate and its temperature evolution. The study of microstructural evolution through phase transformation allows to make the considered solid with better mechanical properties, microstructure, and an accurate geometry. The set up of a high fidelity elasto-plastic model to give a real description of the deformation combined with a phase transformation model, affecting the temperature distribution as well as the stresses and deformation within the solid, is required. This subject has obtained more attention in the pioneering works of Inoue and Raniecki [16], Rammererstorfer et al [17], Fernandes et al [18], Agarwal and Brimacombe [19], Nagasaka et al [20]. The objective was to model the phase transformation in order to improve the accuracy of numerical simulations. In this study, a thermo-elastic-plastic FE numerical phase transformation model developed within the team was employed to give a real description of the solid phase 123 transformation. This model will be coupled with: 1-an adaptive Eulerian framework used to simulate both boiling and evaporation phenomena (cf previous section) occurring at the interface of the heated solid immersed in the quenching bath, 2-an Immersed Volume Method (AIM) (cf Section 5.1.2) allowing to work on two different domains simultaneously, a fluid-solid domain to understand all the physics and the interaction between the fluid and the solid, and a Solid domain to have a follow-up of the solid state. This hybrid model, will be used to simulate 3D industrial real test cases that we will study subsequently.

Quenching of a rectangular solid

As a first validation of the accuracy and efficiency of the developed tools on industrial applications, we consider an air cooling problem of a rectangular solid by forced convection. The solid having a thickness of around 20cm is initially heated up to 150 • and placed inside the quenching chamber (see Figure 5.3). Air at 20 • is pumped into the enclosure at a constant stream of 1m/s from a ventilation injectors positioned at the top cavity wall. The governing equations are solved with no-slip isothermal boundary conditions u = 0 and T = 10 • on ∂Ω, a zero-pressure and a free-slip velocity outflow are imposed in the outlet. The material properties of the air and the solid are updated dynamically in accordance with their temperature. The computational domain, the initial mesh are depicted in Figure 5.4. The parameters used for the simulation are provided in Table 5.1. A comparative study on resolution time and the total number of linear solver iterations is reported in Table 5.2. We identify that the results show again the ability of the proposed framework with adaptive stopping criterion to deal with such challenging test cases. We report in Figure 5.5 the history of temperature at different locations in the workpiece. A good agreement can be noted between the results computed with the fixed stopping criterion and the adaptive one without a need to adjust the physical parameters. We point out that the selected sensors are the most critical ones since two of them are very close to the interface and reflect the accuracy of the heat transfers at that level and the third one in at the heat of the solid.

µ[Pa.s] ρ[kg/m 3 ] C p [J/(kg.K)] K[W/(m.K)] Fluid 0.

Quenching of a solid and comparison to experimental data

We will consider the quenching of a solid into water. We consider a domain of size 4. Finally, Figure 5.10 shows a good agreement for the temperature evolution between the experimental data and the numerical simulation without the need to adjust the physical parameters. Moreover, the strong coupling between the solid, the water and the vapour enables to perform such simulation without use of the heat transfer coefficient at the solid boundaries. The anisotropic boundary layer mesh enables to capture the high thermal gradient, the motion of vapour film as well as the complexity of the flow. In this simulation, as explained in Chapter 2, a single set of equations is solved simultaneously for both fluid and solid domains. The Variational Multiscale approach is employed to stabilize the Navier-Stokes equations. The SUPG method is used to preclude numerical oscillations at the locations of the convection domination and sharp gradients especially in the vicinity of the interface (Chapter 2). The Immersed Volume Method is used to immerse the solid object into the domain and distribute the material properties which are updated dynamically in accordance with their temperature.

In Table 5.4 we reported the simulation results obtained, in terms of iterations and resolution time, using an anisotropic pre-adapted mesh in the narrow layer around the solid interface. As we can see, when we use the proposed adaptive stopping criterion, the total number of Krylov solver decreases as well as the CPU time.

We can also see, that the temperature evolution (see Figure 5.11) using our error error estimator to build the adaptive stopping criterion is in quantitative agreement with its evolution with the fixed criterion. In other words, the used of the proposed stopping criterion does not effect the convergence of the solution because the accuracy is the same between the compared results.

Iterations: F 412190 CPU: F 34522.09 Iterations: A 33551 CPU: A 15095.6 Reduction (%) ≈ 92% Reduction (%) ≈ 56%

Conclusion

In this Chapter, we investigated the applicability of the developed framework and numerical methods presented in previous Chapters to the modelization and the simulation of quenching industrial process in tarpaulins and simplified cold enclosures. The complexity of the process, geometries, the high ratio of physical properties make such a simulation difficult in terms of modelisation and computational cost. The understanding of the involved physical phenomena will allow our industrial partners to control thermal treatment sequences in order to reduce energy consumption and improve the quality of their products in terms of microstructure, flexibility and resistance. We note that in this case the prediction of boiling process is very important for determining the thermal variation and heat transfer at the fluid/solid interface.

In fact a vapour film is formed around the material with quenching properties weaker than those of water and hence if taken into account, they would lead to numerical solution that is closer to the experimental one.

The objective through this study is to demonstrate the applicability of the developed approach to realistic full-scale quenching simulations. The iterative solver coupling Navier-Stokes and Convection-diffusion reaction equation was run with adaptive stopping criterion based on edge-based a posteriori error estimator. The computational time was shown to be significantly reduced while a good quality of the discrete solution with respect to experimental or reference data was maintained: this is promising for improving the numerical accuracy and efficiency of computations and simulations of industrial applications in engineering.

Résumé du chapitre en français

Dans ce chapitre, nous avons étudié l'applicabilité du cadre développé et des méthodes numériques présentées dans les chapitres précédents à la simulation du procédé de trempe industriel. La complexité du processus, les géométries, les propriétés thermo-physique du solide et le transfert de chaleur à l'intérieur des enceintes rendent une telle simulation difficile en termes de modélisation et de coût de calcul. La compréhension des phénomènes physiques intervenant permettra à nos partenaires industriels de contrôler les séquences de traitement thermique afin de réduire la consommation d'énergie et d'améliorer la qualité de leurs produits en termes de microstucture, de flexibilité et de résistance. La prévision du processus d'ébullition est très importante pour déterminer la variation thermique et le transfert de chaleur à l'interface fluide/solide. En effet, un film de vapeur se forme autour de la pièce immergée avec des propriétés de trempe plus faibles que celles de l'eau. Par conséquent leur prise en compte, mèneraient à une solution numérique très proche de l'expérimentale. Le solveur itératif couplant l'équation de la chaleur (CDR) avec les équations de Navier-Stokes a été exécuté avec le critère d'arrêt adaptatif. Les résultats montrent que le temps de calcul est considérablement réduit. La qualité de la solution discrète est en adéquation avec les données expérimentales. Cette étude est une preuve de l'efficacité et de la précision des outils développés sur les applications industrielles de trempe.

Conclusion

This thesis was motivated by the need to improve long time and large-scale simulations' computational cost of complex immersion problems inside quenching baths including fluid-solid interaction phenomena. Such simulations involves turbulent flow, conjugate heat transfer (convection and conduction) in multi-component formulations. These complex physical phenomena simulation is time consuming and expensive and can go on for weeks to months. To overcome this issue, we relied on a posteriori error estimators allowing to reach a precision given by the user and this with a minimal resolution cost.

We started in Chapter 2 by introducing the numerical tools adopted for the modelling and the simulation of conjugate heat transfer and turbulent flow problems. Standard Galerkin formulations of convection-diffusion reaction equation and incompressible Navier-Stokes equations and stabilization techniques to deal with numerical instabilities. The Streamline Upwind Petrov Galerkin (SUPG) to deal with numerical oscillations in convection or diffusion dominated regimes and Variational Multiscale Method (VMS) to deal with the choice of the finite element spaces necessary for flow problems. In Chapter 3 an a posteriori error estimator has been proposed and has been expressed in function of the recovered gradient of the numerical solution. Different methods proposed in the literature to recover the solution's gradient and Hessian from the piecewise linear approximation solution defined at the vertices of the mesh has been implemented and compared in terms of convergence rates. The implemented interpolation operators are of two types: those based on a L 2 projection and those based on Least squares minimisation. This edge-based a posteriori error estimator is used to drive the anisotropic adaptation procedure by building a nodal metric field in order to generate the appropriate well adapted mesh describing the features of the numerical solution.

In Chapter 4 we showed through our different numerical tests, that our adaptive strategy based on early termination of our solvers' iterations (linear and non-linear) effects neither the quality nor the accuracy of the obtained numerical solution. Furthermore, the computational time and solver iterations was shown to be reduced significantly by almost half or more while a good quality of the discrete numerical solution with respect to the quantities of interest was maintained. The validation study was done on different test cases with increasing complexities going from 2 and 3D analytical test cases to fundamental 2 and 3D benchmarks in aerodynamic simulations in order to compare with the literature results.

We conclude our validation study by applying the proposed framework to largescale phenomena, see Chapter 5, which take place in quenching enclosures. Such a process leads to a high mesh resolution and requires an important computational cost even if a turbulence model is used to account for the turbulence effect on the general flow and hence to accelerate the computation. Therefore, accurate and fast algorithms are necessary for the simulation within a reasonable time of the full cooling process. A comparative study on the temperature evolution inside the workpiece obtained numerically together with experimental results reveals a good agreement and this in a reduced computational time. Therefore, even when complex industrial applications are involved a drastic decrease in computational cost can be reached thanks to the adaptive stopping criterion.

Perspectives

Several perspectives on this work emerge naturally. First of all, it is important to prove the usefulness of the proposed framework for more challenging problems such as the simulation of the multiphase flow with phase changes. In fact such simulations, allow to predict and control the microstructure of the workpieces. The objective is to master final geometries of workpieces and production chains in terms of energy consumption and pollutant emissions. The second perspective is linked to the CPU time of error estimation. To avoid an unnecessary burden on computation, it may be interesting to compute error estimation only every n iterations, whenever the residual is less than a given value. Clearly, the framework developed in this work is far from perfect and a number of other considerations have to be taken into account to achieve an accurate, efficient and robust tool usable for HPC-CFD complex problems. Classical methods used to approximate partial differential equations are constrained by the disadvantage of having fixed time steps. However, in many situations, a physical variable (velocity, pressure, temperature, ...) can change quickly a while, then more slowly after that. Consequently, it is advisable to adapt the time step according to the needs, decrease it when the variations are fast and to increase it in the other case, in order to save computational cost; this requires a deep knowledge the the physics of the system. This is the philosophy of adaptive time-steps' algorithm, interesting to adopt alongside adaptive stopping criterion.

For future works, it is important to prove the usefulness of the proposed framework for the simulation of quenching and cooling processes where the effects of heat transfer by thermal radiation and phase changes are taken into account. Indeed, the radiative exchanges play an essential role on the temperature distribution and the physical properties of the flow; phase change allows to predict the material microstructure. On another level, despite the wealth of literature on a posteriori error estimators, little works have been done on contact problems between objects. It would be interesting to apply our methodology to this type of problems. It would also be interesting, later, to extent our study to contact problems between two vibrating membranes. A particular difficulty of biphasic problems is to prove the existence and the uniqueness of the weak solution. This is possible in the context of biphasic flows with one component per phase. For biphasic multicomponent flows it, however, remains very complicated to do. Furthermore, it seems interesting to try formulating a posteriori error estimators for multicomponent biphasic flows.

The problem of computational time can be approached from another angle. Indeed, as conventional parallelization algorithms (domain decomposition into spaces) become less and less effective when they are used in a massively parallel context. The High Performance Computing (HPC) community has been studying for about ten years the possibility to integrate time parallelization as a new approach of parallelization. The Parareal algorithm with spatial coarsening seems to be able to correctly represents physical properties of turbulent flows. It has a good ability to represents all turbulence scales despite the use of a coarse solver not being able to represents the full spectrum of turbulence. Furthermore, it extends parallel efficiency by reducing the volume of inter-process communication.

Résumé du chapitre en français

Cette thèse était motivée par la nécessité de réduire le coût de calcul des simulations en temps long et à grande échelle des problèmes complexes de trempe. Ces simulations sont caractérisés par des écoulements turbulents, un transfert de chaleur conjugué (convection et conduction) dans des formulations à plusieurs composantes. La simulation de ces phénomènes complexes est longue et coûteuse et peut durer des semaines voire des mois. Pour surmonter ce problème, nous nous sommes appuyer sur les estimateurs d'erreurs a posteriori permettant d'atteindre une précision donnée par l'utilisateur et ceci avec un coût de résolution minimal.

Nous avons commencé dans le Chapitre 2 par introduire les outils numériques adoptés pour la modélisation et la simulation des problèmes conjuguées de transfert de chaleur et d'écoulement turbulent. La formulation standard de Galerkin de l'équation de convection-diffusion réaction et les équations incompressible de Navier-Stokes ainsi que les techniques de stabilisation pour traiter les instabilités numériques. La Streamline Upwind Petrov Galerkin (SUPG) ainsi que la méthode Variational Multiscale (VMS) ont été utilisées pour faire face aux oscillations numériques des régimes à convection dominante.

Dans le Chapitre 3, un estimateur d'erreur a posteriori calculé le long des arêtes du maillage a été décrit, il dépend du gradient reconstruit de la solution numérique.

Différentes méthodes de reconstruction du gradient proposées dans la littérature ont été implémentés et comparées en termes, de taux de convergence. Les techniques d'interpolation étudiées sont de deux types: celles basées sur une projection L 2 et celles basées sur une minimisation des moindres carrés. Cet estimateur d'erreur a posteriori calculé le long des arêtes du maillage est utilité pour conduire la procédure d'adaptation de maillage anisotrope en construisant une champ de métrique nodal permettant de générer un maillage approprié bien adapté décrivant les caractéristiques de la solution numérique.

Dans le Chapitre 4, nous avons démontré à travers différents exemples numériques que la stratégie adaptative proposée basée sur un arrêt précoce des itérations de Krylov du solveur utilisé n'affecte ni la qualité ni la précision de la solution numérique obtenue. En outre, il a été démontré que le temps de calcul et les itérations du solveur étaient considérablement réduits (de plus que la moitié) avec un bon maintien de la qualité de la solution numérique par rapport aux quantités d'intérêt. L'étude de validation a été réalisée sur différents cas tests avec des complexités croissantes allant de de cas tests analytique 2 et 3D à des benchmarks fondamentaux en simulations aérodynamiques 2 et 3D afin de comparer les résultats obtenus avec ceux de la littérature.

Nous avons conclu notre étude de validation en appliquant le cadre proposé à des simulations réalistes de trempe, voir Chapitre 5. De telles simulations nécessitent des résolution de maillage élevée et un coût de calcul conséquent même si un modèle de turbulence tes utilisé pour tenir compte de l'effet de la turbulence sur le flux général et donc d'accélérer le calcul. Par conséquent, des algorithmes précis et rapides sont nécessaires pour simuler le processus de refroidissement dans un délais raisonnable. Ainsi, on a pu observer que même lorsque des applications industrielles complexes sont impliquées, une réduction drastique du coût de calcul peut être atteinte grâce au critère d'arrêt adaptatif.

A.0.2 General Minimal Residual method (GMRES)

The GMRES method was introduced by Saad and Schultz [2]. It is based on a least squares minimization of the residual on the Krylov space orthonormalized by Arnoldi method.

The Arnoldi method and Krylov space

We call Krylov space K m (A, r) the vector space generated by the successive products of the vector r by the matrix A until the power m -1 [1].

K m = Span r, Ar, A 2 r, ...., A m-1 r (A.16)

It is in this space that the algorithm GMRES will look for an approximated solution of the linear system Ax = b for r = b -Ax 0 , where x 0 is the initial point of the iterative method. The approximated solution x 0 is always chosen as the null vector, so that the initial residual is equal to the second member.

We obtain an orthonormal basis of this Krylov space with a simple procedure of the Arnoldi orthonormalization, presented in the Algorithm 1 Algorithm 1 The Arnoldi method for K m (A, r) We notice by Hk the Hessenberg matrix of dimension (k + 1) × k composed of coefficient h ij . It is the matrix obtained after k steps of the Arnoldi method. We notice by H k the matrix of size k × K containing the same coefficients. The matrix Hk has one more line than H k , composed only of zeros and a single non-zero term h k+1,k . H k is upper-triangular with a lower diagonal. We notice by V k+1 the column matrix of vector (v j ) j=1,k+1 obtained after k steps of the Arnoldi method. Thus, the fundamental relation of the Arnoldi method is written for i=1,j do 5:

1: v 1 ← r
AV k = V k+1 Hk (A.
h ij = (w j , v i )

6:

w j ← w j -h ij v i 7:

end for This system can be easily solved: The Hessenberg matrix Hm is upper-triangular with a lower diagonal.

This matrix can be made upper-triangular using m Given rotations noticed (Ω i ). Given rotations aim to cancel a term in a matrix, they are describes in [3]. The Ω i rotation cancel the h i+1,i of Hm . The Given rotation matrices are orthonormal, does not therefore change the norm of the minimization problem. Notice Q i = j i=1 Ω i the The algorithm below, presents this GMRES version with m iterations. If, at the end of this iterations, the tolerance on the residual norm is not reached, the algorithm is restarted with taking as a new starting point the approximated solution obtained at the end of iterations. This involves repeating the algorithm bellow by replacing the two first lines by x 0 ← x m and r 0 ← r m . The algorithm can stop before performing the m iterations if the residual norm is smaller that the requested tolerance. In addition, the zero vector is taken as a starting point of the iterative method. The initial residual is then equal to the second member of the linear system.

2. 1 3 △ 4 . 8

 1348 Schematic of turbulent flow scales (left) and three turbulence modelling approaches (right), adopted from[19] [1] . . . . . . . . . . . . 2.2 Temperature evolution at different locations inside the rectangular shape ingot. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Computational domain (left) anisotropic mesh (right). . . . . . . . . . 2.4 Velocity (top) and temperature (bottom) streamlines at different time steps: t = 10s (left), t = 70 (middle) and t = 160s (right) . . . . . . . 3.1 Computation of nodal values for triangular linear finite elements (left) and quadratic (right): ∆ sampling points nodal values determined with the recovery procedure. . . . . . . . . . . . . . . . . . . . . . . . 3.2 The patches S i associated to the node V i for 2D triangular finite elements (left) and 3D triangular finite elements(right). . . . . . . . . 3.External facets associated to the node (left), ▲ patch of cells associated to the node (right) (2D). . . . . . . . . . . . . . . . . . . 3.4 A general FEM mesh (a) and a patch (b). The related patch of a reference node (shown in red). . . . . . . . . . . . . . . . . . . . . . . 3.5 Linear function f (x, y) = x + y: convergence in the L 2 norm of the error in the gradient reconstruction using the three interpolation techniques previously presented. . . . . . . . . . . . . . . . . . . . . . . . 3.6 Quadratic function f (x, y) = x 2 + 2y 2 : convergence in the L 2 norm of the error in the gradient reconstruction using the interpolation techniques previously presented. . . . . . . . . . . . . . . . . . . . . . 3.7 ∥∇u -∇ R u h ∥ L 2 vs ∥∇u -∇u h ∥ L 2 using Zhu-Zienkiewicz, Clément, Arithmetic-mean and Edge-based recovery techniques. . . . . . . . . . 3.8 Convergence of the effectivity index for the two dimensional Laplace problem using: Clément: CL, Arithmetic mean:AM, ER: Edge-Recovery and Zhu-Zienkiewicz gradient recovery: ZZ. . . . . . . . . . 3.9 Irregular mesh used for the Laplace problem. . . . . . . . . . . . . . . 3.10 Laplace problem: the three unstructured meshes used to evaluate recovery techniques performances and effectivity indexes convergence. v 3.11 ∥∇u -∇ R u h ∥ L 2 vs ∥∇u -∇u h ∥ L 2 using Zhu-Zienkiewicz, Clément and Arithmetic-mean recovery techniques. . . . . . . . . . . . . . . . 3.12 Convergence of the effectivity index for the two dimensional Laplace problem using: Clément: CL, Arithmetic mean: AM and Zhu-Zienkiewicz gradient recovery: ZZ. . . . . . . . . . . . . . . . . . . . 3.13 The evolution of error estimator ∥∇ R u h -∇u h ∥ L 2 obtained using different numbers of procs. . . . . . . . . . . . . . . . . . . . . . . . . 3.14 The evolution of efficiency index θ obtained using different numbers of procs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.15 Vertices distribution at the inter-process interface . . . . . . . . . . . 3.16 Edges distribution at the inter-process interface . . . . . . . . . . . . 3.17 Patch associated with the node V i . . . . . . . . . . . . . . . . . . . . 3.18 An example of a collection of edges used to build the length distribution tensor X i . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.19 Mesh adaptation algorithm applied to two immersed bodies: initial mesh (left), Final mesh(right). . . . . . . . . . . . . . . . . . . . . . . 3.20 Evolution of the two immersed solids shapes during the anisotropic mesh adaptation algorithm. . . . . . . . . . . . . . . . . . . . . . . . 3.21 Schematic diagram of the test case. . . . . . . . . . . . . . . . . . . . 3.22 Velocity distribution at different time steps (left) and the corresponding adapted mesh on the velocity field directions (U i ∥U ∥ L 2 ) (right). . . . 3.23 Velocity distribution at different time steps (left) and the corresponding adapted mesh on the normalized velocity field (∥V ∥ L 2 ) (right). . . 3.24 Schematic diagram summarizing the relationship between a posteriori error estimation, anistropic mesh adaptation and adaptive stopping criterion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Comparison of numerical solutions for pressure: error ∥p analyticalp numerical ∥ (left), contours of pressure field (numerical and analytical) (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Comparison of numerical solutions for velocity: error ∥v analyticalv numerical ∥ (left), contours of pressure field (numerical and analytical) (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 2D Taylor-Green vortex: CPU time (left), Krylov iterations (right). . 4.4 2D Taylor-Green vortex: approximation error ∥u -u h ∥ (left), Final residual norm(right). . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 3D Taylor-Green vortex test case: numerical velocity contours (left) and the anisotropic adapted mesh(right). . . . . . . . . . . . . . . . . 4.6 3D Taylor-Green vortex: CPU time (left), Krylov iterations (right). . 4.7 3D Taylor-Green vortex: approximation error ∥u-u h ∥ L 2 (left), Final residual norm(right). . . . . . . . . . . . . . . . . . . . . . . . . . . . vi Schematic representation of natural convection boundary conditions in the square cavity. No-slip conditions are taken on ∂Ω. . . . . . . . 4.9 Natural convection 2D: temperature isotherms at different times, Ra = 10 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Natural convection 2D: velocity streamlines at different times, Ra = 10 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 2D natural convection: anisotropic adaptive mesh on the velocity field at different times, Ra = 10 6 . . . . . . . . . . . . . . . . . . . . . 4.12 Computational time (left) and solver iterations (right) with the fixed vs adaptive stopping criterion . . . . . . . . . . . . . . . . . . . . . . 4.13 Convergence time of the temperature for the fixed (blue) and the adaptive (red) stopping criterion with a structured mesh of 10 5 elements (top left), 10 4 (top right) and 10 5 (bottom). . . . . . . . . . . . 4.14 Thermal validation problem . . . . . . . . . . . . . . . . . . . . . . . 4.15 Zooms on the mesh details inside the 3D cavity. . . . . . . . . . . . . 4.16 Temperature contours when Ra = 10 4 (a, c, e) and Ra = 10 5 (b,d,f)

1 . 1

 11 3D natural convection streamlines for different Rayleigh numbers: Ra=10 6 (top), Ra=10 7 (middle), Ra=10 8 (bottom) on different views: front view(left), side view (middle), oblique view (right) . . . . . . . . 4.18 Computational domain for the flow past a circular cylinder. . . . . . . 4.19 Velocity evolution and the corresponding anisotropic adaptive mesh at different time steps. . . . . . . . . . . . . . . . . . . . . . . . . . . 4.20 Computational time (left) and solver iterations (right) with the fixed vs adaptive stopping criterion. . . . . . . . . . . . . . . . . . . . . . . 4.21 The mean drag coefficient (left) and the final residual norm ∥r∥ l 2 for different values of the Reynolds number. . . . . . . . . . . . . . . . . 5.1 Geometrical representation of the Level Set method: (a)interface representation; (b)Level Set contours for a circular interface. . . . . . . . 5.2 Conductivity distribution in a multi-domain simulation in adapted meshes of 200K nodes based on a gradual refined (left) and on an anisotropic refined meshes (right) . . . . . . . . . . . . . . . . . . . . 5.3 A 2D representation of the set-up of the 3D forced convection (a) and thermocouples' positions (b) . . . . . . . . . . . . . . . . . . . . 5.4 Computational domain, initial mesh and temperature . . . . . . . . . 5.5 Temperature evolution at different locations inside the rectangular ingot. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Simulation results. Contour plot of temperature distribution on the solid domain at different times. . . . . . . . . . . . . . . . . . . . . . Industrial motivation . . . . . . . . . . . . . . . . . . . . . 1.2 Problematic . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Ambition of this work and the state of the art . . . . . . 1.4 Objective of the thesis . . . . . . . . . . . . . . . . . . . . 1.5 Résumé du chapitre en français . . . . . . . . . . . . . . . Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

Figure 2 . 1 :

 21 Figure 2.1: Schematic of turbulent flow scales (left) and three turbulence modelling approaches (right), adopted from [19] [1]

Figure 2 . 2 :

 22 Figure 2.2: Temperature evolution at different locations inside the rectangular shape ingot.

Figure 2 . 3 :

 23 Figure 2.3: Computational domain (left) anisotropic mesh (right).

Figure 2 . 4 :

 24 Figure 2.4: Velocity (top) and temperature (bottom) streamlines at different time steps: t = 10s (left), t = 70 (middle) and t = 160s (right)

3. 5

 5 Numerical validation: error estimator and Effectivity index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .3.5.1 Parallel implementation . . . . . . . . . . . . . . . . . . .3.6 Metric-based anisotropic mesh adaptation . . . . . . . .3.6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3 . 1 :

 31 Figure 3.1: Computation of nodal values for triangular linear finite elements (left) and quadratic (right): ∆ sampling points nodal values determined with the recovery procedure.

Figure 3 . 3 :

 33 Figure 3.3: △ External facets associated to the node (left), ▲ patch of cells associated to the node (right) (2D).

Figure 3 . 4 :

 34 Figure 3.4: A general FEM mesh (a) and a patch (b). The related patch of a reference node (shown in red).

Figure 3 . 5 :

 35 Figure 3.5: Linear function f (x, y) = x + y: convergence in the L 2 norm of the error in the gradient reconstruction using the three interpolation techniques previously presented.

Figure 3 . 6 :

 36 Figure 3.6: Quadratic function f (x, y) = x 2 +2y 2 : convergence in the L 2 norm of the error in the gradient reconstruction using the interpolation techniques previously presented.

  Figure 3.7: ∥∇u-∇ R u h ∥ L 2 vs ∥∇u-∇u h ∥ L 2 using Zhu-Zienkiewicz, Clément, Arithmeticmean and Edge-based recovery techniques.

Figure 3 . 8 :

 38 Figure 3.8: Convergence of the effectivity index for the two dimensional Laplace problem using: Clément: CL, Arithmetic mean:AM, ER: Edge-Recovery and Zhu-Zienkiewicz gradient recovery: ZZ.

Figure 3 . 9 :

 39 Figure 3.9: Irregular mesh used for the Laplace problem.

Figure 3 . 10 :

 310 Figure 3.10: Laplace problem: the three unstructured meshes used to evaluate recovery techniques performances and effectivity indexes convergence.

57

 57

Figure 3 . 14 :

 314 Figure 3.14: The evolution of efficiency index θ obtained using different numbers of procs.

Figure 3 . 15 :Figure 3 . 16 :

 315316 Figure 3.15: Vertices distribution at the inter-process interface

Figure 3 . 17 :

 317 Figure 3.17: Patch associated with the node V i .

Figure 3 . 19 :

 319 Figure 3.19: Mesh adaptation algorithm applied to two immersed bodies: initial mesh (left), Final mesh(right).

Figure 3 . 20 :Figure 3 . 21 :

 320321 Figure 3.20: Evolution of the two immersed solids shapes during the anisotropic mesh adaptation algorithm.

Figure 3 . 22 :

 322 Figure 3.22: Velocity distribution at different time steps (left) and the corresponding adapted mesh on the velocity field directions ( U i ∥U ∥ L 2 ) (right).

Figure 3 . 23 :

 323 Figure 3.23: Velocity distribution at different time steps (left) and the corresponding adapted mesh on the normalized velocity field (∥V ∥ L 2 ) (right).

Figure 3 . 24 :

 324 Figure 3.24: Schematic diagram summarizing the relationship between a posteriori error estimation, anistropic mesh adaptation and adaptive stopping criterion

Figure 4 . 1 :Figure 4 . 2 :Figure 4 . 3 :

 414243 Figure 4.1: Comparison of numerical solutions for pressure: error ∥p analytical -p numerical ∥ (left), contours of pressure field (numerical and analytical) (right)

Figure 4 . 4 :

 44 Figure 4.4: 2D Taylor-Green vortex: approximation error ∥u -u h ∥ (left), Final residual norm(right).

  .22) p(x, y) = 100 + 1 16 [(cos(2x) + cos(2y))(2 + cos(2z)) -2] exp(-π 2 t Re )

Figure 4 . 5 :

 45 Figure 4.5: 3D Taylor-Green vortex test case: numerical velocity contours (left) and the anisotropic adapted mesh(right).

Figure 4 . 6 :

 46 Figure 4.6: 3D Taylor-Green vortex: CPU time (left), Krylov iterations (right).

Figure 4 . 7 :

 47 Figure 4.7: 3D Taylor-Green vortex: approximation error ∥u-u h ∥ L 2 (left), Final residual norm(right).

Figures 4. 2

 2

  and 4.1 allows to compare the solution pair (v, p) obtained numerically (top right) and analytically (bottom right) and the corresponding error (left).Figures 4.

3

 3

  and

Figure 4 . 8 :

 48 Figure 4.8: Schematic representation of natural convection boundary conditions in the square cavity. No-slip conditions are taken on ∂Ω.

Figure 4 . 9 :

 49 Figure 4.9: Natural convection 2D: temperature isotherms at different times, Ra = 10 6 .

Figure 4 . 10 :

 410 Figure 4.10: Natural convection 2D: velocity streamlines at different times, Ra = 10 6 .

Figure 4 . 11 :

 411 Figure 4.11: 2D natural convection: anisotropic adaptive mesh on the velocity field at different times, Ra = 10 6 .

Figure 4 . 12 :Figure 4 . 13 :

 412413 Figure 4.12: Computational time (left) and solver iterations (right) with the fixed vs adaptive stopping criterion

Figure 4 . 14 :

 414 Figure 4.14: Thermal validation problem

Figure 4 . 15 :

 415 Figure 4.15: Zooms on the mesh details inside the 3D cavity.

Figure 4 . 16 :

 416 Figure 4.16: Temperature contours when Ra = 10 4 (a, c, e) and Ra = 10 5 (b,d,f ) obtained by: (a,b) by Li et al.[55]; (c,d) by [54], (e,f ) This study

Figure 4 . 17 :

 417 Figure 4.17: 3D natural convection streamlines for different Rayleigh numbers: Ra=10 6 (top), Ra=10 7 (middle), Ra=10 8 (bottom) on different views: front view(left), side view (middle), oblique view (right) .

Figure 4 . 18 :

 418 Figure 4.18: Computational domain for the flow past a circular cylinder.

Figure 4 . 19 :

 419 Figure 4.19: Velocity evolution and the corresponding anisotropic adaptive mesh at different time steps.

Figure 4 . 20 :

 420 Figure 4.20: Computational time (left) and solver iterations (right) with the fixed vs adaptive stopping criterion.

Figure 4 . 21 :

 421 Figure 4.21: The mean drag coefficient (left) and the final residual norm ∥r∥ l 2 for different values of the Reynolds number.

Figure 5 . 1 :

 51 Figure 5.1: Geometrical representation of the Level Set method: (a)interface representation; (b)Level Set contours for a circular interface.

Figure 5 . 2 :

 52 Figure 5.2: Conductivity distribution in a multi-domain simulation in adapted meshes of 200K nodes based on a gradual refined (left) and on an anisotropic refined meshes (right)

Figure 5 . 3 :

 53 Figure 5.3: A 2D representation of the set-up of the 3D forced convection (a) and thermocouples' positions (b)

Figure 5 . 4 :Figure 5 . 5 :Figure 5 . 6 :

 545556 Figure 5.4: Computational domain, initial mesh and temperature

  7 × 3.1 × 1.75m 3 , three-quarters full of water, in which a circular metallic sample, made of Inconel 718 alloy of radius 0.15m and 0.15m of height, is immersed (see Figure 5.7). The temperature of the sample is T solid = 985 • C and the temperature or water is T water = 25 • C. Several thermocouples are placed at different parts of the workpiece (see Figure 5.8).A free slip boundary condition is prescribed on all the walls: u • n = 0,, where n is the normal to the boundary. The motion of the vapour film is due to the buoyancy force and the surface tension force. The simulations are performed using a boundary layer mesh containing 14472016 cells. The time step is set to ∆t = 0.01s. All the physical parameters related to each phase are presented in Table5.3.

Figure 5 . 7 :

 57 Figure 5.7: Set-up for the 3D industrial quenching(left) and the boundary layer mesh used (right). The solid is represented in red

Figure 5 . 8 :

 58 Figure 5.8: Sensors' position inside the cylindrical workpiece.

Figure 5 . 9 :

 59 Figure 5.9: Industrial quenching. Results at t=1.55s, 6.50 and 12s. Perspective view (left), front view (right). The zero isovalue of the vapour Level Set in blue.

Figure 5 . 10 :Figure 5 . 11 :

 510511 Figure 5.10: Temperature evolution at different locations inside the workpiece: experimental results vs simulation results

  . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Résumé du chapitre en français . . . . . . . . . . . . . . .

1 . 2 .

 12 This choice is motivated by the result of the following proposition. Proposition: Let A ∈ M n (R) as symmetric positive definite matrix and b ∈ R n . Then x ∈ R n is a solution of the linear system[A.4] if and only if∀x ∈ R n , J(x) ≤ J(x) (A.6)where J is the function defined by[A.5].Remark: Remember that the gradient of a functionG : R n -→ R at x ∈ R n verify (∇G(x), z) = lim t→0 G(x + tz) -G(x) tIt follows, considering the function defined by [A.5], that∀z ∈ R n , (∇J(x), z) = (Ax -b, z)and that ∇J(x) = Ax -b We can poof using the following definition of the gradient of a function G :x ∈ R n -→ R at a point a ∈ R n : ∇J(a) = Aa -b Notations: When x ∈ R n verifying [A.6] exists, we say that x is the solution of the following minimization problem inf x∈R n J(x) or x is the solution of the following problemJ( x) = min x∈R n J(x) (A.7)Solving [A.4] is therefore equivalent, according to the above proposition, to find x solution of[A.7].

||r|| 2 : 4 :

 24 for j = 1, m do 3: w j ← Av j for i=1,j do 5:h ij = (w j , v i )

17

 17

) Algorithm 2 1 : v 1 ← r ||r|| 2 :

 212 The modified Arnoldi method for K m (A, r) for j = 1, m do 3:w j ← Av j 4:

min ym∈R m ||βe 1 -

 1 end forThe final residual is as followsr m = b -Ax m = b -AV m y m (A.22)If we take as initial solution x 0 = 0, then the initial residual r 0 used to build the Krylov is equal to r 0 = b. In the Arnoldi method ( Algorithms 1, 2 ), the initial vector v 1 = r 0 is normalized. Let's note β the norm of the initial residual r 0 . Using the Arnoldi relation [A.17] previously introduced notation, we obtainr m = b -Ax m = r 0 -AV m y m = β -V m+1 Hm y m = V m+1 (βe 1 -Hm y m ) (A.23) where e 1 is the first vector of the R m canonical basis. We thus obtain the following relation between the full residual and that obtained in the reduced basis:||r m || 2 = ||βe 1 -Hm y m || 2 (A.24)Indeed, the basis V m+1 is orthonormal, therefore it does not change the L 2 norm of a vector. The residual minimization r m amounts to a simple least-squares minimization problem of size (m + 1) × m: Hm y m || 2 (A.25)

Algorithm 3 8 : 14 : 15 :

 381415 GMRES(m) with restarting1: x 0 ← 0 2: r 0 ← b 3: β ← ∥r 0 ∥ 4: v 1 ← r 0 /β 5: g ← βe 1 6: for j = 1, m do 7:w j ← Av j for i=1,j do 9:h ij = (w j , v i ) 10: w j ← w j -h ij v iComputation of Given rotation Ω j cancelling h j+1,jApplication of Given rotations Ω 1 . . . Ω 1j on H 16:g ← Ω j g 17: end for 18: y m ← H -1 g 19: x m ← x 0 + V m y m 20: r m ← b -Ax m

  and g is a given boundary condition. The scalar function f ∈ L 2 (Ω) is often called source term and can model, depending on the problem, a heat source, a chemical reaction, injection/production wells.

  Zhu-Zienkiewicz, Clément and Arithmetic-mean recovery techniques.

	•10 -2
	1.63
	1.62
	1.61
	1.6
		2	4	6		8
			P rocs
							θ CL
							θ AM
							θ ZZ
		1.04					θ ER
	θ	1.02
		1
		0	0.2	0.4	0.6	0.8	1
				Elements		•10 5
	Figure 3.12: Convergence of the effectivity index for the two dimensional Laplace problem
		using: Clément: CL, Arithmetic mean: AM and Zhu-Zienkiewicz gradient
		recovery: ZZ.

  . It is a periodic analytical solution of the Incompressible Navier-Stokes equations. It describes a set a vortices that amortize over time. The Taylor-Green (TG) vortex was simulated in the unit square, x, y ∈ Ω = [0, 1] 2 with Re = 1000. The Reynolds number Re must be low because chaos and fluctuation can manifest due to vortices instabilities. The initial velocity for the simulation is

		u(x, y) =	-cos(2πx) sin(2πy) sin(2πx) cos(2πy)	(4.19)
	The boundary condition are determined so that the exact solution pair (u, p) is
	described by the following analytical expression
		u(x, y) =	-cos(2πx) sin(2πy) exp(-8νπ2t) sin(2πx) cos(2πy) exp(-8νπt)	(4.20)
	and	p(x, y) = -	(cos(4πx) + cos(4πy)) 4	exp(-8νπt)	(4.21)

Table 4 .

 4 1: Comparison of the Nusselt number obtained in the present work with the literature.

	4	10 5	10 6

Table 4 . 2

 42

	2 -norm and

: 3D natural convection: Results

Table 4 .

 4 3: Convergence results for four different adaptations

H 1 -semi norm as follows:

  Dans ce chapitre, nous proposons un critère d'arrêt adaptatif pour solveurs itératifs dans le cadre des éléments finis adaptatifs anisotropes appliqué aux équations de Convection Diffusion Réaction et Navier-Stokes. Le cadre proposé est conçu à partir de l'indicateur d'erreur a posteriori permettant de raffiner le maillage dans les régions où les variables d'intérêt présentent de fortes variations. L'idée est d'arrêter de résoudre quand l'erreur algébrique est plus petite que l'erreur d'approximation. Les différents exemples numériques présentés et les benchmarks tirés de la littérature illustrent que la stratégie adaptative proposée basée sur un arrêt précoce du solveur itératif n'affecte ni la précision ni la qualité de la solution numérique obtenue. Les résultats illustrent qu'en suivant l'approche adaptative, le nombre total d'itérations nécessaires peut être réduit de plus de 50% sans affecter la précision de la solution.

Table 5 .

 5 1: Initial physical parameters for the 3D forced convection problem.

Table 5 . 2 :

 52 3D forced convection: Resolution time and solver iterations obtained with Fixed vs adaptive stopping criterion 125

Table 5 .

 5 3: Initial physical parameters for the 3D industrial quenching.

	4	2	3
	10	11	14	17
	9	12	15
	1	5	6	7
				37.5
	8	13	16
	15
	50

Table 5 .

 5

4: 3D forced convection: Resolution time and solver iterations obtained with Fixed vs adaptive stopping criterion

Computers & Mathematics with Applications. 68. 10.1016/j.camwa.2014.08.008.

Remerciements

We introduce the set of piecewise constant functions and the approximation set of continuous functions which are piecewise linear.

We denote by S i the support of the basis function ϕ i associated to the node V i .

Parallel implementation

All the interpolation operators described in this chapter have been implemented in the parallel numerical library of the research group using the Message Passing Interface (MPI).

While most computations are related to the finite elements and therefore purely serial on each partition, some quantities introduced in the previous subsections are evaluated on patches of cells incident to a given vertex: in the case of a vertex shared between two or more processors, the patch is split across different partitions. Therefore operations should be performed in a way that information from the patch on one partition should be made available to adjacent processors sharing the patch.

In the first case, the information is cell-based and is only available to one partition; since cells are not shared in the considered framework. For example, in the case of the linear Clément interpolator, the resulting nodal value is a linear expansion of elementwise constant finite element gradients weighed by the cell measure: terms of the linear expansion should be accumulated at the vertex for each sharing processors, as well as the sum of the cell measures used for normalization. This is also the case of Scott-Zhang besides a specific treatment at the domain boundary. A synchronization of the interprocess boundaries using accumulation, similar to the ADD VALUES pattern in PETSc was implemented to allow such computation on a distributed mesh. This is also crucial for the Zhu-Zienkiewicz interpolator where forming the linear system for Least-Squares minimization at each vertex consists of summing elementary matrices evaluated at each sampling point. Given that all the information is available for each adjacent processor, no further synchronization is required (see Figure 3.15).

In the second case, the information is edge-based and therefore available to more than one partition if the edge is shared. While the strategy is similar to the previous case, additional care must be taken to add the contributions only once per edge: only one processor must be responsible for the computation on a shared edge (see Figure 3.16). To this effect, during the construction of the star topology, a function to attribute the ownership of an edge was used: only the owner accumulates the contribution for the given shared edge. This implementation was required by the edge-based recovery method.

To ensure that the parallel implementation is correct, interpolation of gradient of analytical expressions were tested on different numbers of processors and the relative error to the serial solution assessed to be of the order of the truncation error. average representation of the distribution of edges in Γ(i) as in Figure3.18. Note that, X i is symmetric positive definite, and thus invertible, consequently, we get,

Thus, by definition of the trace of a matrix, we have

Then using Equation 3.46 and by identification, we get

M i is an approximation solution of the problem 3.46, this metric allows to transform the mesh edges into unit edges

Mesh adaptation criteria

When we simulate complex physical phenomena involving turbulent flows and conjugate heat transfer, it is highly desirable to accurately capture all the characteristics of the problem, including but not limited to, the flow field, boundary layers, the temperature variations, the flow separation and the fluid/solid interfaces. Two adaptation strategies are generally used: the "a priori " and the "a posteriori " one. The first strategy departs from a pre-adapted mesh around the boundaries and in the wake regions by making an a priori assumptions on the solution behaviour (i.e.

Appendix A Iterative Methods

A.0.1 Conjugate Gradient (CG)

The Conjugate Gradient (CG) method allows to solve linear systems whose matrix is symmetric positive definite. It is a method consisting for a giver vector x 0 to determine at every step a vector p k and a scalar α k allowing to compute x k+1 from

with the aim of minimizing a function. Before introducing the conjugate gradient method to solve a linear system, it is useful to briefly address gradient methods.

A.0.1.1 Gradient methods

We consider the following minimization problem

Where J is a given function from R n to R. In this section, we will chose in [A.1] the gradient of the function to minimize as a descent direction:

Where J : x ∈ R n -→ J(x) ∈ R is the function to minimize. The parameter α k will first be chosen constant than optimal. The aim is to solve the linear system:

where A is a symmetric positive definite matrix of order n and b ∈ R n . We consider the following J function

The gradient method with fixed step size

The fixed step gradient method is defined by the following algorithm

,

where ρ(A) is the spectral radius of A, i.e, the sequence generated by the algorithm [A.8] converges, for every choice of x 0 , toward the unique solution of the linear system:

Remark

In the gradient algorithm described by [A.8,], the parameter r is called the step of the method.

The gradient method with optimal step length

The optimal step length gradient method consists on changing the parameter r in [A.7] at each iteration k by a step length r k , called optimal step length, whose construction is based on the following result. Proposition Let J the function defined on R n by [A.5], x a vector of R n and w a non-zero R n vector. We consider the following problem

If the matrix A involved in [A.5] is symmetric positive definite, then [A.9] has a unique solution given by

Conjugate Gradient method

It is recalled that the optimal step length gradient method is an iterative descent method following the greatest slope in which the passage form x k to x k+1 is done according to the formula

where g k = ∇J(x k ) and r k is the step length at the k ith iteration. The conjugate gradient method consists on taking as a first descent direction w 0 = ∇J(x 0 ) and to chose at the step k + 1 a direction w k+1 such that (w k+1 , Aw k ) = 0 (A.12)

and an optimal step length r k defined as reads

Thus, the descent directions w k and w k+1 are said A-conjugate or orthogonal for the scalar product induced by A.

For a given x k and for all w k ̸ = 0, the parameter r k is defined as follows:

Where g k = Ax k -b is the residual vector at the iteration k or the gradient of J at x k . We use the result of the proposition above to obtain the expression of x k .

Conjugate Gradient algorithm

We initialize the algorithm by choosing a vector x 0 . If g 0 ≡ Ax 0 -b = 0, we stop iterating. If g 0 ̸ = 0, we choose w 0 = g 0 , we take w 0 = g 0 and we execute the following instructions:

(Ag 0 ,g 0 )

The demonstration is done simply by noticing that

We thus obtain

The Arnoli relation is interesting because it shows that the action of the big size matrix A on vectors v i can be reduced to that of a match smaller matrix Hk . The classical size of Krylov space is in the hundred, compared with the matrix A whose order is tens of millions in our case. The Hessenberg matrix H k can be computed only with A and V k :

Indeed, using the Arnoldi relation A.17, and noting I k the identity matrix

The modified Arnoldi method

The previously introduced Arnoldi method performs in just one step the projection of the vector w j on the orthogonal space to the hyperplane of basis V j . We can show, [1,3] that this method is sensitive to numerical errors which can lead to a loss of orthogonality between vectors. The modified form of the Arnoldi method, which performs the projection of the orthogonal of each vector of the basis V j , is more robust to numerical errors. This gives the Algorithm 2.

To further improve the accuracy of the method, we can perform a second orthonormalization after the first one [1]. Finally, a last orthonormalisation method of a set of vectors is based on Householdeer thoughts, is even more accurate than the Arnoldi method, at a higher computational cost [3].

The GMRES algorithm

The GMRES algorithm is based on the minimization, in a least-squares sense, of the residual r = Ax -b on the Krylov space K m (A, r 0 ). A linear combination of basis vectors of the Krylov space is needed. We notice y m ∈ R m its coordinates. The solution x m is as reads

total transformation, rotations product. Thus

where R m is a upper-triangular matrix of size m × m. The same transformation is applied to the second member βe of the minimization problem gives a vector g m+1 of coordinates γ i :

The minimization problem is then rewritten

R m is a upper-triangular matrix, and its invertible (only if the method reach the exact solution, see lucky breakdown notion in [1]). By a triangular ascent method [3], we obtain the solution coordinates

The residual norm is exactly equal to

The result of [A.30] is particularly interesting. Indeed, the solution of the minimization problem can be found without solving the problem. The successive transformation of the second member βe 1 using givens' rotations gives, as a last coordinate, the residual value at the current GMRES iteration, without computing the solution y k = R -1 k g k . For that reason, during GMRES algorithm iterations the residual r and the solution x vectors are not updated before reaching the desired convergence threshold or a restarting point.

ABSTRACT

Quenching is a heat treatment process used to modify the mechanical properties of the forged, moulded or welded metal part. It consists of heating a workpiece to change its microstructure and its properties like hardness, resistance and toughness. The workpiece is then cooled in a medium (oil, water, polymer solution or air). This process is commonly used to harden and reinforce metal parts for the automotive and aerospace sectors such as rings and gears and other transmission parts. It is also used in construction sector to avoid bards' and rods' distortions and in the energy domain (for example, seamless rolled crowns). Nowadays, with the improvement of computing power, the numerical simulation of this process become an essential tool to predict physical phenomena charactering this process such as temperature and cooling velocity, these last two are essential factors allowing to determine final characteristics of the material. Numerical simulation is an excellent tool to understand those results and to optimize them. However, the simulation of such phenomena posed scientific difficulties because their resolution implies long computational times despite the use of important computational resources. In this thesis, we are interested in the resolution of complex long time and large scale problems heat transfer and fluid flow problems. The goal is to offer a general adaptive stopping criteria for each iterative solver used in the aim of reducing the number of iterations and computational time. Those criteria are based on a posteriori error estimators computed on mesh's edges and based on recovery procedures. Those estimators are initially used to lead the anisotropic adaptive process to refine the mesh locally in the regions of interest. They allow to measure the quality of the approximated numerical solution by providing entirely computable upper bounds on the error between the exact solution and the approximated one. Our numerical tests highlight the accuracy of the estimators used and the reduction in terms of iterations' number and computational cost, this reflects the efficiency of our adaptive method. The numerical framework has been validated by confrontations with experimental results provided by our industrial partners.

KEYWORDS

Quenching process, anisotropic mesh adaptation, non-smooth interpolation operators, a posteriori error estimate, adaptive stopping criterion.