Chemical Engineering Science 63 (2008) 4580 -- 4595

Contents lists available at ScienceDirect

Chemical Engineering Science

journal homepage: w w w . e l s e v i e r . c o m / l o c a t e / c e s

Computing the added mass of dispersed particles

M. Simcik, M.C. Ruzicka ∗ , J. Drahoš
Department of Multiphase Reactors, Institute of Chemical Process Fundamentals, Czech Academy of Sciences, Rozvojova 135, 16502 Prague, Czech Republic

A R T I C L E I N F O A B S T R A C T

Article history: This paper investigates the possibility of using computational fluid dynamics (CFD) to evaluate the added
Received 25 February 2008 mass coefficient C of dispersed particles, both individually and in clusters. Using the direct numerical
Received in revised form 6 June 2008 simulation (DNS) approach, the volume of fluid (VOF) model was employed, as implemented in Fluent
Accepted 14 June 2008
software, to move the particles. C was calculated from the initial acceleration of a buoyant particle
Available online 21 June 2008
released from rest. The acceleration data were evaluated from very short initial time intervals (typically
Keywords:
∼ 10−5 s), in which buoyancy and added mass forces are predominant. The numerical parameters of the
Added mass code were optimized in cases of known solutions for C so that other situations could be analyzed. Several
Dispersed particles configurations were calculated for their number of particles, shape and spatial arrangement.
Bubbles © 2008 Elsevier Ltd. All rights reserved.
Multiphase flow
CFD

1. Introduction where C is the added mass coefficient. The problem here is to deter-
mine C.
The concept of added mass1 (denoted by AM), one of the classic Thus far the simplest case has been considered, in which the
achievements of theoretical fluid dynamics, contrasts the accelera- particle is of a simple shape (sphere) and translates without rota-
tion, caused by a given force F, of a particle of mass mp in a vac- tion in an unbounded inviscid fluid. In this case, the quantity C is a
uum and in a fluid. In a vacuum Newton's law states that mp a0 = F, scalar and takes the value of 0.5. When the particle has a lower de-
while in a fluid it holds that mp a1 < F where a1 < a0 . The part of the gree of symmetry and is allowed to rotate, we have a second-order
force corresponding to (a0 − a1 ) is spent on the acceleration of the tensor Ĉ with 6 components (6 × 6 matrix). Matrix Ĉ acts on the
surrounding fluid, which, in effect, is equivalent to increasing parti- vector a of (linear and angular) acceleration to produce vector F of
cle mass by the specific amount ma of fluid mass being accelerated (force + torque), Fa = −(V)Ĉa. Consequently, a and F are generally
together with the particle, (mp + ma )a1 = F. The problem is in the misaligned. The components Cij termed inertia coefficients (analo-
determination of ma . gous to the concept of inertia tensor used in rigid body motion) cor-
The product (ma a1 ) is an unsteady inertial force that opposes the respond to the inertia effects acting in different directions of particle
motion of a particle whenever its velocity changes. When a particle translation (forces) and rotation (torques). They are determined by
moves in a stagnant fluid, this force is usually written as the geometry of the problem: particle shape, particle configuration
and presence of boundaries. Since they are not dependent on the
Fa = −ma (dv/dt). (1.1) physical (material) nature of the particles, these coefficients apply,
without discrimination, to bubbles, drops and solids. Tensor Ĉ may
If the fluid itself has a non-zero velocity u, the force is expressed as
be simplified by the presence of symmetry and by the choice of a
Fa = ma (du/dt − dv/dt). (1.2) suitable coordinate system.
The theoretical approach to AM has been described in a number
The added mass ma is the part C of the fluid mass displaced by the of advanced texts on fluid mechanics (e.g., Batchelor, 1967; Birkhoff,
particle of volume V, 1960; Lamb, 1932; Milne-Thomson, 1968; Robertson, 1965; Yih,
1988), as well as in those focusing on fluid inertia forces (Newman,
ma = C V, (1.3) 1977). In addition, analytical and experimental results for the AM
coefficients of various bodies (e.g. rectangular shapes, rotating ellip-
soids, etc.) can be found in Brennen (1982).
∗
Much has also been written about the practical applications of
Corresponding author. Tel.: +420 220 390 299; fax: +420 220 920 661.
E-mail address: [email protected] (M.C. Ruzicka).
AM. As AM is an important issue when particles move in fluids with
1 names are used for ma and (mp + ma ), e.g., added mass, induced mass, inertial unsteady motion, one area of naval research focuses on the strong
mass, and apparent mass, effective mass, hydrodynamic mass, virtual mass. forces of inertia experienced by floating and underwater objects.

0009-2509/$ - see front matter © 2008 Elsevier Ltd. All rights reserved.
doi:10.1016/j.ces.2008.06.011
 M. Simcik et al. / Chemical Engineering Science 63 (2008) 4580 -- 4595 4581

Multiphase fluid dynamics is another area of research in which However, different codes showed less need for AM (Oey et al., 2003),
dispersed particles are also exposed to these forces (Clift et al., 1978; oscillatory plume being reproducible without factoring in AM. Deen
Michaelides, 2006; Wallis, 1990). Here, the simplest case is a sin- et al. (2001) and Zhang et al. (2006) reported that AM forces had
gle dispersed particle translating in an unbounded domain. Two fac- negligible effect on their E–E simulation results for a centrally aerated
tors complicate this situation; the presence of domain boundaries rectangular bubble column (except at the very top and bottom of the
and the presence of other particles. We need C for the following column). Likewise, Tabib (2007) reported that AM had no significant
typical situations: particles approaching a wall, two and more par- influence on the results they obtained using a 3D transient simulation
ticles in certain configurations (doublets, triplets, arrays, etc.), ex- of a cylindrical bubble column.
panding/collapsing particles, oscillating particles, rotating particles, It is, therefore, clear that AM is an important parameter in the
etc. Some of these simpler situations have already been solved, most description of multiphase systems, and can even play the dominant
of them having been approached analytically from the perspective role in certain specific flow situations taking place in technologi-
of potential flow theory. In addition, experiments have been per- cally important processes. AM particularly applies in those situations
formed on the forces acting on oscillating bodies, where the period where the density ratio / p is large, the motion unsteady, the par-
depends on C. Although C has rarely been obtained using compu- ticle shape and configuration liable to change and the presence of
tational fluid dynamics (CFD), we believe that knowledge about C boundaries influential. Thus, in this paper we use CFD to explore the
is scaleable from the microscale (simple situations: single particle, possibility of determining AM in a variety of important situations.
pairwise interactions, etc.) to the macroscale, where the effective in- In each case, the value of C is obtained from the initial acceleration
ertia forces acting within a multiphase mixture are estimated using of the particle(s); this being a relatively simple and efficient way of
various averaging methods. This forms the basis for modelling iner- evaluating C in situations that are well beyond the power of analyt-
tia coupling between various phases. ical approaches and inaccessible to experimentation. Our approach
In recent decades, in the chemical engineering community the was inspired by the work of Niemann and Laurien (2001) and Laurien
concept of AM has become the subject of more research. Primarily, and Niemann (2004), who, to the best of our knowledge, were the
this is due to the surge of interest in CFD, where it is now apparent first to use this method for the calculation of C.
that AM forces play an important role in the dynamics of dispersed
systems (flow patterns, hydrodynamic stability, etc.). AM forces may 2. Evaluation of C
be relevant in cases where the density of the dispersed (particulate)
phase p is comparable to, or even lower than, that of the contin- Consider a buoyant particle, (p < ) released from rest in a
uous (carrying) phase . This applies to liquid–liquid (drops in liq- stagnant fluid. During a very short initial time interval two forces
uids) and solid–liquid (solids in liquids) systems and, in particular, act upon it: the buoyancy force, Fb = g( − p )V, is a constant and
to gas–liquid systems (bubbles in liquids) where, because the inter- accelerates it upwards; the AM force, Fa = ma a = C Va, is propor-
phase density ratio / p is ∼ 103 , almost all the inertia is concen- tional to the initial acceleration ∼ a and opposes the motion. This
trated in the carrying liquid phase. remains approximately true until other forces develop, namely,
AM also plays an important role in the hydrodynamics of bubble those depending on body speed (e.g., drag) and until the possible
columns, which is our primary research interest. As this is such body deformation is negligible. Upon substitution, Newton's law of
a complex problem, few experimental and analytical studies have force
been conducted, and most of the data in this area have been obtained
using multiphase CFD simulations (see e.g., Jakobsen, 2005; Joshi, mp a = Fb − Fa (2.1)
2001; Ranade, 2002; Prosperetti and Tryggvason, 2007). Delnoij
et al. (1997) carried out a Euler–Lagrange simulation of a partially becomes
aerated bubble column and concluded that AM plays an important
(p + C )Va = g( − p )V. (2.2)
role in the vicinity of the gas distributor. While performing lin-
ear stability analysis of a 1D Euler–Euler (E/E) model of a bubble The formula for the coefficient is thus:
column, Leon-Becerril and Line (2001) found that AM has a signifi-
cant effect on homogenous–heterogeneous flow regime transition. C = g /a − p /  (2.3)
Subsequently, Leon-Becerril et al. (2002), carrying out 3D E/E sim-
ulations of a partially aerated rectangular bubble column, found it where g is the reduced gravity, g( − p )/ , and a the initial bubble
necessary to include AM in the model in order to reproduce bub- acceleration determined by CFD.
ble plume oscillations (data published by Becker et al., 1994). The Our buoyant particles, or bodies, were gas bubbles in liquids. CFD
same authors simulated mixing in a fully aerated rectangular col- simulations of free bubble rise were performed for a very short time
umn, again reporting the importance of AM. Monahan et al. (2005) interval, ∼ 10−5 s. Due to the parabolic nature of the hydrodynamic
observed that AM force is important in stabilizing homogeneous equation (Eq. (3.2)), the disturbance produced by the accelerating
regimes, and cannot be neglected in large columns. In the model for body propagates through the whole flow domain so quickly that
homogeneous–heterogeneous regime transition in bubble columns even such a short time interval is sufficient for the development of
developed by Bhole and Joshi (2005) the AM coefficient was a rele- the flow field. The volume of fluid (VOF) model, as implemented
vant control parameter. AM can also affect the values of the terms in Fluent software, was used to perform the simulations. The bub-
describing the production of turbulent kinetic energy due to the ble speed v was obtained using a user defined function to take ten
presence of gas–liquid interfaces (e.g., Mudde and Simonin, 1999; measurements within each time interval. The value of a was ob-
Joshi, 2001). Chahed et al. (2003) introduced turbulence correlations tained by plotting v against time t and fitting it with a straight
related to AM into their expression for interphase force, and found line.
them to be of great importance in the computation of the void frac- A demonstration example is illustrated in Fig. 1. A single spherical
tion in non-homogeneous bubbly flows. Various experimental stud- bubble was placed in a cylindrical domain and started to rise from
ies have also been conducted (e.g., Cai and Wallis, 1993; Kendoush rest (Fig. 1a). The time courses of the individual forces acting on
et al., 2007; Odar and Hamilton, 1964; Takahashi and Endoh, 1992). the bubble are displayed in Fig. 1b. The buoyancy (Fb ) and added
Other results have also shown the ambiguity of the problem of mass (Fa ) forces are given above. The drag force was evaluated as
AM. Mudde and Simonin (1999) concluded that it is necessary to use Fd = mp a − Fb + Fa , where (mp a) is denoted as the inertia force Fi .
AM in E/E simulations of bubble plume with a given CFD program. Fig. 1c gives details of these forces during the short initial interval of
4582 M. Simcik et al. / Chemical Engineering Science 63 (2008) 4580 -- 4595

