Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Modelling and simulation of drop breakup downstream of an orifice
Riccardo Maniero*1'3, Olivier Masbernat1'3, Eric Climent2'3 and Frederic Risso2'3
1 Laboratoire de Genie Chimique Universite de Toulouse (INPTUPS) et CNRS
4 allee Emile Monso, Toulouse, 31432 cedex 4, France
2 Institute de Mecanique des Fluides Universite de Toulouse, (INPTUPS) et CNRS
Allee du Professeur Camille Soula, Toulouse, 31400, France
3 Federation de recherche FERMAT, CNRS Toulouse, France
Email addresses : Riccardo.Maniero@iunil.ch, olivier.masberat@iensiacet.fr, climent@imft.fr, risso@imft.fr
*Current address:
University de Lausanne, UNIL, Institut de geophysique UNILSorge, Lausanne, 1015, Suisse
Author to whom correspondence should be addressed: climent@imft.fr
Abstract
This work aims at modelling numerically drop breakup in the turbulent inhomogeneous flow that develops downstream of an
orifice inserted in a pipe. Detailed experiments were carried out for investigating the relation between the turbulence
distribution and the breakup events. It has been shown that turbulent fluctuations at the scale of the drop is the major
contribution to the drop deformation and eventually breakup. Results on the drop dynamics are supplemented by PIV
measurements of the carrying fluid flow.
This work reports numerical simulations of these phenomena. The turbulent flow field has been obtained through Direct
Numerical Simulation (DNS) of the NavierStokes equations; the resulting velocity field has been coupled to the Lagrangian
tracking of individual drops and to a deformation model. The drop deformation dynamics are modelled by a harmonic
oscillator forced by the local and instantaneous fluctuations of the turbulent flow (local Weber number). The temporal
evolution of the oscillator response, combined with a critical value of the drop deformation, is used to determine breakup
occurrences. The numerical simulations have been validated through comparisons with flow statistics obtained from PIV
Detailed information on the external stresses experienced by the drop can be achieved by DNS. The contributions of the
velocity fluctuations, pressure gradient and viscous shear stress to the instantaneous Weber number can be compared.
Locations of breakup events have been compared to the experimental results.
Keywords: drop, breakup, turbulence, simulation, experiments
Nomenclature
A2 normalised deformation amplitude
Bzz longitudinal velocity correlation coefficient
Ca capillarity number
d drop diameter (m)
f frequency (s')
D pipe diameter (m)
F force (N)
K prefactor of the external forcing
P pressure (Pa)
r radial direction
Re Reynolds number
7 non dimensional time
t time (s)
u instantaneous velocity (m s')
u' fluctuating velocity (m s')
U mean velocity (m s')
V drop volume (m3)
x axial direction
We Weber number
Greek letters
p damping rate (s ')
g dynamic viscosity (Pa s)
p density (kg m3)
0 angular direction
o surface tension (Nm1)
i damping coefficient
T stress (Nm2)
Subsripts
2 mode 2 of the harmonic oscillator
c continuous phase
cr critical value
p Drop
o orifice
t turbulent
x axial
Paper No
Introduction
Predicting breakup of drops and/or bubbles in
complex flows (i.e. unsteady, heterogeneous, high
concentration) is a major issue in many industrial
applications. Modelling breakup must answer two
questions: the occurrence of breakup and the distribution of
fragments resulting from breakup. These problems are
generally addressed through a critical Weber number, ratio
between the statistically averaged pressure force at the scale
of the particle (drop or bubble) that tends to deform it and
the surface tension force, which tends to restore the particle
shape to the spherical form. This modelling provides a
realistic description of breakup mechanisms in a dispersed
system at a global scale, i.e. in a dispersion volume large
compared to the particle diameter: in average, the greater
the energy input, the greater the breakup rate and the finer
the size distribution (see for instance Galinat et al. 2005).
However such an approach ignores the dynamic response of
the particle to a local instant solicitation of the flow, and
therefore cannot properly describe breakup mechanisms at
small scales, limiting the range of validity of breakup
models based on a critical Weber number in time dependent
or non deterministic flows.
The development of realistic breakage scaling laws in
turbulent flows hence requires the coupling between the
drop or bubble deformation dynamics with the local flow
stress. This type of approach has been successfully applied
to predict bubble deformation statistics in a homogeneous
turbulence (Risso & Fabre 1998) and drop breakup
probability in an inhomogeneous turbulent flow (Galinat et
al. 2007a). In the limit of weak deformations, the drop is
considered as a linear damped oscillator forced by the
external pressure field at the scale of the drop. In this frame,
the criterion for breakup is based upon maximum
amplitude of the drop deformation instead of a maximum of
the instantaneous external forcing (or critical Weber
number). In order to use this model, the characteristic time
scales of the drop (oscillation proper frequency and
damping rate) and the local instantaneous stress field in the
flow must be known. The two first parameters can be
obtained for bubbles from potential flow theory (Lamb
1932) and for drops in liquids from linearized NavierStokes
equations approximation (Miller & Scriven 1968,
Prosperreti 1980). Note that the presence of surface active
species may drastically change these two time scales (Lu &
Apfel 1991, Abi Chebel 2009). The external stress seen by
the drop can be either deduced from flow measurements
(Galinat et al. 2007b) or calculated by numerical
simulations.
In this work, based on this approach we present a numerical
study of breakup probability of oil drops in an
inhomogeneous turbulent flow of an aqueous phase. The
flow configuration considered here is the same as that
experimentally studied by Galinat et al. (2007b). It's a flow
downstream of an orifice in a vertical pipe. In the range of
Reynolds number studied, the flow is turbulent and spatially
inhomogeneous. Drops are injected upstream of the orifice
and their probability of breakup is measured in a finite pipe
volume above the orifice. The flow field is calculated via a
DNS method and drops motion is calculated from a
Lagrangian tracking method (one way coupling). Following
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
the drop trajectory, drop deformation is calculated by
solving a forced dynamic oscillator which response
provides a criterion for breakup occurrence. Statistical
averaging of breakup events leads to a breakup probability.
Numerical results are discussed and compared to
experimental data.
Figure 1: Scheme of the flow test section and domain of
simulation. The locations of the monitor points 1, 2, 3 and 4
are respectively: r/D=0.25 and x/D=0.82 1.09 1.36 and
1.68.
The flow test section investigated by Galinat et al. (2007b)
is a vertical pipe of internal diameter D=0.03m equipped
with a concentric orifice of diameter Do=0.015m. The
restriction generates a turbulent jet flow downstream.
Individual drops are injected in the flow a few pipe
diameters upstream of the orifice (where the flow is laminar
and the drop weakly deformed). The twoliquid phase
system consists of an aqueous solution of glycerin as the
carrying phase and coloured nheptane as the drop phase
(Fig. 1). Properties of the fluids are reported in Table 1.
The velocity field of the continuous phase has been
measured with highspeed PIV in a median plane of the pipe,
with uniform spatial and time resolutions of 0.55mm and
0.5ms respectively. Drop trajectories have been recorded
with the help of a highspeed camera at 500 fps. The
investigated flow region extends from the orifice to 0.06m
downstream. Drop breakup locations in the flow field
downstream of the orifice have been recorded at different
flowrates. The numerical investigation presented here is
restricted to a single case: a drop diameter d=2.4mm and a
mean velocity in the orifice Uo=0.6ms1 corresponding to
Re,=pUoDo/,u=2100.
Table 1: Fluid properties.
Phase p (kgm3) u (Pa.s) o(Nm1)
Waterglycerine 1100 4.7 103
Coloured heptane 683.7 4.5 104 23.6
Due to the inhomogeneous feature of the jet flow
downstream of the orifice, several contributions to the
external forcing applied to a traveling drop have been
considered. Following KolmogorovHinze theory
(Kolomogorov 1949, Hinze 1955), the average turbulent
Weber number is defined as:
We pu2(d)d 1)
a
where u2(d) is the mean square of the flow velocity
difference over a distance equal to the drop diameter d.
Another inertial contribution is provided by the mean flow
deceleration in the flow direction (due to the restriction),
Paper No
Ux/dlx. The Weber number associated with this inertial
stress reads:
(U 2 d3 2)
Weu = Pc Ux d 2)
( ax ) 4a
The jet flow through the orifice is also characterized by the
presence of strong mean velocity gradients in the transverse
plane (see figure 1). They induce a viscous stress which
contribution can be estimated via a capillary number:
Ca = u d 3)
c r 2t7
Both PIV measurements and numerical simulations have
shown that in the flow case studied, the contribution of
turbulence was about 2 orders of magnitude larger than that
of the mean flow (eqs 2 and 3). As a consequence, along
each drop trajectory, the external forcing of the drop
deformation (eq. 10) can be scaled everywhere in the flow
by an instantaneous Weber number:
We(t) pcau2(d) d 4)
a
Flow simulation and turbulence distribution
The flow field of the continuous phase has been simulated
by DNS of the unsteady threedimensional NavierStokes
equations in the same flow geometry and at the same
flowrate as in the experiments and geometry) of the
experiments. The simulation domain which extends towards
x/D=10 has been meshed in cylindrical coordinates with 120
nodes in axial direction x, 80 in radial direction r and 32 in
angular direction 0. The grid is equally spaced in the angular
direction while it has been stretched in radial and axial
directions near the circular orifice. At the inlet x/D=2.5, a
Poiseuille velocity profile is imposed assuming the flow to
be laminar while the outlet is modelled by a numerical
condition of free fluid exit.
66
5 5
3 43
1 1
"1
1 1
0 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3 0.4
r/D r/D
Figure 2: Velocity profiles downstream of the orifice of
the mean axial velocity Ux;
o (x/D=1/3) 0 (x/D=2/3) x (x/D=l)
(x/D=4/3) o (x/D=5/3)
Left: experimental results Right: simulation data
Simulation results have been compared to PIV
measurements (see in figure 1 the window of experimental
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
measurements). PIV measurements have been obtained by
illuminating the flow by a laser sheet crossing the central
part of the pipe. Averages have been computed assuming the
flow to be axisymmetric. In figure 2, the experimental and
numerical profiles of the average axial velocity Ux at
several locations downstream of the orifice have been
compared. Both results are scaled by the average velocity U
flowing through the pipe. These profiles show the gradual
spreading of the jet and the existence of a recirculation zone.
It is important to remind that small discrepancies at large r
lead to strong modifications of the axial flow rate close to
the centre of the pipe.
0.U 2 .2
0 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3
r/D r;D
Figure 3: Profiles of the r.m.s. of the axial
fluctuations 2 downstream of the orifice.
0.4
velocity
o (x/D=1/3) 0 (x/D=2/3) x (x/D=l)
(x/D=4/3) o (x/D=5/3)
Left: experimental results Right: simulation results
General agreement has also been obtained when comparing
the profiles of the r.m.s. of axial velocity fluctuations
(Figure 3). The spatial repartition of the r.m.s. of axial
velocity fluctuations (Figure 4) shows that the DNS is able
to predict the structure of the flow although the intensity is
slightly overestimated in the simulations.
In the experiments, the presence of strong fluctuations close
to the orifice (x/D=1/3 and 2/3) is a possible signature of
selfsustained oscillations of the jet core at low frequency.
Temporal energy spectra of the axial velocity (figures 5ab)
shed some light on this feature. The energy spectra of the
axial velocity have been scaled by the square of the average
axial velocity. The experimental spectra extracted from PIV
measurements have been smoothed applying the Welch
windowing method. The experimental spectrum at point 1
has a clear bump around f4Hz which disappears
downstream (point 2). Due to the exact steady Poiseuille
flow inlet condition, such a largescale oscillation of the jet
does not arise in the simulations. This has an impact on the
velocity fluctuations close to the orifice and consequently
the profile of the average velocity is stiffer in the
experiments than in the simulation. Finally, this peculiar
phenomenon is damped and a good agreement is achieved
downstream x/D= 1.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
0 0.5
0 05
Figure 4: Spatial distribution of the r.m.s. of axial velocity
fluctuations x .
Left: experimental data Right: simulation results
Although the turbulence is inhomogeneous and anisotropic,
the energy spectra from PIV measurements exhibits a clear
inertial subrange characterized by the classical 5/3 power
law. The power law of the inertial subrange in the energy
spectrum is also confirmed in simulations up to the
frequency 2.102s1. For higher frequencies an increase of the
slope is observed and may be related to the progressive
coarsening of the grid filtering the smallest scales. In the
region of interest where we want to compare experiments
and numerical prediction of drop breakup, the spatial and
temporal resolutions are adequate to capture the important
features of the turbulence dynamics.
The drop deformation is not solely the consequence of the
temporal forcing. Spatial stress distribution over scales
comparable to its diameter forces the drop to deform and
might lead to breakup. Therefore comparing spatial
correlations of the velocity fluctuations at the drop scale is
also an important step for DNS validation. The goal here
is to verify that the spatial structure of turbulence
encountered by the drop is well reproduced in the
simulations.
The drop diameter is 2.4mm. At point 2, the average axial
velocity is U,=0.25ms1 and the characteristic frequency at
this scale is the order of U,/d =102s1. In this range, the
inertial subrange is still accurately solved by DNS
simulations. Further downstream, the grid has been
progressively coarsened in the axial direction. In order to
check that accurate statistics of velocity fluctuations at the
drop scale are simulated in the DNS, the structure function
of velocity increments has been computed and compared to
experiments. For an inhomogeneous flow, Risso & Fabre
(1998) showed that the structure function u2(d) in axial
direction could be expressed as:
,(x) [ u,(x+) )
where B U(dxx) is the longitudinal cot2 Bf c:(d)
where B (dx) is the longitudinal coefficient of correlation:
B (dx)
u (x)u (x+ dx)
u ()x 2 +d2x)2
Figure 5: Energy spectra of the axial velocity: simulations
(thinner line) and experiments (ticker line). Top figure: point
1 and bottom figure: point 2 (see figure 1 for the location of
points 1 and 2).
For homogeneous and isotropic turbulence the structure
function & 2(dx) in the inertial subrange follows a 2/3
power law of increasing separation distance dx. Based on
figure 6, it is clear that the drop experiences velocity
fluctuations within the inertial subrange. This is true all the
way through the experimental window (from point 1
towards point 4). Although the geometry produces a
strongly inhomogeneous flow at large scale, the turbulence
at the drop scale follows the universal Kolmogorov energy
cascade. Simulation results closely follow experimental
data at point 1 whereas the 2/3 scaling is preserved at point
4. Point 4 is located at 0.05m from the orifice. The slight
discrepancy is related to the local mesh width which is of
the order of the drop diameter. At point 1, the ratio
between the width of the mesh cell and the drop diameter
is 0.7 whereas it increases to 1.4 at point 4.
Paper No
Paper No
101
213
10, ,
0
o
1 11 , , I 10 _11 , ,
10o 10o
Figure 6: Normalized structure function S5 2(dx d).
(o : PIV experiments ; x : DNS results).
Left : point 1 ; Right : point 4.
In figure 7, we also compared numerical and
spatial distributions of the average turbulent W
We=pcau(d)d/ ()., The agreement is
qualitatively and quantitatively when x/D>
concluded that temporal and spatial resolution
are adequate to reproduce the important fei
flow experienced by deformable drops. 7
increment at the drop size is progressively inc
the orifice section and reaches its maximum a
the experimental data, the presence of i
numbers close to the orifice is related to
oscillations of the jet (see comments on figure
10 1.4
1.2
8
Figure 7: Measured (a) and simulated (b) aver
Weber number. White crosses correspond to loc
drop breakup occurs.
experimental
Teber number
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Simulation of drop deformation and breakup
statistics
Trajectories of individual drops have been simulated in a
Lagrangian framework by solving the force balance (eq. 9)
acting on an isolated drop.
du
p d=
p dt
Forces F experienced by the drop are: the drag force, the
buoyancy, the lift and the addedmass forces and the effect
of the dynamic pressure. In the simulation of drop
trajectories, we neglect twoway coupling (dilute flow)
assuming point drop. The instantaneous velocity u and the
velocity gradients involved in the force balance are provided
through a local interpolation of the DNS flow field at the
drop location. The system of ODEs (eq. 9) has been
efficiently solved with a 4 steps forward RungeKutta
scheme. An example of calculated trajectory is illustrated in
figures 1 and 8.
good both The shape of a drop can be described as a sum of spherical
1. It can be harmonics (see examples in Risso, 2000). Each of them is
s of the DNS associated to a dynamical mode which has its own
atures of the eigenfrequency and damping rate.
The velocity Provided the Weber number is not too large, Risso & Fabre
:reasing after (1998) showed that largeamplitude deformations that are
t x/D=1.5. In responsible for breakup in a turbulent flow essentially imply
moderate We axisymmetric spherical harmonic number 2, which
large scale corresponds to a lengthening or flattening of the drop shape.
4 above). The drop deformation can thus be characterized by the
amplitude A2 of mode 2 and calculated by computing the
response of a linear oscillator to the instantaneous turbulent
forcing. Using the diameter d of the drop as length scale and
14 1 Trru, as time scale, this model reads in dimensionless form:
d2A, dA p&2dd
d 2 +2 +2 = K = KWe(d,x,r,,7t), 10)
12 2 a 2 O"
where the tilde symbol denotes a nondimensional quantity
and =[3 i'rr is the damping coefficient. The lefthand
10 side describes the dynamics of the drop shape. The
frequency f2 and the damping rate 62 for any arbitrary
8 Reynolds numbers were first derived by Miller & Scriven
(1968) and by Prosperetti (1980). Here, we used the
asymptotic expressions proposed by Lu & Apfel (1991) for
6 lowviscosity fluids, which are valid for the present case:
f2=43Hz and Pl=32s1. The righthand of eq. (10) stands for
the turbulent forcing. The Weber number We is computed by
4 DNS from the velocity fluctuation at the scale of the drop,
which has been shown to be accurate from the examination
2 of the structure functions (fig. 5). It is a Lagrangian quantity
since it is considered at the location (x,r ) where the drop is
at the instant t K is an unknown prefactor which is
0.5 introduced to couple the dynamics of the drop shape with
the turbulent forcing. Eq. (10) has been solved numerically
age turbulent by using a fourstep forward RungeKutta scheme and by
nations where considering that each drop was initially spherical (A2=0).
Because eq. (10) is linear the value of K can be set
arbitrarily to unity without loss of generality (A2 is simply
replaced by A2/K)of.
Paper No
Figure 8 shows an example of simulation: the radial drop
location r, the computed deformation A2/K and the
corresponding turbulent forcing We are plotted versus the
axial location x/D. It is worth pointing out the differences
between the turbulent forcing and the deformation that
results from it. The evolution of A2 is smoother than that of
We, which indicates that the shape dynamics preferentially
respond to frequencies close to f2 and filter out higher
turbulent frequencies. Also, the peak of A2 does not coincide
with that of We. As already shown by Galinat et al. (2007b),
a breakup criterion based on the value of the Weber number
is therefore not relevant. Following their recommendations,
we will assume that the breakup occurs when the computed
deformation A2 exceeds a critical value A2,,. It is important
to note that A2,, is the only adjustable parameter of the
model.
Trajectories and deformations of 450 drops have been
simulated. The breakup locations, which have been reported
in figure 6, seem to be in good agreement with the
experiments. Rigorous comparisons require to compute
statistics. We define the probability P(x/D) that a drop
breaks while its axial position belongs to the interval [O,x/D].
P(x/D) starts from zero at the location of the orifice (x/D=0)
and then characterizes the increase of the breakup
probability with the distance from the orifice.
Figure 8 shows axial evolutions of P computed for four
various critical deformations: A2,jK=14, 15, 16 and 17. The
experimental results appear to lie nicely in between the
simulations corresponding to A2,jK =15 and 16. Regarding
the moderated number of samples considered in the
experiments (70), this agreement is remarkable and proves
that the present elementary dynamical model of interface
accounts for the right physical mechanism: the drop shape
responds as an oscillator defined by its frequency and its
damping rate. Coming back to the Weber map shown in fig.
6, we note that some breakup events occur in regions where
the average turbulence is weak, which confirms that the
time response of the interface really needs to be accounted
for.
02
0 0.5 1 15 2
o 0.os 1 15 2
x/D
20
0 '
0 0.5 1 15 2
xflD
0
0 0.5 1 15 2
x/D
Figure 8: Comparison of the simulated drop trajectory (top
figure), its normalised amplitude of oscillation Ai2
(middle figure) and the instantaneous turbulent Weber
number acting on the drop (bottom figure).
From video analysis of the drop shape at the instant of
breakup, Galinat et al. (2007a) estimated that the critical
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
deformation A2,, was approximately unity. Considering that
the best agreement between the experiments and the present
simulations is obtained for A2,/K =16, we find a coupling
factor K=0.06. This estimation is larger than the value K
=0.04 obtained by Galinat et al. (2007a) by solving eq. (10)
in following conditions:
(1) a forcing term We(t) obtained from PIV measurements in
a single radial plane,
(2) drop motions calculated by considering a constant drift
velocity relative to the liquid.
The value of K is thus not accurately determined. However,
once K has been fixed empirically by comparison with
experiments, the present numerical technique allows to
predict the local breakup probability in a strongly
inhomogeneous flow.
05 1 15 2
x/D
Figure 9: Breakup probability versus the distance from
the orifice evaluated from measurements (circles) and
calculated from simulations assuming critical deformation
values between 14 and 17 (continuous lines).
Conclusion
Drop breakup in the turbulent inhomogeneous flow that
develops downstream of an orifice has been numerically
modelled by coupling DNS simulations of the continuous
phase, Lagrangian droplet tracking and a dynamic model of
breakup. This model is an adaptation of KolmogorovHinze
theory of turbulent breakup to the RayleighLamb theory of
drop oscillations. Comparing to PIV measurements, DNS
results have been demonstrated to provide a reliable
prediction of the turbulent flow field and its turbulent
statistics at the drop size scale.
The experimental breakup positions have been correctly
predicted by adjusting only one single parameter in the
numerical model: the value of the critical deformation A2cr
for breakup. The capability of the forced harmonic
oscillation model to provide accurate description of
turbulent drop breakup in a nonhomogeneous turbulent
flow field has been tested and discussed.
These simulations can be used to study the response of
different drop diameters and also different liquid/liquid
systems (i.e. different damping rates and frequencies). It can
also be used in other geometries of inhomogeneous
turbulent flows.
Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
References
Abi Chebel N. Interfacial dynamics and rheology of an
oscillating drop at high frequency, PhD thesis, Univ. of
Toulouse, INPT (2009).
Galinat S. Experimental study of drop breakup in a
turbulent flow. PhD thesis, Univ. of Toulouse, INPT (2005).
Galinat, S., Masbemat, O., Guiraud, P., Dalmazzone, C. &
Noik, C. Drop breakup in turbulent pipe flow downstream
of a restriction. Chem. Eng. Sci., Vol. 60, 65116258 (2005).
Galinat, S., Risso, F., Masbemat, O. & Guiraud, P.
Dynamics of drop breakup in inhomogeneous turbulence at
various volume fractions. J. Fluid Mech., Vol. 578, 8594
(2007 a).
Galinat, S., Garrido Torres, L., Masbemat, O., Guiraud, P.,
Risso, F., Dalmazzone, C. & Noik, C. Breakup of a drop in
a liquidliquid pipe flow through an orifice. AIChE J., Vol.
53 (1), 5668 (2007 b).
Hinze, J.O. Fundamentals of the hydrodynamic mechanism
of splitting in dispersion processes. AIChE J., Vol. 1,
289295 (1955).
Kolmogorov, A.N. On the disintegration of drops in a
turbulent flow. Dokl. Akad. Nauk. SSSR, Vol. 66, 825828
(1949).
Lamb, H. Hydrodynamics. Cambridge University Press
(1932).
Lu, H.L. & Apfel, R.E. Shape oscillations of drops in the
presence of surfactants. J. Fluid Mech., Vol. 222, 351368
(1991).
Miller, C.A. & Scriven, L.E. The oscillations of a fluid
droplet immersed in another fluid. J. Fluid Mech., Vol. 32,
417435 (1968).
Prosperetti A. Normalmode analysis for the oscillations of a
viscous drop in an immiscible liquid. J. de Mecanique, Vol.
19, nl, (1980).
Risso, F. & Fabre, J. Oscillation and breakup of a bubble
immersed in a turbulent field. J. Fluid Mech., Vol. 372,
323355 (1998).
