Francoise J.-P., Naber G.L., Tsun T.S. (editors) Encyclopedia of Mathematical Physics

Подождите немного. Документ загружается.

and a normal vector n

, the interface is given by the

straight line with normal n

such that area beneat h

the line in cell 

is equal to c

. More recently,

parabolic reconstructions of the interface have been

used to gain higher-order accuracy for the surface

tension force (e.g., the ‘‘parabolic reconstruction of

surface tension’’ or PROST algorith m).

Once the interface has been reconstructed, its

motion by the underlying flow field must be

modeled by a suitable advection algorithm. The

key here is that the explicit interface reconstruction

enables fluxes to be developed that exactly conserve

mass and do not diffuse the interface.

Capillary effects may be represented by the

continuous surface stress (Scardovelli and Zaleski

1999),

T ¼ð I  n  nÞjr

cj; F

sing

¼rT ½11

where

c is a smoothed version of the volume

fraction. For the flows in which the capillary force

is the dominant physical mechanism, the PROST

algorithm discussed above can be used to signifi-

cantly reduce spurious currents due to inaccurate

representation of surface tension terms and asso-

ciated pressure jump in normal stress.

The distribution form of the fluid equations [5] is

typically solved using a variant of the projection

method for incompressible single phase flows.

VOF methods are popular and have been used in

commercial multiphase flow codes, in models of

inkjet printers, flows with surfactants and in many

other applications (e.g., see Scardovelli and Zaleski

(1999) and James and Lowengrub (2004) and the

references therein). The principal advantage of VOF

methods is their inherent volume-conserving prop-

erty. Nevertheless, spurious bubbles and drops may

be created. The reconstruction of the interface from

the volume fractions and the computation of

geometric quantities such as curvature are typically

less accurate than other methods discussed here

since the curvature and normal vectors are obtained

by differentiating a nearly discontinuous function

(volume fraction).

Front-tracking methods The basic idea behind the

original front-tracking method is the use of two

grids as illustrated in Figure 2. One is a standard,

Eulerian finite difference mesh that is used to so lve

the fluid equations. The other is a discretized

interface mesh that is used to explicitly track the

interface and compute surface tension force which is

then transferred to the finite difference mesh via a

discrete delta-function. Front tracking was first

proposed by Richtmyer and Morton and further

developed by Glimm and co-workers.

A similar approach was taken by Unverdi and

Tryggvason (see Tryggvason et al. (2001) and Peskin

(2002) for recent reviews), who combined a moving

grid description of the interface with flow computa-

tions on a fixed grid. In this immersed-boundary

approach, all the fluid phases are treated together by

solving a single set of governing equations. This

method has its roots in the original marker-and-cell

(MAC) method, where marker particles are used to

identify each fluid and the immersed-boundary

method of Peskin and McQueen, that was designed

to track moving elastic boundaries in homogeneous

fluids.

The interface is represented discretely by Lagran-

gian markers that are connected to form a front

which lies within and moves through a stationary

Eulerian mesh.

In Tryggvason’s original implementation, the

basic structural unit is a line segment. Since the

interface moves and deforms during the computa-

tion, interface elements must occasionally be added

or deleted to maintain regularity and stability. In the

event of merging/breakup, elements must be relinked

to effect a change in topology.

The interface is represented using an ordered list

of marker particles x

= ((x

)

,(x

)

), 1  k  N.

Fluid 1

Fluid 2

(a)

f,k

(b)

i – 1/2, j + 1

i – 1/2, j

i, j – 1/2

i – 1/2, j

i + 1/2, j + 1

i , j + 1/2

(c)

Figure 2 (a) The basic idea in the front-tracking method is to use two grids – a stationary finite difference mesh and a moving

Lagrangian mesh, which is used to track the interface. (b). Blow-up of the subgrid control volume in (a). (c) Control volume for the

Eulerian mesh, 

i, j þ(1=2)

Interfaces and Multicomponent Fluids 137

The first step in this algorithm is the advection of the

marker particles. A simple bilinear interpolation is

used to find the velocity inside each grid cell (indicated

in Figure 2c). The marker particles are then advected in

a Lagrangian manner. Once the points have been

advected, a list of connected polynomials (p

(s), p

(s))

is constructed using the marker particles. This gives a

parametric representation of the interface, with s

typically an approximation of the arclength. Both

lists are ordered and thus identify the topology of the

interface. In later works, higher-order polynomials

have been used (e.g., cubic splines) and semi-Lagran-

gian evolutions have been implemented where other

tangential velocities have been used.

As the interface evolves, the markers drift along

the interface following tangential velocities an d

more markers may be needed if the interface is

stretched by the flow. Typically, the markers are

redistributed along the interface to maintain an

