LINEBURG

The gas to the left should never take on values above ρ = 1.5, and the

gas to the right should never take on values below ρ = 3, except at the

where Vn = V · N is the local ¬‚uid velocity normal to the interface, and the

smeared-out contact discontinuity, which can produce new extrema for the

superscript “T ” designates the transpose operator. Then the measurements

gas to the left and for the gas to the right. Since both of these gases are well-

are taken in the moving reference frame (speed D) to obtain

behaved gamma-law gases, it turns out that the oscillations in pressure and

« «

velocity can be removed without removing the numerical smearing of the ρ 0

density. However, if one of these gases is replaced with a compressible (and ¬ ρ V T ’ DN T ·

(Vn ’ D) + pN T , (15.3)

sti¬) Tait equation of state for water (that cavitates for densities less than 2

0

ρe + ρ|V ’DN | + p

around 999 kg/m3 ), it becomes rather di¬cult to remove the oscillations 2

in the pressure and velocity while still allowing the rather large dissipative

from which we can de¬ne

errors in the density pro¬le. On the other hand, we will see (below) that it

is rather easy to remove all these errors at once. Fρ = ρ(Vn ’ D), (15.4)

FρV = ρ V T ’ DN T (Vn ’ D) + pN T , (15.5)

15.3 Rankine-Hugoniot Jump Conditions

As can be seen in Figure 15.1, the pressure and velocity are continuous ρ|V ’ DN |2

+ p (Vn ’ D),

FE = ρe + (15.6)

across the contact discontinuity, while the density and entropy are discon- 2

tinuous. If we wish to remove the numerical errors caused by nonphysical

smearing of discontinuous quantities, we need to identify exactly what is as continuous variables across the interface. That is, these quantities should

and what is not continuous across the interface. be continuous across the interface in order to enforce conservation.

We de¬ne the jump in a quantity across the interface as [±] = ±1 ’ ±2 ,

In general, conservation of mass, momentum, and energy can be applied

where ±1 is the value of ± in ¬‚uid 1, and ±2 is the value of ± in ¬‚uid 2.

to an interface in order to abstract continuous variables. One can place a

15.4. Nonconservative Numerical Methods 171 172 15. Two-Phase Compressible Flow

Then we can summarize by stating that [Fρ ] = 0, [FρV ] = 0, and [FE ] = 0 see Puckett et al. [134]. These overshoots can be ignored when they violate

across an interface moving with speed D in the normal direction. conservation, or redistributed in a manner similar to [45] to preserve exact

conservation.

15.4 Nonconservative Numerical Methods 15.5 Capturing Conservation

Traditionally, Eulerian-based numerical methods for compressible ¬‚ow are In summary, conservation of mass, momentum and energy at a disconti-

based on the Lax-Wendro¬ theorem [104], which dictates that numerical nuity tells us which variables are continuous across the interface, although

methods should be fully conservative, and it is well known that noncon- as pointed out by Karni [92] one does not necessarily need exact ¬‚ux-

servative methods produce shocks with incorrect speeds and strengths. di¬erenced conservative form in order to obtain the correct weak solution.

However, Karni [92] advocated nonconservative form at lower-dimensional That is, one can instead obtain the correct weak solution in a di¬erent man-

(e.g., one-dimensional in a two-dimensional calculation) material interfaces ner, for example by solving an associated Riemann problem as in [129].

(contact discontinuities) in order to alleviate the oscillations observed in Therefore, Fedkiw et al. [63] proposed implicitly capturing the Rankine-

[115]. In [92], full conservation was applied away from interfaces, and a Hugoniot jump conditions at the interface in order to avoid the intricate

nonconservative method was applied near the interface without adversely details (e.g., multidimensional solution of Riemann problems) that are in-

a¬ecting the shock speeds or strengths. Since shocks do not move at the volved in explicitly enforcing the conservation at the interface. This idea

local interface velocity, any portion of a shock is in contact only with an of implicitly capturing conservation as opposed to explicitly enforcing it is

interface, and thus the nonconservative discretization employed there, on in the ¬‚avor of level set methods that implicitly capture the location of an

a set of measure zero in space and time, minimizing the accumulation of interface as opposed to explicitly tracking it. Similar to level set methods,

error. this implicit capturing method for enforcing conservation gives rise to a

While it is true that others have used nonconservative discretizations, simple-to-implement numerical algorithm.

Karni [92] is responsible for markedly increasing their popularity in the This implicit approach to capturing conservation is applied by creating a

shock-capturing community, where practitioners usually required conser- set of ¬ctitious ghost cells on each side of the interface corresponding to the

