A Characteristics-Mix Stabilized Finite Element Method for Variable Density Incompressible Navier-Stokes Equations

This paper describes a characteristics-mix finite element method for the computation of incompressible Navier-Stokes equations with variable density. We have introduced a mixed scheme which combines a characteristics finite element scheme for treating the mass conservation equation and a finite element method to deal with the momentum equation and the divergence free constraint. The proposed method has a lot of attractive computational properties: parameter-free, very flexible, and averting the difficulties caused by the original equations. The stability of the method is proved. Finally, several numerical experiments are given to show that this method is efficient for variable density incompressible flows problem.


Introduction
This paper is devoted to the numerical approximation of incompressible viscous flows with variable density.This type of flows is governed by the time-dependent Navier-Stokes equations [1]  where the dependent variables are the density 0 ρ > , the velocity field u , and the pressure p .The constant ν is the dynamic viscosity coefficient and f is a driving external force.In stratified flows we typically have f g ρ = , where g is the gravity field.The fluid occupies a bounded domain Ω in d  (with 2 d = or 3) and a solution to the above problem is sought over a time interval [ ] 0,T .The Navier-Stokes system is supplemented by the following initial and boundary conditions for u and ρ : where Γ = ∂Ω and − Γ is the inflow boundary, which is defined by , with n being the outward unit normal vector.Throughout this paper, we assume that the boundary Γ is impermeable, i.e., 0 u n ⋅ = everywhere on Γ , and − Γ = ∅ .We note that no initial and boundary condition is needed for the pressure p which can be viewed as a Lagrange multiplier whose mathematical role is to enforce the incompressibility condition.
Compared with the constant density incompressible Navier-Stokes equation, the main difficulty for the simulation of the system Equation ( 1)-( 2) is that these equations entangle hyperbolic, parabolic, and elliptic features.Therefore, how to construct stable and efficient numerical schemes for the system Equations ( 1) and ( 2) is challenging.
For developing numerical approximations to this problem, it seems natural to use, as far as possible, the techniques established for the solution of constant density incompressible Navier-Stokes equations, viz., the fractional step projection method of Chorin [4] [5] and Temam [6] [7].The method uses a time splitting, solving separately the transport equation for the density and the momentum for the velocity, the incompressible constraint being treated through a projection method, see [8].This is the methodology followed in [3] [9]- [11].Several algorithms have been developed which extend this idea to the situation that concerns us here, see for example [1]- [3] [12]- [14].Caterina introduces an hybrid scheme which combines a finite volume approach for treating the mass conservation equation and a finite element method to deal with the momentum equation and the divergence free constraint [15].However, we note also that there is a difficulty which arises with the transport equation during the process of calculation.Because the transport equation has the hyperbolic nature, it is not well adapted to a mere treatment by FE methods, but instead requires a specific approach, like Discontinuous Galerkin methods, artificial viscosity, sub-grid stabilization procedure as in [3], see also [16], or the least-square method as used in [17].Becase the characteristic methods have proved their efficiency for this kind of problems, such as convection-dominated problems [18]- [21], a characteristic stabilized finite element scheme is used to deal with the transport equation in this paper.The idea of characteristic methods is to recast the governing equations in terms of the Lagrangian coordinates defined by the particle trajectories (or characteristics) associated with the problem under consideration.The Lagrangian treatment in these methods greatly reduces the time truncation error in Eulerian method [19].In addition, the characteristic methods have been shown to possess remarkable stability properties.
In this paper, we consider the conserved form for variable density incompressible flows which is introduced by Guermond et al. [1] [3].Because in the formulation, the mass conservation and momentum equations are rewritten in a new form that guarantees some control on the 2 L -norm of the density and on the kinetic energy of the fluid [3].We henceforth refer to the conserved form.The complete system of equations for developing unconditionally stable integration schemes for variable density incompressible flows is written in the following system in conserved form ( ) In the above equations Equation (3), the additional term 2 u ρ ∇ ⋅ and ( ) The originality of our work is to use different numerical methods for the transport equation and for evaluating the evolution of the velocity driven by the last two equations in the system Equation (1).To be more specific, we use a time-splitting, solving the first equation for a given velocity by using a characteristic stabilized finite element approach which is efficient when dealing with a pure convection equation, and then, we compute the divergence free solution of the last two equations by exploiting the advantages of FE methods, see [22]- [24].Here, we care to preserve the divergence free constraint between the two steps of the splitting.The results show that the proposed algorithm is stable and efficient.
The paper is organized as follows.In next section, we introduce some notations for this paper.In Section 3, a detailed presentation of the new method is given.In Section 4, the stability of the method is proved.In Section 5, a series of numerical experiments are given.The last section is devoted to concluding remarks.

