On exactly incompressible DG FEM pressure splitting schemes for the Navier-Stokes equation

TL;DR

This paper compares three pressure correction schemes for Navier-Stokes equations using high-order discontinuous Galerkin methods, highlighting their divergence-free properties and computational efficiency.

Contribution

It provides a detailed comparison of IPCS-D, IPCS-A, and SIMPLE schemes, demonstrating the algebraic pressure correction's advantage in mass conservation and efficiency.

Findings

IPCS-A and SIMPLE are exactly mass conserving.

IPCS-D cannot be exactly divergence free due to stabilization.

IPCS-A is more efficient than SIMPLE in convergence.

Abstract

We compare three iterative pressure correction schemes for solving the Navier-Stokes equations with a focus on exactly divergence free solution with higher order discontinuous Galerkin discretisations. The investigated schemes are the incremental pressure correction scheme on the standard differential form (IPCS-D), the same scheme on algebraic form (IPCS-A), and the semi-implicit method for pressure linked equations (SIMPLE). We show algebraically and through numerical examples that the IPCS-A and SIMPLE schemes are exactly mass conserving due to the algebraic pressure correction, while the IPCS-D scheme cannot be exactly divergence free due to the stabilisation terms required in the pressure Poisson equation. The SIMPLE scheme requires a significantly higher number of pressure correction iterations to obtain converged results than the IPCS-A scheme, so for efficient and mass conserving simulation the IPCS-A method is the best option among the three evaluated schemes.

Equations

$$\rho \left(\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{w}\cdot \nabla )\mathbf{u}\right)$$

$$\nabla \cdot \mathbf{u}$$

$$\left[\begin{array}{cc}\mathbf{A}& \mathbf{B}\\ \mathbf{C}& 0\end{array}\right]\left[\begin{array}{c}\mathbf{u}\\ \mathbf{p}\end{array}\right]=\left[\begin{array}{c}\mathbf{d}\\ \mathbf{e}\end{array}\right],$$

$$\displaystyle{\mathcal{{A}}}\,{\mathord{\color[rgb]{0.67,0.67,0.67}\bullet}}\equiv\frac{\rho}{\Delta t}\gamma_{1}{\mathord{\color[rgb]{0.67,0.67,0.67}\bullet}}+\rho({\boldsymbol{{w}}}\cdot\nabla){\mathord{\color[rgb]{0.67,0.67,0.67}\bullet}}-\nabla\cdot\mu\left(\nabla{\mathord{\color[rgb]{0.67,0.67,0.67}\bullet}}+(\nabla{\mathord{\color[rgb]{0.67,0.67,0.67}\bullet}})^{T}\right),$$

$$\rho \left(\frac{1}{\Delta t}({\gamma }_1\mathbf{u}+{\gamma }_2{\mathbf{u}}^n+{\gamma }_3{\mathbf{u}}^{n-1})+(\mathbf{w}\cdot \nabla )\mathbf{u}\right)$$

$$\nabla \cdot \mathbf{u}$$

$$\mathcal{A}(\mathbf{u},\mathbf{v};\mathbf{w})+\mathcal{B}(p,\mathbf{v})$$

$$\mathcal{C}(\mathbf{u},q)$$

$$\mathbf{u}$$

$$\left\{\left\{u\right\}\right\}$$

$$\displaystyle\left\llbracket\,u\,\right\rrbracket$$

$$\displaystyle\left\llbracket\,{\boldsymbol{{u}}}\,\right\rrbracket_{\boldsymbol{{n}}}$$

$${\int }_{\mathcal{T}}\frac{\rho }{\Delta t}({\gamma }_1\mathbf{u}+{\gamma }_2{\mathbf{u}}^n+{\gamma }_3{\mathbf{u}}^{n-1})\mathbf{v}\mathrm{d}\mathbf{x}$$

$$\displaystyle\qquad-\ \int_{\mathcal{T}}{\boldsymbol{{u}}}\cdot\nabla\cdot(\rho{\boldsymbol{{v}}}\otimes{\boldsymbol{{w}}})\,\mathrm{d}{\boldsymbol{{x}}}\ +\ \int_{\mathcal{S}}{\boldsymbol{{w}}}\cdot{\boldsymbol{{n}}}^{+}\,\hat{\boldsymbol{{u}}}^{\boldsymbol{{w}}}\cdot\left\llbracket\,\rho{\boldsymbol{{v}}}\,\right\rrbracket\,\mathrm{d}s$$

$$\displaystyle\qquad+\ \int_{\mathcal{T}}\mu\left(\nabla{\boldsymbol{{u}}}+(\nabla{\boldsymbol{{u}}})^{T}\right):\nabla{\boldsymbol{{v}}}\,\mathrm{d}{\boldsymbol{{x}}}\ +\ \int_{\mathcal{S}_{I}}\kappa_{\mu}\left\llbracket\,{\boldsymbol{{u}}}\,\right\rrbracket\cdot\left\llbracket\,{\boldsymbol{{v}}}\,\right\rrbracket\,\mathrm{d}s$$

