DNS of multiphase flows often result in features much smaller than the dominant flow scales, consisting of very thin films, filaments and drops/bubbles. Here we describe the development of a general strategy to include multi-scale description of small-scale phenomenon in numerical simulations of the dynamics of multiphase systems. Two examples are described in some details. The first is the use of thin film modelling when fluid masses collide and form a thin film. The evolution of the
thickness of the film can be described using lubrication models that allow us to follow the evolution of films much thinner than the grid spacing used to resolve the flow outside the film. By using a thin film model to modify the boundary conditions at the wall, a coarse grid can be used to produce results nearly as accurate as a fully resolved solution on a much finer grid (Thomas et al, 2009). The other example is an attempt to capture thin mass transfer boundary layers in simulations of reactions in
bubble columns by using a boundary layer model to capture the thin mass boundary layer on the liquid site of the interface.

General Note:

The International Conference on Multiphase Flow (ICMF) first was held in Tsukuba, Japan in 1991 and the second ICMF took place in Kyoto, Japan in 1995. During this conference, it was decided to establish an International Governing Board which oversees the major aspects of the conference and makes decisions about future conference locations. Due to the great importance of the field, it was furthermore decided to hold the conference every three years successively in Asia including Australia, Europe including Africa, Russia and the Near East and America. Hence, ICMF 1998 was held in Lyon, France, ICMF 2001 in New Orleans, USA, ICMF 2004 in Yokohama, Japan, and ICMF 2007 in Leipzig, Germany. ICMF-2010 is devoted to all aspects of Multiphase Flow. Researchers from all over the world gathered in order to introduce their recent advances in the field and thereby promote the exchange of new ideas, results and techniques. The conference is a key event in Multiphase Flow and supports the advancement of science in this very important field. The major research topics relevant for the conference are as follows: Bio-Fluid Dynamics; Boiling; Bubbly Flows; Cavitation; Colloidal and Suspension Dynamics; Collision, Agglomeration and Breakup; Computational Techniques for Multiphase Flows; Droplet Flows; Environmental and Geophysical Flows; Experimental Methods for Multiphase Flows; Fluidized and Circulating Fluidized Beds; Fluid Structure Interactions; Granular Media; Industrial Applications; Instabilities; Interfacial Flows; Micro and Nano-Scale Multiphase Flows; Microgravity in Two-Phase Flow; Multiphase Flows with Heat and Mass Transfer; Non-Newtonian Multiphase Flows; Particle-Laden Flows; Particle, Bubble and Drop Dynamics; Reactive Multiphase Flows

DNS of multiphase flows often result in features much smaller than the dominant flow scales, consisting of very thin films,
filaments and drops/bubbles. Here we describe the development of a general strategy to include multi-scale description of
small-scale phenomenon in numerical simulations of the dynamics of multiphase systems. Two examples are described in
some details. The first is the use of thin film modelling when fluid masses collide and form a thin film. The evolution of the
thickness of the film can be described using lubrication models that allow us to follow the evolution of films much thinner than
the grid spacing used to resolve the flow outside the film. By using a thin film model to modify the boundary conditions at the
wall, a coarse grid can be used to produce results nearly as accurate as a fully resolved solution on a much finer grid (Thomas
et al, 2009). The other example is an attempt to capture thin mass transfer boundary layers in simulations of reactions in
bubble columns by using a boundary layer model to capture the thin mass boundary layer on the liquid site of the interface.

Introduction

