Shock Waves (1998) 8: 119 125
Some remarks about the resolution of high velocity flows
near low densities
J.P. Cocchi, R. Saurel, J.C. Loraud
Institute Universitaire des SystŁmes Thermiques Industriels, UMR CNRS 6595, Technop�le de Ch�teau-Gombert, 5, rue Enrico Fermi, 13453 Marseille
Cedex 13, France
Received 28 March 1997 / Accepted 18 June 1997
Abstract. High velocity flows which are exposed to strong HLL Riemann solver proposed by Harten et al. (1983) in the
rarefaction waves and creating low densities regions in it treatment of such flows. This method gives good results on
present difficulties and inaccuracies for numerical resolu- the  123 problem except in the variables related to the in-
tion. In particular, variables related to the internal energy ternal energy. Godunov type schemes using exact Riemann
are wrongly evaluated. Use of classical schemes solving the solvers fail, too.
Euler equations in conservative variables introduces signif- These errors are related to the inaccurate determination
icant errors in the determination of temperature. We rec- of flow velocities. Thus, the internal energy obtained from
ommend to employ a non-conservative formulation of the the total and kinetic energies is wrong. This may cause seri-
energy equation. Results found to be more accurate in using ous consequences when chemical transformations depending
the present internal energy formulation. In order to have the on temperature or internal energy are expected to occur in
formulation available for both shock and strong rarefaction the flow, for example combustion or cavitation. An alter-
waves, we propose a hybrid formulation of conservative and native way to fix this defect is to use a non-conservative
non-conservative ones, depending on a shock indicator. The formulation as investigated by Karni (1993).
results are compared with exact solutions and show a sig- The simple strategy proposed here uses a nearly con-
nificant improvement of the accuracy. The method is then servative formulation. Indeed, this formulation is conser-
extended to two-dimensional cases. vative for shock waves while the total energy equation is
transformed into the internal energy equation away of shock
Key words: Conservative and non-conservative formulation, waves. The calculation of internal energy and temperature is
Strong rarefaction, Hybrid scheme improved by this special treatment. The process of switching
from one formulation to the other is achieved by the help of
a shock indicator.
First, we proceed to a one-dimensional analysis of the
 123 problem (Toro 1997). Various conservative schemes
