Главная Direct Numerical Simulations of Gas-Liquid Multiphase Flows

# Direct Numerical Simulations of Gas-Liquid Multiphase Flows

, ,
Accurately predicting the behaviour of multiphase flows is a problem of immense industrial and scientific interest. Modern computers can now study the dynamics in great detail and these simulations yield unprecedented insight. This book provides a comprehensive introduction to direct numerical simulations of multiphase flows for researchers and graduate students. After a brief overview of the context and history the authors review the governing equations. A particular emphasis is placed on the 'one-fluid' formulation where a single set of equations is used to describe the entire flow field and interface terms are included as singularity distributions. Several applications are discussed, showing how direct numerical simulations have helped researchers advance both our understanding and our ability to make predictions. The final chapter gives an overview of recent studies of flows with relatively complex physics, such as mass transfer and chemical reactions, solidification and boiling, and includes extensive references to current work.
Год: 2011
Издательство: Cambridge University Press
Язык: english
Страниц: 337
ISBN 10: 0521782406
ISBN 13: 9780521782401
Series: Cambridge Monographs on Applied and Computational Mathematics
File: PDF, 7.19 MB

## You may be interested in

### Mathematics of Large Eddy Simulation of Turbulent Flows

Год: 2006
Язык: english
File: PDF, 3.17 MB

### Fluid dynamics: theory, computation, and numerical simulation

Год: 2009
Язык: english
File: PDF, 10.24 MB

## Most frequently terms

You can write a book review and share your experiences. Other readers will always be interested in your opinion of the books you've read. Whether you've loved the book or not, if you give your honest and detailed thoughts then people will find new books that are right for them.
1

### آنها زنده اند

Год: 1354
Язык: persian
File: PDF, 9.05 MB
2

### گذرگاه امواج

