Discontinuous Galerkin method for Navier-Stokes equations
using kinetic flux vector splitting
Praveen Chandrashekar
TIFR Centre for Applicable Mathematics, Bangalore-560065, INDIA

arXiv:1209.4992v1 [cs.NA] 22 Sep 2012

Abstract
Kinetic schemes for compressible flow of gases are constructed by exploiting the connection between Boltzmann
equation and the Navier-Stokes equations. This connection allows us to construct a flux splitting for the NavierStokes equations based on the direction of molecular motion from which a numerical flux can be obtained. The
naive use of such a numerical flux function in a discontinuous Galerkin (DG) discretization leads to an unstable
scheme in the viscous dominated case. Stable schemes are constructed by adding additional terms either in
a symmetric or non-symmetric manner which are motivated by the DG schemes for elliptic equations. The
novelty of the present scheme is the use of kinetic fluxes to construct the stabilization terms. In the symmetric
case, interior penalty terms have to be added for stability and the resulting schemes give optimal convergence
rates in numerical experiments. The non-symmetric schemes lead to a cell energy/entropy inequality but exhibit
sub-optimal convergence rates. These properties are studied by applying the schemes to a scalar convectiondiffusion equation and the 1-D compressible Navier-Stokes equations. In the case of Navier-Stokes equations,
entropy variables are used to construct stable schemes.

1. Introduction
Kinetic schemes for the compressible Navier-Stokes equations are constructed by exploiting the connection
between the Boltzmann equation and the Navier-Stokes equations. The primary unknown in the Boltzmann
equation is the velocity distribution function, which gives the distribution of molecules in velocity and physical
space [1]. The macroscopic behaviour of the flow is obtained by averaging the Boltzmann equation over all the
molecular velocities. In kinetic schemes based on finite volume or discontinuous finite elements, the intercell
fluxes are approximated by exploiting this connection. Kinetic schemes were first developed for the solution
of Euler equations governing inviscid compressible flows by making use of the Maxwell-Boltzmann distribution
function, which is the equilibrium distribution for gases. While there are kinetic schemes based on a particle
model which update the velocity distribution function, our interest in the present work is the use of kinetic
model to construct numerical schemes at the continuum level. The kinetic model is utilized to derive a numerical
flux function by using the upwinding principle with respect to the molecular velocities. The availability of
such kinetic flux based schemes is useful for coupling the continuum solver with a particle-based method like
direct simulation Monte Carlo in the case of rarefied flows where a pure continuum model would not be valid
throughout the flow domain. The use of a DG scheme is expected to be helpful in coupling the continuum solver
with the particle solver. Kinetic schemes at the continuum level based on finite volume discretization have been
constructed in [2, 3]. In the approach of Mandal and Deshpande [4], the two states on either side of a cell face
are used to define two equilibrium distributions from which the inter-cell flux can be computed by taking the
appropriate moments. This gives a numerical flux function which depends smoothly on the two states and can
Email address: praveen@math.tifrbng.res.in (Praveen Chandrashekar)

Preprint submitted to Elsevier

September 25, 2012

be used in a method of lines approach; the resulting set of ODE are solved using any time integration scheme
like the Runge-Kutta method or an implicit scheme. Later, the Chapman-Enskog approximation to the BGK
model [5] of the Boltzmann equation was used by Chou and Baganoff [6] to derive a numerical flux function
for the Navier-Stokes equations. The usual practice in finite volume methods is to use a Godunov scheme or
approximate Riemann solver to construct the numerical flux function for the inviscid part of the flux and use
a central difference-type scheme for the viscous part of the flux. But the use of a kinetic formulation gives a
numerical flux function for the total flux, not just the inviscid flux. This numerical flux depends not only on
the two states at a cell face but also on their derivatives, which are present due to the diffusion terms.
An alternate approach is followed in the gas kinetic finite volume schemes [7, 8, 9] which are based on the
Boltzmann equation with BGK collision model. The exact solution of the BGK model is used to compute
the time integral of the inter-cell flux in a finite volume method. Thus in this approach, the time and space
discretizations are coupled, similar to the Lax-Wendroff scheme. The use of the BGK model automatically gives
a scheme for the Navier-Stokes equations. Ohwada [10] has analyzed the consistency of the gas kinetic scheme,
while a thorough analysis of the stability and accuracy of the BGK scheme for a scalar convection-diffusion
equation has been made in [11].
Discontinuous Galerkin finite element methods were first proposed in [12] for the neutron transport equation.
These methods use a discontinuous representation of the solution which allows them to compute solutions with
steep gradients and shocks. They are useful for solving hyperbolic problems like the Euler equations of inviscid
compressible flows and convection dominated problems like the compressible Navier-Stokes equations at high
Reynolds numbers. The inter-cell discontinuity in the state is resolved by making use of any consistent numerical
flux function for the convective flux. The vast developments in Godunov and approximate Riemann solvers can
be exploited to construct these numerical flux functions. For the case of hyperbolic problems like Euler equations,
the Runge-Kutta discontinuous Galerkin methods have been extensively developed [13, 14] and their high order
of accuracy is demonstrated in numerical experiments. DG methods provide a rational approach to construct
high order accurate schemes for convection-dominated problems and can be considered as a generalization of
finite volume methods, which have been very successful in computing discontinuous solutions and convection
dominated, high Reynolds number flows. They also have many computational advantages like a compact stencil
which aids parallelization, and the ability to easily develop hp-adaptive algorithms.
In the case of Navier-Stokes equations, the viscous fluxes which involve derivatives of the solution, must also
be computed across the element boundaries. However, the derivatives are also possibly discontinuous across
the elements due to the discontinuous numerical approximation. Averaging the derivatives to compute the
inter-element flux leads to an unstable scheme, and is a generic problem for a DG scheme applied to any partial
differential equation with elliptic terms. There are two approaches to construct stable schemes for the diffusion
terms, which are sometimes refered to as the flux formulation and the primal formulation [15]. In the first
approach, additional variables are introduced for the derivatives leading to a mixed finite element formulation.
Some numerical flux function is required for the additional variables whose choice is dictated by the stability of
the resulting scheme. Bassi and Rebay [16] developed such a scheme for compressible Navier-Stokes equations
and applied it to laminar flow problems. Cockburn and Shu [17] developed the local discontinuous Galerkin
(LDG) method for time dependent nonlinear convection-diffusion problems by generalizing the method of [16].
In the second approach, the derivatives are averaged for computing the diffusive fluxes but additional terms are
added to achieve stability of the scheme including interior penalty terms which help to enforce continuity of
the solution across the elements. The resulting schemes, known as interior penalty scheme, may be symmetric
(SIPG) or non-symmetric (NIPG) with respect to the diffusive terms. Arnold [18] proposed a symmetric DG
scheme with interior penalty for the heat equation and showed optimal error estimates in L2 norm. Baumann
and Oden [19, 20] introduced a non-symmetric scheme for diffusion and convection-diffusion problems, which
has been used for compressible Navier-Stokes equations in [21]. Symmetric schemes are preferable since they
are the only ones which have adjoint consistency [15], which is crucial for obtaining optimal error estimates in
L2 norm. A symmetric interior penalty scheme for compressible Navier-Stokes equations has been proposed
in [22] and its optimal accuracy has been shown through numerical experiments. A symmetric space-time DG

2

scheme for Navier-Stokes equations using the appoach of [16] and diffusive Riemann solvers for the viscous terms
has been developed in [23] and its optimal convergence rate has been demonstrated by numerical studies. The
various DG formulations for an elliptic equation have been unified in [15] and their stability, consistency and
adjoint consistency have been analyzed.
The use of a kinetic scheme has the attractive feature that it provides a numerical flux function for the
total flux which can be used in a DG scheme. The numerical flux depends on the derivatives of the solution
which are possibly discontinuous and a naive DG discretization can lead to an unstable scheme, especially in
the viscous dominated situation. Xu [24] has developed a DG discretization for the Navier-Stokes equations
using the BGK schemes, which however gives sub-optimal convergence rates similar to the NIPG scheme. Liu
and Xu [25] proposed a modified version of discontinuous Galerkin BGK scheme for Navier-Stokes equations
which shows improved convergence rates. These BGK-based DG schemes have not yet been analyzed from the
point of view of stability and adjoint consistent as in [15]. In the present work, we follow the interior penalty
approach and construct stable schemes by making use of the kinetic split fluxes to add additional stabilization
terms similar to the NIPG and SIPG schemes. The kinetic model is used to obtain a flux splitting from
which a numerical flux function is constructed even for the diffusion terms. The development of the schemes
is motivated by considering the scalar convection-diffusion equation for which a kinetic formulation can be
given. Through numerical experiments on this model problem, we find the NIPG scheme gives a convergence
rate of O(hk ) when using Pk basis functions which is sub-optimal as compared to the usual Galerkin methods,
while the SIPG scheme gives the optimal O(hk+1 ) convergence rate. Similar convergence rates are observed
through numerical experiments for the 1-D compressible Navier-Stokes equations. In the case of the NavierStokes equations, it is advantageous to use entropy variables since the equations assume a symmetric form in
those variables [26, 27, 28]. It also facilitates the construction of entropy stable schemes which is an important
thermodynamic constraint for compressible flows. Even the pure Galerkin scheme written in entropy variables
satisfies the entropy condition [27] for any degree of the basis functions. For Euler equations, Barth [29, 30] has
studied the entropy condition for discontinuous Galerkin methods and has derived a condition for a numerical
flux function to be entropy consistent, which is a generalization of the E-flux condition [31] to systems of
equations. The numerical flux resulting from kinetic flux vector splitting scheme [4] satisfies this entropy
condition. The entropy variable formulation is also advantageous since the flux jacobians are well behaved in
the limit of incompressible flows [32]; this is not the case with conserved variable formulation where some terms
in the flux jacobian matrix become unbounded or assume indeterminate form.
The rest of the paper is organized as follows. In Section (2) we first consider the scalar convection equation
for which a kinetic formulation is given together with a DG scheme. The energy stability of the DG scheme is
shown for the time continuous case. In Section (3), a kinetic formulation for the convection-diffusion equation is
given together with DG discretization. The instability of the naive DG discretization is shown through numerical
experiments. Stabilized versions of the schemes are constructed and their order of accuracy is studied through
numerical experiments. In Section (4), the Navier-Stokes equations and their entropy variables are discussed.
A novel DG scheme for the Navier-Stokes equations based on kinetic flux vector splitting is presented together
with some numerical results which demonstrate the accuracy of the scheme.
2. Kinetic model for convection equation
The scalar convection equation
ut + cux = 0