are tested in this problem. Then, we present results ob-
1 Introduction
tained with the non-conservative formulation. Some ele-
ments which allow to extend the method to a hybrid formu-
High velocity flows may contain specific zones where most
lation of conservative and non-conservative ones are given.
of total energy is stored as kinetic energy. Methods usually
This method is then generalised to two-dimensional flows
employed in gasdynamics, use the conservative formulation
with the finite volume formulation of a MacCormack (1969)
of the Euler equations. Indeed, this formulation is currently
scheme. The main advantage of the overall method is that
the best for the calculation of discontinuities, such as shock
only a few modifications are necessary for a classical aero-
waves or contact surfaces. In this formulation, the conser-
dynamics code. A two-dimensional comparison between hy-
vative fluxes include the Hugoniot relationships. When the
brid and conservative formulations is achieved for a gas
contribution to total energy is essentially a kinetic one and
flow around an obstacle. Temperatures obtained by these two
the fluid density is very low, this kind of method produces
computations are very different in the wake of a projectile
wrong solutions. The thermodynamic variables related to in-
where the rarefaction wave is very strong.
ternal energy, such as temperature, are wrongly estimated.
This situation takes place in high velocity flows when a
strong rarefaction wave appears, or for the treatment of wall
2 One-dimensional test problem
boundaries. This problem is referenced as  123 problem in
Toro (1997) and is a test case where conservative methods
The simplest test case to evaluate various computational
fail. Einfeldt et al. (1991) have proposed an extension of the
methods is that of a tube containing a gas which has op-
Correspondence to: R. Saurel posite initial velocities. The initial data consists of two con-
120
Fig. 1. Test problem with two rarefaction waves
stant states Wl(�l, ul, Pl), Wr (�r, ur, Pr) separated by Historically, the first scheme based on the Riemann prob-
a discontinuity at position x = x0 where �l = �r, ul = -ur lem is that of Godunov (1959). It is first order accurate and
and Pl = Pr as shown in Fig. 1. �, u and P represent den- uses an exact Riemann solver. This scheme has been ex-
sity, particle velocity and pressure, respectively. This case tended to the second order by using the concept proposed
is a particular Riemann problem in which an exact solution by van Leer (1979). With second order schemes, we can
is available. It will be used to compare the accuracy of the employ either exact or approximate Riemann solvers. A de-
various methods. tailed description of these schemes is given in Hirsch (1988)
The solution of this test case consists of two rare-fac- or Toro (1997). It has already been mentioned that many con-
tion waves propagating into opposite directions. The density servative schemes failed in this test case (Toro 1995). Here,
between these waves decreases when the initial velocity dif- the problems to test are:
ference increases. The situation represented in Fig. 1 is not
S1: First order Godunov scheme with an exact Riemann
very usual in one-dimensional flows. But this configuration
solver;
arises naturally for the treatment of fixed wall boundaries in
S2: Second order Godunov scheme with an exact Riemann
multi-dimensional flows. So it is interesting to investigate
solver;
the one-dimensional behaviour of classical schemes for the
S3: First order Godunov scheme with the HLLE Riemann
resolution of this specific configuration. To impose the fixed
solver; and
wall boundary V.n = 0, where V designates the velocity
S4: Finite difference MacCormack scheme with second or-
vector and n the normal vector to the wall, one assumes a
der artificial viscosity.
mirror condition at a fictitious node which is -V.n. So the
same configuration as represented in Fig. 1 appears along
The scheme S2 uses the MUSCL Hancock strategy
the normal to the solid wall. In the following we show that
(Quirk 1994) for the second order extension. The scheme
some problems may appear with strong rarefaction waves
S4 employs the HLL solver with the right and left wave
when this kind of fixed wall boundary condition is used.
velocities estimates, given by Einfeldt et al.(1991), based
First, we analyse results obtained by usually employed
on Roe (1981) variables. Remind that the HLLE solver was
conservative methods. We distinguish methods based on the
developped specifically for the treatment of this test case.
Riemann problem which are usually the most accurate for
Einfeldt et al. (1991) showed elsewhere that the original
the solution of the Euler equations, in comparison with finite
Roe scheme failed in this test case. The HLLE solver is pre-
difference methods using artificial viscosity. The equations
sented only in the first order method as recommended by the
in conservative formulation are written as follows:
authors. Moreover, we test equally the MacCormack scheme
"� "�u
with artificial viscosity in S4. This scheme is representative
+ =0, (1)
"t "x
of Lax type methods.
"�u "(�u2 + P ) In Fig. 2, we present the velocity, density and pressure
+ =0, (2)
profiles obtained by the solution of the Euler equations writ-
"t "x
ten in conservative formulation (1) to (3) for the cases S1-
1 1
"�(e + u2) "u(�(e + u2) +P )
2 2
S4. The initial conditions are: P0 = 105 Pa, �0 = 1 kg/m3,
+ =0. (3)
"t "x
u0 = 500 m/s ; and ł = 1.4. For these conditions vacuum
The gas is assumed to be the ideal, and the equation of state cannot be created. According to Godunov (1979), a vacuum
is defined as P =(ł - 1)�e where e is the internal energy. (� = 0) appears if the following inequality is verified:
121
Fig. 3a,b. Entropy and temperature profiles obtained with the conservative
formulation and various numerical schemes
efaction waves, we should obtain no entropy variation. How-
ever, all schemes produce entropy variation. This problem is
known elsewhere and is not peculiar in this test case. On the
other hand, one can notice that the temperature which would
decrease through rarefaction waves, diminishes more or less
depending on schemes used, but always increases in a non-
physical way at the centre of the domain. This increase of
temperature may have serious consequences in some prob-
lems involving, for example, combustion or cavitation. The
Fig. 2a c. Velocity, density and pressure profiles obtained with the conser- first order Godunov scheme predicts a wrong decrease of
vative formulation and various numerical schemes
the temperature inside the rarefaction waves. The second
order Godunov scheme predicts better the flow variables in-
side the rarefaction waves, but produces a strong increase
ul - ur of the temperature at the center of the domain. The HLLE
solver tends to give similar results. The best results are ob-
2 łP
where Uvacuum = (ar + al) and a = .
tained by MacCormack scheme which reduces the increase
ł-1 �
Thus, the drawbacks of the various methods in this test case, of temperature by introducing artificial viscosity. The errors
are not due to the creation of vacuum but to the failure of in the calculation of temperature are related to a wrong es-
the tested schemes. timation of the velocity as shown in Fig. 2. On the plateau
The ratio of the internal energy to the initial kinetic en- where the velocity is zero, all the schemes tested produce
ergy is 2. For flow variables represented in Fig. 2, we can some level of oscillations. These inaccuracies have impor-
notice some errors and oscillations in comparison with the tant consequences in the calculation of internal energy and,
exact solution, but this error is usually observed in other test therefore, in the determination of temperature. In order to re-
cases. In Fig. 3, the profiles of entropy and temperature are duce errors related to the calculation of internal energy and
presented. Since this test case consists of two isentropic rar- temperature, we propose to modify the total energy Eq. (3)
122
the predictor - corrector scheme. For Godunov type schemes,
this term is given by the difference of the velocities solu-
tions of the Riemann problem at interfaces between right
and left cells. Figure 4 represents the velocity, density and
pressure profiles calculated from the system (1, 2, 5). This
non-conservative system is here named hybrid formulation.
The HLLE scheme is no longer used in these simulations.
We have seen that it was not particularly efficient in com-
parison with the other schemes tested. Indeed, it has been
built to fix some inconveniences which appear in the Roe
scheme (Roe 1981). We have conserved only the Godunov
schemes which use the exact Riemann solver. As seen in
Fig. 2, the results of Fig. 4 do not present significant devia-
tions in comparison with the exact solution. The results of
Fig. 5 are to compare with those of Fig. 3. We note always
a non-physical entropy creation, but this entropy creation is
reduced by a factor of about 5. All tested schemes are closer
to the exact solution. The first order Godunov scheme (Go-
dunov 1959) and MacCormack scheme (MacCormack 1969)
present no erroneous increase of temperature at the centre
of the domain. The second order Godunov scheme calcu-
lates more accurately the rarefaction waves but produces
some oscillations in the rarefaction tail and a slight increase
of temperature at the centre of the domain. Modifying the
extrapolation of variables on the cell sides and using the
original idea of van Leer (1979) rather than the Hancock
procedure (Quirk 1994), reduce this problem. However, the
formulation (5) must not be applied everywhere in the flow,
since the Rankine-Hugoniot relations are not verified across
shock waves. In order to illustrate this problem, we present
the temperature profiles obtained in a shock tube problem.
The tube contains initially a gas at rest, the left chamber con-
tains a high pressure gas while the right chamber contains
the same gas at atmospheric pressure:
Pl =107Pa, �l = 1kg/m3, ul = 0; and łl =1.4,
Pr =105Pa, �r = 1kg/m3, ur = 0; and łr =1.4.
In Fig. 6, we present the temperature profiles given by the
MacCormack scheme applied to the Euler equations in non-
conservative formulation, in this shock tube problem. As it
was anticipated, the temperature level is wrongly evaluated.
The non-conservative formulation gives an underestimation
of the temperature level. Thus, all other variables linked to
temperature are erroneous as for example the shock prop-
Fig. 4a c. Velocity, density and pressure profiles obtained with the non- agation velocity as seen in Fig. 6. In order to obtain an
conservative formulation and various numerical schemes with minor mod- efficient formulation of both shock and strong rarefaction
ifications
waves, we propose to combine the conservative and the non-
conservative formulations into a unique one, called  hybrid
formulation . The energy equation is written :
into the internal energy Eq. (5), by subtracting the kinetic
1 1
energy equation.
"[�e + �( �u2)] "[�eu + �( �u2 + P )u]
2 2
+ =
"�e "�eu "u "t "x
+ = -P . (5)
"u
"t "x "x
-(1 - �)P , (6)
"x
The system of Eqs. (1, 2, 5) is non-conservative. Its rigorous
solution by a numerical method is not easy. We propose here where � = 1 designates shock waves while � = 0 otherwise.
to treat the right hand side of (5) as a source term. For the
"u
MacCormack scheme (MacCormack 1969), the term is We use, as a shock indicator, the artificial viscosity of
"x
discretised in the same manner as each differencing step of Richtmeyer and Morton (1954), which is defined as follows:
123
Fig. 6a,b. Shock wave indicator profile used with the hybrid formulation
and temperature profiles obtained with the conservative, non-conservative
and hybrid formulations. Results obtained with the MacCormack scheme
Fig. 5a,b. Entropy and temperature profiles obtained with the non-
conservative formulation and various numerical schemes with minor mod-
ifications
�
If " = ui+1 - ui-1 e" 0, then Sshock = "2,
P
otherwise Sshock =0. (7)
If Sshock > 0.01, then we employ the conservative formu-
lation by imposing � = 1 in (6). Some care must be taken
during the use of the hybrid formulation. At the beginning of
each time step, we compute with the shock indicator (7) the
value of � which allows to know what form of the energy
equation must be solved. For this purpose two different for-
Fig. 7. Control volume &!i,j and normal vector to sides
mulations of the flux vector must be considered. Indeed, the
flux differences must be achieved between conservative vec-
tors or non-conservative vectors only. A node where � = 0
may be neighbour to a node where � = 1. The flux differ- 1) At all mesh points, the vector of unknowns is augmented
ence for nodes where � = 0 must be achieved with non- by the internal energy equation. So vector U contains
conservative fluxes while for nodes where � = 1, the flux density, momentum, total energy and internal energy. The
differences must be obtained with conservative fluxes. flux vector is augmented, too, by the flux corresponding
The new formulation of the energy equation allows to to the internal energy equation. Source terms are included
estimate correctly the temperature at shocks and through for the last equation only;
rarefaction waves as seen in Fig. 6. It has the advantage 2) At the beginning of each time step, the shock indicator
to hold for a large domain of density variations and needs a is computed by (7);
few modifications in an aerodynamic code. 3) The conservative system (1-3) and the non-conservative
Several ways to code the algorithm are possible. A sim- Eq. (5) are solved at all mesh points using a classical
ple way is the following: scheme, centred or upwind. The non-conservative term
124
in the internal energy equation is solved as detailed pre-
viously;
4) For mesh points where � = 1, internal energy is obtained
from the total energy (3) as done classically , and for
others mesh points, internal energy is obtained directly
from the solution of (5).
Instead of adding the internal energy equation to the state
Fig. 8. Wave phenomena in a gas flow
vector U, an alternative is to solve (6). In that case, care must
be taken in order that the flux differences must always be
performed with the same conservative or non-conservative
The corrector step is defined by,
formulation.
"t
n
Ui,j = Ui,j - [Śc(Fi,j) +Bi,jŚc(Vi,j)], (11)
&!i,j
3 Two dimensional extention
where for a vector Ai,j, the function Śc is expressed by
As seen in the one-dimensional case, in order to estimate
n n
Śc(Ai,j) = A .Si+1/2,j + A .Si,j+1/2
i,j i,j
correctly the temperature across shock and strong rarefac-
n n
tion waves, we propose a hybrid formulation of conservative
+A .Si-1/2,j + A .Si,j-1/2.
i-1,j i,j-1
and non-conservative forms of the energy equation. Here,
The final expression of vector U at time tn+1 is obtained by:
we generalise this formulation to two dimensions for the fi-
nite volume MacCormack scheme (MacCormack 1969). We
1
n+1
Ui,j = (Ui,j + Ui,j). (12)
have retained the MacCormack scheme not because of its
2
accuracy. Indeed, it is well known that artificial viscosity
At the beginning of each time step, we determine the value
may produce solutions less accurate than high resolution
of � with the shock indicator which allows to know what
scheme. But for specific applications, only centred schemes
form of energy equation is needed.
can be used. For example, when dealing with very compli-
As seen from one-dimension, both conservative and non-
cated equations of state like in liquid cavitation problems,
conservative flux vectors are necessary during the solution
use of other schemes is very difficult. The system of conser-
of the energy equation. Particularly, at node (i, j) where �=
vation laws in hybrid formulation, written in integral form
0, all flux vectors used in (9), at nodes (i, j + 1), (i, j - 1),
is defined as follows:
(i - 1, j), (i +1, j) must be calculated with �= 0, even if �=
"
1 at these neighbouring nodes.
U d&! + F.dS + B(V.dS) =0, (8)
"t
Another alternative of the solution is the one given in
&! s S
the previous section and consists of solving at all nodes the
where
two forms of the conservative and non-conservative energy
1
U =(�, �u, �v, �e + ��(u2 + v2))T , F =(f, g)T , equation.
2
This hybrid formulation is now compared with the con-
servative one, on a problem of gas flow around an obsta-
1
f =(�u, �u2 + P, �uv, u(�e + �(P + �(u2 + v2))))T ,
cle. A perfect gas of ł= 1.4 enters the left part of the
2
domain under study at a velocity u= 1,500 m/s as shown
in Fig. 8. The fluid is initially at atmospheric conditions
1
g =(�v, �uv, �v2 + P, v(�e + �(P + �(u2 + v2))))T ,
(P0 = 105 Pa,T = 300 K) and has an initial velocity
2
V(u = 1500 m/s, v=0). This problem is solved by the conser-
vative and non-conservative formulations. Each formulation
B =(0, 0, 0, (1 - �)P )T and V =(u, v)T .
is solved by the finite volume MacCormack scheme (Mac-
This equation applied to control volume &!i,j writes:
Cormack 1969, see also Hirsch 1988).
The computational mesh has 158 cells in x direction and
"(Ui,j&!i,j)
+ F.S + Bi,j V.S =0, (9)
50 cells in y direction, and the results are presented for
"t
sides sides
the steady state. Figure 8 presents a magnified view of the
temperatures evolution around the obstacle, obtained by the
where S represents the normal vector to the considered side
two latter calculations. The various states of fluid during
as shown in Fig. 7.
changes of the flow direction are shown in Fig. 8.
Equation (9) is solved by the second order accurate Mac-
A strong bow shock is created in front of the obstacle,
Cormack scheme on a quadrilateral composed mesh. Equa-
and at the corner points, two centred rarefaction waves stand
tion (9) is discretised in two steps. The predictor step writes:
to impose the fluid trajectory to be parallel to the wall of
"t
n
the obstacle. Thus, behind the obstacle the fluid pressure is
Ui,j = Ui,j - [Śp(Fi,j) +Bi,jŚp(Vi,j)], (10)
&!i,j
considerably reduced. Downstream in the wake, when the
fluid impinges the axis, a weak shock wave occurs.
where for a vector Ai,j, the function Śp is expressed by
The comparison of the temperatures calculated from con-
Śp(Ai,j) = An .Si+1/2,j + An .Si,j+1/2
i+1,j i,j+1 servative and hybrid formulations shows some differences
+An .Si-1/2,j + An .Si,j-1/2. in the rarefaction waves at points 1 and 2 as seen in Fig. 9.
i,j i+1,j
125
equation in a non-conservative formulation with the internal
energy only. This new formulation gives better results and
improves the calculation of strong rarefaction waves. How-
ever, this non-conservative formulation cannot provide accu-
rate solutions at shock waves. So, we have combined a con-
servative formulation at shock waves and a non-conservative
formulation elsewhere. This combination is achieved by the
help of a shock indicator providing the location of shock
waves and allowing the switching from conservative to non-
conservative formulation. This hybrid formulation has been
tested in one dimension on a shock tube problem. Then it
is extended to two dimensions and tested on a problem of
gas flow around an obstacle. An evident gain of accuracy
appears by the use of the hybrid formulation.
References
Einfeldt B, Munz CD, Roe PL, Sj�green B (1991) On Godunov-type meth-
ods near low densities. J Comp Phys 92: 273 295
Godunov SK (1959) A difference method for numerical calculation of dis-
continuous solutions of hydrodynamics. Matematicheskii Sbornik, 47:
271 306
Godunov SK, Zabrodine A, Ivanov M, Kraiko A, Prokopov G (1979)
R�solution num�rique des problŁmes multidimensionnels de la dy-
namique des gaz. Editions Mir, Moscow
Fig. 9a,b. Temperature contours obtained with the conservative a and non- Harten A, Lax PD, van Leer B (1983) On upstream differencing and
conservative b formulations
Godunov-type schemes for hyperbolic conservation laws. SIAM Re-
view, 25: 35 61
Hirsch C (1988) Numerical computation of internal and external flows.
As pointed out in the one-dimensional case, the computa- Computation Methods for Inviscid and Viscous Flows. Vol.2. John
Wiley and Sons Ltd, England
tion in conservative form presents a non-physical increase
Karni S (1993) Multicomponent flow calculation by a consistent primitive
of temperature in the rarefaction waves. Indeed, a non-
algorithm. J Comp Phys 112: 31 43
physical increase of temperature is observed at the corner
MacCormack RW (1969) The effect of viscosity in hypervelocity impact
points in Fig. 9a, while the calculation in hybrid formula-
cratering. AIAA paper, 69 354
tion Fig. 9b shows a normal decrease of temperature through
Quirk IJ (1994) An alternative to unstructured grids for computing gas
these waves. The conservative calculations do not yield a dynamics flows around arbitrarily complex two dimensional bodies.
Computers and Fluids, 23 : 125 142
significant decrease of the temperature in the wake near the
Richtmyer RD, Morton KW (1967) Difference Methods for Initial-Value
obstacle. On the other hand, the non-conservative calculation
Problems. John Wiley and Sons, New York
produces a more realistic decrease of the temperature. The
Roe PL (1981) Approximate Riemann solvers, parameter vectors and dif-
comparison of these two formulations shows that the hybrid
ference schemes. J Comp Phys 43: 357 372
formulation of the energy equation allows to estimate more
Toro EF (1997) Riemann solvers and upwind methods for fluid dynamics.
correctly the temperature in this kind of two-dimensional Springer Verlag (to appear)
Toro EF (1995) Some IVP s for which conservative methods fail miserably.
problem.
In: Hafez M (ed) Proc 6th Int Symp Comp Fluid Dynamics. Lake
Tahoe, California, USA, September 4 8, 1995
van Leer B (1979) Toward the ultimate conservative difference scheme
4 Conclusion
V. A second order sequel to Godunov s method. J Comp Phys 32:
101 122
We have pointed out some particular problem occurring
in high velocity flows when strong rarefaction waves are
present. Indeed, the solution of the Euler equations in a con-
a
This article was processed by the author using the LTEX style file pljour2
servative formulation with classical Godunov type methods
from Springer-Verlag.
or MacCormack scheme etc., introduces erroneous results in
regions of low density and particularly a wrong determina-
tion of temperature. This problem being due to a wrong es-
timation of velocity, we have written the energy conservative