The three-dimensional simulation of physical phenomena involving moving bodies undergoing large displacements represents a real challenge. These simulations, combine the issues related to instability, meshing and Fluid-Structure Interaction, are generally difficult to perform and very costly in terms of computation time. We present in this research report an overview on moving mesh methods used in the literature. Several mesh deformation methods are proposed, local remeshing operations are explained and finally different approaches that uses PDE are detailed.

Introduction

In this project we are interested in specic metal manufacturing processes. The studied processes involve several heat transfer mechanisms. Studying these heat transfer mechanisms is very important since they condition the temperature distribution of the nal part and its quality. The numerical simulation has proved to be a powerful tool for understanding and controlling these phenomena.

The software THOST is dedicated to the simulation of heat transfer phenomena occurring in furnaces and quenching tanks [START_REF] Hachem | Modeling of heat transfer and turbulent ows inside industrial furnaces[END_REF][START_REF] Hachem | On the stabilized nite element method for steady convection-dominated problems with anisotropic mesh adaptation[END_REF]. In order to make the simulation of complex processes more realistic, THOST software has to be able to simulate the overall process as a single problem and to control both heating and cooling as a unique continuous and integrated process. This is the objective of HECO project: the control of the history of industrial parts along their heat treatment, from HEating to COoling.

The variety of the industrial equipment (turbines, fans,..) has to be taken into consideration but also the displacement of the parts from one equipment to another.

The three-dimensional simulation of these physical phenomena involving moving bodies undergoing large displacements represents a real challenge. These simulations, which combine the diculties related to instability [START_REF] Hachem | A three-eld stabilized nite element method for uid-structure interaction: elastic solid and rigid body limit[END_REF][START_REF] Jannoun | Anisotropic meshing with time-stepping control for unsteady convection-dominated problems[END_REF], meshing and Fluid-Structure Interaction, are generally dicult to perform and very costly in terms of computation time.

In such applications, as the body moves, the mesh may require adaptation at each time step to ensure good representation of the uid-solid interface or the existing mesh has to be allowed to deform to track the moving geometry. The former solution is computationally expensive, especially for 3D problems but also requires an interpolation step to map the solution into the newly generated mesh. This interpolation step is itself time-consuming but also introduces additional errors. The latter option introduces the concept of a moving and deforming mesh.

In the case of a deforming mesh, the mesh deformation has to be accomplished eciently and reliably in order to robustly perform the FSI simulations. Solver stability and accuracy is a primary concern in mesh moving algorithms. As the mesh deforms, elements can become inverted or highly skewed which can result in solver stability issues and errors degrading the accuracy of the model. For this, the algorithm used for moving the mesh has to be ecient and robust and able to handle large deformations. Usually, moving mesh algorithms are coupled to some mesh optimization operations that aim at maintaining the best possible the mesh quality while deforming it using either several local remeshing operations such as vertex addition, vertex collapsing, connectivity change or some mesh smoothing methods based only on vertex displacements. In addition to that, the eciency of the computations in terms of CPU time represents also a great challenge.

The aim of this thesis is to develop a robust moving mesh algorithm able to handle large displacements and ecient in terms of CPU time. In our context, a nite elements scheme and unstructured meshes are used.

A brief bibliographic overview of the moving mesh methods used in the literature will be given.

An outline of this report is given as follows: In section.1 the dierent mesh deformation methods used in the literature will be given. In section.2 some local remeshing operations will be explained and nally in section.3 the dierent approaches used for solving the governing equations are presented.

1

Mesh deformation methods

In moving mesh simulations, the whole mesh has to be allowed to deform in order to follow the displacement of the geometry while preserving the validity of the computational mesh. The whole problem can be formulated as follows: knowing the displacement of the vertices located on the moving body, what displacement must be applied to the rest of the domain vertices in order to respect the above description.

In general, a moving mesh algorithm has to be ecient in terms of CPU time as it is called at each solver time step. It has also to be able to keep a valid mesh so that remeshing operations remain very occasional in order to avoid the induced errors due to the interpolation steps required after any remeshing. It has to be able to preserve a good quality of the mesh so that the solution accuracy is not degraded. The main mesh deformation methods usually used in the literature will be given in this section.

Mechanical Analogy

In the case of the mechanical analogy based methods, two alternatives are generally considered: the spring analogy methods and linear elasticity methods.

Spring analogy method

The method views the uid mesh as a network of springs. In the spring analogy rst presented in [START_REF] Batina | Unsteady Euler airfoil solutions using unstructured dynamic meshes[END_REF], a ctitious tension/compression spring is attached to each edge connecting two vertices i and j, the stiness of the spring k ij is chosen to be inversely proportional to the length of the supporting edge l ij :

K ij ∝ 1 l ij
The above denition of K ij is driven by the fact that if during the mesh movement two vertices tend to get closer, the spring attached to the edge they belong to becomes stier and therefore prevents them from colliding.

This method often results in invalid meshes due to mesh lines crossovers; this is due to the fact that the tension/compression springs prevent two vertices from colliding but does not prevent a vertex from crossing an edge that faces it. This method has been improved by Farhat in [START_REF] Farhat | Torsional springs for two-dimensional dynamic unstructured uid meshes[END_REF] that introduced torsional springs at the mesh vertices to prevent neighboring triangles from interpenetrating each other; for a triangle whose vertices are designed by i, j and k, let θ ijk i denote the angle between two edges ij and ik of this triangle, the stiness of a torsional spring is given by:

C ijk i = 1 sin 2 θ ijk i = l 2 ki l 2 kj 4A ijk 2
where, A ijk is the area of the considered triangle, l ki is as before the length of the edge connecting k and i, and l kj the length of the edge connecting k and j.

The above stiness coecient is meant to prevent both vertex-to-vertex and vertex-to-edge collisions. Actually, whether vertex k moves towards vertex i or j, or whether it moves towards edge i -j, this leads the area A ijk of the triangle to decrease, which in turn results in the increase of the

stiness coecient C ijk i .
This method despite of the attempts to improve it has shown law robustness and tends to deteriorate the quality of the mesh when large displacements are considered [START_REF] Alauzet | A New Changing-Topology ALE Scheme for Moving Mesh Unsteady Simulations[END_REF].

Linear-elasticity method

In the case of the linear elasticity analogy [START_REF] Baker | Adaptive modication for time evolving meshes[END_REF], the motion of the internal nodes is governed by the equations of elasticity.

the movement of the vertices is obtained by solving an elasticity-like equation:

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

σ(ε) = λtrace(ε)I d + 2µε, and ε(σ) = 1+ν E σ -ν E trace(σ)I d
As boundary condition, the motion of the nodes at the interfaces is specied to match the normal velocity of the uid at the interface.

One advantage of this method is that it oers the possibility to introduce a special treatment based on element size, actually, it is possible to adapt the local stiness according to the size of each element so that smaller elements are stiened more than the large ones [START_REF] Tezduyar | Incompressible ow computations with stabilized bilinear and linear equal-order-interpolation velocity-pressure elements[END_REF]. This can be realized by modifying the way the Jacobian transformation from the element domain to the physical domain is accounted for in the FEM matrix assembly [START_REF] Stein | Mesh moving techniques for uid-structure interactions with large displacements[END_REF].

Classically, the P 1 Finite Element formulation of the linear elasticity matrix leads to the evaluation of quantities of the form:

´K s ∂ϕ J ∂x k ∂ϕ I ∂x l dx = s 36|K| ( _ η J ) k ( _ η I ) l
where s is either λ, µ or λ + 2µ and

_ η I = (( _ η I ) x , ( _ η I ) y , ( _ η I ) z )
is the inward non-normalized normal opposite to vertex i in tetrahedron K .In [START_REF] Stein | Mesh moving techniques for uid-structure interactions with large displacements[END_REF], the above quatity is replaced by :

´K s ∂ϕ J ∂x k ∂ϕ I ∂x l | K| |K| χ dx = | K| |K| χ s 36|K| ( _ η J ) k ( _ η I ) l
where K is the reference element. This method stiens each element by a factor proportional to |K| -χ and χ > 0 is the stiening power that determines the degree by which the smaller elements are rendered stier than the larger ones.

But, this method can only handle small incremental deformations and may require much subiteration to complete the full deformation of the mechanical time step. It is usually coupled to regular optimization phases to maintain a quality compatible with the desired accuracy.

Both spring analogy and linear elasticity methods need mesh connectivity information to construct a system of equations of the size of the number of ow nodes n i in order to calculate the displacement of the mesh. Solving this n i × n i matrix repeatedly is computationally intensive.

For the case of the linear elasticity method, in [START_REF] Alauzet | A changing-topology moving mesh technique for large displacements[END_REF], the author tried to nd a way to reduce the CPU time dedicated to this method. Actually, instead of solving the elasticity system at each solver time step, its solution is done once for a large time step ∆t max . After solving the elasticity system, the mesh motion is analyzed. If the mesh is not valid, then the elasticity time step is reduced until validity is achieved: ∆t = ∆tmax 2 n . Some results obtained by this method are shown in the Figure1.

Figure 1: Snapshots of a moving rotating cube using the method in [START_REF] Alauzet | A changing-topology moving mesh technique for large displacements[END_REF] 1

.2 Adaptive MMPDE-based methods

Mesh adaptation plays an indispensable role in the ecient numerical simulation of many physical phenomena since it helps reducing the dedicated CPU time to these simulations while keeping high accuracy [START_REF] Mesri | Automatic coarsening of three dimensional anisotropic unstructured meshes for multigrid applications[END_REF][START_REF] Mesri | Gestion et contrôle des maillages non structurés anisotropes: applications en aérodynamique[END_REF][START_REF] Mesri | On optimal simplicial 3d meshes for minimizing the hessian-based errors[END_REF]. Some authors view the theory of moving meshes as part of an r-adaptive strategy for solving partial dierential equations with moving bodies.

The idea is to move the mesh points and concentrate them in regions where the solution has a special behavior, usually a rapid variation of either the solution or its derivatives [START_REF] Dervieux | Continuous mesh adaptation models for cfd[END_REF][START_REF] Mesri | Hierarchical adaptive multi-mesh partitioning algorithm on heterogeneous systems[END_REF]. The objective is to get the smallest error possible for the number N of mesh nodes. A vector or a scalar monitor function is used in order to control the size, shape and orientation of the elements of the mesh to be generated. The monitor function is usually designed to give an estimate of some measure of the solution error which is then equidistributed over each mesh cell.

