LINEBURG

0.2

0

entropy press

4 5

x 10 x 10

1.7

10

1.6

9

1.5

8

1.4

7

6 1.3

5

1.2

4

1.1

3

1

2

0.6 1

0 0.2 0.8

0.4 0 1

0.4 0.8

0.2 0.6

Figure 15.7. This ¬gure illustrates how the ghost ¬‚uid method converges as the

grid is re¬ned. The solution computed with 400 grid cells is depicted by circles,

while the exact solution is drawn as a solid line.

density around 1 kg/m3 and the other phase is a sti¬ Tait equation of state

Figure 15.8. Schlieren image for an air shock impinging upon a helium bubble

model for water with a density around 1000 kg/m3 . In the ¬gures the air

using the ghost ¬‚uid method to resolve the contact discontinuity. The black circle

is depicted in red and the water is depicted in green. Note that there is no

indicates the initial location of the helium bubble. This image was generated by

numerical smearing of the density at the interface itself which is fortunate,

Tariq Aslam, of Los Alamos National Laboratory.

since water cavitates when it drops to a density slightly above 999 kg/m3 ,

leading to host of nonphysical problems near the interface. Note that the

pressure and velocity are continuous across the interface, although there

are kinks in both of these quantities.

15.9. A Robust Alternative Interpolation 183 184 15. Two-Phase Compressible Flow

15.9 A Robust Alternative Interpolation log(den) vel

3

The interface values of pressure and normal velocity need to be determined 0

using some sort of interpolation technique, where we note that these vari- 2.5

ables are continuous but may possess kinks due to di¬ering equations of

’500

2

state across the interface. Copying these variables into the ghost cells in a

node-by-node fashion, as proposed above (and in [63]) corresponds to one 1.5

choice of interpolation. Using the ¬‚uid on one side of the interface to deter- ’1000

mine the interface pressure and the ¬‚uid on the other side of the interface 1

to determine the interface normal velocity corresponds to another choice. 0.5

Di¬erent interpolation techniques lead to O( x) di¬erences in the interface ’1500

values of pressure and normal velocity, which vanish as the mesh is re¬ned, 0

guaranteeing convergence as the Rankine-Hugoniot jump conditions are

’0.5 ’2000

implicitly captured.

2 10

4 6 0

0 8 2 8

6 10

4

It is not clear exactly which interpolation technique should be used, and

the answer is most likely problem related. For smooth well-behaved prob-

lems with commensurate equations of state, the method proposed above

(and in [63]) is probably superior, while using one ¬‚uid to de¬ne the pres-

entropy

5 press

7

x 10

sure and the other ¬‚uid to de¬ne the normal velocity is probably superior x 10

3

when one ¬‚uid is very sti¬ compared to the other. For example, consider

10

interactions between a sti¬ Tait equation-of-state for water and a gamma-

2.5

law gas model for air as shown, for example, in Figures 15.9, 15.10, 15.11,

8

and 15.12. Since the technique discussed in the last section gives equal

2

weighting to the values of the pressure and normal velocity on both sides

of the interface, any kinks in these values will be smeared out to some ex- 6

1.5

tent, causing small errors in the captured interface values of these variables.

Small errors in the normal velocity of the water create small density errors 4

1

when the equation for conservation of mass in updated. In turn, these small

density errors can lead to large spurious pressure oscillations in the water, 2 0.5

since the Tait equation of state is quite sti¬. While small errors in the ve-

locity of the air cause the same small density errors, these have little e¬ect 0 0

on the gas, since the gamma-law gas equation of state is rather robust.

0 2 8

4 6 10 0 2 6

4 8 10

Again, since the Tait equation of state is rather sti¬, one can expect large

variations in the pressure of the water near the interface, which in turn Figure 15.9. The gamma-law gas is depicted in red, while the sti¬ Tait equation

can lead to poor predictions of the interface pressure. While these errors of state water is depicted in green. Note that the log of the density is shown,

in the interface pressure have a relatively small e¬ect on the heavier water, since the density ratio is approximately 1000 to 1. This calculation uses only 100

they can have a rather large e¬ect on the lighter gas. Conversely, since the grid cells. (See also color ¬gure, Plate 6.)

gamma-law gas equation of state is rather robust, the gas pressure tends

to be smooth near the interface and is therefore a good candidate for the

interface pressure. the total velocity and the entropy are extrapolated into the ghost cells. In

The aforementioned di¬culties can be removed in large part by using the updating the ¬‚uid with the more robust equation of state (in this case the

water to determine the interface normal velocity and the air to determine gamma-law gas air), the normal velocity is still copied over node by node

