Received 18 March 2015; accepted 7 April 2015; published 13 April 2015

ABSTRACT

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.

Keywords:

Characteristic Finite Element, Variable Density Flows, Naiver-Stokes Equations, Stabilized Method

1. 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] [2] :

(1)

where the dependent variables are the density, the velocity field, and the pressure. The constant is the dynamic viscosity coefficient and is a driving external force. In stratified flows we typically have, where is the gravity field. The fluid occupies a bounded domain in (with or 3) and a solution to the above problem is sought over a time interval. The Navier-Stokes system is supplemented by the following initial and boundary conditions for and:

(2)

where and is the inflow boundary, which is defined by, with being the outward unit normal vector. Throughout this paper, we assume that the boundary is impermeable, i.e., 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 incom- pressibility condition.

Compared with the constant density incompressible Navier-Stokes equation, the main difficulty for the simu- lation 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 con- straint 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 trans- port 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 Discon- tinuous 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 -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

(3)

In the above equations Equation (3), the additional term and is consistent since it is

zero if.

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.

2. Notation

In this section, we aim to describe some of the notations which will be frequently used in this paper. We con- sider the time-dependent variable density Navier-Stokes system Equations (1) and (2) on the finite time interval and in an open connected and bounded domain (or 3) with boundary, which we assume to be sufficiently smooth. More precisely, we assume that is such that the Stokes operator possesses the usual regularization properties (see [22] - [24] ).

Let be a time step and let us set for. Let be a normed space equipped with the norm. The space of functions such that, the map is -integrable, is indifferently denoted or. For any time-dependent function, we denote,.

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 that have zero average is denoted. We use the standard Sobolev spaces, for and (see [25] -[27] ). The closure with respect to the norm of the space of -functions compactly supported in is denoted

. To simplify the notation, the Hilbert space (resp.) is denoted (resp.

). 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 and are equipped with their usual scalar product and equivalent norm,

for, here, and denote the usual norm and semi norm of the Sobolev space

or, respectively, for. We define. In particular, there holds

.

We also introduce the following bilinear operator:. Moreover, we define the continuous bilinear forms and on and, respectively, by

and a trilinear form on by

Obviously, the bilinear is continuous and coercive on and the bilinear is continuous on and satisfies the well-known inf-sup condition [23] [24] : there exists a positive constant such that for all

(4)

where.

Henceforth 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.

3. Description of the Numerical Scheme

3.1. The Time Splitting Method

Set, , repeat for:

Step 1. Solve new density field:

Find as the solution of

(5)

Step 2. Solve new velocity and pressure fields:

Find as the solution of

(6)

3.2. 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]

(7)

where is the motion corresponding to the velocity and its reference map. We recall that, according to the standard formalism of continuum mechanics, is the position at time of the material point, while the reference map yields the material point located at position at time. Then

(8)

Let us introduce the characteristic curves, which are simply the the trajectories of the motion associated with the velocity field. Thus, for given the characteristic curve through is defined as the vector function

(9)

which can be obtained by solving the initial value problem

(10)

It represents the trajectory described by a material point that is placed at position at time and is driven by the velocity field. More precisely,.

By using function, we can write an alternative expression for the material time derivative of density field at. Indeed, we have

(11)

For the time variable, the solution will be approximated at times,. Throughout this work, we use the standard notation to denote an approximation of.

In order to discretize the material time derivative in Equation (5), we use the following first order backward Euler formula, namely,

(12)

Moreover, for let be defined by

(13)

So, let us introduce the following characteristics scheme for time semidiscretization of problem Equation (5)

(14)

Now, multiplying Equation (14) by a test function, integrating in we easily get the following weak formulation for the density equation.

Find, such that

(15)

We define finite element space as follows:

(16)

where, is spaces of polynomials with degree 2. Let be a given function defined on, so the standard finite element approximation formulation of Equation (15) is given as follows:

Find, such that

(17)

3.3. 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 varia- ble density Navier-Stokes equation.

Since we aim at using a FE method, it is convenient to write the variational formulation of Equation (6). Let be a given function defined on. We aim at solving the following problem:

Find such that for any

(18)

The domain is approximated by a computational domain, discretized by a conforming and isotropic set of triangles, with mesh-size. To construct a Galerkin approximation of Equation (18), we introduce

FE spaces for the velocity and for the pressure. We choose the pair

of spaces which satisfies a discrete inf-sup condition to discretize the velocity and the pressure. So, we define

where for all, and are spaces of polynomials with degree 2 and 1, respectively. Now, when approximating, given the approximations, we obtain the following standard finite element(FE) approximation formulation of the equations.

Find such that for any

(19)

In the next section, we will prove that the above algorithm is stable.

4. Stability Analysis of the Method

In this section, we recalls some useful Propositions and stability hypothesis assumptions for the characteristics- mix finite element method for the incompressible flow with variational density [34] . The approximation method is used to approximate the mass conservation returns and that this algorithm satisfies the following stability hypothesis:

(4.1)