The mesh adaptation is viewed as a mathematical equivalent of the determination of a coordinate transformation. The adaptive mesh is viewed as the image of a uniform reference mesh under a coordinate transformation (mapping) x = x(ξ) from a regular domain (the computational domain) Ω c = {ξ j , j = 1, ..., N v } to an irregularly-shaped domain (the physical domain)

Ω = {x j , j = 1, ..., N v }. J = dx
dξ is used to denote the Jacobian matrix of the mapping x(ξ).

Several techniques have been developed in order to generate moving adaptive meshes. Most of these approaches can be classied into two groups. The rst group is called the location-based method (or the coordinate-based method) because it controls directly the location of the mesh nodes.

The second group is called the velocity-based approach since it controls the time derivative of the mapping or the mesh velocity. In other words, the mesh equation is formulated for the mesh velocity and the mesh point location is obtained by integrating the velocity eld. In the following paragraphs, good representatives of each group will be given.

Velocity-based approaches:

1.2.1.1 The classical Lagrangian method In uid dynamics, Lagrangian coordinates constitute a moving coordinate system used to follow uid particles. If u(x, t) is the velocity of the uid, ξ the reference coordinate of a uid particle and x(ξ, t) the position of the particle at time t, then the particle and consequently the Lagrangian coordinates envolve with:

∂x ∂t = u (1)
An interesting property of Lagrangian coordinates is that convection terms are eliminated from the governing equations. However, a purely Lagrangian formulation is usually avoided especially in problems involving large deformations since they tend to produce highly skewed meshes resulting in the failure of the numerical calculation.

1.2.1.2 The geometric conservation law method The so-called geometric conservation law (GCL) was used for many years in the engineering community to develop cell-volume-preserving nite volume schemes. In [START_REF] John | Numerical solution of the one-dimensional hydrodynamic equations in an arbitrary time-dependent coordinate system[END_REF] for example, the GCL was used to eliminate oscillations and preserve physical conservation laws for solutions on moving meshes.

In [START_REF] Cao | A moving mesh method based on the geometric conservation law[END_REF][START_REF] Thomas | Geometric conservation law and its application to ow computations on moving grids[END_REF] it is used to transform the algebraic expression specifying the Jacobian into equivalent dierential equation which is the key for the moving mesh strategy. The GCL is given by the following denition:

For any coordinate transformation x(ξ, t) from Ω c to Ω, the mesh speed x t and the time derivative of the Jacobian matrix for the mapping x(ξ) from Ω c to Ω satisfy:

∇.x t = - 1 J DJ Dt (2)
Let's consider the determination of the mapping x(ξ, t) through its time derivative. For a given monitor function ρ(x, t) > 0, we want the cell area to be inversely proportional to ρ(x, t) at each time t. To do so, the mapping x(ξ, t) is required to satisfy:

∇.x t = - 1 ρ Dρ Dt (3)
By comparing equations ( 2) and (3) we have :

D Dt (ρJ) = 0 (4)
Which implies that:

ρJ = constant (5)
,

(3) can be written as:

∇.(ρx t ) + ∂ρ ∂t = 0 (6)
Using only this equation is not sucient to determine the mesh speed x t . Using the Helmholtz Decomposition Theorem for vectors: A continuous and dierentiable vector eld can be resolved into the orthogonal sum of the gradient of a scalar eld and the curl of a vector eld a supplementary equation can be provided. Actually, x t can be calculated using its divergence through [START_REF] Neil | Design and application of a gradient-weighted moving nite element code ii: in two dimensions[END_REF] and its curl. For this x(ξ, t) is required to satisfy :

∇ × ω(x t -u) = 0 (7)
where ω and u are, respectively, a weight function and a background velocity eld to be specied. Dierent choices of ω lead to dierent curl conditions. Equations ( 2) and ( 7) can be formulated as a minimization of the functional:

I GCL = 1 2 ˆΩ |∇.(ρx t ) + ∂ρ ∂t | 2 + ρ ω 2 |∇ × ω(x t -u)| 2 dx (8)
The following boundary condition is considered:

x t .n = 0 (9)
Equations ( 8) and ( 9) constitute an elliptic system for the mesh velocity x t , its solution exists and is smooth. By satisfying the relation ( 5) the Jacobian of the mapping is determined by the monitor function ρ(x, t) and thus, the Jacobian stays positive and the mapping stays itself locally non-singular. However, this method usually tend to produce highly skewed meshes.

1.2.1.3

The moving nite element method The moving nite element method (MFE) developed in [START_REF] Miller | Moving nite elements. I[END_REF][START_REF] Miller | Moving nite elements. II[END_REF] is also a velocity-based method. Given a time dependent physical problem:

∂u ∂t = Lu
where L is a spatial dierential operator. Determining x t using the MFE method, consists in minimizing the L2-norm of the residual involving u and x t over the entire space:

I mf e x t , Du Dt = Ώ Du Dt -∇u.x t -Lu 2 ωdx
In the classical version of MFE [START_REF] Miller | Moving nite elements. I[END_REF][START_REF] Miller | Moving nite elements. II[END_REF] ω is taken equal to 1 and in the gradient weighted MFE (GWMFE) [START_REF] Neil | Design and application of a gradient-weighted moving nite element code i: in one dimension[END_REF][START_REF] Neil | Design and application of a gradient-weighted moving nite element code ii: in two dimensions[END_REF] ω is equal to 1 /(1+|∇u| 2 ).

