The analogue of grad–div stabilization in DG methods for incompressible flows: Limiting behavior and extension to tensor-product meshes

https://doi.org/10.1016/j.cma.2018.07.019Get rights and content

Abstract

grad–div stabilization is a classical remedy in conforming mixed finite element methods for incompressible flow problems, for mitigating velocity errors that are sometimes called poor mass conservation. Such errors arise due to the relaxation of the divergence constraint in classical mixed methods, and are excited whenever the spatial discretization has to deal with comparably large and complicated pressures. In this contribution, an analogue of grad–div stabilization for Discontinuous Galerkin methods is studied. Here, the key is the penalization of the jumps of the normal velocities over facets of the triangulation, which controls the measure-valued part of the distributional divergence of the discrete velocity solution. Our contribution is twofold: first, we characterize the limit for arbitrarily large penalization parameters, which shows that the stabilized nonconforming Discontinuous Galerkin methods remain robust and accurate in this limit; second, we extend these ideas to the case of non-simplicial meshes; here, broken grad–div stabilization must be used in addition to the normal velocity jump penalization, in order to get the desired pressure robustness effect. The analysis is performed for the Stokes equations, and more complex flows and Crouzeix–Raviart elements are considered in numerical examples that also show the relevance of the theory in practical settings.

Introduction

Classical conforming and inf–sup stable mixed finite element methods for the incompressible (Navier–)Stokes equations, such as the mini [1], the Bernardi–Raugel [2] and the Taylor–Hood elements [3], relax the divergence constraint u=0, in order to construct optimally convergent spatial discretizations on regular unstructured triangulations [[4], [5]]. Indeed, the discrete velocity solution uh is not divergence-free, but only discretely divergence-free, i.e., it holds huhπ0(uh)=0, where h denotes the discrete divergence operator and π0 denotes the L2 best approximation in the discrete pressure space.

While relaxing the divergence constraint facilitates the construction of inf–supstable discretizations, it was soon realized that discretely divergence-free is sometimes not good enough. For example, already in 1989 D. Pelletier, A. Fortin and R. Camarero titled their article [6] by the provocative question “Are FEM solutions of incompressible flows really incompressible (or how simple flows can cause headaches!)” and pointed to the problem of poor mass conservation. Poor mass conservation describes velocity errors that are excited when the pressure is comparably large and complicated [7]. This happens in many flow problems, including Boussinesq flows [[8], [6], [9], [10], [11], [12]], potential and generalized Beltrami flows [[13], [14], [7]], quasi-geostrophic flows [[15], [16], [14]], electrophoresis [17], and two-phase flows with surface tension [[18], [19]]. Physical problems where poor mass conservation is not strong are rather limited, and include for example pressure-driven Stokes and pressure-driven, low-Reynolds number Navier–Stokes flows through a channel with zero exterior forcing, where for the momentum balance approximately holds νΔu+p0, i.e., where the pressure gradient is proportional to the friction force.

Phenomenologically, poor mass conservation is often accompanied by comparably large violations of the divergence error uhL2, see for example [12, Table 1]. To address this, researchers since the late 1970’s have enhanced the Navier–Stokes momentum balance of conforming mixed finite element methods by a consistent term [[20], [21]] γgduh,vh,which penalizes large divergence errors, and is nowadays often called grad–div stabilization; here γgd0 denotes the grad–div stabilization parameter. Indeed, grad–div stabilization for conforming mixed finite element methods has recently been investigated in depth, both from a theoretical and computational point of view [[22], [23], [24], [12], [25], [26], [27], [28], [29]]. A better understanding of grad–div stabilization was achieved, when the limit behavior for arbitrarily large stabilization parameters γgd was investigated [[26], [27], [28]]; it turned out that grad–div stabilization is not so much a stabilization, but instead a kind of penalization procedure. On a fixed grid, for γgd the grad–div stabilized discrete velocity solution uhγgd converges to a divergence-free velocity solution uh which is the solution of a divergence-free conforming mixed finite element method with the same discrete velocity space, but with a richer discrete pressure space. Since divergence-free, conforming mixed finite element methods are pressure-robust [[7], [13]], i.e., their velocity error does not depend on the continuous pressure, grad–div stabilized discrete velocities behave in a more robust manner against large and complicated continuous pressures. However, this theoretical understanding also revealed limitations of grad–div stabilization in that large grad–div stabilization parameters can cause classical Poisson locking phenomena, whenever the limiting divergence-free mixed method is not inf–sup stable [28]. On the other hand, on certain mesh families and for certain conforming mixed finite elements, grad–div over-stabilization can be avoided [28, Corollary 1, Case 2].

