A locally conservative LDG method for the incompressible Navier-Stokes equations

By Bernardo Cockburn, Guido Kanschat, and Dominik Schötzau


In this paper a new local discontinuous Galerkin method for the incompressible stationary Navier-Stokes equations is proposed and analyzed. Four important features render this method unique: its stability, its local conservativity, its high-order accuracy, and the exact satisfaction of the incompressibility constraint. Although the method uses completely discontinuous approximations, a globally divergence-free approximate velocity in is obtained by simple, element-by-element post-processing. Optimal error estimates are proven and an iterative procedure used to compute the approximate solution is shown to converge. This procedure is nothing but a discrete version of the classical fixed point iteration used to obtain existence and uniqueness of solutions to the incompressible Navier-Stokes equations by solving a sequence of Oseen problems. Numerical results are shown which verify the theoretical rates of convergence. They also confirm the independence of the number of fixed point iterations with respect to the discretization parameters. Finally, they show that the method works well for a wide range of Reynolds numbers.

1. Introduction

In this paper we propose and analyze a local discontinuous Galerkin (LDG) method for the stationary incompressible Navier-Stokes equations

Here is the kinematic viscosity, the velocity, the pressure, and the external body force. For the sake of simplicity, we take to be a polygonal domain of .

This paper is the fourth in a series (Reference 9, Reference 8, and Reference 7) devoted to the study of the LDG method as applied to incompressible fluid flow problems. In Reference 9, we considered the Stokes equations

and focused on the problem of how to deal with the incompressibility condition. Later, in Reference 8, we considered the Oseen equations

where the convective velocity was taken to be a smooth function, and focused on the problem of how to incorporate the linear convective term. The resulting method was shown to be optimally convergent and robust for a wide range of Reynolds numbers. A succinct review of this work can be found in Reference 7.

In this paper, we continue our study of LDG methods for incompressible flows and consider their application to the Navier-Stokes equations Equation 1.1. Our main concern is to devise LDG methods that are locally conservative and can be proven to be stable. The local conservativity, a property highly valued by practitioners of computational fluid dynamics, is a discrete version of the identity

where is an arbitrary subdomain of with outward normal unit vector . This property can be easily enforced by LDG methods as soon as the equations are written in divergence form. However, to devise an LDG method that can be proven to be stable (that is, that satisfies a discrete version of the stability estimate for the continuous case,

where the constant is the Poincaré constant) is extremely difficult.

The reason for this is that, in order to obtain the stability estimate Equation 1.5, the incompressibility condition must be used. It is well known that, for many numerical methods for the Stokes and the Oseen equations, a weakly enforced incompressibility is enough to guarantee stability. However, this is not so for the Navier-Stokes equations because of the presence of nonlinear convection. Moreover, the now standard solution to this problem, Reference 21, Reference 22, which is based on a suitable modification of the nonlinearity of the Navier-Stokes equations, cannot be used. This happens because such a modification does not have divergence form and hence prevents LDG methods from being locally conservative.

The main contribution of this paper is to show how to overcome this difficulty. In fact, we show that this can be done in two ways. The first one focuses on the convective nonlinearity and is based on a new modification in divergence form of the nonlinearity; it will be explored elsewhere in detail. The second one, which constitutes the main subject of this paper, focuses on the incompressibility constraint and is based on discretizing the Oseen equations Equation 1.3 where the convective velocity is taken to be a projection of the approximate velocity ,

into the space of globally divergence-free functions. This projection is a slight modification of well-known projections with the property

where is an -projection; see Reference 5. Its implementation is very efficient as it can be computed in an element-by-element fashion.

Thus, given the convective velocity , the resulting scheme is nothing but the LDG method Reference 8 already studied for the Oseen equations Equation 1.3. Since that method is stable, high-order accurate and locally conservative, so is the LDG method under consideration. Moreover, the approximation to the velocity given by  has continuous normal components across elements and is globally divergence free in . To the knowledge of the authors, no other numerical scheme for the incompressible Navier-Stokes equations has all these properties.

Of course, the convective velocity depends on the discrete velocity field  through Equation 1.6 and, hence, we need to use an iterative method to compute it. To do so, we note that if is the LDG approximate velocity of the Oseen problem with convective velocity , then the approximate velocity  of the LDG method under consideration is a fixed point of . If is proven to be a contraction, to compute the approximate solution , we can use the fixed point iteration