One positive point of the MFE method is that the resulting mesh has the tendency to follow a path corresponding to the smallest L2-norm of the residual of the discrete equations. The problem is that the mesh resulting from only the minimization of this error can become degenerate or nearly degenerate and requires special regularizing terms (i.e., penalty functions). [START_REF] Winslow | Adaptive-mesh zoning by the equipotential method[END_REF] requires solving a non-linear Poisson-like equation in order to generate the mapping. Suppose that higher mesh concentrations are desired in regions with large values of a given function ω(x) > 0. With the approach of Winslow, the mapping is dened as the minimizer of :

Location-based approaches 1.2.2.1 Winslow's variable diusion functional Winslow's formulation

I W in [ξ] = Ώ 1 ω i (∇ξ i ) T (∇ξ i )dx 1.2.2.2
The functional of Brackbill and Saltzman Brackbill and Saltzman developed a global coordinate approach based on modifying the Winslow algorithm. In [START_REF] Jeremiah | Adaptive zoning for singular problems in two dimensions[END_REF] they tried to combine mesh smoothness, concentration and orthogonality. The functional associated with mesh concentration for example is given by the integral:

I BS [ξ] = Ώω(x)Jdx
One diculty that arises is that there is no theory to guide the relative weights of each integral and they are left to be specied by the user. In addition to that, the various terms are not dimensionally homogenous and must be rescaled according to each specic application. One more thing is that one cannot guarantee the existence of a solution to the previous equation (unless it is a pure smoothness functional).

Harmonic mapping

In [START_REF] Arkady | Adaptive grid generation from harmonic maps on riemannian manifolds[END_REF] Dvinsky developed a method for generating adaptive grids based on harmonic maps. Given an n × n symmetric positive denite matrix G(x) with g = det(G), the mapping is dened as a harmonic mapping between Ω ⊂ R n equipped with metric G and Ω c ⊂ R n equipped with the Euclidean metric. This is equivalent to nding the minimizer of:

I hrm [ξ] = ´√g Ω i (∇ξ i ) T G -1 (∇ξ i ) dx
In [START_REF] Li | Moving mesh methods in multiple dimensions based on harmonic maps[END_REF], the authors used a moving mesh strategy similar to Dvinsly's scheme except the fact that a nite elements approach was used rather than a nite dierence discretization and the fact that the harmonic mapping in this case was constructed by an iterative procedure.

1.2.2.4 Brackbill's direction control functional Based on Winslow method and the harmonic mapping method, Brackbill [START_REF] Brackbill | An adaptive grid with directional control[END_REF] developed a functional which takes into consideration both mesh concentration and mesh alignment. The functional can be written as follows:

I brb [ξ] = Ώ i (∇ξ i ) T G -1 i (∇ξ i ) dx
for user-dened functions G i . Depending on the desired properties of the mesh, they can be chosen in several ways. For example, in the case G i = 1 √ g G, the functional reduces to a harmonic map. If, for instance, the mesh concentration is of primary concern, the natural choice is G = ω(x)I , which leads to Winslow's functional. If the mesh is supposed to align with a certain prescribed

eld(v 1 , v 2 , v 3 ), then G can be chosen such that (∇ξ i ) T G -1 (∇ξ i ) = |v i × ∇ξ i | 2 .
1.2.2.5 Huang's approach The previous methods succeed in concentrating mesh nodes in regions where the monitor function is large. However, the precise relation between the mesh and the monitor function is unclear. One may observe an improvement of the solution accuracy using the previous types of mesh generation but it is not clear whether the obtained mesh is optimal in terms of minimization of the actual error or not. In order to deal with this problem, Huang in [START_REF] Huang | Variational mesh adaptation: isotropy and equidistribution[END_REF] formulates a functional based more directly on error distribution.

In [START_REF] Ngo | A Study on Moving Mesh Finite Element Solution of the Porous Medium Equation[END_REF] the mapping is formulated by minimizing an energy based on the so-called equidistribution

(1) and alignment (2) conditions [START_REF] Huang | Mathematical principles of anisotropic mesh adaptation[END_REF] and the mesh adaptation is controlled through a matrix-valued function.

Equidistribution condition:

|K| det(G K ) = σ h |Kc| |Ωc| , ∀K∈ Ω Alignment condition: 1 d trace((J)G -1 (J) -T ) = det((J)G -1 (J) -T ) 1 d , ∀K ∈ Ω
where, |K| is the volume of K, G is the metric tensor associated with the adaptive mesh, G K is the average of G over K, det(.) and trace(.) denote the determinant and trace of a matrix, respectively, |K c | is the volume of the element K c ∈Ω c corresponding to K, and

σ h = |K| K∈Ω det(G K ), |Ω c | = |K c | Kc∈Ωc
The meshes that closely satisfy these conditions can be obtained by minimizing the energy function:

I h = θ |K| K∈Ω det(G K )trace((J)G -1 K (J) -T ) dp 2 + (1 -2θ)d dp 2 |K| K∈Ω det(G K ) ( |Kc| |K| √ det(G K ) ) p
The above functional is used for generating an adaptive mesh for a specied monitor function.

When one is solving a time dependent problem, and in order to move the mesh in a smooth manner the MMPDE approach is used to dene the moving mesh equation as a gradient system of I h [START_REF] Huang | A high dimensional moving mesh strategy[END_REF][START_REF] Huang | Moving mesh strategy based on a gradient ow equation for two-dimensional problems[END_REF] dξ j dt = -

