Él Éonore Gauci

M Bruno Koobus

Anca Mme

Belme

M Olivier

Allain -Société Lemma

M Frédéric Alauzet

M Alain Dervieux

Marion Merci

Présentée par

Keywords: Metric-based mesh adaptation, anisotropy, adjoint, unsteady simulations, moving mesh, ALE Adaptation de maillage basée métrique, anisotropie, adjoint, simulations instationnaires, mouvement de maillage, ALE Oral Communications

Adaptation de maillage orientée fonctionnelle et basée sur une métrique pour des simulations aérodynamiques en géométrie variable/Goal-oriented metric-based mesh adaptation for unsteady CFD simulations involving moving geometries

Goal-oriented metric-based mesh adaptation for unsteady CFD simulations involving moving geometries Abstract When dealing with CFD problems, mesh adaptation is interesting for its ability to approach the asymptotic convergence and to obtain an accurate prediction for complex flows at a lower cost. Anisotropic mesh adaptation method reduces the number of degrees of freedom required to reach a given solution accuracy, thus impact favorably the CPU time.

Moreover, it reduces the numerical scheme dissipation by automatically taking into account the anisotropy of the physical phenomena inside the mesh. Two main approaches exist in the literature. Feature-based mesh adaptation which is mainly deduced from an interpolation error estimate using the Hessian of the chosen sensor controls the interpolation error of the sensor over the whole computational domain. Such approach is easy to set-up and has a wide range of application, but it does not take into account the considered PDE used to solve the problem.

On the other hand, goal-oriented mesh adaptation, which focuses on a scalar output function, takes into consideration both the solution and the PDE in the error estimation thanks to the adjoint state. But, the design of such error estimate is much more complicated. This thesis presents the results obtained with different CFD methods : the Arbitrary Lagrangian Eulerian (ALE) flow solvers with explicit and implicit schemes are presented and coupled to the moving mesh process, the feature-based unsteady mesh adaptation for moving geometries takes into account the changes of connectivites during the whole simulation, the adjoint state is extended to moving geometries problems and goal-oriented unsteady mesh adaptation for moving meshes is derived from an a priori error estimate. Several numerical examples are considered in the aeronautics sector and the field of civil security.

Convention Simplicial mesh

Let n be the dimension of the physical space and let Ω ⊂ R n be the non-discretized physical domain. Ω is an affine space. The canonical basis of its vectorial space is noted (e 1 , e 2 , . . . , e n ) in general, but notations (e x , e y ) and (e x , e y , e z ) can also be used in two and three dimensions.

Vectors of R n are noted in bold font. The coordinate vector of a point of Ω is generally noted

x = (x 1 , . . . , x n ).
The boundary of Ω, noted ∂Ω, is discretized using simplicial elements the vertices of which are located on ∂Ω. In two dimensions, ∂Ω is discretized with segments while in three dimensions, boundary surfaces ∂Ω are represented by triangles. The discretized boundary is noted ∂Ω h and Ω h denotes the sub-domain of R n having ∂Ω h as boundary.

Building a mesh of Ω h consists in finding a set of simplicial elements -triangles in two dimensions, tetrahedra in three dimensions -noted H, satisfying the following properties:

• Non-degenerescence: Each simplicial element K of H is non-degenerated (no flat elements of null area or volume),

• Covering:

Ω h = K ∈ H K ,
• Non-overlapping: The intersection of the interior of two different elements of H is empty:

Ki ∩ Kj = ∅, ∀ K i , K j ∈ H, i = j .
• Conformity: The intersection of two elements is either a vertex, an edge or a face (in three dimensions) or is empty.

The conformity hypothesis is adopted here, because the meshes will be used with a Finite-Volume solver which requires the enforcement of this constraint. Besides, it facilitates the handling of data structures on the solver side and also enables to save a consequent amount of CPU time.

Finally, a mesh is said to be uniform if all its elements are almost regular (equilateral) and have the same size h.

Notations and orientation conventions

The following notations will be used in this thesis: K is an element of the mesh H, P i is the vertex of H having i as global index, e ij is the edge linking P i to P j . The number of elements, vertices and edges are noted N t , N v and N e , respectively.

The vertices and the edges of an element K are also numbered locally. Vertex numbering inside each element is done in a counter-clockwise (or trigonometric) manner, which enables to compute edges/faces outward normals in a systematic way. This numbering, as well as unit outward normals n and edges/faces orientations are shown for triangles and tetrahedra in Figure 1. Non-normalized normals will be noted η.

Introduction

This manuscript proposes a synthesis of my research during my years of PhD at Inria in Gamma3 and Ecuador teams in the field of anisotropic mesh adaptation applied to unsteady simulations involving moving computational domains.

Often considered as a substitute for modelling systems for which simple closed from analytic solutions are not possible, computer simulations' big advantage is the real ability to generate representive scenarios for which a realisation of all possible states would be prohibitive or impossible in practice. Consequently, numerical simulation has become an integral part of the design process in science and engineering. In sixty years, numerical simulation has come a long way, in conjunction with the development of computing resources. As it became possible to model and simulate more complex problems, industry started to take advantage of its potentiality. However, as the simulation capabilities grow, so do the requirements of both industries and researchers, and a whole set of problems have arisen. In the field of computational fluid dynamics (CFD), the numerical simulations have to deal with ever increasing geometrical and physical complexities. However being able to predict numerically all the features of complex flows around complex geometries remains an unachieved goal. In this sense, scientists and engineers have developed several techniques to make predictions with the maximum of accuracy while mastering the computational efforts in term of CPU time.

In this sense, mesh adaptation methods have been developed to reduce the complexity of the simulations while keeping the same high level of precision. Notably, the increase in computational power allows scientists to not content themselves with steady approximations of the physical phenomena, and to study the effects of time-dependent parameters.

Mesh adaptation for CFD

The computational pipeline for fluid dynamic simulations illustrated by Figure 2 can be summarized as :

Cad → Mesh → Solver → Visualization/Analysis
The mesh is one of the first step of the computation and as the topological base frame, its role may be crucial. Anisotropic mesh adaptation method reduces the number of degrees of freedom required to reach a given solution accuracy and thus impacts favorably the CPU time.

Moreover, it reduces the numerical scheme dissipation by automatically taking into account the anisotropy of the physical phenomena inside the mesh. That's why this technique is nowadays Figure 2: Under the visualization of the numerical solutions of Euler flux solution of f117 in flight on the Left, the mesh of the computational volume around which enables to compute the solution on the Right a great deal of interest. Two main approaches exist in the literature. The first one is the feature-based mesh adaptation which is mainly deduced from an interpolation error estimate using the Hessian of the chosen sensor. It controls the error on the sensor over the whole computational domain. Such approach is easy to set-up and has a wide range of application, but it does not take into account the considered PDE used to solve the problem. The second approach is the goal-oriented mesh adaptation. It focus on a scalar output function and takes into consideration both the solution and the PDE in the error estimation thanks to the adjoint state.

In this thesis, we focus on Computational Fluid Dynamics (CFD), and its applications to vorticity in aeronautic and blast mitigations, although the results presented can be used in a vast field of research and industry. The next section presents its numerical context in details.

Numerical context

The present work was conducted in Gamma3 and Ecuador research groups at Inria. Most numerical choices naturally stem from those made in the group, to benefit from both the considerable experience and the valuable software of its members.

-We consider unstructured meshes made up of triangles (in 2D) or tetrahedra (in 3D). This geometrical choice is initially due to the fact that it is easier to mesh complex geometries with simplexes.

-We consider metric-based mesh adaptation. The theory derives from the classical Riemannian geometry theory and underlies on accurate feature-based or adjoint based error estimates.

-We consider anisotropic mesh-adaptation to fit the natural physical phenomena properties inside the mesh in order to improve their representation.

-Regarding moving-boundary meshes, the body-fitted approach has been preferred in comparison to other methods for its accurate treatment of boundaries.

-From the aspect of solver, the Arbitrary-Lagrangian-Eulerian (ALE) framework has been chosen for its ability to describe physical equations on a mesh moving with an imposed or an abitrary movement.

-Finally, only two-dimensional or three dimensional mono-fluid problems modelled by the inviscid compressible Euler equations are considered in this thesis.

Main contributions

During this thesis, I extend goal-oriented mesh adaptation to time-dependent simulations on moving computational domains.

My work follows the previous works of [Olivier 2011a, Belme 2011, Barral 2015].

• I contributed to the development and implementation of implicit schemes in ALE formulation for unsteady Euler equation involving moving mesh problems and compared it to the explicit time integration simulations.

• I updated the unsteady feature-based ALE metric to handle mesh optimizations during the computational loop with the error analysis from [N. [START_REF][END_REF].

• I extended the adjoint solver of unsteady Euler flows for ALE formulation in two-dimensional (validated) in three-dimensional (partially implemented)

• I added the backward moving mesh to the code.

• I added a new error estimation for goal-oriented mesh adaptation for moving mesh simulations

Organisation and content of thesis

The present thesis is built around two main topics :

• The understanding and improvements of the fluid solver in the ALE framework and the hessian-based mesh adaptation for moving mesh geometries.

• A goal-oriented based anisotropic mesh adaptation for unsteady problems involving moving mesh geometries.

The first three chapters of this work refer to the first point. The remaining parts concentrates on the ALE adjoint solver and the set up of the goal-oriented mesh adaptation.

Chapter 1 The theorical and numerical prerequisites of the next chapters.

Chapter 2 The definition, the implementation and the validation of implicit time integration scheme for the ALE flow solver.

Chapter 3 The completion with respect to several steps for the unsteady mesh adaptation for moving geometries which complement thesis [Barral 2015]. In particular, the importance to update the ALE metric in the flow solver to govern the mesh optimizations during the moving mesh algorithm.

Chapter 4 The development of the unsteady adjoint solver in the ALE framework which requires a backward in time moving mesh algorithm consistent with the forward in time one.

Chapter 5 The time accurate goal-oriented mesh adaptation for moving geometries and its associated error estimate.

Scientific communications

Papers 2015 E. Gauci, F. Alauzet, A. Loseille, A. Dervieux, Towards goal-oriented mesh adaptation for fluid-structure interaction, Coupled Problems in Science and Engineering VI 1197 : 1208, Venice, Italy [E. Gauci 2015] 2017 E. Gauci, F. Alauzet and A. Dervieux, Vortical Flow Prediction of a Moving Aircraft Using Goal-Oriented Anisotropic Mesh Adaptation, Proc. of 23rd AIAA Computational Fluid Dynamics Conference, AIAA-2017-3298, Denver, CO, USA, June 2017Denver, Co, USA Chapter 1

State of the art of previous works

The proposed works in this thesis have been built on three existing theories referred in this chapter. It is quite dense because it corresponds to four thesis : the theoretical foundation have been established in [Loseille 2008] and the three mesh adaptation theories on which this thesis relies have been introduced in [Olivier 2011a] for unsteady feature-based mesh adaptation, in [Belme 2011] for unsteady goal-oriented mesh adaptation and in [Barral 2015] for time-accurate mesh adaptation involving moving geometries.

Continuous mesh framework

Generating an adapted mesh means optimizing the accuracy of some parts of the simulation domain while de-emphazing other parts we don't need precision. To give rise to this type of mesh we must establish a generator for this assignment. And this assignment must be quasi unique depending on the current localization. At first thought and as it comes to geometry, let's look at the available geometrical tools. Considering a triangle of side 1 in an Euclidian affine space of R 2 , its area (i.e. the dot product of one of its pair of vectors) or the angle between these are the same whatever the considered triangle of the space. But generating an unstructured, anisotropic mesh means assigning privileged sizes and local orientations for the elements (triangles or tetrahedra) at each point of the domain. We will see that the use of Riemannian metric spaces is an elegant and efficient way to achieve this goal. This section is organized as follows. We first succinctly recall the metric-based mesh generation in Sections 1.1.1, 1.1.2, then the Section 1.1.3 gives the theoretical framework of anisotropic mesh adaptation based on the continuous mesh concept.

Duality between discrete and continuous element : notion of unit element

In differential geometry, a metric tensor takes as input a pair of tangent vectors at a point of a surface (or any differentiable manifold) and produces a real scalar number in a way that generalizes many of the familiar properties of the dot product of vectors in Euclidian space.

This continuous point of view enables a practical mathematical representation, and the use of powerful tools of differential geometry.

In the same way as a dot product, metric tensors are used to define the length of and angle between tangent vectors. Let's explain the duality between triangle/tetrahedron and metric tensor, in other terms between discrete and continuous element.

Euclidian space

Let us begin by studying this duality in the Euclidian space.

We consider the vector space R n , typically n = 2 or 3 in our case. A scalar product is a Symmetric Positive Definite (SPD) form. This form can be represented by a SPD matrix M = (m ij ) 1≤i,j≤n , which is called a metric tensor or simply metric. The scalar product is then written:

( • , • ) M : R n × R n -→ R + (u, v) -→ (u, v) M = u T M v = n j=1 n i=1 m ij u i v j .
(1.1)

In the simple case where M = I with I the identity matrix, the scalar product is the canonical Euclidian dot product. A real vector space with the scalar product is called an Euclidian space.

The scalar product is useful to compute lengths and angles in the Euclidian space. We can define:

• a distance:

d M : R n × R n -→ R + (P, Q) -→ d M (P, Q) = --→ P Q T M --→ P Q , (1.2)
which induces a metric space structure on the vector space,

• and a norm:

|| • || M : R n -→ R + u -→ ||u|| M = √ u T M u .
(1.3)

From these distance and norm, we deduce the classic geometrical quantities:

• the length of an edge e is given by:

M (e) = √ e T M e , (1.4)

• the angle between two vectors u 1 and u 2 is the unique real number θ ∈ [0, π] such that:

cos θ = (u 1 , u 2 ) M ||u 1 || M ||u 2 || M
(1.5)

• the volume of element K is:

|K| M = √ det M |K| where |K| = |K| I .
(1.6)

A very useful result on metrics is that M is diagonalizable in an orthonormal basis:

M = R Λ R T ,
where

               Λ = diag (λ 1 , . . . , λ n ) is the diagonal matrix made of the eigenvalues of M , R = ( v 1 | v 2 | . . . |v n ) T is the unitary matrix (i.e. R T R = I)
made of the eigenvectors of M .

(1.7)

We will very often use the geometric representation of a metric tensor. In the vicinity V(P ) of point P , the set of points that are at distance ε of P , is given by:

Φ M (ε) = {Q ∈ V(P ) | d M (Q, P ) = ε} .
(1.8)

This relation defines an ellipsoid of center a and whose axes are aligned with the eigen directions v i of M. The set Φ M (1) is called the unit ball of M, and the sizes of its axes are h i = λ -1 2 i . This unit ball is depicted in Figure 1.1 in two and three dimensions. A metric tensor M provides another useful information: the application that maps the unit ball associated with I to the unit ball associated with M. The application Λ 1 2 R where Λ 1 2 = diag(λ 1 2 i ) defines the mapping from the physical space (R n , I) to the Euclidean metric space (R n , M):

h 2 = 1 √ λ 2 r 2 r 1 r 3 h 1 = 1 √ λ 1 h2 = 1 √ λ2 h 3 = 1 √ λ 3 h 1 = 1 √ λ 1
Λ 1 2 R : (R 3 , I) -→ (R 3 , M) x -→ Λ 1 2 R x .
This natural mapping corresponds to the change of basis from the canonical basis to the orthonormal diagonalization basis of M and is depicted in Figure 1.2. The key notion of the discrete-continuous duality is the notion of unit element.

e 2 e 1 h 1 h 2 R 2 , M R 2 , I 2 t R ⇤ 1 2 ⇤ 1 2 R 1 1 v 1 v 2
A tetrahedron K, defined by its list of edges (e i ) i=1..6 , is said to be an unit element with respect to a metric tensor M if the length of all its edges is unit in metric M: ∀i = 1, ..., 6, M (e i ) = 1 with M (e i ) = e i T M e i .

(1.9)

If K is composed only of unit length edges then its volume |K| M in M is constant equal to:

|K| M = √ 2 12 and |K| = √ 2 12 (det(M)) -1 2 , (1.10)
where |K| is its Euclidean volume. It is actually used to define classes of equivalence of discrete elements: let M be a metric tensor, there exists a non-empty infinite set of unit elements with respect to M. Conversely, given an element K such that |K| = 0, then there is a unique metric tensor M for which element K is unit with respect to M.

Consequently, a discrete element can be viewed as a discrete representative of an equivalence class formed by all the unit elements with respect to a metric M. Figure 1.3 depicts some unit elements with respect to a metric tensor, which is geometrically represented by its unit-ball.

M denotes the class of equivalence of all the elements which are unit with respect to M and is called continuous element.

All the discrete representatives of a continuous element M share some common properties, called invariants, which justify the use of this equivalence relation. They connect the geometric properties of the discrete elements to the algebraic properties of metric M. The two main invariants are:

• the edges e i of any unit element K with respect to metric M are unit in M: ∀i = 1, ..., 6, e T i M e i = 1, (1.11)

• conservation of the Euclidean volume for any unit element K with respect to metric

M |K| = √ 3 4 det(M -1 2 ) in 2D and |K| = √ 2 12 det(M -1 2 ) in 3D.
(1.12)

Many other invariants are given in [Loseille 2008].

Riemannian metric space

We've just seen the link between a metric tensor and an unit element. This idea is the starting point to understand how to create a metric-based adapted mesh. This theory, initially introduced in [George 1991], is based on an unit mesh generation in a prescribed Riemannian metric space.

In fact, working in an Euclidian space is no longer sufficient. This is because the scalar product of Euclidian spaces is the same on the whole space, which means that the distance definition is the same for each point of the space (the unit ball is the same everywhere). However, when used to generate adapted meshes with different element sizes in the domain, it is convenient to have a distance that depends on the position of the space. Thus, we now consider a set of SPD tensors M = (M(P )) P ∈Ω , also called metric tensor field, defined on the whole domain Ω ⊂ R n . Locally at point P , M(P ) induces a scalar product on R n × R n . The vector space, with this new structure, is called a Riemannian metric space. In this thesis, we will use the same notation M to speak of the metric field and of the metric tensor at a given point.

Notation M will only be used if the distinction is necessary for pedagogical purposes.

Remark 1. Unlike usual Riemannian spaces, there is no notion of manifold in Riemannian metric spaces. However, Riemannian metric spaces can be assimilated to functions representing Cartesian surfaces, and the metric tensor defined for a point of the space is a scalar product on the tangent plane for that point. There is no global notion of scalar product in a Riemannian metric space, thus no global distance or norm. However, we can extend the notions of length, angle and volume from the Euclidian space case. To do so, we have to consider the variation of the metric in space. The following geometrical quantities are defined:

• the length of edge e = --→ P Q parametrized by γ : t ∈ [0, 1] -→ P + t --→ P Q is computed with the following formula:

M (e) = 1 0 ||γ (t)|| M(γ(t)) dt = 1 0 --→ P Q T M(P + t --→ P Q) --→ P Q dt , (1.13)
• the angle between two vectors u 1 = --→ P Q 1 and u 2 = --→ P Q 2 is the unique real θ ∈ [0, π] such that:

cos θ = (u 1 , u 2 ) M(P )
||u 1 || M(P ) ||u 2 || M(P ) , (1.14)

• the volume of element K with respect to M is more difficult to apprehend. Indeed, due to metric variations, element K, as seen with respect to metric field M is generally curved: it is not a simplex anymore and normally, its volume should be computed with an integration formula:

|K| M = K det M(x) dx .
(1.15)

However, this volume can be approximated at first order:

|K| M ≈ |K| I det M(G K )
, where G K is the barycenter of K.

(1.16)

Duality between discrete and continuous mesh : notion of unit mesh

This local relation of equivalence concerns elements, and needs to be somehow extended to whole meshes. Intuitively, the notion of Riemann metric space is going to play that role. The main difficulty is to take into account the variation of the function x → M(x). The analysis can be simplified if M is rewritten as follows, separating its local and global properties. A

Riemannian metric space M = (M(x)) x∈Ω is locally written :

∀x ∈ Ω, M(x) = d 2 3 (x) R(x)      r -2 3 1 (x) r -2 3 2 (x) r -2 3 3 (x)      R T (x), (1.17)
where

• density d is equal to: d = (λ 1 λ 2 λ 3 ) 1 2 = (h 1 h 2 h 3 ) -1
, with λ i the eigenvalues of M

• anisotropic quotients r i are equal to:

r i = h 3 i (h 1 h 2 h 3 ) -1
• R is the eigenvectors matrix of M representing the orientation.

The density d controls only the local level of accuracy of M: (increasing or decreasing it does not change the anisotropic properties or the orientation), while the anisotropy is given by the anisotropic quotients and the orientation by matrix R. The notion of complexity C of M can also be defined:

C(M) = Ω d(x) dx = Ω det(M(x)) dx. (1.18)
This quantifies the global level of accuracy of (M(x)) x∈Ω .

From this formulation of a Riemannian metric field arises a duality between meshes and Riemannian metric spaces. This duality is justified by the strict analogy between the following discrete and continuous notions: orientation vs. R, stretching vs. r i , size vs. d and number of elements vs. C(M). However, the class of discrete meshes represented by M is complex to describe.

If the notion of unit element in the previous section was quite immediate, things are more complicated when it comes to define an unit mesh. Stricto sensu, an unit mesh is a mesh whose all edges are unit with respect to the prescribed metric field. However, the existence of a mesh made of unit elements is not assured. For example, in R 3 , let's take the simplest case of M(P ) = I(P ) for each point P , i.e. the canonical Euclidian space. It is well known that R 3 cannot be filled with regular tetrahedra (that are unit with respect to the identity metric).

So the constraint on the sizes of the edges has to be relaxed. But this can lead to meshes with very bad quality elements (flat elements). So we have to add a constraint on the volume of the elements, through a quality function.

For this reason, we have to introduce the notion of quasi unit element. A tetrahedron is said to be quasi unit with respect to a metric field M if its edges are close to unit, i.e.

∀i, M (e i ) ∈ [ 1 √ 2 , √ 2].
To avoid elements with a null volume (see [Loseille 2011a]), we have to add a constraint on the volume, which is achieved through a quality function:

Q M (K) = √ 3 216 6 i=1 2 M (e i ) 3 2 |K| M ∈ [1, +∞] . (1.19)
For the regular tetrahedron, the quality function is equal to 1, whereas it tends to +∞ for a null volume tetrahedron. So an element close to a perfectly unit element has a quality close to 1. This leads to the following definition. A tetrahedron K defined by its list of edges (e i ) i=1...6

is said quasi-unit for Riemannian metric space (M(x)) x∈Ω if ∀i ∈ [1, 6], M (e i ) ∈ 1 √ 2 , √ 2 and Q M (K) ∈ [1, α] with α > 1 . (1.20)
The definition of unit mesh consequently becomes : a unit mesh with respect to a Riemannian metric space (M(x)) x∈Ω is a mesh made of quasi-unit elements.

In practice, it is this definition that is used in meshing software. For any kind of desired mesh (uniform, adapted isotropic, adapted anisotropic), the mesh generator will generate a mesh that is unit with respect to the prescribed Riemannian metric space. Thus the resulting mesh is uniform and isotropic in the Riemannian metric space while it is adapted and anisotropic in the canonical Euclidian space. This is illustrated in Figure 1.5. This idea has turned out to be a huge breakthrough in the generation of anisotropic adapted meshes.

Duality between discrete and continuous interpolation error: notion of continuous linear interpolate

The model that we have just defined is used to obtain an analytic expression of the optimal mesh.

The purpose now in this section is to obtain a still valid equivalence between a continuous error estimate for any function on any continuous mesh and a discrete error estimate on a discrete mesh. In mathematical terms, let (M(x)) x∈Ω be a continuous mesh of a domain Ω and let u be a smooth representation This function u is a non-linear function which is assumed to be twice continuously differentiable.

Inputs (H 0 , M i ) i∈H
We define the following approximation space from the standard P 1 FEM approximation space. Let H 1 (Ω) be the Sobolev space defined as the set of applications of L 2 (Ω) such that their first weak derivative exists and also belongs to L 2 (Ω).

H 1 (Ω) = {u ∈ L 2 (Ω) | Du ∈ L 2 (Ω)}.
Let us introduce the following approximation space :

V h = {ϕ ∈ H 1 (Ω) | ϕ h|K is affine ∀K ∈ H}.
The usual P 1 -projector Π h is defined such as Π h u is exact on each vertex of the element K :

Π h : H 1 (Ω) → V h | Π h u(x) = ϕ(x), ∀x vertex of H
In this section, we want to seek a well-posed definition of the continuous linear interpolation error uπ M u L 1 (Ω) related to a continuous mesh (M(x)) x∈Ω which also implies a well-posed definition of a linear continuous interpolate π M u. And obviously, we would like the continuous linear interpolation error to be a reliable mathematical model of u-Π h u L 1 (Ω h ) .

In this sense, the interpolation error is first derived locally (on an element) for a quadratic function, then is extended to a global definition (for any point of the domain).

Local continuous interpolate

Let us first consider a quadratic function u defined on a domain Ω ⊂ R n , and a continuous element M. For all unit discrete elements K with respect to M, the interpolation error of u in L 1 norm does not depend on the element shape and is only a function of the Hessian H u of u and of continuous element M. For more details see [Loseille 2011a].

• In 3D, for all unit elements K for M, the following equality holds:

u -Π h u L 1 (K) = √ 2 240 det(M -1 2 ) trace(M -1 2 H u M -1 2 ).
(1.21)

• In 2D, for all unit elements K for M, the following equality holds:

u -Π h u L 1 (K) = √ 3 64 det(M -1 2 ) trace(M -1 2 H u M -1 2 ).
For all the discrete elements that are unit with respect to M, the interpolation error is the same, and is only based on continuous quantities (metric and Hessian). This last remark is important, since it shows metrics contain all the information required to compute the interpolation error, and are thus well adapted for anisotropic control of this error.

Global continuous interpolate

To define a global continuous linear interpolate, the problem is once again how to move from an expression valid for one continuous element to an expression for the case in which the metric varies point-wise. Let us now suppose that the continuous mesh (M(x)) x∈Ω is varying and that the function u is no more quadratic but only twice continuously differentiable. Equality (1.21) does not hold anymore, but all the terms of the right-hand side M and H are still well defined continuously.

Let's consider a point a ∈ Ω. In the vicinity of a, we denote u Q the quadratic approximation of smooth function u as the truncated second order Taylor expansion of u and (M(x)) x∈Ω reduces to M(a) in the tangent space. There exists a unique function π M such that:

∀a ∈ Ω , |u-π M u|(a) = u Q -Π h u Q L 1 (K) |K| = 1 20 trace M(a) -1 2 |H(a)| M(a) -1 2 , (1.22)
for every K unit element with respect to M(a).

This result underlines another discrete-continuous duality by pointing out a continuous counterpart of the interpolation error. For this reason, the following formalism was adopted:

π M is called continuous linear interpolate and |uπ M u| represents the continuous dual of the interpolation error illustrated in cases of a concave or a convex function on Figure 1.6. From a practical point of view, we deduce the following analogy. Given an unit mesh H of a domain Ω h with respect to a continuous mesh (M(x)) x∈Ω , the global interpolation error is:

u -Π h u L 1 (Ω h ) = K∈H u -Π h u L 1 (K) .
(1.23)

In the continuous case, the discrete summation becomes an integral:

u -π M u L 1 (Ω) = Ω |u -π M u|(x) dx.
(1.24) There is no global guarantee on the reliability of the continuous interpolation error given by Relation (1.24), and in particular, there is no a priori relationship between (1.23) and (1.24).

u concave u convex
The only guarantee is the local equivalence given by Equation (1.22). However, the local guarantee becomes global when the mesh is unit with respect to a constant metric tensor and when the function is quadratic. In the latter case, by neglecting error due to the boundary discretization, we have the equality:

u -Π h u L 1 (Ω h ) = u -π M u L 1 (Ω) , (1.25)
for all unit meshes H with respect to (M(x)) x∈Ω .

Several examples are given in [Loseille 2011b], both analytic and numerical, that confirm the validity of this analogy. They show that the model is accurate and the equivalence (1.23)≈ (1.24) is observed even for non quadratic functions and non-constant continuous meshes, and that the error due to the fact that the mesh generator generates edges with length not strictly equal to one is negligible. In particular, the range for the lengths of the edges given in Section 1.1.2 ensures reliable numerical results.

Summary

We have presented a framework that draws a correspondence between the discrete domain and the continuous domain. This framework is summarized in Table 1.1.

Mesh adaptation 1.2.1 Steady mesh adaptation

Studying steady simulations is a necessary prerequisite for a better understanding of unsteady simulations. Let us first detail the feature-based mesh adaptation algorithm and the attached error analysis for steady problems then we turn to the description of goal-oriented mesh adaptation in the steady case.

Feature-based mesh adaptation for steady flows

The problem can be simplified to the specification of a mesh that is optimal for the interpolation error.

Mesh adaptation consists in finding the mesh that minimizes a certain error on a domain Ω ∈ R n , for a certain sensor function u. The error we consider is the interpolation error, that we control in L p norm -depending on the choice of p, different aspects of the solution are captured, as will be seen later. The problem is stated a priori :

Find H L p having N vertices such that E L p (H L p ) = min H u -Π h u L p (Ω h ) .
(P ) (P ) is too complex to be solved directly, since the unknown, i.e. the mesh, is made up of vertices and their topology, which is far too many degrees of freedom. Besides, several optimal meshes can be found for one sensor function (think about merely swapping an element), so the problem is ill-posed. On the other hand, it is possible to show that the problem moved to the continuous domain is well posed, and can be solved using calculus of variations. What is more, the continuous formulation of the adaptation problem uses a global point of view, while usual methods focus on a local analysis of the error. The reformulated problem is:

DISCRETE CONTINUOUS Element K Continuous element = Metric tensor M Element volume |K| d -1 = det(M -1 ) Mesh H of Ω h Riemannian metric space M(x) = (M(x)) x∈Ω Number of vertices N v Complexity C(M) = Ω d(x) dx P 1 interpolate Π h P 1 -continuous interpolate π M Local element-wise interpolation error Local point-wise interpolation error e h (K) = ||u -Π h u|| K e(x) = (u -π M u) (x) e h (K) = |K| 40 6 i=1 e T i |H u |e i (for u quadratic) e(x) = 1 20 trace M(x) -1 2 |H(x)| M(x) -1 2 Global continuous interpolation error K∈H ||u -Π h u|| K Global interpolation error x∈Ω e(x) dx
Table 1.1: The continuous mesh model draws a correspondence between the discrete domain and the continuous domain.

Find M L p having a complexity of N such that E L p (M L p ) = min M u -π M u L p (Ω) . (1.26)
Using the definition of the linear continuous interpolate π M given by Equation (1.22), the well-posed global optimization problem of finding the optimal continuous mesh minimizing the continuous interpolation error in L p norm can be established:

Find M L p such that E L p (M L p ) = min M Ω u(x) -π M u(x) p dx 1 p (1.27) = min M Ω trace M(x) -1 2 |H u (x)|M(x) -1 2 p dx 1 p , under the constraint C(M) = Ω d(x) dx = N .
The constraint on the complexity notably avoids the trivial solution where all (h i ) i=1,n are zero which provides a null error.

Contrary to a discrete analysis, this problem can be solved globally by using a calculus of variations that is well-defined on the space of continuous meshes. In [Loseille 2011b], it is proved that Problem (1.27) admits a unique solution M L p = (M L p (x)) x∈Ω which is locally defined by:

M L p (x) = N 2 n Ω det(|H u (x)|) p 2p+n dx -2 n det(|H u (x)|) -1 2p+n |H u (x)| . (1.28)
Moreover, it verifies the following properties:

• M L p is unique
• M L p is locally aligned with the eigenvectors basis of H u and has the same anisotropic quotients as H u

• M L p provides an optimal explicit bound of the interpolation error in L p norm:

u -π M L p u L p (Ω) = n N -2 n Ω det (|H u |) p 2p+n 2p+n np
.

(1.29)

• For a sequence of continuous meshes having an increasing complexity with the same orientation and anisotropic quotients (M N L p ) N =1...∞ , the asymptotic order of convergence verifies:

u -π M N L p u L p (Ω) ≤ Cst N 2/n . (1.30)
Relation (1.30) points out a global second order of mesh convergence.

Anisotropic mesh adaptation is a non-linear problem, therefore an iterative procedure is required to solve this problem. For stationary simulations, an adaptive computation is carried out via a mesh adaptation loop inside which an algorithmic convergence of the mesh-solution couple is sought. This mesh adaptation loop is schematized in Algorithm 1 where H, S and M denote respectively meshes, solutions and metric fields.

Algorithm 1 Feature-based mesh adaptation loop for steady flows Input : Initial mesh and solution (H 0 , S 0 0 ) and set targeted complexity N For i = 0, n adap 1. (S i )= Compute solution with the flow solver from pair (H i , S 0 i );

If i = n adap break; 2. (M L p ,i ) = Compute metric (M L p ) according to selected error estimate from (H i , S i ); 3. (H i+1 ) = Generate a new adapted mesh form pair (H i , M L p ,i ); 4. (S 0 i+1 ) = Interpolate new initial solution from (H i+1 , H i , S i );

End For

Before detailing the different steps of the numerical algorithm, let's focus on two major difficulties of it:

• The solution u of the problem is not known. We only can access to the numerical approximation u h .

• As for all numerical simulations, a control of the approximation error uu h is expected.

We demonstrate how this problem can be simplified to the specification of a mesh that is optimal for the interpolation error. We describe how the interpolation theory is applied when only u h , a piecewise linear approximation of the solution, is known. Indeed, in this particular case, the interpolation error estimate is not applied directly to u nor to u h .

Let V 1 h be the space of continuous piecewise linear function associated with a given mesh H of domain Ω h . We denote by R h a reconstruction operator applied to numerical approximation u h . This reconstruction operator can be either a recovery process [Zienkiewicz 1992], a hierarchical basis [Bank 1993], or an operator connected to an a posteriori estimate [Huang 2010a]. We assume that the reconstruction R h u h is better than u h for a given norm . in the sense that:

u -R h u h ≤ α u -u h where 0 ≤ α < 1 .
From the triangle inequality, we deduce:

u -u h ≤ 1 1 -α R h u h -u h .

If reconstruction operator R h has the property: Π

h R h φ h = φ h , ∀φ h ∈ V 1 h
, the approximation error of the solution can be bounded by the interpolation error of reconstructed function R h u h :

u -u h ≤ 1 1 -α R h u h -Π h R h u h .
From previous section, if H L p is an optimal mesh to control the interpolation error in L p norm of R h u h , then the following upper bound of the approximation error can be exhibited:

u -u h L p (Ω h ) ≤ n N -2 n 1 -α Ω det (|H R h u h (x)|) p 2p+n dx 2p+n np .
Remark 2. It is important to note that M L p defined by Relation (1.35) applied to R h u h does not allow us to generate an optimal adapted mesh to control the approximation error uu h .

If all assumptions are verified then we have an upper bound which allows that such generated adapted meshes controls the approximation error.

Now, we details each step of the mesh adaptation algorithm.

Step 1 : The flow solver Wolf . All the numerical simulations of this thesis are concerned with the compressible Euler equations. Results on Navier-Stokes equation are not discussed in this thesis. Assuming that the gas is perfect, inviscid and that there is no thermal diffusion, the 3D compressible Euler equations for mass, momentum and energy conservation read:

∂W ∂t + ∇ • F(W ) = 0, (1.31)
where W is the vector of conservative variables and F is the convection operator

F(W ) = (F 1 (W ), F 2 (W ), F 3 (W )) with: W =          ρ ρu ρv ρw ρE          , F 1 (W ) =          ρu ρu 2 + p ρuv ρuw (ρE + p)u          , F 2 (W ) =          ρv ρuv ρv 2 + p ρvw (ρE + p)v          , F 3 (W ) =          ρw ρuw ρvw ρw 2 + p (ρE + p)w          .
We have denoted by ρ the density, u = (u, v, w) the Cartesian velocity vector, E = T + u 2 2 the total energy and p = (γ -1)ρT the pressure with γ = 1.4 the ratio of specific heats and T the temperature.

To discretize convective terms Wolf uses the mixed-element-volume (MEV) approach [Alauzet 2010a] initiated by [START_REF] Stoufflet | Numerical simulation of 3-D hypersonic Euler flows around space vehicles using adapted finite element[END_REF], Fezoui 1989, Debiez 2000, Cournède 2006]. It is a vertex-centered finite volume scheme applied to tetrahedral unstructured meshes. This scheme uses a particular edge-based formulation with upwind elements.