Moreover, we also need the following results [34] to prove the stability of the discrete method presented.

Propositon 4.1. For, and regular enough and such that we have

(4.2)

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.,.

Lemma 4.2. Let be the solution of (19). Then, there holds that

(4.3)

where.

Proof: Taking and in Equation (19) and using the identity

, we obtain

(4.4)

Using (4.1), implies that

(4.5)

Substituting above inequality into (4.4) and using (4.1), yields

So, we have

(4.6)

Similarly,

(4.7)

Substituting above inequality into (4.4) again, leads to

So, we have

which together with (4.6), we can easily obtain the following result

(4.8)

Next, set, , using inf-sup condition yields

and using the (4.8), a simple derivation leads to

(4.9)

Finally, we obtain the desired stability result.

Theorem 4.3. Let be the solution of (19). Then, there holds that

(4.10)

where.

Proof: Taking and in Equation (19) and using the identity

, we obtain

(4.11)

Using (4.1), implies that

(4.12)

Substituting the above inequality into (4.11) and using (4.1), yields

So, we have

(4.13)

Then using the Poincare inequality, we obtain

Therefore, we can easily get

(4.14)

Next, using the above inequality and Poincare inequality again, we have

(4.15)

Also, combining (4.14), (4.15) with (4.9), we obtain the desired stability result.

5. Numerical Simulations

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

5.1. Rates of Convergence Study

In order to test the accuracy of the algorithm proposed in this paper, we consider a problem with a known analy- tical solution. We solve the variable density Navier-Stokes equations Equations (1) and (2) in the unit square in, having the exact solution

so that the right-hand side to the momentum equation is

We use the approximation for the density, the velocity, and the pressure, respectively. We per- form the accuracy tests with respect to, and Re. The mesh partition of into triangular element.

First, we solve the above mentioned problem for. 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, where is the length of the largest edge of the mesh. We considered. The results are given in Table 1. The accuracy and convergence rate of results are displayed by means of the. 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. 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 are plotted in Table 2. From the Table 2, we can see that the simulation results coincided with the theory.

5.2. 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 com- ments in [35] . We compute the development of a Rayleigh-Taylor instability in the viscous regime as docu- mented in [35] - [38] . This problem consists of two layers of fluid initially at rest in the gravity field. It occupies the domain

which splits into two region with varying density, the heavier fluid superposed to the light one. The interface is slightly smoothed since we set at time:

with, and the initial position 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

;

2) the Reynolds number, defined as

;

where is the dynamic viscosity of the fluid (supposed to be constant in the whole domain) and is the gravitational acceleration. For the system evolves under the action of a vertical downward gravity field of intensity; the source term in the momentum equation is downward and equal to.

The equations are made dimensionless by using the following references: for the density, for lengths, and for time. So, the reference velocity is. 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.

・ A low Atwood number problem: Setting,. 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.

・ A high Atwood number problem: Setting,. For this situation, The time evolution of the interface of the density field is plotted in Figure 3 at times 0.2, 0.3, 0.35, 0.4, 0.45, 0.5. Compared with the above test, we can observe the similar structure and the global characteristics of the flow in the early stage. At the same time, we found that the heavy fluid falls faster compared with the low Atwood number problem. The simulation results coincided with the law of physics and are very close to the results presented in the literature [15] [37] [38] .

・ A very high Atwood number problem: Setting,. As the Atwood value increases, the sensitiveness of the calculation to the numerical instabilities grows. The downward mo- tion of the heavy fluid increases with the density difference. The time evolution of the interface of the den- sity field is plotted in Figure 4. It seems that the evolution of the interface configuration does not change sig- nificantly. But notice that 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.

5.3. 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, where and at the fluid is at rest with density:

where. As in [15] [36] , the equations are made dimensionless by using the follow-

ing references: for density, d for length, for time, then, the reference velocity is. In the di- mensionless equations, the gravity term is and the Reynolds number is defined as in above subsection. In our test, the viscosity of the fluid is supposed to be constant in the whole domain and we have.

The results are displayed in Figure 5. The figure contain snapshots of the fluid interface. The snapshots show how the “droplet” travels up through a heavy fluid and merges with a light fluid above. As the “droplet” rise, its shape remains spherical due to the surface tension and the viscosity. As the droplet hits the interface, it merges with the light fluid above and creates waves on the surface. The simulation results are satisfactory agreement with the results presented in the literature [15] [36] [37] .

5.4. Sloshing Tank

To investigate the capability of our method to work with very large density variations, a two-fluid flow in a sloshing tank is considered next. The setup of the test case follows the description [37] [39] . The domain is a container, where and. The interface separating the two phase is initially given as

The densities of the fluids are and, and the dynamic viscosities are

and. The lighter fluid superposed to the heavy one. No surface tension

is considered here, so the volume force is. Slip-boundary conditions are prescribed along

the walls of the tank, and a zero velocity field is initially assumed. The time-step length is and the situation is observed for

The results are displayed in Figure 6. The evolution of the interface at selected points in time. For confir- mation, 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.

6. 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.

References