$$\displaystyle\qquad-\ \int_{\mathcal{S}}(\left\{\!\left\{\mu\left(\nabla{\boldsymbol{{u}}}+(\nabla{\boldsymbol{{u}}})^{T}\right)\right\}\!\right\}\cdot{\boldsymbol{{n}}}^{+})\cdot\left\llbracket\,{\boldsymbol{{v}}}\,\right\rrbracket\,\mathrm{d}s$$

$$\displaystyle\qquad-\ \int_{\mathcal{S}_{I}}(\left\{\!\left\{\mu\left(\nabla{\boldsymbol{{v}}}+(\nabla{\boldsymbol{{v}}})^{T}\right)\right\}\!\right\}\cdot{\boldsymbol{{n}}}^{+})\cdot\left\llbracket\,{\boldsymbol{{u}}}\,\right\rrbracket\,\mathrm{d}s$$

$$\displaystyle\qquad-\ \int_{\mathcal{T}}p\,\nabla\cdot{\boldsymbol{{v}}}\,\mathrm{d}{\boldsymbol{{x}}}\ +\ \int_{\mathcal{S}}\hat{p}\,{\boldsymbol{{n}}}^{+}\cdot\left\llbracket\,{\boldsymbol{{v}}}\,\right\rrbracket\,\mathrm{d}s+\frac{1}{2}\int_{\mathcal{T}}(\nabla\cdot{\boldsymbol{{w}}}){\boldsymbol{{u}}}\cdot{\boldsymbol{{v}}}\,\mathrm{d}{\boldsymbol{{x}}}=0.$$

$${\kappa }_{\mu }=3\frac{{\mu }_{\max }^2}{{\mu }_{\min }}k(k+1)\underset{K}{\max }\left(\frac{S_K}{V_K}\right),$$

$$\frac{1}{2}{\int }_{\mathcal{T}}(\nabla \cdot \mathbf{w})\mathbf{u}\cdot \mathbf{v}\mathrm{d}\mathbf{x}.$$

$$\displaystyle\int_{\mathcal{S}}\hat{{\boldsymbol{{u}}}}^{p}\cdot{\boldsymbol{{n}}}^{+}\left\llbracket\,q\,\right\rrbracket\,\mathrm{d}s\ -\ \int_{\mathcal{T}}{\boldsymbol{{u}}}\cdot\nabla q\,\mathrm{d}{\boldsymbol{{x}}}\ =\ 0.$$

$$\nabla \cdot \left(\frac{1}{\rho }\nabla p\right)=f.$$

$$-{\int }_{\mathcal{T}}\nabla p\cdot \nabla q\mathrm{d}\mathbf{x}+{\int }_{{\mathcal{S}}_I}q\nabla p\cdot \mathbf{n}\mathrm{d}s$$

$$\displaystyle+\int_{\mathcal{S}_{I}}\left\llbracket\,q\,\right\rrbracket\left\{\!\left\{\nabla p\right\}\!\right\}\cdot{\boldsymbol{{n}}}^{+}\,\mathrm{d}s+\int_{\mathcal{S}_{I}}\left\llbracket\,\nabla p\,\right\rrbracket\left\{\!\left\{q\right\}\!\right\}\cdot{\boldsymbol{{n}}}^{+}\,\mathrm{d}s$$

$$\displaystyle+\int_{\mathcal{S}_{I}}\left\llbracket\,p\,\right\rrbracket\left\{\!\left\{\nabla q\right\}\!\right\}\cdot{\boldsymbol{{n}}}^{+}\,\mathrm{d}s-\int_{\mathcal{S}_{I}}\kappa_{p}\left\llbracket\,p\,\right\rrbracket\left\llbracket\,q\,\right\rrbracket\,\mathrm{d}s=\int_{\mathcal{T}}fq\,\mathrm{d}{\boldsymbol{{x}}},$$

$$\rho \left(\frac{1}{\Delta t}({\gamma }_1{\mathbf{u}}^{\ast }+{\gamma }_2{\mathbf{u}}^n+{\gamma }_3{\mathbf{u}}^{n-1})+(\mathbf{w}\cdot \nabla ){\mathbf{u}}^{\ast }\right)$$

$$(\mathbf{w}\cdot \nabla )(\mathbf{u}-{\mathbf{u}}^{\ast })$$

$$\nabla \cdot \mu \left(\nabla (\mathbf{u}-{\mathbf{u}}^{\ast })+(\nabla (\mathbf{u}-{\mathbf{u}}^{\ast }){)}^T\right)$$

$$\frac{\rho }{\Delta t}({\gamma }_1\mathbf{u}-{\gamma }_1{\mathbf{u}}^{\ast })=-\frac{1}{\rho }\nabla (p-p^{\ast }).$$

$$\nabla \cdot \left(\frac{\Delta t}{{\gamma }_1\rho }\nabla (p-p^{\ast })\right)=\nabla \cdot {\mathbf{u}}^{\ast },$$