For the following, it is important to specify the considered semi-discretization. The system of partial differential equations can be written in semi-discrete form for a vertex P i as:

|C i | dW i dt = R(W i ) ,
where R i is the residual, i.e., the discretization of the convective operator using an approximate Riemann solver and of the boundary conditions, and |C i | is the area/volume of the dual finite volume cell associated with P i .

For an explicit time discretization, the semi-discretized system reads:

|C i | δt n i W n+1 i -W n i = R(W n i )
where W n is the state at iteration n and δt n i is the local time step at iteration n. For implicit time discretization, we have :

|C i | δt n i W n+1 i -W n i = R(W n+1 i )
which is linearized as:

|C i | δt n i I d - ∂R ∂W (W n i ) W n+1 i -W n i = R(W n i )
where ∂R ∂W (W n i ) contributes the i th line of the matrix. We then rewrite the linearized system in compact form for all vertices in the mesh as:

A n δW n = R n where A n = |C| δt n I - ∂R n ∂W and δW n = W n+1 -W n .
Step 2 : The metric module Metrix computes the continuous mesh and performs the metric field gradation [Alauzet 2010b]. The optimal continuous mesh M L p (u) minimizing Problem (1.27) is given in (1.35). The metric gradation field is used to smooth and control the variation of the obtained optimal continuous mesh.

Step 3 : The local adaptive remesher Feflo.a generates an unit mesh with respect to a prescribed metric field. Feflo.a is a classical metric-based local remesher [START_REF] Loseille | [END_REF]] based on a unique cavity operator which is a generalization of standard operators (insertion, collapse, swap of edges and faces, vertex smoothing) [Loseille 2013]. The cavity operator also performs combinations of standard operators in a very natural way. To this end, modifications on the cavity are done to either favor a modification, that would have been rejected with the standard operator, or to improve the final quality by combining automatically many standard operators at once. One important capability is to adapt the volume and the surface mesh in a coupled way so that a valid 3D mesh is always guaranteed on output.

Step 4 : The software component Interpol to project linearly the solution defined on the previous mesh onto the new mesh. The first step is to localize the new vertices in the element of the previous mesh. The second step consists in interpolating linearly the solution using the barycentric coordinates. More details on the interpolation stage can be obtained in [Alauzet 2010b, Alauzet 2016].

Goal-oriented mesh adaptation for steady flows

In the previous paragraph, metric-based anisotropic mesh adaptation method is limited to the minimization of some interpolation errors for some solution fields called sensors and did not take into account the PDE of the considered problem. If, for many applications, this simplifying standpoint is an advantage, there are also many applications where Hessian-based adaptation is far from optimal regarding the way the degrees of freedom are distributed in the computational domain. Indeed, Hessian-based methods aim at controlling the interpolation error but this purpose is not often so close to the objective that consists in obtaining the best solution of a PDE. This is particularly true in many engineering applications where a specific function needs to be accurately evaluated : lift, drag, heat flux, pressure field, etc. Hessian-based adaptation has not been especially designed to address this issue.

The objective in goal-oriented adaptation is to minimize the error on the considered functional j, details can be found in [Loseille 2010b] in the case of compressible Euler equations.

For simplicity's sake boundary terms are not mentioned in the following. The approximation Algorithm 2 Goal-oriented mesh adaptation loop for steady flows Input : Initial mesh and solution (H 0 , S 0 0 ) and set targeted complexity End For error on the functional j(W ) is given by the following error estimate:

N For i = 1, n adap 1. (S i )= Compute
|j(W ) -j(W h )| ≤ Ω h |∇W * | |F(W ) -Π h F(W )|dΩ h (1.32)
where W * is the adjoint state. We observe that this estimate is expressed in term of the interpolation error in L 1 norm of the Euler fluxes applied to the continuous solution W weighted by the gradient of the continuous adjoint state W * .

The goal-oriented mesh adaptation loop is given by Algorithm 2. There is one more step to compute the adjoint state. Let us detail the steps 2 and 3 as all other steps listed earlier stay the same :

Step 2 : The adjoint flow solver Wolf . Let j(W ) be the considered output functional.

Using the same notations as in the previous section, the adjoint state W * is solution of the following linear system:

A * W * = ∂j ∂W (W)
where A * is the transpose of the jacobian matrix:

A * = -∂R n ∂W T .
Step 3 : The metric module Metrix computes the following continuous goal-oriented optimal metric and performs the metric field gradation. Using the same reasoning as for the error estimator for a sensor but applied to the goal-oriented error estimate (1.32), the error estimate is rewritten in a continuous form: (1.33) where M = (M(x)) x∈Ω is a continuous mesh and π M is the continuous linear interpolate. We are now focusing on the following (continuous) mesh optimization problem:

|j(W ) -j(W h )| ≈ E go (M) = Ω |∇W * | • |F(W ) -π M F(W )| dΩ ,
Find M go such that E go (M go ) = min M Ω |∇W * | • |F(W ) -π M F(W )| dΩ = min M Ω trace M -1 2 (x) |H go (x)| M -1 2 (x) dx
under the constraint C(M) = N and where we have introduced the goal-oriented Hessian-metric given in 3D by

|Hgo(x)| = 5 j=1 ∂W * j ∂x (x) • H(F1(Wj))(x) + ∂W * j ∂y (x) • H(F2(Wj))(x) + ∂W * j ∂z (x) • H(F3(Wj))(x) .
(1.34)

The formulation of the optimal goal-oriented continuous mesh M go = (M go (x)) x∈Ω is locally

given by:

M go (x) = N 2 n Ω det(|H go (x)|) 1 2+n dx -2 n det(|H go (x)|) -1 2+n |H go (x)| . (1.35)
Note that we have the same expression as for the feature-based case in L 1 norm where the Hessian of the sensor has been replaced by the goal-oriented Hessian-metric.

Remark 3.

-To minimize the function j we establish a continuous optimisation. But it's today impossible to solve this problem in a continuous way. To overcome this problem, we choose to discretize it. The nodal gradients of the state function and of the adjoint function are recovered from the gradients to the elements using a L 2 local projection and Clément's interpolation operator [Clément 1975].

-The details of the computation for goal-oriented mesh adaptation will be detailed in chapter 5.

-In this document, the metric optimization is defined as continuous and is always discretized to be solved.

Unsteady mesh adaptation

Feature-based mesh adaptation for unsteady flows

In the context of time-dependent problems, the error analysis must control spatial and temporal error. In [Olivier 2011a, Alauzet 2016], time discretization errors are not taken into account but the focus was made on a space-time analysis of the spatial error. In other words, the work seeks for the optimal space-time mesh controlling the space-time spatial discretization error. For the type of considered simulations, the assumption is made that as an explicit time scheme is used for time advancing, then the error in time is controlled by the error in space under CFL condition. This has been demonstrated under specific conditions in [Alauzet 2007]. As long as this hypothesis holds, the spatial interpolation error provides a fair measure of the total space-time error of the discretized unsteady system. Notably in the case of implicit time advancing solvers, the analysis would have to be completed to take into account the temporal discretization error, as is done in [Coupez 2013], which defines the optimal time discretization (i.e., time steps).

Our goal is to solve an unsteady PDE which is set in the computational space-time domain

Q = Ω × [0, T ]
where T is the (positive) maximal time and Ω ⊂ R n is the spatial domain. Let Π h be the usual P 1 projector, we extend it to time-dependent functions:

(Π h ϕ) (t) = Π h (ϕ(t)) , ∀ t ∈ [0, T ] . (1.36)
The considered problem of mesh adaptation consists in finding the space-time mesh H of Ω h × [0, T ] that minimizes the space-time linear interpolation error u-Π h u in L p norm. The problem is thus stated in an a priori way:

Find H L p having N st vertices such that E L p (H L p ) = min H u -Π h u L p (Ω h ×[0,T ]) . (1.37)
This problem is ill-posed and has far too many unknowns to be solved directly, as was explained previously. So it is rewritten in the continuous mesh framework under its continuous form:

Find M L p = (M L p (x, t)) (x,t)∈Q such that E L p (M L p ) = min M u -π M u L p (Ω×[0,T ]) , (1.38)
under the space-time constraint:

C st (M) = T 0 τ (t) -1 Ω d M (x, t) dx dt = N st , (1.39)
where τ (t) is the time step used at time t of interval [0, T ]. Introducing the continuous interpolation error, we recall that we can write the continuous error model as follows:

E L p (M) = T 0 Ω trace M -1 2 (x, t)|H u (x, t)|M -1 2 (x, t) p dx dt 1 p , (1.40)
where H u is the Hessian of sensor u. To find the optimal space-time continuous mesh, Problem (3.3-3.4) is solved in two steps: first, a spatial minimization is done for a fixed t, then a temporal minimization is performed where the time step τ is specified by the user as a function of time t → τ (t). It leads to the expression of the optimal space-time metric M L p for a prescribed time step τ (t) [Alauzet 2016]:

M L p (x, t) = N 2 n st T 0 τ (t) -2p 2p+n K(t)dt -2 n τ (t) 2 2p+n (det |H u (x, t)|) -1 2p+n |H u (x, t)| , (1.41)
where

K(t) = Ω (det |H u (x, t)|) p 2p+n dx .
The computation of the optimal instantaneous continuous mesh given by Relation (3.12) involves a global normalization term which requires the knowledge of quantities over the whole simulation time frame. Thus, the complete simulation must be performed before evaluating any space-time continuous mesh. To solve this issue, we suggest to consider a global fixed-point mesh adaptation algorithm covering the whole time frame [0, T ]. This iterative algorithm is used to converge the non-linear mesh adaptation problem, i.e., converging the mesh-solution couple. This is also a way to predict the solution evolution and to adapt the mesh accordingly.

Moreover, the previous analysis provides the optimal size of the adapted meshes for each time level. Hence, this analysis requires the mesh to be adapted at each flow solver time step which is inconceivable in practical applications. We propose to use a coarse adapted discretization of the time axis. The basic idea consists in splitting the simulation time frame [0, T ] into n adap adaptation sub-intervals:

[0, T ] = [0 = t 0 , t 1 ] ∪ . . . ∪ [t i , t i+1 ] ∪ . . . ∪ [t n adap -1 , t n adap = T ] ,
and to keep the same adapted spatial mesh for each time sub-interval. On each sub-interval, the mesh is adapted to control the solution accuracy from t i to t i+1 . Consequently, the timedependent simulation is performed with n adap different adapted meshes. This drastically reduces the number of remeshings during the simulation, hence the number of solution transfers. The unsteady mesh adaptation algorithm is presented in Algorithm 3.

We now present the modified or new steps involve for time-accurate mesh adaptation.

Step 1 (a) : The software component Interpol to project the solution defined on the previous mesh onto the new mesh between each sub-interval. This stage becomes crucial in the context of unsteady problems as error due to the interpolation step may accumulate throughout the simulation. To minimize the error and to be consistant with the considered PDE (equations of conservation), a P 1 -exact conservative interpolation is considered. More details on the conservative interpolation stage can be obtained in [Alauzet 2010c, Alauzet 2016].

Step 1 (b) : The flow solver Wolf . The main difference between the steady and the unsteady case is that a global time stepping is considered for time-dependent problems. If we consider an explicit first-order time integration scheme, then the semi-discrete unsteady model = Compute all sub-interval metrics according to error estimate from

(C j , {|H L 1 | j i } n adap i=1 ); 4. {H j+1 i } n adap i=1
= Generate all sub-interval adapted meshes from pair ({H j i , M j L p ,i } n adap i=1 );

End For at time t n is:

|C i | W n i -W n-1 i δt n = R(W n-1 i ) .
Step 3 : The metric module Metrix . We extend the previous analysis to the fixed-point mesh adaptation algorithm context -described in Algorithm 3 -where the simulation time interval [0,T] is split into n adap sub-intervals [t i-1 , t i ] for i = 1, .., n adap . Each spatial mesh M i is then kept constant during each sub-interval [t i-1 , t i ]. We could consider this partition as a time discretization of the mesh adaptation problem. Following the previous section analysis, we deduce the optimal continuous mesh M L p = {M i L p } i=1,..,n adap defined locally by:

M i L p (x) = N 2 n st n adap j=1 K j tj tj-1 τ (t) -1 dt 2p 2p+n -2 n ti ti-1 τ (t) -1 dt -2 2p+n (det H i L 1 (x)) -1 2p+n H i L 1 (x) (1.42)
and the associated optimal space-time error: .43) where

E L p (M L p ) = n N -2 n st n adap i=1 K i t i t i-1 τ (t) -1 dt 2p 2p+n 2p+n np , ( 1
K i = Ω det H i L 1 (x)
p 2p+n dx and the Hessian-metric H i L 1 on sub-interval i is the time-average of the Hessian of the sensor u:

H i L 1 (x) = t i t i-1 |H u (x, t)| dt .
(1.44)

Remark 4. The mesh adaptation with unsteady fixed-point method does not converge to the 2nd order for singularities because of the uniform time step. Multi-rate methods as in [Itam 2017] can be a solution to overcome the problem. An other issue is to modify he number of sub-intervals during the adaptive loop [Alauzet 2016].

Goal-oriented mesh adaptation for unsteady flows

We have seen the "in-house" goal-oriented mesh adaptation technique for steady problems and the feature-based mesh adaptation for unsteady flows. Now, let's turn to the goal-oriented mesh adaptation for unsteady flows.

The objective in goal-oriented adaptation for unsteady problems is to minimize the error on the functional j(W ) which is now integrated in time. An a priori error analysis [START_REF] Belme | [END_REF] gives the following error estimate for the compressible Euler equations where we have again neglected the boundary terms for the sake of simplicity :

|j(W ) -j(W h )| = T 0 Ω ∂W * ∂t |W -Π h W |dΩdt + T 0 Ω |∇W * | • |F(W ) -Π h F(W )|dΩdt (1.45)
The Algorithm 4 is further complicated by the fact that the unsteady adjoint solver must be computed backward in time (eg. step 2) after computing the solution over the simulation time frame.

Algorithm 4 Goal-Oriented Mesh Adaptation for Unsteady Flows

Input : Initial mesh and solution (H 0 , S 0 0 ) and set targeted space-time complexity N st # Fixed-point loop to converge the global space-time goal-oriented mesh adaptation problem For j = 1, n ptf x # Adaptive loop to compute forward the solution state in time on time frame [0, T ] = Generate all sub-interval adapted meshes ({H j i ,

1. For i = 1, n adap (a) S j 0,i = Interpolate conservatively next sub-interval initial solution from (H j i-1 , S j i-1 , H j i ); (b) S j i =
M j i } n adap i=1 );
EndFor Now, we describe the new steps of the unsteady goal-oriented mesh adaptation algorithm which are the computation of the unsteady adjoint state and the unsteady goal-oriented error estimate.

Step 2 : The unsteady adjoint flow solver Wolf . Let j(W ) be the considered output functional. using the same notations as previously, the unsteady adjoint state W * is solution of:

- ∂W * ∂t + A * W * = ∂j ∂W (W)
where A * is the transpose of the jacobian matrix:

A * = -∂R n ∂W T .
If we consider an explicit first-order time integration scheme, then the semi-discrete unsteady adjoint model at time t n is:

|C i | W * ,n-1 i -W * ,n i -δt n = ∂j ∂W (W n-1 i ) + W * ,n i T ∂R ∂W (W n-1 i ) .
We notice that the computation of the unsteady adjoint state at time t n-1 requires the knowledge of the adjoint state at time t n and the solution state at time t n-1 . Therefore, the solution state must be known for the whole simulation time frame [0, T ] to be able to compute the adjoint state.

Step 4 : The metric module Metrix . Following the analysis for the steady goal-oriented mesh adaptation and for the unsteady mesh adaptation, the continuous mesh adaptation problem to solve is:

Find M go such that E go (M go ) = min M T 0 Ω ∂W * ∂t |W -π M W | + |∇W * | • |F(W ) -π M F(W )| dΩdt = min M T 0 Ω trace M -1 2 (x, t) |H go (x, t)| M -1 2 (x, t) dxdt
under the space-time constraint:

C st (M) = T 0 τ (t) -1 Ω d M (x, t) dx dt = N st ,
where τ (t) is the time step used at time t of interval [0, T ]. The optimal continuous mesh is defined on the space time domain Q = Ω × [0, T ]: M go = (M go (x, t)) (x,t)∈Q and the unsteady goal-oriented Hessian-metric H go in 3D is

|H go (x, t)| = 5 j=1 ∂W * j ∂t (x, t) • H(W j )(x, t) + ∂W * j ∂x (x, t) • H(F 1 (W j ))(x, t) + ∂W * j ∂y (x, t) • H(F 2 (W j ))(x, t) + ∂W * j ∂z (x, t) • H(F 3 (W j ))(x, t) .
As for the feature-based unsteady mesh adaptation, solving the above optimization problem will provide the optimal instantaneous adapted mesh for each time level. Again, we split the simulation time frame into n adap sub-intervals and we derive an optimal spatial mesh for each sub-interval. The same space-time error analysis provides the expression of the optimal unsteady goal-oriented mesh for each sub-interval M go = M i go i=1,..,n adap defined locally by:

M i go (x) = N 2 n st n adap j=1 K j tj tj-1 τ (t) -1 dt 2p 2p+n -2 n ti ti-1 τ (t) -1 dt -2 2p+n (det H i go,L 1 (x)) -1 2p+n H i go,L 1 (x) (1.46)
where

K i = Ω det H i go,L 1 (x)
p 2p+n dx and the Hessian-metric H i go,L 1 on sub-interval i is the time-average of the unsteady goal-oriented Hessian-metric:

H i go,L 1 (x) = t i t i-1 |H go (x, t)| dt . (1.47)
As for the steady case, we notice that we have the same expression as for the feature-based case where the Hessian of the sensor has been replaced by the unsteady goal-oriented Hessian-metric.

The associated optimal space-time error is:

E go (M go ) = n N -2 n st n adap i=1 K i t i t i-1 τ (t) -1 dt 2p 2p+n 2p+n np
.

(1.48)

Mesh adaptation for moving geometries

Given that feature-based and goal-oriented anisotropic mesh adaptations have proved their efficiency to reduce the CPU time of steady and unsteady simulations while improving their accuracy it has become evident to extend the theory to time-dependent problems with moving geometries. Yet this task is far from straightforward. In [Barral 2015]

Mesh adaptation problem

The feature-based unsteady mesh adaptation cannot directly be extended to moving mesh simulations because it does not take into account the movement of the mesh. Indeed, the error analysis from which the optimal continuous mesh is deduced, assumes that the mesh is constant in time inside a sub-interval, therefore it does not take into account the local deformation of the mesh which necessarily impacts local errors. As regards the discrete representation of the metric field for moving meshes, two options are available. The first approach is an Eulerian approach:

a metric is associated with a fixed position in space, thus the metric field is a "background" field evolving in space and time but independently from the moving mesh. The second approach is a Lagrangian approach, which is considered in this work, in that case a metric is attached to a moving vertex and moves with it, i.e., we have M(x(t), t) which we will write M(x(t)) in what follows.

In the context of dynamic meshes, the optimization problem is still given by Relations (3.3) and (3.4) of the unsteady featured-based mesh adaptation problem. But it is not obvious that for dynamic meshes problem, the space-time error model is given by Relation (1.40).

To complete our reflection, the ALE framework is modelled as follow (see As for the fixed-mesh case, to perform the space-time L p error analysis with time subintervals, we first do a spatial minimization for a given sub-interval. We consider the i th subinterval [t i , t i+1 ]. Given the continuous mesh spatial complexity N i , we seek for the optimal Arbitrary-Lagrangian-Eulerian (ALE) continuous mesh M i,ALE

L p = M i,ALE L p (x(t)) x∈Ω(t)
for the whole sub-interval which is solution of the following problem:

E i,ALE L p (M i,ALE L p ) = min M i t i+1 t i Ω(t) trace (M i ) -1 2 (x(t)) |H u (x(t), t)| (M i ) -1 2 (x(t)) p dx(t) dt such that C(M i ) = N i .
Note that the continuous mesh spatial complexity is constant on this sub-interval. We name it optimal ALE continuous mesh to stress that the continuous mesh is dynamic (i.e., it evolves in time), it represents a moving adapted mesh and it is coupled with an ALE flow solver. And this is the heart of the problem : M i has a dependence in time ! Thus we cannot move the integral over time into the trace to make the time-averaged Hessian-metric appear like in the fixed mesh analysis. The spatial minimization for a given sub-interval can not be done directly.

To solve this issue while preserving the adaptation of the mesh despite the mesh deformation, we are now looking for a mesh at the beginning t i of each sub-interval [t i , t i+1 ] that will be adapted to the solution at each time-step of this sub-interval once moved with the prescribed mesh displacement. The idea which is explained in the next paragraph is to exhibit the optimal instantaneous ALE continuous mesh minimizing the interpolation error in L p norm at a given time t ∈ [t i , t i+1 ]. We shortly call it the ALE metric field because it takes into account the mesh deformation. It is then possible to map this continuous moving mesh of Ω(t), t ∈ [t i , t i+1 ], onto a continuous mesh of Ω(t i ). Thus, the error model can be written only on Ω(t i ) and the spatial minimization can be performed.

1.3.2 Optimal instantaneous ALE continuous mesh minimizing the interpolation error in L p norm: the ALE metric

The computational domain is time-dependent Ω(t) ⊂ R 3 with t ∈ [0, T ]. The problem for one time-step is the following. Let t n and t n+1 ( in red in the Figure 1.7) be two different times.

And let Ω n = Ω(t n ) and Ω n+1 = Ω(t n+1 ) be the spatial domains at t n and t n+1 respectively.

Let's note that generally Ω n = Ω n+1 . We denote by d the given mesh displacement field and

x n = x(t n ). Then, we want to find the optimal continuous mesh M n,ALE

L p = (M n,ALE L p (x n ))
x n ∈Ω n defined on Ω n from which we will generate a mesh at time t n that, once moved with displacement d, will be adapted to a sensor u n+1 at time t n+1 on domain Ω n+1 . We consider that u n+1 is a scalar sensor function, the extension to vector sensor functions being straightforward. The resolution of this problem leads to the optimal instantaneous ALE continuous mesh defined on Ω n minimizing the interpolation error of sensor u n+1 in L p norm at time t n+1 on Ω n+1 after being deformed by displacement d. This ALE metric field involves the gradient of the mesh transformation between t n and t n+1 to take into account the mesh deformation.

To give the result of the analysis involving spatial derivatives, it is very important to state on which domain/mesh these derivatives are computed. To this end, the following notations will be used:

• ∇ n denotes the gradient operator performed on domain Ω n , i.e., computed on mesh H n

• H n+1 denotes the Hessian operator performed on domain Ω n+1 , i.e., computed on mesh

H n+1 .
As operator H n+1 is always applied to sensor u n+1 at time t n+1 , simplified notation

H n+1 u will stand for H n+1 [u n+1 ] • M n+1 L p [u n+1
] denotes the point-wise optimal L p metric for sensor u n+1 computed on domain Ω n+1 , i.e., on mesh H n+1 . Again, simplified notation M n+1 L p will stand for

M n+1 L p [u n+1 ] • M n+1 L p = (M n+1 L p (x)) x∈Ω n+1 is the associated continuous mesh with complexity C(M n+1 L p ) = N n+1 .
Now, let us introduce φ the mapping between domains Ω n and Ω n+1 :

φ : Ω n -→ Ω n+1 x n -→ x n+1 = φ (x n ) ,
(1.49) and d the corresponding mesh displacement field, such that:

x n+1 = φ(x n ) = x n + d(x n ) .
(1.50)

Since φ is a diffeomorphism, we have, for any infinitesimal vector δx n ∈ Ω n :

δx n+1 = ∇ n φ(x n ) T δx n with ∇ n φ(x n ) = I + ∇ n d(x n ) , ∀x n ∈ Ω n . (1.51)
Finally, a . operator is defined which transports a quantity from Ω n+1 to Ω n . We note H n+1 u the Hessian of u n+1 computed on Ω n+1 and transported on domain Ω n . This mathematically writes:

H n+1 u : Ω n -→ R x n -→ H n+1 u (φ(x n )) = H n+1 u (x n ) .
To give the expression of the optimal continuous mesh, the demonstration of the result is given in [N. [START_REF][END_REF], we introduce the ALE Hessian-metric:

|H n,ALE u (x n )| = det ∇ n φ(x n ) 1 p ∇ n φ(x n ) | H n+1 u (x n )| ∇ n φ(x n ) T , (1.52)
and the expression of the optimal continuous mesh in 3D defined on Ω(t n ) is:

M n,ALE L p (x n ) = N n+1 2 3 Ω n det |H n,ALE u (x n )| p 2p+3 dx n -2 3 det |H n,ALE u (x n )| -1 2p+3 |H n,ALE u (x n )| , (1.53)
We find that we have the same expression as the optimal mesh in the stationary case (cf.

Relation (1.35)) the absolute value of the Hessian of the sensor is replaced by the ALE Hessianmetric of the sensor.

Space-time error analysis for dynamic meshes

Section 1.3.2 provides the optimal instantaneous ALE continuous mesh which takes into account the mesh deformation. Now, as precised previously, we can extend the space-time error analysis with time sub-intervals done for fixed meshes to the case of dynamic meshes. The simulation time interval is split into n adap sub-intervals. On each sub-interval, the mesh size (number of vertices) remains constant, but the mesh is deformed to follow the geometry displacement. At each time-step of the sub-interval, we want the moved mesh to be adapted to the current sensor.

The key idea to perform the error analysis is to seek for the optimal dynamic continuous mesh at the beginning of the sub-interval, this continuous mesh being optimal for the whole sub-interval when deformed, instead of seeking for the expression of the optimal continuous mesh at each instant, i.e., as a function of the time.

To perform the spatial minimization for a sub-interval, we consider the i th sub-interval [t i , t i+1 ]. Given the continuous mesh spatial complexity N i , we seek for the optimal ALE continuous mesh

M i,ALE L p = M i,ALE L p (x(t)) x∈Ω(t)
for the whole sub-interval which is solution of the following problem:

E i,ALE L p (M i,ALE L p ) = min M i t i+1 t i Ω(t) trace (M i ) -1 2 (x(t)) |H u (x(t), t)| (M i ) -1 2 (x(t)) p dx(t) dt .
(1.54) such that C(M i ) = N i . The continuous mesh spatial complexity is constant on this sub-interval.

To remove the time dependency of continuous mesh M i , we use the optimal instantaneous ALE continuous mesh. Indeed, continuous mesh M i (x(t)) x∈Ω(t) can be mapped back to Ω(t i ) using

M i,ALE (x(t i )) x∈Ω(t i )
where φ maps Ω(t i ) onto Ω(t). The interpolation error for each time t can be re-written at time t i , see [N. [START_REF][END_REF]]. As we seek for the dynamic continuous mesh at time t i which is optimal to control the interpolation error for the whole sub-interval, we can recast Error Model (1.54) into the following error model where the metric is independent of the time:

E i,ALE L p (M i ) = t i+1 t i Ω(t i ) trace (M i,ALE ) -1 2 (x(t i )) |H i,ALE u (x(t i ))| (M i,ALE ) -1 2 (x(t i )) p dx(t i ) dt ,
where the x(t i ) are in domain Ω(t i ). The time dependency in H i,ALE u is hidden in mapping φ and operator . . The previous expression can now be written on Ω(t i ):

E i,ALE L p (M i ) = Ω(t i ) trace (M i,ALE ) -1 2 (x(t i )) H i,ALE u (x(t i )) (M i,ALE ) -1 2 (x(t i )) p dx(t i ) ,
where the time-average ALE Hessian-metric is:

H i,ALE u (x(t i )) = t i+1 t i |H i,ALE u (x(t i ))| dt.
The expression of the error has exactly the same form as in Problem (1.40), thus the spatial minimization gives the same optimal metric where H i u is replaced by H i,ALE u . Then, the temporal minimization leads to the following optimal space-time ALE continuous mesh

M ALE L p = {M i,ALE L p } i=1,..,n adap : M i,ALE L p (x(t i )) = N 2 3 st n adap j=1 K j,ALE t j+1 t j τ (t) -1 dt 2p 2p+3 -2 3 t i+1 t i τ (t) -1 dt -2 2p+3 (det H i,ALE u (x(t i ))) -1 2p+3 H i,ALE u (x(t i )) ,
where the x(t i ) are in domain Ω(t i ). The ALE continuous mesh dependence in time is hidden in H i,ALE u by means of mapping φ. This way, the mesh generated at time t i is adapted to the solution at any time t > t i within the sub-interval [t i , t i+1 ] once moved with the mesh deformation displacement. We stress that preserving the number of degrees of freedom when moving the mesh is essential in this analysis.

Mesh adaptation algorithm

The optimal space-time ALE continuous mesh is designed to fit in the global fixed-point unsteady mesh adaptation algorithm described in Algorithm 3. Nevertheless, a few things have been modified to extend this algorithm to moving mesh ALE simulations. In fact, first, geometries move so mesh must be moved with all the difficulties that this implies to preserve accuracy and a good mesh quality all along the movement. All of this will be detailed more precisely in this thesis.

Conclusion

This first chapter has provided the basic wide background required by this thesis. It relies on the continuous mesh framework to exhibit an analytical expression of the optimal adapted mesh.

We observe that for the steady and the unsteady cases we have always the same expression of the optimal mesh, the difference in the error estimates resides in the choice of the Hessian-metric which is simply the Hessian of the sensor for the simplest case (eg. steady feature-based mesh adaptation) or more complex combination of Hessian-metric for the other cases.

All the improvements achieved over the theories presented here will be mentioned throughout the dissertation as well as the newest developments regarding goal-oriented adaptation for moving meshes and adjoint solver computations for moving geometries. To this end, the main difficulties of this work are the extension of the unsteady adjoint solver in the context of the moving geometries and the goal-oriented error estimate in this context.

The main contributions of this thesis are:

• the definition, the implementation and the validation of implicit time integration scheme for the ALE flow solver (Chapter 2),

• the completion with respect to several steps for the unsteady mesh adaptation for moving geometries which complement thesis [Barral 2015] (Chapter 3). In particular, the importance to update the ALE metric in the flow solver to govern the mesh optimizations during the moving mesh algorithm,

• the development of the unsteady adjoint solver in the ALE framework which requires a backward in time moving mesh algorithm consistent with the forward in time one (Chapter 4)

• the time accurate goal-oriented mesh adaptation for moving geometries and its associated error estimate (Chapter 5).

Chapter 2

Resolution of Euler equations in the ALE framework

Fluid-structure interaction (FSI) simulations are unavoidable for a wide variety of subjects.

They include simulation of aircraft flutters, ship propellers efficiency [Compère 2010], wind turbines efficiency [START_REF] Bazilevs | [END_REF]], 2D airbag deployment and balloon inflation [START_REF] Saksono | [END_REF]],