P j τ ∂I h ∂ξ j T , j = 1, ....., N v
where the derivative of I h with respect to ξ j , ∂I h /∂ξ j , is considered as a row vector, τ >0 is a parameter used to control the response time of the mesh movement to the change in the metric tensor.

In [START_REF] Tan | Adaptive moving mesh methods for two-dimensional resistive magnetohydrodynamic pde models[END_REF] a specic map is obtained by minimizing a mesh adaptation functional of the following form:

E(ξ, η) = 1 2 ´Ω(∇ξ T G -1 ∇ξ)
This functional is minimized by solving the corresponding Euler-Lagrange equation :

∇.(G -1 ∇ξ) = 0,
This system is solved using Gauss-Seidel iteration. The iteration is continued until there is no signicant change in calculating new grids from one iteration to the next.

Numerical comparison

In [START_REF] Cao | Approaches for generating moving adaptive meshes: location versus velocity[END_REF], the authors tried to present a short numerical comparison of four adaptive moving mesh methods : the GCL method, the winslow's approach, the harmonic map approach and Huang's approach.

The monitor function used is dened as follow:

ρ(x, y, t) = 1 + 10exp -50 x -1 2 -1 4 cos(2πt) 2 + y -1 2 -1 4 sin(2πt) 2 -1 10 2
This equation needs the generation of a moving adaptive mesh from the time-dependent monitor function.

The gures below show the adaptive meshes at dierent times for the four methods.

(a) The adaptive meshes obtained by the GCL method at t=0,0.25,0.5,0.75 [START_REF] Cao | Approaches for generating moving adaptive meshes: location versus velocity[END_REF] (b) The adaptive meshes obtained by Winslow's method at t=0,0.25,0.5,0.75 [START_REF] Cao | Approaches for generating moving adaptive meshes: location versus velocity[END_REF] (c) The adaptive meshes obtained by the harmonic map method at t=0,0.25,0.5,0.75 [START_REF] Cao | Approaches for generating moving adaptive meshes: location versus velocity[END_REF] (d) The adaptive meshes obtained by Huang's method at t=0,0.25,0.5,0.75 [START_REF] Cao | Approaches for generating moving adaptive meshes: location versus velocity[END_REF] Figure 2: Numerical Comparison For winslow's and Huang's methods, the mesh concentration follows closely the evolution of the function ρ. For the harmonic mapping, the cercle with the highest mesh concentration is slightly outside the cercle with the maximum ρ. Moreover, the mesh in the central part of the cercle is much coarser than in the outer region although ρ takes approximately the same valye in both places. This was also noticed in [START_REF] Cao | A study of monitor functions for twodimensional adaptive mesh generation[END_REF]. In the case of the GCL, the mesh becomes highly skewed and the mesh concentration slowly misplaced.

Level-set moving mesh partial dierential equation (LMPDE)

The Level-set deformation method [START_REF] Liao | Level-set-based deformation methods for adaptive grids[END_REF][START_REF] Ong | An hr moving mesh method for onedimensional time-dependent pdes[END_REF] moves the nodes with a proper velocity so that the nodal mapping has the desired Jacobian determinant. The grid in this method is represented in an implicit manner; the intersection points of the level-sets will form the grid nodes at each time.

For example in [START_REF] Ong | An hr moving mesh method for onedimensional time-dependent pdes[END_REF], for 2D simulations, the authors used two level-sets of the following form to represent the grid:

ψ 1 (x, y, 0) = x, ψ 2 (x, y, 0) = y,
The nodes of the mesh are after moved using the usual level-set evolution equation. The velocity vector used in the evolution equation is constructed by solving a Poisson equation determined by the monitor function.

The problem of this method is that it is able to produce only structured grids.

Optimization-Based Rezone strategy

The reference Jacobian optimization-based rezone strategy developed in [START_REF] Knupp | Reference Jacobian Optimization-Based Rezone Strategies for Arbitrary Lagrangian Eulerian Methods[END_REF] aims at producing a mesh of high quality while keeping it as close as possible to the Lagrangian mesh. The essential idea of the RJM method is the recognition that the Lagrangian solution contains sucient information about the ow to constrain our measure of mesh smoothness.

In [START_REF] Knupp | Reference Jacobian Optimization-Based Rezone Strategies for Arbitrary Lagrangian Eulerian Methods[END_REF] the author considered three grids. At the beginning of the computational cycle, the grid of the previous cycle is considered from the previous step and used as an initial grid. This grid is assumed to be of good quality which means that it is unfolded and also well adapted to the ow at this time step. The second grid is the Lagrangian grid; this grid is obtained by simply moving each node of the initial grid with the Lagrangian uid velocity. This Lagrangian grid is not necessarily optimal and may require rezoning . The rezoning strategy aims at improving the quality of the Lagrangian grid while preserving the Lagrangian grid as much as possible. By requiring the rezoned grid to remain as close as possible to the Lagrangian grid, employing locals remappers is justied.

Local remappers are computationally much more ecient than global remappers [START_REF] John | Accurate conservative remapping (rezoning) for arbitrary lagrangian-eulerian computations[END_REF]. One more reason to maintain the rezoned grid close to the Lagrangian one is that, the Lagrangian grid contains important information about the ow; it follows for example the interface between two materials.

