Vertex-centroid finite volume scheme on tetrahedral grids for
conservation laws
Praveen Chandrashekar, Ashish Garg

arXiv:1112.4238v1 [cs.NA] 19 Dec 2011

TIFR Center for Applicable Mathematics, Bangalore-560065, India

Abstract
Vertex-centroid schemes are cell-centered finite volume schemes for conservation laws
which make use of vertex values to construct high resolution schemes. The vertex values
must be obtained through a consistent averaging (interpolation) procedure. A modified
interpolation scheme is proposed which is better than existing schemes in giving positive
weights in the interpolation formula. A simplified reconstruction scheme is also proposed
which is also more accurate and efficient. For scalar conservation laws, we develop limited
versions of the schemes which are stable in maximum norm by constructing suitable
limiters. The schemes are applied to compressible flows governed by the Euler equations
of inviscid gas dynamics.
Keywords: Finite volume method, unstructured grids, reconstruction, maximum
principle, compressible flows
1. Introduction
Finite volume methods on unstructured grids have become very useful in solving conservation laws, especially those arising in fluid dynamics and in which the computational
domain is of a complex shape. This is mainly due to the maturity of grid generation
tools, especially triangular and tetrahedral grids. Also, the local conservation property
of finite volume methods allows them to automatically compute discontinuous solutions.
Finite volume methods can be classified based on the location where the solution values
are stored [2, 14]. In cell-centered methods, the solution is assumed to be stored at the
centroid of the cells, while in vertex-centered methods, the solution is stored at the vertices
of the grid. In the latter case, a cell has to be constructed around each vertex for which
Email addresses: praveen@math.tifrbng.res.in (Praveen Chandrashekar),
ashish_garg31@yahoo.com (Ashish Garg)

Preprint submitted to Elsevier

March 25, 2018

there are several possible alternatives. The most commonly used method is the dual grid
scheme in which the cell is obtained by joining the cell centroids to the face and edge
centroids. Strictly speaking, the unknowns in a finite volume method are cell averages,
which are not associated with any particular spatial location in the cell. However, the
cell average value gives a second order accurate approximation to the solution at the cell
centroid. For constructing second order accurate schemes, the cell average value can be
taken to be the point value at the cell center, which is the approach taken in the present
work.
The construction of second order accurate schemes leads to several possible alternatives. Since our work uses cell centered schemes, we restrict the discussion to this case.
The finite volume method updates the cell average value of the solution, and the detailed
variation of the solution inside the cell is not known. To construct higher order schemes,
the solution inside each cell must be reconstructed by using the cell average values of the
neigbouring cells. The common approach is to reconstruct the solution variation in each
cell by a linear polynomial by making use of the cell average values in the current cell
and some of its neighbours, which forms the stencil of the reconstruction. This can be
achieved by first approximating the derivatives of the solution in the cell by using either
a Green-Gauss procedure or a least squares procedure [2]. The reconstructed solution
may be limited in order to maintain oscillation-free solutions. While in one dimensional
problems, the TVD criterion [8] can be used, for multi-dimensional problems on unstructured grids, the usual criterion is to ensure that the reconstructed solution lies within the
minimum and maximum of the neigbouring cell average values. The two states obtained
on either side of a cell face by the reconstructed solution are used to compute the flux
across the face using some numerical flux function.
The vertex-centroid scheme [4, 5] is a finite volume scheme for conservation laws on
triangles and tetrahedral grids. The basic unknowns are the cell center values but it also
makes use of vertex values in order to construct the second order scheme. The authors
in [4, 5] obtain a formula for the reconstructed value at the center of a cell face by making
use of the geometrical properties of triangles/tetrahedra, and using the cell center and
vertex values of the solution. There is no need to explicitly compute the derivatives of
the solution in each cell as in the usual cell-centered finite volume schemes. The vertex
values are obtained through an interpolation/averaging procedure that is exact for linear
polynomials. This procedure has been refered to as pseudo-laplacian averaging in the
literature. The interpolation formula for each vertex is a linear combination of the cellcenter values adjacent to that vertex and it is desirable from a stability point of view that
the weights in this formula should be positive. Applications to compressible transonic
flows with weak shocks have been given in [4, 5] and the vertex-centroid scheme was able
to compute stable and accurate solutions even without a limiter, demonstrating good
2

stability properties of the scheme. For stronger shocks, the authors recommend using a
minmod limiter though a stability analysis was not shown.
A vertex-centroid scheme which is similar to the schemes in [4, 5] has been used in [10]
for two dimensional compressible flows, but it uses an area averaging procedure which
is not second order accurate on general meshes. However the authors use very regular
meshes and find that the solutions have extremely low levels of numerical dissipation.
They also comment that the vertex-centroid scheme simplified the inter-mesh transfers
of an unstructured mesh multi-grid scheme, in both the fine-to-coarse and coarse-to-fine
directions.
In this paper, we analyze and improve the vertex-centroid schemes in three respects.
We present a modified averaging procedure which reduces the number of negative weights.
With the new scheme the matrices arising in the determination of the weights have better
behaviour. Secondly, we analyze the accuracy of Frink reconstruction procedure and also
propose a simplied procedure which is theoretically more accurate. Thirdly, we develop
limited versions of the schemes which can be shown to be stable in maximum norm when
applied to scalar conservation laws. Under a CFL condition, the solution at the new time
level in each cell is shown to be bounded between the solution at the neigbouring cells
and the average values at the vertices of the cell. This is an important design criterion for
numerical schemes solving hyperbolic conservation laws since they can develop discontinuous solutions, even starting from smooth initial data. The developed schemes are tested
on several problems governed by the inviscid Euler equations modeling compressible flows
including those with discontinuous solutions. The numerical solutions are compared with
exact solutions or with experimental data to demonstrate the performance of the schemes.
The rest of the paper is organized as follows. The basic formulation of the vertexcentroid scheme is explained. The vertex interpolation schemes are discussed and a new
interpolation scheme is proposed. We then study the error in the reconstruction schemes
of Frink and a simpler reconstruction scheme proposed here which we refer to as upwind
reconstruction. Limited versions of both the reconstruction schemes are described and
their maximum stability is analyzed for a scalar conservation law in two and three dimensions. Finally, the schemes are applied to problems governed by Euler equations of gas
dynamics and their performance in computing discontinuous solutions is demonstrated.
2. Vertex-centroid scheme
Consider a system of conservation laws in d space dimensions which can be written as
d