(1)

is the simplest hyperbolic partial differential equation. Throughout this work, we will assume periodic boundary conditions for the above partial differential equation. The exact solution of this equation is obtained by
convecting the initial condition at a constant speed c and the shape of the solution profile is unchanged. We
can obtain this equation as moments of a ”Maxwell-Boltzmann” distribution function
r


β
exp −β(v − c)2
(2)
g(v, u) = u
π
3

u+
j+1/2

j − 1/2

Ij

u−
j+1/2

Ij+1

j + 1/2

j + 3/2

Figure 1: Discontinuous finite element solution

since hgi = u and hvgi = cu, where the angle brackets denote integration over velocity space
Z +∞
h·i =
(·)dv
−∞

The quantity β > 0 is a free parameter that can be chosen arbitrarily1 . The ”Boltzmann equation” for the
distribution function is obtained by differentiating equation (2)
r


β
gt + vgx = (ut + vux )
exp −β(v − c)2
π
and we use equation (1) to replace ut to get
r
gt + vgx = (v − c)ux



β
exp −β(v − c)2 =: Q
π

(3)

Note that the moment of the collision term Q is zero, hQi = 0, so that taking moments of equation (3) gives us
the convection equation (1). The vanishing of the moments of the collision term means that molecular collisions
instantaneously drive the distribution function towards the local equilibrium distribution.
2.1. Kinetic numerical flux: KFVS
In a finite volume or DG method, we need to compute the flux across the inter-element boundaries. The
solution is possibly discontinuous across the elements. This defines two states u+
, u−
on either side of the
j+ 12
j+ 12
cell face xj+ 21 as shown in figure (1), where we have used the notation
u± (x) = lim u(x ∓ )
&0

