<< . .

( 23)

. . >>

ˆ (18.25)
is used to compute the second derivative in the x-direction.
where p does not depend on t or ρ.
ˆ Discretization of equation (18.21) at each cell center leads to a coupled
Boundary conditions can be applied to either the velocity or the pressure. linear system of equations for the unknown pressures (one equation at each
In order to apply boundary conditions to V n+1 , we apply them to V after cell center). This system is symmetric and can be solved in a straightfor-
computing V in equation (18.19) and before solving equation (18.21). ward manner using a preconditioned conjugate gradient (PCG) method
Then in equation (18.21), we set p · N = 0 on the boundary, where N with an incomplete Choleski preconditioner. The interested reader is re-
is the local unit normal to the boundary. Since the ¬‚ow is incompressible, ferred to the comprehensive computational linear algebra text of Golub
the compatibility condition and Van Loan [75]. After the iterative PCG method is used to compute
the pressure at each cell center, the derivatives can be computed at the cell
V ·N =0 (18.26) faces and used to update the velocity ¬eld in equation (18.20).

214 18. Incompressible Flow 18.5. Simulating Smoke for Computer Graphics 215

18.5 Simulating Smoke for Computer Graphics
For computer graphics applications one is less concerned with overall accu-
racy and more concerned with computational e¬ciency, as long as visually
believable results can be obtained. For this reason, computer graphics sim-
ulations are usually undertaken on relatively coarse grids using methods
that allow for large time steps. A popular numerical method in the at-
mospheric sciences community is the semi-Lagrangian method; see, e.g.,
Staniforth and Cote [155] for a review. Semi-Lagrangian methods consist
in backward tracing of characteristic information from grid points to arbi-
trary points in space with subsequent interpolation of data from the grid
points to the backward-traced origins of the characteristics. This method
was ¬rst proposed by Courant et al. [51] and is guaranteed to be stable for
any time step as long as monotone interpolation is used.
In [70], Fedkiw et al. used a ¬rst-order accurate semi-Lagrangian method
to produce visually convincing simulations of smoke. Since this method is
highly dissipative, the viscous terms were ignored in order to arti¬cially
increase the Reynolds number on a coarse grid. To further amplify the
numerical Reynolds number, a vorticity con¬nement term was added as
an arti¬cial forcing function on the right-hand side of the equations. This
Figure 18.1. A warm smoke plume injected from left to right rises under the
vorticity con¬nement method was invented by Steinho¬; see, e.g., [161] for in¬‚uence of buoyancy. (See also color ¬gure, Plate 11.)
more details. Since the air was assumed to have constant density, equa-
tion (18.1) was not needed to model the air. On the other hand, since
equation (18.1) can be used to model any passive scalar, it was used to
independently track the concentration of smoke within the air. Figure 18.1
shows a warm smoke plume injected from left to right rising under the
e¬ects of buoyancy, while Figure 18.2 depicts the ¬‚ow of smoke past a

Figure 18.2. Small-scale eddies are generated as smoke ¬‚ows past a sphere. (See
also color ¬gure, Plate 12.)
218 19. Free Surfaces

