Towards computable flows and robust estimates for inf-sup stable FEM applied to the time-dependent incompressible Navier–Stokes equations

Abstract

Inf-sup stable FEM applied to time-dependent incompressible Navier–Stokes flows are considered. The focus lies on robust estimates for the kinetic and dissipation energies in a twofold sense. Firstly, pressure–robustness ensures the fulfilment of a fundamental invariance principle and velocity error estimates are not corrupted by the pressure approximability. Secondly, Re-semi-robustness means that constants appearing on the right-hand side of kinetic and dissipation energy error estimates (including Gronwall constants) do not explicitly depend on the Reynolds number. Such estimates rely on the essential regularity assumption $\nabla u \in L^1(0,T;L^\infty (\varOmega ))$ which is discussed in detail. In the sense of best practice, we review and establish pressure- and Re-semi-robust estimates for pointwise divergence-free $H^1$-conforming FEM (like Scott–Vogelius pairs or certain isogeometric based FEM) and pointwise divergence-free H(div)-conforming discontinuous Galerkin FEM. For convection-dominated problems, the latter naturally includes an upwind stabilisation for the velocity which is not gradient-based.

Appendix: Computational aspects of \(\mathbf H \)(div)-conforming methods for Navier–Stokes

Improving the efficiency of ${\varvec{H}{{\left( {\mathrm {div}}\right) }}}$-conforming methods In this section we explain how ${\varvec{H}{{\left( {\mathrm {div}}\right) }}}$-conforming FE methods, that are often seen as too complicated and inefficient for real application, can be made efficient. We restrict the discussion here to BDM elements as they are computationally more efficient in the context of incompressible flows compared to RT elements since they have less degrees of freedom (DOFs) for the same velocity approximation.

Choosing the pressure space ${Q}_h$ as the space of (discontinuous) piecewise polynomials of one degree less than the ${\varvec{H}{{\left( {\mathrm {div}}\right) }}}$-conforming velocity space ${\varvec{V}}_h$ renders (16) an equality, that is, ${{\mathrm{\nabla \,\varvec{\cdot }}}}{\varvec{V}}_h = {Q}_h$. A special property of this velocity–pressure pair is that the inf-sup constant is robust in the polynomial degree leading to hp-optimal convergence; cf. [44] for a rigorous analysis in 2D. The strong relation ${{\mathrm{\nabla \,\varvec{\cdot }}}}{\varvec{V}}_h = {Q}_h$ can further be exploited with a smart choice of the basis functions for ${\varvec{V}}_h$ and ${Q}_h$; cf. [63, 67]. The a priori knowledge that the discrete solution will be pointwise divergence-free then allows to remove some DOFs for the velocity and all pressure unknowns except for the piecewise constants; cf. [46, Remark 1] and [45, Section 2.2.4.2]. We make use of this in our numerical experiments.

To account for the tangential discontinuity in the ${\varvec{H}{{\left( {\mathrm {div}}\right) }}}$-conforming FE space, a DG formulation has to be applied. This aspect can be regarded ambivalently. On the one hand, the discontinuous nature of the tangential component offers the possibility of applying an upwind discretisation for the convection, cf. (22), which results in stable discretisations also in the convective limit [35] without adding too much dissipation compared to most convection stabilisations of ${\varvec{H}^{1}{}}$-FEM. On the other hand, the DG formulation results in computationally less attractive features. Due to the break-up of the tangential continuity, several DOFs for the velocity are multiplied compared to ${\varvec{H}^{1}{}}$-conforming methods. Even worse, the number of couplings in a corresponding system matrix increases which results in much higher computational costs for (direct and iterative) solvers of linear systems.

Several measures can be taken to compensate for these costs. To this end, we briefly discuss the concept of hybridisation in the context of DG methods [22]. To reduce the couplings of neighbouring elements, additional unknowns on the facets are introduced (which typically approximate the trace of the unknown field). These additional unknowns are used to replace the direct couplings of neighbouring elements with couplings between element unknowns and the facet unknowns. Due to the lower dimension of the facets, this reduces the overall amount of couplings especially in the higher order case. More importantly, it allows for static condensation, i.e. the elimination of interior unknowns by a local Schur complement strategy which reduces the number of DOFs for which a global linear system needs to be solved.

Fig. 5

Sketch of fourth order FE discretisations with different types of unknowns for velocity and pressure: unknowns that can be remove beforehand if a suitable basis is used, local unknowns that can be eliminated by static condensation and the remaining global unknowns. The space $\mathbf{F}_{3}$ in the ${\varvec{H}{{\left( {\mathrm {div}}\right) }}}$-HDG method is the space of vector-valued functions that are tangential polynomials up to degree three on every facet