Notation
In this section, we aim to describe some of the notations which will be frequently used in this paper.We consider the time-dependent variable density Navier-Stokes system Equations ( 1) and ( 2 No notational distinction is done between scalar or vector-valued functions but spaces of vector-valued functions are identified with bold fonts.The space of functions in ( ) L Ω .We use the standard Sobolev spaces ( ) [25]- [27]).The closure with respect to the norm ( ) The scalar product of ( ) ( ) Ω is denoted , ⋅ ⋅ .We refer readers to [25] [26] for details on these spaces.
For the mathematical setting of problem Equation (1), we introduce the following Hilbert spaces: The spaces W and V are equipped with their usual scalar product and equivalent norm ( ) H Ω , respectively, for 0,1, 2 i = .We define Au u = −∆ .In particular, there holds We also introduce the following bilinear operator: ( ) ( ) ⋅∇ .Moreover, we define the continuous bilinear forms ( ) , a ⋅ ⋅ and ( ) Obviously, the bilinear ( ) , a ⋅ ⋅ is continuous and coercive on V V × and the bilinear ( ) × and satisfies the well-known inf-sup condition [23] [24]: there exists a positive constant 0 0 β > such that for all q M ∈ ( ) where ( ) ( ) Henceforth c denotes a generic constant whose value may change at each occurrence.This constant may depend on the data of the problem and its exact solution, but it does not depend on the discretization parameters or the solution of the numerical scheme.

The Time Splitting Method
as the solution of ( ) Step 2. Solve new velocity and pressure fields: ,

Solving the Density Equation by a Characteristics Finite Element Scheme
The origin of our scheme can be seen by considering the first equation in a space-time framework.First, for density field ρ , we denote by ρ  the material time derivative.It is defined by [28] [29]   ( ) ( ) where X is the motion corresponding to the velocity u and Q its reference map.We recall that, according to the standard formalism of continuum mechanics, ( ) , x X q t = is the position at time t of the material point q , while the reference map ( ) , Q x t yields the material point located at position x at time t .Then grad .u t Let us introduce the characteristic curves, which are simply the the trajectories of the motion associated with the velocity field u .Thus, for given , ; : 0, , , , , which can be obtained by solving the initial value problem .
It represents the trajectory described by a material point that is placed at position x at time t and is driven by the velocity field u .More precisely, ( ) , ; , , x t X Q x t χ τ τ = .By using function χ , we can write an alternative expression for the material time derivative of density field ρ at ( ) q Q x t t x t X q t t x t t For the time variable, the solution will be approximated at times So, let us introduce the following characteristics scheme for time semidiscretization of problem Equation ( 5) Now, multiplying Equation ( 14) by a test function w , integrating in Ω we easily get the following weak formulation for the density equation.Find We define finite element space h W as follows: ( ) where P K is spaces of polynomials with degree 2. Let h u be a given function defined on Ω , so the standard finite element approximation formulation of Equation ( 15) is given as follows: Find

