<< . .

( 23)

. . >>

ENO-RF. More speci¬cally, we use »i0 = f (ui0 ) and »i0 +1 = f (ui0 +1 ) to
spatial terms, we update the entire equation in time with a method-of-lines
decide whether there are sonic points in the vicinity. If »i0 and »i0 +1 agree
approach using, for example, a TVD Runge-Kutta method.
in sign, we use the ENO-Roe discretization where »i0 +1/2 is taken to have
the same sign as »i0 and »i0 +1 . Otherwise, we use the ENO-LLF entropy
¬x discretization given below. Note that ENO-LLF is applied at both ex-
14.5 Systems of Conservation Laws
pansions where »i0 < 0 and »i0 +1 > 0 and at shocks where »i0 > 0 and
»i0 +1 < 0. While this adds extra numerical dissipation at shocks, it is not
In general, a hyperbolic system will simultaneously contain a mixture of
harmful, since shocks are self-sharpening. In fact, this extra dissipation pro-
processes: smooth bulk convection and wave motion, and discontinuous
vides some viscous regularization which is especially desirable in multiple
processes involving contacts, shocks, and rarefactions. For example, if a gas
spatial dimensions. For this reason, authors sometimes use the ENO-LLF
in a tube is initially prepared with a jump in the states (density, velocity,
method everywhere as opposed to mixing in ENO-Roe discretizations where
and temperature) across some surface, as the evolution proceeds in time
the upwind direction is well determined by the eigenvalues ».
these jumps will break up into a combination of shocks, rarefactions, and
The ENO-LLF discretization is formulated as follows. Consider two prim-
itive functions H + and H ’ . We compute a divided di¬erence table for each contacts, in addition to any bulk motion and sound waves that may exist
or develop.
of them, with their ¬rst divided di¬erences being
The hyperbolic systems we encounter in physical problems are written
1 1
Di H ± = in what are e¬ectively the mixed variables where the apparent behavior
f (ui ) ± ±i0 +1/2 ui , (14.27)
2 2 is quite complicated. A transformation is required to decouple them back
into unmixed ¬elds that exhibit the pure contact, shock, and rarefaction
phenomena (as well as bulk convection and waves). In a real system, this
±i0 +1/2 = max (|»i0 |, |»i0 +1 |) (14.28) perfect decoupling is not possible, because the mixing is nonlinear, but it
can be achieved approximately over a small space and time region, and
is our dissipation coe¬cient. The second and third divided di¬erences, this provides the basis for the theoretical understanding of the structure of
Di+1/2 H ± and Di H ± , are then de¬ned in the standard way, like those
2 3
general hyperbolic systems of conservation laws. This is called a transfor-
of H. mation to characteristic variables. As we shall see, this transformation also
For H + , set k = i0 . Then, replacing H with H + everywhere, de¬ne provides the basis for designing appropriate numerical methods.
Q1 (x), Q2 (x), Q3 (x), and ¬nally Fi++1/2 using the ENO-Roe algorithm Consider a simple hyperbolic system of N equations
above. For H , set k = i0 + 1. Then, replacing H with H ’ everywhere,

Ut + [F (U )]x = 0 (14.31)
de¬ne Q1 (x), Q2 (x), Q3 (x), and ¬nally Fi’+1/2 again by using the ENO-
Roe algorithm above. Finally, in one spatial dimension. The basic idea of characteristic numerical schemes
is to transform this nonlinear system to a system of N (nearly) independent
Fi0 +1/2 = Fi++1/2 + Fi’+1/2 (14.29) scalar equations of the form
0 0

is the new numerical ¬‚ux function with added high-order dissipation. ut + »ux = 0 (14.32)
14.5. Systems of Conservation Laws 161 162 14. Hyperbolic Conservation Laws and Compressible Flow