This is nothing but a discrete version of the argument used to prove the existence and uniqueness of the exact solution of Equation 1.1. It ensures the existence and uniqueness of the exact solution of Equation 1.1 under a smallness condition of the type

where only depends on ; see Reference 17, Theorem 10.1.1 and the references therein. We mimic this argument to show that the approximate solution of the LDG method exists and is unique under a similar condition.

Let us point out that exact incompressibility can be achieved trivially if the LDG method has a velocity space that is div-conforming, i.e., it is included in . In this particular case, weak incompressibility implies exact incompressibility, provided that the discrete spaces are matched correctly. As will be discussed below, this approach can be viewed as a particular LDG method for which the operator  is chosen to be the identity. Consequently, all the results of this paper hold true verbatim for methods that are based on div-conforming velocity spaces. Furthermore, we note that, although we have used the LDG method to discretize the terms associated with the viscosity effects, any other DG discretization whose primal form is both coercive and continuous could have been used to that effect; see the discussions in Reference 2 and Reference 19.

The organization of the paper is as follows. In Section 2 we discuss the ideas that motivate the devising of the LDG method we propose in this paper. In Section 3 we present the LDG discretization in detail and verify its local conservativity. In Section 4 we state and discuss the main results, namely, the stability of the method, the convergence of the fixed point iteration, and the a priori error estimates, and in Section 5 we present their proofs. In Section 6 we present numerical experiments verifying the theoretical results. We end in Section 7 with some concluding remarks.

2. Devising the LDG method

In this section we discuss the ideas that led us to the devising of an LDG method that is both stable and locally conservative. To keep the discussion as simple and clear as possible, we do not work with the numerical method. Instead, we work directly with the equations Equation 1.1 and infer, from their structure, the properties of the corresponding LDG method.

2.1. A locally conservative LDG method

Since the incompressible Navier-Stokes equations Equation 1.1 are written in divergence form, a locally conservative LDG method can be easily constructed. However, it is very difficult to prove its stability. Let us illustrate this difficulty by using the equations for the exact solution.

If we multiply the first equation of Equation 1.1 by , integrate by parts and use the boundary conditions, we get

We see that we must use the incompressibility condition to obtain the equation

from which the stability estimate Equation 1.5 immediately follows.

In general, since exact incompressibility is very difficult to achieve after discretization, it is usually only enforced weakly. This weak incompressibility is enough, in a wide variety of cases, to guarantee that the discrete version of the term

is exactly zero, as for most mixed methods for the Stokes and Navier-Stokes equations, or nonnegative, as for the LDG methods considered for the Stokes Reference 9 and Oseen Reference 8 problems. Unfortunately, this is not true for the discrete version of the term

because the square of the modulus of the approximate velocity does not necessarily belong to the space of the approximate pressure.

2.2. The classical modification of the nonlinearity

A solution to this impasse can be obtained by using a now classical technique proposed back in the 1960s; see Reference 21 and Reference 22. From our perspective, it consists in modifying the nonlinearity of the equations as follows:

To see that this solves the problem, multiply the first equation by , integrate by parts, and use the boundary conditions to get

As a consequence, stability can follow from weak incompressibility. This is the case for most mixed methods for the Navier-Stokes equations. It is also the case for the first discontinuous Galerkin method for the incompressible Navier-Stokes equations Reference 13, a method which uses locally divergence-free polynomial approximations of the velocity, and for the more recent discontinuous Galerkin method developed in Reference 10.

The only problem with this approach is that local conservativity cannot be achieved because the first equation is not written in divergence form.

2.3. A new modification of the nonlinearity

To overcome this difficulty, it is enough to take another glance at the first equation in this section to realize that instead of working with the kinematic pressure , we should work with the new variable

If we incorporate this unorthodox pressure into the Navier-Stokes equations, we get

We now see that, since the above modification is in divergence form, locally conservative LDG methods can easily be constructed. Moreover, since we have

we also see that the stability of the LDG method can follow from weak incompressibility. The LDG method obtained with this approach can indeed be proven to have those properties; it is going to be studied thoroughly in a forthcoming paper.

2.4. Enforcing exact incompressibility

As we have seen, suitable modifications of the nonlinearity can be introduced which allow us to obtain stability by using only weakly incompressible approximations to the velocity. The approach we consider in this paper does not rely on a modification of that type. Instead, it is based on enforcing exact incompressibility in the space

