A convergent FEMDG method for
the compressible NavierStokes equations
Abstract.
This paper presents a new numerical method for the compressible NavierStokes equations governing the flow of an ideal isentropic gas. To approximate the continuity equation, the method utilizes a discontinuous Galerkin discretization on piecewise constants and a basic upwind flux. For the momentum equation, the method is a new combined discontinuous Galerkin and finite element method approximating the velocity in the CrouzeixRaviart finite element space. While the diffusion operator is discretized in a standard fashion, the convection and timederivative are discretized using discontinuous Galerkin on the element average velocity and a LaxFriedrich type flux. Our main result is convergence of the method to a global weak solution as discretization parameters go to zero. The convergence analysis constitutes a numerical version of the existence analysis of Lions and Feireisl.
Key words and phrases:
ideal gas, isentropic, Compressible NavierStokes, Convergence, Compactness, Finite elements, Discontinuous Galerkin, Finite Volume, upwind, Weak convergence, Compensated compactness2010 Mathematics Subject Classification:
Primary: 35Q30,74S05; Secondary: 65M12Contents
 1 Introduction
 2 Preliminary material
 3 Numerical method and main result
 4 Energy and stability
 5 Estimates on the numerical operators
 6 Higher integrability on the density
 7 Weak convergence
 8 The discrete Laplace operator
 9 The effective viscous flux
 10 Strong convergence and proof of Theorem 3.5
 A Existence of a numerical solution