∂U X ∂Fi
+
=0
∂t
∂xi
i=1
3

(1)

Vji

Uij−
Ui

Uj

Uij+

Vij

Figure 1: Reconstruction in vertex-centroid scheme

Here U is the vector of conserved variables and Fi , i = 1, . . . , d are the Cartesian components of the flux vector and d is the number of spatial dimensions. We solve this equation
numerically on a triangulation made of triangles in 2-D and tetrahedra in 3-D. We assume
that the cells and vertices are numbered in some way; the i’th cell is denoted by Ci . The
face between the i’th cell and the j’th cell is denoted Sij . The set of indices of the cells
sharing a face with Ci is denoted by N (i). Then the semi-discrete finite volume scheme
for cell Ci is
X
dUi
+
H(Uij+ , Uij− , nij ) = 0
(2)
|Ci |
dt
j∈N (i)

where nij ∈ Rd is the normal to Sij and pointing into cell Cj with |nij | = |Sij |, | · | denotes
the measure of the corresponding geometric objects and H is a numerical flux function.
Note that we have used a mid-point quadrature rule to approximate the flux integral
on the cell faces which is sufficient to obtain second order accuracy. The numerical flux
function must satisfy the following conditions:
• Consistency
H(U, U, n) =

d
X

Fi (U )ni

i=1

• Conservation
H(U1 , U2 , n) + H(U2 , U1 , −n) = 0
In the scalar conservation law case, the numerical flux function is a monotone; in this
case, the numerical flux H(a, b, n) is an increasing function of the first argument and a
decreasing function of the second argument.

4

The values Uij+ and Uij− are the two states at the center of face Sij obtained by some
reconstruction process in the cells Ci and Cj respectively, see figure (1). In any triangle/tetrahedron, the line joining any vertex to the center of the opposite face pasess
through the centroid of the cell. The ratio of the distance between any vertex and the cell
centroid, to the distance between the cell centroid and the center of the opposite face is
2:1 for a triangle and 3:1 for a tetrahedron; these properties are used in the reconstruction
schemes. Let Vij be the value at the vertex of cell Ci which is opposite to the face Sij
and Wij denote the arithmetic average of the vertex values on the face Sij . Note that
Wij = Wji since both the values correspond to the same face, but Vij 6= Vji since they
refer to different vertices. Then the reconstructed values on either side of a face Sij is
given by Frink [4] as
Uij− = Uj + αd (Wij − Vji )

Uij+ = Ui + αd (Wij − Vij ),

(3)

where the subscript d refers to the spatial dimension; for 2-D we have α2 = 1/3 while in
3-D, we have α3 = 1/4. We will refer to this as the Frink scheme.
An alternate reconstruction which does not make use of face averages is given by
Uij− = Uj + βd (Uj − Vji )

Uij+ = Ui + βd (Ui − Vij ),

(4)

where β2 = 1/2 and β3 = 1/3. We will refer to this as the upwind scheme. This scheme
is a simple consequence of extrapolating the solution at the cell center and vertex to the
center of the opposite face. Since this scheme uses states on only one side of the face for the
reconstruction, we refer to this as the upwind reconstruction scheme. Computationally,
the scheme given by equation (4) is more simple and efficient compared to (3) since it
does not require the face average values Wij . Later we will see that the upwind scheme
is also advantageous when limiters are used, since it leads to less restrictions on the
reconstruction to be stable in maximum norm.
3. Vertex interpolation scheme
The basic finite volume scheme used here is based on cell center values which are
updated by the finite volume scheme. In this case, we do not know the solution at the
vertices of the grid, which has to be approximated through some consistent interpolation
procedure. To discuss the interpolation scheme, we take a generic vertex and assume
without loss of generality, that it is located at the origin. We then find all the cells
which contain this vertex; let the centroids of these cells have coordinates (xj , yj , zj ),
j = 1, . . . , n and let rj be the Euclidean distance from the j’th centroid to the vertex.

5

In the literature, the following interpolation procedures have been used, an area/volume
averaging procedure and an inverse distance averaging procedure, given by
P 1
P
j r Uj
j |Cj |Uj
V = P
or
V = P j1
j |Cj |
j rj
where the summations are over all the n cells containing the vertex. The inverse distance
averaging can be viewed as a local Shepard interpolation [18]. On general meshes, both
the above procedures are only first order accurate, i.e., they are exact only for constant
functions. On uniform and isotropic grids, they become second order accurate. Note that
for the finite volume scheme to be second order accurate the interpolation scheme must be
atleast second order accurate. In [10], the area averaging procedure has been used together
with high quality meshes leading to very accurate solutions. While these methods may
not be accurate in general, they are very robust since they are convex combinations of
the cell center values and hence the interpolated value is bounded between the minimum
and maximum of the cell center values. This is a useful property while solving hyperbolic
conservation laws which can have discontinuous solutions.
3.1. Laplacian averaging
In order to construct an interpolation procedure which is exact for linear polynomials, Frink [4] extended a method for 2-D Navier-Stokes solutions given in [9] to three
dimensions. In this approach, we first assume an interpolation formula of the type
P
j wj Uj
(5)
V = P
j wj
The necessary conditions for this formula to be exact for linear polynomials are
X
X
X
wj xj = 0,
wj yj = 0,
wj zj = 0
j