$$\mathbf{A}{\mathbf{u}}^{\ast }=\mathbf{d}-\mathbf{B}{\mathbf{p}}^{\ast }.$$

$$\mathbf{M}(\mathbf{u}-{\mathbf{u}}^{\ast })=-\mathbf{B}(\mathbf{p}-{\mathbf{p}}^{\ast }).$$

$$\mathbf{e}-\mathbf{C}{\mathbf{u}}^{\ast }=-\mathbf{C}{\mathbf{M}}^{-1}\mathbf{B}(\mathbf{p}-{\mathbf{p}}^{\ast }).$$

$$\mathbf{C}{\mathbf{M}}^{-1}\mathbf{B}\mathbf{p}=\mathbf{C}{\mathbf{M}}^{-1}\mathbf{B}{\mathbf{p}}^{\ast }-\mathbf{e}+\mathbf{C}{\mathbf{u}}^{\ast }.$$

$$\mathbf{u}=-{\mathbf{M}}^{-1}\mathbf{B}(\mathbf{p}-{\mathbf{p}}^{\ast }).$$

$$\mathbf{A}{\mathbf{u}}^{\ast }$$

$$\mathbf{C}{\mathbf{u}}^{\ast }$$

Full text

We show algebraically and through numerical examples that the IPCS-A and SIMPLE schemes are exactly mass conserving due to the algebraic pressure correction, while the IPCS-D scheme cannot be exactly divergence free due to the stabilisation terms required in the pressure Poisson equation. The SIMPLE scheme requires a significantly higher number of pressure correction iterations to obtain converged results than the IPCS-A scheme, so for efficient and mass conserving simulation the IPCS-A method is the best option among the three evaluated schemes.

1 Introduction

The incompressible Navier-Stokes equations contain three major obstacles to constructing an efficient solution procedure: the non-linearity of the convective term, the incompressibility criterion, and the intricate pressure-velocity coupling. In this work we will focus on how to maintain exact incompressibility when dealing with the lack of a pressure update equation. The non-linearity will simply be handled by introducing an explicitly extrapolated convecting velocity.

Denoting the unknown velocity vector as ${\boldsymbol{{u}}}$, the explicit convecting velocity as ${\boldsymbol{{w}}}$, the pressure as $p$, the fluid density as $\rho$ and the fluid dynamic viscosity as $\mu$, the incompressible Navier-Stokes equations can be written

[TABLE]

which, after discretisation, leads to a sparse discrete block matrix problem,

[TABLE]

where the sub-matrix ${\mathbf{{A}}}$ is a discrete version of the bilinear operator ${\mathcal{{A}}}$,

[TABLE]

which has been assembled by use of the discontinuous Galerkin finite element (DG FEM) scheme described in section 2. Similar holds for the matrices ${\mathbf{{B}}}$ and ${\mathbf{{C}}}$, the discrete versions of the pressure gradient and velocity divergence operators. The unknowns ${\mathbf{{u}}}$ and ${\mathbf{{p}}}$ (bold, non-italic) are now arrays of unknowns containing the degrees of freedom describing the velocity and pressure fields. The vectors ${\mathbf{{d}}}$ and ${\mathbf{{e}}}$ contain the assembled linear operators from the momentum and continuity equations respectively. Note that due to boundary conditions the vector ${\mathbf{{e}}}$ is not necessarily zero.

There are many methods for solving equation 2. With appropriate boundary and initial conditions the equations are well-formed and can be solved directly on coupled form by an LU decomposition method. Fast parallel sparse LU solvers exist, such as MUMPS, which uses multifrontal LU-factorization (Amestoy et al., 2001, 2006), and SuperLU_DIST, which uses supernodal LU-factorization (Li, 2005; Li et al., 1999; Li and Demmel, 2003). Unfortunately, such direct solution methods only scale up to a relatively small number of parallel processors, but they can deal with the saddle point block matrix problems such as equation 2 without any special preconditioning.

Iterative Krylov subspace methods can scale to thousands of processors, but they require preconditioning for efficient and stable solution of saddle point problems (Hestenes and Stiefel, 1952; Fletcher, 1976; Saad and Schultz, 1986; van der Vorst, 1992). Pressure correction methods are a popular family of methods to deal with the saddle point problem in the Navier-Stokes equations by splitting the equation into individual momentum guess and pressure correction steps, which can then be solved efficiently by sparse iterative Krylov solvers. The earliest pressure correction schemes by Chorin (1968) and Temam (1969) are still popular, often written on iterative form by a method such as the Incremental Pressure Correction Scheme, IPCS, (Goda, 1979). Many other splitting methods exist, such as velocity correction methods (Kawahara and Ohmiya, 1985; Guermond and Shen, 2003), and the more general Schur complement methods (Schur, 1917; Zhang, 2005). The popular PISO (Issa, 1986) and PIMPLE (Weller et al., 1998) algebraic pressure correction techniques should also be mentioned, and we even implemented them in the Ocellaris flow solver that we have used in this work. But, though they are exactly mass conserving, unfortunately we were unable to make them converge properly with the higher order DG discretisation, so they have been left out of this comparison.