Of course, nonconforming mixed methods like the Crouzeix–Raviart finite element method [30] are as much endangered by poor mass conservation as conforming ones, since they are not pressure-robust (their velocity error depends on the continuous pressure). In [[31], [32]] it was recognized that pressure-robustness of a mixed method for incompressible flows does not depend on the fact that the discrete velocity trial functions are divergence-free, but it only depends on the discrete velocity test functions. They have to be divergence-free in the weak sense of Hdiv, in order to be orthogonal to any gradient field in the L2(Ω) scalar product for vector fields [[7], [31]]. Recall that a vector field vL2(Ω) is said to be weakly divergence-free, if its distributional divergence [7] ϕΩvϕdxvanishes for all ϕC0(Ω), i.e., if it is orthogonal in L2(Ω) to all smooth gradient fields [7]. This shows that it is actually a very strong property for a vector field v to be weakly divergence-free, at least compared to being only discretely divergence-free. Indeed, applying the general definition of the distributional divergence to nonconforming finite element methods and Discontinuous Galerkin (DG) methods, it turns out that the distributional divergence of a velocity test function vh is given by ϕKThKϕvhdxFFhiFϕvhnFds,since vh is elementwise polynomial and an integration by parts can be applied. Therefore, the distributional divergence of vh vanishes only if it holds vh=0 elementwise for all KTh and vhnF=0 for all FFhi. Instead, if a vector field vh is discretely divergence-free, this implies for usual, i.e., not pressure-robust [7], discretizations that only ϕhKThKϕhvhdxFFhiFϕhvhnFds,vanishes for all ϕhQh from a finite-dimensional space Qh of pressure test functions. Hence, an analogue for grad–div stabilization for nonconforming methods has to penalize the elementwise defined broken divergence hvh, and the facet jumps vhnF of the normal velocities (i.e., the mass flux). Therefore, in the following a penalization procedure with a penalization parameter γ>0 is employed. This type of penalization for DG methods has been independently proposed in [33] and [34], and the analysis in [33] showed that it provides a benefit to velocity error analogous to what grad–div stabilization provides for conforming elements: it improves the velocity error by reducing the contribution of the pressure discretization error. In the context of the Crouzeix–Raviart finite element method, the importance of penalizing the jumps of the normal velocities was recognized by E. Burman and P. Hansbo in 2005 [35].

Herein, we make two fundamental contributions to advance the study of this type of stabilization. First, we consider the limiting behavior as the mass flux penalization parameter γ, which is important when the viscosity is small and/or the pressure is large. We note that the limiting behavior has not yet been considered in the context of grad–div stabilization, but was considered in [36] in a similar way for constructing a ‘divergence-free’ HDG method applied to a ‘gradient-velocity’ formulation of the Stokes equations. In this limit, the discrete velocity solution uh will be elementwise divergence-free and the jumps of the normal velocities over inner facets will vanish. Often, taking stabilization parameters large can cause over-stabilization, however we prove that the limit solution is the discrete velocity solution of other DG methods, namely of the weakly divergence-free inf–sup stable DG method proposed in [[37], [38]]. Hence the stabilized method is perfectly robust against over-stabilization, and moreover, this suggests that under an assumption of equal computational cost, computing the limit solution directly using the methods of [[37], [38]] would be advantageous over the associated stabilized DG method.

