Continuous, Semi-discrete, and Fully Discretized Navier-Stokes Equations

TL;DR

This paper analyzes the discretization of Navier-Stokes equations, focusing on the numerical difficulties caused by system strangeness and how common time-stepping schemes address these issues.

Contribution

It provides a detailed analysis of the semi-discrete and fully discrete Navier-Stokes equations, highlighting the role of strangeness and the effectiveness of standard time-integration schemes.

Findings

Semi-discrete Navier-Stokes equations form a nonlinear DAE with strangeness.

Fully discrete schemes mitigate strangeness by removing it through time-stepping.

Numerical examples confirm the theoretical analysis.

Abstract

The Navier--Stokes equations are commonly used to model and to simulate flow phenomena. We introduce the basic equations and discuss the standard methods for the spatial and temporal discretization. We analyse the semi-discrete equations -- a semi-explicit nonlinear DAE -- in terms of the strangeness index and quantify the numerical difficulties in the fully discrete schemes, that are induced by the strangeness of the system. By analyzing the Kronecker index of the difference-algebraic equations, that represent commonly and successfully used time stepping schemes for the Navier--Stokes equations, we show that those time-integration schemes factually remove the strangeness. The theoretical considerations are backed and illustrated by numerical examples.

Tables

Table 1. Table 3.1. Selective overview of stable finite element schemes for flow equations.

discretization scheme	discrete function spaces	inf-sup stability
Taylor-Hood	$V_h={\left[{𝒫}_k\left(𝒯\right)\cap 𝒱\right]}^d,k\ge 2$	[TH73, Ver84]
(continuous pressure)	$Q_h=\left[{𝒫}_{k-1}\left(𝒯\right)\cap 𝒱\right]/ℝ$
Crouzeix-Raviart	$V_h={\left[{\mathrm{CR}}_0\left(𝒯\right)\right]}^d$	[CR73, BM11]
(discontinuous velocity)	$Q_h={𝒫}_0\left(𝒯\right)/ℝ$
Bernardi-Raugel	$V_h={\left[{𝒫}_1\left(𝒯\right)\cap 𝒱\right]}^d\oplus ℬ$	[BR85]
(edge-bubble functions)	$Q_h={𝒫}_0\left(𝒯\right)/ℝ$
Rannacher-Turek	on quadrilateral meshes	[RT92]

Equations

$$\mathfrak{v}:\mathcal{I}\times \Omega \to {\mathbb{R}}^3,\mathfrak{p}:\mathcal{I}\times \Omega \to \mathbb{R},\text{and}\rho :\mathcal{I}\times \Omega \to \mathbb{R},$$

$$\frac{\text{d}}{\text{d}t}{\int }_{W}f\text{ d}V={\int }_W\frac{\partial f}{\partial t}+\mathrm{div}(f\mathfrak{v})\text{ d}V.$$

$$\frac{\text{d}}{\text{d}t}{\int }_W\text{ d}V=0\text{if, and only if,}{\int }_W\mathrm{div}\mathfrak{v}\text{ d}V=0.$$

$$\frac{\text{d}}{\text{d}t}{\int }_W\rho \text{ d}V=0\text{if, and only if,}{\int }_W\frac{\partial \rho }{\partial t}+\mathrm{div}(\rho \mathfrak{v})\text{ d}V=0.$$

$$\frac{\text{d}}{\text{d}t}{\int }_W\rho \mathfrak{v}\text{ d}V={\int }_W\rho g\text{ d}V+{\int }_{\partial W}\sigma n\text{ d}S,$$

$${\int }_{\partial W}\sigma n\text{ d}S={\int }_W\mathrm{div}\sigma \text{ d}V,$$

$$\sigma =-\mathfrak{p}I+\tau$$

$$\tau=\mu\bigl{[}\nabla\mathfrak{v}+\nabla\mathfrak{v}^{\mathsf{T}}-\frac{2}{3}(\operatorname{div}\mathfrak{v})I\bigr{]}+\zeta(\operatorname{div}\mathfrak{v})I.$$

$$\frac{\text{d}}{\text{d}t}{\int }_W\rho {\mathfrak{v}}_i\text{ d}V={\int }_W\frac{\partial \rho {\mathfrak{v}}_i}{\partial t}+\mathrm{div}(\rho {\mathfrak{v}}_i\mathfrak{v})\text{ d}V={\int }_W\frac{\partial \rho {\mathfrak{v}}_i}{\partial t}+\rho {\mathfrak{v}}_i\mathrm{div}\mathfrak{v}+\nabla (\rho {\mathfrak{v}}_i)\cdot \mathfrak{v}\text{ d}V$$

$$\frac{\text{d}}{\text{d}t}{\int }_W\rho \mathfrak{v}\text{ d}V={\int }_W\frac{\partial \rho \mathfrak{v}}{\partial t}+\rho \mathfrak{v}\mathrm{div}\mathfrak{v}+(\mathfrak{v}\cdot \nabla )(\rho \mathfrak{v})\text{ d}V.$$

$$\frac{\partial\rho\mathfrak{v}}{\partial t}+\rho\mathfrak{v}\operatorname{div}\mathfrak{v}+(\mathfrak{v}\cdot\nabla)(\rho\mathfrak{v})=\rho g-\nabla\mathfrak{p}+\operatorname{div}\bigl{(}\mu\bigl{[}\nabla\mathfrak{v}+\nabla\mathfrak{v}^{\mathsf{T}}-\frac{2}{3}(\operatorname{div}\mathfrak{v})I\bigr{]}\bigr{)}.$$

$$\displaystyle\rho^{*}\bigl{(}\frac{\partial\mathfrak{v}}{\partial t}+(\mathfrak{v}\cdot\nabla)\mathfrak{v}\bigr{)}+\nabla\mathfrak{p}-\mu^{*}\Delta\mathfrak{v}$$

$$\mathrm{div}\mathfrak{v}$$

$${\mathfrak{v}}^{\prime }=\frac{\mathfrak{v}}{V},{\mathfrak{p}}^{\prime }=\frac{\mathfrak{p}}{{\rho }^{\ast }{U}^2},x^{\prime }=\frac{x}{L},\text{and}{t}^{\prime }=\frac{tV}{L}$$

$$\frac{\partial }{\partial t^{\prime }}{\mathfrak{v}}^{\prime }+({\mathfrak{v}}^{\prime }\cdot \nabla ){\mathfrak{v}}^{\prime }-\frac{1}{\text{Re}}\Delta {\mathfrak{v}}^{\prime }+\nabla {\mathfrak{p}}^{\prime }$$

$$\mathrm{div}{\mathfrak{v}}^{\prime }$$

$$v=g\text{on }\Gamma$$

$$\mathfrak{v}={\mathfrak{v}}_i\text{on }{\Gamma }_{\text{i}},$$

$$\sigma n=0\text{on }{\Gamma }_{\text{o}},$$

$$\frac{1}{\text{Re}}\frac{\partial \mathfrak{v}}{\partial n}-\mathfrak{p}n=0\text{on }{\Gamma }_{\text{o}},$$

$$\frac{\partial \mathfrak{v}}{\partial n}=0\text{and}\frac{\partial \mathfrak{p}}{\partial n}=0\text{on }{\Gamma }_{\text{o}}.$$

$$\mathcal{V}:=[H_0^1(\Omega ){]}^n,\mathcal{H}:=[L^2(\Omega ){]}^n,\text{and }\mathcal{Q}:=L^2(\Omega )/\mathbb{R}.$$

$$\dot{\mathfrak{v}}(t)+\mathcal{K}(\mathfrak{v}(t))-{\mathcal{B}}^{\prime }\mathfrak{p}(t)$$

$$\mathcal{B}\mathfrak{v}(t)$$

$$\mathfrak{v}(0)$$

$$⟨\mathcal{K}(\mathfrak{v}),w⟩={\int }_{\Omega }(\mathfrak{v}\cdot \nabla )\mathfrak{v}\cdot w\text{d}x+\frac{1}{\text{Re}}{\int }_{\Omega }\nabla \mathfrak{v}\cdot \nabla w\text{d}x$$

$$⟨\mathcal{B}\mathfrak{v},q⟩={\int }_{\Omega }(\mathrm{div}\mathfrak{v})q\text{d}x=⟨\mathfrak{v},{\mathcal{B}}^{\prime }q⟩,$$

$$M\dot{v}+K(v)-B^{T}p$$

$$Bv$$