In this paper we will look at the IPCS method on differential form in section 3 and on algebraic form in section 4. The DG SIMPLE method by Klein et al. (2013, 2015, 2016) based on the SIMPLE method by Patankar and Spalding (1972) is described in section 5. The effect of the choice of pressure splitting scheme on the divergence of the resulting velocity field is discussed in section 6, before the methods are compared on several benchmark cases in section 7. The DG FEM spatial discretisation and the temporal scheme used in the numerical experiments are presented in section 2. The findings in this paper are not dependent on the details of the numerical scheme and the results will hold for other methods as long as DG elements are used for the pressure and exact incompressibility can be achieved when the governing equations are solved without pressure-velocity splitting.

2 The DG FEM discretisation

This section contains a brief summary of the DG FEM discretisation of the governing equations (1) used to construct the block matrix in equation 2. More details can be found in Landet et al. (2018) which is based on Cockburn et al. (2005) except for the elliptic term where we have used the symmetric interior penalty (SIP) method by Arnold (1982) and not the local discontinuous Galerkin (LDG) method used by Cockburn et al. The exact details of the presented discretisation are not very important for the focus of this paper. The numerical scheme can be substituted for another numerical scheme as long as (i) the resulting mass matrix is block diagonal, (ii) the pressure is discontinuous and requires continuity-enforcing stabilisation of elliptic operators, and (iii) the velocity field is machine precision divergence free when equation 2 is solved on coupled form by a direct solver.

The Navier-Stokes equations from equation 1 are first discretised in time by use of values from previous time steps, $({\boldsymbol{{u}}}^{n},{\boldsymbol{{u}}}^{n-1})$, along with a suitable choice of time stepping coefficients $(\gamma_{1},\gamma_{2},\gamma_{3})$. The time steps are $t=n\Delta t$ with initial conditions at $t=0$. We have used a second order backwards differencing formulation, BDF2, where the coefficients are $(\gamma_{1},\gamma_{2},\gamma_{3})=(\nicefrac{{3}}{{2}},\,-2,\,\nicefrac{{1}}{{2}})$. If the initial value at $t=-\Delta t$ is not provided then coefficients $(1,\,-1,\,0)$ is used for the first time step. The non-linearity of the convective term is handled semi-implicitly by introducing an explicit convecting velocity, ${\boldsymbol{{w}}}=2{\boldsymbol{{u}}}^{n}-{\boldsymbol{{u}}}^{n-1}$. The spatial differential equations are then

[TABLE]

Let $\mathcal{T}$ be the set of all cells in the mesh and $\mathcal{S}$ the set of all facets. Polynomial function spaces of degree $k$ on each cell $K$ are denoted $P_{k}(K)$. These have no continuity at cell boundaries and no inherent boundary conditions. We use calligraphic typeface to denote operators and sets, bold italic for vectors functions and italic for scalar functions. Equation 4 is cast into the following form: find ${\boldsymbol{{u}}}\in[P_{2}(K)]^{D}$ and $p\in P_{1}(K)$ such that

[TABLE]

in the tessellated domain $\mathcal{T}$ subject to

[TABLE]

The discontinuous Galerkin method works by breaking integrals over the whole domain into a sum of integrals over each mesh cell $K\in\mathcal{T}$, and defining fluxes of the unknown functions between these cells. The average and jump operators across an internal facet between two cells $K^{+}$ and $K^{-}$ are defined as

[TABLE]

The convecting velocity ${\boldsymbol{{w}}}$ is Hdiv-conforming, i.e. the flux is continuous, \left\llbracket\,{\boldsymbol{{w}}}\,\right\rrbracket_{\boldsymbol{{n}}}=0. For exterior facets, let the connected cell be denoted $K^{+}$ such that {\boldsymbol{{n}}}^{+}\!\cdot\left\llbracket\,{\boldsymbol{{u}}}\,\right\rrbracket={\boldsymbol{{n}}}^{+}\!\cdot{\boldsymbol{{u}}}^{+}={\boldsymbol{{n}}}\cdot{\boldsymbol{{u}}}. Take $\left\{\!\left\{u\right\}\!\right\}=u$ and otherwise let all $K^{-}$ values be zero on the boundary facets.

The momentum equation is discretised using the symmetric interior penalty (SIP) method for the elliptic term (Arnold, 1982) and otherwise using the fluxes from Cockburn et al. (2005). The flux of pressure is $\hat{p}=\left\{\!\left\{p\right\}\!\right\}$ and the convective flux $\hat{\boldsymbol{{u}}}^{\boldsymbol{{w}}}$ is a pure upwind flux. After multiplication with ${\boldsymbol{{v}}}$ and integration over $\mathcal{T}$, followed by integration by parts and application of the SIP method to the viscosity; ${\mathcal{{A}}}$ can be found as the bilinear part containing ${\boldsymbol{{u}}}={\boldsymbol{{u}}}^{n+1}$, ${\mathcal{{B}}}$ as the bilinear part containing $p$, and ${\mathcal{{D}}}$ as the linear part of

