Chemical Engineering Science 61 (2006) 6486 6498

www.elsevier.com/locate/ces

Simulation of interfacial mass transfer by droplet dynamics using the

level set method
Kiran B. Deshpande, William B. Zimmerman
Department of Chemical and Process Engineering, University of Shefeld, Newcastle Street, Shefeld S1 3JD, England, UK
Received 9 November 2004; received in revised form 6 June 2006; accepted 7 June 2006
Available online 13 June 2006

Abstract
A unique approach to simulate mass transfer across the moving droplet where mass transport equations and governing equations of the levels
set method are solved separately is proposed in this work. Mass transfer coefcients of the chemical species can be computed by equating the
diffusive ux and the mass transfer ux at the interface, which are found to be of the same order of magnitude as of those obtained using an
empirical correlation. Simulations underestimate mass transfer coefcients by roughly 25% across the range of low Reynolds number studied
systematically. The level set method is used to track the motion of the interface to study droplet dynamics and mass transfer across a moving
droplet because of the ease in dening the local curvature of the interface and in capturing any topological changes. We perform various
numerical simulations by varying the physical properties of the system, in order to analyze the inuence of dimensionless numbers such as
the Reynolds number (Re), the Eotvos number (Eo) and the Morton number (M) on the shape of a buoyancy-driven droplet and compare them
with the various shape regimes of drops and bubbles reported by Clift et al. [1978. Bubbles, Drops and Particles. Academic Press, New York].
It is shown that larger deformation occurs for buoyancy-driven droplets when interfacial forces are considerably greater than viscous forces
(M
1 and Eo > 10) and the droplets are almost undeformed when viscous forces dominate interfacial forces (M > 103 and Eo > 10).
2006 Elsevier Ltd. All rights reserved.
Keywords: Simulation; Mass transfer; Drop

1. Introduction
The phenomenon of mass transfer, which is driven by
the concentration gradient of chemical species, can easily be captured for a homogeneous system by solving the
convectiondiffusion equations for species and product. However, in the case of heterogeneous reactions in laminar ows,
we rst need to track the position of the interface between
the dispersed phase and the continuous phase and then solve
the convectiondiffusion equations in order to study mass
transfer across the interface. There are various computational
methods available in the literature to model multi-phase ows:
the front tracking method (Unverdi and Tryggvason, 1992),
the volume of uid method (Rider and Kothe, 1998), the

Corresponding author. Tel.: +44 114 222 7517; fax: +44 114 222 7501.

E-mail address: [email protected] (W.B. Zimmerman).

0009-2509/$ - see front matter 2006 Elsevier Ltd. All rights reserved.
doi:10.1016/j.ces.2006.06.012

LatticeBoltzman method (Takada et al., 2000), the level set

method (Sethian, 1989; Sussman et al., 1999) and the diffuse interface method (Verschueren et al., 2001). All the above
methods have their advantages and disadvantages. The level
set method has been used in the present work to track the
motion of the interface because of its usefulness in distinctly
representing dispersed and continuous phases. It is computationally very intensive to solve the above equations simultaneously, but relatively simple to code. Hence, it is not well
studied in the literature yet it generates much interest for practical applications. In particular, the coding of re-initialization
(see Section 2.3.2) has been frequently requested from authors.
We rst review the literature available to simulate mass transfer across a moving droplet and then present a relatively simple
methodology.
Petera and Weatherley (2001) modelled mass transfer from
falling drops by calculating mass transfer ux experimentally,
as well as numerically. They solved momentum transport

K.B. Deshpande, W.B. Zimmerman / Chemical Engineering Science 61 (2006) 6486 6498

equations and mass transport equations for the dispersed phase

and the continuous phase separately using a LagrangeGalerkin
nite element method. The shape of the deforming interface
of the droplet was tracked using a particle trajectory analysis that involved the moving element method which could
lead to numerical instability. Remeshing is also required to
locate the new position of the interface and is discussed in
detail in their work. Although they have solved momentum transport equations and mass transport equations for
the dispersed and the continuous phase simultaneously, the
coupling between momentum and mass transport equations
is only one way. Mass transport equations use the velocity eld obtained from momentum transport, but momentum transport equations do not depend on mass transport
equations.
Orthogonal mapping is another potential interface tracking
method. Ryskin and Leal (1984) used the stream functionvorticity formulation in the boundary-tted co-ordinate system
to simulate the steady buoyancy-driven motion of a deformable
droplet for an intermediate Reynolds number. In the orthogonal
mapping method, the physical domains inside and outside the
droplet are transformed into a unit square computational domain. Dandy and Leal (1989) applied the above formulation to
capture the motion of a droplet with considerable deformation.
Mao et al. (2001) further applied a boundary-tted co-ordinate
system to simulate steady and transient mass transfer from
the continuous phase to a single drop dominated by external
resistance. They assumed one way coupling between momentum transport equations and mass transport equations. Hence,
uid ow was rst simulated for a buoyancy-driven droplet
and transient convectiondiffusion equations were solved
using the known uid ow. They proposed a formulation to
evaluate numerically mass transfer coefcients using the diffusion ux. The mass transfer coefcients evaluated were compared with the experiments performed in their laboratory. Li
et al. (2003) extended the numerical work of Mao et al. (2001)
to study the effect of surface active agents on mass transfer of
a solute into single buoyancy-driven droplets. The advantage
of the boundary-tted orthogonal co-ordinated system is that
mass transfer across a moving droplet can be easily achieved
for higher Reynolds number and Peclet number, which is very
difcult using the interface tracking method because of the
discontinuity at the interface. The major disadvantage of the
above method is that the actual deformation of the droplet
at the higher Reynolds number and Peclet number is not
captured. A novel way of inferring the overall mass transfer coefcients has been reported in literature by Deshpande
and Zimmerman (2005a) and is applied to study mass transfer limited characteristics by Deshpande and Zimmerman
(2005b).
The buoyancy-driven droplet rising through quiescent
medium is inuenced by various dimensionless numbers such
as Reynolds number (Re) summarizing the relative strength
of inertial to viscous forces, Eotvos number (Eo) representing
interfacial forces and Morton number (M) indicating the relative strength of viscous to interfacial forces. Clift et al. (1978)
reported various shapes of a moving droplet for different val-

6487