$$M:=[m_{jk}]\in {\mathbb{R}}^{n\times n},m_{jk}:={\int }_{\Omega }{\phi }_j\cdot {\phi }_k\text{d}x.$$

$$K_j(v):={\int }_{\Omega }(\tilde{v}\cdot \nabla )\tilde{v}\cdot {\phi }_j\text{d}x+\frac{1}{\text{Re}}{\int }_{\Omega }\nabla \tilde{v}\cdot \nabla {\phi }_j\text{d}x.$$

$$B:=[b_{ij}]\in {\mathbb{R}}^{m\times n},b_{ij}={\int }_{\Omega }{\psi }_i\mathrm{div}{\phi }_j\text{d}x.$$

$$\dim V_h=n,\dim Q_h=m<n.$$

Full text

**Key words. Navier-Stokes equations, DAEs, strangeness index, difference-algebraic equations

AMS subject classifications. 65L80, 65M12, 35Q30**

1. Introduction

The Navier–Stokes equations (NSE) are a system of nonlinear partial-differential equations that have been commonly used to model fluid flows for more than a century. The NSE are believed to describe all kinds of incompressible flows sufficiently well as long as the setup supports the hypothesis that the fluid is a continuum. Indeed, comparisons of numerical simulations with experiments show arbitrarily good agreement of the model with the observations over a long range from slowly moving flows in small geometries like a pipe up to highly turbulent flows over large spatial extensions like the flow around an airplane or even weather phenomena. Nevertheless, the mere equations and the extent of their applicability have not been fully deciphered by now and a substantial progress in this respect will be eligible for a Clay price111see http://www.claymath.org/millennium-problems/navier–stokes-equation.

Under the assumption of continuity of the observed quantities, the NSE can be derived from fundamental laws of physics; see [LL87] and [CM93]. These considerations are well backed for a macroscopic viewpoint, from which a fluid like water appears as a continuum. On a microscopic level, where discrete molecular structures define the states, the NSE cannot capture the physics right, as it is well-known, e.g., for capillary flows.

On the molecular level, fluids are better described by the Boltzmann equations, which model molecular interactions. This fact seems undisputed the more that the NSE can also be interpreted and derived through a limiting process of the Boltzmann equations in the sense of averaging the microscopic quantities for a macroscopic description [Sai09].

As a mathematical object the NSE have ever been subject to fundamental investigations and led to its own research field and its own subject definitions in the MSC classification scheme222see http://www.ams.org/mathscinet/msc/msc.html?t=35Q30&btn=Current and related . The research on the NSE has focussed on the analysis of the equations and their numerical approximation. Early results on the existence of solutions are due to Leray [Ler33] (weak solutions) and Fujita&Kato [FK64] (smooth solutions). The first textbooks on the functional and on the numerical analysis were written by Ladyzhenskaya [Lad69] and Temam [Tem77], respectively.

For the numerical analysis of the spatial discretization, one may distinguish two lines of development. The mathematical line focusses on Galerkin methods in the realm of variational formulations whereas the engineering orientated line has been advancing finite volume methods (FVM) as they appear well-suited for simulations. On the side of Galerkin methods, and in particular finite element methods (FEM), there have been many efforts in designing stable elements like the famous Taylor–Hood elements [TH73] as well as for general convergence results; see the textbooks [GR86, Lad69, Pir89] for the numerical analysis and [Tur99] for an application oriented overview. On the side of FVM that are the method of choice in most general purpose flow solvers, there have been general developments in view of discretizing conservation laws [LeV92] and particular progress in view of stable approximation of fluid flow [FP02].

A numerical analysis of approximations to the time-dependent NSE with FEM semi-discretizations has been carried out by Heywood&Rannacher [HR90]; see also the textbook [GS00], that covers implementation issues. For time marching schemes for FVM formulations we refer to [FP02]. Strategies for the iterative solution of the arising linear systems can be found in [ESW05] (FEM) or [FP02] (FVM).

In this work we revisit the NSE from a differential-algebraic equations (DAE) perspective. This includes the modelling where the incompressibility is treated as an algebraic constraint in an abstract space and the spatial semi-discretization, that has to be handled with care to respect the incompressibility constraint and to lead to a well-posed classical DAE. It also includes the temporal discretization, in which the DAE properties of the NSE become evident and (hopefully not) problematic. Such a pure DAE perspective has been taken on by Weickert [Wei97] who analysed finite difference approximations and certain time-stepping schemes for the NSE, by Emmrich&Mehrmann [EM13] who provided conditions and solution representations for linearized NSE in abstract spaces in line with linear time-invariant DAEs in finite-dimensional state-spaces. In [Hei14] the nonlinear NSE and its Galerkin approximations have been analysed in view of consistency of reformulations of semi-discrete DAE approximations with the infinite-dimensional model.

The paper is organized as follows. In Section 2 we derive the NSE from first principles and formulate the weak form as an operator DAE. The direct connection to DAEs is then made in Section 3, in which we report on several spatial discretization schemes and that commonly used FEM schemes lead to systems of strangeness index one. In Section 4 we analyse the time approximation schemes in terms of the index of the resulting difference-algebraic equations ($\Delta$AE). We show that the straight-forward temporal discretization leads to a scheme of higher index than that of well-established time-stepping schemes for incompressible flows. We further confirm, that schemes with a $\Delta$AE of lower index can also be obtained from a standard time-discretization applied to a reformulation of the DAE with lower index. In the numerical examples in Section 5, we confirm the superiority of the lower index $\Delta$AE approximations of the NSE over the straight-forward time-discretization. As a benchmark for time-integration schemes for the considered class of nonlinear semi-explicit DAEs of strangeness index 1, we provide reference trajectories and the system coefficients and nonlinear inhomogeneities for direct realization in Python. We conclude this paper with summarizing remarks in Section 6.

2. Continuous model

2.1. Derivation of the Navier–Stokes equations

The Navier–Stokes equations (NSE) provide a model of a flow as a continuum. The basic assumption for their derivation and, thus, the validity of the model is that the flow under consideration forms a continuous entity of flow particles in a spatial domain $\Omega\subset\mathbb{R}^{3}$ and a time interval $\mathcal{I}$ such that the functions

[TABLE]

describing the velocity, the pressure, and the density as measured at the position $x\in\Omega$ at time $t\in\mathcal{I}$ are continuous functions. Here, continuous means that we can apply differential calculus in order to derive basic partial differential equations. Later, when we derive the weak formulation of the NSE, the needed continuity will be specified further.

As mentioned in [CM93, p. 2], these assumptions lead to a model that is believed provide accurate descriptions of common macroscopic flow phenomena. In, e.g., setups of small geometric scales like in capillary flows or under vacuum-like conditions [KBN05], the discrete microscopic molecular structure of the fluid that constitutes the flow has to be taken into account.

The basic assumption of continuity of the matter allows for the consideration of a possibly infinitesimal small control volume $W\subset\Omega$, an open bounded domain in $\mathbb{R}^{3}$ that contains a given agglomerate of fluid particles in the considered flow. Continuity also implies that a fixed $W$ is deformed and convected by the flow but always consists of the same fluid particles.

Under this continuity assumption, one can call on the Reynolds Transport Theorem, that relates the temporal change of an integral quantity over $W$ to the convection velocity $\mathfrak{v}$.

Theorem 2.1** (see [Rey03], Eq.(16) and [CM93], p.10).**

Let $W\colon I\to\mathbb{R}^{3}$ describe a smoothly moving control volume and let $f\colon t\times W(t)\mapsto\mathbb{R}$ be a sufficiently smooth function, then

[TABLE]

Incompressibility and mass conservation

The Reynolds Transport Theorem can be used to show that a flow is incompressible, which means that the volume of any agglomerate $W$ is constant over time, if and only if the velocity field $v$ is divergence free. In fact, if (2.2) applies, then one has that

[TABLE]

With a well-defined density function $\rho$, the mass of a control volume is defined as the integral over the (mass) density $\rho$ and, with the assumption that in the flow there are neither mass sinks nor mass sources, an application of (2.2) gives

[TABLE]

Remark 2.2*.*

The conservation of mass and the incompressibility of a flow are closely related but only equivalent in the case that the density $\rho$ is constant in space and time. Flow models that assume incompressibility and varying density functions $\rho$ are applied, e.g., in oceanography [Tar06].

Remark 2.3*.*

If the volume of a fluid parcel changes over time, the flow is called compressible. In this case, the mass density $\rho$ is modelled as an unknown function and related to the pressure $p$ through constitutive or so-called state equations, like the ideal gas low; see [FKP16].

