7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
A VOFBased Conservative Method for the Simulation of Reactive Mass Transfer
from Rising Bubbles
D. Bothe* M. KrOger* and H.J. Warnecket
Center of Smart Interfaces, Technische Universitit Darmstadt, Darmstadt, Germany
bothe@csi.tudarmstadt.de and kroeger@csi.tudarmstadt.de
t Department of Chemical Engineering, University Paderborn, Germany
wamecke@tc.upb.de
SCorresponding author
Keywords: Reactive mass transfer, fluidic interface, Volume of Fluid
Abstract
In this paper numerical results on reactive mass transfer from single gas bubbles to a surrounding liquid are presented.
The underlying numerical method is based on the solution of the incompressible two phase NavierStokes equations.
The VolumeofFluid method is applied for the description of the liquidgas interface. Within the numerical approach
the concentration of the transfer component is represented by two separate variables, one for each phase. Numerical
results are in good agreement with experimentell data.
1 Introduction
Many chemical reactions are heterogeneous, in which
case the reaction partners are present in different phases,
e.g. in a liquid and in a gas. At least one of these
species has to cross the interface between the two phases
in order to enable the chemical reaction. This process is
called reactive mass transfer, if the physical transfer of
one chemical component from one phase to the other is
followed by a chemical reaction in this phase. A com
mon example is the selective oxidation of cyclohexane
where the gaseous oxygen is first dissolved and then re
acts with the organic liquid to the desired product. This
type of reaction is of high interest, especially for the
chemical process industry. The scaleup from labora
tory sized model apparatuses to large industrial reactors
often turns out to be quite difficult, because several of
the occurring phenomena in such systems are not fully
understood. Direct numerical simulations have proven
to be a useful addition to small scale experiments and
theoretical analysis. Accurate simulations can give in
sight to local details which are otherwise not achievable.
Therefore, a lot of work has been devoted to this topic.
VOFbased simulations of purely physical mass transfer
across deforming interfaces without chemical reaction
have been reported in Davidson and Rudman (2002) and
in Bothe et al. (2004). In the latter paper, transfer of
oxygen from air bubbles rising in water or aqueous so
lutions has been simulated, taking into account the re
alistic jump discontinuity of the oxygen profiles at the
interface. Onea et al. (2009) used a similar approach as
Bothe et al. (2004) to simulate mass transfer in upward
bubble train flow through square and rectangular mini
channels. Darmana et al. (2006) performed 3D simu
lations of mass transfer at rising fluid particles using
the Front Tracking method. There, the transport resis
tance inside the fluid particle is neglected, i.e. a constant
concentration value inside the bubble is assumed. Radl
et al. (2007) performed 2D simulations of deformable
bubbles and bubble swarms with mass transfer in non
Newtonian liquids using a semiLagrangian advection
scheme. To prevent stability problems, a reduced den
sity ratio between gas and liquid is used. Recently, first
papers on numerical simulation of reactive mass trans
fer appeared. In (Khinast et al. 2005), reactive mass
transfer at deformable interfaces is examined using a
2D Front Tracking/Front Capturing hybrid method. In
(Deshpande and Zimmermann 2006a), a Level Set based
method is used to simulate mass transfer across the in
terface of a moving deformable droplet. This method
is extended to reactive mass transfer in (Deshpande and
Zimmermann 2006b) where an instantaneous chemical
reaction occurs inside a moving droplet, which leads
to a quasistationary problem for the mass transfer. In
Radl et al. (2008) 2D simulations are performed using
a FrontTracking method to investigate the effect of dif
ferent Hatta and Schmidt numbers on the catalytic hy
drogenation of nitroarenes for single bubbles and bub
ble clusters. Based on the numerical approach in (Bothe
et al. 2004), in (Bothe et al. 2009) reactive mass trans
fer with parallel consecutive reactions at rising bubbles
is analyzed by using a local selectivity. The approach
from (Bothe et al. 2004) employs a single scalar field for
any transfer component and this scalar quantitiy refers to
a normalized concentration for which the jump disconti
nuity can be removed. This has some disadvantages con
cerning conservativity of the transfer component which
can lead to artificial mass transfer. The latter is caused
by a relative motion of the interface and the concentra
tion discontinuity due to the use of a different advection
algorithm. In (Alke et al. 2009) a new VOFbased ap
proach has been introduced which does not posses these
drawbacks. In the present paper, this approach is ex
tended with a subgrid scale model and is applied to re
active mass transfer from single rising bubbles.
2 Governing Equations
In the mathematical model it is assumed that the sys
tem under consideration consists of two immiscible, in
compressible fluids, separated through a moving and
deformable interface. The interface between the two
phases is presented as a mathematical (sharp) surface of
zero thickness and is denoted by E(t). The local bal
ances for mass and momentum can then be written as
V u 0 (1)
for mass and
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
where K = V ne is the curvature (more precisely,
the sum of the principal curvatures) and expresses the
Laplace pressure jump in stationary cases. It is further
assumed that there is at least one transfer component
with the volumetric molar concentration Ck, being sol
uble in both phases. In case of ideally diluted systems
with small pressure gradients the components have no
influence on the hydrodynamics and can therefore be
considered as passive scalars, transported by the velocity
field. The local balance for the species thus reads as
C+V.(ce u+jk)= R, in Q (t)UQd(t), (6)
at
governing the species transport inside the bulk phases.
At the interface the transmission condition
[jc] ne = 0
holds and the jump condition
[Ip ] 0
is imposed. The diffusive fluxes are given by Fick's law,
ji = Dk Vck
with the constant diffusion coefficient Dh > 0, and
the continuity of the chemical potential is expressed by
Henry's law, i.e.
k/c' = Hk (10)
with a Henry coefficient Hk / 1. The source term Rk
on the righthand side of 6 accounts for chemical reac
tions. Note that the square bracket stands for the interfa
cial jump according to
Vp + pg + V S (2)
for the momentum. In this equation S denotes the vis
cous stress tensor given by
S Vu+ (VuT).
The interface normal unit vector nE points into the con
tinuous phase for the remainder of this paper. This set of
equations represents the onefield formulation of the two
phase NavierStokes equations where the material prop
erties p and refer to the phase dependent density and
viscosity. Additionally, the following jump conditions
are satisfied at the interface:
[u] 0 (4)
[pI S] nr = a K cn,
S(x h hns))
(11)
[4] (xi) = limr ( (xe + hns)
where 0 is an arbitrary scalar.
3 Numerical Method
The set of governing equations in Section 2 is solved
with the inhouse code Free Surface 3D (FS3D), Rieber
(2004). A finite volume discretization is used for the
spatial discretization and an explicit Eulerian scheme for
the time discretization. The volumes are arranged on
a Cartesian, staggered grid, where scalar variables like
pressure or concentration are stored as cell centered val
ues and the velocities are stored on the centers of the cell
faces. The code employs the Volumeoffluid method.
Here, an additional transport equation is solved to keep
track of the location of the different phases and, thereby,
of the interface. This type of method is called a vol
ume tracking scheme, since only the different phases
du
p +p(u. V)u
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
are transported and the interface is geometrically recon
structed from that information. In FS3D, the PLIC al
gorithm (Rider and Kothe 1998) is employed for inter
face reconstruction. The code uses a onefield formula
tion of the twophase NavierStokes equations in which
the volumetric surface tension force is incorporated via
the conservative continuum surface stress (CSS)model
of Lafaurie (Lafaurie et al. 1994). The volume fraction
transport equation reads as
S+ V (u f) 0, (12)
where f is the phase indicator function of the dis
persed phase domain d (t). In the FV discretization
scheme employed here, the cell centered value of f cor
responds to dispersed phase fraction inside a computa
tional cell. Within the VOFmethod, the phase related
material properties are given by
S= fPd + (1 f)pc (13)
and
i = frd + (1 f)rq. (14)
Transport of molar species mass
For the computation of the transport of a transfer species
k, the concentration is represented by two separate scalar
variables according to
.,d { ck(x,t) for xE d(t)
'x, t) 0 for x E Q,(t) (15)
and
(x, t) 0 for xE Qd(t)
'>'= ck(x,t) for x E Qc(t). )
This allows for a representation of the individual one
sided limits of the concentration at the interface. Within
the FV discretization, these scalars are related to the cell
volume V, i.e. the cell centered values are given as
(d 1
(M = I
SIVI
(1 = f
V JvnoW)
Ck(t)dV
Ck (t)dV.
The computation is carried out with a directional split
ting. This means, that for each of the three dimensions
a one dimensional transport step is calculated. The or
der of the directional steps is changed in every time step
to reduce systematical errors. Inside the bulk phases the
convective transport of the volume fraction f is calcu
lated with a simple first order upwind scheme. Since
f has only one discrete value in the bulk phase, this is
accurate. The concentration of a chemical species, how
ever, can take arbitrary, non negative values. A first or
der upwind scheme applied in that case would lead to
unacceptable numerical diffusion. Therefore the convec
tive transfer is based on the limiter scheme of Van Leer
(Van Leer 1979). In this scheme the concentration in
side a computational cell is approximated by a linear
function instead of a constant value. The slope of this
linear function is then used to extrapolate the concentra
tion onto the cell face, where the convective fluxes are
calculated. This is restricted in a manner that no new
maxmia or minima are created in the solution. The lim
iter scheme reduces numerical diffusion in areas where
the concentration has steep gradients, especially in the
wake of the bubble. In areas with flat gradients it falls
back to the first order upwind scheme. This allows us to
use the advantage of a higher order scheme without tak
ing the drawback of over or undershoots, respectively,
which are often created by higher order schemes. At
the interface the same PLIC algorithm is used as for the
ffield. Since both onesided concentration limits are
known, a clean distinction between the convective fluxes
for the two phases can be made. This ensures that no ar
tificial mass transfer due to convection occurs. Diffusive
fluxes inside the bulk phases are obtained by a standard
central differencing scheme. In interfacial cells, differ
ent cases have to be distinguished. Diffusion between
an interfacial cell and a bulk phase cell is calculated for
the appropriate scalar only, using the according diffusion
coefficient for that phase. For the diffusion between two
interface cells two diffusive fluxes are calculated, one
for each phase with the respective diffusion coefficients.
The effective area in that case is approximated by the
cell face area multiplied by the fraction of the respec
tive phase. To ensure conservativness inside cells with
very small values of f, the diffusive fluxes are restricted.
Since diffusion can only take place until the concentra
tion gradient is zero, this can be used as a limit for the
diffusion in one time step.
Mass transfer source term
The mass transfer between the phases takes place in the
cells carrying a part of the interface and is accounted for
by means of source terms for the two scalars represent
ing the full concentration field. In the following, a brief
description of two possible ways to calculate the inter
facial source term is given. For more details please we
refer to (Kroeger et al.). The calculation of the source
terms is restricted to interfacial cells. Inside those cells,
the interface is assumed to be planar. Further, it can
be assumed that inside an interfacial cell, the concen
tration in the dispersed phase is homogeneous due to
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
the high diffusion coefficient. In a one dimensional ap
proach, one can then calculate the onesided interfacial
limit value of the concentration for the dispersed phase
by simply extrapolating the cell centered value onto the
interface. The limit value for the continuous phase can
be obtained by using Henry's law (10). Together with
one neighboring value from the continuous phase, a lin
ear approximation for the gradient of the concentration
at the interface can be obtained from this value. Insert
ing this gradient into Fick's law ( 9) yields the molar flux
of the transfer species across the interface for the contin
uous phase. Since (7) states that both fluxes have to be
equal, this also gives directly the flux for the dispersed
phase, leading to a conservative algorithm. In a dimen
sional splitting scheme, the fluxes are calculated for each
direction separately. The second method is quite similar
to the linear gradient method. Here, instead of a linear
approximation a different function is used to estimate the
gradient. The type of function used in this subgrid scale
model is obtained analytical from the mass transfer for
an overflown plane (Bird et al. 2002) and has the form
o. 25
J
S0 15
o,.z
0
8
5
Sc = 1000
dial distance om paicle ce in m
radial distance from particle race in m
Sc 10000
radial distance rom particle surface in m
Ck(x,y) = ck(
erf))
69 n
8y 2 D Y
V
It can be shown in simulations that in case of
fine grids, both methods yield the same results. I
ever, with the subgrid scale model the result can air
be obtained with a coarser resolution. Therefore
method is to be preferred. The only drawback he
that in this form it can be only used in situations, w
mass transfer resistance in the dispersed phase is si
All numerical results in the next sections are obta
employing the subgrid scale model.
Figure 1: Comparison of exact (lines) and with the
(19) VOFmethod obtained concentration profiles at the
equator of a rising fluid particle with Re = 0.284, Sc =
1000 (top) and Re=0.284, Sc=10000 (bottom).
The Reynolds number is 0.284 for these cases. The pro
very files show very good agreement in this case.
low As a measure for the overall mass transfer the integral
eadv Sherwood number
this
re is
There
mall.
lined
4 Validation
In this section, comparisons between theoretical and nu
merical as well as between experimental and numeri
cal result are used to validate the described numerical
approach. The first comparison employs the flow field
obtained by Hadamard (1911) and Rybczynski (1911)
which is valid for Reynolds number below 0.3. In Fig
ure 1 comparisons between the theoretical and simulated
concentration profiles at the equator of a rising bubble
are shown. In the numerical setup the bubble equivalent
diameter is set to dc = 4 mm. The computational do
main for this 3D simulation has the dimension of 4x2x2
dp and is discretized with 64 cells per diameter. The vis
cosity of the continuous phase is set to 460 mPa s. Due
to the high viscosity, the Schmidt number is quite large.
SL
Sh = (20)
can be used. There are different theoretical correlations
for the Sherwood number for certain special cases:
Sh= 1 + (1 + 0.546Pe2) (21)
which is assumed to be valid for Re 0 and arbitrary
Peclet numbers Pe = ReSc (Clift et al. 1978),
pel1.72
Sh = 2 + 0.651 1.22 (22)
which is assumed to be valid for Re 0 and Sc 
oc (Oellrich et al. 1973). For physical mass transfer
the numerical results compare reasonably well with the
experimental data as can be seen in Figure 2. Since the
setup for the experiment differs slightly from that of the
simulation, the pictures are only qualitatively compared.
5 Results
The method described above can be applied to simu
late reactive mass transfer with simple and complex re
o. OO
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Re 0.284 Re 0.284
Sc 1000 Sc 10000
Numerical Result 18.5 64.1
Eq. (21) 12.0 35.1
Eq. (22) 13.0 36.7
Table 1: Comparison of integral Sherwood numbers
Figure 2: Distribution of oxygen in the wake of a rising
bubble. LiF experiment (left) dp 1.6 mm, numerical
3D Simulation (right) dp 1.8 mm.
actions inside one of the phases. As an example, the
simulation of the metal catalized oxidation of sulfite is
considered in this section. The overall reaction scheme
for this reaction is given as
HS03 + 0, i HS04 (23)
with k* given as the gross reaction rate assumed con
stant for this reaction. Note that the complete reaction
scheme is far more complex, hence there is a radical re
action mechanism involved. For simplification, the gross
reaction is taken into account. In a more abstract nota
tion this results in a reaction scheme of the type
A + B k* P. (24)
In this scheme A denotes the transfer component (oxy
gen), B is the dissolved component in the liquid phase
(the hydrogen sulfite ion) and P is the desired product
(the hydrogen sulfate ion). In the left part of Figure 3 a
snapshot of an LiFexperiment for the sulfite oxidation
is shown. In the center of the wake it can be seen that the
oxygen is depleted because of the reaction, thus creating
a gap in the wake. The numerical result reproduces this
gap very well, where so far only a qualitatively compar
ison is possible.
Figure 3: Concentration distribution of oxygen
6 Conclusions and Outlook
A numerical approach for the simulation of reactive
mass transfer is introduced in this contribution. It is vali
dated for the case of physical mass transfer. The method
is applied to the case of a simple reaction and the numer
ical results are compared with experimental data. Vari
ation of the gross reaction rate constant can be used to
simulate the effect of rising temperature on this reaction,
which we are currently doing. Since it is quite difficult
to measure the velocity at high resolution in an exper
iment, it is also planned to combine the simulated ve
locity fields with measured concentration fields in order
to obtain integral Sherwood numbers from the experi
ments.
7 Acknowledgements
The authors gratefully acknowledge financial sup
port from the Deutsche Forschungsgemeinschaft (DFG)
within the DFGproject "Reactive mass transfer from
rising gas bubbles" (PAK119). We also thank Prof.
Michael Schiiter (TH HamburgHarburg, Germany) and
Prof. Norbert Rabiger (University Bremen, Germany)
for providing the experimental results.
References
A. Alke, D. Bothe, M. Kroeger, and H.J. Warnecke.
VOFbased simulation of conjugate mass transfer from
freely moving fluid particles. In Mammoli, AA and
Brebbia, CA, editor, COMPUTATIONAL METHODS IN
MULTIPHASE FLOW V, WIT Transactions on Engi
neering Sciences, pages 157168, 2009.
R. B. Bird, E. S. Warren, and E. N. Lightfood. Trans
port Phenomena, 2nd ed. John Wiley & Sons, Inc.,
New York, Chichester, Weinheim, Brisane, Singapore,
Toronto, 2002.
D. Bothe, M. Koebe, K. Wielage, J. Prtiss, and H.J.
Warnecke. Direct numerical simulation of mass trans
fer between rising gas bubbles and water. In M. Som
merfeld, editor, Bubbly Flows: Analysis, Modelling, and
Calculation, Heat and Mass Tranfer, pages 159 174.
Springer, Berlin, 2004.
D. Bothe, M. Kriger, A. Alke, and H.J. Warnecke. Vof
based simulation of reactive mass transfer across de
formable interfaces. PCFD, 9:325 331, 2009.
J.R. Clift, M.E. Grace, and M.E. Weber. Bubbles, Drops,
and Particles. Academic Press, New York, 1978.
D. Darmana, N.G. Deen, and J.A.M. Kuipers. Detailed
3d modelling of mass transfer processes in twophase
flows with dynamic interfaces. Chem. Eng. Technol., 29:
1027 1033,2006.
M.R. Davidson and M.J. Rudman. Volumeoffluid cal
culation of heat or mass transfer across deforming inter
faces in twofluid flow. Numer Heat Transfer Part B,
41:291 308, 2002.
K.B. Deshpande and W.B. Zimmermann. Simulation
of interfacial mass transfer by droplet dynamics using
the levelset method. Chem. Eng. Sci., 61:6486 6498,
2006a.
K.B. Deshpande and W.B. Zimmermann. Simulation
of mass transfer limited reaction in a moving droplet to
study transport limited characteristics. Chem. Eng. Sci.,
61:6424 6441, 2006b.
J.S. Hadamard. Mouvement permanent lent d'une
sphere liquid et visquese dans un liquid visqueux. C.R.
Acad. Sci. Paris, 152:1735 1738, 1911.
J.G. Khinast, A. Koynov, and G. Tryggvason. Mass
transfer and chemical reactions in bubble swarms with
dynamic interfaces. AIChE J., 51:2786 2800, 2005.
M. Kroeger, A. Alke, and D. Bothe. A VOFbased
method for mass transfer processes at fluidic particles.
in preparation.
B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, and
G. Zanetti. Modelling merging and fragmentation in
multiphase flows with surfer. J. Comp. Phys., 113:134 
147, 1994.
H. Oellrich, H. SchmidtTraub, and H. Brauer.
Theoretische Berechnung des Stofftransports in der
Umgebung einer Einzelblase. Chem. Eng. Sci., 28:711 
721, 1973.
A. Onea, M. Wirner, and D.G. Cacuci. A qualitative
computational study of mass transfer in upward bubble
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
train flow through square and rectangular minichannels.
Chem. Eng. Sci., 64:1416 1435, 2009.
S. Radl, G. Tryggvason, and J.G. Khinast. Flow and
mass transfer of fully resolved bubbles in nonnewtonian
fluids. AIChE. J., 53:1861 1878, 2007.
S. Radl, A. Koynov, G. Tryggvason, and J.G. Khinast.
DNSbased prediction of the selectivity of fast multi
phase reactions: hydrogenations of nitroarenes. Chem.
Eng. Sci., 63:3279 3291, 2008.
W.J. Rider and D.B. Kothe. Reconstructing volume
tracking. J. Comput. Phys., 2:112 152, 1998.
M. Rieber. Numerische Modellierung der Dynamik
freier Go .. .'l. I. .i in Zweiphasenstroimungen, Ph.D.
Thesis University of Stuttgart. VDI Verlag, Dusseldorf,
2004.
W. Rybczynski. Uber die fortschreitende Bewegung
einer fliissigen Kugel in einem zahen Medium. Bull.
Inst. Acad. Sci. Cracovie, A, pages 40 46, 1911.
B Van Leer. Towards the ultimate conservative differ
ence scheme v. a second order sequel to godunov's
method. Journal of Computational Physics, 1979.
