systems based on the use of bicharacteristics
Maria Lukacova-Medvid'ova
, Jitka Saibertova y
Abstract
In this paper we present recent results for the bicharacteristic based nite vol-
ume schemes, theso-called nitevolume evolution Galerkin (FVEG) schemes. These
methodswere proposedto solve multi-dimensionalhyperbolicconservationlaws. They
combine the usually conicting design objectives of using the conservation form and
following the characteristics, or bicharacteristics. This is realized by combining the
nite volume formulation with approximate evolution operators, which use bicharac-
teristics ofmulti-dimensional hyperbolicsystem. In thisway all ofthe innitelymany
directions of wave propagationare taken into account. The maingoal of thispaper is
to present a self contained overview on the recent results. We study the L 1
-stability
of thenitevolumeschemesobtainedbydierentapproximationsoftheuxintegrals.
Severalnumericalexperimentspresentedinthelastsectionconrmrobustnessandcor-
rect multi-dimensional behaviourof theFVEG methods.
KEYWORDS: multidimensionalnite volumemethods, bicharacteristics, hyperbolicsystems,
waveequation,Eulerequations
1 Introduction
In general numerical solution of truly multi-dimensional systems of conservation laws is a
challengingtask. Themainreasonisthat even forsmallinitialdatawedonothaveexistence
and qualitativeresultsforthesolutionofmulti-dimensionalRiemannproblems. Inprinciple,
we have two nite volume approaches to overcome this fact. Firstly, the are the so-called
centralnite volumemethods (FVM),whichdoes not use the Riemannproblem,see e.g. [4]
and the references therein. However if no characteristic information is taken into account
they may not provide a satisfactory resolution when small time steps are enforced by the
stability condition. Note that for multi-dimensionalproblems the central schemes have the
CFL stability restrictions strongly less than 1.
Thesecondapproachis basedonaquasi-dimensionalsplittingandonthe useof approxi-
mate solutionto the one-dimensionalRiemannproblem, which structureis well-understood.
If main features, that we want to approximate are just one-dimensional, this approach can
produce good qualitative results. But for complex genuinely multi-dimensional structures,
Departmentofmathematics, HamburgUniversityofTechnology,Schwarzenbergstrae95,Germany,E-
mail: lukacova@tu-harburg.de
y
Institute ofmathematics, Universityof Technology Brno,Technicka2, Czech Republic,E-mail: saiber-
tova@um.fme.vutbr.cz
rious localwave structure resolutions.
Looking back to the literature of the last decade we nd several new genuinely multi-
dimensional methods. For example, the wave propagation algorithm of LeVeque [5], the
methodoftransport(MoT)ofFey[2]anditssimpliedversionofNoelle[13],orthemultistate
FVM of Brio et al. [1]. The latter approach is based on the use of the Kircho formulae
for the waveequationorthe linearized Eulerequations inordertocorrectmulti-dimensional
contributions incornersofthe computationalcells. Infact,our approachissimilartothatof
Brio. However, insteadofthe Kirchhoformulaewhichare explicitintimebut singularover
the soniccircle,weuseadierentmethod. Weworkwithageneraltheoryofbicharacteristics
for linear hyperbolic systems of rst order and derive the so-called approximate evolution
operators. This is the most involved part of the derivation ofour schemes.
The basic idea of the evolution Galerkin schemes (EG), introduced by Morton, see e.g.
[12],isthefollowing. Transportquantitiesareshiftedalongcharacteristicsandthenprojected
onto a nite elementspace. Using the results of Ostkamp [14] we have derived in[7] several
new evolution Galerkin methods for the linear system of the wave equation, which have
better stability properties as well as global accuracy. Their generalization to the second
order EGmethodwasdonein[10]forlineartwo-dimensionalsystems. In[6]wehave studied
two-dimensional Riemann problem for the wave equation system and demonstrated good
accuracyofthe EGschemesaswellascorrectmulti-dimensionalresolutionofobliqueshocks.
In their dissertation [15] Ostkamp proposed a generalization of the EG method to the
nonlinear Euler equations. However, in order to implement her scheme quite tedious cal-
culations of three-dimensional integrals had to be done. It was barely feasible for practical
applications, such as the shallow water equations and the Euler equations, especially for
higher order methods. The decisive step was to take a dierent approach, which yield us
to the nite volume framework. In the so-called nite volume evolution Galerkin (FVEG)
methodstheapproximateevolutionsareusedonlyoncellinterfacestoevaluatethenumerical
uxes.
Note that there is a connection between the FVEG method and the interface centered
MoT ofNoelle. Bothmethodsuse multi-dimensionalevolutiononlyoncell-interfacesinstead
of onwhole computationalcells. This leads tothe crucialsimplication of originalmethods
within the FV framework.
The aim of this paper is to give a self-contained overview on the recent results of the
evolutionGalerkinschemes. Wewanttopresentbasicideasofthetheoryofbicharacteristics,
which are used for the derivation of exact and approximate evolution operators for general
multi-dimensionalhyperbolic systems. We illustrate the application of these techniques on
the Euler equations system. Further interesting applications are, for example, the wave
equation system, the shallow water equations, the magneto-hydrodynamic equations or the
equations of nonlinear elasticity.
Thepaperisorganizedasfollows. InSection2the formulationof thenitevolumeevolu-
tion Galerkin scheme ispresented. The exact integral representation forthe linearized Euler
equations are given in Section 3. In order to apply the approximate evolution operator to
fully nonlinear systems, such as the Euler equations of gas dynamics or the shallow water
equations,rstasuitablelinearizationisneeded. ItisdonebyfreezingtheJacobianmatrices
locallyaroundsuitableconstantstates. Despitethe linearizationproceduretheFVEGmeth-
ods satisfy the entropy condition on sonic rarefaction waves and no entropy x is needed,
see [8], [9]for numericalexperiments. Approximation ofthe so-calledmantleintegralsisdis-
cussed inSection 4. We describe the EG3 approximate evolution operator [7] and following
the linesof ourrecentpaper[9]wepresent newEG5approximate evolutionoperators,which
order scheme obtained by means of a bilinear recovery in space is introduced. We present
here L 1
-stability analysis of the FVM obtained by dierent cellinterface ux integration.
The error analysis of the FVEG methods was presented in [9] and [16]. For linear or
linearized systems it was proved that if a bilinear recovery is used the method is of second
order inspaceand time. In Section6 weillustratethrough numericalexperimentsbehaviour
of the scheme onvarioustest examples.
2 Finite volume evolution Galerkin method
Letbeour two-dimensionalcomputationaldomaincovered by theregularsquare mesh
cells
ij
i 1
2
h;
i+ 1
2
h
j 1
2
h;
j+ 1
2
h
= h
x
i 1
2
;x
i+
1
2 i
h
y
j 1
2
;y
j+
1
2 i
;
where i;j 2Z, and h>0is the mesh size parameter.
In nite volume schemes typically the one-dimensional Riemannproblems in normal di-
rection to the cell interfaces are used to approximate uxes on cell interfaces. Instead of
this dimensional-splittingtechnique weuse in our scheme agenuinely multi-dimensionalap-
proach. Inordertocomputeuxesoncellinterfacesthe valueofapproximatesolutionwillbe
determinedbymeansofasuitableapproximateevolutionoperator. Inthis wayalldirections
of wave propagation are taken into accountexplicitly.
Asanexampleofgeneralhyperbolicconservationlaws letusconsidertheEuler equations
written inthe conservative variables
u
t +f
1 (u)
x +f
2 (u)
y
=0: (2.1)
Here the vector of conservative variables isu :=(;u;v;e) T
and the uxes are
f
1 (u):=
0
B
B
@ u
u 2
+p
uv
(e+p)u 1
C
C
A
; f
2
(u):=
0
B
B
@ v
uv
v 2
+p
(e+p)v 1
C
C
A
;
where e stands for the total energy, i.e. e=p=( 1)+(u 2
+v 2
)=2.
Letusintegrate(2.1)overameshcell
ij
andtime interval[t
n
;t
n+1
]:Applying theGauss
theorem for the uxintegralsyields the followingequality
Z
ij
u(x;y;t)dxd y
tn+1
tn +
Z
tn+1
t
n Z
@
ij (f
1 (u)n
x +f
2 (u)n
y
)dSdt=0; (2.2)
where (n
x
;n
y
) is the unit outer normal to interface of control volume @
ij
. Let U n+1
;U n
denote the piecewise constant functions obtained by the integral averages evaluated attime
t
n+1 or t
n
,respectively. More precisely,
U n
ij
= 1
h 2
Z
ij
u(x;y;t
n
)dxd y:
The crucialpointofthe FVEGschemesisthe use oftheapproximateevolution operatorsfor
the evaluationof uxes at cellinterfaces.
Starting from some initialdata U at time t = 0 and taking into account the fact that
the regular rectangular mesh is used, the nite volume evolution Galerkin method(FVEG)
is recursively dened by means of
U n+1
=U n
1
h Z
t
n+1
tn F
1
(U(t))dt 1
h Z
t
n+1
tn F
2
(U(t))dt: (2.3)
Here F
k
(U(t)) represent the approximations of the physical uxes f
k
(u);k = 1;2, on the
cellinterfaces. Forthetime integrationof uxesthe midpointruleisused. Theintermediate
value of the solution atcell interfaces are calculated by means of the approximate evolution
operatorE
4t=2
,whichhasbeenderivedusingthecharacteristicstheoryofhyperbolicsystems,
cf. Section 3. The approximate uxes (e.g. on vertical edges) have a following form
F
k (U
n+
1
2
)= 1
h Z
h
0 f
k (E
4t=2 R
h U
n
)dS
y
: k=1;2; (2.4)
Analogous formulaholds for horizontal edges. Here R
h
is a recovery operator, which trans-
forms piecewise constant function U toa piecewise bilinear functionR
h
U, cf. Section5.
Ifnorecoveryin(2.4)isusedthewholemethodisofrstorder. Inthiscaseweevaluateall
spaceintegralsexactly. Forthehigherorderschemestheinterfaceintegralsareapproximated
byasuitablenumericalquadrature. Weshouldalsonotethat,if theinterfaceintegralsinthe
rst order scheme are computed exactly we actually compute the cell interface uxes (2.4)
in the followingway
F
k (U
n+
1
2
)=f
k
1
h Z
h
0 E
t=2 U
n
dS
: (2.5)
3 Linearized Euler equations and evolution operator
Inordertoderiveanevolutionoperatorforthe Eulerequationsitissuitabletoworkwith
this system in primitive variables
v
t +A
1 (v)v
x +A
2 (v)v
y
=0; x=(x;y) T
; (3.1)
where v :=(;u;v;p) T
is thevector of primitivevariablesand the JacobianmatricesA
1 (v),
A
2
(v) are given by
A
1 (v):=
0
B
B
@
u 0 0
0 u 0 1
0 0 u 0
0 p 0 u 1
C
C
A
; A
2 (v):=
0
B
B
@
v 0 0
0 v 0 0
0 0 v 1
0 0 p v 1
C
C
A :
Here denotes the density, u and v components of velocity, p pressure and isentropic
exponent; =1:4fordryair. Inwhatfollowswebrieydescribemaintechniqueforderiving
an integral representation (or anexact evolution operator). First, we linearize system (3.1)
by freezingthe Jacobianmatricesat asuitablepoint
~
P =(~x;y;~
~
t). Denote by v~=(;~ u;~ ~v;p)~
the local variables at the point
~
P and by ~c the local speed of sound there, i.e. ~c = q
p~
~ .
Thus, the linearized system (3.1) with frozen constant coeÆcient matrices has the form
v
t +A
1 (~v)v
x +A
2 (~v)v
y
=0; x=(x;y) T
: (3.2)
The eigenvalues of the matrix pencil v)= A
1 v)n
x +A
2 v)n
y
, where n =n () =
(n
x
;n
y )
T
=(cos;sin) T
2R 2
is a unit vector, are
1
=u~cos +v~sin c;~
2
=
3
=u~ cos +v~sin ;
4
=u~cos +v~sin +c;~
and the corresponding linearly independent righteigenvectors are
r
1
=
~
~ c
;cos;sin; ~~c
T
; r
3
=(0;sin; cos;0) T
;
r
2
=(1;0;0;0) T
; r
4
=
~
~ c
;cos;sin;~~c
T
:
Let R(v)~ be the matrix of the right eigenvectors and R 1
(~v) its inverse. Denote by w
the vector of characteristicvariables
w=R 1
(~v)v= 0
B
B
@ 1
2 (
p
~
~c
+ucos+vsin)
p
~ c 2
usin vcos
1
2 (
p
~
~c
+ucos+vsin) 1
C
C
A :
Multiplying system (3.2)by R 1
(~
v) from the left we obtainthe following system written in
characteristicvariables
w
t +B
1 (~v)w
x +B
2 (~v)w
y
=0; (3.3)
where B
k
(~v)=R 1
(~v)A
k
(~v)R(v),~ k =1;2, are transformed Jacobian matrices.
Beinginone spacedimension the system (3.3)reduces to adiagonalsystemconsisting of
separated advection equations. Their exact evolutionoperator reads
w
`
(x;t)=w
`
(x
`
t;0); `=1;:::;4: (3.4)
In the multi-dimensionalcase the system (3.3) will reduce to a diagonal one only if the
Jacobian matrices A
1
; A
2
commute, which is not the case of the two-dimensional Euler
equations. Thus, we rewritethe system (3.3)in the formof the followingquasi-diagonalized
system
w
t +
0
B
@
~
u c~cos 0 0 0
0 u~ 0 0
0 0 u~ 0
0 0 0 u~+~ccos 1
C
A w
x +
0
B
@
~
v c~sin 0 0 0
0 v~ 0 0
0 0 v~ 0
0 0 0 ~v+c~sin 1
C
A w
y
=S
(3.5)
with
S = 0
B
B
B
@
1
2
~ c(sin
@w
3
@x
cos
@w
3
@y )
0
~ csin(
@w1
@x
@w4
@x
) ~ccos(
@w1
@y
@w4
@y )
1
2
~
c( sin
@w3
@x
+cos
@w3
@y )
1
C
C
C
A :
Each characteristic variable w
`
; ` = 1;:::;4; is evolved in time along the corresponding
bicharacteristic curve x
`
dened by
dx
`
dt
=b
``
(n):=(b 1
``
;b 2
``
) T
; (3.6)
t
x y
P 0
Q
l ()
Figure1: Bicharacteristics along the Mach cone through P and Q
` ().
where B
1
=(b 1
jk )
j;k=1;::;4 , B
2
=(b 2
jk )
j;k=1;::;4
. The set of all bicharacteristicscreates amantle
of the so-called Mach cone, see Figure 1. In order to obtain the exact evolution of each
characteristic variable w
`
we integrate the `-th equation of the system (3.5) from the apex
P =(x;y;t+t)down tothe corresponding footpointQ
` ():
Q
1
() = (x (~u c~cos)t;y (~v ~csin)t;t);
Q
2
= Q
3
=(x u4t;~ y vt;~ t);
Q
4
() = (x (~u+~ccos)t;y (~v+~csin)t;t):
Multiplyingtheresultingsystemfromtheleftbythematrix R(~
v)yieldstheexactintegral
equations,i.e. theexactevolutionoperators,fortheprimitivevariablesofthelinearizedEuler
equations (3.2).
(P) = (Q
2 )
p(Q
2 )
~ c 2
+ 1
2 Z
2
0 h
p(Q
1 )
~ c 2
~
~ c
u(Q
1 )cos
~
~ c
v(Q
1 )sin
i
d
~
~ c
1
2 Z
2
0 Z
t+4t
t
S(;
~
t;)d
~
td; (3.7)
u(P) = 1
2 Z
2
0 h
p(Q
1 )
~
~c
cos+u(Q
1 ) cos
2
+v(Q
1
)sincos i
d
+ 1
2 Z
2
0 Z
t+4t
t
cosS(;
~
t;)d
~
t d+ 1
2 u(Q
2 )
1
2~ Z
t+4t
t p
x (Q
2 (
~
t))d
~
t;
(3.8)
v(P) = 1
2 Z
2
0 h
p(Q
1 )
~
~c
sin+u(Q
1
)cossin+v(Q
1 )sin
2
i
d
+ 1
2 Z
2
0 Z
t+4t
t
sinS(;
~
t;)d
~
td+ 1
2 v(Q
2 )
1
2~
Z
t+4t
t p
y (Q
2 (
~
t))d
~
t;
(3.9)
p(P) = 1
2 Z
2
0 [p(Q
1
) ~~cu(Q
1
)cos ~~cv(Q
1
)sin]d
~
~c 1
2 2
0
t+4t
t
S(;
~
t;)d
~
td; (3.10)
where =(x (~u ~ccos)(t+4t
~
t);y (~v c~sin)(t+4t
~
t));and the so-calledsource
term S is given inthe following form
S(x;t;):=~c[u
x
(x;t;)sin 2
(u
y
(x;t;)+v
x
(x;t;))sincos+v
y
(x;t;)cos 2
]:
4 Approximate evolution operators
Theexactintegralequations(3.7)-(3.10)forthesolutiontothe linearizedEuler equations
(3.2) willbeabasis for ourfurther numericalapproximations. Theabove exact integralrep-
resentation (3.7)-(3.10) is implicit in time. In order to derive a numerical scheme, which is
explicit intime, timeintegralsofthe source terms,the so-calledmantleintegrals,havetobe
approximated with suitablenumerical quadratures.
4.1 Approximate evolution operator EG3
As in [7] the integrals of the source term with respect to time will be approximated by the
rectangle rule. Thus, we would need to evaluate derivatives of the velocity components at
time t. However, it was proved in [7, Lemma 2.1], that the integrals of the source S can be
simplied through integration by parts, which yields
t Z
2
0
S(t;)d= Z
2
0 [u
Q
cos+v
Q
sin]d: (4.1)
Analogouslywecan derive
t Z
2
0
S(t;)sind= Z
2
0 [2u
Q
sincos+v
Q (2sin
2
1)]d (4.2)
and
t Z
2
0
S(t;)cosd= Z
2
0 [u
Q (2cos
2
1)+2v
Q
sincos]d : (4.3)
Furthertheintegralsin(3.8)and(3.9)involvingp
x andp
y
needtobereplacedbyintegrals
over the cone mantle. This is done by using the Taylor expansion
p
x (P
0
)=p
x
(Q)+O(jP 0
Qj)=p
x
(Q)+O(t):
Then we use the rectangle rule in time and the Gauss theorem on a sonic circle O =
f(x;y);x 2
+y 2
~c 2
t 2
g
1
2~ Z
t+4t
t p
x (P
0
(
~
t))d
~
t =
4t
2~
p
x (P
0
)+O(t 2
)
=
1
2~~c 2
4t Z
O p
x
(Q)d xd y+O(t 2
)
=
1
2~~c 2
4t I
p(Q) d y+O(t 2
)
=
1
2~~c Z
2
0
p(Q) cos d+O(t 2
): (4.4)
y
1
2~ Z
t+4t
t p
y (P
0
(
~
t)) d
~
t = 1
2~~c Z
2
0
p(Q) sin d:
We have generated the approximate evolution operator for the Euler equations, which we,
analogously asfor the wave equation in [7], callthe EG3operator.
Approximate evolution operator EG3
(P) = (Q
2 )+
1
2 Z
2
0 h
(Q
1 )
2
~
~ c
u(Q
1
)cos+v(Q
1 )sin
i
d+O(t 2
); (4.5)
u(P) = 1
2 u(Q
2 )+
1
2 Z
2
0 h
p(Q
1 )
~
~c
cos+u(Q
1
) 3cos 2
1
+3v(Q
1
)sincos i
d
+O(t 2
); (4.6)
v(P) = 1
2 Z
2
0 h
p(Q
1 )
c
sin+3u(Q
1
)sincos+v(Q
1
) 3sin 2
1
i
d
+O(t 2
); (4.7)
p(P) = 1
2 Z
2
0 h
p(Q
1
) 2~~c
u(Q
1
)cos+v(Q
1 )sin
i
d+ O(t 2
): (4.8)
Note that other integral rules for the time approximation of the source terms can be
used as well. These lead to the approximate evolution operators EG1, EG2 and EG4. The
notation is used in analogy to [7], [17]. We decided to present here only the operator EG3,
which yieldsthe most accurate numericalapproximation among EG1- EG4.
Onthe otherhandtheaboveFVEGmethodssueredfromtherestrictivestabilitylimits.
Forexample,for the FVEG3scheme a typicalCFL stability limitwas 0.63 and 0.56for rst
and second order schemes, respectively. The CFL number is dene as CFL = maxfjuj+
c;jvj+cg 4t
h
. Note that integrals around the sonic circle, i.e.
R
2
0
d, are evaluated for the
piecewise constant or piecewise bilinear data exactly. Therefore, obviously the only step
where the stability could be reduced was the time approximation of the mantle integral.
For example, in the rst order EG scheme we work with piecewise constant data, in which
case a discontinuity cuts through the cone mantle. Naturally, classicalquadratures, such as
the rectangle or the trapezoidal rule, which were used for the EG1 - EG4 schemes, cannot
correctly reproduce integration of discontinuous data.
4.2 Approximate evolution operator EG5
Since the approximate evolution operators EG1-EG4 did not provide full stability up to
the CFL number 1 a more appropriate numerical quadratures for time integration along
the mantle of the Mach cone have been constructed in [9]. For one-dimensional advection
equation it is known that any scheme, which is stable up to CFL = 1 has to reproduce
the exact solution at CFL = 1. In [9] the following design principle was proposed for the
wave equation system. Consider plane wave data which are parallel to one of the spatial
axis. Fora rst order scheme these are taken tobepiecewise constant, i.e. one-dimensional
reproduce the exact solutionatthe apex ofthe bicharacteristiccone centered atthe original
discontinuity. When considering slopes for second order schemes we derive approximate
evolutionoperatorsfor the slopesthat againreproduce the solutionfor continuous piecewise
linear data exactly atapex of the bicharacteristic cone centered atthe kink of such data.
In the dissertation of the second author [16] the rigorous derivation of the approximate
evolution operator for the Euler equations using the above design principles has been done.
Wehavedecidednottopresenttheseratherlengthycalculationshere,sinceitwillenlargethe
size of the papersubstantially. Instead we onlypresent the formulationsof the approximate
evolution operator EG5 for both piecewise constant and piecewise bilinear data. Numerical
experimentspresentedin[9],[11],[16]demonstratethattheFVEG5methodisapproximately
three times moreaccurate than FVEG3methodand has the CFLlimitclose to1. Theoreti-
calerroranalysisoftheFVEG5schemeforlinearandlinearizedsystemshasbeendonein[16].
Approximate evolution operator E const
for piecewise constant data
(P) = (Q
2 )
p(Q
2 )
~ c 2
+ 1
2 Z
2
0
p(Q
1 )
~ c 2
~
~ c
u(Q
1
)sgn(cos)
~
~ c
v(Q
1
)sgn(sin)
d
(4.9)
u(P) = 1
2 Z
2
0
p(Q
1 )
~
~c
sgn(cos)+u(Q
1 )
1
2 +cos
2
+v(Q
1
) sin cos
d(4.10)
v(P) = 1
2 Z
2
0
p(Q
1 )
~
~c
sgn(sin)+u(Q
1
) cos sin+v(Q
1 )
1
2 +sin
2
d (4.11)
p(P) = 1
2 Z
2
0 h
p(Q
1
) ~~cu(Q
1
)sgn(cos) ~~cv(Q
1
)sgn (sin) i
d (4.12)
Approximate evolution operator E bil in
for piecewise bilinear data
(P) = (Q
2 )+
1
4 Z
2
0 1
~ c 2
[p(Q
1
) p(Q
2 )]d
1
Z
2
0
~
~ c
[u(Q
1
)cos+v(Q
1
)sin]d
+O(t 2
); (4.13)
u(P) = u(Q
2 )
1
Z
2
0
p(Q
1 )
~
~c
cosd
+ 1
4 Z
2
0 h
3(u(Q
1
)cos+v(Q
1
)sin)cos u(Q
1 )
1
2 u(Q
2 )
i
d +O(t 2
);
(4.14)
v(P) = v(Q
2 )
1
Z
2
0
p(Q
1 )
a
sind
+ 1
4 Z
2
0 h
3(u(Q
1
)cos+v(Q
1
)sin) sin v(Q
1 )
1
2 v(Q
2 )
i
d+O(t 2
);
p(P) = p(Q
2 )+
1
4 Z
2
0
[p(Q
1
) p(Q
2 )]d
1
Z
2
0
~
~c[u(Q
1
)cos+v(Q
1
)sin]d
+O(t 2
): (4.16)
5 Second order schemes
In the rst order schemes no reconstruction is used and the exact solution u(x;y;t
n ) is
approximated by the piecewise constant function U n
. The integrals along the Mach cone
and overeach cellinterface are computed exactly, i.e. no approximation is used. Ifwe want
to obtain higher order nite volume EGschemes, wehave to replace the piecewise constant
functionU n
byahigherdegreepolynomial. Onepossibilitytoobtainthesecondorderscheme
is to apply oneach mesh cell
ij
; i;j 2Z the following discontinuous bilinearrecovery
R D
h U
n
ij
=
1+
(x x
i )
h
x
2
y Æ
x +
(y y
j )
h
2
x
y Æ
y +
(x x
i
)(y y
j )
h 2
x
y Æ
x Æ
y
U n
ij
;
where
x
v(x)= 1
2
v
x+ h
2
+v
x h
2
,Æ
x
v(x)=v
x+ h
2
v
x h
2
,ananalogous
notation is used for y-direction. Note, that this bilinear recovery is constructed in such a
way, that itis conservative, i.e.
1
h 2
Z
ij R
h U
n
dxd y = 1
h 2
Z
ij U
n
dxd y:
Another possibility is to use continuous, but non-conservative bilinear recovery. In this
case wehave oneach mesh cell
R C
h U
n
ij
=
2
x
2
y +
(x x
i )
h
x
2
y Æ
x +
(y y
j )
h
2
x
y Æ
y +
(x x
i
)(y y
j )
h 2
x
y Æ
x Æ
y
U n
ij :
(5.1)
Inour numericalexperimentswehaveused forthe FVEG1-FVEG4schemesthe conser-
vativebilinearrecoveryR D
h
. ThesecondorderFVEG5schemeisconstructedasacombination
of theapproximateevolutionoperatorE bil in
4
,whichevolvesthe continuousbilineardata,and
the approximate evolutionoperatorE const
4
, which corrects the evolutionof the constantpart
in order to preserve conservativity of the cell interface values U n+
1
2
, cf. also [9]. Thus for
FVEG5 scheme weuse aformula
U n+
1
2
=E const
(1
2
x
2
y )U
n
+E bil in
R
C
h U
n
:
It is a well-known fact that the higher order methods can suer from oscillations near
discontinuities. In order to avoid developing oscillationsin solution we controllgradients of
the recoveredfunctionsbyalimiter. Therearemanypossibilitiestochooseasuitablelimiter.
In our numericalcomputationswehave used the so-called minmodlimiter, see [5], [16].
As we mentioned above in the second order scheme we compute exactly only the Mach
cone integrals. Integralovercellinterfacesarecomputed bysome suitablenumericalquadra-
tures. In the next subsection we study the L 1
-stability for these quadrature rules.
5.1 L -stability analysis of the cell interface ux integrals
Natural quadrature points on cell interface are vertices, used for the trapezoidal rule and
midpoints used in the midpoint rule. Combining vertices and midpoints yields Simpson's
rule. In what follows we consider these three quadrature rules for approximation of the cell
interface integrals.
Letus consider a scalar linear advection equationin two-dimensions
U
t +uU~
x +~vU
y
=0; (5.2)
where u,~ ~v > 0 are positive constant velocities. We can write nite volume scheme for the
equation (5.2) in the following way
h 2
U n+1
ij U
n
ij
+ht
F n+
1
2
i+
1
2 j
F n+
1
2
i 1
2 j
+G n+
1
2
ij+
1
2 G
n+
1
2
ij 1
2
=0; (5.3)
where h> 0 is the mesh size parameter. The exact evaluation of the uxes at time t
n+
1
2
=
t
n +
1
2
t gives
F n+
1
2
i+
1
2 j
= u~
"
(~uth 1
2
~ u~vt
2
)U n
ij +
1
2
~ u~vt
2
U n
i 1j
~ uth
#
= u~
1
~ vt
2h
U n
ij +
~ vt
2h U
n
i 1j
; (5.4)
F n+
1
2
i 1
2 j
= u~
1
~ vt
2h
U n
i 1j +
~ vt
2h U
n
i 1j 1
; (5.5)
G n+
1
2
ij+
1
2
= v~
1
~ ut
2h
U n
ij +
~ ut
2h U
n
ij 1
; (5.6)
G n+
1
2
ij 1
2
= v~
1
~ ut
2h
U n
ij 1 +
~ ut
2h U
n
i 1j 1
: (5.7)
Then the nitevolume scheme (5.3) becomes
U n+1
ij
= U n
ij
~ ut
h
x
~ u~vt
2
2h 2
x
y
U n
ij
~ vt
h
y
~ u~vt
2
2h 2
x
y
U n
ij
=
1
~ ut
h
x
1
~ vt
h
y
U n
ij
; (5.8)
where
x U
ij ,
y U
ij
, are the backward dierences in x-, y-direction, respectively, i.e. e.g.
x U
ij
= U
ij U
i 1j
. This scheme is a tensor product of the one-dimensional upwind
schemes and it is well-known that it is monotone and stable up to CFL 1, where we
dene CFL:=maxf
~ ut
h
;
~ vt
h
g. More precisely, under the above CFL condition let us show
the L 1
-stability of the scheme, i.e. kU n+1
k
1 kU
n
k
1
. For the discrete grid function U n
we
dene the discrete L 1
-norm by
kU n
k
1
=h 2
X
i;j2Z jU
n
ij
j: (5.9)
We can rewrite equation(5.8) in the following way
U n+1
ij
=
1
~ ut
h
1
~ vt
h
U n
ij +
~ vt
h
1
~ ut
h
U n
ij 1
+
~ ut
h
1
~ vt
h
U n
i 1j +
~ u~vt
2
h 2
U n
i 1j 1 :
From the denition of the L -norm we have
kU n+1
k
1
= h 2
X
i;j jU
n+1
ij j
h
"
X
i;j
1
~ ut
h
1
~ vt
h
U n
ij
+
X
i;j
~ vt
h
1
~ ut
h
U n
ij 1
+ X
i;j
~ ut
h
1
~ vt
h
U n
i 1j
+
X
i;j
~ u~vt
2
h 2
U n
i 1j 1
#
:
Thus, if the followingconditions hold
1
~ ut
h
1
~ vt
h
0
~ vt
h
1
~ ut
h
0
~ ut
h
1
~ vt
h
0
~ u~vt
2
h 2
0; (5.10)
then the scheme (5.8) is stable, since wecan take terms from (5.10) infront of the absolute
value. We get
kU n+1
k
1
h
2
"
1
~ ut
h
1
~ vt
h
X
i;j
U n
ij
+
~ vt
h
1
~ ut
h
X
i;j
U n
ij 1
+
~ ut
h
1
~ vt
h
X
i;j
U n
i 1j
+
~ u~vt
2
h 2
X
i;j
U n
i 1j 1
#
=
1
~ ut
h
1
~ vt
h
kU n
k
1 +
~ vt
h
1
~ ut
h
kU n
k
1
+
~ ut
h
1
~ vt
h
kU n
k
1 +
~ u~vt
2
h 2
kU n
k
1
= kU n
k
1 :
Note that conditions (5.10) can be rewritten inthe following form
0
~ ut
h
1 and 0
~ vt
h 1;
which corresponds tothe CFL condition given above.
Now we would like to nd stability conditions for the trapezoidal rule approximation
interface integrals of uxes. If we replace the exact uxes (5.4) - (5.7) by the trapezoidal
rule, we get
F n+
1
2
i+
1
2 j
=
~ u
2 U
n
ij +U
n
ij 1
; F
n+
1
2
i 1
2 j
=
~ u
2 U
n
i 1j +U
n
i 1j 1
;
G n+
1
2
ij+
1
2
=
~ v
2 U
n
ij +U
n
i 1j
; G
n+
1
2
ij 1
2
=
~ v
2 U
n
ij 1 +U
n
i 1j 1
:
Then the nitevolume scheme (5.3) has the form
U n+1
ij
= U n
ij
~ ut
2h [
x +
x
x
y ]U
n
ij
~ vt
2h [
y +
y
x
y ]U
n
ij
=
1
~ ut
h
x
1
~ vt
h
y
U n
ij +
~ ut
2h
~ u~vt
2
h 2
+
~ vt
2h
x
y U
n
ij :
(5.11)
Letus showthat this scheme isstable only if
h
=
h
. Wecan rewrite(5.11) as
U n+1
ij
= 1
2
2
~ ut
h
~ vt
h
U n
ij +
1
2
~ ut
h +
~ vt
h
U n
ij 1
+ 1
2
~ ut
h
~ vt
h
U n
i 1j +
1
2
~ ut
h +
~ vt
h
U n
i 1j 1
and hencethe L 1
-normcan be bounded from abovein the followingway
kU n+1
k
1
= h 2
X
i;j jU
n+1
ij j
h
2
2
"
X
i;j
2
~ ut
h
~ vt
h
U n
ij
+
X
i;j
~ ut
h +
~ vt
h
U n
ij 1
+ X
i;j
~ ut
h
~ vt
h
U n
i 1j
+
X
i;j
~ ut
h +
~ vt
h
U n
i 1j 1
#
: (5.12)
In order toobtainstabilityof the scheme the followingfour conditionshaveto holdsimulta-
neously
~ ut
h +
~ vt
h 2;
~ ut
h
~ vt
h 0
~ ut
h
~ vt
h 0;
~ ut
h +
~ vt
h 0:
Thusweneedtorequirethat
~ ut
h
=
~ vt
h ,
~ ut
h
1and
~ vt
h
1. Ifweputtheseconditions
into(5.12), we obtain
kU n+1
k
1
h 2
2
"
2
~ ut
h
~ vt
h
X
i;j
U n
ij
+
~ ut
h +
~ vt
h
X
i;j
U n
ij 1
+
~ ut
h
~ vt
h
X
i;j
U n
i 1j
+
~ ut
h +
~ vt
h
X
i;j
U n
i 1j 1
#
= 1
2
2 2
~ ut
h
kU n
k
1 +
2
~ ut
h
kU n
k
1
= kU n
k
1 :
We have shown, that the trapezoidalrule yieldsa stable nite volume scheme onlyif u~=v~
and hence itissuited onlyforspecial situations, e.g. for theapproximationof edgeintegrals
for the wave equation without advection, e.g. u~ = 0= ~v. For other systems with arbitrary
advections,such as the Euler equations,this quadrature rule yieldsunconditionaly unstable
scheme.
In an analogous way wecan derive stability condition forthe FV scheme using the mid-
pointrule for ux integrals, which reads
~ ut
h +
~ vt
h
1; 0
~ ut
h
; 0
~ vt
h :
If Simpson's rule is used the followingstability conditions have to be satised
1
5
~ vt
h
~ ut
h 5
~ vt
h
and 0
~ ut
h +
~ vt
h
6
5 :
~
u;v~ Simpson's rule leads to in general the most stable discretization of the ux integrals.
This is the numerical quadrature that we use in our numerical experiments in Section 6.
Note that this L 1
-stability analysis holdfor general two-dimensional FVschemes.
6 Numerical experiments
In this section, through numerical simulations, we illustrate performance of the FVEG
methods for the linear wave equation system and for fully nonlinear Euler equations of gas
dynamics.
6.1 Water waves propagation
It is known that the approximation of circular waves on rectangular meshes can cause diÆ-
culties. Particularly,if thedimensionalsplittingapproachisusedthespuriousmeshoriented
structures can be developed, see e.g. [5], [7], [13]and the references therein.
We consider the wave equationsystem
v
t +A
1 (v)v
x +A
2 (v)v
y
=0; x=(x;y) T
; (6.1)
where v :=(;u;v) T
. The Jacobian matricesA
1 (v);A
2
(v) are given by
A
1 (v):=
0
@
0 c 0
c 0 0
0 0 0 1
A
; A
2 (v):=
0
@
0 0 c
0 0 0
c 0 0 1
A
;
where c is a sound speed; see [9] for approximate evolution operators of (6.1). Consider the
following initialdata modeling a pointwise disturbance
(x;0)= cexp ( 15x 2
15y 2
);
u(x;0)=0=v(x;0):
The computationaldomain [ 3;3][ 3;3] is divided into 100100 cells. The solution
obtained by the second order FVEG scheme at dierent times from T = 0:2 until T = 8:0
is shown inFigure 2. We can notice well resolved symmetric circularwave. As time evolves
the wave propagates and is being reected from the left boundary. The Mach stem which
is evolving behind the main wave can be recognized in time t = 8:0. This problem can be
considered asa model for apointwise disturbance of water surface, e.g. asit occurs when a
stone is thrown into alake.
As mentioned above, we set reected boundary conditions on the left vertical boundary
and absorbing boundary conditions elsewhere. In numerical experiments presented in this
paper we have implemented absorbing boundary conditions by linear extrapolation of all
quantities to the so-called ghost cells, which are adjacent to the boundary of the computa-
tional domain. Thus, we have, e.g. for the pressure wave ,
1j
=
0j
;
2j
=
1j
; j 2Z;
where
1j ,
2j
and
0j ,
1j
are the ghost cells and the cells belongingto , respectively.
Note that due to the secondorder method two layers of the ghost cells are needed. The
reected boundary conditions are easily modeled by reecting the interior data across the
boundary
u
1j
= u
0j
; u
2j
= u
1j
; j 2Z:
Other quantities, i.e. ;v, are extrapolated.
6.2 Interaction between circular shocks and reected waves
In this example weconsider again the wave equationsystem (6.1) with the following discon-
tinuous initialdata
=1; u=0; v =0; kxk<0:4
=0; u=0; v =0; else.
Thecomputationaldomain[ 1;1][ 1;1]isdividedinto200200cells. Weimplemented
reected boundaryconditionsonthe vertical boundaries andabsorbing boundaryconditions
on the horizontal ones. In Figure 3 the pressure wave distribution at dierent times T =
0:3;1:0;1:3isdepicted. Wecan notice awell-resolved circularshock traveling away fromthe
centerofthecomputationaldomain. Astimeevolvestheshockreachestheverticalboundaries
and isreectedintothecomputationaldomain. Duetothelinearmodelinteractionsbetween
the linearcircular shock and reected waves can be observed verywell.
6.3 Static disc problem for the Euler equations
This is a two-dimensional problem with a circular symmetry. We consider nonlinear hyper-
bolic systems of Euler equations (3.1). The computational domain is [ 1;1][ 1;1] and
the boundary conditions are periodical. Here the initialconditions are following
(x;y;0) =
3 if (x 2
+y 2
)0:5
1; otherwise;
u(x;y;0) = 0;
v(x;y;0) = 0;
p(x;y;0) = 1:
Thus wehave alongthe disc boundary contact discontinuity. Theexact solutionisthe same
asthe initialcondition,i.e. disc withradius0.5. Wecomputethis problemup totime t=10
on dierent meshes (1616, 3232, 6464, 128128). The initialconditions for cells
lyingontheboundaryof thecirculardiscare implementedbyapproximateweightedintegral
averagesof initialdata. InFigure4weshowisolinesofdensity computedondierentmeshes
by the second order nite volume schemes EG3 and EG5. In middle column the weighted
initialdata, which weused, are plotted. This exampledemonstrates thatthe steady contact
discontinuity is resolved in a correct way preserving the multi-dimensional phenomena as
well. Note that Kroger and Noellereported in[3], that their MoT-ICE method failson this
example producing totallyincorrect resolution.
Acknowledgements
This researchhas been supported underthe VolkswagenStiftung grantI 76859,bythe grant
of the Czech grant Agency GA
CR 201/03/0570, by the MSM 262100001, by the Deutsche
by the EC as contract HPRN-CT-2002-00282. Authors gratefully acknowledge these sup-
ports.
References
[1] M. Brio, A.R. Zakharian,G.M. Webb, Two-dimensional Riemann solver for Euler equa-
tions of gas dynamics, J. Comput. Phys.167(1), (2001) 177-195.
[2] M.Fey,Multidimensionalupwinding,PartII. Decomposition ofthe Euler equations into
advection equations, J.Comp. Phys. 143, (1998) 181-199.
[3] T. Kroger, S. Noelle Numericalcomparison of the Method of Transport to a standard
scheme. IGPM Preprint no. 229, submited toComp. Fluids,2003
[4] A.Kurganov,S.Noelle,G.Petrova,Semidiscretecentral-unpwindschemesforhyperbolic
conservation laws and Hamilton-Jacobi equations, SIAM J. Sci. Comput. 23(3), (2001)
707-740.
[5] R.J. LeVeque, Wave propagation algorithms for multi-dimensional hyperbolic systems,
J.Comp.Phys. 131, (1997)327-353.
[6] J. Li, M. Lukacova-Medvid'ova, G.Warnecke, Evolution Galerkin schemes for the two-
dimensional Riemmanproblems, Discrete and Continuous Dynamical Systems (A), 9(3),
(2003) 559{578.
[7] M.Lukacova - Medvid'ova, K.W. Morton,G. Warnecke, EvolutionGalerkin methodsfor
hyperbolic systems in two space dimensions, MathComp.69, (2000) 1355{1384.
[8] M. Lukacova-Medvid'ova, K.W. Morton, G.Warnecke.Finite volume evolutionGalerkin
methods for Euler equations of gas dynamics, Int. J. Numer. Meth. Fluids 40, (2002)
425-434.
[9] M. Lukacova-Medvid'ova, K.W. Morton, G.Warnecke.Finite volume evolutionGalerkin
(FVEG) methods for hyperbolicproblems, accepted to SISC, (2003).
[10] M. Lukacova - Medvid'ova, G. Warnecke, Lax-Wendro type second order evolution
Galerkin methods for multidimesnional hyperbolic systems, East-West Journal 8(2),
(2000) 127{152.
[11] M.Lukacova-Medvid'ova,G.Warnecke,Y.Zahaykah,OntheStabilityoftheEvolution
Galerkin Schemes Applied to a Two-dimensional Wave Equation System, submitted to
Numer. Math. (2004).
[12] K.W.Morton,Ontheanalysisofnitevolumemethodsforevolutionaryproblems,SIAM
J. Numer. Anal.35(6), (1998)2195-2222.
[13] S.Noelle, The MOT-ICE: anew high-resolutionwave-propagationalgorithmfor multi-
dimensionalsystems of conservativelawsbased onFey'smethodof transport, J.Comput.
Phys. 164, (2000) 283{334.
[14] S.Ostkamp, MultidimensionalcharacterisiticGalerkin schemesand evolution operators
for hyperbolic systems, Math. Meth. Appl. Sci. 20, (1997)1111-1125.
for hyperbolic systems, PhD thesis (University Hannover, 1995).
[16] J. Saibertova, Genuinely multi-dimensional nite volume schemes for systems of con-
servation laws, PhD thesis (Technical University Brno, 2003).
[17] Y.Zahaykah, Evolution Galerkin schemesand discrete boundary conditions formultidi-
mensional rst order systems, PhD thesis (University of Magdeburg, 2002).
−0.5 0 0.5
−0.5 0 0.5
FV EG3
16x16
−0.5 0 0.5
−0.5 0 0.5
32x32
−0.5 0 0.5
−0.5 0 0.5
64x64
−0.5 0 0.5
−0.5 0 0.5
128x128
−0.5 0 0.5
−0.5 0 0.5
Initial data
−0.5 0 0.5
−0.5 0 0.5
−0.5 0 0.5
−0.5 0 0.5
−0.5 0 0.5
−0.5 0 0.5
density
−0.5 0 0.5
−0.5 0 0.5
FV EG5
−0.5 0 0.5
−0.5 0 0.5
−0.5 0 0.5
−0.5 0 0.5
−0.5 0 0.5
−0.5 0 0.5
Figure 4: Static disc problem;plots of density isolineson dierent meshes.