ues of these dimensionless numbers. In this work, we perform

numerical simulations for various operating conditions, in
order to study the inuence of various dimensionless numbers on the deformation of a moving droplet and interfacial
transport.
All the work we have reviewed so far in this section is mainly
focussed on the transport of a single species and cannot be
readily extended to multiple species. Our aim is to infer transfer rates of multiple species, which could potentially result in
cross-over due to asymmetric transfer rates. In particular, the
condition of asymmetric transfer rates is common in microuidic chemical reaction (Zimmerman, 2004). The details about
existence of cross-over phenomenon in a moving droplet can
be found in our next communication (Deshpande and Zimmerman, 2006). Recently, Mchedlov-Petrossyan et al. (2003a,b)
have developed a new theory of heterogeneous reaction that
is mass transfer limited in dispersed phases which illustrates
cross-over. Hence, we propose a methodology to evaluate mass
transfer coefcients computationally in this work as a building
block for the full heterogeneous reaction simulation. This article is organised as follows: the level set method is discussed in
the next section by illustrating the simulation of a buoyancydriven droplet. Simulations are performed for different values
for Eotvos number and Morton number to capture the deformation of a moving droplet in droplet dynamics section. The level
set method is then applied to simulate mass transfer across a
moving droplet. A numerical scheme is then proposed to evaluate mass transfer coefcients without using any empirical correlations.
2. Level set method
The level set approach was originally introduced by Osher
and Sethian (1988) to simulate the motion of an incompressible
two phase ow. In this front tracking technique, the interface
is represented by a characteristic level set function, which is
smooth in nature.
2.1. Characteristics of the level set function
The level set function used here is the signed distance function referring to the shortest distance to the front. The free surface in the two phase ow is represented as the zero level set
of a smooth function. The level set function used as an initial
condition in our simulations can be written as




2
2

(t = 0) = min (3 y), x + (y 2) 0.3 .
(1)
The level set function, as dened by Eq. (1), is applied and
solved over the entire computational domain. The contour plot
of the above-mentioned function is plotted for different level
sets as shown in Fig. 1. A particular characteristic of level set
function is that it carries a value equal to zero at the interface;
it is always positive in the continuous phase and is always
negative in the dispersed phase. These characteristics can be
represented schematically as shown in Fig. 2. Thus, the free

6488

K.B. Deshpande, W.B. Zimmerman / Chemical Engineering Science 61 (2006) 6486 6498

2.2. Governing equations

The motion of the interface is captured by solving the advection equation of the level set function, represented as
j

+ u
= 0.
jt

(2)

The velocity vector, u, required for the advection of the

level set function, is obtained by solving the incompressible
NavierStokes equations


ju
(u + (u)T ) + (u )u + p = F
jt

(3)

and
u = 0,

(4)

where  is density,  is viscosity, F is a body force term and p

is pressure.
The physical properties such as density, , and viscosity, ,
in the above equation are represented in terms of a Heaviside
function as,
 = c H (
) + d (1 H (
))

(5)

and
Fig. 1. A contour plot of the initial condition of the level set function for
different level sets.

 = c H (
) + d (1 H (
)).

(6)

The subscripts c and d denote the continuous phase and the dispersed phase, respectively. To avoid numerical artefacts, H (
)
is approximated as a smoothed Heaviside function represented
in a polynomial form and carries a value approaching 0 when


0 and approaching 1 when
> 0.
In the level set method, surface tension stress and gravitational force are incorporated as a body force term, F, in the
incompressible NavierStokes equation (Eq. (3)). Since gravitational force acts only in the vertical direction, the two components of the body force term can be written as

>0

=0

Fx = (
)nx
<0

(7)

and
Fy = (
)ny + g,

Fig. 2. Characteristics of the level set function indicating that the level set
function carries positive value in the continuous phase, negative value in the
dispersed phase, and is equal to zero at the interface.

surface is implicitly represented by the set of points for which

the level set function is always zero. A computational domain
is represented by , which is often used in this work to evaluate
the average value of a particular variable in the dispersed phase
or in the continuous phase using the characteristics of the level
set function.

(8)

where  is the surface tension, g is acceleration due to gravity,

 is local interfacial curvature, and nx and ny are the two components of the unit normal vector at the interface pointing from
the dispersed phase to the continuous phase. The surface tension term includes a (smoothed) delta function, (
) to ensure
that it is applied only at the interface, which is determined by
the position of the zero level set, which have as many uiduid
interfaces as necessary to demarcate the dispersed phase.
The unit normal vector at the interface pointing from the
dispersed phase to the continuous phase, n, can be written in
terms of the level set function as
n=

.
|
|

(9)

K.B. Deshpande, W.B. Zimmerman / Chemical Engineering Science 61 (2006) 6486 6498

The local interfacial curvature, , can be represented in terms

of the level set function as
=

.
|
|

(10)

The above equation (Eq. (10)) can be further expanded in terms

of cartesian components:
=

boundary conditions is shown in Fig. 3. With these boundary

conditions, it is ensured that there is no ow entering into or
leaving the computational domain. Thus, the motion of the
interface can be attributed solely to buoyancy force. We have
considered a two phase system, where density of the dispersed
phase, d , is 0.8 g cm3 , density of the continuous phase, c ,
is 1 g cm3 , viscosity of both the dispersed phase, d , and the

(j2
/jx 2 )(j
/jy)2 2(j
/jx)(j
/jy)(j
/jxjy) + (j2
/jy 2 )(j
/jx)2

((j
/jx)2 + (j
/jy)2 )3/2

Since the surface tension term and the local interfacial