5.E-06

4.E-06

3.E-06 Fb

force [N]
2.5 Fa
10
2.E-06 Fd
Fi

1.E-06

0.5
0.E+00
2.5
0 0.01 0.02 0.03 0.04 0.05
time [s]

5.0E-06
1.5E-04
4.0E-06 v = 14.39 t
Rxy = 0.9999
1.0E-04
force [N]

3.0E-06 Fb
v [m/s]

Fa
2.0E-06 Fd
Fi 5.0E-05
1.0E-06

0.0E+00 0.0E+00
0.E+00 5.E-06 1.E-05 0.E+00 5.E-06 1.E-05
t [s] t [s]

Fig. 1. Evaluation of C by CFD: a test case. (a) A bubble in cylindrical domain, released from rest, dimensions in mm. (b) Development of the forces on the bubble:
Fb —buoyancy; Fa —added mass; Fd —drag and other possible forces; Fi —`inertia'. Time span 0.05 s. (c) Detail of the forces on a short time interval ∼ 10−5 s. (d) Initial
acceleration of the bubble a = dv/dt and then added mass C by Eq. (2.3).

10−5 s. The AM force was roughly constant over this initial interval 3. Numerical simulations
and the drag was negligible, being less than 3% of the AM force.
The bubble speed is plotted in Fig. 1d, its slope representing the Unless otherwise stated, all simulations were performed using
acceleration a inserted into Eq. (2.3). Fluent 6.1 software (Fluent 6.1 User Guide 2003). To simulate bub-
Over such small time scales the bubble shape remained spherical. ble motion, we used the VOF model in which the evolution of the
For instance, the period of the first mode of bubble shape oscillations bubble interface is described by an advection equation. This model
is T1 ≈ 0. 007 s, which can serve as a deformation time scale. For is designed for two or more immiscible fluids where the position
completeness, the values of certain other parameters are given: Re = of the interface between these fluids is of interest. A single set of
120, We=0. 197, Mo=2. 30×10−11 , Eo=0. 121. Another deformation continuity and momentum equations is shared by the fluids:
time scale T2 can be obtained from We = (V 2 )/(/L). Assuming that
j/ jt + ∇ · (u) = 0. (3.1)
V = L/T2 , we have, T2 = (L3 / We)1/2 , which in our case gives T2 ≈
0. 008 s. Consequently, as both T1 and T2 are larger than 10−5 s, j(u)/ jt + ∇ · (uu) = − ∇p + ∇ · [(∇u + ∇uT )]
bubble deformation is insignificant. + g + S. (3.2)
The tensorial nature of C was not addressed in this study. In-
stead, we focused on the value of C along the dominant direction For each phase, volume fraction q of phase q is introduced and
of motion. Nevertheless, in few cases, the values of the other com- serves as a marker function in the computational cell. In each control
ponents of Cij were also calculated to estimate their magnitude. volume, the volume fractions of all phases combined sum to unity.
For example, the lateral acceleration of the bubble in Fig. 7 was The volume fraction q for a particular phase is: unity when this
found to be 2–3 orders of magnitude smaller than the main ac- phase occupies the whole cell; zero when no phase q is present in the
celeration. And the acceleration perpendicular to gravity sketched cell; 0 < q < 1 at the gas–liquid interface. All variables and properties
in Fig. 8a was less than 6% of the main acceleration (h/r = 1. 1, are shared by the phases and represent volume-averaged values.
 = 45◦ ). The interface motion is inferred indirectly from the motion of
With the arrays, which consisted of multiple bubbles, all bub- phases and is described by the advection equation for the volume
bles were released from rest in stagnant fluid simultaneously and fraction q :
the individual accelerations determined. The values of C obtained in
jq / jt + u · ∇ q = 0. (3.3)
these cases may not correspond to those actually occurring in the
complex flow situations of multiphase mixtures, but are considered Ideally, volume fraction changes discontinuously from 1 to 0 at the
good enough to serve as a first estimate. interface. The discretizations and approximations used in the VOF
 M. Simcik et al. / Chemical Engineering Science 63 (2008) 4580 -- 4595 4583

Table 1
Computational parameters and physical properties of gas/liquid phases used in simulations

Case CPD in bubble region Timestep (s) L (kg m−3 ) G (kg m−3 ) L (kg (m s)−1 ) G (kg (m s)−1 )  (N m−1 )
(dimensionless)

A 80 10−6 –5 × 10−5 1000 100 0.001003 0.001003 0.073

B 80 1.E−06 998.2 1.225 0.01003 1. 7894 × 10−5 0.00073–0.073
C 80 1.E−06 998.2 1.225 0.01003 1. 7894 × 10−5 0.0073
D 40 1.E−06 998.2 1.225 0.01003 1. 7894 × 10−5 0.00073
E 80 1.E−06 998.2 1.225 0.01003 1. 7894 × 10−5 0.00073
F 20–40 per 2 mm 1.E−06 998.2 1.225 0.01003 1. 7894 × 10−5 0.00073
G 40 1.E−06 998.2 1.225 0.001003 1. 7894 × 10−5 0.00073
H 20–80 1.E−06 998.2 1.225 0.001003 1. 7894 × 10−5 0
I 80 (across b) 1.E−06 998.2 1.225 inviscid inviscid 0
Nine different cases were employed (cases from A to I). CPB: number of grid `cells per bubble'; numerical time step; L : liquid density, G : gas density, L : liquid viscosity,
G : gas viscosity, : surface tension.