Balance of momentum

Under the continuum assumption, the momentum of a fluid agglomerate $W$ can be expressed as the integral of the mass density times velocity and the temporal change equated with volume and surface forces on $W$:

[TABLE]

where $g\colon\mathcal{I}\times\Omega\to\mathbb{R}^{3}$ is the density of a body force, and where $\sigma\colon\mathcal{I}\times\Omega\to\mathbb{R}^{3,3}$ is a tensor such that $\sigma n$, where $n$ is the normal field on the boundary $\partial W$ of $W$, represents the density of the forces acting on the surface. Note that the expression in (2.5) is vector valued and that integration and differentiation is performed componentwise.

Applying the Divergence Theorem componentwise, one can write the term with the surface forces as a volume integral

[TABLE]

with the divergence of a tensor defined accordingly.

So far, all assumptions have based on first principles. For the mathematical modelling of the tensor $\sigma$, however, ad hoc assumptions and heuristics are employed. First of all, it is assumed that $\sigma$ can be written as

[TABLE]

where $\mathfrak{p}$ – the pressure – is a smooth scalar function and $\tau\colon\mathcal{I}\times\Omega\to\mathbb{R}^{3,3}$ is the tensor of shear stresses. It is assumed that $\tau(t,x)$ is symmetric and invariant under rigid body rotation. Furthermore, $\tau$ is a linear function of the velocity gradient \nabla\mathfrak{v}:=\bigl{[}\frac{\partial\mathfrak{v}_{j}}{\partial x_{i}}\bigr{]}_{i,j=1,2,3}.

Remark 2.4*.*

The separation of $\sigma$ into a pressure and shear stress part is motivated by requirement to recover the equations of hydrostatic or ideal flows, where there are no shear stresses because of the absence either of motion or of viscosity. The symmetry of the tensor $\tau$ can be derived from the requirement that the angular momentum of a flow agglomerate is conserved. The invariance of $\tau$ with respect to rigid body rotations is derived from the general assumption that the considered fluid is isotropic, i.e., the physical properties of the fluid are the same in all spatial directions. Finally, the assumption that $\tau$ depends linearly on $\nabla u$ is well grounded for many fluids of interest. In fact, the validity of this assumption classifies a fluid as Newtonian fluid.

Under the given assumptions, $\tau$ is defined by the velocity field and two further parameters. It is commonly written as

[TABLE]

The parameter $\mu$ is the (first coefficient of the) viscosity and has been experimentally determined and tabulated for many gases and fluids. The parameter $\zeta$ is called the bulk viscosity and, in line with the Stokes Hypothesis, often set to zero.

As for the left hand side in (2.5), one proceeds as follows. Let $\mathfrak{v}_{i}$ denote the $i$-th component of $\mathfrak{v}$, $i=1,2,3$. Then, an application of (2.2) gives that

[TABLE]

where basic vector calculus has been applied. If one considers $\nabla$ as the formal column vector of the three space derivatives, relation (2.9) for all components of $v$ can be written in compact form as

[TABLE]

The Navier–Stokes equations

In the preceding derivations, integral quantities over a flow agglomerate $W$ were considered. If the underlying continuity assumption includes that all relations hold on arbitrary (small) control volumes $W$, instead of equating the integrals, one can equate the integrands pointwise in space, which leads to partial-differential equations.

Thus, putting together all assumptions and derivations, the balance of momentum (2.5) for an isotropic Newtonian fluid under the Stokes Hypothesis defines the (NSE) as

[TABLE]

In the case that the flow is incompressible with a constant density $\rho\equiv\rho^{*}$ and a constant viscosity $\mu\equiv\mu^{*}$, the combination of (2.4) and (2.11) results in the system

[TABLE]

Remark 2.5*.*

Note that in the derivation, the pressure $\mathfrak{p}$ is not a variable but a function which is assumed to be known and to describe the normal forces on a fluid element, cf. (2.7) and [CM93, Ch. 1.1.ii]. Thus, the NSE (2.11) has only the velocity $\mathfrak{v}$ as an unknown and the divergence free formulation (2.12) can be seen as an abstract ODE for $\mathfrak{v}$ with an invariant [HLW06]. Nonetheless, since it turns out that the divergence constraint defines the function $\mathfrak{p}$, system (2.12) is commonly considered as an abstract differential-algebraic equation with $\mathfrak{v}$ and $\mathfrak{p}$ as unknowns.

For the numerical treatment and for similarity considerations, one relates all dependent and independent variables to a characteristic length $L$ and characteristic velocity $V$ via

[TABLE]

With this and with $f^{\prime}:=\frac{U^{2}}{L}g$, equation (2.12) can be rewritten in dimensionless form

[TABLE]

Thus, the system is completely parameterized by only one parameter $\text{Re}:=\frac{UL\rho^{*}}{\mu^{*}}$, the Reynolds number. In what follows, we will always consider the dimensionless NSE (2.14) but drop the dashes of the dimensionless variables.

Boundary conditions

For the considered domain $\Omega$, let $\Gamma$ denote the boundary in the abstract and in the physical sense. If the domain is bounded by nonpermeable walls, then the no-slip condition

[TABLE]

applies, where $g$ is the velocity of the wall. The no-slip conditions align well with experiments and macroscopic considerations; see [Lay08, Ch. 5.3] for references but also for examples where no-slip conditions seem insufficient to describe the flow at walls.

If the domain is not fully bounded by walls, typically because the computational domain needs to be bounded whereas the physical domain of the flow is unbounded, artificial boundary conditions are needed. We assume that the boundary is composed of a part $\Gamma_{\text{w}}$ that is associated with a wall, a part $\Gamma_{\text{i}}$ where an (incoming) velocity profile is known, i.e.,

[TABLE]

and a part $\Gamma_{\text{o}}$ that is often associated with an outflow.

There is no immediate physical insight into what should be the mathematical conditions at an artificial boundary condition $\Gamma_{\text{o}}$ that models the part where the fluid leaves the domain. Accordingly, outflow boundary conditions are motivated as being useful for the implementation of downstream boundary conditions (see, e.g., [Glo03, Pir89]) and equipped with the advice to put the outlet sufficiently far away from the region of interest; see [FP02, Ch.8.10.2].

The most common outflow boundary conditions are the no stress conditions:

[TABLE]

where $\sigma$ is the stress tensor as defined in its general form in (2.7), the do-nothing conditions:

[TABLE]

that are formulated for the nondimensional Navier–Stokes equations (2.14), or the no gradient conditions:

[TABLE]

that are also applied at symmetry planes.

Remark 2.6*.*

The do-nothing conditions have been extended; see [BM14], to the case of backflow, i.e., when some, possibly spurious, inflow occurs at the boundary $\Gamma_{\text{o}}$.

2.2. Formulation as operator DAE

In this subsection we provide yet another formulation of the dimensionless NSE (2.14), namely in the weak form as an operator DAE. This is a DAE in an abstract setting, where the solution is an element of a Sobolev space instead of a vector. We consider here homogeneous Dirichlet boundary conditions. More general boundary conditions with in- and outflow may modeled in a similar way. We define the spaces

[TABLE]

By $\mathcal{V}^{\prime}$ we denote the dual space of $\mathcal{V}$. Note that the spaces $\mathcal{V}$, $\mathcal{H}$, $\mathcal{V}^{\prime}$ form a Gelfand or evolution triple [Zei90, Ch. 23.4]. Furthermore, we define $\mathcal{W}(0,T)$ as the space of functions $u\in L^{2}(0,T;\mathcal{V})$, which contain a weak time derivative $\dot{u}\in L^{2}(0,T;\mathcal{V}^{\prime})$. Well-known embedding results then imply $u\in C([0,T];\mathcal{H})$, cf. [Rou05, Lem. 7.3].

We now consider the weak formulation of (2.14) in operator form. This means that for given right-hand sides $\mathcal{F}\in L^{2}(0,T;\mathcal{V}^{\prime})$, $\mathcal{G}\in L^{2}(0,T;\mathcal{Q}^{\prime})$ and an initial condition $a\in\mathcal{H}$, we seek for a pair $(\mathfrak{v},\mathfrak{p})\in\mathcal{W}(0,T)\times L^{2}(0,T;\mathcal{Q})$ satisfying

[TABLE]