curvature term are easily represented in terms of the level set
function as a body force, the level set method can be used
to compute any changes in topology induced due to surface
tension-driven ows. Surface tension and local curvature terms
are implemented by second order interpolation of Lagrange
quadratic elements.
2.3. Simulations of a buoyancy-driven droplet
The level set method discussed in the previous section is
applied to simulate the motion of a single drop rising due to
buoyancy. The computations are performed using a nite element method-based commercial software FEMLAB. The
advection equation of the level set function and the incompressible NavierStokes equations are coupled with each other
and are solved using a utility called Multiphysics available
in FEMLAB. In multiphysics, any number of equations can
be coupled using different application modes. In the present
simulations, an in-built application mode has been chosen for
the incompressible NavierStokes equations, and the advection
equation has been dened using the convectiondiffusion application mode. The advection equation utilises the velocity eld
calculated using the incompressible NavierStokes equations.
The dependence of the level set function in the incompressible
NavierStokes equations comes through the physical properties
as represented by Eqs. (5) and (6).
All the computations in the present work are performed in
a two-dimensional (2-D) computational domain. We have considered a rectangular cross-section of 2 cm 3 cm. The droplet
of radius 0.3 cm is introduced in the computational domain using an initial condition of the level set function,
, as dened
by Eq. (1). It should be noted that two side boundaries of the
computational domain are considered as a wall and are more
than three times radius of the droplet distance away from the
centre of the droplet. This distance is found to be sufcient so
that the velocity of the droplet is not affected due to the presence of the wall and hence, wall effects are not considered in
the present computations.
The boundary conditions initially used for numerical testing of the level set method are no slip (i.e. both components
of velocity vector are zero) on all the boundaries for the
incompressible NavierStokes application mode, and insulation (i.e. no ux across the boundary) on all the boundaries
of the computational domain for the advection of the level
set function application mode. The schematic of the above

6489

(11)

continuous phase, c , is varied from 0.01 cm1 s1 to

1 g cm1 s1 , acceleration due to gravity, g, is 980 cm s2 and
interfacial tension, , is varied from 1 to 10 dyn cm1 . The
above physical properties are chosen considering oilwater
system. The CGS system of units is used for convenience. The
droplet is expected to rise due to buoyancy force. The position
of interface at various time steps is tracked.
The level set method is found to capture the deformation of
the droplet as it rises in a column. The accuracy of the present
method can be analysed by checking for area conservation of
the dispersed phase.
2.3.1. Area conservation
In the present simulations, since there is no inow or outow
from the computational domain, correctness of the level set
method can be veried using area (volume in 3-D) conservation
analysis for a single drop rising due to buoyancy. Since it is an

Fig. 3. Schematic of boundary conditions used for simulation of the moving

droplet using the level set methodology.

6490

K.B. Deshpande, W.B. Zimmerman / Chemical Engineering Science 61 (2006) 6486 6498

of the level set function (i.e. after every iteration of the advection of
). The re-initialization equation can be represented as



2
2
j

j

j

,
(13)
= S(
) 1
+
jt
jx
jy

Cross-sectional area of a droplet

0.4
without reinitialization
with reinitialization

where t is a temporal representation of iteration, not the physical time. Clearly, t yields the distance formula differentially. S(
) is a smoothed sign function which can be written
in terms of a smoothed Heaviside function as

0.2

S(
) = 2H (
) 1.
0

0.1

0.2

0.3
0.4
Time, sec

0.5

0.6

0.7

Fig. 4. Comparison of the cross-sectional area evaluated with and without

re-initialization of the level set function, for a droplet rising due to buoyancy,
at various time steps using the level set method.

inherent disadvantage of the level set method that it does not

conserve area (volume) (Sussman et al., 1994, 1999; Lakehal et
al., 2002), area conservation analysis is particularly important.
One of the major advantages of the level set method is that
the dispersed phase and the continuous phase can be easily
distinguished using the unique characteristics of the level set
function, as shown in Fig. 2. The average cross-sectional area
of the dispersed phase can be calculated using this denition of
the level set function (Pillapakkam and Singh, 2001) as follows:

Area =
d.
(12)


0

The above integration is restricted over a domain where

0,

i.e. the dispersed phase. This can be implemented in FEMLAB
using logical operation (

0) which evaluates to 1 when true
and 0 when false. FEMLAB can compute any function of
,
including a logical function. We use this logical operation extensively to evaluate the average quantity of a particular variable in the dispersed phase, in the continuous phase or at the
interface. Using the above formulation, the cross-sectional area
of the dispersed phase can easily be calculated at various time
steps.
The change in the cross-sectional area of the dispersed phase
with time is shown in Fig. 4. It can be seen that the crosssectional area of the dispersed phase is continuously decreasing
with time, indicating that area is not conserved in the present
simulation. In the level set method, the level set function is a distance function when dened as an initial condition. But, it does
not necessarily remain a distance function at later times. Hence,
we need to force the level set function to be a distance function
after every iteration in order to ensure area conservation.
2.3.2. Re-initialization of

We can force the level set function to be a distance function
by solving a re-initialization equation after every advancement

(14)

The re-initialization equation is subject to the initial condition

(t = 0) =
(t) where the level set function,
, is obtained
from the level set method. We need to advance the solution for
the incompressible NavierStokes equations and the advection
equation of the level set function by one time step, and then
solve the re-initialization equation using the values of the level
set function obtained from the previous step as an initial condition. After solving (13), the replacement
(t)
(t )
is made.
Ideally, we need to iterate the solution of the re-initialization
equation until the level set function becomes a distance function. Since we are solving the re-initialization equation after every iteration of the level set method, 10 iterations are generally
enough to force the level set function to be a distance function.
Thus, the pseudo-time variable involved in the re-initialization
step (Eq. (13)) is considered as a dummy variable, as it is not
dependent on time steps involved in the governing equations of
the level set method.
The level set method has been employed after implementing
the re-initialization of the level set function,
, to simulate the
motion of a single droplet rising due to buoyancy. The surface
plot of y-direction velocity is overlaid by the contour plot of
the level set function at
= 0 for various time steps, as shown
in Fig. 5.
After incorporating re-initialization of
, the cross-sectional
area of the dispersed phase is found to be conserved with time
more than that without re-initialization of
, as shown in Fig. 4.
There is still a 9% loss in the cross-sectional area of the droplet
after incorporating re-initialization of
, but is better than the
32% loss in the same without re-initialization of
. The loss in
cross-sectional area of the dispersed droplet can be overcome
at the expense of computational time. This would have some
effect in estimating the average concentration of the species in
the dispersed phase and in the continuous phase which is not
accounted for in the present work. An inherent drawback of the
level set method of mass non-conservation due to the behaviour
of the level set function at later times is thus partially overcome
by incorporating re-initialization of the level set function.
The level set method, developed after incorporating reinitialization of
, seems to be robust enough to model multiphase ows. We perform the convergence test using the area
conservation analysis over the computational domain on a
coarse mesh and a rened mesh. The computational results for
the cross-sectional area of the dispersed phase evaluated for