In spite of the enormous information and understanding that
DNS are providing for relatively complex flows Prosperetti
& Tryggvason, 2007; Lu, Fernandez & Tryggvason, 2005;
van Sint Annaland, Deen & Kuipers, 2005; Hua, Stene, Jan
& Lin, 2008) real systems provide challenges that still limit
the range of situations that can be simulated, even when we
limit our studies to systems well described by continuum
theories. Generally the smallest length scale of the flow has
to be resolved by 0(10) grid points or so, and as the number
of grid points available increases, the range of scales that
can be resolved goes up. Current computer power makes it
possible to simulate two fluid systems in domains resolved
by several hundred grid points in each spatial direction (for
thousands of time steps) relatively routinely and simulations
of systems resolved by over 10003, or billion grid points, are
becoming possible. Simulations of this size do obviously
offer the opportunities to span scales whose ratios span over
two orders to magnitude. This is, however, often not enough.
Starting with simulations where what we might call the
"dominant small-scales" (see below) are fully resolved, it is
frequently found that multiphase flows also can generate
features much smaller than the dominant flow scales,
consisting of very thin films, filaments, and drops.
Frequently there is a clear separation of scales between
these "features," usually inertia effects are relatively small
for the local evolution, and in isolation these features are
often well described by analytical models. While these
features can, in principle, be captured by local, adaptive,
grid refinement, doing so increases the complexity of the
computations significantly and usually results in greatly

increased computational time.
The "natural" or a dominant small-scale in many multiphase
fluid-fluid systems is set by the balance of surface tension
and viscosity and inertia. For the breakup of a jet, for
example, this scale determines the average droplet size. In
most cases, this scale corresponds to roughly where the
appropriately chosen nondimensional numbers, such as
Weber, Capillary, Ohnsorge, and Reynolds, are 0(1) (and
the key word here is obviously "appropriately").
Unfortunately, however, this is often not the only length
scale that exists and frequently the flow also has much
smaller features. Consider, for example, the relatively tame
problem of the collision of two droplets that are a few
hundred micrometers in diameters. As the drops collide,
they deform and trap air in a thin film between them. The
air drains out of the film and if the drops stay in contact
long enough, the film ruptures. It is, however, generally
believed that the film thickness must get down to just about
few hundred Angstroms before it ruptures. 100 A = 10nm =
0.01 p m, so if the drops are 500 i m in diameter, the
range of scales, from the thickness of the film to the
diameter of the drop is 104. It is, at least in principle,
possible to capture this motion for very simple situations
(such as the head-on collision of two drops, where the
motion could be resolved by unevenly spaced grid points
using axisymmetry to reduce the computations effort).
However, for more complex systems of several drops
moving arbitrarily, such a simulation is essentially
impossible at the present time.
The possibility of representing small-scale multiphase flow
features in simplified form is perhaps best demonstrated by
the point particle approximations. For very small drops, that

are smaller than the smallest flow scales (and whose Stokes
number is smaller than unity) several investigators have
carried out simulations where the flow away from the
particles is fully resolved but the particles are approximated
as points. Although the flow is not resolved by the
computational grid, there are good reasons to believe that
the accuracy of these models can be very high and when the
Reynolds number of the particle is low enough, analytical
results for the forces between the particle and the fluid
eliminates any empirical adjustments. Another example of a
multiscale approach in multiphase flow simulations is the
near field correction used in simulations of multi particle
systems in the Stokes flow limit, where either a lubrication
approximation (Bossis & Brady, 1984, and others) or an
analytical solution (Ingber et al, 2008) is used to resolve the
flow between nearly touching solids.

Multiscale DNS-General strategies