1. Introduction
In this paper, we will construct a convergent numerical method for the compressible NavierStokes equations:
(1.1)  
(1.2) 
where is an open, bounded, domain with Lipschitz boundary , and is a given finite final time. The unknowns are the fluid density and vector velocity , while the operator denotes the tensor product of two vectors (). The mechanism driving the flow is the pressure which is assumed to be that of an ideal isentropic gas (constant entropy);
To close the system of equations (1.1)  (1.2), standard noslip boundary condition is assumed
(1.3) 
together with initial data
(1.4) 
From the point of view of applications, the system (1.1)  (1.2) is the simplest form of the equations governing the flow of an ideal viscous and isentropic gas [1, 22]. In the available engineering literature, the reader can find a variaty of flows for which the assumption of constant entropy (reversibility) is a reasonable approximation. However, while viscosity in (1.2) is modeled by the Laplace operator (), practical applications will make use of a more appropriate stress tensor, the simplest being that of a Newtonian fluid with constant coefficients
Note that our diffusion is a special case of . Indeed, reduces to the Laplace operator when and . The analysis in this paper can be generalized to all relevant cases of constant nonvanishing and . However, this comes at the cost of more terms as the chosen finite element space (CrouzeixRaviart) does not satisfy Korn’s inequality [19], but with a negligible gain in terms of ideas and novelty. The more physically relevant case where and are functions of the density of the form is not included in the available existence theory (cf. [4, 13]) and there does not seem to be any obvious way of including it here either.
In terms of physical applicability of the forthcoming results, a more pressing issue is the equation of state for the pressure. For purely technical reasons, we will be forced to require that is greater than the spatial dimension;
This is a severe restriction on with no physical applications (to the author’s knowledge). Kinetic theory predicts a value of depending on the specific gas in question: for monoatomic gas (e.g helium), for diatomic gas (e.g air), and creeping towards one for more complex molecules and/or higher temperatures. It should be noted that global existence is only known for and hence only in the case of monoatomic gas (see [10] and the references therein). In this paper, the restriction on seems absolutely necessary to prove convergence of the method, but is not required for stability. In fact, the strict condition on is related to the numerical diffusion introduced by the method and is in perfect analogy to the problems encountered in a vanishing diffusion limit for the continuity equation (1.1) (see for instance [10]). That being said, it might be calming that the upcoming analysis can be easily extended to pressures of the form
where may be chosen arbitrarily small. Hence, for all practical purposes, the numerical method presented herein is convergent for cases very close to the physically relevant ones.
The numerical literature contains a vast body of methods for compressible fluid flows. Many of them are widely applied and constitute an indispensable tool in a variety of disciplines such as engineering, meteorology, or astrophysics. While the complexity and range of applicability of numerical methods for compressible viscous fluids is increasing, the convergence properties of any of these methods have thus far remained unknown. In fact, prior to this paper, there have not been reported any convergence results for numerical approximations of compressible viscous flow in more than one spatial dimension. In one dimension, the available results are due to David Hoff and collaborators [5, 26, 27, 28] (see also [13]) and concerns the equations posed in Lagrangian variables and relies on the 1D existence theory which have yet to be extended to multiple dimensions. That being said, in the recent years there have appeared a number of convergent methods for simplified versions of (1.1)(1.2). In [8, 9, 11], convergent finite volume and finite element methods for the stationary compressible Stokes equations was developed. Simultaneously, in [14, 15, 16], K. Karlsen and the author developed convergent finite element methods for the nonstationary compressible Stokes equations. In the upcoming analysis, we will utilize ideas from all of these recent papers.
Let us now discuss the choice of numerical method. For the approximation of the continuity equation, we will use the standard upwind finite volume method. However, we will formulate this method as a discontinuous Galerkin method where the density is approximated by piecewise constants. For the velocity we will use the CrouzeixRaviart finite element space. The method and some its properties was originally inspired by [25, 18] (cf. [14]), but variants can be found several places in the literature (see for instance [12]). For the momentum equation, we are leaning on the knowledge gained through [8, 9, 11, 14, 15, 16], from which it is evident that proving convergence for any numerical method is a hard task. In particular, to develop a numerical analogy of the continuous existence theory, it seems necessary that the numerical discretization of the diffusion and pressure respects orthogonal Hodge decompositions (see Section 8 for an explanation) in some form.
The distinctively new and completely original feature of the upcoming method is the discretization of the material derivative in the momentum equation. Our discretization will be of the form
where denotes the projection of onto piecewise constants. The operator is the specific upwind flux which we shall use (see Section 3 for precision). Hence, the material derivative is discretized using discontinuous Galerkin with the same polynomial order as the continuity discretization. This stands in contrast to the diffusion and pressure terms which are solved with the CrouzeixRaviart element space and hence in particular with piecewise constant divergence matching the density (and pressure) space. In the language of finite differences, this means that the pressure and diffusion is solved using staggered grid while the material derivative are solved at the same nodal values as the density. The great benefit of this approach is that it solves the longlasting problem of incorporating both the hyperbolic nature of the material derivative and the nature of the diffusionpressure coupling. In particular, our main result yields stability and convergence for all Mach and Reynolds numbers.
By proving convergence of a numerical method we will in effect also give an alternative existence proof for global weak solutions to the equations (1.1)(1.2). While the first global existence result for incompressible NavierStokes was achieved by Leray more than 80 years ago, a similar result for compressible flow was obtained by PL. Lions in the mid 90s. In the celebrated book [17], Lions obtains weak solutions of (1.1)(1.2) as the a.e weak limit of a sequence of approximate solutions. Consequently, from the point of view of analysis, we will in this paper perform similar analysis to that of [17], but for a numerical method. That is, we will develop a numerical analogy of the continuous existence theory. The key difficulty when passing to a limit is presented by the nonlinear pressure . Specifically, compactness of is needed, while the available estimates provides no form of continuity of . In all current proofs of existence, the necessary compactness is derived using renormalized solutions together with a remarkable sequential continuity result for the quantity . The result provides a.e convergence of the density, but gives no insight into continuity properties of . In the original proof, Lions needed that . The existence theory was further developed by Feireisl and the requirement lowered to . This seems to be optimal in the absence of pointwise bounds on the density. However, it is interesting that this still does not include the case of air in three dimensions. The reader is strongly encouraged to consult [23] for a thorough and wellwritten introduction to the mathematical theory of solutions to (1.1)(1.2).
Organization of the paper:
In the next section, we will go through some preliminary knowledge where we attempt to make clear the solution concept, basic compactness results, and the finite element theory used in the analysis. Then, in Section 3, we present the numerical method, give some basic properties of the method, and state the main convergence result (Theorem 3.5). We will then move on to Section 4 in which we establish stability of the method and draw conclusions in terms of uniform integrability. The following section, Section 5, is a fundamental section where we will estimate the weak error (in a weak norm) of the transport operators in the discretization. The material contained in this section will be used ubiquitously in the convergence analysis. In Section 6 we prove higher integrability of the density. That is, the density enjoys more integrability than the energy estimate provides. Then, in Section 7, we will pass to the limit in the method and conclude that the limit is almost a global weak solution. It will then only remain to prove strong convergence of the density approximation. In Section 9, we establish the fundamental ingredient in the proof of density compactness;  weak sequential continuity of the effective viscous flux. Finally, in Section 10 we prove strong convergence of the density and conclude the main result (Theorem 3.5). The paper ends with an appendix section containing the proof of wellposedness for the numerical method.
2. Preliminary material
The purpose of this section is to state some results that will be needed in the upcoming convergence analysis.
2.1. Weak and renormalized solutions
Definition 2.1.
Definition 2.2 (Renormalized solutions).
We shall need the following wellknown lemma [17] stating that squareintegrable weak solutions are also renormalized solutions.
2.2. Compactness results
In the analysis, we shall need a number of compactness results.
Lemma 2.4.
Let be a separable Banach space, and suppose , , is a sequence for which , for some constant independent of . Suppose the sequence , , is equicontinuous for every that belongs to a dense subset of . Then belongs to for every , and there exists a function such that along a subsequence as there holds in .
To obtain strong compactness of the density approximation, we will utilize the following lemma.
Lemma 2.5.
Let be a bounded open subset of , . Suppose is a lower semicontinuous convex function and is a sequence of functions on for which in , for each , in . Then a.e. on , , and . If, in addition, is strictly convex on an open interval and a.e. on , then, passing to a subsequence if necessary, for a.e. .
In what follows, we will often obtain a priori estimates for a sequence that we write as “” for some functional space . What this really means is that we have a bound on that is independent of . The following lemma is taken from [14].
Lemma 2.6.
Given and a small number , write with and . Let , be two sequences such that:

the mappings and are constant on each interval .

and converge weakly to and in and , respectively, where and

the discrete time derivative satisfies

as , uniformly in .
Then, in the sense of distributions on .
2.3. Finite element spaces and some basic properties
Let denote a shape regular tetrahedral mesh of . Let denote the set of faces in . We will approximate the density in the space of piecewise constants on and denote this space by . For the approximation of the velocity we use the Crouzeix–Raviart element space [6]:
(2.2) 
where denotes the jump across a face . To incorporate the boundary condition, we let the degrees of freedom of vanish at the boundary. Consequently, the finite element method is nonconforming in the sense that the velocity approximation space is not a subspace of the corresponding continuous space, .
For the purpose of analysis, we shall also need the divconforming Nédélec finite element space of first order and kind [20, 21]
(2.3) 
We introduce the canonical interpolation operators
(2.4) 
defined by
(2.5) 
Then, by virtue of (2.5) and Stokes’ theorem,
(2.6) 
for all . Here, and denote the curl and divergence operators, respectively, taken inside each element.
Let us now state some basic properties of the finite element spaces. We start by recalling from [3, 6, 20] a few interpolation error estimates.
Lemma 2.7.
There exists a constant , depending only on the shape regularity of , such that for any ,
for all and .
By scaling arguments, the trace theorem, and the Poincaré inequality, we obtain
Lemma 2.8.
For any and , we have

trace inequality,

Poincaré inequality,
In both estimates, is the diameter of the element .
Scaling arguments and the equivalence of finite dimensional norms yields the classical inverse estimate (cf. [2]):
Lemma 2.9.
There exists a positive constant , depending only on the shape regularity of , such that for and ,
for any and all polynomial functions , .
Since the CrouzeixRaviart element space is not a subspace of , it is not a priori clear that functions in this space are compact in , . However, from the Sobolev inequality on each element and an interpolation argument we obtain the needed result.
Lemma 2.10.
For , let satisfy and determine such that
The following space translation estimate holds
where is independent of and .
Proof.
From the work of Stummel [24], we have that
(2.7) 
The standard Sobolev inequality gives
(2.8) 
Hence, the proof follows from interpolation between (2.7) and (2.8).
∎
Finally, we recall the following wellknown property of the CrouzeixRaviart element space.
Lemma 2.11.
For any and ,
(2.9) 
Proof.
By direct calculation, we obtain
Now, since is linear on each element and is constant. Moreover, since the normal vector is constant on each face of the element, we have that
by definition of the interpolation operator . Hence, we have proved (2.9).
∎
3. Numerical method and main result
For a given timestep , we divide the time interval in terms of the points , , where we assume . To discretize space, we let be a shaperegular family of tetrahedral meshes of , where is the maximal diameter. It will be a standing assumption that and are related like for some constant . We also let denote the set of faces in .
The functions that are piecewise constant with respect to the elements of a mesh are denoted by and by we denote the Crouzeix–Raviart finite element space (2.2) formed on . To incorporate the boundary condition, we let the degrees of freedom of vanish at the boundary:
We will need some additional notation to accommodate discontinuous Galerkin discretization. Related to the boundary of an element , we write for the trace of the function taken within the element , and for the trace of from the outside. Related to a face shared between two elements and , we will write for the trace of within , and for the trace of within . Here and are defined such that points from to , where is fixed as one of the two possible normal components on . The jump of across the face is denoted .
To pose the method, and in the convergence analysis, we shall need the canonical interpolation operators (2.4). In fact, we shall need the operators and to such an extent that we introduce the convenient notation
(3.1) 
To discretize the convective operator in the continuity equation (1.1), we will utilize a standard upwind method in the degrees of freedom of . The upwind discretization is defined as follows
(3.2) 
where we have used the notation (3.1).
For the convective operator in the momentum equation (1.2) we will use the following LaxFriedrich type upwind flux
(3.3) 
Observe that this operator is posed for the average value of over each element. This is nonstandard and has to the author’s knowledge not been studied previously. We are now ready to pose the method.
  the mesh.  
  an element in the mesh.  
  the boundary of .  
  a face in the mesh.  
  all faces in the mesh.  
  the space of piece constant scalars on .  
  the CrouzeixRaviart vector space on .  
  the div conforming Nédélec space  