K.B. Deshpande, W.B. Zimmerman / Chemical Engineering Science 61 (2006) 6486 6498

6491

Fig. 5. The position of the droplet rising due to buoyancy captured for various time steps, using the level set method after implementing re-initialization of

. (a) t = 0.025 s; (b) t = 0.25 s; (c) t = 0.6 s.

Fig. 6. The internal circulation within the droplet rising due to buoyancy for various time steps, in the frame of reference of the moving droplet: (a) t = 0.025 s;
(b) t = 0.25 s; (c) t = 0.06 s.

both the coarse mesh domain and the rened mesh domain are
found to be negligibly different (Deshpande, 2004).
2.3.3. Internal recirculation
The surface plot of y-direction velocity, as shown in Fig. 5
for various time steps, indicate that the velocity of the droplet
is higher than that of the continuous phase. However, the above
simulations are performed in the Eulerian frame of reference
and hence, internal recirculation inside the droplet is not readily apparent. In order to capture internal circulation within the
moving droplet, the velocity eld needs to be resolved in the
frame of reference of the moving droplet. This can be achieved
by subtracting the average droplet x- and y-direction velocity
from the x- and y-direction velocity, respectively, while plot-

ting the streamlines. The streamlines obtained in the frame of

reference of the moving droplet are shown in Fig. 6.
It can be seen that internal recirculation with the moving
droplet can be illustrated using the level set method in moving
frame of reference of the droplet. Now, we apply the level set
method to study the droplet dynamics by performing simulations for various dimensionless numbers and then extend the
same to study mass transfer across moving droplets.
3. Droplet dynamics
We have described the interface tracking method called the
level set method in the previous section by illustrating the simulations of a buoyancy-driven droplet. In this section, we ex-

6492

K.B. Deshpande, W.B. Zimmerman / Chemical Engineering Science 61 (2006) 6486 6498

tend the level set methodology further to study the variation in

the shape of drops rising due to buoyancy for different physical
parameters. Clift et al. (1978) broadly classied freely moving
bubbles and drops under the inuence of gravity as:

105
LOG M
-14
-13

104

-12
SPHERICAL-CAP

-11

103

102

REYNOLDS NUMBER, Re

(1) Spherical: The drops and bubbles are considered to be

spherical if interfacial and viscous forces are dominating
over inertial forces.
(2) Ellipsoidal: The drops and bubbles are considered to be
ellipsoidal if they are oblate with a convex surface.
(3) Spherical cap or ellipsoidal cap: The larger drops or
bubbles which tend to be at; dimpled or skirted at the
rear end fall under the category of spherical or ellipsoidal
cap.

WOBBLING

-10
-9
-8
-7
-6

Chapter 5

Chapter 7

-5

Chapter 8

-4
-3
SKIRTED

-2

Photographs of drops and bubbles under different regimes are

shown in Clift et al. (1978). The different shapes of drops
and bubbles can be better characterised by three dimensionless
numbers: Eotvos number, Morton number and Reynolds number and are dened as

ELLIPSOIDAL
-1

10

DIMPLED
ELLIPSOIDAL-CAB

SPHERICAL

Eo = g
de2 /,

(15)

3
4
Chapter 3

= g4c
/2c 3 ,

(16)

10

-1

10-2

10-1

1
10
EOTVOS NUMBER, Eo

102

103

and
Re = c de U/c ,

(17)

respectively, where g is acceleration due gravity,
 is the

density difference between the continuous phase and the dispersed phase, c is density of the continuous phase, c is viscosity of the continuous phase,  is interfacial tension between
the two phases and de is diameter of the volume-equivalent
sphere.
The effect of these dimensionless numbers on the shape of
drops and bubbles moving through the liquid due to gravity
can be very well represented by a plot of Reynolds number
against Eotvos number for different values of Morton number
(ReEoM plot), as shown in Fig. 7, from Clift et al. (1978).
Here, we perform simulations for a buoyancy-driven droplet in
laminar ow (Re is of the order of 10) by varying interfacial
force and viscous force.
We performed several numerical simulations by varying the
physical properties of the dispersed phase and the continuous
phase, in order to capture the various shapes of the moving
droplet. Numerical results for various values of the Eotvos number and the Morton number are shown in Fig. 8. The dimensionless numbers reported here are evaluated at the initial stage
using the droplet diameter. It can be seen that there is very little
deformation of the droplet at high Morton numbers, as shown
in Figs. 8(a) and (b). At high Morton numbers, such as 104
and 107 , corresponding to the cases (a) and (b), the droplet is
at the bottom right of the ReEoM plot, shown in Fig. 7. As
the Morton number is decreased further, a dimpled ellipsoidal
cap, as shown in Fig. 8(c), occurs. Morton number is further
decreased to move up on the ReEoM plot. The droplet is

Fig. 7. The shape of droplet can be well characterized in terms of dimensionless numbers such as Reynolds number, Eotvos number and Morton number. A plot of Reynolds number against Eotvos number for different Morton
number is represented along with corresponding shapes of the droplet. (Clift
et al. (1978), with permission).

skirted, as shown in Figs. 8(d) and (e), at Morton numbers as

low as 1 and 104 . Simulations performed by decreasing both
the Eotvos number and the Morton number, indicate that there
is little deformation and the droplet acquires a nearly spherical
shape (Fig. 8(f)) and is at the bottom part of the ReEoM plot.
Larger deformation occurs for buoyancy-driven droplets when
interfacial forces are considerably greater than viscous forces
(M
1 and Eo > 10) and the droplets are almost undeformed
when viscous forces dominate interfacial forces (M > 103 and
Eo > 10).
Thus, the various shapes of the droplet, including a significant deformation of the droplet as shown by the skirted drop
(Fig. 8(e)), for various dimensionless numbers, as reported by
Clift et al. (1978) are captured by this simulation scheme.
All the simulations performed hereafter in order to investigate
mass transfer across the moving droplet, deal with the case (c),
(Fig. 8(c)), with a considerable deformation of the droplet. The
purpose of these deformation studies is to provide visual proof
that the simulation technique captures the gross effects of interfacial and hydrodynamics, since quantitative measures of the
experimental studies are not readily available in the literature.
This provides plausibility for extending our simulations to estimating interfacial mass transfer in dispersed phase ows in
the next section.

