7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Lagrangian Simulations of Solid Particle Motion in A Curved Duct
Sassan Etemad*
Cooling Analysis Dept 26463, Volvo 3P, Volvo Truck Corporation
Gothenburg SE40508, Sweden
sassan.etemad@volvo.com
Keywords: Lagrangian simulations, curved duct, secondary flow, turbulence model
Abstract
Lagrangian computational simulations are carried out for an air flow in a superellipse sectioned curved duct with solid particles
at a Reynolds number of 105620 and a Dean number of 61500. The boundary layer impact on the onset of the secondary flow has
been investigated. Due to the importance of the role of the boundary layer, alternative wall treatment approaches were used. The
performance of eddy viscosity turbulence models with wall function and damping function, i.e. high Reynolds number model
and low Reynolds number models, respectively, have been investigated. The secondary flow and its impact on the solid particle
motion are studied. The gained knowledge can be used for simple, costefficient purification of the fluid with low pressure drop
penalty by removing the particles using a customized particle trap. It is believed that CFD can be used for similar purification
studied of confined contaminated flows in curved conduits.
Introduction
In most of the fluid purification devices, due to the
geometrical constraints, routing of the fluid imposes the
usage of curved ducts with associated pressure drop. The
curvature leads to the centrifugal force which is in
equilibrium with pressure force in general, except in the
boundary layer where the circumferential velocity decreases.
This in turn leads to a reduction of the centrifugal force which
is linearly linked to the circumferential velocity. The
opposing pressure force, however, remains unchanged. It
initiates a secondary motion which has the same orientation
as the pressure force in the near wall region, i.e. towards the
centre of curvature. Consequently, a pair of secondary
vortices is generated near the side walls of the duct (not inner
wall near the centre of curvature, nor the outer wall far away
from the centre). The motion of the fluid and its particle
content is affected by parameters such as duct shape,
Reynolds number, Dean number, particle properties and
particle concentration. By studying such mechanisms, in
particular by using CFD simulations, one can learn more
about particle behaviour in such flow fields. In the search for
compact design, cost reduction and energy saving, the
existing secondary flow can be favourably used for efficient
decontamination.
Nomenclature
D Crosssection diameter[m], here 0.08 m, see Fig. 1
Ewf wall roughness coefficient [], here 9.0, see eq. (17)
H height of the duct crosssection [m], here 0.0889 m
k turbulent kinetic energy [m2/s2], = 0.5 u,u,
P turbulent production [m2/s3], see Table 3
R Crosssection radius[m], here 0.04 m,=D/2, see Fig. 1
Ri inner radius [m], see Fig. 1
Ro outer radius [m], see Fig. 1
Re Reynolds number = UmH/V
R, turbulent Reynolds number =k2/(vE)
R, turbulent Reynolds number = k2 (vt)
R, bend midline radius [m], =1.5D=0.5(Ri+ Rio)
S strain invariant [],= r(2 Sij Sj)5
S strain invariant [],= (k/E) (0.5 Sij S)0o5
Si strain rate [s1], = 0.5 (dUi,/dj + dUj/,di)
t time [s]
U, time averaged velocity in xi direction [m/s]
Um mean flow velocity in the duct [m/s]
ui fluctuating velocity in xi direction [m/s]
u+ dimensionless streamwise velocity [], =Ui/u,
u, friction velocity [m/s], =( l/p)0O.
W vorticity invariant [], = r(2 Wi Wj).5
V vorticity invariant [], = (kl/) (0.5 Wi Wy)5
W, vorticity rate [s1], = 0.5 (dU,/ix dUj,/d,)
X, streamwise coordinates or X, Y Z [m]
xi Cartesian coordinates, or x, y, z [m]
Y+ dimensionless wall distance [], = X2ul/v
Greek Letters
cij Kronecker delta []
e turbulent dissipation rate [m2/s3], = vaui, juj
T the isotropic part of e, E = E2v(k/laxi)2
i von Kirmin's constant [], here 0.4153
p/ molecular viscosity [Pa s]
/4 eddy viscosity [Pa s]
v kinematic viscosity [mn/s]
p density [kg/m3]
r time scale [s] = k/c
Wall shear stress [Pa] = (,UdU1/d2) aty=O
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Problem Statement
The geometry concerned is super ellipse sectioned 1800
Ubend ducts. The shape of the crosssection (perpendicular
to the plane of the Ubend midline) is described as a
combination of a rectangular core with two halfcircles on
each side, as shown in Figure 1. The half circle diameter is
D=0.08 m, which is equal to the rectangle height, while the
rectangle width is only D/4. The bend midline radius, r, was
1.5D, as illustrated in Figure. 1. The Reynolds number based
on D was ReD=Un/v=56000 where Un,=20 m/s is the mean
velocity and the Dean number was
(ReDUv) I _' ))0.5=61500. To avoid confusion, it should
be noted that in these definitions, the height, D, is used
instead of the hydraulic diameter Dh. The bend is provided
with a 12.5D long straight inlet part and a 2.5D straight outlet
part as shown in Figure 1. Based on earlier work by Etemad
et al. [1] it is known that the inlet length is not sufficiently
long for the flow to reach fully developed status. It was,
however, found that 12.5D is sufficiently long for the
boundary layer to develop to an extent which has significant
impact on the flow in the bend. Two different crosssection
grid types were made. Model A, made for the LowRe model
has a finer grid near the wall. Therefore, for this case an
expanding grid was made for the near wall region, up to a
distance D/4 away from the wall. In the offwall region, i.e.
from a distance D/4 from the wall, a coarser, non expanding
(equidistant) grid was used as shown in Figure 2. For case B,
made for the HighRe model a coarser, non expanding
(equidistant) grid was used throughout the crosssection, i.e.
even in the nearwall region. The grid details are specified in
Table 1. Due to the symmetry of the geometry all the grids
were made only in halfmodels and symmetryplane
boundary conditions were applied on the geometry midplane.
At the inlet of the upstream straight duct a uniform velocity
profile with the inlet mean velocity of U,, =20m/s was
applied as the inlet boundary condition. The I' value raged
between 0.1 and 1.3 for the LowRe model grid and between
24 and 99 for the HighRe model grid.
Figure 2: Grids for the halfmodels. Fine crosssection
grid for case A and coarse crosssection grid for case B.
Table 1. Grid specifications
Case Cells in crosssection. Mean 7+
(inlet+bend+outlet) Y+ range
A 3040 (184+90+80) 0.88 0.11.3
B 1000 (184+90+80) 72.7 2499
Figure 1: Geometry description.
Turbulence Models
The equations for the conservation of mass and momentum
can be expressed as given by, e.g., Versteeg and
Malalasekera [2]. The turbulence models used in this work
are Chen's linear highRe ke model [3], and Suga's cubic
lowRe ke model [45]. The latter employs a damping
function to resolve the boundary layer down to the solid
wall in contrast to the highRe model which uses the
conventional wall function. Both models are briefly
explained below.
It is well known that linear eddy viscosity models can not
account for the nonisotropic nature of the turbulence. This
deficiency can be avoided by adding series expansions of
functional terms to the expression for linear eddy viscosity
model to obtain nonlinear eddy viscosity models.
Depending on the desired theoretical accuracy level, several
functional terms, referred to as quadratic and cubic terms,
are added so that the general expression for the Reynolds
stresses becomes as presented in eq. (1).
2
puu = . 2,u,S
+4Cflzu(SS, S,,S, )+4C2,uz(WS, +WS )+4Cyzur(WW,1 WkW6S, )
+. ?(SkW +s WI, )S + '*. ,(SkSk WklWk (1)
The first line in eq. (1) contains only the linear version while
by adding the second line, which contains the quadratic
products of strain and vorticity rate, S, and Wy, the quadratic
model is obtained. Likewise, by adding the third line,
containing the cubic products of S, and W,, the cubic model
is obtained. This means that in the present study, for Chen's
model, only the first line and for Suga's model all the three
lines are used. Note that r= k/e is the turbulent dissipation
time scale. The coefficients C1Cs, the turbulent viscosity, /,,
and the damping function (in Suga's model), f,, for the used
models, are expressed below. For Chen's model f, is 1.
SCpk2
, = f, 
Cj = CNLJC2
Cj= CNLJ
C 0.3 [1exp(
+ 0.3577
for J=4,5
for J=1,2,3
0.36exp(0.7517))]
S= max(S,W) (6)
f lexp( (R,i/ ,R,' /n') (7)
The values for the empirical coefficients for C, and Cj are
given in Table 2.
Table 2: Values for the coefficients in C, & Cj
Ao A, A2 A3 CNLI CNL2 CNL3 CNL4 CNL5
.667 1.25 1. 0.9 0.1 0.1 0.26 10. 5.
Both of the used turbulence models in this work are ke
model variants, hence, two transport equations, one for the
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
turbulent kinetic energy, k, and another one for its
dissipation rate, c, should also be solved. These two
equations can be expressed in a condensed general form as
presented below.
8(pk) 8 ,}k
a U k)= +a xak +,Ppe (8)
at Dxj a L ) k)j
a(pe) \ e A C, E2
ux e  + +Cel Ce2+
at acj aJ i 0x )O'jd k k
U2 2
+ Ce + E + Y, (9)
p k
As can be seen the kequation has its conventional form
while the three last terms in the eequation are additional
terms. The first of these terms belongs to Chen's model and
the last two terms, E and Y, are used only in Suga's model.
In addition, for Suga's model, in the cequation, e is
replaced by E where is the isotropic part of e as
defined below. T is zero at the wall.
= 2v 2 (10)
( ,
The coefficients in Suga's model are C9= 0.09, C1=1.44,
C,2=1.92, Cr=1.44 (0 for P<0) and C4=0.33. P, fi, Cs,
E and Y, are expressed as given in Table 3 and below,
depending on the choice of turbulence model.
Table 3. Expressions used in k & e equations
Model P fi Cs E Y,
Chen S 1 0.25 0 0
0U,
U 0 E Y
Suga puu, k, f, 0 E Y,
where
f, = 0.3exp ( R)
E = 0.0022 SK, k [ u
Y = max k.83 k ,3/2 y 2 ]
Y, = max 0.83 \ k 1 I ,0
L [2.5y I[2.5Ey k
The differences in the turbulence models explained above
give them their particular characters which are expected to
affect the flow they predict.
In Chen's model, in addition to the dissipation time scale,
k/e, which is used in the standard ke model, the production
time scale, k/P, which appears in eq. (9), is also used. This
extra time scale is claimed to allow the energy transfer
mechanism of turbulence to respond to the mean strain rate
more efficiently. This causes the extra constant, C,5, in the
eequation which is expected to improve jet flows, elliptic
flows and confined swirling flows. This model uses the
conventional logarithmic wall function, as expressed below,
for the walladjacent cells.
1
u = n( E Y+) (14)
K
where Ewf =9.0 is the wall roughness coefficient and icis the
von Kirmin's constant and here i=0.4153.
Suga's model is nonlinear and it has both high and lowRe
versions. Only the latter is used here. This model uses the
stress invariant as a wall detector in the nearwall region. In
addition, the role of Yc is to reduce the otherwise excessive
levels of length scales in the nearwall region in separated
and stagnating flows [45].
For the temperature field the simple eddy diffusivity model,
SED, was used. The turbulent heat flux is expressed as
puh /, Sh
CJ JDxJ
0,h =0.9
Numerical Solution Procedure
The governing equations are solved by an implicit finite
volume nonstaggered collocatedd) solver. The SIMPLE
algorithm was used for the pressurevelocity coupling.
Earlier work by Etemad et al. [1], showed that MARS
(second order) performed better than the Upwind scheme
(first order) while the QUICK scheme (third order) did not
give any significant improvement compared to MARS.
U/Urn
0.6
0.4
0.2
Ri Ro
Figure 3: Normalized streamwise velocity profiles on
the symmetry plane at 10D (upstream the bend) at the 0
(bend inlet) and 900 (in midbend). Ri and Ri are the inner
and outer radius respectively.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Hence, MARS was used throughout this work.The fluid was
air and constant physical properties were assumed. The
particles were assumed to be spherical with a diameter of
25pum and a density of 2200 kg/m3. Lagrangian coupled
simulation was carried out with perfectly rebounding
(reflecting) walls. Maximum particle tracking time was 2s.
The fluid/particle mass fraction was 186220. The gravity
impact was not accounted for. All calculations were
performed by using the code STARCD [6].
Results and Discussion
The results obtained from the simulations are presented in
two forms. In order to quantify the results, profiles of
normalized streamwise velocity on the symmetry plane at
various streamwise positions are presented in Figure 3. The
profiles show results from 10D (upstream the bend), 0 and
900. These positions are selected due to the importance of the
velocity profile and the boundary layer thickness on the onset
of the secondary flow in the curved duct.
Figure 4: Velocity magnitude distribution in the bend
for case A (left) and B (right). The velocity is normalized
with the mean velocity U,n.
To give a qualitative picture of the flow in the bend, Figures
46 are provided. They show sequences of plots of velocity
magnitude contours and velocity vectors at the streamwise
positions 10D, 0 and 900 and in the entire bend. The plots
are viewed by an observer facing the duct inlet and travelling
streamwise in the duct with the sidewall on his right and the
symmetry plane on his left. The upper and lower edges of the
plots indicate the outer and inner radius, respectively. A
review of the velocity profiles in Figure 3 shows that Suga's
model and Chen's model gave different results.
Case B
" 10D
e 0 degrees
x 90 degrees
A
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Case A
Case A
Case B
Case B
R,
at O0D
Velocity/Un,
1.3
1.2
1.1
1.0
0.9
0.8
0.7 R
0.6 Ro
0.5
0.4
0.3
0.2
0.1
0.0
Ri
Ri :
Ro
at 00
at 900
' : t  ':: :
R i
Ro
at 1800
RI
Figure 5: Velocity magnitude distribution across
positions at 10D (upstream the bend), 0 (bend inlet), 900
(midbend) and 1800 (bend outlet). Ri and Ri are the inner
and outer radius respectively.
It is obvious that the flow develops differently in the bend in
these two cases. The high speed flow core comes closer to the
duct inner radius in both cases. This feature is, however,
more pronounced for case B than for case A. This can also
clearly be seen in the velocity magnitude distribution plots in
Figure 5. Figure 6 which shows the velocity vectors at he
same crosssections as the ones in Figure 5 contribute to the
understanding of the occurring feature in the bend. It appears
that the high speed core, prior to entering the bend, comes
closer to the inner radius. However, somewhere between 0
and 900, this tendency reverses and the by 1800 it becomes
more pronounced.
Figure6: Velocity vectors and illustration of the
secondary flow across positions at 10D, 0, 900 and 180.
To assist understanding of the motion descriptive arrows
have been added to the plots.
As Ito [7] describes, a flow with uniform streamwise velocity
profile approaching the curved duct has, behaves similar to
the free vortex flow, where the velocity near the inner radius
is higher than near the outer radius. This is the case here
because due to the rather short inlet section the approaching
flow is still rather uniform and it has developed a thin
boundary layer. In the bend, the pressure gradient, which is
oriented towards the bend centre, is in equilibrium with the
centrifugal force with the opposite orientation. Near the
sidewall (opposite from the symmetry plane), due to the
friction, the centrifugal forces decrease, and thereby, the
pressure forces prevail. Therefore they start to displace fluid
from the sidewall towards the inner wall. This is what can be
observed as a secondary motion at positions 90 and 180.
This makes the flow to move alongside the inner wall
towards the symmetry plane to eventually depart from the
inner wall and move towards the duct core. This mechanism
at 10D
Ri
Ro
at 00
at 900
at 1800
Iif~
;;;+
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
referred to as the pressuredriven secondary motion generates
a main longitudinal vortex accompanied by a number of
smaller satellite vortices. The thickness of the boundary layer
has a strong influence on the strength of this motion. The
plots in Figures 56 show this motion. This effect gets .
stronger with increasing inlet length leading the more high
speed fluid to move closer to the outer wall. Here, due to the .
rather short inlet section prior to the bend, this phenomenon J
is rather weak. Nevertheless, it exists, and it has just managed
to move away the high speed core from the inner radius and
bring some low speed fluid between the high speed core and Case A ,
the inner wall. This phenomenon can be observed well in 
Figure 4 as well. It is clear that the simulation using Suga's 
LowRe model for case A predict a stronger secondary '
motion that the one using Chen's HighRe model for case B. .
This is directly linked to the fact that Suga's model predicts a '
thicker boundary layer in the bend than Chen's model, as can A
be observed in Figure 3. Consequently the impact of the two
different flows on the particle motion in and after the bend
will be significant. I i
Figure 7 shows particle tracks in the computational domain
for both cases. The bend region has been enlarged and view
at a slightly different angle. Due to the very low particle/fluid
mass fraction the impact of the particle motion on the fluid /
motion was found to be quite negligible. .
In case A, the particles are under the influence of both the. ,,
centrifugal force due to the bend curvature and also the
strong secondary vortex in the second half of the bend.
Therefore, the particles gather rather well towards the outer
wall of the bend in the first half of the bend. In the second
half, however, due to the increasing strength of the secondary
vortex, the particles start to obey the rotational motion of the
vortex. This in turn gives a secondary centrifugal force to the
particles that pushes them near the wall and the symmetry. ":"
plane, leaving the mid region rather clean. Most particles can ;
be found in a region about 0.05D away from the wall and the
symmetry plane at the outlet of the model. Only 4% of the. :
particles are in the mid part of the crosssection at the outlet.
In case B, the particles are seemed to be influenced by only Case B
the centrifugal force due to the bend curvature. Very weak 
signs the influence of the secondary vortex is observed. # t. i.
Hence, the particles gather in a concentrated stoke on the .
outer wall leaving almost all of the duct crosssection clean.
From the point of view of purification objectives, this gives a 3. .
very optimistic picture leading to the conclusion that
providing that the concentrated particle stroke is removed,
most of the fluid can be purified at a very low expense.
However, the results from case A show that the secondary
vortex in the bend prevents the particles from gathering as
nicely as in case B. On the other hand, due to its rotational
motion it throws away the particles towards the wall (and the
symmetry plane). Over 96% of the particles were in the
nearwall region. This gives a great purification potential. By
removing the nearwall particles a high degree of purification 
can be achieved. It should be noted that such a solution is
efficient for both case A and case B.
Figure 7: Particle tracks in the entire computational
domain and in the bend section.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
The turbulent flow with solid particles in a super ellipse
sectioned Ubend duct was investigated using two different
variants of the kemodel, Chen's (highRe) model and Suga's
cubic lowRe model. It was found that the boundary layer
upstream the bend has a significant impact on the flow in the
bend and consequently on the particle motion. Suga's model
predicted a thicker boundary layer and thereby a resulting
secondary vortex that significantly impacts the particle
motion. It is concluded that great care should be taken when
choosing turbulence model to achieve high accuracy in
prediction of the boundary layer. Suga's cubic lowRe model
is preferred for such analysis. It requires, however, a more
demanding computational grid and needs much more
computational power. The flow was found to display a
complex motion with a number of rotating vortices
accompanying the main streamwise flow. Over 96% of the
particles were kept in the nearwall region when using Suga's
model. With Chen's model, the particle concentration near
the wall was even greater. This enables efficient particle
decontamination with limited penalty.
Acknowledgements
The author wishes to express his thanks to Reimer Ryrholm
for valuable discussions about purification and to Volvo 3P
for financing this work.
[1] Etemad, S., Sund6n, B. and Daunius, O., 2006, "Turbulent
Flow and Heat Transfer in a SquareSectioned UBend",
Progress in Computational Fluid Dynamics, Vol. 6, Nos.
1/2/3, pp. 89100.
[2] Versteeg, H.K. & Malalasekera, W., 1995, "An
introduction to computational fluid dynamicsThe finite
volume method". Longman Group Ltd.
[3] Chen, Y.S. & Kim, S.W., 1987, "Computation of turbulent
flows using an extended ke turbulence closure model".
NASA CR179204.
[4] Craft, T.J., Launder, B.E., & Suga, K., 1993, "Extending
the applicability of eddy viscosity models through the use
of deformation invariants and nonlinear elements", Proc.
5th Int. Symp. on Refined Flow Modelling and
Turbulence Measurements, pp. 125132.
[5] Craft, T.J., Launder, B.E., & Suga, K., 1996, Development
and application of a cubic eddy viscosity models of
turbulence", Int. J. Heat and Fluid Flow, 17, pp. 108115.
[6] STARCD, "'l ,t..,...1.,, STARCD VERSION 3.2",
2003, CD Adapco Group.
[7] Ito. H., Flow In Curved Pipes, JSME Intl. Journal, Vol.
30, No. 262. pp. 543552, 1987.
Conclusions
References