The rezoning strategy consists of two components: A sequence of local optimizations and a global optimization. The sequence of local optimizations aims at dening a set of reference Jacobians which incorporate the denition of mesh quality at each grid point. Actually, at each node a local patch is formed from the adjacent cells of the Lagrangian grid and a local realization of the Winslow smoothness functional [START_REF] Winslow | Numerical solution of the quasilinear poisson equation in a nonuniform triangle mesh[END_REF] on this patch is constructed. The minimization of this functional with respect to the position of the central node denes its virtual location. By connecting this virtually moved node to its stationary neighbors, a reference Jacobian (RJM) that represents the best locally achievable geometric grid quality is dened.

A local realization of the Winslow function is used to determine the new virtual position of the central node. The optimized position is obtained by minimizing the following functional [START_REF] Jeremiah | Adaptive zoning for singular problems in two dimensions[END_REF] with two variables while keeping all other coordinates xed at their Lagragian positions.

F (x, y) = ´1 0 ´1 0 [(x ξ ) 2 +(y ξ ) 2 ]+[(xη) 2 +(yη) 2 ]
|J| dξdη where J =

x ξ x η y ξ y η is the Jacobian matrix of map x(ξ, η), y(ξ, η), and |J|is its determinant.

It can be noticed that this functional has a barrier function, which means that the value of the functional approaches ininity when the unknown map approaches any degenerate map where

|J| = 0.
The RJMs are constructed by connecting these virtual positions.

It is not possible to use directly the RJMs themselves to determine the new positions of the vertices. Actually, each vertex belongs to four cells: a simple rezone strategy applied individually to each cell will lead in general to four incompatible specications of the rezoned position of any vertex.

For this, a global optimization step is performed in order to resolve incompatibilities of the locally dened RJMs. Actually, the nal rezoned grid results from a minimization of a global objective function that measures the distance between the RJMs and the Jacobian of the rezoned grid (the solution of the global optimization step). This function does not guarantee that the rezoned grid is unfolded, for this, a barrier function that penalizes grids whose cells are inverted is included. This function is minimized using a line search procedure coupled with a conjugate gradient algorithm. The optimization procedure is iterative and every time that it results in a folded mesh, the mesh is rejected and the objective function is set to a large number.In [START_REF] Knupp | Reference Jacobian Optimization-Based Rezone Strategies for Arbitrary Lagrangian Eulerian Methods[END_REF] the method was used for structured grids but it can also be used for unstructured grids.

F (x(ξ, η), y(ξ, η)) = ´1 0 ´1 0 J-J ref 2 F |J| / | J ref | dξdη

Interpolation methods

Interpolation methods treat mesh deformation as a problem of interpolating deformations between the moving body and the boundaries of the domain of calculation. These methods are able to handle large deformations and arbitrary mesh topologies and can be easily implemented in parallel.

Radial Basis Function (RBF)

One main interpolating method is the Radial Basis Function (RBF).The RBF method [START_REF] Rendall | Ecient mesh motion using radial basis functions with data reduction algorithms[END_REF][START_REF] Nail | Fast radial basis function interpolation via preconditioned krylov iteration[END_REF]consists in interpolating the displacement of the moving body to all the ow mesh using an interpolation function:

s(x) = n b j=1 α j φ(||X -X b j ||) + p(X) (10)
where X b j = [x b j , y b j , z b j ] are the centers in which the values are known, in this case the boundary nodes of the moving body, p a polynomial, n b the number of boundary nodes and φ is a given basis function with respect to the Euclidean distance ||X||.

The coecients {α} and the polynomial p(X) are chosen to satisfy the n b tting conditions:

s(X b j ) = d b j
with d b containing the discrete known values of the displacements at the boundary nodes and the constraints :

n b j=1 α j q(X b j ) = 0
for all polynomials q with a degree less or equal than that of polynomial p.

The values of the coecients and the linear polynomial can be obtained by solving the system:

d b 0 = M b
i b j = φ( X b i -X bj i ).
Solving this system can be done in an iterative manner. RBF is an ecient method to achieve high quality mesh and to handle large mesh deformations caused by translations and rotations for both 2D and 3D meshes. However, the performance depends on the used RBF [START_REF] De Boer | Mesh deformation based on radial basis function interpolation[END_REF].

The main drawback of this method is that the n b × n b matrix is often dense and ill-conditioned and solving it especially for large meshes is expensive due to the large number of surface points. This method can be optimized using specialized preconditioners, however, they add more complexity to the implementation of the method. Some other methods were developed to vastly reduce the number of surface points used, such as the Greedy point selection used in [START_REF] Rendall | Reduced surface point selection options for ecient mesh deformation using radial basis functions[END_REF]. The method uses a chosen error function on the surface mesh to select a reduced subset of surface points; this subset contains a suciently small number of points so as to make the volume deformation fast.

Inverse Distance Weighting interpolation method (IDW)

One other interesting method is based on the Inverse Distance Weighting interpolation method (IDW) . The IDW [START_REF] Luke | A fast mesh deformation method using explicit interpolation[END_REF] is a weighted average interpolation technique for multivariate interpolation of scattered data points.

The interpolated value is an average of the known values at the data points weighted by the inverse of the distance to the unsampled point. For mesh deformation, that means that the inuence of boundary node displacements on the displacement of an internal mesh node is inversely proportional to the distance between the two points.

The moving boundaries nodal displacement eld is dened as s(x b i ), where x b i is the position vector of the boundary node i. Then, the displacement eld s in the volume mesh is then described through a weighted average of all moving body boundary nodes displacement eld as given by:

s(x) = m i=1 w i (x)s(x b i ) m i=1 w i (x) (11)
where m is the total number of boundary nodes and w i (x) is a two-exponent weighting function of the reciprocal distance

w i (x) = A i L ref x-x b i a + αL ref x-x b i b ,
where A i is the average area of all moving boundary faces containing the node i, this is used so that the mesh renement of a region does not increase its inuence in the interpolation method, L ref is the distance between the mesh centroid and the farthest point of the domain, a is the basic IDW exponent, α is the maximal distance that nodes near to the moving body have to a boundary, and b is the IDW exponent which rigidies nodes in a near-body region. α is given by : α = 5

L ref m max i=1 s(x b i ) -s , where s = m i=1 .s(x b i ) and A i = A i m j=1 A j
The weighting function is designed such that it preserves a rigid body motion of nodes close to the boundaries, in addition to ensuring a smooth deformation transition through the mesh.

In order to increase the robustness and lower computational cost, this method was coupled in [START_REF] Landry | Robust moving mesh algorithms for hybrid stretched meshes: Application to moving boundaries problems[END_REF] to the so-called Moving Submesh Approach (MSA) proposed in [START_REF] Lefrançois | A simple mesh deformation technique for uid-structure interaction based on a submesh approach[END_REF]. It consists in interpolating the displacement to the ne computational mesh from the solved motion of a coarse mesh which is obtained by the IDW method.

In contrast to RBF interpolation, IDW method is an explicit mesh deformation method which does not require solving a system of equations for deforming the volume mesh. In fact, it results in an algebraic expression for the nodes displacements as function of the boundary deformation.

This explicit evaluation reduces the computational costs and simplies the implementation and the parallelization of the mesh deformation routines and can handle translations, rotations and arbitrary mesh topologies.

2 Post-treatment using local remeshing operations While moving, the quality of the mesh decreases. Rather than generating a new mesh, some local remeshing operations can be performed in order to maintain a good mesh quality.

The standard mesh motion algorithms are able to move the mesh points. Then some kind of mesh smoothing can be applied to the mesh in order to maintain a reasonable mesh quality. Using this kind of operations, the topology of the mesh is not modied, which means that the underlying graph of the mesh remains unchanged. However, this approach fails for even simple motions such as large rigid body translation and rotation. Moreover, although the mesh topology may be preserved for simple mesh motions, such a procedure gives little control of the mesh size, hence, some cells may be compressed or stretched undesirably due to the mesh motion and smoothing, which results in a large error in the solution.

H-adaptation methods can be used to improve the quality of the mesh by modifying its topology [START_REF] Compère | A mesh adaptation framework for dealing with large deforming meshes[END_REF].

They consist in locally coarsening or rening the mesh by the inclusion or deletion of mesh nodes or in changing the connectivities between the existing nodes. But, since these operations change the topology of the mesh, once performed they have to be followed by an interpolation step.

If no attention is paid, local remeshing operations may invalidate the mesh. Consequently, adequate checks of the mesh quality have to be performed. In [START_REF] Alauzet | A New Changing-Topology ALE Scheme for Moving Mesh Unsteady Simulations[END_REF], after moving the mesh swap operations are performed to optimize the mesh.

Edge swap operation (Figure3a) plays an important role in improving the mesh quality and acts only on the connectivities of the mesh leaving the vertices' positions unchanged. In 2D, edge swapping is a simple operation which consists in ipping an edge shared by two triangles. In 2D, this operation changes neither the number of triangles nor the number of edges of the mesh. Face swapping (Figure3b) consists in suppressing face F of a tetrahedron K.

(a) edge swap operation [START_REF] Compere | A mesh adaptation framework for dealing with large deforming meshes[END_REF] (b) face swap operation [START_REF] Compere | A mesh adaptation framework for dealing with large deforming meshes[END_REF] In [START_REF] Dapogny | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF] the authors used, in addition to the swapping operations, edge collapse operation (Figure3d)

which is used to remove one vertex from the mesh. Collapsing an edge pq consists in merging its two end points to a single one. They also used edge split (Figure3c) that consists in splitting one edge pq by introducting a new point m in the mesh.

(c) edge split operation [START_REF] Compere | A mesh adaptation framework for dealing with large deforming meshes[END_REF] (d) edge collapse operation [START_REF] Compere | A mesh adaptation framework for dealing with large deforming meshes[END_REF] Usually, swap operators (edge swaps and face swaps) aim at improving locally the quality of the elements. Splits and collapses are there to control the mesh size: long edges are split whereas one of the two vertices of short edges is collapsed.

These operations are more ecient than the xed-topology optimization methods in maintaining good mesh quality but, once used, local interpolations steps have to be performed. These interpolations are time consuming and induce extra errors.

Resolution of equations in moving mesh framework

A fundamentally important consideration when developing a computer code for simulating timedependent physical or mechanical problems occurring in evolving domains is the choice of an appropriate kinematical description of the continuum. Actually, this choice determines the relationship between the deforming continuum and the mesh. In the case of FSI problems two main classes of methods are usually used depending on whether an implicit or explicit representation of the considered evolving domain is used.

As parts of the rst class, Eulerian methods use a xed mesh of a computational box. A Level-Set function is used according to which the evolving domain is dened in the sense of implicit function.