is fully discretized, we multiply the entire system by L’1 = R0 to return
and discretize each scalar equation independently in an upwind biased 0
fashion using the characteristic velocity ». Then transform the discretized to the original variables.
system back into the original variables. In terms of our original equation (14.31), our procedure for discretizing
at a point x0 is simply to multiply the entire system by the left eigenvector
matrix L0 ,
14.5.1 The Eigensystem
[L0 U ]t + [L0 F (U )]x = 0, (14.38)
In a smooth region of the ¬‚ow, we can get a better understanding of the
and discretize the N scalar components of this system, indexed by p,
structure of the system by expanding out the derivative as
Ut + J Ux = 0, (14.33) [(L0 U )p ]t + [(L0 F (U ))p ]x = 0 (14.39)
where J = ‚ F /‚ U is the Jacobian matrix of the convective ¬‚ux function. independently, using upwind biased di¬erencing with the upwind direction
for the pth equation determined by the sign of »p . We then multiply the
If J were a diagonal matrix with real diagonal elements, this system would
decouple into N independent scalar equations, as desired. In general, J resulting spatially discretized system of equations by R0 to recover the
is not of this form, but we can transform this system to that form by spatially discretized ¬‚uxes for the original variables
multiplying through by a matrix that diagonalizes J. If the system is indeed
Ut + R0 ∆(L0 F (U )) = 0, (14.40)
hyperbolic, J will have N real eigenvalues »p , p = 1, . . . , N , and N linearly
independent right eigenvectors. If we use these as columns of a matrix R, where ∆ stands for the upwind biased discretization operator, i.e., either
this is expressed by the matrix equation the ENO-RF or ENO-LLF discretization.
We call »p the pth characteristic velocity or speed, (L0 U )p = Lp · U
JR = RΛ, (14.34) 0
the pth characteristic state or ¬eld (here Lp denotes the pth row of L,
where Λ is a diagonal matrix with the elements »p , p = 1, . . . , N , on the di- i.e., the pth left eigenvector of J), and (L0 F (U ))p = Lp · F (U ) the pth
agonal. Similarly, there are N linearly independent left eigenvectors. When characteristic ¬‚ux. According to the local linearization, it is approximately
these are used as the rows of a matrix L, this is expressed by the matrix true that the pth characteristic ¬eld rigidly translates in space at the pth
equation characteristic velocity. Thus this decomposition corresponds to the local
physical propagation of independent waves or signals.
LJ = ΛL, (14.35)
where L and R can be chosen to be inverses of each other, LR = RL =
14.5.2 Discretization
I. These matrices transform to a system of coordinates in which J is
At a speci¬c ¬‚ux location xi0 +1/2 midway between two grid nodes, we wish
LJR = Λ, (14.36) to ¬nd the vector numerical ¬‚ux function Fi0 +1/2 . First we evaluate the
eigensystem at the point xi0 +1/2 using the standard average Ui0 +1/2 =
as desired.
(Ui + Ui+1 )/2. Note that there are more advanced ways to evaluate the
Suppose we want to discretize our equation at the node x0 , where L and R
eigensystem, as detailed by Donat and Marquina [55]; see also Fedkiw et al.
have values L0 and R0 . To get a locally diagonalized form, we multiply our
[65]. Then, in the pth characteristic ¬eld we have an eigenvalue »p (Ui0 +1/2 ),
system equation by the constant matrix L0 that nearly diagonalizes J over
the region near x0 . We require a constant matrix so that we can move it left eigenvector Lp (Ui0 +1/2 ), and right eigenvector Rp (Ui0 +1/2 ). We put U
inside all derivatives to obtain
values and F (U ) values into the pth characteristic ¬eld by taking the dot
product with the left eigenvector,
[L0 U ]t + L0 JR0 [L0 U ]x = 0, (14.37)
u = Lp (Ui0 +1/2 ) · U
where we have inserted I = R0 L0 to put the equation in a more recognizable (14.41)
form. The spatially varying matrix L0 JR0 is exactly diagonalized at the
f (u) = Lp (Ui0 +1/2 ) · F (U ) (14.42)
point x0 , with eigenvalues »p , and it is nearly diagonalized at nearby points.
Thus the equations are su¬ciently decoupled for us to apply upwind biased where u and f (u) are scalars. Once in the characteristic ¬eld we perform a
discretizations independently to each component, with »p determining the scalar version of the conservative ENO scheme, obtaining a scalar numer-
ical ¬‚ux function Fi0 +1/2 in the scalar ¬eld. We take this ¬‚ux out of the
upwind biased direction for the pth component equation. Once this system
14.6. Compressible Flow Equations 163 164 14. Hyperbolic Conservation Laws and Compressible Flow