The second major contribution of this work is to extend the study of this stabilization to DG methods on tensor-product meshes. The key idea here is to use the stabilization above for the mass flux and an elementwise divergence penalization. In [39] an elementwise divergence penalization has been used (seemingly) for the first time. It was shown that even this incomplete stabilization method can improve poor mass conservation for both inf–sup stable QkdcQk1dc and equal-order DG methods QkdcQkdc on tensor-product meshes in practical applications of non-isothermal flows, with the superscript ‘dc’ denoting that these elements are piecewise discontinuous. Further, it was recognized independently in the recent works [34] and [40] that in DG methods, the mass balance across interior facets has to be accounted for, in addition to the classical (broken) grad–div stabilization. For laminar and turbulent flows, the authors of Ref. [34] –improving earlier results [40] –use a velocity-correction time integration with equal-order QkdcQkdc elements on tensor-product meshes and an implementation in deal.ii. The additional terms they favor are slightly different from the stabilization that is considered in this work, see [34, 4.2.3 and 4.2.4]. Further, the authors do not rigorously justify their approach by numerical analysis and do not consider a possible lack of robustness against over-stabilization.

Organization of the article: In Section 2 we provide some notation and mathematical preliminaries to allow for a cleaner analysis to follow. Section 3 considers the case of the mass flux penalization for DG on simplices, Section 4 considers tensor product meshes, and Section 5 considers the enhancement in Crouzeix–Raviart elements. Several numerical tests of concept are given in Sections 3 Mass flux penalization applied to inf–sup stable DG methods on simplicial meshes, 4 Stabilization of inf–sup stable DG methods on tensor-product meshes, and in Section 6, we consider applications of the penalization outside of the Stokes setting, to Navier–Stokes equations and to Boussinesq equations. In all numerical tests the (sometimes dramatic) improvement offered by the penalization is clear. Finally, conclusions are drawn in Section 7, and future research directions are discussed.

Section snippets

Stokes problem and DG setting

We consider a domain ΩRd, d = 2,3, to be a simply connected set with smooth boundary, or a convex polygon. For KΩ we use the standard Sobolev spaces Wm,pK for scalar-valued functions with associated norms Wm,pK and seminorms ||Wm,pK for 0mR and p1. Spaces and norms for vector- and tensor-valued functions are indicated with bold letters. We use the Lebesgue space LpK=W0,pK and the Hilbert space HmK=Wm,2K. Additionally, the closed subspaces H01K consisting of H1K-functions with vanishing

Mass flux penalization applied to inf–sup stable DG methods on simplicial meshes

We consider in this section analysis of DG methods on simplicial meshes, with the mass flux penalization. We note that much of Section 3.1 has essentially been done in [33] where the spatial convergence of steady Navier–Stokes with the same DG discretization and mass flux penalization is considered. We include the reduction of their analysis to the Stokes case herein because it defines the framework for our study of the limiting behavior in Section 3.4, and it gives a starting point from which

Stabilization of inf–sup stable DG methods on tensor-product meshes

As mentioned in the previous section, on tensor-product elements (quadrilateral and hexagons) using QkdcQk1dc, the situation is slightly more involved since hQkdcQk1dc. In order to demonstrate the difference, we repeat the no-flow test from the previous section with QkdcQk1dc elements on a structured mesh consisting of squares. The results (obtained with NGSolve) are shown in Table 3.

Interestingly, neither mass flux penalization alone (γ>0, γgd=0), nor broken grad–div stabilization

Improving pressure-robustness in Crouzeix–Raviart approximations

The mass flux penalization can also be applied, with similar results, to the Crouzeix–Raviart (CR) element, and we include some results for completeness. In a sense, it is somewhat easier to analyze than the DG case considered above, and we find a similar fundamental result: penalization of the mass flux reduces the effect of the pressure error on the velocity error. Such a result is proven implicitly for CR elements in a recent work of Burman and Hansbo [35] for the Darcy–Stokes problem, where

