- Email: [email protected]
Flux difference splittings and limiters for the resolution of contact discontinuities
Flux Difference Splittings and Limiters Resolution of Contact Discontinuities
for the
Stephen F. Davis Department of Mathematics and Statistics NSF Engineering Research Center for Computational Field Simulation Mississippi State University Mississippi State, Mississippi 39762
Transmitted
by Melvin Scott
ABSTRACT Harten,
Lax, and van Leer developed
tions and used it to derive
numerical
In the first part of this paper, compare
various
compare
the resolution
perfectly,
standard subcell
the differences
and show why it is important
1.
of approximate
for systems this theory
formulae.
resolution
method
solulaws.
so that it can be used to
by various
we show how to
methods
can resolve an isolated
for comparison.
Riemann
of conservation
In particular
discontinuities
Since the Roe method
it is a good
We discuss
splitting
of contact
we show that Harten’s iter.
we reformulate
flux difference
of the Roe method.
a theory
methods
In the second
with that
steady
contact
part of the paper,
is the ultimate
compressive
lim-
between
compressive
and noncompressive
limiters
to choose
the limiter
to suit the problem
at hand.
INTRODUCTION
This paper is concerned with the resolution of contact discontinuities in numerical solutions of systems of conservation laws,
In the above, u and f are m-vector functions and x and t are scalars. Such systems are hyperbolic if the Jacobian matrix [af/&] has real eigenvalues and a complete set of linearly independent eigenvectors. The simplest methods for approximating systems of the form (1.1) are APPLIED @
MATHEMATICS
AND
COMPUTATION
65:3-18
(1994)
3
Elsevier Science Inc., 1994
655 Avenue of the Americas,
New York, NY 10010
0096-3003/94/$7.00
4
STEPHEN
the three
point,
explicit
schemes.
U~+‘=u,“-~
In conservation U:,U~+,)
F. DAVIS
form, these look like
-F(U:_r,UF)],
(1.2)
where X = At/Ax and the function F(Ui,U,) is known as the numerical flux. The scheme (1.2) is said to be consistent with the conservation laws (1.1) if the numerical flux satisfies the condition F(u, Additional conditions following, only stable One reason merically numerical
that
u) = f(u).
(1.3)
must be satisfied for the scheme to be stable. schemes are considered.
it is a challenge
to solve hyperbolic
systems
(1.1)
nu-
is that such systems admit discontinuous solutions. Standard methods tend to approximate a discontinuity either by a smooth
profile spread over many grid points or by a less smeared quence of spurious oscillations which pollute the solution. and 80’s considerable resolve discontinuities smearing.
In the
effort was devoted without producing
jump and a seDuring the 70’s
to finding methods spurious oscillations
This effort has been quite successful.
A number
which could or excessive
of methods
de-
veloped then are now routinely applied to very complex practical problems. Despite this success, there is still some debate over the relative merits of various methods. The reason for this debate is that most methods seem to require a trade off between resolution and robustness. Those schemes which require the fewest grid points to resolve a discontinuity tend to produce spurious solutions under extreme circumstances. Those which are more robust tend to produce too much smearing of contact discontinuities and boundary layers. A recent paper by Quirk [l] catalogs some of the more spectacular failures that occur in gasdynamic flow calculations. Van Leer et al. [2] demonstrate some difficulties that occur in viscous flows. This debate is usually carried out by constructing a test case for which some method fails and showing that some other method is not subject to this particular failure. Later somebody shows that the second method fails for a case that the first method can handle and the debate continues. The purpose of the first part of this paper is to change the way that this debate is conducted. Since a large part of the debate concerns the behavior of various approximate Riemann solutions, we propose a rational way to compare approximate Riemann solutions. This permits the focus of the debate to change from a discussion of particular examples to a discussion of the characteristics of various approximations. This information might make it possible to correlate specific properties of the different schemes
Resolution
of Contact
Discontinuities
5
with their behavior and adapt the schemes to the problem at hand. This is a refinement of the Quirk proposal to use different methods in different parts of the domain. As a first step towards carrying out this program, we show that the theory of approximate Riemann solutions developed by Harten, Lax, and van Leer [3] can be used to compare different approximate This is carried out in the next section of this paper.
Riemann solutions. We note that this
method cannot be used to study methods such as Osher’s [4] method or flux vector splitting [5, 61, which are not based on the approximate solution of a Riemann problem. Perhaps an improved theory which includes these methods will be developed in the future. The second part of this paper is a study of limiters.
These
were intro-
duced into numerical methods independently by Boris and Book [7] and van Leer [8] to suppress spurious oscillations in shock calculations. Harten [9] and Roe [lo] showed that limiters can also be used to prevent the spreading of contact discontinuities. Numerical experiments have shown that some limiters (referred to as “overly compressive”) approximate even smooth data by a near discontinuity while others cause actual discontinuities to spread. The behavior of specific limiters is known to practitioners but more needs to be known if this behavior is to be controlled. In this section, we show that Harten’s subcell resolution method is the ultimate compressive fects the performance
limiter. We discuss how varying the parameters afof commonly used limiters and we show that more
work is needed to characterize 2.
THEORY
limiters
OF APPROXIMATE
Many of the more successful
and to characterize
RIEMANN methods
problems.
SOLVERS
for solving
systems
of conserva-
tion laws are descendants of the Godunov [ll] method. The idea behind this method is to approximate the solution by piecewise constants. The discontinuities between the different constant states are taken to be the initial conditions for a local Riemann problem which is solved exactly. The solution to the Riemann problem depends only on the states to the left and right of the discontinuity, UL and un and the ratio x/t. All signal velocities are finite. For time increments small enough so that solutions to neighboring Riemann problems do not interact, a new piecewise constant solution is constructed by averaging the exact solution. It can be shown [3] that this procedure is equivalent to a scheme of the form (1.2) with numerical flux
(2.1)
6
STEPHEN
F. DAVIS
The notation on the right hand side of (2.1) means that the exact defined by (1.1) is evaluated at the point u which is the solution the ray (zr - xj + r)/t = 0 of the Riemann UL = Uj and un = Uj+r. Although
the method
described
problem
in the preceding
with initial paragraph
flux f along
conditions resolves sta-
tionary discontinuities as well as possible, it smears moving discontinuities and is only first order accurate. These shortcomings have led to the development of higher order extensions which will be discussed later in this paper. A second disadvantage of this method is that solutions of the Riemann problem cannot usually be obtained in closed form. Some iterative procedure must be used to compute the solution. Since this solution is required on each cell boundary, the method requires considerably more computation than conventional finite difference methods. Since only the solution at x/t = 0 is required to compute the flux, most of the details of the exact solution are never used. It is obvious that considerable computation could be avoided if the exact solution were replaced by some approximation that was easy to compute. An industry has sprung up to provide approximations and a debate rages over which approximation is “best.” Most approximate Riemann solutions are derived by simplifying the process used to find the exact solution or by simplifying the equations. The discussion of these methods usually concentrate on how much physics is included in a particular approximation. is lost when the solution is averaged. A different
approach
was taken
Unfortunately,
by Harten,
much of the physics
Lax, and van Leer [3]. They
consider what conditions need to be satisfied by an approximate Riemann solution if it is to produce a Godunov-type method in conservation form which is consistent and satisfies an entropy inequality. The approximate Riemann solution depends only on UL, un, and x/t. It turns out that all that is required is that the approximate Riemann solution itself satisfy the conservation laws and the entropy inequality in an averaged sense. They outline how an approximation could be designed which satisfies these minimal requirements and how additional physics could be added to improve the resolution of selected discontinuities. In the following we carry this one step farther and use this theory to compare different approximate Riemann solutions. The approximate derived by assuming
Riemann solution of Harten, Lax, and van Leer is that an initial discontinuity in the solution breaks up
into K waves, each moving at a different velocity aj, j = 1,2, . , k. This situation is illustrated in Figure 1. If we integrate the conservation laws
Resolution
of Contact Discontinuities
I
a2
FIG. 1. Approximate
(1.1) over the rectangle f(UR)_f(UL)
(-h/2,
h/2)
ak
ak-l
Riemann solution model.
x (t, t + At), the result is
=.,(,;-~~)+.,(U;-U;)+...+a,,(u~-u;_l),
where u; =
a,,lAt
1 aj+lAt
- ajAt
s a,At
and u$ = UL and u; = UR. Formula (2.2) defines a flux d#erence can be split in a similar manner. ,,-u,=(u;-u;)+(u;-u;)+-..+(u;-u;_,).
u(x,
(2.2)
t + At) dx
splitting (c.f.
[3]). The solution
(2.3)
If we assume that all wave speeds aj for j = 1, . . . , k and the jumps (~5 + 1 ~5) for interior waves j = 1,. . . , k - 2 are known, we can solve (2.2) and (2.3) for the two remaining jumps.
(u;-UL)
=-,fR-fL-+,-u,)
+(a-oi)(u;-u;_:,]/(w -~1 (2.4)
8
STEPHEN
(UR-u;_1)
=
F. DAVIS
[fR-fL-n,(u,-UL)
-.,)(uf
_g
Any of the following
formulae
F(UL,UR)
= fL +
-U;-#Uk
-al).
can be used to compute
c
nj(u;
(2.5)
a numerical
flux.
-U,-1)
Cl a, < 01 c
= fR -
(Lj(U; - Ll_l)
(2.6)
{jl% >O}
=$L+fR)-&(U;-U;_l). j=l
The last term in the last equation
is usually called the artificial
viscosity.
It determines how fast a particular approximate discontinuity will spread. Some modern methods which are based on wave decompositions do not decompose the solution jump into single valued waves. Examples include various flux vector splitting methods and Osher’s method. More work is needed to analyze these schemes. On the other hand, the method of Roe fits this framework very well and thus serves as a standard for comparison. Briefly, Roe has shown that for certain systems of equations it is possible to construct a mean value matrix things, the condition
f(UR)
The numerical
A(uL,uR) which
- f(UL) = A(UL,UR)(UR
flux for this method
satisfies, among other
- UL).
(2.7)
is defined by taking
where Xj is an eigenvalue and Ri is the corresponding eigenvector of the matrix A(uL,uR).cxj is the projection of (uR - UL) onto the vector Ri. To compare various methods, (fu-fL) and uR-uL in (2.4) and (2.5) are expanded in terms of eigenvalues and eigenvectors of the Roe matrix and substituted into (2.6). After some manipulation, we obtain the following
Resolution
of Contact Discontinuities
result
F(q,,
uR)
= FRo~(UL,UR)
+ 2 j=l
(u&;;) CV!J3k - (u;“i,) QIR1. U+X-
-
In the above a+ = max(O, u) and a- = min(O, u). The first extra term is the difference in numerical flux due to the difference between uj and Xj; the second term is due to the difference between UT -u,” _ r and the corresponding Roe projection. The last two terms are nonzero if the first and last wave speeds have different signs from the corresponding Roe eigenvalues. can occur with the Rusanov and Lax-Friedrichs methods. To illustrate the use of (2.8), we apply it to the Euler equations. is a system
This This
of the form (1.1) with u = [PI P% EIT f(u)
(2.9)
= [P? W2 + P, (E + P)U] T p = (y - 1) [E - 1/2p??].
The exact
solution
to a Riemann
problem
(2.10)
for these equations
would con-
tain three waves but under some circumstances an approximate Riemann solution containing only two waves can give good results. To determine when this approximation is acceptable, we set uz - UT = 0 and u2 = 0. For illustration purposes assume that al, Xi < 0 and X2, us, X3 > 0. Then k = 3 and (2.8)
F(uL,
becomes
UR) = FRO.S(~L> uR) + a1 (a3 (a3 -
X2) al)
+
a3(a1 (a3 -
a2R.2.
-
Xl) al)
WRI
+
m(a3 - X3) (a3 _ ulj ~3R3
(2.11)
In most cases, a; z Xi and u$ z X3 and the second and third terms in (2.11) are small. By contrast, the coefficient of the last term is approximately
10
STEPHEN
F. DAVIS
equal to (u - c)/2, where u is the flow velocity
and c is the speed of sound.
This adds a large
amount of additional artificial viscosity to the contact discontinuity. When u = 0 the Roe scheme resolves the contact perfectly and the effect of the additional dissipation is most severe. This explains the results of van Leer et al. [2] which show that some sort of contact correction is needed to help resolve viscous boundary layers. For inviscid problems with moving contact discontinuities, no first order method gives acceptable results. To illustrate this we compute the solution to a Riemann problem with initial data
Harten [la] attributes this problem to Lax. It is a challenging problem because the solution contains a very large contact discontinuity. Figure 2 shows a solution after 85 time steps, computed using the Godunov method with a Courant number of 0.8. At this time the contact discontinuity has spread to the point where it has reduced the size of the jump across the shock. This is about the best result that can be obtained using a first order method. Figure 3 shows the same result computed using a two wave numerical flux proposed in [13]. As expected, the contact discontinuity is spread more than in Figure 2 but both results are unacceptable. After considering the results shown above, one might ask why any method other than Roe’s is considered at all? One reason is that Roe’s method has difficulty computing slow moving strong shocks. This problem is discussed in [l] and [14] w h ere it is shown that other methods do not share this difficulty. A second reason is that it is very difficult to find a Roe matrix for more general conservation laws. The discovery of the mean value matrix for the Euler equations is a testament to Philip Roe’s ingenuity. Others have extended this work to related problems and the existence of a Roe matrix has been proved for a large class of problems but a general procedure for constructing a Roe matrix is not known. By contrast the procedure of Harten, Lax, and van Leer was developed for general conservation laws. The two wave flux gives good results under many circumstances and simple approximations to the jumps across selected interior waves can be computed using the method of characteristics.
Resolution
of Contact
Discontinuities
11
0.0 0
20
40
60
60
100
X FIG. 2. Solution of Lax problem after 85 steps of the Godunov method with a Courant number of 0.8.
0
20
40
60
00
100
X FIG. 3. Solution of Lax problem after 85 steps of the upwind approximate method with a Courant number of 0.8.
flux
12
3.
STEPHEN
LIMITERS
AND ARTIFICIAL
F. DAVIS
COMPRESSION
It has long been known that it requires too many grid points to compute accurate solutions to practical problems with first order methods. On the other
hand,
when higher
order methods
are used to solve problems
with
discontinuous solutions, they produce spurious oscillations which pollute the entire solution. This problem was solved in the early 70’s when Boris and Book [7] and van Leer [S] independently introduced limiters to suppress oscillations. It was not until much later that Roe [lo] observed that certain limiters also preserve moving contact discontinuities. Since then, we have looked for some simple condition that characterizes nonspreading limiters, so far without success. Numerical experiments have identified particular limiters as spreading or nonspreading but analytic conditions which distinguish between them do not yet exist. It is important to distinguish between limiters because numerical experiments have shown that the nonspreading limiters tend to turn a smooth profile into a discontinuous one. Zalesak [15] has created the problem of a moving ellipse which demonstrates how poorly these limiters can perform. We need a better understanding of what limiters do so that we can control their behavior. In the remainder of this paper we summarize what is currently known about nonspreading limiters and we carry out some numerical experiments to determine how the performance of a basic class of limiters changes with changes in parameter values. The ultimate nonspreading limiter is the Harten method [16]. To show that this is a limiter method, linear
advection
eouation
Subcell consider
g+g=o.
Resolution the simple
(3.1)
For the numerical methods considered in this paper, numerical solution values represent approximate cell averages of the actual solution. Thus a numerical solution given by two constant regions separated by a single intermediate value as shown in Figure 4(a) can be reconstructed as a single discontinuity as shown in Figure 4(b). From the definition of cell averages the location of the discontinuity is 5 = ax(uA4
-
uR)/(uL
-
where the notation comes from Figure 4. The discontinuity moves to the right with unit speed. cell M provided
(3.2)
UR),
It will remain
in
that
at <
AX
-
X
<
aX(uL
-
u,f,f)/(uL - UR).
(3.3)
Resolution
of
Contact
Discontinuities
13
(a) 1.)
3
OUPA
0.5 -
UR 0 0.0 L 0.0
0
0
0
0.5
c: 1 .o
X
(b)
“L
UR 0
FIG. 4. Arrangement subcell resolution method.
!Ax
X
of cell averages
discontinuity
(b) for
If, under these circumstances, the discontinuous profile is advected tance At to the right, the new cell average values are
a dis-
n+l_
UL
uz
n+l_
-u~-~(u;4-u~).
n+l
=
UM UR If condition
-
(a) and reconstructed
(3.3) is violated
(3.4)
Un,
R
and the discontinuity
moves into the next cell
14
STEPHEN
to the right, the new cell averages n+l
=Un
n+l
=un
UL UM
n+i
UR
F. DAVIS
are
L (3.5)
L
=un .-g-(u;-u;l)
-
(UZ-ug.
In each of these cases the pattern of two constant regions separated by a single intermediate point is preserved. This is usually called a perfect representation of the discontinuity. Conditions (3.3), (3.4), and (3.5) equations. These equations are n+l
UM
=Un
can be combined
into a single set of
M+min[v(u”,-u~)(uT;-U~)]
~~+l=u~+max[O,V(u~-7&)-(7$-U&)], where u = &j&r. After additional
manipulation,
(3.6)
these equations
can be combined
into the form studied by Sweby [17]. F or advection
and put
to the right, this form is (3.7a)
where (3.7b) and Au,-1
(3.7c)
FT.
Space (3.6)
limitations
prohibit
us from carrying
out the manipulations
but
can be put into the form (3.7) with 4 defined by $(T) = min
[&?I.
(3.3)
This result establishes a relationship between subcell resolution and limiter methods. An interesting aspect of (3.8) is that it represents the outer boundary of the region where nonoscillatory solutions are possible. This condition was first discovered by van Leer [18] and has been rediscovered since by many people in many different ways.
Resolution of Contact Discontinuities
15
2.5
20
0
40
60
80
100
X FIG. 5. Computed solution for a periodic contact discontinuity after every 1200 steps until step 12000 using the Godunov-MUSCL method, a Courant number of 0.25, and the Fromm limiter.
Consider
limiters of the form 4(r)
= Kmin[l,r].
(3.9)
Methods based on these limiters are only first order accurate but these limiters contain only one parameter. It is well known that K = 1 causes square wave data to spread and that K = 2 is a nonspreading limiter. Our numerical experiments indicate that K = 1.5 does not spread. These tests were conducted using square wave initial data and periodic boundary conditions. The solution was sampled each time that the wave returned to its initial position. Figures 5 and 6 show typical spreading and nonspreading results. We could not find an exact spreading boundary and we could not detect any significant dependence on the Courant Two commonly used limiters are $(7-) = min[2,2r,
$l
number u.
+ 7-)]
(3.10)
max(1, r)]
(3.11)
and 4(r)
= min[2,2r,
when T > 0 and $ = 0 otherwise.
The limiter (3.10) reduces to the Fromm
16
STEPHEN
0.01 0
20
40
60
80
F. DAVIS
100
X
FIG. 6. Computed solution for a periodic contact discontinuity after every 1200 steps until step 12000 using the Godunov-MUSCL method, a Cow-ant number of 0.25, and the superbee limiter.
method when the solution is smooth and the limiter (3.11) is Roe’s superbee limiter. Both reduce to (3.9) with K = 2 near the corners of an approximate discontinuity but (3.10) spreads and (3.11) does not. Numerical experiments indicate that any,limiter that reduces to a strictly convex combination of 1 and T near T = 1 will spread. Such convex combinations represent linearly stable second order methods. The max(r, 1) term in superbee is less dissipative than any such method. Thus our numerical experiments indicate that this reduced dissipation near r = 1 is necessary to prevent spreading but we do not know how to prove this. We conclude this section by showing some results for the Zalesak [15] ellipse problem. This problem makes all limiter methods look bad. In general, spreading limiters give better results for smooth solutions and nonspreading limiters give better results when discontinuities are present. The Zalesak problem falls somewhere in between. The discontinuous slope and steep profile activate the limiters. Nonspreading limiters try to turn the original profile into a square wave and spreading limiters eventually smear it out beyond recognition. Figures 7 and 8 shown that the Fromm limiter gives more attractive results than superbee but neither solution is acceptable. Zalasak suggests that this problem be solved with very high order methods. This may or may not be necessary but additional work is
Resolution
of Contact Discontinuities
0.0
0.0
20.0
17
40.0
60.0
60.0
100.0
X FIG. 7. Computed solution for a periodic Fromm limiter and a Courant number of 0.5.
half ellipse after 2000 steps using the
2.0
1.5
3
1.0
0.5
0.0
’
0.0
20.0
40.0
60.0
80.0
1
X FIG. 8. Computed solution for a periodic Superbee limiter and a Courant number of 0.5.
half ellipse afer 2000 steps using the
REFERENCES 1
J. J. Quirk,
A Contribution
to the Great
Riemann
Solver
Debate,
ICASE
18
2
3
4 5
STEPHEN F. DAVIS
Report No. 92-64 (1992). B. van Leer, J. L. Thomas, P. L. Roe, and R. W. Newsome, A Comparison of Numerical Flux Formulas for the Euler and Navier-Stokes Equations, 8th AIAA CFD Conference, Honolulu, Hawaii, 1987. A. Harten, P. D. Lax, and B. van Leer, On upstream differencing and Godunov-type schemes for hyperbolic conservation laws, SIAM Rev. 25:3561 (1983). S. Osher and F. Solomon, Upwind difference schemes for hyperbolic systems of conservation laws, Math Comput. 38:339-374 (1982). J. Steger and R. F. Warming, Flux vector splitting of the inviscid gasdynamic equations with applications to finite difference methods, J. Comput. Phys.
6 7 8
9 10
11
12
13
40:263-293
Statist.
14 15
16 17
Comput.
9:445-473
(1988).
T. W. Roberts, The behavior if flux difference splitting schemes near slowly moving shock waves, J. Comput. Phys. 90:141-160 (1990). S. T. Zalesak, A preliminary comparison of modern shock-capturing schemes: Linear advection, in Advances in Computer Methods for Partial Differential Equations, VI (R. Vichnevetsky, R. S. Stepleman, Eds.) IMACS, New Brunswick, N.J. 1987. A. Harten, EN0 schemes with subcell resolution, J. Comput. Phys. 83:148184 (1989). P. K. Sweby, High resolution servation
18
(1981).
B. van Leer, Flux-vector splitting for the Euler equations, in Lecture Notes in Physics, Springer, 1982, pp. 507-512. J. P. Boris and D. L. Book, Flux-corrected transport. I. SHASTA, A fluid transport algorithm that works, J. Comput. Phys. 11:38-69 (1973). B. van Leer, Towards the ultimate conservative difference scheme. I. The quest for monotonicity, in Lecture Notes in Physics, Springer, Berlin, 1973, pp. 163-168. A. Harten, The artificial compression method for computation of shocks and contact discontinuities., Comm. Pure Appl. Math. 30:611-638 (1977). P. L. Roe and M. J. Baines, Asymptotic behavior of some non-linear schemes for linear advection, in Proceedings 5th GAMM Conference on Numerical Methods in Fluid Mechanics (M. Pandolfi, R. Piva, Eds.) 1983. S. K. Godunov, A finite difference method for the numerical computation of discontinuous solutions of the equations of fluid dynamics, Mat. Zb. 47:271290 (1959). A. Harten, High resolution schemes for hyperbolic conservation laws, J. Comput. Phys. 49:357-393 (1983). S. F. Davis, Simplified second-order Godunov-type methods, SIAM J. Sci.
laws,
B. van Leer, tonicity Phys.
SIAM
Towards
the ultimate
and conservation 14(4):361-370
schemes
using flux limiters
J. Numer. Anal.
(1974).
combined
21:995-1011
conservative in a second
for hyperbolic
con-
(1984).
difference order
scheme.
scheme,
II. MonoJ.
Comput.
for the
Stephen F. Davis Department of Mathematics and Statistics NSF Engineering Research Center for Computational Field Simulation Mississippi State University Mississippi State, Mississippi 39762
Transmitted
by Melvin Scott
ABSTRACT Harten,
Lax, and van Leer developed
tions and used it to derive
numerical
In the first part of this paper, compare
various
compare
the resolution
perfectly,
standard subcell
the differences
and show why it is important
1.
of approximate
for systems this theory
formulae.
resolution
method
solulaws.
so that it can be used to
by various
we show how to
methods
can resolve an isolated
for comparison.
Riemann
of conservation
In particular
discontinuities
Since the Roe method
it is a good
We discuss
splitting
of contact
we show that Harten’s iter.
we reformulate
flux difference
of the Roe method.
a theory
methods
In the second
with that
steady
contact
part of the paper,
is the ultimate
compressive
lim-
between
compressive
and noncompressive
limiters
to choose
the limiter
to suit the problem
at hand.
INTRODUCTION
This paper is concerned with the resolution of contact discontinuities in numerical solutions of systems of conservation laws,
In the above, u and f are m-vector functions and x and t are scalars. Such systems are hyperbolic if the Jacobian matrix [af/&] has real eigenvalues and a complete set of linearly independent eigenvectors. The simplest methods for approximating systems of the form (1.1) are APPLIED @
MATHEMATICS
AND
COMPUTATION
65:3-18
(1994)
3
Elsevier Science Inc., 1994
655 Avenue of the Americas,
New York, NY 10010
0096-3003/94/$7.00
4
STEPHEN
the three
point,
explicit
schemes.
U~+‘=u,“-~
In conservation U:,U~+,)
F. DAVIS
form, these look like
-F(U:_r,UF)],
(1.2)
where X = At/Ax and the function F(Ui,U,) is known as the numerical flux. The scheme (1.2) is said to be consistent with the conservation laws (1.1) if the numerical flux satisfies the condition F(u, Additional conditions following, only stable One reason merically numerical
that
u) = f(u).
(1.3)
must be satisfied for the scheme to be stable. schemes are considered.
it is a challenge
to solve hyperbolic
systems
(1.1)
nu-
is that such systems admit discontinuous solutions. Standard methods tend to approximate a discontinuity either by a smooth
profile spread over many grid points or by a less smeared quence of spurious oscillations which pollute the solution. and 80’s considerable resolve discontinuities smearing.
In the
effort was devoted without producing
jump and a seDuring the 70’s
to finding methods spurious oscillations
This effort has been quite successful.
A number
which could or excessive
of methods
de-
veloped then are now routinely applied to very complex practical problems. Despite this success, there is still some debate over the relative merits of various methods. The reason for this debate is that most methods seem to require a trade off between resolution and robustness. Those schemes which require the fewest grid points to resolve a discontinuity tend to produce spurious solutions under extreme circumstances. Those which are more robust tend to produce too much smearing of contact discontinuities and boundary layers. A recent paper by Quirk [l] catalogs some of the more spectacular failures that occur in gasdynamic flow calculations. Van Leer et al. [2] demonstrate some difficulties that occur in viscous flows. This debate is usually carried out by constructing a test case for which some method fails and showing that some other method is not subject to this particular failure. Later somebody shows that the second method fails for a case that the first method can handle and the debate continues. The purpose of the first part of this paper is to change the way that this debate is conducted. Since a large part of the debate concerns the behavior of various approximate Riemann solutions, we propose a rational way to compare approximate Riemann solutions. This permits the focus of the debate to change from a discussion of particular examples to a discussion of the characteristics of various approximations. This information might make it possible to correlate specific properties of the different schemes
Resolution
of Contact
Discontinuities
5
with their behavior and adapt the schemes to the problem at hand. This is a refinement of the Quirk proposal to use different methods in different parts of the domain. As a first step towards carrying out this program, we show that the theory of approximate Riemann solutions developed by Harten, Lax, and van Leer [3] can be used to compare different approximate This is carried out in the next section of this paper.
Riemann solutions. We note that this
method cannot be used to study methods such as Osher’s [4] method or flux vector splitting [5, 61, which are not based on the approximate solution of a Riemann problem. Perhaps an improved theory which includes these methods will be developed in the future. The second part of this paper is a study of limiters.
These
were intro-
duced into numerical methods independently by Boris and Book [7] and van Leer [8] to suppress spurious oscillations in shock calculations. Harten [9] and Roe [lo] showed that limiters can also be used to prevent the spreading of contact discontinuities. Numerical experiments have shown that some limiters (referred to as “overly compressive”) approximate even smooth data by a near discontinuity while others cause actual discontinuities to spread. The behavior of specific limiters is known to practitioners but more needs to be known if this behavior is to be controlled. In this section, we show that Harten’s subcell resolution method is the ultimate compressive fects the performance
limiter. We discuss how varying the parameters afof commonly used limiters and we show that more
work is needed to characterize 2.
THEORY
limiters
OF APPROXIMATE
Many of the more successful
and to characterize
RIEMANN methods
problems.
SOLVERS
for solving
systems
of conserva-
tion laws are descendants of the Godunov [ll] method. The idea behind this method is to approximate the solution by piecewise constants. The discontinuities between the different constant states are taken to be the initial conditions for a local Riemann problem which is solved exactly. The solution to the Riemann problem depends only on the states to the left and right of the discontinuity, UL and un and the ratio x/t. All signal velocities are finite. For time increments small enough so that solutions to neighboring Riemann problems do not interact, a new piecewise constant solution is constructed by averaging the exact solution. It can be shown [3] that this procedure is equivalent to a scheme of the form (1.2) with numerical flux
(2.1)
6
STEPHEN
F. DAVIS
The notation on the right hand side of (2.1) means that the exact defined by (1.1) is evaluated at the point u which is the solution the ray (zr - xj + r)/t = 0 of the Riemann UL = Uj and un = Uj+r. Although
the method
described
problem
in the preceding
with initial paragraph
flux f along
conditions resolves sta-
tionary discontinuities as well as possible, it smears moving discontinuities and is only first order accurate. These shortcomings have led to the development of higher order extensions which will be discussed later in this paper. A second disadvantage of this method is that solutions of the Riemann problem cannot usually be obtained in closed form. Some iterative procedure must be used to compute the solution. Since this solution is required on each cell boundary, the method requires considerably more computation than conventional finite difference methods. Since only the solution at x/t = 0 is required to compute the flux, most of the details of the exact solution are never used. It is obvious that considerable computation could be avoided if the exact solution were replaced by some approximation that was easy to compute. An industry has sprung up to provide approximations and a debate rages over which approximation is “best.” Most approximate Riemann solutions are derived by simplifying the process used to find the exact solution or by simplifying the equations. The discussion of these methods usually concentrate on how much physics is included in a particular approximation. is lost when the solution is averaged. A different
approach
was taken
Unfortunately,
by Harten,
much of the physics
Lax, and van Leer [3]. They
consider what conditions need to be satisfied by an approximate Riemann solution if it is to produce a Godunov-type method in conservation form which is consistent and satisfies an entropy inequality. The approximate Riemann solution depends only on UL, un, and x/t. It turns out that all that is required is that the approximate Riemann solution itself satisfy the conservation laws and the entropy inequality in an averaged sense. They outline how an approximation could be designed which satisfies these minimal requirements and how additional physics could be added to improve the resolution of selected discontinuities. In the following we carry this one step farther and use this theory to compare different approximate Riemann solutions. The approximate derived by assuming
Riemann solution of Harten, Lax, and van Leer is that an initial discontinuity in the solution breaks up
into K waves, each moving at a different velocity aj, j = 1,2, . , k. This situation is illustrated in Figure 1. If we integrate the conservation laws
Resolution
of Contact Discontinuities
I
a2
FIG. 1. Approximate
(1.1) over the rectangle f(UR)_f(UL)
(-h/2,
h/2)
ak
ak-l
Riemann solution model.
x (t, t + At), the result is
=.,(,;-~~)+.,(U;-U;)+...+a,,(u~-u;_l),
where u; =
a,,lAt
1 aj+lAt
- ajAt
s a,At
and u$ = UL and u; = UR. Formula (2.2) defines a flux d#erence can be split in a similar manner. ,,-u,=(u;-u;)+(u;-u;)+-..+(u;-u;_,).
u(x,
(2.2)
t + At) dx
splitting (c.f.
[3]). The solution
(2.3)
If we assume that all wave speeds aj for j = 1, . . . , k and the jumps (~5 + 1 ~5) for interior waves j = 1,. . . , k - 2 are known, we can solve (2.2) and (2.3) for the two remaining jumps.
(u;-UL)
=-,fR-fL-+,-u,)
+(a-oi)(u;-u;_:,]/(w -~1 (2.4)
8
STEPHEN
(UR-u;_1)
=
F. DAVIS
[fR-fL-n,(u,-UL)
-.,)(uf
_g
Any of the following
formulae
F(UL,UR)
= fL +
-U;-#Uk
-al).
can be used to compute
c
nj(u;
(2.5)
a numerical
flux.
-U,-1)
Cl a, < 01 c
= fR -
(Lj(U; - Ll_l)
(2.6)
{jl% >O}
=$L+fR)-&(U;-U;_l). j=l
The last term in the last equation
is usually called the artificial
viscosity.
It determines how fast a particular approximate discontinuity will spread. Some modern methods which are based on wave decompositions do not decompose the solution jump into single valued waves. Examples include various flux vector splitting methods and Osher’s method. More work is needed to analyze these schemes. On the other hand, the method of Roe fits this framework very well and thus serves as a standard for comparison. Briefly, Roe has shown that for certain systems of equations it is possible to construct a mean value matrix things, the condition
f(UR)
The numerical
A(uL,uR) which
- f(UL) = A(UL,UR)(UR
flux for this method
satisfies, among other
- UL).
(2.7)
is defined by taking
where Xj is an eigenvalue and Ri is the corresponding eigenvector of the matrix A(uL,uR).cxj is the projection of (uR - UL) onto the vector Ri. To compare various methods, (fu-fL) and uR-uL in (2.4) and (2.5) are expanded in terms of eigenvalues and eigenvectors of the Roe matrix and substituted into (2.6). After some manipulation, we obtain the following
Resolution
of Contact Discontinuities
result
F(q,,
uR)
= FRo~(UL,UR)
+ 2 j=l
(u&;;) CV!J3k - (u;“i,) QIR1. U+X-
-
In the above a+ = max(O, u) and a- = min(O, u). The first extra term is the difference in numerical flux due to the difference between uj and Xj; the second term is due to the difference between UT -u,” _ r and the corresponding Roe projection. The last two terms are nonzero if the first and last wave speeds have different signs from the corresponding Roe eigenvalues. can occur with the Rusanov and Lax-Friedrichs methods. To illustrate the use of (2.8), we apply it to the Euler equations. is a system
This This
of the form (1.1) with u = [PI P% EIT f(u)
(2.9)
= [P? W2 + P, (E + P)U] T p = (y - 1) [E - 1/2p??].
The exact
solution
to a Riemann
problem
(2.10)
for these equations
would con-
tain three waves but under some circumstances an approximate Riemann solution containing only two waves can give good results. To determine when this approximation is acceptable, we set uz - UT = 0 and u2 = 0. For illustration purposes assume that al, Xi < 0 and X2, us, X3 > 0. Then k = 3 and (2.8)
F(uL,
becomes
UR) = FRO.S(~L> uR) + a1 (a3 (a3 -
X2) al)
+
a3(a1 (a3 -
a2R.2.
-
Xl) al)
WRI
+
m(a3 - X3) (a3 _ ulj ~3R3
(2.11)
In most cases, a; z Xi and u$ z X3 and the second and third terms in (2.11) are small. By contrast, the coefficient of the last term is approximately
10
STEPHEN
F. DAVIS
equal to (u - c)/2, where u is the flow velocity
and c is the speed of sound.
This adds a large
amount of additional artificial viscosity to the contact discontinuity. When u = 0 the Roe scheme resolves the contact perfectly and the effect of the additional dissipation is most severe. This explains the results of van Leer et al. [2] which show that some sort of contact correction is needed to help resolve viscous boundary layers. For inviscid problems with moving contact discontinuities, no first order method gives acceptable results. To illustrate this we compute the solution to a Riemann problem with initial data
Harten [la] attributes this problem to Lax. It is a challenging problem because the solution contains a very large contact discontinuity. Figure 2 shows a solution after 85 time steps, computed using the Godunov method with a Courant number of 0.8. At this time the contact discontinuity has spread to the point where it has reduced the size of the jump across the shock. This is about the best result that can be obtained using a first order method. Figure 3 shows the same result computed using a two wave numerical flux proposed in [13]. As expected, the contact discontinuity is spread more than in Figure 2 but both results are unacceptable. After considering the results shown above, one might ask why any method other than Roe’s is considered at all? One reason is that Roe’s method has difficulty computing slow moving strong shocks. This problem is discussed in [l] and [14] w h ere it is shown that other methods do not share this difficulty. A second reason is that it is very difficult to find a Roe matrix for more general conservation laws. The discovery of the mean value matrix for the Euler equations is a testament to Philip Roe’s ingenuity. Others have extended this work to related problems and the existence of a Roe matrix has been proved for a large class of problems but a general procedure for constructing a Roe matrix is not known. By contrast the procedure of Harten, Lax, and van Leer was developed for general conservation laws. The two wave flux gives good results under many circumstances and simple approximations to the jumps across selected interior waves can be computed using the method of characteristics.
Resolution
of Contact
Discontinuities
11
0.0 0
20
40
60
60
100
X FIG. 2. Solution of Lax problem after 85 steps of the Godunov method with a Courant number of 0.8.
0
20
40
60
00
100
X FIG. 3. Solution of Lax problem after 85 steps of the upwind approximate method with a Courant number of 0.8.
flux
12
3.
STEPHEN
LIMITERS
AND ARTIFICIAL
F. DAVIS
COMPRESSION
It has long been known that it requires too many grid points to compute accurate solutions to practical problems with first order methods. On the other
hand,
when higher
order methods
are used to solve problems
with
discontinuous solutions, they produce spurious oscillations which pollute the entire solution. This problem was solved in the early 70’s when Boris and Book [7] and van Leer [S] independently introduced limiters to suppress oscillations. It was not until much later that Roe [lo] observed that certain limiters also preserve moving contact discontinuities. Since then, we have looked for some simple condition that characterizes nonspreading limiters, so far without success. Numerical experiments have identified particular limiters as spreading or nonspreading but analytic conditions which distinguish between them do not yet exist. It is important to distinguish between limiters because numerical experiments have shown that the nonspreading limiters tend to turn a smooth profile into a discontinuous one. Zalesak [15] has created the problem of a moving ellipse which demonstrates how poorly these limiters can perform. We need a better understanding of what limiters do so that we can control their behavior. In the remainder of this paper we summarize what is currently known about nonspreading limiters and we carry out some numerical experiments to determine how the performance of a basic class of limiters changes with changes in parameter values. The ultimate nonspreading limiter is the Harten method [16]. To show that this is a limiter method, linear
advection
eouation
Subcell consider
g+g=o.
Resolution the simple
(3.1)
For the numerical methods considered in this paper, numerical solution values represent approximate cell averages of the actual solution. Thus a numerical solution given by two constant regions separated by a single intermediate value as shown in Figure 4(a) can be reconstructed as a single discontinuity as shown in Figure 4(b). From the definition of cell averages the location of the discontinuity is 5 = ax(uA4
-
uR)/(uL
-
where the notation comes from Figure 4. The discontinuity moves to the right with unit speed. cell M provided
(3.2)
UR),
It will remain
in
that
at <
AX
-
X
<
aX(uL
-
u,f,f)/(uL - UR).
(3.3)
Resolution
of
Contact
Discontinuities
13
(a) 1.)
3
OUPA
0.5 -
UR 0 0.0 L 0.0
0
0
0
0.5
c: 1 .o
X
(b)
“L
UR 0
FIG. 4. Arrangement subcell resolution method.
!Ax
X
of cell averages
discontinuity
(b) for
If, under these circumstances, the discontinuous profile is advected tance At to the right, the new cell average values are
a dis-
n+l_
UL
uz
n+l_
-u~-~(u;4-u~).
n+l
=
UM UR If condition
-
(a) and reconstructed
(3.3) is violated
(3.4)
Un,
R
and the discontinuity
moves into the next cell
14
STEPHEN
to the right, the new cell averages n+l
=Un
n+l
=un
UL UM
n+i
UR
F. DAVIS
are
L (3.5)
L
=un .-g-(u;-u;l)
-
(UZ-ug.
In each of these cases the pattern of two constant regions separated by a single intermediate point is preserved. This is usually called a perfect representation of the discontinuity. Conditions (3.3), (3.4), and (3.5) equations. These equations are n+l
UM
=Un
can be combined
into a single set of
M+min[v(u”,-u~)(uT;-U~)]
~~+l=u~+max[O,V(u~-7&)-(7$-U&)], where u = &j&r. After additional
manipulation,
(3.6)
these equations
can be combined
into the form studied by Sweby [17]. F or advection
and put
to the right, this form is (3.7a)
where (3.7b) and Au,-1
(3.7c)
FT.
Space (3.6)
limitations
prohibit
us from carrying
out the manipulations
but
can be put into the form (3.7) with 4 defined by $(T) = min
[&?I.
(3.3)
This result establishes a relationship between subcell resolution and limiter methods. An interesting aspect of (3.8) is that it represents the outer boundary of the region where nonoscillatory solutions are possible. This condition was first discovered by van Leer [18] and has been rediscovered since by many people in many different ways.
Resolution of Contact Discontinuities
15
2.5
20
0
40
60
80
100
X FIG. 5. Computed solution for a periodic contact discontinuity after every 1200 steps until step 12000 using the Godunov-MUSCL method, a Courant number of 0.25, and the Fromm limiter.
Consider
limiters of the form 4(r)
= Kmin[l,r].
(3.9)
Methods based on these limiters are only first order accurate but these limiters contain only one parameter. It is well known that K = 1 causes square wave data to spread and that K = 2 is a nonspreading limiter. Our numerical experiments indicate that K = 1.5 does not spread. These tests were conducted using square wave initial data and periodic boundary conditions. The solution was sampled each time that the wave returned to its initial position. Figures 5 and 6 show typical spreading and nonspreading results. We could not find an exact spreading boundary and we could not detect any significant dependence on the Courant Two commonly used limiters are $(7-) = min[2,2r,
$l
number u.
+ 7-)]
(3.10)
max(1, r)]
(3.11)
and 4(r)
= min[2,2r,
when T > 0 and $ = 0 otherwise.
The limiter (3.10) reduces to the Fromm
16
STEPHEN
0.01 0
20
40
60
80
F. DAVIS
100
X
FIG. 6. Computed solution for a periodic contact discontinuity after every 1200 steps until step 12000 using the Godunov-MUSCL method, a Cow-ant number of 0.25, and the superbee limiter.
method when the solution is smooth and the limiter (3.11) is Roe’s superbee limiter. Both reduce to (3.9) with K = 2 near the corners of an approximate discontinuity but (3.10) spreads and (3.11) does not. Numerical experiments indicate that any,limiter that reduces to a strictly convex combination of 1 and T near T = 1 will spread. Such convex combinations represent linearly stable second order methods. The max(r, 1) term in superbee is less dissipative than any such method. Thus our numerical experiments indicate that this reduced dissipation near r = 1 is necessary to prevent spreading but we do not know how to prove this. We conclude this section by showing some results for the Zalesak [15] ellipse problem. This problem makes all limiter methods look bad. In general, spreading limiters give better results for smooth solutions and nonspreading limiters give better results when discontinuities are present. The Zalesak problem falls somewhere in between. The discontinuous slope and steep profile activate the limiters. Nonspreading limiters try to turn the original profile into a square wave and spreading limiters eventually smear it out beyond recognition. Figures 7 and 8 shown that the Fromm limiter gives more attractive results than superbee but neither solution is acceptable. Zalasak suggests that this problem be solved with very high order methods. This may or may not be necessary but additional work is
Resolution
of Contact Discontinuities
0.0
0.0
20.0
17
40.0
60.0
60.0
100.0
X FIG. 7. Computed solution for a periodic Fromm limiter and a Courant number of 0.5.
half ellipse after 2000 steps using the
2.0
1.5
3
1.0
0.5
0.0
’
0.0
20.0
40.0
60.0
80.0
1
X FIG. 8. Computed solution for a periodic Superbee limiter and a Courant number of 0.5.
half ellipse afer 2000 steps using the
REFERENCES 1
J. J. Quirk,
A Contribution
to the Great
Riemann
Solver
Debate,
ICASE
18
2
3
4 5
STEPHEN F. DAVIS
Report No. 92-64 (1992). B. van Leer, J. L. Thomas, P. L. Roe, and R. W. Newsome, A Comparison of Numerical Flux Formulas for the Euler and Navier-Stokes Equations, 8th AIAA CFD Conference, Honolulu, Hawaii, 1987. A. Harten, P. D. Lax, and B. van Leer, On upstream differencing and Godunov-type schemes for hyperbolic conservation laws, SIAM Rev. 25:3561 (1983). S. Osher and F. Solomon, Upwind difference schemes for hyperbolic systems of conservation laws, Math Comput. 38:339-374 (1982). J. Steger and R. F. Warming, Flux vector splitting of the inviscid gasdynamic equations with applications to finite difference methods, J. Comput. Phys.
6 7 8
9 10
11
12
13
40:263-293
Statist.
14 15
16 17
Comput.
9:445-473
(1988).
T. W. Roberts, The behavior if flux difference splitting schemes near slowly moving shock waves, J. Comput. Phys. 90:141-160 (1990). S. T. Zalesak, A preliminary comparison of modern shock-capturing schemes: Linear advection, in Advances in Computer Methods for Partial Differential Equations, VI (R. Vichnevetsky, R. S. Stepleman, Eds.) IMACS, New Brunswick, N.J. 1987. A. Harten, EN0 schemes with subcell resolution, J. Comput. Phys. 83:148184 (1989). P. K. Sweby, High resolution servation
18
(1981).
B. van Leer, Flux-vector splitting for the Euler equations, in Lecture Notes in Physics, Springer, 1982, pp. 507-512. J. P. Boris and D. L. Book, Flux-corrected transport. I. SHASTA, A fluid transport algorithm that works, J. Comput. Phys. 11:38-69 (1973). B. van Leer, Towards the ultimate conservative difference scheme. I. The quest for monotonicity, in Lecture Notes in Physics, Springer, Berlin, 1973, pp. 163-168. A. Harten, The artificial compression method for computation of shocks and contact discontinuities., Comm. Pure Appl. Math. 30:611-638 (1977). P. L. Roe and M. J. Baines, Asymptotic behavior of some non-linear schemes for linear advection, in Proceedings 5th GAMM Conference on Numerical Methods in Fluid Mechanics (M. Pandolfi, R. Piva, Eds.) 1983. S. K. Godunov, A finite difference method for the numerical computation of discontinuous solutions of the equations of fluid dynamics, Mat. Zb. 47:271290 (1959). A. Harten, High resolution schemes for hyperbolic conservation laws, J. Comput. Phys. 49:357-393 (1983). S. F. Davis, Simplified second-order Godunov-type methods, SIAM J. Sci.
laws,
B. van Leer, tonicity Phys.
SIAM
Towards
the ultimate
and conservation 14(4):361-370
schemes
using flux limiters
J. Numer. Anal.
(1974).
combined
21:995-1011
conservative in a second
for hyperbolic
con-
(1984).
difference order
scheme.
scheme,
II. MonoJ.
Comput.