releasing of a missile [START_REF] Murman | [END_REF][START_REF] Hassan | [END_REF]].

In the acoustic field one can cite tuning forks [START_REF] Froehle | [END_REF]].

And in the field of biology, simulation of aortic valves [START_REF] Astorino | [END_REF]], cardiovascular systems [Formaggia 2009, Gerbeau 2014] or the simulation of jellyfish [Etienne 2010] also use FSI simulation. In the case of fluid-structure interaction, the framework is a movable solid whether stiff or deformable and a fluid whether gaseous or liquid. But other problems such as icing [START_REF] Tong | [END_REF], Pendenza 2014] or blast studies [START_REF] Baum | [END_REF] do not involve FSI but share a lot of common problems with FSI simulations. And all of these applications require appropriate computational methodologies customized to the specific context.

Three main types of methods are possible to treat the deformation of the domain : the multiphasic method, the mesh deformation method or the fictional domain method. The multiphasic approach initially refers to a multiphase flow then any fluid flow consisting of more than one phase or component. It is limited in the case of FSI simulation to the cases where the solid and the fluid can be described by the same equations with variable physical parameters assigned to each phase and advected during the movement of the interface. In general, the Lagrangian approach for a fluid simulation is conceivable and gives good results for some multiphase issues for example in the industry of nuclear energy [Causin 2005], but when the flow becomes complex, the interfaces displacement induces the deformation of the mesh and imposes a remeshing of the complete domain. In a realistic state, this recall occurs in very short times.

An Eulerian formulation is so better suited to the simulation of a flow fluid. On the other hand, in the context of the fluid-structure interaction, a formulation purely eulerian does not allow to effectively monitor the phenomena at the fluid-solid interface and poses the question of the law of behavior to be imposed in the cells cut by this interface.

The mesh deformation approach with an Arbitrary Lagrangian Eulerian (ALE) approach is built to avoid these drawbacks. This method allows to ensure a smooth transition between these two methods of modelling because the flow is calculated on a domain that is deformed to follow the movement of the interface (Lagrangian close to the solid) the rate of deformation does not follow that of the interior field. To our knowledge, this technique was introduced in the 70s in [Hirt 1974[START_REF] Hughes | [END_REF][START_REF] Donea | [END_REF]. Since then, so many developments have been made in that field that a complete list of them would not fit in this thesis. However, one may in particular refer to [START_REF] Nkonga | [END_REF], Baum 1994[START_REF] Farhat | [END_REF], Formaggia 2004, Mavriplis 2006[START_REF] Hassan | [END_REF][START_REF] Hay | High-order temporal accuracy for 3D finite-element ALE flow simulations[END_REF]], which mainly focus on improving temporal schemes for ALE simulations. When the displacement of the geometry is small enough, slightly deforming the original mesh [Batina 1990, Degand 2002] can generally be acceptable. However, difficulties arise when the displacement of the structure is too important: the mesh quickly becomes distorted and the numerical error due to this distortion quickly becomes too great, until the elements of the mesh finally become invalid, and the simulation has to be stopped. This is the major difficulty in this method. The computational domain of the fluid varies depending on the time according to the displacement of the solid. And even if the great development of computing capacities has made it possible to run increasingly complex simulations, the engineers are still far from performing such simulations, largely due to the difficulty of handling the moving meshes induced by large deformations of the moving geometries. Specific strategies need to be developed to deal with the displacement of moving boundary problems. For now, this method is interesting for large displacements problems if the fluid domain is quite often remeshed or coupled with mesh adaptation.

It's to get around the question of remeshing that the fictional domain methods have been developed also called embedded methods. These embedded boundary approach [Bruchon 2009, Löhner 2001] uses meshes that are not body-fitted at all : the bodies are embedded in a fixed grid, and techniques such as level-sets are used to recover their moving boundaries. The inconvenient of this method especially is the low CFL condition for unstructured vertex-centered method.

All three approaches have their own strengths and weaknesses. The chosen strategy in this thesis to handle moving geometries is the body-fitted framework where the geometry is explicitly represented inside the mesh. This choice is motivated by its compliance with anisotropic mesh adaptation and by the fact that body-fitted approaches are the more accurate methods to simulate viscous flows (even if only inviscid flow will be considered in this thesis). For body-fitted moving mesh simulations, the whole mesh must be deformed to follow the moving boundaries. And for this strategy to be successful, it is essential to take into account this move into the numerical scheme. The ALE framework is the chosen strategy in this work and it has been studied using explicit time integration in [N. [START_REF][END_REF]. However, the use of an explicit time integration schemes leads to an important restriction of the admissible time step because it is specified by the size of small elements (usually the smallest element) in the mesh.

Explicit schemes are widely used for blast problems (violent flows), LES, aeroacoustic, ... to control efficiently dissipation and dispersion errors that can spoil the overall solution accuracy [J. [START_REF] Berland | [END_REF], Löhner 2001, N. Gourdain 2009].

To reduce the time step restriction due to the CFL condition, improved time advancing schemes exist such as multi-rate approaches [Constantinescu 2007, B. Seny 2014].However, the efficiency these approaches is highly reduced when coupled with mesh adaptation because they increase the time steps in regions where mesh adaptation generates a coarse mesh.

The other strategy to reduce the time step restriction due to the CFL condition is to use implicit time integration schemes that are theoretically unconditionally stable (in practice, implicit schemes have a larger region of stability compared to explicit schemes, leading to larger time steps during time marching). These schemes are very attractive for quasi static flows, aeronautical smooth flow simulations (vortex shedding, ...), viscous flows, ... But, to see a gain in CPU time, the considered implicit time-step should be sufficiently large to overcome the over cost of the implicit time integration without impacting the solution accuracy.

Section 2.1 presents a reminder of moving mesh strategy [Alauzet 2014]. Then, the ALE flow solver Wolf is described, Section 2.2. Section 2.2.3 presents the different implicit time integration approaches that have been developed in this thesis. And finally, comparaison between explicit and implicit time integration schemes are provided in the section 2.4.

Mesh-connectivity-change moving mesh strategy

To handle moving boundaries, we adopt a body-fitted approach, with a single mesh: the inner vertices of the mesh are moved following the moving boundaries to preserve the validity of mesh (i.e. to prevent the mesh from getting tangled). Our strategy involves two main parts:

• Computing the mesh deformation: inner vertices are assigned a trajectory depending on the displacement of the boundaries, and thus a position for future time steps.

• Optimizing the mesh: the trajectories computed in the mesh deformation phase are corrected, and the connectivity of the mesh is modified to preserve the quality of the mesh.

This strategy, recalled below, has proven to be very powerful in 3D [Alauzet 2014], since large displacement of complex geometries can be performed while preserving a good mesh quality without any global remeshing (i.e. without ever generating a whole new mesh).

To our knowledge, very few examples of ALE solvers coupled with connectivity-change moving mesh techniques can be found in the literature. In [Kucharik 2008] a conservative interpolation is proposed to handle the swaps. In [START_REF] Guardone | [END_REF], Olivier 2011b] an ALE formulation of the swap operator is built. However, these studies are limited to 2D. In [Barral 2015] a method for computations in 3D for moving geometries coupled with connectivity-change is done and the numerical solver used is an explicit time integration method.

Linear elasticity mesh deformation method

During the mesh deformation step, a displacement field is computed for the whole computational domain, given the displacement of its boundaries. Trajectories can thus be assigned to inner vertices, or in other words, positions at a future solver time step. Several techniques can be found to compute this displacement field:

• implicit or direct interpolation [de Boer 2007, Luke 2012], in which the displacement of the inner vertices is a weighted average of the displacement of the boundary vertices,

• solving PDEs -the most common of which being Laplacian smoothing [Löhner 1998], a spring analogy [Degand 2002] and a linear elasticity analogy [START_REF] Baker | [END_REF]].

It is this last method that we selected, due to its robustness in 3D [START_REF][END_REF]]. The computational domain is assimilated to a soft elastic material, which is deformed by the displacement of its boundaries.

The inner vertices movement is obtained by solving an elasticity-like equation with a P 1 Finite Element Method (FEM):

div(σ(E)) = 0 , with E = ∇d + T ∇d 2 , (2.1)
where σ and E are respectively the Cauchy stress and strain tensors, and d is the Lagrangian displacement of the vertices. The Cauchy stress tensor follows Hooke's law for isotropic homogeneous medium, where ν is the Poisson ratio, E the Young modulus of the material and λ, µ are the Lamé coefficients:

σ(E) = λ trace(E) I d + 2 µ E or E(σ) = 1 + ν E σ - ν E trace(σ) I d .
In our context, ν is typically chosen of the order of 0.48, which corresponds to a nearly incompressible material1 . Dirichlet boundary conditions are used and the displacement of vertices located on the domain boundary is strongly enforced in the linear system. The linear system is solved by a Conjugate Gradient algorithm coupled with an LU-SGS pre-conditioner. An advantage of elasticity-like methods is the opportunity they offer to adapt the local material properties of the mesh, especially its stiffness, according to the distortion and efforts borne by each element. In particular, the stiffness of the elements is increased for small elements, in order to limit their distortion. More details can be found in [Alauzet 2014].

Improving mesh deformation algorithm efficiency

The computation of the mesh deformation -here the solution of a linear elasticity problemis known to be an expensive part of dynamic mesh simulations, and the fact that it is usually performed at every solver time step makes it all the more so.

We propose to combine several techniques to improve the time efficiency of this step. Some regions are rigidified, more specifically a few layers around tiny complex details of the moving bodies, with very small elements. They are moved with exactly the same rigid displacement as the corresponding body, thus avoiding very stiff elements in the elasticity matrix. On the other hand, the elasticity can be solved only on a reduced region, if the domain is big compared to the displacement. A coarse mesh can also be used to solve the elasticity problem, the displacement of the vertices then being interpolated on the computational mesh.

The major improvement we proposed is to reduce the number of mesh deformation computations: the elasticity problem is solved for a large time frame of length ∆t instead of doing it at each solver time step δt. While there is a risk of a less effective mesh displacement solution, it is a worthwhile strategy if our methodology is able to handle large displacements while preserving the mesh quality. Solving the previously described mesh deformation problem once for large time frame could be problematic in the case of: (i) curved trajectories of the boundary vertices and (ii) accelerating bodies. To enhance the mesh deformation prescription, accelerated-velocity curved, i.e., high-order, vertex trajectories are computed.

The paths of inner vertices can be improved if a constant acceleration a is provided to each vertex in addition to its speed, which results in an accelerated and curved trajectory. During time frame [t, t + ∆t], the position and the velocity of a vertex are updated as follows:

x(t + δt) = x(t) + δt v(t) + δt 2 2 a v(t + δt) = v(t) + δt a .
Prescribing a velocity vector and an acceleration vector to each vertex requires solving two elasticity systems. For both systems, the same matrix, thus the same pre-conditioner, is considered. Only boundary conditions change. If inner vertex displacement is sought for time frame [t, t + ∆t], boundary conditions are imposed by the location of the body at time t + ∆t/2 and t + ∆t. These locations are computed using body velocity and acceleration. Note that solving the second linear system is cheaper than solving the first one as a good prediction of the expected solution can be obtained from the solution of the first linear system. Now, to define the trajectory of each vertex, the velocity and acceleration are deduced from evaluated middle and final positions:

∆t v(t) = -3 x(t) + 4 x(t + ∆t/2) -x(t + ∆t) ∆t 2 2 a = 2 x(t) -4 x(t + ∆t/2) + 2x(t + ∆t) .
In this context, it is mandatory to make sure that the mesh remains valid for the whole time frame [t, t + ∆t], which is done by computing the sign of the volume of the elements all along their path [Alauzet 2014].

Local mesh optimization

In order to preserve the mesh quality between two mesh deformation computations, it has been proposed [Alauzet 2014] to couple mesh deformation with local mesh optimization using smoothing and generalized swapping to efficiently achieve large displacement in moving mesh applications. Connectivity changes are really effective in handling shear and removing highly skewed elements. Here, we briefly recall the mesh optimization procedure.

Mesh optimizations are performed regularly to preserve the quality of the mesh. The optimization procedure consists in one or several passes of vertex smoothing and one pass of generalized swapping. For 3D meshes, the quality of an element is measured in terms of the element shape by the quality function:

Q(K) = √ 3 216 6 i=1 2 M (e i ) 3 2 |K| M ∈ [1, +∞] , (2.2)
where M (e) and |K| M are edge length and element volume in metric M. Q(K) = 1 corresponds to a perfectly regular element and Q(K) < 2 corresponds to excellent quality elements, while a high value of Q(K) indicates a nearly degenerated element. For non-adapted meshes, the identity matrix I 3 is chosen as metric tensor.

Mesh smoothing. The first mesh optimization tool is vertex smoothing which consists in relocating each vertex inside its ball of elements, i.e., the set of elements having P i as their vertex. For each tetrahedron K j of the ball of P i , a new optimal position P opt j for P i can be proposed to form a regular tetrahedron:

P opt j = G j + 2 3 n j (n j ) ,
where F j is the face of K j opposite vertex P i , G j is the center of gravity of F j , n j is the inward normal to F j and (n j ) the length of n j . The final optimal position P opt i is computed as a weighted average of all these optimal positions {P opt j } K j ⊃P i , the weight coefficients being the quality of K j . This way, an element of the ball is all the more dominant if its quality in the original mesh is bad. Finally, the new position is analyzed: if it improves the worst quality of the ball, the vertex is directly moved to its new position.

Generalized edge/face swapping. The second mesh optimization tool to improve mesh quality is generalized swapping/local-reconnection. Let α and β be the two tetrahedra vertices opposite the common face P 1 P 2 P 3 . Face swapping consists of suppressing this face and creating the edge e = αβ. In this case, the two original tetrahedra are deleted and three new tetrahedra are created. This swap is called 2 → 3. The reverse operator can also be defined by deleting three tetrahedra sharing such a common edge αβ and creating two new tetrahedra sharing face

P 1 P 2 P 3 . This swap is called 3 → 2.
A generalization of this operation exists and acts on shells of tetrahedra [Alauzet 2014[START_REF] Frey | [END_REF]. For an internal edge e = αβ, the shell of e is the set of tetrahedra having e as common edge. The different edge swaps are generally denoted n → m where n is the size of the shell and m is the number of new tetrahedra. In this work, edge swaps 3 → 2, 4 → 4, 5 → 6, 6 → 8 and 7 → 10 have been implemented. In our algorithm, swaps are only performed if they improve the quality of the mesh. These operations are well-know in the field of mesh generation [J. F. Thompson 1999[START_REF] Frey | [END_REF]], but are not necessarily efficient in the moving mesh context. Notably, performing too many of them results in slow codes, whereas the use of bad quality functions results in poor quality meshes. The interest of the method used here lies in how and when optimizations are performed. The mesh optimizations are performed entity by entity, and only when they are needed. Smoothing is performed for every vertex above a prescribed quality threshold ( the quality of a vertex is the maximal element quality of the elements in its ball 2 ), provided the new position increases the quality of the worst element in its ball. Swaps are analyzed element by element if and only if the element quality is above a prescribed quality threshold. Then, swaps are performed if the quality of the considered elements decrease after the connectivity change. To avoid multiple identical checks, tetrahedra are treated in quality order from the worst to the best one. In the moving mesh context, the optimization operation is performed only if it verifies quality criteria on the current position of the mesh and on the final position given by the mesh deformation. This avoid to perform a mesh optimization which will degrade the quality of the mesh in the future. A key to performing efficient swaps in the moving mesh context is to allow a slight quality degradation in the future. Details on this optimization step can be found in [Alauzet 2014].

Handling of boundaries

The mesh of the boundaries is moved rigidly, and the vertices are not generally moved on the surface (no displacement in the tangential directions). However, in some cases, such as when a body is moving very close to the bounding box of the domain, it can be useful to move the vertices of the bounding box as well. In this case, we can allow tangential displacement on the boundary. For curved boundary, vertices are re-projected on the surface after the displacement.

To do so, the displacements along the tangential axes are simply considered as new degrees of freedom. For instance, for a plane (x,y), the displacements along the x-axis and the y-axis are considered as degrees of freedom and are added to the elasticity system. The displacement along the z-axis is still set to 0, and thus is not added to the system. This can be generalized for any tangent plane.

Moving mesh algorithm

The overall connectivity-change moving mesh algorithm is described in Algorithm 5, where the different phases described above are put together. When coupled with a flow solver (see Sections 2.2), the flow solver is called after the optimization phase. In this algorithm, H stands for meshes, S for solutions, Q for quality (see Relation (2.2)), d |∂Ω h for the displacement on the boundary, and v and a for speed and acceleration. ∆t and δt are time steps whose meaning is detailed below.

2 The elements ball of a vertex is the set of elements sharing that vertex. Else T els = t + ∆t 3. Moving mesh stage, with mesh optimizations and solver solution In this algorithm, three time steps appear: a large one ∆t for the mesh deformation computation, a smaller one δt opt corresponding to the steps where the mesh is optimized, and the solver time step δt solver . The first one, ∆t, is currently set manually at the beginning of the computation. After each mesh deformation solution, the quality of the mesh in the future is analyzed: if the quality is too low, the mesh deformation is problem is solved again with a smaller ∆t (cf step 2.(d)). Moreover, if the mesh quality degrades, a new mesh deformation solution is computed (cf step 3.(g)). The second one is computed automatically, using the CF L geom parameter as described below. Determining the third one will be discussed in Section 2.2. If the solver time step is greater than the optimization time step, then the solver time step is truncated to follow the optimizations. If the solver time step is smaller than the optimization time step -which is almost always the case -, several (many) iterations of the flow solver are performed between two optimization steps.

While (t < T els ) (a) δt opt = Get optimization time step H k , v, CF L geom (b) δt solver = Get solver time step H k , S k , v, CF L (c) δt = min(δt opt , δt solver ) (d) If t > t opt i. H k = Swaps optimization H k , Q swap target ii. v opt = Vertex smoothing H k , Q smoothing target , Q max iii. t opt = t + δt opt (e) S k+1 = Solve Equation of State H k , S k , δt, v, v opt

Moving mesh time steps

A good restriction to be imposed on the mesh movement in order to limit the apparition of flat or inverted elements is that vertices cannot cross too many elements on a single move between two mesh optimizations. Therefore, a geometric parameter CF L geom is introduced to control the number of stages used to perform the mesh displacement between t and t + ∆t. If CF L geom is greater than one, the mesh is authorized to cross more than one element in a single move. In practice, CF L geom is usually set between 1 and 8. The moving geometric time step is given by:

δt opt = CF L geom max P i h i v i , (2.3)
where h i is the smallest height of all the elements in the ball of vertex P i and v i the mesh velocity at vertex P i . In practice, when coupled with a flow solver, the actual time step is the minimum between the flow solver time step and the geometric one.

Arbitrary Lagrangian Eulerian flow solver

This section discusses in details the implementation of the Arbitrary Lagrangian Eulerian (ALE) flow solver Wolf that has been coupled to the moving mesh process described in Algorithm 5.

After presenting the compressible Euler equations in the ALE framework, we describe the second order spatial discretization. As regards the temporal discretization, we first discuss the geometric conservation law (GCL) and then we present the explicit and implicit time integration schemes. This thesis has focused on the deriving and implementing implicit schemes in the framework of Mavriplis and Yang [Mavriplis 2006]. They are compared to explicit schemes in the last section.

Euler equations in the ALE framework

We consider the 3D unsteady compressible Euler equations for a Newtonian fluid in their ALE formulation set in the space-time computational domain

Q = Ω(t) × [0, T ],
where T is the (positive) maximal time. The ALE formulation allows the equations to take arbitrary motion of the (domain) mesh into account. Assuming that the gas is perfect, inviscid and that there is no thermal diffusion, the ALE formulation of the equations is written, for any arbitrary closed volume C(t) of boundary ∂C(t) moved with mesh velocity w:

d dt C(t) W dx + ∂C(t) (F(W ) -W ⊗ w) • n ds = C(t) F ext dx ⇐⇒ d dt C(t) W dx + ∂C(t) (F (W ) -W (w • n)) ds = C(t) F ext dx ,
where

             W = (ρ, ρu, ρE) T is the conservative variables vector F(W ) = (ρu, ρuu + pe 1 , ρvu + pe 2 , ρwu + pe 3 , (ρE + p)u) is the flux tensor F (W ) = F(W ) • n = (ρη, ρuη + pn 1 , ρvη + pn 2 , ρwη + pn 3 , (ρE + p)η) T F ext = (0, ρ f ext , ρu • f ext )
T is the contribution of the external forces, and we have noted ρ the density of the fluid, p the pressure, u = (u, v, w) its Eulerian velocity, η = u • n , q = u , ε the internal energy per unit mass, E = 1/2 q 2 + ε the total energy per unit mass, e = (e 1 , e 2 , e 3 ) the canonical basis, f ext the resultant of the volumic external forces applied on the particle and n = (n 1 , n 2 , n 3 ) the outward normal to interface ∂C(t) of C(t).

Spatial discretization

The spatial discretization of the governing equations is based on a vertex-centered finite volume formulation on unstructured meshes. This formulation consists in associating with each vertex of the mesh a control volume (or Finite-Volume cell) built by the rule of medians. Wolf proposes several choices of approximate Riemann solvers to compute numerical fluxes. Second-order space accuracy is achieved through a piecewise-linear extrapolation based on the Monotonic Upwind Scheme for Conservation Law (MUSCL) procedure with a particular edge-based formulation with upwind elements. Specific slope limiters are employed to damp or eliminate spurious oscillations that may occur in the vicinity of discontinuities. The main difference when translating these schemes from the standard formulation to the ALE formulation is the addition of the mesh velocities in the wave speeds of the Riemann solver problem.

Dual mesh construction

The spatial domain Ω(t) is discretized by a tetrahedral unstructured mesh H. The vertexcentered finite volume formulation consists in associating a control volume denoted C i (t) with each vertex P i of the mesh and at each time t. The dual finite volume cell mesh is built by the rule of medians. Discretized domain Ω h (t) can be written as the union of the elements or the union of the finite volume cells:

Ω h (t) = N K i=1 K i (t) = N V i=1 C i (t) ,
where N K is the number of elements and N V the number of vertices. Note that the dual mesh (composed of cells) is built in a preprocessing step. Consequently, only a simplicial mesh is needed in the input.

In 2D, this standard method consists in building cells bounded by segments of medians.

Each triangle is split into three quadrangles (one associated with each one of its three vertices).

The quadrangle's four vertices of the triangle associated with vertex P i are:

• M i , M j : the middle of the two edges incident to P i ,

• G: the gravity center of the triangle,

• P i : the considered vertex.

Median cell C i (t) of vertex P i is the union of all quadrangles of the triangles surrounding P i .

An exemple of median cell is represented in Figure 2.2 (left). Consequently, at each edge P i P j of the mesh is associated a cell interface ∂C ij (t) which is common boundary between the two

neighboring cells C i (t) and C j (t), i.e., ∂C ij (t) = ∂C i (t)∩∂C j (t). Note that cell interface ∂C ij (t)
is composed of two bi-segments (or one bi-segment for a boundary edge).

In 3D each tetrahedra is subdivided into four hexahedra (one associated with each one of its four vertices). The eight vertices of the hexahedron associated with a point P i are given by:

• M i , M j , M k : the middle of the three edges incident to P i

• Gf i , Gf j , Gf k : the gravity centers of the three faces containing P i

• G: the gravity center of the tetrahedra

• P i : the considered vertex.
Such an hexahedron is represented in Figure 2.2 (right). Median cell C i (t) of vertex P i is the union of all hexahedra of the tetrahedra surrounding P i . Consequently, at each edge P i P j of the mesh is associated a cell interface ∂C ij (t) which is common boundary between the two neighboring cells C i (t) and C j (t), i.e., ∂C ij (t) = ∂C i (t) ∩ ∂C j (t). In three dimensions this cell interface ∂C ij (t) is composed of many triangular facets (two triangular facts per tetrahedra in the shell of edge P i P j ).

Edge-based Finite-Volume discretization

Assuming there is no contribution of the external forces, based on a finite volume formulation, the compressible Euler equations are integrated on each finite volume cell C i (t) (using the Green formula):

d dt (|C i (t)|W i (t)) + F(W i (t)) = 0 , (2.4)
where W i (t) is the mean value of the solution W on cell C i (t) and F i (W (t)) is the numerical convective flux term:

F(W i (t)) = ∂C i (t) (F(W i (t)) -W i (t) ⊗ w i ) • n i ds , (2.5)
where n i is the outer normal to the finite volume cell surface ∂C i (t) as depicted in Figure 2.3, F is the convective term flux functions, and w i is the mesh velocity at vertex P i . The integration of convective fluxes F of Equation (2.5) is done by decomposing the cell boundary into many interfaces ∂C ij composed of facets:

F(W i (t)) = P j ∈V(P i ) F ij (t) • ∂C ij (t) n i ds -W ij (t) ∂C ij (t) (w i • n i ) ds ,
where V(P i ) is the set of all neighboring vertices linked by an edge to P i , and W ij (t) and F ij (t)

represents the constant value of W and F(W ) at interface ∂C ij (t). The flow is calculated using a numerical flux function to approximate the flux at cell interface ∂C ij (t), denoted by Φ ij :

Φ ij (t) = Φ ij (W i (t), W j (t), n ij (t), σ ij (t)) = F ij (t) • n ij (t) -W ij (t)σ ij (t) ,
where

• n ij (t) = 1 |∂C ij (t)| ∂C ij (t) n i dγ is the outward normalized normal (with respect to cell C i ) of cell interface ∂C ij , • σ ij (t) = 1 |∂C ij (t)| ∂C ij (t) w ij (t) • n ij (t)ds is the normal velocity of cell interface ∂C ij (t).
We finally get the following semi-discretization at P i :

d dt (|C i (t)|W i (t)) + P j ∈V(P i ) |∂C ij (t)| Φ ij (W i (t), W j (t), n ij (t), σ ij (t)) = 0 , (2.6)
The numerical flux function approximates the hyperbolic terms on the common boundary ∂C ij (t). We notice that the computation of the convective fluxes is performed monodimensionally in the direction normal to the boundary of the finite volume cell. Therefore, the numerical calculation of the flux function Φ ij (t) at the interface ∂C ij (t) can be achieved by solving, at each time step, a one-dimensional Riemann problem in the direction of the normal n ij by means of an approximate Riemann solver.

P j P i M i M j downwind triangle K j upwind triangle K i n 1 ij n ij ( t) n 2 ij P i P j W ij W ji Figure 2
.3: Illustration of finite volume cells configuration around an edge P i P j in 2D: two neighboring cells C i and C j , their common interface ∂C ij and the upwind triangles K i and K j associated with edge P i P j .

Three approximate Riemann solvers are available in Wolf to compute the numerical flux Φ:

Rusanov, Roe and HLLC. In the following, we describe their extensions to the ALE framework.

For sake of clarity, the dependency in t of the variable is omitted.

Numerical flux computation

ALE Rusanov approximate Riemann solver. In the Rusanov's approach [Rusanov 1961], the upwinding term is defined thanks to the maximal spectral radius of the Jacobian matrix

∂(F(W ) • n) ∂W
given by: λ = |u • n| + c. The flux resolution uses the following formulation of Φ:

Φ Rusanov (W i , W j , n ij ) = F(W i ) + F(W j ) 2 • n ij + |λ ij | W i -W j 2 ,
where

|λ ij | = max(|u i .n ij | + c i , |u j .n ij | + c j ) and c = γp ρ
is the speed of sound. For ALE simulation, the flux formulation takes into account the mesh displacement and is given by:

Φ Rusanov ale (W i , W j , n ij , σ ij ) = F(W i ) • n ij -σ ij W i + F(W j ) • n ij -σ ij W j 2 + | λ ij | W i -W j 2 = F i + F j 2 -σ ij W i + W j 2 + | λ ij | W i -W j 2 where | λ ij | = max(|u i .n ij | + c i -σ ij , |u j .n ij | + c j -σ ij ) and the notation F k = F(W k ) • n ij .
ALE Roe approximate Riemann solver. In the Roe's approach [Roe 1981], the upwinding term is defined by the Jacobian matrix

A(W ) = ∂(F(W ) • n) ∂W . The eigenvalues of A(W ) are real and given by |u • n|, |u • n| + c and |u • n| -c, thus A(W ) is diagonalisable.
In the context of the Euler equations, the hyperbolic flux is homogeneous of order one leading to the property:

F(W ) • n = A(W ) W ,
that allows non-oscillatory conservative schemes to be built. The flux resolution uses the formulation of Φ introduced by Roe:

Φ Roe (W i , W j , n ij ) = F(W i ) + F(W j ) 2 • n ij + | Ã(W i , W j )| W i -W j 2 ,
where à is the Jacobian matrix evaluated for the Roe's average variables and for all diagonalisable matrix A = P ΛP -1 we have denoted |A| = P |Λ|P -1 . Let W i and W j be the two states, then the Roe's average variables are given by:

ρ = √ ρ i ρ j ũ = √ ρ l u i + √ ρ r u j √ ρ i + √ ρ j ṽ = √ ρ l v i + √ ρ r v j √ ρ i + √ ρ j w = √ ρ l w i + √ ρ r w j √ ρ i + √ ρ j h = √ ρ l h i + √ ρ r h j √ ρ i + √ ρ j
where h is the enthalpy per mass unit: h = ρE+p ρ , and from which we get c2 = (γ -1) ( h -1 2 q2 ) where q2 = ũ2 + ṽ2 + w2

For ALE simulation, the flux formulation takes into account the mesh displacement:

Φ Roe ale (W i , W j , n ij , σ ij ) = F(W i ) • n ij -σ ij W i + F(W j ) • n ij -σ ij W j 2 + | Ã(W i , W j ) -σ ij I| W i -W j 2 = F i + F j 2 -σ ij W i + W j 2 + | Ã(W i , W j ) -σ ij I| W i -W j 2
where I is the identity matrix and with the notation

F k = F(W k ) • n ij .
HLLC numerical flux. The idea of the HLLC flow solver is to consider locally a simplified Riemann problem with two intermediate states depending on the local left and right states [Batten 1997]. The simplified solution to the Riemann problem consists of a contact wave with a velocity S M and two acoustic waves, which may be either shocks or expansion fans. The acoustic waves have the smallest and the largest velocities (S i and S j , respectively) of all the waves present in the exact solution. If S i > 0 then the flow is supersonic from left to right and the upwind flux is simply defined from F (W i ) where W i is the state to the left of the discontinuity. Similarly, if S j < 0 then the flow is supersonic from right to left and the flux is defined from F (W j ) where W j is the state to the right of the discontinuity. In the more difficult subsonic case when S i < 0 < S j we have to calculate F (W i ) or F (W j ). Consequently, the HLLC flux is given by:

Φ HLLC (W i , W j , n ij ) =                F(W i ) • n ij if S i > 0 F(W i ) • n ij if S i ≤ 0 < S M F(W j ) • n ij if S M ≤ 0 ≤ S j F(W j ) • n ij if S j < 0 .
where W i and W j are evaluated as follows. We denote by

η k = u k • n ij . Assuming that η = η i = η j = S M
, the following evaluations are proposed [Batten 1997]:

W i =                                      ρ i = ρ i S i -η i S i -S M p i = p = ρ i (η i -S i ) (η i -S M ) + p i (ρu) i = (S i -η i ) ρu i + (p -p i ) n ij S i -S M (ρE) i = (S i -η i ) ρE i -p i η i + p S M S i -S M , W j =                                        ρ j = ρ j S j -η j S j -S M p j = p = ρ j (η j -S j ) (η j -S M ) + p j (ρu) j = (S j -η j ) ρu j + (p -p j ) n ij S j -S M (ρE) j = (S j -η j ) ρE j -p j η j + p S M S j -S M ,
A key feature of this solver is in the definition of the three waves velocity. For the contact wave we consider:

S M = ρ j η j (S j -η j ) -ρ i η i (S i -η i ) + p i -p j ρ j (S j -η j ) -ρ i (S i -η i ) ,
and the acoustic wave speeds based on the Roe average variable as defined above:

S i = min(η i -c i , η -c) and S j = max(η j + c j , η + c)
where . are Roe average variables [Roe 1981]. With such waves velocities, the approximate HLLC Riemann solver has the following properties. It automatically (i) satisfies the entropy inequality, (ii) resolves isolated contacts exactly, (iii) resolves isolated shocks exactly, and (iv) preserves positivity.

The methodology provided by Batten [Batten 1997] can be extended to the Euler equations in their ALE formulation:

Φ HLLC ale (W i , W j , n ij , σ ij ) =                F (W i ) • n ij -σ ij W i if S i > σ ij F (W i ) • n ij -σ ij W i if S i ≤ σ ij < S M F (W j ) • n ij -σ ij W j if S M ≤ σ ij ≤ S j F (W j ) • n ij -σ ij W j if S j < σ ij
It is important to note that the considered waves speed (S i , S j and S M ) and the star states are the classical one and that the waves speed are compared in the moving frame.

Second-order accurate version

The previous formulation reaches at best a first-order spatial accuracy. A MUSCL type reconstruction method has been designed to increase the order of accuracy of the scheme. This method was introduced by Van Leer in a series of papers, see for instance [Leer 1972]. The idea is to use extrapolated values W ij and W ji instead of W i and W j at the interface ∂C ij to evaluate the flux:

Φ ij = Φ ij (W ij , W ji , n ij , σ ij ) ,
and to achieve second-order accuracy.

In the implementation, we do not extrapolate directly the conservative variable W . In fact, the primitive variables U = (ρ, u, v, w, p) are extrapolated to guarantee the positivity of the density and the pressure, then the conservative variables are reconstructed from these values.

The extrapolation steps are the following:

In the following, we present how the MUSCL gradients are computed to design secondorder accurate scheme with a fourth-order numerical dissipation. In contrast to the original Algorithm 6 MUSCL extrapolation 1. Get primitive variables U i and U j from conservatives ones W i and W j 2. Compute primitive variable MUSCL gradients ∇U ij and ∇U ji 3. Extrapolates the primitive variables:

U ij = U i + 1 2 (∇U ) ij • --→ P i P j and U ji = U j + 1 2 (∇U ) ji • --→ P j P i .
4. Get conservative variables W ij and W ji from primitive ones U ij and U ji MUSCL approach, the approximate "slopes" (∇U ) ij and (∇U ) ji are defined for any edge and are obtained using a combination of centered, upwind and nodal gradients.

The centered gradient related to edge --→ P i P j , is defined implicitly along edge --→ P i P j via relation:

(∇U ) C ij --→ P i P j = U j -U i .
(2.7)

Upwind and downwind gradients, which are also related to edge --→ P i P j , are computed using the upstream and downstream triangles/tetrahedra associated with this edge. These triangles/tetrahedra are respectively denoted K ij and K ji . K ij (resp. K ji ) is the unique triangle/tetrahedron of the ball of P i (resp. P j ) whose opposite face is crossed by the straight line prolongating edge --→ P i P j , see Figures 2.3 (2D) and 2.4 (3D). Upwind and downwind gradients of edge --→ P i P j are then defined as:

(∇U ) U ij = (∇U ) | K ij and (∇U ) D ij = (∇U ) | K ji ,
where

∇U | K = P i ∈K (∇φ P i ⊗ U i )
is the P 1 -Galerkin gradient on element K and φ P i is the usual basis function associated with P i . Parametrized nodal gradients are built by introducing the β-scheme:

(∇U ) ij --→ P i P j = (1 -β) (∇U ) C ij --→ P i P j + β (∇U ) U ij --→ P i P j , (∇U ) ji --→ P i P j = (1 -β) (∇U ) C ij --→ P i P j + β (∇U ) D ij --→ P i P j ,
where β ∈ [0, 1] is a parameter controlling the amount of upwinding. For instance, the scheme is centered for β = 0 and fully upwind for β = 1.

The most accurate β-scheme is obtained for β = 1/3, also called the V4-scheme. This scheme is third-order for the two-dimensional linear advection problem on structured triangular

M j M i P i P j K ij K ji Figure 2
.4: Downstream K ij and upstream K ji tetrahedra associated with edge --→ P i P j .

meshes. In our case, for the non-linear Euler equations on unstructured meshes, a second-order scheme with a fourth-order numerical dissipation is obtained [Debiez 2000]. The high-order gradients are given by:

(∇U ) V 4 ij --→ P i P j = 2 3 (∇U ) C ij --→ P i P j + 1 3 (∇U ) U ij --→ P i P j , (∇U ) V 4 ji --→ P i P j = 2 3 (∇U ) C ij --→ P i P j + 1 3 (∇U ) D ij --→ P i P j .

Limiter function

The aforementioned MUSCL scheme is not monotone and can be a source of spurious oscillations. These oscillations can affect the accuracy of the final solution or simply end the computation because (for instance) of negative pressures. A widely used technique for addressing this issue is to guarantee the TVD property in 1D [Harten 1983] or the LED property in 2D/3D of the scheme, which ensures that the extrapolated values W ij and W ji are not invalid.

To guarantee the TVD or the LED properties, limiting functions are coupled with the previous high-order gradient evaluations. The gradient is substituted by a limited gradient denoted (∇W ) lim ij . The choice of the limiting function is crucial as it directly affects the convergence of the simulation. Two limiters are possible with the V4-scheme.

Piperno Limiter. This limiter has been proposed by [START_REF] Piperno | Criteria for the design of limiters yielding efficient high resolution TVD schemes[END_REF]]. Here, we present slightly modified version which is less dissipative. This limiter is expressed in a factorized form,