The idea that allows this to happen is based on two observations. The first is that we can rewrite the Navier-Stokes equations as the Oseen problem

where, of course, . If we multiply the first equation by , integrate by parts, and use the boundary conditions, we get

This suggests that we consider an LDG method with two different (but strongly related) approximations to the velocity: one approximation for and another for . The stability for the LDG method would then be achieved if the approximation to  is weakly incompressible and if the approximation to is exactly incompressible. Further, local conservativity can be readily achieved for such an LDG method. Indeed, the fact that the equations are not written in conservative form can be compensated for by the fact that the approximation to is globally divergence free.

The second observation is that, if the approximation to given by an LDG method is weakly incompressible, it is possible to compute, in an element-by-element fashion, another approximation, , which is exactly incompressible. As we shall see, the post-processing operator is a slight modification of the well-known Brezzi-Douglas-Marini (BDM) interpolation operator; see Reference 4.

We note that the post-processing procedure can be omitted in the particular case where the velocity space is div-conforming. Indeed, in this case the fact that the approximation to the velocity is weakly incompressible does imply that it is exactly incompressible, provided the pressure space is chosen suitably. Hence, we can take .

We are now ready to describe the LDG method in detail.

3. The LDG method

In this section, we introduce a locally conservative LDG discretization for the Navier-Stokes equations Equation 1.1.

3.1. Meshes and trace operators

We begin by introducing some notation. We denote by a regular and shape-regular triangulation of mesh-size of the domain  into triangles . We further denote by  the set of all interior edges of  and by  the set of all boundary edges. We set .

Next we introduce notation associated with traces. Let and be two adjacent elements of . Let be an arbitrary point of the interior edge . Let be a piecewise smooth scalar-, vector-, or matrix-valued function and let us denote by the traces of on taken from within the interior of . Then, we define the mean value at as

Further, for a generic multiplication operator , we define the jump at as

Here, denotes an outward unit normal vector on the boundary of element . On boundary edges, we set accordingly , and , with denoting the outward unit normal vector on .

3.2. The LDG method for the Oseen equations

We now recall the LDG method for the Oseen equations Equation 1.3. We assume that the convective velocity field  is in the space

We begin by introducing the auxiliary variable and rewriting the Oseen equations as

Next, we introduce the space where

for an approximation order . Here denotes the space of polynomials of total degree at most  on . For simplicity we consider here only so-called mixed-order elements where the approximation degree in the pressure is of one order lower than the one in the velocity.

Finally, we define the approximate solution by requesting that for each ,

for all test functions . Each of the above equations is enforced locally, that is, element by element, due to the appearance of the so-called numerical fluxes , , , , and .

Thanks to this structure of the LDG method we immediately get that

and since is globally divergence free, we obtain a discrete version of the property of local conservativity Equation 1.4, namely,

In other words, the LDG method is locally conservative.

To ensure that the method is also stable (and high-order accurate), the numerical fluxes, which are nothing but discrete approximations to the traces on the boundary of the elements, must be defined carefully. As we shall prove, the numerical fluxes that define the LDG method for the Oseen equations Reference 8 do ensure stability. For the sake of clarity we consider the fluxes in their simplest form.

The convective numerical flux

For the convective flux in Equation 3.3, we take the standard upwind flux introduced in Reference 15Reference 18. For an element , we set

where is the inflow part of given by

The diffusive numerical fluxes

If a face lies inside the domain , we take

and, if lies on the boundary, we take

As will be shown later, the role of the parameter is to ensure the stability of the method; see also Reference 6.

The numerical fluxes related to the incompressibility constraint

The numerical fluxes associated with the incompressibility constraint, and , are defined by using an analogous recipe. If the face lies on the interior of , we take

On the boundary, we set

This completes the definition of the LDG method for the Oseen problem in Equation 3.1.

Remark 3.1.

Notice that one can take the div-conforming velocity space

while keeping the other spaces and definitions above unchanged.

3.3. The post-processing operator

To complete the definition of the LDG method for the Navier-Stokes equations Equation 1.1, it only remains to introduce what we refer to as the post-processing operator and to set in the approximation Equation 3.2Equation 3.4.

For a piecewise smooth velocity field , we define the operator by

where is the numerical flux Equation 3.8Equation 3.9 related to the incompressibility constraint. For each element, the local operator is given via the moments