A splitting method for incompressible flows with variable density based on a pressure Poisson equation
Introduction
This paper deals with the numerical approximation of incompressible viscous flows with variable density. This type of flows are governed by the time-dependent Navier–Stokes equations: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 (with d = 2 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 Γ− = {x ∈ Γ:u(x) · n < 0}, with n being the outward unit normal vector. Throughout this paper we assume that the boundary Γ is impermeable, i.e., u · n = 0 everywhere on Γ, and Γ− = ∅.
The mathematical theory of existence and uniqueness for (1.1), (1.2) is quite involved and we refer to Lions [22] for further details. The difficulty comes from the fact that these equations entangle hyperbolic, parabolic, and elliptic features. Approximating (1.1), (1.2) efficiently is a challenging task. A testimony of the difficulty is that, so far, very few papers have been dedicated to the mathematical analysis of the approximation of (1.1), (1.2). We refer to Liu and Walkington [23] for one of the few attempts in this direction.
Approximating (1.1), (1.2) can be done by solving the coupled system (1.1), but this approach is computer intensive due to the elliptic character induced by incompressibility. Alternative, more efficient, approaches advocated in the literature consist of using fractional time-stepping and exploiting, as far as possible, techniques already established for the solution of constant density incompressible fluid flows. The starting point of most fractional time-stepping algorithms consists of decoupling the incompressibility constraint and diffusion in the spirit of Chorin’s [5] and Temam’s [28] projection method. Several algorithms have been developed which extend this idea to the situation that concerns us here, see for example [2], [1], [18], [24]. To the best of our knowledge, Guermond and Quartapelle [18] gave the first stability proof of a projection method for variable density flows. The algorithm proposed in [18] is somewhat expensive since it is composed of two time-consuming projections. An alternative algorithm composed of only one projection per time step was proposed in [24] and proved to be stable. It seems that so far [18], [24] are the only papers where projection methods for variable density flows have been proved to be stable, the best available results being that of Pyo and Shen [24].
The common feature of all the projection-like methods referred to above is that at each time step, say tn+1, the pressure or some related scalar quantity, say Φ, is determined by solving an equation of the following form:where ρn+1 is an approximation of the density at time tn+1 and Ψ is some right-hand side that varies at each time step. The problem (1.3) is far more complicated to solve than just a Poisson equation. It is time consuming since it requires assembling and pre-conditioning a variable-coefficient stiffness matrix at each time step. Note also that, it is necessary to have a uniform lower bound on the value of the density for (1.3) to be solvable. This condition is a key to the method that we propose and it seems that it has not been given enough attention in the literature.
The objective of the present work is to introduce a fractional time-stepping method for solving variable density flows that involves solving only one Poisson problem per time step instead of problems like (1.3). The proposed algorithm is proved to be stable and numerically illustrated.
The paper is organized as follows: in Section 2 we introduce the non-incremental version of our method and prove its stability. The incremental version of the method is introduced in Section 3. First-order Euler time stepping is used in Sections 2 Non-incremental projection method, 3 Incremental projection method. The most accurate version of the method using second-order Backward Second Difference (BDF2) is presented in Section 4. Finally, Section 5 contains some numerical experiments that demonstrate the performance of the method.
Section snippets
Non-incremental projection method
To introduce the main characteristics of the method, we first focus our attention on its simplest form, which, using the terminology from Guermond et al. [13], we henceforth refer to as the non-incremental version. More elaborate versions of this method are introduced in the subsequent sections.
Incremental projection method
It is now well established that the non-incremental pressure correction method is low-order accurate. More precisely, the error is for the velocity in the L2-norm and for the velocity in the H1-norm and the pressure in the L2-norm; see e.g. [25], [26], [17]. We now introduce an incremental version of the method to overcome this accuracy deficiency. We use the same notation as in the previous section and, under the same assumptions on the density, we prove that this method is stable.
BDF second-order rotational projection method
The method presented in the previous section has a second-order splitting error (see [11, Theorem 5.1] for the constant density version of this argument). But, to obtain a scheme of formal second-order accuracy in time it is necessary to replace the two-level semi-implicit time discretization used in the mass conservation and momentum equations with a second-order accurate time stepping method. The purpose of this section is to rewrite the rotational version of the incremental
Convergence tests
In order to test the algorithm (4.1), (4.2), (4.3), (4.4), we consider a problem with a known analytical solution. We solve the variable density Navier–Stokes equations in the unit diskhaving the exact solutionso that the right-hand side to the momentum equation is
Using the second-order scheme
Acknowledgments
The authors are supported by the National Science Foundation Grants DMS-0510650 and DMS-0713829.
References (32)
- et al.
A conservative adaptive projection method for the variable density incompressible Navier–Stokes equations
J. Comput. Phys.
(1998) - et al.
A second-order projection method for variable-density flows
J. Comput. Phys.
(1992) - et al.
Accurate projection methods for the incompressible Navier–Stokes equations
J. Comput. Phys.
(2001) - et al.
An overview of projection methods for incompressible flows
Comput. Methods Appl. Mech. Eng.
(2006) - et al.
Entropy-based nonlinear viscosity for fourier approximations of conservation laws
C.R. Math. Acad. Sci. Paris
(2008) - et al.
Calculation of incompressible viscous flows by an unconditionally stable projection FEM
J. Comput. Phys.
(1997) - et al.
A projection FEM for variable density incompressible flows
J. Comput. Phys.
(2000) - et al.
A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/Least-Squares method for advection-diffusive equations
Comput. Methods Appl. Mech. Eng.
(1989) - et al.
Finite element methods for linear hyperbolic equations
Comput. Methods Appl. Mech. Eng.
(1984) - et al.
Gauge-Uzawa methods for incompressible flows with variable density
J. Comput. Phys.
(2007)
Numerical simulation of Rayleigh–Taylor instability
J. Comput. Phys.
Mixed and Hybrid Finite Element Methods
Numerical solution of the Navier–Stokes equations
Math. Comput.
Numerical methods for convection-dominated diffusion problems based on combining the method of characteristics with finite element or finite difference procedures
SIAM J. Numer. Anal.
Theory and Practice of Finite Elements
Approximation of variable density incompressible flows by means of finite elements and finite volumes
Commun. Numer. Methods Eng.
Cited by (0)
- 1
On long leave from LIMSI (CNRS-UPR 3251), BP 133, 91403, Orsay, France.