LINEBURG

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

sphere.

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

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0 0.1 0.8

0 0.7

0.6

0.5 0.9

0.2 0.4 1

0.3

0.1 0.8

0 0.7

0.6

0.5 0.9

0.2 0.4 1

0.3

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

δ(φ)σκN

, (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

1

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

10

5

0

Figure 21.2. A water drop falls through the air into the water. Surface tension

’5

forces cause the spherically shaped region at the top of the water jet in the last

frame. (See also color ¬gure, Plate 20.)

1

’10

21.2 Jump Conditions

0.5

3

Applying conservation to the two-phase incompressible ¬‚ow interface

0

2.5

results in the jump conditions

2

®« «

1.5 ’0.5

N σκ

1

° T1 (pI ’ „ )N T » = 0 , (21.5)

0.5

’1

0

0

T2

1

· ( ρ 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-

1

[ ρ 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

mass.

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

T

· „)

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

Copyright Design by: Sunlight webdesign