Multiphase flows possess a multitude of scales and the
multiscale modeling discussed here can be put into a
broader context by considering the classification of E &
Enquist (21' 1). They divide multiscale problems into two
broad categories:
Type A Problems: Dealing with Isolated Defects
Type B Problems: Constitutive Modeling Based on
Microscopic Models
Figure 1 shows the classification schematically for bubbly
flow in an inclined channel. Problem B is the classical
averaging encountered in turbulence and multiphase flow
modelling where we seek a model that accurately captures
the aspects of the flow that are of relevance to us, without
having to compute the intricate details of the flow. Using
DNS data for multiphase flow modeling has already shown
that this is a powerful approach (Lu, Biswas & Tryggvason,
2006, for example). While much is yet to be done and
considerable progress continues to be made, problem B is
not the focus of the present discussion. Here we address
Problem A.
In direct numerical simulations of multiphase flows, all
scales are fully resolved. Most of the flow is resolved in
standard ways, where each term in the governing equations
is approximated using a finite difference or a finite volume
technique. To resolve certain small scale features, however,
such as thin films and threads, thin boundary/reaction layers,
and small drops and bubbles, we can take advantage of the
fact that these features often have a relatively simple
structure-since high surface tension imposes a simple
geometry and high viscosity leads to relatively simple
flow-and can therefore be described accurately using
analytical or quasi-analytical theories. Given a simple
description (which may have to be solved numerically) of a
feature not resolved on the fluid grid, the challenges are to
identify when and where the small scale description should
be used and how to pass information between the
numerically computed flow and the small-scale description.
Obviously we have to introduce new computational objects
to capture the small-scale features. In front tracking
methods (Unverdi & Tryggvason, 1992; Tryggvason et al,
2001), where connected marker points are used to represent
the interface between different fluids, those already exist
and the general data structure for other things. Information
about the resolved flow (such as strain-rate) is easily

Thin film f
model-A
Figure 1. A schematic of multiscale modeling of
bubbly flow in an inclined channel, showing where "A"
and "B" type problems arise.

interpolated from the fixed fluid grid. Passing information
back, from the small-scale description to the flow grid is
more complex, and dependent on exactly what is being
transferred. Some quantities, such as forces, can be directly
smoothed onto the fixed grid, but in other cases we may
need to transfer fluxes, rather then the quantity itself (for
heat and mass flow, for example). We note that for the most
part computational efficiency of how the models are
evolved and incorporated is a relatively small concern. In
front tracking simulations of multiphase flows the interfaces
are of on dimension lower than the flow field so any
operations performed on the interface generally take much
less time than the flow solver. For methods based on
following the fluid interface by advecting a marker function
directly, such as VOF and level set methods, it is of course
necessary to introduce additional computational elements to
follow the evolution of the small scale processes. Below we
discuss two examples of the use of an analytical model to
follow the evolution of a scale processes.

Thin films and topology changes

Multiphase flows often develop thin films when two fluid
masses collide and trap another fluid between then and
narrow threads when a fluid mass breaks up. Both processes
can lead to topology changes. The thin films can rupture
and the thin threads can snap. Thin threads appear to
be-by far-the more easier to deal with. There are good
reasons to believe that the Navier-Stokes equations describe
how their diameter becomes zero in a finite time and
simulations suggest that the overall results are relatively
insensitive to how they are treated. In computations we
often simply do nothing and are left with a very thin string
of points, representing a filament of essentially zero
thickness. For very accurate treatment of this filament, we
should keep track of how it breaks into small drops and here
we will develop models that do that, representing the
resulting drops as point particles.
The tracking of the draining of thin films and their rupture is
much more critical to the overall accuracy of the
calculations. For clean interfaces that are fully mobile, we
generally find that we capture well their evolution even if
they are poorly resolved (not surprising since to the rest of

Figure 2. Simulations of a drop falling onto a sloping
wall. The solution represented by the thick line is
computed on a grid fully resolving the film between the
drop and the wall.

the flow the film is simply a membrane), with the exception
of its thickness. Thickness is critical for predictions of
rupture and the main goal of tracking the film accurately is
to be able to predict when it can rupture. In the front
tracking method that we use, where the fluid interface is
explicitly marked by connected markers, the film never
ruptures by default but in simulations where the different
fluids are followed by the advection of a marker function,
such as in VOF or level set methods, the film would always
rupture once the film is thinner than a grid spacing or so.
Both methods can, of course, be modified to behave
differently and we have conducted a number of simulations
where we rupture films either when they are thin enough or
at given times. Not capturing rupture correctly can prevent
results computed using marker functions from converging
since finer grids postpone the rupture. The accurate
capturing of the rupture of a thin film requires three steps.
First of all, the film must be identified. This can be done in
many ways, although some ways are more efficient then
others. In the example shown below, this is accomplished in
a relatively simple way by checking the distance of front
points from the bottom wall. Second, the film must be
modeled. To do so we use thin film theories where we solve
equations for the thickness of the film, using the velocity
field computed on the fixed grid as input. Although our
ultimate plan is to be able to deal with films appearing
anywhere in the flow, we have started by looking at the

Figure 3. Drop sliding down a sloping wall computed
on a coarse grid (thin line) and on a coarse grid using a
simple multiscale model (dashed line) along with a fully
converged solution (thick line) at a later time than in
figure 2.

simpler case of films next to walls. Our work on modelling
such films has been described in detail Thomas, Esmaeeli &
Tryggvason (2009) so here we only give a brief review of
the main results.
In figure 2 we show four frames from simulations of a drop
falling onto a sloping wall and then sliding down under the
action of gravity. The computations have been done using
two different grid resolutions, both of which are sufficient
to resolve the motion in the absence of a wall. The thin line
shows results for a uniform grid where the film between the
drop and the wall is under resolved. The thick line shows
the results when the grid has been refined near the wall and
the film is fully resolved. Obviously there are significant
differences at the later times. In figure 3 we plot the drop
shape for the last time shown in figure 2 for both resolutions
as well as a simulation carried out on the coarser grid but
with a thin film model to describe the flow in the film
between a drop and the wall. The film model assumes a
liner velocity profile, resulting in a two-equation
(conservation of mass and momentum) model:
Sh 1 9
-+-- (hU,) = 0;
dt 2 dx
2 2h dp
-(hU,)+ (hUW) .
dt 3 dx p dx ),
Outside the film we apply the usual no-slip boundary
conditions when computing the flow field, but where the
film is, the shear stress, as found from the film model, is
used. Obviously the use of the model produces results that
are relatively close to the fully converged solution (Thomas,
Esmaeeli & Tryggvason, 2009). The film model used in the
simulations in figure 2 and 3 is fairly simple but
nevertheless improves the agreement with the fully resolved
solution significantly. More sophisticated models, similar to
those of Yoon et al (2007), Dai et al (2008), Baldessari &
Leal (2007), for example, should result in an even better
agreement. As of yet we have not incorporated a rupture
model into simulations with a film model.
Studies of thin films have a long history and model
equations have been developed for a wide range of
conditions. Such models have been coupled with flow
outside the film, such as by Davis, Schonberg, & Rallison
(1989) who examined the collision of drops in the zero
Reynolds number limit. In addition to the film between