K.B. Deshpande, W.B. Zimmerman / Chemical Engineering Science 61 (2006) 6486 6498

6493

Fig. 8. Spherical, ellipsoidal, dimpled and skirted shapes of the droplet are captured numerically for various values of the Eotvos number and the Morton
number using the level set methodology. (a) Eo = 360, M = 107 ; (b) Eo = 36, M = 104 ; (c) Eo = 72, M = 200; (d) Eo = 36, M = 1; (e) Eo = 36, M = 104 ;
(f) Eo = 3.6, M = 10.

4. Mass transfer across a moving droplet

Mass transfer phenomenon is generally quantied in terms
of the mass transfer coefcient, which is obtained using empirical correlations. These correlations are specic to the system
considered and are typically parameterized by a tted power
law and proportionality constant. One of the major objectives
of the present work is to simulate mass transfer phenomenon
across a moving droplet to evaluate mass transfer coefcients
without using any empirical correlations.
4.1. Physical description of the system

dispersed phase. The system that we have considered here can

be schematically represented as shown in Fig. 9.
4.2. Governing equations
Transport of the species from the continuous phase
to the dispersed phase can be captured by solving the
convectiondiffusion equation which is written for species A
and B, respectively, as follows:
jCA
+ u CA = DA 2 CA
jt

(18)

and
We are interested in a system which has two species A and
B, which are dissolved in the continuous phase but do not react
with each other in the continuous phase. The two species diffuse
into the dispersed phase. We are particularly interested in the
transport of the two species from the continuous phase to the

jCB
+ u CB = DB 2 CB ,
jt

(19)

where CA and CB are the concentrations of species A and

B, respectively, DA and DB are the diffusion coefcients of

6494

K.B. Deshpande, W.B. Zimmerman / Chemical Engineering Science 61 (2006) 6486 6498

Dispersed
Phase

Continuous
Phase

Fig. 9. Schematic representation of the system considered for studying mass

transfer across moving droplets where the two species transport from the
continuous phase to the dispersed phase.

species A and B, respectively. Species A and B are convected

by the velocity vector, u, which is obtained by solving the
governing equations of the level set methodology. The mass
transport equation contains a time-dependent term to capture
the change in the concentration of species with time, a ux
term due to convection of the species and a ux term due
to diffusion of the species. The limitations of this model to
molecular diffusion and convective mass transfer alone should
be noted. Additional auxiliary equations for non-equilibrium
effects (see Zimmerman et al., 2005) and liquidliquid phase
equilibrium (Zimmerman et al., 2006) if the thermodynamics
of phase transport are important can be modelled. Similarly, ion
transport mechanisms could be included. Since the focus here
is to develop the simplest model for interphase transport that
predicts mass transfer coefcients, the interfacial phenomena
are limited to those modelled in (18) and (19).
Ideally, we would have liked to solve simultaneously the governing equations of the level set method and the mass transport
equations for the two species, in order to simulate mass transfer across a moving droplet. Since simulation using the level
set method incorporating re-initialization of the level set function in itself is computationally intensive, solving mass transport equations simultaneously along with the level set method
turns out to be memory limited as well. We did try to solve the
simultaneous equations of mass transfer and level set method,
but without convergence given the mesh/memory limitations.
We propose an alternative numerical strategy to simulate mass
transfer of species across a moving droplet, and this is discussed
in the next section.
4.3. Numerical scheme
In the present simulations of mass transfer across a moving
droplet, we assume that there is a little effect of mass transfer
on the surface dynamics of a single moving droplet. In other
words, we are not incorporating the change in interfacial tension
due to mass transfer, i.e. surface excess concentration leading
to Marangoni stresses.

Since there is no dependence of the mass transport equations

on the level set method that is used in the present work, coupling between the mass transport equations and the governing
equations of the level set method is only in one direction. The
mass transport equation is dependent on the level set method,
but the converse is not true, unlike the level set method and
the re-initialization equation where coupling was in both directions. Thus, it is possible to decouple the mass transport
equations from the level set method. We propose a numerical
scheme to simulate mass transfer across a moving droplet that
exploits this one way coupling.
In the proposed scheme, we solve the level set method in the
rst stage to track the position of the interface which is changing with time, as discussed in the previous section and store the
solution. In the second stage, we solve the mass transport equations (Eqs. (18) and (19)) separately using the hydrodynamics
solved in the previous stage using table look up and interpolation. We also assume that the species A and B are convected
by the same velocity eld that is obtained from the level set
method. Thus, the velocity vector, u, in the convection term of
mass transport equations, is imported from the level set method.
We can now solve the mass transport equations separately
with the stored velocity eld (u and v) and the level set function
eld (
) obtained from the level set method. The velocity eld
is required to dene the ux due to convection in the mass
transport equations. Clearly, this staged approach is storage
limited which relieves memory and computational limitations
for moderate length simulations. Long simulations are storage
intensive and run up against mass storage device limits.
4.4. Numerical results
The level set method is applied to capture the motion of a
droplet rising due to buoyancy. The physical properties of the
system considered here are exactly the same as mentioned in
Section 2.3. The only difference in the present simulations is
the boundary conditions.
The inlet boundary condition is used at the bottom boundary
to dene the velocity of the uid entering into the computational domain. The y-component (vertical direction) of the inlet velocity is 1 cm s1 in the present case and is later varied to
obtain different hydrodynamic conditions. The top boundary is
dened as the outow boundary condition with pressure at the
outlet equal to zero. The schematic of the above boundary conditions is shown in Fig. 10. Thus, instead of a closed domain
(no-slip on all boundaries) as described in Section 2.3, uid
enters into and leaves the computational domain in this case.
The diffusion coefcients of the species A and B, DA
and DB , are chosen to be 1 and 5 cm2 s1 , respectively, to
induce kinetic asymmetry of the species, which plays a crucial role in a phenomenon called cross-over (seeMchedlovPetrossyan et al. (2003a?))), which is elaborately discussed in
our next communication (Deshpande and Zimmerman, 2006).
Since we are interested in transport of the species from the
continuous phase to the dispersed phase, asymmetric diffusion
coefcients DA and DB are applied in the continuous phase
and diffusion coefcients of both the species in the dispersed

Average relative droplet velocity, cm/sec