the velocity are needed on the air side of the interface. These are de¬ned
by extrapolating the velocity ¬eld across the interface using equation (8.1).
Dirichlet pressure boundary conditions can be applied at the free surface by
setting the pressure to atmospheric pressure at every cell center located in
19 the air. This use of Dirichlet boundary conditions on the pressure reduces
the need to enforce the compatibility condition.
The atmospheric pressure boundary condition should be applied directly
Free Surfaces on the free surface, not at cell centers that are a ¬nite distance away
from the surface. Setting the pressure to atmospheric pressure at every
cell center in the air causes an overprediction of the pressure at the in-
terface itself. Raad et al. [136] reduces this problem to some degree by
using micro cells near the interface, lowering the size of the error constant
in this ¬rst-order accurate approximation of the boundary condition. Re-
cently, Cheng et al. [44] devised a fully second-order accurate method for
enforcing the atmospheric pressure boundary conditions on the free surface.
This is accomplished by setting the pressure at cells in the air to specially
calculated values that are lower than atmospheric pressure. This implicitly
enforces the atmospheric pressure boundary condition exactly on the free
surface to second-order accuracy. Notably, this method does not disturb
the symmetric nature of the coe¬cient matrix that needs to be inverted to
¬nd the pressures.
19.1 Description of the Model
Consider a two-phase incompressible ¬‚ow consisting of water and air. The
19.2 Simulating Water for Computer Graphics
two phase interface separating these ¬‚uids has rather complicated dynamics
that we will discuss later (in Chapter 21). In many situations the behavior of
One of the di¬culties associated with using level set methods to simulate
the heavier ¬‚uid signi¬cantly dominates the behavior of the lighter ¬‚uid. In
free surfaces (and ¬‚uids in general) is that level set methods tend to su¬er
these situations, the model can be simpli¬ed by treating the air as a simple
from mass loss in underresolved regions of the ¬‚ow. Foster and Fedkiw [71]
constant-pressure ¬‚uid that exerts no other stress (i.e., except pressure
addressed this problem in the context of splashing water by devising a new
forces) on the interface. This removes any relevant dynamic e¬ects from
numerical approach that hybridizes the particle method and the level set
the air, leaving only a simple uniform pressure force normal to the interface
method using the local interface curvature as a diagnostic. The curvature
independent of the interface topology or motion.
was used to monitor the interface topology by classifying regions of high
Harlow and Welch [79] used marker particles to identify which grid cells
curvature as underresolved. In these underresolved regions, the particles are
contained water and which grid cells contained air. Later, Raad et al. [136]
used to rebuild the level set function, reducing mass loss to a large degree.
improved the treatment of the pressure boundary conditions at the in-
On the other hand, in regions of relatively low curvature the high-order
terface, introducing a subcell treatment of the pressure. Chen et al. [39]
accurate level set solution is preferred. A sample calculation of a splash
improved the velocity boundary conditions at the free surface to obtain a
generated by a sphere is shown in Figure 19.1. In some regions of the ¬‚ow
more accurate ¬‚ow ¬eld. Furthermore, Chen et al. [40] reduced the need
the grid is too coarse to capture the splashing behavior, even with the aid
to resolve the three-dimensional volume with particles by introducing a
of the particles, and some particles will inevitably escape from the water
method that required particles only near the free surface itself.
side of the interface and appear on the air side. These escaped particles can
Particle locations give only a rough indication of the location of the
be used to generate an interesting spray e¬ect, as shown in Figure 19.2,
free surface. We instead prefer to use the level set function φ to more
where a slightly submerged ellipse slips through the water.
accurately model the free surface. Then the Navier-Stokes equations are
Although the local interface curvature is a good diagnostic of potential
solved on the water side of the interface only. In order to discretize both
mass loss in level set methods, the approach in [71] is still somewhat ad hoc
the level set equation and the Navier-Stokes equations, ghost cell values of
19.2. Simulating Water for Computer Graphics 219 220 19. Free Surfaces

Figure 19.1. A splash is generated as a sphere is thrown into the water. (See also
color ¬gure, Plate 13.)

in nature. A more thorough approach to hybridizing particle methods and
level set methods (the particle level set method) was presented in Chap-
ter 9. Figure 19.3 shows the highly realistic thin water sheet generated as
a sphere splashes into the water. This highlights the ability of the particle
level set method to accurately capture thin interfacial regions that may be
underresolved by the Cartesian grid. Figure 19.4 shows the impressive re-
sults generated using this method to model water ¬‚owing into a cylindrical
glass. A close-up view of the bottom of the glass is shown in Figure 19.5. Figure 19.2. An interesting spray e¬ect is generated as a slightly submerged ellipse
slips through the water. (See also color ¬gure, Plate 14.)
19.2. Simulating Water for Computer Graphics 221 222 19. Free Surfaces

Figure 19.3. A thin water sheet is generated by a sphere thrown into the water.
(See also color ¬gure, Plate 15.)

Figure 19.5. Pouring water into a cylindrical glass using the particle level set
method. (See also color ¬gure, Plate 17.)

Figure 19.4. Pouring water into a cylindrical glass using the particle level set
method. (See also color ¬gure, Plate 16.)
224 20. Liquid-Gas Interactions