colliding drops, we expect thin, relatively flat films also
between drops and bubbles sliding along a solid wall, and
when a liquid is extended in a strain field where
compressive strain is applied in one dimension and
stretching in the other two. While the incorporation of
models for thin films in simulations using the Navier-Stokes
equations is rare, it is not completely non-existing. Son &
Dhir (1998) in their level-set simulations of film boiling
developed a simple axisymmetric model for the evolution of
the liquid micro layer left behind on the wall as the vapor
bubble expands. The heat transfer in this layer is believed to
be important for the total evaporation rate. By modeling the
micro layer using the usual thin film approximations, Son &
Dhir developed an equation for its evolution that provided
boundary conditions for the otherwise fully resolved vapor
bubble. We are currently using a similar approach in our
own simulations of boiling and have experimented both
with the model of Son & Dhir as well as with a simpler
model that we have developed. For a simple multiscale
model of boiling in a film between drops colliding with hot
solid particles, see Ge & Fan (2006).
As we go to smaller and smaller scales, we eventually reach
a point where it is no longer fully justified to assume that
the usual continuum hypothesis is accurate. It is then
necessary to either change the modeling approach
completely by, for example, using molecular simulations or
possibly something like the dissipative particle dynamics
approach, or work with modified continuum formulations
designed to account for small scale effects. See
Koumoutsakos (2005) and Nie, Chen, E & Robbins (2004)
for recent attempts to couple molecular dynamics with
continuum simulations. The second approach leads to phase
field models where a thermodynamically consistent model
is derived for the transition zone. Phase field models have
been widely used for simulations of solidification where a
properly selected energy function leads to the proper surface
tension anisotrophy and kinetic effects. For multifluid
systems Jacqmin (1999) has used a phase field method to
study the dynamics of a moving contact line and Lowengrub
et al. (1999) has examined the draining and rupturing of thin
films. Such modeling should integrate well with the
approach taken here.