of first order and kind.  
  the projection operator onto .  
  the canonical interpolation operator onto .  
  the canonical interpolation operator onto .  
  (the piecewise constant projection).  
  (the Nédélec interpolation).  
  .  
  .  
  the trace of taken from within .  
  the trace of taken from outside .  
  the trace of taken against the normal vector .  
  the trace of taken in the direction of .  
  .  
  .  
  .  
  .  
  . 
Definition 3.1 (Numerical method).
Let and be given initial data and assume that is a given finite final time. Define the numerical initial data
(3.4) 
where is a small positive number. Determine sequentially
satisfying, for all ,
(3.5) 
and for all ,
(3.6) 
where should be chosen as small as possible.
For the purpose of analysis, we will need to extend the pointwiseintime numerical solution , , to a piecewise constant in time. For this purpose, we will use the following definition
(3.7) 
Remark 3.2.
3.1. The method is welldefined
Since the numerical method is nonlinear and implicit, it is not trivial that it is actually welldefined. In addition, the transport operators in the momentum equation is posed for the element average velocity and hence does not provide a full set of equations in themselves. In fact, it is only due to the diffusion that the system has a sufficient number of equations for the degrees of freedom of . We shall prove the existence of a numerical solution through a topological degree argument. The proof is very similar to that of [12, 14] and is for this reason deferred to the appendix.
Proposition 3.3.
In the upcoming analysis, we will need that the density solution is positive. For this purpose, we recall the following result from [14] (see also [12]):
Lemma 3.4.
Fix any and suppose , are given bounded functions. Then the solution of the discontinuous Galerkin method (3.5) satisfies
Consequently, if , then .
3.2. Main result
Our main result is that the numerical method converges to a weak solution of the compressible NavierStokes equations (1.1)  (1.4).
Theorem 3.5.
Suppose , is a given finite final time, and that the initial data satisfies
Let be a sequence of numerical solutions constructed according to Definition 3.1 and (3.7) with . Along a subsequence as ,
where is a weak solution of the isentropic compressible NavierStokes equations (1.1)  (1.2) in the sense of Definition 2.1.
4. Energy and stability
In this section we will prove that our method is stable and satisfies a numerical analogy of the continuous energy inequality (2.1). Both the stability estimate and large parts of the subsequent convergence analysis relies on a renormalized type identity derived from the continuity scheme (3.5). The proof of this identity can be found in [14, 15]. We shall utilize the following form:
Lemma 4.1.
Let solve the continuity scheme (3.5). Then, for any with , there holds
where and are some numbers in the range and , respectively.
We now prove our main stability result.
Proposition 4.2.
Proof.
Let in the momentum scheme (3.6), to obtain
(4.2) 
From Lemma 4.1 with , we see that the pressure term
(4.3) 
Next, we turn our attention to the first term after the equality in (4.2). From the definition of , we have that
(4.4)  
By setting in the continuity scheme (3.5), we see that the first term after the equality in (4) appears
(4.5)  
To see the contribution of the second term in (4), we first recall that
since the normal vector has opposite signs. Using this, we obtain
By applying this together with (4) in (4), we discover
(4.6) 