vation at all cost. It is interesting to note that many front-tracking and real ¬‚uid on the other side. These ghost cells are populated with a specially

volume-of-¬‚uid schemes are actually nonconservative; i.e., they do not sat- chosen (ghost) ¬‚uid that implicitly captures the Rankine-Hugoniot jump

isfy the strict ¬‚ux-di¬erencing conservation form usually thought to be conditions across the interface. In [63], this method was referred to as the

required by the Lax-Wendro¬ theorem. In this sense, many of these schemes ghost ¬‚uid method (GFM).

share similar properties with the ideology proposed in [92]. For example,

consider the front-tracking approach of Pember et al. [129], where a high-

order Godunov method is used to obtain a nonconservative update near the

15.6 A Degree of Freedom

tracked interface and a fully conservative update away from the tracked in-

terface. All ¬‚ow features including shock speeds and strengths as well as

A contact discontinuity (or material interface) has speed D = Vn . Rewriting

the speed of the tracked front are correctly determined, as is ensured by the

[Fρ ] = 0 as ρ1 (Vn ’ D) = ρ2 (Vn ’ D) using equation (15.4) and choosing

1 2

solutions of the appropriate Riemann problems. Then the authors go one 1 2 1 2

either D = Vn or D = Vn leads directly to Vn = Vn or [Vn ] = 0. This

step further and correct the lack of conservation at the interface using a re-

is our ¬rst jump condition, [Vn ] = 0, implying that the normal velocity is

distribution procedure due to Chern and Colella [45] that is (presumably) 1 2

continuous across the interface. Setting D = Vn = Vn in [FρV ] = 0 leads to

not necessary for obtaining a grid-resolved solution, but is used only to

[p]N T = 0, using equation (15.5) and the fact that the normal is the same

maintain exact conservation. In fact, the nature of this redistribution pro-

on both sides of the interface, i.e., [N ] = 0. Multiplying [p]N T = 0 by N

cedure does not allow strict application of the Lax-Wendro¬ theorem, and

leads to our second jump condition, [p] = 0, implying that the pressure is

one must assume that the correct solutions are obtained because the numer-

continuous across the interface. Note that multiplying [p]N T = 0 by any

ical method is fully conservative except at the lower-dimensional tracked

interface, which is updated correctly based on solutions of Riemann prob- tangent vector (there are two of these in three spatial dimensions) leads to

lems. Similar loss of exact conservation occurs in volume-of-¬‚uid methods, 0 = 0 and the failure to produce a jump condition or a continuous variable.

1 2

where nonphysical overshoots may occur in the volume fraction equation; Plugging D = Vn = Vn into [FE ] = 0 also leads to 0 = 0 and again a failure

15.7. Isobaric Fix 173 174 15. Two-Phase Compressible Flow

to produce a jump condition or a continuous variable. Thus, at a contact arti¬cial heat conduction to the numerical method in a form similar to ar-

discontinuity in three spatial dimensions we have two jump conditions, ti¬cial viscosity. Later, Donat and Marquina [55] proposed a ¬‚ux-splitting

[Vn ] = 0 and [p] = 0, along with three degrees of freedom corresponding to method with a built-in heat conduction mechanism that dissipates these

the 0 = 0 trivially satis¬ed jump conditions. errors throughout the ¬‚uid.

There are ¬ve eigenvalues present in the equations of compressible ¬‚ow For the one-dimensional Euler equations, the Rankine-Hugoniot jump

in three spatial dimensions. Two of these correspond to the genuinely conditions for the solid wall moving at the local ¬‚ow velocity D = VN are

nonlinear ¬eld associated with sound waves, and the other three of these [VN ] = 0, [p] = 0, and 0 = 0. These describe the relationship between the

correspond to the linearly degenerate particle velocity. Since a contact dis- external ¬‚ow ¬eld and the internal one; i.e., both the normal velocity and

continuity moves with the linearly degenerate particle velocity, nothing the pressure must be continuous across the solid wall boundary extending

represented by the linearly degenerate ¬elds can cross the interface, mean- into the ghost cells. Since these jump conditions are inherently part of

ing that these quantities, the two-dimensional tangential velocity and the the equations and thus part of any consistent numerical method, jumps

entropy, do not cross the interface and may be discontinuous there. In fact, in pressure and velocity are hard to maintain for any duration of time at

since these quantities are uncoupled across the interface, they are usually a solid wall boundary; i.e., jumps between the ¬‚uid values and the ghost

discontinuous there. We can write cell values are quickly dissipated. In this sense, one can think of pressure

and velocity equilibration at a solid wall boundary as an intrinsic action of