[TABLE]

On Dirichlet boundaries $\hat{\boldsymbol{{u}}}^{\boldsymbol{{w}}}$ is the upwind value. Depending on flow direction this is either ${\boldsymbol{{u}}}$ or ${\boldsymbol{{u}}}_{D}$. On all boundaries take $\hat{p}=p$. The SIP penalty parameter is selected based on Epshteyn and Rivière (2007); Shahbazi et al. (2007),

[TABLE]

where the $k$ is the order of the approximating polynomials, $S_{K}$ is the cell surface area and $V_{K}$ the cell volume. The penalty parameter is doubled on external facets.

Since the IPCS-D scheme is not exactly divergence free, the standard skew symmetric term (Gresho, 1991) is included in the weak form used for all the methods. The term is the last integral in equation 11,

[TABLE]

The continuity equation is multiplied by $q$ and integrated over $\mathcal{T}$. The flux of velocity is here $\hat{{\boldsymbol{{u}}}}^{p}=\left\{\!\left\{{\boldsymbol{{u}}}\right\}\!\right\}$. After integration by parts ${\mathcal{{C}}}({\boldsymbol{{u}}},q)$ and ${\mathcal{{E}}}(q)$ can be found from

[TABLE]

The non-zero ${\mathcal{{E}}}(q)$ results from using the boundary conditions in the flux, $\hat{\boldsymbol{{u}}}^{p}={\boldsymbol{{u}}}_{D}$ on Dirichlet boundaries.

After solving equations 5 and 6 by one of the presented pressure splitting schemes, the resulting velocity field is projected into a space where it is exactly divergence free and where the normal velocities across internal facets are continuous. This projection from weak to strong incompressibility requires only that the discrete version of the weak incompressibility criterion on the form shown in equation 14 is satisfied exactly. The projection is described in Cockburn et al. (2005). It is not necessary to perform the projection inside the pressure iteration loops, it is run once as a post processing operation at the end of each time step.

2.1 Differential Poisson equation for the pressure

In the IPCS-D method, which will be described in section 3, a Poisson equation must be solved for the pressure,

[TABLE]

This equation is discretised with the same SIP penalty method used for the viscous term in the momentum equation and the same penalty parameter as in equation 12 is used, with $\mu$ replaced by $\rho^{-1}$. The applied boundary condition is $\nabla p\cdot{\boldsymbol{{n}}}=0$ on all boundaries. After multiplying by the test function $q\in P_{1}[K]$, integrating over the domain, performing integration by parts, and including the SIP stabilisation terms, the result is

[TABLE]

where $\mathcal{S}_{I}$ is the set of all internal facets. The null space due to using pure Neumann boundary conditions is removed in the Krylov solver by providing the basis of the null space along with the system matrix (Balay et al., 2014). The right-hand side will be the divergence of a known velocity field. This will be discretised as shown in equation 14 to minimise the differences between the differential and the algebraic methods.

The main difference between solving for the pressure this way, by a forming an elliptic differential equation, instead of performing the elliptic-like algebraic operation used in the IPCS-A method described in section 3, is the introduction of the penalty parameter. The need for stabilisation of elliptic operators in DG methods means that the resulting discrete elliptic matrices are different between the algebraic and the differential methods. The consequences of this are explored in section 6 and shown numerically in section 7.

3 The IPCS-D method

The Incremental Pressure Correction Scheme on Differential form (IPCS-D) is an iterative version of the classic Chorin-Temam pressure correction scheme (Chorin, 1968; Temam, 1969; Goda, 1979). The method begins with solving the momentum equation for an approximate velocity field ${\boldsymbol{{u}}}^{*}$ by use of a guessed pressure field, $p^{*}$,

[TABLE]

A splitting error is now introduced. Make the assumption that the true velocity ${\boldsymbol{{u}}}$ and the velocity estimate ${\boldsymbol{{u}}}^{*}$ are close enough so that

[TABLE]

By subtracting equation 17 from the true momentum equation 4 and using the assumptions in equations 18 and 19 the result is

[TABLE]

The velocity field should be divergence free, but the velocity estimate ${\boldsymbol{{u}}}^{*}$ is not due to having used an approximate pressure field $p^{*}$, and not the true velocity field $p$. The continuity equation, $\nabla\cdot{\boldsymbol{{u}}}=0$, is now used to eliminate the unknown velocity field ${\boldsymbol{{u}}}$ from equation 20 by taking the divergence,

[TABLE]

which gives a Poisson equation for the unknown pressure $p$ with the guessed pressure $p^{*}$ and the velocity estimate ${\boldsymbol{{u}}}^{*}$ as known coefficients. When the pressure $p$ is found it can be put back into equation 20 to find the updated velocity field, ${\boldsymbol{{u}}}$.

