7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Lagrangian Concentration Differential Equation for Multiphase Flows
Sida Wang and Eric Loth
Department of Aerospace Engineering, University of Illinois UrbanaChampaign, 104 South Wright Street, Urbana, IL 61801,
USA
Department of Mechanical & Aerospace Engineering, University of Virginia, 122 Engineer's Way,
Charlottesville, VA 229044746
swang46@illinois.edu and loth@virginia.edu
Keywords: Lagrangian, Eulerian, particle tracking, concentration, impact efficiency
Abstract
A new technique is proposed for computing Eulerian particle concentration and fluxes with Lagrangian trajectories. In
particular, a Lagrangian concentration differential equation is solved along a particle path using Eulerian derivatives for the
particle velocity divergence field. This is referred to as the Lagrangian Concentration Differential Equation (LCDE)
technique and is shown to provide efficient and accurate results for steady flows. It is compared to other Lagrangian
techniques using higherorder trajectory simulations for particle with a linear drag law in simple flows. In particular, steady
twodimensional flows in a comer and past a cylinder are considered. The flux fields are predicted for tracer particles and
for particle of various inertias, i.e. Stokes numbers. The LCDE approach combines the accuracy of areabased Lagrangian
methods (which are limited to steady flows) and the unsteady capabilities of binbased Lagrangian methods. In future work,
this approach may be investigated with respect to unsteady and threedimensional flows.
Introduction
A natural approach when simulating individual
particle interactions (with the fluid, with another particle,
or with a wall) is to use particle trajectories, i.e. the
Lagrangian dispersedphase representation. In this case,
particle pathlines are defined based on the center of mass
of the particle (or the center of mass for a cloud of
particles) along which Ordinary Differential Equations are
used to update particle position, velocity, mass,
temperature, etc.. This provides high accuracy with
respect along a particle path, including particlewall and
particleparticle reflections. However, determination of
volume or areaaveraged quantities, such as particle
concentration or mass flux per unit area, on an Eulerian
grid requires summation or interpolation from the
Lagrangian trajectories. This can lead to uncertainty
and/or inaccuracy due to spatial averaging, especially for
unsteady flows. Numerical errors associated with
concentration predictions can lead to problems for
Lagrangian particles when twoway coupling is important
(Sivier et al. 1996).
To avoid this problem, Eulerian dispersedphase
representations are often used when particle concentrations
are important, especially when twoway coupling
phenomenons are important. Such an approach employs
Partial Differential equations (PDEs) for the particle
characteristics (e.g. velocity and concentration) which are
discretized as cellaverages on an Eulerian grid. These
fields are integrated in time based on discrete spatial and
temporal gradients and assumptions of continuity.
However, this can lead to nonphysical diffusion of particle
mass and momentum during advection. Furthermore,
reflecting boundary conditions are difficult to treat within
an Eulerian approach.
As such, it is desirable to develop a method which can
incorporate the trajectory accuracy of the Lagrangian
approach and the concentration accuracy of the Eulerian
approach. A recent hybrid concept termed the Full
Lagrangian method was developed by Healey & Young
(2005) based on a linear drag force and a Jacobian system
of equations to be solved. The present approach seeks to
use elements of this approach while preserving the
simplicity of conventional Lagrangian approaches, e.g.
allowing for generalized surface forces and trajectory
ODE's. The latter is accomplished by using Eulerian
spatial derivatives to estimate the particle velocity
divergence. The present technique is evaluated by
comparison to conventional techniques for two
fundamental flowfields.
Nomenclature
Area
particlepath derivative
Stokes drag correction factor
gravity acceleration vector
m particle mass
in mass flux
M mass flux per unit time
N number of particles
n outward normal vector
Re Reynolds number
St Stokes number
t time
u fluid velocity
v particle velocity
w particle relative velocity
x horizontal particle location in Cartesian
y vertical particle location in Cartesian
F Force
Greek letters
a particle concentration
+ velocity potential
S discretization shape function
11 vertical direction in natural coordinates
1' dynamic viscosity
P density
T response time
S volume
horizontal direction in natural coordinates
Subscripts
x horizontal direction in Cartesian system
y vertical direction in Cartesian system
turb turbulent
p particle
f fluid
D Drag
Numerical Schemes
The overall particle translational equation of motion
equates the rate of change of the particle's linear
momentum to the net sum of the forces acting on the
particle of a given mass. This mass (mp) can be
expressed in terms of the particle volume (Vp) and density
(p,), or in terms of the particle diameter (d):
mp = ppp 6
A = ~6
If one assumes negligible contact interactions (i.e. particle
concentrations are small enough such that particleparticle
collisions are not significant), the Lagrangian equation for
particle velocity (v) is given by
mp (dr / d) = mpg +Fsurf 2
The left hand side (LHS) includes the particle path
temporal derivative of the particle velocity while the right
hand side (RHS) includes the body force based on gravity
(g) and fluid dynamic surface forces (Fsf). Similarly, an
ordinary differential equation (ODE) can be written for
particle position (Xp) and mass as
dxp /df =v 3a
dinp / f = rp 3b
As in equation (Eq.) 2, the LHS of equations (Eqs.) 3a and
3b include particlepath temporal derivatives. The RHS
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
of Eq. 3b is the mass transfer from the surrounding fluid to
the particle, which is a function of the difference between
vapor species fractions on the particle surface and that in
the farfield in the case of droplet evaporation.
The RHS of Eqs. 2 and 3b include F,,f and
mp which are the surface integrated fluid dynamic force
and mass flux to the particle. One may simulate the
detailed flow and species around the particle to obtain
these, but this is typically computationally prohibitive if
there are many small particles. The pointforce approach
avoids this problem by employing theoretical and
empirical models for Fs. and p based on relative
differences between particle characteristics and that of the
surrounding fluid.
The dominant contribution of the surface force is
typically the drag component which is nonzero in the
presence of finite viscosity of the surrounding fluid (Lf).
For a spherical solid particle surrounded by a uniform
continuum flow of with constant density (pf), this force can
be represented by the Stokes drag and a Stokes correction
(f) as:
FD = 37tdfgtfW 4a
Rep = 4b
Lf
W (t) v(t) up (t) 4c
The particle relative velocity (w) is based on the
"unhindered" continuousfluid velocity (u@p),
extrapolated to the particle centroid (x,). The
"unhindered" velocity neglects local flow disturbances
caused by the particle itself, i.e. it is the velocity that
would occur if the particle was not present.
The Stokes correction factor is unity for a particle
Reynolds number which is small (Rep< 1) and otherwise
can be modeled with the SchillerNaumann (1935)
empirical expression for finite Reynolds numbers:
f = 1+0.15Re 687 for Rep < 1000 5
If the particle density is large compared to that of the
surrounding fluid (pppf, e.g. water droplet in air), the
virtual mass is negligible and the effective mass is equal to
the particle mass. In this case, temporal response to the
restoring drag force can be used to define a particle
response time as
Smjwj pd2 6
p 6
S FD 18 lff
In addition, the high density ratio limit indicates that the
other surface force components are generally small
compared to the drag as long as particle spin and fluid
shear effects are weak (Maxey & Riley, 1983) so that
Fs F 7
Herein, this simple surface force will be employed with
Eqs. 4a and 5, but the proposed particle concentration
approach is general enough to include any surface force
model, i.e. can be extended to include lift, history force,
fluid stress and added mass force terms using the
techniques discussed by Loth & Dorgan (2009).
Lagrangian Methods for Particle Trajectories
The momentum and position ODE's (Eqs. 2 and 3a)
may be integrated in time for a discrete time increment (At)
assuming constant zp and u@p method proposed by Barton
(1996) which assumes quadratic variation of the fluid
velocity in time:
n+l ne At/zp + U l e At/,
Vp = Vpe +u@p e
+B{At2 2TAt+2yfl2r eAt/z
+ ( Atn I At/,(l,. +u~pt
8a
xn+l n n n nt +un At,8b
p p p @p P p 8b
3At2 ,n nz
+ABJ rAt.ixill e
+B{ At2 + 2 )2 At 2 (,)3 (1 e At
Where the vectors A and B are represented as
un+ n1 (/i n
A @p @p 9a
2At a )
un+ 2un +un1 2U n
2At2 2 t29
The above method will be referred as Scheme 3 and is
secondorder in time with respect to acceleration
(firstorder in time with respect to velocity). Barton
(1996) also introduced two lowerorder schemes that can
be obtained from Eqs. 8 and 9. Scheme 1 sets the vectors
A and B to zero and uses a predictorcorrect step to obtain
the fluid velocity at the n+l time level and average with
the fluid velocity at the present time. Scheme 2 sets the
vector B to zero and averages the flow velocity at the
present and n+l time level over the change in time for A.
The accuracy benefit of Scheme 3 lies in the ability to
incorporate changes in fluid velocity along the particle
path (given by B).
Ensemble Average Method
Perhaps the simplest and most common approach to
determine Eulerian concentration from Lagrangian particle
trajectories is the ensembleaverage bin method. At a
given instant of time, a summation over a computational
cell volume (Fig. 1) can be used to express concentration
as:
SN,,
a = V, 10
VA,1 jl
In this expression, VA,1 is the volume associated with
node x, which contains NpA number of particles which
volumes and interphase forces are given by Vpj. This
can be extended to parcels which have Npp particles per
parcel and NpA parcels per cell as:
a1 _1 Np 11
a,1 j=i
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
The parcel volume is assumed to be less than the
computational cell volume.
Figure 1: Schematic of a twodimensional Eulerian
continuousphase grid showing summation of particle
volumes in a computational volume to compute volume
fraction associated with a node x,.
Parcels are more computationally efficient than
tracking individual particles if there are many particles in
the domain. However, both Eqs. 10 and 11 are subjected
to statistical averaging if the control volume is independent
of the particle spacing. In particular, the statistical
uncertainty of the particle concentration (Lbt ) for a
random correlation with individual particles is inversely
proportional to the square root of the sampling frequency:
iLx 1
12
Based on this, 100 particles per cell translates into an
uncertainty of 10%, while 104 particles would be needed to
ensure an uncertainty of 1%. Thus, a large number of
particles per cell is needed for accurate predictions of
Eulerian concentration gradients, i.e.:
NpA > 1 for accurate cell concentration 13
As a result, the number of particles required for
cellaveraged concentration convergence greatly exceeds
the number of fluid nodes in a domain (DeAngelis et al.
1997; Loth, 2000).
To obtain particle mass flux per unit time (M ) for a
given discrete area (AA,1) with the bin approach, a flux
summation through a discrete computational area over a
time average can be used:
N,,
I Np 14
p,_1 mp,3 14
avg ]=1
In this expression, Tavg is the time used for temporal
averaging and NpA, is total number of trajectories crossing
the discrete area in this averaging time. Often the
discrete mass flux per area is of interest, which may then
be obtained for bin i as
SA,1 AI t avg j
In either case, statistical uncertainty requires that the
number of particles (or parcels) used in the summation be
large per Eq. 15. For a 2D flow with a vertical surface,
the incremental area along the yaxis is based on the
location of the bin centroid (y,) and the neighboring bin
centroids:
Ay, =(y, 1+y,) (y, +y 1) 16
Thus, particles are only counted in a bin if
(y, +yi)
For uniform spacing, this simplifies to a constant value for
all bins. A representation of a 1D bin is shown in Fig. 2a.
1.2I
1

V V V
1.5 1 0.5 0 0.5
(yy.)/Ay
Figure 2: Shape function along a
Ensemble Averaging
1 1.5
line segment for
The timeaveraged particle concentration by transport
theorem as can be obtained as
1Kp = (app)(v.n) dA 18
cs
Defining n as the normal to the surface area, the dot
product effectively modifies the area integral to be over an
area perpendicular to the particle velocity. The mass flux
for onedimensional flow is thus conserved if there is no
mass transfer between the particle and fluid phases (as is
assumed herein). For the particle injection flux, this
yields:
Mlnj = anPpVln An 19
Using Eqs. 14 and 19, the mean volume fraction at an
Eulerian flux node (x,) can be expressed as a flux sum of
the number of particles through a discrete area as:
S(app)(vn) dAA
1 NpA
1 "^
Ym jl
avg Z=1
Assuming constant density for each particle the integral
can be expressed as:
(pp) AA v.n = m 21
Tavg j=1
Averaging the velocity for each particle, the mean
concentration becomes:
NpA,,
jpal 1
a = F)PAA^ir N 22
P p A A 1T a g V.
J1
In the above expression AA,1, is the surface area per bin
number (i). A similar summation can be used for NPA
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
parcels for a parcel approach. With the knowledge of the
concentration, the change in concentration for the
ensemble binaverage method can be expressed as:
aC 23
a,nj
For a steady flow with Np,, particles, the particles are
injected with constant mass over an area AmJ at uniform
v n,1. Hence, this ratio simplifies to:
AA,lnjNpA, ( n) 24
a = 24
NpA,I
AA,1 (V n,)
j=1
If the particles are injected over a vertical plane whose
normal is in the xdirection, then vj ni = vx,j and if the
particles are tracers which follow the fluid streamlines (i.e.
v=u for St=0), then the concentration density is constant so
that a*=l. Again, the number of particles (or parcels)
used in the summation should be large for accurate results.
Weighted Average Method
A variant of the above ensemble bin method is to use a
shape function to weight the contribution from each
particle (or each parcel). The shape functions may be
quadratic or higherorder, but linear variations are typically
used since this is generally consistent with the treatment of
the fluid properties. For a 2D flow, the discrete flux area
per unit depth is simply a line segment. If one considers
a vertical surface discretized into segments of length Ay
with junctions y,, the linear shape function over two
neighboring segments can be described as:
l ypyJ/Ay if ypYl Ay
ti,= 25
0 if ypy, > Ay
A 2D linear shape function with particles impacting the
collection plane is shown in Fig. 3. For the
weightedaverage bin approach, the mass flux through an
incremental area (between two neighboring junctions)
becomes
SNpA4,
, Z (() jmpj) 26
Tavg j1
For a vertical surface, the incremental area along the yaxis
is centered at the node junction and can be expressed using
Eq. 17. Note that the junction values are not the same as
the centroid values used in the ensembleaverage bin
method. For the discrete mass flux per area
NpAI
Zp = 1 ((m)l 27
SAJA1 AA'Tlavg j=1
Note that the incremental areas overlap between adjoining
nodes so that the mass flux for a particle between two
nodes is distributed to the associated discrete area for each
node. The mean volume fraction at an Eulerian junction
node can be expressed as a flux sum through a discrete
area over the number of particles:
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
NpA,I
Z (mDl,,mp1)
N NA
pa^,l 28
Pa A Tavg NZ (I
J1
If we consider a steady flow with Np,, particles of constant
particle mass injected over an area the concentration ratio
for the weighted method is:
NpA ,
Aln, (v n)l NpAI (
a _ j 29
lnj AA,1 NpA'I
Npnj YZ ((),jVj "n1
J=1
To apply these methods, particles must be injected at an
upstream inflow plane and then collected at the flux plane
of interest. If the flow is also steady and the inflow
particle velocity is uniform, the mass flux may be
computed by releasing a single set of particles at the inflow
plane with this spacing and determining the mass fluxes
over discrete downstream areas based on avg=Atm,.
1.2
1
0.5
 .5
1 0.5
y
0
(yyi)/Ay
0.5
0.5
Figure 3: Shape function along a
Weighted Averaging
1 1.5
line segment for
Areabased Method
The areabased technique employs the steady uniform
inflow assumptions needed for Eq. 19, to relate the change
in the streamtube area to the change in mass flux
perpendicular to that area. This technique has been
developed by NASA (Ruff and Berkowitz, 1990) and is
efficient for simple geometries, e.g. where trajectories do
not cross. It notes that the discrete mass flux through a
particle streamtube is constant:
Mpnj = Mpflux
Since the streamtube is defined by a set of particle paths,
the 2D streamtube and trajectories and the associated
areas can be represent as in Fig. 4. For 2D, the area
between trajectories maybe calculated as:
A, (X x)2 +(y y)2 31
From this, the mass flux per unit area can be expressed for
2D flows as:
i M
AA ,
(mp1 +mp,)
2AA, Tavg
Note that 3D flows require three or four paths. The
volume fraction for a streamtube between trajectories can
be expressed as a flux through the discrete area:
NpA,,
SiZ(mp,)
1 =l _
l1
PPAA'lZavg I (Vi .ni)
Jl
Note that this is the same as the ensemble average equation.
However, the area method has a variable discrete area
which depends on the a preselected NpA,l (=2 in 2D).
Because it captures the stream tube exactly, the areabased
methods and is not subjected to statistical uncertainty
associated with the binbased the ensemble average
methods. The concentration ratio for the area ratio
method is given by:
Anj NpA,1(v )
a "n 34
NpA,
A ,Z( ,n,)
j1
Since the number of particles required for the streamtube
for each discrete region is much less than the ones for the
ensemble method and its variant, the computational cost
for the areabased method is substantially lower.
rtidletrajectories 
Ai / v.n
Figure 4. Stream tube for the areabased technique
Concentration Differential Equation
The present method is based on "full Lagrangian"
method of Healy & Young (2005) but avoids the numerical
integration of a Jacobian system and instead samples the
local particle velocity gradient field for a streamtube
(based on two particle trajectories that are close to each
other). It uses the Eulerian particle concentration PDE
(assuming constant particle density):
'8a a (avj) ai ip
+_ 35
5t Ox Vppp
This can be formulated as:
Sln(ca)
dt
mli
V.V+
The differential equation above is integrated along the
particle path along with the particle position and velocity.
Assuming no added mass, the relationship can be
expressed in exponential finite difference form for a
discrete time increment (At) to yield:
an+l = a exp(AtV v)
For each n+1 timestep, the previous time concentration is
needed to for the computation. Hence, the concentration
is computed along with particle trajectory and velocity
when marching the particle trajectory forward in time
using the third order accurate exponential methods in by
Barton (1996).
The key to the present method is to utilize the two
particle trajectories to compute the velocity gradient. It
also avoids the assumption that the drag force is linear.
Hence, the concentration ratio at different locations is the
same as Eq. 23.
Discretization of the Particle Velocity Divergence
Particle velocity divergence is of pinnacle importance
for the numerical treatment of the concentration
differential equation. A streamtube approach is used for
its computation with isoparametric formulations that
express the element position and element velocity in the
form of interpolations using the natural coordinate system.
For the purpose of validation only 2D quadrilateral
elements are being considered.
Four points for the covers of a quadrilateral element
are computed using two trajectories for both velocity and
position. These values are then being used to compute the
velocity divergence in the middle of the element using the
following formulations.
The general twodimensional element coordinate
interpolations are:
4
x = ,x,1 38
i1i
4
y = ,y, 38b
Where x and y are the coordinates at any point of the
element, and x1, y,, i = 1,...,4, are the coordinates of the 4
element nodes. (1 is the interpolation function in the
natural coordinate system, which has variables 4 and fr that
each vary from +1 to 1. These variables are the direction
coordinates in natural coordinates that are equivalent to x
and y in the Cartesian coordinate system. The same
interpolation can be applied to element velocities. To
evaluate the displacement derivatives, the following can be
used:
9 39
The first term on the RHS is the Jacobian. Then from Eq.
39 inverting the Jacobian will results in
? 1 
S40
Apply Eq. 40 to the element displacement:
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
~v~ v av,, NvY
c9x c9x 9xV
9v 9v 2v, 41v
_y 9r y 9ly
Using Eq. 41 the divergence between the particle
streamtube can be obtained and by the following equation:
V.v + 42
9x 0y
The above treatment of divergence can be simplified for
uniform Cartesian grids on simple flowfields by finite
difference. The Lagrangian Concentration Different
Equation (LCDE) method uses the same number of
particles as the areabased method for a 2D fluid.
Although the LCDE method is more computationally
expensive than the areabased method, it handles
concentration effectively and is applicable to complex
flows.
Flow Fields and Test Conditions
Corer Flow
An inviscid comer flow can be described by the
velocity potential:
A(y2 x2) 43
=^f  43
2
In this expression, A is a positive constant. The velocity
components for the fluid flowfield are then:
u Ax 44
Ox
u =c9=Ay 44b
ay
For computational simplicity the constant A is set to be 1.
The stagnationpoint is at the origin.
The particle trajectories are obtained analytically by
solving the following differential equations of position
with the initial conditions x(0) = 1, vx(0)= 1, y(0) =yo and
Vy(0)=yo, where yo is the initial y position of the particle
released:
92x Ox
S+ +2x = 0 45
C 2zy = 0 45b
at 2 t
The solution to Eq. 45 is then used as the analytical
trajectory in later studies.
Cylinder Flow
An inviscid, incompressible twodimensional flow
over a cylinder is being using for the simple cylinder flow
case. This case is unphysical because there are no
boundary layers, but it is still a useful test because the
existence of an analytical solution for the velocity field.
The fluid flowfield is defined by the velocity potential:
4f = Ux x 1+ 46
In this expression, U is the freestream velocity of the flow
In this expression, U. is the freestream velocity of the flow
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
field, and R is the radius of the cylinder. The velocity
components for the fluid flowfield are consequently as:
1 (y2x2R2
u, U, 1+ 47
S(x2 +y)2 )2
8df 2U.xyR2 47b
U y (x2 +2)2
For computational simplicity the constant U~ and R are
both set to 1. They cylinder diameter is D=2R.
Stokes Number and Trajectory Accuracy
The Stokes number (St) is also called the Modified
Inertia Parameters and is the ratio of the particle response
time (T ) and the continuousfluid velocity changes (rf)
St 2 48
Tf 18TLfrf
This expression assumes linear drag and the fluid response
time to be 1 for all of the studies.
Results and Discussion
Particle Equation of Motion Accuracy
The analytical trajectories of the comer and cylinder flows
can be used to evaluate the accuracy of the particle
trajectories for different numerical methods, and the max
error along a particle trajectory (e) in the horizontal of the
velocity is defined as:
e = max a v 49
Vx,a
for 1=1,max
Where x,a is the analytical particle x velocity, vx,n is the
numerical particle horizontal velocity from the different
numerical schemes, and i is the time step index.
The accuracy of the trajectories was first tested using
the analytical solutions at variable time steps ranging from
101 to 1025 which is nondimensionalized by the comer
length and velocity. With the collection plane at x = 0
and releasing from x = 1, a number of particles were
released for Stokes number of 0.1 using three exponential
methods from Barton (1996) to compute their trajectories.
The results can be seen in Fig. 5, which shows that the
velocity error decreases linearly with decreasing timestep
which is consistent with the accuracy expected from these
methods. The difference of values of error were
negligible at low time steps for all three methods.
e scheme 1
scheme 2
0 scheme 3
10_4
103 102 / 101
time step
Figure 5. Prediction of xvelocity errors for a comer flow
with St = 0.1
The cylinder flow is also tested where particles are
released from x/D = 1.5 and then collected at x/D = 1.5
for the zeroinertia case where the particle velocity is equal
to the flow velocity. An analytical solution can be obtained
from the zeroinertia. The maximum error is plotted
against different time steps in Fig. 6. As the time step
decreases, the error for Scheme 3 decreases consistent with
firstorder accuracy for particle velocity even when the
fluid velocity along the particle path is nonlinear. In
contrast, Schemes 1 and 2 do not show improvement in
error as the timestep decreases. Hence, Scheme 3 was
used for all subsequent studies.
101
102
scheme 1
scheme 2
scheme 3
104
103
102
time step
Figure 6. Prediction of xvelocity errors for trajectories in
a cylinder flow with St = 0.
Particle Concentration in a Comer Flowfield
The LCDE, ensemble bin, weighted averaging, and
areabased methods are investigated to evaluate the
concentration ratio (a*) for a simple steady
twodimensional comer flow in the second quadrant where
the fluid velocity for the x and y direction are described in
Eq. 44. The particles are released from x=2, and y= [0
1]. The comer flow is evaluated first because it consists
of simple fluid velocity fields and would give smooth
variation for the velocity and position that induces finite
gradient. Fig. 7 shows particle trajectories in a comer
flowfield for different Stokes numbers. In each case, the
particle initial velocities are set equal to that to of the
surrounding fluid. For St=0, the particle velocities are
simply set to the flow velocities so the trajectories would
follow the fluid streamlines (Fig. 7a) where it can be seen
that no particles hit the vertical wall. For St=0.5 (Fig. 7b),
the particle trajectories indicate less turning due to their
inertia. For St=10, (Fig. 7c) there is little turning and the
particles continue along at their injected velocity angles
and impact the wall accordingly.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
weighted methods The concentration ratio (a*) results
are plotted against the final vertical location (y/ymax@m) of
the particles or bins in Figs. 810.
From Fig. 8 shows that increasing the number of
particles hitting each bin results in a* converging to the
analytical result. For the case with only 5 particles
(distributed over 20 bins), a* has substantial errors since
the particle number per bin is low and some bins do not get
any contribution from the particles. When the particle
number released is increased to 50 particles the errors in
a* have decreased and are within 20% of the analytical
result. When the particle number reaches 500, a* is close
to the analytical solution. This is consistent with Eq. 13,
i.e. the bin based method requires a large number of
particles to provide consistently accurate result for a*.
2 Np = 5
 Np=5
Np = 50
1Np = 500
a*
Analytical
5
5
2 1.5
1 0.5
1 0.5
1 1.5 2 2.5
Y/YmaxOinj
Figure 8. Effect of Np= 5, 50, 500 on Concentration ratio
for St=0 for the ensemble method (different symbols for
5,50 and 500) with analytical results.
Similar results as the ensemble method can be seen for
the weighted averaging method in Fig. 9. As the number of
particles per weighted bin increase the concentration ratio
distributions became smoother and converge to that of the
numerical analytical solution. For the 5 particle case, a*
behavior similar to that for the ensemble method As Np
goes to 50, a* has improved greatly and is within 10% of
the analytical result. The effect of shape function
inclusion allows a smoother and more diffuse profile. For
the 500 particles case, the profile is quite similar to that for
the ensemble bin method, consistent again with Eq. 13.
2,n
(c)
Figure 7. Particle trajectories for Np=10 at a) St=0, b) St=
0.5, and c) St = 10
The ensemble average, weightedfunction average,
and LCDE methods are selected to see the effect of particle
resolutions on their concentration ratio. The particles are
released from the x=2 and then collected at x=1 (marked
in Fig. 7a) with the same Stokes number (St = 0) for
zeroinertia. The same release range of y = 0 to y = 1 and
bin size of O.lymax@mi are used for both ensemble and
a*
0.5
0 0.5 1 1.5 2 2.5
Y/Ymax@inj
Figure 9. Effect of Np= 5, 50, 500 on concentration ratio
for St=0 for weighted method with analytical results.
i ;: I: :' \I
S~S~
In Fig 10, the LCDE method is being tested for the
same conditions as the previous two numerical methods. It
can be observed that the number of particles gives no effect
on the results for different particle resolutions. The results
for a* are very similar for 5, 50, and 500 particles
collected. This is consistent with the formulation of the
LCDE method based on a* does not rely statistical
convergence and only depends on the trajectory accuracy
of the particles. Since the areabased method also uses a
two particle trajectories, it provides the same level of
accuracy.
2
18
1 6
1 4
1 2
%* 1.H LL:L VW '
0 05 1 15 2 25
Y/Ymax@inj
Figure 10. Effect of Np= 5, 50, 500 on Concentration ratio
for St=0 for LCDE method with analytical results
Based on the above particle resolution study, Np=50
gives a reasonable result for the bin based methods and is
therefore used for the remaining studies. For the
following cases, the particles are injected at x=2, and
collected at x=1 with a bin size of 0. lymax@m and with the
same injection range from y = 0 to y = 1. The
concentration ratio results are plotted against the final
vertical location of the particles (bins for the ensemble and
weighted method). Concentration ratios are plotted with
the analytical solution for St=0 (no change in concentration
consistent with incompressibility for the fluid flow) in
Figure 11. It can be observed that all four methods have
overall good agreement with the analytical result. However,
the bin based methods have significant fluctuations
compared to the other methods. This is particularly
noticeable for the ensembleaverage method because few
particles hit per bin. For the weightedaverage method,
the shape functions distribute the contributions resulting in
less errors. There are errors for the first and last bin of
the weightedaverage method and ensemble bin because
the bins are not getting consistent contributions. For the
areabased and LCDE method, the values of the
concentration are almost identical as the analytical
solution.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
1.5
 weighted
0 analytical
0 0.5 1 1.5 2
y/y
Ymax@inj
Figure 11. Normalized concentration at vertical wall of
comer flow forAx=0.1, N,=50, a) St=0
The numerical methods are further investigates for St
= 1 in Fig. 12 with otherwise the same conditions used for
Fig. 11. Similar trends for all four methods were
observed as the zeroinertia case. The ensembleaverage
bin method has the largest errors because of variations in
the number of particles hitting each bin. In this particular
case, the weighted averaging method gives relatively small
errors compared to the ensemble method by its use of the
linear shape function. The areabased and LCDE
methods give a constant distribution at the collection plane
and provides good agreement with the weighted averaging
method.
1.51
0 0.5 1
y/y .
Y/Ymax@inj
Figure 12. Normalized concentration
comer flow forAx=0.1, Np=50, St = 1
1.5 2
at vertical wall of
The next case considers St = 10 where the particle
trajectories generally follow the initial trajectory angles
and the concentration ratios are shown in Fig. 13. The
ensembleaverage bin method gives various fluctuations
because of the low number of particles hitting per bin. The
weighted averaging method reduces the fluctuation by the
use of the shape functions but have spikes at the initial and
final bins. The areabased and LCDE methods give
constant concentration ratios that agree well with the
weighted method.
1.5
1
0.5
0
0
1 1.5 2
y/ymaxinj
maxfinj
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
The concentration ratios are plotted against the
vertical collection plane location in Fig. 15 for the
zeroinertia case (St = 0). Theoretically, a*=l where the
trajectories hit, but is otherwise zero. The various
methods are approximately consistent with this result.
Nonphysical fluctuations again exist for the ensemble bin
and weighted averaging methods because each bin does
not have a high number of particles hitting to give accurate
results. The areabased method and LCDE method gives
a smooth concentration ratio that closely follows unity
despite a small number of trajectories.
2 LCDE
D area
ensemble
Figure 13. Normalized concentration
comer flow for Ax=0.1, Np=50, St = 10
at vertical wall of
1
To consider the potential performance on angles
surfaces, a collection plane set as a simple ramped surface
with a uniform slope of 1 was placed around the comer to
test the behavior of various methods. For the ensemble
and weightedaverage bin methods, 500 particles are
injected with a collection bin size of 0.02ymax@mj to ensure
accuracy. For the areabased and LCDE methods, 20
particles were used. For St=0 and St = 1, the trajectories
are plotted in Fig. 14. For the St=0 case, the particle
trajectories follow the flow field. For the finite inertia
case, the particle trajectories are not affected by the comer.
Hence the trajectories are relatively straight.
1 I
weighted
analytical
u0 0.5 1 1.5 2
Y/Ymax@inj
Figure 15. Concentration ratio for a ramped surface with
ramp slope = 1 and St = 0 to ensure low numerical error
500 particles were used for the ensemble and weighted bin
methods, though only 20 particles were needed for the
areabased and LCDE methods.
A finiteinertia case is also considered for the ramped
comer flow. The concentration ratios are plotted against
the vertical collection plane location in Fig. 16 for St = 1.
All four methods agree well with each other. The
binbased methods (ensemble and weighted) only show
minor differences of concentration along the collection
plane because of there are more particles hitting each bin
and the accuracy of a* is increased. The areabased and
LCDE methods do not show such problems.
a
:.: [ _[. # aA k A .E ,
0.5 1 1.5 2
max@inj
Figure 16. Concentration ratio for St=1
number of particles are used as Figure 15.
with the same
2 1.5 1
x
(b)
Figure 14. Trajectories for a ramped surface with ramp
slope = 1 a) St = 0 and b) St = 1
Particle Concentration in Cylinder Flow
For the incompressible, inviscid twodimensional flow
described by Eq. 47, the various methods are investigated
to evaluate concentration ratio (a*) for vertical collection
0.5
planes. Trajectories for particles entering at the left
boundary of the computational domain (x=1.5D) and
exiting through the right boundary of the domain (x= 1.5D)
are plotted for different Stokes numbers in Fig. 17. The
initial condition of the particles are the same as the fluid
and the particles march with a time step of 102. For the
zeroinertia (St=0), all of the particles go around the
cylinder so the particle velocity is the same as the flowfield
velocity. For St = 1, the particles around the edge of the
cylinder do not return to the centerline leaving a void wake.
For St = 10, the particles are not significantly influenced by
the flowfield and exhibit little turning.
1.5
0.5
y/D 0
0.5
1
1.5
1.5
1
0.5
y/D 0
0.5
1 0 1
x/D
(a)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
The zeroinertia case can be used as a reference to test
the accuracy of different methods since this case has zero
divergence, so the concentration ratio should be uniform
throughout the flow field. In Figure 18, different bin
sizes were investigated with the same number of particles
(Np=101) and flow conditions as Fig. 10 where the
particles are collected at the plane where x/D = 0. When
the bin side is y/D = 0.05, the ensemble and weighted
average methods give a* with substantial spatial diffusions.
In contrast, the areabased and LCDE methods are
consistent with the analytical result except near the edge of
the cylinder. This is due to the time step induced
numerical error in the position at the collection plane. If
smaller time steps are used this higher a* can be
eliminated. When the binsize is halved, the areabased
and LCDE methods are unaffected but again show an error
near the discontinuity at y/D = 1/2. This problem does
not appear in the binbased methods though nonphysical
fluctuations remain indicating that such methods are more
accurate as resolution decreases. .
2
1.5
a*
1
0.5
00
0
1 0
x/D
(b)
a 1
0.5
0
0
0.
y/D
E * LCDE
E area
E ensemble
S weighted
analytical
0.2 0.4 0.6 0.8 1
0.2 0.4 0.6 0.8 1
Figure 18. Normalized concentration bin
for x/D = 0, Np = 101, St = 0, and Ay/D
0.025
resolution study
= a) 0.05 and b)
1.5
1 0 1
x/D
(c)
Figure 17. Particle trajectories for St=l and two flux
planes at x/D=0 and x/D=1.5 at a) St=0, b) St = 1, andb) St
= 10 all for ymax@,=D
Next, a study for particle resolution was conducted for
the ensemble and weighted average methods with
zeroinertia and uniform bin size of Ay/D=0.05. The
particles are collection at the top of the cylinder (x=0) to
compute a*. For the ensemble method, Fig. 19 shows
that for a low number of particles released the numerical
prediction is greatly lower than the analytical solution.
As the number of particles released from the injection
plane increases, the prediction using the ensemble method
converges to that of the analytical solution. This result is
consistent with the particle resolution study for the corer
flow and shows that the ensemble bin method need a high
number of particles per bin to obtain accurate results.
Np = 5
Np =50
1 Np = 500
* analytical
0.5
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
concentration is high just above the surface as the particles
trajectories become quite close. All of the methods are
consistent away from the cylinder.
2
a*
1
"VO~ 0.2 0.44 0.6 0.8 1 i
y/D
Figure 19. Normalized concentration particle resolution
study for x/D = 0 Ay/D = 0.05, St = 0, and Np = 5, 50, 500
for the ensemble bin method.
For the weighted average method, a* predictions give
similar trends to that of the ensembleaverage method but
with improved convergence. In Fig. 20, when a low
number of particles is released, a* is significantly lower
than the analytical solution. As the number of particles
released increases, the numerical results of a* approach the
analytical solution value.
0.2 0.4 0.6 0.8 1
Figure 21. Normalized concentration for Ay/D
= 50 at x/D=0, and St = 1
0.05, Np
For St=10, the concentration ratio at the collection
plane tends back to unity (Fig. 22) as there is little
trajectory curvatures. For binbased methods, even fewer
particles hit the collection plane which results in greater
statistical error. For the weighted averaging method,
because of the shape function only a few bins getting
concentration at the collection plane so there is significant
diffusion. The areabased method once again closely
agrees with the LCDE method and yields almost constant
concentration values at the impact plane.
2 i
a*
1
0.5
0.2 0.4 0.6 0.8 1
u u.2 u.4 U.o U.8 1
Figure 20. Normalized concentration particle resolution
study for x/D = 0 Ay/D = 0.05, St = 0, and Np = 5, 50, 500
for the weighted method
For the case St = 1, the a* results are shown in Fig. 21.
Because some of the particles hit the surface of the cylinder,
only a fraction of the released particles impact the
collection plane. This explains the fluctuations in a* for
the ensemble method. The weighted method gives a more
smooth behavior because of the shape function weighting.
The areabased method gives results almost identical to
that of the LCDE method with and captures the sharp
discontinuity at y/D=0.5. At this location, the
Figure 22. Normalized concentration for Ay/D
= 50 at x/D=0, and St = 10
0.05, Np
The collection plane is next moved to the wake of the
cylinder (x/D=1.5). The zeroinertia case is first tested with
the analytical solution. The results are displayed in Fig.
23. For bins close to y/D = 0, the ensembleaverage and
weighted methods give fluctuations and errors near the
centerline but trend to the analytical solution away from
the centerline of the cylinder. The deviations of
concentration for the area and LCDE methods near the
centerline of the cylinder are an artifact of the errors
introduced by the time step in the trajectories. As the
impact locations get away from the centerline, these
method are almost identical to the analytical solution.
t
LCDE
Area
ensemble
weighted
analytical
0.2 0.4 0.6 0.8 1
y/D
Figure 23. Normalized concentration
= 50 at x/D=0, and St = 0
for Ay/D = 0.05, Np
For the St =10 case the results of a* plotted against
the collection location normalized by the cylinder diameter
can be seen in Fig. 24. The ensembleaverage bin method
has only a few values due to a large amount of particles
hitting the surface of the airfoil and this gives a big
variation in a* compared to the areabased and LCDE
methods. The weighted averaging method has the same
behavior as the ensemble method but the shape function
smooths out the results. The areabased and LCDE
method shows good agreement with each other indicating
nearly unity concentration which tends to the analytical
solution of unity for infinite Stokes number.
1.5
a
1
0.5
0 0.2 0.4 0.6 0.8 1
Figure 24. Normalized
Np= 50 at x/D=0, and St:
y/D
(14c)
concentration for Ay/D
S10
0.05,
Conclusions
Simulations were conducted for particle trajectories
and concentrations for two flows and several particle
Stokes numbers. For linear flowfield of a corner, the
error in velocity does not affect the integrity of the
numerical methods. For the nonlinear flowfield around a
cylinder, velocity error is substantial and higherorder
methods with small time step are needed for high accuracy.
Of the tow statistical methods, the weighted average
method is most accurate especially when a high number of
particles is hitting the collection bin. It requires less
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
number of particles per bin to obtain accurate results than
the ensembleaverage method since the shape function will
smooth out the distribution but can sometimes give
inaccurate results at the edge bins due numerical diffusion.
The ensembleaverage bin method is accurate when a high
number of particles is hitting each bin. It has the
advantage of being simple and applicable to complex flow
but requires the most number of particles to be accurate.
The areabased method is the most efficient method
in computing the particle concentration ratio for the
present flows. It uses calculations that require only a few
particles' trajectories at the present time step and is not as
sensitive to the numerical trajectory errors as the LCDE
method. However, the areabased method is not
applicable for unsteady flows and also yields increased
errors near discontinuities.
Like the areabased method, the LCDE method
computes the particle concentration ratio with the use of
only a few particles for flows with either linear or
nonlinear flow velocity variations. This method requires
high accuracy in trajectory computation otherwise the
numerical errors will be introduced in the computation for
the concentration ratio. It also requires solving an
additional differential equation along the particle path and
thus needs more memory. Importantly, this method can
be extended to unsteady flows. Further study for
complex flowfields (unsteady, with added mass effects,
etc.) is recommended to test the integrity of all the
methods.
Acknowledgements
This work was funded by NASA Glenn Research Center
icing branch for their Icing Research Tunnel, with Colin
Bidwell as technical manager.
References
Barton, I.E., ExponentialLagrangian Tracking
Schemes Applied to Stokes Law, J. Fluids Eng, Vol. 118, pp.
8589 (1996)
DeAngelis, B., Loth, E., Lankford, D. and Bartlett,
C.S. "Computations of turbulent droplet dispersion for wind
tunnel icing tests," AIAA Journal of Aircraft, Vol. 34, No.2,
pp. 213219, MarchApril (1997)
Dorgan, A. & Loth, E. "Dispersion of bubbles in a
turbulent boundary layer," ASME Fluids Engineering
Division Summer Conference (2009)
Healy, D.P and Young. J.B. "Full Lagrangian Methods
for Calculating Particle Concentration Fields in Dilute
GasParticle Flows", Proc. R. Soc. A, 2005, Vol. 461, pp.
21972225 (2005)
Loth, E. "Numerical approaches for motion of
dispersed particles, bubbles, and droplets," Progress in
Energy and Combustion Sciences, Vol. 26, pp. 161223.
(2000)
Maxey, M.R. and Riley, J.J. "Equation of motion for
a small rigid sphere in a nonzeroinertia," Phys. Fluids, Vol.
26, No. 4, pp. 883889 (1983)
Osiptsov, A.N. 1998 Modified Lagrangian method for
calculating the particle concentration in dustygas flows
with intersecting particle trajectories. In Proceedings of
Third International Conference on Multiphase Flow
ICMF98, Lyon, France, Paper 236 (1998)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
G. A. Ruff and B. M. Berkowitz, "Users manual for the
NASA Lewice Ice Accretion Prediction Code (LEWICE),"
NASA CR 185129 (1990)
Schiller L., Naumann Z. "A drag coefficient
correlation" Z. Ver. Deutsch. Ing., 77318 (1935)
Sivier, S.A., Loth, E., and Baum, J.D. "Dusty shock
flow with unstructured adaptive finite elements and
parcels," AIAA Journal, Vol. 34, No. 5, pp. 10781080
(1996)