a.e. on $(0,T)$. The derivative of $\mathfrak{v}$ should be understood in the weak sense. The operators $\mathcal{K}\colon\mathcal{V}\to\mathcal{V}^{\prime}$ and $\mathcal{B}\colon\mathcal{V}\to\mathcal{Q}^{\prime}$ are defined via

[TABLE]

and

[TABLE]

respectively, for a given $\mathfrak{v}\in\mathcal{V}$ and for all test functions $w\in\mathcal{V}$ and $q\in\mathcal{Q}$. We emphasize that system (2.18) not only covers the NSE but also more general flow equations, since we have introduced an inhomogeneity $\mathcal{G}$.

For results on the existence solutions to the weak formulation of the NSE, we refer to [Tar06]. For a compact summary that also considers nonzero $\mathcal{G}$ in (2.18b) see [AH15]. The differential-algebraic structure of (2.18) and the possible decoupling of differential and algebraic parts and variables has been discussed in [Hei14].

3. Semi-discrete equations

In this section, we consider the DAE, which results from a spatial discretization of system (2.14) or (2.18). Using finite elements, finite differences, or finite volumes, we obtain a DAE of the form

[TABLE]

Here, $v$ and $p$ denote the finite-dimensional approximations of $\mathfrak{v}$ and $\mathfrak{p}$, respectively. In addition, we assume a consistent initial condition of the form $v(0)=v_{0}$, i.e., we demand $Bv_{0}=g(0)$.

There are many reasons to include a right-hand side $g$ in equation (3.1b) rather than a zero as it seems naturally for incompressible flows. Firstly, a nonzero $g$ in the continuous continuity equations (2.14b) may appear also in generalizations of the NSE model to fluid-structure interactions (see, e.g., [Ray10]) or in optimal control setups; see, e.g., [Hin00]. Secondly, a nonzero $g$ maybe a numerical artifact from the semidiscretization of flows with nonzero Dirichlet boundary conditions as in the example we will provide below. Finally, in view of analysing the DAE, this inclusion of a nonzero $g$ leads to a better understanding as one can track where the derivatives of the right-hand sides appear within the solution. This is of importance also in the case $g=0$, since inexact solves lead to errors in the (3.1b) that act like a nonsmooth inhomogeneity.

3.1. Spatial discretization by finite elements

We consider a discretization by finite elements, i.e., we construct finite-dimensional (sub)spaces $V_{h}$ and $Q_{h}$ of $\mathcal{V}$ and $\mathcal{Q}$, respectively, based on a (shape) regular triangulation $\mathcal{T}$ of the polygonal Lipschitz domain $\Omega$, cf. [Cia78]. Given a basis $\{\varphi_{1},\dots,\varphi_{n}\}$ of $V_{h}$, we can identify the finite-dimensional approximation of the velocity $\mathfrak{v}(t)$ by the coefficient vector $v(t)\in\mathbb{R}^{n}$. The discrete representative of the pressure $\mathfrak{p}$(t) is denoted by $p(t)\in\mathbb{R}^{m}$ and corresponds to a basis $\{\psi_{1},\dots,\psi_{m}\}$ of $Q_{h}$.

With the basis functions of $V_{h}$ we define the symmetric and positive definite mass matrix $M\in\mathbb{R}^{n,n}$ by

[TABLE]

For specific discretization schemes the mass matrix may even equal the identity matrix. The discretization of the nonlinearity of the Navier-Stokes equation, i.e., the discretization of the operator $\mathcal{K}$ in Section 2.2, is denoted by $K\colon\mathbb{R}^{n}\to\mathbb{R}^{n}$. Note that the given model also includes linearizations of the Navier-Stokes equation such as the unsteady Stokes or Oseen equation. In this case, $K$ can be written as a $n\times n$ matrix. In view of the index analysis of system (3.1), however, this is not of importance; see Appendix A. We define for $v\in\mathbb{R}^{n}$ and its representative $\tilde{v}=\sum_{j=1}^{n}v_{j}\varphi_{j}\in V_{h}$,

[TABLE]

Finally, we define the matrix $B$, which corresponds to the divergence operator, i.e.,

[TABLE]

Since $\operatorname{div}$ and $-\nabla$ are dual operators in the continuous setting, in the semi-discrete equations we get a matrix $B$ and (minus) its transpose $B^{T}$. Note that this yields the saddle point structure of system (3.1).

Especially for the stable approximation of the pressure, it is necessary that the chosen finite element spaces are compatible [BF91, Ch. VI.3]. Let $V_{h}$ and $Q_{h}$ denote again finite-dimensional (sub)spaces of $\mathcal{V}$ and $\mathcal{Q}$, respectively, with

[TABLE]

The spaces $V_{h}$ and $Q_{h}$ are compatible if they satisfy a so-called inf-sup or Ladyzhenskaya-Babu${\check{\text{s}}}$ka-Brezzi condition [BF91, Ch. VI.3]. This means that there exists a constant $\beta>0$, independent of the chosen mesh size, such that

[TABLE]

The inf-sup condition, with $\beta$ independent of $h$, is a necessary condition for the convergence of the FEM; see, e.g., [Lad69] or [GR86]. For a fixed spatial discretization, this condition implies that the matrix $B$ resulting from the discretization scheme $V_{h}$, $Q_{h}$ is of full rank. One scheme which is known not to satisfy this condition (in dimension $d=2$) is given by

[TABLE]

Here $\mathcal{P}_{k}(\mathcal{T})$ denotes the space of piecewise polynomials of degree $k$. Thus, the space $V_{h}$ involves continuous finite element functions of polynomial degree one, whereas $Q_{h}$ consists of piecewise constant functions. Furthermore, $V_{h}$ includes homogeneous Dirichlet boundary conditions, since they are part of $\mathcal{V}$. It has been shown that that the velocity space is too small and thus $\beta\not>0$ [BF91, Ex. VI.3.1]. This flaw can be fixed in various ways. One possibility is given by introducing so-called nonconforming finite element spaces. Crouzeix and Raviart [CR73] proposed to use piecewise polynomials of degree one, which are only continuous in the midpoints of the edges of the triangulation. Together with homogeneous Dirichlet boundary conditions the ansatz space is denoted by $\operatorname{CR}_{0}(\mathcal{T})$.

Another strategy is to enrich the space $V_{h}$ by so-called bubble functions, namely $\mathcal{B}$, cf. [BR85]. With this we introduce ansatz functions of polynomial degree two. A small collection of stable finite element schemes is given in Table 3.1. Further stable schemes are addressed in [GR86, Ch. II] as well as in [GS00, Ch. 3].

Remark 3.1*.*

Factoring out $\mathbb{R}$ in the definition of the pressure approximation spaces $\mathcal{Q}_{h}$ is needed to eliminate the kernel of the gradient operator or $B^{T}$, which are the constant functions both in finite and infinite dimensions; see [GR86]. This reflects the fact that in internal flows the pressure is defined up to a constant. In fact, the standard results for inf-sup stability always assume the case of internal flow, i.e. zero Dirichlet conditions everywhere at the boundary. There have been but a few attempts to derive inf-sup conditions that apply to nonzero boundary conditions [GH92].

3.2. Finite volumes and finite differences

In this section we briefly touch the spatial discretization of the incompressible NSE by the methods of finite volumes (FVM) and finite differences (FDM). If applied to (2.14) in a straight forward manner, both approaches lead to a DAE of type (3.1).

Finite volume approximations base on the integral formulation of the conservation laws as it reads for the momentum equation

[TABLE]

cf. (2.5) and (2.7). Balances, like (3.3) in particular, hold for the whole domain of computation. For the discretization, the domain of the flow is subdivided into small volumes often referred to as cells. The velocities and the pressure are assumed to be, say, constant over the cells, and their approximated values are determined by evaluating and equating the volume and surface integrals associated with every cell; see [FP02, Ch. 8.6].

Because of its flexibility in the discretization, because of its variants that provide unconditional stability properties, and since typical models for the effective treatment of turbulence are formulated as conservation laws too, the FVM is the method of choice in most general purpose solvers. As for the discussion in the DAE context we note that, typically, the DAE (3.1) is never assembled but rather decoupled during the time discretization; see [FP02, Ch. 7] and Section 4.3.2 below.

Approximations of the NSE via FDM are not widely used because of the known flaws of finite difference approximations like high regularity requirements and confinement to regular grids. Nonetheless, FDM discretizations have been successfully used in flow discretizations and are probably still in use in certain specified, say, single-purpose codes. Also, FDM are often well suited for teaching the fundamentals of flow simulation of laminar and turbulent setups; see [GDN97].