K.B. Deshpande, W.B. Zimmerman / Chemical Engineering Science 61 (2006) 6486 6498

6495

10
Re=2.88
Re=2.75
Re=2.61
Re=2.49
Re=2.36
Re=2.22
Re=2.1
Re=1.97
Re=1.83
Re=1.56

9
8
7
6
5
4
3
2
1
0

0.05

0.1
0.15
Time, sec

0.2

0.25

Fig. 11. Transient behaviour of the relative velocity of a moving droplet for
various operating conditions.

Fig. 10. Schematic of boundary conditions used for simulation of mass

transfer across the moving droplet.

phase are chosen to be relatively high, 50 cm2 s1 , to ensure

uniform concentration inside the dispersed phase. Thus, diffusion term in the present work is dened for reactant A as
DA (
> 0) + 50(

0). Similarly, diffusion term for reactant
B can be written as, DB (
> 0) + 50(

0).
The initial concentration of species A and B is chosen to 1
and 0.4 M, respectively. The boundary conditions chosen for
both the species A and B are the same, inlet boundary at the
bottom, convective outow at the top and insulation boundary
at the side walls. The inlet concentration of reactant A is 1 M
and that of B is 0.4 M.
The concentration proles of species A and B are obtained for
various time steps, as shown in Figs. 10 and 11. The position of
the interface is found to be exactly the same for the various time
steps as that obtained for the simulations of buoyancy-driven
droplet, indicating that mass transfer simulations are indeed
using the same velocity and level set function elds obtained
from the level set method.
As the transport rate of species B is higher than that of
species A and as there is no consumption of any of the species
inside the dispersed phase, the dispersed phase is equilibrated
with species B quicker than that with species A.
Now, having simulated mass transfer across a moving
droplet, we will explore the possibility of inferring mass transfer coefcients using the present numerical scheme.

phase and the continuous phase using the level set function
characteristics shown in Fig. 2. Given the average concentration of species in dispersed and continuous phases, we can infer mass transfer coefcients, and this is discussed in detail in
this section.
5.1. Formulation
In simulations of mass transfer across a moving droplet, two
species with different transport properties are used. We use
the reactant (reactant B) with a higher transport rate (diffusion
coefcient) to evaluate the mass transfer coefcient.
Average concentration of reactant B in the continuous phase,
also termed as bulk concentration, can be represented in terms
of the level set function as

>0 CB d
CB,b =

.
(20)

>0 d
Average concentration of reactant B in the dispersed phase, also
termed as surface concentration, can be represented in terms of
the level set function as

0 CB d
CB,s =

.
(21)


0 d
Concentration of reactant B is averaged over the cross-sectional
area of the dispersed phase and the continuous phase to evaluate
CB,s and CB,b , respectively.
Diffusion ux, j, across the interface can be written in terms
of the mass transfer coefcient as

5. Evaluation of mass transfer coefcient

j = (DB CB )|
=0 = B (CB,b CB,s ).

One of the signicant features of the nite element method

(FEM) is that we can estimate the concentration of species at
every node in the computational domain. From this we can calculate the average concentration of the species in the dispersed

Although it is quite common practice to equate the diffusion

ux to the mass transfer ux, it should be noted that the above
equation is not comparable as the right-hand side of the above
equation is a scalar quantity, whereas the left-hand side is a

(22)

6496

K.B. Deshpande, W.B. Zimmerman / Chemical Engineering Science 61 (2006) 6486 6498

vector quantity. Hence, we need to take a dot product of the

left-hand side with a normal vector. The above equation can be
written to evaluate the mass transfer coefcient as
B =

n.(DB CB )|
=0
.
CB,b CB,s

(23)

A similar denition of mass transfer coefcient was found in

the literature for bubbles and bubble swarms by Bailey and Ollis (1986). Mao et al. (2001) reported a similar denition for
numerical simulation of steady and transient mass transfer to
a single drop, dominated by external resistance. As discussed
earlier, they used a boundary-tted orthogonal coordinate system, which essentially does not capture the deformation of the
interface in highly deforming systems. Also, the droplet is not
physically moving in their study, and hence true motion of a
moving droplet cannot be captured. The formulation discussed
in this section is applied for simulations of mass transfer across
a moving droplet for different operating conditions, to evaluate the mass transfer coefcient, which is then compared with
empirical correlations described in the next section.
5.2. Empirical correlation
Mass transfer phenomenon is often quantied in terms of
a mass transfer coefcient which is obtained using empirical correlations. These correlations are generally represented
in terms of dimensionless numbers like the Sherwood number
(Sh = B d/DB ), the Reynolds number (Re = dUR c /c ) and
the Schmidt number (Sc = c /c DB ), where c and c are
density and viscosity of the continuous phase and d is the diameter of the droplet. The relative velocity, UR , used in these
correlations is the velocity of the rising or falling droplet or
bubble, relative to the velocity of the continuous phase. In the
present simulations, since we are dealing with 2-D simulations,
the empirical correlation for a cylinder (McCabe et al., 2000)
could be more appropriate and is represented as
Sh = 0.61Re1/2 Sc1/3 .

(24)

Mass transfer coefcient obtained using empirical correlation

is later compared with the numerically evaluated mass transfer
coefcient, for various operating conditions.
5.3. Numerical results
One of the major advantages of the present simulations is
that we can numerically estimate the mass transfer coefcient
by comparing the diffusion ux and the mass transfer driving
force, as discussed in the previous section. In this section, we
compare numerically evaluated mass transfer coefcients with
those obtained using empirical correlation. The mass transfer
coefcient is often represented in a dimensionless form in terms
of the Sherwood number, Sh, which is dened as a function of
other dimensionless numbers such as the Reynolds number, Re,
and the Schmidt number, Sc. The relative velocity of the droplet
is used in the above-mentioned dimensionless numbers. The
effect of hydrodynamic conditions on mass transfer coefcients

can be brought out by changing the relative velocity of species.

Therefore, we perform simulations of mass transfer across a
moving droplet by changing the hydrodynamic conditions. We
rst discuss how we dene the relative velocity in the present
simulations.
5.3.1. Relative velocity
In the numerical simulations of a droplet rising due to buoyancy, both the continuous phase and the dispersed phase are
moving in the computational domain. The effective velocity of
a moving droplet depends on average velocity of a continuous
phase. Hence, we dene the average relative velocity as the
average velocity of a moving droplet relative to the average velocity of the continuous phase. Average velocity of a droplet is
dened in terms of the level set function as