the interface pressure, producing a more robust version of the interpolation. in the ghost region, while the pressure, entropy, and tangential velocity are

When the sti¬er ¬‚uid (in this case the Tait equation of state water) is extrapolated into the ghost cells. This new robust interpolation technique

updated, pressure is still copied over node by node in the ghost region, while was ¬rst proposed by Fedkiw in [62]. Numerical results have shown that

15.9. A Robust Alternative Interpolation 185 186 15. Two-Phase Compressible Flow

and x = 8. The robustnesss of the new interpolation technique is illustrated

log(den) vel

by the high quality of the solution obtained with as few as 100 grid cells,

3

as shown in Figure 15.9.

0

In Figures 15.11 and 15.12 interfaces separate gas on the outside of the

2.5

domain from water on the inside of the domain. A solid-wall boundary is

’500

2 enforced on the left, and an out¬‚ow boundary condition is enforced on the

right. Initially, all the ¬‚uids are moving to the right at 500m/s causing

1.5

’1000

a rarefaction wave to start at the solid wall on the left. This rarefaction

wave propagates to the right, slowing down the ¬‚uids. Note that it is much

1

easier to slow down the lighter gas as opposed to the heavier water. The

’1500

0.5

¬gures show the steep pressure pro¬le that forms in the water and acts to

to slow the water down. One of the di¬culties encountered in [63] was a

0

’2000

nonphysical pressure overshoot in the water near the interface on the left.

’0.5

This new robust interpolation technique removes the overshoot, producing

2 10

4 6 0

0 8 2 8

6 10

4 a monotone pressure pro¬le near the interface even in the coarse 100-grid-

cell solution in Figure 15.11.

entropy

5 press

7

x 10 x 10

12

6

10 5

8 4

6 3

4 2

2 1

0 0

0 2 8

4 6 10 4 6 8

2 10

0

Figure 15.10. This is the same calculation as in Figure 15.9, except that 500 grid

cells are used. (See also color ¬gure, Plate 7.)

this new method behaves in a fashion similar to the original method, except

for the increased interface dissipation, which leads to greater stability.

In Figures 15.9 and 15.10 an interface separates gas on the left from

water on the right. Solid-wall boundary conditions are enforced on both

sides of the domain. Initially, a right-going shock wave is located in the

gas, and a left-going shock wave is located in the water. These shock waves

propagate toward the interface, producing a complex wave interaction. In

the ¬gures one can see re¬‚ected shock waves traveling outward near x = 1

15.9. A Robust Alternative Interpolation 187 188 15. Two-Phase Compressible Flow

vel

log(den)

vel

log(den)

500

3

500

3

400

2.5 400

2.5

300

300

2

2

200

200

1.5

1.5

100

100

1

1

0

0

4 6 8

0 10

2 8 10

6

2

0 4

4 6 8

0 10

2 8 10

6

2

0 4

entropy press

4 6

x 10 x 10

entropy press

4 6

x 10 x 10

10

10

10

10

9

8 8

8 8

7

6

6 6

6

5

4

4

4

4

3

2

2

2

2

1

0

0

0 4 8

2 6 0

10 2 6 8 10

4 0

0 4 8

2 6 0

10 2 6 8 10

4

Figure 15.11. In this calculation two interfaces are present, since the air surrounds

Figure 15.12. This is the same calculation as in Figure 15.11, except that 500

the water on both sides. This calculation uses only 100 grid cells. (See also color

grid cells are used. (See also color ¬gure, Plate 9.)

¬gure, Plate 8.)

190 16. Shocks, Detonations, and De¬‚agrations

16.2 Computing the Velocity of the Discontinuity

For a contact discontinuity, we use equation (15.1) to update the interface

location, but a more general discontinuity moving at speed D in the normal

direction is governed by

16 φt + W · φ = 0, (16.1)

Shocks, Detonations, and De¬‚agrations where W = DN is the velocity of the discontinuity. The authors of [64]

proposed a capturing method to compute W without explicitly computing

it on the interface. Suppose U (1) and U (2) represent states on di¬erent sides

of the interface. Then, in general, the velocity of the interface is de¬ned

(1) (2)

by W = W (Uint , Uint ) where the “int” subscript designates a variable

that has been interpolated to the interface in a one-sided fashion. Gen-

erally, W is a continuous function of its arguments, and application of

W = W (U (1) , U (2) ) in a node-by-node fashion will implicitly capture the