characteristic ¬eld by multiplying by the right eigenvector, internal energy and the kinetic energy,
E = ρe + ρ(u2 + v 2 + w2 )/2,
Fip +1/2 = Fi0 +1/2 Rp (Ui0 +1/2 ), (14.48)

where e is the internal energy per unit mass. The two-dimensional Euler
where Fip +1/2 is the portion of the numerical ¬‚ux function Fi0 +1/2 from
equations are obtained by setting w = 0, while the one-dimensional Euler
the pth ¬eld. Once we have evaluated the contribution to the numerical equations are obtained by setting both v = 0 and w = 0.
¬‚ux function from each ¬eld, we get the total numerical ¬‚ux by summing The pressure can be written as a function of density and internal energy,
the contributions from each ¬eld, p = p(ρ, e). The speed of sound is de¬ned by
Fip +1/2 ,
Fi0 +1/2 = (14.44) ppe
c= pρ + , (14.49)

completing the evaluation of our numerical ¬‚ux function at the point where pρ and pe are partial derivatives of the pressure with respect to the
xi0 +1/2 . density and internal energy, respectively.

14.6.1 Ideal Gas Equation of State
14.6 Compressible Flow Equations For an ideal gas we have p = ρRT where R = Ru /M is the speci¬c gas
constant, with Ru ≈ 8.31451 J/(mol K) the universal gas constant and M
The equations for one-phase compressible ¬‚ow are a general system of
the molecular weight of the gas. Also valid for an ideal gas is cp ’ cv = R,
convection-di¬usion-reaction conservation equations in up to three spatial
where cp is the speci¬c heat at constant pressure and cv is the speci¬c heat
dimensions. For example, in two spatial dimensions, the equations are of
at constant volume. The ratio of speci¬c heats is given by γ = cp /cv . For
the form
an ideal gas, one can write
Ut + F (U )x + G(U )y = Fd ( U )x + Gd ( U )y + S(U ), (14.45)
de = cv dT, (14.50)
where U is the vector of conserved variables, F (U ) and G(U ) are the vectors and assuming that cv does not depend on temperature (calorically perfect
of convective ¬‚uxes, Fd ( U ) and Gd ( U ) are the vectors of di¬usive ¬‚uxes, gas), integration yields
and S(U ) is the vector of reaction terms. Again, for high-speed ¬‚ow with
e = eo + cv T, (14.51)
shocks, one can usually ignore the di¬use ¬‚uxes. We choose to ignore the
where eo is not uniquely determined, and one could choose any value for e
source terms (e.g., the e¬ects of chemical reaction) here as well. For more
at 0 K (although one needs to use caution when dealing with more than
details on the di¬use terms and the source terms, see, for example, Fedkiw
one material to be sure that integration constants are consistent with the
et al. [68, 67, 69].
heat release in any chemical reactions that occur). For more details, see,
The inviscid Euler equations for one-phase compressible ¬‚ow in the
e.g., Atkins [7]. Note that
absence of chemical reactions are then
Ut + F (U )x + G(U )y + H(U )z = 0, (14.46) ρ(e ’ eo ) = (γ ’ 1)ρ(e ’ eo ),
p = ρRT = (14.52)
which can be written in detail as
and equation (14.51) are used frequently with eo = 0 arbitrarily for
«  «  «  « 
ρ ρu ρv ρw simplicity.
¬ ρu · ¬ ρu + p · ¬ · ¬ ·
ρuv ρuw
¬ · ¬ · ¬ · ¬ ·
¬ ρv · + ¬ · + ¬ ρv 2 + p · +¬ · =0
ρuv ρvw
¬ · ¬ · ¬ · ¬ ·
14.6.2 Eigensystem
 ρw       