0 v d
,
(25)
Udrop =

0 d
where v is the y-component of velocity eld. Average velocity
of the continuous phase is dened in terms of the level set
function as

>0 v d
Ubulk =

.
(26)

>0 d
Averaging is done over the cross-sectional area in a similar
fashion as was done to evaluate average bulk and surface concentration, as shown by Eqs. (20) and (21), respectively. The
relative velocity, UR can be written as
UR = Udrop Ubulk .

(27)

The relative velocity, as dened by Eq. (27), is plotted against

time as shown in Fig. 11, which indicates the transient behaviour of the relative velocity for various operating conditions.
The inlet velocity is changed, ranging from 1 to 5 cm s1 . The
Reynolds number is calculated on the basis of the diameter of
the undeformed droplet, which is the case at the initial stage. It
is clearly seen that relative velocity of a moving droplet changes
with time and eventually attains a steady state.
Since relative velocity of the moving droplet and the shape
of the droplet change with time, hydrodynamic conditions are
changing with time for a buoyancy-driven droplet. However,
the deformation of the droplet is not considered in the empirical
correlation, which are proposed for a rigid particle. The dimensionless numbers calculated in order to evaluate mass transfer
coefcients using empirical correlation are based on the diameter of the undeformed droplet. The assumption of a rigid particle is reasonable at the initial stage (t = 0.025 s) and hence
mass transfer coefcients obtained by the level set methodology are compared with empirical correlations only at the initial
stage.
5.3.2. Mass transfer coefcient
The relative velocity, discussed in the previous section, is
used to calculate Pe and Re, and is varied by changing the inlet

K.B. Deshpande, W.B. Zimmerman / Chemical Engineering Science 61 (2006) 6486 6498
Table 1
Comparison of mass transfer coefcients evaluated using the various empirical
correlation and those obtained using the level set methodology, for various
operating conditions
Inlet
velocity

Relative
velocity

Reynolds
number

Peclet
number

Empirical
MTC

Simulation
MTC

0.5
1
1.5
2
2.5
3
3.5
4
4.5
5.0

2.625
3.05
3.275
3.5
3.7
3.925
4.15
4.35
4.575
4.8

1.575
1.83
1.965
2.1
2.22
2.355
2.49
2.61
2.745
2.88

0.315
0.366
0.393
0.42
0.44
0.471
0.498
0.522
0.549
0.576

3.731
4.02
4.17
4.31
4.43
4.56
4.69
4.8
4.93
5.04

2.5
2.57
2.55
2.68
2.66
2.65
4.19
4.0
3.68
3.85

velocity of species. The mass transfer coefcients evaluated using empirical correlation and using simulations are represented
in Table 1 for various operating conditions and the corresponding non-dimensional numbers.
The level set methodology is applied to the simulations of
mass transfer across a moving droplet for various operating
conditions to evaluate the mass transfer coefcient. The mass
transfer coefcient can be calculated by using Eq. (23), where
we need to calculate the diffusive ux at the interface and
the average bulk and surface concentrations of reactant B. The
average bulk and surface concentrations are calculated using
Eqs. (20) and (21), respectively. The diffusive ux at the interface can be calculated using the characteristics of the level
set function shown in Fig. 2. To avoid any arbitrary numerical noise, we consider the interface to be a region where
0.015


0.015. It is ensured that the interface region is
indeed very small relative to the droplet diameter.
Mass transfer coefcients at the initial stage (t = 0.025 s),
where there is no deformation of the droplet, are evaluated using the level set methodology and are tabulated in Table 1 for
various operating conditions. It was observed that the numerically evaluated mass transfer coefcient decreased initially and
then increased at later time steps time for all the operating
conditions considered in the present study (results not reported
here). The initial decrease can be attributed to the deformation
of the droplet, but this could be attributed to the decreasing
diffusion ux at the interface, whereas increase in mass transfer coefcients at later stages is due to average concentration
gradient (difference in the average bulk and surface concentrations), which is a decreasing due to transport of species from
the continuous phase and accumulating in the dispersed phase.
We could not obtain appropriate empirical correlations for the
transient behaviour of mass transfer coefcients and hence,
numerically observed transient behaviour of the mass transfer
coefcient for a buoyancy-driven droplet is not reported here.
Since the mass transfer coefcient obtained using the present
numerical scheme is changing with time, and also since there
is a considerable deformation of the rising droplet (moderate
Eo and M) at later stages, we will consider mass transfer coefcients calculated at the initial stages for comparison with

6497

empirical correlations which do not take deformability into account.

The mass transfer coefcients obtained using empirical correlation and that obtained using the level set methodology are
tabulated against in Table 1, along with various operating conditions. Mass transfer coefcients evaluated using the level set
methodology are found to be of the same order of magnitude
as those evaluated using the empirical correlation. Roughly,
they underestimate mass transfer coefcients by 25% across
the range of low Reynolds number studied systematically. This
suggests that accurate estimates are possible by taking into account the underestimate.
Thus, mass transfer coefcients can be evaluated using the
present formulation without using any empirical correlations
and have been compared with the empirical correlation to check
the correctness of the methodology. The transient nature of
mass transfer coefcients still remains an unexplored eld and
requires further attention.

6. Conclusion
A numerical scheme, called the level set method, is discussed in this article to track the motion of the interface
for a buoyancy-driven droplet. The inherent disadvantage
(mass non-conservation) of the level set method of mass
non-conservation is addressed and is improved upon by implementing re-initialization of the level set function at every
time step. After incorporating re-initialization, 91% of area
(volume) has been conserved, which is far better than without
re-initialization where 32% loss has been observed.
The shape of a buoyancy-driven droplet is inuenced by
dimensionless numbers such as the Reynolds number (Re), the
Eotvos number (Eo) and the Morton number (M). We performed
various numerical simulations for different sets of the above
mentioned dimensionless numbers and analysed the shapes of
the droplets by comparing them with the shape regimes of drops
and bubbles reported by Clift et al. (1978). It is shown that the
larger deformation occurs for buoyancy-driven droplets when
interfacial forces are considerably greater than viscous forces
(M
1 and Eo > 10) and the droplets are almost undeformed
when viscous forces dominate interfacial forces (M > 103 and
Eo > 10).
Mass transfer across a moving droplet is simulated by
adopting a two stage approach where the convectiondiffusion
equations for mass transfer are decoupled from the governing
equations of the level set methodology, which is otherwise
computationally very intensive. The major result of the present
formulation is that we can accurately infer mass transfer coefcients without using any empirical correlations. Mass transfer
coefcients obtained using simulations are found to be of the
same order of magnitude of those obtained using empirical
correlations.
The simulation of dispersed phase mass transfer developed
here can be reliably extended to predict more general features of
heterogeneous reactions and interfacial dynamics of multiphase
ow (Deshpande and Zimmerman, 2005b,2006).

