Thi Thanh

Mai Ta

Pascal Frey

M Florian De Vuyst

M Pascal Frey

M Yvon

M Bijan Mohammadi

M Franck Pigeonneau

Saint-Gobain Recherche

PRIVAT Invité M Yannick

M Pierre Saramito

Contents

Introduction

We present here two problems related to numerical simulation of viscous incompressible fluids that will be described and analyzed in the two main parts of this thesis:

-Part 1: Numerical simulation of the two-phase flow problem.

-Part 2: An example of shape optimization in fluid mechanics. This general introduction is aim at enlightening the structure of this manuscript, it contains neither technical details, nor references on the aforementioned topics (which are left for the detailed introduction of each part). We sketch out here the main problems, hence leading to define the content of each chapter.

The first concern of this thesis is the problem of two fluids flow or two-phase flow, i.e, we are interested in the simulation of the evolution of an interface (or a free surface) between two immiscible viscous fluids or two phases of a fluid. The main challenge in solving these problems is to define and to handle an accurate representation of the interface that separates the different fluids, even in the event of geometry or topology changes. In strongly deforming flows, this difficulty is related to maintaining a high-quality interface while conserving the total mass and allowing some extreme changes like folding, merging and breaking (see Figure 1). Another difficulty is the accurate computation or estimation of some algebraic quantities like the normal vector and the local curvature along the interface. Furthermore, a surface tension force must be considered in the model and accurately evaluated.

It seems obvious that any proposed approach for solving this problem involves two major problems, related respectively to: i. The resolution of the geometry and the evolution of the interface.

ii. The resolution of the equation governing the "movement" of the flows in a non-homogeneous medium. Among the numerous methods that have been developed for the simulation of free surface flows, the level set method has been extensively used in the past few years due to its simplicity and efficiency. More specifically, we propose here a general scheme for solving two fluids flow or twophase flows which takes advantage of the flexibility of the level set method for capturing evolution of the interfaces, including topological changes. Unlike similar approaches that solve the flow problem and the transport equation related to the evolution of the interface on Cartesian grids, our approach relies on an adaptive unstructured mesh to carry out these computations and enjoys an exact and accurate description of the interface. The explicit representation of the manifold separating the two fluids will be extracted to compute approximately the surface tension as well as some algebraic quantities like the normal vector and the curvature at the interface.

Figure 2: Illustration of a computational domain Ω(t) composed of two phases Ω 1 (t) and Ω 2 (t) having different material properties and separated by an interface Γ(t) which may evolve in time.

More precise description, the moving interface Γ(t) between two immiscible fluids occupying two sub-domain Ω i (t)(i = 1, 2) (see Figure 2) is represented as zero-level set of a level set function φ(x, t): Γ(t) = {x ∈ Ω(t) : φ(x, t) = 0}. Hence, the evolution of the interface is governed by the advection of the level set function φ: ∂φ ∂t + u.∇φ = 0

where the advection vector u is the velocity of the flows governed by the incompressible Navier-Stokes equations solving:

       ρ i ∂u i ∂t + (u i • ∇)u i -µ i ∆u i + ∇p i = ρ i f i in Ω i (i = 1, 2
)
divu i = 0 in Ω i (i = 1, 2)
with the interfacial condition on Γ :

u 1 -u 2 = 0 (σ 1 -σ 2 ).n 1 = -γκn 1
The surface tension is taken into account as the term of a localized force -γκn 1 at the interface.

In a nutshell, the resolution of a two-fluid problem is summarized by the steps involves the following ingredients:

-a numerical method for solving incompressible Navier-Stokes equations (Chapter 1); -a geometrical treatment to evaluate the surface tension (Chapter 2); -a numerical method for solving the level set advection (Chapter 2);

-the techniques of mesh adaptation (Chapter 3).

Let us finish the introduction of our first part by specifying the role of mesh adaptation in this scheme. It is obvious that no numerical method is completely exact in solving the PDE problem at hand, hence, we need a discretized computational domain. However, the accuracy of numerical solutions or the mass loss/gain can generally be improved with mesh refinement. The question that arises is related to where and how to refine the mesh. At each time, our mesh adaptation produces the adapted mesh based on the geometric properties of the interface and the physical properties of the fluid, simply speaking, only one adapted mesh at each time step to assume both the resolution of Navier-Stokes and the advection equations. It answers to the need for an accurate representation of the interface and an accurate approximation of the velocity of fluids with a minimal number of elements, then decreasing the amount of computational time.

CONTENTS

The second part of this thesis is related to shape optimization in fluid mechanics. Broadly speaking, this problem can be formulated as the minimization of an objective function J(Ω) of the domain variable Ω, via the solution u Ω of a mechanical problems (Stokes system in our context). By the classical way, the study of the derivative of J with respect to the domain makes it possible to compute a descent direction V for J from a given shape Ω, then Ω is updated to obtained the new shape Ω(t). However, there are a lot of difficulties in implementing this idea:

-The practical computation of the descent direction V usually involves the resolution of one, or several PDE systems posed on Ω (Stokes system and adjoint system). At this point, we need to construct the adjoint system which is also solved on Ω and leads to the resolution of V .

-The evolution of the shape Ω along the velocity field V is fairly straightforward in the theoretical framework, but unfortunately much harder in the numerical practice. In particular, it inherently depends on how Ω is parametrized. For instance, if Ω is described by a mesh, the evolution of Ω is naturally translated to the moving of the associated vertices in the direction of V . This sometimes produces an ill-shaped (or even invalid) mesh for the new shape (see Figure 3). In general, mesh evolution is a difficult issue, especially in three space dimensions.

Figure 3: Ilustration for the cases of ill-shaped mesh (top) and invalid mesh (bottom): (left) A velocity field V, defined at the vertices of a mesh; (right) deformed mesh.

So, in order to satisfy all the requirements of the computation of a descent direction for J and of the description of the domain evolution, several authors proposed to combine the aforementioned techniques of shape sensitivity analysis with the level set method. In these approaches, the general idea is enclosing all the possible shapes in a fixed, large computational domain D equipped with a fixed mesh D and to describe any shape Ω from an implicit point of view, via a scalar level set function φ : D → R which fulfills the following properties (see Figure 4): ∀x ∈ D,

             φ(x) < 0 if x ∈ Ω φ(x) = 0 if x ∈ ∂Ω φ(x) > 0 if x ∈ D \ Ω
The evolution of Ω(t) along the velocity field V is reformulated in terms of an associated level set function φ(., t) as the following level set advection equation:

∂φ ∂t + V.∇φ = 0
which can be solved on D, using its mesh D.

The computation of V is however not so easy in this context: we indeed evoked the fact that it requires solving PDE systems posed on Ω -a mesh of which is not available. These systems must then be approximated as PDE systems posed on the whole domain D (in the context of linearized elasticity, this is generally achieved by using the Ersatz material approach), and solved on the fixed mesh D. This operation may turn out to be difficult in the study of mechanical models which require a high accuracy in the description of the boundaries of shapes.

CONTENTS

Our approach arises from the observation that a slight modification in this methodology allows us to retain its great versatility when it comes to tracking the evolution of shapes, while benefiting from an exact description of any considered shape Ω ⊂ D.

Indeed, since a mesh of a shape Ω is needed at each time, one could imagine to modify D in such a way that an explicit discretization of Ω appears in it (see Figure 4). Hence, the computation of the descent direction V from Ω would become straightforward, and would not involve any approximation of the considered mechanical problem. Carrying out this idea inherently requires to be able to perform local mesh operations on D, hence to work with fully unstructured meshes; in our case, we shall use simplicial meshes (i.e triangulation in two dimensions, tetrahedralization in three dimensions).

Moreover, it implies the following ingredients:

-a numerical method for solving the mechanical problem (Stokes system and adjoint system). This is the reduced system (not containing the advection operation) of the Navier-Stokes equations that have been solved in Chapter 1.

-a numerical method for solving the level set advection equation on a simplicial computational mesh of D. This is one of the purposes of the work in Chapter 2 (Section 2.3) -a meshing technique for discretizing explicitly a shape known via an associated level set function, on a mesh of D. This is described in Chapter 3.

Let us now outline the description of the different chapters of this thesis.

The plan of the thesis

Chapter 1: Numerical resolution of Navier-Stokes equations

Chapter 1 focuses on the numerical solving of the incompressible Navier-Stokes equations for one viscous fluid. We rely on the Lagrange-Galerkin formulation for treating the convection operator ∂ ∂t + u.∇. By this way, at each time step the Navier-Stokes problem becomes an unsteady Stokes problem that can be solved using the finite element method. The resolution of the linear system from the finite element discretization is based on the Uzawa method or the Penalty method using Conjugate Gradient algorithm and a preconditionned Conjugate Gradient algorithm. A classical numerical test Lid-driven cavity is investigated in order to valid and assess our code with considered options in two and three dimensions spaces.

Chapter 2: The scheme for two-phase flow simulation

The second chapter begins with the introduction of the model of two-fluid flows with the interface.

Our model is based on the incompressible Navier-Stokes equation for the velocity field of the fluids, and the level set advection equation for the evolution of the interface. At each time step, the level function is transported by the velocity field of the fluid to define a new position of the interface. The resolution of this equation is carried out using the method of characteristics. Here, the characteristics are approximated by piecewise affine sequences using a Euler or 4 th -order Runge-Kutta schemes. Due to the explicit discretization of the interface in the computational mesh, the surface tension can be estimated geometrically. As known, a drawback of the level set method is the difficulty in preserving the regularity of the interface and the conservation of the mass. To overcome these problems, a redistancing procedure and a technique of mass correction are supplemented to our scheme. A numerical scheme based on mesh adaptation for two-phase flow simulation is proposed at the end of this chapter.

Chapter 3: Mesh adaptation

In this chapter we recall the major concepts and aspects of two mesh adaptation processes that have been used successfully in the general schemes of two-fluid flows as well as for shape optimization.

-Anisotropic mesh adaptation: the adapted mesh is generated based on the information stored in a metric tensor. This metric tensor is constructed from the requirements of the accuracy of geometrical approximation as well as from a numerical approximation of the concerned problems.

-Discrete three dimensional domain remeshing: The adapted mesh is generated based on a local size map. The local size map of the adapted mesh is computed to satisfy a given geometric approximation as well as to comply with numerical approximation constraints.

In these mesh adaptation processes, the final adapted mesh should be "well-shaped" mesh with respect to some quality criterions.

Chapter 4: Numerical examples

This chapter presents several examples of numerical resolutions of incompressible Navier-Stokes equations as well as of two-fluid flows involving an interface. Each case test will be considered in two and three dimensions. The first example is the Lid-driven cavity. This test has been presented in Chapter 1 to valid all implements in our code. In this chapter, examples are exhaustively investigated CONTENTS and compared with the results gathered in the references. Not only taking care of the number of the vortices, the streamlines, we compare also the position of the primary vortex and the velocity profile as well as the pressure of the fluid.

The efficiency and reliability of our numerical scheme for two-fluid flows are examined throughout several examples for which the surface tension is taken into account, namely, the rising bubble, the Rayleigh-Taylor instability, the coalescence of two rising bubbles. The results obtained are consistent with those found in the literature and show that the proposed scheme is able to manage complex interfacial movements, include topology changes.

Chapter 5: Shape optimization in fluid mechanics

The chapter begins with the shape sensitivity analysis based on Hadamard's boundary variation method in which the computation of shape derivative as well as the adjoint system are considered in the context of Stokes system. From the inferring of the descent direction, we play out a global algorithm for shape optimization in fluid mechanics. This algorithm is set in the context of level set method on unstructured meshes and mesh evolution to obtain a numerical scheme for shape optimization in fluid mechanics. A numerical example with the objective function of energy dissipation is presented to demonstrate the efficiency of the proposed scheme.

Introduction

Over the last decades, tremendous progress has been achieved on the development and the analysis of numerical methods for one-phase incompressible Stokes and Navier-Stokes flows, as emphasized by the vast literature on this topic. Open source and commercial software packages are now readily available and can be used as black-box solvers for a large class of industrial problems. However, challenging topics still require further investigation to reach the same level of maturity. For instance, the work that has been done on numerical methods for incompressible Navier-Stokes equations for one fluid is already a good starting point for dealing with two fluids or two-phase flow problems.

Actually, research on this topic has started in the last few years. However, there are several specific issues relevant to two-fluid problems that are not present in one-phase or one-fluid incompressible flow problems. In this context, the numerical treatment of the interface between two immiscible fluids is certainly a difficult challenge that raises many related and yet mostly unresolved problems:

the coupling between the fluid dynamics and the interface evolution, the conservation of mass, the treatment of singularities (geometrical and topological) and the approximation of surface tension forces. Not surprisingly given the level of difficulty of mathematical analysis, most research results in this topic have been published in the engineering literature, at the noticeable exception of the monograph [START_REF] Gross | Numerical methods for two-phase incompressible flows[END_REF]; see for instance [START_REF] Anderson | Diffuse-interface methods in fluid mechanics[END_REF] or [TBE + 01] for an overview of numerical methods for the simulation of multiphase flows.

Well-posedness results for the general weak formulation of the Navier-Stokes problem for twophase flows including the interface condition V Γ = u • n (see Section 2.1) have been analyzed only for special cases (e.g. unbounded domain Ω = R 3 and lim |x|→∞ u(x, t) = 0) and when the initial interface Γ(0) is a closed manifold [START_REF] Denisova | Problem of the motion of two viscous incompressible fluids separated by a closed free interface[END_REF]. The case of a bounded domain Ω for arbitrary time intervals [0, T ], T > 0 is treated in [START_REF] Tanaka | Global existence of two-phase nonhomogeneous viscous incompressible fluid flow[END_REF]; it provides a well-posedness result for the Navier-Stokes problem in a weak formulation. Most analysis apply to cases with sufficiently smooth data and do not apply when the regularity of the interface drops down, like when bubbles collide, for example. In such cases, curvature is no longer well-defined and weak alternatives have to be considered. These alternatives involve different representations of the interface which in turn induce relevant numerical techniques for the simulation of two-fluid flows.

Let us briefly recall the two most important methods for interface representation. Broadly speaking, these methods can be classified as Lagrangian ODE techniques and Eulerian PDE techniques.

Authors often introduce the terminology interface tracking versus interface capturing to characterize the treatment of the interface.

1. Interface tracking: When the interface regularity is sufficient, then normal, curvature and immiscibility condition like V Γ are well-defined. Given a velocity field u ∈ V , where V a suitable functional space, the trace u| Γ is well defined and the interface evolution can be described in Lagrangian coordinates. Each and every infinitely small particle in the domain Ω is transported (advected) by the flow field u(x, t) and we define the characteristic X(t) = X(x 0 , t 0 ; t) as the path of this particle with initial position x 0 . This trajectory is described by the following set of ordinary differential equations:

     dX(x 0 , t 0 ; t) dt = u(X(x 0 , t 0 ; t), t) t ≥ 0 X(x 0 , t 0 ; t 0 ) = x 0 .
(1)

For u(x, t) Lipschitz-continuous (with respect to x), this system has a unique solution and the regularity of X is related to the regularity of u. For T > 0 sufficiently small, there is a one to one mapping between (x, t) ∈ Ω × [0, T ] and a unique x 0 such that x = X(x 0 , t 0 ; t). Following the flow backwards in time starting from (x, t) yields the equation:

x = x 0 + t t 0 u(X(x 0 , t 0 ; t))dt , (2)
which represents a transformation from Eulerian to Lagrangian coordinates. The Navier-Stokes problem can be transformed accordingly into a non stationary problem with a stationary interface Γ(0) [START_REF] Tanaka | Global existence of two-phase nonhomogeneous viscous incompressible fluid flow[END_REF]. Likewise, the evolution of the interface Γ(t) can be described by using the Lagrangian coordinates and Γ(t) is simply characterized as the set of x ∈ Γ(t) satisfying

(2). This class of method is called interface tracking.

Practically, the Navier-Stokes problem is solved on a fixed grid or an unstructured mesh using an Eulerian approach and a Lagrangian approach is used to solve the evolution of the interface.

Marker points are equidistributed along the interface Γ(t 0 ) at time t 0 and then advected by the flow field u over a time period ∆t. Their final location mark the position of the interface at time t + ∆t. After several time steps, the equidistribution property is usually lost and marker points have to be redistributed along the new interface Γ(t + ∆t). Markers are usually connected to define a piecewise affine interpolation of the interface and can coincide with the set of vertices of a triangulation of Γ(t 0 ). In addition, this method requires to transfer information between the interface and the fixed grid once the interface has moved. Obviously, this approach is not very well-suited for dealing with topology changes or severe displacements (distortion) of the interface between two time steps. More details about this method and its implementation can be found in the references [UT92, ET98, ET99, HAC97].

To partially overcome these problems, hybrid approaches like the arbitrary Lagrangian-Eulerian (ALE) method have been proposed and revealed especially efficient for solving fluid-structure interaction [Bae01, Beh01, GT09, GT94]. They consist in solving the interface using a grid or a mesh and then the later is moved according to the flow velocity u. The mesh velocity in the interior of the domain generally differs from the flow velocity field, in order to avoid strong distortions.

In Volume of Fluid methods (VOF), the treatment of the interface is based on a weak formulation of the advection equation

∂χ 1 ∂t + u • ∇χ 1 = 0 (3)
where χ 1 (•, t) is the characteristic function for the subdomain Ω 1 (t), one of the two subdomains delimited by the interface Γ(t) hereby defined as

Γ(t) = ∂Ω 1 (t) = ∂supp(χ 1 (•, t)).
The function χ 1 is discontinuous across the interface and the transport equation requires a specific treatment.

Given an arbitrary small elementary volume of fluid W and integrating leads to the following equation:

∂χ 1 ∂t W χ 1 dx + ∂W χ 1 u • n ds = 0 , (4)
that can be interpreted as a weak formulation of (3) corresponding to volume conservation. The method typically involves two steps: at first, the reconstruction of the interface (approximation of the characteristic function) and then, the advection of the volume fraction function.

This VOF approach is widely used for the numerical simulation of two-fluid flows, mainly

φ(x) < 0, x ∈ Ω 1 (0), φ(x) > 0, x ∈ Ω 2 (0), φ(x) = 0, x ∈ Γ(0) . (5)
Assuming the velocity field u(x, t) is sufficiently smooth, then for t > 0 the level set values φ(x, t) are defined by considering the values constant along characteristics, namely by writing:

φ(X(x 0 , t 0 ; t)) = φ(x 0 ) , x 0 ∈ Ω, t ≥ 0 , (6)
and when differating this advection equation with respect to t it comes:

∂φ ∂t + u∇φ = 0 , for all 1 2 x ∈ Ω , t ≥ 0 . (7)
This equation (similar to (3)) is well-defined in its current fomulation given the velocity field u is Lipschitz-continuous with respect to x. Furthermore, the interface Γ(t) can be defined by the values of the auxiliary function φ at any time t:

Γ(t) = {x ∈ Ω ; φ(x) = 0} . (8)
There is no uniqueness of a solution for a general continuous velocity field u in this strong formulation. However, the notion of viscosity solutions of transport equations with a continuous velocity field eventually applies here, which yields to sub-or supersolutions [START_REF] Crandall | User's guide to viscosity solutions of second order partial differential equations[END_REF].

The level set (7) is not only used for the mathematical analysis of well-posedness of two-fluid flows but has also very attractive numerical features for representing and handling the interface Γ(t). We consider the initialization of φ with the signed distance function φ = d(x, Γ(0)) to the interface (equation (8)) at t = 0. The velocity field is the result of Navier-Stokes equation and the transport equation ( 7) is discretized using suitable numerical methods in space and time (see Sections 2.1 and 2.3 for more details of our implementation). As the level set function is continuous, its discretization is more accurate than that of the characteristic function considered in VOF methods. As the iterations in time increase, so does the discrepancy between the numerical solution φ h (x, t) and the signed distance function. A reinitialization of the level set function is then carried out when for example ∇φ h (x, t) 2 exceeds some given tolerance value.

Due to its simplicity and its ability to efficiently deal with topology changes, the level set method has been extensively used in engineering applications. Additionally, the extension from two to three dimensions can be achieved easily. Like the ALE method described above, variants of the level set method can be considered, in which the interface is explicitely discretized (hence the terminology interface capturing) by a mesh and this mesh is moved with the flow velocity.

As will be seen, intrinsic properties of the interface can be accurately computed using the level set function [START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF]. In the remainder of this manuscript, we will restrict our numerical investigation to only one of these, the level set representation.

In this part, we will describe in details our work and propose a general strategy for solving two-phase flows which takes advantage of the flexibility of the level set method for capturing and tracking evolution of the interfaces, including topological changes, and enjoys an exact and accurate description of the interface using a conforming unstructured mesh. The numerical resolution of incompressible Navier-Stokes equations is presented in Chapter 1, our approach is based on the Lagrange-Galerkin scheme was first introduced by [START_REF] Benqué | A finite element method for the Navier-Stokes equations[END_REF] and studied in [START_REF] Pironneau | On the transport-diffusion algorithm and its applications to the Navier-Stokes equations[END_REF]. Chapter 2 presents the combining an implicit method for dealing with the domain evolution and an explicit representation of the manifold separating the two fluids. This idea is obviously not new but unlike similar approaches that solve the flow problem and the transport equation on Cartesian grids [ALTP98, GTBD06, EFFM02, MG07, BCdV14], our approach relies on an adaptive unstructured mesh to carry out these computations. Chapter 3 is devoted to the techniques of mesh adaptation as well as the construction of adapted mesh process. Many numerical tests are investigated in Chapter 4 in two and three dimensions in order to assess the efficiency and reliability of our proposed numerical schemes. The purpose of this chapter is to present the numerical methods we have designed to resolve Navier-Stokes equations. Our approach is based on the combination of the method of characteristics with a finite element formulation. This method was first introduced by Benqué [BIKL80] and analyzed in [START_REF] Pironneau | On the transport-diffusion algorithm and its applications to the Navier-Stokes equations[END_REF], and is also known as the Lagrange-Galerkin method. The main idea behind this method is that the convection operator (the non linear term) of the Navier-Stokes equation can be turned into a total derivative (also called material derivative) by using a Lagrangian formulation.

The outline of this chapter is the following: After introducing the incompressible Navier-Stokes equations for a single viscous fluid in Section 1.1, we express the operator of velocity ∂ ∂t +u.∇ by material derivative d dt which allows us to discretize in time this problem in Section 1.2. In this approach, the non-linear part of Navier-Stokes equations is associated with a Cauchy problem, and this leads us to solve an unsteady Stokes problem at each time step. The following sections are concerned with the numerical resolution of a generalized Stokes problem by the finite element method: section 1.3 presents the variational formulation for different boundary conditions. Section 1.4 presents the discrete formulation corresponding to Lagrange-Galerkin finite elements of P1bubble/P1 and P2/P1. Section 1.5 describes the matrix formulation and the resolution of the resulting linear system discretized by finite elements spaces. In the last section 1.6 we detail a numerical test, the Lid-driven cavity in order to validate and to compare the numerical resolution of the incompressible Navier-Stokes equations in two and three dimensions with different parameters and settings related to: finite elements, the order of the scheme for approximating of the characteristic lines and the resolution of linear systems.