3.3. Index of the DAE

As argued above, in (3.1) the coefficient matrix $B$ may be rank-deficient and, thus, make the differentiation index not well-defined [Wei96]. However, the strangeness index [KM94, KM01, KM06] that applies to over- and underdetermined DAE systems and, thus, also can be determined in the case of a rank-deficient coefficient matrix $B$. A rough index analysis has been realized in [Wei97] with the result that under general and reasonable assumptions, system (3.1) has strangeness index 1. See also Appendix A, where the strangeness index has been determined in a rigorous way. This matches the, say, observations in [HW96, Ch. VII.1] that the system is of (differentiation) index 2 if the matrix $B$ is of full rank.

In the sequel we always assume that the discretization scheme is chosen in such a way that the matrix $B$ is of full rank. Note that this can always be realized, choosing linearly independent finite element basis functions with respect to the space $\mathcal{Q}$, i.e., keeping in mind that constant functions are in the same equivalence class as the zero-function. If the ansatz functions form a partition of unity, then this means nothing else than eliminating one of these ansatz functions. Numerical schemes used in practice usually satisfy this condition as the discrete ansatz space for the pressure $Q_{h}$ is chosen appropriately.

4. Fully discrete approximation schemes

The final step of a numerical approximation of the infinite-dimensional NSE (2.14) is the numerical time integration of the semi-discrete approximation (3.1). For this, one discretizes the time interval $(0,T)$ via the grid

[TABLE]

and computes a sequence $(v^{k},p^{k})_{k\in\mathbb{N}}$ of values that are supposed to approximate the solution of (3.1) at the discrete time instances, i.e., $v^{k}\approx v(t^{k})$ and $p^{k}\approx p(t^{k})$.

For example, a single-step time integration scheme applied to an ODE $\dot{x}=\phi(t,x)$ leads to an approximating sequence defined through

[TABLE]

where $\tau^{k}:=t^{k+1}-t^{k}$ and $\Phi$ is the increment function of the single-step scheme.

4.1. Index and causality of discrete systems

For the analysis, we will restrict our considerations to a linear time-invariant setup. This is no restriction, since in the particular semi-linear semi-explicit form, the DAE structure is not affected by the nonlinearity the more that in the presented methods the momentum equation is typically discretized with an explicit scheme. Also we assume that the time grid is equidistantly spaced and of size $\tau=t^{k+1}-t^{k}$.

We will cast the fully-discrete schemes into the standard form

[TABLE]

with coefficient matrices $\mathcal{E}$ and $\mathcal{A}$ and $x^{k}$ containing all variables like $x^{k}=[v^{k};p^{k}]$. We call such a fully discrete approximation scheme (4.2) a difference-algebraic equation ($\Delta$AE).

Remark 4.1*.*

This recast of the scheme, which will basically amount to a shift of the time indices for some variables, is needed to avoid ambiguities. In fact, for continuous time DAEs with delays, one can deliberately apply the operations of time-shift and differentiation to produce equivalent formulations of different indices [Ha15].

In what follows, we argue that the matrix pair \bigl{(}\mathcal{E},\mathcal{A}\bigr{)} can be assumed to be regular (cf. [KM06, Def. 2.5]) such that it can be brought into a particular Kronecker form

[TABLE]

where $J$ is a matrix in Jordan form and $N$ is a nilpotent matrix, cf. [KM06, Thm. 2.7].

Definition 4.2**.**

The index of nilpotency – that integer $\nu$ for which $N^{\nu}=0$ while $N^{\nu-1}\neq 0$ – is called the index of the matrix pair \bigl{(}\mathcal{E},\mathcal{A}\bigr{)}; see [KM06, Def. 2.9], or the Kronecker index of the $\Delta$AE (4.2), cf. [Rei06, p. 39].

The excursion to the discrete Kronecker index is necessary, since the common time stepping schemes considered in Section 4.3 below can not be interpreted as a time-discretization of time-continuous DAE. Equivalently, we can can classify the discrete schemes using the notion of causality that attributes systems with states depending only on the past or current inputs.

Definition 4.3** ([Dai89, Def. 8-1.1]).**

The sequence $x^{k}$, $k=1,2,\dots$, defined through (4.2) is called causal, if $x^{k}$ is determined completely by an initial condition $x^{0}$ and former (and the current) inputs $f^{0},f^{1},\dotsc,f^{k}$.

Conversely, in a noncausal system, the current state $x^{k}$ depends on future inputs like $f^{k+1}$. We will analyse time-stepping schemes for causality, which in the considered case is equivalent to being of Kronecker index $1$ [Dai89, Thm. 8-1.1] and discuss how and why a noncausal system poses difficulties in the numerical approximation.

To introduce the procedure, to fix the notation, and to have a benchmark for further comparisons, we start with analysing an implicit-explicit Euler discretization of a linearized version of (3.1), namely

[TABLE]

This defines the difference equations for $v^{k}$ and $p^{k}$, which approximate the velocity and pressure at the discrete time instances $t^{k}$, $k=1,2,\dots$. Similarly, $f^{k}$ and $g^{k}$ stand for the approximations of $f$ and $g$ at the time instances.

The difference scheme is of the standard form (4.2) with $x^{k}=[v^{k};p^{k}]$, $h^{k}=[f^{k};g^{k}]$, a shift of the index in (4.4b), and

[TABLE]

In Appendix B it is shown that under standard conditions and for $\tau$ sufficiently small, the pair $(\mathcal{E},\mathcal{A})$ in (4.5) is regular and equivalent to

[TABLE]

This means that the difference scheme based on an implicit-explicit Euler discretization is of Kronecker index 2. In fact, a straight forward calculation reveals that in (4.4) the state $p^{k}$ depends on $g^{k+1}$.

Remark 4.4*.*

The index $k$ for the pressure $p$ in (4.4) is consistent with $f^{k}$, as can be directly derived from the case that $A=0$ and $g=0$. A fully implicit scheme, i.e., considering $f^{k+1}$ and $p^{k+1}$ in (4.4a), can not be brought into the standard form (4.2).

4.2. Inherent instabilities of $\Delta$AEs of higher index

In this section, we illustrate a mechanism that leads to a numerical instability and, thus, possibly to divergence of the approximation of a dynamical system through a time discretization, i.e. a $\Delta$AE, with an index greater than $1$. Consider a DAE in Kronecker form

[TABLE]

with $N\neq 0$ and $N^{2}=0$ and it’s time discrete approximation through an Euler scheme:

[TABLE]

The $\Delta$AE (4.6) is of Kronecker index 2, according to Definition 4.2 and as it can be read off after a premultiplication by $(I-N)$, and it has the solution

[TABLE]

as it follows from an adaption of the arguments in [KM06, Lem. 2.8] to the discrete case. From this solution representation one can conclude, that the solution to a $\Delta$AE of higher index that discretizes a DAE may depend on numerical differentiations and that any error in the computation may be amplified by the factor $\tau^{-1}$. Note that for index-1 $\Delta$AEs, where $N=0$, this derivative is not present and that for even higher indices higher (numerical) derivatives will appear in the solution. Also note, that for systems that are not in Kronecker form such as (4.4), these derivations will be realized implicitly; see [AH15].

4.3. Common time-stepping schemes as index-1 $\Delta$AEs

In this subsection, we discuss the different strategies, which are used in practice, to solve the spatially discretized NSE (3.1). In practical applications one typically uses schemes that decouple pressure and velocity computations. Although, as we have argued from a DAE perspective [AH15], this may lead to instabilities, the advantages that

• •

one has to solve two smaller systems rather than one large and

• •

one basically solves Poisson equations and convection-diffusion rather than saddle-point problems

seem to prevail. The second point is important in so far as for the automated solution of the coupled problem there are no generally well performing solvers for the arising linear systems. For example, in the cases of practical relevance namely that of high Reynolds numbers or highly non-normal coefficient matrices, there are no efficient preconditioners for solving the resulting saddle-point systems iteratively. In fact, recent developments in this direction were tested in low to moderate Re-setups and the reported numbers indicate an unfavourable scaling of the iteration numbers with Re; see, e.g., [ESW05, Ch. 8] or [HV16]. Sometimes, in particular when it comes to parallelization, Gauss–Seidel-type methods are used, which iteratively and locally update the variables. As far as convection and diffusion is concerned such a formulation, that the quantities of interests in a cell depend on the local neighbors, is physically well backed.