(∇U ) Lim ij --→ P i P j = (∇U ) C ij --→ P i P j ψ P I (∇U ) C ij --→ P i P j (∇U ) V 4 ij --→ P i P j , with
ψ P I (R) = ( 1 3 + 2 3 R)        3 1 R 2 -6 1 R + 19 1 R 3 -3 1 R + 18 if R < 1 1 + ( 3 2 1 R + 1)( 1 R -1) 3 if R ≥ 1 , where R = (∇U ) C ij --→ P i P j (∇U ) V 4 ij --→ P i P j .
Koren-Dervieux Limiter. The Koren-Dervieux limiter [Cournède 2006, Koren 1993] is a three-entry limiter which is a generalization of the SuperBee limiter for the V4-scheme. It reads:

(∇U ) lim ij --→ P i P j = ψ KD ∇U D ij --→ P i P j , ∇U C ij --→ P i P j , ∇U V 4 ij --→ P i P j (2.8) with ψ KD (a, b, c) =    0 if ab ≤ 0 sign(a) min(2 |a|, 2 |b|, |c|) otherwise .
The operator ψ KD defined above is applied component by component.

Boundary conditions

Boundary conditions are computed vertex-wise. Several conditions (inflow, outflow, symmetry, ...) are available in Wolf. We refer to [Alauzet 2010a] for their description. As, in this thesis, we are only concern with inviscid flows, the only boundary condition which is modified in the ALE framework is the slipping boundary condition applied to moving bodies. This boundary condition has to take into account the displacement of the body. To this end, we impose weakly 3

u i • n i = σ i , (2.9)
where n i is the unitary boundary face normal and σ i is the boundary face velocity. The standard ALE slipping boundary flux of vertex P i reduces to:

Φ slip ale (W i , n i , σ i ) =      0 -p i n i n i -p i σ i      , (2.10)
where p i is the vertex pressure and n i =

F j P i |F j |n F j F j P i |F j |
is the mean outward normal of the boundary interface. The F j are the mesh boundary faces (edges/triangles). This flux is integrated over the boundary interface area

|∂C i | Γ = F j P i 1 3 |F j |, see Section 2.2.3.
However, when the second-order numerical schemes is considered, such a boundary condition creates oscillations in the density and the pressure when shock waves impact normally the boundary, see Figure 2.5. We thus prefer considering a mirror state and apply an approximate Riemann solver to diminish these oscillations. We thus have to evaluate the flux between the state on the boundary W and the ALE mirror state W :

W i =      ρ i ρu i ρE i      and W i =      ρ i ρu i -2 ρ i (u i • n i -σ i ) n i ρE i -2 ρ i σ i (u i • n i -σ i )     
as the mirror state verifies

p i = p i , c i = c i , u i • n i = 2 σ i -u i • n i and h i = h i -2 σ i (u i • n i -σ i ) .
To evaluate the boundary flux, we consider the HLLC approximate Riemann solver between the state and the mirror state:

Φ slip ale (W i , n i , σ i ) = Φ HLLC ale W i , W i , n i , σ i . (2.11)
Note that by definition we have

Φ HLLC ale (W i , W i , n i , σ i ) = F(W i ) -σ i W i .
Thus, if Condition (2.9) is satisfied, then W i = W i and the flux Φ slip ale (W i , n i , σ i ) simplifies to the form in Relation (2.10). In general, this condition is not satisfied, so we use Relation (2.11).

Figure 2.5: Exemple of a piston problem where a shock wave impacts the pistion. The plot represents the pressure evolution in time at a sensor on the piston for the first-order scheme (blue curve) and the second-order scheme using the classical slip boundary condition (green curve) and the Riemann slip boundary condition (red curve).

Time discretization

Temporal discretization is a more complex matter. In this section, we first recall the Geometric Conservation Law (GCL), and then, we described explicit time integration schemes. Finally, we presents the implicit time integration schemes that have been developed in this thesis. All the time discretizations are compliant with the Discrete Geometric Conservation Law (DGCL), which can be used to rigorously determine when the geometric parameters that appear in the fluxes should be computed.

The Geometric Conservation Law

We need to make sure that the movement of the mesh is not responsible for any artificial alteration of the physical phenomena involved, or at least, to make our best from a numerical point of view for the mesh movement to introduce an error of the same order as the one introduced by the numerical scheme. If System (2.4) is written for a constant state, assuming F ext = 0, we get, for any arbitrary closed volume C = C(t):

d (|C(t)|) dt - ∂C(t) (w • n) ds = 0 .
(2.12)

As the constant state is a solution of the Euler equations, if boundaries transmit the flux towards the outside as it comes, we find a purely geometrical relation inherent to the continuous problem.

For any arbitrary closed volume C = C(t) of boundary ∂C(t), Relation (2.12) is integrated into:

|C(t + δt)| -|C(t)| = t+δt t ∂C(t) (w • n) dsdt, with t and t + δt ∈ [0, T ] , (2.13)
which is usually known as the Geometric Conservation Law (GCL). From a geometric point of view, this relation states that the algebraic variation of the volume of C between two instants equals the algebraic volume swept by its boundary.

The role of the GCL in ALE simulations has been analyzed in [START_REF][END_REF]]. It has been

shown that the GCL is not a necessary condition to preserve time accuracy. However, violating it can lead to numerical oscillation [Mavriplis 2006]. In [START_REF] Farhat | [END_REF]] the authors show that compliance with the GCL guarantees an accuracy of at least the first order in some conditions.

Therefore, most would agree that the GCL should be enforced at the discrete level for a large majority of cases.

Discrete GCL enforcement

A new approach to enforcing the Discrete GCL was proposed in [Mavriplis 2006, Yang 2005[START_REF][END_REF]], in which the authors proposed a framework to build ALE high order temporal schemes that reach approximately the design order of accuracy. The originality of this approach consists in precisely defining which ALE parameters are true degrees of freedom and which are not. In contrast to other approaches [Koobus 1999, Lesoinne 1996, Nkonga 2000], they consider that the times and configurations at which the fluxes are evaluated do not constitute a new degree of freedom to be set thanks to the ALE scheme. To maintain the design accuracy of the fixed-mesh temporal integration, the moment at which the geometric parameters, such as the cells' interfaces' normals or the upwind/downwind tetrahedra must be computed, is entirely determined by the intermediate configurations involved in the chosen temporal scheme. The only degree of freedom to be set by enforcing the GCL at the discrete level is σ. Incidentally, it is implicitly stated that w is never involved alone but only hidden in the term σ n which represents the instantaneous algebraic volume swept.

In practical terms, the interfaces normal speeds are found by simply rewriting the scheme for a constant discrete solution, which leads to a small linear system that is easily invertible by hand.

This procedure is detailed in the next sections for one explicit Runge-Kutta scheme and several implicit schemes. Any fixed-mesh time-integration scheme can be extended to the case of moving meshes thanks to this methodology, and the resulting temporal scheme is naturally DGCL. Even if this has not been proven theoretically, the expected temporal order of convergence has also been observed numerically for several schemes designed using this method [Yang 2005].

Explicit Runge-Kutta schemes

The system of partial differential equations -given by Relation (2.6) -can be written in semidiscrete form for a vertex P i as:

d (|C i (t)|W i (t)) dt = F(W i (t)) ,
where F (W i (t)) is the residual, i.e., the discretization of the convective operator using an approximate Riemann solver and of the boundary conditions, and |C i | is the area/volume of the dual finite volume cell associated with P i . For a first-order explicit time discretization, the semi-discretized system reads:

|C n+1 i |W n+1 i -|C n i |W n i δt n i = F(W n i ) ,
where W n is the state at iteration n and δt n i is the local time step at iteration n. For high-order scheme in time Runge-Kutta scheme are considered.

Runge-Kutta (RK) methods are famous multi-stage methods to integrate ODEs. In the numerical solution of hyperbolic PDEs, notably the Euler equations, the favorite schemes among the huge family of Runge-Kutta schemes are those satisfying the Strong Stability Preserving (SSP) property [Shu 1988, Spiteri 2002]. In what follows, we denote by SSPRK(S,P ) the Sstage RK scheme of order P . We adopt the following notations:

F(W s i ) = n i j=1 |∂C s ij |Φ ij (W s i , W s j , n s ij , σ s ij ) = n i j=1 Φ ij (W s i , W s j , η s ij , σ s ij )
where

• n i is the number of vertices P j (edges e ij ) in the ball of vertex P i , i.e., we have P j ∈ V(P i ),

• η ij is the outward non-normalized normal to the portion of the interface of cell C s i around edge e ij , i.e., ∂C ij , and we have:

η ij = |∂C ij |n ij ,
• σ ij is the normal speed of the interface around edge e ij of cell C s i .

Superscript notation X s indicates that the quantity considered is the X obtained at stage s of the Runge-Kutta process. For instance, C s i is the cell associated with vertex P i when the mesh has been moved to its s th Runge-Kutta configuration.

In the following, coefficients (c s ) 0≤s≤S indicate the relative position in time of the current Runge-Kutta configuration: t s = t n + c s δt n with δt n = t n+1t n . Finally, we denote by A s ij the volume swept by the interface around edge e ij of cell C i between the initial Runge-Kutta configuration and the s th configuration.

Application to the SSPRK(4,3) scheme. This approach was used, for example, to build the 3rd order 4-step Runge-Kutta scheme [Olivier 2011b], whose Butcher and Shu-Osher representations are given in Table 2.1. For this scheme to be DGCL, it must preserve a constant solution

Butcher representation

Shu-Osher representation

t s = t n + c s δt n Y i 0 = Y i n Y i 0 = Y i n c 0 = 0, t 0 = t n Y i 1 = Y i 0 + δt n 2 F 0 i Y i 1 = Y i 0 + δt n 2 F 0 i c 1 = 1 2 , t 1 = t n + δt n Y i 2 = Y i 0 + δt n 2 F 0 i + F 1 i Y i 2 = Y i 1 + δt n 2 F 1 i c 2 = 1, t 2 = t n+1 Y i 3 = Y i 0 + δt n 6 F 0 i + F 1 i + F 2 i Y i 3 = 2 3 Y i 0 + 1 3 Y i 2 + δt n 6 F 2 i c 3 = 1 2 , t 3 = t n + δt n Y i 4 = Y i 0 + δt n 2 1 3 F 0 i + 1 3 F 1 i + 1 3 F 2 i + F 3 i Y i 4 = Y i 3 + δt n 2 F 3 i c 4 = 1, t 4 = t n+1
Table 2.1: Butcher and Shu-Osher representations of the 3rd order 4-step Runge-Kutta scheme (SSPRK(4,3).

W i = W 0 , as stated above. In this specific case, our conservative variable is Y i = |C i |W 0 and the purely physical fluxes vanish, leading to

F(W s i ) = -W 0 n i j=1 η s ij σ s ij .
Therefore, the scheme reads:

|C 0 i | = |C n i | |C 1 i | -|C 0 i | = n i j=1 A 1 ij = δt n 2 n i j=1 η 0 ij σ 0 ij |C 2 i | -|C 0 i | = n i j=1 A 2 ij = δt n 2 n i j=1 η 0 ij σ 0 ij + η 1 ij σ 1 ij |C 3 i | -|C 0 i | = n i j=1 A 3 ij = δt n 6 n i j=1 η 0 ij σ 0 ij + η 1 ij σ 1 ij + η 2 ij σ 2 ij |C 4 i | -|C 0 i | = n i j=1 A 4 ij = δt n 2 n i j=1 1 3 η 0 ij σ 0 ij + 1 3 η 1 ij σ 1 ij + 1 3 η 2 ij σ 2 ij + η 3 ij σ 3 ij ,
where A s ij is the volume swept by the interface around edge e ij of cell C i between the initial and the s th Runge-Kutta configuration. A natural necessary condition for the above relations to be satisfied is to have, for each interface around edge e ij of each finite volume cell C i :

        A 1 ij A 2 ij A 3 ij A 4 ij         = δt n         1 2 0 0 0 1 2 1 2 0 0 1 6 1 6 1 6 0 1 6 1 6 1 6 1 2                 η 0 ij σ 0 ij η 1 ij σ 1 ij η 2 ij σ 2 ij η 3 ij σ 3 ij         ⇐⇒         η 0 ij σ 0 ij η 1 ij σ 1 ij η 2 ij σ 2 ij η 3 ij σ 3 ij         = 1 δt n         2 0 0 0 -2 2 0 0 0 -2 6 0 0 0 -2 2                 A 1 ij A 2 ij A 3 ij A 4 ij         .
Therefore, the normal speed of the interface around edge e ij of cell C i must be updated in the Runge-Kutta process as follows :

σ 0 ij = 2A 1 ij δt n η 0 ij , σ 1 ij = -2A 1 ij + 2A 2 ij δt n η 1 ij , σ 2 ij = -2A 2 ij + 6A 3 ij δt n η 2 ij , and
σ 3 ij = -2A 3 ij + 2A 4 ij δt n η 3 ij , (2.14)
and the η s ij and A s ij are computed on the mesh once it has been moved to the s th Runge-Kutta configuration.

Implicit time discretization

For numerical scheme stability, explicit time discretization have strong constraint on the allowable time-step defined by the CFL condition. Such explicit schemes are widely used for blast problems (violent flows), LES, aeroacoustic, ... to control efficiently dissipation and dispersion errors that can spoil the overall solution accuracy. However, for other applications, implicit scheme are of main interest to reduce the time-step restriction. But, implicit schemes are more difficult to implement as they require the solution of a linear system and their cost per iterations is generally very large w.r.t explicit schemes.

As above, the system of partial differential equations can be written in semi-discrete form for a vertex P i as:

d (|C i (t)|W i (t)) dt = F(W i (t)) , (2.15)
where F(W i (t)) is the residual and |C i (t)| is the area/volume of the dual finite volume cell associated with P i at time t. To introduce the notations, we first assume we are not in the ALE framework and the finite volume cells have no dependency in time:

|C i | d (W i (t)) dt = F(W i (t)) .
Then, for a first order implicit time discretization (BDF1, i.e., backward differentiation formula),

we have :

|C i | δt n i (|W n+1 i -W n i ) = F(W n+1 i )
which is linearized as:

|C i | δt n i I d - ∂F ∂W (W n i ) W n+1 i -W n i = F(W n i )
where

∂F ∂W (W n i )
is the Jacobian contributing to the i th line of the matrix. We then rewrite the linearized system in compact form for all vertices in the mesh as:

A n δW n = R n where A n = |C| δt n I - ∂R n ∂W and δW n = W n+1 -W n .
In this section, we present three implicit integration schemes in the ALE context: BDF1, BDF2, Crank-Nicolson, and the Defect Correction methods (DeC) to increase the spatial accuracy. Again, we follow the Mavriplis and Yang approach [Mavriplis 2006] and use the same notations as for the explicit case:

• n i is the number of vertices P j (edges e ij ) in the ball of vertex P i , i.e., we have P j ∈ V(P i ),

• η ij is the outward non-normalized normal to the portion of the interface of cell C s i around edge e ij , i.e., ∂C ij , and we have:

η ij = |∂C ij |n ij and η ij = |∂C ij | is the area of the portion of the cell interface, • n ij = η ij η ij the normalized cell interface normal
• ω ij the cell interface velocity vector

• σ ij = ω ij • n ij is the normal velocity of the cell interface around edge e ij of cell C s i , thus η ij σ ij = ω ij • η ij is the interface normal integrated velocity.
For each implicit time integration scheme, we give the normal velocity for each cell interface (as its computation depends on the considered scheme) and the linear system to be solved. ALE BDF1 scheme. The first-order backward differentiation formula (BDF1) scheme reads:

|C n+1 |W n+1 -|C n |W n δt n = R(W n+1 , t n+1 ) , (2.16)
the mesh coordinates and velocities in the residual should only be evaluated at the new time step t n+1 . Thus, we have:

R(W n+1 , t n+1 ) = R(W n+1 , x n+1 , η n+1 σ n+1 ) .
As the mesh coordinates and the cell-face data should be taken at t n+1 , the only remaining degree of freedom is the mesh velocity ω n+1 (or σ n+1 ) and it is necessary to devise an implicit scheme which obeys the GCL. In fact, verifying the GCL defines completely the mesh velocity.

Computing σ. To compute the cells interface normal velocity, the GCL states that a constant state is preserved, thus for the BDF1 scheme the GCL reads:

|C n+1 i | -|C n i | = δt n n i j=1 η n+1 ij σ n+1 ij ,
The left-handside, which can be computed exactly (following [Nkonga 1993]),but it may also be written as a linear combination of the incremental changes in volume between the two time levels, which themselves may be decomposed into the elemental volumes swept by each control volume boundary face between the two time levels. The equality is then required to hold for each moving control volume boundary face ∂C ij , leading to:

σ n+1 ij = |∂C n+1 ij | -|∂C n ij | η n+1 ij δt n .
(2.17)

BDF1 linear system. As regards the BDF1 linear system, we linearize the residual of System (2.16) with respect to W n and let x and σ fixed :

|C n+1 |W n+1 -|C n |W n = δt n R(W n+1 , x n+1 , σ n+1 ) , ⇐⇒ |C n+1 |δW n + (|C n+1 | -|C n |)W n = δt n R(W n , x n+1 , σ n+1 ) + δt n ∂R ∂W (W n , x n+1 , σ n+1 ) δW n ,
to obtain the following BDF1 linear system

|C n+1 | δt n I - ∂R ∂W (W n , x n+1 , σ n+1 ) δW n = R(W n , x n+1 , σ n+1 ) - |C n+1 | -|C n | δt n W n . (2.18)
We notice that the Jacobian (in fact the whole matrix) and the residual are evaluated on the mesh at t n+1 but we use the solution at t n .

In the unsteady case, the implicit scheme given by Equation (2.18) is A-stable but it is only first order accurate in time and in space. To increase the time accuracy to second order, two approaches are possible: the use of second-order backward difference formulae (BDF2) or the Crank-Nicolson method. To increase the space accuracy to second order a defect-correction (DeC) method is considered. ALE BDF2 scheme. The second-order backward differentiation formula (BDF2) scheme is L-stable and uses three-points to integrate in time:

λ n+1 |C n+1 |W n+1 + λ n |C n |W n + λ n-1 |C n-1 |W n-1 δt n = R(W n+1 , t n+1 ) , (2.19)
where

λ n+1 = 3 2 , λ n = -2 , λ n-1 = 1 2 ,
if the time-step δt is constant. If the time-step changes at each iteration, then we have [Denner 2014, Koobus 1999]:

λ n+1 = 1 + 2r n 1 + r n , λ n = -1 -r n , λ n-1 = (r n ) 2 1 + r n with r n = δt n δt n-1 .
The mesh coordinates and velocities in the residual should only be evaluated at the new time step t n+1 . Thus, we have:

R(W n+1 , t n+1 ) = R(W n+1 , x n+1 , η n+1 σ n+1 ) .
As the mesh coordinates and the cell-face data should be taken at t n+1 , the only remaining degree of freedom is the mesh velocity ω n+1 (or σ n+1 ) and it is necessary to devise an implicit scheme which obeys the GCL. In fact, verifying the GCL defines completely the mesh velocity.

Computing σ. The GCL states that a constant state is preserved, thus for the BDF2 scheme the GCL reads:

λ n+1 |C n+1 | + λ n |C n | + λ n-1 |C n-1 | = δt n n i j=1 η n+1 ij σ n+1 ij ,
The left-handside, which can be computed exactly (following [Nkonga 1993]), may also be written as a linear combination of the incremental changes in volume between the various time levels, which themselves may be decomposed into the elemental volumes swept by each control volume boundary face between time levels. The equality is then required to hold for each moving control volume boundary face, leading to:

δt n η n+1 ij σ n+1 ij = γ n |∂C n+1 ij | -|∂C n ij | + γ n-1 |∂C n ij | -|∂C n-1 ij | =⇒ λ n+1 |C n+1 | + λ n |C n | + λ n-1 |C n-1 | = ni j=1 γ n |∂C n+1 ij | -|∂C n ij | + γ n-1 |∂C n ij | -|∂C n-1 ij | ,
leading to γ n = λ n+1 and γ n-1 = -λ n-1 . And, we verify that:

λ n = γ n-1 -γ n = -λ n+1 -λ n-1 = -1 -2r n -(r n ) 2 1 + r n = -(1 + r n ) = λ n .
Consequently, the face normal integrated velocity is defined by the relation:

η n+1 f σ n+1 ij = 1 δt n λ n+1 |∂C n+1 ij | -|∂C n ij | -λ n-1 |∂C n ij | -|∂C n-1 ij | ,
and thus the cell face normal velocity is

σ n+1 ij = λ n+1 |∂C n+1 ij | -|∂C n ij | -λ n-1 |∂C n ij | -|∂C n-1 ij | η n+1 ij δt n .
(2.20)

BDF2 linear system. To simplify the notation, the residual is now written R(W n+1 , x n+1 , σ n+1 ) meaning that it applies to the state at t n+1 , uses the mesh at t n+1 to compute all geometric quantities and considers σ n+1 defined above. The BDF2 scheme

λ n+1 |C n+1 |W n+1 + λ n |C n |W n + λ n-1 |C n-1 |W n-1 = δt n R(W n+1 , x n+1 , σ n+1 )
is linearized as follows:

λ n+1 |C n+1 |δW n + (λ n+1 |C n+1 | + λ n |C n |)W n + λ n-1 |C n-1 |W n-1 = δt n R(W n , x n+1 , σ n+1 ) + δt n ∂R ∂W (W n , x n+1 , σ n+1 ) δW n ,
to obtain the BDF2 linear system:

λ n+1 |C n+1 | δt n I - ∂R ∂W (W n , x n+1 , σ n+1 ) δW n = R(W n , x n+1 , σ n+1 ) (2.21) - 1 δt n (λ n+1 |C n+1 | + λ n |C n |)W n + λ n-1 |C n-1 |W n-1 .
We notice that the Jacobian (in fact the whole matrix) and the residual are evaluated on the mesh at t n+1 but we use the solution at t n . ALE Crank-Nicolson method. The Crank-Nicolson method (CN2) consists in using a general trapezoidal scheme:

|C n+1 |W n+1 -|C n |W n δt n = βR(W n+1 , t n+1 ) + (1 -β)R(W n , t n ) , (2.22)
with β = 1 2 (for β = 0 we recover the explicit scheme and for β = 1 we recover the Euler implicit time integration). This implicit scheme is A-stable and second order accurate in time.

The mesh coordinates and velocities in the residual should be evaluated at the current time step t n and the new time step t n+1 .

Computing σ. The GCL states that a constant state is preserved, thus for the CN2 scheme the GCL reads:

|C n+1 | -|C n | = δt n 2 n i j=1 η n+1 ij σ n+1 ij + η n ij σ n ij
Unfortunately here, we have two unknowns σ n ij and σ n+1 ij for one equation. We make the following choice to solve this issue (as a choice has to be made):

|C n+1 | -|C n | = δt n n i j=1 η n ij σ n ij and |C n+1 | -|C n | = δt n n i j=1 η n+1 ij σ n+1 ij
Note that σ n ij and σ n+1 ij will not be identical because the cell interface normal varies between t n and t n+1 . Following, the results of the BDF1 scheme, we get:

σ n ij = |∂C n+1 ij | -|∂C n ij | η n ij δt n and σ n+1 ij = |∂C n+1 ij | -|∂C n ij | η n+1 ij δt n .
(2.23)

Crank-Nicolson linear system. The CN2 scheme:

|C n+1 |W n+1 -|C n |W n = δt n 2 R(W n+1 , x n+1 , σ n+1 ) + δt n 2 R(W n , x n , σ n ) ,
is linearized as follows:

|C n+1 |δW n +(|C n+1 |-|C n |)W n = δt n 2 R(W n , x n+1 , σ n+1 )+ δt n 2 ∂R ∂W (W n , x n+1 , σ n+1 ) δW n + δt n 2 R(W n , x n , σ n ) ,
to obtain our linear system

|C n+1 | δt n I - 1 2 ∂R ∂W (W n , x n+1 , σ n+1 ) δW n = 1 2 R(W n , x n+1 , σ n+1 ) + 1 2 R(W n , x n , σ n ) - |C n+1 | -|C n | δt n W n .
(2.24)

We note that in comparison to the fixed mesh case, the CN2 scheme require two residual evaluations, one on the mesh at t n and the other one on the mesh at t n+1 .

Defect correction method. The BDF2 and the CN2 time integration schemes provide second order accurate schemes in time, but only first order accurate scheme in space. To recover the second order space convergence of the numerical scheme, the linearized linear system is solved by an iterative method proposed in [Martin 1996] and called Defect Correction (DeC). In the case of fixed meshes, the solution W n+1 is obtained after α max DeC iterations:

Set W n+1,1 = W n For α = 1, . . . , α max solve |C| δt n I - ∂R n ∂W (W n+1,α+1 -W n+1,α ) = R n - |C| δt n I - ∂R n ∂W (W n+1,α -W n ) = R n+1,α - |C| δt n W n+1,α -W n EndFor Set W n+1 = W n+1,αmax
To recover the second order in time and in space, it is sufficient to couple the above methods. Now, we consider the case with dynamic meshes and we write the coupling of the DeC scheme with the BDF2 and CN2 schemes.

Coupling the DeC method with BDF2 leads to solve the following sub-iteration linear system:

λ n+1 |C n+1 | δt n I - ∂R ∂W (W n , x n+1 , σ n+1 ) W n+1,α+1 -W n+1,α = R(W n+1,α , x n+1 , σ n+1 ) -λ n+1 |C n+1 | δt n W n+1,α -W n - 1 δt n (λ n+1 |C n+1 | + λ n |C n |)W n + λ n-1 |C n-1 |W n-1 .
Coupling the DeC method with Crank-Nicolson leads to solve the following sub-iteration linear system:

|C n+1 | δt n I - 1 2 ∂R ∂W (W n , x n+1 , σ n+1 ) (W n+1,α+1 -W n+1,α ) = 1 2 R(W n , x n , σ n ) + 1 2 R(W n , x n+1 , σ n+1 ) - |C n+1 | -|C n | δt n W n - |C n+1 | δt n I - 1 2 ∂R n ∂W (W n , x n+1 , σ n+1 ) (W n+1,α -W n ) = 1 2 R(W n , x n , σ n ) + 1 2 R(W n+1,α , x n+1 , σ n+1 ) - |C n+1 | -|C n | δt n W n - |C n+1 | δt n (W n+1,α -W n ) = 1 2 R(W n , x n , σ n ) + 1 2 R(W n+1,α , x n+1 , σ n+1 ) - 1 δt n |C n+1 |W n+1,α -|C n |W n .
Practical implementation. In practice, the following steps are done in the implementation:

1. We compute the volume swept by the cell faces using Nkonga formula, see next Section (with the cell-face velocity ω ij and the cell-face non-normalized pseudo-normal ηij ), 2. We compute the cell-face normal velocity σ ij using Equation (2.17) or (2.20) or (2.23) 3. Numerical fluxes can be evaluated using σ ij .

Practical computation of the volumes swept

The interface of a finite volume cell is made up of several triangles, connecting the middle of an edge to the center of gravity of a face and the center of gravity of that tetrahedron, see Figure 2.6. The two triangles of the interface sharing one edge within a tetrahedron are coplanar (i.e. the middle of an edge, the center of gravity of the two faces neighboring this edge and the center of gravity of tetrahedron are coplanar). The union of these two triangles is called the facet associated to the edge and the tetrahedron. At configuration t s = t n + c s δt n , (with t n+1 = t n + δt n ), the outward non-normalized pseudo normal ηs ij and the volume swept A s ij are computed as described in [START_REF] Nkonga | [END_REF]]. As the cell interface is made up of several facets, the total swept volume is the sum of the volumes swept by each facet. Let us assume that the facet considered is associated with edge e ij = P i P j and belongs to tetrahedron K = (P 0 , P 1 , P 2 , P 3 ). In what follows, i = j = l = m ∈ [[0, 3]], G denotes the center of gravity of tetrahedron K, M m denotes the gravity center of face F m = (P i , P j , P l ) of tetrahedron K and M l the center of gravity of face F l = (P i , P j , P m ). The outward nonnormalized pseudo normal of the facet is given by:

ηs ij,K = 1 4 G n M n m ∧ G s M s l + 1 4 G s M s m ∧ G n M n l + 1 2 G n M n m ∧ G n M n l + 1 2 G s M s m ∧ G s M s l , (2.25)
where

G n M n m = 1 12 P n i + P n j -3P n m + P n l , G s M s m = 1 12 P s i + P s j -3P s m + P s l , G n M n l = 1 12 P n i + P n j + P n m -3P n l , G s M s l = 1 12 P s i + P s j + P s m -3P s l .
The volume swept by the facet is:

A s ij,K = c s δt n (w G ) s ij,K • ηs ij,K , (2.26)
with the mean velocity written as:

(w G ) s ij,K = 1 36 c s δt n 13w s i + 13w s j + 5w s m + 5w s l with w s ξ = P s ξ -P n ξ .
Finally, the total volume swept by the interface around edge e ij of cell C i is obtained by summing over the shell of the edge, i.e. all the tetrahedra sharing the edge:

A s ij = K∈Shell(i,j) A s ij,K . (2.27)
It is important to understand that normals ηs ij,K are pseudo-normals which are used only to compute the volumes swept by the facets. They must not be mistaken for normals to facets taken at t s , η s ij,K , which are used for the computation of ALE fluxes.

Volumes swept by boundary interfaces

The pseudo-normals and swept volumes of boundary faces are computed in a similar way. Let K = (P 0 , P 1 , P 2 ) be a boundary triangle, as in Figure 2.6. Let M 0 , M 1 and M 2 be the middles of the edges and G the center of gravity of K. The triangle is made up of three quadrangle finite volume interfaces: (P 0 , M 2 , G , M 1 ) associated with cell C 0 , (P 1 , M 0 , G , M 2 ) with C 1 and (P 2 , M 1 , G , M 0 ) with C 2 . Each boundary quadrangle interface I i is made of two sub-triangles, noted T ij and T ik with j, k = i.

The volume swept by interface I i between the initial and current configurations is the sum of the volumes swept by its two sub-triangles:

A s i,K = A s T ij + A s T ik = c s δt n w s G ij • ηs ij + c s δt n w s G ik • ηs ik , (2.28)
where w s G ij is the velocity of the center of gravity G ij of triangle T ij and ηs ij is the pseudo-normal associated with T ij , computed between the initial and current configurations. The six triangles T ij are coplanar, so their pseudo-normals have the same direction. Moreover, as median cells are used, their pseudo-normals also have the same norm, which is equal to one sixth of the norm of the pseudo-normal to triangle K. This common pseudo-normal is therefore equal to:

ηs = 1 6 ηs K = 1 6 • 1 3 1 4 P n 0 P n 1 ∧ P s 0 P s 2 + 1 4 P s 0 P s 1 ∧ P n 0 P n 2 + 1 2 P n 0 P n 1 ∧ P n 0 P n 2 + 1 2 P s 0 P s 1 ∧ P s 0 P s 2 , (2.29) thus A s i,K = c s δt n w s G ij + w s G ik • ηs , (2.30)
where :

w s G ij = 1 18 c s δt n 11w s i + 5w s j + 2w s k with w s ξ = P s ξ -P n ξ .
Finally, the total volume swept by the boundary interface is:

A s i = K∈Ball(i)
A s i,K .

(2.31)

MUSCL approach and temporal schemes

Regarding spatial accuracy, we have seen that the order of accuracy can be enhanced using the MUSCL-type reconstruction with upwinding. However, in the ALE context, one must determine how and when upwind/downwind elements should be evaluated to compute upwind/downwind gradients which are necessary for the β-schemes. This question is neither answered in the literature and generally approximations are carried out. For instance, some papers propose to use upwind and downwind elements at t n for the whole Runge-Kutta process. However, this choice should be consistent with the considered time integration scheme. Following the framework of [Yang 2005], it is clear that preserving the expected order of accuracy in time imposes that the upwind/downwind elements and the gradients are computed on the current configuration,

i.e., on the mesh at t s for the Runge-Kutta stage or the mesh specified in the residual for implicit schemes. Therefore, similarly to geometric parameters, the upwind/downwind elements and the gradients should be re-evaluated at each step of the Runge-Kutta stage or at each time step for implicite schemes.

Computation of the time step

For explicit time integration : The maximal allowable time step for the numerical scheme is:

δt n (P i ) = h(P i ) c i + u i -w i (2.32)
where h(P i ) is the smallest height in the ball of vertex P i , c i is the speed of sound, u i is the Eulerian speed (the speed of the fluid, computed by the solver) and w i is the speed of the mesh vertex. The global time step is then given by δt n = CF L min P i (δt n (P i )).

For implicit time integration : Time step is up to the user.

Handling the swaps

An ALE formulation of the swap operator, satisfying the DGCL, was proposed in 2D

in [Olivier 2011b], but its extension to 3D is very delicate because it requires to handle 4D geometry, and it has not been carried out yet. Instead of considering an ALE scheme for the swap operator, our choice in this work is to perform the swaps between two solver iterations, i.e.

at a fixed time t n . This consequently means that during the swap phase, the mesh vertices do not move, and thus the swaps do not impact the ALE parameters η and w, unlike [Olivier 2011b] where the swaps are performed during the solver iteration, i.e. between t n and t n+1 . After each swap, the solution should be updated on the new configuration. Two interpolation methods are considered here.

The first, and simplest, one is to perform a linear interpolation to recover the solution. As only the connectivity changes and not the vertices positions, the solution at the vertices does not change, i.e., nothing has to be done. This interpolation is DGCL compliant, since the constant state is preserved (in fact, any linear state is preserved), but it does not conserve the mass (i.e., it does not conserve the integral of the conservative variable) which is problematic for conservative equations when discontinuities are involved in the flow.

The second method is the P 1 -exact conservative interpolation following [Alauzet 2010b, Alauzet 2016]. It is a simplified version of the latter because the cavity of the swap configuration is fixed. The mass conservation property of the interpolation operator is achieved by elementelement intersections. The idea is to find, for each element of the new configuration, its geometric intersection with all the elements of the previous configuration it overlaps and to mesh this geometric intersection with simplices. We are then able to use a Gauss quadrature formula to exactly compute the mass which has been locally transferred. High-order accuracy is obtained through the reconstruction of the gradient of the solution from the discrete data and the use of some Taylor formulae. Unfortunately, this high-order interpolation can lead to a loss of monotonicity. The maximum principle is recovered by correcting the interpolated solution in a conservative manner, using a limiter strategy very similar to the one used for Finite-Volume solvers. Finally, the solution values at vertices are reconstructed from this piecewise linear by element discontinuous representation of the solution. The algorithm is summarized in Algorithm 9 where m K stands for the integral of any conservative quantities (density, momentum and energy) on the considered element. This method is also compliant with the DGCL.

Algorithm 7 Conservative Interpolation Process 1. For all elements K back of the original cavity, compute solution mass m K back and gradient ∇ K back 2. For all elements K new of the new cavity, recover solution mass m Knew and gradient ∇ Knew :

(a) compute the intersection of K new with all K i back it overlaps (b) mesh the intersection polygon/polyhedra of each couple of elements (K new , K i back ) (c) compute m Knew and ∇ Knew using Gauss quadrature formulae =⇒ a piecewise linear discontinuous representation of the mass on the new cavity is obtained 3. Correct the gradient to enforce the maximum principle 4. Set the solution values to vertices by an averaging procedure Moreover, after each swap, the data of the finite volume cells (volume and interface normals) are updated, together with the topology of the mesh (edges and tetrahedra). This requires to have a flow solver with dynamic data.

FSI coupling

The moving boundaries can have an imposed motion, or be driven by fluid-structure interaction.

A simple solid mechanics solver is coupled to the flow solver described previously. The chosen approach is the 6-DOF (6 Degrees of Freedom) approach for rigid bodies.

In this work, the bodies are assumed to be rigid, of constant mass and homogeneous, i.e., their mass is uniformly distributed in their volume. The bodies we consider will never break into different parts. Each rigid body B is fully described by: Physical quantities: its boundary ∂B and its associated inward normal n, its mass m assumed to be constant, its d×d matrix of inertia J G computed at G which is symmetric and depends only on the shape and physical nature of the solid object4 .

Kinematic quantities: the position of its center of gravity x G = (x(t), y(t), z(t)), its angular displacement vector θ = θ(t) and its angular speed vector ω = dθ/dt.

F ext denotes the resultant vector of the external forces applied on B, M G (F ext ) the kinetic moment of the external forces applied on B computed at G and g the gravity vector. We assume that the bodies are only submitted to forces of gravity and fluid pressure. The equations for solid dynamics in an inertial frame then read:

       m d 2 x G dt 2 = F ext = ∂B p(s)n(s)ds + mg J G d 2 θ dt 2 = M G (F ext ) = ∂B (s -x G ) × p(s)n(s) ds .
(2.33)

The equation governing the position of the center of gravity of the body is easy to solve since it is linear. Its discretization is straightforward. However, the discretization of the second equation, which controls the orientation of the body, is more delicate. Since the matrix of inertia J G depends on θ, it is a non-linear second order ODE. The chosen discretization is extensively detailed in [Olivier 2011a] and is based on rewriting of the equations in the frame of the moving body.

As the geometry must be moved in accordance with the fluid computation, the same time integration scheme has been taken to integrate the fluid and the solid equations. Therefore, time-advancing of the rigid bodies ODE System is performed using the same RKSSP or implicit scheme as the one used to advance the fluid numerical solution. The coupling is loose and explicit as the external forces and moments acting on rigid objects are computed on the current configuration.

Numerical experiments

Naca0012 Pitching

The formation of vortices in the wake of an aircraft is essentially the result of the separation of flow and the vorticity created by the lifting surfaces. Air below the wing is at higher pressure than the air pressure above. This pressure difference causes air to flow from the lower surface of the wing, around the wingtip towards the upper surface of the wing. This phenomenon creates a vortex which is lift-generating. But wingtip vortices in lift generation are also associated to wake turbulence and increase the drag of the aircraft. This could be particularly dangerous during take-off, landing or when the aircraft operates a high-angle of attack. Vortical flows, if not properly understood and controlled, may be undesirable aerodynamic effects and can severely influence the flight stability or speed up the degradation of the aircraft structure. Consequently, there is a continuous need for better understanding of vortex flow characteristics to predict the flow behavior and to help the designers to develop more effective means.

Mesh

In a circular domain of radius 20 is put a NACA0012 of length 1. The entire mesh includes 3291 vertices.

Initial conditions

To generate similar conditions of take-off, the initial Mach number chosen is 0.4.

Moving mesh

The moving of the NACA0012 is imposed. In phase 1, the coordinates of the geometry do not moved. The flow around the aircraft is established. Let t 1 = 0.1 be the time of the beginning of the movement and θ max = 20 o the maximal angle reached at first time t θ max . The movement of the wing profile is defined in phase 2 by calculating at each time step its angle of rotation α = 0.5 * θ max * (1cos(π * (tt 1 )/(t θ maxt 0 ))). The coordinates of the new points are then

x new = x old * cos α + y old * sin α and y new = -x old * sin α + y old * cos α. A phase 3 follows where the NACA0012 stays in upward position before the descent computed in the same way as phase 2.

Numerical results and comparison between implicit and explicit solvers

The simulation is done both in explicit Euler SSPRK(2,5) and BDF2 DeC. The results are shown in Figures 2.7 and 2.8.

On the upper surface of the airfoil, counter-clockwise vortices are observed with both temporal schemes. The initial formation of this bubble can be observed and is moved on the airfoil upper surface towards the leading edge as angle of attack increases. Once the maximal angle reached, the flow separated at the leading edge of the airfoil reattaches at the half of the chord. The formation of the new counterclockwise bubble causes this reattachment just behind the midchord. The bubble length is approximately two times smaller than the first one length and the second separation is observed about the last quarter chord location. The trailing edge vortex in counter-clockwise direction now causes once more the flow to attach at the end of the chord.

This bubble is well captured in the BDF-2 scheme simulation but the vortex in the wake is more accurate in the SSPRK(2,5) scheme. We know that the SSPRK(2,5) and the BDF-2 are schemes of order 2. And with the numerical results, the two schemes for concluding have quite the same results in this 2D test case.

Explicit Implicit But the CFL for the SSPRK(2,5) explicit scheme cannot exceed 3which is the faster secondorder in time explicit scheme we have, when CFL 500 is used for the BDF-2 scheme. The CPU times of both simulations are summarized in the

FSI Case : 2D Bomb drop

Mesh

To compute this case, we first define a 2D mesh as a cross-section view of an initial 3D mesh represented in Figure 2.10. We first keep only the triangles involved in the intersection of the initial mesh surface and the plane z = 3.2. Then for each of these triangles, we compute some points at the intersection and necessarily the points at the extremities. As the wing profile looks like a NACA0012, we decide to use it instead of my section view and to keep the ratios considering the NACA0012 is 1-length. Finally we choose a circular domain of radius 20 and center (0.5, 0). The final 2D mesh (see Figure 2.11) has 25779 vertices. It is generated by AFLR [START_REF] Marcum | Efficient Generation of High-Quality Unstructured Surface and Volume Grids[END_REF]]. AFLR is an unstructured mesh generation software based on an advancing-front/local-reconnection method.

Initial conditions

-Aerodynamic conditions : Mach number = 0.8, Angle of Attack = 2 • -External forces To eject the bomb an ejection force is applied to provide it to move up due to the external forces and strike the aircraft. It's supposed to be F = (0.; -0.02469).

-Gravity vector Considering the aircraft is flying below 20 km above ground so the gravity vector stands for g = (0.; -9.74). -Bomb charcteristics Obviously, its gravity center must be located in the head of the bomb to fall. We consider it at (0.06, -0.21). To be realistic, its mass is 500 kgs and its matrix of inertia is defined as for a cylinder of radius 0.032 and height 0.6 i.e. :

     15.125 0 0 0 15.125 0 0 0 0.251     

Numerical results

The simulation is done with a BDF-2 scheme with DeC in FSI motion.

Observations and Interpretations

At initial time the plane is flying and carrying the bomb. At the moment the bomb is released, it has the same velocity as the aircraft. At time t = 0, an attached oblique shockwave is located The chosen case concerns vortex flow created by the nosing up/nosing down of a subsonic F-117 aircraft. As a first step the F-117 aircraft geometry nosing up, that creates a vortical wake.

3D Nosing-up F117

An inflow of air at Mach 0.4 arrives in front of the aircraft, initially in horizontal position, that noses up, stays up for a while, then noses down. In this example, the aircraft rotates around its center of gravity. Let T = 1 be the characteristic time of the movement and θ max = 20 o the maximal angle reached, the movement is defined by its angle of rotation, of which the evolution is divided in 7 phases:

θ(t) = θ max                                  0 if 0 ≤ t ≤ T /2 (i) 2(t-T 2 ) 2 T 2 if T 2 < t ≤ T (ii) 1 2 + 2 ( 2 (t-T 2 ) T -1) -1 2 ( 4 (t-T 2 ) 2 T 2 -1) if T < t ≤ 3T 2 (iii) 1 if 3T 2 < t ≤ 7T 2 (iv) 1 - 2 (t-7T 2 ) 2 T 2 if 7T 2 < t ≤ 4T (v) 1 2 -2 ( 2 (t-7T 2 ) T -1) + 1 2 ( 4 (t-7T 2 ) 2 T 2 -1) if 4T < t ≤ 9T 2 (vi) 0 if 9T 2 t ≤ 5T (vii)
Phase (i) is an initialization phase, during which the flow around the aircraft is established.

Phases (ii) and (iii) are respectively phases of accelerated and decelerated ascension. Vortices start to grow behind the aircraft, and they expand during phase (iv), where the aircraft stays in upward position. Phases (v) and (vi) are phases of accelerated and decelerated descent, the vortices start to move away and they slowly disappear in phase (vii). Free-stream conditions are imposed on the faces of the surrounding box, and slipping conditions on the aircraft. Table of CPU time for the three simulations 2.3. The first column stands for the explicit RK2 scheme computation, the second for the BDF1 scheme, the third for the BDF2 scheme :

Comparison between implicit and explicit solvers

SSPRK2 BDF1 BDF2
CPU Time 3h2min13s 0h39min38s 0h44min8s The difference of time between BDF1 and BDF2 can be explained by the additional step computed for the second order scheme.

Conclusion

In this chapter, we presented the 3D ALE solver for the compressible Euler equations that was implemented in two and three dimensions. As regards time discretization, the multi-step Runge-Kutta schemes were already implemented in [Barral 2015] and the contribution of this thesis is the implementation of implicit methods. Details on the computation of the geometric parameters are provided. A 6-DOF mechanical model for rigid bodies was considered to add FSI coupling to the fluid solver. This solver was implemented in 3D in the same code as the moving mesh algorithm, and was plugged into it as explained in Algorithm 5. We also present numerous cases using this ALE solver and the connectivity-change moving mesh method.

In this chapter, we have seen the two aspects of simulations with moving boundaries: the handling of the moving mesh, and its consequences on the solver side. This algorithm allows us to run large displacement complex moving mesh simulations without remeshing.

Following from [Barral 2015], a novelty was presented in this chapter :

• A new integration in time was introduced. Three implicit time integration schemes ALE solver were implemented : BDF1, BDF2, Crank Nicolson schemes and parallelized in 3D.

The formulas for the computation of the geometrical parameters were clarified.

• The two BDF implicit solver were validated on test cases in 2D and 3D.

• ALE simulations were run in 2D and 3D, with imposed motion or Fluid-Structure Interaction.

• These considered implicit solvers are efficient to speed-up low-mach simulations, in particular, in the presented test cases, the subsonic aeronautic cases.

Perspectives arising from this work concern both the moving mesh aspect and the solver side.

• A strategy for moving boundary layer meshes around deformable geometries naturally arises from what is presented. It needs to be developed and applied on various cases.

• More generally, handling boundary layer meshes would make it possible to run Navier-Stokes simulations.

Now, to improve the accuracy of these moving mesh simulations without increasing unreasonably their cost, we use metric-based mesh adaptation. We will see it in the next chapter.

Chapter 3

Hessian-based mesh adaptation for moving geometries

Simulating complex moving geometries evolving in unsteady flows in three dimensions is more and more required by industry. But it still remains a challenge. The first reason of this challenge is the modelling as we already seen but this is not the only explanation. Simulating moving geometries is also very costful. In fact preserving a good accuracy is often synonym of fine mesh and great number of vertices. To increase the accuracy of these simulations while preserving the same number of vertices and a similar CPU time, anisotropic metric-based mesh adaptation is one solution. Generating an adapted mesh means optimizing the accuracy of some parts of the simulation domain while de-emphazing other parts we don't need precision. Mesh adaptation allows computation of a sufficiently precise solution without enormous requirements for memory and CPU time. Anisotropic mesh adaptation is one of these mesh adaptations and is based on the control of the interpolation error of a piecewise polynomial interpolation of the exact solution on the triangular grid. The mesh adaptation criterion checks the behaviour of the second order derivatives of the solution of the considered problem. Even if in this thesis we only use it with the solver previously developed in Chapter 2, we can mention that the additional advantage of it is that it can be used without any modification for arbitrary boundary value problem and arbitrary numerical method (finite element method, finite volume method, discontinuous Galerkin method, etc.). Anisotropic mesh adaptation has already proved its efficiency for steady problems, and appears as a good perspective for unsteady problems. However, its extension to the unsteady case involving moving geometries is absolutely not straightforward. These simulations combine the difficulties arising from unsteadiness and geometrical complexity :

global time step driven by the mesh smallest altitude, evolution of the phenomena in the whole domain, interpolation spoiling, but also three-dimensional meshing and remeshing issues with an imposed discretized surface. In fact, the introduction of moving geometries in this process raises new difficulties : handling of the mesh movement without deterioration of its quality, handling the specific numerical schemes imposed by moving mesh schemes and their restrictions, as well as fluid/structure coupling.

Three different approaches dealing with time-dependent mesh adaptation in the literature can be distinguished : re-meshing methods, fixed-point algorithms and adaptative moving meshes.

First, an isotropic mesh is adapted frequently in order to maintain the solution within refined regions and a safety area around critical regions is introduced [Löhner 1990, Löhner 1992, Rausch 1992, Speares 1997]. Another approach is to use an unsteady mesh adaptation algorithm [de Sampaio 1993[START_REF] Wu | [END_REF]] based on the estimation of the error every n flow solver iterations and global remeshing techniques. If the error is greater than a prescribed threshold, the mesh is re-adapted. A similar strategy is adopted in [Picasso 2003, Picasso 2009],

where, at each re-adaptation step, a nearly optimal metric is found by scaling the metrics in every directions depending on directional error factors. More recently, local adaptive remeshing enabling the construction of anisotropic meshes has been considered. In this case [START_REF] Pain | Tetrahedral mesh optimisation and adaptivity for steady-state and transient finite element calculations[END_REF], Remacle 2005, Compère 2008], the mesh is frequently adapted in order to guarantee that the solution always evolves in refined regions. All these approaches involve a large number of adaptations while introducing unquantified errors due to the transfer of the solution from the old mesh to the new one. Moreover, none of them considers the inherent non-linear nature of the mesh adaptation problem: the convergence of the mesh adaptation process is never addressed and therefore obtaining the optimal mesh cannot be expected.

A first answer to these issues has already been proposed in [Alauzet 2007]. It relies on the assumption that the temporal error is always controlled by the spatial error, which is indeed the case when solving a linear advection problem under a CFL condition. A fixedpoint algorithm is used to predict the solution and converge the couple mesh/solution. The simulation interval is split into fixed-size sub-intervals on which an unique adapted mesh is used. More recently, mesh adaptation and adaptive time-stepping methods were combined for incompressible unsteady simulations [Coupez 2013]. In that approach, moving boundaries are handled via level-set methods.

In adaptive moving meshes, also called r-adaptation, the mesh is moved at each time step into a mesh adapted to the solution via a mesh equation. A large range of methods share this common principle, the differences between each lying in the adaptation criterion, the way it influences the mesh movement and the way the mesh movement is taken into account reciprocally in the physical equations. In general, they strongly couple the mesh movement and the physical equation being solved, resulting in methods with no spoiling interpolation errors. Some methods are velocity-based, meaning the mesh position is found through the vertices velocity. They include Lagrangian methods, where the velocity of the mesh is the actual Lagrangian velocity of the fluid [Loubère 2010], Moving Finite Elements, where the mesh celerity is the solution of a minimization problem [START_REF] Miller | [END_REF], Baines 1994],

and methods based on the conservation of the Geometric Conservation Law [START_REF] Cao | [END_REF]].

These methods tend to quickly produce tangled meshes. In location-based methods, the mesh equation directly gives the position of the vertices. The adapted meshes are seen as images of a reference mesh by a coordinate transformation that is the solution of a minimization problem. This problem usually involves a monitor function to specify the sizes of the meshes, which is either heuristic or derived from error estimates. The problems can be based on Laplace [Winslow 1963] or Poisson [START_REF] Thompson | [END_REF]] equations, or harmonic maps [Dvinsky 1991]. The application of the equidistribution principle leads to so called Moving-Mesh PDEs (MMPDEs) [Huang 2010b]. More recently, Monge-Ampère equation was also used to determine the mesh movement [Chacón 2011[START_REF] Browne | [END_REF]]. These approaches seem very time-consuming for now, since the PDEs are solved at every solver time step.

Moreover, the solution of the Monge-Ampère equation is known to be very difficult in three dimensions.

This section is the continuation of the work of [Olivier 2011a], [Barral 2015], [N. [START_REF][END_REF].

In those works, the authors proposed a new error estimation for simulations with moving geometries and a new ALE metric. But, mesh optimizations were not adapted consistently during the process. The contribution made in this thesis is the maintain of a consistent ALE metric in time-accurate anisotropic mesh adaptation for time-dependent problems with moving geometries.

This chapter is divided into two axes. The first sections detail the theory of error estimation and optimal metric for unsteady simulations and the numerical implementation in practice.

The last sections focus on the update for dynamic meshes adaptation and the numerical choices made.

Error estimate for unsteady simulations

Let us first detail the error analysis that leads to the adaptation algorithm for unsteady simulations. This analysis is performed in the context of the continuous mesh model: first, an error model is proposed, then this error is minimized to derive an optimal mesh.

In the context of time-dependent problems, the error analysis for the steady case recalled in Section 1.2.1 is not sufficient, since it controls only spatial errors, whereas temporal errors need to be addressed too. In this work, we do not account for time discretization errors but we focus on a space-time analysis of the spatial error. In other words, we seek for the optimal space-time mesh controlling the space-time spatial discretization error. This assumption makes sense for the type of simulations that are considered in this work considering that as an explicit time scheme is used for time advancing, then the error in time is controlled by the error in space under CFL condition. This has been demonstrated under specific conditions in [Alauzet 2007].

As long as this hypothesis holds, the spatial interpolation error provides a fair measure of the total space-time error of the discretized unsteady system. But it will be no longer valid in the case of future implicit time advancing solvers. The analysis would have to be completed to take into account the temporal discretization error.

Error model

Our goal is to solve an unsteady PDE which is set in the computational space-time domain

Q = Ω × [0, T ]
where T is the (positive) maximal time and Ω ⊂ R 3 is the spatial domain. Let Π h be the usual P 1 projector, we extend it to time-dependent functions:

(Π h ϕ) (t) = Π h (ϕ(t)) , ∀ t ∈ [0, T ] . (3.1)
The considered problem of mesh adaptation consists in finding the space-time mesh H of Ω that minimizes the space-time linear interpolation error u -Π h u in L p norm. The problem is thus stated in an a priori way:

Find H opt having N st vertices such that E L p (H opt ) = min H u -Π h u L p (Ω h ×[0,T ]) . (3.2)
This problem is ill-posed and has far too many unknowns to be solved directly, as was explained previously. So it is rewritten in the continuous mesh framework under its continuous form:

Find M L p = (M L p (x, t)) (x,t)∈Q such that E L p (M L p ) = min M u -π M u L p (Ω×[0,T ]) , (3.3)
under the space-time constraint:

C st (M) = T 0 τ (t) -1 Ω d M (x, t) dx dt = N st . (3.4)
where τ (t) is the time step used at time t of interval [0, T ]. Introducing the continuous interpolation error, we recall that we can write the continuous error model as follows:

E L p (M) = T 0 Ω trace M -1 2 (x, t)|H u (x, t)|M -1 2 (x, t) p dx dt 1 p . (3.5)
where H u is the Hessian of sensor u. To find the optimal space-time continuous mesh, ) is solved in two steps:

(i) First, a spatial minimization is done for a fixed t.

(ii) Second, a temporal minimization is performed.

Note that both minimizations are performed formally.

Spatial minimization for a fixed t

Let us assume that at time t, we seek for the optimal continuous mesh M L p (t) which minimizes the instantaneous error, i.e., the spatial error for a fixed time t:

ẼL p (M(t)) = Ω trace M -1 2 (x, t) |H u (x, t)| M -1 2 (x, t) p dx
under the constraint that the number of vertices is prescribed to

C(M(t)) = Ω d M(t) (x, t) dx = N (t). (3.6)
Solving the optimality conditions provides the optimal instantaneous continuous mesh in L pnorm M L p (t) = (M L p (x, t)) x∈Ω at time t defined by:

M L p (x, t) = N (t) 2 3 M L p ,1 (x, t) , (3.7)
where M L p ,1 is the optimum for C(M(t)) = 1:

M L p ,1 (x, t) = Ω (det |H u (x, t)|) p 2p+3 dx -2 3 (det |H u (x, t)|) -1 2p+3 |H u (x, t)|. (3.8)
The corresponding optimal instantaneous error at time t writes:

ẼL p (M L p (t)) = 3 p N (t) -2p 3 Ω (det |H u (x, t)|) p 2p+3 dx 2p+3 3 = 3 p N (t) -2p 3 K(t) . (3.9)
Throughout this thesis we denote:

K(t) = Ω (det |H u (x, t)|) p 2p+3 dx 2p+3 3
.

Temporal minimization

To complete the resolution of optimization Problem (3.3-3.4), a temporal minimization is performed in order to get the optimal space-time continuous mesh. In other words, we need to find the optimal time law t → N (t) for the instantaneous mesh size. In this work, we consider the case where the time step τ is specified by the user as a function of time t → τ (t). The analysis can be extended to the case of an explicit time advancing solver subject to Courant time step condition, but the resulting optimal mesh is too complex to be used in practice.

Let the time step τ be specified by a function of time t → τ (t). After the spatial optimization, the space-time error writes:

E L p (M L p ) = T 0 ẼL p (M L p (t)) dt 1 p = 3 T 0 N (t) -2p 3 K(t) dt 1 p (3.10)
and we aim at minimizing it under the following space-time complexity constraint:

T 0 τ (t) -1 N (t) dt = N st . (3.11)
We want to find the optimal distribution of N (t) when the space-time total number of nodes N st is prescribed. We can apply the following one-to-one change of variables:

Ñ (t) = N (t)τ (t) -1 and K(t) = τ (t) -2p 3 K(t) .
Then, the temporal optimization problem becomes:

min M E L p (M) p = 3 p T 0 Ñ (t) -2p 3 K(t) dt under constraint T 0 Ñ (t) dt = N st .
The solution of this problem is given by:

Ñopt (t) -2p+3 3 K(t) = const ⇒ N opt (t) = C(N st ) (τ (t) K(t)) 3 2p+3 .
Here, constant C(N st ) can be obtained by introducing the above expression in space-time complexity Constraint (3.11), leading to:

C(N st ) = T 0 τ (t) -2p 2p+3 K(t) 3 2p+3 dt -1 N st ,
which completes the description of the optimal space-time metric for a prescribed time step. Using Relations (3.7) and (3.8), the analytic expression of the optimal space-time metric in L p -norm M L p writes:

M L p (x, t) = N 2 3 st T 0 τ (t) -2p 2p+3 Ω (det |H u (x, t)|) p 2p+3 dx dt -2 3 τ (t) 2 2p+3 (det |H u (x, t)|) -1 2p+3 |H u (x, t)| .
(3.12)

We get the following optimal error:

E L p (M L p ) = 3 N -2 3 st T 0 τ (t) -2p 2p+3 Ω (det |H u (x, t)|) p 2p+3 dx dt 2p+3 3p
.

(3.13)

Error analysis for time sub-intervals

The previous analysis provides the optimal size of the adapted meshes for each time level. Hence, this analysis requires the mesh to be adapted at each flow solver time step. In practice this approach involves a very large number of remeshings which is very CPU consuming and spoils solution accuracy due to many solution transfers. In consequence, an adaptive strategy has been proposed in [Alauzet 2007, Alauzet 2011b] where the number of remeshings is controlled (thus drastically reduced) by considering a coarse adapted discretization of the time axis, and generating adapted meshes for several solver time steps.

The idea is to split the simulation time interval into n adap sub-intervals [t i-1 , t i ] for i = 1, .., n adap . Each spatial mesh M i is then kept constant during each sub-interval [t i-1 , t i ]. We could consider this partition as a time discretization of the mesh adaptation problem. In other words, the number of nodes N i of the i th adapted mesh M i on sub-interval [t i-1 , t i ] should for example be taken equal to:

N i = t i t i-1 N opt (t)τ (t) -1 dt t i t i-1 τ (t) -1 dt
.

Here, we propose a different option in which we get an optimal discrete answer.

Spatial minimization on a sub-interval. Given the continuous mesh complexity N i for the single adapted mesh used during time sub-interval [t i-1 , t i ], we seek for the optimal continuous mesh M i L p solution of the following problem: (3.14) where matrix H i u on the sub-interval can be defined by either using an L 1 or an L ∞ norm:

min M i E i L p (M i ) = Ω trace (M i ) -1 2 (x) H i u (x) (M i ) -1 2 (x) p dx such that C(M i ) = N i ,
H i L 1 (x) = t i t i-1 |H u (x, t)| dt or H i L ∞ (x) = ∆t i max t∈[t i-1 ,t i ] |H u (x, t)| , with ∆t i = t i -t i-1 .
Processing as previously, we get the spatial optimality condition:

M i L p (x) = (N i ) 2 3 M i L p ,1 (x)
with

M i L p ,1 (x) = Ω (det H i u (x)) p 2p+3 dx -2 3 (det H i u (x)) -1 2p+3 H i u (x).
The corresponding optimal error E i (M i L p ) writes:

E i L p (M i L p ) = 3 p (N i ) -2p 3 Ω (det H i u (x)) p 2p+3 dx 2p+3 3 = 3 p (N i ) -2p 3 K i . (3.15)
where

K i = Ω det H i u (x) p 2p+3 dx 2p+3 3
.

To complete our analysis, we shall perform a temporal minimization. Again, we consider the case where the time step τ is specified by a function of time and.

Temporal minimization for specified τ . After the spatial minimization, the temporal optimization problem reads:

min M E L p (M) p = n adap i=1 E i L p (M i L p ) = 3 p n adap i=1 (N i ) -2p 3 K i
under the constraint:

n adap i=1 N i t i t i-1 τ (t) -1 dt = N st .
We set the one-to-one mapping:

Ñ i = N i t i t i-1 τ (t) -1 dt and Ki = K i t i t i-1 τ (t) -1 dt 2p 3
, then the optimization problem reduces to:

min M n adap i=1 ( Ñ i ) -2p 3 Ki such that n adap i=1 Ñ i = N st .
The solution is:

Ñ i opt = C(N st ) ( Ki ) 3 2p+3 with C(N st ) = N st n adap i=1 ( Ki ) 3 2p+3 -1 ⇒ N i = N st   n adap i=1 (K i ) 3 2p+3 t i t i-1 τ (t) -1 dt 2p 2p+3   -1 (K i ) 3 2p+3 t i t i-1 τ (t) -1 dt -3 2p+3 .
and we deduce the following optimal continuous mesh M L p = {M i L p } i=1,..,n adap and error:

M i L p (x) = N 2 3 st n adap i=1 (K i ) 3 2p+3 ti ti-1 τ (t) -1 dt 2p 2p+3 -2 3 ti ti-1 τ (t) -1 dt -2 2p+3 (det H i u (x)) -1 2p+3 H i u (x) , (3.16)
E L p (M L p ) = 3 N -2 3 st n adap i=1 (K i ) 3 2p+3 t i t i-1 τ (t) -1 dt 2p 2p+3 2p+3 3p , (3.17)
with

K i = Ω (det H i u (x)) p 2p+3 dx 2p+3 3
.

Global fixed-point mesh adaptation algorithm for unsteady simulation

Finally, the unsteady adaptation algorithm can be derived from this error analysis. Three main ideas govern this algorithm:

• It is based on splitting the simulation into sub-intervals.

• It is an iterative fixed point algorithm.

• It is a global algorithm.

The methodology consists in splitting the simulation time frame [0, T ] into n adap adaptation sub-intervals:

[0, T ] = [0 = t 0 , t 1 ] ∪ . . . ∪ [t i , t i+1 ] ∪ . . . ∪ [t n adap -1 , t n adap ] ,
and to keep the same adapted mesh for each time sub-interval. On each sub-interval, the mesh is adapted to control the solution accuracy from t i to t i+1 . Consequently, the time-dependent simulation is performed with n adap different adapted meshes. This drastically reduces the number of remeshing during the simulation, hence the number of solution transfers. This can been seen as a coarse adapted discretization of the time axis, the spatial mesh being kept constant for each sub-interval when the global space-time mesh is visualized, thus providing a first answer to the adaptation of the whole space-time mesh.

We have seen in Chapter 1 that mesh adaptation is a non-linear problem, and that this is addressed with iterative algorithms to converge the mesh/solution couple. To that end we propose a fixed-point mesh adaptation algorithm. This is also a way to predict the solution evolution and to adapt the mesh accordingly.

Previously [Alauzet 2007], the optimal metric of a sub-interval could be computed directly once the simulation on the sub-interval had been run. However, the computation of the optimal continuous mesh for sub-intervals given by Relation (3.16) involves a global normalization term which requires the knowledge of quantities over the whole simulation time frame. In this case, the normalization term is:

N 2 3 st T 0 τ (t) -2p 2p+3 Ω (det |H u (x, t)|) p 2p+3 dx dt -2 3 ,
which requires to know all the time steps τ (t) and Hessians H u (x, t) over time frame [0, T ].

Thus, the complete simulation must be performed before evaluating any continuous mesh.

Against, considering a global fixed-point mesh adaptation algorithm covering the whole time frame [0, T ] enables to compare the normalization term. All the solutions and Hessian-metrics are computed, and only then can the global normalization term and thus the metrics for each sub-interval be computed. This algorithm is schematized in Algorithm 8 where H, S and M denote respectively meshes, solutions and metrics, and H is the Hessian-metric.

Algorithm 8 Mesh Adaptation Loop for Unsteady Flows Initial mesh and solution(H 0 , S 0 0 ) and set targeted space-time complexity N st // Fixed-point loop to converge the global space-time mesh adaptation problem For j = 1, n ptf x 1. // Adaptive loop to advance the solution in time on time frame [0, T ]

For i = 1, n adap (a) S j 0,i = Interpolate conservatively next sub-interval initial sol. from (H j i-1 , S j i-1 , H j i ); (b) S j i = Compute solution on sub-interval from pair (S j 0,i , H j i ); (c) |H| j i = Compute sub-interval Hessian-metric from sol. sample (H j i , {S j i (k)} k=1,nk );
EndFor 2. C j = Compute space-time complexity from all Hessian-metrics ({|H| j i } i=1,n adap );

3. {M j i } i=1,n adap = Compute all sub-interval unsteady metrics (C j , {|H max | j i } i=1,n adap ); 4. {H j+1 i } i=1,n adap = Generate all sub-interval adapted meshes ({H j i , M j i } i=1,n adap );
EndFor 3.2 From theory to practice

Computation of the Hessian metric

Since Hessians are the basis of the computation of metrics, it is important to compute them with enough accuracy.

Gradient and Hessian recovery techniques

From a piecewise linear function u h , the chosen method must be able to recover a smooth Hessian field. Two main approaches coexist: a double L 2 projection method, and a weighted least square method. The second method is detailed in [Menier 2015]. The first one is the one generally used in this thesis, and works as follows.

Nodal gradients. Let K be an element and (P i ) i=0...3 its vertices. Let x ∈ K. u h is written:

u h (x) = 3 j=0 u h (P j ) ϕ j (x) ,
where (ϕ j ) j are the P 1 shape functions, given by:

∇ x ϕ 0 (x) = 1 6|K| η 0 , ∇ x ϕ 1 (x) = 1 6|K| η 1 , ∇ x ϕ 2 (x) = 1 6|K| η 2 , ∇ x ϕ 3 (x) = 1 6|K| η 3 .
The expression is differentiated term by term:

∇u h | K = n j=0 u h (P j ) ∇ x ϕ j (x) .
As ∇u h is not defined at the vertices of the mesh, the nodal gradients are recovered from the gradients to the elements using a Clément's L 2 local projection.

Let P i be a vertex of H. The stencil of the shape function ϕ i is the ball of vertices around P i , Ball(P i ). The following spaces are introduced:

V 0 h = {v ∈ L 2 (Ω) v |K ∈ P 0 ∀K ∈ H} V 1 h = {v ∈ C 0 (Ω) v |K ∈ P 1 ∀K ∈ H} .
where P 0 and P 1 are the set of polynomials of degree respectively 0 (constant polynomials) and 1.

For v ∈ L 2 (Ω), we set Π L 2 v ∈ V 0 h : ∀ Ball(P i ) ⊂ H ,    (Π L 2 v) |Ball(P i ) ∈ P 0 Ball(P i ) (Π L 2 v -v) w = 0 , ∀w ∈ P 0 . Clement's operator Π c : V 0 h -→ V 1 h is defined as follows: Π c v = n i=0 Π L 2 v(P i ) ϕ i .
With Clément's operator, we can recover nodal gradients from the ∇u h ∈ P 0 . For each Ball(P i ) ⊂ H we can in particular take v = 1 ∈ P 0 as test function:

Ball(P i ) (Π L 2 (∇u h ) -∇u h ) = 0 ⇐⇒ Ball(P i ) Π L 2 (∇u h ) = Ball(P i ) ∇u h ⇐⇒ |Ball(P i )| Π L 2 (∇u h ) |Ball(P i ) = K j ∈Ball(P i ) K j ∇u h ⇐⇒ Π L 2 (∇u h ) |Ball(P i ) = K j ∈Ball(P i ) |K j | ∇u h|K j K j ∈Ball(P i ) |K j |
where |K| is the volume of element K . For each vertex P i , we can write:

∇ R u h (P i ) = K j ∈Ball(P i ) |K j | ∇u h|K j K j ∈Ball(P i ) |K j |
This procedure comes to recovering the gradient as an average of the gradients at the elements, weighted by the volume of the elements.

We thus have defined values to the vertices for the gradient, and consequently, thanks to Clément's operator, a piecewise P 1 gradient on the mesh.

Hessian recovery. To recover the Hessian matrix by means a double L 2 projection, the same procedure is applied once more, but using the gradient as input.

The optimal L p metric involves an averaged Hessian-metric H i u on sub-interval i, but it still remains to know how to compute it practically, i.e., how it is discretized. The strategy adopted [Alauzet 2007] is to sample the solution on the time sub-interval. More precisely, n k solutions equally distributed on the sub-interval time frame are saved, including the initial solution at t i-1 and the final solution at t i . Positive Hessian |H u (x, t k )| is evaluated for each sample. If samples are distributed regularly in time, the time elapsed between two samples is ∆t i n k -1 . In practice, 20 samples per sub-intervals are usually used.

Once these Hessians are computed, one need to average them. In the previous algorithm [Alauzet 2007], the following discretization was used:

H i L ∞ (x) = ∆t i n k k=1 |H u (x, t k )| = ∆t i |H i max (x)| ,
where ∩ has to be understand as the metric intersection in time of all samples. This corresponds to an integration in time in L ∞ norm of the Hessians.

However, the new error analysis leads to write H i as the integral over time of the Hessian matrices, and thus the following discretization is preferred:

H i L 1 (x) = 1 2 ∆t i n k -1 |H u (x, t i-1 )| + ∆t i n k -1 n k -1 k=2 |H u (x, t k )| + 1 2 ∆t i n k -1 |H u (x, t i )| = ∆t i |H i avg (x)| ,
where ∆t i = t it i-1 is the sub-interval time length and t k = t i-1 + k-1 n k -1 ∆t i . This corresponds to an integration in time in L 1 norm of the Hessians.

A comparative study of both discretizations was performed and will be detailed in section 3.2.2.

Choice of the mean Hessian-metric

All the unsteady adaptation algorithm relies on the computation of a mean metric for each sub-intervals. Instead of averaging real metrics, we average the Hessians on the go during the solver iterations.

There is a choice to make concerning the way this average is done, between the historical intersection or the sum dictated by theory. The intersection of metrics reads:

H i L ∞ (x) = ∆t i n k k=1 |H u (x, t k )| = ∆t i |H i max (x)| ,
It comes to making an average in L ∞ norm, and has the advantage of having a clear geometric interpretation: the intersection of two metrics is the largest metric included in the two metrics.

In the context of unsteady adaptation, this means that the mean Hessian-metric has the smallest sizes met over the sub-interval. The sum of metrics comes to making an average in L 1 norm, and reads:

H i L 1 (x) = 1 2 ∆t i n k -1 |H u (x, t i-1 )| + ∆t i n k -1 n k -1 k=2 |H u (x, t k )| + 1 2 ∆t i n k -1 |H u (x, t i )| = ∆t i |H i avg (x)| ,
It seems to be more coherent with theory. Indeed, when writing the minimization problem (3.14)

with sub-intervals, we can formally derive the error from the unsteady error without subintervals (3.5) by moving the integral over time into the trace. However, it is only a formal manipulation, and the geometrical interpretation of the sum is much less clear: a term-to-term sum of metrics is clearly still a metric (a symmetric positive definite form), but the sizes of the axes of the resulting metric are not the sum of the sizes of the two input metrics. This justifies numerical experimentation to select the best approach in practice.

Choice of the optimal continuous mesh

The optimal adapted mesh for each sub-interval is generated according to analysis of Section 3.1.4. For the numerical results presented below, we select the optimal mesh given by Relation (3.16) and the following particular choice has been made:

• the Hessian-metric for sub-interval i is discretized in time in L 1 norm.

• function τ : t → τ (t) is constant and equal to 1

• all sub-intervals have the same time length ∆t.

Moreover, it is clear that integral

t i t i-1 τ (t)
-1 dt corresponds to the number of time steps (iterations) performed by the flow solver during the i th sub-interval. In practice using the number of iterations of the flow solver for each sub-interval to define the continuous mesh is an issue because it highly depends on the discrete representation of the continuous mesh, i.e., the generated discrete mesh. Thus, the number of iterations of a sub-interval may substantially varies between two fixed-point iterations. To avoid this issue, we prefer considering the flow solver time step constant on each sub-interval, and consider the time step τ (t) constant equal to ∆t (because of the global normalization term). The integral

t i t i-1 τ (t) -1 dt then reduces to 1.
With these choices, the optimal continuous mesh M L p = {M i L p } i=1,..,n adap simplifies to:

M i L p (x) = N 2 3 st   n adap j=1 Ω (det |H j L 1 (x)|) p 2p+3 dx   -2 3 det |H i L 1 (x)| -1 2p+3 |H i L 1 (x)| . (3.18)
In that case, as we assume that theoretically one time step is done by sub-interval, N st /n adap represents the average spatial complexity by sub-interval and thus N st is the total spatial complexity by summing the sub-intervals average complexity. We do not prescribe the temporal complexity, i.e., we do not control the number of time steps done at each sub-interval.

Specification of the space-time complexity. The optimal continuous mesh from Equation (3.16) is obtained for a given space-time complexity that is the number of vertices times the number of solver time steps. Such a complexity should thus be prescribed for each adaptative simulation. However, in practice, it depends on the time discretization of the sub-intervals.

Moreover, as explained above, the number of solver time steps is highly dependent on the sizes of the elements of the mesh, which is not perfectly controlled throughout the algorithm. For these reasons, we prefer to prescribe an average number of vertices per sub-interval, which we sometimes refer to as "simplified space-time complexity" in what follows [Alauzet 2012].

Matrix-free P1-exact conservative solution transfer

For the unsteady adaptation algorithm, the interpolation step is crucial. Indeed, when moving from a sub-interval to the following, the solution has to be transferred from the mesh of the current sub-interval to the mesh of the following sub-interval, and this is done as many times as there are sub-intervals. It is obvious that if information is lost during the transfer process, then the accuracy of the global solution is greatly affected. The following properties need to be satisfied in the interpolation method: (i) mass conservation, (ii) P 1 exactness preserving the second order of the adaptive strategy and (iii) verify the maximum principle. The interpolation scheme used is detailed below.

Interpolation step

The classical approach involves two steps: first the localization of the vertices of the new mesh in the background mesh, and then the application of an interpolation scheme. The localization step is described in several references [START_REF] Frey | [END_REF], Alauzet 2010b]. As concerns interpolation schemes, the easiest is the classical P 1 interpolation, which is written:

u(p) = 3 i=0 β i (p) u(q i ) , (3.19)
where p is a vertex of the new mesh that has been located in tetrahedron K = [q 0 , q 1 , q 2 , q 3 ] of the background mesh and β i are the barycentric coordinates of p with respect to K. Higher order interpolation schemes can be devised, for instance using P 2 Lagrange shape functions.

However these schemes do not conserve the mass.

A new conservation scheme was devised to match these criteria [Alauzet 2010b]. The mass conservation property of the interpolation operator is achieved by local mesh intersections, i.e., intersections are performed at the element level. The use of mesh intersection to build a conservative interpolation process seems natural for unconnected meshes. The idea is to find, for each element of the new mesh, its geometric intersection with all the elements of the background mesh it overlaps and to mesh this geometric intersection with simplices. We are then able to use a Gauss quadrature formula to exactly compute the mass which has been locally transferred.

High-order accuracy is obtained through the reconstruction of the gradient of the solution from the discrete data and the use of some Taylor formulas. Unfortunately, this high-order interpolation can lead to a loss of monotonicity. The maximum principle is recovered by correcting the interpolated solution in a conservative manner. Finally, the solution values at vertices are reconstructed from this piecewise linear by element discontinuous representation of the solution.

The approach is summarized in Algorithm 9:

Algorithm 9 Conservative Interpolation Process Piecewise linear (continuous or discontinuous) representation of the solution on H back 1. For all elements K back ∈ H back , compute solution mass m K back and gradient ∇ K back 2. For all elements K new ∈ H new , recover solution mass m Knew and gradient ∇ Knew :

(a) compute the intersection of K new with all K i back ∈ H back it overlaps (b) mesh the intersection polygon/polyhedron of each couple of elements (K new , K i back ) (c) compute m Knew and ∇ Knew using Gauss quadrature formulas =⇒ a piecewise linear discontinuous representation of the mass on H new is obtained 3. Correct the gradient to enforce the maximum principle 4. Set the solution values to vertices by an averaging procedure.

The remeshing step

The remeshing step is also crucial in the adaptation process, as a poorly adapted mesh or a mesh with bad quality elements will impact the accuracy of the solution and spoil the efforts on all other parts of the process. In this paper, the remesher Feflo.a [Loseille 2013] was used.

Feflo.a belongs to the class of 3D anisotropic local remeshers that aims at generating a unit mesh with respect to a prescribed metric field. Its main particularity is to adapt the volume and the surface mesh in a coupled way so that a valid 3D mesh is always guaranteed on output. It uses a unique cavity-based operator (that generalizes standard operators: point insertion, edge removal, edges swapping, point smoothing). This operator is governed by dedicated metricbased quality functions. This allows to reach a good balance between high level of anisotropy, necessary to capture the physical features of the solution, and mesh quality, necessary to ensure the stability and accuracy of the numerical scheme. Feflo.a has several enhancements designed for CFD computations, including explicit control and optimization of tetrahedra to ensure a maximal time step for unsteady simulations, and surface mesh re-projection based either on CAD or on background discrete meshes.

Software used

Our implementation of the mesh adaptation algorithm described above requires successively four different software components:

• Wolf the second order Finite-Volume flow solver, whose ALE version was described in Chapter 2 -when fixed-mesh simulations are run, a standard version is used, but the numerical schemes are the same as the ALE ones, without the extra ALE velocity terms -• Metrix to compute the continuous space-time mesh and perform the metric fields gradation,

• Feflo.a the local adaptive remesher based on the cavity operator, described in Section 3.2.5,

• Interpol for the P 1 -conservative solution transfer, described in Section 3.2.4.

Unsteady mesh adaptation for dynamic meshes

The strategy adopted to move the meshes is the one described in Section 2.1 : one body-fitted mesh is deformed to follow the moving boundaries, and an Arbitrary-Lagrangian-Eulerian (ALE) formulation of the equations is solved on this moving mesh. This moving mesh strategy preserves the number of vertices, which is a requirement in our error previous analysis framework. An ALE metric, which brings a first answer to this problem, was proposed in [Olivier 2011b] and was improved in [N. [START_REF][END_REF].

Section 1.3.2 provides the optimal instantaneous ALE continuous mesh which takes into account the mesh deformation. Now, as stated previously, we can extend the space-time error analysis with time sub-intervals done for fixed meshes to the case of dynamic meshes. The simulation time interval is still split into n adap sub-intervals. On each sub-interval, the mesh size (number of vertices) remains constant, but the mesh is deformed to follow the geometry displacement. At each time-step of the sub-interval, we want the moved mesh to be adapted to the current sensor. The key idea to perform the error analysis is to seek for the optimal dynamic continuous mesh at the beginning of the sub-interval, this continuous mesh being optimal for the whole sub-interval when deformed, instead of seeking for the expression of the optimal continuous mesh at each instant, i.e., as a function of the time.

The optimal space-time ALE continuous mesh is designed to fit in the global fixed-point unsteady mesh adaptation algorithm described in Algorithm 10. Nevertheless, a few things have been modified to extend this algorithm to moving mesh ALE simulations. In fact, first, geometries move so mesh must be moved with all the difficulties that this implies to preserve accuracy and a good mesh quality all along the movement.

Algorithm 10 Feature-based mesh adaptation for unsteady flows with moving geometries Input : Initial mesh and solution (H 0 , S 0 0 ) and set targeted space-time complexity N st # Fixed-point loop to converge the global space-time goal-oriented mesh adaptation problem For j = 1, n ptf x # Adaptive loop to advance the solution in time on time frame [0, T ] 1. For i = 0, n adap (a) (S j 0,i )= Interpolate conservatively next sub-interval initial solution from pair (H j i-1 , S j i-1 , H j i ); (b) For n = 0, n sample i. (S j i )= Compute solution with the flow solver on sub-interval from pair (H j i , S 0 i ); ii. |H| j i = Compute sub-interval Hessian-metric from solution sample (H j i , = Compute all sub-interval metrics according to error estimate from (C j , {|H| j i }

{S j i (k)} nk k=1 ); iii. M ALE,n L p = Compute ALE metric;
n adap i=1 ); 4. {H j+1 i } n adap i=1
= Generate all sub-interval adapted meshes from pair ({H j i , M j L p ,i } n adap i=1 );

End For

Handling the swaps : update of the adaptation algorithm

The mesh optimization procedure of the connectivity-change moving mesh strategy requires a proper metric field to evaluate geometric quantities and elements qualities of anisotropic adapted meshes. But, at a given fixed-point iteration j and at a given time sub-interval i, the input is the optimal ALE continuous mesh M i,ALE L p defined on Ω(t i ). Therefore, when the mesh is deformed during the sub-interval due to the geometry displacement, the input metric is no more compatible with the moved mesh at a given time t. The use of an incorrect metric field will drive the smoothing to move vertices to a wrong location and the swaps to break the anisotropy, thus spoiling the adaptation of the mesh. As the mesh is evolving in time, the input metric field (i.e., the input continuous mesh) must also evolve in time. Consequently, we have to apply the deformation of the adapted mesh to the continuous mesh in order to maintain the consistency between the discrete and the continuous meshes.

Metric for optimizations

But, in the theory developed in Chapter 1, the mesh deformation correction is applied to the Hessian-metric |H i,ALE u | and not directly to the continuous mesh. Therefore, from the optimal ALE continuous mesh

M i,ALE L p (x(t i )) x∈Ω(t i )
at the beginning of sub-interval [t i , t i+1 ],

we cannot directly find the deformed optimal ALE continuous mesh M i,ALE L p (x(t))

x∈Ω(t) for t ∈ [t i , t i+1 ].
Two possibilities arise from these remarks. First, several ALE metric field samples denoted

M i,ALE L p (x(t k )) x∈Ω(t k ) k
can be pre-calculated for each sub-interval at the previous fixedpoint iteration, and when an ALE metric field is required at time t, it is linearly interpolated in time on the current mesh from two of these pre-calculated ALE metric fields. This method has two drawbacks: it requires computing and saving all these ALE metric fields and it requires a metric interpolation stage at each mesh optimization which leads to consequent memory, I/Os and CPU time overhead. The second possibility is to modify the input ALE metric field according to the current mesh deformation to get a consistant metric field at time t. To deform the continuous mesh defined at the beginning of the sub-interval, we adopt the same reasoning as the one that leads to the optimal instantaneous ALE continuous mesh in Section 1.3.2, but on a reverse time frame.

Let us consider t k a time in the sub-interval [t i , t i+1 ] and d the displacement field of the mesh between t i and t k . Here, the problem is reversed, we seek for the continuous mesh at time t k that will result in the continuous mesh at time t i once moved with displacement -d. At time t i , the continuous mesh is M i,ALE L p (the adapted mesh at time t i is generated directly using this metric field).

Displacement d is the displacement of vertices between t i and t k : x k = x i + d(x i ) and we write d (x k ) = x ix k = -d(x i ). We define :

φ : Ω k -→ Ω i x k -→ x i = φ x k = x k + d (x k ) .
(3.20)

We start from an adapted mesh which is unit for metric field M i,ALE L p at time t i and we search for the expression of the metric field at time t k for which the deformed adapted mesh at time t k is unit. Following Section 1.3.2, we consider an edge e k at time t k having edge e i as image by φ at time t i . They verify:

1 = e i (x i ) T M i,ALE L p x i e i (x i ) = e i φ (x k ) T M i,ALE L p x i e i φ (x k ) = ∇ k φ (x k ) T e k (x k ) T M i,ALE L p (x i ) ∇ k φ (x k ) T e k (x k ) = e k (x k ) T ∇ k φ (x k ) M i,ALE L p (x i ) ∇ k φ (x k ) T e k (x k ) .
In other words, the metric field M i,optim L p we are looking for at time t k is:

M i,optim L p (x(t k )) = ∇ k φ x(t k ) M i,ALE L p x(t i ) ∇ k φ x(t k ) T . (3.21)
This metric field is easy to compute on the fly, because the central term is the input metric field M i,ALE L p

x(t i ) which is known, and the gradients of φ are easy to compute, since the displacement considered in φ is the opposite of the current displacement at each vertex. Note that the gradients of φ should be computed on Ω(t k ) that is to say on the current mesh.

Modification of the algorithm

The overall global fixed-point mesh adaptation algorithm remains unchanged from [Barral 2015].

the simulation time frame is still divided into sub-intervals, a global fixed-point strategy is used to converge the meshes and the solutions within a sub-interval and the mesh is moved using the connectivity-change moving mesh algorithm presented On each sub-interval, the mesh number of vertices remains the same. One or several mesh deformation steps are performed during a sub-interval. The connectivity-change moving mesh algorithm uses mesh optimizations to make the mesh deformation efficient and robust. These mesh optimizations, smoothing and swapping operators, use the dynamic (corrected) metric field to be compliant with the current dynamic adapted mesh.

The difference concerns the computation of the ALE Hessian-metric for the next adaptation step. ALE Hessian-metrics |H i,ALE u | need to be sampled in time. To compute each sample, we use the displacement given by the mesh deformation step without taking into account the smoothing optimization. In numerical experiments, we did not see any significant difference on the solution accuracy and the quality of the adapted meshes by considering or not the smoothing in the ALE Hessian-metrics computations. This is due mainly to the fact that the mesh smoothing has only a slight influence in 3D, the mesh displacement is predominantly governed by the mesh deformation [Alauzet 2014]. The mesh smoothing is only a slight perturbation around the mesh deformation displacement to get better shaped elements. Nevertheless, we prefer not to consider the smoothing displacement when computing ALE Hessian-metrics because:

• the geometry displacement is clearly defined, thus the mesh deformation is converging toward a unique mesh deformation and therefore we converge toward a unique optimal ALE continuous mesh

• the smoothing correction for one node may change from one fixed-point iteration to the other, thus for each node it does not converge to a fixed correction.

The gradient of the mesh deformation displacement field is computed on the mesh position at the beginning of the sub-interval. The Hessian of the sensor is computed on the mesh current position and is mapped back to the mesh position at the beginning of the sub-interval. Then, as in the standard algorithm, the ALE Hessian-metrics are computed at each vertex.

Note that the connectivity-change moving mesh algorithm allows the mesh vertices to move to their final position given by the mesh deformation solution while it would be impossible to reach such positions without skewing mesh elements with a classical method. Thanks to the robustness of this algorithm, we are always able to calculate the ALE hessian-metrics.

⌦ i ⌦ k u k r c H k H k M i,ALE d H i,ALE adap H k adap ⌦ k d 0 , 0 r 0 d H i,ALE adap H k adap M i,ALE
Mk,ALE ⌦ i

Procedure

Let H n 0 be an uniform mesh, d a displacement field between two times t n and t n+1 , and u n+1 a sensor at the final time. The goal is to generate a mesh H n at the initial time that, once moved into H n+1 ALE = d(H n ) will be adapted to the sensor. At the same time, a reference adapted mesh is computed, that is directly adapted to sensor u n+1 . To do so, the classical steady-state procedure described in Algorithm 8. Since the sensor is analytic and to be fair with the ALE version, only one iteration of the algorithm is performed. The metric used is normalized to look like the regular optimal metric for unsteady simulations, obtained from the hessian of the sensor. In what follows, we note H n+1 ALE the mesh obtained with the ALE procedure and H n+1 ref the mesh directly adapted.

∇ n φ = ComputeTransformationGradient(H n 0 , d) H n+1 0 = MoveMesh(H n 0 , d) u n+1 = ComputeTargetSensor(H n+1 0 ) M n,ALE L p = ComputeALEMetric(H n+1 0 , u n+1 , ∇ n φ) H n = AdaptMesh(H n , M n,ALE L p ) H n+1 ALE = MoveMesh(H n , d) M n+1 L p = ComputeMetric(H n 0 , u n+1 , H u n+1 ) H n+1 ref = AdaptMesh(H n 0 , M n+1 L p ) .
According to the above developments, we expect the mesh H n+1 , obtained by moving the vertices of H n with d, to be optimal for the control of the interpolation error of u n+1 in L 1 norm.

To study the adapted characteristic of the resulting meshes, two quantities are considered:

• the mesh quality, computed using the metric M n+1 L p of the reference adapted mesh (i.e. the one that was adapted directly, .... The quality measure is the one recalled in Section 1. 1.2,Equation (1.19).

• the error comitted when the sensor function is projected on the mesh remains the best way to evaluate if the mesh is really adapted to a sensor. Here we consider the interpolation error Π h uu for sensor u, where Π h u is the piecewise P 1 function with Π h u(P i ) = u(P i ) for any vertex P i of the mesh.

The displacement considered are straight-line displacements of somewhat large amplitude, whereas infinitesimal displacements have been considered in the theoretical analysis. However, this is closer to real life simulations, with large displacements even within a sub-interval. Note that the ALE metric does not guarantee that the mesh moved at time t n+1 is still valid. In practice, mesh optimizations deal with this issue and ensure the validity of the mesh together with its quality. In this study, no mesh optimization is performed on the adapted mesh, or during the adaptation or moving process, to avoid altering the results of the quality analysis.

For the same reason, unlike in the actual adaptation loop, no gradation is performed.

Analytic function considered

The function here is taken from [Olivier 2011a] in order to confirm and complete the analysis proposed in the present work. Its sensor presents specific features (smoothness, discontinuities, features of different scales...) that are as many difficulties for the adaptation process, so that handling them correctly shows the robustness of the approach. As for the displacements, they are set so that the displacement is null on the boundaries.

d(x, y) =               -0.3 (x + 1) y 2 -1 exp -5x 2 , if x ≥ 0 0.3 (x -1) y 2 -1 exp -5x 2 , if x < 0    -0.3 x 2 -1 (y + 1) exp -5y 2 , if y ≥ 0 0.3 x 2 -1 (y -1) exp -5y 2 , if y < 0           
.

Results

In Figure 3.4, the mesh on the left is the one that was adapted directly quality, computed using the metric M n+1 L p of the reference adapted mesh (i.e. the one that was adapted directly maximal angle reached, the movement is defined by its angle of rotation, of which the evolution is divided in 7 phases:

θ(t) = θ max                                  0 if 0 ≤ t ≤ T /2 (i) 2(t-T 2 ) 2 T 2 if T 2 < t ≤ T (ii) 1 2 + 2 ( 2 (t-T 2 ) T -1) -1 2 ( 4 (t-T 2 ) 2 T 2 -1) if T < t ≤ 3T 2 (iii) 1 if 3T 2 < t ≤ 7T 2 (iv) 1 - 2 (t-7T 2 ) 2 T 2 if 7T 2 < t ≤ 4T (v) 1 2 -2 ( 2 (t-7T 2 ) T -1) + 1 2 ( 4 (t-7T 2 ) 2 T 2 -1) if 4T < t ≤ 9T 2 (vi) 0 if 9T 2 t ≤ 5T (vii)
Phase (i) is an initialization phase, during which the flow around the aircraft is established.

Phases (ii) and (iii) are respectively phases of accelerated and decelerated ascension. Vortices start to grow behind the aircraft, and they expand during phase (iv), where the aircraft stays in upward position. Phases (v) and (vi) are phases of accelerated and decelerated descent, the vortices start to move away and they slowly disappear in phase (vii). Free-stream conditions are imposed on the faces of the surrounding box, and slipping conditions on the aircraft.

The time simulation interval [0, 5] was divided into 96 sub-intervals, and 5 fixed-point iterations are performed. The adaptation sensor is the local Mach number. A space-time complexity of 48 million was prescribed leading to a sub-interval average spatial complexity N avg = 500, 000.

The adapted meshes at the last fixed-point iteration have between 100, 000 and 1, 400, 000 vertices, for an average number of vertices equal to 841, 000. A total of n iter = 38, 980 time-steps have been performed in the last fixed-point iteration. The total number of space-time vertices N st at the last fixed-point iteration is equal to 35 billion. The total CPU time of this simulation is 95 hours. 31 hours being spent in the last fixed-point iteration (error estimate, metric gradation, adaptive mesh generation, solution interpolation, and flow solver).

Views of the adapted meshes and the corresponding solutions are displayed in Fig. 3.6. The vortical wake is propagated far from the aircraft, and the patterns of the vortices are highly resolved. The observation of the adapted meshes shows that they are actually refined only in the vicinity of the wake, and with such a precision that one can see the vortices being created, evolving and vanishing just by looking at the meshes. A view at several meshes within the same sub-interval shows that the mesh evolves continuously following the physical phenomena.

In the wake far from the aircraft, the elements are highly anisotropic, whereas they have to be isotropic in the area of vorticity due to the characteristic of the physical phenomena.

Conclusion

This chapter addresses metric-based mesh adaptation for time dependent and dynamic meshes simulations. In the first sections, a brief survey of mesh adaptation is presented, and the metric-based adaptation concept is detailed, notably the error analysis and choose of metric optimization. This concept is used to derive an improved version of the previous unsteady fixedpoint mesh adaptation algorithm. Notably, a global normalization term arises, which results in a global fixed-point algorithm. The method and recovery to compute the Hessian metric was considered. An efficient parallelization of the algorithm was carried out, which makes it possible to run a number of unsteady 3D simulations with an accuracy unthinkable a few years ago. The ALE metrics proposed in [Olivier 2011a] then in [Barral 2015] allow to generate meshes that will be adapted once moved with a certain displacement. Analytic examples validated this metric, and it was used within a new error analysis, to derive an ALE metric for adaptive sub-intervals.

The regular time-dependent algorithm was modified to use this metric, and coupled with the moving mesh algorithm and example is done on a 2D analytic example.

The main contributions of this thesis with respect to pre-existing work in the GAMMA group is:

• An update strategy was proposed for the metric at mesh optimization times, to avoid spoiling the adaptation with optimizations.

• A study of the error on analytic case for the ALE metric that confirmed its adapted characteristic.

• A new version of the global fixed-point unsteady mesh adaptation algorithm was used.

Several perspective arise from this work, notably:

• To improve the adaptation algorithm, the next step is to control the temporal error. This would in particular lead to automatic adaptation of the size of the sub-intervals. This appears to be all the more important with moving meshes, where the size of the mesh is not even constant throughout a sub-interval. This improvement is necessary to deal with the ALE implicit scheme presented in 2.

• The convergence of this unsteady adaptative method needs to be studied, depending on the number of vertices and/or the size of the sub-intervals, as well as the influence of the initial solution on the final solution, and the ability of the algorithm to recover features of the solution lost due to coarse initial meshes.

• The metric used for optimizations during the adaptation loop should be analyzed in detail.

Indeed, this step is crucial in the adaptation process, since it controls where the elements actually go.

• Vertices should be allowed to move on non plane surfaces. However, this is very complicated, since it requires reconstructing the surfaces with a good geometric accuracy to avoid introducing errors due to bad surface discretization.

• At longer term, adaptative moving mesh methods mentioned in [?] appear as an interesting alternative to prescribe the displacement of the vertices, and possible interactions with metric based unsteady adaptation should be investigated.

Chapter 4

Adjoint Methods

Adjoint techniques have gained widespread use over past decades as simulation capabilities have progressed. The ability to compute not only the simulation but to focus also on a functional output of interest has become increasingly recognized as a significant capability for enabling accuracy. We count nowadays many applications of the adjoint approach. The most classical ones are in control theory [Lions 1971], optimization [Allaire 2007], sensitivity analysis [Allaire 2015],

inverse problems [A. [START_REF] Tikhonov | [END_REF] and data assimilation [F.X. Le Dimet 1986]. The aim of this chapter is not to be exhaustive but to report on the universality nature of adjoint methods.

Moreover understanding its use in other domains can lead to good ideas relevant to the implementation of our subject. In CFD, their utility is also reflected in numerous applications, mainly in design optimization [Giles 2000, Mohammadi 1997, Mavriplis 2000, A.Stück 2009] and lesscommonly in mesh adaptation [START_REF] Venditti | [END_REF], Loseille 2010c[START_REF] Belme | [END_REF], Fidkowski 2017].

In the first domain, the adjoint method is the only way to compute a gradient and perform efficient optimization due to the large number of design variables. In the domain of mesh adaptation, the use of adjoint computation is not so common. It is therefore appropriate : if adjoint methods can be used to compute sensitivities of a functional then they can be used to determine regions of the computational domain that are most sensitive to a given functional and provide an error estimation in a specific quantity of interest, rather than an estimate of total simulation error. In this case, error estimation purposes, which in turned can be used to drive adaptive error control schemes in the interest of producing more accurate, reliable and predictive simulation outcomes. Among linearisation methods, adjoint methods demonstrate their efficiency. For many engineering problems, this is an important advantage as it allows to focus computational resources on optimizing the accuracy of the specific simulation objective while de-emphazing other aspects of the simulation which are less important for the objective of interest. The adjoint-based error method distinguishes from other approaches that may have little relation to errors in integral outputs of primary concern to the engineer. As an example, one can consider the wake behind an aircraft. In order to resolve the wake with accuracy, a finer grid is locally required. But it is often the case that the computed wake, located a chord or two downstream of the wing, passes into a region represented by a rather coarse grid. In this case, a mesh adaptation based on an usual error estimation, as seen on the previous chapter, would probably induce further refinement in this region. However, the effort of refinement coud be useless if errors in this region has a very small contribution to computed lift and drag.

Probably, refinement in region closer to the wing, near the leading and trailing edges would lead to much greater reduction in the lift and drag errors.

The use of adjoint methods for estimating errors in simulation output objectives is now fairly well-established, particularly for spatial discretization errors in steady-state aerodynamics problems [START_REF] Venditti | [END_REF], Loseille 2010c]. Additionally adjoint methods for unsteady problems have been pursued over the last decade by various contributors [START_REF] Belme | [END_REF], Fidkowski 2017].

But in spite of theses advances, adjoint methods have remained confined to relatively small subset of simulation problems and there is a need to extend these methods to more complex moving geometries problems to justify their interest.

State of the art of adjoint method

In this section, let's take at first a back step to consider the mathematical environment to know where the method comes from.

Mathematical point of view

In mathematics, the term "adjoint" applies in several situations and one can encounter this term in linear algebra, in functional analysis, but also in group theory or further in category theory. Their common particularity is that there is always a link to the principle of duality which establishes a full meaning of the term "adjoint". And duality is an old and complex principle because it enables to view a same problem with two different points of view. The interest of this rewrite is that the two problems are closely related and takes to the great advantage that for various reasons, it may sometimes be easier to solve the dual problem than the primal problem.

But this is also which makes this own major difficulty : one have to also consider the primal problem with an other perspective : the dual problem and leads to two spaces, two equations, two sets of variables, etc.

Let's make an aside to illustrate the commodity but, in the same time, the complexity of duality with the linear duality in plane geometry. Anyone knows the assumption : "Two different points can be joined by an unique line". The dual of this assumption is : "Two different lines meet in one point". But the vision of classical Euclidian geometry does not allow to check this last assumption in every cases. In fact two parallel lines never join. The dual space associated to this dual assumption needs an addition of a point at infinity and is called the projective plane.

Now the relationship between points and lines is perfectly symmetrical in this space.

It can be stated that in mathematical optimization theory, the adjoint method has its origin in the theory of Lagrange multipliers. The first use was probably by the mathematician and one of the founder of the variational calculus J.L Lagrange. He introduced it for Lagrangian functions in dynamical systems. The method of Lagrange multipliers is a strategy for finding the local maxima or minima of a function under equality constraints. The great advantage of it is that the optimization method can be solved without explicit parametrization in terms of the constraints. We recall the method as follows.

Consider the problem : Let's f and g have continuous first partial derivatives, maximize f (x)

under the condition g(x) = 0.

To solve this problem we introduce a new variable λ called a Lagrange multiplier and study the Lagrangian function defined by :

L(x, λ) = f (x) -λ • g(x).
If f (x 0 ) is a maximum of f for the constrained problem then there exists λ 0 such that (x 0 , λ 0 ) is a stationary point for L (i.e its partial derivatives at this point are zero).

Interpretation

The Lagrangian multipliers measure the sensitivity of the optimal value by relation to the variation of the constraint. One way of looking at them is that they give the influence of an arbitrary source term on the functional of interest.

The adjoint space is chosen to simplify the physical interpretation of equation constraints. It may stands as a Hilbert space. In the field of functional analysis, each bounded linear operator on a Hilbert space has a corresponding adjoint operators.

Numerical adjoint approaches

The basic derivation for adjoint error estimation can be performed in several manners in the literature. And depending on the sphere of use and underlying motivations the adjoint could be considered as "continuous" or "discrete". Let's make in this section a quick overview of the use of adjoint in numerical computations before returning to the scope of output-based mesh adaptation.

Continuous adjoint approach : In the continuous formulation, the continuous governing equations are first derived using analytic differentiation and integration by parts and then linearized.

In functional analysis each bounded linear operator on a Hilbert space has a corresponding adjoint operator. Thus from functional analysis theory, a continuous system can be well-posed and considered in these terms and we could say that a continuous adjoint system can be derived for any fuctional output. These however doesn't mean that all functionals lead to properly defined adjoint boundary conditions. Several works in litterature [START_REF] Arian | [END_REF], Castro 2007, A. Bueno-Orovio 2007, Anderson 1997] illustrate this problem and propose solutions usually by adding auxiliary boundary terms. In [Castro 2007], the authors conclude that for Euler flows, only output functional depending on pressure only can be considered. This choice remains a deep challenge when working with continuous adjoint formulation.

Discrete adjoint approach : In the discrete formulation, the governing primal equations are discretized then linearized and transposed for the adjoint problem. The discrete approach formulates directly the adjoint of the forward numerical scheme. This approach is very attractive since these dicrete adjoints can be in general generated automatically.

For discrete adjoint systems, the functionals area choice does not seem to be a problem and as far as the discrete adjoint system is well-posed, any admissible functional (in terms of non-null integral) can be considered.

For consistent discretizations, these two approaches converge to the same result while increasing the resolution of the discretization. However, the consistency of the discretization is not always guaranteed and needs to be performed. For some substantial works in the literature about this consistency for various schemes one can refer to [A. Sei 1995, Giles 2002].

Automatic Differentiation : Automatic Differentiation (AD) is a technique which automatically evaluates the derivatives of a function defined by a computer program. The AD exploits the fact that every computer program executes a sequence of elementary arithmetic operations and functions. To differentiate the code, it applies the chain rule repeatedly to these. Derivatives of arbitrary order can be then be computed automatically. [START_REF] Giles | [END_REF], N.A. Pierce 2004] showed the efficiency of adjoint systems for a better prediction of an output functional associated with a given system of PDE. Finally, let's cite some of AD tools : ADIFOR [Bischof 1998],

TAPENADE [Hascoet 2004], ODYSSEE [C. Faure 1998], OpenAD/F [J. [START_REF] Utke | [END_REF]] and more recently package for Python programs as PyCppAD or PyADOL-C [Walter ].

Our approach : Although continuous or discrete adjoint approaches can be applied in the case of mesh adaptation, in this thesis and to match with the discrete formal solver, we focus on the discrete one for implementation. However, none of the Automatic Differentiation tools is used because the adjoint state has to be highly correlated to the mesh even more with a moving mesh. It is computed in the in-house code solver Wolf.

Computation of adjoint Euler system in ALE formulation

In this section we detail the discrete adjoint system associated to unsteady inviscid Euler flows computed with an ALE formulation [E. [START_REF] Gauci | [END_REF].

Unsteady adjoint Euler equations in ALE formulation

We consider the 3D unsteady compressible Euler equations for a Newtonian fluid in their ALE formulation set in the space-time computational domain Q = Ω(t) × [0, T ], where T is the (positive) maximal time. Assuming that the gas is perfect, inviscid and that there is no thermal diffusion, the continuous unsteady Euler formulation of the equations is written for ALE formulation :

Ψ(W ) = W t + ∇ • F ale (W ) = 0 (4.1)        W = (ρ, ρu, ρE) T is the conservative variables vector F ale (W ) = F(W ) -W ⊗ w
where w is the mesh velocity F(W ) = (ρu, ρuu + pe 1 , ρvu + pe 2 , ρwu + pe 3 , (ρE + p)u) is the flux tensor and we have noted ρ the density of the fluid, p the pressure, u = (u, v, w) its Eulerian velocity, E = 1/2 q 2 + ε the total energy per unit mass, q = u and ε the internal energy per unit mass.

Discrete adjoint state in ALE formulation

Although any consistent approximation of the continuous adjoint system could be built by discretizing the continuous adjoint equation, we choose the option to build the discrete adjoint system from the discrete state system in order to be numerically closer to the true error. In this thesis we use the same notations as in [START_REF] Belme | [END_REF]].

Explicit time integration scheme

Consider the following semi-discrete unsteady compressible Euler model (explicit RK1 time integration):

Ψ n h (W n h , W n-1 h ) = W n h -W n-1 h δt n + Φ ale,h (W n-1 h ) = 0 for n = 1, ..., N. (4.2)
The time-dependent functional is discretized as follows. The functional doesn't depend on the mesh displacement :

j h (W h ) = N n=1 δt n j n-1 h (W n-1 h ).
For the sake of simplicity, we restrict to the case g T = 0 for the functional output defined. The problem of minimizing the error committed on the target functional j(W h ) = (g, W h ), subject to the Euler system (4.2), can be transformed into an unconstrained problem for the following Lagrangian functional [Giles 2000]:

L(W h , W * h ) = N n=1 δt n j n-1 h (W n-1 h ) - N n=1 δt n (W * ,n h ) T Ψ n h (W n h , W n-1 h ) ,
where W * ,n h are the N vectors of the Lagrange multipliers (which are the time-dependent adjoint states). The conditions for an extremum are: The first condition is clearly verified from Relation (4.2). Thus the Lagrangian multipliers W * ,n h must be chosen such that the second condition of extrema is verified. This provides the unsteady discrete adjoint system:

     W * ,N h = 0 W * ,n-1 h = W * ,n h + δt n ∂j n-1 h ∂W n-1 h (W n-1 h ) -δt n (W * ,n h ) T ∂Φ ale,h ∂W n-1 h (W n-1 h ) , (4.3)
or equivalently, the semi-discrete unsteady adjoint model reads:

Ψ * ,n h (W * ,n h , W * ,n-1 h , W n-1 h ) = W * ,n-1 h -W * ,n h -δt n + Φ * ale,h (W * ,n h , W n-1 h ) = 0 for n = 1, ..., N with Φ * ale,h (W * ,n h , W n-1 h ) = ∂j n-1 h ∂W n-1 h (W n-1 h ) -(W * ,n h ) T ∂Φ ale,h ∂W n-1 h (W n-1 h ) .

Implicit scheme

We present in this paragraph the implicit scheme for unsteady adjoint Euler equation. Only the discrete adjoint solution which satisfies the implicit equation system is presented. The numerical implementation has not been done in the solver.

Consider the following unsteady Euler model using implicit time integration :

Ψ n (W n+1 , W n ) = W n+1 -W n δt n + Φ(W n+1 ) = 0
, where Φ is the Euler flux functions, and a time-dependent functional :

j(W ) = N n=1 δt n j n (W n )
The problem of minimizing the error committed on the target functional J n (W n ), subject to the Euler system, can be solved using the method of Lagrange multipliers.

In this sense, we introduce the adjoint state W * which stands as a Lagrange multiplier and the Lagrangian defined by :

L(W, W * ) = N n=1 δt n j n (W n ) - N n=1 δt n (W * ,n ) T Ψ n (W n+1 , W n ),
and study : ∂L ∂W * ,n (W, W * ) = 0 and ∂L ∂W n (W, W * ) = 0. The first condition is clearly verified. The adjoint state is thus chosen such that the second condition is verified.

∂L ∂W n (W, W * ) = δt n ∂j n ∂W n (W n ) -δt n (W * ,n ) T ∂Ψ n ∂W n (W n , W n+1 ) -δt n-1 (W * ,n-1 ) T ∂Ψ n-1 ∂W n (W n-1 , W n ) = δt n ∂j n ∂W n (W n ) + W * ,n -W * ,n-1 -δt n-1 (W * ,n-1 ) T ∂Φ ∂W n (W n ) = 0.
This equation implies :

W * ,n-1 -W * ,n δt n = ∂j n ∂W n (W n ) - δt n-1 δt n (W * ,n-1 ) T ∂Φ ∂W n (W n ) = ∂j n ∂W n (W n ) - δt n-1 δt n W * ,n-1 -W * ,n T ∂Φ ∂W n (W n ) - δt n-1 δt n (W * ,n ) T ∂Φ ∂W n (W n ).
Let us divide by δt n-1 δt n , we get :

Id δt n-1 + ∂Φ ∂W n (W n ) T W * ,n-1 -W * ,n = δt n δt n-1 ∂j n ∂W n (W n ) -(W * ,n ) T ∂Φ ∂W n (W n ).
This provides the unsteady discrete adjoint model :

Ψ * ,n (W * ,n , W * ,n-1 , W n ) = Id -δt n-1 - ∂Φ ∂W n (W n ) T W * ,n-1 -W * ,n + Φ * (W * ,n , W n ) = 0, with Φ * (W * ,n , W n ) = δt n δt n-1 ∂j n ∂W n (W n ) -(W * ,n ) T ∂Φ ∂W n (W n ).
In the following of this thesis, we only consider the explicit time integration scheme.

Connectivity-change moving mesh strategy

The overall connectivity-change moving mesh algorithm in the primal state was described in chapter 2 section 2.1 in Algorithm 5. When coupled with a flow solver (see Sections 2.2), the flow solver is called after the optimization phase. To compute the adjoint state and to be consistent with the connectivity-change moving mesh algorithm in the primal state, a choice must be done for the backward connectivity-change moving mesh algorithm. THe choice made was :

• Store all the coordinates of the vertices at each iteration during the forward loop in order to not recompute the deformation during the backward loop

• Not store the connectivities but remake connectivity-change optimizations during the backward loop.

The first point is very costly in term of memory as we already have to store the states W n-1 h n=1,N to solve the adjoint states. But in this way, the vertices make exactly the reverse movement. However, our choice has been to not store all the connectivities even if it would have been possible by storing a char by edge. Thus during the backward loop optimizations have remade over time. For this reason, we can have different connectivities between forward and backward meshes at the same corresponding time. Let's detail the process.

During the computation of the equation state, the forward meshes are deformed. Each iteration a connectivity-change occurs a tag d * ,opt is used. Therefore, during the computation of the adjoint equation state, the backward meshes are moved. The initial mesh in backward computation is the last mesh of the forward computation. It is moved with the opposite displacement and if an iteration is tagged, an optimization of the mesh is performed.

The backward meshes are then topologically equivalent to the forward meshes but not identical.

An example of forward mesh and backward mesh at the same time iteration is given in

The Discrete Geometric Conservation Law

As said in Section 2.2.3 for the forward mesh computation, we need to ensure that the movement of the mesh backward is not responsible for any artificial alteration of the adjoint computation, or at least, to make our best from a numerical point of view for the mesh movement to introduce an error of the same order as the one introduced by the numerical scheme.

The continuous adjoint state W * , integrated on any arbitrary closed volume C = C(t), writes :

-

C(t)

Backward Implementation

The properties of the backward implementation moving mesh has been detailed in the previous section. In the adjoint solver point of view, it is then necessary to introduce mesh velocities in the ALE formulation. In this section we detail the different Riemann solvers.

The implementation of the adjoint ALE solver was done in the in-house code Wolf.

Numerical flux computation

In the formula 4.3, the functional output is not affected by the ALE formulation. Its computation does not change from [START_REF] Belme | [END_REF]]. The great difference focuses on the Jacobian solver computation. Let's detail some of them.

Jacobian in ALE Rusanov Riemann solver

∂Φ Rusanov ale (W i , W j , n ij , σ ij ) ∂W j = 1 2 ∂F (W i ) • n ij ∂W j -σ ij Id -|λ ij |Id . (4.7)
where

|λ ij | = max(|W i • n ij | + c i -σ ij , |W j • n ij | + c j -σ ij ).
Jacobian in ALE Roe Riemann solver The HLLC flux depends on the waves speed, so does the associated Jacobian.

∂Φ Roe ale (W i , W j , n ij , σ ij ) ∂W j = 1 2 ∂(F (W j ) • n ij ) ∂W j -σ ij Id -| Ã(W i , W j ) -σ ij Id| . ( 4
if

S L > σ lr ∂Φ HLLC ale (W l ) ∂W l = ∂F (W l ) ∂W l -σ lr Id and ∂Φ HLLC ale (W l ) ∂W r = 0 if S L ≤ σ lr ≤ S M ∂Φ HLLC ale (W * l ) ∂W l = ∂F (W * l ) ∂W l -σ lr Id and ∂Φ HLLC ale (W * l ) ∂W r = ∂F (W * l ) ∂W r -σ lr Id if σ lr > S R ∂Φ HLLC ale (W r ) ∂W r = ∂F (W r ) ∂W r -σ lr Id and ∂Φ HLLC ale (W r ) ∂W l = 0 Symmetrically, if S L ≤ σ lr ≤ S M ∂Φ HLLC ale (W * r ) ∂W r = ∂F (W * r ) ∂W r -σ lr Id and ∂Φ HLLC ale (W * r ) ∂W l = ∂F (W * r ) ∂W l -σ lr Id

Memory management

As the adjoint system runs in reverse time, the first expression in the adjoint System (4.3) is referred to as adjoint "initialization".

Solve adjoint state backward: Ψ * (W, W * ) = 0

Solve state foreward: Ψ(W ) = 0

Computing W * ,n-1 h at time t n-1 requires the knowledge of state W n-1 h and adjoint state W * ,n h . Therefore, the knowledge of all states W n-1 h n=1,N is needed to compute backward the adjoint state from time T to 0 which involves large memory storage effort. For instance, if we consider a 3D simulation with a mesh composed of one million vertices then we need to store at each iteration five millions solution data (we have 5 conservative variables). If we perform 1000 iterations, then the memory effort to store all states is 37.25 Gb for double-type data storage (or 18.62 for float-type data storage).

Dealing with memory 4.6.1 Packing and interpolation

After each swap step, a linear interpolation is performed to recover the solution. In other words, the solution at the vertices does not change, which is not conservative. The data of the finite volume cells (volume and interface normals) are then updated, together with the connectivity of the mesh (edges and elements). This approach however is DGCL compliant, since the constant state is preserved (in fact, any linear state is preserved). A better approach would be to use conservative interpolation [Alauzet 2010c] which is conservative and DGCL compliant and the best would be to develop the ALE formulation of the generalized swap in 3D from the formulation in 2D [Olivier 2011a].

Parallelization

Most parts of the code are parallelized for shared memory multi-core computers with a p-threads approach using a semi-automatic p-thread parallelization library [Maréchal 2016] coupled with a space filling curve renumbering strategy [Alauzet 2009].

The automatic parallelization library cuts the data into small chunks that are compact in terms of memory (because of the renumbering), then uses a dynamic scheduler to allocate the chunks to the threads to limit concurrent memory accesses. In the case of a loop performing the same operation on each entry of a table, the loop is split into many sub-loops. Each sub-loop will perform the same operation (symmetric parallelism) on equally-sized portions of the main table and will be concurrently executed. It is the scheduler's job to make sure that two threads do not write on the same memory location simultaneously. When indirect memory accesses occur in loops, memory write conflicts can still arise. To deal with this issue, an asynchronous parallelization is considered rather than of a classic gather/scatter technique, i.e., each thread writes in its own working array and then the data are merged.

However, with this approach, a subtle management of cache misses and cache-line overwrites is required to avoid drops in scaling factors. Space filling curves exhibit clustering properties very helpful for renumbering algorithms [Sagan 1994]. The Hilbert space filling curve based renumbering strategy is used to map mesh geometric entities, such as vertices, edges, triangles and tetrahedra, into a one dimensional interval. In numerical applications, it can be viewed as a mapping from the computational domain onto the memory of a computer. The local property of the Hilbert space filling curve implies that entities which are neighbors on the memory 1D interval are also neighbors in the computational domain. Therefore, such a renumbering strategy has a significant impact on the efficiency of a code. We can state the following benefits: it reduces the number of cache misses during indirect addressing loops, and it reduces cache-line overwrites during indirect addressing loops, which is fundamental for shared memory parallelism.

Validation of the adjoint code

Once the numerical adjoint scheme coded, the validation made is to check if the Jacobian fluxes in ALE formulation are equal to the variation rate of the fluxes in ALE formulation. For that we perform a validation test. We compute on one hand : This example was already presented in Chapter 2. We study the movement of a wing profile NACA0012 induced by an imposed of rotational displacement.

Mesh

In a circular domain of radius 20 is put a NACA0012 of length 1. The entire mesh includes 3291 vertices.

Initial conditions

To generate similar conditions of take-off, the initial Mach number chosen is 0.4.

Moving mesh

The moving of the NACA0012 is imposed. In phase 1, the coordinates of the geometry do not moved. The flow around the aircraft is established. Let t 1 = 0.1 be the time of the beginning of the movement and θ max = 20 o the maximal angle reached at first time t θ max . The movement of the wing profile is defined in phase 2 by calculating at each time step its angle of rotation α = 0.5 * θ max * (1cos(π * (tt 1 )/(t θ maxt 0 ))). The coordinates of the new points are then

x new = x old * cos α + y old * sin α and y new = -x old * sin α + y old * cos α. A phase 3 follows where the NACA0012 stays in upward position before the descent computed in the same way as phase 2.

Cost function

The considered functional output here is the drag of the surface S of the aircraft.

FSI Case : 2D Shock-ball interaction in a tube

This example is a variation of the classic Sod shock-tube problem. We consider a ball in a shock tube.

Mesh and initial conditions

The dimensions of the tube are: [0, 1] × [0, 0.2]. At initial time, the gas is split in two states:

for x ≤ 0.08 the state is (ρ lef t , u lef t , p lef t ) = (1, 0, 4) whereas for x > 0.08 the state is (ρ right , u right , p right ) = (0, 0, 2.5). A ball of radius r = 0.05 is immersed in the gas, its center being in position (0.5, 0.1). The tube is fixed. The simulation is run until final time t f = 0.5.

This initial state can be produced by having a diaphragm in the middle of the tube. The gas to the left and right of the diaphragm is initially at rest. The pressure and density are discontinuous across the diaphragm.

Cost function

The cost fonction considered is the pressure on the surface S of the ball. The output functional of interest is the quadratic deviation from ambient pressure on S :

j(W ) = T 0 S 1 2 (p(t) -p air ) 2 dSdt.
The observation is the pressure on this circle

Numerical results

At t = 0, the diaphragm is broken. The gas is then left to evolve freely and the classic rarefaction wave, contact discontinuity and shock wave appear. After a while, the ball is touched by the shock and the contact discontinuity, creating complex patterns both in front of it and in its wake.

The primal state and corresponding adjoint state at different time steps are shown in Fig- The bump is extruded in a bump-like cylinder.

Initial conditions

The half-circular blast is located on (1.2, 0). In 3D, it is a half-cylinder blast of symmetrical axis (1.2, 0, z). At initial time, the state is (ρ blast , u blast , p blast ) = (10, 0, 25) whereas for the rest of the domain the state is (ρ domain , u domain , p domain ) = (1, 0, 2.5).

Moving mesh

The moving mesh function is :

A = 0.

Conclusion

The discrete adjoint system associated to unsteady inviscid Euler flows computed with an ALE formulation was presented for both explicit and implicit time integration schemes. Considering the adjoint state context rises new difficulties. Indeed, the unsteady adjoint state has to be computed on a moving mesh and the mesh has to be moved backward in time. We need to introduce mesh velocities in the ALE adjoint state formulation. Therefore, time properties for backward moving mesh simulation were detailed. The opposite mesh velocity of forward moving mesh computation is considered to make backward meshes consistent with forward meshes. Thus, we have to address DGCL-like condition for the adjoint state introducing mesh velocities in the ALE adjoint state formulation. And numerical Jacobians for Riemann solvers have been presented.

In order to handle connectivity-changes a choice has been made to preserve memory. Each mesh optimization in the forward computation is tagged in order to make an optimization of the mesh at the same iteration in the backward moving-mesh computation. A validation of the adjoint has been performed to check the consistance of the jacobian fluxes.

Finally numerical cases were presented in 2D for an imposed displacement or for a fluid-structure interaction case. A 3D case of validation with comparison between moving mesh case and not moving mesh case have also been presented.

Chapter 5

Goal-oriented mesh adaptation for moving mesh simulations

When dealing with CFD problems, mesh adaptation is interesting for its ability to approach the asymptotic convergence and to obtain an accurate prediction for complex flows at a lower cost. Anisotropic mesh adaptation method reduces the number of degrees of freedom required to reach a given solution accuracy, thus impact favorably the CPU time. Moreover, it reduces the numerical scheme dissipation by automatically taking into account the anisotropy of the physical phenomena inside the mesh. Two main approaches of mesh adaptations exist in the literature. We have already seen the feature-based mesh adaptation in Chapter 2.

Feature-based mesh adaptation is mainly deduced from an interpolation error estimate using the Hessian of the chosen sensor. It controls the interpolation error of the sensor over the whole computational domain. Such approach is easy to set-up and has a wide range of application, but it does not take into account the considered PDE used to solve the problem. If we want to focus on a particular output functional g, as the drag or the lift for example, we need to use a method of mesh optimisation which takes g into account. This is the role of the goal-oriented mesh adaptation. However, considering a simulation with moving geometries in a goal-oriented mesh adaptation context rises new difficulties. Indeed, the unsteady adjoint state has to be computed on a moving mesh and the mesh has to be moved backward in time when computing the unsteady adjoint state. For more explanations on these, see Chapter 4.

Goal-oriented methods result from a series of papers dealing with a posteriori estimates as [Becker 1996, Giles 2002[START_REF] Wintzer | [END_REF]]. But, extracting informations concerning a posteriori estimate is a difficult task. Starting from a priori estimates, Loseille et al. proposed a fully anisotropic goal-oriented mesh adaptation technique for steady problems [Loseille 2010b].

The authors combine goal-oriented rationale and application of Hessian-based analysis to truncation error. Then the study is extended to unsteady problems in [START_REF] Belme | [END_REF]]. The error does not take into account the temporal error but concentrate only on spatial error for unsteady simulations.