In the kinetic model, molecules carry information due to their streaming motion. This allows us to use the
upwind principle to approximate the velocity distribution function at xj+ 12 based on the sign of the molecular
velocity as
(
g(v, u+
) v>0
j+ 21
gj+ 21 =
−
g(v, uj+ 1 ) v < 0
2

1 In the case of compressible gases for which kinetic theory holds, there is no free parameter in the Maxwell-Boltzmann velocity
distribution function. A parameter like β is present but is related to the absolute temperature.

4

The numerical flux function is obtained by integrating the above distribution function over the velocity space
leading to
, u−
) = hvgj+ 21 i = F + (u+
) + F − (u−
)
Fj+ 12 = F (u+
j+ 1
j+ 1
j+ 1
j+ 1
2

2

2

2

The split fluxes are obtained by integrating over half velocity spaces, and are given by
F ± (u) = cuA± + uB ± ,

A± =

1
[1 ± erf(s)],
2

1
B± = ± √
exp(−s2 ),
2 πβ

p
s=c β

The flux F + is due to all molecules moving to the right (v > 0) while the flux F − is due to all molecules moving
to the left (v < 0); this is the upwind principle for kinetic schemes. The numerical flux function F (u− , u+ ) is
consistent, i.e., F (u, u) = cu and is obviously a smooth function of u. The split flux F + (u) is an increasing
function of u while F − (u) is a decreasing function. Hence the numerical flux F (u− , u+ ) is a monotone flux. We
can also write the numerical flux function as a central flux with dissipation
F (u+ , u− ) =

1
1
(cu+ + cu− ) + D(u+ − u− ),
2
2

1
exp(−s2 )
D = c erf(s) + √
πβ

(4)

and it is easy to check that D > 0. In the limit of β → ∞, the split fluxes become
(
(
cu, c > 0
0,
c>0
+
−
F (u) =
,
F (u) =
0,
c<0
cu, c < 0
and the numerical flux function of equation (4) becomes the standard upwind flux
(
cu+ , c > 0
F (u+ , u− ) =
cu− , c < 0
We will see that the quantity D controls the rate of dissipation of energy and it is the smallest for the upwind
scheme. In all the computations we use a value of β = 1.
2.2. DG scheme
Consider a discretization of the one dimensional domain Ω by N elements Ij = [xj− 12 , xj+ 12 ], j = 1, . . . , N .
The DG method uses the broken space of polynomials of degree k

Vhk = φ ∈ L2 (Ω) : φ|Ij ∈ Pk (Ij )
The DG scheme is obtained by multiplying equation (1) by a test function φh from the broken space Vhk and
integrating over element Ij
Z
Z
φh ∂t uh dx −
cuh ∂x φh dx + F (uh )j+ 12 φh (x+
) − F (uh )j− 21 φh (x−
)=0
j+ 1
j− 1
Ij

2

Ij

2

The same equation would have been obtained if we had started from equation (3) by multiplying it with a
test function and perform the integrals over velocity space and physical space; the collision term would have
vanished due to the integration over velocity space. Note that we have used the numerical flux function to
resolve the discontinuity at cell boundaries, i.e.,
, u−
)
F (uh )j+ 21 = F (u+
j+ 1
j+ 1
2

2

If we take φh = uh and sum over all elements, we obtain the energy equation
N Z
N
X
X
1 d
2
kuh k −
cuh ∂x uh dx +
F (uh )j+ 12 Juh Kj+ 1 = 0
2
2 dt
j=1 Ij
j=1

5

(5)

where we have defined the inter-element jump in the solution by
Juh Kj+ 1 := uh (x+
) − uh (x−
)
j+ 1
j+ 1
2

2

2

After some simple calculations, the energy equation reduces to
N

1 d
1 X
2
kuh k2 + D
Juh Kj+ 1 = 0
2
2 dt
2 j=1
The jumps at the inter-element boundaries lead to dissipation of energy and the scheme is stable. Note that
even if D = 0, i.e., if we use a central flux, the scheme preserves the energy for any degree of the basis functions,
which is consistent with the exact solution of the partial differential equation.
3. Kinetic model for convection-diffusion equation
Consider the linear convection-diffusion equation
ut + cux = µuxx ,

µ>0

(6)

We can give a kinetic formulation for this equation using the Boltzmann equation with BGK model [11]
ft + vfx =

g−f
τR

(7)

where g is the same Maxwell-Botzmann distribution as used in the previous section. The basic idea of the BGK
model is that the distribution function f is close to the local equilibrium distribution g, and due to molecular
collisions, it relaxes towards g over a time interval τR which is the relaxation time. We can look for solutions to
the BGK model by using the Chapman-Enskog expansion in the relaxation time τR
f = f0 + τR f1 + τR2 f2 + . . .

(8)

Substituting this in equation (7) and collecting terms with the same coefficient of τR , we get
f0 = g,

f1 = −(gt + vgx )

as the first two terms in the expansion. The zeroth order term in the Chapman-Enskog expansion leads to the
convection equation as seen in the previous sections. The first order term becomes
r


β
exp −β(v − c)2
f1 = −(v − c)ux
π
where we have used the inviscid equation ut = −cux to replace the time derivative term since f1 should depend
only on the lower order terms in the Chapman-Enskog expansion. The first two terms constitute the ChapmanEnskog distribution function
r


β
ce
f = f0 + τR f1 = [u − τR (v − c)ux ]
exp −β(v − c)2
(9)
π
Upon taking moments of this distribution function, we obtain
hf ce i = u,

hvf ce i = cu −
6

τR
ux
2β

In order to recover the correct flux of equation (6), we have to choose τR = 2βµ so that hvf ce i = cu − µux .
Since f ce is a truncated solution in the Chapman-Enskog expansion it does not satisfy equation (7). We can
derive an equation for f ce by differentiating equation (9) as follows.
r


β
ce
ce
ft + vfx = [ut + vux − τR (v − c)(uxt + vuxx )]
exp −β(v − c)2
π
We now use the equation (6) to replace ut and uxt to obtain
r
ftce

+

vfxce



= (v − c)ux + (µ − τR (v − c)2 )uxx − τR µ(v − c)uxxx



β
exp −β(v − c)2 =: Q
π

Upon taking moments of this equation, the right hand side vanishes, i.e., hQi = 0, and we obtain the convectiondiffusion equation (6).
3.1. Flux vector splitting
The split fluxes for the viscous problem are obtained by integrating the Chapman-Enskog distribution
function f ce over the half velocity spaces, −∞ < v < 0 and 0 < v < +∞. The expressions for the kinetic split
fluxes can be written as
F ± = (cu − µux )A± + uB ±
The numerical flux function F (u+ , u− ) = F + (u+ )+F − (u− ) can be written as a sum of convective and dissipative
fluxes
−
F = Fc (u+ , u− ) + Fd (u+
x , ux )

where the convective flux Fc is identical to equation (4), while the dissipative flux is given by
−
+ +
− −
Fd (u+
x , ux ) = −µux A − µux A

Note that the diffusive flux can also be written in split form as
+ +
− −
−
Fd (u+
x , ux ) = Fd (ux ) + Fd (ux ),

Fd± (ux ) = −µux A±

In the case of the convection-diffusion equation the viscous fluxes depend only on the derivatives of the solution,
but for the Navier-Stokes equations, the viscous fluxes will depend on the solution also.
3.2. Test case
We will consider the convection-diffusion equation (6) in the domain (−1, +1) together with the initial
condition
u(x, 0) = − sin(πx)
(10)
and periodic boundary conditions. The exact solution to this problem is given by
u(x, t) = − exp(−µπ 2 t) sin(π(x − ct))

(11)

In all the tests, we set the convection speed to be c = 1. We define two test cases based on different values
of the viscosity coefficient µ as given in table (1) with tf being the final time. Test case 1 corresponds to a
situation in which convection dominates diffusion, while in Test case 2, diffusion dominates convection effects.
The time integration is performed using the 3-stage, strong stability preserving Runge-Kutta scheme of Shu
and Osher [33]. The time-step is chosen based on a CFL condition of the form ∆t = CFL · max(h/|c|, h2 /µ).
In the case of convection equation solved with a DG scheme, the CFL number depends on the degree of the
polynomial basis functions [34].
7

Test Case
µ
tf

1
0.001
30

2
1
0.5

Table 1: Parameters for the test case based on convection-diffusion equation (6)

3.3. DG scheme
The DG discretization of equation (6) is obtained by multiplying it by a test function φh and integrating
over any element Ij
Z
Z
Z
) − F (uh )j− 21 φh (x−
)=0
(12)
φh ∂t uh dx −
cuh ∂x φh dx + µ
∂x uh ∂x φh dx + F (uh )j+ 12 φh (x+
j+ 1
j− 1
Ij

Ij

2

Ij

2

In the above equation, the inter-element flux is approximated by the numerical flux which is defined as
−
+
−
F (uh ) = Fc (u+
h , uh ) + Fd (∂x uh , ∂x uh )

and it includes both convective and diffusive fluxes. Substituting φh = uh and summing over all elements we
obtain
N
N Z
N Z
X
X
X
1 d
F (uh )j+ 21 Juh Kj+ 1 = 0
(∂x uh )2 dx +
cuh ∂x uh dx + µ
kuh k2 −
2
2 dt
j=1
j=1 Ij
j=1 Ij
Following similar steps as in the case of linear convection equation, we obtain the following energy equation
N
N
N Z
X
X
1 d
1 X
2
2
Fd (∂x uh )j+ 21 Juh Kj+ 1 = 0
Juh Kj+ 1 + µ
(∂x uh )2 dx +
kuh k + D
2
2
2 dt
2 j=1
j=1
j=1 Ij

(13)

The second and third terms are dissipative since they lead to a decrease of the energy. But we do not know
the sign of the last term, which can be either positive or negative. We apply this DG scheme to the two test
problems using P1 basis functions and the results are shown in figure (2). When µ = 0.001 the results look
correct but this is not so when µ = 1 which corresponds to a viscous dominated case. The scheme is not able
to dissipate energy at the correct rate. We see that a naive discretization of convection-diffusion equation by
the discontinuous Galerkin method leads to unstable schemes when viscosity is dominant.
3.4. DG scheme with stabilization
Motivated by the energy analysis, we propose the following stabilized DG scheme
Z
Z
Z
∂x uh ∂x φh dx + F (uh )j+ 12 φh (x+
φh ∂t uh dx −
cuh ∂x φh dx + µ
) − F (uh )j− 21 φh (x−
)
j+ 1
j− 1
Ij

−

Ij
Fd+ (∂x φ+
h )j+ 12

2

Ij

Juh Kj+ 1 −
2

Fd− (∂x φ−
h )j− 21

2

(14)

Juh Kj− 1 = 0
2

Note that the last two terms in the above equation are constructed from the kinetic split viscous fluxes. Choosing
φh = uh and summing up over all elements gives us
N
N Z
X
1 d
1 X
2
kuh k2 + D
Juh Kj+ 1 + µ
(∂x uh )2 dx = 0
2
2 dt
2 j=1
I
j
j=1

(15)

which shows that the scheme is stable. The choice of the additional stabilization terms in equation (14) was
just enough to eliminate the problematic term in the energy equation (13). The stabilized scheme is applied to
8

1
Exact
KFVS

Exact
KFVS
0.5

0.5

0

0

-0.5

-1
-1

-0.5

-0.5

0

1

0.5

-1

0

-0.5

(a)

0.5

1

(b)

Figure 2: Convection-diffusion problem using equation (12) and P1 basis functions. (a) Test case 1, (b) Test case 2

convection-diffusion problem with µ = 1 and the results are shown in figure 3. The stabilized scheme is able to
dissipate the energy better than the unstabilized one, but it exhibits a significant amount of phase error for P1
basis functions; the energy dissipation is still seen to be incorrect with P1 basis functions. However, the use of
P2 basis functions is seen to give excellent results without any phase error. In the next sections, we see that
adding a penalty term gives correct results even with P1 basis functions.
3.5. Cell energy analysis
In the previous sections, we performed a global energy analysis; here we make a cell energy analysis. If we
multiply the convection-diffusion equation by u and integrate over element Ij , we obtain

j+ 12
Z
1
d
j+ 1
kuk2Ij + cu2
+ [−µuux ]j− 21 = −µ
(ux )2 dx ≤ 0
2
dt
2
1
Ij
j−

(16)

2

We will now try to mimic this property for the DG scheme for which we will have to define consistent numerical
fluxes for the second and third terms on the left hand side of the above equation. Substituting φh = uh in
equation (14), we obtain
i
d
ch 2 +
+
kuh k2Ij −
uh (xj+ 1 ) − u2h (x−
) − F (uh )j− 12 uh (x−
)
1 ) + F (uh )j+ 1 uh (x
j− 2
j+ 12
j− 12
2
2
dt
2
Z
−
−
− Fd+ (∂x u+
h )j+ 12 Juh Kj+ 1 − Fd (∂x uh )j− 12 Juh Kj− 1 = −µ
2

2

(∂x uh )2 dx

Ij

The numerical flux F (uh )j+ 21 consists of the convective and viscous fluxes. Define the average value
{{uh }}j+ 1 :=
2

1
(uh (x+
) + u(x−
))
j+ 21
j+ 12
2

and convective and diffusive numerical fluxes
F̄c (uh )j+ 12

=

F̄d (uh )j+ 12

=

c  2
uh j+ 1
2
2
+
+
−
−
+
Fd (∂x uh )j+ 12 uh (xj+ 1 ) + Fd (∂x u−
h )j+ 21 uh (xj+ 1 )
Fc (uh )j+ 12 {{uh }}j+ 1 −
2

2

9

2

0.4
Exact
KFVS

0.4

Exact
KFVS
0.2

0.2

0

0

-0.2
-0.2

-0.4
-1

-0.5

0

-0.4
-1

1

0.5

0

-0.5

(a)

0.5

1

(b)

Figure 3: Results for Test case 2 using equation (14): (a) P1 and (b) P2 basis functions

The last two equations provide consistent numerical fluxes for 12 cu2 and −µuux , respectively. Then we obtain
the cell energy equation for the DG scheme as
Z
 j+ 1  j+ 1
1
1
d
2
2
kuk2Ij + F̄c j− 21 + F̄d j− 21 = −µ
(∂x uh )2 dx − D Juh Kj− 1 − D Juh Kj+ 1 ≤ 0
2
2
2
2
dt
4
4
IJ
which is consistent with the exact energy equation (16). The inter-element jumps in the solution add additional
dissipation in the numerical solution.
3.6. DG scheme with interior penalty
The above examples have shown the importance of energy stability in constructing good DG scheme for the
convection-diffusion equation. The schemes were written for each element but it is more instructive to write
a weak formulation for the whole computational domain since it helps to understand some structure of the
equations. By summing up the equation (14) for each element, we write the DG scheme for convection-diffusion
equation as a weak statement: Find uh (t) ∈ Vhk such that for all φh ∈ Vhk , the following equation is satisfied
Z
Z
Z
X
φh ∂t uh dx −
cuh ∂x φh dx + µ
∂x uh ∂x φh dx +
Fc (uh )j+ 12 Jφh Kj+ 1
Ω

Ω

+

X

2

Ω

j

Fd (∂x uh )j+ 12 Jφh Kj+ 1 + 

X

2

j

Fd (∂x φh )j+ 21 Juh Kj+ 1 +

X

2

j

(17)

δh (uh )j+ 21 Jφh Kj+ 1 = 0
2

j

where the last term is an additional term that has been added, called the interior penalty term, and is of the
form [18, 15]
µ
δh (u) = Cip JuK
h
where Cip is a non-dimensional positive number. Two classes of schemes are obtained based on whether  = −1
or +1. Define the bilinear form consisting of all diffusive terms
X
X
X
a (vh , wh ) =
Fd (∂x vh )j+ 12 Jwh Kj+ 1 + 
Fd (∂x wh )j+ 12 Jvh Kj+ 1 +
δh (vh )j+ 21 Jwh Kj+ 1
2

j

2

j

2

j

10

1. If we choose  = −1 in (17) and omit the interior penalty term, then we obtain the scheme of equation (14).
The scheme is non-symmetric since a−1 (vh , wh ) 6= a−1 (wh , vh ). As we have seen, this scheme satisfies a
cell energy inequality. The interior penalty terms are necessary to make the scheme stable since
Z
a−1 (vh , vh ) = µ

(∂x vh )2 dx +

Ω

N
X
j=1

Cip

µ
2
Jvh K
h

and the right handside is a norm on the space of discontinuous functions only in the presence of the penalty
term for any Cip > 0. Without the penalty term (Cip = 0), the right hand side is only a semi-norm and
the scheme is only weakly stable [15]. In this case, the discretization of diffusive terms is similar to that
of Baumann and Oden [20] for elliptic equation, who also show sub-optimal convergence in the case k ≥ 2
only. There is no proof of convergence for k = 1 case which combined with only weak stability explains
the poor performance of the scheme shown in figure (3a). With the interior penalty terms, the global
energy equation is
N
N Z
N
X
X
1 X
1 d
µ
2
2
2
kuh k + D
Juh Kj+ 1 + µ
(∂x uh )2 dx +
Cip Juh K = 0
2
2 dt
2 j=1
h
j=1 Ij
j=1

which shows that interior penalty adds additional stabilization to the scheme which is important when
viscosity is dominant. The scheme with interior penalty is known as NIPG and is stable for any k ≥ 1.
2. If  = +1, the discretization of the diffusive terms becomes symmetric, i.e., a1 (vh , wh ) = a1 (wh , vh ), and
the scheme is known as SIPG scheme. However, the SIPG scheme does not satisy a cell energy inequality,
though in the case of elliptic equation it can be shown to be globally stable [15]. This scheme is stable
only with the addition of a penalty term δh , with a sufficiently large value of Cip . The determination
of a good value of Cip is an open problem; in practice a value of Cip between 5 to 10 seems to yield
satisfactory results. A larger value of Cip reduces the allowable time step for explicit time integration
schemes, requiring the necessity of implicit time integration schemes.
The jump terms in the above DG scheme do not destroy the consistency of the scheme; if we substitute the
exact solution of the convection-diffusion equation which is continuous in equation (17), it exactly satisfies the
equation since the jump terms vanish.
3.7. Numerical results using interior penalty scheme
In this section, we present results for the two test cases, using the interior penalty scheme given by equation (17). We first present results using 20 cells. The schemes NIPG and SIPG are tested with a penalty
parameter Cip = 10. Figure (4) shows the results for Test case 1 which is convection dominated situation.
Figure (5) shows results for Test case 2 which is a diffusion dominated problem. We see that the addition of
the interior penalty leads to better solutions with P1 basis functions, as compared to figure (3).
In order to compute the convergence rates of the various schemes, we perform a grid refinement study.
Errors in the solution and its derivative are computed in L2 norm and shown in tables (2-9) for P1 and P2 basis
functions. We observe the following trends:
• SIPG gives more accurate results than NIPG. This can be seen by comparing the error values for the two
schemes given in the tables.
• NIPG converges at sub-optimal rates; both the solution and its derivative converge at O(hk ).
• SIPG converges at optimal rates; the solution converges at O(hk+1 ) and the derivative converges at O(hk ).

11

1

1
Exact
KFVS

Exact
KFVS

0.5

0.5

0

0

-0.5

-0.5

-1
-1

-0.5

0

-1
-1

1

0.5

-0.5

NIPG, P1

0

1

1
Exact
KFVS

Exact
KFVS

0.5

0.5

0

0

-0.5

-0.5

-1
-1

1

0.5

NIPG, P2

-0.5

0

0.5

-1
-1

1

SIPG, P1

-0.5

0

0.5

SIPG, P2

Figure 4: Results for Test case 1 using interior penalty scheme

12

1

0.4

0.4
Exact
KFVS

Exact
KFVS

0.2

0.2

0

0

-0.2

-0.2

-0.4
-1

-0.5

0

-0.4
-1

1

0.5

-0.5

NIPG, P1

0

0.4

0.4
Exact
KFVS

Exact
KFVS

0.2

0.2

0

0

-0.2

-0.2

-0.4
-1

1

0.5

NIPG, P2

-0.5

0

0.5

-0.4
-1

1

SIPG, P1

-0.5

0

0.5

SIPG, P2

Figure 5: Results for Test case 2 using interior penalty scheme

13

1

cells
20
40
80
160
320

dofs
40
80
160
320
640

L2
1.259e-02
3.535e-03
1.245e-03
3.511e-04
3.986e-05

1.83
1.51
1.83
3.14

H1
2.757e-01
1.356e-01
6.639e-02
3.257e-02
1.614e-02

1.02
1.03
1.03
1.01

Table 2: Test case 1: NIPG, P1 basis

cells
20
40
80
160
320

dofs
60
120
240
480
960

L2
3.376e-04
7.361e-05
1.660e-05
3.433e-06
6.892e-07

2.20
2.15
2.27
2.32

H1
1.416e-02
3.612e-03
9.186e-04
2.243e-04
5.103e-05

1.97
1.98
2.03
2.14

Table 3: Test case 1: NIPG, P2 basis

Actually both the schemes seem to exhibit slightly higher convergence rates on Test case 1 which is convection
dominated. In Test case 2 which is viscosity dominated, the convergence rates are very close to the above
mentioned values. These convergence rates are consistent with what is observed in the case of DG schemes
applied to elliptic equation [15]. The NIPG scheme for elliptic problems is numerically found to show optimal
convergence rate for odd degree polynomials but we have not observed this phenomenon for the convectiondiffusion equation with the current NIPG scheme.
3.8. Test case 3
Torrilhon and Xu [11] have analyzed the BGK finite volume scheme for scalar convection-diffusion equation
and find that the BGK scheme exhibits only first order accuracy in the asymptotic limit, though on coarse grids
it shows third order accuracy. They find a non-monotone convergence of the error with grid size. In order to
study the behaviour of the error convergence, we apply the present DG scheme on a test case similar to the one
given in [11] with β = 1, µ = 0.005, c = 1, using P1 basis functions. The grid is refined starting with 50 cells
and increased by 50 cells in each refinement step. The initial condition is taken to be
u(x, 0) = 4 +

8
16
sin(πx/2) +
sin(3πx/2)
π
3π

for which an exact solution is available [11] for the case of periodic boundary conditions using Fourier approach.
The convergence of the L2 error in the solution as a function of the number of cells is shown in figure (6a)
cells
20
40
80
160
320

dofs
40
80
160
320
640

L2
1.245e-02
1.843e-03
3.170e-04
6.181e-05
1.351e-05

2.76
2.54
2.36
2.19

H1
2.866e-01
1.403e-01
6.807e-02
3.293e-02
1.615e-02

Table 4: Test case 1: SIPG, P1 basis

14

1.03
1.04
1.05
1.03

cells
20
40
80
160
320

dofs
60
120
240
480
960

L2
1.354e-04
1.622e-05
1.900e-06
2.153e-07
2.487e-08

3.06
3.09
3.14
3.11

H1
1.367e-02
3.296e-03
7.692e-04
1.748e-04
4.324e-05

2.05
2.10
2.14
2.02

Table 5: Test case 1: SIPG, P2 basis

cells
20
40
80
160
320

dofs
40
80
160
320
640

L2
1.592e-02
8.018e-03
4.031e-03
2.022e-03
1.013e-03

0.99
0.99
1.00
1.00

H1
1.088e-01
5.451e-02
2.727e-02
1.364e-02
6.820e-03

1.00
1.00
1.00
1.00

Table 6: Test case 2: NIPG, P1 basis

cells
20
40
80
160
320

dofs
60
120
240
480
960

L2
7.467e-04
1.860e-04
4.648e-05
1.162e-05
2.906e-06

2.01
2.00
2.00
2.00

H1
4.778e-03
1.195e-03
2.989e-04
7.472e-05
1.868e-05

2.00
2.00
2.00
2.00

Table 7: Test case 2: NIPG, P2 basis

cells
20
40
80
160
320

dofs
40
80
160
320
640

L2
2.255e-03
5.653e-04
1.415e-04
3.538e-05
8.846e-06

2.00
2.00
2.00
2.00

H1
1.124e-01
5.638e-02
2.823e-02
1.412e-02
7.063e-03

1.00
1.00
1.00
1.00

Table 8: Test case 2: SIPG, P1 basis

cells
20
40
80
160
320

dofs
60
120
240
480
960

L2
1.355e-04
1.713e-05
2.158e-06
2.709e-07
3.394e-08

2.98
2.99
2.99
3.00

H1
9.936e-03
2.497e-03
6.268e-04
1.571e-04
3.932e-05

Table 9: Test case 2: SIPG, P2 basis

15

1.99
1.99
2.00
2.00

−2

NIPG
SIPG

10

NIPG
SIPG
−4

10

||u−u ||

Slope = 1

h

||u−uh||

−3

10

−5

10

Slope = 2

−6

10
−4

10

−7

10

Slope = 2
2

3

10

10
log(N)

Slope = 3
2

3

10

10
log(N)

(a)

(b)

Figure 6: Test case 3; convergence of error for (a) P1 and (b) P2 basis functions

for P1 basis functions; this plot must be compared with figure (7) in [11]. It is clear that both the NIPG and
SIPG schemes converge monotonically with grid refinement. Asymptotically, the NIPG scheme converges at
first order while the SIPG scheme converges at second order. Surprisingly, the NIPG scheme exhibits second
order convergence on coarse grids, a superconvergence phenomenon that is also observed in the BGK finite
volume scheme [11]. The convergence of the error for P2 basis functions is shown in figure (6b) which shows
second and third order convergence for NIPG and SIPG schemes, respectively. On coarser grids, the NIPG
scheme again shows optimal convergence rate. These convergence rates are consistent with what is observed in
the previous two test cases, and also in the 1-D Navier-Stokes equations given in the later sections.
4. Navier-Stokes equations
The Navier-Stokes equations which represent conservation laws for mass, momentum and energy for a fluid
can be written as a system of partial differential equations
∂t U + ∂i Fi (U ) + ∂i Gi (U, ∇U ) = 0
where U is the vector of conserved variables corresponding to density of mass, momentum, and energy, while
Fi , i = 1, 2, 3 are the inviscid or convective fluxes and Gi are the viscous fluxes. The conserved variables are
given by

>
U = ρ ρu1 ρu2 ρu3 ρe
where ρ is the density, u = (u1 , u2 , u3 ) is the fluid velocity and e is the energy per unit mass, which can be
written in terms of the internal energy ε and the kinetic energy as e = ε + 21 u2 . An equation of state ε = ε(ρ, p)
is required to close the system of equations, where p is the pressure; for a polytropic ideal gas, the equation
C
p
of state is given by ε = ρ(γ−1)
where γ = Cvp is the ratio of specific heats at constant pressure and constant

16

volume [35]. The flux vectors are given by





0
 −τi1




−τ
Gi = 
i2


 −τi3

−τij uj + qi


ρui
pδi1 + ρui u1 



Fi = 
pδi2 + ρui u2  ,
pδi3 + ρui u3 
(ρe + p)ui

(18)

The shear stress tensor τ is related to the strain rate by Newton’s constitutive law while the heat flux is related
to the gradient of absolute temperature T by Fourier’s law
2
τij = µ(∂i uj + ∂j ui ) − µ(∂k uk )δij ,
3

qi = −κ∂i T

(19)

The coefficient of dynamic viscosity µ and the coefficient of heat conduction κ are related through the Prandtl
µC
number Pr = κ p . We will use the notation F = [F1 , F2 , F3 ] and G = [G1 , G2, G3 ]. Moreover, for any vector
n ∈ R3 , we denote F · n = Fi ni and G · n = Gi ni for the fluxes in the direction n. Since the numerical
experiments in this paper are conducted for the 1-D Navier-Stokes equations, the corresponding equations are
given in Appendix A.
4.1. Euler equations, entropy variables and symmetric hyperbolic form
The system of equations containing only the convective terms
∂t U + ∂i Fi = 0

(20)

are called the Euler equations, and are hyperbolic in nature. This means that for any vector n = (n1 , n2 , n3 ),
the matrices
A(U, n) = Ai (U )ni ,
Ai (U ) = Fi0 (U )
have all real eigenvalues and a complete set of eigenvectors. The Euler equations admit an entropy function
η(U ) with entropy fluxes θi (U ), i = 1, 2, 3 such that2
θi0 (U ) = η 0 (U )Fi0 (U )
This implies that for smooth solutions, an additional conservation law is satisfied, i.e., ∂t η + ∂i θi = 0, while for
discontinuous solutions the inequality ∂t η + ∂i θi ≤ 0 is satisfied in the weak sense. The entropy and entropy
fluxes can be related to the physical entropy s = ln(p/ργ ) + s0 by
η=−

ρs
,
γ−1

θi = −

ρsui
γ−1

(21)

Actually any function of the form η = ρh(s) where h is convex is an entropy for the Euler equations, but the
above linear choice is the only one which is consistent with the second law of thermodynamics in the presence of
heat conduction [27]. Due to the mass conservation equation, the constant term s0 can be ignored and moreover
it can be shown that η(U ) is a strictly convex function. We can then define the entropy variables by
V > = η 0 (U )

(22)

which can be inverted in a unique manner to obtain U = U (V ). Let us define the dual variables by a Legendre
transform
η ∗ (V ) = V · U (V ) − η(U (V )),
θi∗ (V ) = V · Fi (U (V )) − θi (U (V ))
2 For

a scalar function of a vector like η(U ) we consider η 0 (U ) to be a row vector.

17

0

0

Differentiating these expressions it is easy to check that η ∗ (V ) = U > and θi∗ (V ) = Fi (U (V ))> . It then follows
00
00
that the matrices U 0 (V ) = η ∗ (V ) and Fi0 (U (V ))U 0 (V ) = θi∗ (V ) are symmetric. Moreover the matrix
U 0 (V ) = η 00 (U (V ))−1
is positive-definite. Making a change of variables from U to V in equation (20), we get
Ã0 ∂t V + Ãi ∂i V = 0
where we have defined

00

00

Ã0 = U 0 (V ) = η ∗ (V ),

Ãi = Ai Ã0 = θi∗ (V )

the matrix Ã0 is symmetric positive definite, while the matrices Ãi are symmetric.
4.2. Kinetic description of Euler equations
The Boltzmann equation [1] gives a molecular statistical description of gas dynamics in the form of the
velocity distribution function f (v, x, t) which can be written as
∂t f + vi ∂i f = Q(f, f )
where the right hand side represents the effect of intermolecular collisions. The meaning of the distribution
function is that the quantity f (v, x, t)dxdv gives the average number of particles found in the spatial volume
dx around the point x at time t and having velocities in the range dv around the velocity v. Thus f gives the
number density of molecules in phase space and by a suitable scaling, it can be considered as providing the
mass density. The macrosopic fluid motion is the average effect of all molecular motions, so that the conserved
variables are obtained by integrating the distribution function over the velocity space
Z
U=
ψ(v)f (v, x, t)dv = hψf i
R3

where
ψ(v) =



1

v1

v2

v3

1
2
2 |v|

>

are the microscopic collisional invariants. Since collisions do not destroy mass, momentum and energy, we have
hψQi = 0. In fact, it can bePshown [1] that any function φ(v) of the molecular velocities for which hφQi = 0
must be of the form φ(v) = ai ψi .
A gas in thermodynamic equilibrium is characterized by the fact that there are no spatial or temporal
variations in the distribution function, i.e., Q(f, f ) = 0. It can be shown that the unique solution g = f of this
equation is of the form
g = exp(V · ψ)

for some constants V ∈ R5 . This is the usual Maxwell-Boltzmann distribution function and can be written in
the more familiar form
 3/2


β
1
g(v) = ρ
exp −β|v − u|2 ,
β=
π
2RT
The Maxwell distribution can also be characterized as the solution of the following minimization problem [36]
min {H[f ] : hψf i = U } ,

H[f ] = hf ln f i

and the quantities V appear as the Lagrange multipliers to enforce the constraint. If we define
η ∗ (V ) = hf i = hexp(V · ψ)i
then

0

η ∗ (V ) = hψ > exp(V · ψ)i = U >

(23)

and hence the Lagrange multipliers V are precisely the entropy variables corresponding to the macroscopic state
U . The modifications in the case of polyatomic gases are briefly discussed in Appendix B.
18

K′
n

Vh+

Vh−

K

Figure 7: Definition of edge normal and trace values

4.3. From Maxwell distribution to Euler equations
The Maxwell distribution was obtained under the condition of thermodynamic equilibrium corresponding to
a gas which has uniform properties in space and time. The appropriate moments of the distribution function
yield the conserved variables and the inviscid fluxes, i.e., U = hψ(v)gi and Fi = hvi ψ(v)gi. Due to this
correspondence, it is common to use the collision-less Boltzmann equation
∂t g + vi ∂i g = 0
to construct numerical schemes for the Euler equations. When the gas is not uniform, then the assumption is
that locally, it is still in thermodynamic equilibrium, so that one has a local Maxwell distribution
g(v, x, t) = exp(V (x, t) · ψ)
defined at every point of space. This distribution however does not satisfy the collision-less Boltzmann equation.
We can derive an equation for g by differentiating the above equation
h
i
∂t g + vi ∂i g = g [∂t V + vi ∂i V ] · ψ = g −Ã−1
Ã
+
v
I
∂i V · ψ
i
i
0
Taking moments on both sides, we obtain
hψ(∂t g + vi ∂i g)i = −hψ ⊗ ψgiÃ−1
0 Ãi ∂i V + hvi ψ ⊗ ψgi∂i V

(24)

Now from equation (23)
00

η ∗ (V ) = hψ ⊗ ψgi = U 0 (V ) = Ã0

Moreover, since Fi = hvi ψgi, we get Fi0 (V ) = hvi ψ ⊗ ψgi = Ãi , which shows that the right hand side of equation
(24) vanishes, and we obtain the Euler equations. Due to this reason the right hand side is ignored and we
use the collision-less Boltzmann equation with Maxwell-Boltzmann distribution function to derive numerical
schemes for Euler equations.
4.4. DG scheme for Euler equations
Consider a triangulation Th = {K} of the domain Ω by polygonal, non-overlapping elements. We assume
that the edges of the triangulation are oriented, so that each edge has a unit normal vector n with trace values

19

of the discontinuous functions on the edge defined as in figure (7). Let us also define the jump and average
operators as
1
J·K = (·)+ − (·)− ,
{{·}} = [(·)+ + (·)− ]
2
respectively. Note that we do not involve the normal vector n in the definition of these operators since it is used
in the definition of the numerical fluxes. The DG scheme for an interior element K is given by multiplying the
Euler equation by a test function Wh and doing some integration by parts to obtain
Z
Z
Z
Wh · ∂t U (Vh )dx −
Fi (Vh ) · ∂i Wh dx +
Hc (Vh+ , Vh− , n) · Wh+ dσ = 0
(25)
K

K

∂K

R
The integral
on the element boundary ∂K is to be interpreted as the sum of the edge integrals, i.e., K (·)dσ =
R
P
e∈∂K e (·)dσ. Note that the flux across the element boundaries is approximated by a numerical flux function.
The numerical flux is given by the concept of upwinding at the kinetic level; since molecules carry mass,
momentum and energy, the flux across any edge is determined as a contribution due to molecules crossing the
edge in both directions, which leads to the following numerical flux function
Hc (Vh+ , Vh− , n) = (v · n)+ ψg(Vh+ ) + (v · n)− ψg(Vh− )

(26)

where we have used the notation

1
[(v · n) ± |v · n|]
2
Note that each of the integrals in the two terms of the numerical flux are effectively taken over half the velocity
space, i.e, v · n > 0 and v · n < 0, respectively. The kinetic flux can also be written as
(v · n)± =



1
1
1
(v · n)ψg(Vh+ ) +
(v · n)ψg(Vh− ) +
|v · n|ψ g(Vh+ ) − g(Vh− )
2
2
2

Hc (Vh+ , Vh− , n) =
and using
g(Vh+ ) − g(Vh− ) =

1

Z
0

d
g(V̄ (s))ds =
ds

Z
0

1

g(V̄ (s))ψ · JVh K ds,

V̄ (s) = sVh+ + (1 − s)Vh−

we can write the numerical flux as
Hc (Vh+ , Vh− , n)


1
1
F (Vh+ ) + F (Vh− ) · n +
=
2
2

Z
0

1

|v · n|ψ ⊗ ψg(V̄ (s)) · JVh K ds

This is precisely the kinetic mean-value (KMV) flux of Barth [37, 30]. Hence for the Euler equations, the
KMV flux is identical to the KFVS flux, which can be written in terms of explicit formulae but involving error
functions and exponentials [4]. The expressions for the one dimensional split fluxes are given in the Appendix
C.
4.4.1. Entropy stability of DG scheme
To study entropy stability of the scheme, we follow Barth and substitute Wh = Vh in equation (25) to get
Z
Z


∂t η(Vh )dx = −
Θ(Vh+ , Vh− , n) + D(Vh+ , Vh− , n) dσ
K

∂K

where
Θ(V + , V − , n) = {{V }} · Hc (V + , V − , n) − {{θ∗ · n}}

20

is a consistent and conservative numerical flux for the entropy flux θ · n, and
Z


 1 1
1
+
−
+
−
∗
D(V , V , n) =
JV K · Hc (V + , V − , n) − F (V̄ (s)) · n ds
JV K · Hc (V , V , n) − Jθ · nK =
2
2 0
For entropy stability, we have to show that D is positive; a sufficient condition is given by [14, 30]


JV K · Hc (V + , V − , n) − F (V̄ (s)) · n ≥ 0,
∀s ∈ [0, 1]

(27)

which is a generalization to systems of equations, of Osher’s E-flux condition [31] for scalar conservations laws.
Barth [30] has shown that the KMV flux and hence also the KFVS flux, satisfies the E-flux condition; thus the
KFVS-based DG scheme for the Euler equations satisfies the entropy condition for any degree of basis functions.
4.5. Chapman-Enskog distribution function
Consider the Boltzmann equation with the BGK model for the collision term
∂t f + vi ∂i f =

g−f
τR

(28)

If we write the solution as a series in the relaxation time τR as in equation (8), we obtain the following first two
terms of the solution
f0 = g,
f1 = −(∂t g + vi ∂i g)

The distribution function g depends on the macroscopic variables ρ, u, T ; we substitute the inviscid equations
governing these quantities into the expression for f1 to obtain the following Chapman-Enskog velocity distribution function [6]



ρqi ci
2
ρ
1 − β|c|2
(29)
f ce = g 1 − 2 τij ci cj − 2
2p
p
5
where the relaxation time is given by τR = µ/p and ci = vi − ui is the peculiar velocity of the molecules. It can
be checked that U = hψf ce i and Fi + Gi = hvi ψf ce i, i.e., the moments of the Chapman-Enskog distribution
yield the conserved variables and fluxes for the Navier-Stokes equations. Since f ce is a truncated solution in
the Chapman-Enskog expansion, it does not satisfy the Boltzmann equation (28); we can derive an equation for
f ce as ∂t f ce + vi ∂i f ce = Q and it can be shown through explicit but lengthy computations that hψQi = 0. The
Chapman-Enskog distribution is derived for a monatic gas for which Prandtl number is one; however written
in the form of equation (29), the effect of Prandtl number is accounted in the computation of the heat flux q
and we take the above distribution function to be valid for polyatomic gases also [6, 38].
4.6. DG scheme for Navier-Stokes equations
In terms of the entropy variables, the Navier-Stokes equations can be written as
Ã0 ∂t V + Ãi ∂i V − ∂i (Kij ∂j V ) = 0,

Gi = −Kij ∂j V

and the matrices Kij are symmetric positive semi-definite [39, 27]. Taking dot product with entropy variables,
we obtain
∂t η + ∂i θi + ∂i (V · Gi ) = −(∂i V )> Kij ∂j V ≤ 0
(30)
which is essentially the entropy condition. In order to mimic the above property, it is useful to write the DG
scheme for Navier-Stokes equations using entropy variables. For an interior element K, the DG scheme which
is the analoque of the scheme given by equation (14) is
Z
Z
Wh · ∂t U (Vh )dx −
[Fi (Vh ) + Gi (Vh , ∇Vh )] · ∂i Wh dx
K
K
Z
Z
(31)
+
H(Vh+ , ∇Vh+ , Vh− , ∇Vh− , n) · Wh+ dσ −
Hd+ (Vh+ , ∇Wh+ , n) · JVh K dσ = 0
∂K

∂K

21

where n is taken to be the outward normal to the element K. The numerical flux function H which consists of
both inviscid and viscous fluxes is given by kinetic splitting as
H(V + , ∇V + , V − , ∇V − , n) = Hc (V + , V − , n) + Hd (V + , ∇V + , V − , ∇V − , n)
Here Hc is the KFVS numerical flux for Euler equations as given by equation (26) while Hd is the KFVS
numerical flux for the viscous fluxes obtained from Chapman-Enskog distribution. This flux can be written as
Hd (V + , ∇V + , V − , ∇V − , n) = Hd+ (V + , ∇V + , n) + Hd− (V − , ∇V − , n)
Note that equation (31) is a non-symmetric DG scheme for the NS equation at the element level.
4.6.1. Entropy stability
In order to study entropy stability of the scheme given by equation (31) we substitute Wh = Vh to obtain
Z
Z


∂t η(Vh )dx +
Θ(Vh+ , Vh− , n) + D(Vh+ , Vh− , n) dσ
K
Z∂K
Z
 + +

+
−
−
−
−
+
+
Hd (Vh , ∇Vh , n) · Vh + Hd (Vh , ∇Vh , n) · Vh dσ −
Gi (Vh , ∇Vh ) · ∂i V dx = 0
∂K

K

The first two integrals are common with the Euler equations, while the last two terms are the contributions
from NS equations. We observe that
Hd+ (Vh+ , ∇Vh+ , n) · Vh− + Hd− (Vh− , ∇Vh− , n) · Vh+
is a consistent and conservative numerical flux for the viscous entropy flux V · (Gi ni ), while from the NavierStokes equations, we have
Gi (Vh , ∇Vh ) · ∂i Vh = −(∂i Vh )> Kij ∂j Vh ≤ 0
which leads to a cell entropy inequality for the DG scheme which mimics equation (30). The non-symmetric
nature of the scheme (31) is required to obtain the above cell entropy inequality. This type of local or cell
entropy condition does not seem to be valid in the case of the symmetric version of the scheme.
4.6.2. DG scheme for NS equation
We now state the DG scheme for NS equation in a global formulation including penalty terms which is
analogous to the scheme given in equation (17). Let ΓI denote the set of all interior edges/faces of the finite
element grid and let Γ denote the boundary faces. The time continuous DG scheme for Navier-Stokes equations
is given by
Z
Z
Z
Wh · ∂t U (Vh )dx −
[Fi (Vh ) + Gi (Vh , ∇Vh )] · ∂i Wh dx +
H(Vh+ , ∇Vh+ , Vh− , ∇Vh− , n) · JWh K dσ
Ω
Ω
ΓI
Z
Z
(32)
+
+
−
−
+
Hd (Vh , ∇Wh , Vh , ∇Wh , n) · JVh K dσ +
δh (Vh ) · JWh K dσ + NΓ (Vh , Wh ) = 0
ΓI

ΓI

Choosing  = −1 yields the NIPG scheme while  = +1 yields the SIPG scheme. The interior penalty term is
of a similar form as in the scalar problem and given by
δh (Vh ) = Cip

k2 ν
JU (Vh )K
h

22

where we use the coefficient of kinematic viscosity in order to have dimensional consistency. The boundary
conditions are imposed through the boundary terms NΓ in a weak form [22] and given by
Z
NΓ (Vh , Wh ) = H(ΓVh+ , Γ∇Vh+ , ΓVh− , Γ∇Vh− , n) · Wh+ dσ
Γ
Z
Z
+
+
−
−
+
+
+ Hd (ΓVh , Γ∇Wh , ΓVh , Γ∇Wh , n) · [Vh − ΓVh ]dσ + δΓ (Vh ) · Wh+ dσ
Γ

Γ

The operator Γ modifies the boundary trace of the solution to apply the boundary conditions; for example, in
the case of a no-slip boundary

>
ΓVh+ = ΓVh− = V1+ , 0, 0, 0, V5+ ,

Γ∇Vh+ = Γ∇Vh− = ∇Vh+

For an isothermal wall, V5+ is modified by using the specified wall temperature. For an adiabatic wall, ∇Vh+
is modified so that the heat flux is zero while at a farfield boundary with conditions V∞ , we take ΓVh+ = Vh+ ,
ΓVh− = V∞ and Γ∇Vh+ = Γ∇Vh− = ∇Vh+ . The boundary penalty term δΓ is similar to the interior penalty and
is of the form
k2 ν
[U (Vh+ ) − U (ΓVh+ )]
δΓ (Vh ) = Cip
h
4.7. Test case: Order of accuracy
The NS equations do not have 1-D analytical solutions which could be used to compute the convergence
rate of the numerical schemes. Hence we construct exact solutions by the method of manufactured solutions.
We assume the following form for the exact solution
ρ(x) = 1 +

1
cos(2πx),
2

u(x) = 10x2 (1 − x)2 sin(2πx),

T (x) = 1 + 2x2 (1 − x)2

(33)

which solves the following stationary 1-D NS equations with source terms
(ρu)x = f1 ,

(p + ρu2 )x − τx = f2 ,

((ρe + p)u)x − (τ u)x + qx = f3

(34)

together with the boundary conditions u(0) = u(1) = 0 and q(0) = q(1) = 0, which are the adiabatic boundary
conditions corresponding to zero heat flux. The source terms f1 , f2 , f3 are computed by substituting the exact
solutions into the above equations. The viscosity coefficient is µ = 0.01 or µ = 1 which corresponds to convection
dominated and viscous dominated cases, while the ratio of specific heats γ = 7/5 corresponding to polyatomic
gas like air. Computations are performed on a sequence of four uniform grids with N = 40, 80, 160 and 320
elements using P1 , P2 and P3 basis functions. The L2 norm of the error in velocity and temperature and their
derivatives are computed from the exact solution by performing a high order quadrature. The convergence
of the error in velocity and temperature and their derivatives are plotted in figure (8) and (9) for the case of
µ = 0.01 and µ = 1 respectively. The short lines are drawn with slopes indicated beside them. From these
figures we can observe that for the NIPG scheme, the solution and its derivative converge approximately at the
rate of O(hk ). For the SIPG scheme, the solution converges at the rate O(hk+1 ) while the derivative converges
at the rate O(hk ). The SIPG scheme is slightly better in terms of the error in the solution while both schemes
yield nearly similar solution for the derivatives.
4.8. Test case: 1-D shock structure
We compute the shock structure using the Navier-Stokes equations. While it is known that the NS equations
are not valid inside the shock, we use this as a verification test case. The exact solution of the NS equations
for a stationary shock are obtained using the method of [40] which is implemented as a matlab code in [9].
This solution is not truly exact since it uses an ODE solver, and hence we are not able to compute errors in
23

0

−2

10

10

1
1

−4

−6

10

Error in velocity derivative

Error in velocity

10

2

−8

10

NIPG P1
NIPG P2
NIPG P3
SIPG P1
SIPG P2
SIPG P3

−10

10

−12

10

1

3

2
−4

10

NIPG P1
NIPG P2
NIPG P3
SIPG P1
SIPG P2
SIPG P3

−6

10

4
−8

2

10

−2

10

10

3

10
Number of cells

1

2

10

10

3

3

10
Number of cells

(a)

10

(b)
0

−2

10

10

1
−4

−6

10

2

3

−8

10

NIPG P1
NIPG P2
NIPG P3
SIPG P1
SIPG P2
SIPG P3

−10

10

−12

10

Error in velocity derivative

Error in temperature

10

1

10

1

−2

10

−4

10

2

NIPG P1
NIPG P2
NIPG P3
SIPG P1
SIPG P2
SIPG P3

−6

10

4

−8

2

10
Number of cells

3

10

(c)

10

1

10

3

2

10
Number of cells

3

10

(d)

Figure 8: Order of accuracy study for 1-D NS, µ = 0.01: (a) Velocity (b) Velocity derivative (c) Temperature (d) Temperature
derivative

24

0

−2

10

10

1
1

−4

−6

10

Error in velocity derivative

Error in velocity

10

2

−8

10

NIPG P1
NIPG P2
NIPG P3
SIPG P1
SIPG P2
SIPG P3

−10

10

−12

10

1

3

2
−4

10

NIPG P1
NIPG P2
NIPG P3
SIPG P1
SIPG P2
SIPG P3

−6

10

4
−8

2

10

−2

10

10

3

10
Number of cells

1

2

10

10

3

3

10
Number of cells

(a)

10

(b)
0

−2

10

10

Error in temperature derivative

1
−4

Error in temperature

10

2
−6

10

3

−8

10

NIPG P1
NIPG P2
NIPG P3
SIPG P1
SIPG P2
SIPG P3

−10

10

−12

10

1

10

4

1
−2

10

2

−4

10

NIPG P1
NIPG P2
NIPG P3
SIPG P1
SIPG P2
SIPG P3

−6

10

−8

2

10
Number of cells

3

10

(c)

10

1

10

3

2

10
Number of cells

3

10

(d)

Figure 9: Order of accuracy study for 1-D NS, µ = 1.0: (a) Velocity (b) Velocity derivative (c) Temperature (d) Temperature
derivative

25

N=100, P2
1.8

1.7

1.7

1.6

1.6

1.5

1.5
Density

Density

N=100, P1
1.8

1.4
1.3
1.2

1.3
1.2

1.1

1.1
NIPG
SIPG
Exact

1
0.9
-0.14

1.4

-0.135

-0.13

-0.125
x

-0.12

-0.115

NIPG
SIPG
Exact

1
0.9
-0.14

-0.11

-0.135

1.8

1.7

1.7

1.6

1.6

1.5

1.5

1.4
1.3
1.2

-0.12

-0.115

-0.11

1.4
1.3
1.2

1.1

1.1
NIPG
SIPG
Exact

1
0.9
-0.14

-0.125
x

N=200, P2

1.8

Density

Density

N=200, P1

-0.13

-0.135

-0.13

-0.125
x

-0.12

-0.115

NIPG
SIPG
Exact

1
-0.11

0.9
-0.14

-0.135

-0.13

-0.125
x

-0.12

-0.115

-0.11

Figure 10: Density for shock structure problem using entropy variables

the DG solution. The parameters defining the problem are: Mach number ahead of the shock is M1 = 1.5,
γ = 5/3, Pr = 2/3 while the viscosity law is given by µ = µ1 (T /T1 )0.8 where the subscript “1” denotes pre-shock
conditions and µ1 = 0.0005. The time discretization is made with backward Euler implicit scheme and the time
h
step is chosen based on the convective speeds as ∆t = CFL · max(|u|+c)
with CFL = 5. Figure (10) shows the
results for the density using NIPG/SIPG schemes and P1 /P2 basis functions on grids with N = 100, h = 1/400
and N = 200, h = 1/800 elements, while figures (11) and (12) show the corresponding results for shear stress
and heat flux. It is clear that the numerical solutions are able to compute the shock structure quite accurately.
It is more important to resolve the variation of shear stress and heat flux since these exhibit an extremum inside
the shock. In the case of P1 basis functions, the solution derivative is constant inside each cell and so we plot
the constant value of shear stress and heat flux at the center of each cell with a symbol. It is clear from these
plots that both the schemes give accurate solutions and are able to resolve the shock structure well, even for the
shear stress and heat flux. The DG solution compares very well with the analytical solution of the NS equation.
With P2 basis functions, we see that the DG solution is able to capture the variation of shear stress τ and heat
flux q inside each cell very accurately; the accuracy is seen to improve with grid refinement.
5. Summary and conclusions
We have proposed a DG scheme for Navier-Stokes equations which makes use of kinetic flux vector splitting
to construct a numerical flux function for both the convective and diffusive fluxes. These schemes are motivated
by applying them to scalar convection-diffusion problem and studying their accuracy. Energy/entropy stable
26

N=100, P1

N=100, P2

0.005

0.005

0

0

-0.005

-0.005
-0.01
Shear stress

Shear stress

-0.01
-0.015
-0.02
-0.025

-0.015
-0.02
-0.025

-0.03

-0.03

-0.035

-0.035
NIPG
SIPG
Exact

-0.04

NIPG
SIPG
Exact

-0.04

-0.045

-0.045
-0.14

-0.135

-0.13

-0.125
x

-0.12

-0.115

-0.11

-0.14

-0.135

N=200, P1

-0.125
x

-0.12

-0.115

-0.11

N=200, P2

0.005

0.005

0

0

-0.005

-0.005
-0.01
Shear stress

-0.01
Shear stress

-0.13

-0.015
-0.02
-0.025

-0.015
-0.02
-0.025

-0.03

-0.03

-0.035

-0.035
NIPG
SIPG
Exact

-0.04

NIPG
SIPG
Exact

-0.04

-0.045

-0.045
-0.14

-0.135

-0.13

-0.125
x

-0.12

-0.115

-0.11

-0.14

-0.135

-0.13

-0.125
x

-0.12

Figure 11: Shear stress for shock structure problem using entropy variables

27

-0.115

-0.11

N=100, P1

N=100, P2
0.005

0

0

-0.005

-0.005

-0.01

-0.01
Heat flux

Heat flux

0.005

-0.015
-0.02
-0.025

-0.015
-0.02
-0.025

-0.03

-0.03
NIPG
SIPG
Exact

-0.035

NIPG
SIPG
Exact

-0.035

-0.04

-0.04
-0.14

-0.135

-0.13

-0.125
x

-0.12

-0.115

-0.11

-0.14

-0.135

N=200, P1

-0.125
x

-0.12

-0.115

-0.11

N=200, P2

0.005

0.005

0

0

-0.005

-0.005
-0.01
Heat flux

-0.01
Heat flux

-0.13

-0.015
-0.02
-0.025

-0.015
-0.02
-0.025

-0.03

-0.03
NIPG
SIPG
Exact

-0.035

NIPG
SIPG
Exact

-0.035

-0.04

-0.04
-0.14

-0.135

-0.13

-0.125
x

-0.12

-0.115

-0.11

-0.14

-0.135

-0.13

-0.125
x

-0.12

Figure 12: Heat flux for shock structure problem using using entropy variables

28

-0.115

-0.11

schemes are obtained by adding additional stabilization terms which are constructed from the kinetic numerical
flux and interior penalty terms. The resulting schemes which can be either symmetric or non-symmetric in
the discretization of diffusive terms, exhibit good accuracy properties. The symmetric schemes show optimal
convergence rates in numerical tests on the scalar convection-diffusion equation. These good properties are
observed also for Navier-Stokes equations based on entropy variables. Future work with these schemes will be
focused on applying them for two dimensional viscous flow problems to study if the observed convergence rates
carry over to higher dimensions, and trying to couple them with particle methods for rarefied flow situations.
Appendix A. 1-D Navier-Stokes equations
In conservation form the 1-D Navier-Stokes equations can be written as
∂F
∂G
∂U
+
+
=0
∂t
∂x
∂x
where




ρ
U =  ρu  ,
ρe




ρu
F =  p + ρu2  ,
(ρe + p)u




0

−τ
G=
−τ u + q

and the shear stress and heat flux are given by
τ=

4 ∂u
µ ,
3 ∂x

The coefficient of heat conduction is given byκ =
Prandtl number.

q = −κ

µCp
Pr

=

∂T
∂x

µγR
(γ−1) Pr

where R is the gas constant and Pr is the

Appendix B. Polyatomic gases
The Maxwell-Boltzmann distribution function discussed above corresponds to a monatic gas for which the
ratio of specific heats is γ = 53 . Real gases like air which are mostly composed of diatomic molecules N2 and O2
have a value of γ = 1.4 and effectively behave like a diatomic gas. The internal energy of such gases is different
since they have contributions from rotational degrees of freedom. In order to get the correct value of γ and hence
the internal energy, additional internal energy variables can be introduced. In the approach of Deshpande [28],
the collisional invariant corresponding to energy is modified to 21 |v|2 + I while in the BGK schemes [9], it is
2
modified as 12 |v|2 + ξ 2 , with ξ 2 = ξ12 + . . . + ξK
, the value of K depending on γ. A third approach is to use
1
2
δ
|v|
+
I
as
introduced
by
Perthame
[36];
this
has the advantage of preserving the form of the Boltzmann
2
entropy function H = f ln f [36, 30] and will be adopted in the present work. Note that the kinetic split fluxes
resulting from any of the above approaches are identical since there is no splitting with respect to the additional
degrees of freedom. With the above choice of the collisional invariants, the Maxwell-Boltzmann distribution
function is given by [30]
Z
Z
δ
1
ρ
1
( 12 |v−u|2 +I δ )
− 21 |v|2
− RT
g(ρ, u, T ; v, I) =
,
α
=
e
dv
·
e−I dI,
δ= 1
e
d
α(RT )1/(γ−1)
d
+
R
R
γ−1 − 2
The moments of the above distribution function yield the correct expressions for the energy and energy flux for a
polyatomic gas while the mass and momentum equations are unaffected. The entropy variables for a polytropic
gas corresponding to the convex entropy (21) are given by [30]

V = −

s
|u|2
u
1
−
,
, −
γ − 1 2RT RT
RT
29

>

Any constant terms in the entropy variables, e.g., in the definition of the physical entropy s, can be ignored
since they do not change the entropy stability properties of the NS equations or of the DG schemes since only
the jump in the entropy variables or their derivatives appear in the equations.
Appendix C. KFVS fluxes for NS equation in 1-D
The inviscid kinetic split fluxes are given by

F±


ρuA± + ρB ±

(p + ρu2 )A± + ρuB ±
=
(ρe + p)uA± + (ρe + p/2)B ±

where
p
s = u β,

A± =

1
(1 ± erf(s)),
2

1
B± = ± √
exp(−s2 )
2 πβ

The viscous kinetic split fluxes are given by





− τ + 54 βuq βB ±

−τ A± + 45 βqB ±
G± = 

(−uτ + q)A± − 32 τ + 25 βuq B ±

where τ and q are the shear stress and heat flux.
References
References
[1] C. Cercignani, The Boltzmann equation and its applications, Applied mathematical sciences, SpringerVerlag, 1988.
[2] D. I. Pullin, Direct simulation methods for compressible inviscid ideal-gas flow, Journal of Computational
Physics 34 (2) (1980) 231 – 244. doi:10.1016/0021-9991(80)90107-2.
[3] S. M. Deshpande, A second-order accurate kinetic theory based method for inviscid compressible flows,
Tech. Rep. NASA TP-2613, NASA (1986).
[4] J. C. Mandal, S. M. Deshpande, Kinetic flux vector splitting for Euler equations, Computers and Fluids
23 (2) (1994) 447 – 478. doi:DOI:10.1016/0045-7930(94)90050-7.
[5] P. L. Bhatnagar, E. P. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude
processes in charged and neutral one-component systems, Phys. Rev. 94 (1954) 511–525. doi:10.1103/
PhysRev.94.511.
[6] S. Y. Chou, D. Baganoff, Kinetic flux-vector splitting for the Navier-Stokes equations, J. Comput. Phys.
130 (1997) 217–230. doi:10.1006/jcph.1996.5579.
[7] K. H. Prendergast, K. Xu, Numerical hydrodynamics from gas-kinetic theory, J. Comput. Phys. 109 (1993)
53–66. doi:10.1006/jcph.1993.1198.
[8] K. Xu, L. Martinelli, A. Jameson, Gas-kinetic finite volume methods, in: S. M. Deshpande, S. Desai,
R. Narasimha (Eds.), Fourteenth International Conference on Numerical Methods in Fluid Dynamics, Vol.
453 of Lecture Notes in Physics, Springer, Berlin, 1995, pp. 106–111.
30

[9] K. Xu, A gas-kinetic BGK scheme for the Navier-Stokes equations and its connection with artificial dissipation and Godunov method, Journal of Computational Physics 171 (1) (2001) 289 – 335.
doi:DOI:10.1006/jcph.2001.6790.
[10] T. Ohwada, On the construction of kinetic schemes, Journal of Computational Physics 177 (1) (2002) 156
– 175.
[11] M. Torrilhon, K. Xu, Stability and consistency of kinetic upwinding for advection–diffusion equations, IMA
Journal of Numerical Analysis 26 (4) (2006) 686–722.
[12] W. H. Reed, T. R. Hill, Triangular mesh methods for the neutron transport equation, Tech. Rep. LA-UR73-476, Los Alamos Scientific Laboratory (1973).
[13] B. Cockburn, S.-Y. Lin, C.-W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite
element method for conservation laws III: One-dimensional systems, J. Comput. Phys. 84 (1989) 90–113.
[14] B. Cockburn, C.-W. Shu, The Runge-Kutta discontinuous Galerkin method for conservation laws V: Multidimensional systems, J. Comput. Phys. 141 (1998) 199–224.
[15] D. N. Arnold, F. Brezzi, B. Cockburn, L. D. Marini, Unified analysis of discontinuous Galerkin methods
for elliptic problems, SIAM Journal on Numerical Analysis 39 (5) (2002) pp. 1749–1779.
[16] F. Bassi, S. Rebay, A high-order accurate discontinuous finite element method for the numerical solution
of the compressible Navier-Stokes equations, Journal of Computational Physics 131 (2) (1997) 267 – 279.
doi:10.1006/jcph.1996.5572.
[17] B. Cockburn, C.-W. Shu, The local discontinuous Galerkin method for time-dependent convection-diffusion
systems, SIAM J. Numer. Anal. 35 (1998) 2440–2463.
[18] D. N. Arnold, An interior penalty finite element method with discontinuous elements, SIAM Journal on
Numerical Analysis 19 (4) (1982) pp. 742–760.
[19] J. Oden, I. Babuska, C. E. Baumann, A discontinuous hp finite element method for diffusion problems,
Journal of Computational Physics 146 (2) (1998) 491 – 519.
[20] C. E. Baumann, J. T. Oden, A discontinuous hp finite element method for convection-diffusion problems,
Computer Methods in Applied Mechanics and Engineering 175 (3-4) (1999) 311 – 341.
[21] F. Bassi, S. Rebay, Numerical evaluation of two discontinuous Galerkin methods for the compressible
Navier–Stokes equations, International Journal for Numerical Methods in Fluids 40 (1-2) (2002) 197–207.
doi:10.1002/fld.338.
[22] R. Hartmann, P. Houston, An optimal order interior penalty discontinuous Galerkin discretization of the
compressible Navier-Stokes equations, Journal of Computational Physics 227 (22) (2008) 9670 – 9685.
doi:10.1016/j.jcp.2008.07.015.
[23] G. Gassner, F. Lörcher, C. D. Munz, A discontinuous Galerkin scheme based on a space-time expansion
II. Viscous flow equations in multi dimensions, J. Sci. Comput. 34 (3) (2008) 260–286. doi:10.1007/
s10915-007-9169-1.
URL http://dx.doi.org/10.1007/s10915-007-9169-1
[24] K. Xu, Discontinuous Galerkin BGK method for viscous flow equations: One-dimensional systems, SIAM
J. Sci. Comput. 25 (2004) 1941–1963. doi:10.1137/S1064827502416113.

31

[25] H. Liu, K. Xu, A gas-kinetic discontinuous Galerkin method for viscous flow equations, Journal of Mechanical Science and Technology 21 (2007) 1344–1351, 10.1007/BF03177419.
[26] A. Harten, On the symmetric form of systems of conservation laws with entropy, Journal of Computational
Physics 49 (1) (1983) 151 – 164. doi:10.1016/0021-9991(83)90118-3.
[27] F. Shakib, T. J. Hughes, Z. Johan, A new finite element formulation for computational fluid dynamics:
X. The compressible Euler and Navier-Stokes equations, Computer Methods in Applied Mechanics and
Engineering 89 (1-3) (1991) 141 – 219. doi:10.1016/0045-7825(91)90041-4.
[28] S. M. Deshpande, On the Maxwellian distribution, symmetric form and entropy conservation for the Euler
equations, Tech. Rep. TP-2583, NASA Langley, Hampton, VA (1986).
[29] T. Barth, Numerical methods for gasdynamic systems on unstructured grids, in: Kroner, Ohlberger, Rohde
(Eds.), An introduction to recent developments in theory and numerics for conservation laws, Vol. 5 of
Lecture notes in computational science and engineering, Springer-Verlag, 1998, pp. 198–285.
[30] T. Barth, On discontinuous Galerkin approximations of Boltzmann moment systems with Levermore
closure, Computer Methods in Applied Mechanics and Engineering 195 (25-28) (2006) 3311 – 3330.
doi:10.1016/j.cma.2005.06.016.
[31] S. Osher, Convergence of generalized MUSCL schemes, SIAM Journal on Numerical Analysis 22 (5) (1985)
947–961. doi:10.1137/0722057.
[32] G. Hauke, T. J. R. Hughes, A comparative study of different sets of variables for solving compressible
and incompressible flows, Computer Methods in Applied Mechanics and Engineering 153 (1) (1998) 1–44.
doi:doi:10.1016/S0045-7825(97)00043-1.
[33] C.-W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, Journal of Computational Physics 77 (2) (1988) 439 – 471. doi:10.1016/0021-9991(88)90177-5.
[34] B. Cockburn, C.-W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method
for conservation laws II: General framework, Mathematics of Computation 52 (186) (1989) pp. 411–435.
[35] H. W. Liepmann, A. Roshko, Elements of Gasdynamics, Dover, 2002.
[36] B. Perthame, Boltzmann type schemes for gas dynamics and the entropy property, SIAM Journal on
Numerical Analysis 27 (6) (1990) pp. 1405–1421.
[37] T. J. Barth, P. Charrier, Energy stable flux formulas for the discontinuous Galerkin discretization of firstorder nonlinear conservation laws, Tech. Rep. NAS-01-001, NASA (2001).
[38] G. May, B. Srinivasan, A. Jameson, An improved gas-kinetic BGK finite-volume method for threedimensional transonic flow, J. Comput. Phys. 220 (2007) 856–878. doi:10.1016/j.jcp.2006.05.027.
[39] P. Dutt, Stable boundary conditions and difference schemes for Navier-Stokes equations, SIAM Journal on
Numerical Analysis 25 (2) (1988) pp. 245–267.
[40] D. Gilbarg, D. Paolucci, The structure of shock waves in the continuum theory of fluids, Journal of Rational
Mechanics and Analysis 2 (1953) 617.

32