accurate interface repre sentation.

Next, we compute the surface tension force,

sing

ðx; tÞ¼

ðtÞ



ð x  x

ðsÞÞn

ds ½12

where the subscript f means values evaluated at the

interface (t) and s is arclength. The discrete

numerical implementation of this distribution onto

the fixed grid is in the form of a sum over interface

elements, x

f , k

ðxÞ¼

ð x  x

f ;k

Þs

½13

where s

is the average of the straight line

distances from the point x

f , k

to the two neighboring

points x

f , kþ1

and x

f , k1

as indicated by the subgrid

control volume shown in Figures 2a and 2b. The

delta-function is typically taken to be Peskin’s

discrete Dirac delta-function:

ðx x

f ;k

i¼1

1 þcos

½x

ðx

f ;i





if jx x

f ;k

j2h

0 otherwise [14]

Other higher -order alternative forms of the regular-

ized delt a-function using the product formula have

recently been proposed.

Using the Frenet relation, the surface tension force

on a short segment of the front is given by



ds ¼



ds ¼ ðt

 t

Þ½15

where A and B are the segment endpoints that lie

on the boundary of the subgrid control volume

(Figures 2a and 2b), and t

is a tangent vector

computed by fitting a polynomial to the endpoints

of each element.

In the case of flows with varying density and/or

viscosity between the fluid components, there is a

need to calculate the phase indicator function I(x, t)

(defined by interface geometry and position), which

has the value 0 in fluid 1 and 1 in fluid 2. The

indication function can be determined via the

solution of the equation

Iðx; tÞ¼r

ðtÞ

ð x  x

ðs; tÞÞds ½16

This equation is discretized on the Eulerian mesh

and a discrete delta-function (e.g., eqn [14]) is used.

The fluid properties such as density and viscosity are

determined via the indicator function, that is,

(x, t) = 

þ (

 

)I(x, t), etc.

As in the volume of fluid algorithm, the distribu-

tion form of the Navier–Stokes equations [5] are

typically solved using a version of Chorin’s projec-

tion method.

An alternative flow solver that can be used to

integrate the flow equations in the presence of an

interface is the immersed-interface method (IIM).

The IIM was developed by Leveque and Li (see the

review Li 2003), and can be used together with

front-tracking as well as level-set methods.

The IIM directly incorporates jump conditions for

the normal stress into the finite difference stencil. The

key idea of this method is to use the jump conditions

in Taylor series expansions of pressure and velocity

near interfaces to derive difference equations that

achieve pointwise second-order accuracy.

The principal advantage of front -tracking algo-

rithms is their inherent accuracy, due in part to the

ability to use a large number of grid points on the

interface. Front-tracking methods can be compli-

cated to implement, particularly in 3D, but give the

precise location and geometry of the interface. In

addition, explicit front tracking permit s more than

one interface to be present in a single computational

cell without coalescence, which can be important in

dense bubbly flows, emulsions, etc. One of major

handicaps of front-tracking methods is the difficulty

in modeling topolog ical changes of the interface

such as breakup and coalescence without ad hoc cut-

and-connect and reconnecting parameterized inter-

face (particularly, difficulties in 3D).

Interface-Capturing Methods

Level-set method Level-set methods, introduced by

Osher and Sethian (see the recent review papers

(Osher and Fedkiw 2001, Sethian and Smereka

138 Interfaces and Multicomponent Fluids

2003) and the recent texts (Osher and Fedkiw

2002, Sethian 1999)), are popular computational

techniques for tracking moving interfaces. These

methods rely on an implicit representation of the

interface as the zero set of an auxiliary function

(level-set function). The application of these meth-

ods to incompressible, multiphase flows started with

the work of Osher, M erriman, Sussman, Smereka,

Hou, and their collaborators.

In the level-set method, the level-set function

(x, t) is defined as follows (see Figure 3):

ðx; tÞ

> 0ifx 2 fluid 1

¼ 0ifx 2  ðthe interface between fluidsÞ

< 0ifx 2 fluid 2