Mass transfer and chemical reactions

Mass transfer and reaction (at least simple ones)
superficially add little complexity to DNS of multiphase
flows. The species concentration is governed by a linear
advection-diffusion equation and the reaction constitutes a
(nonlinear) source term. In practice, however, the slow
diffusion of mass, compared to the evolution of the flow
scales (and sometimes temperature) and the rapidity by
which the reactions take place pose formidable stiffness and
resolution requirements. For simulations of both mass
transfer and reactions, but for relatively simple systems
where we have selected the parameters in such a way that
the systems can be computed relatively easily using
two-dimensional flow and slow reactions (see Koynov et al.,
2005; Radl, Koynov, Tryggvason & Khinast, 2008). In the
present project we intend to deal with the problem in a
different way. The key observation is that in spite of the
numerical difficulties, the solution is-in most
cases-relatively simple. Since mass diffusion is slow, the

Figure 4. The scalar diffusing from a bubble.
Coarse grid on the left, fully converged solution on
the right. The middle frame is the coarse grid with
the boundary layer model.

0.14

0.12

0.3 0.4
Time

Figure 5. The total amount of scalar in the liquid as a
function of time for three different ways of capturing
the diffusion from the bubble. The use of the
boundary layer model on the coarse grid results in
improved agreement with the fully converged

competition with advection usually results in thin boundary
layers (that require fine grid if we resolve them fully). Thus,
the concentration of a gas dissolving from a rising bubble
will change rapidly next to the bubble, but since it is the
gradient at the bubble surface that determines the rate at
which the gas leaves the bubble it must be determined
accurately. The thickness of the mass-boundary layer is
usually much thinner than the fluid boundary layer and the
fluid motion "seen" by the mass-boundary layer changes
slowly on scales comparable to its thickness. The obvious
way to accommodate the mass diffusion next to the bubble,
without fully resolving the layer by refining the grid, is
therefore to use a boundary layer approximation to capture
it. We solve a simplified model equation for the boundary
layer and use the solution of the simplified model as a
boundary condition for the region outside the boundary
layer. If the distance of a grid point from the bubble surface
is larger than the thickness of the boundary layer then that
grid point does not "see" the scalar being diffused from the
bubble, but if the layer thickness is larger then the value
predicted by the boundary layer model is used to set the grid
point value. Grid points further away from the bubble