20.2 Treating the Interface
Since the interface is a contact discontinuity moving with the local ¬‚uid ve-
locity, the Rankine-Hugoniot jump conditions imply that both the pressure
and the normal velocity are continuous across the interface. An incompress-
ible liquid can be thought of as the limiting case obtained by increasing
the sound speed of a compressible liquid to in¬nity. In this sense, an in-
compressible ¬‚uid can be thought of as a very sti¬ compressible ¬‚uid, in
20 fact, the sti¬est. The interface separating the compressible ¬‚ow from the
incompressible ¬‚ow is treated using the robust interpolation procedure out-
lined in Section 15.9. That is, the compressible gas is used to determine
Liquid-Gas Interactions the interface pressure, while the incompressible liquid is used to determine
the interface normal velocity.
Advancing the solution forward in time consists of four steps. First, the
entire incompressible velocity ¬eld is extrapolated across the interface. The
ghost cell values are used to ¬nd the intermediate incompressible velocity
¬eld V . Second, the entire compressible state vector is extrapolated across
the interface, and the extrapolated tangential velocity is combined with the
incompressible normal velocity to obtain a ghost cell velocity for the com-
pressible ¬‚uid. Then the compressible gas is updated in time. Third, the
level set function is advanced forward in time using the incompressible ve-
locity ¬eld only, since the interface velocity is de¬ned by the incompressible
¬‚ow. The extrapolated ghost cell values of the incompressible velocity ¬eld
are useful in this step. Fourth, the intermediate incompressible ¬‚ow veloc-
ity V is projected into a divergence-free state using the updated level set
location and the updated values of the compressible pressure as Dirichlet
20.1 Modeling
boundary conditions at the interface. This last step accounts for the in-
terface forces imposed by the pressure of the compressible ¬‚uid. Surface
Liquid-gas interactions, e.g., the vaporization and subsequent combustion
tension e¬ects are easily included in this last step using Dirichlet pressure
of liquid fuel droplets or the shock-induced mixing of liquids, are rather
boundary conditions of p = pc + σκ, where pc is the compressible pressure,
di¬cult problems in computational ¬‚uid dynamics. These problems ad-
σ is a constant, and κ is the local interface curvature. This accounts for
dress the interaction of liquid droplets with a compressible gas medium.
the jump in pressure due to surface tension forces, i.e., [p] = σκ.
There are three classical approaches to such problems: Both phases can be
Figure 20.1 shows an incompressible liquid droplet moving from left to
treated as compressible ¬‚uids (as we did in Chapter 15), both phases can
right in a compressible gas ¬‚ow. Notice the lead shock in the compressible
be treated as incompressible ¬‚uids (as we did in Chapter 21), or the gas can
gas. Figure 20.2 shows a shock wave impinging on an incompressible liquid
be treated as a compressible ¬‚uid while the liquid is treated as an incom-
droplet. A re¬‚ected wave can be seen to the left, and a faint transmitted
pressible ¬‚uid. A completely incompressible treatment can be ruled out any
wave can be seen to the right.
time one is interested in shock waves or other compressible phenomena. A
completely compressible treatment is not desirable, since a relatively high
sound speed in the liquid phase can impose a restrictive (and ine¬cient)
CFL condition. Moreover, a completely compressible approach is limited to
liquids for which there are acceptable models for their compressible evolu-
tion. To overcome these di¬culties, Caiden et al. [23] modeled the gas as a
compressible ¬‚uid and the liquid as an incompressible ¬‚uid. They coupled
a high-resolution shock-capturing scheme for the compressible gas ¬‚ow to
a standard incompressible ¬‚ow solver for the liquid phase.
20.2. Treating the Interface 225 226 20. Liquid-Gas Interactions

velocity field
velocity field 1










0 0.1 0.8
0 0.7
0.5 0.9
0.2 0.4 1
0.1 0.8
0 0.7
0.5 0.9
0.2 0.4 1

Figure 20.2. A shock wave impinging on an incompressible droplet producing
Figure 20.1. An incompressible droplet traveling to the right in a compressible
a re¬‚ected wave and a (very) weak transmitted wave. (See also color ¬gure,
gas ¬‚ow. Note the lead shock wave. (See also color ¬gure, Plate 18.)
Plate 19.)
228 21. Two-Phase Incompressible Flow

the jump in pressure due to surface tension, [p] = σκ, where σ is a constant
coe¬cient and κ is the local curvature of the interface. Using the immersed
boundary method to smear out the pressure across the interface leads to

21 continuity of the pressure, [p] = 0, and loss of all surface-tension e¬ects.
This was remedied by Brackbill et al. [18] in the context of volume of ¬‚uid
(VOF) methods and by Unverdi and Tryggvason [168] in the context of
Two-Phase Incompressible Flow front-tracking methods by adding a new forcing term to the right-hand
side of the momentum equations. In the context of level set methods (see
[160]) this new forcing term takes the form