Depending on the problem at hand there are many ways to make use of hybridisation. For an overview we refer to the review article [18]. For Stokes and Navier–Stokes discretisations many variants have been considered; see, for instance, [19,20,21]. Exactly divergence-free HDG methods have also been considered in [17, 57, 58] where additional facet unknowns can be used to enforce normal continuity on a standard DG space which circumvents the construction of ${\varvec{H}{{\left( {\mathrm {div}}\right) }}}$-conforming FE spaces. Here, we use the formulation presented in [46] where, additionally to an ${\varvec{H}{{\left( {\mathrm {div}}\right) }}}$-conforming FE space ${\varvec{V}}_h$ for the velocity and a discontinuous pressure space ${Q}_h$, facet unknowns are introduced only for the tangential component of the velocity. The DG terms in the variational formulation are then adjusted correspondingly. Finally, the element unknowns of the ${\varvec{H}{{\left( {\mathrm {div}}\right) }}}$-conforming FE space couple with neighbour elements only through facet unknowns. These facet unknowns are either the DOFs for the normal continuity of ${\varvec{V}}_h$ or the additional facet unknowns. All remaining velocity unknowns, as well as the pressure unknowns, have only element local couplings such that these—except for the mean element pressure—can be eliminated during static condensation; cf. Fig. 5 for a sketch.

In the viscosity dominated case hybridisation can be optimised further so that only facet unknowns of one degree less need to be considered; cf. [46, 55, 56]. A similar optimisation can also be made for the unknowns for the normal continuity by relaxing the ${\varvec{H}{{\left( {\mathrm {div}}\right) }}}$-conformity slightly. We do not treat this here but instead refer to [42]. To make use of these superconvergence properties of HDG methods we apply—as suggested in [46]—an operator splitting time integration method where the convection operator is treated only explicitly while the remaining time-independent operators are treated implicitly. Note that such an operator splitting is not only desirable for hybrid DG discretisations. Several time integration methods allow for such a splitting; cf. [46, Section 3]. For the experiments in Sect. 5.3 a second-order implicit-explicit BDF2 method has been used.

Some performance comparisons for the numerical study in Sect. 5.3 In Sect. 5.3 the errors for Taylor–Hood, Scott–Vogelius, BDM and the hybridised BDM FE discretisation on two different meshes are compared. At this point, this study shall be complemented with information on the computational costs of the methods. The results are shown in Table 1 where we make this comparison only in terms of the following measures. Firstly, the numbers of DOFs for velocity and pressure ($\#{\left\{ {\varvec{u}}\,{\mathrm {DOFs}} \right\} }$, $\#{\left\{ p\,{\mathrm {DOFs}} \right\} }$, $\#{\left\{ {\mathrm {DOFs}} \right\} }$) are compared. Secondly, we consider the same numbers that remain in a global linear system after static condensation and a potential reduction of the basis (in brackets). Thirdly, the non-zero entries in the global matrix $M^*$ before ($\#{\left\{ \text {nz}(M^*) \right\} }$) and after reduction and static condensation (in brackets) are considered. Note that these numbers can only give an indication of the computational efficiency of the methods. Many different practically relevant aspects, as for example parallelisability or the availability and performance of suitable preconditioners, are not reflected in these numbers.

Regarding static condensation in the Taylor–Hood method, independent of the grad-div stabilisation, we can eliminate all interior unknowns for velocity and pressure. On general meshes, the pressure unknowns for the Scott–Vogelius element cannot be eliminated and hence, static condensation is only applied with respect to the interior velocity DOFs. We note that on barycentre refined meshes static condensation can also be applied for the pressure unknowns; cf. [24]. In case of a DG formulation with BDM elements we utilise the special basis introduced in [63, 67] to eliminate some velocity unknowns and all pressure unknowns except for the mean element pressure. However, static condensation cannot be applied to any additional DOFs due to the DG couplings. Note that this could potentially be improved slightly by choosing a nodal basis similar to the one in [37] where interior unknowns only couple with the boundary nodes of neighbouring elements. To the best of the authors’ knowledge, such a basis has not yet been proposed for an ${\varvec{H}{{\left( {\mathrm {div}}\right) }}}$-conforming FE space. For the hybridised DG method we can apply the reduction of the basis for the ${\varvec{H}{{\left( {\mathrm {div}}\right) }}}$-conforming FE space as well as static condensation. Note that in this work, the formulation from [46] is used which only involves tangential facet unknowns of degree 7. The results are shown in Table 1.

Table 1 Overview of meshes, DOFs and non-zero entries of $M^*$

We observe that the effect of the basis reduction and especially the hybridisation reduces the computational costs of the ${\varvec{H}{{\left( {\mathrm {div}}\right) }}}$-conforming methods drastically, thereby rendering them competitive not only in terms of accuracy, cf. Sect. 5.3, but also in terms of computing time; see also the benchmark results in [46, Section 4.5].