Solving the Velocity Equation by a FE Method
In the numerical simulation of the Navier-Stokes Equation ( 6), a major difficulty is that the velocity and the pressure are coupled by the incompressibility constraint.Many researchers have done a lot of work about Navier-Stokes system with constant density, for example [30]- [33].Here, we can try to use these methods to the variable density Navier-Stokes equation.
Since we aim at using a FE method, it is convenient to write the variational formulation of Equation ( 6).Let 1 n ρ + be a given function defined on Ω .We aim at solving the following problem: , such that for any ( ) ) The domain Ω is approximated by a computational domain h Ω , discretized by a conforming and isotropic set of triangles h  , with mesh-size h .To construct a Galerkin approximation of Equation ( 18), we introduce FE spaces ( ) ( ) for the velocity h u and ( ) for the pressure h p .We choose the pair of spaces ( ) V M which satisfies a discrete inf-sup condition to discretize the velocity and the pressure.So, we define , and 0 , , , where for all ( ) P K are spaces of polynomials with degree 2 and 1, respectively.Now, when ρ + , given the approximations n u , we obtain the following standard finite element(FE) approximation formulation of the equations. , such that for any ( ) .
In the next section, we will prove that the above algorithm is stable.

Stability Analysis of the Method
In this section, we recalls some useful Propositions and stability hypothesis assumptions for the characteristicsmix finite element method for the incompressible flow with variational density [34].The approximation method is used to approximate the mass conservation returns 1 k h ρ + and that this algorithm satisfies the following stability hypothesis: ( ) ( ) Moreover, we also need the following results [34] to prove the stability of the discrete method presented.
Propositon 4.1.For φ , u and v regular enough and v such that Next, we start with the stability proof of the system.To avoid irrelevant technicalities, we assume that there is no external driving force, i.e., 0 f = .Lemma 4.2.Let ( ) be the solution of (19).Then, there holds that ( ) where k k h σ ρ = .
Proof: Taking  19) and using the identity Using (4.1), implies that ( ) Substituting the above inequality into (4.11) and using (4.1), yields ( ) Then using the Poincare inequality, we obtain ( ) Therefore, we can easily get Next, using the above inequality and Poincare inequality again, we have Also, combining (4.14), (4.15) with (4.9), we obtain the desired stability result.

Numerical Simulations
In this section, we present four series of numerical results to illustrate the theoretical analysis of the algorithm proposed in this paper.