surface are then updated using a discrete version of the
advection-diffusion equation, computed using a scheme that
is stable for the operating parameters being used. If we
parabolize the scalar advection-diffusion equation, assume a
linear profile in the boundary layer and ignore the curvature
of the boundary, then the evolution of the square of the
boundary layer thickness is given by:
d62 du d62
4D 22 --u- (2)
dt ds ds
where the velocity is interpolated from the grid used to
solve for the fluid flow, D is the mass diffusion coefficient
and s is a coordinate along the bubble surface.
In figure 4 we show the scalar field around a single bubble,
computed using a 114 x 338 grid for the Navier-Stokes
equations, which results in a fully converged solution for
the flow field. In frame (a) we use the same grid for the
advection/diffusion equation. A cantered difference scheme
for the advection term results in a very oscillatory field (due
to the high cell Reynolds number) so we have used a second
order ENO scheme for the advection terms. In frame (c),
furthest to the right we have used a much finer grid, 338 x
1010 grid points, for the scalar equation and simulations on
an ever finer grid (for early times) suggest that this is a fully
converged solution. The fluid mechanics is computed in
exactly the same way in all three simulations so the
difference is due to the treatment of the scalar. Obviously
the low grid resolution results in a much larger amount of
scalar in the liquid---as we expect. This is shown even more
clearly in figure 5 where the integrated scalar in the liquid is
plotted versus time. The coarser grid results in almost twice
as fast transfer of scalar into the liquid than the fine grid.
The middle frame in figure 4 and the blue line in figure 9
show results using the coarse scalar grid but a model for the
scalar boundary layer. Although these results are
preliminary and both the model and the coupling between
the boundary layer and the rest of the liquid are still
relatively primitive, it is clear that the scalar distribution and
the amount of scalar in the liquid are now much closer to
the fully converged solution than on the coarse grid without
the model. We note that once the scalar makes it out of the
very thin and under-resolved boundary layer, it is relatively
well represented by the coarser grid.
While we have not started to model the reactions yet, we
note that Chang, Dahm & Tryggvason (1991), Tryggvason
& Dahm (1992) and Dahm, Tryggvason & Zhuang (1996)
developed subscale models for strained diffusion flames
using an approach similar to the one we are using for the
mass diffusion. Since the reaction terms are nonlinear it is
necessary to make specific assumptions about the shape of
the species profiles to be able to integrate the species
conservation equation across the boundary layer. In the
references above it was shown that even relatively simple
profiles were sufficient to describe reasonably complex
reactions (reduced four step mechanism of Peeters & Kee,
1986). Not only was the total amount of products and heat
predicted accurately but flameout as well.
Since the multiscale models used here are intended to
reduce the resolution requirement but do not alter the
physics that is being modeled, validations are relatively
simple: we pick simple cases that we can fully resolve with
a large number of grid points (in general those will be
two-dimensional or axisymmetric) and then we coarsen the

Figure 6. The scalar diffusing from three freely
moving bubbles, computed using the model
described in the text. The top frame shows the scalar
distribution at one time and the bottom graph shows
the total scalar versus time for a fully resolved
simulation and a coarse grid simulation and a coarse
grid simulation using the model.

grid and introduce the models (see figure 3.) Figure 6 shows
the same information for a simulation of three freely and
interacting bubbles.

Conclusions

Multiphase flows play a major role in the working of Nature
and the enterprise of man. They are important for energy
production, manufacturing and chemical processes,
wastewater treatment, agriculture, and many others.
Specific problems range from the cooling of high energy
density electronic component to the production of synthetic
fuel, for example. As the use of DNS for multiphase flow
studies has increased and it has been applied to an increased
number of problems, it has become clear that in many
situations the formation of small-scale features such as thin
films or drops require excessive resolution. Here we
described a strategy to allow us to use simple analytical
models to capture features that are not resolved by a grid

appropriate for the rest of the flow.
Approximating boundary layers by known profiles is a very
old idea and although this approach has fallen somewhat out
of favor over the last couple of decades (we now generally
prefer adaptive grid refinement) it is alive and well as
wall-functions in turbulence modeling. The fundamental
idea of boundary layer modeling is, however, a natural way
of accounting for already understood phenomenon in
multiscale simulations. When something is already
understood, why recomputing it?
Developing methodologies for problems where it is
necessary to go beyond the continuum scale is a rapidly
emerging discipline and to many investigators "multiscale"
implies precisely such hybrid systems. As we have
discussed above however, even within the realms of the
continuum assumptions there are scales that need to be
bridges and we use multiscale more generally to imply any
large range of scales that needs to be accommodated and
treated differently.

Acknowledgements

This research was funded in part by NSF grant
CTS-0522581 "Multiscale simulations of multiphase
systems."