Another common feature of the schemes used in practice is their explicit approach to the momentum equation which avoids the repeated assembling of Jacobians. The error is then either controlled through a small time-step or through some fixed point iterations.

In this subsection, we consider the schemes Projection as it was described in [Gre90], SIMPLE as in [FP02], and artificial compressibility [FP02]. We will use a description and formulation general enough, to also accommodate numerous variants of the methods.

4.3.1. Projection

The principle of these methods is to solve for an intermediate velocity approximation that does not need to be divergence-free, and project it onto the divergence-free constraint in a second step. As far as the velocity approximation is concerned, Projection methods can be formulated both in infinite and finite dimensions. Since the first work by Chorin [Cho68], a number of variants have been developed mainly proposing different approaches to the approximation of the pressure tackling or circumventing the need of solving a Poisson equation for the pressure, which requires certain regularity assumptions (cf. [GS00, p. 642]). Another problem is the requirement of boundary conditions for the pressure update that do not have a physical motivation and may cause inaccuracies close to the boundary. Nonetheless, these schemes have been extensively studied and certain heuristics ensure satisfactory convergence behaviour; see, e.g, [Gre90, GS00].

As an example for a Projection scheme, we present the variant proposed in [Gre90]:

1. (1)

Solve for intermediate velocity with the old pressure

[TABLE] 2. (2)

Determine the new velocity $v^{k+1}$ as the projection of the intermediate velocity onto $\ker B$ by solving

[TABLE] 3. (3)

Update the pressure via

[TABLE]

Remark 4.5*.*

Instead of solving the saddle-point problem (4.8) as a whole and in order to avoid the division of a numerically computed quantity by $\tau$ in (4.9), one can decouple the system. Then, one solves for $\tilde{\phi}^{k+1}:=\frac{2}{\tau}\phi^{k+1}$ through

[TABLE]

and obtains the updates via

[TABLE]

and

[TABLE]

Remark 4.6*.*

The formula for the update of the pressure is derived from the relation

[TABLE]

cf. [Gre90, p. 595], that holds under certain regularity assumptions. Note also that, if formulated for the space-continuous problem, this Projection 2 algorithm requires a sophisticated treatment of the boundary conditions. The presented variant is a simplification for spatially discretized equations; see [GS00, Ch. 3.16.6c].

As a single system, this projection scheme defines a $\Delta$AE in the form of (4.2) with $x^{k}:=[\tilde{v}^{k};\phi^{k};v^{k};p^{k}]$ and

[TABLE]

which is a matrix pair of Kronecker index 1.

4.3.2. SIMPLE scheme – implicit pressure correction

The SIMPLE scheme and its variants are based on the decomposition

[TABLE]

where $\tilde{v}^{k+1}$ is the tentative velocity computed by means of the old pressure. We present the basic variant, in which the velocity correction $v_{\Delta}^{k+1}$ is discarded when solving for $p_{\Delta}^{k+1}$:

1. (1)

Solve for the intermediate velocity $\tilde{v}^{k+1}$ with the old pressure as in (4.7). 2. (2)

Compute $p_{\Delta}^{k+1}$ through

[TABLE]

as the correction to $p^{k}$ such that the 3. (3)

updates of the velocity and pressure defined through

[TABLE]

jointly fulfill the time discrete momentum and the continuity equation.

Remark 4.7*.*

Step (2) of the presented SIMPLE algorithm computes $p_{\Delta}^{k+1}$ under the temporary assumption that $v_{\Delta}^{k+1}=0$. This is hardly justified and probably a reason for slow convergence, when applied in a fixed-point iteration within fully implicit schemes. Certain variants of the SIMPLE scheme try to approximate $v_{\Delta}^{k+1}$ at this step, e.g., through interpolation; see [FP02, Ch. 7.3.4].

As a single system, the SIMPLE scheme defines a $\Delta$AE in the form of (4.2) with $x^{k}:=[\tilde{v}^{k};p_{\Delta}^{k};v^{k};p^{k}]$ and

[TABLE]

which is a matrix pair of Kronecker index 1; see Appendix B.3.

4.3.3. Artificial compressibility

In this class of methods, the divergence-free constraint (2.14b) is relaxed by adding a scaled time-derivative of the pressure, i.e.,

[TABLE]

The parameter is assumed to satisfy $\beta\gg 1$. The corresponding spatial semi-discretization (3.1) then reads

[TABLE]

with $M_{p}$ denoting the mass matrix of the pressure approximation.

In theory, system (4.16) is an ODE and could be solved by standard time-stepping schemes. In practice, however, the common solution approaches to System (4.16) decouple pressure and velocity computations through introducing auxiliary quantities and, thus, a DAE structure. A possible method, that defines a pressure update similar to the SIMPLE scheme (4.13), can be interpreted and implemented as follows, cf. [FP02, Ch. 7.4.3]. As in the SIMPLE approach a first velocity approximation $\tilde{v}^{k+1}$ is computed via (4.7) on the base of the old pressure value $p^{k}$.

Next, one substracts the contribution of the old pressure,

[TABLE]

and defines the new velocity $v^{k+1}$ through the linear expansion of $\bar{v}^{k+1}$ in terms of the pressure gradient, i.e.,

[TABLE]

Note that one uses (4.17) to determine the Jacobian $\frac{\partial(\bar{v}^{k+1})}{\partial(B^{T}p)}$ and that we have defined $p_{\Delta}^{k+1}=p^{k+1}-p^{k}$. Expression (4.18) is then inserted in a, e.g., implicit Euler discretization of (4.16b), which gives an equation for $p_{\Delta}^{k+1}$. Accordingly, the fully discrete scheme using artificial compressibility reads

[TABLE]

The corresponding $\Delta$AEequals the $\Delta$AEof the SIMPLE scheme up to a slight modification of the equation for the pressure update such that we can use the arguments laid out in Appendix 4.4 to conclude that it is of index 1.

Remark 4.8*.*

Note that in practice there may be nonlinear solves to determine, e.g., $\tilde{v}^{k+1}$ such that, e.g., the correction $\bar{v}^{k+1}$ is possibly different from $\tilde{v}^{k+1}$ computed with $p^{k}=0$.

4.4. Time-stepping schemes resulting from index reduction

Besides the presented schemes in Section 4.3, one may also apply an index reduction to system (3.1) and then discretize in time. This then also leads to matrix pairs $(\mathcal{E},\mathcal{A})$ of Kronecker index 1.

4.4.1. Penalty methods

Similar to Section 4.3.3, we may reduce the index of the DAE (3.1) by relaxing the divergence-free constraint or, in other words, add a penalty term [She95]. With the penalty parameter $\beta\gg 1$ we replace the incompressibility condition by

[TABLE]

In the semi-discrete case this corresponds to the constraint equation $M_{p}p+Bv=g$. It can be shown that this then leads to a DAE of (differentiation) index 1. However, this approach changes the solution of the system. In order to keep this difference of reasonable size, $\beta$ should be choosen relatively large. On the other hand, the condition number of the involved matrices increase with $\beta$ and thus, lead to numerical difficulties. The difficulty of a reasonable choice of $\beta$ is one drawback of this approach. A second disadvantage is that small velocities of order $\beta^{-1}$ or less cannot be resolved [HV95]. Further modifications of the penalty method, which are also applicable for slightly compressible fluids (Newtonian fluids)), are discussed, e.g., in [HV95, Sect. 5].

4.4.2. Derivative of the constraint

Another very simple possibility to reduce the index of the given DAE is to replace the constraint by its derivative, the so-called hidden constraint. Instead of (3.1) we then consider the system

[TABLE]

The initial condition remains unchanged, i.e., $\tilde{v}(0)=v_{0}$. It is well-known that this system has (differentiation) index 1. Nevertheless, numerical simulations, which rely on this formulation, show a linear drift from the solution manifold given by $Bv=g$. This can be seen as follows.

Although the two systems are equivalent, numerical errors are integrated over time and thus amplified. Solving the constraints only up to a small error, i.e.,

[TABLE]

for small and constant $\varepsilon$ and $\eta$, we calculate for $\tilde{v}$ that

[TABLE]

Note that the last step holds because of the assumed consistency condition $Bv_{0}=g(0)$. This shows that a constant error in the constraint of $\tilde{v}$ leads to an error which grows linearly in time. Because of this, the method of replacing the constraint by its derivative is – although it is of index 1 – not advisable.

4.4.3. Minimal extension