The IPCS-D algorithm

Guess $p^{*}$ and solve equation 17 for ${\boldsymbol{{u}}}^{*}$. 2. 2.

Solve equation 21 for $p$. 3. 3.

Compute the updated velocity ${\boldsymbol{{u}}}$ from equation 20. 4. 4.

Check for convergence by computing the residual $r_{u}=||{\mathbf{{u}}}^{*}-{\mathbf{{u}}}||$ and optionally go to step 1 using $p$ as a new guess for the pressure if $r_{u}$ is not sufficiently small.

4 The IPCS-A method

The Incremental Pressure Correction Scheme on Algebraic form (IPCS-A) starts from the Navier-Stokes equations on block matrix form, as can be seen in equation 2. The momentum equation can be solved for an approximate velocity field ${\mathbf{{u}}}^{*}$ by use of a guessed pressure field ${\mathbf{{p}}}^{*}$,

[TABLE]

A splitting error is now introduced. First let ${\mathbf{{A}}}={\mathbf{{M}}}+{\mathbf{{R}}}$ where ${\mathbf{{M}}}$ is a scaled mass matrix resulting from the assembly of the first term in equation 3, and ${\mathbf{{R}}}$ contains the convective and diffusive operators. Then make the assumption that ${\mathbf{{R}}}({\mathbf{{u}}}-{\mathbf{{u}}}^{*})\approx 0$ and use this when subtracting equation 22 from the first line of equation 2, giving

[TABLE]

The matrix ${\mathbf{{M}}}$ is block diagonal and can easily be inverted. Use this property and the divergence free criterion, ${\mathbf{{C}}}{\mathbf{{u}}}={\mathbf{{e}}}$, to remove the unknown ${\mathbf{{u}}}$ from equation 23

[TABLE]

Reorganising equation 24 results in an equation for ${\mathbf{{p}}}$

[TABLE]

The velocity ${\mathbf{{u}}}$ can be recovered without solving a linear system, simply by substituting the pressure ${\mathbf{{p}}}$ from the solution of equation 25 into equation 23,

[TABLE]

The IPCS-A algorithm

Guess ${\mathbf{{p}}}^{*}$ and solve equation 22 for ${\mathbf{{u}}}^{*}$. 2. 2.

Solve equation 25 for ${\mathbf{{p}}}$. 3. 3.

Compute the updated velocity ${\mathbf{{u}}}$ from equation 26. 4. 4.

Check for convergence by computing the residual $r_{u}=||{\mathbf{{u}}}^{*}-{\mathbf{{u}}}||$ and optionally go to step 1 using ${\mathbf{{p}}}$ as a new guess for the pressure if $r_{u}$ is not sufficiently small.

5 The SIMPLE method

This description is a summary of the method presented in Klein et al. (2013), the Semi-Implicit Method for Pressure-Linked Equations, SIMPLE. The method starts with the Navier-Stokes equations on algebraic block matrix form, equation 2. A guess ${\mathbf{{p}}}^{*}$ is made and inserted into the governing equations which can be solved for an estimated velocity field, ${\mathbf{{u}}}^{*}$,

[TABLE]

where ${\mathbf{{e}}}^{*}$ is not zero when ${\mathbf{{p}}}^{*}$ is not a perfect guess due to $\nabla\cdot{\boldsymbol{{u}}}^{*}\neq 0$. Subtracting these equations from equation 2 and defining corrections ${\hat{\mathbf{{u}}}}={\mathbf{{u}}}-{\mathbf{{u}}}^{*}$ and ${\hat{\mathbf{{p}}}}={\mathbf{{p}}}-{\mathbf{{p}}}^{*}$ gives

[TABLE]

A splitting error is now introduced. Construct a matrix $\tilde{{\mathbf{{A}}}}\approx{\mathbf{{A}}}$ where $\tilde{{\mathbf{{A}}}}$ approximates ${\mathbf{{A}}}$, but is much easier to invert. A diagonal or block diagonal version of ${\mathbf{{A}}}$ is used as an approximation in the numerical examples below. Note that this approximation allows the time derivative to be removed for steady-state problems, unlike in the IPCS-A method where we must require ${\mathbf{{M}}}\neq 0$. The velocity correction in equation 29 can then be approximated as

[TABLE]

and this approximation can be used to solve for ${\hat{\mathbf{{p}}}}$ by substitution into equation 30

[TABLE]

The SIMPLE method is not guaranteed to converge without under-relaxation in the updates of ${\mathbf{{p}}}$ and ${\mathbf{{u}}}$ due to the approximate $\tilde{{\mathbf{{A}}}}$. The under-relaxation of the pressure is performed explicitly,

[TABLE]

while an implicit scheme is used for the velocity,

[TABLE]

with $0<(\alpha_{p},\alpha_{u})<1$.

The SIMPLE algorithm

Solve for ${\mathbf{{u}}}^{*}$ using equation 34 starting with guesses ${\mathbf{{u}}}^{*}_{\text{prev}}$ and ${\mathbf{{p}}}^{*}$ 2. 2.