, (21.3)
where δ(φ) is the smeared-out delta function given by equation (1.23). In
the spirit of the immersed boundary method, [18] referred to this as the
continuum surface force (CSF) method.
In the interest of solving for the pressure jump directly, Liu et
al. [106] devised a new boundary-condition-capturing approach for the
variable-coe¬cient Poisson equation to solve problems of the form
p = f, (21.4)
21.1 Introduction
where the jump conditions [p] = g and [(1/ρ) p · N ] = h are given. Here,
The earliest real success in the coupling of the level set method to problems ρ can be discontinuous across the interface as well. Figure 21.1 shows a
involving external physics came in computing two-phase incompressible typical example of the discontinuous solutions obtained using this method.
¬‚ow, in particular see the work of Sussman et al. [160] and Chang et al. Note that both the pressure and its derivatives are clearly discontinuous
[38]. In two-phase incompressible ¬‚ow, the Navier-Stokes equations are used across the interface. Kang et al. [91] applied this technique to multiphase
to model the ¬‚uids on both sides of the interface. Generally, the ¬‚uids incompressible ¬‚ow, illustrating the ability to solve these equations without
have di¬erent densities and viscosities, so these quantities are discontinuous smearing out the density, the viscosity, or the pressure across the interface.
across the interface. Both the discontinuous viscosity and surface tension Moreover, the δ(φ)σκN /ρ forcing term was not needed, since the pressure
forces cause the pressure to be discontinuous across the interface as well. jump was modeled directly. Figure 21.2 shows a water drop falling through
In addition, a discontinuous viscosity leads to kinks in the velocity ¬eld the air into the water. Here, surface tension forces cause the spherically
across the interface. shaped region at the top of the resulting water jet in the last frame of the
In [132], Peskin introduced the “immersed boundary” method for simu- ¬gure.
lating an elastic membrane immersed in an incompressible ¬‚uid ¬‚ow. This LeVeque and Li [102] proposed a second-order accurate sharp interface
method uses a δ-function formulation to smear out the numerical solution method to solve equation (21.4). In general, one needs to solve a linear sys-
in a thin region about the immersed interface. This concept has been used tem of equations of the form Ap = b, where p are the unknown pressures, A
by a variety of authors to solve a number of interface-related problems. For is the coe¬cient matrix, and b is the right-hand side. Unfortunately, the dis-
example, [160] de¬ned a numerically smeared-out density and viscosity as cretization in [102] leads to a complicated asymmetric coe¬cient matrix,
functions of the level set function, making this linear system di¬cult to solve. So far, this method has not
been applied to two-phase incompressible ¬‚ow. In contrast, the discretiza-
ρ(φ) = ρ’ + ρ+ ’ ρ’ H(φ), (21.1)
tion proposed in [106] leads to a symmetric coe¬cient matrix identical
µ(φ) = µ’ + µ+ ’ µ’ H(φ), (21.2) to the standard one obtained when both the pressure and its derivatives
are continuous across the interface. Adding the jump conditions only re-
where H(φ) is the numerically smeared-out Heaviside function de¬ned by
quires modi¬cation of the right-hand side, b. This allows the use of standard
equation (1.22). This removes all discontinuities across the interface, except
21.1. Introduction 229 230 21. Two-Phase Incompressible Flow




Figure 21.2. A water drop falls through the air into the water. Surface tension
forces cause the spherically shaped region at the top of the water jet in the last
frame. (See also color ¬gure, Plate 20.)

21.2 Jump Conditions

Applying conservation to the two-phase incompressible ¬‚ow interface
results in the jump conditions
®«  « 
1.5 ’0.5
N σκ
° T1  (pI ’ „ )N T » =  0  , (21.5)
· ( ρ p) = f (x, y), with [p] = g(x, y) and
Figure 21.1. Solution of
where T1 and T2 are orthogonal unit tangent vectors, I is the identity ma-
[ ρ p · N ] = h(x, y). Note the sharp resolution of the discontinuity. trix, and „ is the viscous stress tensor; see equation (18.6). Equation (21.5)
states that the net stress on the interface must be zero, since it has no
Since the ¬‚ow is viscous, the velocities and their tangential derivatives
black-box linear system solvers even in the presence of complicated jump are continuous across the interface, i.e.,
conditions. Recently, Li and Lai [103] extended [106] by adding a second-
[u] = [v] = [w] = 0, (21.6)
order accurate correction term to the method. This correction term is valid
in the presence of an immersed interface, raising the order of accuracy from [ u · T1 ] = [ v · T1 ] = [ w · T1 ] = 0, (21.7)
one to two in that instance. Unfortunately, a correction term for two-phase
[ u · T2 ] = [ v · T2 ] = [ w · T2 ] = 0. (21.8)
incompressible ¬‚ow does not yet exist.
21.2. Jump Conditions 231 232 21. Two-Phase Incompressible Flow

This leads to the jump condition This can be combined with equation (21.12) to obtain
· „)
p (
u · N, v · N, w · N · N = 0, (21.9) = , (21.14)
ρ ρ
which states that the normal derivative of the normal component of the which is equivalent to

<< . .

( 23)

. . >>

Copyright Design by: Sunlight webdesign