Finally, we present the index reduction technique of minimal extension [KM06, Ch. 6.4]. A general framework for an index reduction based on derivative arrays is given in [KM06, Ch. 6]; see also [Cam87]. Because of the special saddle point structure of system (3.1), in which the constraint is explicitly given, this procedure can be simplified by using so-called dummy variables, cf. [MS93, KM04].

Since we assume that $B$ is of full rank, there exists an invertible matrix $Q\in\mathbb{R}^{n,n}$ such that $BQ$ has the block structure $BQ=[B_{1},\ B_{2}]$ with an invertible matrix $B_{2}\in\mathbb{R}^{m,m}$. Note that the choice of $Q$ is not unique. We then use this transformation to partition the variable $v$, namely

[TABLE]

with according dimensions $v_{1}\in\mathbb{R}^{n-m}$ and $v_{2}\in\mathbb{R}^{m}$. With this, the constraint and its derivative may be written as

[TABLE]

Together with the differential equation (3.1a) this yields an overdetermined system. Instead of an (expensive) search for projectors [KM06, Ch. 6.2] or selectors [Ste06], which would bring the system back to its original size, we introduce a dummy variable $w_{2}:=\dot{v}_{2}$ to get rid of the redundancy. Note, however, that this leads to a slightly bigger system. More precisely, we extend the system dimensions from $n+m$ to $n+2m$, which is still moderate, since we usually have $m\ll n$. The extended system is again square and has the form

[TABLE]

This DAE is of index $1$ and has the same solution set as the original system, cf. [AH15].

The drawback of this method is the need of a transformation matrix $Q$. With a suitable reordering of the variables, however, we can choose $Q$ to be the identity matrix. In this case, the needed variable transformation is just a permutation and thus, all variables keep their physical meaning. For some specific finite element schemes an algorithm to find such a permutation (which leads to an invertible $B_{2}$ block) is given in [AH15]. This paper also considers the implicit-explicit Euler scheme applied to the minimally extended system, which leads to a pair $(\mathcal{E},\mathcal{A})$ of index 1.

5. Numerical experiments

To illustrate the performance and particular issues of the time-stepping schemes for the NSE, we consider the numerical simulation of the flow passing a cylinder in two space dimensions. This problem, also known as cylinder wake, is a popular flow benchmark problem and a test field for flow control [Wil96, NAM*+*03, BH15].

We consider the incompressible NSE (2.14) on the domain as illustrated in Figure 5.1 with boundary $\Gamma$ and boundary conditions as follows. At the inflow $\Gamma_{\text{i}}$, we prescribe a parabolic velocity profile through the function

[TABLE]

At the outflow $\Gamma_{\text{o}}$, we impose do-nothing conditions, cf. (2.17). At the upper and the lower wall of the channel and at the cylinder periphery, we employ no-slip, i.e., zero Dirichlet conditions. Thus, we set

[TABLE]

We set $\text{Re}=60$ and consider a $\mathcal{P}_{2}$–$\mathcal{P}_{1}$ (Taylor–Hood) finite element discretization of (2.14) on the grid depicted in Figure 5.1 with $9356$ degrees of freedom in the velocity and $1289$ degrees of freedom in the pressure approximation. The result of the semi-discretization is a DAE of the form (3.1), which we write as

[TABLE]

on the time interval $(0,1]$. Here, $A$ is the Laplacian or the linear part of $K$ as defined in (3.2) and $N$ denotes the convection part. The right-hand sides $f$ and $g$ account for the (static) boundary conditions. As initial value we take the corresponding steady-state Stokes solution $v_{\text{S}}$, which is part of the solution to

[TABLE]

For the schemes that need an initial value for the pressure, we provide it as $p_{\text{NS}}$ solving (5.1) at $t=0$ with $v(0)=v_{\text{S}}$ and $\dot{v}(0)=0$. Note that this gives a consistent initial pressure, since $g$ is constant, and thus, $B\dot{v}=0$.

We consider the time-discretization of (5.1) by means of the implicit-explicit Euler scheme and compare it to the time-discretization via the SIMPLE scheme as described in Section 4.3.2. Both schemes treat the nonlinearity explicitly. Comparing the computed approximations to a reference, obtained by the time-discretization via the implicit trapezoidal rule on a fine grid, we show that both schemes are convergent of order 1 as long as the resulting linear systems are solved with sufficient precision. For inexact solves, we show that the pressure approximation in the implicit-explicit Euler scheme diverges unlike for the SIMPLE approximation. For snapshots of the approximate velocity solution, see Figure 5.2.

The finite element implementation uses FEniCS, Version 2017.2 [LORW12]. For the iterative solutions of the linear system, we employ Krypy [Gau17]. The code used for the numerical investigations is freely available for reproducing the reported results and as a benchmark for further developments in the time integration of semi-explicit DAEs of strangeness-index 1; see Figure 5.3 for a stable link to the online repository.

5.1. Time integration with implicit-explicit Euler

With the implicit-explicit Euler time discretization, at every time step the linear system

[TABLE]

is solved for the velocity and the pressure approximations. For the approximate solution of the linear systems, we use MinRes iterations, which are stopped as soon as the relative residual drops below a given tolerance tol, i.e.,

[TABLE]

where we use the norm induced by the inverses of the mass matrix $M$ and of the mass matrix of the pressure approximation space $M_{\mathcal{Q}}$.

It can be seen from the error plots in Figure 5.4, for exact solves, that the implicit-explicit Euler converges linearly in the velocity and the pressure approximation, which is in line with the theory [HLR89, Tab. 2.3]. For inexact solves, however, for smaller time-steps, the pressure approximation diverges linearly. As we have shown in [AH15] this is an inherent instability of the index-2 formulation.

The reported residuals and errors are defined as follows:

[TABLE]

where the subscript ref denotes the reference solution, the subscripts $\tau;\texttt{tol}$ denote the approximation that is computed on the grid of size $\tau$ with the linear equations solved with tolerance tol, and where $\operatorname{trp}(s)$ denotes the approximation to the integral $\int_{0}^{t_{N}}s(t)\text{~{}{d}}t$ by means of the piecewise trapezoidal rule with step size $\tau$. Analogously, we define the error in the pressure approximation

[TABLE]

and the integrated residuals in the momentum equation $r_{\tau;\texttt{tol}}^{M}$ as the integral of the function

[TABLE]

and in the continuum equation as

[TABLE]

5.2. Time integration with SIMPLE

For this particular time discretization, we have to solve three linear systems, namely

[TABLE]

to compute the updates as $v^{k+1}=\tilde{v}^{k+1}+v_{\Delta}^{k+1}$ and $p^{k+1}=p^{k}+p_{\Delta}^{k+1}$. For the approximate solution, we use CG iterations until the relative residuals, measured in the $M^{-1}$ norm (or $M_{\mathcal{Q}}^{-1}$ for (5.4b)), drop below a given tolerance tol.

As one can see from the error plots in Figure 5.5, the inexact solves affect the approximation only for a rough tolerance $\texttt{tol}=10^{-4}$, which is only one order of magnitude smaller than the actual approximation error. Thus the breakdown in the convergence observed in Figure 5.5 is probably due to the accumulation of errors that every single step scheme suffers from. For smaller tolerances, the approximations almost achieve the accuracy of the direct solves. Interestingly, the continuity equation (3.1b), which is only an implicit part of the SIMPLE scheme, is fulfilled at much better accuracy than the choices of the tolerances suggest; see Figure 5.5.

6. Conclusion

In this paper, we have discussed the incompressible NSE from a DAE point of view, in which the incompressibility is interpreted as an algebraic constraint. Thus, a spatial discretization leads to a DAE. If the time is discretized as well, a $\Delta$AE is obtained – a sequence of equations, that define the numerical approximations.

We have discussed suitable approaches to well-posed semi-discrete approximations and investigated the $\Delta$AE stemming from different time-discretizations in terms of the Kronecker index. It turned out that commonly and successfully used time integration schemes like the SIMPLE algorithm define a $\Delta$AE of index 1, whereas a time-discretization of the semi-discrete NSE by an implicit-explicit Euler scheme leads to a $\Delta$AE of index 2. Alternatively, one may first apply an index reduction on the semi-discrete level and then apply time marching schemes.

The advantage of the discrete index-1 formulations is that they avoid implicit derivations that amplify computational errors. This mechanism will likely lead to larger errors, in particular if the equations are solved with limited accuracy. We have illustrated the origin of this behavior in a small analytical example and verified it in a numerical simulation. The code of the numerical test case, an implementation of the 2D cylinder wake, is provided to serve as a benchmark for future developments of time integration schemes for DAEs of Navier–Stokes type.

