<< . .

( 23)

. . >>

E+p 0
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
ρ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

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
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
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
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
not completely eliminate all of the density and temperature errors, but does
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
15.8 Ghost Fluid Method
In [63], Fedkiw et al. pointed out that a two-phase contact discontinuity
press temp
could be discretized with techniques similar to those used for the solid- x 10
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
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
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
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
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,
VN , is put into a vector of length three, VN N , and then the tangential
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
and the tangential component of velocity, V ’ VN N , from the extrapo-
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
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
x 10 ¬‚uid 1 at nodes to the right, including node i+1. For each of these nodes we
de¬ne the ghost ¬‚uid values by combining ¬‚uid 2™s pressure and velocity at
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
the entropy of ¬‚uid 2 from node i + 1. In order to apply the isobaric ¬x
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
be equal to the entropy at node i + 2.
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
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

5 120
Fluid 2
i +2
i+1 i+3
V=velocity 80
S =entropy 3.5


S 2 20

i ’2 i ’1 i 0.6 1
0 0.2 0.8
0 0.2 0.8

Fluid 1 Ghost Cells

entropy press
4 5
x 10 x 10
Figure 15.4. Fluid 1™s ghost ¬‚uid values are constructed by combining ¬‚uid 2™s
pressure and velocity with the entropy of ¬‚uid 1 from node i.
Fluid 2
6 1.3
i +1 i +2 i+3
P=pressure 5
S=entropy 1.1
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



2 20

0 0.4 0.8
0.6 1 1

<< . .

( 23)

. . >>

Copyright Design by: Sunlight webdesign