(

and the evolution of  is given by



þ u r ¼ 0 ½17

which means that the interface moves with fluid.

To keep the interface geometry well resolved, the

level-set function  should be a distance function near

the interface. However, under the evolution [17] it

will not necessarily remain as such. We note that

special velocity extensions v off the interface (i.e.,

v = u at the interface, v 6¼ u away from interface)

have been recently developed to better maintain  as

a distance function (e.g., Sethian and Smereka (2003)

and Macklin and Lowengrub (2005)). Typically, a

reinitialization step (solving a Hamilton–Jacobi type

equation, eqn [18]) below, is performed to keep  as

a distance function near the interface while keeping

original zero-level set unchanged. More specifically,

given a level-set function, , at time t, the contours

are redistributed according to the steady-state solu-

tion of the equation

@

¼ S



ðÞð1 jrdjÞ; dðx; 0 Þ¼ðxÞ½18

where S



is the smoothed sign function defined as



ðÞ¼



ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ



þ 

½19

where  is usually is one or two grid lengths. After

solving eqn [18] to steady state (x, t) is then

replaced by d(x, 

steady

). Note that d(x, 

steady

)is

typically a good approximation of the signed

distance function.

The density and viscosity are defin ed as

ðÞ¼

þð

 

ÞH



ðÞ

and

ð Þ¼

þð

 

ÞH



ðÞ½20

where H



() is the smoothed Heaviside function

given by



ðÞ¼

0if<

1 þ







sinð=Þ



if jj

if >

The mollified delta-function is 



() = dH



=d. The

surface tension force is given as

sing

¼r

r

jrj







ðÞ

r

jrj

½21

The fluid equations [5] are solved using projecti on

methods, the IIM or the ghost-fluid (GF) method

(e.g., Osher and Fedkiw (2001, 2002) and Fedkiw

et al. (2003)). The GF method is similar to the IIM

in that jump discontinuities are incorporated in the

finite difference stencil. In the GF algorithm, subcell

resolution is used to mark the interface position and

the values of discontinuous quantities are artificially

extended to grid points neighboring the interface via

extrapolation. A fully second order accurate GF

method for moving interfaces has recently been

developed (Macklin and Lowengrub 2005).

Applications of the level-set method include

multiphase flows, viscoelastic fluid flows and fluid–

structure interactions (e.g., see the reviews Osher and

Fedkiw (2001, 2002), Sethian (1999),andSethian

and Smereka (2003)).

Advantages of the level-set algorithm include the

simplicity with which it can be implemented, the

ability to capture merging and breakup of interfaces

automatically, and the ease with which the interface

geometry can be described using the level-set

function. A disadvantage of the level-set method is

that mass is not conserved.

Accurate numer ical simulations of multiphase

flow and topology transitions require the computa-

tional mesh to resolve both the macroscales (e.g.,

droplet size, flow geometry) and the microscales to

accurately capture local interface geometries near

contact region, van der Waals forces, surfactant

distribution, and Marangoni stresses. Adaptive mesh

–2 –1.5 –1 –0.5 0 0.5 1 1.5 2

–2

–1.5

–1

–0.5

0.5

1.5

φ < 0

> 0

φ = 0

Fluid 1

Fluid 2

–2

–1

–2

–1.5

–1

–0.5

0.5

(a) (b)

Figure 3 (a) Zero contour of  representing the interface .

(b) Surface of  with zero contour.

Interfaces and Multicomponent Fluids 139

algorithms have recently been used greatly to

increase accuracy and computational efficiency in

level-set methods. Typically, the methods involve

Cartesian adaptive mesh refinement. Problems

tackled using this approach include droplet forma-

tion in inkjet printers and wake development behind

a ship. Another approach, recently developed, is to

use adaptive unstructured mesh refinement (Zheng

et al. 2005), as shown in Figure 4, in which the

impact of a drop onto a fluid interface is captured.

Hybrid Methods

More recently, a number of hybrid methods, which

combine good featu res of each algorithm, have been

developed. These include coupled level-set volume-

of-fluid (CLSVOF) algorithms, particle level-set

methods, marker-VOF methods and level-contour

front-tracking methods.

Level-set and VOF methods have recently been

combined. The volume fraction is used to maintain

volume conservation, while the level-set function is

used to describe the interface geometry. After every

time step, the volume-fraction function and level-set

function are made compatible. The coupling

between the level-set function  and the VOF

function c occurs through the normal of the

reconstructed interface and through the fact that

the level-set function is reset to the exact signed

normal distance to the reconstructed interface

(where the area below the reconstructed interface is

given by the volume-fraction function).

In the particle leve l-set method, Lagrangian

disconnected marker particles are rand omly posi-

tioned near the interface and are passively advected

by the flow in or der to rebuild the level-set function

in under-resolved zones, such as high-curvature

regions and near filaments. In these regions, the

standard nonadaptive level-set method regularizes

excessively the interface structure and mass is lost.

The use of marker particles significantly ameliorates

these difficulties.

Recently, a hybrid method has been developed,

which uses both marker particles, to reconstruct and

move the interface, and the volum e-fraction function

to conserve volume. In this approach, a smooth

motion of the interface, typical of marker methods is

obtained together with volume conservation, as in

standard VOF methods. This work improves both

the accuracy of interface tracking, when compared

to standard VOF methods, and the conservation of

mass, with respect to the original marker method.

Finally, a hybrid method that combines a level

contour reconstruction technique with front-tracking

methods has recently been developed to auto-

matically model the merging and breakup of inter-

faces in three-dimensional flows.

Phase-Field Method

Phase-field, or diffuse-interface, model s are an

increasingly popular choice for modeling the motion

of multiphase fluids (see Anderson et al. (1998) for a

recent review). In the phase-field model, sharp fluid

interfaces are replaced by thin but nonzero thickness

transition regions where the interfacial forces are

smoothly distributed. The basic idea is to introduce

a conserved order parameter (e.g., mass concentra-

tion) that varies continuously over thin interfacial

layers and is mostly uniform in the bulk phases (see

Figure 5).

For density-matched binary liquids (let  = 1

for simplicity), the coupling of the convective

Cahn–Hilliard equation for the mass concentration

with a modified momentum equation that includes a

phase-field-dependent surface force is known as

Model H (Hohenberg and Halperin 1977). In the

case of fluids with different densities a phase-field

model has been proposed by Lowengrub and

Truskinovsky. Complex flow morphologies and

topological transitions such as coalescence and

interface breakup can be captured naturally and in

a mass-conservative and energy-dissipative fashion

since there is an associated free energy functional.

–1

–2

–3

–4

–5

–4 –2 0 2 4

–3.2

–3.25

–3.3

–3.35

–3.4

–

0.75 – 0.7 – 0.65 – 0.6 – 0.55

–1 –0.5 0 0.5 1

–2

–2.25

–2.5

–2.75

–3

–3.25

–3.5

–3.75

–4

–3.275

–3.28

–3.285

–3.29

–0.655 –

0.65 – 0.645 – 0.64

Figure 4 Each of the first three figures has a boxed region that is magnified in the next figure. The rates of magnification are 5, 10,

40/3, respectively. The meshes in the figure are used to simulate the drop-impacting interface problem. Source: Zheng X, Anderson A,

Lowengrub JS, and Cristini V (unpublished).

140 Interfaces and Multicomponent Fluids

The phase field is governed by the following

advective Cahn–Hilliard equation:

þ u rc ¼rðMðcÞrÞ½22

 ¼ F

ðcÞ

c ½23

where M(c) = c(1  c) is the mobility, F(c) =

(1=4)c

(1  c)

is a Hel mholtz free energy that

describe the coexistence of immiscible phases, and

 is a measure of interface thickness and    (see

Figure 5). It can be shown that in the sharp interface

limit  !0, the classi cal Navier–Stokes system

equations and jump conditions are recovered.

The singular surface tension force is F

sing

6

ﬃﬃﬃ

r(rc rc), where  is the surface ten-

sion coefficient. An alternative surface tension force

formulation based on the CSF is F

sing

= 6

ﬃﬃﬃ

r

(rc=jrcj)jrcjrc.

Recently, very efficient nonli near multigrid meth-

ods have been developed to solve implicit discretiza-

tions of the Cahn–Hilliard equation (e.g., Kim et al.

(2004)). These schemes have been combined with

projection methods to solve the Navier–Stokes

equations to perform simulations of multiphase

flows.

An example of simulation of liquid thread breakup

using a phase-field method is shown in Figure 6.

A long cylindrical thread of a viscous fluid 1 is in an

infinite mass of another viscous fluid 2. If the thread

becomes varicose with wavelength , the equilibrium

of the column is unstable, provided  exceeds the

circumference of the cylinder. This is the Rayleigh

capillary instability that results in surface-tension-

driven breakup of the thread.

An advantage of the phase-field approach is that it

is straightforward to include more complex physical

effects. For example, the binary model can be

straightforwardly extended to describe three-

component flows as follows.

Consider a ternary mixture and denote the

composition of components 1, 2, and 3, expressed

as mass fractions, by c

, c

, and c

, respectively.

Therefore,

i¼1

¼ 1; 0  c

 1 ½24

The composition of a ternary mixture (A, B, and C)

can be mapped onto an equilateral triangle (the

Gibbs triangle (Porter and Easterling 1993)) whose

corners represent 100% concentration of A, B, or C

as shown in Figure 7a. Mixtures with components

lying on lines parallel to

BC contain the same

percentage of A, those with lines parallel to

AC have

the same percentage of B concentration, and

analogously for the C concentration. In Figure 7a,

the mixture at the position marked ‘’ contains 60% A,

10% B, and 30% C. Because the concentrations sum

to unity, only two of them need to be determined,

say c

, c

The evolution of c

and c

is governed by the

following advective ternary Cahn–Hilliard equation:

þ u rc

¼rðMðc

; c

Þr

Þ½25

–2 –1.5 –1 –0.5 0 0.5 1 1.5 2

–0.2

0.2

0.4

0.6

0.8

Figure 5 A concentration prome across an interface with

interface thickness, .

Figure 6 Time evolution leading to multiple pinch-offs. The

evolution is from top to bottom and left to right. The domain is

axisymmetric, the initial velocities are zero everywhere, and the

concentration field is given by c(r, z) = 0.5(1  tanh ((r  0.5

0.05 cos ( z))=(2

ﬃﬃﬃ

))) on =(0, )  (0, 2). Densities are

matched and viscosity ratio is 0.5.

(a) (b)

Figure 7 (a) Gibbs triangle. (b) Contour plot of the free energy

F (c

, c

) on the Gibbs triangle.

Interfaces and Multicomponent Fluids 141

þ u rc

¼rðMðc

; c

Þr

Þ½26



@Fðc

; c

 

c

 0:5

c

½27



@Fðc

; c

 0:5

c

 

c

½28

where M(c

, c

) =

i<j

is the mobility and

F(c

, c

) is the Helmholtz free energy that can be

used to model the miscibility of the components. An

example of a free energy (used in the simulation

shown in Figure 8 below) for which fluids 1 and 3

are immiscible and fluid 2 is preferentially miscible

with fluid 3 is:

Fðc

; c

Þ¼2c

1  c

 c

ðÞ

þ c

þ 0:2ðÞc

 0:2ðÞ

þ 1:2  c

 c

ðÞ

 0:4

ðÞ

The contours of F on the Gibbs triangle are shown

in Figure 7b.

The singular surface tension force is F

sing

6

ﬃﬃﬃ



i = 1



r(rc

rc

), where the physical

surface tension coefficients 

between two fluids i

and j are decom posed into the phase-specific surface

tensions 

such that 

= 

þ 

As a demonstration of the evolution possible in

partially miscible liquid systems, we present an

example in which there is a gravity-driven

(Rayleigh–Taylor) instability that enhances the

transfer of a preferentially miscible contaminant

from one immiscible fluid to another in 2D. In this

system, the ternary Cahn–Hilliard system is solved

using nonlinear multigrid methods and a projection

method (Kim and Lowengrub (in press)) is used to

solve the flow equations [5].

In Figure 8 (first column), the top half of the domain

initially consists of a mixture of fluids 1 and 2,

and the bottom half consists of fluid 3, which is

immiscible with fluid 1. The contours of c

, c

,andc

are visualized in gray-scale where darker regions

denote larger values of c

, c

,andc

,respectively.

In the top row, the contours of fluid 1 are shown, the

middle and bottom rows correspond to fluids 2 and 3,

respectively.

Fluid 2 is preferentially miscible with fluid 3.

Fluid 1 is assumed to be the lightest and fluid 2 the

heaviest. The density of the 1/2 mixture is heavier

than that of fluid 3, so the density gradient induces

the Rayleigh–Taylor instability.

The evolution of the three phases is shown in

Figure 8. As the simulation begins, the 1/2 mixture

falls and fluid 2 diffuses into fluid 3. A characteristic

Rayleigh–Taylor (inverted) mushroom forms, the

Figure 8 Evolution of concentration of fluid 1 (top row), 2 (middle row), and 3 (bottom row). The contours of c

, c

, and c

are

visualized in gray-scale where darker regions denote larger values of c

, c

, and c

, respectively.

142 Interfaces and Multicomponent Fluids

surface area of the 1/3 interface increases, and

vorticity is generated and shed into the bulk.

As fluid 2 is diffused from fluid 1, the pure fluid

1 rises to the top as shown in Figure 8. Imagining

that fluid 2 is a contaminant in fluid 1, this

configuration provides an efficient m eans of cleans-

ing fluid 1 since the buoyancy-driven flow enhances

the diffusional transfer of fluid 2 from fluid 1 to

fluid 3.

The advantages of the phase-field method are:

(1) topology changes are automatically described;

(2) the composition field c has a physical meaning

not only near interface but also in the bulk phases;

(3) complex physics can easily be incorporated into

the framework, the methods can be straightforwardly

extended to multicomponent systems, and miscible,

immiscible, partially miscible, and lamellar phases

can be modeled.

Associated with diffuse interfaces is a small scale

, proportional to the width of the interface. In real

physical systems describing immiscible fluids,  can

be vanishingly small. However, for numerical

accuracy  must be at least a few grid lengths in

size. This can make computations expensive. One

way of ameliorating this problem is to adaptively

refine the grid only near the transition layer. Such

methods are under development by various research

groups.

Phase-field methods have been used to model

viscoelastic flow, therm ocapillary flow, spinodal

decomposition, the mixing and inte rfac ial stret ch-

ing, in a shear flow, droplet breakup process,

wave-breaking and sloshing, the fluid motion near

a moving contact line, and the nucleation and

annihilation of an equilibrium droplet (see the

references in the review paper Anderson et al.

(1998)).

Conclusions and Future Directions

In this paper we have reviewed the basic ideas of

interface-tracking and interface-capturing methods

that are critical in simulating the motion of inter-

faces in multicomponent fluid flows. The differences

between these various formulations lie in the

representation and the reconstruction of interfaces.

The advantages and disadvantages of the algorithms

have been discussed. While there has been much

progress on the development of robust multifluid

solvers, there is much more work to be done.

Promising future directions for research include the

incorporation of adaptive mesh refinement into the

algorithms and the development of efficient hybrid

schemes that combine the best features of individual

methods.

See also: Breaking Water Waves; Capillary Surfaces;

Fluid Mechanics: Numerical Methods; Incompressible

Euler Equations: Mathematical Theory; Inviscid

Flows; Non-Newtonian Fluids; Partial Differential

Equations: Some Examples; Viscous Incompressible

Fluids: Mathematical Theory; Vortex Dynamics.

Further Reading

Anderson DM, McFadden GB, and Wheeler AA (1998) Diffuse-

interface methods in fluid mechanics. Ann. Rev. Fluid Mech.

30: 139–165.

Cristini V, Blawzdziewicz J, and Loewenberg M (2001) An

adaptive mesh algorithm for evolving surfaces: simulations of

drop breakup and coalescence. Journal of Computational

Physics 168: 445–463.

Fedkiw RP, Sapirop G, and Shu C-W (2003) Shock capturing,

level sets and PDE based methods in computer vision and

image processing: a review of Osher’s contributions. Journal

of Computational Physics 185: 309–341.

Hohenberg PC and Halperin BI (1977) Theory of dynamic critical

phenomena. Reviews of Modern Physics 49: 435–479.

Hou TY, Lowengrub JS, and Shelley MJ (2001) Boundary integral

methods for multicomponent fluids and multiphase materials.

Journal of Computational Physics 169: 302–362.

James AJ and Lowengr ub J (2004) A sur fa ctan t-c onse rving

volume-of-fluid method for interfacial flows with insoluble

surfactant. Journal of Computational Physics 201:

685–722.

Kim JS, Kang KK, and Lowengrub JS (2004) Conservative

multigrid methods for Cahn–Hilliard fluids. Journal of

Computational Physics 193: 511–543.

Kim JS and Lowengrub JS Phase field modeling and simulation of

three-phase flows, Int. Free Bound (in press).

Li Z (2003) An overview of the immersed interface method and

its applications. Taiwanese J. Math. 7: 1–49.

Macklin P and Lowengrub JS (2005) Evolving interfaces via

gradients of geometry-dependent interior Poisson problems:

application to tumor growth. Journal of Computational

Physics 203: 191–220.

Osher S and Fedkiw RP (2001) Level set methods: an overview

and some recent results. Journal of Computational Physics

169: 463–502.

Osher S and Fedkiw RP (2002) Level Set Methods and Dynamic

Implicit Surfaces. Springer.

Peskin CS (2002) The immersed boundary method. Acta

Numerica 1–39.

Porter DA and Easterling KE (1993) Phase Transformations in

Metals and Alloys. Van Nostrand Reinhold.

Pozrikidis C (2001) Interfacial dynamics for Stokes flow. Journal

of Computational Physics 169: 250–301.

Scardovelli R and Zaleski S (1999) Direct numerical simulation of

free-surface and interfacial flow. Annu. Rev. Fluid Mech. 31:

567–603.

Sethian JA (1999) Level-Set Methods and Fast Marching

Methods: Evolving Interfaces in Computational Geometry,

Fluid Mechanics, Computer Vision and Materials Science.

Cambridge, MA: Cambridge University Press.

Sethian JA and Smereka P (2003) Level set methods for fluid

interfaces. Annu. Rev. Fluid Mech. 35: 341–372.

Interfaces and Multicomponent Fluids 143

TryggvasonG,BunnerB,EsmaeeliA,JuricD,Al-RawahiN

et al. (2001) A front-tracking method for the computations

of multiphase flow. Journal of Computational Physics 169:

708–759.

Zheng X, Lowengrub J, Anderson A, and Cristini V (2005)

Adaptive unstructured volume remeshing II. Application to

two- and three-dimensional level-set simulations of multiphase

flow. Journal of Computational Physics 208: 626–650.

Intermittency in Turbulence

J Jime´ nez, Universidad Politecnica de Madrid,

Madrid, Spain

Introduction

Intermittency has several meanings in turbulence.

The oldest one, now most often label ed ‘‘external’’

or ‘‘large-scale’’ intermittency, refers to the coex-

istence of turbulent and laminar regions in inho-

mogeneous turbulent flows, such as in boundary

layers or in free shear layers. In those cases, the

interface between lami nar irrotational flow and

turbulent vortical fluid is typically sharp and

corrugated. An observer sitting near the edge of

the layer is immersed in turbulent fluid only part of

the time.

The intermittency coefficient  measures the

fraction of turbulent fluid over the sampling

universe over which the statistics are taken. For

example, in a boundary layer such as that in

Figure 1, the intermittency coefficient as a function

of wall distance measures the fraction of turbulent

fluid at a given distance from the wall. External

intermittency is important in an y attempt to model

realistic turbulent flows, which are almost always

inhomogeneous. Consider, for example, the classical

homogeneous relation in eqn [1] between the mean

kinetic energy K of the turbulent fluctuations and

the energy dissipation rate " :

" ¼ C

3=2

½1

where L is the length scale of the largest eddies, and

C  0.1 is an experimentally determined constant.

Such relations are often implicit in turbulent models,

and they have to be modified to account for

intermittency. Equation [1] only holds within the

turbulent regions where the energy and the dissipa-

tion rates are K

and "

, while the overall mean

values used in the modeling conservation equations

are K = K

and " = "

. The true overall relation

should therefore be

" ¼ C

1=2

3=2

½2

which may differ substantially from eqn [1],

especially near the edge of the layer. Experimental

values and rough theoretical estimates for the

distribution of the intermittency coefficient are

available for most practical turbulent flows.

Internal Intermittency

While the external intermittency just described is

probably the most important one from the point of

view of applications, it is not the most interesting

from the theoretical point of view. Turbulence is a

multiscale phenomenon which is inhomogeneous

at all length scales, from the largest ones to the

inner viscous cutoff (see Turbulence Theories).

Moreover, this inhomogeneity goes beyond what

could be expected just from the s tatistics of a

random process. Consider, for example, the velo-

city difference u between two points separated

by a distance r. The original Kolmogorov formula-

tion of the energy c ascade assumes that the

probability density function (PDF), p(u), is a

universal function in the inertial range of scales,

whose only parameter is a velocity scale depending

on r. It then follows from Kolmogorov’s analysis

that

pðuÞ¼F u= ð



"rÞ

1=3

½3

where



" is the average energy transfer rate across

scales per unit mass, and the average

( ) is taken

either over the whole flow or over a suitably designed

ensemble of experiments. In an equilibrium system,

Flow

Figure 1 Sketch of a turbulent boundary layer, and of the

associated intermittency factor. An observer such as A, at a

distance y from the wall, only sees turbulent flow for a fraction 

of the time.

144 Intermittency in Turbulence

global energy conservation implies that



" is equal to

the average viscous dissipation per unit mass:



" ¼ 

jruj

½4

In eqn [4], the kinematic viscosity of the fluid is , and

jruj is the L

-norm of the velocity gradient tensor.

Equation [3] is valid as long as the separation r is

much larger than the Kolmogorov viscous cutoff

 = (



1=4

, and much smaller than the integral

scale of the largest eddies L

= u



", where u

is the

root-mean-square value of the fluctuations of one

velocity component. The extent of this iner tial range

is a function of the Reynolds number Re

= u

= :

= ¼ Re

3=4

½5

The strict similarity hypothesis in eqn [3] is not well

satisfied by experiments. While the velocity distribu-

tion at a given point is approximately Gaussian,

Figure 2a shows that the v elocity increments become

increasingly non-Gaussian as the spatial separation

is made much smaller than L

. It was also soon

noted that the dependence of eqn [3] on a single

parameter such as



" was theoretically suspect, since

it is difficult to see how the PDFs of a whole set of

local properties, such as the u for different

intervals, could depend only on a single global

property. Kolmogorov himself sought to bypass that

difficulty by substituting eqn [3] by a ‘‘refined

similarity’’ hypothesis,

pðuÞ¼F u= ð"

rÞ

1=3

½6

where "

is no longer a global average, but the mean

value of the dissipation over a ball of radius of order

r centered at the midpoint of the interval. This

refined similarity is better sati sfied by experiments

(see Figure 2b), although, from the practical poin t of

view, it just transfers the problem of characterizing

u to that of characterizing the statistics of "

It has become customary to measure the behavior

of p(u) in terms of its structure funct ions,

SðnÞ¼

1

u

pðuÞdu ½7

which can be normalized as generalized flatness

factors,

ðnÞ¼SðnÞ=Sð2Þ

n=2

½8

It follows from the strict similarity hypothesis [3] that

SðnÞr

n=3

½9

and that all the  (n) should be independent of the

separation.

For example, the fourth-order flatness of a

Gaussian distribution is (4) = 3. Figure 3 shows

that this is not true. The flatness increases as the

separation decreases, and it only levels off at lengths

of the order of the Kolmogorov viscous scale. For

separations in that viscous range the flow is smooth,

u  (@

u)r,and

ðnÞ

ð@

uÞ

=ð@

uÞ

n=2

½10

It follows from eqn [10] and from Figure 3 that the

velocity gradients become increasingly non-Gaussian

as L

and  separate at high Reynolds numbers. The

velocity differences across intervals which are large

with respect to  also become very non-Gaussian

when r  L

Because the velocity difference between two

points which are not too close to each other can be

expressed as the sum of velocity differences over

subintervals, a loose application of the central limit

–10 –5

5 10

–1

–2

–3

–4

Δu/(εr )

1/3

PDF

(a)

–10

0510

–1

–2

–3

–4

Δu/(ε

r )

1/3

PDF

(b)

–5

Figure 2 PDFs of the differences of the velocity component in the direction of the separation (for separations in the inertial range of

scales). r =L

= 0:020:36, increasing by factors of 2; equivalent to r = = 1803000: Nominally isotropic turbulence at Reynolds

number Re

= 10

.. (a) u is normalized with the global energy dissipation rate



"; distributions are wider as the separation decreases.

(b) u is scaled with the locally averaged dissipation over the separation interval. Data courtesy of H Willaime and P Tabeling.

Intermittency in Turbulence 145

theorem would suggest that its PDF should be

roughly Gaussian. The key conditions for that to

happen are that the summands should be mutually

independent, that their magnitudes should be com-

parable, and that each of them has a probability

distribution with a finite variance. The first of those

three conditions is probably a good approximation

if the separation is much longer than the viscous

cutoff, but the second one depends on the structure

of the flow. The experimental non-Gaussian beha-

vior suggests the existence of occasional very strong

velocity jumps. In the viscous range of scales, those

structures have been identified both experimentally

and numerically as very strong linear vortices, in

whose neighborhoods the strongest gradients are

generated. An example of a tangle of such structures

is shown in Figure 4.

In another example, the vorticity in decaying

two-dimen sional turbulence concentrates very

quickly into relatively few strong compact vortices,

which are stable except when they interact with

each other. The velocity field is dominated by them,

and the flatness of the velocity increme nts reaches

values of the order of (4) 50–100, even at

moderate Reynolds numbers. That case is interest-

ing because somethi ng can be said about the

probability distribution of the velocity g radients.

We have noted that the PDF of a sum of mutually

comparabl e independent random variables with

finite variances tends to Gaussian when the number

of summands is la rge. This well-known theorem is a

particular case of a more general result about sums

of random variables whose incomplete second

moments diverge as



ðsÞ¼

s

pðxÞdx  s

2

when s !1 ½11

When 0 < 2, the sums of such variables tend

to a family of ‘‘stable’’ distributions parametrized by

. The Gaussian case is the limit of that family when

 = 2. In the case of two-dimensional vortices with

very small cores, the velocity gradients at a distance

R from the center of the vortex behave as 1=R

.If

we take s in eqn [11] to be one of those velocity

derivatives, its probability distribution is propor-

tional to the area covered by gradients with a given

magnitude, and



ðsÞ

1=2

4

2R dR  s

1

½12

The velocity derivatives at any point, which are

the sums of the velocity derivatives induced by all

the randomly distributed neighboring vortices,

should therefore be distributed accordi ng to the

stable distribution with  = 1, which is Cauchy’s

pðsÞ¼

ðc

þ s

½13

This distribution has no moments for n >1. Its

tails decay as s

2

, and the distribution of the

gradients e ssentially reflects the properties of the

closest vortex. In real two-dimensional turbulent

flows, the distribution [13] is followed fairly well,

but its e xtreme tails only reach to the maximum

values of the velocity gradient found within the

viscous vortex cores, which are not exactly point

vortices.

–5

–4

–3

–2

–1

r /L

σ(4)

–0.12

Gaussian

Figure 3 Fourth-order flatness of the differences of the

velocity component in the direction of the separation, for

separations in the inertial range of scales, r=L

= 0:5to

r= = 2.. The Reynolds numbers of the different flows range

from Re

= 1800 to 10

.. Data in part courtesy of H Willaime,

P Tabeling, and R A Antonia.

Figure 4 Intense vortex tangle in the logarithmic layer of a

turbulent channel. The vortex diameters are of the order of 10,

and the size of the bounding box is of the order of the channel

width. Reproduced with permission of J C del A

lamo.

146 Intermittency in Turbulence