correct value of W at the interface. This computation of W = W (U (1) , U (2) )

in a node-by-node fashion requires values of U (1) and U (2) at every node.

We accomplish this by extrapolating U (1) across the interface into the re-

gion occupied by U (2) (using equation (8.1)), and likewise extrapolating

U (2) across the interface into the region occupied by U (1) . Of course, we

only need to extrapolate values in a thin band of grid cells near the in-

(1) (2) (1) (2)

terface. When W = DN , i.e., W = W (Uint , Uint ) = D(Uint , Uint )N , we

16.1 Introduction

can compute D in a node-by-node fashion and obtain W by multiplying D

by N .

For a contact discontinuity we separated the variables into two sets based on

When using the ghost ¬‚uid method for general discontinuities, we need

their continuity across the interface. The continuous variables were copied

to accurately determine the interface speed D. For shock waves and deto-

into the ghost ¬‚uid in a node-by-node fashion, capturing the correct in-

nation waves, D can be found by solving an appropriate Riemann problem

terface values, while the discontinuous variables were extrapolated in a

in a node-by-node fashion [64]. In fact, there is no reason one cannot solve

one-sided fashion to avoid numerical dissipation errors. In order to apply

a Riemann problem in the case of a contact discontinuity as well, using

this idea to a general interface moving at speed D in the normal direction,

W = DN in equation (16.1) to update the level set function as opposed

we need to correctly determine the continuous and discontinuous variables

to using equation (15.1) with the less-accurate local ¬‚uid velocity. In fact,

across the interface. For example, consider a shock wave where all variables

a combination of ghost cells and Riemann problems is commonly used in

are discontinuous, and extrapolation of all variables for both the preshock

front tracking algorithms; see, e.g., Glimm et al. [73, 72], where a Rie-

and postshock ¬‚uids obviously gives the wrong answer, since the physical

mann problem is solved at the interface and the results are extrapolated

coupling is ignored. We generally state, For each degree of freedom that is

into ghost cells. The di¬erence between the ghost ¬‚uid method and typi-

coupled across a discontinuity, one can de¬ne a variable that is continu-

cal front-tracking algorithms is in the order of operations. Front-tracking

ous across the discontinuity, and all remaining degrees of freedom can be

algorithms ¬rst solve a Riemann problem using ¬‚ow variables interpolated

expressed as discontinuous variables that can be extrapolated across the dis-

to the (possibly) multidimensional interface and then extrapolate the re-

continuity in a one-sided fashion, as the key to extending the GFM. For

sults into ghost cells, while the ghost ¬‚uid method ¬rst extrapolates into

the Euler equations, conservation of mass, momentum, and energy can be

ghost cells and then solves the Riemann problem in a node-by-node fashion,

applied to any discontinuity in order to abstract continuous variables; i.e.,

removing complications due to interface geometry.

the Rankine-Hugoniot jump conditions always dictate the coupling between

For a de¬‚agration wave, the Riemann problem is not well posed unless

the prediscontinuity and postdiscontinuity ¬‚uids. This idea was proposed

the speed of the de¬‚agration (D) is given. However, the G-equation for

by Fedkiw et al. [64] to create sharp pro¬les for shock waves, detonation

¬‚ame discontinuities, see Markstein [110], represents a ¬‚ame front as a dis-

waves, and de¬‚agration waves.

16.3. Limitations of the Level Set Representation 191 192 16. Shocks, Detonations, and De¬‚agrations

continuity in the same fashion as the level set method, so that one can easily and to rewrite equation (16.4) as

consult the abundant literature on the G-equation to obtain de¬‚agration

ρ(Vn ’ D)2

wave speeds. + p (Vn ’ D)

FE = ρe + (16.6)

2

using the fact that N = ±1 in one spatial dimension.

16.3 Limitations of the Level Set Representation Our goal is to implicitly capture the Rankine-Hugoniot jump conditions

by implicitly enforcing continuity of Fρ , FρVn , and FE . These quantities

can be evaluated at each real grid node by plugging the local values of the

The level set function is designed to represent interfaces where the interface

R

conserved variables into equations (16.2), (16.5) and (16.6) to obtain Fρ ,

crosses material at most once due to an entropy condition; see Sethian [147]

R R

FρVn , and FE , respectively, where the “R” superscript designates a real

and Osher and Seian [126]. Contact discontinuities move with the local

grid node value. In order to implicitly capture the values of these variables

material velocity, and thus never cross over material. On the other hand, if

with ghost node values, we want the ghost node values to result exactly in

Copyright Design by: Sunlight webdesign