Язык: persian
File: PDF, 15.77 MB
```DIRECT NUMERICAL SIMULATIONS
OF GAS–LIQUID MULTIPHASE FLOWS
Accurately predicting the behavior of multiphase flows is a problem of immense
industrial and scientific interest. Using modern computers, researchers can now
study the dynamics in great detail, and computer simulations are yielding unprecedented insight. This book provides a comprehensive introduction to direct numerical simulations of multiphase flows for researchers and graduate students.
After a brief overview of the context and history, the authors review the governing equations. A particular emphasis is placed on the “one-fluid” formulation,
where a single set of equations is used to describe the entire flow field and interface terms are included as singularity distributions. Several applications are
discussed, such as atomization, droplet impact, breakup and collision, and bubbly
flows, showing how direct numerical simulations have helped researchers advance
both our understanding and our ability to make predictions. The final chapter
gives an overview of recent studies of flows with relatively complex physics, such
as mass transfer and chemical reactions, solidification, and boiling, and includes
extensive references to current work.
G R É T A R T R Y G G V A S O N is the Viola D. Hank Professor of Aerospace and
Mechanical Engineering at the University of Notre Dame, Indiana.

is an Associate Professor in the Dipartimento di Ingegneria Energetica, Nucleare e del Controllo Ambientale (DIENCA) of the Università degli Studi di Bologna.
RUBEN SCARDOVELLI

S T É P H A N E Z A L E S K I is a Professor of Mechanics at the Université Pierre et
Marie Curie (UPMC) in Paris and Head of the Jean Le Rond d’Alembert Institute
(CNRS UMR 7190).

DIRECT NUMERICAL SIMULATIONS
OF GAS–LIQUID MULTIPHASE FLOWS
GR ÉTAR TRYGGVASON
University of Notre Dame, Indiana

RUBEN SCARDOVELLI
Università degli Studi di Bologna

ST ÉPHANE ZALESKI
Université Pierre et Marie Curie, Paris 6

CAMBRIDGE UNIVERSITY PRESS

Cambridge, New York, Melbourne, Madrid, Cape Tow; n,
Singapore, São Paulo, Delhi, Tokyo, Mexico City
Cambridge University Press
The Edinburgh Building, Cambridge CB2 8RU, UK
Published in the United States of America by Cambridge University Press, New York
www.cambridge.org
Information on this title: www.cambridge.org/9780521782401

C

Cambridge University Press 2011

This publication is in copyright. Subject to statutory exception
and to the provisions of relevant collective licensing agreements,
no reproduction of any part may take place without the written
permission of Cambridge University Press.
First published 2011
Printed in the United Kingdom at the University Press, Cambridge
A catalogue record for this publication is available from the British Library
ISBN 978-0-521-78240-1 Hardback

Cambridge University Press has no responsibility for the persistence or
accuracy of URLs for external or third-party internet websites referred to
in this publication, and does not guarantee that any content on such
websites is, or will remain, accurate or appropriate.

Contents

Preface

page ix

1

Introduction
1.1 Examples of multiphase flows
1.2 Computational modeling

1
3
7
18

2

Fluid mechanics with interfaces
2.1 General principles
2.2 Basic equations
2.3 Interfaces: description and definitions
2.4 Fluid mechanics with interfaces
2.5 Fluid mechanics with interfaces: the one-fluid formulation
2.6 Nondimensional numbers
2.7 Thin films, intermolecular forces, and contact lines
2.8 Notes

21
21
22
30
36
41
42
44
47

3

Numerical solutions of the Navier–Stokes equations
3.1 Time integration
3.2 Spatial discretization
3.3 Discretization of the advection terms
3.4 The viscous terms
3.5 The pressure equation
3.6 Velocity boundary conditions
3.7 Outflow boundary conditions
3.9 Summary
3.10 Postscript: conservative versus non-conservative form

50
51
55
59
61
64
69
70
71
72
73

4

4.1 Notations

75
76
v

vi

Contents

4.2
4.3
4.4
4.5
4.6
4.7
4.8

The volume-of-fluid (VOF) method
Front tracking
The level-set method
Phase-field methods
The CIP method
Summary

77
81
84
87
90
91
93

5

The volume-of-fluid method
5.1 Basic properties
5.2 Interface reconstruction
5.3 Tests of reconstruction methods
5.5 Tests of reconstruction and advection methods
5.6 Hybrid methods

95
95
98
106
108
122
128

6

6.1 The structure of the front
6.2 Restructuring the fronts
6.3 The front-grid communications
6.5 Constructing the marker function
6.6 Changes in the front topology
6.7 Notes

133
134
143
145
150
152
158
160

7

Surface tension
7.1 Computing surface tension from marker functions
7.2 Computing the surface tension of a tracked front
7.3 Testing the surface tension methods
7.4 More sophisticated surface tension methods
7.5 Conclusion on numerical methods

161
161
168
177
181
186

8

Disperse bubbly flows
8.1 Introduction
8.2 Homogeneous bubbly flows
8.3 Bubbly flows in vertical channels
8.4 Discussion

187
187
189
194
201

9

Atomization and breakup
9.1 Introduction
9.2 Thread, sheet, and rim breakup
9.3 High-speed jets
9.4 Atomization simulations

204
204
205
214
219

Contents

vii

10

Droplet collision, impact, and splashing
10.1 Introduction
10.2 Early simulations
10.3 Low-velocity impacts and collisions
10.4 More complex slow impacts
10.5 Corolla, crowns, and splashing impacts

228
228
229
229
232
235

11

Extensions
11.1 Additional fields and surface physics
11.2 Imbedded boundaries
11.3 Multiscale issues
11.4 Summary

243
243
256
266
269

Appendix A Interfaces: description and definitions
A.1 Two-dimensional geometry
A.2 Three-dimensional geometry
A.3 Axisymmetric geometry
A.4 Differentiation and integration on surfaces

270
270
272
274
275

Appendix B Distributions concentrated on the interface
B.1 A simple example

279
281

Appendix C Cube-chopping algorithm
C.1 Two-dimensional problem
C.2 Three-dimensional problem

284
285
286

Appendix D The dynamics of liquid sheets: linearized theory
D.1 Flow configuration
D.2 Inviscid results
D.3 Viscous theory for the Kelvin–Helmholtz instability
References
Index

288
288
288
293
295
322

Preface
Progress is usually a sequence of events where advances in one field open up new
opportunities in another, which in turn makes it possible to push yet another field
forward, and so on. Thus, the development of fast and powerful computers has
led to the development of new numerical methods for direct numerical simulations (DNS) of multiphase flows that have produced detailed studies and improved
knowledge of multiphase flows. While the origin of DNS of multiphase flows goes
back to the beginning of computational fluid dynamics in the early sixties, it is
only in the last decade and a half that the field has taken off. We, the authors of
this book, have had the privilege of being among the pioneers in the development
of these methods and among the first researchers to apply DNS to study relatively
complex multiphase flows. We have also had the opportunity to follow the progress
of others closely, as participants in numerous meetings, as visitors to many laboratories, and as editors of scientific journals such as the Journal of Computational
Physics and the International Journal of Multiphase Flows. To us, the state of the
art can be summarized by two observations:
• Even though there are superficial differences between the various approaches
being pursued for DNS of multiphase flows, the similarities and commonalities
of the approaches are considerably greater than the differences.
• As methods become more sophisticated and the problems of interest become
more complex, the barrier that must be overcome by a new investigator wishing
to do DNS of multiphase flows keeps increasing.
This book is an attempt to address both issues.
The development of numerical methods for flows containing a sharp interface,
as fluids consisting of two or more immiscible components inherently do, is currently a “hot” topic and significant progress has been made by a number of groups.
Indeed, for a while there was hardly an issue of the Journal of Computational
Physics that did not contain one or more papers describing such methods. In the
present book we have elected to focus mostly on two specific classes of methods:
volume of fluid (VOF) and front-tracking methods. This choice reflects our own
background, as well as the fact that both types of method have been very successful
and are responsible for some of the most significant new insights into multiphase
flow dynamics that DNS has revealed. Furthermore, as emphasized by the first
bullet point, the similarities in the different approaches are sufficiently great that
ix

x

Preface

a reader of the present book would most likely find it relatively easy to switch to
other methods capable of capturing the interface, such as level set and phase field.
The goal of DNS of multiphase flows is the understanding of the behavior and
properties of such flows. We believe that while the development of numerical
methods is important, it is in their applications where the most significant rewards
are to be found. Thus, we include in the book several chapters where we describe
the use of DNS to understand specific systems and what has been learned up to now.
This is inherently a somewhat biased sample (since we elected to talk about studies
that we know well – our own!), but we feel that the importance of these chapters
goes beyond the specific topics treated. We furthermore firmly believe that the
methods that we describe here have now reached sufficient maturity so that they
can be used to probe the mysteries of a large number of complex flows. Therefore,
the application of existing methods to problems that they are suited for and the
development of new numerical methods for more complex flows, such as those
described in the final chapter, are among the most exciting immediate directions
for DNS of multiphase flows.
Our work has benefitted from the efforts of many colleagues and friends. First
and foremost we thank our students, postdoctoral researchers and visitors for the
many and significant contributions they have made to the work presented here. GT
would like to thank his students, Drs. D. Yu, M. Song, S.O. Unverdi, E. Ervin,
M.R. Nobari, C.H.H. Chang, Y.-J. Jan, S. Nas, M. Saeed, A. Esmaeeli, F. Tounsi,
D. Juric, N.C. Suresh, J. Han, J. Che, B. Bunner, N. Al-Rawahi, W. Tauber, M.
Stock, S. Biswas, and S. Thomas, as well as the following visitors and postdoctoral
researchers: S. Homma, J. Wells, A. Fernandez, and J. Lu. RS would like to thank
his students, Drs. E. Aulisa, L. Campioni, A. Cervone, and V. Marra, as well as
his collaborators S. Manservisi, P. Yecko, and G. Zanetti. SZ would like to thank
his students, Drs. B. Lafaurie, F.-X. Keller, J. Li, D. Gueyffier, S. Popinet, A.
Leboissetier, L. Duchemin, O. Devauchelle, A. Bagué, and G. Agbaglah, as well
as his collaborators, visitors, and postdoctoral researchers, G. Zanetti, A. Nadim,
J.-M. Fullana, C. Josserand, P. Yecko, M. and Y. Renardy, E. Lopez-Pages, T.
Boeck, P. Ray, D. Fuster, G. Tomar, and J. Hoepffner. SZ would also like to thank
his trusted friends and mentors, Y. Pomeau, D.H. Rothman, and E.A. Spiegel, for
their invaluable advice. We also thank Ms. Victoria Tsengué Ingoba for reading the
complete book twice and pointing out numerous typos and mistakes. Any errors,
omissions, and ambiguities are, of course, the fault of the authors alone.
And last, but certainly not least, we would like to thank our families for unerring
support, acceptance of long working hours, and tolerance of (what sometimes must
have seemed) obscure priorities.
Grétar Tryggvason
Ruben Scardovelli
Stéphane Zaleski

1
Introduction

Gas–liquid multiphase flows play an essential role in the workings of Nature and
the enterprises of mankind. Our everyday encounter with liquids is nearly always
at a free surface, such as when drinking, washing, rinsing, and cooking. Similarly,
such flows are in abundance in industrial applications: heat transfer by boiling is
the preferred mode in both conventional and nuclear power plants, and bubbledriven circulation systems are used in metal processing operations such as steel
making, ladle metallurgy, and the secondary refining of aluminum and copper. A
significant fraction of the energy needs of mankind is met by burning liquid fuel,
and a liquid must evaporate before it burns. In almost all cases the liquid is therefore atomized to generate a large number of small droplets and, hence, a large
surface area. Indeed, except for drag (including pressure drops in pipes) and mixing of gaseous fuels, we would not be far off to assert that nearly all industrial
applications of fluids involve a multiphase flow of one sort or another. Sometimes,
one of the phases is a solid, such as in slurries and fluidized beds, but in a large
number of applications one phase is a liquid and the other is a gas. Of natural gas–
liquid multiphase flows, rain is perhaps the experience that first comes to mind, but
bubbles and droplets play a major role in the exchange of heat and mass between
the oceans and the atmosphere and in volcanic explosions. Living organisms are
essentially large and complex multiphase systems.
Understanding the dynamics of gas–liquid multiphase flows is of critical engineering and scientific importance and the literature is extensive. From a mathematical point of view, multiphase flow problems are notoriously difficult and much
of what we know has been obtained by experimentation and scaling analysis. Not
only are the equations, governing the fluid flow in both phases, highly nonlinear,
but the position of the phase boundary must generally be found as a part of the solution. Exact analytical solutions, therefore, exist only for the simplest problems,
such as the steady-state motion of bubbles and droplets in Stokes flow, linear inviscid waves, and small oscillations of bubbles and droplets. Experimental studies of
1

2

Introduction

Fig. 1.1. A picture of many buoyant bubbles rising in an otherwise quiescent liquid pool.
The average bubble diameter is about 2.2 mm and the void fraction is approximately 0.75%.
From Bröder and Sommerfeld (2007). Reproduced with permission.

multiphase flows are not easy either. For many flows of practical interest the length
scales are small, the time scales are short, and optical access to much of the flow is
limited. The need for numerical solutions of the governing equations has, therefore,
been felt by the multiphase research community since the origin of computational
fluid dynamics, in the late fifties and early sixties. Although much has been accomplished, simulations of multiphase flows have remained far behind homogeneous
flows where direct numerical simulations (DNS) have become a standard tool in
turbulence research.
In this book we use DNS to mean simulations of unsteady flow containing a
non-trivial range of scales, where the governing equations are solved using sufficiently fine grids so that all continuum time- and length-scales are fully resolved.
We believe that this conforms reasonably well with commonly accepted usage,
although we recognize that there are exceptions. Some authors feel that DNS
refers exclusively to fully resolved simulations of turbulent flows, while others
seem to use DNS for any computation of fluid flow that does not include a turbulence model. Our definition falls somewhere in the middle. We also note that some
authors, especially in the field of atomization – which is of some importance in this
book – refer to unresolved simulations without a turbulence model as LES. We also

1.1 Examples of multiphase flows

3

Fig. 1.2. A photograph of an atomization experiment performed with coaxial water and
air jets reproduced from Villermaux et al. (2004). Reproduced with permission. Copyright
American Physical Society.

prefer to call such computations DNS, especially as a continuous effort is made in
such simulations to check the results as the grid is refined. While it is not surprising that DNS of multiphase flows lags behind homogeneous flows, considering the
added difficulty, the situation is certainly not due to lack of effort. However, in the
last decade and a half or so, these efforts have started to pay off and rather significant progress has been accomplished on many fronts. It is now possible to do DNS
for a large number of fairly complex systems and DNS are starting to yield information that are likely to be unobtainable in any other way. This book is an effort
to assess the state of the art, to review how we came to where we are, and to provide the foundation for further progress, involving even more complex multiphase
flows.
1.1 Examples of multiphase flows
Since this is a book about numerical simulations, it seems appropriate to start by
showing a few “real” systems. The following examples are picked somewhat randomly, but give some insight into the kind of systems that can be examined by
direct numerical simulations.
Bubbles are found in a large number of industrial applications. For example,
they carry vapor away from hot surfaces in boiling heat transfer, disperse gases
and provide stirring in various chemical processing systems, and also affect the
propagation of sound in the ocean. To design systems that involve bubbly flows it
is necessary to understand how the collective rise velocity of many bubbles depends

4

Introduction

Fig. 1.3. The splash generated when a droplet hits a free surface. From A. Davidhazy.
Rochester Institute of Technology. Reproduced with permission.

on the void fraction and the bubble size distribution, how bubbles disperse and how
they stir up the fluid. Figure 1.1 is a picture of air bubbles rising through water in
a small bubble column. The average bubble diameter is about 2.2 mm and the
void fraction is approximately 0.75%. At these parameters the bubbles rise with an
average velocity of roughly 0.27 m/s, but since the bubbles are not all of the same
size they will generally rise with different velocities.
To generate sprays for combustion, coating and painting, irrigation, humidification, and a large number of other applications, a liquid jet must be atomized.
Predicting the rate of atomization and the resulting droplet size distribution, as
well as droplet velocity, is critical to the successful design of such processes. In
Fig. 1.2, a liquid jet is ejected from a nozzle of diameter 8 mm with a velocity of
0.6 m/s. To accelerate its breakup, the jet is injected into a co-flowing air stream,
with a velocity of 35 m/s. Initially, the shear between the air and the liquid leads
to large axisymmetric waves, but as the waves move downstream the air pulls long
filaments from the crest of the wave. The filaments then break into droplets by a
capillary instability. See Marmottant and Villermaux (2001) and Villermaux et al.
(2004) for details.
Droplets impacting solid or liquid surfaces generally splash, often disrupting the
surface significantly. Rain droplets falling on the ground often result in soil erosion,

1.1 Examples of multiphase flows

5

Fig. 1.4. Massive cavitation near the maximum thickness of an airfoil. The flow is from
the left to right. In addition to volume change due to phase change, the compressibility of
the bubbles is often important (see Section 11.1.5). From Kermeen (1956). Reproduced
with permission.

for example. But droplet impact can also help to increase the heat transfer, such as
in quenching and spray cooling, and rain often greatly enhances the mixing at the
ocean surface. Figure 1.3 shows the splash created when a droplet of a diameter of
about 3 mm, released from nearly 0.5 m above the surface, impacts a liquid layer a
little over a droplet-diameter deep. The impact of the droplet creates a liquid crater
and a rim that often breaks into droplets. As the crater collapses, air bubbles are
sometimes trapped in the liquid.
While bubbles are often generated by air injection into a pool of liquid or are
formed by entrainment at a free surface, such as when waves break, they also frequently form when a liquid changes phase into vapor. Such a phase change is
often nucleated at a solid surface and can take place either by heating the liquid
above the saturation temperature, as in boiling, or by lowering the pressure below
the vapor pressure, as in cavitation. Figure 1.4 shows massive cavitation near the
maximum thickness of an airfoil submerged in water. The chord of the airfoil is
7.6 cm, the flow speed is 13.7 m/s from left to right, and the increase in the liquid
velocity as it passes over the leading edge of the airfoil leads to a drop in pressure
that is sufficiently large so that the liquid “boils.” As the vapor bubbles move into
regions of higher pressure at the back of the airfoil, they collapse. However, residual gases, dissolved in the liquid, diffuse into the bubbles during their existence,
leaving traces that are visible after the vapor has condensed.

6

Introduction

Fig. 1.5. Microstructure of an aluminum–silicon alloy. From D. Apelian, Worcester Polytechnic Institute. Reproduced with permission.

In many multiphase systems one phase is a solid. Suspensions of solid particles
in liquids or gases are common and the definition of multiphase flows is sometimes
extended to cover flows through or over complex stationary solids, such as packed
beds, porous media, forests, and cities. The main difference between gas–liquid
multiphase flows and solid–gas or solid–liquid multiphase flows is usually that the
interface maintains its shape in the latter cases, even though the location of the
solid may change. In some instances, however, that is not the case. Flexible solids
can change their shape in response to fluid flow, and during solidification or erosion the boundary can evolve, sometimes into shapes that are just as convoluted
as encountered for gas–liquid systems. When a metal alloy solidifies, the solute is
initially rejected by the solid phase. This leads to constitutional undercooling and
an instability of the solidification front. The solute-rich phase eventually solidifies,
but with a very different composition than the material that first became solid. The
size, shape, and composition of the resulting microstructures determine the properties of the material, and those are usually sensitively dependent on the various
process parameters. A representative micrograph of an Al–Si alloy prepared by
metallographic techniques and etching to reveal phase boundaries and interfaces is
shown in Fig. 1.5. The light gray phase is almost pure aluminum and solidifies
first, but constitutional undercooling leads to dendritic structures of a size on the
order of a few tens of micrometers.
Living systems provide an abundance of multiphase flow examples. Suspended
blood cells and aerosol in pulmonary flow are obvious examples at the “body”

1.2 Computational modeling

7

Fig. 1.6. A school of yellow-tailed goatfish (Mulloidichthys flavolineatus) near the Northwest Hawaiian Islands. Since self-propelled bodies develop a thrust-producing wake, their
collective dynamics is likely to differ significantly from rising or falling bodies. From the
NOAA Photo Library.

scale, as are the motion of organs and even complete individuals. But even more
complex systems, such as the motion of a flock of birds through air and a school
of fish through water, are also multiphase flows. Figure 1.6 show a large number
of yellow-tailed goatfish swimming together and coordinating their movement. An
understanding of the motion of both a single fish and the collective motion of a
large school may have implication for population control and harvesting, as well
as the construction of mechanical swimming and flying devices.
1.2 Computational modeling
Computations of multifluid (two different fluids) and multiphase (same fluid, different phases) flows are nearly as old as computations of constant-density flows.
As for such flows, a number of different approaches have been tried and a number
of simplifications used. In this section we will attempt to give a brief but comprehensive overview of the major efforts to simulate multi-fluid flows. We make no
attempt to cite every paper, but hope to mention all major developments.
1.2.1 Simple flows (Re = 0 and Re = ∞)
In the limit of either very large or very small viscosity (as measured by the Reynolds
number, see Section 2.2.6), it is sometimes possible to simplify considerably the

8

Introduction

flow description by either ignoring inertia completely (Stokes flow) or by ignoring
viscous effects completely (inviscid flow). For inviscid flows it is usually further
necessary to assume that the flow is irrotational, except at fluid interfaces. Most
success has been achieved for disperse flows of undeformable spheres where, in
both these limits, it is possible to reduce the governing equations to a system of coupled ordinary differential equations (ODEs) for the particle positions. For Stokes
flow the main developer was Brady and his collaborators (see Brady and Bossis
(1988) for a review of early work) who have investigated extensively the properties of suspensions of particles in shear flows, among other problems. For inviscid
flows, Sangani and Didwania (1993) and Smereka (1993) simulated the motion
of spherical bubbles in a periodic box and observed that the bubbles tended to
form horizontal clusters, particularly when the variance of the bubble velocity was
small.
For both Stokes flows and inviscid potential flows, problems with deformable
boundaries can be simulated with boundary integral techniques. One of the earliest attempts was due to Birkhoff (1954), where the evolution of the interface between a heavy fluid initially on top of a lighter one (the Rayleigh–Taylor instability)
was followed by a method tracking the interface between two inviscid and irrotational fluids. Both the method and the problem later became a staple of multiphase
flow simulations. A boundary integral method for water waves was presented by
Longuet-Higgins and Cokelet (1976) and used to examine breaking waves. This
paper had enormous influence and was followed by a large number of very successful extensions and applications, particularly for water waves (e.g. Vinje and
Brevig, 1981; Baker et al., 1982; Schultz et al., 1994). Other applications include
the evolution of the Rayleigh–Taylor instability (Baker et al., 1980), the growth and
collapse of cavitation bubbles (Blake and Gibson, 1981; Robinson et al., 2001), the
generation of bubbles and droplets due to the coalescence of bubbles with a free
surface (Oguz and Prosperetti, 1990; Boulton-Stone and Blake, 1993), the formation of bubbles and droplets from an orifice (Oguz and Prosperetti, 1993), and the
interactions of vortical flows with a free surface (Yu and Tryggvason, 1990), just
to name a few. All boundary integral (or boundary element, when the integration is
element based) methods for inviscid flows are based on following the evolution of
the strength of surface singularities in time by integrating a Bernoulli-type equation. The surface singularities give one velocity component and Green’s second
theorem yields the other, thus allowing the position of the surface to be advanced
in time. Different surface singularities allow for a large number of different methods (some that can only deal with a free surface and others that are suited for
two-fluid problems), and different implementations multiply the possibilities even
further. For an extensive discussion and recent progress, see Hou et al. (2001).
Although continuous improvements are being made and new applications continue

1.2 Computational modeling

9

Fig. 1.7. A Stokes flow simulation of the breakup of a droplet in a linear shear flow. The
barely visible line behind the numerical results is the outline of a drop traced from an
experimental photograph. Reprinted with permission from Cristini et al. (1998). Copyright
2005, American Institute of Physics.

to appear, two-dimensional boundary integral techniques for inviscid flows are by
now – more than 30 years after the publication of the paper by Longuet-Higgins
and Cokelet – a fairly mature technology. Fully three-dimensional computations
are, however, still rare. Chahine and Duraiswami (1992) computed the interactions of a few inviscid cavitation bubbles and Xue et al. (2001) have simulated a
three-dimensional breaking wave. While the potential flow assumption has led to
many spectacular successes, particularly for short-time transient flows, its inherent
limitations are many. The lack of a small-scale dissipative mechanism makes those
models susceptible to singularity formation and the absence of dissipation usually makes them unsuitable for the predictions of the long-time evolution of any
system.
The key to the reformulation of inviscid interface problems with irrotational
flow in terms of a boundary integral is the linearity of the potential equation.
In the opposite limit, where inertia effects can be ignored and the flow is dominated by viscous dissipation, the Navier–Stokes equations become linear (the socalled Stokes flow limit) and it is also possible to recast the governing equations
as an integral equation on a moving surface. Boundary integral simulations of unsteady two-fluid Stokes problems originated with Youngren and Acrivos (1976)
and Rallison and Acrivos (1978), who simulated the deformation of a bubble and
a droplet, respectively, in an extensional flow. Subsequently, several authors have

10

Introduction

investigated a number of problems. Pozrikidis and collaborators have examined
several aspects of suspensions of droplets, starting with a study by
Zhou and Pozrikidis (1993) of the suspension of a few two-dimensional droplets in
a channel. Simulations of fully three-dimensional suspensions have been done by
Loewenberg and Hinch (1996) and Zinchenko and Davis (2000). The method has
been described in detail in the book by Pozrikidis (1992), and Pozrikidis (2001)
gives a very complete summary of the various applications. An example of a computation of the breakup of a very viscous droplet in a linear shear flow, using a
method that adaptively refines the surface grid as the droplet deforms, is shown in
Fig. 1.7.
In addition to inviscid flows and Stokes flows, boundary integral methods have
been used by a number of authors to examine two-dimensional, two-fluid flows in
Hele–Shaw cells. Although the flow is completely viscous, away from the interface
it is a potential flow. The interface can be represented by the singularities used for
inviscid flows (de Josselin de Jong, 1960), but the evolution equation for the singularity strength is different. This was used by Tryggvason and Aref (1983, 1985) to
examine the Saffman–Taylor instability, where an interface separating two fluids of
different viscosity deforms if the less viscous fluid is displacing the more viscous
one. They used a fixed grid to solve for the normal velocity component (instead of
Green’s theorem), but Green’s theorem was subsequently used by several authors
to develop boundary integral methods for interfaces in Hele–Shaw cells. See, for
example, DeGregoria and Schwartz (1985), Meiburg and Homsy (1988), and the
review by Hou et al. (2001).
Under the heading of simple flows we should also mention simulations of the
motion of solid particles, in the limit where the fluid motion can be neglected and
the dynamics is governed only by the inertia of the particles. Several authors have
followed the motion of a large number of particles that interact only when they
collide with each other. Here, it is also sufficient to solve a system of ODEs for
the particle motion. Simulations of this kind are usually called “granular dynamics.” For an early discussion, see Louge (1994); a more recent one can be found
in Pöschel and Schwage (2005), for example. While these methods have been
enormously successful in simulating certain types of solid–gas multiphase flows,
they are limited to a very small class of problems. One could, however, argue
that simulations of the motion of particles interacting through a potential, such as
simulations of the gravitational interactions of planets or galaxies and molecular
dynamics, also fall into this class. Discussing such methods and their applications
would enlarge the scope of the present work enormously, and so we will confine
our coverage by simply suggesting that the interested reader consults the appropriate references, such as Schlick (2002) for molecular simulations and Hockney and
Eastwood (1981) for astrophysical and other systems.

1.2 Computational modeling

11

Fig. 1.8. The beginning of computational studies of multiphase flows. The evolution of
the nonlinear Rayleigh–Taylor instability, computed using the two-fluid MAC method.
Reprinted with permission from Daly (1969b). Copyright 2005, American Institute of
Physics.

1.2.2 Finite Reynolds number flows
For intermediate Reynolds numbers it is necessary to solve the full Navier–Stokes
equations. Nearly 10 years after Birkhoff’s effort to simulate the inviscid Rayleigh–
Taylor problem by a boundary integral technique, the marker-and-cell (MAC)
method was developed at Los Alamos by Harlow and collaborators. In the MAC
method the fluid is identified by marker particles distributed throughout the fluid
region and the governing equations solved on a regular grid that covers both the
fluid-filled and the empty part of the domain. The method was introduced in Harlow and Welch (1965) and two sample computations of the so-called dam breaking
problem were shown in that first paper. Several papers quickly followed: Harlow
and Welch (1966) examined the Rayleigh–Taylor problem (Fig. 1.8) and Harlow
and Shannon (1967) studied the splash when a droplet hits a liquid surface. As
originally implemented, the MAC method assumed a free surface, so there was
only one fluid involved. This required boundary conditions to be applied at the
surface and the fluid in the rest of the domain to be completely passive. The Los
Alamos group realized, however, that the same methodology could be applied to
two-fluid problems. Daly (1969b) computed the evolution of the Rayleigh–Taylor
instability for finite density ratios and Daly and Pracht (1968) examined the initial motion of density currents. Surface tension was then added by Daly (1969a)
and the method again used to examine the Rayleigh–Taylor instability. The MAC
method quickly attracted a small group of followers that used it to study several
problems: Chan and Street (1970) applied it to free-surface waves, Foote (1973)

12

Introduction

and Foote (1975) simulated the oscillations of an axisymmetric droplet and the collision of a droplet with a rigid wall, respectively, and Chapman and Plesset (1972)
and Mitchell and Hammitt (1973) followed the collapse of a cavitation bubble.
While the Los Alamos group did a number of computations of various problems
in the sixties and early seventies and Harlow described the basic idea in a Scientific American article (Harlow and Fromm, 1965), the enormous potential of this
newfound tool did not, for the most part, capture the fancy of the fluid mechanics research community. Although the MAC method was designed specifically for
multifluid problems (hence the M for markers!), it was also the first method to
successfully solve the Navier–Stokes equation using the primitive variables (velocity and pressure). The staggered grid used was a novelty, and today it is a common practice to refer to any method using a projection-based time integration on a
staggered grid as a MAC method (see Chapter 3).
The next generation of methods for multifluid flows evolved gradually from the
MAC method. It was already clear in the Harlow and Welch (1965) paper that the
marker particles could cause inaccuracies, and of the many algorithmic ideas explored by the Los Alamos group, the replacement of the particles by a marker function soon became the most popular alternative. Thus, the volume-of-fluid (VOF)
method was born. VOF was first discussed in a refereed journal article by Hirt and
Nichols (1981), but the method originated earlier (DeBar, 1974; Noh and Woodward, 1976). The basic problem with advecting a marker function is the numerical
diffusion resulting from working with a cell-averaged marker function (see Chapter 4). To prevent the marker function from continuing to diffuse, the interface is
“reconstructed” in the VOF method in such a way that the marker does not start to
flow into a new cell until the current cell is full. The one-dimensional implementation of this idea is essentially trivial, and in the early implementation of VOF the
interface in each cell was simply assumed to be a vertical plane for advection in
the horizontal direction and a horizontal plane for advection in the vertical direction. This relatively crude reconstruction often led to large amount of “floatsam
and jetsam” (small unphysical droplets that break away from the interface) that
degraded the accuracy of the computation. To improve the representation, Youngs
(1982), Ashgriz and Poo (1991), and others introduced more complex reconstructions of the interface, representing it with a line (two dimensions) or a plane (three
dimensions) that could be oriented arbitrarily in such a way as to best fit the interface. This increased the complexity of the method considerably, but resulted in
greatly improved advection of the marker function. Even with higher order representation of the fluid interface in each cell, the accurate computation of surface
tension remained a major problem. In his simulations of surface tension effects
on the Rayleigh–Taylor instability, using the MAC method, Daly (1969b) introduced explicit surface markers for this purpose. However, the premise behind the

1.2 Computational modeling

13

Fig. 1.9. Computation of a splashing drop using an advanced VOF method. Reprinted
from Rieber and Frohn (1999) with permission from Elsevier.

development of the VOF method was to get away from using any kind of surface
marker so that the surface tension had to be obtained from the marker function instead. This was achieved by Brackbill et al. (1992), who showed that the curvature
(and hence surface tension) could be computed by taking the discrete divergence
of the marker function. A “conservative” version of this “continuum surface force”
method was developed by Lafaurie et al. (1994). The VOF method has been extended in various ways by a number of authors. In addition to better ways to
reconstruct the interface (Rider and Kothe, 1998; Scardovelli and Zaleski, 2000;
Aulisa et al., 2007) and compute the surface tension (Renardy and Renardy, 2002;
better solvers for the pressure equation have been introduced (see Rudman (1997),
for example). Other refinements include the use of sub-cells to keep the interface
as sharp as possible (Chen et al., 1997a). VOF methods are in widespread use today, and many commercial codes include VOF to track interfaces and free surfaces.
Figure 1.9 shows one example of a computation of the splash made when a liquid
droplet hits a free surface, done by a modern VOF method. We will discuss the use
of VOF extensively in later chapters.
The basic ideas behind the MAC and the VOF methods gave rise to several
new approaches in the early nineties. Unverdi and Tryggvason (1992) introduced
a front-tracking method for multifluid flows where the interface was marked by
connected marker points. The markers are used to advect the material properties
(such as density and viscosity) and to compute surface tension, but the rest of the
computations are done on a fixed grid as in the VOF method. Although using

14

Introduction

connected markers to update the material function was new, marker particles had
already been used by Daly (1969a), who used them to evaluate surface tension in
simulations with the MAC method, in the immersed-boundary method of Peskin
(1977) for one-dimensional elastic fibers in homogeneous viscous fluids and in the
vortex-in-cell method of Tryggvason and Aref (1983) for two-fluid interfaces in
a Hele–Shaw cell, for example. The front-tracking method of Unverdi and Tryggvason (1992) has been very successful for simulations of finite Reynolds number
flows of immiscible fluids and Tryggvason and collaborators have used it to explore
a large number of problems.
The early nineties also saw the introduction of the level-set, the CIP, and the
phase-field methods to track fluid interfaces on stationary grids. The level-set
method was introduced by Osher and Sethian (1988), but its first use to track fluid
interfaces appears to be in the work of Sussman et al. (1994) and Chang et al.
(1996), who used it to simulate the rise of bubbles and the fall of droplets in two
dimensions. An axisymmetric version was used subsequently by Sussman and
Smereka (1997) to examine the behavior of bubbles and droplets. Unlike the VOF
method, where a discontinuous marker function is advected with the flow, in the
level-set method a continuous level-set function is used. The interface is then identified with the zero contour of the level-set function. To reconstruct the material
properties of the flow (density and viscosity, for example) a marker function is constructed from the level-set function. The marker function is given a smooth transition zone from one fluid to the next, thus increasing the regularity of the interface
over the VOF method where the interface is confined to only one grid space. However, this mapping from the level-set function to the marker function requires the
level-set function to maintain the same shape near the interface and to deal with this
problem, Sussman et al. (1994) introduced a reinitialization procedure where the
level-set function is adjusted in such a way that its value is equal to the shortest distance to the interface at all times. This step was critical in making level-sets work
for fluid-dynamics simulations. Surface tension is found in the same way as in the
continuous surface force technique introduced for VOF methods by Brackbill et al.
(1992). The early implementation of the level-set method did not conserve mass
very well, and a number of improvements and extensions followed its original introduction. Sussman et al. (1998) and Sussman and Fatemi (1999) introduced ways
to improved mass conservation, Sussman et al. (1999) coupled level-set tracking
with adaptive grid refinement and a hybrid VOF/level-set method was developed
by Sussman and Puckett (2000), for example.
The constrained interpolated propagation (CIP) method introduced by Takewaki
et al. (1985) has been particularly popular with Japanese authors, who have applied
it to a wide variety of multiphase problems. In the CIP method, the transition from
one fluid to another is described by a cubic polynomial. Both the marker function

1.2 Computational modeling

15

and its derivative are then updated to advect the interface. In addition to simulating two-fluid problems, the method has been used for a number of more complex
applications, such as those involving floating solids; see Yabe et al. (2001).
In the phase-field method the governing equations are modified in such a way
that the dynamics of the smoothed region between the different fluids is described
in a thermodynamically consistent way. In actual implementations the thickness
of the transition is, however, much larger than it is in real systems and the net
effect of the modification is to provide an “antidiffusive” term that keeps the interface reasonably sharp. While superficially there are considerable similarities between phase-field and level-set methods, the fundamental ideas behind the methods
are very different. In the level-set method the smoothness of the phase boundary
is completely artificial and introduced for numerical reasons only. In phase-field
methods, on the other hand, the transition zone is real, although it is made much
thicker than it should be for numerical reasons. It is not clear, at the time of this
writing, whether keeping the correct thermodynamic conditions in an artificially
interface. The key drawback seems to be that since the propagation and properties of the interface depend sensitively on the dynamics in the transition zone, it
must be well resolved. For the motion of two immiscible fluids, that are well described by assuming a sharp interface, this adds a resolution requirement that is
more stringent than for other “one-fluid” methods. The phase-field approach was
originally introduced to model solidification (see Kobayashi (1992, 1993)) and has
found widespread use in such simulations. With the exception of the modeling of
solidification in the presence of flows (Beckermann et al., 1999; Tonhardt and Amberg, 1998), its use for fluid dynamic simulations is relatively limited (Jacqmin,
1999; Jamet et al., 2001). The main appeal of the phase-field methods appears to
be for problems where small-scale physics must be accounted for and it is difficult
to do so in the sharp interface limit.
In the “one-fluid” methods described above, where a single set of governing
equations is used to describe the fluid motion in both fluids, the fluid motion is
mostly computed on regular structured grids and the main difference between the
various methods is how a marker function is advected (and how surface tension is
found). The thickness of the interface varies from one cell in VOF methods to a few
cells in level-set and front-tracking methods, but once the marker function has been
found, the specific scheme for the interface advection is essentially irrelevant for
the rest of the computations. While these methods have been enormously successful, their accuracy is generally somewhat limited. There have, therefore, recently
been several attempts to generate methods that retain most of the advantages of
these methods but treat the interface as “fully sharp.” The origin of these attempts
can be traced to the work of Glimm and collaborators (Glimm et al., 1981; Glimm

16

Introduction

and McBryan, 1985; Chern et al., 1986), who used grids that were modified locally
near an interface in such a way that the interface coincided with a cell boundary,
and more recent “cut-cell” methods for the inclusion of complex bodies in simulations of inviscid flows (Quirk, 1994; Powell, 1998). In their modern incarnation,
sharp interface methods include the ghost fluid method, the immersed-interface
method and the method of Udaykumar et al. (2001). In the “ghost fluid” method
introduced by Fedkiw et al. (1999) the interface is marked by advecting a levelset function, but to find numerical approximations to derivatives near the interface,
where the finite difference stencil involves values from the other side of the interface, fictitious values are assigned to those grid points. The values are obtained by
extrapolation, and a few different possibilities for doing so are discussed by Glimm
et al. (2001), for example. The “immersed-interface” method of Lee and LeVeque
(2003), on the other hand, is based on modifying the numerical approximations
near the interface by explicitly incorporating the jump across the interface into the
finite difference equations. While this is easily done for relatively simple jump conditions, it becomes more involved for complex situations. Lee and LeVeque (2003)
thus found it necessary to limit their development to fluids with the same viscosity.
In the method of Udaykumar et al. (2001), complex solid boundaries are represented on a regular grid by merging cells near the interface and using polynomial
fitting to find field values at the interface. This method, which is related to the “cutcell” methods used for inviscid compressible flows (Powell, 1998), has so far only
been implemented for solids and fluids, including solidification (Yang and Udaykumar, 2005), but there seems to be no reason why the method cannot be used for multifluid problems. For an extension to three dimensions, see Marella et al. (2005).
While the original “one-fluid” methods require essentially no modification of
the flow solver near the interface (except allowing for variable density and viscosity), the sharp interface methods all require localized modifications of the basic
scheme. This results in considerably more complex numerical schemes, but is also
likely to improve the accuracy. That may be important for extreme values of the
governing parameters, such as large differences between the material properties of
the different fluids and low viscosities. The sharp interface approach may also be
required for flows with very complex interface physics. However, methods based
on a straightforward implementation of the “one-fluid” formulation of the governing equations, coupled with advanced schemes to advect the interface (or marker
function), have already demonstrated their usefulness for a large range of problems, and it is likely that their simplicity will ensure that they will continue to be
widely used.
In addition to the development of more accurate implementations of the “onefluid” approach, many investigators have pursued extension of the basic schemes
to problems that are more complex than the flow of two immiscible liquids. More

1.2 Computational modeling

17

Fig. 1.10. The interaction of two falling spheres. The spheres are shown at four different
times, going from left to right. Reprinted from Hu et al. (2001), with permission from
Elsevier.

complex physics has been incorporated to simulate contaminated interfaces, mass
transfer and chemical reactions, electrorheological effects, boiling, solidification,
and the interaction of solid bodies with a free surface or a fluid interface. We
will briefly review such advanced applications at the very end of the book, in
Chapter 11.
While methods based on the “one-fluid” approach were being developed, other
techniques were also explored. Hirt et al. (1970) describe one of the earliest use of
structured, boundary-fitted Lagrangian grids. In this approach a logically rectangular structured grid is used, but the grid points move with the fluid velocity, thus deforming the grid. This approach is particularly well suited when the interface topology is relatively simple and no unexpected interface configurations develop. In a
related approach, a grid line is aligned with the fluid interface, but the grid away
from the interface is generated using standard grid-generation techniques, such as
conformal mapping, or other more advanced elliptic grid-generation schemes. The
method was used by Ryskin and Leal (1984) to compute the steady rise of buoyant,
deformable, axisymmetric bubbles. They assumed that the fluid inside the bubble
could be neglected, but Dandy and Leal (1989) and Kang and Leal (1987) extended the method to two-fluid problems and unsteady flows. Several authors have
used this approach to examine relatively simple problems, such as the steady-state
motion of single particles or moderate deformation of free surfaces. Fully threedimensional simulations are relatively rare (see, though, Takagi et al. (1997)), and
it is probably fair to say that it is unlikely that this approach will be the method of
choice for very complex problems such as the three-dimensional unsteady motion
of several particles.

18

Introduction

A much more general approach to continuously represent a fluid interface by a
grid line is to use fully unstructured grids. This allows grid points to be inserted
and deleted as needed and distorted grid cells to be reshaped. While the grid was
moved with the fluid velocity in some of the early applications of this method, the
more modern approach is either to move only the interface points or to move the
interior nodes with a velocity different from the fluid velocity, in such a way that
the grid distortion is reduced but adequate resolution is still maintained. A large
number of methods have been developed that fall into this general category, but we
will only reference a few examples. Oran and Boris (1987) simulated the breakup
of a two-dimensional drop; Shopov et al. (1990) examined the initial deformation
of a buoyant bubble; Feng et al. (1994, 1995) and Hu (1996) computed the unsteady
two-dimensional motion of several particles and Fukai et al. (1995) followed the
collision of a single axisymmetric droplet with a wall. Although this appears to
be a fairly complex approach, Johnson and Tezduyar (1997) and Hu et al. (2001)
have produced very impressive results for the three-dimensional unsteady motion
of many spherical particles. Figure 1.10 shows an example of a simulation done
using the arbitrary-Lagrangian–Eulerian method of Hu et al. (2001). Here, two
solid spheres are initially falling in-line (left frame). Since the trailing sphere is
sheltered from the flow by the leading one, it catches up and “kisses” the leading
one. The in-line configuration is unstable and the spheres “tumble” (two middle
frames). After tumbling, the spheres drift apart (right frame).
The most recent addition to the collection of methods to simulate finite Reynolds
number multiphase flows is the lattice–Boltzmann method (LBM). It is now clear
that LBM can be used to obtain results of accuracy comparable to more conventional methods. It is still not clear, however, whether the LBM is significantly
faster or simpler than other methods (as sometimes claimed), but most likely these
methods are here to stay. For a discussion see, for example, Shan and Chen (1993)
and Sankaranarayanan et al. (2002). A comparison of results obtained by the LBM
method and the front-tracking method of Unverdi and Tryggvason (1992) can be
found in Sankaranarayanan et al. (2003). We will not discuss LBM in this book,
but refer the reader to Rothman and Zaleski (1997) and Chapter 6 in Prosperetti
and Tryggvason (2007).

Direct numerical simulations of multiphase flows have come a long way in the last
decade and a half or so. It is now possible to simulate accurately the evolution of
disperse flows of several hundred bubbles, droplets, and particles for sufficiently
long times so that reliable values can be obtained for various statistical quantities. Similarly, major progress has been achieved in the development of methods

19

for more complex flows, including those where a liquid solidifies or evaporates.
Simulations of large systems undergoing boiling and solidification are therefore
within reach.
Much remains to be done, however, and it is probably fair to say that the use
of direct numerical simulations of multiphase flows for research and design is still
in the embryonic state. The possibility of computing the evolution of complex
multiphase flows – such as churn-turbulent bubbly flow undergoing boiling, or the
breakup of a jet into evaporating droplets – will transform our understanding of
flows of enormous economic significance. Currently, control of most multiphase
flow processes is fairly rudimentary and almost exclusively based on intuition and
empirical observations. Industries that deal primarily with multiphase flows are,
however, multibillion dollar operations, and the savings realized if atomizers for
spray generation, bubble injectors in bubble columns, and inserts into pipes to
break up droplets, just to name a few examples, could be improved by just a little
bit would add up to a substantial amount of money. Reliable predictions would
also reduce the design cost significantly for situations such as space vehicles and
habitats where experimental investigations are expensive. And, as the possibilities
of manipulating flows at the very smallest scales by either stationary or free flowing microelectromechanical devices become more realistic, the need to predict the
effect of such manipulations becomes critical.
While speculating about the long-term impact of any new technology is a dangerous thing – and we will simply state that the impact of direct numerical simulations of multiphase flows will without doubt be significant – it is easier to predict
the near future. Apart from the obvious prediction that computers will continue
to become faster and more available, we expect that the development of numerical methods will focus mainly on flows with complex physics. Although some
progress has already been achieved for flows with variable surface tension, flows
coupled to temperature and electric fields, and flows with phase change, simulations of such systems are still far from being commonplace. In addition to the need
to solve a large number of equations, coupled systems generally possess much
larger ranges of length and time scales than simple two-fluid systems. Thus, the incorporation of implicit time-integrators for stiff systems and adaptive gridding will
become even more important. It is also likely, as more and more complex problems
are dealt with, that the differences between direct numerical simulations – where
everything is resolved fully – and simulations where the smallest scales are modeled will become blurred. Simulations of atomization where the evolution of thin
films are computed by “subgrid” models and very small droplets are included as
point particles are relatively obvious examples of such simulations (for a discussion
of the point-particle approximation, see Chapter 9 in Prosperetti and Tryggvason
(2007), for example). Other examples include possible couplings of continuum

20

Introduction

approaches such as those described in this book with microscopic simulations of
moving contact lines, kinetics effects at a solidifying interface, and reactions in
thin flames. Simulations of non-Newtonian fluids, where the microstructure has
to be modeled in such a way that the molecular structure is accounted for in some
way, also fall under this category.
In addition to the development of more powerful numerical methods, it is increasingly critical to deal with the “human” aspect of large-scale numerical simulations. The physical problems that we must deal with and the computational
tools that are available are rapidly becoming very complex. The difficulty of developing fully parallelized software to solve the continuum equations (fluid flow,
mass and heat transfer, etc.), where three-dimensional interfaces must be handled
and the grids must be dynamically adapted, is putting such simulations beyond the
reach of a single graduate student. In the future these simulations may even be
beyond the capacity of small research groups. It is becoming very difficult for a
graduate student to learn everything that they need to know and make significant
new progress in 4 to 5 years. Lowering the “knowledge barrier” and ensuring that
new investigators can enter the field of direct numerical simulations of multiphase
flow may well become as important as improving the efficiency and accuracy of
the numerical methods. The present book is an attempt to ease the entry of new
researchers into this field.

2
Fluid mechanics with interfaces

The equations governing multiphase flows, where a sharp interface separates immiscible fluids or phases, are presented in this chapter. We first derive the equations
for flows without interfaces, in a relatively standard manner. Then we discuss the
mathematical representation of a moving interface and the appropriate jump conditions needed to couple the equations across the interfaces. Finally, we introduce
the so-called “one-fluid” approach, where the interface is introduced as a singular
distribution in equations written for the whole flow field. The “one-fluid” form of
the equations plays a fundamental rôle for the numerical methods discussed in the
rest of the book.

2.1 General principles
The derivation of the governing equations is based on three general principles:
the continuum hypothesis, the hypothesis of sharp interfaces, and the neglect of
intermolecular forces. The assumption that fluids can be treated as a continuum
is usually an excellent approximation. Real fluids are, of course, made of atoms
or molecules. To understand the continuum hypothesis, consider the density or
amount of mass per unit volume. If this amount were measured in a box of sufficiently small dimensions , it would be a wildly fluctuating quantity (see Batchelor
(1970), for a detailed discussion). However, as the box side  increases, the density
becomes ever smoother, until it is well approximated by a smooth function ρ . For
liquids in ambient conditions this happens for  above a few tens of nanometers
(1 nm = 10−9 m). In some cases, such as in dilute gases, the discrete nature of
matter may be felt over much larger length scales. For dilute gases, the average
distance between molecular collisions, or the mean free path mfp , is the important length scale. The gas obeys the Navier–Stokes equations for scales   mfp .
Molecular simulations, where the motion of many individual molecules is followed
for sufficiently long times so that meaningful averages can be computed, show that
21

22

Fluid mechanics with interfaces

the fluid behaves as a continuum for a surprisingly small number of molecules. Koplik et al. (1988) found, for example, that under realistic pressure and temperature
a few hundred molecules in a channel resulted in a Poiseuille flow that agreed with
the predictions of continuum theory.
Beyond the continuum hypothesis, for multiphase flows we shall make the assumption of sharp interfaces. Interfaces separate different fluids, such as air and
water, oil and vinegar, or any other pair of immiscible fluids and different thermodynamic phases, such as solid and liquid or vapor and liquid. The properties of the
fluids, including their equation of state, density, viscosity and heat conductivity,
generally change across the interface. The transition from one phase to another occurs on very small scales, as described above. For continuum scales we may safely
assume that interfaces have vanishing thickness.
We also impose certain restrictions on the type of forces that are taken into account. Long-range forces between fluid particles, such as electromagnetic forces
in charged fluids, shall not be considered. Intermolecular forces, such as van der
Waals forces that play an important rôle in interface physics, are modelled by retaining their most important effect: capillarity. This effect, also called surface
tension, amounts to a stress concentrated at the sharp interfaces.
The three assumptions above also reflect the fact that it would be nearly impossible, with the current state of the art, to describe complex droplet and bubble
interactions while keeping the microscopic physics. For instance, simulating physical phenomena from the nanometer to the centimeter scale would require 107 grid
points in every direction, an extravagant requirement for any type of computation, even with the use of cleverly employed adaptive mesh refinement, at least at
present.
Beyond the three assumptions above, we mostly deal with incompressible flows
in this book, although in the present chapter we derive the equations initially for
general flow situations.
2.2 Basic equations
Expressing the basic principles of conservation of mass, momentum, and energy
mathematically leads to the governing equations for fluid flow. In addition to the
general conservation principles, we also need constitutive assumptions about the
specific nature of each fluid. Here we will work only with Newtonian fluids.
2.2.1 Mass conservation
The principle of conservation of mass states that mass cannot be created or destroyed. Therefore, if we consider a volume V , fixed in space, then the mass inside
this volume can only change if mass flows in or out through its boundary S. The

2.2 Basic equations

23

Fig. 2.1. A stationary control volume V . The surface is denoted by S.

flow out of V , through a surface element ds, is ρ u · n ds where n is the outward
normal, ρ is the density, and u is the velocity. The notation is shown in Fig. 2.1.
Stated in integral form, the principle of mass conservation is
d
dt


V

ρ dv = −


S

ρ u · n ds.

(2.1)

Here, the left-hand side is the rate of change of mass in the volume V and the
right-hand side represents the net flow through its boundary S. Since the volume is
fixed in space, we can take the derivative inside the integral and, by applying the
divergence theorem to the integral of the fluxes through the boundary, we have

 
∂ρ
(2.2)
+ ∇ · (ρ u) dv = 0.
V ∂t
This relation must hold for any arbitrary volume, no matter how small, and that
can only be true if the quantity inside the square brackets is zero. The partial
differential equation expressing conservation of mass is therefore

∂ρ
+ ∇ · (ρ u) = 0.
∂t

(2.3)

By using the definition of the substantial derivative
D() ∂ ()
=
+ u · ∇(),
Dt
∂t

(2.4)

and expanding the divergence, ∇ · (ρ u) = u · ∇ρ + ρ ∇ · u, the continuity equation
can be rewritten in convective form as
Dρ
= −ρ ∇ · u,
Dt

(2.5)

24

Fluid mechanics with interfaces

emphasizing that the density of a material particle can only change if the fluid is
compressed or expanded (∇ · u = 0).
2.2.2 Momentum conservation
The equation of motion is derived by using the momentum-conservation principle,
stating that the rate of change of fluid momentum in the fixed volume V is the
difference in momentum flux across the boundary S plus the net forces acting on
the volume. Therefore,








d
ρ u dv = − ρ u(u · n) ds + f dv + n · T ds.
(2.6)
dt V
S
V
S
The first term on the right-hand side is the momentum flux through the boundary
of V and the next term is the total body force on V . Frequently, the force per unit
volume f is only the gravitational force, f = ρ g. The last term is the total surface
force. Here, the tensor T is a symmetric stress tensor constructed in such a way
that n · T ds is the force on a surface element ds with a normal n.
By the same argument as applied to the mass conservation equation, Equation
(2.6) must be valid at every point in the fluid, so that

∂ ρu
(2.7)
= −∇ · (ρ uu) + f + ∇ · T.
∂t
Here we denote with ab, whose i jth component is ai b j , the dyadic product of the
two vectors a and b; hence, uu has components ui u j . The nonlinear advection term
can be written as
∇ · (ρ uu) = ρ u · ∇u + u∇ · (ρ u),

(2.8)

and using the definition of the substantial derivative and the continuity equation we
can rewrite Equation (2.7) as
Du
= f + ∇ · T.
(2.9)
Dt
This is Cauchy’s equation of motion and is valid for any continuous medium. For
fluids like water, oil, and air (as well as many others that are generally referred to
as Newtonian fluids) the stress may be assumed to be a linear function of the rate
of strain

ρ

T = (−p + λ ∇ · u)I + 2µ S.

(2.10)

Here, I is the unit tensor, p the pressure, µ the viscosity, and S = 12 (∇u + ∇uT ) is
the rate-of-strain or deformation tensor whose components are
1  ∂ ui ∂ u j 
+
.
(2.11)
Si j =
2 ∂ x j ∂ xi

2.2 Basic equations

25

λ is the second coefficient of viscosity, and if Stokes’ hypothesis is assumed to
hold, then λ = −(2/3)µ .1 Substituting the expression for the stress tensor into
Cauchy’s equation of motion results in
Du
= f − ∇p + ∇(λ ∇ · u) + ∇ · (2µ S),
Dt
which is the Navier–Stokes equation for fluid flow.

ρ

(2.12)

2.2.3 Energy conservation
In integral form the conservation of energy principle, as applied to a control volume
V fixed in space, is



1 2
d
(2.13)
ρ e + u dv
dt V
2

 



1 2
= − ρ e + u u · n ds + u · f dv + n · (u · T) ds − q · n ds.
2
S
V
S
S
Here, u2 = u · u and e is the total internal energy per unit mass. The left-hand side
is the rate of change of the internal and kinetic energy, the first term on the righthand side is the flow of internal and kinetic energy across the boundary, the second
term represents the work done by body forces, the third term is the work done by
the stresses at the boundary (pressure and viscous shear), and q in the fourth term
is the heat-flux vector. This equation can be simplified by using the momentum
equation. Taking the dot product of the velocity with Equation (2.9) gives

∂ u2 /2
= −ρ u · ∇u2 /2 + u · f + u · (∇ · T)
(2.14)
∂t
for the mechanical energy. After using this equation to cancel terms in (2.13)
and applying the same arguments as before, we obtain the convective form of the
energy equation:
De
ρ
− T : ∇u + ∇ · q = 0.
(2.15)
Dt
Here we denote with A : B = ∑i ∑ j Ai j B ji the scalar product of the two tensors A
and B.
This equation needs to be supplemented by constitutive equations for the specific
fluids we are considering. We will assume that the flux of heat is proportional to
the gradient of the temperature T , so that
ρ

q = −k∇T.
1

(2.16)

In several texts the discussion is based on the bulk viscosity κ = 2/3 µ + λ . Stokes’ hypothesis is that the bulk
viscosity vanishes.

26

Fluid mechanics with interfaces

This is called Fourier’s law and k is the thermal conductivity. Using Fourier’s law,
assuming a Newtonian fluid, so that the stress tensor is given by Equation (2.10),
and that radiative heat transfer is negligible, then the energy equation takes the
form
De
ρ
+ p∇ · u = Φ + ∇ · k∇T,
(2.17)
Dt
where Φ = λ (∇ · u)2 + 2µ S : S is called the dissipation function. It can be shown
that Φ, which represents the rate at which work is converted into heat, is always
greater or equal to zero.
In general we also need an equation of state giving, say, pressure as a function
of density and internal energy
p = p(e, ρ ),

(2.18)

as well as equations for the transport coefficients µ , λ , and k as functions of the
state of the fluid.
The governing equations are summarized in Panels 2.1 to 2.3. The integral form
obtained by applying the conservation principles directly to a small control volume
is shown in Panel 2.1. While usually not very convenient for analytical work, the
integral form is the starting point for finite-volume numerical methods. In Panel
2.2 we show the differential form obtained directly by assuming that the integral
laws hold at a point. This form can be used for finite difference numerical methods
and usually leads to discretizations that are essentially identical to those obtained
by the finite-volume method. In Panel 2.3 we have regrouped the terms to obtain
the convective or the non-conservative form of the equations. This is the form
usually shown in textbooks and is often used as a starting point for finite difference
methods.

2.2.4 Incompressible flow
For an important class of flows the density of each fluid particle does not change
as it moves. This is generally the case when the maximum or characteristic flow
velocity U is much smaller than the velocity of sound cs , or equivalently when
the Mach number Ma = U/cs is much smaller than unity. The equation for the
evolution of the density is then
Dρ
∂ρ
=
+ u · ∇ρ = 0,
Dt
∂t

(2.19)

and the mass conservation equation, Equation (2.5), becomes
∇ · u = 0.

(2.20)

2.2 Basic equations
d
dt
d
dt
d
dt


V


V

ρ u dv =

ρ dv +


S



ρ u · n ds = 0,

 

f dv +
V

27

S


n · T − ρ u(u · n) ds,








1 2
1 2
ρ e + u dv = u · f dv + n · u · T − ρ e + u u − q ds.
2
2
V
V
S





Panel 2.1. The equations of fluid motion in integral form. The flux terms have been moved
under the same surface integral as the stress terms.

∂ρ
+ ∇ · (ρ u) = 0,
∂t
∂ ρu
= f + ∇ · (T − ρ uu),
∂t


 


∂
1
1
ρ e + u2 = −∇ · ρ e + u2 u − T · u + q + u · f.
∂t
2
2
Panel 2.2. The equations of fluid motion in conservative form.
Dρ
+ ρ ∇ · u = 0,
Dt

ρ

ρ

Du
= f + ∇ · T,
Dt

De
= T : ∇u − ∇ · q.
Dt

Panel 2.3. The equations of fluid motion in convective or non-conservative form.

Equation (2.20) states that the volume of any fluid element cannot be changed and
these flows are therefore referred to as incompressible flows. If we integrate this
equation over a finite volume V with boundary S and use the divergence theorem
(or expand the divergence in (2.2) and use Equation (2.3)), we find that the integral

28

Fluid mechanics with interfaces

form of the mass conservation equation for incompressible flows is


S

u · n ds = 0,

(2.21)

stating that inflow balances outflow.
Notice that there is no requirement that the density is the same everywhere, for
incompressible flows. The density of a material particle can vary from one particle
to the next one, but the density of each particle must stay constant. When the
density is not the same everywhere its value at any given point in space can change
with time as material particles of different density are advected with the flow. In
this case the density field must be updated using Equation (2.19). If the density is
constant everywhere, then this is, of course, not necessary.
The pressure plays a special rôle for incompressible flows. Instead of being a
thermodynamic function of (say) density and temperature, it is determined solely
by the velocity field and will take on whatever value is necessary to make the flow
divergence free. It is sometimes convenient – and we will use this extensively –
to think of the pressure as projecting the velocity field into the space of incompressible functions. That is, we imagine the velocity first being predicted by (2.12)
without the pressure, then we find the pressure necessary to enforce incompressibility and correct the velocity field. Notice that for incompressible flow it is not
necessary to solve the energy equation to find the velocity and the pressure, unless the material properties are functions of the temperature. Thus, the flow field
is found by solving the momentum equation, Equation (2.12), along with the incompressibility condition, Equation (2.20) or (2.21). The governing equations for
incompressible, Newtonian flows in convective form are summarized in Panel 2.4.
The total derivative in the momentum equation in Panel 2.4 can be written in several different ways, all of which are equivalent analytically but that generally lead
to slightly different numerical approximations. The most common ones are

∂u
∂u 1
Du ∂ u
=
+ ∇ · (uu) =
+ u · ∇u =
+ ∇(u · u) − u × (∇ × u). (2.22)
Dt
∂t
∂t
∂t 2
We will usually work with the second form, but we note that the third one is very
commonly used as well. For compressible flows the fully conservative form of
the momentum equations is usually used, but as discussed at the end of Chapter 3,
doing so for incompressible flow with a density that changes abruptly can lead to
certain numerical difficulties.
The special case of incompressible fluids with constant density and viscosity is
of considerable importance. By writing the deformation tensor in component form
it is easily shown that ∇ · (∇u + ∇uT ) = ∇2 u and the momentum equation becomes
∂u
∇p f
+ + ν ∇2 u,
+ u · ∇u = −
∂t
ρ
ρ

(2.23)

2.2 Basic equations

29

∇ · u = 0,

ρ

Du
= −∇p + f + ∇ · µ (∇u + ∇uT ).
Dt

Panel 2.4. The equations of fluid motion for incompressible Newtonian flow in convective
form.

where ν = µ /ρ is the kinematic viscosity.
2.2.5 Boundary conditions
One of the major difficulties in numerical simulations of fluid flows is the correct
implementation of the boundary conditions. In principle the conditions at boundaries are well defined. For viscous, incompressible fluids we require the fluid to
stick to the wall so the fluid velocity there is equal to the wall velocity:
u = Uwall .

(2.24)

This equation includes both the normal and the tangential components of the velocity. For inviscid flows, where viscous stresses are absent and the fluid can slip
freely at the wall, only the normal velocity is equal to that of the wall. In some
cases, for instance for flow in very small channels, it is useful to introduce a slip
boundary condition for the tangential velocity component ut :
ut −Uwall = β

∂ ut
.
∂n

(2.25)

Here, β is a slip coefficient, with the dimensions of a (small) length, and ∂ /∂ n is
the derivative in the direction normal to the wall. The slip length is the distance at
which the velocity would vanish if extrapolated inside the wall (Fig. 2.2).
Frequently, we are interested in simulating a fluid domain with in- and out-flow
boundaries, or we are interested in simulating only a part of a larger fluid-filled domain. In those cases it is necessary to specify in- and out-flow boundary conditions.
For inflow boundaries the velocity is usually given. Realistic outflow boundaries,
on the other hand, pose a challenge that unfortunately does not have a simple solution. Similar difficulties are faced by the experimentalist who must carefully design
their wind tunnel to provide uniform inlet velocity and outlet conditions that have
minimal influence on the upstream flow. While it is probably easier to implement
uniform inlet conditions computationally than in experiments, the outlet flow is
more problematic, since the solution often available to the experimentalist, simply
making the wind tunnel long enough, usually requires an extensive number of grid

30

Fluid mechanics with interfaces
n

b

solid

Fig. 2.2. In some models the fluid is allowed to slip on the solid surface. The velocity
vanishes at some distance β inside the solid, called the slip length. This model is useful
when dealing with triple lines (contact lines) involving two fluids and a solid.

points. The goal is generally to truncate the domain using boundary conditions
that make it “appear” longer. Conditions on the rate of change of the velocity, such
as imposing a zero gradient on the velocity field, usually can be used if the outflow boundary is sufficiently far from the region of interest. Once the boundary
conditions for the velocity have been specified, conditions for the pressure can be
derived at both in- and out-flow boundaries.
In many cases computations are done in domains that are periodic in one or
more coordinate directions. If we assume that the x direction is periodic, then
the boundary conditions for the velocity are found from u(x, y, z) = u(x + L, y, z),
where L is the length of the period. Boundary conditions for pressure and other
variables are found in the same way, and similar formulae apply when the domain
is periodic in y and/or z. For theoretical work, periodic boundaries are often very
attractive, since it is not necessary to deal with the often difficult task of specifying
appropriate in- and out-flow boundary conditions. The only new consideration is
that if the domain is periodic in the direction of gravity, then it is necessary to
add a body force equal to ρav g, where ρav is the average density of the whole
domain, to prevent a uniform acceleration of the fluid. This force gives rise to a
hydrostatic pressure gradient, but since it results in the conservation of momentum
in the domain – not the average velocity – it does not correspond exactly to a
container with a rigid bottom.
2.3 Interfaces: description and definitions
Following the motion of a deformable interface separating different fluids or phases
is at the center of accurate predictions of multiphase flows. Describing the interface
location and how it moves can be accomplished in several ways. The simplest case

2.3 Interfaces: description and definitions

31

y

h(x,t)
x

Fig. 2.3. A simple way to parameterize an interface is through the equation y = h(x). However, overhangs in the curve may make the height function h multivalued. Other definitions
then become necessary.

is when the location of the interface can be described by a single-valued function
of one (in two dimensions) or two (in three dimensions) coordinates. For a thin
liquid layer, we can, for example, write
y = h(x) (two dimensions);

z = h(x, y) (three dimensions).

(2.26)

This situation is sketched in Fig. 2.3, which also shows the limitation of this approach: if part of the interface “overhangs,” h becomes a multivalued function.
While there are many situations, particularly for thin films, where Equation (2.26)
is the most convenient description, we generally need a more sophisticated way to
describe the interface.
To handle interfaces of arbitrary shape we can parameterize the interface by
introducing a coordinate u in two dimensions, such that the location of the interface
is given by


x(u) = x(u), y(u) .
(2.27)
Figure 2.4 shows a closed contour separating fluid “1” from fluid “2,” described
in this way. For a curve in two dimensions, described by Equation (2.27), dx/du
yields a vector tangent to the curve. Although u can be any parameterization of
the curve, the simplest case is when u is taken to be the arc length s, where ds2 =
dx2 + dy2 . The tangent vector t = dx/ds is then the unit tangent and |dx/ds| = 1.
The unit normal n is perpendicular to the curve, such that t · n = 0. The orientation
of the normal is arbitrary, and one should take the most convenient one depending
on the problem. For a closed curve it is customary to let the normal point outwards.

32

Fluid mechanics with interfaces
Fluid 2
t

Fluid 1

y
u

n

x(u)=(x(u),y(u))

x

Fig. 2.4. A parameterization of the interface. The coordinate u follows the interface.

However, later on we will use the convention that the normal points from one of
the fluids to the other. For instance, for an air–water mixture, water may be taken
as the reference fluid, air being the other fluid. The normal then points from water
to air. We usually label the reference fluid as fluid 1 and the other as fluid 2, and
the normal points from 1 to 2.
Consider the case when the normal points outwards and the parameterization follows the curve in the trigonometric (anticlockwise) direction as in Fig. 2.4. Then,
if the tangent vector is given by t = (t1 ,t2 ), the normal is given by n = (t2 , −t1 ).
The length of the unit tangent is by definition constant; so, as we move along the
curve, the change of the unit tangent vector must be in the normal direction. The
magnitude of the change of the tangent vector, with respect to the arc length s, is
the curvature of the interface. We therefore have
dt
= κ n.
ds

(2.28)

dn
= −κ t.
ds

(2.29)

Similarly, it can be shown that

The two previous equations are known as an instance of the Frenet–Serret formulae. The sign of the curvature may be a vexing question. It depends on one’s choice
of normal orientation, as shown on Fig. 2.5. With our previous conventions, for a
closed circle of radius R, the curvature is negative and equal to −1/R.
By taking the dot product of these equations with the normal and the tangent,
respectively, the curvature is found to be

κ = n·

dt
dn
= −t · .
ds
ds

(2.30)

2.3 Interfaces: description and definitions

33

n
n
k>0

k<0

Fig. 2.5. The sign of the curvature depends on one’s choice for the orientation of the normal. If the curve folds towards the normal the curvature is positive.
xu

n
dA

xv
t
p

u = const

v=const
x(u,v)=(x(u,v), y(u,v), z(u,v)

Fig. 2.6. A small surface element. The surface is parameterized by u and v. p is a vector
perpendicular to the edge of the element but tangent to the surface.

In Appendix A we examine the properties of a two-dimensional curve in more
detail and work out the above expressions in component form.
For an interface in three-dimensional space we need two independent coordinates, u and v, and the interface is described by


x(u, v) = x(u, v), y(u, v), z(u, v) ,
(2.31)
as shown in Fig. 2.6. The geometry of three-dimensional surfaces is an elaborate
subject and we will only need a few specific results. To avoid cluttering the discussion unnecessarily we have gathered some of the more pertinent elements of the
theory in Appendix A and here we simply quote the results that we need.
In addition to knowing the location of the interface, we frequently need to compute the normal to the interface, tangent vectors, and the mean curvature. Differentiating x(u, v) with respect to the surface coordinates yields two vectors tangent
to the surface. We will use the shorthand
xu =

∂x
∂u

and

xv =

∂x
.
∂v

(2.32)

With xu and xv both tangent to the interface, we can find the unit normal by taking

34

Fluid mechanics with interfaces

the cross product of the tangent vectors and dividing by the length of the product:
n=

xu × xv
.
|xu × xv |

(2.33)

In Appendix A it is shown that the mean curvature can be found by taking the
surface divergence of the normal vector:

κ = −∇S · n.

(2.34)

Here, ∇S denotes the surface divergence, defined in Appendix A. However, as also
shown in the appendix, if we define a normal field that extends off the interface,
then the curvature can be found by taking the usual divergence of this extended
normal field:

κ = −∇ · n.

(2.35)

We will also use the fact that the curvature can be found as the integral, along the
edges of an infinitesimal surface element of area δ A, of the interface tangent vector
p perpendicular to the edge of the surface element (see Fig. 2.6 and a derivation in
Appendix A):
1
δ A→0 δ A

κ n = lim



p dl.

(2.36)

Here, the vector p = t × n is tangent to the surface and perpendicular to both the
normal vector n and the vector t tangent to the edge of the surface element. Thus,
p “pulls” on the boundary of the surface element.
Instead of identifying the interface by explicitly specifying the location of every
point on the interface, the interface can also be given by a marker function defined
in the whole domain. Such marker functions can take many forms. We can, for
example, use the characteristic function, defined as a discontinuous function by
H(x) =

1 if inside a closed interface;
0 if outside a closed interface.

(2.37)

Alternatively, we may decide that H = 1 for fluid “1” and H = 0 for fluid “2.” With
H given, the interface is identified with the sharp change from one value to the
other (Fig. 2.7). H can be constructed in different ways, but for our purpose it is
convenient to express it in terms of an integral over the product of one-dimensional
δ -functions. For a two-dimensional domain:


H(x, y) =
V

δ (x − x )δ (y − y ) dv .

(2.38)

2.3 Interfaces: description and definitions

35

Fig. 2.7. A step function identifying the region occupied by a specific fluid.

Here, the integration is over the region bounded by the contour S and dv = dx dy
(Fig. 2.7). Obviously, H = 1 if the point (x, y) is inside the contour, where it will
coincide with (x , y ) during the integration, and zero if it is outside and will never
be equal to (x , y ) .
The properties of the step function are discussed in more detail in Appendix B,
where we show that the gradient of H is given by
∇H(x, y) = −


S

δ (x − x )δ (y − y )n ds = −δ (n)n,

(2.39)

where n is the coordinate normal to the interface in a local coordinate system
aligned with the interface. We can also define a surface distribution δS (x) = δ (n),
and write
∇H = −δS n.

(2.40)

In addition, the interface can be described by a smooth function F, if the interface is identified with a particular value of the function, say F = 0 (Fig. 2.8).
Obviously, F < 0 on one side of the interface and F > 0 on the other. One usually
takes the reference fluid (fluid “1”) to be in the F > 0 region. Taking F > 0 inside
the closed surface, the normal points outwards and is found by
n=−

∇F
|∇F|

(2.41)

and the curvature is given by

κ = −∇ · n = ∇ ·

∇F
.
|∇F|

(2.42)

36

Fluid mechanics with interfaces

Fig. 2.8. The identification of an interface by a level-set function.

The representation of the interface as a contour with a specific value is used in
level-set methods.
The motion of the interface S is determined by the normal velocity V (u, v,t), for
each point on S. This velocity may be that of the fluid itself, or may be different
due to, for example, evaporation or condensation. In parametric form:

∂
x(u, v,t).
(2.43)
∂t
To update the location of S, we note that there is no need to specify a tangential velocity, although from a physical point of view the interface is made of fluid particles
that travel at some tangential as well as normal velocity. The fact that the motion of
the interface is determined only by the normal velocity at the interface can be seen
by considering the interface as a level-set function F. Since the function moves
with the fluid velocity, its motion is described by
V (u, v,t) = n ·

∂F
DF
=
+ u · ∇F = 0.
(2.44)
Dt
∂t
The normal to the interface is given by (2.41), so we can rewrite the above equation
as
∂F
∂F
− u · n|∇F| =
−V |∇F| = 0,
(2.45)
∂t
∂t
which involves only the normal velocity V = u · n.
2.4 Fluid mechanics with interfaces
Consider the incompressible flow of two immiscible fluids filling a given domain.
The domain may be decomposed into any number of subdomains filled with the

2.4 Fluid mechanics with interfaces
n

37

V
u2

dV
dS
u1

S

Fig. 2.9. A thin control volume δ V with boundary δ S including a portion of the interface
S. The thickness of the control volume is taken to be zero, so no accumulation takes place.

individual phases. In any subdomain the usual Navier–Stokes equations hold, and
here we derive the interface conditions that allow us to couple the fluid motion in
the different fluids, across the interface. These conditions are, for the most part,
derived using mass and momentum conservation. In one case (for surface tension)
it is necessary to introduce additional physics.

2.4.1 Mass conservation and velocity conditions
To derive the jump conditions for the normal velocity, we apply the conservation
principle used in Section 2.2.1 to the control volume shown in Fig. 2.9. Since its
thickness is taken to be zero, there can be no accumulation of mass inside it, so the
mass flux into the control volume must be equal to the flow out. If the velocity on
one side of the interface is u1 and on the other it is u2 and the normal velocity of the
interface is V , then the inflow is u1 · n −V and the outflow is u2 · n −V . Denoting
the mass flow across the interface by ṁ, we have

ρ1 (u1 · n −V ) = ρ2 (u2 · n −V ) = ṁ .

(2.46)

This is simply the Rankine–Hugoniot condition. If there is no change of phase,
ṁ = 0. For incompressible flows this condition must hold for arbitrary density
ratios and we therefore have
V = u1 · n = u2 · n.

(2.47)

Mass conservation places no restriction on the tangential velocity components.
Indeed, inviscid fluids are usually assumed to slip at the interface. For very low
viscosity and rarefied gases it is possible that finite slip is the appropriate interface
condition, but for viscous fluid under normal operating conditions it is an experimentally observed fact – like the no-slip boundary conditions at solid walls – that
no slip takes place. It is also possible to argue that if the sharp interface is assumed
to come about by taking the limit of a finite interface zone as its thickness goes

38

Fluid mechanics with interfaces

to zero, then anything but no slip would result in infinitely high stresses if the interface zone consists of a fluid with a finite viscosity. Thus, if there is no phase
change and the fluids are incompressible, the interfacial condition for viscous fluid
is simply u1 = u2 , or
[u]S = 0,

(2.48)

where we have used the jump notation, i.e. the notation [x]S = x2 − x1 .
2.4.2 Surface tension
From a molecular point of view, surface tension arises because the interface is not
an optimal region thermodynamically: the molecules “prefer” to be at the gas or
the liquid density, which minimizes the free energy. The non-optimal conditions
near the interface result in an excess energy
dE σ = σ ds,

(2.49)

where ds is an infinitesimal interface or surface area and σ is a material property, usually referred to as the surface tension coefficient or simply surface tension.
From a mechanical point of view, however, surface tension is simply a force per
unit length acting perpendicularly on any line segment in the surface. If p is a
vector perpendicular to the line segment, the “pull” is simply σ p per unit length.
These points of view are equivalent: to stretch the surface one has to pull on it.
Suppose that we pull in the x direction and that the surface has extension L in the
other direction. Stretching it increases the area by ds = L dx, at the expense of a
work σ L dx = σ ds increasing as much the interfacial free energy.
By arguments that are exactly analogous to the introduction of the stress tensor
for a three-dimensional continuum (see Aris (1962), for example), it can be shown
that the force on the edge of a small surface element can be written as TσS · p, where
TσS is the surface tension tensor and p is a vector normal to the edge of the element,
in the tangent plane to the surface. Since
TσS · p = σ p = σ IS · p

(2.50)

we find that TσS = σ IS , where IS is the surface identity tensor.
The surface tension tensor was defined above in terms of the surface identity
tensor. Frequently, it is more convenient to work with the full three-dimensional
operator and use the fact that IS is the tangential projection of the three-dimensional
identity tensor I. Thus, IS = (I − nn) and
TσS = σ (I − nn).

(2.51)

Using that n = t1 × t2 , where t1 and t2 are orthonormal tangent vectors to the

2.4 Fluid mechanics with interfaces

39

surface, the surface tension tensor can also be written as
TσS = σ (t1 t1 + t2 t2 ).

(2.52)

The force on a surface element of area S, bounded by a contour C, is the integral of
the “pull” on its edges:

δ Fσ =


C

TσS · p dl =


S

∇S · TσS ds.

(2.53)

Taking the limit of Equation (2.53) as the surface area shrinks to a point, we define
the surface force per unit area as
fσ = ∇ · TσS = ∇ · σ IS = σ ∇S · IS + IS · ∇S σ .

(2.54)

The first term can be shown to be σ κ n and the second term is simply ∇S σ . Thus,
fσ = σ κ n + ∇S σ ,

(2.55)

where the last term, the surface gradient of σ , is obviously zero for constant surface
tension.

2.4.3 Momentum conservation with interfaces
Applying the conservation of momentum principle to the control volume in Fig.
2.9, moving with the interface S, we obtain
0=−


δS

ρ u(u · n −V )ds +


δS

n · T ds +


S

fσ ds.

(2.56)

Notice that the integration is around the edges of the control volume δ V for the
first two terms and along the interface for the last term. The first term on the righthand side is zero for incompressible flows in the absence of phase change, since the
control volume is moving with the fluid velocity. As the thickness of the control
volume approaches zero, the boundary coincides with the interface, and by integrating the stresses first on one side and then on the other, the second term yields
the jump in T across the interface. Since the results hold for any control volume
that includes the interface, they must be valid at any interface point, resulting in
− [T]S · n = σ κ n + ∇S σ .

(2.57)

This condition may be split into a normal and tangential stress condition, using
Equation (2.55) for fσ , yielding the boxed jump conditions in Panel 2.5.

40

Fluid mechanics with interfaces
[u]S = 0,
− [−p + 2µ n · S · n]S = σ κ ,
− 2µ t(k) · S · n

S

= t(k) · ∇S σ ,

∇ · u = 0,

ρ

Du
= −∇p + f + µ ∇2 u.
Dt

Panel 2.5. The jump conditions and the partial differential equations for a Newtonian,
incompressible fluid, without evaporation or condensation. The last two equations apply
in the bulk of each phase. The first three are jump conditions on the interface S. There,
the t(k) are two unit tangent vectors. In the two-dimensional flow case there is a single
tangent vector. Together with boundary conditions on walls or entry and exit conditions,
these equations form a complete set for flows with interfaces.

(−p + 2µ n · S · n)|S = −pfree + σ κ ,


2µ t(k) · S · n = t(k) · ∇S σ ,
S

∇ · u = 0,

ρ

Du
= −∇p + f + µ ∇2 u.
Dt

Panel 2.6. The boundary conditions and the partial differential equations for a Newtonian,
incompressible fluid with a free surface. In this example the pressure just outside the free
surface is set to pfree .

2.4.4 Free-surface flow
Free-surface flow is a limiting case of flows with interfaces, in which the treatment
of one of the fluids is simplified. For air–water flow, for instance, we can sometimes
assume that the pressure pfree in the air is a constant, or depends only on time, and
that the viscous stresses in the air are negligible. The jump conditions then become
boundary conditions for the liquid domain, as shown in Panel 2.6. In the case of
constant σ , the tangential stresses vanish and


t(k) · S · n = 0.
S

(2.58)

2.5 One-fluid formulation

41

This is a purely kinematic condition: it does not involve material properties at all.
For this reason it has interesting consequences. The vorticity vector is defined in
general as ∇ × u and its study reveals interesting aspects of the dynamics of the
flow. In two dimensions, we note ω = (∇ × u) · ez , its only non-zero component.
Then on the free surface

ω = 2κ q,

(2.59)

where q = (u · u)1/2 is the norm of the velocity, as shown in Batchelor (1970).
This implies that a large vorticity will be seen on the free surface in regions of
high velocity and high curvature. This is were the departure from the hypothesis
of irrotational flow will be the most pronounced. Moreover, the vorticity is also
seen in these regions for two-phase flows with a strong density and viscosity contrast, since these flows resemble free-surface flows. Numerical simulations show
characteristic patches of vorticity in the regions of high curvature.

2.5 Fluid mechanics with interfaces: the one-fluid formulation
In contrast to the approach described in Section 2.4, where we wrote the governing
equations separately for each phase and used jump conditions to couple the solutions at the fluid interface, it is possible to write one set of governing equations
for the whole flow domain occupied by the various phases, without resorting to
jump conditions. The various phases are treated as one fluid with variable material
properties that change abruptly at the phase boundary. To account for the “extra”
forces at the phase boundary, however, it is necessary to add singular terms (δ functions) to the equations. These singular terms are the counterpart of the jump
conditions of the preceding section, and it can be shown that both formulations are
equivalent.
This form of the equations is often referred to as the “one-fluid” approach. Since
the solution can change discontinuously across the interface, we must either interpret the governing equations in a weak sense where they are satisfied only in an
integral sense, or admit solutions that include generalized functions, such as delta
functions and step functions. The latter approach is taken here. The “one-fluid”
formulation is the starting point for several numerical methods based on the use of
fixed grids, as we will see in later chapters.
The derivation of the “one-fluid” equations is exactly the same as before (Section 2.2), except that we need to add the surface tension as a body force to the
momentum equation. The mass conservation equation for the total mass has no
source terms. For reacting flows, the conservation equation for each species may
have localized source terms, and for flows with phase change (Chapter 11), we
shall see that a localized volume source emerges naturally at the phase boundary.

42

Fluid mechanics with interfaces

For a control volume including an interface, the surface tension force is given by
fσ integrated over S, the part of the surface enclosed in the control volume. Just as
we use the divergence theorem to convert surface integrals to volume integrals, we
use Equation (B.9) to transform the surface integral




S

fσ ds =

V

fσ δS dv,

(2.60)

where δS = δS (x − xs ). Adding this force to the integral form of the momentum
equation, Equation (2.6), and applying the same arguments as before – that the integral can only be zero for all possible control volumes if the integrand is zero – we
obtain the “one-fluid” version of the Navier–Stokes equation for incompressible,
Newtonian flows with sharp interfaces:
∂u
ρ
(2.61)
+ ρ ∇ · uu = −∇p + f + ∇ · µ [∇u + (∇u)T ] + fσ δS .
∂t
We emphasize that this equation holds for the whole flow field, even if the density field ρ and the viscosity field µ change discontinuously. For constant surface
tension, the surface force is
fσ δS = σ κ nδS .

(2.62)

For completeness, we observe that the last term can be written in several different, but equivalent forms:
f σ δS =



fσ δ (x − xs ) dA

and

σ κ n δS =



σ κ nδ (x − xs ) dA,

(2.63)

where we have used that δ (x − xs ) = δ (x − xs )δ (y − ys )δ (z − zs ). Notice also that
it is possible to add an arbitrary constant to the surface tension term, replacing σ κ
by (σ κ + ∆po ). The constant ∆po is balanced exactly by a change in the pressure
jump across the interface, and since the flow on either side of the interface is driven
by the pressure gradient – not the absolute value of the pressure – adding the constant has no effect on the flow. Jan (1994) used this observation to subtract the
average pressure inside a bubble to reduce the pressure jump across the interface
in his computations of the motion of clean and contaminated axisymmetric and
two-dimensional bubbles.
2.6 Nondimensional numbers
The dynamics of multifluid and multiphase flows is governed by a variety of nondimensional numbers, depending on the exact situation driving the flow. Those include the Reynolds number, usually defined based on the properties of one of the
fluids, using an external length scale L and velocity U:
Re =

ρ LU
LU
=
.
µ
ν

(2.64)

2.6 Nondimensional numbers

43

Re represents
the ratio of inertial to viscous forces. If we construct a velocity scale

by (∆ρ /ρ )gL, then the square of the Reynolds number becomes the Galileo
number:
gρ ∆ρ L 3
.
(2.65)
N=
µ2
When this number is large, gravity is balanced by inertia and turbulent dissipation;
when it is small, gravity and viscous forces balance. The ratio of inertia to surface
tension gives the Weber number, defined by

ρ LU 2
,
σ

We =

(2.66)

and the ratio of viscous to capillary stresses yields the capillary n```
{"ft_lang_label":"__label__en","ft_lang_prob":0.91982484,"math_prob":0.9418609,"size":104554,"snap":"2020-10-2020-16","text_gpt3_token_len":23620,"char_repetition_ratio":0.16968915,"word_repetition_ratio":0.038786985,"special_character_ratio":0.22197142,"punctuation_ratio":0.10743468,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.9780824,"pos_list":[0],"im_url_duplicate_count":[null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2020-03-30T20:02:42Z\",\"WARC-Record-ID\":\"<urn:uuid:00de5136-d6d5-4e02-8d87-53b489702253>\",\"Content-Length\":\"162753\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:a577018a-be68-4d4d-992f-1c97a0a6b0c6>\",\"WARC-Concurrent-To\":\"<urn:uuid:fe7c4670-de2e-4b13-8cb0-4dfc8b1f0af3>\",\"WARC-IP-Address\":\"179.43.147.124\",\"WARC-Target-URI\":\"https://ru.b-ok.org/book/1256304/da2840\",\"WARC-Payload-Digest\":\"sha1:535GRA2V55BI5E275PHGSIBGVF2XX7NR\",\"WARC-Block-Digest\":\"sha1:DVIBSZCJ3H5H2U3PCFTWL2VBJ6GZYZFE\",\"WARC-Identified-Payload-Type\":\"text/html\",\"warc_filename\":\"/cc_download/warc_2020/CC-MAIN-2020-16/CC-MAIN-2020-16_segments_1585370497301.29_warc_CC-MAIN-20200330181842-20200330211842-00072.warc.gz\"}"}