7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Bubble induced agitation at a particle Reynolds number around 10
A. Cartellier*, M. Andreotti and Ph. Sechet*
LEGI, Laboratoire des Ecoulements Geophysiques et Industriels,
CNRS/UJF/GrenobleINP, UMR 5519, BP 53, 38041, Grenoble, France
alain.cartellier@hmg.inpg.fr ; philippe.sechet@hmg.inpg.fr
Sygma Motors, SAo Jose dos Campos, Brazil
marceloandreotti @sigmamotors.com.br
Keywords: dispersed flow, induced agitation, pseudoturbulence, pair density, screening mechanisms, hybrid model
Abstract
We investigated the induced agitation in monodispersed bubbly flows that are spatially homogeneous at large scale.
Experiments performed at a fixed moderate (\10) particle Reynolds number indicate that the velocity variance exhibits a
change in the scaling with the concentration. This new feature allows to explain some apparent inconsistencies between
published results. We have also shown that the transition between the two regimes is linked with the role of the particle
diffusivity in the screening mechanisms [Phys. Rev. E 80, 065301(R) (2009)]. We present here an attempt to reproduce these
observations using a simplified version of the hybrid model and accounting only for the average twoinclusions interaction. A
numerical procedure has been adapted to solve the equations coupling the microstructure and the flow disturbance. It is shown
that a buoyancy screening mechanism is recovered in the simulation in which the pair density distribution is governed by the
relative motion of the trailing particle in the wake of the test inclusion. Various scalings observed in the experiments, such as
those governing the screening distance, are recovered with this approach. Yet, significant quantitative disagreements are found,
in particular on the magnitude of the induced agitation, and ways to circumvent such limitations are discussed.
Introduction
Inclusions (bubbles, droplets, solid particles) left to rise or
fall in an inert continuous phase induce velocity fluctuations.
The later are helpful to enhance mixing. They also give rise
to the socalled pseudoturbulence or agitation stress tensor
A=, where v' is the continuous phase velocity
fluctuation relative to its unconditional average v, X, is the
continuous phase indicator function and <> denotes the
ensemble average over all the possible positions and
velocities of the inclusions. That tensor acts as a Reynolds
stress in the continuous phase momentum balance, and it
tends to uniformize the spatial distribution of inclusions (see
for ex. [1,2]). Many investigations have been focused on the
investigation of the mechanisms controlling such velocity
fluctuations. It is clear that these mechanisms evolve with
the particle Reynolds number (see among others Caflisch &
Luke 1985, Koch & Shaqfeh 1991, Ladd 1996, Brenner
1999 for suspensions, Koch 1993, Kumaran & Koch 1993,
Cartellier & Riviere 2001 for small to moderate Rep,
Biesheuvel & Van Wijngaarden 1984, Lance & Bataille
1991, Hunt & Eames 2002, Risso et al. 2008 for large Rep].
Yet, they are still not fully understood, nor are the
corresponding scaling laws. This is especially true at
moderate particle Reynolds numbers (Rep is defined here on
the inclusion radius a, its relative velocity Ur and on carrier
fluid kinematic viscosity vC). Indeed, as shown in Fig.l,
available experiments seem to exhibit different evolutions of
the velocity variance with the dispersed phase concentration
a for similar Rep values. This is especially true for Rep
around 10 where all the experimental data shown
correspond to spherical and nearly monodispersed bubbles.
Also, in the range of small to moderate Rep, direct
numerical simulations fail so far to reproduce the observed
velocity variance evolution with the void fraction (Fig. 1).
1000 '
(v' UZ/u /
C
0 Carteller & Rviere (2001)
S Marnez Merado et al (2007)
Lance & Bataille (1991)
f S CQameretal (2001)
S Lare de Toumemlne (2001)
Chment&Maxey(2003)
E Bunner & Tyggvason (2002)
001 01
S
10 100 1000
Figure 1: Experimental data on the axial velocity variance
in bubbly flows scaled by the square of the relative velocity
Ur and by the void fraction as a function of Rep. Some direct
simulation results are also shown (Climent & Maxey data
correspond to the axial dispersed phase velocity variance).
Another way to handle these questions is to exploit
statistical twofluid models [4, 20, 21]. In such models, A
can be expressed as the sum over all test particle positions
x each position x being weighted by the corresponding
number density 4(x) of the average velocity perturbation
v*(xlxo) induced at a location x in the carrier phase by a test
particle centered at x.
A(x) = J(x) v*v*(xxo)dxo+O(a2) (1)
Here, the average disturbance v*(xlxo)=v(xlxo)v(x) is the
difference between the continuous phase velocity v(xlx) at
x conditioned by the presence of an inclusion at xo, and the
unconditional continuous phase velocity v(x). Higher order
terms in eq.(1) are in proportion of the square of the
dispersed phase concentration a. That equation has been
exploited a number of times assuming an uniform
probability for test particle positions, and that the average
disturbance equals the perturbation by an isolated particle.
The ensuing shortcomings of that approach applied to finite
Rep situations are well known [3, 7]. To circumvent these
difficulties, one or both of the above assumptions have to be
relaxed. In particular, Koch & Shaqfeh [4] have shown how
a non uniform probability for particle position can arise due
to multiple body interactions, and how that feature leads to a
screening of velocity disturbances and therefore to a finite
velocity variance. In addition, experimental evidence
demonstrates that for moderate Rep both assumptions are
incorrect (Cartellier & Riviere 2001).
To progress in our understanding of the mechanisms,
experiments have been undertaken at a fixed particle
Reynolds number close to 10, and over a wide range of void
fractions, from 0.2% up to 10%. The main experimental
results are summarised in the next section. In particular they
demonstrate that the mechanisms are also evolving with the
concentration. In parallel, we attempted to predict the
induced agitation by exploiting the eq.(1) in the framework
of the hybrid twofluid model (Achard & Cartellier, 2000,
2001). This approach as well as the simplifications we
introduced will be presented. Then, the preliminary results
will be discussed and compared with the experiments.
Summary of experimental results
Experiments have been performed in a gaslift setup. Data
were collected in the middle of a tube where both the mean
phasic velocities and the concentration were uniform. Care
was taken to ensure a narrow size distribution and clean
bubble interfaces, so that the nondimensional parameters in
that problem reduce to Rep which has been maintained
almost constant, and to the local phase fraction a. The main
results are summarized here below (more information is
available in [22]). The probability density functions (pdf) of
the axial liquid velocity fluctuations in the continuous phase
collapse once normalized by the velocity root mean square.
These single parameter distributions are neatly asymmetric,
positive realisations (i.e. velocities directed upwards)
extending up to 10 times the RMS.
The velocity variance itself, scaled by the mean relative
velocity, was found to increase with the concentration. Two
regimes were identified: a close to linear growth was
observed for a concentration above a critical void fraction
(those value is about 1.2%) while below that threshold, the
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
increase was nonlinear and much less steep. The existence
of these two different regimes is an important new result. In
particular, they allow to reinterpret available experiments at
moderate Re, and to solve the apparent contradictions
between previous results. Another aspect to be noticed in
Figure 3 is that the increase seems slightly steeper than
linear for concentrations close to 10%. This possible trend
needs to be checked further by experiments at larger void
fractions.
0.01
0.001
0.0001
a=0.16%
. a=0.59%
x a=1.2%
 =4.5%
A a=8.1%
A
x x' (w 1 10
5 W'(W)1/2 10
w (
Figure 2: Pdf of the axial velocity fluctuations in the
carrier phase normalised by the root mean square
(8
In all the conditions investigated, experiments indicate that
the velocity disturbances are screened at a finite distance
behind the test bubble. Also, a deficit of neighbour is always
registered in the rear of test bubbles. The axial extents of the
velocity deficit Xv and of the screening zone Xs are
nearly the same: they decrease with the concentration.
These quantities are always much smaller than the tube
diameter, so that the magnitude of the induced agitation is
indeed insensitive to the vessel size. The induced agitation
is therefore liable to a local closure.
(v, / U )2
1.I
S
0.7
001
0001 001
Figure 3: Evolution of the axial liquid velocity variance
scaled by the mean relative velocity with the void fraction
(8
In [22], we have also shown that the linear regime
corresponds to a weak diffusivity of the leading particle,
that is to say its trajectory is almost unaffected by the
presence of the trailing inclusion. On the opposite, the
nonlinear regime corresponds to a higher diffusivity. In that
case, it is the lateral motion of the test inclusion instead of
a (abs)
the carrier phase viscosity alone that controls the wake
spreading. Such a mechanism, that has been proposed by
Brenner (1999) for suspensions, happens to be also effective
at moderate Rep. For the socalled linear regime, we have
proposed a simplified analysis that predicts a screening
length X/a evolving as Rep/a (see [22]). The 1/a behaviour
was indeed recovered in the experiments for void fractions
above the critical value, but the linear increase of X with
Rep was not tested for lack of data.
Let us now introduce the simplified version of the hybrid
model used to compute the velocity variance.
Hybrid modelling and the induced agitation
In order to compute the stress tensor in the liquid phase A,
we exploited an hybrid Euler/Euler twofluid formulation
developed for dispersed twophase flows (Achard &
Cartellier, 2000, 2001). That approach relies on a statistical
description of the dispersed phase where the configurations
correspond to the state of the particles, coupled with a
continuous description of the carrier flow. The ensuing
equations are structured as a hierarchy of coupled phasic
balance equations. The first level of the hierarchy governs
unconditional variables, such as the average inclusion
velocity, the average phasic velocity for the carrier phase...
Next levels bring in more and more refined flow
descriptions using variables defined over subsets of the
outcomes where one, two, etc... test inclusions are fixed in
space. In such nonlinear multibody problems, higher order
dynamics must be uncoupled at some stage in order to
obtain a closed set of equations. The truncation is mainly
based on diluteness.
We are interested here in monodispersed inclusions at small
to moderate particle Reynolds numbers, say between 1 and
100. Our bet is that, in that Rep range, the induced agitation
should be rather correctly captured by accounting for two
bodies interactions.
Unconditional flow:
In the experimental situation considered, the mean dispersed
flow is one dimensional, steady, almostfully developed and
remains uniform at large scales (i.e. at scales much larger
than a). In the model, we thus considered an infinite,
uniform base flow. Let us recall that our experiments do
show that the induced agitation is independent on the vessel
size provided that later is large enough. In these conditions,
the unconditional momentum equation for the dispersed
phase takes a simple form, which in the limit of a vanishing
inclusion density writes:
_dx ax ax
(+R1) I] FV= (2)
where 6 is the volume of an inclusion, and Kp is the
dynamic viscosity of the carrier phase. Due to the scale
separation between the inclusion radius a and the vessel
dimension L, the contribution from the operator R1, which is
of order O(a/L)2 (see [21]), vanishes. Eq.(2) simply states
that the resisting force F* acting on a test inclusion
equilibrates the generalized Archimedian force induced by
the unconditional velocity (v) and pressure (p) fields in the
continuous phase. In the problem considered, the inclusion
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
size and the mean relative velocity Ur are given quantities,
so that F* is also given and fixed. The question now is to
evaluate the average perturbation v* in the continuous
phase.
Equation for the average disturbance in the carrier phase:
The averaged first order disturbance in the continuous phase
is the one evaluated for a fixed inclusion at xo while all
others inclusions scan the available phase space. These
velocity v* and pressure p* disturbances at first order obey
continuity and momentum equations that can be simplified
in the dilute limit. In particular, in that limit, v* remains
solenoidal. Since the dispersed flow we consider is uniform
at large scale, unconditional quantities relative to velocity,
pressure and concentration are slowly varying over
distances of the order of a. This not so for the disturbances,
and thus there is no longer any scale separation at this order.
As a first approach, we also neglected the operator R1 in the
momentum equation for the disturbance. Its simplified form
at x given a test inclusion centered at xo writes:
P V*+(v*Ur ) + Lx PCAv*
lat dx dx (3)
= [1 S(x I xo)]OF*(x) S(x I x)0F**(x I x)
In the l.h.s., one recognizes the NavierStokes equation for
the flow around an isolated inclusion traveling at a relative
velocity Ur. Yet eq.(3) differs from the one associated with
an isolated particle dynamics because of the presence of
source terms in the r.h.s. that do not vanish for finite Rep.
The term (F* arises from the interfacial momentum
exchange between unconditional fields (note that a= ( 6 ).
The remaining terms arise from the interfacial momentum
exchange that enters the momentum balance for the
continuous phase velocity field at a location x given an
inclusion at xo. The equation (3) is supplemented with the
appropriate boundary conditions at the surface of the
inclusion (depending on its nature, solid, clean bubble...),
and with the condition that the disturbance decays to zero
far from the test particle, namely: v*(xlxo)0 as xxol .
Clearly, eq.(3) cannot be solved without the knowledge of
the quantities S and F**. Here, S is the pair density
".''\,x) normalized by the number densities ( at x and at
x. F**(xlxO) represents the extra force acting on the particle
at x due to the presence of a test particle at xo. In other
words, the total force on the inclusion at x given another one
at x is the sum F*(x) + F**(xlx). Even if the pair
probability for inclusion positions is uniform, that is S=I
everywhere, a source terms remains in the r.h.s. of eq.(3)
since F** does not vanish whenever some amount of inertia
is present (Rep,0). Let us now introduce the equations
driving the quantities S and F**.
Pair density equation:
The equation for the pair density comes out from the
continuity equation at the second order for the dispersed
phase. For the homogeneous situation at large scale
considered here, and distinguishing between quick (i.e.
varying over a scale a) and slow (i.e. varying over a scale L)
variables, that equation simplifies into:
*{(Au +Au*) S(x,r)}= 0 (4)
r
Here, the trailing particle center is at x and the test inclusion
center is located at x+r. The quantity Au=u(x+rlx)u(xlx+r)
represents the difference in the two inclusions mean
velocities. Since all inclusions possess the same mean
relative velocity and since the continuous flow field is
uniform, Au=0 here. The quantity Au*=u*(x+rlx)u*(xlx+r)
is similarly defined, but it is based on the disturbances of the
inclusions velocities. As for v*, the later are defined as
u*(xlxo)=u(xlxo)u(x). In these first computations, we
assume that the leading particle motion is unaffected by the
trailing one, so that u*(x+rlx) z 0. The remaining term
u*(xlx+r) that represents the influence of the wake of the
leading inclusion on the motion of the trailing one is kept.
Its determination is based on its definition, namely one can
write: u*(xlx+r) = u,*(xx+r) + v*(xlx+r). The disturbance
v* is determined from the momentum equation (3): this is a
first coupling between eq.(3) and (4). The determination of
the quantity Ur* that represents the relative velocity of the
trailing inclusion with respect to the disturbed field v*, will
be detailed here after.
A first boundary condition associated with eq.(2) is that the
positions of widely separated particle are uncorrelated, i.e.
S1 as xx>oo. The second condition is an integral one:
it states that, over the domain, deficit and excess of S should
compensate each over.
Expression of the extra force on the trailing inclusion:
The extra force F** can be evaluated from the next level of
the hierarchy for the continuous phase, namely the two
bodies problem. Yet, that step can be bypassed by
exploiting the momentum balance for the dispersed phase at
first order. Indeed for light particles (vanishing density
limit), the later reduces to:
p[ ~ ax F**
(l+R) PcAv* v* = (5)
a x ax ax (1
That equation states that the force disturbance F** is
equilibrated by the generalized Archimedian force induced
by the disturbed field v*, p* in the continuous phase. Here,
the scale separation is no longer effective, and the operator
R1 should be kept. Yet, as a first approach, we choose to
neglect its influence, that is to say that the finite size of the
trailing particle will not be accounted for.
Relative disturbed motion of the trailing inclusion:
The last step required to close the above set of equations is
an expression for the relative velocity Ur*. We assume
(subject to a posteriori verification) that the relative velocity
ur* between the trailing inclusion and the disturbed flow
field v* remains small (in particular ur* << Ur). This is
equivalent to postulate a vanishing Reynolds number for the
relative motion between the two inclusions. In such
conditions, the resisting force F** is linearly related with
the relative velocity ur*, namely:
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
using either a Stokes (for solid particles, k=6) or an
HadamardRybczynski (for clean bubbles, k=4) formula.
This leads to a second source of coupling between eq.(3)
and (4). Now, the equations (3) and (4) supplemented by the
expressions for F** and for Ur*, and by the boundary
conditions do form a closed set. Finaly, let us note that in
this simplified model proposed here, there is no equation
accounting for the perturbation of the motion of the test
particle. Therefore, that model will not be able to reproduce
the observations at the smallest void fractions for which the
particle diffusivity is an important ingredient of the pair
interaction.
Numerical Scheme
The coupled equations (3) and (4) governing the average
disturbance v*, p* and the normalized pair density S
respectively were solved by way of an iteration procedure.
Starting with S=0 everywhere, a first disturbance field is
evaluated, which is then used to predict a first solution for
the pair density. The later is reinjected in the momentum
balance (3) to update the disturbance fields, and so on until
convergence. Such a convergence was obtained without
difficulties except for very dilute conditions. The origin of
such a defect is not fully understood: it could be due to
some inaccurate numerical treatment because the source
terms are indeed very weak in that case. It may also be
related with some consistency of the set of equations, since
in that approach we do not account for the test particle
diffusivity while experiments do show that the later
becomes central at low void fractions.
The JADIM code from IMFT, based on a finite volume
method with secondorder accuracy, has been adapted to
solve these coupled equations. The twodimensional polar
grid used was composed of 72 by 72 cells. Note that the
numerical solution was first tested for the far wake of a
single inclusion at Rep=l and 35, and a good prediction of
the velocity in the wake was ensured up to a distance about
50 to 60 radii downstream.
Results and Discussion
Simulations were mainly performed for light solid
inclusions. Let us first qualitatively discuss the disturbance
force F**. As shown in Figures 4 and 5, where the test
particle is centred at (0,0), the modulus of F** is about two
orders of magnitude below the average resisting force F*.
This result implies that the disturbed relative velocity ur* is
much smaller than Ur, and thus that the linear drag
relationship introduced in eq.(6) is correct.
Concerning now the orientation of F**, let us first note that
the force driving the motion of the trailing bubble in the
disturbed flow field is F**. Along a vertical (i.e. the
direction x, the gravity is oriented from left to right), that
force tends to slow down the motion of the approaching
bubble on the axis (red zone near the axis in figure 4), and
to accelerate it on the side (blue region in figure 4).
F** = k7rcau,
oF
500x10'
695x10"
297x10 "
238x1011
5 84x10n
SOaxio D
2 97x105
3 87x10 u
5 13x10'
1 21x101
5 50x 10
326x1 0"
1 01x10i2
500x1012
negative
1 1 I I I
10 5 0 5 10 15 20 25 30
x/a
Figure 4: Axial component of F**/F* in the vicinity of
the test inclusion, and for a=0.042 and Rep=10.
30 r
9 13x10
3 98x100"
1 38x100
9.83x1014
5.26x104
2.38x1014
1.71X1014
1.19x10"
9.10X10i0
4.41x10'
8.10x10 6
1.36x10
2.56x10"
1.26x10
4 53x10
10 5 0 5 10 15 20 25 30
x/a
Figure 5: Radial component of F**/IF* in the vicinity of
the test inclusion, and for a=0.042 and Rep=10.
The main effect of F** is along the radial direction since it
is the only active force in that direction (F* is aligned with
x). For a trailing particle located on the axis behind the test
inclusion, that force deviates its trajectory away from the
axis (positive red region near the axis in figure 5). The
opposite holds when the trailing particle is off axis (blue
region outside the wake in figure 5). Therefore, the average
trajectory of a trailing particle approaching a test particle
from behind corresponds to an axial slow down combined
with an ejection away from the axis, followed then by an
axial acceleration until it catches up the test particle. Such a
binary interaction is sketched in Figure 6.

Figure 6: Average relative trajectory
of inclusions.
According to this "average" relative
M 15
of an interacting pair
trajectory of the test
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
particle, one expects a deficit of neighbours in the near
wake of any test inclusion. This is indeed what the model
predicts as shown by the distribution of the pair density in
the vicinity of the test inclusion (see blue and green zones in
Figure 7). That deficit can be quite strong: for the flow
conditions of Figure 7, S is as small as 0.1 in the near wake.
Another feature is the increased probability for horizontally
aligned pairs (see the red zone on the side of the test particle
in Figure 7). Indeed, once the trailing particle catches up
with the leading one, it enters a region where the axial force
is reversed while the radial one is directed towards to the
inclusion: such a force reversal traps the trailing inclusion
and explains the increase in the pair density above unity. In
addition, the increased pair probability (S>1) extends along
the outer wake region and forms like shell. Such features are
fully consistent with the experiments: they demonstrate that
the force F** is truly the key for the formation of the
microstructure in monodispersed systems.
10
a a
am
4
I . . I . . L . . L . . L .
0 5 ie 15 20
x/a
Figure 7: Computed spatial organisation of S around the
test particle (a=0.021, Rep=10).
Before comparing experiments and model predictions in a
quantitative terms, it is useful to consider now how the
forcing term in eq.(3), that is 4{(1S)F*SF**}, is
distributed around the test inclusion and how its affects the
velocity disturbance in the wake. Again, in the radial
direction, the F** term is the only one active: near the
symmetry axis, it tends to spread the wake so that the later
would be enlarged compared with the one formed behind an
isolated inclusion. This is agreement (at least qualitatively)
with the observations made in [9]. In the axial direction,
since the modulus of F** is quite small compared with F*I,
that term can contribute to the forcing only when S is near
unity. Comparing figure 7 with figures 4 and 5, it is clear the
axial contribution of F** to the forcing is marginal. The
main axial forcing comes therefore from (1S)F*, and it is
essentially active in the pair density deficit region where S
strongly differs from unity. That force density is directed
downward (toward positive x in the figures) that is it is
opposite to v*. Therefore, that forcing hampers the
momentum flux in the wake, and contributes to its quick
decay. Such an extinction of the velocity disturbance is
illustrated Figure 8 for various concentrations. Clearly, for a
fixed Rep (here a fixed F*), the screening effect is more
pronounced as the concentration increases since the forcing
term is proportional to the number density 0.
0001
x/a
Figure 8: Screening of the disturbance velocity along the
axis behind the test particle for various concentrations (light
solid inclusions, Rep=10). The solid line corresponds to an
isolated inclusion at the same Rep.
The coupled mechanism leading to disturbance screening in
monodispersed systems at moderate Rep is therefore the
following. The force modification F** due to the wake of a
nearby test inclusion induces, for a trailing inclusion, a
deviation from the rectilinear trajectory and thus it produces
a pair density deficit in the wake. In its turn, that deficit
induces a forcing of the disturbance flow that is in
proportion of the mean drag force F*. This forcing
extinguishes the velocity perturbation at a finite distance.
Hence, no momentum escapes from that celllike structure
surrounding the test inclusion, and the integral involved in
the expression of the induced agitation (eq.1) will be then
finite, as we shall see below.
It is here interesting to notice that the forcing, which is in
proportion of the mean drag force F* multiplied by the
difference in concentration between the region of deficit and
the outer uniform zone of the flow, is exactly the
buoyancy screening effect identified by Koch (1993). Yet,
the ejection mechanism responsible for the deficit of
neighbours in the near wake drastically differs from the one
exhibited by Koch for polydispersed systems. In the later
case, there exists a mean relative velocity between two
inclusions (i.e. the term Au in eq.(4) differs from zero)
which is controlled by the difference in size or density
between the two particles, and the interaction of that relative
velocity with the wake velocity field gives rise to a lift force
that ejects the trailing particle outside the wake. In
monodispersed conditions, Au is zero as we have already
seen, and the relative velocity ur* of the trailing bubble with
respect to the wake velocity field is only due to the
generalised Archimedian force at this order (eq.5). Note that
since ur* is much smaller than Ur, lift effects which can
possibly be included in the analysis, are negligible here.
Such differences between monoand polydispersed systems
become manifest when one compares the screening lengths.
In polydispersed systems, the screening length is predicted
by Koch to evolve as X/a oc P/(aRep). Naturally, its
expression involves the polydispersity P of the dispersed
phase. In monodispersed conditions, the equilibrium
between the momentum flux due to the test inclusion and
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
the opposite action of the source ((1S)F* we discussed
above has been exploited in [22] to derive an expression for
the screening length. The later is found to evolve as X/a oc
Rep/a. The proportionality of X with Ur found here arises
from the fact that it is the mean buoyancy force (or drag F*)
that drives the disturbance field v*, p* and thus, by way of
the equilibrium (5), the force F** acting on the trailing
particle and that generates the microstructure. Another
consequence of that scaling is a longer interaction time and
thus a longer extent of the pair density deficit
monodispersed systems compared with polydispersed ones,
the ratio of which being of order Rep2 at a fixed
concentration.
Evolution of the screening distance with a and Rep
We exploited the simulations to examine how the screening
distance evolves with the two control parameters, namely a
and Rep. For a fixed Rep, the screening distances defined
either on S (X,) or on v* (Xv) are plotted versus the
concentration in figure 9. Most of the data follow the 1/a
behaviour expected from the above mentioned analysis and
found in experiments (Note that, as discussed in [22], the
experimental data for the two smallest void fractions in
Figure 9 pertain to the range where the particle diffusivity is
important and should not be considered in the discussion).
Yet, two simulated data for X, do not follow that trend most
probably because the correct solutions should have
exceeded the extent of the computational domain (The
effect of the limited domain used in these preliminary
computations is also visible in Figure 10). An odd feature in
Figure 9 is the quantitative difference between X, and Xv;
such a difference was not observed in the experiments.
Another noticeable difference is the much larger extend of
the simulated screening lengths compared with experiments.
This is probably due to the fact that the finite size of the
trailing inclusion was neglected in that simplified version of
the model. Indeed, with regards to the actual extents of X,
which are typically less than 10 radii, such an assumption is
only very marginally acceptable.
100
 Extent of the comptation domain^
U
10
X X/ a simulated
SX / a simulated
X/a experiment
0001 001
01 a(abs)
Figure 9: Evolution of the screening length with the
concentration (simulations are for light solid inclusions at
Rep= 10).
We also investigated the evolution of the predicted
screening distance with Rep. As shown in figure 10, the
screening both in terms of extend and of magnitude is
enhanced as the particle Reynolds number increases. This is
in qualitative agreement with the available observations [9,
h i
* *`*
t
t
1
i
I
U
22]. For the smallest concentration considered (a=0.021),
the increase of X, is almost linear in Rep, as expected from
the proposed scaling. Yet as the concentration increases, the
screening length still increases with Rep but in a less steep
manner (typically X,/a evolves as Rep044 for a =0.126).
Such deviations from the linearity may possibly correspond
to the limit of the assumed diluteness used to simplify the
model. We also considered the lateral extent of the deficit
region although we have no experimental data to compare
with. It happens that the difference between X, and Xv is
even more marked in the radial than in the axial direction. A
close examination of the simulations indicates that the
region of pair density deficit is strongly squeezed near the
axis of symmetry, while the velocity disturbances are much
more laterally spread. Such effects are almost unaffected by
the concentration but are sensitive to the particulate
Reynolds number. Again, we suspect that the point particle
approximation used for the trailing particle could be
responsible for the enhanced confinement of the pair density
deficit near the axis. Such a confinement is compensated by
an extra elongation of the structure in the axial direction, a
feature that has already been identified in figure 9.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Reynolds numbers. As already said, the computations were
restricted at the smallest concentrations because of
convergence difficulties. The first striking feature is the
linear increase of v'2/U,2 with a for all the particulate
Reynolds numbers. The linearity holds almost over the
entire concentration range considered. Yet, clear departures
from linearity appear at the largest concentrations, close to
10%, a value for which the diluteness assumption is
probably no longer valid. For Rep=l, the experiments (see
[9]) indicate that the linearity holds from 0.1% up to 8% in
concentration: the simulation results are therefore quite
consistent with these findings. For Rep=10, we have seen in
Figure 3 that the linearity starts at a critical concentration
about 1% and holds up to about 10%, with a possible
steeper increase at the largest void fractions. These trends
are recovered in simulations except for the change in the
regime evoked in the previous section. As already
mentioned, that change is associated with the influence of
the particle diffusivity. Since the later is not accounted for in
the model, it is not surprising that the change in the scaling
is not captured in these simulations. Interestingly, by
neglecting the diffusivity, the agitation remains linear with
the concentration down to very small void fractions.
Still, despite all these qualitative agreements between
experiments and simulations, the amount of induced
agitation is significantly underestimated by the model. This
is even more critical when considering clean bubbles instead
of solid particles (see Figure 11).
v' /U
 Rep=1 (solid particles)
* Rep=5 (solid particles)
* Rep=10 (solid particles)
SRep=10 (bubbles)
Rep=15 (solid particles)
S Bubbly flow experiments at Rep=10
x/a
Figure 10: Evolution of the normalised pair density S
behind a test inclusion as a function of Rep for a fixed
concentration (a=0.042, solid particles).
Therefore, a number of encouraging features concerning the
screening of velocity disturbances and the microstructure
come out from the simplified version of the hybrid model.
Yet, the quantitative agreement with the experiments is not
satisfactory. We nevertheless attempted to estimate the
induced agitation from these preliminary simulations.
Axial velocity variance:
The axial velocity variance v' in the carrier phase was
computed using eq.(1). The integration was performed over
the whole computational domain since the velocity
disturbance is mostly extinguished at the domain boundaries.
Yet, near the axis does some amount of perturbation remains
and flows out of the domain, and this effect is most
noticeable for the largest concentrations. Despite that
limitation, let us first consider the evolution with the two
independent parameters. Figure 11 provides the axial
velocity variance in the carrier phase, scaled by as a
function of the concentration and for various particulate
a (abs)
0 001
0 001
Figure 11: Simulated axial induced agitation v'2/Ur2 versus
the concentration for different particle Reynolds numbers.
Concerning now the effect of Rep, Figure 11 exhibits a neat
sensitivity of the predicted induced agitation to that
parameter. This is expected according to the discussion on
the screening mechanisms. Let us recall that for the
experiments available at Rep close to unity (for Rep between
0.3 and 3.5, see [9]) and for those presented in figure 3 (for
8< Rep <12), the axial velocity variance scaled by Ur2
correlates with the ratio a/Rep (at least for a/Rep > 103 that
corresponds to nearly linear increase of the agitation with
the concentration). We tested whether such a correlation
also holds for the simulations, and, as shown in figure 12,
this is indeed so. Moreover, the slopes for the experimental
and for the simulated results are quite similar (discarding
again the simulated data at the largest concentrations). This
is quite encouraging, especially since an a/Rep behaviour for
the quantity v'2/Ur2 was derived by Koch for buoyancy
screening. Yet, and contrary to experimental results, the
simulations for the various Rep do not collapse on a single
curve when the agitation v'2/Ur2 is plotted against a/Rep.
This extra sensitivity to the particle Reynolds number that
remains in the present simulations most probably originates
from the deficiencies in the prediction of the microstructure
we have discussed above. In particular, the simulated
screening distances are far too large and they are associated
with pair density deficits that are concentrated in the close
vicinity of the axis of symmetry. It is interesting to compare
these results with those obtained for an imposed
microstructure. In [23, 24], we considered a thick pair
density deficit region in the rear of a test particle, with an
extent about 4 radii coming from experiments. Eq.(3) was
solved (the contribution from F** was neglected) and the
induced agitation was found to be quite close to
experimental results as shown in Figure 12. This feature
demonstrates the sensitivity of the induced agitation to the
microstructure. Besides, it supports our conclusion
concerning the origins of the limitations of the simulations
presented here and which couple the disturbance and the
pair density equations. More investigations are required to
check whether a part of the observed shortcomings arise
from some numerical problems or if they are only due to the
treatment of the trailing inclusion as a point particle. Note
that this assumption has been used at various steps in the
derivation of the simplified model. First, it is required for
the closure on F**. Second, it has also been exploited to
discard the operator R1 both in the equilibrium eq.(5) and in
the evaluation of the forcing terms in the r.h.s. of eq.(3).
10
 Rep=1 (solid particles)
V2 2 "" Rep=5 (solid particles)
v /U Rep=10 (solid particles)
r Rep=10 (bubbles)
Rep=15 (solid particles)
Experiments at Rep=10 [Ref 22]
S Experments at Rep=0(1) [Ref 9]
1 Simulations with an imposed 0
pair density [Ref 23, 24] 0
01
^ I/ '. .
o .. o .o
001 "
o/Re
o p
0001
00001 0001 001 0 1
Figure 12: Comparison between simulations and
experiments on a v'2/Ur2 versus a/Rep plot. The socalled
"linear" regime is located at the right hand side of the
double bar.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Conclusions
Experiments and simulations based on the hybrid model
have been presented for monodispersed uniform flows at
small to moderate particulate Reynolds numbers. Results
indicate that the active screening mechanism in that case is
due to particleparticle interactions mediated by the
interstitial dense carrier phase. More specifically, the motion
of a particle trailing behind a test inclusion is deflected
outside its wake by the action of the continuous phase
disturbance field. In its turn, the pair density deficit thus
created in the rear of the test inclusion drives the forcing
terms and alters that disturbance. The coupling between
these two effects leads to the extinction of the disturbance at
a finite distance, which has been shown to be of order X/a
a/Rep. Interactions between inclusions of the same size and
density have a longer duration, and lead to an enhanced pair
density deficit compared with polydispersed conditions.
The resulting induced agitation is also much (typically 10
times) stronger in mono than in polydispersed systems.
So far, the hybrid modeling adapted to handle simplified
binary interactions has been able to capture the key features
of the above mechanism. However, the predicted amount of
fluctuations is still significantly underestimated and it
exhibits a sensitivity to the particle Reynolds number that
differs from the one detected is experiments. The origin of
these defects is most probably related with the assumption
of a point trailing particle we introduced to drastically
simplify the system of equations. This point needs to be
tested. Another interesting extension of that approach would
be to take into account the disturbance induced on the
motion of the test particle. This would permit to investigate
the regime that we have exhibited at small concentrations
and for which the particle diffusivity tends to become a key
ingredient in the pair interaction.
Acknowledgements
The authors would like to thanks Dominique Legendre from
IMFT and Christophe Corre from LEGI for their guidance
and assistance in the adaptation of the JADIM code.
References
[1] Achard J.L., Cartellier A., Local characteristics of upward
laminar bubbly flows. PhysicoChemical Hydrodynamics, 6, 5/6,
841852 (1985).
[2] Ho T., S6chet P. & Cartellier A., Sensitivity analysis of the
transverse distribution of bubbles in laminar flows in the limit of
small particle Reynolds numbers. Proc. CHISA Meeting Praha,
Czech Rep. (2002).
[3] Caflisch R., Luke J., Variance in the sedimentation speed of a
suspension, Phys. Fluids 28, 759760 (1985).
[4] Koch D., Shaqfeh E., Screening in sedimenting suspensions, J.
FluidMech. 224,275 303 (1991).
[5] Ladd A., Hydrodynamic screening in sedimenting suspensions
of nonbrownian spheres, PRL 76, 13921395 (1996).
[6] Brenner M.P., Screening mechanisms in sedimentation, Phys
Fluids 11 (4) 754772 (1999).
[7] Koch D., Hydrodynamic diffusion in dilute sedimenting
suspensions at moderate Reynolds numbers. Phys. Fluids 5, 1141
(1993).
[8] Kumaran V, Koch D., The effect of hydrodynamic interactions
on the average properties of a bidisperse suspension of high
Reynolds number, low Weber number bubbles, Phys. Fluids, A, 5,
11231134 (1993).
[9] Cartellier A. & Riviere R., Bubbleinduced agitation and
microstructure in uniform bubbly flows at small to moderate
particle Reynolds numbers. Phys. Fluids 13, 2165 (2001).
[10] Biesheuvel A. & Van Wijngaarden L., Twophase flow
equations for a dilute dispersion of gas bubbles in liquid, J. Fluid
Mech. 148, 301 318 (1984).
[11] Lance M. & Bataille J., Turbulence in the liquid phase of a
uniform bubbly airwater flow. J. FluidMech., 222, 95118 (1991).
[12] Hunt J. C. R. & Eames I., The disappearance of laminar and
turbulent wakes in complex flows. J. FluidMech. 457, 111132
(2002)
[13] Risso F., Roig V, Amoura Z. & Billet AM., Wake attenuation
in large Reynolds number dispersed twophase flows. Phil. Trans.
R. Soc. A, 366, p.2177 (2008).
[14] MartinezMercado J., PalaciosMorales C., and Zenit R.,
Measurement of pseudoturbulence intensity in monodispersed
bubbly liquids for 10
[15] Gamier C., Lance M. and Marie J.L., Measurement of local
flow characteristics in buoyancydriven bubbly flow at high void
fraction, Exp. Thermal & Fluid Sc., 26 (67) 811815 (2002).
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
[16] Tournemine A., Etude exp6rimentale de l'effet du taux de vide
en 6coulements diphasiques a bulles. Ph.D. INPT, Toulouse (2001).
[17] Bunner B. & Tryggvason G., Dynamics of homogeneous
bubbly flows, part 1: Rise velocity and microstructure of the
bubbles p.17, part 2: Velocity fluctuations p.53, J. Fluid. Mech.
466 (2002).
[18] Esmaeli, A. & Tryggvason, G. 1999, Direct numerical
simulation of bubbly flows. Part2. Moderate Reynolds number
arrays. J. FluidMech., 385, 325358 (1999).
[19] Climent E. & Maxey M., Numerical simulations of random
suspensions at finite Reynolds number. IJMF29, p.579, (2003).
[20] Zhang D. Prosperetti A., Ensemble phaseaverage equations
for bubbly flows, Phys. Fluids, 6 (9), 9 29562970 (1994).
[21]Achard J.L. & Cartellier A., Laminar dispersed twophase
flows at low concentration. Arch. Mech. Part 1: Generalized
system of equations 52, p.25 (2000). Part 2: Disturbance equations
52, p.275 (2000). Part 3: Pseudoturbulence 53, p.123 (2001).
[22] Cartellier A., Andreotti M. and Sechet Ph., Induced agitation
in homogeneous bubbly flows at moderate particle Reynolds
number, Phys. Rev. E. 80, 065301(R) (2009) DOI:
10.1103/PhysRevE.80.065301