j

(6)

j

In order to determine the weights wj , the following minimization problem is solved with
respect to the weights together with the constraints given by equations (6)
min
{wj }

1X
(wj − 1)2
2 j

(7)

This approach tries to keep the weights as close to unity as possible while still satisfying the consistency conditions. In order to solve the constrained minimization problem,
6

we introduce Lagrange multipliers λx , λy , λz ; the constrained minimization problem is
converted into an unconstrained problem in which we have to now minimize the function
X
X
X
1X
(wj − 1)2 + λx
w j xj + λ y
w j y j + λz
wj zj
2 j
j
j
j
with respect to the weights {wj } and the Lagrange multipliers λx , λy , λz as independent
variables. The solution is obtained by setting the derivative of the above function with
respect to each of these variables to zero; the weights are then given by
wj = 1 + λx xj + λy yj + λz zj
where the Lagrange multipliers are obtained by solving the following system of equations
 
P 
 P 2 P
P
λx
P xj
Pxj y2 j P xj zj
P xj




 xj yj
λy = − P y j 
(8)
P
P yj
Pyj z2 j
zj
λz
xj zj
yj zj
zj
In [4] the above set of equations is solved explictly using Cramers rule.
3.2. Consistent Shepard interpolation
We propose a modification of the Shepard interpolation procedure which makes it
exact for linear polynomials. In this case, we begin by assuming an interpolation formula
of the form
P wj
j r Uj
V = P jwj
(9)
j rj

where the weights should now satisfy the consistency conditions
X xj
X yj
X zj
wj = 0,
wj = 0,
wj = 0
r
r
rj
j
j
j
j
j

(10)

In order to determine the weights wj , the minimization probem (7) is solved with respect
to the weights together with the constraints given by equations (10). In order to solve the
constrained minimization, we again introduce Lagrange multipliers λx , λy , λz in terms of
which the solution can be written as
wj = 1 + λx

xj
yj
zj
+ λy + λz
rj
rj
rj

7

where the Lagrange multipliers are obtained by solving the following system of equations


 P  
P  xj 2
P  xj   yj  P  xj   zj 
xj


rj
rj
rj
rj  λ
rj 
  rj  



 2
    x
P x
P
P
P


y
yj
yj
yj
zj  
j
j


λy = − 
(11)
rj 

rj
rj
rj
rj
rj 
P  


P     P    

2
λ
zj
P
z
xj
rj

zj
rj

yj
rj

zj
rj

zj
rj

rj

We solve the above matrix equation explicitly using Cramer’s rule.
Remarks. A first observation in comparing the above two methods is that all the formulae
in the modified Shepard interpolation are independent of the spatial scales. If h is a
measure of the local grid size, then the determinant of the matrix in equation (8) is O(h6 )
and the numerical value of this determinant can become very small. For example, if
h = 10−2 , then the determinant could be O(10−12 ). In the examples, we find that the
determinant can reach the levels of machine precision even for inviscid grids. Due to
this reason, the solution of equation (8) will also suffer from round-off errors. On the
other hand, the determinant of the matrix in equation (11) is independent of h and the
determinant is found to be well behaved, as given in the later examples.
Both of the above methods do not guarantee that the weights wj will be positive. A
necessary condition for equations (6) or (10) to be satisfied with positive weights is that
the cell centers must be distributed on all sides of the vertex. This means that there must
be atleast one cell with xj < 0 and atleast one cell with xj > 0, with similar conditions
for the other two coordinates. Also, these conditions must be satisfied for any orientation
of the coordinate axes. In equations (6) the magnitude of the different terms depends on
the spatial distribution of the cell centers around the vertex while in (10) the magnitudes
depend only on the angular distribution of the cell centers. We conjecture that due to this
reason, the modified Shepard interpolation scheme will be able to preserve the positivity
of the weights better than the Laplacian averaging scheme. We indeed find this to be
the case in our numerical tests given in the later sections. In many cases, the modified
scheme gives all positive weights while the Laplacian scheme yields some negative weights.
In other cases, the number of vertices with negative weights is considerably reduced with
the modified scheme.
Due to the presence of the inverse distance factor in equation (9), the modified Shepard
interpolation scheme gives more weightage to the cells which are closer to the given vertex.
This might be beneficial for anisotropic grids which would arise if we perform anisotropic
grid adaptation, or in the case of boundary layer grids. However in the present work, we
have only considered inviscid problems for which the grids are not highly anistropic.

8

V1
V3

U0

U0

V2

U+

U+

V3

V1
V2

V4

Figure 2: 2-D and 3-D reconstruction in vertex-centroid scheme

4. Accuracy of reconstruction schemes
The derivation of the reconstruction scheme (3) makes use of the geometrical properties
of triangles and tetrahedra, leading to second order accuracy. This scheme makes use of
all the vertex values including the cell center value. A simpler scheme which is also second
order accurate is given by equation (4). This is based on a one dimensional extrapolation
of the vertex and cell center values to the mid-point of the opposite face. In this section
we study the accuracy of these two schemes by assuming that the vertex values are exact.
Assume that U is twice continuously differentiable function of the space coordinates.
We denote by D2 U the symmetric matrix of second derivatives. In any cell C, we have
kD2 U k ≤ K

(12)

for some positive constant K, where the matrix norm k · k is induced by the l2 norm for
vectors [11].
4.1. 2-D reconstruction
Consider a triangle C whose vertex values are V1 , V2 , V3 and whose cell center value is
U0 , see figure (2). Without loss of generality, assume that the cell center is at the origin;
then the position vector of the vertices r1 , r2 , r3 ∈ R2 are such that r1 + r2 + r3 = 0. Let
us look at the reconstructed value on the face formed by vertices 2 and 3. The position
vector of the center of this face is 21 (r2 + r3 ) = − 12 r1 . The exact value at the face center
9