Rates of Convergence Study
In order to test the accuracy of the algorithm proposed in this paper, we consider a problem with a known analytical solution.We solve the variable density Navier-Stokes equations Equations ( 1) and (2) in the unit square so that the right-hand side to the momentum equation is ( We use the ( ) , , P P P approximation for the density, the velocity, and the pressure, respectively.We perform the accuracy tests with respect to τ , h and Re.The mesh partition of Ω into triangular element.First, we solve the above mentioned problem for 0.5 T = . The time step is chosen small enough so that the error from the discretization in time is negligible compared to the space error.We give our results for different mesh size max h h = , where max h is the length of the largest edge of the mesh.We considered Re 1000 = .The results are given in Table 1.The accuracy and convergence rate of results are displayed by means of the h .From the Table 1, we can see that we obtained a better convergence rates compared with the results presented in the literature [15].
Secondly, computation are made on a fixed mesh size for different Reynolds number (Re = 1000, 3000, 5000, 8000, 10000).Taking , the results is presented in Figure 1.From the Figure 1, we can see that the stability still keeps well when the Reynolds number increases.These demonstrate that our method is very effective for high Reynolds number.
Next, computation are made on a fixed mesh size and a fixed Reynolds number with different time steps.The computation has been performed for 0 1 t ≤ ≤ .The mesh size is chosen small enough so that the error from the discretization in space is negligible compared to the time stepping error.The convergence results with respect to t ∆ are plotted in Table 2. From the Table 2, we can see that the simulation results coincided with the theory.

Rayleigh-Taylor Instability
In this Subsection we illustrate the performance of the method on a realistic problem, namely we investigate a Rayleigh-Taylor instability.The problem has been considered in [1] [3] [15] starting from the results and comments in [35].We compute the development of a Rayleigh-Taylor instability in the viscous regime as documented in [35]- [38].This problem consists of two layers of fluid initially at rest in the gravity field.It occupies the domain   where 0 µ > is the dynamic viscosity of the fluid (supposed to be constant in the whole domain) and g is the gravitational acceleration.For 0 t > the system evolves under the action of a vertical downward gravity field of intensity g ; the source term in the momentum equation is downward and equal to g ρ .The equations are made dimensionless by using the following references: m ρ for the density, d for lengths, and 1 2 1 2 d g for time.So, the reference velocity is 1 2 1 2 d g .We assume that the symmetry of the ini-tial condition is maintained during the time evolution.The no-slip condition is enforced at the bottom and top walls and symmetry is imposed on the two vertical sides.
Next, we compare the solutions obtained at different Atwood numbers.As the Atwood value increases, the sensitiveness of the calculation to the numerical instabilities grows.The downward motion of the heavy fluid increases with the density difference.The time evolution of the interface of the density field is plotted in Figure 4.It seems that the evolution of the interface configuration does not change significantly.But notice that at 0.9 At = it is very difficult to continue the simulation in the literature.A series of numerical experiments are given to show that this method is highly efficient.

Rising Bubble Test
To investigate the capability of our method to work with larger density variations, we give the computational results for rising bubble test.This simulation is inspired from [15] [36]- [38].A light "droplet" rise through a heavy fluid and impacts into the plane surface of the heavy fluid in a cavity.The computational domain is   The results are displayed in Figure 6.The evolution of the interface at selected points in time.For confirmation, these patterns may be compared to the respective patterns displayed in Figure 9 in [37] and Figure 15 in [39], the numerical results are in very good agreement.

Conclusions
In this paper, we proposed a characteristics-mix finite element method to the case of incompressible viscous flows with variable density.The originality of our approach is to use different numerical methods for the transport equation and evaluating the evolution of the velocity pressure.The new method uses a time splitting, solving separately the transport equation and the momentum equation.To be more specific, we solve the first equation for a given velocity by using a characteristic stabilized finite element approach which is efficient when dealing with a pure convection equation, and then, we compute the divergence free solution of the last two equations by exploiting the advantages of FE methods.The stability proof of the method we proposed for variable density flows was given in the paper.
To verify the correctness of the method, it has been applied to the test cases previously considered in the literature.The spatial approximation is performed by means of Lagrangian finite elements with P2 interpolation for density and velocity and P1 interpolation for pressure.First, the rates of convergence of the method were proved to be in accordance with the theoretical expected ones, leading so to an accurate solver.Then, the simulation of the viscous Rayleigh-Taylor instability was also investigated.We obtained very good results, even for rather high Atwood numbers.Finally, we considered the rising bubble test and sloshing tank to investigate the robustness property of the scheme with regard to high density ratios.The simulation results coincided with the law of physics are very close to the results presented in the literature.Compared with some established methods, the numerical results show that the new method exhibits good stability behavior even large time steps or the high Atwood numbers are used in computation.
of C ∞ -functions compactly supported in Ω is denoted

Figure 1 .
Figure 1.Effect of varying h at different Reynolds number.(a) L 2 error for the density; (b) L 2 error for the velocity; (c) L 2 error for the pressure.
two region with varying density, the heavier fluid superposed to the light one.The interface is slightly smoothed since we set at time 0 t = : of the perturbed interface.The difficulty of the problem essentially depends on:1) the density ratio between the light and the heavy fluid, which is measured by the so-called Atwood number = .The time evolution of the interface of the density field is displayed in Figure 2 at times 0.0, 0.2, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.8.The results are very close to those in [15] [37] [38].Coming to the comparison of the structure, there is satisfactory agreement of the global characteristics of the flow in the early stage and we can observe some slight difference only at large times of the calculation.

Table 1 .
Rates of convergence and error with different mesh size.

Table 2 .
Rates of convergence and error in time.