In this present thesis, we combine fully anisotropic goal-oriented method of [START_REF] Belme | [END_REF]] and the fixed point advances of [N. [START_REF][END_REF] for time-accurate mesh adaptation involving moving geometries.

We start this chapter with a formal description of the error analysis in its most general expression, then the application to unsteady compressible Euler flows is presented. In Section 5.1, we introduce the optimal adjoint-based metric definition and all its relative issues, then Section 5.2 is dedicated to the implementation of the algorithm mesh optimization. This chapter ends with Section 5.3 and numerical experiments for blast wave problems.

Optimal adjoint-based ALE metric

We are interested in a scalar functional. The functions W, W h ∈ L 2 (Ω) are assumed to be the solution of the nonlinear continuous and discrete unsteady Euler PDE :

Ψ(W) = 0 ; Ψ h (W h ) = 0
where :

Ψ(W ) = W t + ∇ • F ale (W ) = 0.
(5.1)

Functional output of interest

Let j be a smooth linear functional applied to W into the scalar number : j(W ) = (g, W ) L 2 (Ω) ; j(W h ) = (g, W h ) L 2 (Ω) .

where (g, W ) holds for the following rather most general functional output formulation :

(g, W ) = T 0 Ω (g Ω , W )dΩ + Γ (g Γ , W )dΓ dt
The function g is a given L 2 (Ω) function. And we are interested in minimizing the approximation error committed on the evaluation δj of j : δj = j(W )j(W h )

Continuous adjoint state

Let us introduce the continuous and discrete systems in their variational formulation :

Find W ∈ V | ∀ϕ ∈ V, (Ψ(W ), ϕ) = 0, Find W h ∈ V h | ∀ϕ h ∈ V h , (Ψ h (W h ), ϕ h ) = 0.
The continuous and discrete adjoint systems in their variational formuation are:

Find W * ∈ V | ∀ψ ∈ V, ∂Ψ ∂W (W )ψ, W * = (g, ψ),

Find W * h ∈ V h | ∀ψ h ∈ V h , ∂Ψ h ∂W (W h )ψ h , W * h = (g, ψ h ).
An integration by parts gives :

∂Ψ ∂W (W )ψ, W * = Ω(t)
(ψ(0)W * (0)ψ(T )W * (T )) dΩ

+ T 0 Ω(t)
W * (T ) = g T .
The adjoint Euler equations is also a system of advection equations, where the temporal integration goes backward.

Formal error analysis

The idea is now to compute the difference of variational residual for a discrete test function:

δj = j(W ) -j(W h ) = (g, W ) -(g, W h ), (Ψ h (W ), ϕ h ) -(Ψ h (W h ), ϕ h ) = (Ψ h (W ) -Ψ(W ), ϕ h ) .
Then assuming that W * , Π h W * and W * h and their gradients are close to each other: δj ≈ (Ψ h (W ) -Ψ(W ), W * ) ≈ (Ψ h (W ) -Ψ(W ), Π h W * ) .

(5.4)

The term Ψ h (W ) -Ψ(W ) is an a posteriori local error which we now evaluate.

Local error analysis

We replace in Estimation (5.4) operators Ψ and Ψ h by their expressions given by Relation (5.1).

We also discard the error committed when imposing the initial condition. We finally get the following simplified error model:

δj ≈ n=nmax n=1 1 t n+1 -t n Ω t n+1 Π h W * ,n+1 (Π h W n+1 -W n+1 ) dΩ - n=nmax n=1 1 t n+1 -t n Ω t n Π h W * ,n (Π h W n -W n ) dΩ + n=nmax n=1 Ω t n Π h W * ,n ∇ • (Π h F ale (W n ) -F ale (W n )) dΩ - n=nmax n=1 ∂Ω t n Π h W * ,n (Π h Fale (W n ) -Fale (W n ))
.n dΓ.

(5.5)

Integrating by parts leads to:

δj ≈ n=nmax n=1 1 t n+1 -t n Ω t n+1 Π h W * ,n+1 (Π h W n+1 -W n+1 ) dΩ - n=nmax n=1 1 t n+1 -t n Ω t n Π h W * ,n (Π h W n -W n ) dΩ - n=nmax n=1 Ω t n ∇Π h W * ,n • (Π h F ale (W n ) -F ale (W n )) dΩ - n=nmax n=1 ∂Ω t n Π h W * ,n (Π h Fale (W n ) -Fale (W n ))
.n dΓ.

(5.6)

with Fale = Fale -F ale . We observe that this estimate of δj is expressed in terms of interpolation errors of the Euler fluxes and of the time derivative weighted by continuous functions Π h W * ≈ W * and ∇Π h W * ≈ ∇W * . The integrands in Error Estimation (5.6) contain positive and negative parts which can compensate for some particular meshes. In our strategy, we prefer not to rely on these parasitic effects and to slightly over-estimate the error. To this end, all integrands are bounded by their absolute values:

δj ≤ n=nmax n=1 1 t n+1 -t n Ω t n+1 |Π h W * ,n+1 ||Π h W n+1 -W n+1 | dΩ + n=nmax n=1 1 t n+1 -t n Ω t n |Π h W * ,n ||Π h W n -W n | dΩ + n=nmax n=1 Ω t n |∇Π h W * ,n | |Π h F ale (W n ) -F ale (W n )| dΩ + n=nmax n=1 ∂Ω t n |Π h W * ,n | |(Π h Fale (W n ) -Fale (W n ))
.n| dΓ.

(5.7)

Continuous error model

Working in the continuous framework enables to write Estimate (5.7) in a spatially-continuous form, in which the Π h are discarded, and in which the interpolation error Id -Π h is replaced by its continuous Idπ M . Then, we are interested in minimizing the following error functional: (5.8)

E(M) = n=nmax n=1 1 t n+1 -t n Ω t n+1 |W * ,n+1 ||π M W n+1 -W n+1 | dΩ + n=nmax n=1 1 t n+1 -t n
We observe that the fourth term introduces a dependency of the error with respect to the boundary surface mesh. In the present paper, we discard this term and refer to [Loseille 2010a] for a discussion on its influence. The first term can be transformed as follows without introducing a large error:

Ω t n+1 |W * ,n+1 ||π M W n+1 -W n+1 | dΩ = Ω t n |J n+1 n | -1 | Ŵ * ,n+1 ||π M Ŵ n+1 -Ŵ n+1 | dΩ where |J n+1 n
| is the determinant of the transformation from Ω n to Ω n+1 , J n+1 n = det ∂φ(., t n+1 ) • φ -1 (., t n ) ∂x (5.9) and Ŵ * ,n+1 resp. Ŵ n+1 the functions of Ω n obtained by reverse transportation from Ω n+1 :

∀x ∈ Ω t n Ŵ n+1 (x) = W n+1 φ(φ -1 (x, t n ), t n+1 , ∀x ∈ Ω t n Ŵ * ,n+1 (x) = W * ,n+1 φ(φ -1 (x, t n ), t n+1 .

(5.10)

Then, introducing the continuous interpolation error, we can write the simplified error model as follows: is defined on Ω t n . Here, W * j denotes the j th component of the adjoint vector W * , H(F ale,i (W j )) the Hessian of the j th component of the vector F ale,i (W ), and H(W j,t ) the Hessian of the j th component of the time derivative of W .

E(M) = n=nmax n=1 Ω trace M -1 2 (x,
We observe that the first line of RHS is a time derivative, since : where:

∀x ∈ Ω t n Ŵ n (x) = W n φ(φ -1 (x, t n ), t n , ∀x ∈ Ω t n Ŵ * ,n (x) = W * ,n φ(φ -1 (x, t n ), t n .

(5.13)

5.2 Implementation of goal-oriented mesh adaptation for unsteady problems involving moving geometries

Choice of the goal-oriented metric

The optimal adapted meshes for each sub-interval are generated according to the previous analysis. In this work, the following particular choices have been made :

• the Hessian metric for sub-interval i is based on a control of the temporal error in L 1 norm

• function τ : t → τ (t) is constant and equal to 1

• all sub-intervals have the same time length ∆t.

The optimal goal-oriented metric simplifies to : Remark 5. We notice that we obtain a similar expression of the optimal metric to those presented in Chapter 1.

l M i,ALE

Memory management

The main difference with the feature-based fixed-point mesh adaptation lies in the backward computation of the unsteady adjoint state which requires to know the solution on the whole simulation interval [0,T].

As seen in Chapter 4, the memory effort can be reduced by out-of-core storage of checkpoints as shown in the picture below. First the state-simulation is performed to store checkpoints.

Second, when computing backward the adjoint, we first recompute all states from the checkpoint and store them in memory and then we compute the unsteady adjoint until the checkpoint physical time. This method implies a recomputing effort of the state W .

As the fixed-point mesh adaptation algorithm split the time frame into n adap sub-intervals for mesh adaptation, this sub-intervals are used as checkpoints for the adjoint computation.

If the computation of the adjoint states need to much memory storage, the other strategy consists in storing solution states in memory only each m solver iterations. When the unsteady adjoint is solved, solution states between two savings are linearly interpolated. So the process is the following : at first, no interpolation is done, then if the system complains about a lack of memory, one solution out of two is deleted and interpolations are done and so on : if the memory capacity is still insufficient, only one solution out of four is kept and interpolation is made etc.

This method leads inevitably to a loss of accuracy for the unsteady adjoint computation.

Solve state once to get checkpoints Ψ(W ) = 0 Ψ * (W, W * ) = 0 Ψ(W ) = 0

Solve state and backward adjoint state from checkpoints

Algorithm

The Algorithm 11 is complicated by the fact that the unsteady adjoint solver must be computed backward in time (eg. step 2) after computing the solution over the simulation time frame.

Algorithm 11 Goal-Oriented Mesh Adaptation for Unsteady Flows

Input : Initial mesh and solution (H 0 , S 0 0 ) and set targeted space-time complexity N st # Fixed-point loop to converge the global space-time goal-oriented mesh adaptation problem For j = 1, n ptf x # Adaptive loop to compute forward the solution state in time on time frame [0, T ] 1. For i = 1, n adap (a) S j 0,i = Interpolate conservatively next sub-interval initial solution from (H j i-1 , S j i-1 , H j i ); (b) S j i = Compute solution on sub-interval from pair (H j i , S j 0,i ); = Generate all sub-interval adapted meshes ({H j i , M j i } n adap i=1 ); EndFor

Numerical test cases

We present in this section some preliminary results of our goal-oriented mesh adaptation algorithm for an unsteady simulation with moving domain.

Numerical test case : a Lohner blast-like on a 2D Bump

The test case presented is a Löhner semi-circular blast-like simulation in a domain [0, 5]×[0, 0.5] with a bump centered in (1.875, 0) of diameter 0.4. The circular blast is located on (1.2, 0). At initial time, the state is (ρ blast , u blast , p blast ) = (10, 0, 25) whereas for the rest of the domain the state is (ρ domain , u domain , p domain ) = (1, 0, 2.5).

Moving mesh

The moving mesh function is : The observation is on this area In goal-oriented mesh adaptation the mesh is adapted only on the flow interesting for the functional output while in hess adaptation, the remeshing is not quite dense.

A = 0.
Figure 5.3: Goal-oriented mesh adaptation

Conclusion

This chapter addresses goal-oriented mesh adaptation for time dependent and dynamic meshes simulations. In the first section, the error analysis is widely detailed. In the second section, the choice of metric optimization is presented. This concept is used to extend the previous unsteady goal-oriented mesh adaptation of [Belme 2011] to moving geometries.

The algorithm follows the work of [Olivier 2011a] and [Barral 2015] to generate meshes that are adapted once moved with a certain displacement. Analytic examples validated this metric, and it was used within a new error analysis, to derive an ALE metric for adaptive sub-intervals.

The regular time-dependent algorithm was modified to use this metric, and coupled with the moving mesh algorithm and example is done on a 2D analytic example.

The main contributions of this thesis with respect to pre-existing work in the GAMMA group is:

• An update strategy was proposed for the metric at mesh optimization times, to avoid spoiling the adaptation with optimizations.

• A study of the error on analytic case for the ALE goal-oriented metric that confirmed its adapted characteristic.

• A new version of the global fixed-point unsteady mesh adaptation algorithm was used.

Several perspectives arise from this work, notably:

• To improve the adaptation algorithm, the next step is to control the temporal error. This would in particular lead to automatic adaptation of the size of the sub-intervals. This appears to be all the more important with moving meshes, where the size of the mesh is not even constant throughout a sub-interval and time step. This improvement is necessary to deal with the ALE implicit scheme presented in Chapter 2.

• The convergence of this unsteady adaptive method needs to be studied, depending on the number of vertices and/or the size of the sub-intervals, as well as the influence of the initial solution on the final solution, and the ability of the algorithm to recover features of the solution lost due to coarse initial meshes.

Figure 1 . 1 :

 11 Figure 1.1: Unit balls associated with metric M = R Λ R T in two and three dimensions.

Figure 1 . 2 :

 12 Figure 1.2: Natural mapping Λ 1 2 R associated with metric M in two dimensions. It sends the unit ball of I 2 onto the unit ball of M.

Figure 1 . 3 :

 13 Figure 1.3: Several unit elements with respect to a metric tensor in 3D.

  Figure 1.4 gives an example of a Cartesian surface associated with a Riemannian metric space. This Riemannian metric space is pictured by drawing the unit ball of the metric at some points of the domain. A link with differential geometry is proposed in [Alauzet 2011a].

Figure 1 . 4 :

 14 Figure 1.4: Left, example of a Cartesian surface embedded in R 3 . Right, geometric visualization of a Riemannian metric space (M(x)) x∈[0,1]×[0,1] associated with this surface. At some points x of the domain, the unit ball of M(x) represented by ellipses is drawn.

Figure 1

 1 Figure 1.5: Metric-based mesh generation. Left, specified Riemannian metric space. Right, unit mesh in the prescribed Riemannian metric space which becomes adapted anisotropic in the Euclidean space.

Figure 1 . 6 :

 16 Figure 1.6: Representations of discrete (red) and continuous (blue) projectors for u concave and u convex for an unspecified mesh (first row) and for an unit mesh (second row).

  Figure 1.7: Notations and mapping in ALE framework.

Figure 2 . 1 :

 21 Figure 2.1: Top left, the swap operation in two dimensions. Top right, edge swap of type 3 → 2 and face swap 2 → 3. Bottom left, the five possible triangulations of the pseudo-polygon for a shell having five elements. Bottom right, an example of 5 → 6 edge swap. For all these figures, shells are in black, old edges are in red, new edges in green and the pseudo-polygon is in blue.

  (f) H k+1 = Move the mesh and update vertex speed H k , δt, v, v opt , a (g) Check mesh quality. If too distorted: solve new deformation problem or stop.

  (h) t = t + δt EndWhile EndWhile

Figure 2 . 2 :

 22 Figure 2.2: Left, 2D finite volume median cells and interfaces made up of facets and, right, hexahedron associated with vertex P i considered to build the finite volume cell of P i .

Figure 2

 2 Figure 2.6: Interface of a finite volume cell made up of facets (top) and boundary interface (bottom).

Figure 2 . 7 :

 27 Figure 2.7: Nosing up NACA0012 test case with explicit (left) and implicit schemes (right) Mach isolines at different time steps (t = 0, 30, 45).

Figure 2

 2 Figure 2.9: B61-12LEP dropping from a F16 (left). Corresponding 2D-computation (right)

Figure 2 .

 2 Figure 2.10: Initial 3D Mesh of the wing profile with attached bomb.

Figure 2 .

 2 Figure 2.11: 2D Mesh : naca0012 and sectional view of the bomb.

Figure 2 .

 2 Figure 2.12: Results of the bomb drop computation at different times (left). Corresponding zoom-in (right).

Figure 2 .

 2 Figure 2.14: F117 flight (left). Corresponding 3D-computation (right).

Figure 2 .

 2 Figure 2.15: Mach isolines of nosing up f117 test case in SSPRK2 explicit scheme (left)/ implicit BDF2 scheme with defect correction (right) at different time steps [ (view from the top) t = 0, 100, 200, 300 ].

Figure 2 .

 2 Figure 2.16: Mach isolines of nosing up f117 test case in SSPRK2 explicit scheme (left)/ implicit BDF2 scheme with defect correction (right) at different time steps [ (view from the top) t = 0, 100, 200, 300 ].

  j = Compute space-time complexity from all Hessian metrics {|H| j i }

Figure 3

 3 Figure 3.1: Fom left to right : Classical metric, ALE metric

In 2D .

 2D Figure 3.2: Initial sensor

Figure 3 . 3 :

 33 Figure 3.3: Initial Mesh is the uniform unstructured square mesh here is the representation of the displacement at t n+1

Figure 3

 3 Figure 3.4: Fom left to right : Hessian metric at t n+1 , ALE Hessian metric at t n , Moved and Updated ALE Hessian metric at t n+1

Figure 3 . 5 :

 35 Figure 3.5: Fom left to right : Updated Metric , ALE Metric

  for n = 1, ..., N.

Figure 4

 4 Figure 4.1. The case test used is the FSI Case 2D shock-ball interaction in a tube detailed in Section 4.8. The positions of the vertices are the same but the connectivity is different.

. 8 )

 8 Jacobian in ALE HLLC Riemann solverThe Figure4.2 illustrates the definition of our version of HLLC flux presented in section 2.2.2.

Figure 4

 4 Figure 4.2: Region of integration containing the Riemann fan in Godunov's method

Φ

  ale (a + da) -Φ ale (a), and on the other hand :∂Φ ale ∂W (a) • da.The validation for 7 iterations of adjoint of the case test used is the FSI Case 2D shock-ball interaction in a tube is illustrated in Figure4.3.

Figure 4

 4 Figure 4.3: Residual |Φ ale (a + da) -Φ ale (a) -∂Φ ale ∂W (a) • da|.

Figure 4

 4 Figure 4.4: Nosing up NACA0012 test case: (view from the top to the bottom) Pressure at different time steps( t = 165, 135, 120, 105, 90, 75, 60, 45).

Figure 4 . 5 :

 45 Figure 4.5: Ball in a shock tube test case: (view from the top to the bottom) Pressure

  at different time steps ( t = 0, 0.025, 0.05, 0.075, 0.1, 0.125, 0.150, 0.175, 0.2, 0.225 ) and corresponding adjoint state.

4. 8 . 3

 83 Numerical qualitative verification in 3D : Lohner blast-like on a BumpMeshes and initial conditions This 2D case is a Löhner semi-circular blast-like simulation in a domain [0, 5] × [0, 0.5] with a sinus-like bump of diameter 0.45 and center 1.775. The equation of this bump is similar to the Nasa 2D Bump-in-channel. The bumps extends from 1.55 to 2. Its maximum height is h = 04 (π(x -1.55)/0.45) , if 1.55 < x < 2 0 , else .It is then extruded into a 3D mesh of 261918 vertices and dimensions [0, 5] × [0, 0.5] × [0, 0.8].

  15 is the amplitude of the oscillation ω = π/0.7 is the frequence of the oscillationd(x, y, z) cos (2(x -2.3)π + π)) * A * sin (ωt) + x , if 2.29 < x < 2.8 A * sin(ωt) + x , if 2.8 < x < 3.8 0.5(1 + cos(2(x -2.3)π + π)) * A * sin(ωt) + x , if 3.8 < x < 4moved region is illustrated in Figure 4.6.

Figure 4 . 6 :

 46 Figure 4.6: Bump extruded mesh with colored moved region.

Figure 4 . 7 :Figure 4

 474 Figure 4.7: Surface area of interest.

  Ω t n |W * ,n ||π M W n -W n | dΩ + n=nmax n=1 Ω t n |∇W * ,n | |π M F ale (W n ) -F ale (W n )| dΩ + n=nmax n=1 ∂Ω t n |W * ,n | |(π M Fale (W n ) -Fale (W n )).n| dΓ.

  t) • H( Ŵ n+1 j )(x) -|J n n | -1 Ŵ * ,n j (x) • H( Ŵ n j )(x, t) (5.12)

EndFor#=

  Adaptive loop to compute backward the adjoint state in time on time frame [T, 0] 2. For i = n adap , 1 (a) (S * ) j i = Interpolate previous sub-interval final adjoint state from (H j i , H j i+1 , (S * 0 ) j i+1 ); (b) (S * 0 ) j i = Compute backward adjoint state on sub-interval from (H j i , S j i , (S * ) j i ); (c) |H go,L 1 | j i j = Compute space-time complexity from all goal-oriented Hessian-metrics {|H go,L 1 | j i } Compute all sub-interval unsteady metrics (C j , {|H go,L 1 | j i }

  15 is the amplitude of the oscillation ω = π/0.7 is the frequence of the oscillationd(x, y) cos (2(x -2.3)π + π)) * A * sin (ωt) + x , if 2.29 < x < 2.8 A * sin(ωt) + x , if 2.8 < x < 3.8 0.5(1 + cos(2(x -2.3)π + π)) * A * sin(ωt) + x , if 3.8 < x < 4is illustrated inFigure 5.1 (Top). It is represented by the red triangles. The blue ones are fixed.

Figure 5

 5 Figure 5.1: (Top)The region in blue is fixed while the red region moves. (Bottom) Illustration of the output functional observation area

Figure 5

 5 Figure 5.2: Hessian-based mesh adaptation

  solution with the flow solver from pair (H i , S 0

	i );
	If i = 1 break;
	2. (S *
	3. (M go,i ) = Compute metric (M go ) according to selected error estimate from (H i , S i , S * i );
	4. (H i+1 ) = Generate a new adapted mesh form pair (H i , M go,i );
	5. (S 0 i+1 ) = Interpolate new initial solution from (H i+1 , H i , S i );

i )= Compute adjoint solution with the flow solver from (H i , S i , j i );

  Algorithm 5 Connectivity-Change Moving Mesh Algorithm with Curved Trajectories Input: H 0 , S 0 , ∆t 0 Compute body vertex displacement from current translation speed v body , rotation speed θ body and acceleration a body for [t, t + ∆t/2] d(t + ∆t/2) = Solve elasticity system d body |∂Ω h (t + ∆t/2), ∆t/2 (b) d body |∂Ω h (t + ∆t) = Compute body vertex displacement from current translation speed v body , rotation speed θ body and acceleration a body for [t, t + ∆t] d(t + ∆t) = Solve elasticity system d body |∂Ω h (t + ∆t), ∆t (c) {v, a} = Deduce inner vertex speed and acceleration from both displacements {d(t + ∆t/2), d(t + ∆t)} (d) If predicted mesh motion is invalid then ∆t = ∆t/2 and goto 2.

	While (t < T end )
	1. ∆t = ∆t 0
	2. Solve mesh deformation: compute vertices trajectories
	(a) d body |∂Ω h (t + ∆t/2) =

Table 2 .

 2 Table 2.2. We only divided the computation time by 2 in this two-dimensional case.

	RK2	BDF2
	CPU Time 55min42s 29min24s

Explicit Implicit Figure 2.8: Nosing up NACA0012 test case with explicit (left) and implicit schemes (right) Mach isolines at different time steps (t = 60, 75, 90, 105, 147). 2: Nosing up NACA0012 test case. Final CPU time.

Table 2 .

 2 3: Nosing up F117 test case. Final CPU time.

  For a code computing Ψ h (W h )

	it can compute it as
		∂Ψ h ∂W h	U the forward mode (stands for forward accumulation)
	or	∂Ψ h ∂W h	T	V the reverse mode (stands for reverse accumulation).
		Automatic Differentiation generally takes as input a computer source program, plus a
	request for differentiation and then returns the differentiated source program that evalu-
	ates the required derivatives. Several authors use Automatic Differentiation to generate
	adjoint CFD codes as [M. Martinelli 2010]. For numerical applications, Becker and Ran-
	nacher [R. Becker 2001] and Giles and Pierce

  As the constant adjoint state is a solution of the adjoint Euler equations, if boundaries transmit the flux towards the backward outside as it comes, we find a purely geometrical relation inherent to the continuous problem. For any arbitrary closed volume C = C(t) of boundary ∂C(t), reverse time of the usual Geometric Conservation Law (GCL). From a geometric point of view, this relation states that the algebraic variation of the volume of C between two instants in backward equals the algebraic volume swept by its boundary.

	Relation (4.5) is integrated into:
		t+δt
	-|C(t + δt)| + |C(t)| =	t		∂C(t)	(w • n) dsdt, with t and t + δt ∈ [0, T ] ,	(4.6)
	which is a
	W * t -	∂F ale ∂W	*	∇W * dx -	∂C(t)	∂F ale ∂W	*	W * .n ds = 0	(4.4)
	If Equation (4.4) is written for a constant adjoint state, we obtain :
			-	d (|C(t)|) dt	-	∂C(t)	(w • n) ds = 0 .	(4.5)

  The local drag coefficient C d leads to the global drag coefficient for a given body :

	Primal State		Adjoint State
	S	1 S ref	C d (x, t)dx
	Numerical Results
	Views of the two states : the primal state and the dual state for the considered functional
	output are displayed in Figure 4.4.

  Consequently, the continuous adjoint state W * must be such that:

	ψ -W * t -	∂F ale ∂W	*	∇W * dΩ dt
	+ -W * T 0 Γ(t) ψ ∂F ale ∂W ∂W t -∂F ale

* W * .n -∂ Fale ∂W * W * .n dΓ dt . * ∇W * = g Ω(t) in Ω(t),

(5.2)

with the associated adjoint boundary conditions:

∂F ale ∂W * W * .n -∂ Fale ∂W * W * .n = g Γ(t) on Γ(t),

(5.3)

and the final adjoint state condition:

  t n ) H ALE go (x, t n ) M -1 2 (x, t n ) dΩ dt with H ALE go (x, t n ) = 5 j=1 H GO ale,j (x, t n ), in which

	H ALE j,go (x, t n ) = + +	1 t n+1 -t n |J n+1 n 1 j (x) • H(W n | -1 Ŵ * ,n+1 j (x, t) • H( Ŵ n+1 j j )(x, t) t n+1 -t n W * ,n ∂W * ,n j ∂x (x) • H(F ale,1 (W n j ))(x) + ∂W * ,n )(x) j ∂y (x) • H(F ale,2 (W n j ))(x)
	+	∂W * ,n j ∂z	(x) • H(F ale,3 (W n j ))(x)
			(5.11)

Note that the closer to 0.5 ν is, the harder to solve the system is.

The moment of inertia relative to an axis of direction a passing through G where a is an arbitrary unit vector is given by: J (Ga) = a T JGa

Remerciements

Conclusions and

Perspectives Conclusions

This thesis was carried out at Inria within the teams GAMMA3 (Inria Saclay) and ECUADOR (Inria Sophia-Antipolis). It has been part of the continuous works on mesh adaptation of these two teams and contributes to the developpement of the numerical in-house fluid solver Wolf.

This thesis has presented several novelties regarding the extension of both hessian-based and goal-oriented metric-based anisoltropic mesh adaptation to unsteady simulations for moving computational domains.

The problematics raised by mesh adaptation, mesh movement and adjoint method have been studied and analyzed at all the levels of the numerical resolution : the primal and dual solver phases, the forward and backward moving mesh algorithms and the mesh adaptation.

The work presented in this thesis follows from [Olivier 2011a, Belme 2011, Barral 2015] and presented updates for most stages of the resolution loop. Aside from the important theorical and technical prerequisites required for this work, the challenge was triple. To compute mesh adaptation for moving computational domains, the moving mesh must be performed with precision, then intelligent solutions must be performed to compute the adjoint solver backward in ALE formulation on a moving mesh and make it consistent with the flowfield state. And the last difficulty is to preserve the mesh quality after remeshing. For that reason, the mesh adaptation algorithm is divided into sub-intervals and coupled with an optimization algortihm of connectivity changes.

The main contribution on this thesis concerns the adjoint solver and the goal-oriented mesh adaptation. A new space-time goal-oriented error estimator has been designed to handle moving mesh simulations. It has been derived from the unsteady Euler equation in the ALE formulation. This ALE goal-oriented metric takes the movement of the mesh into account in the error estimate. Moreover the consideration of the mesh optimization made after the remeshing and during the connectivity change loop ensures the accuracy and the mesh quality of the ALE-mesh adaptation both for hessian-based mesh adaptation presented in Chapter 3 and for goal-oriented mesh adaptation in Chapter 5. The ALE solver has been adapted to handle this new backward moving mesh framework. In fact, this work also presents an implementation of adjoint solver for time dependent simulations with moving meshes in Chapter 4. This solver allows in particular mesh connectivity changes. Therefore, the backward changing-connectivity ALE scheme has been developed and validated in two dimensions and introduced in three dimensions. Finally and in the perspective of a space-time mesh adaptation, implicit solvers in ALE formulation for moving geometries problems have been developed and implemented in three-dimensional moving mesh simulations as illustrated in Chapter 2.

It is also already effective in terms of CPU time as compared to explicit time integration solvers.

Perspectives

A wide variety of perspectives emerge from this work. We list there below :

• The first natural one is the validation of 3D adjoint solver in ALE formulation and computation for more complex simulations. Then extend the goal-oriented mesh adaptation for 3D simulations.

• Since only rigid bodies are considered in this thesis, one perspective could be to consider deformable bodies.

• In the field of mesh deformation, other kinds of moving mesh equations can be investigated.

A premice of this investigation can be find in [Barral 2015].

• Future work will certainly focus on improving the space-time error analysis with moving meshes, notably by introducing a better control of the temporal error. This would notably result in an automatic adaptation of the sub-intervals sizes of the fixed-point algorithm.

Combined with the implicit time-integration solver in ALE formulation presented in this thesis could improve the simulation time significantly.

• In the future, extending this work to goal-oriented mesh adaptation for Navier-Stokes equations will lead to a greatest interest for the method to industrials. Already, a lot of ALE schemes are available for viscous simulations [START_REF] Hay | High-order temporal accuracy for 3D finite-element ALE flow simulations[END_REF]], which facilitate the task, together with many works on adaptation for viscous flows [Menier 2015, Frazza 2018].

• Time consumption is a huge limitation in the size and the complexity of the problems considered. Mesh adaptation appears as an answer to this problem, by helping to reduce greatly the size of the meshes while preserving the solution accuracy. However, it is only really efficient if it is coupled to an efficient parallelization of the codes. In this thesis, only shared memory parallelization has been considered. In the medium term, it seems that considering other kinds of parallelization will be unavoidable: reflections on parallelization of meshing software on distributed memory were started in [Loseille 2013], and the parallelization of the fluid solver on GPU was started in [Barral 2015]. Both approaches may have to be combined in order to take advantage of recent hybrid computation clusters.

• The FSI aspect can also be improved, by replacing the 6-DOF model by continuous elastic materials. The loose coupling could be replaced by more sophisticated strong coupling schemes. In this way, complex problems from aeroelasticity to biology could be run, in the interest of both industries and researchers.

• Another promising research perspective for the future is the extension of metric-based mesh adaptation to curved meshes. In fact, high-order elements have been used for a long while, mostly in Finite Element Methods, but only by increasing the order of the approximation polynomials. It has not been followed by meshing software, who have been generating flat pseudo high-order elements until very recently, i.e. standard P 1 elements with more nodes on the edges. The generation of really high-order meshes is only in its earliest infancy. Lately, progress has been made on the generation of P 2 curved meshes [Johnen 2013[START_REF] Abgrall | [END_REF]].

However, a smart use of curved meshes seems to be promising, notably by enabling a much more accurate representation of curved boundaries.