ρw2 + p
ρuw ρvw
E (E + p)u x (E + p)v (E + p)w For brevity we consider only the two-dimensional eigensystem here. The
t y z
(14.47) two-dimensional Euler equations can be obtained by setting w = 0 so that
where ρ is the density, V = u, v, w are the velocities, E is the total energy both the fourth equation in equations (14.47) and the entire H(U )z term
per unit volume, and p is the pressure. The total energy is the sum of the vanish.
14.6. Compressible Flow Equations 165 166 14. Hyperbolic Conservation Laws and Compressible Flow

The eigenvalues and eigenvectors for the Jacobian matrix of F (U ) are den vel
obtained by setting A = 1 and B = 0 in the following formulas, while those 80
for the Jacobian of G(U ) are obtained with A = 0 and B = 1.
The eigenvalues are
»1 = u ’ c, »4 = u + c,
»2 = »3 = u,
ˆ ˆ ˆ (14.53)
and the eigenvectors are
b2 u
ˆ b1 u A b1 v B b1
L1 = + ,’ ’ ,’ ’, , (14.54) 30
2 2c 2 2c 2 2c 2
L = (1 ’ b2 , b1 u, b1 v, ’b1 ) , (14.55)
L3 = (ˆ, B, ’A, 0) ,
v (14.56)
b2 u
ˆ b1 u A b1 v B b1 1.5
’ ,’ + ,’
L= +, , (14.57)
2 2c 2 2c 2 2c 2 0.4
0 0.6 1
0.8 0 0.8 1
0.2 0.6
«  « 
1 1
¬ u ’ Ac · ¬ ·
R1 = ¬ ·, R2 = ¬ ·, (14.58)
 v ’ Bc   
temp press
x 10
H ’ uc H ’ 1/b1
«  « 
0 1 2
¬B· ¬ u + Ac ·
R3 = ¬ ·, R4 = ¬ · 325
 v + Bc  , (14.59)
 ’A  1.8
’ˆv H + ucˆ
315 1.6
2 2 2
v = Av ’ Bu,
q =u +v , u = Au + Bv,
ˆ ˆ (14.60) 1.4

“p 300
“ = pe /ρ, c= pρ + , H = (E + p) /ρ, (14.61) 1.2
ρ 295

b1 = “/c2 , b2 = 1 + b1 q 2 ’ b1 H. (14.62) 1

The choice of eigenvectors one and four is unique (up to scalar multi- 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
ples), but the choice for eigenvectors two and three is not unique. Any two
Figure 14.1. Standard shock tube test case. The solution computed with 100 grid
independent vectors from the span of eigenvectors two and three could be
cells is depicted by circles, while the exact solution is drawn as a solid line.
used instead. In fact, the numerical method designed by Fedkiw et al. [67]
exploits this fact.
The eigensystem for the one-dimensional Euler equations is obtained by be independently approximated using either an ENO-RF or ENO-LLF dis-
setting v = 0. cretization scheme. A sample calculation in one spatial dimension is shown
in Figure 14.1. The initial data consisting of two constant states form (from
left to right) a rarefaction wave, a contact discontinuity, and a shock wave.
14.6.3 Numerical Approach The solution computed with 100 grid cells is depicted by circles, while the
exact solution is drawn as a solid line. The ¬gures show solutions at a later
Since the three-dimensional Euler equations are a system of conservation
time for density, velocity, temperature, and pressure.
laws, the methods outlined earlier in this chapter can be applied in a
straight-foward fashion. That is, each of F (U )x , G(U )y , and H(U )z can
168 15. Two-Phase Compressible Flow

den vel


15 3


Two-Phase Compressible Flow 2 20


0.6 1
0 0.2 0.8
0.4 0.2
0 0.4 0.8
0.6 1

entropy press
4 5
x 10 x 10