Application to more complex flows

We show here that the pressure-robust approach above can have a considerable impact on more complicated problems than the steady incompressible Stokes equations. We consider first a numerical test for steady incompressible Navier–Stokes equations (NSE), followed by a numerical test for non-isothermal flow.

Conclusions and future directions

We have analyzed and tested a penalization(s) for nonconforming methods that has an effect on solutions analogous to that of grad–div stabilization for conforming methods. The penalization is a mass flux penalization, and also a broken divergence stabilization in elements where the divergence of the velocity space is not contained in the pressure space.

Our new theoretical results include error estimates that reduce the scaling of the velocity error by the pressure from ν1 to γ12ν12,

Acknowledgments

The authors would especially like to thank Christoph Lehrenfeld for several related fruitful discussions on stabilization and hybridization and the invaluable help he provided in using the finite element library NGSolve in the context of this work. Mine Akbas acknowledges support from the German Academic Exchange Service (DAAD) with the program “Research Grants for Doctoral Candidates and Young Academics and Scientists”, 2017/18 (57299291). The third author was supported by National Science

References (52)

  • LinkeA.

    A divergence-free velocity reconstruction for incompressible flows

    C. R. Math. Acad. Sci. Paris

    (2012)
  • KrankB. et al.

    A high-order semi-explicit discontinuous Galerkin solver for 3D incompressible flow with application to DNS and LES of turbulent channel flow

    J. Comput. Phys.

    (2017)
  • SchroederP.W. et al.

    Stabilised dG-FEM for incompressible natural convection flows with boundary and moving interior layers on non-adapted meshes

    J. Comput. Phys.

    (2017)
  • JoshiS.M. et al.

    A post-processing technique for stabilizing the discontinuous pressure projection operator in marginally-resolved incompressible inviscid flow

    Comput.& Fluids

    (2016)
  • LehrenfeldC. et al.

    High order exactly divergence-free hybrid discontinuous Galerkin methods for unsteady incompressible flows

    Comput. Methods Appl. Mech. Engrg.

    (2016)
  • ArnoldD.N. et al.

    A stable finite element for the Stokes equations

    Calcolo

    (1984)
  • BernardiC. et al.

    Analysis of some finite elements for the Stokes problem

    Math. Comp.

    (1985)
  • VerfürthR.

    Error estimates for a mixed finite element approximation of the Stokes equation

    RAIRO Anal. Numer.

    (1984)
  • GiraultV. et al.

    Finite Element Methods for Navier–Stokes Equations

  • BoffiD. et al.

    Mixed Finite Element Methods and Applications

    (2013)
  • PelletierD. et al.

    Are FEM solutions of incompressible flows really incompressible? (or how simple flows can cause headaches!)

    Internat. J. Numer. Methods Fluids

    (1989)
  • JohnV. et al.

    On the divergence constraint in mixed finite element methods for incompressible flows

    SIAM Rev.

    (2017)
  • GreshoP.M. et al.

    A new finite element for incompressible or Boussinesq fluids

  • DorokO. et al.

    Aspects of finite element discretizations for solving the Boussinesq approximation of the Navier–Stokes Equations

  • GerbeauJ.F. et al.

    Spurious velocities in the steady flow of an incompressible fluid subjected to external forces

    Internat. J. Numer. Methods Fluids

    (1997)
  • FrolkovicP.

    Consistent velocity approximation for density driven flow and transport

  • Cited by (25)

    • An equal-order hybridized discontinuous Galerkin method with a small pressure penalty parameter for the Stokes equations

      2021, Computers and Mathematics with Applications
      Citation Excerpt :

      However, taking grad-div stabilization parameter too large will lead to an over-stabilizing effect in conforming finite element methods [17,18]. For the equal-order DG method, both grad-div stabilization and mass flux penalization are usually adopted simultaneously by following the idea of [19]. In addition, the DG methods are known to be computationally expensive.

    View all citing articles on Scopus
    View full text