St + V · S = 0, (15.7)

the boundary conditions. There is no such condition for the temperature

where S is the entropy, to indicate that entropy is advected along stream- or the density. In the case of a complete equation of state (see Davis [54])

lines. Then since S and the contact discontinuity move at the same speed only one variable in the linearly degenerate ¬eld need be de¬ned, and all

in the normal direction, information associated with S cannot cross the other variables can be determined from the equation of state relations. In

interface. (As a disclaimer we note that equation (15.7) is not valid for this sense, there is no boundary condition for the linearly degenerate ¬eld,

streamlines that cross shock waves, i.e., entropy jumps across a shock wave. as is emphasized by the trivially satis¬ed jump condition 0 = 0. Since a

However, shock waves do not move at the linearly degenerate speed, so this solid wall boundary is an initial boundary value problem, the value of the

equation is true except for a lower-dimensional subset of space and time.) temperature at the wall must come from the initial data, as one can see from

Note that in the case of the full viscous Navier-Stokes equations, equation (15.7) which states that entropy is advected along streamlines of

the physical viscosity imposes continuity of the tangential velocities and the ¬‚uid. This implies that the entropy near the wall stays near the wall,

thermal conductivity imposes continuity of the temperature. since the wall moves with the local ¬‚uid velocity.

In [66], equation (15.7) was used to develop the isobaric ¬x, which is a

boundary condition type of treatment for the linearly degenerate ¬eld at a

15.7 Isobaric Fix solid wall boundary. The isobaric ¬x modi¬es the linearly degenerate ¬eld

at a solid wall without changing the values of the pressure or the normal

In [66], Fedkiw et al. exploited this degree of freedom in order to signi¬- velocity. Noting that entropy is advected along streamlines and that the

cantly reduce errors caused by nonphysical wall heating. The well-known entropy within a ¬‚uid is usually continuous, we see that the entropy errors

“overheating e¬ect” occurs when a shock re¬‚ects o¬ of a solid wall bound- at the wall are repaired using new values of entropy extrapolated from

ary, causing overshoots in the temperature and density, while the pressure the surrounding ¬‚ow. For example, replacing the entropy at the wall with

and velocity remain constant. In one spatial dimension, a solid wall bound- the entropy of the neighboring cell gives a ¬rst-order accurate value of

ary condition can be applied with the aid of ghost cells by constructing a the entropy at the wall for smooth entropy pro¬les. Higher-order accurate

symmetric pressure and density re¬‚ection and an asymmetric normal ve- extrapolation can be used as well, but this has been found to be quite

locity re¬‚ection about the solid wall. Then a shock wave impinging on the dangerous in practice due to the presence of discontinuous shock waves

wall will collide with a shock in the ghost cells that has equal strength that can cause large overshoots when one extrapolates with higher than

traveling in the opposite direction, producing the desired shock re¬‚ection. ¬rst-order accuracy. In multiple spatial dimensions the solid wall can be

Meniko¬ [113] and Noh [122] showed that overheating errors are a symptom represented as the zero isocontour of a level set function and moved rigidly

of smeared-out shock pro¬les and that sharper shocks usually produce less using the level set equation (3.2) where the velocity is chosen as the rigid-

overheating. They also showed that the pressure and velocity equilibrate wall velocity (or even the velocity of a deforming wall). Then the isobaric

quickly, while errors in the temperature and density persist. In order to ¬x can be applied using extrapolation of entropy in the normal direction,

dissipate these errors in temperature and density, [122] proposed adding as discussed in Chapter 8. Note that one does not have to deal with the

15.8. Ghost Fluid Method 175 176 15. Two-Phase Compressible Flow

entropy directly, which can sometimes be di¬cult to compute for general den vel

equations of state, but one can choose any variable corresponding to the

45

2000

equation of state degree of freedom in the linearly degenerate ¬eld, e.g.,

density or temperature. 40

Although the overheating e¬ect is traditionally discussed in the context

1500

35

of shock re¬‚ection from stationary walls, more signi¬cant cumulative over-

heating e¬ects are generated by moving walls. Figure 15.2 shows the errors 30

generated in density and temperature by a solid wall moving from left to 1000

25

right. The wall, initially at rest, is instantaneously accelerated to a velocity

of 1000 m/s forming the shock wave seen in the ¬gure. The isobaric ¬x does 20

500

not completely eliminate all of the density and temperature errors, but does

15

reduce them by at least an order of magnitude, as shown in Figure 15.3.

10 0

0.4 0.45

0.55 0.5

0.45 0.5 0.55

0.4

15.8 Ghost Fluid Method