The Navier-Stokes equations for incompressible fluid

Let Ω be a bounded domain in R d (d ∈ {2, 3}) and let denote Σ, its boundary. We consider the Navier-Stokes equations governing unsteady incompressible flows in general:

       ρ ∂u ∂t + (u • ∇)u -µ∆u + ∇p = ρf in Ω × [0, T ] divu = 0 in Ω × [0, T ] (1.1)
where:

u = u(x, t) is velocity of a fluid particle occupying the position x at the time t and p = p(x, t) is its pressure;

ρ and µ are the constant density and dynamic viscosity, respectively;

f = f (x, t) is an internal force exerted on the fluid per unit mass;

-and

(u • ∇)u = Σ d i=1 u i ∂ i u, u i is the i th -component of velocity field u, ∂ i u is the partial derivative according to the variable x i ; divu = Σ d i=1 ∂ i u i .
The first equation of (1.1) expresses the conservation of motion of the fluid while the second one expresses the conservation of mass for an incompressible fluid. By posing ν = µ ρ ; p = p ρ we obtain the formula of the incompressible Navier-Stokes equations as follows:

       ∂u ∂t + (u • ∇)u -ν∆u + ∇p = f in Ω × [0, T ] divu = 0 in Ω × [0, T ] (1.2)
In these equations, ν represents the constant kinematic viscosity of the fluid. By giving L a characteristic length of the fluid, we can define the Reynolds number Re as:

Re = |u|.L ν (1.3)
This number represents the ratio between the inertial force and the viscous force. In case of large value of ν, the viscous force dominates and the convection term in the Navier-Stokes equations can be neglected. This leads to consider the unsteady Stokes equations:

     ∂u ∂t -ν∆u + ∇p = f in Ω × [0, T ] divu = 0 in Ω × [0, T ] (1.4)
In order to ensure the well-posedness of the problem, the equations (1.2) need to be supplemented with a set of adequate boundary conditions (see Section 1.3). They are also endowed with an initial condition: a divergence-free velocity field u 0 (x) is specified over the domain Ω at time t = 0, i.e.

u(x, 0) = u 0 (x), divu 0 = 0.

The method of characteristics

In this section we focus on the time discretization of the Navier-Stokes equations by the method of characteristics. This method, also known as the Lagrange-Galerkin method was introduced by Benqué [START_REF] Benqué | A finite element method for the Navier-Stokes equations[END_REF]. It provides a time discretization without the necessity of restrictive stability conditions as in a classical discretization by an explicit scheme.

The Lagrange-Galerkin formulation

We are now introducing the Lagrange-Galerkin formulation to recast the convective operator of the Navier-Stokes equations. Let X(x, s; t) be a characteristic curve associated with the velocity field u, solution of the ordinary differential equations:

     dX(x, s; t) dt = u(X(x, s; t), t) X(x, s; s) = x (1.5)
where the point X(x, s; t) denotes the position at time t of a fluid particle located at position x at the time s.

We write the derivative of the velocity field u along the characteristic lines as: It can be observed that the treatment of the nonlinear convection term is therefore reduced to a problem of searching the characteristic foot X(x, s; t), i.e the position of the particle at the previous time. This approach allows us to avoid theoretically the constraint of CFL (Courant-Friedrichs-Levy) on the time step (for example in [START_REF] Pironneau | Characteristic-Galerkin and Galerkin/leastsquares space-time formulations for the advectiondiffusion equation with timedenpendent domains[END_REF] the good choice for time step is proposed as ∆t ≈ 1.5h/|u|).