Find ${\hat{\mathbf{{p}}}}$ using equation 32 3. 3.

Update ${\mathbf{{p}}}$ using equation 33 4. 4.

Update ${\mathbf{{u}}}$ using equation 31 5. 5.

Check for convergence by computing the residual $r_{u}=||{\mathbf{{u}}}^{*}-{\mathbf{{u}}}||$ and optionally go to step 1 with new guesses for ${\mathbf{{u}}}^{*}$ and ${\mathbf{{p}}}^{*}$

6 Exact mass conservation

The numerical representation of the velocity field will be divergence free to the precision of the discrete operators if ${\mathbf{{C}}}{\mathbf{{u}}}-{\mathbf{{e}}}=0$. In the IPCS-A method this is ensured on the algebraic level,

[TABLE]

The IPCS-D pressure Poisson equation 21 is not equivalent to equation 24 due to the stabilisation jump terms which will, even if the boundary conditions applied to equation 21 are in perfect agreement with equation 25, introduce a residual divergence in the velocity field since the assembled matrix from the elliptic operator with stabilisation does not satisfy equation 24,

[TABLE]

Other elliptic DG FEM discretisations like LDG and NIPG will also contain stabilisation terms (Arnold et al., 2002), and will hence have the same problem. A similar argument is shown in Klein et al. (2015). They also show that the SIMPLE method is exactly divergence free with an argument similar to equation 35.

7 Numerical experiments

The weak form described in section 2 is implemented with quadratic discontinuous Galerkin elements for the velocity (DG2) and linear DG elements for the pressure (DG1). The numerical experiments are run with iterative solvers from PETSc 3.8 (Balay et al., 2014). All results are from simulations running in parallel with MPI on the Abel computational cluster at the University of Oslo. The GMRES iterative solver has been used for both velocities and pressures and both the relative and absolute error criteria are $1\text{\times}{10}^{-15}$ with a maximum of $100$ Krylov iterations. Convergence in the Krylov solver is typically achieved in less than $100$ Krylov iterations for all but the first pressure correction iteration in the first time step.

An Additive Schwarz preconditioner (ASM) is used for the velocities and for the pressure an Algebraic Multigrid preconditioner (HYPRE BoomerAMG) is applied. The solvers are set up to use the previous solution as the initial guess. The Additive Schwarz preconditioner uses PETSc’s default settings, which is incomplete LU factorization on each CPU. The BoomerAMG preconditioner also uses PETSc’s default settings.

All the results shown below have been produced by Ocellaris (Landet, 2019a), a single and multi-phase DG FEM Navier-Stokes solver. Input files for Ocellaris and post-processing scripts, which can be used to reproduce all figures shown in this paper, can be found in Landet (2019b).

7.1 Taylor-Green 2D flow

The Taylor-Green vortex is an analytical solution to the 2D incompressible Navier-Stokes equations. The solution is

[TABLE]

Spatial convergence of the pressure correction methods is examined on a regular grid. The domain, $(x,y)\in[0,2]\times[0,2]$, is divided into $N\times N$ rectangles that are then subdivided into two triangles each. Grid sizes $N=\{8,16,24,32\}$ are used to study the spatial convergence rate. The time step is $\Delta t=0.01$ and the numerical results are compared to the analytical expressions at $t=1.0$. Exactly 160 pressure correction iterations are performed for each time step, which is sufficient for all methods, see figure 1. The physical parameters applied are $\rho=1.0$ and $\nu=0.005$.

In the SIMPLE method an under-relaxation factor of $\alpha_{u}=0.7$ is used for the velocity and $\alpha_{p}=1.0$ is used for the pressure. This is in line with Klein et al. (2013). Setting both under-relaxation factors to $1.0$ (no under-relaxation) as was done in Klein et al. (2015) is unstable for the current simulations and causes the iteration procedure to blow up. A diagonal $\tilde{{\mathbf{{A}}}}$ matrix is used with no lumping.

The number of pressure correction iterations required to reach convergence is shown in figure 1. The SIMPLE method requires far more iterations in order to converge than the IPCS methods. The influence of some of the tunable parameters in the SIMPLE method are also shown in the figure. Some sensitivity is found both in regard to how the approximate inverse of ${\mathbf{{A}}}$ is computed and in the choice of under-relaxation parameters. The effect of these choices are significantly smaller than the effect of changing the pressure correction method.

Spatial convergence rates of the velocity and pressure can be seen in figures 2(a) and 2(b). The velocity converges with the expected order while the pressure super-converges due to the regular mesh. It is normal to observe super-convergence on highly regular meshes (Guillén-gonzález and Tierra, 2012). The divergence of the velocity field is shown in figure 2(c). The divergence has been projected to a piecewise constant (DG0) function space before computing the $L_{2}$ norm. The IPCS-D method gives a velocity field with approximately ten orders of magnitude higher divergence than IPCS-A and SIMPLE. The divergence is also computed in the space of the velocity, piecewise quadratics (DG2), and this is shown in figure 2(d). The same behaviour is shown here, though the error is higher.