15.1 Introduction 6
Chronologically, the ¬rst attempt to use the level set method for ¬‚ows in- 1.2
volving external physics was in the area of two-phase inviscid compressible 1.1
¬‚ow. Mulder et al. [115] appended the level set equation 3
φt + V · 2
φ=0 (15.1)
0 0.4 0.8
0.6 1 0.8 1
0 0.6
0.2 0.4
to the standard equations for one-phase compressible ¬‚ow, equation (14.47).
Figure 15.1. Spurious oscillations in pressure and velocity obtained using the
Here V is taken to be the velocity of the compressible ¬‚ow ¬eld, so that the
method proposed by Mulder et al. [115]. The solution computed with 100 grid
zero level set of φ corresponds to particle velocities and can be used to track
cells is depicted by circles, while the exact solution is drawn as a solid line.
an interface separating two di¬erent compressible ¬‚uids. The sign of φ is
used to identify which gas occupied which region, i.e., to determine the local
equation of state. In [115], only gamma law gas models were considered,
15.2 Errors at Discontinuities
with γ = γ1 for φ > 0 and γ = γ2 for φ ¤ 0. Later, Karni [93] pointed
out that this method su¬ered from spurious oscillations at the interface.
The exact solution in Figure 15.1 clearly shows that the pressure and ve-
Figure 15.1 shows a sample calculation using the method proposed in [115].
locity are continuous across the contact discontinuity (in fact, they are
Here a right going shock wave impinges upon the interface, producing both
constant in this case), while the density and entropy are discontinuous.
re¬‚ected and transmitted shock waves. Note the spurious oscillations in
Since discontinuous quantities indicate the absence of spatial derivatives
both the pressure and the velocity pro¬les near the centrally located contact
needed in equation (14.47), one should be suspicious of the behavior of nu-
discontinuity that separates the two di¬erent gamma-law gases.
15.3. Rankine-Hugoniot Jump Conditions 169 170 15. Two-Phase Compressible Flow

merical methods in that region. In fact, even the supposedly well-behaved ¬‚ux on the interface-oriented tangent to the interface so that material that
solution shown in Figure 14.1 (page 166) is not entirely adequate near the passes through this ¬‚ux passes through the interface. If the interface is
discontinuities. While the rarefaction wave (to the left) is continuous in na- moving with speed D in the normal direction, this ¬‚ux will also move with
ture, both the contact discontinuity and the shock wave should have jump speed D. From conservation, the mass, momentum, and energy that ¬‚ow
discontinuities in the computed solution. However, the computed solution into this ¬‚ux from one side of the interface must ¬‚ow back out the other side
smears out these discontinuities over a number of grid cells, leading to of the interface. Otherwise, there would be a mass, momentum, or energy
O(1) errors. Turning our attention to Figure 15.1 we note that the density sink (or source) at the interface, and conservation would be violated. This
pro¬le near the contact discontinuity (near the center) should only have tells us that the mass, momentum, and energy ¬‚ux in this moving reference
values near 1.5 on the left and values near 4.75 on the right; i.e., none of frame (moving at speed D) are continuous variables. We denote the mass,
the intermediate values should be present. Intermediate values, such as the momentum, and energy ¬‚ux in this moving reference frame by Fρ , FρV ,
one near 2.5, represent a rather signi¬cant O(1) error. In light of this, the and FE , respectively. The statement that these variables are continuous
oscillations in the pressure and velocity shown in Figure 15.1 are no worse is equivalent to the Rankine-Hugoniot jump conditions for an interface
than should be expected given the signi¬cant errors in the density pro¬le. moving with speed D in the normal direction.
The only di¬erence is that the density errors are dissipative in nature, while To de¬ne Fρ , FρV , and FE , we write the equations in conservation form
the pressure and velocity errors are dispersive in nature. for mass, momentum, and energy as in equation (14.47). The ¬‚uxes for
Of course, one could argue that dispersive errors are worse than dissi- these variables are then rewritten in the reference frame of a ¬‚ux that is
pative errors, since dispersive errors produce new extrema, changing the tangent to the interface by simply taking the dot product with the normal
monotonicity of the solution, while dissipative errors only connect two ad- direction,
missible states, producing no new extrema. While this argument is valid for «  « 
ρ 0
the shock and the contact discontinuity in Figure 14.1 and valid for both
F (U ), G(U ), H(U ) · N =  ρV T  Vn +  pN T  , (15.2)
of the shocks in Figure 15.1, it is not valid for the contact discontinuity.

<< . .

( 23)

. . >>

Copyright Design by: Sunlight webdesign