A splitting method for incompressible flows with variable density based on a pressure Poisson equation

https://doi.org/10.1016/j.jcp.2008.12.036Get rights and content

Abstract

A new fractional time-stepping technique for solving incompressible flows with variable density is proposed. The main feature of this method is that, as opposed to other known algorithms, the pressure is determined by just solving one Poisson equation per time step, which greatly reduces the computational cost. The stability of the method is proved and the performance of the method is numerically illustrated.

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:ρt+·(ρu)=0,ρ(ut+u·u)+p-μΔu=f,·u=0,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 Rd (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 ρ:ρ(x,0)=ρ0(x),ρ(x,t)|Γ-=a(x,t),u(x,0)=u0(x),u(x,t)|Γ=b(x,t),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:-·1ρn+1Φ=Ψ,nΦ|Γ=0,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 O(Δt) for the velocity in the L2-norm and OΔt12 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 diskΩ={(x,y)R2:x2+y2<1},having the exact solutionρ(x,y,t)=2+xcos(sin(t))+ysin(sin(t)),u(x,y,t)=-ycos(t)xcos(t),p(x,y,t)=sin(x)sin(y)sin(t),so that the right-hand side to the momentum equation isf=ρ(x,y,t)(ysin(t)-xcos2(t))+cos(x)sin(y)sin(t)-ρ(x,y,t)(xsin(t)+ycos2(t))+sin(x)cos(y)sin(t).

Using the second-order scheme

Acknowledgments

The authors are supported by the National Science Foundation Grants DMS-0510650 and DMS-0713829.

References (32)

  • G. Tryggvason

    Numerical simulation of Rayleigh–Taylor instability

    J. Comput. Phys.

    (1988)
  • F. Brezzi et al.

    Mixed and Hybrid Finite Element Methods

    (1991)
  • A.J. Chorin

    Numerical solution of the Navier–Stokes equations

    Math. Comput.

    (1968)
  • J. Douglas et al.

    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.

    (1982)
  • A. Ern et al.

    Theory and Practice of Finite Elements

    (2004)
  • Y. Fraigneau et al.

    Approximation of variable density incompressible flows by means of finite elements and finite volumes

    Commun. Numer. Methods Eng.

    (2001)
  • Cited by (0)

    1

    On long leave from LIMSI (CNRS-UPR 3251), BP 133, 91403, Orsay, France.

    View full text