Studying temporal convergence for the selected Taylor-Green flow requires a very fine spatial discretisation since the error in the time stepping routine is much smaller than the spatial error. The results obtained with a very fine grid where $N=200$ is shown in figure 3. The time step is gradually refined with $\Delta t=\{2.00,1.00,0.50,0.25\}$ and the numerical solution is compared with the analytical at $t=6.0$. The number of pressure correction iterations per time step is $200$ and the number of Krylov solver iterations is also set to $200$. This is a necessary increase to achieve the theoretical convergence rate for the velocity. The pressure will converge at a rate slightly above $2.0$ with fewer iterations, but even with the fine spatial discretisation, the convergence rate of the velocity is slightly below $2.0$ and the Krylov solver is never fully converged even after $200$ GMRES iterations per pressure iteration. Between refinements $\Delta t=0.5$ and $\Delta t=0.25$ the observed convergence rate of the velocity is $1.97$.

7.2 Ethier-Steinman 3D flow

An analytical solution to the 3D incompressible Navier-Stokes equations is given by Ethier and Steinman (1994, equation 15) as

[TABLE]

A constant correction for the pressure, $\bar{p}$, is included in equation 38 to ensure that the average pressure is zero. This is enforced in the linear solver to remove the pressure kernel from the system. On a domain $(x,y,z)\in[-1,1]\times[-1,1]\times[-1,1]$ this gives

[TABLE]

Dirichlet boundary conditions are applied to the velocity. The domain is divided into a regular mesh with $N\times N\times N$ cubes which are further subdivided into six tetrahedra in each cube. The time step is $\Delta t=0.001$ and the numerical solution is compared to the analytical at time $t=0.1$.

The spatial convergence of the methods in Ethier-Steinman flow can be seen in figure 4. The results for the IPCS methods are presented with 100 pressure correction iterations per time step, but the results differ only in the third significant digit from the results obtained with only 5 iterations. The SIMPLE method has been run with $\alpha_{u}=\alpha_{p}=0.5$ to get stable iterations at the finer grid resolutions where using $\alpha_{u}=0.7$ and $\alpha_{p}=1.0$ would blow up. The SIMPLE results are produced using 500 pressure correction iterations per time step. This gives an optimal convergence rate for the pressure at all resolutions and for the velocity at all but the finest resolution. With 100 pressure correction iterations only two points would fall on the line in figure 4(a). The divergence error norms for the algebraic methods are again close to machine precision. The IPCS-D divergence error is much higher, but getting smaller with increasing mesh resolution, as expected.

7.3 Efficiency

Producing a fair wall clock running time comparison of the pressure correction methods is more tricky than reporting convergence results. For the convergence rates it only matters that the Krylov iterations inside each pressure correction iteration converge in the end, not how many iterations are needed for convergence. The total running time is totally dominated by the time spent in the momentum prediction and pressure correction steps, and except for matrix assembly, which takes around $1$% of the wall clock time, this is all time spent in the Krylov solvers. Tweaking the preconditioner parameters used for the Krylov solvers in each method may improve one of the methods more than the other, causing this comparison—which uses the same preconditioner parameters for all three methods—to not be very relevant. With that caveat in mind, the time spent in the two steps can be seen in figure 5.

The total time per pressure correction iteration is less in the SIMPLE method than in the IPCS methods, which is to be expected due to the under-relaxation procedure causing the unknown vectors to be closer to the initial guesses. The IPCS-D method spends more time in the pressure correction than the IPCS-A method. The reason may be that the penalty parameter, which ensures coercivity, also makes the resulting matrix more and more ill-conditioned as the penalty parameter increases, although the parameter we have used, equation 12, should be close to optimal. The momentum prediction step takes significantly more time in IPCS-A than in IPCS-D. The reason may be related to the skew symmetric convective term, equation 13, making the IPCS-D ${\mathbf{{A}}}$ matrix more diagonal dominant. Unfavourable right hand sides may also require more iterations.

8 Conclusions

For both the 2D (Taylor-Green, figure 2) and 3D (Ethier-Steinman, figure 4) numerical experiments, the convergence results show that the algebraic pressure splitting methods are exactly divergence free while the differential IPCS-D method is not. This is in line with the expected results from the discussion in section 6. For both 2D and 3D tests the IPCS methods converge perfectly with 5 pressure correction iterations per time step, while the SIMPLE method requires between 10 and 100 times more iterations to produce converged results. The SIMPLE method spends slightly less time per iteration (figure 5), but is not fast enough to outweigh the vastly increased number of iterations required to obtain converged results. For mass conserving pressure correction iterations the IPCS-A method is most efficient, and hence the recommended option among the three tested methods.

Acknowledgements

The simulations were performed on resources provided by UNINETT Sigma2, the National Infrastructure for High Performance Computing and Data Storage in Norway.