calculations lead to interface smearing, but this can be minimized ical results. The variance was only greater in a few exceptional cases
using various techniques. (two or more bubbles separated by a short distance; bubbles almost
Surface tension effects are modelled using the continuous surface touching either each other or the wall). Where we were unable to
force model (Brackbill et al., 1992), in which the surface force is rep- find a corresponding value in the literature, we developed a sim-
resented by a smoothly varying volumetric force (S in Eq. (3.2)) that ple correlation based on our numerical data. Errors between these
acts on all fluid elements in the interface region, S = 2q ∇ q /(q + correlations and our data are discussed in relation to the particular
p ). In this model, the surface curvature is q =∇ ·n∗ , the unit normal cases to which they apply.
is n∗ = n/|n| and the normal is n = ∇ q . In the next section, we present our results for C and discuss the
The PRESTO! discretization scheme was used for pressure, effects of the various numerical and other parameters.
the QUICK algorithm for momentum and the PISO algorithm for
pressure–velocity coupling. Further details about the software and 4. Results and discussion
the implementation of these particular procedures can be found in
the reference manual (Fluent 6.1, User Guide 2003). 4.1. Single bubble in an unbounded domain
The following boundary conditions were used. The domain
boundaries were treated either as free-slip walls or as periodic. For 4.1.1. Spherical bubble
the cylindrical coordinates, the axis-of-symmetry boundary condi- The AM for a single spherical bubble in an infinite domain was
tion was used. The difference between free-slip and no-slip condi- obtained for a sufficiently large finite domain, approaching C0 = 0. 5
tions was tested and found to be small (∼ 10−2 relatively) in the as the distance grew from the walls or the other bubbles. This was
short time scales employed (e.g., 1.3% for a 2 mm spherical bubble subsequently used as a reference value. Thus, for example, in the case
in 2.5 mm spherical shell). For the periodic boundaries a pressure of a bubble in a spherical domain (Section 4.6.2), C = 0. 500 when the
drop of p = gHmix (H is domain height, mix is average density domain-to-bubble diameter ratio D/d = 8. Or, for an infinite row of
in the domain) was prescribed between the periodic boundaries. bubbles in a tube (Section 4.5), C =0. 49 for the largest tube diameter
In terms of the initial conditions, bubbles of the prescribed shape and bubble spacing considered.
and configuration were placed into computational domains and all In agreement with the literature, a value of 0.5 was analytically
velocities set to zero. The values of key parameters used in our derived for potential flow as a classical result (e.g., Lamb, 1932).
simulations are shown in Table 1. Laurien and Niemann (2004) reproduced this value by CFD simula-
Computations were performed on a standard PC (Intel Pentium tion, calculating C from the initial acceleration.
4, 2.8 GHz, 1 GB RAM). Depending on the setup, the time needed for
one simulation run ranged from several seconds to 1 h, but in most
4.1.2. Ellipsoidal bubble
cases, was completed within a few minutes. The effect of the domain
The definition sketch for this arrangement is shown in Fig. 2a and
size on the resulting value of C was tested and was generally found
the results are given in Tables 2 and 3 and Fig. 2b. The 2D axisym-
to be small due to the short-range nature of the AM forces. When the
metric cylindrical coordinates were used. The relevant parameters
domain size was enlarged by a factor of 2 or 4, the results typically
are specified under Case F in Table 1.
differed by ∼ 10−2 (relatively). Bubble size was normally set to
The domain size was 20 × 20 mm. Grid independence was tested,
∼ 1 mm. Larger bubbles and hence, domains, are more demanding
with both coarse ( = 0. 1 mm) and fine (0.05 mm) grids giving rea-
on computer resources, while smaller bubbles are generally more
sonable results (Table 2). The effect of surface tension on the fine grid
difficult to calculate because of the problems associated with surface
was investigated, and found to be negligible because the deforma-
forces at large curvatures.
tion time scale was longer than the CFD run (∼ 10−5 s) (see Table 3).
As a result of the symmetry of potential flow, motion in both
The bubble aspect ratio, E = a/b, ranged from 1:4 (oblate) to 4:1 (pro-
the forward and backward directions gives the same value of C.
late). The main result was that C increases as the bubble becomes
Generally only one direction (indicated by the arrow in the definition
more flat (Fig. 2b—marks), recovering C0 at spherical shape E = 1.
sketches) was considered in our calculations. The initial acceleration
The following analytical solution exists for an ellipsoid in infinite
was determined with good precision (see Rxy ≈ 1 in Fig. 1d).
fluid moving along the direction of one of its axes (Lamb, 1932;
Just as with a deterministic machine, each input generates a cor-
Milne-Thomson, 1968, see also Kendoush, 2007):
responding unique output. Each of our results was thus obtained
from a single computational run. We compared our results with  ∞
d
those in the literature. Typically, the variance between our simula- Ce = Z/(2 − Z), Z = abc and
0 (a2 + )k
tion results and the published theoretical solutions was less than 3%
with the simulations generally slightly underpredicting the theoret- k2 = (a2 + )(b2 + )(c2 + ). (4.1.1)
4584 M. Simcik et al. / Chemical Engineering Science 63 (2008) 4580 -- 4595

2.5

2.0

1.5

C [-]
1.0

2b
0.5

2a
0.0
0 1 2 3 4
E = a/b [-]

Fig. 2. Ellipsoidal bubble in unbounded domain (Section 4.1.2). (a) Definition sketch; horizontal semi-axes coincide (b = c), aspect ratio is E = a/b. (b) Result: added mass C
decreases with bubble aspect ratio E; CFD simulation (marks), theoretical solution by Eq. (4.1.1) (line).

Table 2
Ellipsoidal bubble in unbounded domain (Section 4.1.2) r
a (mm) b (mm) C Coarse grid C Fine grid C Analyt.

0.5 2 2.341 2.374 2.374

0.7 1.5 1.193 1.2 1.204
0.9 1.1 0.636 0.635 0.634
1.1 0.9 0.397 0.393 0.392
1.5 0.7 0.196 0.193 0.192 Fig. 3. Spherical cap bubble in unbounded domain (Section 4.1.3). Definition sketch;
2 0.5 0.085 0.083 0.082 bubble of semi-spherical shape.
Effect of computational grid on added mass C, simulations vs. analytical solution by
Eq. (4.1.1) coarse grid  = 0. 1 mm, fine grid  = 0. 05 mm (see Fig. 2a).

Table 4
Spherical cap bubble in unbounded domain (Section 4.1.3)
Table 3
Ellipsoidal bubble in unbounded domain (Section 4.1.2) C (dimensionless)
 (N m−1 ) C (dimensionless) CFD 1 0.852
CFD 2 1.24
0.00073 2.374
CFD 3 1.22
0 2.363
Kendoush, 2003 2.04
Analytical solution 2.374
Three different CFD methods (CFD1,2,3) for calculating added mass C. Approximate
Effect of surface tension  on added mass C, a = 0. 5 mm, b = 0. 5 mm, fine grid
solution by Kendoush (2003) (see text and Fig. 3).
 = 0. 05 mm, analytical solution by Eq. (4.1.1) (see Fig. 2a).

Lamb (1932) gives the closed formulas: Table 5

Spherical cap bubble in unbounded domain (Section 4.1.3)

E cos−1 E − 1 − E2 CPD (dimensionless) C (dimensionless)
C= 
1 − E2 − E cos−1 E
E2 20 0.824
oblate ellipsoid (E < 1) (4.1.2) 40 0.839
80 0.852
and Effect of numerical grid on added mass C (method CFD1) (see text and Fig. 3).
 
E ln(E + E2 − 1) − E2 − 1
C=  
E2 E2 − 1 − E ln(E + E2 − 1) Table 6
prolate ellipsoid (E > 1) (4.1.3) The 2D bubble: effect of computational grid (CPB) on value of C (see Section 4.1.4)