is given by a Taylor series with remainder [16] as
1
1
(13)
Ue = U0 − r1T G0 + r1T H01 r1
2
8
where we denote G0 = ∇U0 and H01 = D2 U (−sr1 /2) for some s ∈ [0, 1]. Using another
Taylor series with remainder term, we can write the vertex values in terms of the cell
center values as
1
Vj = U0 + rjT G0 + rjT Hj rj
(14)
2
where Hj = D2 U (sj rj ) ∈ R2×2 for some sj ∈ [0, 1] and j = 1, 2, 3.
4.1.1. Frink scheme
The reconstructed value on the face formed by vertices 2, 3 is according to equation (3)


1 1
+
U = U0 +
(V2 + V3 ) − V1
(15)
3 2
which upon using equation (14) leads to
1
1
1
1
U + = U0 − r1T G0 + r2T H2 r2 + r3T H3 r3 − r1T H1 r1
2
12
12
6
Using equation (13, the error in the interpolated value can be bounded as
7
1
1
11
Kkr1 k2 + Kkr2 k2 + Kkr3 k3 ≤ Kh2
24
12
12
24
where the diameter of the cell is defined as
|U + − Ue | ≤

h = max krj k
1≤j≤3

4.1.2. Upwind scheme
The reconstructed value is given by equation (4) as
U + = U0 +

1
[U0 − V1 ]
2

(16)

which upon using equation (14) leads to
1
1
U + = U0 − r1T G0 − r1T H1 r1
2
4
The error in this approximation can be bounded as
3
3
|U + − Ue | ≤ Kkr1 k2 ≤ Kh2
8
8
We see that the upwind scheme given by equation (16) has smaller error constant than
the Frink scheme given by equation (15).
10

4.2. 3-D reconstruction
Consider a tetrahedron C whose vertex values are V1 , V2 , V3 , V4 and whose cell center
value is U0 , see figure (2). Without loss of generality, assume that the cell center is
at the origin; then the position vectors of the vertices r1 , r2 , r3 , r4 ∈ R3 are such that
r1 + r2 + r3 + r4 = 0. Let us look at the reconstructed value on the face formed by vertices
2, 3 and 4. The position vector of the center of this face is 31 (r2 + r3 + r4 ) = − 31 r1 . The
exact value at the face center is given by a Taylor expansion with remainder term as
1
1
Ue = U0 − r1T G0 + r1T H01 r1
3
18

(17)

where we denote G0 = ∇U0 and H01 = D2 U (−sr1 /3) for some s ∈ [0, 1]. Using another
Taylor series with remainder term, we can write the vertex values in terms of the cell
center values as
1
Vj = U0 + rjT G0 + rjT Hj rj
(18)
2
where Hj = D2 U (sj rj ) ∈ R3×3 for some sj ∈ [0, 1] and j = 1, . . . , 4.
4.2.1. Frink scheme
The reconstructed value on the face formed by vertices 2, 3 and 4 according to equation (3) is


1 1
+
U = U0 +
(V2 + V3 + V4 ) − V1
(19)
4 3
which upon using equation (18) leads to
1
1
1
1
1
U + = U0 − r1T G0 + r2T H2 r2 + r3T H3 r3 + r4T H4 r4 − r1T H1 r1
3
24
24
24
8
Using equation (17), the error in the interpolated value can be bounded as
|U + − Ue | ≤

1
1
1
11
13
Kkr1 k2 + Kkr2 k2 + Kkr3 k3 + Kkr4 k2 ≤ Kh2
72
24
24
24
36

where the cell diameter is defined as
h = max krj k
1≤j≤4

11

(20)

4.2.2. Upwind scheme
The reconstructed value is given by equation (4) as
U + = U0 +

1
[U0 − V1 ]
3

(21)

which upon using equation (18) leads to
1
1
U + = U0 − r1T G0 − r1T H1 r1
3
6
Using equation (17), the error in this approximation can be bounded as
2
2
|U + − Ue | ≤ Kkr1 k2 ≤ Kh2
9
9
We see that the upwind scheme given by equation (21) has smaller error constant than
the Frink scheme given by equation (19).
Remark. This analysis shows that the simpler upwind scheme is actually slightly more
accurate compared to the Frink scheme. However both methods are second order accurate.
The difference in the error constant does not translate into any appreciable difference in
the numerical solutions obtained by the two methods, atleast for the test cases that we
have studied. The purpose of this analysis was to show that there is no advantage in
using the more elaborate scheme given by equation (3) and the more simpler scheme of
equation (4) is more accurate in theory. In the next section where we construct schemes
stable in the maximum norm, we find that the upwind scheme is more beneficial since it
requires less restrictions to satisfy the maximum stability.
5. Maximum stability for scalar conservation law
Solutions of nonlinear hyperbolic PDEs can develop discontinuous even if the initial
condition is smooth. While first order accurate finite volume methods on unstructured
grids are capable of computing discontinuous solutions in a stable manner [3], the higher
order numerical schemes might produce spurious oscillations near discontinuities. The
classic approach is to construct total variation diminishing schemes which will prevent
the generation of oscillations. The TVD concept however is difficult to extend to multiple
dimensions and unstructured grids. Moreover, there is a negative result by Goodman and
LeVeque [7] that a TVD scheme can be at most first order accurate in more than one
spatial dimension. Due to these reasons, one constructs schemes which are stable in the
maximum norm. This prevents the generation of new extrema and hence oscillations in
12

the solution will be avoided. The first order accurate finite volume schemes combined
with a monotone flux are stable in the maximum norm. The basic idea to achieve this
stability for higher order schemes is to reduce the accuracy of reconstruction to first order
wherever the solution has a discontinuity and might lead to oscillations.
In the case of structured grids, there are well established conditions that can be used
to check the positivity of the schemes. For one dimensional probems, we can write the
second order semi-discrete scheme as
h

dUi
= Ai (Ui+1 − Ui ) + Bi (Ui−1 − Ui )
dt

If all the coefficients A, B are positive, then the scheme will be stable in the maximum
norm under a CFL condition. In particular, if Ui is a local maximum, then it will not
increase in time, and similarly, if Ui is a local minimum, it will not decrease in time. For
unstructured grids, the cells are not arranged along grid lines which makes it difficult to
write the scheme in the above form. However, for the vertex-centroid scheme, this can
be achieved by using the vertex values along with the cell-center values, which is the key
step in the present construction of limited schemes.
We propose to write the semi-discrete vertex-centroid finite volume scheme in the
following form, which not only contains the difference of cell average values but also the
difference between vertex value and cell average value.
|Ci |

X
dUi
=
[Aij (Uj − Ui ) + Bij (Vij − Ui )]
dt

(22)

j∈N (i)

We discretize the above set of ordinary differential equations using Euler time stepping
scheme which leads to the following fully discrete update equation


∆t X
∆t X
∆t X
Uin+1 = 1 −
(Aij + Bij ) Uin +
Aij Uj +
Bij Vij
(23)
|Ci |
|Ci |
|Ci |
j∈N (i)

j∈N (i)

j∈N (i)

If all the coefficients A, B in the above equation are non-negative, then the resulting
scheme satisfies a maximum principle under a CFL condition, i.e.,
min {Uin , Ujn , Vijn } ≤ Uin+1 ≤ max {Uin , Ujn , Vijn }
j∈N (i)

j∈N (i)

(24)

The vertex values are obtained by an interpolation formula of the type given in equation
(5) or (9). If all the weights in the interpolation formula are positive, then the vertex
values are bounded by the cell average values. The higher order accurate scheme must
13

be constructed to ensure the positivity of the coefficients in equation (22) so that spurious oscillations near shocks are prevented. If the second order reconstruction given by
equation (3) or equation (4) is written as
Uij− = Uj + ∆Uji

Uij+ = Ui + ∆Uij ,

(25)

then the corresponding limited reconstruction scheme is defined to be
Uij− = Uj + θji ∆Uji

Uij+ = Ui + θij ∆Uij ,

(26)

where the limiter θ ∈ [0, 1] should be chosen to satisy the positivity conditions. When
θ = 1, the scheme is second order accurate while if θ = 0, it becomes first order accurate.
Consider a numerical flux function which can be written as
H(a, b, n) = H + (a, n) + H − (b, n)

(27)

with H + being a non-decreasing function and H − being a non-increasing function. Then
from the consistency of the numerical flux function,
X
H(Ui , Ui , nij ) = 0
j∈N (i)

we can separate the positive and negative fluxes and write equation (2) as
|Ci |

X
dUi
= −
[H(Uij+ , Uij− , nij ) − H(Ui , Ui , nij )]
dt
j∈N (i)

X H + (Uij+ , nij ) − H + (Ui , nij )
= −
θij ∆Uij
Uij+ − Ui
j∈N (i) |
{z
}
−

X H
j∈N (i) |

−

Pij ≥0

(Uij− , nij )
Uij−

− H − (Ui , nij )
(Uj + θji ∆Uji − Ui )
− Ui
{z
}

Qij ≤0

Due to the properties of the flux function, we know the sign of the terms defined as Pij
and Qij . We next derive conditions on the limiter function θ for the two reconstruction
schemes so that when they are written in the form of equation (22), they have non-negative
coefficients.

14

5.1. Frink scheme
For the reconstruction scheme given by equation (3), i.e., ∆Uij = αd (Wij − Vij ), the
semi-discrete scheme can be written as



X  θij ∆Uij
dUi
θji ∆Uji
|Ci |
=−
Pij
(Vij − Ui ) + Qij 1 +
(Uj − Ui )
(28)
dt
Vij − Ui
Uj − Ui
j∈N (i)

Defining
rij =

∆Uij
αd (Wij − Vij )
= 1
− Ui )
(Uj − Ui )
2