du(X(
Moreover, it has been shown in [START_REF] Pironneau | On the transport-diffusion algorithm and its applications to the Navier-Stokes equations[END_REF] that if the characteristic trajectory is calculated exactly and the second member is also integrated exactly, the resulting scheme is unconditionally stable.

The discretization in time

We assume that the interval [0, T ] is divided into a number of subintervals ∆t and we denote t n = n∆t. For the numerical solution, we only compute u(x, t n ) for all n, so from now on, we will denote u(x, t n ) by u n (x).

Using the following approximation of the total derivative along the characteristic curves:

du(X(x, t n ; t), t) dt u(x, t n ) -u(X n-1 (x), t n-1 ) ∆t
where x is the position of the particle at current time t n and X n-1 (x) is its foot at t n-1 on the characteristic curve (i.e X n-1 (x) stands for X(x, t n ; t n-1 )), we find the discretization formulation of the Navier-Stokes equations (1.2):

     u n (x) -u n-1 (X n-1 (x)) ∆t -ν∆u n (x) + ∇p n (x) = f n divu n (x) = 0
which are equivalent to:

     u n (x) ∆t -ν∆u n (x) + ∇p n (x) = f n + u n-1 (X n-1 (x)) ∆t divu n (x) = 0 (1.8)
where u n-1 (X n-1 (x)) represents the velocity at the point X n-1 (x) at the time t n-1 .

Hence, at each time step, the Navier-Stokes problem reduces to an unsteady Stokes problem plus a transport of the previous solution along the characteristic lines under a velocity field u.

At each time steps, the discretization in time of the Navier-Stokes equations is carried on in two steps :

i. Resolution of the problem (1.5) in the sub-interval

[t n-1 , t n ].
ii. Resolution of the generalized Stokes equations defined by:

     ku -ν∆u + ∇p = g divu = 0
(1.9)

where k is denoting 1 ∆t and g is denoting all terms in the right-hand side and is updated in time according to formulation (1.8).

Remark 1.

1. The advantage of the discretization of Navier-Stokes equations using this scheme lies in the fact that at each iteration, we only need to update the right-hand side, i.e the left-hand side remains constant during the resolution.

2. Concerning the problem of the discretisation of the characteristic lines, or more precisely the approximation of X n-1 (x), we forecast some difficulties: i. The first problem is how to approximate the velocity field u(t) in the sub-interval time [t n-1 , t n ].

In solving Navier-Stokes equations, u n is unknown, and thus X n-1 (x) is associated with u n-1 .

ii. The second difficulty is that the characteristic line may cross some elements of the domain or even go out of the computational domain, hence it is necessary to locate the element which contains the last point of the advection or to find the cases in which the latter hits a boundary edge. In numerical practice, the backward tracking the characteristic points are performed by Euler's scheme and Runge-Kutta 4. These scheme will be detailed in Section 2.3 where the method of characteristics is again applied to resolve the advection equation of the interface. See We know that in practice the exact calculation of the characteristic curves or the second member is impossible, it must be done using numerical integration formulas. Unfortunately, the numerical approaches prevent the unconditional stability of the Lagrange-Galerkin scheme, either of order 1 or order 2. For more exhaustive studies about the stability of Lagrange-Galerkin method, we refer to [START_REF] Morton | Stability of the Lagrange-Galerkin method with non-exact integration[END_REF] or [START_REF] Fourestey | Stabilité des méthodes de Lagrange-Galerkin du premier et du second ordre[END_REF] and references therein. We turn to the numerical practice where the stopping condition is performed when the residual of solution is small enough, one say that the solution reaches a steady-state (of course if it exists). More precisely, given a computational mesh, we construct a global scheme for solving Navier-Stokes equations: while ((n ≤ nmax) and (res > ε)) do Integrate X n-1 (x) (equation (1.5)) for each node x of mesh;

Compute u n-1 (X n-1 (x)) at each node x of mesh;

Update the right hand side g n := f n + u n-1 (X n-1 (x));

Solve Stokes equations (1.9) to obtain (u n , p n );

Compute residual at n th -step: res := u nu n-1 0 ; n:=n+1; end where ε is a given small value introduced to define the steady-state solution in each test case.

The variational formulation

As we discussed in the previous section, at each time step, we have to solve an unsteady Stokes problem (1.9). The numerical resolution of this problem by a finite element method evidently involves a variational formulation.

Let us firstly introduce the following Sobolev spaces which will be used hereafter:

H 1 0 (Ω) d = {v ∈ H 1 (Ω) d : v| Σ = 0} H 1 Σ D (Ω) d = {v ∈ H 1 (Ω) d : v| Σ D = 0} L 2 0 (Ω) = {q ∈ L 2 (Ω) : Ω q = 0} 1.3.

Homogeneous Dirichlet boundary conditions

Denoting V and M for the functional spaces of velocity and pressure, respectively. We consider the variational formulation for the general Stokes problem in case of homogeneous Dirichlet boundary condition, i.e u = 0 on Σ. In this case, V will be chosen to coincide with H 1 0 (Ω) d and we take M = L 2 0 (Ω) for the pressure. Given v ∈ V and using Green's formula we have:

-

Ω ν∆uvdx = Ω ν∇u : ∇vdx - ∂Ω ν∇un.vds Ω ∇p.vdx = - Ω pdivvdx + ∂Ω pn.vds
Using these equalities and taking the inner product in L 2 (Ω) d of the first equation of systems (1.9), we come to the following equation:

k Ω u.vdx + Ω ν∇u : ∇vdx - Ω pdivvdx = Ω g.vdx + ∂Ω (ν∇u -pI d )n.vds (1.10)
Since v = 0 on ∂Ω, we have:

k Ω u.vdx + Ω ν∇u : ∇vdx - Ω pdivvdx = Ω g.vdx
Choose the test function q ∈ M , taking the inner product in L 2 (Ω) of the second equation of system

(1.9), we obtain finally the variational formulation of the homogeneous problem: Given the function g ∈ H -1 (Ω) d , constant k and viscosity ν; find u ∈ V and p ∈ M solution of:

       k Ω u.vdx + Ω ν∇u : ∇vdx - Ω pdivvdx = Ω g.vdx ∀v ∈ V Ω qdivudx = 0 ∀q ∈ M CHAPTER 1. NUMERICAL RESOLUTION OF NAVIER-STOKES EQUATIONS
This problem can be written in the equivalent weak formulation: Find (u, p) ∈ (V × M ) such that:

     ∀v ∈ V a(u, v) + b(v, p) = l(v) ∀q ∈ M b(u, q) = 0 (1.11)
where l(•) is the continuous linear form defined on V:

l(v) = Ω g.vdx
and a(u, v), b(v, p) are the continuous bilinear forms defined on V × V and V × M respectively as:

       a(u, v) = k Ω u.vdx + Ω ν∇u : ∇vdx b(v, p) = - Ω pdivvdx (1.12)
The existence and the uniqueness of a solution for the weak formulation of the generalized Stokes problem has been proven (see [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations[END_REF] or [START_REF] Ern | Theory and Practice of Finite Elements[END_REF]). This proof involves: i. the ellipticity of the form a(., .) results from the Poincare-Friedrichs inequality; ii. the compatibility of the spaces of velocity and pressure results the satisfying Babuska-Brezzi condition, also called inf-sup condition on the form b(., .), i.e. existing a positive constant C such that:

inf q∈M sup v∈V b(v, q) v 1 q 0 ≥ C > 0 (1.13)
where

v 1 = Σ d i=1 v i 1 2 1/2
and . 1 , . 0 are standard notations of norms in the Sobolev spaces

H 1 (Ω), L 2 (Ω) respectively.

Case of general boundary conditions

Let us now consider the variational formulation in the case of more complicated boundary conditions. We limit ourselves to consider somme classical boundary conditions (Dirichlet, Newmann and Slip) which are used in most numerical tests. Suppose that the boundary Σ of Ω is partitioned into a finite number of components corresponding to different types of boundary condition:

Σ = Σ D ∪ Σ N ∪ Σ S . i. Dirichlet conditions: u = u D on Σ D .
ii. Newmann conditions: σn = ϕ N on Σ N .

iii. Slip conditions: u.n = 0, αu.τ + τ.σn = 0 on Σ S .

where σ is defined as: σ = ν∇u -pI d ; I d is the identity matrix and τ, n are the unit tangent vector and unit outward normal vector to Σ, respectively; α ≥ 0 is a constant of friction coefficient. In the case of Dirichlet boundary condition, V is usually chosen in such a way that the test function vanishes on Σ D :

V = H 1 Σ D d = {v ∈ H 1 (Ω) d : v| Σ D = 0}
Hence,

Σ D σn.vds = 0
It can be observed that the condition Dirichlet affects only on the selection of space functions for the velocity, it doesn't affect the variational formulation, so it is called essential boundary condition.

When the non-homogeneous Dirichlet is imposed, i.e u D = 0, u D is supposed to be a sufficiently smooth on Σ D such that there exists a lifting function ūD with divergence-free property. Using the expression: u = ūD + ū where ū ∈ V . The weak formulation becomes:

Find (ū, p) ∈ (V × M ) such that:      ∀v ∈ V a(ū, v) + b(v, p) = l(v) -a( ūD , v) ∀q ∈ M b(ū, q) = 0 (1.14)
In practice the Dirichlet condition is simply imposed by multiplying a very large value on the righthand side and the diagonal matrix of matrix formulation (will be presented in the next section).

Notice that V will coincide with (H 1 0 (Ω)) d if Σ D ≡ Σ that means only boundary condition of Dirichlet type is imposed. In this case, the pressure only appears through its gradient, so, if (u * , p * ) is a solution of Navier-Stokes equation then for any constant c, (u * , p * + c) will be also a solution.

We chose

M = L 2 0 (Ω) = {p ∈ L 2 (Ω) : Ω p = 0}
to avoid this indeterminacy.

The term corresponding to a Newmann condition can be written as:

Σ N σn.vds = Σ N ϕ N .vds
So, Newmann condition needs only adding this boundary term in the right-hand side of the variational formulation. This is one kind of natural boundary condition. In the case ϕ N = 0, the subset Σ N is called a free-out flow part, this integral term disappears from variational formulation. Notice that when Σ N is not empty, this means that we have condition for pressure on Σ N , the problem of pressure indeterminacy previously in the case only Dirichlet condition, no longer exists, so we can take M = L 2 (Ω).

In the case of slip boundary condition, the test function v is chosen in V and satisfies v.n| Σ S = 0.

The term of this condition can be written as:

The spatial discretization

By using a Galerkin finite element approximation, the discrete problem corresponding to the formulation (1.11) writes as follows:

Find (u h , p h ) ∈ V h × M h s.t      ∀v h ∈ V h a h (u h , v h ) + b h (v h , p h ) = l h (v h ) ∀q h ∈ M h b(u h , q h ) = 0 (1.15) with V h ⊂ V and M h ⊂ M represent two families of finite dimensional subspaces and a h (u h , v h ), b h (v h , p h ), l h (v h ) are defined on V h × V h , V h × M h and V h respectively as follows:                      a h (u h , v h ) = K∈T h k K u h .v h dx + K∈T h K ν∇u h : ∇v h dx b h (v h , p h ) = K∈T h K -p h divv h l h (v h ) = K∈T h K g h .v h dx (1.16)
It is well-known that the approximate problem also requires the discrete compatibility condition, meaning that the discrete spaces of velocity needs to be "rich" enough comparing with the one of pressure, i.e there exists a positive constant C h s.t:

inf q h ∈M h sup v h ∈V h b h (v h , p h ) v h 1 q h 0 ≥ C h > 0 (1.17)
which is also known as the discrete inf-sup condition.

There are many possible choices for the discrete spaces which satisfy this condition. A description of some compatible couple of spaces for finite element spaces can be found in [START_REF] Pironneau | The finite element methods for fluids[END_REF], [BBD + 06] or [START_REF] Quarteroni | Numerical Models for Differential Problems[END_REF]. We have decided to retain mini elements (P 1 bubble/P 1 ) and Taylor-Hood elements (P2/P1) have been implemented. We briefly recall here these elements.

The mini elements: All polynomial approximation of same order for both velocity and pressure are ill-posed. Suppose that T h is composed of elements K = (P 1 , ..., P d+1 ), from the approximation of piecewise linear P1/P1 finite elements, in order to satisfy the inf-sup condition (1.17), the spaces of velocity can be enriched by adding one degree of freedom b K on each elements where b K is socalled bubble function defined as: b K = (d + 1) d+1 d+1 i=1 ϕ i with ϕ i is P 1 Lagrange basic function corresponds to vertex P i of elements K. We have following discretizations of the spaces of velocity and pressure in this case (see figure 1.2 (a)):

V h = {v h ∈ (C 0 ( Ω)) d : v h | K ∈ (P 1 (K) ⊕ b K ) d ∀K ∈ T h } ∩ V M h = {q h ∈ C 0 ( Ω) : q h | K ∈ P 1 (K)∀K ∈ T h } ∩ M
The Taylor-Hood elements: Generally, this choice is to use piecewise polynomials of degree k ≥ 2 for the velocity space V h and of k -1 for pressure space M h . In particular, we use piecewise quadratic P 2 finite elements for each velocity component, and piecewise linear P 1 finite elements for pressure. This is the lowest degree representative of the family of so-called Taylor-Hood elements P k /P k-1 , k ≥ 2 (continuous velocities and continuous pressure), that are inf-sup stable (see figure 1.2 (b)):

V h = {v h ∈ (C 0 ( Ω)) d : v h | K ∈ (P 2 (K)) d ∀K ∈ T h } ∩ V M h = {q h ∈ C 0 ( Ω) : q h | K ∈ P 1 (K)∀K ∈ T h } ∩ M
See figure 1.2 for ilustrations of these elements in two dimensions with the denoting of symbol • for degrees of freedom of the velocity and symbol • for ones of the pressure.

The discrete linear systems

Construction of the matrix formulation

Let {ϕ i } i=1,...,nv and {φ j } j=1,...,np are basis functions of V h and M h , respectively with nv = dimV h and np = dimM h . Given u h ∈ V h , p h ∈ M h , we have the following decompositions:

u h = Σ nv i=1 u i ϕ i p h = Σ np j=1 p j φ j (1.18)
We can introduce the following matrix coefficients:

A i,j = a h (ϕ i , ϕ j ) B i,j = b h (ϕ i , φ j ) G i = l h (ϕ i )
where the matrices A, B correspond to the bilinear forms a h , b h , respectively and the vector G corresponds to the linear form l h according to expression (1.16). Denoting the vectors U = (u i ) i=1,...,nv and P = (p j ) j=1,...,np , the problem (1.15) is equivalent to solving the square linear system:

  A B t B 0     U P   =   G 0   (1.19)
The system (1.19) is sparse, symmetric but not definite and its size is dimV h + dimM h .

Matrix assembly

In practice, the assembling matrices A, B and G in section 1.5.1 are accomplished by following the mesh, or in other words, from programming point of view, by using a loop on the indices of mesh elements. All the contributions on each element K will be computed and stored in the local matrices We present here the contributions of these local matrices in two dimensions. In three dimensions, the procedure is almost identical.

denoted by [Ae] K , [Be] K , [Ge] K . Matrices A,

The local matrix of matrix A

Recall the bilinear form of a(u h , v h ):

a h (u h , v h ) = K∈T h k K u h .v h dx + K∈T h K ν∇u h : ∇v h dx
We have the bilinear form a h (u h , v h ) on each element K of mesh as follows:

a h (u h , v h )| K = k K u h .v h dx + K ν∇u h : ∇v h dx (1.20)
Using the matrical notations in the transformation of bilinear form a h (u h , v h ):

k K u h .v h dx = k K   v 1 v 2   t   u 1 u 2   dx = k K   v 1 | K v 2 | K   t   ϕ| K 0 0 ϕ| K   t   ϕ| K 0 0 ϕ| K     u 1 | K u 2 | K   dx (1.21) K ν∇u h : ∇v h dx = K         ∂ 1 v 1 ∂ 2 v 1 ∂ 1 v 2 ∂ 2 v 2         t         ν 0 0 0 0 ν 0 0 0 0 ν 0 0 0 0 ν                 ∂ 1 u 1 ∂ 2 u 1 ∂ 1 u 2 ∂ 2 u 2         dx = K   v 1 | K v 2 | K   t   Dϕ| K 0 0 Dϕ| K   t         ν 0 0 0 0 ν 0 0 0 0 ν 0 0 0 0 ν           Dϕ| K 0 0 Dϕ| K     u 1 | K u 2 | K   dx (1.22)
where:

[u i | K ] t = (u ϕ 1 i , . . . , u ϕnev i )(i = 1, 2); [ϕ| K ] = (ϕ 1 , . . . , ϕ nev ); [Dϕ| K ] =   ∂ 1 ϕ 1 . . . ∂ 1 ϕ nev ∂ 2 ϕ 1 . . . ∂ 2 ϕ nev  
with ϕ 1 , . . . , ϕ nev are Lagrange basis functions of V h on element K, nev is the number of degree of freedom of velocity on each element and u ϕ 1 i , . . . , u ϕnev i are nodal values of velocity components correspond to ϕ 1 , . . . , ϕ nev ; respectively.

Hence, the contributions associated with element K of global matrix A will be stored in the local matrix [Ae] K defined by:

[Ae] K = k K   ϕ| K 0 0 ϕ| K   t   ϕ| K 0 0 ϕ| K   dx + K   Dϕ| K 0 0 Dϕ| K   t         ν 0 0 0 0 ν 0 0 0 0 ν 0 0 0 0 ν           Dϕ| K 0 0 Dϕ| K   dx (1.23)
This integral is computed based on the quadrature formulation, we have:

[Ae] K = k np p=1   ϕ| K 0 0 ϕ| K   t   ϕ| K 0 0 ϕ| K   xp + np p=1   Dϕ| K 0 0 Dϕ| K   t         ν 0 0 0 0 ν 0 0 0 0 ν 0 0 0 0 ν           Dϕ| K 0 0 Dϕ| K   xp (1.24)
where x p are Gaussian points and n p is the number of Gaussian point in element K. It can be seen that the size of local matrix [Ae] K is 2nev × 2nev. (nev = 4 in case of P1-bubble and nev = 6 in case of P2)

The local matrix of matrix B

Similarly, recall the bilinear form of b h (v h , p h ):

b h (v h , p h ) = K∈T h K -p h divv h dx (1.25)
and using the matrical notations we have:

b h (u h , q h )| K = - K q h divu h dx = K   -q -q   t   ∂ 1 u 1 ∂ 2 u 2   dx = K   q| K q| K   t   -φ| K 0 0 -φ| K   t   ∂ 1 ϕ| K 0 0 ∂ 2 ϕ| K     u 1 | K u 2 | K   dx (1.26)
where:

[q| K ] t = (q φ 1 , . . . , q φnep ), [φ| K ] = (φ 1 , . . . , φ nep ), [∂ 1 ϕ| K ] = [∂ 1 ϕ 1 , . . . , ∂ 1 ϕ nev ], [∂ 2 ϕ| K ] = [∂ 2 ϕ 1 , . . . , ∂ 2 ϕ nev ]
with φ 1 , . . . , φ nep are Lagrange basis functions of M h on element K, nep is the number of degree of freedom of pressure on each mesh element and q φ 1 , . . . , q φnep are nodal values of test function q h correspond to φ 1 , . . . , φ nep , respectively.

The local matrix [Be] K is then determined as follows:

[Be] K = K   -φ| K 0 0 -φ| K   t   ∂ 1 ϕ| K 0 0 ∂ 2 ϕ| K   dx (1.27) We have [Be] K is size of 2nep × 2nev.
This integral is also computed by a quadrature formulation as in the computation of [Ae] K .

The local matrix of matrix G

The last term is linear form l h :

l h (v h ) = K∈T h K g h .v h dx Hence, l h (v h )| K = K g h .v h dx = K   v 1 v 2   t   g 1 g 2   dx = K   v 1 | K v 2 | K   t   ϕ| K 0 0 ϕ| K   t   ϕ| K 0 0 ϕ| K     g 1 | K g 2 | K   dx (1.28)
where:

[q i | K ] t = (q ϕ 1 i , . . . , q ϕnev i )(i = 1, 2); [ϕ| K ] = (ϕ 1 , . . . , ϕ nev ).
We have the local vector on the right-hand side:

[Ge] K = K   ϕ| K 0 0 ϕ| K   t   ϕ| K 0 0 ϕ| K     g 1 | K g 2 | K   dx (1.29)
is size of 2nev and computed by using a quadrature formulation.

The overall assembly algorithm of matrix A is given as following:

Algorithm 2: Assembly of matrix A, case 2D Initialization A = 0.

Loop on k = 1, 2, . . . , ne. Let K be the k th element.

Loop on local indices

il = 1, 2, . . . , 2nev of [Ae] K Find global index ig of matrix A corresponding to local index il of [Ae] K : ig := 2 * [pt → v[il%nev]] + il/nev; (where pt → v[il%nev] is global index corresponding to (il%nev) th point of element K).
Loop on local indices jl = il, . . . , 2nev.

Find global index jg corresponding to local index jl of [Ae] K :

jg := 2 * [pt → v[jl%nev]] + jl/nev; (where pt → v[jl%nev] is global index corresponding to (jl%nev) th point of element K).
If (jg <= ig) then A jg,ig := +Ae jl,il ;

else A ig,jg := +Ae il,jl ;

The assembly matrices B, G are implemented similarly.

Solving the linear systems

Numerous methods have been proposed to solve the linear system (1.19). A classical alternative is Uzawa method [AHU58] which relies on coupling (1.19) in to two subsystems, one for the unknown velocity U and other for pressure P , which leads to solving successively:

BA -1 B t P = BA -1 G, AU = G -B t P. (1.30)
Attentively, the so-called Penalty method can be used to solving U and P in an equivalent linear system as follows:

  A B t B I d     U P   =   G 0   (1.31)
where I d is identity matrix, its size is dimM h and values about 10 -6 to 10 -4 .

So, at each time step a generalized Stokes problem leads to solve one or severals linear systems.

Our code based upon solving linear system by a Conjugate Gradient (CG) algorithm and a SSORpreconditionned CG algorithms. These techniques have been applied popularly in solving sparse linear system, see [She94, Saa03, YJ07] for details.

Numerical examples

We have developed one solver Navier-Stokes in two dimension (2D) and three dimensions (3D)

with the following options:

1. The discretisation in time is based on the method of characteristics where the characteristic curves is approximated by Euler's scheme or scheme Runge-Kutta 4.

2. The disretisation in space is programmed the by P 1 bubble/P 1 or P 2 /P1.

3. The resolution of linear system is implemented by method of Uzawa or Penalty method.

In order to valid the solver with the previous options, we perform the numerical test of Lid-driven cavity 2D with Re = 1000. The detailed results of this test will be investigated in Chapter 4

(Section 4.1) in comparison with the other well-known references. This section aims to compare the convergence of all the scheme in our solver. The first one (the scheme of element P1bP1, Euler's scheme and Uzawa's method) is chosen as the reference. In case of this Reynolds number, all the experimental and numerical observations showed that the solution reach to a steady-state. With the same stopping condition when the residual between the solutions reaches to 10 -6 , it has been observed that the solutions are exactly the same for all schemes. The required CPU time of these resolutions to obtain the steady solution are detailed on Table 1.1 in comparison with the reference one.

From the details of the In the theorical analysis of method of characteristics for the Navier-Stokes equations, the error estimation has been determined of the form O(h m + ∆t m + h m+1 /∆t) in [START_REF] Pironneau | On the transport-diffusion algorithm and its applications to the Navier-Stokes equations[END_REF] (a scheme of order (h + ∆t + h 2 /∆t) for Crouzeix-Raviant element and a scheme of order (h 2 + ∆t 2 + h 3 /∆t) for Taylor-Hood element). It has seen that the error estimate is controlled by either characteristic mesh element h, ∆t or the relation of h and ∆t. An optimal error estimates for the Lagrange-Galerkin mixed finite element approximation of the Navier-Stokes equations was given by Suli [START_REF] Suli | Convergence and nonlinear stability of the Lagrange-Galerkin method for the Navier-Stokes equations[END_REF], see also [START_REF] Zhi | Second order modified method of characteristics mixed defect-correction finite element method for time dependent Navier-Stokes problems[END_REF] and references therein. More exhaustive studies about the error estimations for each finite elements lie far beyond the scope of this thesis. However, by numerical experience we propose here the choice of time step to obtain the best solution with the present scheme. In fact, the choice of time step is very important: with the different time steps we obtain the different behaviour of the steady solutions. For this illustration, we perform the test on the mesh carre7 (see figure 1 The present solver can handle with the cases of high Reynold number. In fact, we have investigated the test of cavity with the number of Reynolds up to 10000 and the resultats are in good agreement with other references, (see section 4.1).

Conclusion

We have presented in detail all the features in the resolution of viscous incompressible Navier-Stokes equations in two and three dimensions. The Navier-Stokes solver is capable of dealing with high Reynold number problems. Moreover, it will be still valuable for solving bifluid problems. In fact, the resolution of Navier-Stokes equation for two fluids is one part in our simulation of bifluid flows with interface will be detailed in the next chapter.

Regarding the perspectives for this chapter, we would like to:

-Improve the Lagrange-Galerkin scheme of second order to comparer with the one of first order.

-Consider the P1P1 elements in combination with a stable scheme that is hoped that dramatically reduce the computational time of solving linear systems, either solving by Uzawa or Penalty method.

-Develope the present solver with a technique of mesh adaptation to treat the more complex test of Navier-Stokes problem.

Chapter 2

Two-phase flow simulations In the previous chapter, all the resolution of Navier-Stokes for a single fluid has been described.

This chapter continues with the numerical resolution of bifluid flow. In section 2.1 we introduce the model of bifluid flow with the continuous condition on the interface between two fluids. This condition induces effect of surface tension as a localized force added to the right-hand side of the momentum equation. This model for a bifluid incompressible flow problem is often used in the literature, see [START_REF] Girault | One time-step finite element discretization of the equation of motion of two-fluid flows[END_REF][START_REF] Grob | A finite element based level set method for twophase incompressible flows[END_REF] and many references therein. Hence, the motion of two fluids can be achieved by the solver presented in Chapter 1 to obtain the velocity field of the flows.

The appearance of localized surface tension force and its geometrical approximation are described in section 2.2. In our approach, the interface is implicitly captured as zero-level set of a levelset function. The advection of level-set functions by the velocity of the flows (may be extended and regularized) obeys a Hamilton-Jacobi equation will be handled numerically in section 2.3. As known before, in the problem involving incompressible flows, the level-set function does not assure the regularity during the time of evolution leading the diffusion on the interface and difficulties in preserving the mass conservation. In order to overcome these problems, some numerical techniques are proposed in section 2.4 to ameliorate the numerical resolution. A numerical scheme for two-phase flow is proposed in the last section 2.5.

Introduction the model of two-phase flow

The equation of motion of bifluid flow

We consider the problem of two immiscible fluids illustrated by figure 2.1 (in 2D). At each time and the interface between the two fluids is Γ(t) = Ω 1 (t) ∩ Ω 2 (t), and we suppose that:

Ω 1 (t) ∪ Ω 2 (t) = Ω, Ω 1 (t) ∩ Ω 2 (t) = ∅
Supposing ρ(x, t) = ρ(x) and µ(x, t) = µ(x), meaning that the physical characteristics of each fluid remain constant during each time step of the resolution. For the sake of simplicity, we omit the dependence on t, write simply Ω i and Γ instead of Ω i (t) and Γ(t), at each time t ∈ [0, T ], the motion of bifluid flow governed by incompressible Navier-Stokes equations can be written as follows:

       ρ i ∂u i ∂t + (u i • ∇)u i -µ i ∆u i + ∇p i = ρ i f i in Ω i (i = 1, 2) divu i = 0 in Ω i (i = 1, 2) (2.1)
We introduce the functions µ and ρ for the viscosity and the density of the fluid, respectively defined on the whole domain Ω as follows:

µ = χ 1 µ 1 + χ 2 µ 2 , ρ = χ 1 ρ 1 + χ 2 ρ 2
where χ i is the characteristic function of the subdomain Ω i and µ i ,ρ i are the constant dynamic viscosities and densities of each fluid (i=1,2), respectively.

Using these notations and based on the definition ν of the constant kinematic viscocity of the fluid for all subdomains, equations (2.1) are rewritten then in reduced form:

       ∂u ∂t + (u • ∇)u -ν∆u + ∇p = f in each Ω i (i = 1, 2) divu = 0 in Ω (2.2)
where u = u i , p = p i and f = f i in Ω i (i=1,2).

These equations need to be endowed by adequate boundary and initial conditions which have been discussed in the previous chapter for single fluid. Besides, they are also endowed with interface conditions imposing: the continuity of the velocity and the balance of the normal stress with the surface tension across the interface Γ(t), namely:

[u] Γ = 0, [σ] Γ .n 1 = -γκn 1 (2.3)
where σ here is the stress tensor defined as:

σ = ν ∇u + t ∇u -pI d , n 1 is the unit normal vector to Γ(t), exterior to Ω 1 (t); [.] Γ denotes the jump of quality across Γ in the normal direction of n 1 , i.e. [.] Γ = .| Ω 1 -.| Ω 2
and γ > 0 is the constant surface tension coefficient along the interface; κ is the algebraic mean curvature of the interface, being positive if the interface curve/surface bends towards Ω 1 and negative otherwise.

Remark 2. The effect of the surface tension here is expressed in term of a localized force -γκn 1 at the interface which is also mentioned as the continuum surface force (CSF) technique in references [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF][START_REF] Chang | A level set formulation of Eulerian interface capturing methods for incompressible fluid flows[END_REF]. We will show that this localized force of surface tension is added to the righthand side of the variational formulation in case of bifluid flow (in Section 2.2).

The interface capturing by level set method

The level set method, introduced by Osher and Sethian in [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations-jacobi formulations[END_REF], has become very popular in the past few years for describing boundaries of domains or evolution of free surfaces and interfaces.

This approach appeared in Computational Fluid Dynamics with the study of the motion of two compressible gases, separated by a sharp interface in [START_REF] Mulder | Computing interface motion in compressible gas dynamics[END_REF] and was used in [START_REF] Sussman | A level set approach for computing solutions to in compressible two-phase flow[END_REF] for describing the interface between two immiscible incompressible fluids, driven by the Navier-Stokes equations.

Since these seminal works, level set method has been extensively used to treat problems involving free surfaces or interfaces between the fluids and strongly promoted in two books by Sethian [START_REF] Sethian | Level Set Methods[END_REF] and [START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF].

Following [START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF], we introduce the signed distance function to the interface Γ(t) as follows:

φ(x, t) = ±d(x, Γ(t)), ∀x ∈ Ω(t) (2.4)
where d(x, Γ(t)) denotes the usual Euclidean distance function from x to Γ(t) and:

             φ(x, t) > 0 if x ∈ Ω 1 (t) φ(x, t) < 0 if x ∈ Ω 2 (t) φ(x, t) = 0 if x ∈ Γ(t). (2.5)
See figures 2.2, 2.3 for illustrations of computational domain with the level-set function.

The function φ is initialized by the signed distance function to Γ 0 :

φ 0 (x) = ±d(x, Γ 0 ), ∀x ∈ Ω(0)
where Γ 0 is the initial position of the interface and indicates the initial shape of each subdomain Ω i (0). We use here the approach suggested in [START_REF] Sussman | A level set approach for computing solutions to in compressible two-phase flow[END_REF] for the incompressible two-phase flows. The interface between two fluids at each time t is then captured as zero-level set of φ(x, t) which is solution to advection equation:

     ∂φ ∂t (x, t) + u(x, t)∇φ(x, t) = 0 ∀(x, t) ∈ Ω × R + φ(x, 0) = φ 0 (x) ∀x ∈ Ω (2.6)
where u is the solution of (2.2) and (2.3). Notice that in these equations, the density and the viscosity are constant in each fluid and take different values depending on the characteristic function χ i of each subdomain Ω i (hence, depending on the sign of the level set function): Furthermore, in the level set framework, the unit normal vector n to the interface and the mean curvature κ at the interface are classically computed via the function φ as:

µ = χ 1 µ 1 + χ 2 µ 2 , ρ = χ 1 ρ 1 + χ 2 ρ 2
n = ∇φ |∇φ| | φ=0 ; κ = divn = div ∇φ |∇φ| | φ=0
However, we will introduce below in section 2.2 another geometric technique to approximate n and κ, that revealed less sensitive to numerical artifacts when the shape of Γ(t) becomes more complex.

The extension and regularity of flow fields

One should note that the evolution of the interface depends strongly only on the flow field in its vicinity. Moreover, in case of complex deformation, the large velocity may cause uncontrolled oscillation and jeopardize the numerical stability. This is the reasons we propose to extend and regularize the velocity of flows u (defined only on the interface Γ) into vector field u (defined on the whole Ω) for advection of the interface. There are many approachs for extension and regularity velocity, for example from suggest in [START_REF] Burger | A framework for the construction of level-set methods for shape optimization and reconstruction[END_REF] the vector u can be seached in a regular space V by solving the following problem:

             -α∆ u + u = 0 inΩ u = 0 onΣ u = u onΓ (2.7)
where α > 0 can be interpreted as a regularization lengthscale. This leads us to the variational problem: Find u ∈ V such that:

α Ω ∇ u.∇v + Ω uv = Γ uv ∀v ∈ V (2.8)
It can be seen that the left-hand side of (2.8) is a coercive bilinear form on V which is close to I (so u is expected close to u). This equation can be solved easily by a finite element method.

Description of the numerical scheme

The bifluid model is represented by the equations (2.2)-(2.7). In summary, the general scheme to solve this system is the following:

1. Initialization of the level set function φ 0 and velocity u 0 .

At each time step:

-update the interface Γ, as well as the densities ρ and the viscosities µ of fluids;

-solve the Navier-Stokes equations to find u, p;

-extend and regularize the velocity u to obtain a vector field for solving the advection equation;

-solve the advection equation to obtain φ in the next time.

3. Repeat step 2 until the final time.

The reduced scheme is hiding some complex and necessary numerical routines like: solution interpolation, error estimate and mesh adaptation. All these numerical techniques will be detailed in chapter 3. However, let us say few words about mesh adaptation, an important routine. Mesh adaptation aims at improving the accuracy of numerical solution and saving the memory resource.

In our approach, the adapted mesh will be generated in order to assure the minimal error estimation of velocity and geometrical approximation of the interface, it plays the role of a computational mesh for both the Navier-Stokes and advection equation in each time step. A more precise scheme of the whole three parts (Navier-Stokes, advection, mesh adaptation) will be presented in section 2.5.

The surface tension

Surface tension is very important in the simulation of fluid flows with free surface or interface.

In the case of interfaces between two fluids, the surface tension contributes a surface pressure that is a normal force per unit interfacial area. Actually, this surface tension balances the normal stress at the interface (see equation (2.3)). Intuitively, the surface tension balances the interaction between two fluids, if the interface is concave in one domain then the surface tension tends to reduce this concavity.

For an illustration of the effect of the surface tension, we consider the evolution of a rising bubble in a fluid. The bubble with lower density than surrounding fluid rises and final at the top of the domain. It has been observed the difference behavior of shape bubble when the surface tension term effects or not in figure 2.4. In the first case, without surface tension (γ = 0), the bubble rapidly becomes deformed in time while in the second case, under the effect of surface tension, the bubble is almost circular during the evolution.

As mentioned in the previous section, due to the condition on the interface (2.3) the surface tension appears in the right-hand side of the momentum equation as a localized force concentrated at the interface. This section will consider the effect of this term in the equation of bifluid flow (2.2).

After, we will treat the geometric approximation of this term using the explicit discretization of the interface in the computational mesh. (bottom) with surface tension.

The surface tension term in the variational formulation

In solving Navier-Stokes equations for two fluids flow, the discretization in time the leads us solve the following generalized Stokes equations in each time:

     ku -ν∆u + ∇p = g in each Ω i (i = 1, 2) divu = 0 in Ω (2.9)
where u = u i , p = p i and f = f i in Ω i (i=1,2) and k is denoting of 1/∆t.

By introducing D(u) = ∇u + t ∇u 2
, the symmetric gradient of u, also called the rate of deformation tensor, and thanks to the condition of incompressibility we have: div(2D(u)) = ∆u.

Using Green's formula we have:

-

Ω i div(2µD(u i ))v i dx = Ω i 2µD(u i ) : ∇v i dx - ∂Ω i 2µD(u i )n i .v i ds = Ω i 2µD(u i ) : D(v i )dx - ∂Ω i 2µD(u i )n i .v i ds Ω i ∇p i .v i dx = - Ω i p i divv i dx + ∂Ω i p i n i .v i ds
Given a test function v, taking the inner product in L 2 (Ω) d of the first equation of system (2.9) and summing on i, we have following equation:

k Ω u.vdx + Ω 2µD(u) : D(v)dx - Ω pdivvdx = Ω g.vdx + 2 i=1 ∂Ω i 2µD(u i )n i .v i ds - 2 i=1 ∂Ω i p i n i .v i ds (2.10)
Due to the interface condition (2.3) and the geometric condition n 2 = -n 1 on Γ and noticing also that v = 0 on Σ we have:

2 i=1 ∂Ω i 2µD(u i )n i .v i ds - 2 i=1 ∂Ω i p i n i .v i ds = 2 i=1 Γ (2µD(u i ) -p i I)n i v i ds = Γ (σ 1 n 1 v 1 -σ 2 n 1 v 2 )ds = Γ -γκn 1 .vds
This term corresponds to the effect of the surface tension on the interface in the variational formulation. As can be seen that the surface tension is taken into account as a localized force -γκn 1 added to the right-hand side of the momentum equation and proportioned to the local mean curvature of the interface.

The discretisation of the surface tension term

In the context of level set method, the level-set function allows to calculate the unit normal vector and the mean curvature at the interface by following formulas:

n = ∇φ |∇φ| | φ=0 ; and κ = divn = div ∇φ |∇φ| | φ=0 (2.11)
Introducing δ Γ is a Dirac δ-function with support on Γ acting on a smooth test function ϕ as follows:

Ω δ Γ ϕdx = Γ ϕds (2.12)
The surface tension force is then represented by, see also in [START_REF] Sussman | A level set approach for computing solutions to in compressible two-phase flow[END_REF] :

l Γ (v) : = - Γ γκn 1 .vds = - Ω γκδ Γ (φ)∇φ.vdx
Hence, we have the expression of the surface tension in the discrete formulation:

l Γ h (v h ) = K∈T h - K γκδ(φ)∇φ.v h dx (2.13)
However, this approach requires an approximation of the gradient of φ which is known that errorprone on unstructured triangulations. In addition, the main challenge is to compute the mean curvature κ which is related to a second order derivative of the level-set function.

Another approach has been used to compute the effect of surface tension based on the so-called

Laplace-Beltrami operator. The Laplace-Beltrami operator ∆ Γ is a function of the tangential gradient ∇ Γ and is defined for a given function f by:

∆ Γ f = ∇ Γ .∇ Γ f (2.14)
with ∇ Γ is the tangential derivative (along Γ):

∇ Γ f = ∇f -n Γ n Γ t ∇f = (I -n Γ n Γ t )∇f (2.15)
then the second order derivatives of curvature can be eliminated by the following formulation (see [START_REF] Gross | Numerical methods for two-phase incompressible flows[END_REF][START_REF] Grob | A finite element based level set method for twophase incompressible flows[END_REF][START_REF] Hysing | A new implicit surface tension implementation for interfacial flows[END_REF]):

-∆ Γ id Γ (x) = κ(x)n(x), x ∈ Γ (2.16)
Using this result, we see that the surface tension term can be rewritten as follows:

l Γ (v) = -γ Γ ∇ Γ id Γ .∇ Γ v (2.17)
Given an approximation Γ h of Γ, the localized force term l Γ (v h ) can be approximated by:

l Γ h (v h ) = -γ Γ h ∇ Γ h id Γ h .∇ Γ h v h (2.18)
The discretization of the surface tension by this approach has been detailed in [START_REF] Gross | Numerical methods for two-phase incompressible flows[END_REF], however the computation is rather complex as requires also the approximation of the normal vector according to the first formula of (2.11).

In our approach the interface is explicitly discretized in the triangulation T h , this representation

gives us an alternate technique to approximate the interface via a set of connected segments (faces in three dimensions). Denoting (x i ) 1≤i≤ns the set of ordered vertices along the discrete curve Γ h such that x i-1 , x i , x i+1 represent its three consecutive points and x 0 ≡ x ns , x 1 ≡ x ns+1 in case Γ h is closed curve, see figure 2.5. Using quadrature formula along each edge E of Γ h , the surface tension term can be rewritten as follows, for all v h ∈ V h :

l Γ h (v h ) : = - Γ h γκn 1 h .v h ds = E⊂Γ h |E| 2 x i ∈E γκ(x i )n 1 h (x i ).v h (x i ) = x i ∈Γ h γκ(x i )n 1 h (x i ).v h (x i ) E x i |E| 2 (2.19)
where the unit normal vector n 1 is computed from the approximation of the unit tangent vector

τ = (τ 1 , τ 2 ) t at each vertex x i ∈ Γ h : τ (x i ) = ------→ x i+1 x i-1 / ------→ x i+1 x i-1 , hence n(x i ) = (τ 2 (x i ), -τ 1 (x i )) t .
The mean curvature κ(x i ) is obtained as the inverse of the radius r(x i ) which can be computed via the following approximation [START_REF] Frey | Mesh generation. Application to finite elements[END_REF]:

r(x i ) = 1 4 ----→ x i x i-1 , ----→ x i x i-1 -n 1 (x i ), ----→ x i x i-1 + ----→ x i x i-1 , ----→ x i x i+1 -n 1 (x i ), ----→ x i x i+1
This technique can be extended straightforwardly to three dimensions, where the unit normal is then taken as the weighted average value of the unit normals of all triangles sharing vertex x i .

Other formula can be used to approximate r(x) or κ(x) due to explicit discretization of geometric of interface, see for instances [TBE + 01] or [START_REF] Girault | One time-step finite element discretization of the equation of motion of two-fluid flows[END_REF].

Numerical resolution of the level set advection equation

We sketch out in this section numerical resolution of the evolution of the interface. We recall that, at each time step, the moving of the interface between two fluids is represented by the zero isocontour of the level set function φ that is convected by the velocity field to find the new position of the interface. The advection of the level set function by the flow velocity u (may be extended and regularized) obeys the so-called advection equation:

∂φ ∂t (x, t) + u(x, t).∇φ(x, t) = 0 ∀ (x, t) ∈ Ω × [O, T ].
(2.20)

The discretisation in time by the method of characteristics

Consider the Cauchy problem for advection equation has been introduced in (2.20) : Given an initial function φ 0 (x) : Ω → R and a velocity field u(x, t) : Ω → R d defined on Ω, find φ(x, t) :

Ω × [O, T ] → R such that:      ∂φ ∂t (x, t) + u(x, t)∇φ(x, t) = 0 ∀(x, t) ∈ Ω × (0, T ) φ(x, 0) = φ 0 (x) ∀x ∈ Ω (2.21)
This problem is solved by following the characteristic curves of the fluid particules. Given a particule

x ∈ Ω at the time s, its curve is described by the following equations:

     dX(x, s; t) dt = u(X(x, s; t), t) ∀t ∈ (0, T ) X(x, s; s) = x (2.22)
where X(x, s; t) is the position of x at the time t.

The first equation of (2.21) implies that φ(x, t) is constant along the characteristic lines X(x, s; t).

Hence the solution of the Cauchy problem (2.21) can be expressed as:

φ(x, t) = φ 0 (X(x, t; 0), 0) ∀(x, t) ∈ Ω × [0, T ] (2.23)
Assume that the interval [0, T ] is divided into a finite number intervals ∆t of the form (t n-1 , t n ) with t n = n∆t, we have the discretization in time of the equations (2.22) for all n:

     dX(x, t n ; t) dt = u(X(x, t n ; t), t) ∀t ∈ (t n-1 , t n ) X(x, t n ; t n ) = x (2.24)
We only compute φ(x, t n ) for all n, so if denote φ(x, t n ) by φ n (x), by substituting the time interval

[t n-1 , t n ] into (2.
23), we obtain the following result:

φ n (x) = φ n-1 (X(x, t n ; t n-1 )) ∀x ∈ Ω (2.25)
where X(x, t n ; t n-1 ) is the position at the time t n-1 of the characteristic emerging from x at the time t n .

The expression (2.25) follows us to calculate φ n from φ n-1 . In numerical solution, this problem can be solved by a Galerkin numerical scheme which involves the resolution of a linear system with using the quadrature formulas for approximating the integrals in the right-hand side (see [START_REF] Bermudez | Numerical analysis of convectiondiffusion-reaction problems with higher order characteristics/finite elements. part ii : fully discretized scheme and quadrature formulas[END_REF] for example). Another alternative is simply combination with a Lagrange interpolation, i.e φ n (x) is computed by values of φ n-1 (x) at the degrees of freedom of element K which locates X(x, t n ; t n-1 ).

In fact, the first approach is more expensive than the second one. Moreover, the second approach brings an estimation to control the geometric error of the interface by the interpolation error (see section 3.3) that is very efficient in generation of mesh adaptation, so, its computation is applied for present work.

Algorithm 3: Numerical scheme for advection equation n=0, start with the initializations (φ 0 , Γ 0 ); for n=1,...,N do 1. Solve the ODE (2.24) to search X n-1 (x) for each node x of mesh;

2. Compute φ n-1 (X n-1 (x)) for each node x of mesh;

3. Update the new level set values φ n (x) according to (2.25). end

Spatial approximation of the characteristic points

Let us now explain the approximation of the characteristic lines, or more precisely the approximation of X(x, t n ; t n-1 ) in space. We notice that in numerical resolution only an approximation u h of u is known at each degree of freedom of computational mesh T h of the domain, and we need to compute an approximation of X h (x, t n ; t n-1 ) where X h (x, t n ; t n-1 ) is solution at the time t n-1 of the approximated characteristic curve:

     dX h (x, t n ; t) dt = u h (X h (x, t n ; t), t) X h (x, t n ; t n ) = x
(2.26) that implies the "formal" expression:

X h (x, t n ; t) = x - t n t u h (X h (x, t n ; t), t)dt (2.27)
There are many approaches to approximate X h (x, t n ; t n-1 ). For example the simplest algorithm for computing the characteristics is obtained directly from (2.27) as:

X h (x, t n ; t n-1 ) = x -∆tu h (x) (2.28)
In this way, the characteristic curve is considered to be a straight line from the starting point x to its foot X h (x, t n ; t n-1 ). But in fact, this algorithm should be used only when the time step is small enough, one should not travel too many elements mesh at once to ensure the numerical stability.

So, we prefer to use Runge Kutta 4 or Euler's scheme, especially in case the time step is relatively large. Given a substep δt of each interval [t n-1 , t n ] and suppose that existing an integer number M such that: ∆t = M δt, we have the backing forward schemes as follows:

1. Euler's scheme:

             X h (x, t n ; t n ) = x X h (x, t n ; t n -mδt) = X h (x, t n ; t n -(m -1)δt) -δtu h (X h (x, t n ; t n -mδt)); (m = 1, ..., M ) (2.29)
2. Runge-Kutta 4 scheme: See figure 2.6 for illustration of characteristic curve approximation.

                                   X h (x, t n ; t n ) = x X h (x, t n ; t n -mδt) = X h (x, t n ; t n -(m -1)δt) - δt 6 (v 1 + 2v 2 + 2v 3 + v 4 ), with : v 1 = u h (X h (x, t n ; t n -mδt)) v 2 = u h (X h (x, t n ; t n -mδt) - δt 2 v 1 ) v 3 = u h (X h (x, t n ; t n -mδt) - δt 2 v 2 ) v 4 = u h (X h (x, t n ; t n -mδt) -δtv 3 ); (m = 1, ...M ) (2.30)
Remark 3.

1. At now, two different time steps have been used. The first one, ∆t = t n -t n-1 , is a "large" time step on a generic period of time [0, T ], that is the period of time for which we assume the velocity and physical properties of the fluids are unchanged. Time step ∆t is then used for both resolution of Navier-Stokes and advection equation. The second one, δt < ∆t is a sub-time step involved the integration of the characteristic curves. In fact, the approximation of characteristic curves with the hold time step ∆t means that the curves (x, X h (x, t n ; t n-1 )) is approximated by straight line from the starting point x to its foot X h (x, t n ; t n-1 ) and directed only by u(x), this may lead to inaccuracy dramatically. By using δt we achieve a polygonal approximation with the velocity field is updated more frequently, so we can improve the approximation of the characteristic curve.

2. The approximation by these schemes needs the information of the velocity at the immediate points as u h (X h (x, t n ; t n -mδt)), ... that always requires the exact location of these points.

Normally, given a point x located at element K, the immediate points is somewhere in K or adjacents of K. Therefore, the efficiency of the algorithm is related to the fast identification of the triangle adjacent to a given triangle. To this end, we use a specific data structure in which the three adjacent triangles are associated with the informations of each triangle (the so-called adjacency matrix [START_REF] Frey | Mesh generation. Application to finite elements[END_REF]).

3. Some time, the procedure of searching the characteristic point can not implemented for all the interval (t n-1 , t n ), especially with the vertices next to the boundary of domain. This implies that during the procedure, one of immediate points is go out to domain (see figure 2.7). Although, in that case these vertices need to be traveled with the maximal time that possible, may be this leads to the last point hit the boundary. In fact, this requires one algorithm of calculating the time dt that is necessary for one point to hit an edge or go out the element contains it, see 1. Calcul the barycentric λ i(1≤i≤3) of point P .

2. Calcul the barycentric η i(1≤i≤3) of point P -u(P ).

3. for i=1,...3 do If (λ i > η i and dt > λ i λ i -η i ) then dt := λ i λ i -η i end 2.

Redistancing procedure and conversation the mass

It has been known that, in the context of level set method, the level set function must usually satisfy:

|∇φ| = 1 (2.31)
Unfortunately, when φ is transported by a physical velocity field, this property is not preserved as the isolines do not travel at the same speed and the level set does not remain a distance function.

A natural choice to reinitialize the level set function is the signed distance function to the interface:

φ(x) =          d(x, Γ) if x ∈ Ω 1 0 if x ∈ Γ -d(x, Γ) if x ∈ Ω 2
In our scheme, this signed distance function is computed approximatively by redistancing procedure studied in [START_REF] Dapogny | Computation of the signed distance function to a discrete contour on adapted triangulation[END_REF], we only briefly recall here the implement of this procedure:

-Step 1: Initialization φ 0 of φ: denoting T Γ the set of mesh elements intersected by the interface, i.e. T Γ = {K ∈ T h : K ∩ Γ = ∅}, φ 0 (x) is defined as:

φ 0 (x) =          φ(x) if x ∈ T Γ +∞ if x ∈ Ω 1 \ T Γ -∞ if x ∈ Ω 2 \ T Γ -Step 2:
Calculation numerically φ as steady solution of so-called Eikonal equation:

     ∂φ ∂t (x, t) + sgn(φ 0 )(|∇φ| -1) = 0 ∀(x, t) ∈ Ω × (0, T ) φ(x, 0) = φ 0 (x) ∀x ∈ Ω (2.32)
Its purposes are to remain the behaviour in "vicinity" of the zero isoline, i.e. the position of interface Γ is not modified and to ensure the constraint (2.31).

As pointed out, other drawbacks of this method in case of free surface problems involving incompressible flows is the lack of mass conservation. In this work no particular treatment for mass conservation during the simulation, this problem is remedied by replacing zero-level set by -level set with is appropriate negative constant and depends strongly on quality of computational mesh.

It has been observed that without this correction of level set function, the derivation of mass is not notable when our mesh adaptation asks a lot of elements at each iteration, however, for the efficiency of computational time especially intending to the problem in 3D, this cure is applied in our approach.

The proposed scheme

In this section we will describe the general scheme for all the time [0, T ]. Suppose that [0, T ]

is divided by N subintervals [t n-1 , t n ].
At n-iteration, we have to solve the Navier-Stokes equation and the advection equation approximatively, emphasize that they are solved on the same mesh T n to obtain the approximative solution (u n , p n ) and φ n .

Our scheme is an iterative procedure based on the mesh adaptation process. At each time t n , given mesh T n , initialization u n 0 for velocity and φ n 0 for level-set function. After solving Navier-Stokes problem we have the new velocity field u n , this velocity (may be extended and regularized) is used for solving the advection equation to obtain a new level-set φ n , the new interface Γ n is captured as the zero-level of φ n . The mesh T n is adapted with u n , φ n 0 , φ n in order to create a new mesh T n+1 . The generation of adapted mesh at each time step will be detailed in chapter 3. The projection of u n and φ n on T n+1 are implemented to have initialisation u n+1 0 and φ n+1 0 for the (n + 1)-iteration.

The overall algorithm is given as following:

Numerical scheme for bifluid flow over [0, T ] 1. Start with mesh T 1 and initialization u 1 0 , φ 1 0 2. For n = 1, ..., N do:

Problem Input Output 2.1 Solving Navier Stokes (T n ,u n 0 ) u n 2.2 Solving Advection (T n , u n , φ n 0 ) φ n

Extension and regularity

(T n , u n | Γ ) u n 2.4 Adaptation (T n , φ n , φ n 0 , u n ) T n+1 2.5 Redistancing (T n+1 , φ n ) φ n 2.6 Projection (u n , φ n , T n+1 ) u n+1 0 , φ n+1 0 3. Return (u N , φ N , T N )
In comparison with previous study in [START_REF] Bui | A coupling strategy based on anisotropic mesh adaptation for solving two-fluid flows[END_REF] this scheme has been much reduced. In the before approach mesh adaptation needed two mesh independent, one for the resolution of the fluid and other for the advection of the interface, so, the interpolations is always necessary to correspond to each solving of the problem. The simplicity of present scheme is due to the requirement only one adapted mesh at each time of the simulation, yet all difficulties have been driven for mesh adaptation. On the one hand, the adapted mesh has to respond the capturing of the interface with high quality, especially in the vicinity of interface. On the other hand, the solving Navier-Stokes equation requires the accuracy of numerical solution of velocity of fluids and the efficiency of the numerical scheme depends strongly on the quality of the computational mesh. However, we are only interested in the vicinity of the interface, so mesh adaptation is only much concerned in this "important" region of the computational mesh, resulting actually our mesh adaptation is implemented rapidly and efficiently.

Conclusion

We have presented all the numerical scheme for bifluid flows simulation. Our scheme is based on unstructured adapted mesh at each time step. This strategy has several assets; on the one hand, no projection is needed between different meshes. In addition, the method does not present any theoretical difficulty to the extension from the two-dimensional case to the three-dimensional case.

This is a tremendous feature insofar as mesh adaptation is concerned; in fact mesh generation is known to be more difficult to deal with in three dimensions. Most of the difficulty related to meshing is recast as a robust remeshing problem. We go to the next chapter to detail the mesh adaptation process for the previous scheme.

Chapter 3

Mesh adaptation This chapter is devoted to mesh adaptation. It is known that in numerical simulation based on finite element methods, the accuracy of the solution depends strongly on the quality of the computational mesh. The purpose of mesh adaptation is to increase the efficiency of the numerical schemes and improve the accuracy of numerical solutions. This process relies on: 65 CHAPTER 3. MESH ADAPTATION -the control of a priori/ posteriori error approximations (interpolation error, geometry error,...) of numerical solutions.

-the quality assessment of mesh elements respect to an appropriate quality function.

-the concentration of more degrees of freedom in important regions (such as boundary/interface or regions of large solution variations) than in any other region of the computational domain.

Consequently, the adapted mesh enjoys a nice quality with fewer elements that ultimately results in a diminution of the computational cost to achieve the desired accuracy.

We will introduce hereafter two mesh adaptation procedures which have been used in our numerical experiment.

1. Anisotropic mesh adaptation. This is the inheritance of a mesh adaptation technique that has been developed for many years, see [AF03, FA05, Fre08] and applied in many previous works for the model of level set advection equations, see [START_REF] Claisse | Levelsets and anisotropic mesh adaptation[END_REF][START_REF] Bui | An accurate anisotropic adaptation method for solving the level set advection equation[END_REF] and the model of Stokes bifuid flow in 2D [START_REF] Bui | A coupling strategy based on anisotropic mesh adaptation for solving two-fluid flows[END_REF]. In this thesis, this anisotropic mesh adaptation is ameliorated to solve the model of Navier-Stokes bifuid flow in two dimensions.

2. Discrete three dimensional domain remeshing. This is a tremendous developing of mesh adaptation since mesh generation is known to be more difficult to deal with in three dimensions.

Most of the difficulty related to meshing is recast as a robust remeshing problem. Thanks to the results of this mesh adaptation process (see [START_REF] Dapogny | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF]) several test cases of bifluid flows in three dimensions have been implemented efficiently, as will be seen later in chapter 4.

The outline of this chapter is the following: In section 3.1 we recall some basic notations and definitions of mesh adaptation. Section 3.2 and section 3.3 describe the two previous mesh adaptation procedures and at the end of each section we will explain how to construct the adapted mesh for two fluids flow problem in chapter 2.

Basic notations and definitions

Isotropic and anisotropic mesh adaptation

Definition 1. Let Ω be a open, bounded, polygonal domain in R d (d=2,3) and given its simplicial mesh T (i.e triangulation in two dimensions, tetrahedralization in three dimensions). We assume that each element K in T is a closed d-simplex (triangle in two dimensions, tetrahedra in three dimensions) and Ω = K∈T K. T is called a conforming mesh if satisfies the so-called conformity condition such that the intersection of any two different simplices, if not empty, is reduced to:

-either a point, or a common edge in two dimensions, -either a point, or a common edge or a common (triangular) face in three dimensions.

See figure 3.1 for examples of conforming and non conforming meshes.

Given T is a conforming mesh, then T is: -a uniform mesh, if all its simplices are equally sized and regular (equilateral), -a quasi-uniform mesh, if the variation of its simplices size is bounded and there exists a constant c such that h K ρ K ≤ c for all elements K in T with h K is diameter of longest edge of K and ρ K is radius of its inscribed circle or sphere. -Isotropic mesh adaptation generates the mesh where almost elements are regular (equilateral) and only adjusted in size based on an error estimate.

-Anisotropic mesh adaptation based on a metric tensor consists in controlling the size, the shape and the orientation of the elements altogether. -Both anisotropic and isotropic mesh adaptation have been applied in practice. Traditionally, isotropic mesh adaptation has received much attention, however, in regions of large solution gradient, adaptive isotropic meshes usually contain too many elements. Moreover, in several cases elements with high aspect ratio are needed for better represent of the solution variations or boundary layers in fluid dynamics, thus it is desirable to adjust the element size as well as the element shape, resulting in an anisotropic mesh.

-From the geometrical point of view, isotropic mesh adaptation can be considered as a special case of anisotropic one with a special metric tensor (see below in section 3.2.1).

The quality of a mesh

In addition to the requirement of adapting mesh to approximate the continuous domain Ω, we also need to define criteria to evaluate the quality of an element K of a simplicial mesh T .

There are many definitions of the quality of a simplex K ∈ T . For instance one can rely on the ratio ρ K h K , see [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]. This measure only depends on the shape of the considered simplex and not on its size. From practical point of view, this implies that the accuracy of a finite element computation performed on T is of course influenced by the size of its elements, but also by their "well-shapeness".

However, we will rather rely on the following quality function Q(K) of a d-simplex K:

Q(K) := α V (K) (Σ na i=1 l(e i )) d/2 (3.1)
where V (K) is the area/volume of K, na = d(d + 1)/2 is the number of edges of K denoted by e 1 , . . . , e na and l(e i ) stands for the length of e i . The quality function in (3.1) retains the same theoretical meaning as ρ K h K , but shows a better numerical ability when it comes to discriminating "good" from "average", or "bad" elements. One can show that, for any tetrahedron K, Q(K) ≤ 1, and equality holds if and only if K is regular.

Anisotropic mesh adaptation

Research on anisotropic adaptation is based on the idea of metric-based mesh adaptation: the local desired size, shape and orientation of the elements are prescribed using a metric tensor field.

Usually, an anisotropic tensor field can be defined from an error indicator or an error estimate that relates the approximation error for a size, shape and orientation prescriptions, see [BA76, VHM91, Ape99, FP01, DVB + 02, Hua05] and the references therein for examples. Theoretical analysis of this approach in our program of mesh adaptation have been detailed in [START_REF] Alauzet | Estimateur d'erreur géométrique et métriques anisotropes pour l'adaptation de maillage. partie i: aspects théoriques[END_REF][START_REF] Frey | A differential approach to mesh generation, volume 9 of Differential Geometry: Theory and Applications in[END_REF][START_REF] Frey | Mesh generation. Application to finite elements[END_REF]. We only recall here some basic notations related to adaptation process in this thesis.

Metric tensor fields

l M (γ) = 1 0 < M (γ(t))γ (t), γ (t) > dt V M (K) = K det(M (x))dx d M (x, y) = inf γ∈C 1 ([0,1],R d ) γ(0)=x,γ(1)=y l M (γ)
From the geometrical viewpoint, the metric defining an element can be represented by an ellipsoid whose volume, ratios between the lengths of semi-axes and the principal axis vectors are associated respectively with notations of the size, the shape and the orientation of the element. Indeed, given a metric tensor M , let T a unit mesh with respect to M . Given x 0 is a vertex of T such that M (x) = M is almost constant around x 0 , then every simplex K of T lying in the ball B(x 0 ) (implies that K sharing x 0 ) is inscribed in the ellipsoid B M (x 0 )) defined as:

B M (x 0 ) = {x ∈ R d : d M (x, x 0 ) = 1} = {x = Σ d i=1 x i e i ∈ R d : λ 1 x 1 2 + . . . + λ d x d 2 = 1} (3.2)
where e 1 , . . . , e d denote the normalized eigenvectors of M (x 0 ) and λ 1 , . . . , λ d are the associated eigenvalues.

The set in (3.2) represents an ellipse (in 2D case) or an ellipsoid (in 3D) centered at x 0 , its axes are the eigenvectors e i and the length of the semi-axes are 1/ √ λ i (see figure 3.4). Notice that when all the λ i are equal, M (x 0 ) = λI d , where I d is the identity matrix, then the ellipsoid B M (x 0 ) becomes a sphere with radius 1/ √ λ, thus it accounts for isotropic mesh. On the contrary, in case of anisotropic mesh, size, shape and orientation notions of an element are associated with the volume, the ratios between the lengths of semi-axes and the principal axis vectors of ellipsoids, respectively. In cases of very stretched elements may be desired which do not fulfill the standard quality requirements follow definition 3.1, these quality functions must be traded for their anisotropic counterparts, obtained by using the same expressions except that the distance and volume notions are those supplied by metric tensor M , see quality function in expression (3.5).

Metric definition based on interpolation error

We present in this section the determination of metric tensor such that the according anisotropic mesh guarantees a given interpolation error. This determination will be used to construct a metric tensor for adapted mesh in the scheme of two fluids flow (see section 3.2.6).

Classically, the approximation error is bounded from above by the interpolation error, thus it is natural to get an upper bound of the interpolation error. Let us now recall the L ∞ error estimate for the Lagrange finite element P 1 -interpolation in [START_REF] Frey | Anisotropic mesh adaptation for CFD computations[END_REF]. For a function ϕ : R d → R, its linear interpolation is denoted by Πϕ, for all mesh element K we have following error estimate:

||ϕ -Πϕ|| L ∞ (K) ≤ c d max x∈K max -→ v ⊂K < - → v , |H(ϕ)(x)| - → v > (3.3) ≤ c d max x∈K max -→ e ⊂E K < - → e , |H(ϕ)(x)| - → e >
where H(ϕ) is the Hessian matrix of function ϕ, E K is the set of edges of the element K and c d is a constant depending on the dimension d. This estimate indicates that controlling the size of element edges allows to control the interpolation error. Given ε > 0 is desired tolerance, we aim at defining a metric tensor such that: ||ϕ -Πϕ|| L ∞ (K) ≤ ε for every element K in the adapted mesh.

As suggested in [START_REF] Alauzet | Estimateur d'erreur géométrique et métriques anisotropes pour l'adaptation de maillage. partie i: aspects théoriques[END_REF], we define an anisotropic metric as follows:

M ϕ = P t ΛP ; Λ =      λ 1 • • • 0 . . . . . . . . . 0 • • • λ d      (3.4) with λ i = min max c d |λ i | , 1 h 2 max , 1 h 2 min
where P is the eigenvector matrix and the coefficients λ i are the eigenvalues of the Hessian matrix H(ϕ), h min and h max are the prescribed minimal and maximal edge size of the mesh. This metric construction is related to the Hessian of the exact solution ϕ, that is however not known. Therefore, a procedure for approximating the Hessian of the solution from the discrete solution is necessary, see [START_REF] Alauzet | Estimateur d'erreur géométrique et métriques anisotropes pour l'adaptation de maillage. partie i: aspects théoriques[END_REF].

Metric intersection

In several applications, we need to adapt a mesh at the same time to several priori independent informations, supplied by two (or more) metric tensor fields M 1 , M 2 (for examples the mesh adaptation for two fluids flow, see below in section 3.2.6). This is classically achieved by a so-called metric intersection procedure: operating on the simultaneous reductions of M 1 (x) and M 2 (x) at any point

x ∈ R d . Suppose that M 1 (x) = t P (x)      λ 1 (x) . . . 0 . . . . . . . . . 0 . . . λ d (x)      P (x) and M 2 (x) = t P (x)      ν 1 (x) . . . 0 . . . . . . . . . 0 . . . ν d (x)      P (x)
where P (x) is an invertible matrix; λ i > 0, ν i > 0(i = 1, . . . , d) then the intersected metric is determinated by:

M 1 ∩ M 2 (x) := t P (x)      max(λ 1 (x), ν 1 (x)) . . . 0 . . . . . . . . . 0 . . . max(λ d (x), ν d (x))      P (x)
Geometrically, this intersected metric implies that the ellipsoid B M 1 ∩M 2 (x) is the largest one in-

The generation of anisotropic mesh

We assume that a metric tensor M is defined at the mesh vertices of a given simplicial T (note that in practice, M (x) is defined only at the nodes of any background structure and interpolated from these values, see [START_REF] Frey | Mesh generation. Application to finite elements[END_REF]). Several techniques have been devised for generating anisotropic meshes according to a metric tensor, such as Delaunay based triangulation methods (see [START_REF] Shewchuk | Lecture notes on Delaunay mesh generation[END_REF] for instances) and advancing-front methods where the construction of mesh elements is from the surface mesh of its boundary, see [START_REF] Frederick | Two-dimensional automatic mesh generation for structural analysis[END_REF][START_REF] George | Computer Implementation of the Finite Element Method[END_REF]. Another approach is the local mesh modification method which starts from an existing non-adapted mesh and adapt it so that it fits at best conditions, see [START_REF] Dobrzynski | Anisotropic Delaunay mesh adaptation for unsteady simulations[END_REF] for a description of mesh procedure. In our strategy, the generation of an adapted mesh is obtained using a Delaunay-based local mesh modification procedure (described in [START_REF] Dobrzynski | Anisotropic Delaunay mesh adaptation for unsteady simulations[END_REF]) include following classical operations:

-splitting of mesh edges, -collapsing of mesh edge endpoints, -swapping of mesh edges or faces, -vertex relocation.

Anisotropic mesh adaptation aims at modifying T iteratively by local operations in order to get a uniform (or a quasi-uniform) mesh with respect to this metric: in the desired mesh, all its simplices have edges lengths equal to 1 (respectively, lying in [l min ,l max ]) in the sense of Definition 3 and the anisotropic quality measure:

Q M (K) := α V M (K) (Σ na i=1 l M (e i )) d/2 (3.5)
of its simplices are as close to 1 as possible.

Explicit discretization of the interface

In the context of anisotropic level set adaptation, an interface Γ is approximated by a piecewise affine interface Γ T such that d(Γ, Γ T ) < ε on the adapted mesh T . We recall that Γ is defined basing on zero level set of function φ: Γ φ = {x ∈ Ω : φ(x) = 0}. For practical reasons, we need to modify T in order to obtain a new triangulation T in which Γ T is explicitly discretized. This is achieved by inserting the set of intersection points of Γ φ with the mesh edges/surfaces of T into the set of mesh vertices of T . The level set metric is updated on the resulting mesh and this mesh is then modified as previously to obtain final mesh T with respect to the updated metric.

Anisotropic mesh adaptation process for two-phase flow

We are now in position to construct an adapted mesh based on anisotropic mesh adaptation for the problem two fluid in section 2.5. In each time step, the initial mesh T have to solve two following problems:

1. Solving Navier-Stokes equation to have the vector field u. This resolution is based on method of characteristic which is known that unconditional stable with the following error estimate, see [START_REF] Pironneau | On the transport-diffusion algorithm and its applications to the Navier-Stokes equations[END_REF]:

u -u h L 2 (Ω) ≤ c(h m + ∆t m + h m+1 /∆t) (3.6) 2.
Find the new position of the interface Γ as zero iso-contour of level set function, solution to advection equation:

     ∂φ ∂t (x, t) + u(x, t)∇φ(x, t) = 0 ∀(x, t) ∈ Ω × R + φ(x, 0) = φ 0 (x) ∀x ∈ Ω
We aim to alter the initial mesh T to obtain a new mesh T such that:

i. The interface Γ is explicitly discretized the mesh T .

ii. The geometry error of the zero-level set and its approximation in the context of Hausdorff distance is small enough, i.e d H (Γ, Γ T ) < ε with ε is the desired error and:

d H (Γ, Γ T ) := max sup x∈Γ inf y∈Γ T |x -y|, sup y∈Γ T inf x∈Γ |x -y|
iii. The mesh T is nice quality according to quality function defined in (3.5) and responds the resolution of the Navier-Stokes equation at the next iteration with small error approximation.

It can be interpreted that:

i. The first requirement is assured following the lines in section 3.2.5.

ii. The second requirement is implemented by the construction a metric tensor as the following:

We have show in [START_REF] Bui | An accurate anisotropic adaptation method for solving the level set advection equation[END_REF] the following error estimates analysis for advection equation:

-The interpolation error of level set function and its numerical solution:

φ -φ h L ∞ (Ω) ≤ φ -Π h φ L ∞ (Ω) + φ 0 -Π h φ 0 L ∞ (Ω) + c 1 u -u h L ∞ (Ω) + c 2 e δt δt. (3.7)
where c 1 , c 2 are constants depending on initial datas φ 0 and velocity u.

-The geometry error of zero-level set and its approximation in the context of Hausdorff distance:

d H (Γ, Γ h ) ≤ sup   sup x∈Ω ∇φ(x) inf x∈Ω ∇φ(x) 2 , sup K∈T h ∇φ h (x)| K inf K∈T h ∇φ h (x)| K 2   φ -φ h L ∞ (Ω) ( 3
-Π h φ L ∞ (Ω) , φ 0 -Π h φ 0 L ∞ (Ω)
and the approximation error u -

u h L ∞ (Ω)
. In other words, the initial mesh is adapted such that these errors are as small as necessary. Otherwise, it is well-known that in case of approximation by Galerkin finite element, thanks to Cea's lemma the approximation error u-u h L ∞ (Ω) is bounded by the interpolation error u -πu L ∞ (Ω) . In summary, in order to perform a good discrete approximate Γ h of Γ the desired mesh should be adapted to interpolation of functions: φ, φ 0 and u.

The last term in the right-hand side of (3.7) is related to substep δt, this means that the geometry error can be controlled by the adjusting of this substep.

In the section 3.2.3 we have detailed how to construct a metric tensor in such a way that the generated mesh satisfies the given interpolation error. In our scheme, we intend to construct only one mesh adapted to both Lagrange interpolation of function φ 0 , φ and u. Denoting M φ 0 , M φ and M u are respectively the resulting metric tensor according to these interpolations.

By using intersection metric (see section 3.2.4) two times we obtain a tensor metric M = (M φ 0 ∩ M φ ) ∩ M u . The desired mesh is built according to metric tensor M .

iii. The third condition is responded rely on the estimation (3.6). It relates to characteristic mesh elements size h and the time step ∆t. This implies that the approximation error of Navier-Stokes problem can be decreased by diminution h as well as chosen the appropriate time step ∆t.

Hence, the final adapted mesh is an anisotropic mesh generated according to metric tensor M following the lines in section 3.2.4 with the appropriate characteristic element size h adjusted by [l min , l max ]. See figure 3.6 for the generation of adapted mesh at each time step. This illustration is one step of the simulation of the rising bubble in 2D (section 4.2).

Discrete three dimensional domain remeshing

We now address the problem of isotropic three-dimensional domain remeshing which has been used in mesh adaptation of some numerical tests in 3D. As in the previous section, we keep always the notations Ω of continue domain and its three-dimensional simplicial mesh T associated with "surface triangulation" Γ T which is an approximation of surface part Γ of Ω. The aim is to improve mesh T to obtain a resulting mesh T which is expected to be nice approximation of domain Ω follows some criterions, for examples:

-The surface triangulation Γ T is the better approximation of Γ.

-T is well-shaped according an evaluating of quality mesh, for instances in this strategy we use the following nomalized quality mesh function element for any tetrahedron K ⊂ R 3 :

Q(K) := α V (K) (Σ 6 i=1 l(e i )) 3/2 , with α = 144 √ 3 (3.9)
where e i , i = 1, . . . , 6 are the edges of K and l(e i ) stands for the length of e i . Explicit discretization of zero-level set of φ n , the obtained mesh is not well-shaped, (d): High-quality adapted mesh in which the zero-level set of φ n is explicitly discretized.

.

The local size map

The remeshing process of T is based on the classical local remeshing operators (edge split, edge collapse, edge swap, and node relocation) which enjoy two forms, depending on whether they are applied to a surface configuration or to a purely internal one, see [START_REF] Frey | Mesh generation. Application to finite elements[END_REF][START_REF] Dapogny | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF]. Although, we need a global vision to drive the remeshing strategy, that is to identify (or classify) those edges of T that should be split, collapsed, or swapped. Since [START_REF] Vallet | Anisotropic control of mesh generation based upon a Voronoï type method[END_REF][START_REF] Frey | About surface remeshing[END_REF], a very convenient means to encode such information has been rely on a size map h : Ω → R such that for any x ∈ Ω, h(x)

accounts for the local desired size for the edges of T surrounding x. The final aim of the process is then to produce a new mesh T of Ω, whose edges pq have (as far as possible) unit length l h (pq) respect to h, where:

l h (pq) = 1 0 |pq| h(p + t(q -p)) dt (3.10)
In numerical practice, h is defined and stored at the vertices of T then interpolated from these data whenever needed elsewhere.

Map size determination adapted to the geometric approximation

This section will present the definition of size map h such that the adapted mesh according to h satisfies the desired geometry error of a given surface.

Consider T is a conforming tetrahedral mesh of a continue domain Ω with the interface surface Γ.

Let denote Γ T is approximation of Γ. Our attempt to modify T into a new mesh T in which the approximating interface surface Γ T of surface Γ is demanded should be close to Γ up to a given tolerance ε, that is:

d H (Γ T , Γ) ≤ ε (3.11)
where d H (Γ T , Γ) is Hausdorff distance between Γ T and Γ.

In practice, T will be generated rely on the information of the size map h which is defined as follows:

i. at a surface vertex x ∈ Γ T , the size prescription h(x) in a neighborhood of x based on the following theorem (see [START_REF] Dapogny | Shape optimization, level set methods on unstructured meshes and mesh evolution[END_REF] for a proof, Theorem 8.2):

Theorem 1. Let Ω ∈ R d a domain, and T a mesh, whose associated surface mesh Γ T is "close" from Γ. Denote as d Γ the signed distance function to Γ, and H(d Γ ) its Hessian matrix. Then

d H (Γ, Γ T ) ≤ 1 2 d -1 d 2 max K∈Γ T max x∈K max -→ e ⊂E K - → e , |H(d Γ )(x)| - → e (3.12)
Hence, to satisfy (3.11) it is enough to ask that, for each triangle K ∈ Γ T , one has:

2 9 max x∈K max -→ e ⊂E K - → e , |H(d Γ )(x)| - → e (3.13)
Now, let K ∈ Γ T be a triangle which is "very close" to a point x ∈ Γ, so that we can assume H(d Γ ) is nearly constant around K. An approximation of the latter sufficient condition is then:

∀ - → e ⊂ E K , 2 9 - → e , |H(d Γ )(x)| - → e ≤ ε
On the other hand, introducing κ i (x), (i = 1, 2) the two principal curvatures of Γ at x and e i (x) the two associated principal directions, the matrix H(d Γ )(x) writes in the orthonormal basis (e 1 (x), e 2 (x), n(x)) of R d :

H(d Γ )(x) =      κ 1 (x) 0 0 0 κ 2 (x) 0 0 0 0      (3.14)
Then, we have the upper bound

9ε 2max(|κ 1 (x)|,|κ 2 (x)|
) of any edges of K. We are led to the definition of metric size map h on the surface as

h : Γ → R ∀x ∈ Γ, h(x) = min h max , max h min , 9ε 2max(|κ 1 (x)|, |κ 2 (x)|) (3.15)
where h min and h max are respectively lower and upper bounds on the authorized lengths of edges in T .

ii. at an internal point x, no constraint from the geometry, i.e we only need apply the maximal authorized size h(x) = h max for an edge of the resulting mesh T .

Remark 5. The size map h defined above only accounts for the local size feature associated to the geometric approximation of surface Γ. In many case, the local size feature need to adapted another size map k : Ω → R arise from a posteriori error analysis associated to some numerical resolution of a problem (see the requirement of mesh adaptation in chapter 2 of this thesis or the problem of [START_REF] Ainsworth | A posteriori error estimation in finite element analysis[END_REF] for instances). In the problem of two fluids flows, the local size map need to adapted size map k arise from the error analysis associated to numerical resolution of Navier-Stokes equation problem.

In such a case, this additional information is taken into account by trading h for a new size map h defined as: ∀x ∈ Ω : h(x) = min(h(x), k(x)).

Gradation of the size map

Unfortunately, conforming to the size prescription discussed above is not sufficient in itself to guarantee the resulting mesh will enjoy a nice mesh quality. Indeed, the computed size map h may vary very sharply from one point to one of its neighbors, because the computation of h suffered from noise one the input data Γ T , or because the ideal surface Γ itself shows sharp variations of curvatures. This shock of size prescriptions between close areas may urge the formation of undesired ill-shaped elements during the remeshing process. For this reason, it may be desirable to drive the remeshing operators in such a way that the resulting triangulation Γ T from the process shows smooth variations in edge lengths. More accurately, we expect the resulting mesh to be such that two edges ap and bp sharing a common vertex p have Euclidean lengths satisfying:

1 r ≤ |b -p| |a -p| ≤ r (3.16)
where r is a user-defined bound (typically, we use r = 1.1, 1.2, 1.3, 1.4, 1.5).

To achieve this, the seminal work [START_REF] Borouchaki | Mesh gradation control[END_REF] proposed to enforce a gradation on the size map h. More precisely, once the size map h has been computed in such a way that:

∀ edge pq ∈ T, |h(p) -h(q)| |p -q| ≤ h grad (3.17)
This criterion is imposed on the values of h stored at the vertices of T by traveling repeatedly the edges pq ∈ T and decreasing the largest of the two values h(p) and h(q) so that (3.17) holds.

Explicit discretization of the zero level set

As mentioned, in our approach we need an explicit discretization of the zero level set of function φ: Γ φ = {x ∈ Ω : φ(x) = 0} into mesh T . This is achieved through the following marching tetrahedra procedure, which is a well-known variation of the marching cubes algorithm, see [LC87, DK91, FB96]: i. Identify the set K of elements K ∈ T intersecting Γ φ : a tetrahedron K = (P 0 P 1 P 2 P 3 ) belongs to K if and only if there exists i = j ∈ {0, 1, 2, 3} with φ(P i ) ≤ 0, and φ(P j ) ≥ 0.

ii. For an element K = (P 0 P 1 P 2 P 3 ) ∈ K, the intersection of Γ φ ∩ K is a plane portion of surface.

Identify the edges P i P j of K which intersect Γ φ (i.e. such that φ(P i ) and φ(P j ) have different signs), and compute the coordinates of the associated intersection points M ij .

iii. Travel all elements K ∈ K, and split them, introducing the pre-computed points M ij then using splitting patterns. Up to permutations, there are four possible configurations, depending on the relative signs of the φ(P i )

This procedure creates a new conforming mesh T of T but it is likely very ill-shaped since the intersections of the elements of T with Γ φ are quite arbitrary.

3.3.5

The complete adapted mesh strategy for two-phase flows in 3D

Now, starting from an initial simplicial mesh T , with associated surface triangulation Γ T , and given the four parameters h hausd , h min , h max and h grad , with h min , h max , h grad have the same meaning as in previous subsections and h hausd is the desired geometry error between the interface and its approximation in the context of Hausdoff distance. The proposed remeshing algorithm reads as follows:

i. The explicit discretisation of the interface is completed following the lines in the section 3.3.4, hence we obtain a new mesh T 1 may be still of poor quality.

ii. Construction of the size map. Although it may still be of poor quality, the new mesh T 1 accounts for a suitable discretisation of the geometry of Γ. A size map h : Ω → R dedicated to the geometric approximation of within tolerance h hausd in terms of Hausdorff distance is computed as the section 3.3.2 and stored on T 1 . This size map is then graded as detailed in section 3.3.3.

iii. Mesh modifications with respect to the size map. In this step, mesh modifications aims at producing another intermediate mesh T 2 of Ω which is a nice geometric approximation, with respect to the prescribed tolerance: d H (Γ T 2 , Γ) < h hausd . Furthermore, we rely now on lengths measured with respect to map size h by formula (3.10). Aiming at getting a new mesh whose edges have length 1, we impose that all the edges of the mesh should lie in [l min , l max ] where l min , l max are rough bounds around the target size 1.

iv. "Fine" mesh modifications with respect to the size map. Mesh T 2 should now be good in terms of geometric approximation of Γ and of better quality, and in this last stage, we perform delicately driven operations so as to get the final mesh T . Lengths of edges are still evaluated with respect to h, except we now impose T have no edge with length lying outside a sharper interval as in step iii. We are also even stricter as far as the authorized degradation in mesh quality entailed by our operators is concerned. We eventually use the vertex relocation operator to help improving the quality of the mesh.

It is clear that the final resulting T responds all the requirements of an adapted mesh in each time step for resolution of two-phase flows in section 2.5:

-The explicit discretization of the interface in T .

-The nice geometrical approximation for the interface as well as the desired error approximation of the velocity.

-The quality assessment of mesh elements.

For illustration, we show in the figure 3.7 the generation of adapted mesh in one iteration of the simulation of Rayleigh Taylor instability in 3D (detailed in section 4.

Conclusion

We have presented in this chapter some basic features of mesh adaptation. Together with chapter 2 our approach for the bifluid flows simulation with interface is fully described.

This chapter completes the theoretical presentation of the first part which contains two main problems:

Chapter 4

Numerical examples In this chapter, we present several numerical results obtained with our method. Firstly, the Navier-Stokes solver for monofluid has been valided by the Lid-driven cavity test. Next, the test of a rising bubble and an instability Rayleigh-Taylor are investigated to evaluate the ability of the scheme of two fluids problem. The results in both two dimensions (2D) and three dimensions (3D) of these simulations are given in comparison with some results in other references.

The Lid-driven cavity problem

The Lid-driven cavity problem is known as standard benchmark for Navier-Stokes solver in numerical methods and there is also a great deal of reference to compare with. The problem corresponds to the flow confined in the unit domain Ω = [0, 1] d (d = 2; 3) (a domain in three dimension is given by extending the two dimensional one in z-direction with a unit width) and the homogeneous Dirichlet boundary conditions are imposed on all boundaries: zero-velocity everywhere except on the upper one. The fluid motion is then generated by the upper lid that moves in the x -direction with a constant velocity u x = 1m/s. The viscosity is adjusted to obtain the desired Reynolds number.

Two-dimensional lid-driven cavity

In two dimensions we investigated the simulations for Reynolds number from 100 up to 10000. , in [START_REF] Gorazd | NSIKE-an incompressible Navier-Stokes solver for unstructured meshes[END_REF], and the benchmark result for cavity flow in [START_REF] Erturk | Numerical solutions of 2D steady incompressible driven cavity flow at high Reynolds numbers[END_REF] obtained with a fine uniform grid mesh of 601 × 601.

We have also obtained a good agreement of the pressure solution with the result showed in [HRK + 10], see figure 4.3. tetrahedrals for Re = 1000. In many references the results in 3D is examined by plane in 2D, so we simulate 3D problem with the same number of Reynolds in 2D: Re = 100, Re = 400, Re = 1000.

Lid-driven cavity in 3D

As expected, we obtain the streamlines in the each plane z = const are correspondent to those in 2D, see the figure 4.4 for example.

We would like to compare our numerical solutions with the results in [START_REF] Ku | A pseudospectral method for solution of the three-dimensional incompressible Navier-Stokes equations[END_REF] because the Reynolds numbers are exactly the same in our tests, but the data of velocity profiles has not been detailed there, so the comparison is on the images, see figure 4.5 for a good agreements of our velocity profiles in 2D, 3D and those obtained in given reference. x Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 Re=1000 2D 3D

Rising bubble

We consider the rising and the deformation of single bubble under gravity in fluid contained in a vertical, rectangular domain. The bubble with lower density than the surrounding fluid rises and final at the top of the domain.

Rising bubble in 2D

The initial configuration consists of a circular bubble of radius r = 0.5 centered at [2, 1.5] in a (also known as Eotvos number) defined as follows:

Re = ρ 1 √ g(2r) 3/2 µ 1 , Bo = 4ρ 1 gr 2 γ (4.1)
The problem has been set up with the constant densities and viscosities are:

ρ 1 = 100kg.m -3 , µ 1 = 0.1kg.m -1 .s -1 , ρ 2 = 1.0kg.m -3 , µ 2 = 0.01kg.m -1 .s -1 .
The gravity is: g = 9.81e -3m.s -2 . The evolution of the bubble with the coefficient of surface tension γ = 6.e -3N.m -1 and adaptive meshes with h min = 0.02, h max = 1.0 is showed in the figure 4.6. It can be seen that the bubble shape deforms during rising and the terminal bubble shape is slightly dimpled at the bottom.

In order to impress the effect of surface tension, we investigate this simulation with different coefficients of surface tension. It can be seen in the figure 4.8 when surface tension is rather small here (γ = 6e -5), the bottom of the bubble becomes more dimpled while it is flat in cases of more important coefficient (γ = 2.5e -2) and the bubble remains almost circular with the further increase in this coefficient (γ = 9.e -2). This result of bubble shapes is in good agreement with figure 6 We have known that during the long simulation, the mass of bubble can't be guaranteed. The advection of bubble front and the interpolation in the each iteration may lead to the loss of mass.

We measure the variation of mass in each time step and show in the figure 4.9. The variation of mass can be seen as about 2.E -3(0, 1%) in each time step.

Rising bubble in 3D

Extending to simulate bubble rising under gravity in 3D. We consider the problem with the same conditions in 2D: a bubble with diameter of 0.5m initialized at [0.75, 0.75, 1.0] in the domain As the simulation with low Reynolds number, we observe that the shape of bubble deforms slowly from the beginning, it becomes dimpled ellipsoidal and more distorted in time. This result is similar to those in many references with corresponding Reynolds and Bond number, see [START_REF] Lee | Single bubble rising dynamics for moderate Reynolds number using Lattice Boltzmann method[END_REF], [START_REF] Hua | Numerical simulation of 3D bubbles rising in viscous liquids using a front tracking method[END_REF] for examples. In figure 4.10 we represent the evolution of bubble from t = 0 to t = 10. It can be seen that when the bubble is very close to upper wall of the domain, its shape is rapidly distorted.

Rayleigh-Taylor instability

In this section, we carry out the simulation of more interesting problem, named Rayleigh-Taylor instability. This instability would occur along the interface of two layers where the heavy fluid is superposed the light one (ρ 2 > ρ 1 ) in the gravity field g. Starting from the results concerning the inviscid flow of [START_REF] Tryggvason | Numerical simulation of the Rayleigh-Taylor instability[END_REF], this problem has been considered by many references for viscous regime, see [START_REF] Fraigneau | Approximation of variable density incompressible flows by means of finite elements and finite volumes[END_REF] and [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF], for examples.

In our knowledge, the difficulty of the problem concerning viscous fluid depends on many criterions, among those we impress two numbers:

-The density difference represented by the Atwood number defined as:

At = ρ 2 -ρ 1 ρ 2 + ρ 1
-The Reynolds number is defined by: where d is the width of the computational domain.

Re = ρ 1 d 3 2 g 1 2 µ t = 0s t = 2.0s t = 4.0s t = 6.0s t = 8.0s t = 10.0s

Rayleigh-Taylor instability in 2D

Firstly, we set up the problem in the rectangular domain with a width of d = 1 and a height of 4d. The no-slip conditions are imposed on the upper and lower boundaries while free-slip conditions are enforced on the vertical sides. The initial interface is set as in by:

tanh y -2 -0.1cos2πx 0.01 √ 2 = 0
We obtain the results is in good symmetry during the time evolution with the Atwood number is 0.3. The results display in figure 4.12 can be compared with figure 8 of [START_REF] Lee | On the long time simulation of the Rayleigh-Taylor instability[END_REF], notice that the mushroom shapes are not identical because the given comparison showed the simulation with At = 0.5 and higher Reynolds number. We will investigate this case of the test in reduced domain in which the computational time is much cheaper and we can also verify boundary conditions for this type of domain as in [START_REF] Fraigneau | Approximation of variable density incompressible flows by means of finite elements and finite volumes[END_REF]. We can observe that the mushroom shape is more roll-up in the increasing of Atwood number. This results suggest the effect of this number on the ratio of the width of bubble and spike fluid is in good agreement with the measuring played in the figure 9 in [START_REF] Lee | On the long time simulation of the Rayleigh-Taylor instability[END_REF]. Secondly, assuming the symmetry of the initial condition is maintained during the time evolution, we consider the problem with the reduced domain with a width of d/2 and a height of 2d. We see again the same configurations with different Atwood numbers on the reduced domain in figure 4.14

As mentioned, we also validate our code by the same value of the test in [START_REF] Fraigneau | Approximation of variable density incompressible flows by means of finite elements and finite volumes[END_REF]. The evolution of the interface is plotted in figure 4.15 in time scale according to [START_REF] Tryggvason | Numerical simulation of the Rayleigh-Taylor instability[END_REF] which is related to ours by t ref = t √ d.At.g. We can see that these results are qualitatively close to those of figure 1 in [START_REF] Fraigneau | Approximation of variable density incompressible flows by means of finite elements and finite volumes[END_REF], also in [START_REF] Calgaro | An hybrid finite volume-finite element method for variable density incompressible flows[END_REF], there are in good agreement of the global characteristic of the flow in the early stage and some slight discrepancies at the large times of the evolution.

Rayleigh-Taylor instability in 3D

We consider the 3D computation of Rayleigh -Taylor instability problem. All the computational conditions are exactly the same as the 2D case and 3D instability has been well captured by the proposed scheme. The results performed in the figure 4.18 are very close to the results in figure 7 of [START_REF] Li | Numerical study for the three-dimensional Rayleigh-Taylor instability through the TVD/AC scheme and parallel computation[END_REF], or figure 18 of [START_REF] Li | A compressible Navier-Stokes flow solver with scalar transport[END_REF]. The slight difference of geometric configurations is the reason that our fluids are incompressible while the comparisons are compressible. On the other hand, the density ratio is not identical, so more detailed comparison with these results may be not interesting.

In this simulation, we begin with the mesh of 9.901 vertices (50429 tetrahedras) and the final

The coalescence of two rising bubbles

Finally, we consider the rising of two bubbles in a fluid under the gravity and the coalescence during the rising. The aim of this simulation is to show the capacity of managing topology changes of the present scheme.

The coalescence of two rising bubbles in 2D

The computational domain is

Ω = [0, 1] × [0, 4].
The under circular bubble is initially centered at (0.5, 1), with a radius of 0.25 and the upper one is initially centered at (0.5, 1.75), with a radius of 0.15. The density and the viscosity of the fluid and two bubbles are respectively:

ν 1 = 10 -1 kg/(m.s), ρ 1 = 100kg/m 3 ν 2 = 10 -2 kg/(m.s), ρ 2 = 1kg/m 3
The surface tension coefficient is set as γ = 6×10 -3 . The gravity is g = 9.81m/s 2 . In this simulation, the problem is endowed with the following boundary conditions:

-homogeneous Dirichlet conditions, u = 0, on the horizontal walls; -free-slip condition on the vertical walls: τ.σn = 0 and u.n=0;

The time step has been set to ∆t = 0.1. It is observed that the rising and the coalescence process can be divided into three consecutive stages:

-The rising of two separate bubbles.

-The coalescence of two bubbles to the one.

-The rising of aggregate bubble to the top boundary.

Initially, the shape of two bubbles are slightly distorted. We see that the shape of two bubbles become rapidly distorted and the deformations become even more pronounced until the coalescence occurs (see figure 4.20). This simulation is implemented with 200 iterations of the algorithm in

The coalescence of two rising bubbles in 3D

The computational domain is

Ω = [0, 1] × [0, 4].
The under bubble is initially centered at (0.5, 0.5, 1.0), with a radius of 0.25 and the upper one is initially centered at (0.5, 0.5, 1.6), with a radius of 0.15. The physical parameters of the fluid (ν 1 , ρ 1 ) and two bubbles (ν 2 , ρ 2 ) are similar in the 2D case:

ν 1 = 10 -1 kg/(m.s), ρ 1 = 100kg/m 3 ν 2 = 10 -2 kg/(m.s), ρ 2 = 1kg/m 3
We impose also the homogeneous Dirichlet conditions: u = 0 on the horizontal walls and the free-slip condition on the vertical walls: τ.σn = 0 and u.n=0. This simulation is implemented with 200 iterations. The adapted mesh is about 30000 vertices at each time with the parameters: h min = 0.005, h max = 0.2, h grad = 1.1, h hausd = 0.005. See figure 4.27 for more detail images of initial mesh and adapted one at the coalescence. The computational time takes about 2.5 minutes for each iteration.

Conclusion

The numerical results obtained in this chapter are in good accordance with those in many references which have been cited in each test case. It can be seen that the proposed scheme for bi-fluid flows, that is based on mesh adaptation, is capable of managing complex movements of the interface include topology change. Furthermore, the computational adapted triangulations have a modest number of nodes, while maintaining a relatively high accuracy near the interface, and the computational time is acceptable.

This chapter finishes also the first part of this thesis. We go to the second part where the considered problem is related to the optimization of computational domain, i.e the computational domain is variable during the time evolution. Part II

Shape optimization in fluid mechanics

In the first part of this thesis, we have mentioned the following problems:

-The numerical resolution of Navier-Stokes equations for two fluids flow with the interface.

-The numerical resolution of the evolution of the interface (free surface) by a formulation of the level set method.

-Mesh adaptation techniques, especially for controlling the geometrical discrepancy between the free surface and its approximation.

All these problems and their numerical counterparts allow us to resolve any stationary or non stationary problems for two phase flows in a fixed domain. However, in many aspects of the application another challenge is related to the optimization of the geometry of the computational domain with respect to a given functional. This part of computational mechanics has soared over the last decades and is devoted to finding the optimal size, shape or topology of a mechanical part, with respect to a given mechanical criterion -such as, for instance, the work of external loads, or the internal stress.

These problems strongly depend on the physics at play (elasticity, thermoelasticity, fluid mechanics ...), and on the constraints that should be fulfilled by the shapes.

In the context of structural optimization, size, shape and topology optimization problems address different aspects of the structural design problem. In a typical size problem the goal may be to find the optimal thickness distribution of a linearly elastic plate or the optimal member areas in a truss structure. On the other hand, in a shape optimization problem the goal is to search the optimum shape of this domain, that is, the shape is now the design variable. Unlike the case of classical -Last but not least, shape optimization in fluid mechanics has received a large amount of attention from both engineers and mathematicians since the pioneering works of Pironneau [START_REF] Pironneau | On optimum profiles in Stokes flow[END_REF][START_REF] Pironneau | On optimum design in fluid mechanics[END_REF]. Many applications of shape optimization concepts in this domain, to name a few: the optimal shape design (OSD) of airplane wing which induces minimal drag, i.e. minimal reaction from the surrounding fluid in aeronautic industry, see [START_REF] Mohammadi | Shape optimization in fluid mechanics[END_REF][START_REF] Painchaud-Oullet | Airfoil shape optimization using a nonuniform rational b-spline parametrization under thickness constraint[END_REF]; the optimal control approaches to shape optimization of aorto-coronaric bypass anastomoses are applied in [START_REF] Quarteroni | Optimal control and shape optimization of aortocoronaric bypass anastomoses[END_REF][START_REF] Agoshkov | A mathematical approach in the design of arterial bypass using unsteady Stokes equations[END_REF][START_REF] Agoshkov | Shape design in aorto-coronaric bypass anastomoses using perturbation theory[END_REF]; the problem of optimal swimming of microorganisms at low Reynolds number is considered in [START_REF] Alouges | Numerical strategies for stroke optimization of axisymmetric microswimmers[END_REF], more specifically, the authors are in search for the shape of such microorganism that allows to reach a prescribed displacement while undergoing minimal stress from the surrounding fluid, the shape of the microorganisms is parametrized by means of a few physical parameters, and the sensitivity of the stress exerted by the fluid on the microorganisms with respect to perturbations on these parameters is computed; an interesting application of the topology optimization method for fluid flow can be found in [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF] for minimum power dissipation Stokes flow problems in two dimensions. Since then, alternative parameterizations have been suggested for the Stokes flow problem [START_REF] Guest | Topology optimization of creeping fluid flows using a Darcy-Stokes finite element[END_REF] and the Navier Stokes flow problem [START_REF] Evgrafov | Topology optimization of slightly compressible fluids[END_REF]. Topology optimization for Stokes flow has been extended to large-scale problems [START_REF] Aage | Topology optimization of large scale Stokes flow problems[END_REF] and has been applied to fluid flows with regions of Darcy and

Stokes flow [START_REF] Wiker | Topology optimization of regions of Darcy and Stokes flow[END_REF] and with low-to-moderate Reynolds numbers [START_REF] Olesen | A high-level programming-language implementation of topology optimization applied to steady-state Navier-Stokes flow[END_REF][START_REF] Duan | Shape-topology optimization for Navier-Stokes problem using variational level set method[END_REF][START_REF] Zhou | A variational level set method for the topology optimization of steady-state Navier-Stokes flow[END_REF].

For more works and many other applications of shape or topology optimization concepts in fluid mechanics, we refer to [BS03, Gun03, MP04, MP10] and references therein.

In this part our aim is to construct a numerical scheme for the problem of shape optimization in

     -ν∆u + ∇p = f in Ω divu = 0 in Ω (1)
with some adequate boundary conditions. The problem we consider is to minimize an objective functional denoted by J(Ω) which depends on the domain Ω via the solution u Ω of (1) where the variable shape Ω belongs to a set of admissible shapes U ad . Formally speaking, our model of shape optimization is:

inf

Ω∈U ad J(Ω)
This problem has also been considered in [START_REF] Borrvall | Topology optimization of fluids in Stokes flow[END_REF][START_REF] Duan | Shape-topology optimization of Stokes flow via variational level set method[END_REF][START_REF] Duan | Shape-topology optimization for Navier-Stokes problem using variational level set method[END_REF]. Before talking about our approach, let us mention some techniques of shape optimization that have been introduced in computational mechanics since the 1960s. Regarding numerous methods to handle the structural optimization problems in the literatures, one can see that these methods are mainly distinguished by the way they represent the shapes and the way they compute the sensitivity of the objective criterion with respect to the design.

Describing shapes is a complicated problem since it needs to satisfy two requirements. Firstly, shape optimization techniques must be able to perform mechanical computations on the considered shapes, e.g. by means of finite differences, finite element or finite volume methods, and not all kind of representations lend themselves to such computations. Secondly, the adopted representation must be versatile enough to allow for a robust account of shapes' deformations. Some typical approaches for descriptions of shapes can be named as:

-explicit methods: in this case, shapes are explicitly accounted by a set of degrees of freedom, see for instance [START_REF] Braibant | Shape optimal design using B-splines[END_REF] around these issues, or the review article [START_REF] Haftka | Structural shape optimization -a survey[END_REF], and monograph [START_REF] Pironneau | Optimal Shape Design for Elliptic Systems[END_REF].

-homogenization methods and their variants (power-law method or solid isotropic microstructure with penalization (SIMP) method), may be considered quite classical in view of the number of publications see [START_REF] Bendsøe | Generating optimal topologies in structural design using a homogenization method[END_REF][START_REF] Bendsøe | Methods for optimization of structural topology[END_REF][START_REF] Bendsøe | Topology Optimization, Theory, Methods and Applications[END_REF][START_REF] Allaire | Shape optimization by the homogenization method[END_REF][START_REF] Allaire | Shape optimization by the homogenization method[END_REF] for examples. These methods belong to the class of density methods. They have made a notable change in perspectives in structural optimization by using a density function θ : D → [0, 1] defined over a computational domain D, where θ is close to 0 if there is almost only void (or a very "soft" material, mimicking void), and where θ is close to 1 if there is almost only the shape. The problem of finding the "best"

shape is transformed into that of the optimal distribution of a mixture of material and void in a large computational domain. Therefore, the shape optimization problem ends rephrased as a parametric optimization problem.

-implicit methods: have been used broadly for interface-tracking or interface-capturing in various fields and this method became very popular in shape-topology optimization, the most proposed new evolution equations for the level set function in order to achieve a smooth evolution without reinitialization. Zhou and Li [START_REF] Zhou | A variational level set method for the topology optimization of steady-state Navier-Stokes flow[END_REF] have also used the level set method for the topology optimization of steady-state Navier Stokes flows in both two and three dimensions.

They stated that the computational costs incurred by remeshing the fluid domain, extending the velocity and reinitializing the level set function were limitations of their model.

It is obvious that the performing of perturbations of the considered shapes is closely related to the problem of shape description. Our approach for shape description belongs to the last class:

we apply a conventional level set approach [WWG03, AJT04] where a precise description of the boundary is kept under the form of an auxiliary function (level set function). In our approach the deformation of shape Ω is reformulated in terms of an associated level set function φ(., t) as the following level set advection equation over a larger domain D :

∂φ ∂t + V.∇φ = 0
where the advection vector V is the steepest descent direction for optimization corresponding to the shape derivative of the given objective function. We also compute a shape derivative by using an adjoint problem, however unlike the approach in [AJT04] using the so-called "ersatz material" which amounts to fill the holes by a weak phase, we only need impose the Stokes and adjoint problems on a sub-mesh T of the shape domain. This can be achieved thanks to the explicit discretization of generated shape during the optimization process.

It is worth noting that although this method is not specifically designed for topology optimization, it can easily handle topology changes. For this reason, if the optimum design of the structure includes information on the topology and shape of the structure, the level set models allow to deal with these two problems simultaneously. We shall describe the details of our proposed approach and its numerical implementation in the chapter 5 hereafter.

Chapter 5

Shape optimization in fluid mechanics As mentioned in the introduction, our numerical scheme is based on a combination of classical shape derivative, level set method and mesh adaptation. The motivation of this approach related to the following components:

-In the context of level set method, the variational shape Ω is classically represented as the negative subdomain of a level set function φ of "larger" computational domain D. In the first part of this thesis, this method has been approached for the evolution of the interface Γ where 115 Γ is represented as the zero-level set of φ. Hence, if we consider the boundary ∂Ω of variable shape Ω as the "interface" defined by the zero-level set: ∂Ω = {x ∈ D : φ(x) = 0}, then the problem of shape deformation will be considered as the problem of domain evolution which involves only to the advection of ∂Ω. In this way, we can apply the numerical resolution of advection equation (is detailed in Chapter 2) for the evolution of boundary of Ω, and hence for the deformation of Ω.

-Consequently, during the shape deformation the mesh adaptation process is constructed by the same approach (Chapter 3) to guarantee the accuracy of geometrical approximation of the continuous shape Ω as well as the quality of the adapted mesh for the resolution of mechanical problems.

-Level set method has been used broadly in structural optimization. In [START_REF] Dapogny | Shape optimization, level set methods on unstructured meshes and mesh evolution[END_REF], the combination of level set methods on unstructured meshes and mesh evolution has led to successful results in shape optimization of elastic structures.

The main difficulty is related to the problem of introducing perturbations of the considered shapes, and thus on how to compute the sensitivity of the objective function J with respect to the shape. We overcome this problem in section 5.1, and show that the computation of the sensitivity of objective function is related to the Hadamard's boundary variation method and the shape derivatives is computed by Céa's formal method. They are the key ingredients to construct a general scheme for shape optimization in fluid mechanics (section 5.2). Section 5.3 will detail our numerical scheme using level set method and mesh adaptation. Last, in section 5.4 a numerical example with the objective function of energy dissipation is presented to assess the efficiency and the reliability of the proposed scheme.

Shape sensitivity analysis using Hadamard's boundary variation method

In this section we focus on one particular method for describing variations of a shape, namely

Hadamard's boundary variation method, as well as the notions of differentiation with respect to the domain. At first, we introduce the basic ideas of Hadamard's method and set some notations, that we shall use in the rest of this part. Then, we present the derived notions of differentiation with respect to a domain: the notion of shape derivative of a scalar function of the domain, or of a function which is itself defined on the domain is introduced. Some classical results is recalled to

SHAPE SENSITIVITY ANALYSIS USING HADAMARD'S BOUNDARY VARIATION METHOD117

prepare for the computation of shape derivative. Finally, we apply the stress on a specific context: the optimization of mechanics fluid in context of Stokes system.

Hadamard's boundary variation method

The central idea of Hadamard's boundary variation is proposed in its seminal paper [START_REF] Hadamard | Mémoire sur le problème d'analyse relatif à l'équilibre des plaque élastiques encastrées[END_REF] (see also [START_REF] Sokolowski | Introduction to Shape Optimization; Shape Sensitivity Analysis[END_REF]). It was then exploited in depth in [START_REF] Murat | Etudes de problèmes d'optimal design[END_REF]. Let θ be a displacement field of "small" amplitude, we consider variations of a given reference shape Ω of the form (I + θ)(Ω) (see figure 5.1).

Thus, all the considered transformations (I + θ)(Ω) are homeomorphisms "close" to the identity; in particular, all the variations of Ω achieved in this way share the same topology.

Note that in the framework of Hadamard's method, the variations Ω θ of Ω only depend on the values taken by θ on ∂Ω, and we have following consequence of Picard's fixed point theorem (see lemma 6.13 in [START_REF] Allaire | Conception optimale de structures[END_REF] for a proof):

Lemma 1. For every deformation field θ ∈ W 1,∞ (R d , R d ) such that θ W 1,∞ (R d ,R d ) < 1, the appli- cation (I + θ) : R d → R d is a Lipschitz homeomorphism with Lispchitz inverse.
To the end of this part Ω ⊂ R d stands for a fixed domain, which enjoys at least Lispchitz regularity. For any θ

∈ W 1,∞ (R d , R d ), θ W 1,∞ (R d ,R d ) < 1, we denote as Ω θ := (I + θ)(Ω) the deformed shape with respect to θ. Note that W 1,∞ (R d , R d ) ⊂ L ∞ (R d ) d
is the Banach space of bounded functions θ : R d → R d , endowed with the natural norm:

∀θ ∈ W 1,∞ (R d , R d ), θ W 1,∞ (R d ,R d ) := θ L ∞ (R d ) d + ∇θ L ∞ (R d ) d×d
Hence, variations of a given shape Ω end up parametrized by means of an open subset of a Banach space. As we shall see in the next section, this allows among other things to introduce a notion a differentiability with respect to the shape by rewriting any operation performed on shapes close to Ω in terms of the underlying deformation field θ. Definition 4. Let J(Ω) a functional of the domain. J is shape differentiable at Ω if the underlying application:

Shape differentiability

W 1,∞ (R d , R d ) -→ R θ -→ J(Ω θ )
is Fréchet differentiable at θ = 0. The associated Fréchet differential; denoted as J (Ω) is called the shape derivative of J at Ω. We have following expansion in the vicinity of 0

∈ W 1,∞ (R d , R d ): J(Ω θ ) = J(Ω) + J (Ω)(θ) + o(θ), with lim θ→0 |o(θ)| θ = 0 (5.1) Remark 6. 1. J (Ω) is a continuous linear form on W 1,∞ (R d , R d ).
2. In the case of a functional which also depends on other variables than Ω, the partial Fréchet differential with respect to the domain is denoted as ∂ ∂Ω .

We recall a classical result which states that the directional derivative J (Ω)(θ) depends only the normal trace θ.n on the boundary ∂Ω, see [START_REF] Allaire | Structural optimization using shape sensitivity analysis and a level-set method[END_REF].

Lemma 2. Let Ω be a smooth bounded open set and J(Ω) a differentiable function at

Ω, if θ 1 , θ 2 ∈ W 1,∞ (R d , R d ) are such that: θ 2 -θ 1 ∈ C(R d , R d ) and θ 1 .n = θ 2 .
n, then we have:

J (Ω)(θ 1 ) = J (Ω)(θ 2 )
We also recall here two results on shape derivative that will be useful in the computation of shape derivatives in the section 5.1.3 (see [START_REF] Henrot | Variation et optimisation de formes, une analyse géométrique[END_REF] or [START_REF] Allaire | Structural optimization using shape sensitivity analysis and a level-set method[END_REF] for the proofs).

SHAPE SENSITIVITY ANALYSIS USING HADAMARD'S BOUNDARY VARIATION METHOD119

Lemma 3. Let Ω be a bounded Lispchitz domain and let f (x) ∈ W 1,1 (R d ). Define: J(Ω) = Ω f (x)dx then J(Ω) a differentiable function at Ω, and we have:

J (Ω)(θ) = Ω div(θf )dx = ∂Ω f θ.nds for any θ ∈ W 1,∞ (R d , R d ) Lemma 4.
Let Ω be a bounded Lispchitz domain in class C 2 and g(x) ∈ W 2,1 (R d ). Define: J(Ω) = ∂Ω g(x)ds then J(Ω) a differentiable function at Ω, and:

J (Ω)(θ) = ∂Ω ( ∂g ∂n + κg)θ.nds for any θ ∈ W 1,∞ (R d , R d ),
where κ is the mean curvature of ∂Ω defined by κ = divn. Furthermore, this result still holds true if one replaces ∂Ω by Γ, a smooth open subset of ∂Ω, and assumes that h = 0 on the surface boundary of Γ.

In particular, Lemma 3 is useful in order to compute the shape derivative of a volume constraint V (Ω) = C. Indeed, we have:

V (Ω) = Ω dx and V (Ω)(θ) = ∂Ω θ(x).n(x)ds. (5.2)
Similarly, Lemma 4 is useful in order to compute the shape derivative of a perimeter constraint P (Ω) = C since we have:

P (Ω) = ∂ Ω
ds and P (Ω)(θ) = ∂Ω κθ(x).n(x)ds.

5.1.3 Lagrange's approach for the computation of shape derivatives in the context of Stokes system

In this section, we present the computation of the shape derivative for objective functionals J(Ω) which depend on the domain Ω via the solutions u Ω to the Stokes system:

     -∆u + ∇p = f in Ω divu = 0 in Ω (5.3)
endowed with some adequate boundary conditions.

The rigorous presentation of this topic is not an easy task. In practice, we rely on a formal method to obtain the expressions of these shape derivatives, namely Céa's fast derivation method, introduced in [START_REF] Céa | Conception optimale ou identification de formes, calcul rapide de la dérivée directionnelle de la fonction coût[END_REF]. It is formal in the sense that all the data at hand: objective functions, domains, deformations are assumed smooth enough and the state solution u Ω , p Ω is differentiable with respect to the shape. Consider the following cases of boundary conditions:

Newmann boundary condition

We consider the objective function of type:

J(Ω) = Ω j(u Ω )
and the state equations (5.3) are endowed by Newmann condition as follows:

             -∆u + ∇p = f in Ω divu = 0 in Ω σ(u, p)n = ϕ N on ∂Ω (5.4)
In this case, we introduce the Lagrange functional defined for (u, p, v, q) as follows:

L(Ω, u, p, v, q) = Ω j(u)dx + Ω (-∆u + ∇p -f ).vdx - Ω qdivu Taking v ∈ (H 1 (R d ) d
) and q ∈ L 2 (R d ) and using Green formula:

L(Ω, u, p, v, q) = Ω j(u)dx + Ω (∇u.∇v -pdivv -f.v -qdivu)dx - ∂Ω ϕ N .v
We have the partial differential of L with respect to the each variable as follows:

∂L ∂v (Ω, u, p, v, q), ū = Ω (∇u.∇ū -pdivū -f.ū)dx - ∂Ω ϕ N .ū (5.5)
∂L ∂q (Ω, u, p, v, q), p = Ω -pdivu (5.6)

∂L ∂u (Ω, u, p, v, q), v = Ω j (u).vdx + Ω (∇v.∇v -qdivv)dx (5.7) ∂L ∂p (Ω, u, p, v, q), q = Ω -qdivv (5.8)
We can see that (5.5), (5.6) lead the variational formulation of the state system (5.4). It has been known that in the weak form this problem has unique solution (u Ω , p Ω ), and we can easily verify that:

J(Ω) = L(Ω, u Ω , p Ω , v, q), ∀ (v, q) ∈ (H 1 (R d ) d , L 2 (R d ))
Taking the partial differential with respect to domain Ω we have:

J (Ω)(θ) = ∂L ∂Ω (Ω, u Ω , p Ω , v, q)(θ) + ∂L ∂u (Ω, u Ω , p Ω , v, q), u (Ω)(θ) + ∂L ∂p (Ω, u Ω , p Ω , v, q), p (Ω)(θ) (5.9)
Hence, suppose that (v Ω , q Ω ) is solution of the system:

     ∂L ∂u (Ω, u Ω , p Ω , v, q) = 0 ∂L ∂p (Ω, u Ω , p Ω , v, q) = 0
(5.10)

SHAPE SENSITIVITY ANALYSIS USING HADAMARD'S BOUNDARY VARIATION METHOD121

then J (Ω)(θ) can be written with the reduced formula:

J (Ω)(θ) = ∂L ∂Ω (Ω, u Ω , p Ω , v Ω , q Ω )(θ) (5.11)
The transformations (5.7), (5.8) lead us to formulate the system (5.10), also called the adjoint equations as:

             -∆v + ∇q = -j (u Ω ) in Ω divv = 0 in Ω σ(v, q)n = 0 on ∂Ω (5.12)
Applying Lemma 3, Lemma 4 and taking into account the expression (5.11), the shape derivative of objective function can be computed as follows:

J (Ω)(θ) = ∂Ω θ.n j(u) + ∇u.∇v -pdivv -f.v -qdivu - ∂(ϕ N .v) ∂n -κϕ N .v ds (5.13)

Dirichlet boundary condition

Let us now consider the case of Dirichlet condition, i.e the state equations are:

             -∆u + ∇p = f in Ω divu = 0 in Ω u = u D on ∂Ω (5.14)
We construct the Lagrange function as:

L(Ω, u, p, v, q, λ) = Ω j(u)dx+ Ω (-∆u+∇p-f ).vdx- Ω qdivudx+ ∂Ω λ.(u-u D )ds (5.15)
The partial differential of L respect to the each variable:

∂L ∂v (Ω, u, p, v, q, λ), ū = Ω (-∆u + ∇p -f ).ūdx (5.16) ∂L ∂q (Ω, u, p, v, q, λ), p = Ω -pdivu (5.17) ∂L ∂λ (Ω, u, p, v, q, λ), λ = ∂Ω λ.(u -u D ) (5.18) ∂L ∂u (Ω, u, p, v, q, λ), v = Ω j (u).vdx + Ω -∆v.vdx - Ω qdivvdx + ∂Ω λ.vds = Ω j (u).vdx + Ω -∆v.vdx + Ω ∇q.vdx + ∂Ω (∂ n v -qn + λ).vds - ∂Ω v.∂ n vds = Ω (j (u) -∆v + ∇q).vdx + ∂Ω (σ(v, q)n + λ) .vds - ∂Ω v.∂ n vds (5.19) ∂L ∂p (Ω, u, p, v, q, λ), q = Ω ∇q.v = - Ω qdivv + ∂Ω qn.v (5.20)
Introducing the condition σ(v, q)n + λ = 0 on ∂Ω, the transformation (5.19) and(5.20) lead us to the adjoint system: It has been observed that the equilibrium of ∂L ∂v (Ω, u, p, v, q) and ∂L ∂q (Ω, u, p, v, q) yields the variational formulation associated to the state system (5.24). There exists unique solution (u Ω , p Ω ) of the weak form of the state system and we can easily verify that:

             -∆v + ∇q = -j (u) in Ω divv = 0 in Ω v =
J(Ω) = L(Ω, u Ω , p Ω , v, q) ∀ (v, q) ∈ (H 1 Γ D (R d ) d , L 2 (R d ))
where w Ω is a certain scalar defined over ∂Ω (depending on the Stokes and the adjoint systems), which makes the computation of the possible deformation field for J(Ω): θ t = -tw Ω n (5.35) Indeed, from the expansion (5.1), for t > 0 small enough, we have:

J(Ω tθ ) = J(Ω) -t ∂Ω w 2 Ω ds + o(t) < J(Ω),
hence, the new shape Ω tθ is obtained from initial form Ω as: Ω tθ = (I + θ t )(Ω).

The global algorithm

We are now in position to introduce a shape gradient algorithm for the considered problem (5.3).

Algorithm 5: Shape gradient algorithm n=0, start with an initial guess Ω 0 ;

for n = 1, . . . , untill convergence do 1. Compute the solution (u Ω n , p Ω n ) to the Stokes system on Ω n , and the solution (v Ω n , q Ω n ) to the adjoint system (section 5.1.3).

2. Infer a descent direction θ n for the objective functional following the lines in section 5.2.1.

3. Choose a descent step τ n > 0 small enough so that J(Ω n τ n θ n ) < J(Ω n ). 4. The new shape is obtained as Ω n+1 = Ω n τ n θ n end 5.3 Shape optimization basing on level set method and mesh adaptation

Description of the level set method for shape optimization

The idea of this method is to combine an implicit domain evolution method with an explicit type of shape representation which have been used in many references: for examples, in the twodimensional work [START_REF] Yamasaki | A level set-based topology optimization method targeting metallic waveguide design problems[END_REF], the evolution of shapes is tracked on a triangular mesh D of a working domain D owing to the level set method, and at each iteration of the process, an exact mesh of the current shape Ω is obtained by relocating vertices of D onto ∂Ω. In [XSLW12], a similar strategy is applied for dealing with the motion of shapes; a computational mesh for any shape arising during the process is moreover constructed by first identifying the intersection points of the implicitly-defined boundary ∂Ω with the edges of the mesh D of D, then using them as an input for a Delaunay-based mesh generation algorithm. In our approach, a computational domain D is defined, and is equipped with an unstructured mesh which is modified at each iteration of the algorithm, in such a way that all the shapes Ω 0 , . . . , Ω n produced during the shape optimization process are explicitly discretized in this mesh, i.e all elements (vertices, edges, faces) of mesh of shape are also as elements of the mesh of D (see figure 5.2, left). The direction of the deformation of each shape Ω n is inferred from computing the shape gradient J (Ω n ) (which depends on solutions of the state and the adjoint systems on Ω n , see section 5.1.3). The determination and the deformation of the shape Ω n is related to the level set function φ n defined on D: at each step, φ n is advected by descent direction of Ω n to obtain the new level set function φ n+1 . This function leads to define the new shape Ω n+1 as the negative subdomain of D (see figure 5.2, right). This strategy presents several advantages: problem only on a submesh (the negative part). hence the computational time can be reduced, since we know that normally we need to resolve both state and adjoint problem which always are challenges in numerical resolution.

-As a shape is known as a submesh of the computational mesh of D, no projection between different meshes is involved in the computation of a descent direction for the considered objective function of the domain.

-The proposed method does not pose any theoretical difficulty to be extended to the threedimensional case. so that the set Θ ad of admissible variations is:

Θ ad = {θ ∈ W 1,∞ (R d , R d ) : θ = 0 on Γ D ∪ Γ N }.
Given a objective function J(Ω), in fact the minimization problem: where U ad is defined by (5.38) and l > 0 is a positive Lagrange multiplier.

Other regularized variants of (5.39) have been proposed in [START_REF] Chenais | On the existence of a solution in a domain identification problem[END_REF][START_REF] Chambolle | A density result in two-dimensional linearized elasticity and applications[END_REF] for such existence theories or more efficient algorithms for handling constraints such as the Method of Moving Asymptotes described in [START_REF] Svanberg | The method of moving asymptotes -a new method for structural optimization[END_REF], or the Method of Feasible Directions in [START_REF] Vanderplaats | Structural optimization by methods of feasible directions[END_REF], but they are not dealt in the scope of this thesis. This implies that for the optimization problem under a volume constraint, the computation of the scalar function w Ω in (5.35) requires adding the Lagrange multiplier l. In order to simplify the notation, the notation J(Ω) of the objective function will hereafter have the meaning of (J(Ω) + lV (Ω)) while w Ω stands for (w Ω + l).

Computation of a descent direction

From a given shape Ω ∈ U ad , a descent direction V ∈ Θ ad for the considered objective function J(Ω) is computed on a whole mesh D Ω of D which encloses an explicit discretization of Ω.

The generic expression (5.35) for the shape derivative of J suggests the immediate choice:

∀x ∈ ∂Ω, V (x) = -w Ω (x)n(x).

(5.41)

As we have seen, w Ω depends on the solution of Stokes and adjoint systems posed on Ω, which can be accurately solved on the submesh T Ω of D Ω , using the finite element method. Unfortunately, the choice (5.35) for a descent direction turns out to be hazardous for two independent reasons:

SHAPE OPTIMIZATION BASING ON LEVEL SET METHOD AND MESH ADAPTATION129

-Formula (5.35) is only defined on ∂Ω, whereas we are in search for a velocity field V which is defined at least in a vicinity of ∂Ω, as well to comply with the requirement V ∈ Θ ad as to use

V in the context of the level set method.

-The scalar field w Ω usually depends on the derivatives of the Stokes solution and the adjoint solution, which may be quite irregular -in the theoretical framework as well as when it comes to their numerical approximation. This may jeopardize the numerical stability of the process.

As advocated by [START_REF] Burger | A framework for the construction of level-set methods for shape optimization and reconstruction[END_REF][START_REF] De Gournay | Velocity extension for the level-set method and multiple eigenvalues in shape optimization[END_REF] (see also in chapter 2) an efficient way to address both problems at the same time consists in using the gradient of J as a descent direction associated with a different scalar product instead of canonical one of L 2 (Γ).

More accurately, α > 0 being a small "extension -regularization" parameter, let us consider the functional space:

H 1 Γ D ∪Γ N (D) = {w ∈ H 1 (D), w = 0 on Γ D ∪ Γ N }
and let w ∈ H 1 Γ D ∪Γ N be the unique solution to the variational problem (see [START_REF] Burger | A framework for the construction of level-set methods for shape optimization and reconstruction[END_REF] for alternative choices): where n is an extension to D of the normal vector field to ∂Ω. Using the asymptotic expansion (5.1) we can show that V is still a descent direction for J (for t small enough: J(Ω t V ) = J(Ω)t Γ w 2 w Ω nds + o(t)). However, V intrinsically enjoys more regularity than V owing to the classical regularity theory for elliptic equations, and is inherently defined on the whole domain D.

In the numerical setting, w is easily computed by solving (5.42) with the classical finite element method carried out on mesh D Ω , after computing w Ω (the discretization of the right-hand side being straightforward since the computational mesh D Ω encloses an explicit discretization of ∂Ω).

The proposed algorithm

We are now in position to outline the proposed general strategy for handling mesh evolution in the context of shape optimization.

Starting with an initial shape Ω 0 , and a simplicial mesh D Ω 0 of D in which Ω 0 is explicitly discretized.

For n = 0, ... untill convergence, the current shape Ω n is known via a submesh T Ω n of D Ω n where T Ω n is a mesh of Ω n .

1. Generate the level set function φ n (as the signed distance function) to Ω n on the whole mesh

D Ω n of D.
2. Compute the value of the scalar field w Ω n appearing in the shape derivative of the considered functional in (5.35). This may involve one, or several finite element analyses for solving the Stokes and adjoint systems, to be held on the part T Ω n of the mesh D Ω n corresponding to Ω n .

The quantity w Ω n is defined only on ∂Ω n , i.e. in the numerically setting, on the discretization of ∂Ω n which explicitly appears in both T Ω n and D Ω n .

3. Extend w Ω n to a vector field V n defined on the whole mesh D Ω n of D, following the lines in section 5. 6. Evaluate J(Ω n+1 ). If J(Ω n+1 ) < J(Ω n ), Ω n+1 is retained as the new shape; else Ω n+1 = Ω n .

Then, go back to step 4, decreasing the chosen value for the time step.

Numerical examples

We consider the situation depicted in figure 5.3: a flow with the constant viscosity ν travels throughout a domain of elbow Ω with three parts of the boundary:

∂Ω = Γ ∪ Γ D ∪ Γ N
We endow a parabolic velocity profile u D on the inlet boundary Γ D and assume that the flow comes out at the boundary Γ N with the free-out condition. This domain Ω is the optimization variable of the problem where Γ D , Γ N are the subset of fixed part of the boundary and only the part Γ is variable during the optimization process. We are interested in the minimization the energy dissipation in fluid, then the object functional J(Ω) defined as: w Ω = ν(∇u.∇u -∇u.∇v)

J(Ω) =
where u is the solution of Stokes systems (5.37) and v is solution of adjoint systems constructed by (5.31) w Ω is extended to obtain the descent direction V defined on the whole mesh D. As explained in the section 5.3.2, we will consider this minimization problem with a volume constraint incorporated under the form of a penalization by a fixed Lagrange multiplier l. In this case test, we set l = 2.5

and 300 iterations of the algorithm of section 5. We have reduced considerably the energy dissipation:

J(Ω f inal ) -J(Ω 0 ) J(Ω 0 ) × 100 40%

It can be observed that in this case test, the final shape is the domain connected by the "almost straight" lines between input and output parts. This result is in good agreement with the one obtained in [BP03, DMZ08b, Pri08, CG09].

Conclusion

This chapter has given a numerical scheme for shape-topology optimization in fluid mechanics. -Extend the proposed scheme in the context of Navier-Stokes problem like in [START_REF] Duan | Shape-topology optimization for Navier-Stokes problem using variational level set method[END_REF] and may be combine with the different objective functions as the suggest of maximum permeability in [START_REF] Zhou | A variational level set method for the topology optimization of steady-state Navier-Stokes flow[END_REF]. -Our proposed scheme which is based on a combination of classical shape derivative, level set method and mesh adaptation process has been utilized successfully in the elastic structure [START_REF] Dapogny | Shape optimization, level set methods on unstructured meshes and mesh evolution[END_REF]. This study has demonstrated its efficiency and its reliability in fluid mechanics problem. We could expect a promising application of the present scheme for shape-topology optimization in other physical problems (thermal structure,... ).

IntroductionI

  Numerical simulation of the problem two-fluid flows Introduction 1 Numerical resolution of Navier-Stokes equations 1.1 The Navier-Stokes equations for incompressible fluid . . . . . . . . . . . . . . . . . . 1.2 The method of characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 The Lagrange-Galerkin formulation . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 The discretization in time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 The variational formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Homogeneous Dirichlet boundary conditions . . . . . . . . . . . . . . . . . . . 1.3.2 Case of general boundary conditions . . . . . . . . . . . . . . . . . . . . . . . 1.4 The spatial discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 The discrete linear systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5.1 Construction of the matrix formulation . . . . . . . . . . . . . . . . . . . . . . 1.5.2 Matrix assembly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5.3 Solving the linear systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Numerical examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Two-phase flow simulations 2.1 Introduction the model of two-phase flow . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 The equation of motion of bifluid flow . . . . . . . . . . . . . . . . . . . . . . 3

7 CONTENTSFigure 1 :

 71 Figure 1: Extreme changes of the interface: folding (top, left), breaking (bottom, left) and merging (right).

Figure 4 :

 4 Figure 4: (Left) A shape Ω ⊂ D, (Center) representation of shape Ω by associated level set function and (Right) triangular mesh of D enclosing a mesh of Ω (yellow elements).

  Navier-Stokes equations for incompressible fluid . . . . . . . . . . . 24 1.2 The method of characteristics . . . . . . . . . . . . . . . . . . . . . . . . . 25 1.2.1 The Lagrange-Galerkin formulation . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 The discretization in time . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 The variational formulation . . . . . . . . . . . . . . . . . . . . . . . . . . 29 1.3.1 Homogeneous Dirichlet boundary conditions . . . . . . . . . . . . . . . . . . 1.3.2 Case of general boundary conditions . . . . . . . . . . . . . . . . . . . . . . 1.4 The spatial discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 1.5 The discrete linear systems . . . . . . . . . . . . . . . . . . . . . . . . . . 34 1.5.1 Construction of the matrix formulation . . . . . . . . . . . . . . . . . . . . 1.5.2 Matrix assembly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5.3 Solving the linear systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Numerical examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 1.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

  figure 1.1 for illustrations the travel of given point x in domain Ω with substep time δt ∆t.

Figure 1 . 1 :

 11 Figure 1.1: (left) travel characteristic point with time δt and (right) characteristic point is go out to domain.

Algorithm 1 :

 1 Numerical scheme for solving Navier-Stokes equation over [0, T ] n=0, initializations (u 0 , p 0 ); Given time step ∆t; nmax = [T /∆t];

Figure 1 . 2 :

 12 Figure 1.2: Degree of freedom: (a) mini element (P 1 bubble/P 1 ); (b)Taylor-Hood element (P2/P1)

  B, G are initialized by zero and add them up to the appropriate entries by finding index global corresponding to index local of [Ae] K , [Be] K , [Ge] K . See algorithm 2 for example of the assembly matrix A.

. 4 ,Figure 1 . 3 :Figure 1 . 4 :

 41314 Figure 1.3: Cavity in 2D: velocity profile for u x and u y in case of Re = 1000 with the different time steps.

Figure 1 . 5 :

 15 Figure 1.5: Cavity in 3D: the streamtraces for the test case Re = 1000.

Figure 1 . 6 :

 16 Figure 1.6: Cavity in 3D: (left) velocity of the flow and (right) streamlines in the plane (z = 0.5) for Re = 1000.

Figure 2 . 1 :

 21 Figure 2.1: Two examples of computational domains for bifluid flow simulations; left: disconnected components, right: connected components.

Figure 2 . 2 :

 22 Figure 2.2: Examples of computational domains with level-set function in 2D.

Figure 2 . 3 :

 23 Figure 2.3: (left) Cutting plane through a computational domains with level-set function and (right) several iso-surfaces in 3D.

  t = 0s t = 2.5s t = 5.0 t = 7.5s t = 10.0s

Figure 2 . 4 :

 24 Figure 2.4: Rising bubble in 2D: evolution of the interface in time: (top) without surface tension;

Figure 2 . 5 :

 25 Figure 2.5: Ilustration of geometric discretization of interface in 2D.

Figure 2 . 6 :

 26 Figure 2.6: Ilustration of characteristic curve approximation.

Figure 2 . 7 :Figure 2 . 8 :

 2728 Figure 2.7: Ilustration of characteristic points is go out to domain.

Figure 3 . 1 :

 31 Figure 3.1: Examples of non conforming mesh (left) and conforming mesh (right).

Figure 3 . 2 :Definition 2 .

 322 Figure 3.2: Example of an element in isotropic (left) and anisotropic (right) mesh adaptation with circumscribed circle and ellipse.

See figure 3 .

 3 2 and figure 3.3 for illustrations of isotropic and anisotropic mesh adaptation.

Figure 3 . 3 :

 33 Figure 3.3: Example of mesh adaptation about 7000 vertices: isotropic mesh (top, left) and anisotropic mesh (top, right); zoom in the vicinity of the interface triangulations (bottom).

Definition 3 .

 3 Let M be a metric tensor over R d , i.e at each pointx ∈ R d , M (x) is a symmetric, positive definite d × d matrix. The length of a curve γ : [0, 1] → R d , the volume V M (K) of a simplexK, and the distance d M (x, y) between two points x, y ∈ R d with respect to M an be defined as follows, respectively:

Figure 3 . 4 :

 34 Figure 3.4: Geometrical illustration of a unit ellipsoid B M (x 0 ) with an inscribed element K associated to a metric tensor in 3D.

Figure 3 . 5 :

 35 Figure 3.5: Intersection of two metric tensors.

  .8) where: d H (Γ, Γ h ) := max sup x∈Γ inf y∈Γ h |x -y|, sup y∈Γ h inf x∈Γ |x -y| . The estimations (3.7) and (3.8) permit us to control the Hausdoff distance d H (Γ, Γ h ) by the interpolation errors φ

Figure 3 . 6 :

 36 Figure 3.6: Illustration for mesh adaptation in n th iteration for two fluids flow. (a): Initial mesh T n (used to obtain u n ), (b): The new level set fuction φ n with the red line is the zero-level set , (c):

Figure 3 . 7 :

 37 Figure 3.7: Illustration for mesh adaptation in 3D (h min = 1.e -2, h max = 0.2, h grad = 1.3, h hausd = 5e -3). (a): isosurfaces of level set function, (b) and (c): Final adapted mesh where zero-surfaces is explicitly discretized , (d): a cut of final mesh.

.

Contents 4 . 4 . 4

 444 1 The Lid-driven cavity problem . . . . . . . . . . . . . . . . . . . . . . . . 4.1.1 Two-dimensional lid-driven cavity . . . . . . . . . . . . . . . . . . . . . . . 4.1.2 Lid-driven cavity in 3D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Rising bubble . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Rising bubble in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.2 Rising bubble in 3D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Rayleigh-Taylor instability . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 Rayleigh-Taylor instability in 2D . . . . . . . . . . . . . . . . . . . . . . . . 4.3.2 Rayleigh-Taylor instability in 3D . . . . . . . . . . . . . . . . . . . . . . . . The coalescence of two rising bubbles . . . . . . . . . . . . . . . . . . . . 4.4.1 The coalescence of two rising bubbles in 2D . . . . . . . . . . . . . . . . . . 4.4.2 The coalescence of two rising bubbles in 3D . . . . . . . . . . . . . . . . . . 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Four 4 . 2 .

 42 figure 4.2. Our numerical computations in these cases are compared to the results obtained in the very well-known references[START_REF] Shin | High-resolution for incompressible flow using the Navier-Stokes equations and a multigrid method[END_REF], in[START_REF] Gorazd | NSIKE-an incompressible Navier-Stokes solver for unstructured meshes[END_REF], and the benchmark result for cavity flow in[START_REF] Erturk | Numerical solutions of 2D steady incompressible driven cavity flow at high Reynolds numbers[END_REF]

  Two mesh are employed to simulate 3D problem: the coarse one consists of 11037 vertices, 56244 tetrahedrals for the cases of Re = 100, Re = 400 and the other one consists of 35723 vertices, 516 x = 0.515 x = 0.520 x = 0.521 x = 0.531 x = 0.545 y = 0.569 y = 0.564 y = 0.570 y = 0.570 y = 0.562 y = 0.560 Table 4.1: Cavity in 2D: comparison of the positions of the main vortex for different Reynolds. Re = 400 Re = 1000 Re = 5000 Re = 10000

Figure 4 . 1 :

 41 Figure 4.1: Cavity in 2D: Streamlines for different Reynolds number.

Figure 4 . 2 :

 42 Figure 4.2: Cavity in 2D: velocity profile for u x and u y in cases of Re=100 and Re =1000.

Figure 4 . 3 :

 43 Figure 4.3: Cavity in 2D: Isolines of pressure with Re = 10000. Left: result in [VG04]. Center: result in [HRK + 10]. Right: present result.

Figure 4 . 4 :

 44 Figure 4.4: Cavity in 3D: from left to right, streamlines in 2D and in the plane (z = 0.5) of 3D for Re = 400 (top), Re = 1000 (below).

Figure 4 . 5 :

 45 Figure 4.5: Cavity in 3D: velocity profiles on vertical centerline. The first and third rows: present results. The second and fourth rows: results obtained in [KHT87].

[ 4 ,Figure 4 . 6 :

 446 Figure 4.6: Rising bubble in 2D: initialisation (left) and evolution in time: t = 5.0s (center), t = 10.0s (right).

  of [HSL07] in cases of low Reynolds number and Bond number is from 10 to 200 extracted in figure 4.8, corresponding with our differences of coefficient surface tension γ. (see figure (4.1))

Figure 4 . 7 :

 47 Figure 4.7: Rising bubble in 2D: interaction of surface tension on the final bubble (at time t=10s) with different tension coefficients, from left to right: γ = 6e -5; γ = 6e -3; γ = 2.5e -2; γ = 9e -2.

Figure 4 . 8 :

 48 Figure 4.8: Rising bubble in 2D: interaction of surface tension in case of Re = 5 and from left to right: Bo = 10; Bo = 20 ; Bo = 50; Bo = 100; Bo = 200 in [HSL07].

Figure 4 . 9 :

 49 Figure 4.9: Rising bubble in 2D: variation of mass correction with time evolution.

Figure 4 . 10 :

 410 Figure 4.10: Rising bubble in 3D: evolution of the interface in time.

Figure 4 .

 4 Figure 4.11: Rising bubble in 3D: zoom of and adapted mesh at time t = 6.0s.

Figure 4 .

 4 Figure 4.13 shows the results according to different Atwood numbers at the same Reynolds number.

Figure 4 .

 4 Figure 4.12: Rayleigh-Taylor instability in 2D: evolution of the interface in time with At = 0.3.

At = 0. 3 5 Figure 4 .

 354 Figure 4.13: Rayleigh-Taylor instability in 2D: evolution of the interface with different Atwood numbers.

Figure 4 .

 4 Figure 4.14: Rayleigh-Taylor instability in 2D: evolution of the interface with different Atwood numbers and reduced domain.

Figure 4 .

 4 Figure 4.15: Rayleigh-Taylor instability in 2D: evolution of the interface with At = 0.5, Re = 1000.

Figure 4 .

 4 Figure 4.16: Rayleigh-Taylor instability in 2D: extracted results of evolution of the interface with At = 0.5, Re = 1000 in [FGQ01] (top) and [CCG08] (bottom).

Figure 4 .

 4 Figure 4.17: Rayleigh-Taylor instability in 2D: zoom of adapted mesh in the vicinity of the interface, left: anisotropic mesh, right: isotropic mesh.

Figure 4 .

 4 Figure 4.18: Rayleigh-Taylor instability in 3D, evolution of the interface in time with At = 0.5.

Figure 4 .

 4 Figure 4.19: Rayleigh-Taylor instability in 3D: results obtained in [LJG96] (left) and [LFX05] (right).

  t = 0s t = 5.0s t = 10.0s t = 12.6s t = 12.7s t = 13.5s t = 16.0s t = 20.0s

Figure 4 .

 4 Figure 4.20: Coalescence of two bubbles: evolution of the interface in time. The coalescence occurs at t = 12.7s.

Figure 4 .

 4 Figure 4.21: Zoom of mesh and velocity of 2 bubbles at t = 12.6s and t = 12.7s.

Figure 4 .

 4 Figure 4.22 shows the evolution of the interfaces during the rising with the coalescence of two bubbles.

Figure 4 .

 4 Figure 4.22: Coalescence of two bubbles in 3D: evolution of the interfaces in time. The coalescence occurs at t = 4.6s.

Figure 4 .

 4 Figure 4.23: Coalescence of two bubbles in 3D: results obtained in [NGK74].

Figure 4 . 24 :

 424 Figure 4.24: Coalescence of two bubbles in 3D: results obtained in [CL98].

Figure 4 .

 4 Figure 4.25: Coalescence of two bubbles in 3D: results obtained in [dSMN + 04]: (a) 0.0 s; (b) 0.03 s; (c) 0.06 s; (d) 0.09 s; (e) 0.12 s; (f ) 0.15 s.

Figure 4 .Figure 4 . 27 :

 4427 Figure 4.26: Coalescence of two bubbles in 3D: Velocity fields, cutting in the plan y = 0.5.

Figure 4 . 28 :

 428 Figure 4.28: Three categories of structural optimization. (a) Size optimization of a truss structure, (b) shape optimization and (c) topology optimization. The initial problems are shown at the left hand side and the optimal solutions are shown at the right. (from [BS03])

  shape optimization where only the boundary of the domain is optimized, in topology optimization the structure of the domain (the number of holes and the connectivity of the domain for instance) may change during the optimization process. See figure 4.28 for the illustration of three categories of structural optimization. Altogether different problems coming from physics, mechanics, biology, etc... can be cast into the different frameworks of structural optimization. To name a few: -The broad approaches on the design optimization of heat conduction problems can be found in [SK96, Mer98, YPG01, HMY03]. Shape optimization methods for thermal structure have also been utilized for multiobjective function in [PCL03, NPR06], for practical uses. See also in [LSQX99] an extended algorithm for evolutionary structural optimization (ESO) to shape and topology design problems subjected to steady heat conduction. -Many structure topology optimization approaches have been introduced to the field of elastic structure design in recent years, see [SW00, WWG03] for examples. The requiring of structural optimization for a model of linear or nonlinear elasticity with different objective functions have been investigated in [AJT02, AJT04].

  fluid mechanics: Let Ω ∈ R d (d = 2 or 3) be a bounded open set occupied by the fluid governed by Stokes equations:

  famous of them being the level set method. For examples, Wang et al. developed a topology optimization method for linearly elastic structures using an implicit moving boundary [WWG03]. They represented the structural boundary using the level set model embedded in a higher-dimensional scalar function. Allaire et al. proposed a structural optimization method based on a combination of shape derivative and the level set method for front propagation in [AJT04]. In their formulation, the compliance and target displacement objectives are considered. Optimization problems for heat conduction have been solved by level setbased optimization and extended to nonlinear problems by Ha and Cho [HC05, HC08]. Kim et al. developed a level set-based topology optimization method for heat conduction problems utilizing topological derivatives as a nucleation criterion [KHC09]. Recently Duan et al. have applied variational level set methods to shape and topology optimization of fluid flows in [DMZ08b, DMZ08a]. They demonstrated their approach with two-dimensional examples and

Figure 5 . 1 :

 51 Figure 5.1: Variation I + θ of a reference shape Ω.

Following

  the approach of Murat and Simon [MS76, Sim80], starting from a smooth reference open set Ω, we consider domains of the type Ω θ = (I + θ)(Ω) with θ ∈ W 1,∞ (R d , R d ) as been defined in the previous subsection. Let us present some classical definition and results about the differentiation with respect to the domain of functionals of type: Ω → J(Ω) ∈ R.

3 .

 3 0 on ∂Ω (5.21) Hence if we consider the Lagrange function at the equilibrium points of the state and the adjoint problems with the condition σn + λ = 0 on ∂Ω, we obtain the formula:J (Ω)(θ) = ∂L ∂Ω (Ω, u Ω , p Ω , v Ω , q Ω , λ)(θ) (5.22)Applying Lemma 3 and Lemma 4 to this expression leads us to write the shape derivative of the objective function as:J (Ω)(θ) = ∂Ω θ.n j(u) + (-∇u + ∆p -f ).v -qdivu + ∂λ(u -u D ) ∂n + κλ(u -u D ) ds (5.23) Mixed boundary conditionsLet us come to the more general case in which the state system is imposed both Dirichlet and5.1. SHAPE SENSITIVITY ANALYSIS USING HADAMARD'S BOUNDARY VARIATION METHOD123Newmann conditions (∂Ω = Γ D ∪ Γ N ). We consider here the state system:∇p = f in Ω divu = 0 in Ω σ(u, p)n = ϕ N on Γ N u = 0 on Γ D (5.24)and the objective function we consider is the following:J(Ω) = Ω j(u Ω (x))dx + Γ N l(u Ω (x))dsThe Lagrange function can be constructed as:L(Ω, u, p, v, q) = Ω j(u)dx + Γ N l(u)ds + Ω (∇u.∇v -p.divv -q.divu -f.v)dx -Γ D (σ(v, q)n.u + σ(u, p)n.v) -Γ N ϕ N .vdsThe partial differential of L with respect to the each variable read:∂L ∂v (Ω, u, p, v, q), ū = Ω (∇u.∇ū -p.divū -f.ū)dx -Γ D (∂ n ū.u + σ(u, p)n.ū) -+ ∇p -f ).ūdx + Γ N (σ(u, p)n -ϕ N ).ūds -∇v -qdivv)dx -Γ D (σ(v, q)n.v + ∂ n v.v)= Ω (-∆v + ∇q + j (u)).vdx + Γ N (σ(v, q)n + l (u)).vds -

Figure 5 . 2 :

 52 Figure 5.2: Left: explicit discretisation of shape Ω (yellow zone) in mesh of computational domain D. Right: representation of shape Ω as negative subdomain of a level set function defined on D

  smoothness or geometrical or topological constraints, for example, one could require the volume V (Ω) := Ω dx or the perimeter P (Ω) := ∂Ω ds of an admissible shape Ω, to be equal (or lower or equal) to a prescribed upper bound. Here, we limit ourselves to consider the problem with the constraint on the volume of shapes: a variant of (5.39) turns out to be a minimization problem: inf Ω∈U ad (J(Ω) + lV (Ω))(5.40)

  Note that: Result (5.2) of Lemma (3) and formulation (5.34) lead to write:J (Ω) + lV (Ω) (θ) = ∂Ω (w Ω + l)θ.n

  ∀w ∈ H 1 Γ D ∪Γ N (D) , D ( ww + α∇ w.∇wdx) = Γ w Ω wds = J (Ω)(wn)(5.42)Consider now the choice:∀x ∈ D, V (x) = wn(x)

3 . 3 4.∂ 5 .

 335 Choose a descent step τ n > 0, and solve the following level set advection equation onD Ω n : t φ(x, t) + V n (x).∇φ(x, t) for (x, t) ∈ D × (0, τ n ) φ(x, 0) = φ n (x) for x ∈ D (5.43)This produces a new level set function φ n+1 := φ(x, τ n ) associated to the new shape Ω n+1 . Generate the meshed representation of Ω n+1 from the set of data D Ω n , φ n+1 . A new mesh D Ω n+1 of D is produced through an intermedia "ill-shaped" mesh D Ω n+1 in which Ω n+1 is explicitly discretized. D Ω n+1 is remeshed by local modification operators to obtain the resulting "well-shaped" mesh D Ω n+1 for better quality and geometric approximation (see Chapter 3).

Finally, we obtain

  the following result:J (Ω) = Γ (ν∇u.∇u -ν∇u.∇v)θ.nds (5.47)So, the scalar function w Ω (on Γ) in step 1 is:

  3.4 are performed. Each mesh D Ω of D arising in the course of the process has approximately 1200 vertices, and the whole computation takes about 10 minutes. The detail implements in each iteration are displayed on figure 5.4 and the convergence history for the aggregated objective functional (J(Ω) + lV (Ω)) is reported on figure 5.5.

  One numerical example in 2D for Stokes problems with the minimizing energy dissipation has been presented. Regarding the perspectives of this work, we would like to mention a few options: -Investigate the different test cases with Stokes problem and dissipated object function in 2D and 3D, for examples the test cases such as in [BP03, DMZ08b, CG09].

Figure 5 . 4 :

 54 Figure 5.4: Implements in n th -iteration. (a): Explicit discretization of shape Ω n in mesh D Ω n , (b): Velocity field of Stokes problem on sub-mesh T Ω n of shape Ω n , (c): Shape-domain Ω n in mesh D Ω n corresponding level set fuction φ n , (d): Velocity field V n on D Ω n for advection of φ n , (e): Explicit discretization of zero-level set of φ n+1 , the obtained mesh D Ω n+1 is very ill-shaped, (f ): High-quality mesh D Ω n+1 in which the new shape Ω n+1 is explicitly discretized..

Figure 5 . 5 :

 55 Figure 5.5: Deformation process of domain in the case test of elbow. From top to bottom, left to right: 1 th , 50 th , 70 th , 85 th , 100 th and 150 th iterations of the case test.

.

Figure 5 . 6 :

 56 Figure 5.6: Convergence history for the case test to 300 th iteration.

Table 1

 1

	Finite element Scheme of	Method Number of CPU time
		characteristic		iterations
	P 1 bP 1	Euler	Uzawa	70	1
	P 1 bP 1	Euler	Penalty	76	1,2
	P 1 bP 1	Runge Kutta 4 Uzawa	70	1,79
	P 1 bP 1	Runge Kutta 4 Penalty	76	2,13
	P 2 P 1	Euler	Uzawa	76	4,58
	P 2 P 1	Runge Kutta 4 Uzawa	76	7.97

Table 1.1 we see that: the implement of scheme Euler or Runge-Kutta 4 .1: Cavity in 2D: Details on CPU time to reach the steady-state with Re = 1000, uniform mesh of 676 vertices, 1250 elements. (RK4) does not affect the number of iteration to reach the steady-state, on the contrary, the RK4 increases the CPU time in comparison with Euler's scheme. This increasing is explained by the reason that in RK4, at each sub-step δt we have to pass 4 intermediate points to calculate one characteristic point while it requires only one in Euler's scheme (see section 2.3.2 for the details).

Table 1 .

 1

	Finite element Scheme of	Method Number of CPU time
		characteristic		iterations
	P 1 bP 1	Euler	Uzawa	41	1
	P 1 bP 1	Euler	Penalty	44	0,72
	P 1 bP 1	Runge Kutta 4 Uzawa	45	1,16
	P 1 bP 1	Runge Kutta 4 Penalty	47	0.88

2: Cavity in 3D: Details on CPU time to reach the steady-state with Re = 400, uniform mesh of 4978 vertices, 23837 elements.

  Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	3.3.5 The complete adapted mesh strategy for two-phase flows in 3D . . . . . . .
	3.4

Contents 3.1 Basic notations and definitions . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Isotropic and anisotropic mesh adaptation . . . . . . . . . . . . . . . . . . . 3.1.2 The quality of a mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.2 Anisotropic mesh adaptation . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Metric tensor fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Metric definition based on interpolation error . . . . . . . . . . . . . . . . . 3.2.3 Metric intersection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.4 The generation of anisotropic mesh . . . . . . . . . . . . . . . . . . . . . . . 3.2.5 Explicit discretization of the interface . . . . . . . . . . . . . . . . . . . . . 3.2.6 Anisotropic mesh adaptation process for two-phase flow . . . . . . . . . . . 3.3 Discrete three dimensional domain remeshing . . . . . . . . . . . . . . . 3.3.1 The local size map . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Map size determination adapted to the geometric approximation . . . . . . 3.3.3 Gradation of the size map . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.4 Explicit discretization of the zero level set . . . . . . . . . . . . . . . . . . .

  Model of fluid mechanics in the context of Stokes system . . . . . . . . . . 127 5.3.3 Computation of a descent direction . . . . . . . . . . . . . . . . . . . . . . . 128 5.3.4 The proposed algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 5.4 Numerical examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132

Contents 5.1 Shape sensitivity analysis using Hadamard's boundary variation method116 5.1.1 Hadamard's boundary variation method . . . . . . . . . . . . . . . . . . . . 117 5.1.2 Shape differentiability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 5.1.3 Lagrange's approach for the computation of shape derivatives in the context of Stokes system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 5.2 The generic shape optimization algorithm . . . . . . . . . . . . . . . . . 124 5.2.1 The gradient method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 5.2.2 The global algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 5.3 Shape optimization basing on level set method and mesh adaptation . 125 5.3.1 Description of the level set method for shape optimization . . . . . . . . . . 125 5.3.2

Ω ν|∇u| 2 dx

By taking the partial derivative of two parts with respect to domain Ω, we have:

(Ω, u Ω , p Ω , v, q)(θ) (5.29)

(Ω, u Ω , p Ω , v, q)u Ω (θ) + ∂L ∂p (Ω, u Ω , p Ω , v, q)p Ω (θ) (5.30)

Hence, if (v Ω , q Ω ) are solutions of the system of ∂L ∂u (Ω, u Ω , p Ω , v, q) = 0 and ∂L ∂p (Ω, u Ω , p Ω , v, q) which also called adjoint system written as follows:

then the partial differential of objective function with respect to domain Ω can be formulated in a rather simple way:

Applying Lemma 3, Lemma 4 and noting that u = v = 0/Γ D , we have the shape derivative of the objective function in this case:

We have computed the shape derivative for objective functionals. Next, we turn to construct an algorithm for shape optimization problem in the context of Stokes system.

5.2

The generic shape optimization algorithm

The gradient method

As suggested by the name, this method relies on the knowledge of the first-order shape derivative J (Ω) to produce a descent direction θ for J. We have stated that only the normal component of the deformations is considered and play a role in the shape derivative of objective functions. Actually, the shape derivatives of all the functionals we shall get interested in enjoy a structure of the form:

Let us now recast the algorithm 5 in the context of level set method:

Algorithm 6: Shape gradient algorithm (in combination with level set method) n=0, start with an initial guess Ω 0 ;

for n = 1, . . . , untill convergence do 1. Compute the solution (u Ω n , p Ω n ) to the Stokes system on Ω n , and the solution (v Ω n , q Ω n ) to the adjoint state which is constructed in section (5.1.3)

2. Infer a descent direction θ n for the objective functional following the lines of section (5.2.1)

3. Choose a descent step τ n > 0,and solve the Hamilton-Jacobi equation:

to obtain the new level set function φ n .

4. The new shape is obtained as

Model of fluid mechanics in the context of Stokes system

Let us now come back to the problem (5.3). Suppose that the boundary of Ω is decomposed by three disjoint parts:

where Γ = ∅ is the variable part during the optimization process, Γ D is a fixed part on where we impose an inlet (a constant velocity or a profile) and Γ N is also a subset of fixed part where we apply one condition for the outlet (one type of Newmann condition, see chapter 1). As consequence, the velocity u of the fluid is the solution of the following systems:

According to the previous requirements, the considered set U ad of admissible shapes is defined as: where u is the solution of the following the Stokes systems:

The domain Ω is enclosed in a fixed working domain D. Applying the results of shape derivative in the section 5.1.3 for j(u) = ν∇u.∇u, f = 0 and note that only the part Γ of ∂Ω is variable during the optimization process, we have the computation the shape derivative of J(Ω):

The second integral can be reduced as follows:

Hence, we have:

List of Figures