In [63], Fedkiw et al. pointed out that a two-phase contact discontinuity

press temp

6

could be discretized with techniques similar to those used for the solid- x 10

1400

15

wall boundary, except that they are applied twice, i.e., once for each ¬‚uid.

Conceptually, each grid point corresponds to one ¬‚uid or the other, and

1200

ghost cells can be de¬ned at every point in the computational domain so

that each grid point contains the mass, momentum, and energy for the

10 1000

real ¬‚uid that exists at that point (according to the sign of the level set

function) and a ghost mass, momentum, and energy for the other ¬‚uid that

does not really exist at that grid point (the ¬‚uid from the other side of the 800

interface). Once the ghost cells are de¬ned, standard one-phase numerical

5

methods can be used on the entire domain for each ¬‚uid; i.e., we now have 600

two separate single-¬‚uid problems. After each ¬‚uid is advanced in time,

the level set function is updated using equation (3.2) to advect the level 400

set with the local ¬‚uid velocity V , and the sign of the level set function

0

is used to determine the appropriate real ¬‚uid values at each grid point. 200

0.4 0.45 0.5 0.55 0.4 0.45 0.5 0.55

While ghost cells are de¬ned everywhere for the sake of exposition, only a

band of 3 to 5 ghost cells is actually needed in practice. Figure 15.2. Overheating errors in density and temperature generated by a piston

Contact discontinuities move at the local ¬‚uid velocity, and the Rankine- moving to the right.

Hugoniot jump conditions are [VN ] = 0, [p] = 0, and 0 = 0 three times. Only

the pressure and normal velocities can be determined from the boundary

conditions, while the entropy and both tangential velocities remain unde- locity, the ghost ¬‚uid values can be set equal to the real ¬‚uid values at each

termined. Since certain properties are discontinuous across the interface, grid point, implicitly capturing the correct interface values of these vari-

one should be careful in applying ¬nite di¬erence methods across the inter- ables. This is the key mechanism in coupling the two distinct sets of Euler

face, since di¬erencing discontinuous quantities leads erroneously to terms equations. On the other hand, the discontinuous variables move with the

of the form 1/ x that increase without bound as the grid is re¬ned. There- speed of the interface, and information in these variables does not cross the

fore, the layer of ghost cells should be introduced so that there is continuity interface and is not coupled to the corresponding information on the other

with the neighboring ¬‚uid that needs to be discretized. For variables that side of the interface. Moreover, in order to avoid numerical smearing or

are already continuous across the interface, e.g., pressure and normal ve- spurious oscillations, these discontinuous variables should not be nonphys-

15.8. Ghost Fluid Method 177 178 15. Two-Phase Compressible Flow

In order to to extrapolate the tangential velocity, we ¬rst extrapolate

den vel

the entire velocity ¬eld V . Then, at every cell in the ghost region we have

45

2000 two separate velocity ¬elds, one from the real ¬‚uid and one from the ex-

trapolated ¬‚uid. For each velocity ¬eld, the normal component of velocity,

40

VN , is put into a vector of length three, VN N , and then the tangential

1500

35

velocity ¬eld is de¬ned by another vector of length three, V ’ VN N . We

add together the normal component of velocity, VN N , from the real ¬‚uid

30

and the tangential component of velocity, V ’ VN N , from the extrapo-

1000

25 lated ¬‚uid. This new velocity is our ghost ¬‚uid velocity. For the full viscous

Navier-Stokes equations, the physical viscosity imposes continuity of the

20

500 tangential velocities, so that the entire velocity ¬eld is continuous. In this

15 case the entire velocity ¬eld can be copied over into the ghost cells in a

node-by-node fashion. Similar statements hold for the temperature in the

10 0

presence of a nonzero thermal conductivity.

Next, we describe the one-dimensional method in more detail. Suppose

0.4 0.45

0.55 0.5

0.4 0.45 0.5 0.55

that the interface lies between nodes i and i + 1. Then ¬‚uid 1 is de¬ned to

the left and includes node i, while ¬‚uid 2 is de¬ned to the right and includes

node i+1. In order to update ¬‚uid 1, we need to de¬ne ghost ¬‚uid values of

press temp

6

x 10 ¬‚uid 1 at nodes to the right, including node i+1. For each of these nodes we

15

de¬ne the ghost ¬‚uid values by combining ¬‚uid 2™s pressure and velocity at

1100

each node with the entropy of ¬‚uid 1 from node i; see Figure 15.4. Likewise,

1000 we create a ghost ¬‚uid for ¬‚uid 2 in the region to the left, including node i.

This is done by combining ¬‚uid 1™s pressure and velocity at each node with