The so-called level set approach enables to represent the interface between two domains as a zero level of a smooth function [START_REF] Hachem | Stabilized nite element solution to handle complex heat and uid ows in industrial furnaces using the immersed volume method[END_REF]. In practice, a signed distance function is used to localize the interface.

Let W i1 , W i2 and G i represent, respectively, the rst domain, the second domain and the interface.

They verify:

Ω i1 ∪ Ω i2 = Ω and Ω i1 ∩ Ω i2 = Γ i
For each node of the computational domain W, the level set function a which is the signed distance from the interface reads:

α(x) =      > 0 x ∈ Ω i1 0 x ∈ Γ i < 0 x ∈ Ω i2
Thanks to this implicit representation, the domain evolution can easily be accounted using the following advection equation [START_REF] Khallou | High delity anisotropic adaptive variational multiscale method for multiphase ows with surface tension[END_REF]:

∂α ∂t + v.∇α = 0 α(t = 0, x) = α 0 (x)
In this case, any modication of the mesh that aims at following the movement of the moving body has to be followed by an interpolation step.

In our case, we use the so-called immersed Volume method . This approach is of a great exibility since it allows us to solve one set of equations in the whole domain. It consists in generating one unique volume mesh in which the dierent parts of the complex geometry are immersed and treated as dierent materials using Level-Set function and an anisotropic mesh adaptation at the interface between two materials. In order to ensure a continuous variation of the material properties across an interface, the so-called mixing laws are used [START_REF] Schmid | Unied formulation for modeling heat and uid ow in complex real industrial equipment[END_REF].

On the opposite, Arbitrary-Lagrangian-Eulerian (ALE) methods [START_REF] Donea | Arbitrary Lagrangian Eulerian Methods[END_REF] use an explicit discretization of the considered evolving domain. So, an exact and conformal mesh to the domain has to be generated.

In the ALE framework the nodes of the domain can move with the continuum and their velocity is directly integrated in the conservation equations. In fact, when the mesh is deformed while keeping its topology xed, the equations are solved in a fully ALE manner so no interpolation step is needed.

But, due to the xed-topology constraint imposed by the classical Arbitrary-Lagrangian-Eulerian formulation, once the connectivity of the mesh changes between two time steps, local or global interpolation operations have to be performed in order to get the solution on the new mesh. In order to get rid of these interpolation operations, some attempts to get a true changing topology ALE scheme have been done, but only in 2D [START_REF] Alauzet | A New Changing-Topology ALE Scheme for Moving Mesh Unsteady Simulations[END_REF].

In [START_REF] Dapogny | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF] the authors tried to combine the advantages of both methods; on the one hand they have considered a conformal mesh to the interface (so an explicit representation of the interface) when it comes to solving the mechanical problem and on the other hand, they considered an implicit representation of the evolving geometry when it comes to compute its evolution using the Level-set function described above.

Conclusion and perspectives

We tried in this report to give a brief bibliographic overview of the main existing moving mesh methods.

We started our study with the mechanical analogy approach which includes mainly the spring and the elastic methods. In both approaches, the stiness of the springs or elastic material is varied locally based on mesh cell size or proximity to the moving body. The choice of the parameters of these methods is left to the user and is usually based on ad hoc rezoning about specic quality objectives. These methods cannot handle all types of meshes and can have problems especially when high aspect ratio meshes are used. They may also require many sub-iterations to complete the full deformation of a mechanical time step which makes them time-consuming. These schemes are also typically dicult to implement.

We presented in the following paragraph the adaptive MMPDE-based methods that can also be divided into two groups: the velocity-based and the location-based methods. One main dierence between the two methods is that with the velocity-based methods the adaptive mesh is the result of time integration of the mesh velocity elds. Thus, the mesh history for previous time steps will inuence the mesh behavior at the current instant and the mesh can become highly skew. For the location-based methods , the adaptive mesh at a given time is determined by the monitor function at that time.

We also presented an original method called the LMPDE method. This method is quite interesting but, as mentioned before, it is still not clear how it may be extented to the use of unstructured meshes.

An optimization method was also explained. The main idea is to generate at each time step a mesh as close as possible to the Lagrangian mesh but with higher quality. This method is exible, generates high quality meshes and is also able to handle all types of meshes. One main drawback is that it is time-consuming.

Finally, we presented two interpolation methods that treat the problem of mesh deformation as a problem of interpolating the deformation between the moving geometry and the boundaries of the computational domain. The RBF method seems to be an interesting method in terms of robustness.

The problem is that it is a quite costly method in terms of CPU time especially, when compared to the IDW method. This latter seems to be a very ecient method. Actually, it is an explicit mesh deformation method that does not require solving a system of equations. This explicit formulation reduces the computational cost and makes the parallelization easier. This method is also robust especially when coupled to some mesh smoothing methods.

As a perspective, since the IDW seems to be an ecient method, it will be implemented and tested and probably coupled to some smoothing or local remeshing operations. In the case of using some local remeshing operations, the use of a conservative interpolation step is also required [Phd Thesis Chahrazade BAHBAH].

where .F

  is the Euclidean norm of the matrix, |.|the determinant of the matrix, J the Jacobian of the rezoned grid and J ref the Jacobian of the rezoned grid. The term |J| |Jref| in the denominator is the barrier function used to penalize degenerated solutions.

  α the coecients α j , β the coecients of the linear polynomial p, M b,b an n b × n b matrix containing the evaluation of the basis fucntion φ b