7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Influence of the history term in a Lagrangian method for oilwater separation
D.F. van Eijkeren* and H.W.M. Hoeijmakers*
Faculty of Engineering Technology, University of Twente, PO Box 217, 7500 AE, Enschede, The Netherlands
d.f.vaneijkeren@utwente.nl and h.w.m.hoeijmakers@utwente.nl
Keywords: Lagrangian simulations, history term, swirling flow, oilwater separation
Abstract
In the oil industry, both settling tanks and swirling flow separators are used to provide separation of phases. In the
present paper the motion of oil droplets is investigated numerically for the case of swirling flow fields as well as for
the case of settling in quiescent fluids. During phase separation the oil droplets are subjected to various forces, i.e.
drag force, buoyancy force, pressure gradient force, virtual mass force, Saffman and Magnus lift force and the history
or Basset force. In oilwater separation, the density of the continuous phase is of the same order of magnitude as that
of the dispersed phase. This implies that, for a relatively dense continuous phase, the history term cannot be neglected.
Several experimental and numerical data, e.g. Mordant and Pinton (2000); Mei et al. (1991); Loth and Dorgan (2009),
show that using the classical Basset kernel for the history term overestimates this unsteady force significantly.
In the present research, Lagrangian simulations have been performed using, for the history term, variants of the kernel
proposed by Mei et al. (1991). Obtained results are compared with experimental and numerical data. The influence
of the history force using various kernels has also been investigated for a generic swirling flow field, for the moment
without accounting for effects of turbulence. The results show that inclusion of the history term leads to an increase
of the time required for particles to reach the axis of rotation. This implies that in numerical methods for swirling
flow separators the history term should be included. For the case of solidbody rotation, present results have been
compared with experimental results by Candelier et al. (2005). This comparison shows that, in order to account for
the forces adequately, additional effects have to be taken into account for solidbody rotation flow fields, as also
suggested by Candelier (2008), who also accounts for nonBasset type of history effects.
Nomenclature
Roman symbols
C coefficient for force abbreviated in subscript ()
c empirical parameters in history term ()
D particle or droplet diameter (m)
F force(kgms 2)
f force per unit mass (m s 2)
f drag factor ()
fH history factor ()
g gravitational constant (m s 2)
g gravitational acceleration vector (m s 2)
I moment of inertia (kg mn2)
K dimensionless kernel in integration ()
m mass (kg)
R radius (m)
Re Reynolds number based on D()
t time (s)
T dimensionless time ()
T Torque (kg n2s2)
u flow velocity vector (m s 1)
v particle velocity vector (ms 1)
w v u (ms 1)
x coordinate in space (m)
Greek symbols
a dimensionless particle rotation ()
3 dimensionless flow rotation ()
F Circulation (m2 s 1)
p dynamic viscosity (kg m 1 s 1)
v kinematic viscosity (m2 s 1)
S arbitrary vector variable
p density (kg m 3)
p density ratio ()
(a stress tensor (N m 2)
T characteristic and integration time (s)
4 arbitrary forcing vector
w rate of rotation (s1)
Subscripts
c continuous
d diffusive
D nondimensional based on particle diameter
H history
v referring to viscous decay time
p particle, either solid or fluid
R nondimensional based on initial radius
t terminal
V referring to relaxation time
w referring to rotation Reynolds number
0 at initial time to
Superscipts
T transposed
t at discrete time
Introduction
In oil industry, water is used to maintain pressure in oil
wells during the process of extraction of oil. Therefore, a
mixture of oil and water will be produced during the ex
traction process. It is required to separate this oilwater
mixture, not only to obtain the desired oil, but also to
be able to re inject the water into the oil field, or to dis
pense it in the environment. Separation of phases can
be divided in primary, or bulk, separation and secondary
separation. A mixture with a very small amount of oil,
as well as a mixture with a very small amount of wa
ter are produced during the initial bulk separation. The
secondary separation of oil from produced water, is also
known as produced water treatment. Tiny oil droplets
that remain in the water after bulk separation will need
to be removed in order to obtain sufficiently clean water.
This requires more advanced separation methods.
Separation of oil and water is most often based on the
density difference between oil and water. A widely used
method based on the density difference is gravity set
tling. A mixture is put in a settling tank and after some
time oil will float on top of the water and can be ex
tracted from the original mixture. The method of sepa
ration using gravity settling is most often implemented
as a batch process. Moreover, the settling times during
produced water treatment would be days or even longer.
Therefore, more advanced methods have been developed
to obtain separation of phases.
An important quantity to estimate the speed of set
tling, and thus to identify improvements for separation,
is the terminal velocity vt of droplets. As a first approx
imation the droplet will be approximated by a spherical
particle. Another approximation for an initial estimate
of vt is that the flow around the particle will be Stokes'
flow. Therefore the drag will be Stokes' drag. The ter
minal velocity vt is then:
Dg (1
18V,
g (1
P
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
with the particle relaxation time TV and density ratio p
defined as
=pD2
18V,
and p
respectively.
The most important parameter influencing the settling
time is the droplet diameter D. Coalescence of droplets
causes a huge improvement in settling time. However,
in dilute mixtures with very small droplets, coalescence
of droplets is not occurring very often. Another param
eter in equation (1) that can be influenced is the gravita
tional acceleration. Centrifugal forces create an artificial
gravitational field, stronger than the natural gravitational
field. This is used in swirling flow separators. Swirl is
generated by either an internal swirl element in the pipe
or a tangential inlet to the pipe. Due to the centrifugal
force the oil droplets will move towards the axis of rota
tion. A major advantage with respect to settling tanks is
that a continuous, rather than a batch, separation process
is possible.
To control the efficiency of swirling flow separators,
it is not only necessary to be able to compute the flow
field, but also to be able to compute the droplet trajec
tories. A small change in inward velocity of a particle
could imply, for example, that the separator should be
much longer than foreseen. Therefore, the equations of
motion for a particle have been investigated in order to
be able to predict particle behaviour adequately. These
studies involve both particle behaviour in settling tank
conditions and particle behaviour in a generic swirling
flow field. The equations of motion involve various flow
phenomena resulting in forces on particles.
One of the forces involved in particle tracking cap
tures the effects of the history of the particle movement.
This history term is often neglected, not only because the
term is often deemed very small, but more importantly
because calculation of the force involves an integral over
the time from the start of the motion. Such an integral
is computationally very expensive, especially if the time
steps are very small. However, present research shows
that the history term cannot be neglected in the case of
oil droplets in water. This force will generally increase
the predicted settling time, both in conventional settling
tanks and in swirling flow separators. This will result in
either more, or larger, settling tanks, or longer swirling
flow separators.
First the equations of motion for a particle will be de
scribed. Subsequently solutions for these equations of
motion will be compared with experimental data from
Mordant and Pinton (2000). The data will be used to
compare the effect of several kernels that are used in the
history term. Empirical coefficients that appear in the
kernels have been varied to obtain a proper behaviour of
the history term. Subsequently results will be presented
of a comparison of numerical results with numerical and
experimental data obtained by Candelier et al. (2005);
Candelier (2008). Finally results will be shown for a
generic swirling flow field and conclusions will be pre
sented about the history kernel and its importance for
water treatment using swirling flow fields.
Equations of motion for a particle
The motion of a particle is dictated by the velocity due
to forces acting on that particle. Newton's second law
states that d = F. The change in time of the posi
tion of a particle is then the velocity of the particle. The
torque on a particle is related to the rotation of the parti
cle, coupled by its moment of inertia: dIP T. Mak
ing the assumption that no mass transfer is occurring and
the moment of inertia is constant, the differential equa
tions for the position, velocity and rotation of a particle
can be written as:
dv
Spdt =m volumee + fVstress d+ rag + fSaffman
+fMagnus + faddedmass + history + unknown) (3b)
dwp T
dt LI
(3c)
where the forcing term per unit mass is split into several
additive contributions. The last forcing term captures
all forcing terms currently not taken into account. Us
ing known relations for the forcing terms the equations
of motion can be rewritten. The unknown forces will be
neglected. However, unknown forces can have a major
impact on the motion of particles. This is for example
observed for a particle settling along the axis of rotation
in a solid body rotation, e.g. see Candelier et al. (2005);
Candelier (2008).
For the derivation of the forces, it will be assumed that
oil droplets can be approximated by solid spherical par
ticles. As the particle Reynolds number will generally
be low for settling conditions, deformations due to parti
cle induced flow will be negligibly small. The particles
used in settling experiments are solid spheres, therefore
no assumptions have to be made with respect to geome
try or substance.
The volumetric forces such as gravity are written as:
volume = f (4)
The forcing due to the divergence of the stress of the
undisturbed flow field at the location of the particle V o
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
has been rewritten using the NavierStokes equations for
oneway coupling in nonconservation formulation:
Vao 1Du
fVstress f, (5)
PP Dt
with the material derivative defined as:
D a
D +u.V. (6)
Dt Ot
The effect of the drag is written as a function of the
drag factor f, the particle relaxation time TV and the rel
ative velocity w.
fdrag  (7)
TV
A larger particle relaxation time implies more slip and
therefore smaller influence of the drag on the particle
acceleration. The drag factor is related to the drag co
efficient by Stokes' drag. It shows the influence of non
potential flow effects on the drag:
Re C'
f = D 2 6 fStokes 1, (8)
24 CD, Stokes
with the Reynolds number Re:
Dlwl
Re (9)
Vc
In the present research the relation for the drag factor
found by Clift and Gauvin (1970) is used, although the
relation for the drag factor found by Schiller and Nau
mann (1933), a good approximation for Re < 800 ac
cording to Crowe et al. (1998), can be used for most
cases considered. The drag factor, from the drag coeffi
cient found by Clift and Gauvin (1970), is:
f 1 + 0.15Re0 68
0.0175Re
1 + 4.25 x 104Re 116 (10)
Figure 1 shows four common expressions for the drag
coefficient. It can be noted that Stokes drag underesti
mates the actual drag for all Reynolds numbers, there
fore f > 1. However, for Re < 1 Stokes' drag is a very
good approximation. The drag factor can be corrected
for the case of a nonsolid particle, i.e. a droplet or a
bubble. For the present research such a correction is not
used, as the viability of such a correction for the non
Stokes regime does not appear to be proven. Moreover,
in case of contaminated bubbles such a correction would
underestimate the drag force on a bubble significantly.
Also, in the case of oil droplets in water, the viscosity of
the droplet is generally larger than the viscosity of water
and therefore the correction is quite small.
The Saffman and Magnus lift forces capture phenom
ena related to rotation of particle and flow field. Figure
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
10
101
101 0'5
100 10
Re
Figure 1: The drag coefficient CD for several expres
sions found in literature.
FL,S
U
Figure 2: Saffman force effect(left) and Magnus force
effect(right).
2 shows the two flow phenomena taken into account by
these forces. The Saffman lift captures the effect of a
higher fluid velocity on one side of the particle with re
spect to the other side of the particle due to a flow field
velocity gradient. This effect causes a pressure differ
ence over the particle resulting in a lift force that de
pends on the direction of the rotation and the direction
of the relative velocity. The Magnus lift force represents
the effect of the flow field being curved by a rotating
particle. Rotation of the particle causes a change in an
gle between the flow ahead of the particle and the flow
in the wake of the particle. This is caused by a force
of the particle on the flow and vice versa. Crowe et al.
(1998) use the relative rotation for the Magnus lift, while
Loth and Dorgan (2009) use the absolute rotation of the
particle. For the cases studied, both methods have been
investigated, and differences turn out to be small in cur
rent research. An approach using the absolute rotation
will therefore be used.
The Saffman lift force per unit mass is:
9.66 n^
fSaffman LS x w6, (11)
sann Di) V~ I W 0
with the rotation of the flow field:
W = VX U.
The Saffman lift coefficient CL can be written as, see
Mei (1992); Crowe et al. (1998):
L,S o0.33143 (1
CL, 0524 Re
0.0524 7Re
Ro+ \
Re
+e 10
Re < 40,
Re > 40,
where
D I
= 0.005 < 3 < 0.4, (14)
Iw
This is a correction for the lift force which in the limit
of low Reynolds numbers, for which the influence of
B becomes negligible, reduces to the relation found by
Saffman (1965), and for high Reynolds numbers be
comes the relation from numerical results by Dandy and
Dwyer (1990).
The Magnus lift force per unit mass is:
31
fMagnus 4 3 CL,MWp X W. (15)
The Magnus lift coefficient CL,M is a Reynolds number
and spin correction on the Magnus lift as derived by Ru
binow and Keller (1961). Loth and Dorgan (2009) use
a relation that gives the lift coefficient in the same order
of magnitude as measurements. This relation is:
CL = 1 +0.15) tanh (0.18Re)
0.15tanh(0.28(a 2))tanh (0.18vRe), (16)
with
c =D p (17)
w
When a particle is moving with respect to the flow,
a similar volume of fluid has to be moved in opposite
direction. In case of an acceleration of either the contin
uous phase, or the dispersed phase, a resulting force will
remain. This effect is written as:
C 1Du
added mass Cam
p Dt
dv
dt)
The expression for the added mass force coefficient pro
posed by Auton et al. (1988) is Cam = 1, and is ac
curate for a fairly large range of Reynolds numbers, see
Loth and Dorgan (2009). There has been some debate
whether or not the use of the material derivative is cor
rect, but general consensus appears to be that the ma
terial derivative should be used. However, in case the
3
+ L,M
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
particle velocity does not differ too much from the fluid
velocity, the derivative of the flow velocity along the par
ticle trajectory is a good approximation of the material
derivative.
The history force will be presented in more detail after
introducing the equations of motion. It can be written as
an integral over time, having the dimensionless kernel
K (t, T), with T the integration time variable:
fh T K dT.
C
Using the relations for the forcing terms, neglecting
the unknown forces, and rewriting parts for convenience,
the equation of motion for a particle used in present re
search is written as:
dv
(p+ Cam)
dt
SPf
(P 1) f w
TV
Du 3
+ (1 + Ca,) Du+ 3CL,Mp X W
9.66 L
DCLS1 W X Ww
7nD V FW o
t
TV
din
K (t, T) dw (20)
dT
The history term is related to Stokes' second problem
for the boundary layer on a flat plate oscillating in the
direction tangential to the plate. The drag force in equa
tion (7) is based on a developed flow field around the
particle. However, the boundary layer will only be fully
developed after a certain time span. To correct for this
phenomenon the history term will have to be included,
describing the noninstantaneous effects. The kernel that
was originally derived for this force, also called the Bas
set force, will be called the Basset kernel and can be
written as:
d1 Td 1
KBasset (t, T) = (21)
with the diffusive timescale:
The kernel is based on the solution for the case of
constant acceleration in the low Reynolds limit. This
solution involves the derivative of an error function. For
a variable acceleration infinitesimal steps can be com
bined to the resulting kernel with a 1 behaviour.
It was observed that the decay of influence of the ac
celeration over time was more rapid for the long time
100
(tt)/d []
Figure 3: Reproduction of figure in Loth and Dorgan
(2009) of the history kernel with cl = 2.5 and
C2 0.2 for three Reynolds numbers. The
history time TH is as defined in equation (25).
behaviour than is predicted by the Basset kernel. There
fore alternative kernels have been proposed. An alter
native kernel, proposed by Mei et al. (1991); Mei and
Adrian (1992), has a square root decay for short time be
haviour, and a more rapid quadratic decay for long time
behaviour:
K (t, T)
1 +
KBasset (t, T)
( Re
KBa ,e(t,T))
1] cl (23)
TW c
with the history factor fH as:
fH = (0.75 + caRe)3,
and cl and c2 dimensionless empirical coefficients that
have been taken cl = 2.5 and c2 0.2 by Loth and
Dorgan (2009), using several sources of settling experi
ments. However, Mei and Adrian (1992) originally pro
posed cl = 2 and c2 0.105 on semianalytical basis.
In the line of Loth and Dorgan (2009) a history time is
now defined to indicate the distance in time from time t
at which the influence of the square root behaviour is
equal to the influence of the quadratic behaviour, i.e.
when the denominator in equation (23) is equal to 2.
Therefore:
H fH Re2 (25)
7 Re
The history time is only influenced by the choice of the
parameter c2, the Reynolds number Re and the particle
diffusive timescale Td. In figure 3 the history kernel is
presented including the history time. For the limiting
case of Re = 0 the history time would be infinite, as
no transition to the quadratic decay regime will occur.
The parameter cl does not influence TH as defined in
equation (25). However, cl does influence the size of
the transition region of the curve between the square root
and quadratic decay regime, i.e. a larger cl results in a
larger transition region.
The kernel was subsequently refined by Lawrence and
Mei (1995), leading to a drag factor depending kernel,
with history factor fH:
3
fH (t, T) 3 (f (t) + Re (t) f' (t)) f' () (26)
4
For low Reynolds numbers Lawrence and Mei (1995)
used the drag factor
f =1 + Re. (27)
16
It can be derived that using this expression for the drag
coefficient in equation (26), the drag coefficient based
history factor and the empirically determined history
factor in equation (24) only slightly deviate in values up
to order 0 (Re).
where f' is the derivative of f to Re. Lawrence and
Mei (1995) also found that the quadratic time compo
nent can be described by an integral of the velocity from
T to t. In present research, after comparison of numer
ical results with experimental data, the original expres
sion for the history kernel is here adapted with an added
drag factor correction leading to:
K (t, T) 4fH (t, T)
7T
x Re (t)j(t, + T (t,r) (28)
with the history factor fH as defined in equation (26)
and dimensionless time T as:
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
there will be some history prior to to, the time at which
initial information is known. Therefore the integral is
split into a part for which data can be obtained, either
exactly or by interpolation, and a part for which data is
not available:
t
dw
00O
to t
L dw f dw
K dT dr+ K dT. (30)
dT J dT
CO tO
Subsequently the first part of the integral can be
rewritten using integration by parts, as
to
K dT
j dr
00
to
f dK
[Kw] o, J w dT.
The integral in the right hand side is still unknown.
However, the derivative of the kernel has a larger rate
of decay than the kernel itself, while relative velocities
can be expected to be near zero. The kernel times the
velocity at minus infinity will be zero, and therefore the
integral over the entire time domain can be written as if
the particle obtains its velocity instantanuously at t = to:
t
f dw
K dr
00
SKdw dT + Kowo.
S dr
The torque on the particle is also required for solu
tions of the equations of motion of the particle. The
torque on a particle tends to force the rate of rotation
of the particle to be the same as half the rate of ro
tation of the flow. In Crowe et al. (1998) a relation
for the torque vector is given, and subsequently a rela
tion for the torque in a quiescent fluid depending on the
Reynolds number is provided. It will be assumed that
the Reynolds number correction will be approximately
valid for the relation with a fluid flow, leading to:
SRe (') ,
Tdf (T) d .
The expressions in equations (28),(29) differ from the
equation proposed by Lawrence and Mei (1995) due to
added drag coefficient corrections and a smooth switch
from the square root to the quadratic decay regime.
These corrections are based on a comparison of numeri
cal results with experimental results of Mordant and Pin
ton (2000), this will be shown in more detail in the next
section.
It must be noted that the lower limit of integration for
the history term is minus infinity. If the particle and fluid
are initially, at time to, at rest, this is not a problem as
the contribution prior to to is zero. If this is not the case,
T = 2.01/D3 cwe
where:
1 )
p (1 + 0.201 Re,
p D2
4)P 7we[ D2
Analytic solution of the equations of motion is only
possible in very specific situations and often the assump
tion of Stokes flow. For a more generic case a numerical
integration scheme is required to obtain an approximate
particle trajectory. In the present research two numerical
methods have been used. An explicit algorithm using the
classical RungeKutta scheme, which is globally fourth
order accurate in time, is used for most calculations. It
is especially useful for the comparison of the numerical
solution with experimental data and fully resolved flow
solutions. The RungeKutta scheme has good capabil
ities for adaptive timestepping which is used for some
of the calculations. For particle tracking in a generic
swirling flow field, an implicit CrankNicholson scheme
is used, which is globally secondorder accurate in time.
The major advantage is that in this scheme there is no
strict limitation to step size for a stable solution. How
ever, if flow time scales are small, solutions generally do
not show these high frequency effects. A major disad
vantage of an implicit scheme is the requirement to solve
a system of equations. In the case of the present equa
tions, this will involve a nonlinear system of equations.
However, iterative methods can generally approximate
solutions quite efficiently with only a few iterations.
The history term requires a numerical integration
scheme. A Simpson quadrature method is used to ob
tain a proper order of accuracy. However, one of the
major problems for a standard integration scheme is the
singularity at T t in the kernel for both the Basset ker
nel and the alternative kernel. One of the solutions for
this problem, which is only possible for the Basset force,
is to write the entire integral as a fractional derivative,
see e.g. Bombardelli et al. (2008); Coimbra and Rangel
(1998). The fractional derivative makes use of the divi
sion by the square root, and therefore cannot be used for
the kernel involving quadratic decay. Another solution,
also only applicable to the Basset kernel, is to assume a
constant acceleration in the final time step and use the
analytic integral for the Basset kernel for this situation.
For the alternative kernel it can be assumed that in the
limit of very small time steps (t < t/ which
implies that the kernel reduces to the Basset kernel in the
final step. This makes the second method, in which an
analytic solution for the Basset kernel is used valid for
the calculation of the history term. Therefore the method
used in the present research for calculation of the history
term will have an altered integral:
dw
/K d dr
j dr
CO
t t
f dwi
lirn Kdw dT
StRo J dT
to
dw
dt
+ d  .KBasst drT+ Kowo. (35)
2 t6t
Settling sphere
In 2000 Mordant and Pinton (2000) published experi
mental results of experiments of spheres settling in qui
escent water. These experiments have been used to adapt
the kernel in the history term in such a way that par
ticle behaviour is predicted fairly accurate for a large
range of Reynolds numbers. The results show that pre
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
dictions without a history term overestimate the initial
settling velocity and underestimate the time required for
the sphere to reach its terminal velocity. On the other
hand, predictions using the Basset kernel significantly
underestimate the long term velocity and overestimate
the time required for the sphere to reach terminal veloc
ity. Both predictions sometimes overestimate the veloc
ity of the sphere in part of the time domain for which
the Reynolds number is high. Therefore a kernel should
not only show a more rapid decay than the Basset kernel,
but it should also result in larger values than the Basset
kernel in specific situations.
Mordant and Pinton (2000) calculated the history
force for a specific case. The calculated force is the re
sult of the force determined by the acceleration of the
particle, with known forces such as the drag and buoy
ancy substracted. Figure 4 shows the resulting history
force. It includes the history force caclulated by Mor
dant and Pinton (2000) using the Basset kernel. It was
observed that while the peak force is of the proper or
der of magnitude, however, the long term decay of the
Basset kernel is much too slow. This case has been
used to compare the numerically obtained history force
with experimental data. It has been assumed in the nu
merical simulations that the experiments have been per
formed in Lyon. Using the altitude and latitude of Lyon
the gravitational acceleration has been calculated to be
g z 9.807 [ms 2]. Mordant and Pinton (2000) de
termined water conditions to be conditions at 25 [ C].
At that temperature, density and kinematic viscosity are
pc 997 [kgm 3] and v = 8.9 10 [m2 s ], re
spectively.
The numerically determined history force using the
kernel as defined in equation (28) is also shown in fig
ure 4 for a visual comparison. The general trend of the
force is predicted correctly, although the peak value of
the force is a bit lower than that of the experimentally
determined force. In figure 5 the history force predicted
using several kernels are shown. While the Basset ker
nel overestimates the force significantly, the kernels pro
posed by Loth and Dorgan (2009) and the one proposed
by Mei and Adrian (1992) predict a lower history force
in the part of the domain where accellerations are rela
tively large.
Another quantity that can be used to compare numeri
cal results with experimental results is the velocity. Usu
ally the velocity time history is provided in experimental
data found in literature. Mordant and Pinton (2000) pro
vide velocity histories for several cases. The first case
considered here is the case for which the history force
was determined. Figure 6 shows the velocity history
for this case including the velocity predicted using the
kernel of equation (28). While the numerical prediction
slightly overestimates the velocity reaching terminal ve
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
z ........4 .. ..... . ......" ............t......
0 0.1 0.2 03 04 0.5 0.6
time (s)
Figure 4: Experimentally determined history force, ir
regular solid black line, and predicted history
force using the Basset kernel, dashed line, by
Mordant and Pinton (2000). Prediction using
the kernel of equation (28) is added for com
parison, solid blue line. pp 7850 [kg m3],
D 1 103 [m].
x106
7
Basset
Mel et al
SLoth and Dorgan
6 . Equation 28
4
S3
2
00 01 02 03 04 05 06
t[s]
Figure 5: Numerical results of the history force for sev
eral kernels. pp 7850 [kg m3], D
1 103 [m].
600
Figure 6: Experimental obtained velocity profile, solid
black line, numerical predictions without his
tory force, dashed line, and with history force
using the Basset kernel by Mordant and Pin
ton (2000). Prediction using the kernel of
equation (28) is added for comparison, solid
blue line. pp 7850 [kg m3], D 1 1
103 [m].
locity, and slightly underpredicts the initial velocity, re
sults are quite satisfactory.
Another case from the experiments has been selected
to present a case for which a sphere attains a high
Reynolds number during the settling process. The high
est Reynolds number that was observed in the experi
ments was Re = 7700. This Reynolds number implies
that the drag coefficient is in the Newton regime, the
regime with a near constant drag coefficient. Figure 7
shows the result of the experiment. It can be noted that
the experimentally obtained velocity is lower for part of
the domain than the calculation both without and with
the history force, using the Basset kernel. Numerical
prediction using the kernel of equation (28) does not
show the initial overprediction of the velocity. Due to
the drag correction in the kernel, the history force can
obtain larger values than when using the Basset kernel.
However, the current kernel fails to correctly predict the
particle velocity near its terminal velocity and further in
vestigation is required.
One of the causes for different particle behaviour
could be the kinematic viscosity used in the numeri
cal calculations. The Reynolds number listed by Mor
dant and Pinton (2000) can be obtained using either a
kinematic viscosity of v, 9.02 10 7 [in2 s1] or
a diameter of D 5.6 103 [m], while the kine
matic viscosity used in numerical calculations was set
to v, 8.92 10 7 [m2s1] and the diameter to
D 6 103 [m]. Although such a difference in kine
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
0 200 400 600 800
time ( ms )
Figure 7: Experimental obtained velocity profile, solid
black line, numerical predictions without his
tory force, dashed line, and with history force
using the Basset kernel by Mordant and Pin
ton (2000). Prediction using the kernel of
equation (28) is added for comparison, solid
blue line. pp 7750 [kg m3], D 6 6
103[m].
matic viscosity or diameter seems minor, it has a sig
nificant impact on the obtained Reynolds number and
subsequently on the obtained drag factor, equation (10),
and history factor, equation (26). As the flow is in the
Newton regime, a near constant drag coefficient implies
that the drag factor increases almost linearly with the
Reynolds number.
The above two cases are two important cases, one
for a relatively low Reynolds number, and one for a
very large Reynolds number. The kernel introduced in
the present research is able to predict particle behaviour
more accurately than the Basset kernel and agrees with
the experimentally determined history force. In the re
mainder of this paper the obtained kernel, equation (28),
will be used in the prediction of particle behaviour in
rotating and swirling flow fields.
Swirling flow field
One of the properties of swirling flow, is the rotating
flow field. Therefore, a flow field consisting of purely
solid body rotation has been studied to further improve
the prediction of particle trajectories. In Candelier et al.
(2005) and Candelier (2008) such a case of solid body
rotation has been presented, and predictions of particle
behaviour using the present model will be compared to
results obtained by Candelier (2008).
A solid body rotation is described with a single pa
rameter w for the rotation. Using this parameter the flow
can be described by:
S(
0WY
0o
Important properties of such a flow field for particle
tracking include the material derivative, the rotation of
the flow field and the Laplacian of the flow field:
DuW 2 Y ;wc 0 ;V2U =0,
0 2w
(37)
respectively. The material derivative is used in both the
added mass force and the description of the pressure gra
dient and stress divergence forces. The rotation is re
quired for the lift forces. When the Laplacian of the
flow field is not zero, corrections have to be added to
the forces, known as the Faxen force corrections, see e.g.
Crowe et al. (1998). As the Laplacian of a solid body ro
tation is zero, such corrections are not required and will
therefore simplify the equations of motion to the form
presented in equation (20). It will also be important to
notice that u x w = 2.D
Candelier (2008) studied two situations for the case of
solid body rotation. The first case is a particle released
in the flow field with gravity working in the direction of
the axis of rotation. The second case, only numerically
studied, is a particle released in the flow field with grav
ity working in the direction perpendicular to the axis of
rotation. The first case shows that the rate of rotation has
an influence on the particle settling velocity. In equa
tion (20) all terms involving rotation only act in a plane
perpendicular to the axis of rotation, and therefore the
currently used equation fails to describe this particle be
haviour.
In the second case gravity is acting in the plane per
pendicular to the axis of rotation. There will be no move
ment in a direction parallel to the axis of rotation, and
therefore an influence of the rotation to the drag in the
direction parallel to the axis of rotation will not influence
results. If a particle is released at the axis of rotation in
such a flow field, the particle will be forced away from
the axis of rotation due to gravity. However, centrifu
gal forces, caused by the flow field, will force the parti
cle back to the axis of rotation. In the case studied, the
forces acting on the particle will cause the particle to at
tain equilibrium at a specific position with respect to the
axis of rotation. The location of that position will first be
predicted using the assumption of steady state. The par
ticle velocity and its derivatives will then be zero at this
equilibrium point. The steady state of the particle also
implies that rotation of the particle will be equal to half
the rotation rate of the flow field, because the torque on
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
the particle will than be zero. The reduced equation of
motion for the particle in terms of x, y, w, g, with gravity
working in negative ydirection will then become:
0(1
w2[A ( Y
XA 
B(x
)13
(38)
with the dimensionless factor A (x):
A (x) W T
A=TV
and dimensionless factor B:
3 9.66 v
B(x)= 1 +Ca ,M "2 9 L,S c. (40)
4 7D 2w
Equation (38) provides a relation for the equilibrium
location of the particle:
Depending on the sign of A, B and 1 p, the point
where steady state is attained is within a specific quad
rant with respect to the axis of rotation. A can never
be zero, and therefore the equilibrium point can only be
on the line x = 0 if the density ratio = 1. This im
plies that the equilibrium point can only be at the origin
if the equilibrium point is located at x 0. Candelier
(2008) used as parameters D = 3 103 [m], = 0.1,
v, = 2 104 [m2 s ] and w 15 [rads ]. Gravity
is chosen to be the same as in the numerical simulations
for a settling sphere. Using an iterative approach, the
equilibrium location for this situation is then obtained
as:
X x ( 0.45
XD D 0.067 ) (42)
However, the location of the equilibrium position shown
by Candelier (2008) has a positive value for the y coor
dinate. The sign of B determines the sign of the y coor
dinate, and if lift forces are taken into account in such a
way that B > 0 the ycoordinate will become positive.
In figure 8 the particle trajectory predicted by the
equations of motion in equation (20) is presented. The
trajectory converges to the equilibrium point predicted in
equation (42). In figure 9 the result of a simulation tak
ing the history force into account has been compared to
the result without taking the history force into account.
The history force decreases the maximum excursion of
the particle in x and ydirection. Therefore, the particle
trajectory without history force will initially be a curve
with a larger radius than the trajectory with history force.
After some time however, the trajectory will converge to
05 04
XD[
00
Figure 8: Prediction of a particle trajectory in a solid
body rotation flow field with gravity acting
in a plane perpendicular to the axis of rota
tion. D = 3 103 [m], 0.1, vc = 2
104 [m2 s1] and w 15 [rads ]. x and
y coordinates have been divided by the par
ticle diameter D to obtain nondimensional
variables.
08 07 06 05 04
XD[
00
Figure 9: Predictions of particle trajectories in the solid
body rotation flow field also used in figure
8 with and without taking the history force
into account. x and y coordinates have been
divided by the particle diameter D to obtain
nondimensional variables.
_ ( A (x)) 1
B (x) W 2 (A2 (x) + B2 (X)) 9.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Table 1: Summary of fluid, flow field and particle prop
erties used in computations for swirling flow.
Variable Value Units
U 2 ms 1
F 2 m2s 1
TV 500 s
Pc 997 kgm 3
vc 8.9 107 m2s 1
pp 800 kg m3
D 1 i 104 m
Ro 5 102 [m]
the equilibrium position faster than the trajectory pre
dicted taking history force into account.
Forces in the direction parallel to the axis of rotation
caused by rotation still need to be investigated. How
ever, present research first focuses on the effect of the
history force for particles in a swirling flow field, i.e. a
flow field with both rotation and an axial velocity. To in
vestigate the influence of the history force for such flow
fields, a generic swirling flow field is described as:
U (
(1
(1
2rrr2 i
2 )
2 )
where F is a parameter determining the circulation of
the vortical flow with dimension [m2s1], and T, is a
parameter determining the decay in time of the circula
tion away from the axis of rotation with dimension [s].
The rotating part of the velocity is a LambOseen vortex,
however to provide a steady flow field the decay in time
has been frozen by taking T, constant. The parameters
for the performed simulations are summarized in table
1.
The parameters for the flow field and the particle have
been chosen to represent conditions similar to swirling
flow produced water treatment. During such a separation
process oil droplets up to 1 104 [m] are found. As
a reduction in diameter significantly increases the time
for an oil droplet to reach the axis of rotation, a droplet
size has been chosen close to the upper limit of often
occurring droplet diameters.
Figures 10 and 11 show the resulting trajectory of a
droplet, predicted taking the history force into account
as well as the trajectory of a droplet without taking the
history force into account. The droplet is released at
x Ro with the velocity of the flow field at that lo
cation. Distances have been divided by the initial radius
with respect to the axis of rotation to obtain dimension
less distances. The trajectories follow the rotating flow
25  I
S With history
I Without history
20
20   
00 01 02 03 04 05 06 07 08 09 10
R1
Figure 10: Dimensionless distance travelled in axial di
rection versus dimensionless radial position
of the particle trajectory for a swirling flow
field. Radius and axial distance have been
divided by the initial radius Ro to obtain
nondimensional variables. Table 1 summa
rizes the flow conditions. Density ratio of
particle and fluid p 0.8
1 00 L
1 00
NxlO
080 060 0 40 0 20 0 00 020 040 060 080 1 00
X1
Figure 11: Particle trajectories projected on the plane
perpendicular to the axis of rotation. x and y
position have been divided by the initial ra
dius Ro to obtain nondimensional variables.
Table 1 summarizes the flow conditions
and gradually move to the axis of rotation. Initially the
effect of the history force is not very large. However,
the distance required for the particle to move within one
tenth of the initial radius shows a difference of about
two times the initial radius. This case shows that a sig
nificant influence of the history force for swirling flow
separation is to be expected.
Conclusions
To provide an adequate description of the trajectories of
particles in separation processes a Lagrangian particle
tracking model is used. First the settling of particles in
quiescent fluids is investigated, numerical predictions
are compared with available experimental results. The
results showed the importance of a history force acting
on the particle. A square root decay, resulting from the
Basset kernel, is not able to describe the quadratic decay
found for the long time behaviour of the history force.
Therefore, a kernel initially proposed by Lawrence and
Mei (1995) is adapted to provide a kernel for the history
term with a smooth transition from square root decay to
a quadratic decay. The kernel used in the present re
search depends among others on the drag factor, and is
therefore able to adequately predict particle behaviour
for a large range of Reynolds numbers However, in
some cases predictions of the particle reaching terminal
velocity do not agree completely with experimental data.
The influence of the history term is subsequently in
vestigated in a solid body rotation flow field. A numeri
cal study performed by Candelier (2008) is used to pro
vide for a case of a solid rotation flow field. For the case
of solid body rotation flow with gravity acting parallel
to the axis of rotation a significant influence of rotation
rate on the particle trajectory can be found. This cannot
be predicted by equations of motion, used in present re
search, that include only Saffman and Magnus lift forces
as rotation dependant forces. However, the observed ef
fect will be relatively small for swirling flow, as axial
velocity differences will generally be small.
The second case studied is the case of solid body ro
tation flow with gravity acting perpendicular to the axis
of rotation. Simplification of the equations of motion for
a steady state situation provides a prediction of the point
of equilibrium for a steady state situation. Both the par
ticle trajectory taking history force into account and the
particle trajectory without taking history force into ac
count were found to converge to the predicted equilib
rium point. However, history force significantly influ
ences the maximum excursion of the particle as well as
the time in which the particle converges.
Finally the case of swirling flow is investigated. The
history force shows a significant influence under condi
tions occurring during produced water treatment. Dif
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
ferences in the order of pipe diameter can be observed
for the distance it takes to reach a radius will outside the
axis of rotation.
Acknowledgements
The research presented is made possible by FMC Tech
nologies/CDS Separation Systems.
References
T. R. Auton, J. C. R. Hunt, and M. Prud'Homme.
The force exerted on a body in inviscid unsteady non
uniform rotational flow. Journal of Fluid Mechanics,
197:241257, December 1988.
Fabian A. Bombardelli, Andrea E. Gonzalez, and
Yarko I. Nino. Computation of the particle basset force
with a fractionalderivative approach. Journal of Hy
draulic Engineering, 134(10):15131520, 2008.
F. Candelier. Timedependent force acting on a particle
moving arbitrarily in a rotating flow, at small Reynolds
and Taylor numbers. Journal of Fluid Mechanics, 608:
319336, August 2008.
F Candelier, J.R. Angilella, and M. Souhar. On the ef
fect of inertia and history forces on the slow motion of a
spherical solid or gaseous inclusion in a solidbody ro
tation flow. Journal of Fluid Mechanics, 545:113139,
December 2005.
R. Clift and W.H. Gauvin. The motion of particles in
turbulent gas streams. Proceedings of Chemeca, 1:14
28, August 1970.
C. F M. Coimbra and R. H. Rangel. General solution
of the particle momentum equation in unsteady Stokes
flows. Journal ofFluidMechanics, 370:5372, Septem
ber 1998.
Clayton T Crowe, Martin Sommerfeld, and Yutaka
Tsuji. Multiphase flows with droplets and particles.
CRC Press, 6000 Broken Sound Parkway, NW, (Suite
300) Boca Raton, Florida 33487, USA, 1998. URL
www.crcpress.com.
David S. Dandy and Harry A. Dwyer. A sphere in shear
flow at finite reynolds number: effect of shear on particle
lift, drag, and heat transfer. Journal ofFluidMechanics,
216:381410, July 1990.
Christopher J. Lawrence and Renwei Mei. Longtime
behaviour of the drag on a body in impulsive mo
tion. Journal ofFluidMechanics, 283:307327, January
1995.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
E. Loth and A. J. Dorgan. An equation of motion for
particles of finite Reynolds number and size. Environ
mental FluidMechanics, 9(2):187206, April 2009.
Renwei Mei. An approximate expression for the shear
lift force on a spherical particle at finite Reynolds num
ber. International Journal of Multiphase Flow, 18(1):
145147, January 1992.
Renwei Mei and Ronald J. Adrian. Flow past a sphere
with an oscillation in the freestream velocity and un
stready drag at finite Reynolds number. Journal ofFluid
Mechanics, 237:323341, April 1992.
Renwei Mei, Christopher J. Lawrence, and Ronald J.
Adrian. Unsteady drag on a sphere at finite Reynolds
number with small fluctuations in the freestream veloc
ity. Journal ofFluid Mechanics, 233:613631, Decem
ber 1991.
N. Mordant and J.F. Pinton. Velocity measurement of
a settling sphere. European Physical Journal B, 18(2):
343352, November 2000.
S. I. Rubinow and Joseph B. Keller. The transverse force
on a spinning sphere moving in a viscous fluid. Journal
ofFluid Mechanics, 11(03):447459, November 1961.
P. G. Saffman. The lift on a small sphere in a slow shear
flow. Journal ofFluid Mechanics, 22(2):385400, June
1965.
L. Schiller and A. Naumann. Uber die grundlegenden
berechungen bei der schwerkraftaufbereitung. Vereines
DeutscherIngenieure, 77:318320, 1933.