900

10

the entropy of ¬‚uid 2 from node i + 1. In order to apply the isobaric ¬x

800

technique, we change the entropy at node i to be equal to the entropy at

700 node i ’ 1 without modifying the values of the pressure and velocity at

node i; see Figure 15.5. Likewise, we change the entropy at node i + 1 to

600

5

be equal to the entropy at node i + 2.

500

An important aspect of this method is its simplicity. We do not need to

400 solve a Riemann problem, consider the Rankine-Hugoniot jump conditions,

or solve an initial boundary value problem at the interface. We capture the

300

0 appropriate interface conditions by de¬ning a ¬‚uid that has the pressure

0.4 0.45 0.5 0.55 0.4 0.45 0.5 0.55

and velocity of the real ¬‚uid at each point, but the entropy of some other

¬‚uid. Consider the case of air and water. In order to solve for the air,

Figure 15.3. The isobaric ¬x signi¬cantly reduces the overheating errors in both

we replace the water with ghost air that acts like the water in every way

density and temperature.

(pressure and velocity) but appears to be air (entropy). In order to solve

for the water, we replace the air with ghost water that acts like the air in

every way (pressure and velocity) but appears to be water (entropy). Since

ically coupled together or forced to be continuous across the interface. The

the ghost ¬‚uids behave in a fashion consistent with the real ¬‚uids that they

most obvious way of de¬ning the discontinuous variables in the ghost cells

are replacing, the appropriate boundary conditions are captured. Since the

is by extrapolating that information from the neighboring real ¬‚uid nodes;

ghost ¬‚uids have the same entropy as the real ¬‚uid that is not replaced, we

e.g., the entropy can be extrapolated into the ghost cells using extrapo-

are solving a one-phase problem.

latation in the normal equation in the same fashion as it was in applying

Figure 15.6 shows how the ghost ¬‚uid method removes the spurious os-

the isobaric ¬x. Again, as with to the isobaric ¬x, one does not have to

ciallations in the pressure and velocity obtained using the method proposed

deal with the entropy directly, but can choose any variable in the linearly

in [115] as shown in Figure 15.1. Note that the density pro¬le remains sharp

degenerate ¬eld, e.g., density or temperature.

15.8. Ghost Fluid Method 179 180 15. Two-Phase Compressible Flow

vel

den

Interface

5 120

Fluid 2

4.5

100

i +2

i+1 i+3

P=pressure

4

V=velocity 80

S =entropy 3.5

60

3

PV PV PV 40

2.5

S 2 20

1.5

0

1

i ’2 i ’1 i 0.6 1

0 0.2 0.8

0.4

1

0 0.2 0.8

0.6

0.4

Fluid 1 Ghost Cells

entropy press

4 5

x 10 x 10

1.7

Figure 15.4. Fluid 1™s ghost ¬‚uid values are constructed by combining ¬‚uid 2™s

10

pressure and velocity with the entropy of ¬‚uid 1 from node i.

1.6

9

1.5

8

Interface

1.4

7

Fluid 2

6 1.3

i +1 i +2 i+3

P=pressure 5

1.2

V=velocity

4

S=entropy 1.1

3

1

PV PV PV 2

0.2

0 0.6

0.4 0.8 1 0.2

0 0.4 0.8

0.6 1

Figure 15.6. The spurious oscillations are removed using the ghost ¬‚uid method.

S Moreover, the density pro¬e remains sharp at the contact discontinuity. The

solution computed with 100 grid cells is depicted by circles, while the exact

solution is drawn as a solid line.

i ’2 i ’1 i

across the contact discontinuity. While Figure 15.6 is computed using only

Fluid 1 Ghost Cells

100 grid cells, Figure 15.7 is computed with 400 grid cells to illustrate the

behavior of the method under mesh re¬nement. A two-dimensional exam-

ple of an air shock hitting a helium bubble is shown in Figure 15.8. The

Figure 15.5. In applying the isobaric ¬x in conjunction with the ghost ¬‚uid

black circle indicates the initial location of the Helium bubble before it was

method, the entropy from node i ’ 1 is used instead of the entropy from node i.

hit by the air shock. Figures 15.9, 15.10, 15.11, and 15.12 show two phase-

¬‚ow calculations where one phase is a gamma-law gas model of air with a

15.8. Ghost Fluid Method 181 182 15. Two-Phase Compressible Flow

den vel

5 120

4.5

100

4

80

3.5

60

3

40

2.5

2 20

1.5

0

1

0.2

0 0.4 0.8

0.6 1 1

0.8

0.6

Copyright Design by: Sunlight webdesign