6498

K.B. Deshpande, W.B. Zimmerman / Chemical Engineering Science 61 (2006) 6486 6498

Acknowledgements
KBD would like to thank Dr. B. Hewakandamby for very
useful discussions. W.B. Z. would like to thank the EPSRC
(GR/A01435, GR/S67845, GR/R72754 and GR/N20676).
References
Bailey, J.E., Ollis, D.F., 1986. Biochemical Engineering Fundamentals.
McGraw-Hill, New York.
Clift, R., Grace, J.R., Weber, M.E., 1978. Bubbles, Drops and Particles.
Academic Press, New York.
Dandy, D.S., Leal, L.G., 1989. Buoyancy-driven motion of a deformable drop
through a quiescence liquid at intermediate Reynolds numbers. Journal of
Fluid Mechanics 208, 161192.
Deshpande, K.B., 2004. Study of transport limited heterogeneous reaction in
the dispersed phase. Ph.D. Thesis, University of Shefeld, UK.
Deshpande, K.B., Zimmerman, W.B., 2005a. Experimental study of mass
transfer limited reaction Part I: use of bre optic spectrometry to infer
asymmetric mass transfer coefcients. Chemical Engineering Science 60,
28792893.
Deshpande, K.B., Zimmerman, W.B., 2005b. Experimental study of mass
transfer limited reaction Part II: existence of cross-over phenomenon.
Chemical Engineering Science 60, 41474156.
Deshpande, K.B., Zimmerman, W.B., 2006. Simulations of mass
transfer limited reaction in a moving droplet to study transport
limited characteristics. Chemical Engineering Science, in press, doi:
10.1016/j.ces.2006.06.013.
Lakehal, D., Meier, M., Fulgosi, M., 2002. Interface tracking towards the direct
simulation of heat and mass transfer in multiphase ows. International
Journal of Heat and Fluid Flow 23, 242257.
Li, X., Mao, Z.S., Fei, W., 2003. Effects of surface-active agents on mass
transfer of a solute into single buoyancy driven drops in solvent extraction
systems. Chemical Engineering Science 58, 37933806.
Mao, Z.S., Li, T., Chen, J., 2001. Numerical simulation of steady and
transient mass transfer to a single drop dominated by external resistance.
International Journal of Heat and Mass Transfer 44, 12351247.
McCabe, W., Smith, J., Harriott, P., 2000. Unit Operations of Chemical
Engineering. McGraw-Hill, New York.
Mchedlov-Petrossyan, P.O., Khomenko, G., Zimmerman, W.B., 2003a. Nearly
irreversible, fast hetrogeneous reactions in premixed ow. Chemical
Engineering Science 58, 30053023.

Mchedlov-Petrossyan, P.O., Zimmerman, W.B., Khomenko, G., 2003b. Fast

binary reactions in a heterogeneous catalyic batch reactor. Chemical
Engineering Science 58, 26912703.
Osher, S., Sethian, J.A., 1988. Front propagating with curvature-dependent
speed: algorithm based on HamiltonJacobi formulations. Journal of
Computational Physics 79, 1249.
Petera, J., Weatherley, L.R., 2001. Modelling of mass transfer from falling
drops. Chemical Engineering Science 56, 49294947.
Pillapakkam, S.B., Singh, P., 2001. A level set method for computing solutions
to viscoelastic two-phase ow. Journal of Computational Physics 174,
552578.
Rider, W.J., Kothe, D.B., 1998. Reconstructing volume tracking. Journal of
Computational Physics 141, 112152.
Ryskin, G., Leal, L.G., 1984. Numerical solution of free-boundary problems
in uid mechanics. Part 2. Buoyancy-driven motion of a gas bubble through
a quiescence liquid. Journal of Fluid Mechanics 148, 1935.
Sethian, J.A., 1989. Level Set Methods and Fast Marching Methods.
Cambridge University Press, Cambridge.
Sussman, M., Smereka, P., Osher, S., 1994. A level set approach for computing
solutions to incompressible two-phase ow. Journal of Computational
Physics 114, 146159.
Sussman, M., Almgren, A.S., Bell, J.B., Colella, P., Howell, L.H., Welcome,
M.L., 1999. An adaptive level set approach for incompressible two-phase
ows. Journal of Computational Physics 149, 81124.
Takada, N., Misawa, M., Tomiyama, A., Fujiwara, S., 2000. Numerical
simulation of two and three dimensional two phase uid motion by Lattice
Boltzmann method. Computer Physics Communications 129, 233246.
Unverdi, S.O., Tryggvason, G., 1992. A front tracking method for viscous,
incompressible, multi-uid ows. Journal of Computational Physics 100,
2537.
Verschueren, M., van de Vosse, F.N., Meijer, H.E.H., 2001. Diffuse-interface
modelling of thermocapillary instabilities in a HeleShaw cell. Journal of
Fluid Mechanics 434, 153166.
Zimmerman, W.B., 2004. Heterogeneous catalysis in microchannel reactors.
ESDA Manchester, ASME Transactions 58441.
Zimmerman, W.B., Mchedlov-Petrossyan, P.O., Khomenko, G.A., 2005.
Nonequilibrium effects on fast binary reactions in a heterogeneous catalytic
batch reactor. Chemical Engineering Science 60, 30613076.
Zimmerman, W.B., Mchedlov-Petrossyan, P.O., Khomenko, G.A., 2006.
Effects of transport and phase equilibrium on fast, nearly
irreversible reactive extraction. Chemical Engineering Science, in press,
doi:10.1016/j.ces.2006.04.011.