C (dimensionless)
Our calculations agreed with this analytical solution (Fig. 2b—marks
vs. line). CPD = 20 0.983
CPD = 40 0.998
Analytical .
4.1.3. Spherical cap bubble
The definition sketch is in Fig. 3 and the results are in Tables 4
and 5.
Three simulations were performed for a semi-spherical bubble case H in Table 1) and gave C = 0. 852 (Table 4). The results were
with a radius r of 1 mm. CFD1 involved using VOF to accelerate a relatively grid independent (Table 5). Because a large surface force S
fluid particle (domain size 8 × 8 mm; parameters corresponding to (last term in Eq. (3.2)) was generated by the curvature at the bubble
 M. Simcik et al. / Chemical Engineering Science 63 (2008) 4580 -- 4595 4585

0.70

0.65

0.60
r

C [-]
0.55

h 0.50

wall 0.45
1 2 3 4 5
h/r

Fig. 4. Spherical bubble perpendicular to wall (Section 4.2.1.1). (a) Definition sketch; bubble size fixed to r = 1 mm, wall distance h varies. (b) Result: added mass C decreases
with wall distance h; CFD simulation (marks), theoretical solution by Eq. (4.2.1) (line).

edges, a very low surface tension had to be used ( = 0) to secure Table 7

Spherical bubble perpendicular to wall (Section 4.2.1.1)
grid independence.
CFD2 was a stationary rigid semi-spherical bubble measured in CBD (dimensionless) h/r = 1. 1 h/r = 1. 5 h/r = 2. 0 h/r = 5. 0
an accelerating flow field with a domain size of 16 × 32 mm. The AM 20 0.614 0.536
force was calculated, yielding C = 1. 24. 40 0.624 0.545
CFD3 was a rigid body accelerated from rest in stagnant fluid 80 0.634 0.553 0.52 0.499
Analytical 0.641 0.556 0.523 0.502
(domain size 16×32 mm). The AM force was then calculated, yielding
C = 1. 22. Effect of grid (CPB) and wall distance (h) on value of C, analytical solution by Eq.
(4.2.1) (see Fig. 4a).
Although we are not currently able to explain the variance
between these three results, the fact that a difference exists si-
multaneously serves as both a warning and a challenge for further
studies. All three methods yielded the same result, C = 0. 5, for a
fully spherical particle.
For an approximate inviscid solution (potential flow ahead of the
bubble, Hill's vortex behind it), Kendoush (2003) reported that AM
increases rapidly with bubble flatness: C = ∞ for a flat disk; C = 2. 04 2b
for a half-sphere.

4.1.4. 2D bubble 2a
Simulations for a 2D circular bubble (infinite cylinder) in un-
bounded fluid were performed using 2D Cartesian coordinates. A
h
single bubble (d = 2 mm) was placed in the middle of a 20 × 20 mm
square domain (effectively `unbounded'), the parameters corre-
sponding to case D in Table 1. Two grids were tested and C = 1
was obtained for each. These results agree with various analytical
Fig. 5. Ellipsoidal bubble perpendicular to wall (Section 4.2.1.1). Definition sketch:
solutions (e.g. Keulegan and Carpenter, 1958; Newman, 1977) (see
horizontal semi-axes coincide (b = c), aspect ratio is E = a/b, wall distance h varies.
Table 6).

4.2. Single bubble and wall radii), which is typical for the inertial force. At h/r = 3, C differs from
C0 by only 1.5%.
4.2.1. Perpendicular motion Milne-Thompson's (1968) approximate analytical solution (first
4.2.1.1. Spherical bubble. The definition sketch is in Fig. 4a, and the order method of images) states that
results are in Table 7 and Fig. 4b. The 2D axisymmetric cylindrical
coordinates were used. The relevant parameters correspond to case Cpe = C0 (1 + (3/8)r3 /h3 ). (4.2.1)
E in Table 1.
The domain size was 20 × 20 mm, the bubble radius r = 1 mm and Our calculations agreed with this leading order solution
the wall distance h/r ranged from 1.1 to 5.0. A relatively fine grid was (Fig. 4b—marks vs. line).
used (Table 7). The surface tension was maintained at the low level This case is equivalent to two bubbles moving along a line of
of 0.00073 N/m. The effect of the domain size on C was small (errors centres due to the plane of symmetry. Further terms in the two-
less than 1%). The main result observed was that the value of C fell bubble expansion have been published by Kok (1993), see also
quickly as the wall distance grew (Fig. 4b—marks), the influence of Legendre et al. (2005): 1 + (3/8)r3 /h3 + (3/64)r6 /h6 + (9/256)r8 /h8 +
the wall only being significant at a very short distance (few bubble (3/512)r9 /h9 +, . . . , which applies also to Cpe .
4586 M. Simcik et al. / Chemical Engineering Science 63 (2008) 4580 -- 4595

1.7 5

1.6
4
1.5
3
C [-] 1.4

C [-]
1.3 2
1.2
1
1.1

1.0 0
1 3 5 7 9 0.2 0.3 0.4 0.5 0.6 0.7
h/a or h/z [-] E = a/b [-]

Fig. 6. Ellipsoidal bubble perpendicular to wall (Section 4.2.1.2). Result: added mass C decreases with both wall distance h and bubble aspect ratio E; CFD simulation (marks),
empirical formula by Eq. (4.2.2) (line). (a) h is scaled either by a (white marks) or by z (black marks), E = 0. 5. (b) h is scaled both on a (dashed line) and on z (full line),
wall distance is h/a = 2.

4.2.1.2. Ellipsoidal bubble. Definition sketch is in Fig. 5 and the re-

sults are given in Fig. 6. In this case Fluent ver. 6.3 was used. 2D
axisymmetric cylindrical coordinates were used, the parameters cor-
responding to case I in Table 1.
The domain size was 16 × 16 mm. One semi-axis was fixed (b =
1 mm), while the other semi-axis (a) varied in length. The bubble
distance h/a ranged from 1.1 to 10 at fixed aspect ratio E = a/b = 1/2. h
E ranged from 1/5 to 2/3 at a fixed distance h/a = 2. The effect of the r
domain size on C was small. The main result being that the value of
C fell with both h and E (Fig. 6—marks).
As we were unable to find an analytical solution for this case
in the literature, we suggest the following empirical formula, which
preserves the structure of Eq. (4.2.1): Wall

C = Ce (1 + (3/8)(z3 /h3 )) (4.2.2) Fig. 7. Spherical bubble parallel to wall (Section 4.2.2.2). Definition sketch; bubble
size fixed to r = 1 mm, wall distance fixed at h/r = 1. 1 (close to wall).
In this formula, the spherical coefficient C0 is replaced by the ellip-
soidal one Ce given by Eq. (4.1.1) in which r was replaced by another
length scale defined as z = bh/(h − a + b). This resulted in a good data of approach  ranged from 0◦ (perpendicular) to 90◦ (parallel) by
fit with Eq. (4.2.2). Note that z approaches a on contact (h = a) and changing the direction of gravity. The computational grid consisted
recovers Eq. (4.2.1) at a = b = r for spherical shapes. For a flat bubble of several zones, the finest grid being near the bubble and the coarser
(E = 1/5, h/a = 2) and close spacing (h/a = 1. 1, E = 1/2), Eq. (4.2.2) grids further away (Table 8). The effect of the domain size on C
fits the data with an error of about 10%, which rapidly fell as both h was small. The main result observed was the dependence of C on 
and E increased. (Fig. 8b).
As we could not find an analytical solution for this case in the
4.2.2. Parallel motion of spherical bubble literature, we propose the following empirical formula:
Definition sketch is in Fig. 7. The 3D Cartesian simulations were
C = A Cos(2) + B, A = (Cpe − Cpa )/2,
performed, the parameters corresponding to case G in Table 1.
The wall distance was fixed at h/r =1. 1. The bubble radius was r = B = (Cpe + Cpa )/2, (4.2.4)
1 mm. The resulting value of C was 0.558, which is in good agreement where Cpe and Cpa are the limits of perpendicular motion ( = 0◦ )
with Milne-Thompson's (1968) analytical solution, where and parallel motion ( =90◦ ) given by Eqs. (4.2.1) and (4.2.3), respec-
tively. At these limits, the numerical results agree with the analytical
Cpa = C0 (1 + (3/16)(r/h)3 ) (4.2.3)
solution within 3%. This case is equivalent to two bubbles mirrored
yields C = 0. 570. by the plane of symmetry.
This case is equivalent to two bubbles moving perpendicular to
a line of centres, due to symmetry. 4.3. Two spherical bubbles in unbounded domain

4.2.3. Oblique motion of spherical bubble Definition sketch is in Fig. 9, results in Figs. 10–12. The 3D Carte-
Definition sketch is in Fig. 8a and the results are given in Fig. 8b. sian simulations were performed, the parameters corresponding to
The 3D Cartesian simulations were performed, the parameters cor- case G in Table 1.
responding to case G in Table 1. Two spherical equal-sized bubbles of fixed diameter (d = 2 mm)
A cubed domain was used (16 × 16 × 16 mm3 ). The bubble, with were positioned inside a cubed domain (27. 2 × 27. 2 × 27. 2 mm3 ).
a fixed radius r = 1 mm, was positioned near one of the cube walls. The bubble distance D/d ranged from 1.1 to 5.0 and the angle of
The wall distance was fixed at h/r = 1. 1 (close spacing). The angle alignment  ranged from 0◦ (vertical) to 90◦ (horizontal). To save
 M. Simcik et al. / Chemical Engineering Science 63 (2008) 4580 -- 4595 4587

0.70

r 0.65
h

C [-]
0.60

0.55

0.50
0 10 20 30 40 50 60 70 80 90
α [°]

Fig. 8. Spherical bubble at inclined wall (Section 4.2.3). (a) Definition sketch; bubble size fixed to r = 1 mm, wall distance fixed to h/r = 1. 1 (close to wall), approach angle
 varies. (b) Result: added mass C decreases with ; CFD simulation (black marks), empirical formula by Eq. (4.2.4) (line), two limits  = 0◦ and 90◦ correspond to Cpe and
Cpa by Eqs. 4.2.1) and (4.2.3) (white marks).

Table 8 4.3.1. Vertical alignment

Computational grid (CPB) for 3D simulations near bubble (Section 4.3)
The AM of two bubbles in-line ( = 0◦ ) decreases because this
Range x CBD (dimensionless) doublet is more streamlined than a single sphere (Fig. 10a). However,
due to the short-range character of inertia forces, this decrease is
x < 1. 5r 40
1. 5r < x < 2. 0r 20 only evident with close spacing (a few radii). The following empirical
2r < x < 2. 4r 10 formula was developed to fit our CFD data:
2. 4r < x < 3. 0r 5
x > 3. 0r 2.5 Cv = 0. 35 + 0. 15(D/d − 1)/(D/d − 0. 546), (4.3.1)

giving 0.35 for touching bubbles and 0.5 where there was no contact.
Our results are in relatively good agreement with those of Helfinstine
(1974), who used the method of images (Fig. 10a). Note that this
arrangement is different from a bubble perpendicular to the wall in
Section 4.2.1 (cf. Figs. 4b and 10a).

4.3.2. Horizontal alignment

The AM of two bubbles side by side ( = 90◦ ) is greater (Fig. 10b)
because this doublet is less streamlined than a single sphere. Again,
this effect only applies at short range. At a D/d of about 2–3, AM
D approaches a constant value, but slightly underpredicts the reference
value C0 expected for wide spacing. The following empirical formula
was developed to fit our CFD data:
α
d Ch = 0. 61 − 0. 11(D/d − 1)/(D/d − 0. 862), (4.3.2)

giving 0.61 for touching bubbles and 0.5 where there was no con-
tact. Again, our results compare relatively well with the approxi-
mate solution of Helfinstine (1974), who used the method of images
(Fig. 10b).
Note that this arrangement is equivalent to a bubble rising along
a wall (Section 4.2.2). Consequently, the case where h/r = 1. 1 and
Fig. 9. Two bubbles in general position (Section 4.3). Definition sketch; bubble size C = 0. 558 in Section (4.2.2) corresponds to D/d = 1. 1 and C = 0. 558
fixed to d = 2 mm; bubble distance D and alignement angle  vary.
in Fig. 10b (black marks). Ideally, Cpa = Ch and therefore, the
Eqs. (4.2.3) and (4.3.2) are equivalent. In our simulations, we en-
countered some discrepancies, the magnitude of which can be used
processing time and power, the simulations were initially performed to assess possible errors connected with these types of calculation.
for one bubble only in order to determine a suitable grid and soft-
ware configuration. A fine grid was used up to 1.5 bubble radii 4.3.3. General position
from the bubble surface. Beyond this, coarser grids were used (see The effect of the inclination angle  is demonstrated in Fig. 11
Table 8). Forty CPB ensured acceptable grid independence. The effect (marks) for fixed bubble spacing D/d = 1. 5. The AM gradually in-
of the domain size on C was small. Due to the fore–aft symmetry of creases from the limit of the vertical pair (C ≈ 0. 44) to the limit of
potential flow, the same value of C was obtained for both bubbles. the horizontal pair (C ≈ 0. 52), passing through the reference value
4588 M. Simcik et al. / Chemical Engineering Science 63 (2008) 4580 -- 4595

0.50
0.60

0.45

0.56

C [-]
C [-]

0.40

0.52
0.35

0.30 0.48
1 2 3 4 5 1 2 3 4 5
D/d [-] D/d [-]

Fig. 10. Two bubbles in vertical and horizontal position (Sections 4.3.1 and 4.3.2). (a) In-line arrangement ( = 0◦ ), added mass C increases with distance D. (b) Side by side
arrangement ( = 90◦ ), added mass C decreases with distance D. CFD simulation (black marks), empirical formulas by Eqs. (4.3.1) and (4.3.2) (full lines), Helfinstine (1974)
theoretical results (white marks), a bubble parallel to wall by Eq. (4.2.3) (dashed line) and C = 0. 558 (white square) from Section (4.2.2).

4.4. Three spherical bubbles

0.54
Definition sketch is shown in Fig. 13 and the results in Fig. 14.
The 3D Cartesian simulations were performed with parameter values
corresponding to case G in Table 1.
The computational domain was a 27.2 mm cube. Three bubbles
(d = 2 mm) forming an equilateral triangle were placed in the cen-
0.50 tre of the domain in one of two different configurations, horizon-
tal (Fig. 13a) or vertical (Fig. 13b). The horizontal bubbles behaved
C [-]

identically due to symmetry (finite domain shape was assumed to

have no effect), whereas the vertical cluster consisted of one lead-
ing bubble (1) and two identical trailing bubbles (2, 3). The bub-
ble distance D/d ranged from 1.1 to 2.0. The computational grid
0.46 was finer closer to the bubbles and coarser further away from (see
Table 8). In both cases, the effect of the domain size on C was
small.
In the horizontal triangle, each bubble had the same value of C,
which was larger than both C0 and Ch (Fig. 14a). The flat triangle
was much less streamlined. This result clearly demonstrates that
0.42
the presence of other bubbles can strongly affect inertial coupling
0 10 20 30 40 50 60 70 80 90
based only on pairwise interactions. As such multiple interactions
α [°] are difficult to obtain both analytically and experimentally, and the
CFD approach can be of great help.
Fig. 11. Two bubbles in general position: effect of alignment angle  (Section 4.3.3).
Added mass C increases with . CFD simulation (black marks), empirical formula by
The results obtained for the vertical triangle can be summarized
Eq. (4.3.3) (dashed line), Kok (1993) theoretical result (full line). D/d = 1. 5. as C1 < C2 = C3 < C0 (Fig. 14b—marks). The value of C1 is close to that
of Cv obtained for the vertical doublet. The C2,3 values fall between
the values of Cv and Ch obtained for the vertical and horizontal
C0 = 0. 5 at about  ≈ 60◦ . The following empirical formula was de- doublets, and were close to a doublet inclined at  ≈ 50◦ . Part of the
veloped to fit our CFD data: vertical triangle (bubbles 1 and 2) can be considered as a doublet
inclined at  ≈ 30◦ . For comparison purposes, Fig. 14b also plots
C = A sin(2( − 45◦ )) + B, A = 12 (Ch − Cv ), the corresponding line obtained using Eq. (4.3.3), this line passing
between the lines for C1 and C2 .
B = 12 (Cv + Ch ), (4.3.3) We believe that horizontal and vertical limits on bubbles in tri-
angular configurations set bounds on the value of C between these
where Cv and Ch , the aforementioned vertical and horizontal coef-
two limits. Our results for three bubbles also indicate the relevance
ficients, are dependent on distance D/d (Fig. 11—dashed line). Our
of multiple inertial coupling in bubbly mixtures.
results are in relatively good agreement with the approximate so-
lution published by Kok (1993), who considered the first four terms
in the method of images (Fig. 11—full line). The differences between 4.5. The 1D array of spherical bubbles in a tube
our results and Kok's are more pronounced at low  and D/d. For
completeness, Fig. 12 shows the influence of both D/d and  on C. Definition sketch is shown in Fig. 15 and the results are given, in
Their combined effect on AM can be clearly seen, causing both a de- Tables 9–12 and Fig. 16. The 2D axisymmetric cylindrical coordinates
crease and an increase in the value of C. Similar behaviour has been were used with the parameters corresponding to case C in Table 1
reported for drag and lift (Legendre et al., 2003). (main results).
 M. Simcik et al. / Chemical Engineering Science 63 (2008) 4580 -- 4595 4589

0.60 0.60

0.50 0.50

C [-]
C [-] alpha = 0°
alpha = 15°
alpha = 30°
0.40 alpha = 45° 0.40 D/d = 1.1
alpha = 60° D/d = 1.5
alpha = 75° D/d = 2
alpha = 90° D/d = 5
0.5 0.5
0.30 0.30
1 2 3 4 5 0 10 20 30 40 50 60 70 80 90
D/d [-] α [°]

Fig. 12. Two bubbles in general position: combined effect of bubble spacing D and alignment angle  (Section 4.3.3). (a) Effect of D, parameterized by . (b) Effect of ,
parameterized by D. CFD simulation (marks).

2 3
D D

Fig. 13. Three bubbles in triangle (Section 4.4). Definition sketch; bubble size fixed to d = 2 mm, bubble spacing D varies. (a) Horizontal triangle; bubbles equivalent. (b)
Vertical triangle; leading bubble (1), trailing bubbles (2,3).

0.60
0.65

0.60
0.50
C [-]
C [-]

0.55

0.40 C1
C2 = C3
0.50 CV
CH
C30
C50
C=0.5
0.45 0.30
1 1.5 2 1 1.5 2
D/d [-] D/d [-]

Fig. 14. Three bubbles in triangle (Section 4.4). (a) Horizontal triplet; added mass C decreases with bubble spacing D (marks). For comparison Ch is shown for horizontal
doublet (full line, Eq. (4.3.2)). (b) Vertical triplet; added mass C increases with bubble spacing D (marks), C1 < C2 = C3 . For comparison, other results are shown by lines; from
below: vertical doublet (bold full line, Cv by Eq. (4.3.1)), doublet aligned at  = 30◦ and 50◦ (grey line and thin full line, C30,50 by Eq. (4.3.3)), horizontal doublet (dashed
line, Ch by Eq. (4.3.2)).

The free-slip boundary condition was used at the tube wall; the to 20; H being correspondingly enlarged as the number of bubbles
periodic boundary condition along the tube axis. The bubble radius in the domain increased (N = 1, 2, etc.). Four different grids were
r was fixed at 1 mm. The tube radius R/r ranged from 1.1 to 20. By used within a range of 20–160 CPB. The grid test showed that 80 CPB
altering the value of H, the bubble spacing H/r=2H/d ranged from 1.1 was acceptable (Table 9). Changing the viscosity (by four orders)
4590 M. Simcik et al. / Chemical Engineering Science 63 (2008) 4580 -- 4595

Table 12
Bubble array in tube (Section 4.5)

CBD (dimensionless) No. of bubbles (N)

1 2 10
R 20 0.463 0.476 0.477
40 0.501 0.505 0.505
80 0.508 0.508 0.508
160 0.509 0.51

Effect of number of bubbles N in periodic domain on resulting value of C (H = 2 mm,

2H R = 1. 63 mm, r = 1 mm, CPB = 40) (see Fig. 15).

C decreases. To provide a comparison, the data for the unbounded

vertical doublet Cv (from the simulation depicted in Fig. 10a) are also
shown (white marks). At wider spacing H, the rise in the value of C
above 0.5 is due to the confining tube wall. At fixed bubble spacing
H, C decreases with tube size R, which is similar to a bubble moving
r
parallel to the wall (Fig. 16b—black marks). In the tube, the value of
C is higher due to the presence of the wall. At high values of both H
Fig. 15. One-dimensional bubble array in tube (Section 4.5). Definition sketch; bubble
size fixed to r = 1 mm, bubble spacing H and tube radius R vary. and R, the reference value C0 is recovered.
The important consequence of these results is that C is dependent
on the aspect ratio of the domain. For narrow and tall domains, the
wall effect prevails and AM increases (C > 0. 5), whereas for wide and
Table 9 short domains, the array effect prevails and AM decreases (C < 0. 5).
Bubble array in tube (Section 4.5) Obviously, these effects mutually cancel themselves out at a certain
CBD (dimensionless) H (mm) domain proportion (C ≈ 0. 5).
H = 1. 1 1.5 2 3 6 7.5 12 15 20
Our simulation results were compared with Cai and Wallis (1992)
approximate solution (Fig. 16—black lines). For potential flow past
20 0.141 0.308 0.431 0.538 0.655 0.716 0.735 0.748 0.773
an infinite 1D array of equidistant spheres in a circular open-ended
40 0.149 0.333 0.45 0.563 0.686 0.754 0.769 0.782 0.806
80 0.154 0.338 0.451 0.566 0.692 0.761 0.776 0.787 0.812 tube, they obtained
160 0.155 0.339 0.452 0.566 0.692 0.761
Cai and Wallis 0.228 0.34 0.434 0.548 0.701 C = ( − 1)HR2 /((2/3) ) − 1. (4.5.1)
Grid refinement study: effect of CPB and bubble spacing H on resulting value of C
For 1 < H/r < 6 and 1 < R/r < 2. 5, the parameter was correlated as
(r = 1 mm, R = 1. 63 mm) (see Fig. 15).

= 1 + a/(H/r)b ,

Table 10
a(R/r) = exp(1. 3916 − 7. 2445 ln(R/r) + 7. 049(ln(R/r))2
Bubble array in tube (Section 4.5) − 3. 3251(ln(R/r))3 ,
L (kg/m s) C (dimensionless)
−1
b(R/r) = 0. 99832(R/r)−0.084924 . (4.5.2)
1. 003 × 10 0.509
1. 003 × 10−2 0.502 Good agreement was found between their results and ours. An array
1. 003 × 10−3 0.500
in a tube can be used as an `elementary cell' in dispersed flow, possi-
1. 003 × 10−4 0.501
1. 003 × 10−5 0.501 bly approximating to an ensemble of `bubble chains' with a spacing
of 2H in the vertical direction and of 2R in the horizontal direction.
Effect of liquid viscosity L on resulting value of C (H = 2 mm, R = 1. 63 mm, r = 1 mm,
CPB = 40) (see Fig. 15). Qualitatively similar trends also apply to arrays in channels with
non-circular cross sections.

4.6. The 3D bubble arrays

Table 11
Bubble array in tube (Section 4.5)
If real dispersions are to be correctly modelled in general and in-
 (N/m) C (dimensionless)
terphase inertial coupling to be determined in particular, the effects
7. 3 × 10−2 0.501 of AM must be clearly understood. An important part of this involves
7. 3 × 10−3 0.501 understanding how C is dependent on bubble volume fraction and
7. 3 × 10−4 0.502
configuration. Certain important details of bubble arrangements are
7. 3 × 10−5 0.501
not always resolved, often being considered `random' or `regular'.
Effect of liquid surface tension  on resulting value of C (H = 2 mm, R = 1. 63 mm,
Typically AM increases with bubble concentration c. Some results
r = 1 mm, CPB = 40) (see Fig. 15).
have been obtained using simple geometric concepts. For example,
Zuber (1964) applied the classical formula for a spherical body in a
spherical domain (e.g., Lamb, 1932),
and surface tension (by three orders) had virtually no effect on the
C = C0 (1 + 2(d/D)3 )/(1 − (d/D)3 ), (4.6.1)
computation results (Tables 10 and 11). As the number of bubbles
in the domain did not affect the results (Table 12), N = 1 was used and considering the domain as an elementary cell in the dispersion
in all calculations. of concentration c = (d/D)3 , obtained
At fixed tube size R, C increases with bubble spacing H
(Fig. 16a—black marks). As a dense array becomes more streamlined, C = C0 (1 + 2c)/(1 − c). (4.6.2)
 M. Simcik et al. / Chemical Engineering Science 63 (2008) 4580 -- 4595 4591

3.0
0.8

2.5
0.6
2.0

C [-]
C [-] 1.5
0.4

1.0
0.2
0.5

0 0.0
1 2 3 4 5 6 1.0 1.5 2.0 2.5
H/r [-] R/r [-]

Fig. 16. One-dimensional bubble array in tube (Section 4.5). (a) Added mass C increases with bubble spacing H (black marks), R/r = 1. 63; theoretical solution of Cai and
Wallis (1993) (line), for comparison vertical doublet Cv by Eq. (4.3.1) (white marks). (b) Added mass C decreases with tube radius R (black marks), H/r = 3(•) and 6();
theoretical solution of Cai and Wallis (1993) (full lines); for comparison bubble parallel to wall Cpa by Eq. (4.2.3) (dashed line).

0.6

0.4
C [-]

d
Cg2
Cg3
0.2 Cg1
CT
g1 g2 Eq. 4.6.7
Eq. 4.6.8
g3 Eq. 4.6.9
D
0
0 10 20 30 40
c [%]

Fig. 17. Bubble in periodic box (Section 4.6.1). (a) Definition sketch; bubble size fixed to d = 2 mm, domain size D varies, three directions of gravity used g1, 2, 3. Situation
corresponds to 3D regular array with bubble concentration c = ( /6)(d/D)3 . (b) Result: added mass C decreases with bubble concentration c due to the `wide' domain (cf.
Figs. 17 and 20). Three direction of gravity: Cg1 , Cg2 , Cg3 (marks, coincide); empirical formula by Eq. (4.6.9) (line). For comparison: CT , bubble in tube, same distance and
cross-section (Section 4.5, Eq. (4.5.1)); Cai and Wallis (1994) result for spherical shell and pressure-compliance condition (Eqs. (4.6.7) and (4.6.8)).

This cell model is relatively robust, works for all values of c and gives Re, they evaluated C from the increase in bubble speed due to in-
reasonable results. At the dilute limit, the expansion of Eq. (4.6.1) creased gravity, and obtained k=4 as a compromise for most of their
yields data. Laurien and Niemann (2004) and Niemann and Laurien (2001)
employed the DNS approach and calculated C from initial bubble ac-
C = C0 (1 + 3c) + O(c2 ). (4.6.3) celeration. For a single bubble, a value of 0.5 was correctly obtained
Other results have been obtained on theoretical grounds using de- for C0 . Using a periodic box, for a regular array they found that
manding analytical efforts based on averaging techniques or energy
C = C0 (1 + 3. 26c + 7. 70c2 ), (4.6.5)
concepts. These formulas usually take the form

C = C0 (1 + kc) + O(higher terms), (4.6.4) where k = 3. 26, which is in agreement with the analytical results.
Alternatively, the generalized cell model used by Cai and Wallis
where the value of k is typically close to 3, but varies according to (1994),
the particular approach used, the specific assumptions made and the
bubble configurations, velocity distributions, etc. employed. For in- C = [(1 − c) + k(0. 5 + c)]/[(2 + c) + k(1 − c)], (4.6.6)
stance, Wijngaarden (1976) obtained k = 2. 78 while both Kok (1988)
and Biesheuvel and Spoelstra (1989) reported k = 3. 32. Spelt and predicts both an increase and decrease in AM with bubble concentra-
Sangani (1998) used the averaged equations of motion of bubbly tion c depending on the value of parameter k. For k = ∞, their model
mixtures at large Re and small We, solved them numerically, and recovers the increasing Zuber result (Eq. (4.6.2)). For k = 0, it gives a
reported k = 3. Sankaranarayanan et al. (2002) used a two-continua decreasing result:
model within the Lattice Boltzmann approach to simulate bubble ar-
rays in a periodic box. For spherical and distorted bubbles at finite C = (1 − c)/(2 + c), (4.6.7)
4592 M. Simcik et al. / Chemical Engineering Science 63 (2008) 4580 -- 4595

with the dilute limit corresponds to a movable or penetrable outer shell (ideal pressure
release surface, Eq. (4.6.7)). Parameter k can be seen as the ratio of
C = C0 (1 − 1. 5c) + O(c2 ). (4.6.8) (outer)/(inner) fluid density.

The increase in AM at large k corresponds to a fixed outer shell

4.6.1. Spherical bubble in a periodic box
(stationary rigid wall, Eq. (4.6.2)). The decrease in AM at low k
Definition sketch is shown in Fig. 17a and the results in Fig. 17b.
The 3D Cartesian simulations were performed and the parameters
corresponding to case G is given in Table 1.
The bubble size was fixed at d = 2 mm. The periodic box size
D/d ranged from 1.1 to 6 by altering the value of D. The finer grid
d (CPB = 40) was closer to the bubble, the coarser grid further away
(Table 8). Three directions of bubble motion (gravity) were tested:
D box face (g1), edge (g2) and corner (g3). The direction of gravity had
no effect on the CFD results (Fig. 17b). This situation corresponds to

6
Fig. 18. Spherical bubble in spherical domain (Section 4.6.2). Definition sketch;
bubble size fixed to d = 2 mm, domain size D varies. Situation corresponds to 3D
random array with bubble concentration c = (d/D)3 .
5

Table 13
Spherical bubble in spherical domain (Section 4.6.2)
4
CPD D (mm)

2.2 2.5 3 4 5 6 8 16
C [-]

20 4.805 2.04 1.061 0.667 0.595 0.567 0.504 0.491 3

40 4.98 2.088 1.144 0.704 0.598 0.546
80 5.094 2.098 1.118 0.702 0.599 0.552 0.52 0.5
160 5.139 1.989 1.126 0.599
Analytical 5.032 2.074 1.132 0.714 0.603 0.558 0.524 0.503
2
Grid refinement study: effect of CPB and domain size D on resulting value of C,
analytical solution by Eq. (4.6.1) (see Fig. 18).

1
Table 14
Spherical bubble in spherical domain (Section 4.6.2)

L (kg/m s) C (dimensionless) 0
0.01003 0.702
0.0 0.2 0.4 0.6 0.8
0.001 0.7 c [-]
0.0001 0.7
Inviscid 0.698
Fig. 20. Spherical bubble in spherical domain (Section 4.6.2). Result: added mass C
Analytical solution 0.714
increases with bubble concentration c (marks) (cf. Figs. 20 and 17). For comparison:
Effect of liquid viscosity L on resulting value of C (D = 4 mm, d = 2 mm, CPB = 80), Zuber cell model by Eq. (4.6.2) (bold line), its linearization by Eq. (4.6.3) (thin line),
inviscid solution for L = 0, analytical solution by Eq. (4.6.1). (see Fig. 18). Laurien and Niemann result for periodic box by Eq. (4.6.5) (dashed line).

1.2 1.2

1.0 1.0

0.8 0.8
C [-]
C [-]

0.6 0.6

0.4 0.4
surf. t.= 0.073N/m D = 4 mm
surf. t.= 0.0073N/m D = 8mm
0.2 surf. t.= 0.00073N/m 0.2 D = 40 mm
Analytical solution Analytical solution
0.0 0.0
0.0 0.1 0.1 0.00 0.05 0.10
δ [mm] δ [mm]

Fig. 19. Spherical bubble in spherical domain (Section 4.6.2). Numerical tests. (a) Effect of grid size  and surface tension  (d = 2 mm, D = 2. 2 mm). (b) Effect of grid size 
and domain size D ( = 0. 073 N/m). Analytical solution by Eq. (4.6.1) (line).
 M. Simcik et al. / Chemical Engineering Science 63 (2008) 4580 -- 4595 4593

3.0

R
2.5

2.0
2b

C [-]
1.5
2a
1.0

0.5

0.0
0.0 0.1 0.2 0.3 0.4
c [-]

Fig. 21. Ellipsoidal bubble in spherical domain (Section 4.6.2). (a) Definition sketch; bubble size fixed to a = 0. 5 mm, b = 1 mm, domain size R varies. Situation corresponds
to 3D random array with bubble concentration c = (ab2 /D3 ). (b) Result: added mass C increases with bubble concentration c (marks), empirical formula by Eq. (4.6.10) (full
line). For comparison: Zuber cell model for spherical bubble by Eq. (4.6.2) (dashed line).

an infinite 3D regular bubble array where c = ( /6)(d/D)3 and ranged centration (Fig. 20), thereby, as expected, conforming to the trend
from 0.002 to 0.393. The periodic boundary condition mimicked the to increase. Agreement with the Zuber cell model was good.
pressure release surface, rather than the fixed wall.
Contrary to general expectation, we found that C decreases with 4.6.2.2. Ellipsoidal bubble. Definition sketch and results are pre-
bubble concentration (Fig. 17b). The reason for this is as follows. sented in Fig. 21. The 2D axisymmetric coordinates were used with
Firstly, the periodic box is close to pressure release case, so that, in the parameters corresponding to case I in Table 1.
principle, AM can decrease (see Eq. (4.6.7)). Secondly, the periodic The bubble size was fixed: a = 0. 5 mm; b = 1 mm; aspect E =
box H × D × D is equivalent to a 1D bubble array of spacing H in a a/b = 0. 5. By altering R the shell size R/a ranged from 1.5 to 8. The
rectangular channel of cross-section D×D. We know that C increases numerical parameters were the same as in the case of the spherical
with H (Fig. 16a) and decreases with D (Fig. 16b). Thus for a tall bubble described above. These simulations were performed using
domain C will increase, while for a wide domain it will decrease. Fluent ver. 6.3. This situation corresponds to that of an infinite 3D
With our isometric domain H = D, the decreasing tendency prevails. random bubble array where c = ab2 /R3 and ranged from ≈ 0. 002 to
Fig. 17b also plots Cai and Wallis's (1994) results (Eqs. (4.6.7) and ≈ 0. 3.
(4.6.8) and shows that their data agrees with ours in this case. The main result was that the value of C was higher than for the
The following empirical formula was developed to fit our CFD spherical bubble, increasing more quickly with bubble concentra-
data (see Fig. 17b): tion c (Fig. 21b). We were unable to find an analytical solution for
this case. The following empirical formula based on the cell model
C = C0 (1 − 1. 048c − 0. 588c2 ). (4.6.9) (Eq. (4.6.2)) was developed to fit our CFD data:

Laurien and Niemann (2004) performed the same calculations, but C = Ce (1 + 2c)/(1 − c) (4.6.10)
obtained an increase (Eq. (4.6.5)). The reason for this discrepancy
may be their use of a tall anisometric domain where H = 5D (see which can be used for moderate values of c. When E=0. 5 and c  0. 2,
their Figs. 8 and 9). errors of approximately only ≈ 1% are obtained. For flatter bubbles,
the error is larger.
4.6.2. Bubble in a spherical domain
4.6.2.1. Spherical bubble. Definition sketch is shown in Fig. 18 and 5. Conclusions
the results are given in Tables 13 and 14 and Figs. 19 and 20. 2D ax-
isymmetric coordinates were used, with the parameters correspond- The added mass coefficient C of dispersed particles was calcu-
ing to case B in Table 1 (80 CPB are used for the main results; CPB lated using computational fluid dynamics (CFD). As our interest is
20–160 are used for the grid independence study). primarily in gas–liquid systems, the calculations were performed for
The free-slip boundary condition was used at the wall. The bubble various configurations of gas bubbles in liquids. The value of C was
size d was fixed at 2 mm. By modifying D the shell size D/d ranged obtained from initial bubble acceleration for very short time inter-
from 1.1 to 8. The grid independence study showed that where vals of ∼ 10−5 s, where motion is controlled by buoyancy and in-
CPB = 80 the results can be considered grid independent; the differ- ertia. In such cases, the effects of other phenomena are negligible
ence between 80 and 160 CPB being less than 1% except for D/d = 2. 5 (e.g., velocity dependent forces, bubble deformation, boundary layer
(about 5%) (Table 13). The influence of liquid viscosity on the simu- development, vorticity production, etc). Under these conditions, the
lations was negligible (Table 14). Due to the formulation of the vol- difference between slip and no-slip boundary conditions is small,
umetric surface force in VOF (Fig. 19), a realistic surface tension of and the flow is effectively inviscid and potential. Consequently, our
0.073 N/m caused problems for the small domains (D < 5 mm). It was, results obtained for gas bubbles in liquids also apply to drops and
therefore, reduced by a factor of 10–100 to secure grid independence solids in fluids.
without affecting the output. This situation corresponds to that of an Several flow situations were considered with one or more bubbles
infinite 3D random bubble array where c = (d/D)3 and ranged from in 1D, 2D and 3D configurations. In each case, the effects of both
0.002 to 0.75. The main result was that C increases with bubble con- numerical and physico-chemical parameters were tested. Although
4594 M. Simcik et al. / Chemical Engineering Science 63 (2008) 4580 -- 4595

the value of C should not be dependent on the material properties Others

of the phases (e.g., viscosity and surface tension), these quantities
AM added mass
can impact the numerical procedures included in the volume of fluid
CPB number of grid cells per bubble diameter
(VOF) model. For example, in some cases it is necessary to reduce
Eo Eotwos number
the surface tension value in order to obtain grid independence and
Mo Morton number
reasonable output. Also, in certain situations, our CFD results were
Re Reynolds number
at variance with solutions obtained by other authors. We intend to
We Weber number
investigate the reasons for this in our future work.
G gas
Our results show that CFD is a useful tool for providing informa-
L liquid
tion about added mass in dispersions. With relative ease, it enables
one to compute C for situations inaccessible to analytical methods
and experimental measurements. Where we were unable to find re- Acknowledgements
sults for certain cases in the literature, we developed simple empiri-
cal formulas and fitted them to our CFD data. These simple formulas We gratefully acknowledge the financial support of GACR (Grant
are suitable for further use. Particularly interesting were the results nos.: 104/06/1418, 104/07/0111, 104/08/H055), GAAV (Grant nos.:
that showed that the presence of a third bubble can change the pat- IAAX 00130702, IAA200720801), MSMT CR (Project no. KONTAKT -
tern of pair-wise bubble interactions that form the starting point for ME 952) and AVCR CZ (Project no. 11 - 20213). We also wish to thank
averaging efforts aimed at acquiring mean-field macroscopic equa- our colleagues J. Havlica and J. Vejrazka for reading the manuscript
tions for multiphase flows. and providing valuable feedback.

References
Notation
Batchelor, G.K., 1967. An Introduction to Fluid Dynamics. Cambridge University Press,
Cambridge.
a body acceleration (vacuum a0 ,fluid a1 ),m s−2 ;ellipsoid Becker, S., et al., 1994. Gas–liquid flow in bubble columns and loop reactors: Part II.
semiaxis,m Comparison of detailed experiments and flow simulations. Chemical Engineering
b ellipsoid semiaxis,m Science 49 (24B), 5747–5762.
Bhole, M.R., Joshi, J.B., 2005. Stability analysis of bubble columns: predictions for
c bubble volume fraction in the domain,dimensionless regime transition. Chemical Engineering Science 60, 4493–4507.
C added mass coefficient,dimensionless Biesheuvel, A., Spoelstra, S., 1989. The added mass coefficient of a dispersion of
d bubble diameter,m spherical bubbles in liquid. International Journal of Multiphase Flow 15 (6),
911–924.
D distance,domain diameter or width,m
Birkhoff, G., 1960. Hydrodynamics: a Study in Logic, Fact and Similitude. Princeton
E ellipsoid bubble aspect ratio a/b,dimensionless University Press, Princeton, NJ.
F force on body in fluid,N Brackbill, J.U., et al., 1992. A continuum method for modeling surface-tension. Journal
g gravity (g reduced),m s−2 of Computational Physics 100 (2), 335–354.
Brennen, C.E., 1982. A review of added mass and fluid inertial forces. Report no. CR
h distance,spacing,m 82.010, Naval Civil Engineering Laboratory, Port Hueneme, California.
H distance,spacing,domain height,m Cai, X.L., Wallis, G.B., 1992. Potential flow around a row of spheres in a circular
L length scale,m tube. Physics of Fluids a-Fluid Dynamics 4 (5), 904–912.
Cai, X.L., Wallis, G.B., 1993. The added mass coefficient for rows and arrays of
m body (bubble) mass,kg spheres oscillating along the axes of tubes. Physics of Fluids a-Fluid Dynamics
p pressure,Pa 5 (7), 1614–1629.
r bubble radius,m Cai, X.L., Wallis, G.B., 1994. A more general cell model for added mass in two-phase
flow. Chemical Engineering Science 49 (10), 1631–1638.
R domain radius,m Chahed, J.V., et al., 2003. Eulerian–Eulerian two-fluid model for turbulent gas–liquid
t time,s bubbly flows. International Journal of Multiphase Flow 29 (1), 23–49.
T deformation time scale,s Clift, R., et al., 1978. Bubbles, Drops, and Particles. Dover Publications, New York.
Deen, N.G., et al., 2001. Large eddy simulation of the gas–liquid flow in a
u fluid velocity,m s−1 square cross-sectioned bubble column. Chemical Engineering Science 59 (8–9),
v body (bubble,particle) velocity,m s−1 6341–6349.
V body (bubble) volume,m3 ;velocity scale,m s−1 Delnoij, E., et al., 1997. Dynamic simulation of dispersed gas–liquid two-phase flow
using a discrete bubble model. Chemical Engineering Science 52 (9), 1429–1458.
z characteristic bubble dimension,m Fluent_Inc. 2003. Fluent 6.1 User's Guide.
Helfinstine, R.A., 1974. Unsteady potential flow past a group of spheres. Computer
Greek letters and Fluids 2, 99–112.
Jakobsen, H.A., 2005. Modeling of bubble column reactors: progress and limitations.
 angle, ◦ Industrial and Engineering Chemistry Research 44, 5107–5151.
 grid cell size,m Joshi, J.B., 2001. Computational flow modelling and design of bubble column reactors.
q phase q volume fraction,dimensionless Chemical Engineering Science 56, 5893–5933.
Kendoush, A.A., 2003. The virtual mass of a spherical-cap bubble. Physics of Fluids
 dynamic viscosity,Pa s 15 (9), 2782–2785.
 fluid density,kg m−3 Kendoush, A.A., 2007. The virtual mass of an oblate ellipsoidal bubble. Physics Letters
p body (bubble,particle) density,kg m−3 A 366, 253–255.
Kendoush, A.A., et al., 2007. Experimental evaluation of the virtual mass of two
 surface tension, N m−1 solid spheres accelerating in fluids. Experimental Thermal and Fluid Science 31,
813–823.
Subscripts Keulegan, H., Carpenter, H., 1958. Forces on cylinders and plates in an oscillating
fluid. Journal of Research of the National Bureau of Standards, 423–440.
a added mass Kok, J.B.W., 1988. Kinetic-energy and added mass of hydrodynamically interacting
b buoyancy, bubble gas-bubbles in liquid. Physica A 148 (1–2), 240–252.
Kok, J.B.W., 1993. Dynamics of a pair of gas-bubbles moving through liquid. Part 1.
d drag Theory. European Journal of Mechanics B-Fluids 12 (4), 515–540.
e ellipsoid Lamb, H., 1932. Hydrodynamics. Cambridge University Press, New York.
h/v horizontal/vertical Laurien, E., Niemann, J., 2004. Determination of the Virtual Mass Coefficient for
Dense Bubbly Flows by Direct Numerical Simulation. In: The fifth International
p particle, body
Conference on Multiphase Flow, ICMF'04, May 30–June 4, Yokohama, Japan.
pe/pa perpendicular/parallel Legendre, D., et al., 2003. Hydrodynamic interactions between two spherical bubble
0 vacuum; reference value rising side by side in a viscous liquid. Journal of Fluid Mechanics 497, 133–166.
 M. Simcik et al. / Chemical Engineering Science 63 (2008) 4580 -- 4595 4595

Legendre, D., Daniel, C., Guiraud, P., 2005. Experimental study of a drop bouncing Ranade, V.V., 2002. Computational Flow Modeling for Chemical Reactor Engineering.
on a wall in a liquid. Physics of Fluids 17, 097105. Academy Press, New York.
Leon-Becerril, E., Line, A., 2001. Stability analysis of a bubble column. Chemical Robertson, J.M., 1965. Hydrodynamics in Theory and Application. Prentice-Hall,
Engineering Science 56, 6135–6141. Englewood Cliffs, NJ.
Leon-Becerril, E., et al., 2002. Effect of bubble deformation on stability and mixing Sankaranarayanan, K., et al., 2002. Analysis of drag and virtual mass forces in bubbly
in bubble columns. Chemical Engineering Science 57 (16), 3283–3297. suspensions using an implicit formulation of the lattice Boltzmann method.
Michaelides, E.E., 2006. Particles, Bubbles and Drops: Their Motion, Heat and Mass Journal of Fluid Mechanics 452, 61–96.
Transfer. World Scientific, Singapore. Spelt, P.D.M., Sangani, A.S., 1998. Properties and Averaged Equations for Flows of
Milne-Thomson, L.M., 1968. Theoretical Hydrodynamics. Dover Publications, New Bubbly Liquids. Applied Scientific Research 58, 337–386.
York. Tabib, M.V., 2007. CFD simulation of bubble column—an analysis of interphase forces
Monahan, S.M., et al., 2005. CFD predictions for flow-regime transitions in bubble and turbulent models. Chemical Engineering Journal.
columns. A.I.Ch.E. Journal 51 (7), 1897–1923. Takahashi, K., Endoh, K., 1992. A virtual mass and drag coefficients for an oscillating
Mudde, R.F., Simonin, O., 1999. Two- and three-dimensional simulations of a particle. Journal of Chemical Engineering of Japan 25 (6), 683–685.
bubble plume using a two-fluid model. Chemical Engineering Science 54 (21), Wallis, G.B., 1990. Inertial coupling in two-phase flow: macroscopic properties of
5061–5069. suspensions in an inviscid fluid. Advances in Multiphase Science and Technology
Newman, J.N., 1977. Marine Hydrodynamics. The MIT Press, London. 5, 239–361.
Niemann, J., Laurien, E., 2001. Computing virtual mass by direct numerical simulation. Wijngaarden, L.V., 1976. Hydrodynamic interaction between gas-bubbles in liquid.
Zeitschrift Fur Angewandte Mathematik Und Mechanik 81, S555–S556. Journal of Fluid Mechanics 77 (9), 27–44.
Odar, F., Hamilton, W.S., 1964. Forces on a sphere accelerating in a viscous fluid. Yih, C.S. 1988. Fluid Mechanics. West River Press, CT, USA.
Journal of Fluid Mechanics 18 (2), 302–314. Zhang, D., et al., 2006. Numerical simulation of the dynamic flow behavior in a
Oey, R.S., et al., 2003. Sensitivity study on interfacial closure laws in two-fluid bubbly bubble column: a study of closures for turbulence and interface forces. Chemical
flow simulations. A.I.Ch.E. Journal 49 (7), 1621–1636. Engineering Science 61 (23), 7593–7608.
Prosperetti, A., Tryggvason, G., 2007. Computational Methods for Multiphase Flow. Zuber, N., 1964. On the dispersed 2-phase flow in the laminar flow regime. Chemical
Cambridge University Press, New York. Engineering Science 19 (11), 897–917.