LINEBURG

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

1

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

where

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

0

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-

0

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

0

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

diagonalized,

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.

0

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-

0

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)

(14.43)

0

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

0

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)

0

ρ2

p

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

R

Ut + F (U )x + G(U )y + H(U )z = 0, (14.46) ρ(e ’ eo ) = (γ ’ 1)ρ(e ’ eo ),

p = ρRT = (14.52)

cv

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 · ¬ · ¬ ·

2

ρ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

3

for the Jacobian of G(U ) are obtained with A = 0 and B = 1.

70

The eigenvalues are

60

»1 = u ’ c, »4 = u + c,

»2 = »3 = u,

ˆ ˆ ˆ (14.53)

50

2.5

and the eigenvectors are

40

b2 u

ˆ b1 u A b1 v B b1

L1 = + ,’ ’ ,’ ’, , (14.54) 30

2 2c 2 2c 2 2c 2

2

20

2

L = (1 ’ b2 , b1 u, b1 v, ’b1 ) , (14.55)

10

L3 = (ˆ, B, ’A, 0) ,

v (14.56)

0

b2 u

ˆ b1 u A b1 v B b1 1.5

4

’ ,’ + ,’

L= +, , (14.57)

2 2c 2 2c 2 2c 2 0.4

0.2

0 0.6 1

0.8 0 0.8 1

0.2 0.6

0.4

« «

1 1

¬ u ’ Ac · ¬ ·

u

R1 = ¬ ·, R2 = ¬ ·, (14.58)

v ’ Bc

v

temp press

5

x 10

H ’ uc H ’ 1/b1

ˆ

« «

0 1 2

330

¬B· ¬ u + Ac ·

R3 = ¬ ·, R4 = ¬ · 325

v + Bc , (14.59)

’A 1.8

320

’ˆv H + ucˆ

315 1.6

where

310

2 2 2

v = Av ’ Bu,

q =u +v , u = Au + Bv,

ˆ ˆ (14.60) 1.4

305

“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

290

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

5

120

4.5

100

4

80

3.5

60

15 3

40

2.5

Two-Phase Compressible Flow 2 20

1.5

0

1

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

1.8

10

1.7

9

1.6

8

1.5

7

1.4

15.1 Introduction 6

1.3

5

Chronologically, the ¬rst attempt to use the level set method for ¬‚ows in- 1.2

4

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

1

φt + V · 2

φ=0 (15.1)

0.2

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.

Copyright Design by: Sunlight webdesign