Appendix A Strangeness index of equation (3.1)

We analyse in detail the strangeness-index of the DAE (3.1). Note that we do not ask for any assumptions on the nonlinearity $K$ and that we allow the matrix $B$ to be rank-deficient. This then also implies that the d-index of (3.1) equals 2 if it is well-defined, i.e., if $B$ is of full rank.

A.1. Linear case

Considering any linearization of the Navier-Stokes equations, i.e., $K(u)=Ku$ in (3.1), we deal with the matrix pair

[TABLE]

Following [KM06, Th. 3.11], we can construct a to $(E,A)$ (globally) equivalent pair $(\tilde{E},\tilde{A})$ of the form

[TABLE]

with $A_{13}\in\mathbb{R}^{b,m}$ being of full rank. Thus, the original system (3.1) is equivalent to a system of the form

[TABLE]

with dimensions $x_{1}(t)\in\mathbb{R}^{b}$, $x_{2}(t)\in\mathbb{R}^{n-b}$, $x_{3}(t)\in\mathbb{R}^{m}$. Since we have a differential and an algebraic equation for $x_{1}$ (this causes the ’strangeness’), we use the derivative of $0=x_{1}+f_{3}$ in order to eliminate $\dot{x}_{1}$ in the first equation. Hence, we consider the pair $(E_{\text{mod}},A_{\text{mod}})$ with

[TABLE]

Since $A_{13}$ is of full rank, one can show that system $(E_{\text{mod}},A_{\text{mod}})$ is strangeness-free, cf.  the calculation in [KM06, Th. 3.7]. Since we have obtained a strangeness-free system with only one differentiation, system (3.1) has strangeness-index one.

A.2. Nonlinear case

The general form of a nonlinear DAE is given by

[TABLE]

In regard of system (3.1) we set $x:=[q^{T},\ p^{T}]^{T}$ and define

[TABLE]

with

[TABLE]

In the sequel we show that (3.1) has strangeness-index 1 also in the nonlinear case. For this, we assume that $B$ has full rank such that there are no vanishing equations and the pressure variable is uniquely defined. In the case $\operatorname{rank}B=b<m$, we consider the following transformation.

Let $C_{0}\in\mathbb{R}^{m,m-b}$ be the matrix of full rank satisfying $B^{T}C_{0}=0$. Furthermore, $C^{\prime}\in\mathbb{R}^{m,b}$ defines any matrix such that $C=[C_{0}\ C^{\prime}]\in\mathbb{R}^{m,m}$ is invertible. With this, we obtain the relation

[TABLE]

with $\tilde{B}\in\mathbb{R}^{b,n}$ having full rank. With the matrix $C$ in hand, we first introduce the new pressure variable $\tilde{p}:=C^{-1}p$. Thus, we consider the the pair $z:=[q^{T},\tilde{p}^{T}]^{T}$. As a second step, we multiply equation (3.1) by the block-diagonal matrix $\operatorname{diag}(I_{n},C^{T})$ from the left. In total, this yields the equivalent DAE

[TABLE]

Note that the constraint contains $(m-b)$ consistency equations of the form $0=g_{1}$. Assuming that system (3.1) is solvable, we suppose that these are in fact vanishing equations. Thus, they have no influence on the index of the system. Furthermore, the first $(m-b)$ components of the transformed pressure $\tilde{p}$ do not influence the system. These components are underdetermined and may be omitted, again without changing the index. Leaving out the underdetermined parts as well as the vanishing equations, we obtain a system of the form (3.1) with a full rank matrix $B$.

In the sequel, we assume that $\operatorname{rank}B=m$ and show that [KM06, Hyp. 4.2] is satisfied for $\mu=1$. Note that this hypothesis is not satisfied for $\mu=0$, i.e., the system is not strangeness-free. We define the matrices

[TABLE]

We now pass through the list of points of the hypothesis in [KM06, Hyp. 4.2]:

1. (1)

First, we note that the rank of $M_{1}$ equals $2n$ and we set $a:=2(n+m)-2n=2m$. Thus, the system contains $2m$ algebraic variables (the pressure and the part of $q$, which is not divergence-free). Furthermore, we define $Z_{2}\in\mathbb{R}^{2(n+m),2m}$ by $Z_{2}^{T}M_{1}=0$, i.e.,

[TABLE] 2. (2)

As a second step we define $\hat{A}_{2}:=Z_{2}^{T}N_{1}[I_{n+m},\ 0]^{T}$, which yields

[TABLE]

This matrix has rank $2m$, since the full-rank property of $B$ implies that $BM^{-1}B^{T}$ is invertible. We define $d:=n-m$ as the number of differential variables and $T_{2}\in\mathbb{R}^{n+m,n-m}$ by $\hat{A}_{2}T_{2}=0$. Let $C\in\mathbb{R}^{n,n-m}$ be a matrix of full rank with $BC=0$ and $C_{2}:=-(BM^{-1}B^{T})^{-1}BM^{-1}K_{u}C\in\mathbb{R}^{m,n-m}$. Then, we set

[TABLE] 3. (3)

Finally, we compute the rank of $ET_{2}$. Since $C$ has full rank, this equals $\operatorname{rank}MC=n-m=d$. The matrix $Z_{1}^{T}:=[C^{T}\ 0]\in\mathbb{R}^{n-m,n+m}$ satisfies

[TABLE]

Thus, the hypothesis in [KM06, Hyp. 4.2] is satisfied for $\mu=1$, which implies that the nonlinear DAE (3.1) has strangeness-index one.

Appendix B Difference-algebraic equation index of the considered systems

In this appendix, we derive the Kronecker index for the discrete schemes considered in Section  4.3.

B.1. IMEX Euler

We start with the implicit-explicit Euler discretization, that gives a scheme $\mathcal{E}x^{k+1}=\mathcal{A}^{k}x^{k}+h^{k}$ with the matrix pair

[TABLE]

as in (4.5). For sufficiently small $\tau$, due to the definiteness of $M$ and the full-rank property of $B$, the matrix $\mathcal{A}$ is invertible and thus, the pair $(\mathcal{E},\mathcal{A})$ is regular. Let $S$ denote the matrix $BM^{-1}B^{T}$. If one applies

[TABLE]

from the left and the right, one finds that $(\mathcal{E},\mathcal{A})$ is similar to

[TABLE]

Since $B$ is of full rank, there exists an orthogonal matrix $Q$ and an invertible matrix $R$ such that BM^{-\frac{1}{2}}Q=\big{[}0\ \ R\big{]} and, in particular,

[TABLE]

Thus, the corresponding similarity transformation transforms $(\mathcal{E},\mathcal{A})$ into

[TABLE]

where $*$ stands for the block matrix entries that need not be specified further. With another few regular row and column transformations, one can eliminate the entry $*_{2}$ and read off the Kronecker index of \bigl{(}\mathcal{E},\mathcal{A}\bigr{)} as the index of nilpotency of $\begin{bmatrix}0&0\\ \tfrac{1}{\tau}R&0\end{bmatrix}$ which is $2$.

B.2. Projection scheme

The matrix coefficient pair of the Projection scheme (4.10) reads

[TABLE]

If we define $S:=BM^{-1}B^{T}$, if we move the second row and column to the left and bottom, respectively, and if we rescale certain rows and columns, we find that the pair is equivalent to

[TABLE]

where the $*$’s stand for unspecified but possibly nonzero entries. Since, in particular, $S$ is invertible, one can eliminate the entries $*_{1}$–$*_{4}$ by regular row and column manipulations without affecting the invertibility of the left upper $3\times 3$ block in the transformed $\mathcal{E}$ and read off the Kronecker index of (4.10) as the index of nilpotency of [math] which is $1$.

B.3. SIMPLE

The matrix coefficient pair of the SIMPLE scheme (4.14) reads

[TABLE]

If we define $S_{A}:=B(\frac{1}{\tau}M^{-1}+A)^{-1}B^{T}$, move the second row and column to the left and bottom, respectively, and rescale certain rows and columns, then we find that the pair is equivalent to

[TABLE]

where the $*$’s stand for unspecified but possibly nonzero entries. Since $S_{A}$ is invertible for sufficiently small $\tau$, we find that this matrix pair has the very same structure as the one of the projection scheme (see Appendix B.2) and, thus, is of index-1.