1
(Uj
2

The above scheme can be written as



X 1
dUi
Uj − Ui
1
|Ci |
=−
Pij θij rij
(Vij − Ui ) + Qij 1 − θji rji (Uj − Ui )
dt
2
Vij − Ui
2

(29)

j∈N (i)

The first term contains the reconstruction in cell Ci while the second term is due to the
reconstruction in cell Cj . From the coefficients of the terms containing Pij and Qji , this
scheme will have positive coefficients provided the following two conditions are satisfied
by the limiter function θij ,
θij rij

Uj − Ui
≤ 0,
Vij − Ui

1
1 − θij rij ≥ 0
2

The above two conditions will be satisfied if we choose the limiter as follows
(
0
(Uj − Ui )(Ui − Vij ) ≤ 0
θij =
θ(rij ) otherwise

(30)

where the function θ(r) must satisfy the following conditions:
1. θ(r) = 0 for all r < 0
2. 0 ≤ θ(r) ≤ 1 for all r ∈ [0, ∞)
3. θ(r) ≤ 2r for all r > 0
These conditions are satisfied if we choose
θ(r) = max(0, min(1, 2/r))

(31)

The function θ(r) is sketched in figure (3). The first condition in equation (30) checks
whether Ui lies between the values Vij and Uj ; if it does not lie in this range, the limiter
is set to zero and the order of the reconstruction is reduced to first order.
15

2
2/r
θ(r)

1.8
1.6
1.4
1.2
1
0.8
0.6
0.4
0.2
0
0

0.5

1

1.5

2
r

2.5

3

3.5

4

Figure 3: Limiter function θ(r)

5.2. Upwind scheme
For the reconstruction scheme given by equation (4), i.e., ∆Uij = βd (Ui − Vij ), the
semi-discrete scheme can be written as



X 
dUi
θji ∆Uji
=
Pij θij βd (Vij − Ui ) − Qij 1 +
(Uj − Ui )
(32)
|Ci |
dt
Uj − Ui
j∈N (i)

In this case, the first term inside the curly braces already satisfies the positivity condition.
The coefficient in the second term is positive if we choose θij = θ(rij ), with θ(r) being
given by equation (31) and
rij =

βd (Wij − Vij )
∆Uij
= 1
− Ui )
(Uj − Ui )
2

1
(Uj
2

We note in this case that there are less conditions on the limiter function θij as compared
to the previous scheme since we do not have to enforce the first condition in equation (30).
This is due to the fact that the reconstruction scheme does not make use of face average
values Wij but only the vertex and cell values.
6. Numerical results for Euler equations
We apply the schemes developed here to inviscid compressible flows governed by Euler equations of gas dynamics. These are a system of hyperbolic PDE representing the
16

conservation of mass, momentum and energy. In three space dimensions, the vector of
unknowns U and the Cartesian components of the flux vector are given by
 


ρ
ρui
ρu1 
pδi1 + ρu1 ui 
 


 , Fi = pδi2 + ρu2 ui  , i = 1, . . . , 3
ρu
U =
(33)
2
 


ρu3 
pδi3 + ρu3 ui 
E
(E + p)ui
where ρ is the density, (u1 , u2 , u3 ) are the Cartesian components of the velocity, p is the
pressure and E is the total energy per unit volume given by
1
E = ρε + ρ|u|2
2

(34)

which is the sum of internal and kinetic energy. For an ideal gas, the pressure is related
to the density and internal energy per unit mass ε by
p = (γ − 1)ρε

(35)

where γ is the ratio of specific heats and is a constant for a given gas. In all the numerical
examples, we take γ = 1.4 which corresponds to air under normal conditions.
The Euler equations are solved on tetrahedral grids using a matrix-free LUSGS scheme [17].
In the case of unsteady problems, the 3-stage strong stability preserving Runge-Kutta
scheme is used [19]. As in [4], we reconstruct the primitive variables (ρ, u1 , u2 , u3 , p)
rather than the conserved variables. This is beneficial in maintaining the positivity of
density and pressure. The numerical flux function is either based on the Roe scheme [15]
or the KFVS scheme [13]. The latter scheme is known to be entropy consistent [12] and
hence very robust; we use it for the high Mach number problems which contain strong
shocks. Most of the grids used in this work were generated using the open source tool
GMSH [6].
We also compare the current schemes with the limited schemes of Jameson [10] where
a central limiter is used. The reconstructed values on any face Sij are given by
1
Uij+ = Ui + L(Ui − Vij , Vji − Uj ),
3

1
Uij− = Uj − L(Ui − Vij , Vji − Uj )
3

(36)

where L(a, b) is a limited average which is taken to be
1
L(a, b) = (a + b)[1 − R(a, b)],
2

q

R(a, b) =
17

a−b
,
max(|a| + |b|, h3/2 )

q ∈ [1, 3] (37)

Case
Test 1
Test 2

ρl
1
1

ul
0
-2

pl
1
0.4

ρr
ur
0.125 0
1
2

pr
0.1
0.4

Table 1: Initial conditions for shock tube problem

If q = 1, then the above limiter reduces to the minmod limiter. Note that larger values of
q lead to less limiting. For this scheme, we are not able to prove the maximum stability
property. Note that the term Jameson scheme is used to refer to the reconstruction
scheme given by equations (36) and (37).
All the results are obtained using the modified Shepard interpolation. For some problems, a few of the vertices have negative weights which we do not modify. In each case, we
indicate the number of negative weights using the original interpolation scheme of Frink
and the modified scheme.
6.1. Shock tube problem
The standard shock tube problem is solved on a three dimensional grid. The computational domain is in the form of a channel with a square cross-section and whose axis
is along the x-axis. The side walls of the channel are treated as slip walls. The computational grid contains about 100 points along the axis of the channel. The solution is
taken along the center line of the channel for comparison with the exact solutions of the
Riemann problem. Two standard Riemann problems are considered, the Sod problem
and a problem which involves low densities in the solution. In both cases, the initial
discontinuity is placed at x = 0.5. The initial conditions defining the Riemann problem
are shown in table (1). For Test 1 we use the Roe flux function while for Test 2, we
use the KFVS flux. Limited versions of the reconstruction schemes are used since the
solutions have shock and contact discontinuities. With the original interpolation scheme,
91 vertices have some negative weights with the smallest weight being -0.03, while the
smallest determinant is 1.44 × 10−15 . With the modified interpolation scheme, there is
only one vertex with a negative weight of −1.5 × 10−5 while the smallest determinant is
1.18.
The scheme given by equation (36) did not work for Test 1. Figure (4) shows the
density and velocity obtained for Test 1; all the other schemes are able to give oscillation free solutions. There is no appreciable difference between the Frink and upwind
reconstructions.
Figure (5) shows the density for Test 2 obtained from all the three schemes. Due to the
expansion wave, the density reaches a very low value at the middle of the computational
domain as shown in the zoomed view in figure (5-b). All the schemes are able to preserve
18

1

Exact
Frink
Upwind

Exact
Frink
Upwind

0.6
Velocity

Density

0.8

0.5

0.4

0.2

0

0

0.2

0.4

0.6

0.8

1

0

0

0.2

0.4

0.6

0.8

1

Figure 4: Solution for Test case 1 at time t = 0.2

the positivity of density. The upwind scheme has the least overshoot of density while
Jameson scheme has the highest. The zoomed view in figure (5-c) show the tip of the
expansion wave where the solution has a corner shape. We see that the two new limited
reconstructions proposed here give the best results in terms of preserving the corner shape.
6.2. Transonic flow over Onera M6 wing
This is a very standard test case for compressible flows using finite volume methods
for which experimental results are available. The problem involves transonic flow over
the Onera M6 wing at a free stream Mach number of 0.839 and an angle of attack of
3.06 degrees [1]. The solution consists of shocks on the upper surface of the wing with
a lambda structure on the wing. The tetrahedral grid used in the computations consists
of 341797 cells. With the original vertex averaging procedure, there are 2044 vertices
with negative weights and the minimum determinant is 1.15E-15, while with the modified
averaging procedure, there are only 86 vertices with negative weights and the minimum
determinant is 0.098.
The pressure coefficient at different spanwise stations are plotted and compared with
experimental results in figure (6). The Frink and upwind reconstruction schemes yield
very similar results, while the Jameson scheme is slightly dissipative at the shocks. The
comparison with experimental results is remarkably good for all the schemes considering
that the grid used is not very fine. The pressure variation on the wing surface is shown
in figure (7) for the three schemes with the same contour levels; the lambda shock is well
captured by all of them. However, the Jameson reconstruction can be again seen to give
slightly more diffused shock compared to the other two schemes.

19

0.05

1

Exact
Frink
Jameson
Upwind

1.1

Exact
Frink
Jameson
Upwind

0.04

Exact
Frink
Jameson
Upwind

1.05

0.8

Density

Density

Density

0.03

0.6

1

0.02

0.4

0.95

0.01

0.2

0

0

0.2

0.4

(a)

0.6

0.8

1

0
0.4

0.45

0.5

0.55

0.9
0.8

(b)

0.85

0.9

0.95

(c)

Figure 5: Solution for Test case 2 at time t = 0.15

6.3. Supersonic flow over hemi-sphere
This problem consists of a Mach 3 flow over a hemi-sphere which leads to the formation
of a bow shock in front of the sphere. We use the KFVS scheme for the numerical flux
function, which is known to be entropy consistent and gives physically correct solutions
even for strong shocks. The radius of the hemi-sphere is one and the problem is solved
on two grids: grid G1 has a cell size of 0.05 on the surface of the hemi-sphere with
908330 tetrahedra while grid G2 has a cell size of 0.025 with 1934791 tetrahedra. All the
averaging weights were positive using both methods. In the original procedure of Frink,
the smallest determinants were 2.9e-11 and 3.5e-13. In the case of the modified averaging
procedure, the smallest determinants for the two grids were 1.54 and 1.98 respectively.
The variation of pressure along the stagnation stream-line is shown in figure (8). Along
this line, the shock is normal to the line and we see that there are no oscillations in the
pressure distribution. This is true of the density and velocity also but the plots are not
shown here. The shock width is seen to decrease when we go from the coarser grid G1
to the finer grid G2. Figure (9) shows the pressure contours obtained on G2 using the
three reconstruction methods. It is clear that all of them give similar results without any
oscillations.
6.4. Blast wave
In order to test the limiters under strong shocks, we consider a problem with a large
temperature/pressure discontinuity. The initial condition consists of a small region of
high temperature. The solution develops into an expanding shock wave. In the initial
condition, the density is uniform with a value of 1.228 while the velocity is zero. The
temperature inside a spherical core of radius 5 meters is 8.1 × 107 K while outside the core
20

1

20% span

44% span

1

Exp
Frink
Jameson
Upwind

0.8

Exp
Frink
Jameson
Upwind

1
0.8

0.6
0.6
0.4

0.2

−Cp

−Cp

0.4

0

0.2
0

−0.2

−0.2

−0.4

−0.4

−0.6

−0.6

−0.8
0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

x/c
65% span

0.6

0.8

1

x/c
80% span
1.2
Exp
Frink
Jameson
Upwind

1

0.8

0.6

0.6

0.4

0.4
−Cp

−Cp

0.8

0.2

0.2

0

0

−0.2

−0.2

−0.4

−0.4

−0.6

−0.6

0

0.2

0.4

0.6

0.8

Exp
Frink
Jameson
Upwind

1

1

0

0.2

0.4

x/c

0.6
x/c

Figure 6: Pressure coefficient for Onera M6 wing

21

0.8

1

Frink

Jameson

Upwind

Figure 7: Pressure contours for Onera M6 wing

Frink
Upwind
Jameson

Frink
Upwind
Jameson

10

Pressure

Pressure

10

5

0

5

1

1.2

1.4

(a) Grid G1

0

1

1.2

(b) Grid G2

Figure 8: Pressure variation along stagnation stream-line for flow over hemi-sphere

22

1.4

(a) Frink

(b) Upwind

(c) Jameson

Figure 9: Pressure contour for flow over hemi-sphere using grid G2

it is 298 K. The numerical flux function is based on KFVS scheme. The computational
domain is a cube whose sides are 81 meters and the grid consists of 3182815 tetrahedra.
With the original averaging procedure, there was one vertex with negative weight and the
minimum determinant was 8.7 × 10−4 . With the modified averaging procedure, there was
no negative weight and the minimum determinant was 1.5.
The computations are performed in a time accurate manner using 3-stage RungeKutta scheme and a constant time step of 6 × 10−8 for which the CFL number is less than
one. The Frink and Jameson reconstruction schemes lead to loss of positivity of pressure
in the first iteration itself. Only the upwind scheme is able to give stable solutions for this
problem. In figure (10) we compare the radial pressure variation for the first order scheme
and the second order scheme with upwind reconstruction at two different times; note that
the pressure jump across the shock has a large magnitude in the higher order scheme.
It is seen that the second order scheme is also free of oscillations. Moreover, the second
order scheme is able to resolve more variations in the solution as compared to first order
scheme, especially at the earlier times. This indicates that the limiter is not completely
reducing the scheme to first order. A similarity solution for the radius of the blast wave
is given by Taylor [20] according to which the quantity R5/2 t−1 is a constant, where R
is the radius of the blast wave. From the numerical solution of the pressure variation
obtained with the upwind reconstructed scheme, we determine the blast wave radius and
the results are plotted in figure (11) along with the analytical results. Except for the
23

8

10

8

x 10

2.5

First order
Upwind

9

x 10

First order
Upwind

8

2

6
Pressure

Pressure

7

5
4

1.5

1

3
2

0.5

1
0
0

5

10

15

20
r

25

30

35

40

(a) 25 time steps

0
0

5

10

15

20
r

25

30

35

40

(b) 100 time steps

Figure 10: Pressure variation along a radial line for blast wave problem

early times, the evolution of the radius is well predicted by the current computations as
compared to theoretical results.
7. Conclusions
Vertex-centroid finite volume schemes on tetrahedral grids have a simple structure
even for second order schemes since they make use of vertex values for reconstruction.
Maintaining the positivity of the interpolation weights is important for stability of the
scheme. The new averaging procedure given here is shown to be more successful in maintaining positivity of the interpolation weights and has good properties in terms of having
O(1) determinants. The simplified reconstruction scheme (upwind scheme) proposed here
is shown to be theoretically more accurate as compared to the Frink scheme, and gives
comparably similar results on test cases. By writing the scheme as a difference of vertex
values and cell-center values, limited versions of the scheme are developed and shown to
be stable in maximum norm for scalar conservation laws. Upwind scheme is seen to be
less restrictive in terms of the conditions that the limiter has to satisfy. These schemes
give oscillation free solutions when applied to Euler equations. For the blast wave problem which has a very strong shock, only the upwind limited scheme was stable while the
other schemes failed due to loss of positivity of density/pressure, demonstrating the good
stability properties of the proposed reconstruction scheme.

24

1.6
Upwind
Taylor

1.55

log10(Radius m)

1.5
1.45
1.4
1.35
1.3
1.25

−4

−3.8

−3.6

−3.4
−3.2
log10(Time sec)

Figure 11: Variation of blast radius with time.
R5/2 t−1 =const.

−3

−2.8

−2.6

The straight line corresponds to the relation

References
[1] AGARD. Experimental database for computer program assessment. Technical Report AR-138, AGARD, 1979.
[2] T. J. Barth. Aspects of unstructured grids and finite volume solvers for the Euler
and Navier-Stokes equations. In Lecture Series in Computational Fluid Dynamics,
No. 1994–04. von Karman Institute for Fluid Dynamics, 1994.
[3] T. J. Barth and D. C. Jespersen. The design and analysis of upwind schemes on
unstructured meshes. Number AIAA-89-0366, 1989.
[4] N. T. Frink. Recent progress towards a three-dimensional unstructured Navier-Stokes
flow solver. In 32nd Aerospace Sciences Meeting, number AIAA-94-0061, Reno, NV,
January 1994.
[5] N. T. Frink and S. Z. Pirzadeh. Tetrahedral finite-volume solutions to the NavierStokes equations on complex configurations. Technical Report TM-1998-208961, Langley Research Center, 1998.
[6] C. Geuzaine and J.-F. Remacle. Gmsh: A 3-d finite element mesh generator with
built-in pre- and post-processing facilities. International Journal for Numerical Methods in Engineering, 79(11):1309–1331, 2009.

25

[7] J. B. Goodman and R. J. LeVeque. On the accuracy of stable schemes for 2D scalar
conservation laws. Mathematics of Computation, 45(171):15–21, 1985.
[8] A. Harten. On a class of high resolution total-variation-stable finite-difference
schemes. SIAM Journal on Numerical Analysis, 21(1):1–23, 1984.
[9] D. G. Holmes and S. D. Connell. Solution of the 2D Navier-Stokes equations on unstructured adaptive grids. In AIAA 9th Computational Fluid Dynamics Conference,
June 1989.
[10] A. Jameson and J. C. Vassberg. A vertex-centroid (v-c) scheme for the gas-dynamics
equations. In N. Satofuka, editor, Computational Fluid Dynamics (2000). Springer,
July 2000.
[11] E. Kreyszig. Introductory Functional Analysis with Applications. Wiley, 1989.
[12] S. H. Lui and K. Xu. Entropy analysis of kinetic flux vector splitting schemes for the
compressible Euler equations. Zeitschrift für Angewandte Mathematik und Physik
(ZAMP), 52:62–78, 2001. 10.1007/PL00001540.
[13] J. C. Mandal and S. M. Deshpande. Kinetic flux vector splitting for Euler equations.
Computers and Fluids, 23(2):447 – 478, 1994.
[14] D. J. Mavriplis. Unstructured mesh discretizations and solvers for computational
aerodynamics. In 18th AIAA Computational Fluid Dynamics Conference, 2007.
[15] P. L. Roe. Approximate riemann solvers, parameter vectors, and difference schemes.
Journal of Computational Physics, 43(2):357 – 372, 1981.
[16] W. Rudin. Principles of Mathematical Analysis. McGraw-Hill Book Company, 1976.
[17] D. Sharov and K. Nakahashi. Low speed preconditioning and LU-SGS scheme for
3-D viscous flow computations. AIAA Paper 98-0614, 1998.
[18] D. Shepard. A two-dimensional interpolation function for irregularly-spaced data.
In Proceedings of the 1968 23rd ACM national conference, ACM ’68, pages 517–524,
New York, NY, USA, 1968. ACM.
[19] C.-W. Shu and S. Osher. Efficient implementation of essentially non-oscillatory shockcapturing schemes. Journal of Computational Physics, 77(2):439 – 471, 1988.
[20] G. I. Taylor. The formation of a blast wave by a very intense explosion. I. Theoretical
discussion. Royal Society of London Proceedings Series A, 201:159–174, Mar. 1950.
26

