The effect of particle density in turbulent channel flow laden with finite size
particles in semi-dilute conditions
W. Fornari, A. Formenti, F. Picano, and L. Brandt
Citation: Physics of Fluids 28, 033301 (2016); doi: 10.1063/1.4942518
View online: http://dx.doi.org/10.1063/1.4942518
View Table of Contents: http://scitation.aip.org/content/aip/journal/pof2/28/3?ver=pdfcov
Published by the AIP Publishing
Articles you may be interested in
On the relative rotational motion between rigid fibers and fluid in turbulent channel flow
Phys. Fluids 28, 013301 (2016); 10.1063/1.4937757
Slip velocity of rigid fibers in turbulent channel flow
Phys. Fluids 26, 063302 (2014); 10.1063/1.4881942
The effect of neutrally buoyant finite-size particles on channel flows in the laminar-turbulent
transition regime
Phys. Fluids 25, 123304 (2013); 10.1063/1.4848856
Anisotropy in pair dispersion of inertial particles in turbulent channel flow
Phys. Fluids 24, 073305 (2012); 10.1063/1.4737655
Numerical simulation of particle-laden turbulent channel flow
Phys. Fluids 13, 2957 (2001); 10.1063/1.1396846
Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 185.5.62.153
On: Tue, 01 Mar 2016 17:58:53
PHYSICS OF FLUIDS 28, 033301 (2016)
The effect of particle density in turbulent channel flow
laden with finite size particles in semi-dilute conditions
W. Fornari,
1
A. Formenti,
2
F. Picano,
2
and L. Brandt
1
1
SeRC and Linné FLOW Centre, KTH Mechanics, SE-100 44 Stockholm, Sweden
2
Department of Industrial Engineering, University of Padova,
Via Venezia 1, 35131 Padova, Italy
(Received 14 October 2015; accepted 5 February 2016; published online 1 March 2016)
We study the effect of varying the mass and volume fraction of a suspension of rigid
spheres dispersed in a turbulent channel flow. We performed several direct numerical
simulations using an immersed boundary method for finite-size particles changing
the solid to fluid density ratio R, the mass fraction χ, and the volume fraction φ. We
find that varying the density ratio R between 1 and 10 at constant volume fraction
does not alter the flow statistics as much as when varying the volume fraction φ at
constant R and at constant mass fraction. Interestingly, the increase in overall drag
found when varying the volume fraction is considerably higher than that obtained for
increasing density ratios at same volume fraction. The main effect at density ratios
R of the order of 10 is a strong shear-induced migration towards the centerline of
the channel. When the density ratio R is further increased up to 1000, the particle
dynamics decouple from that of the fluid. The solid phase behaves as a dense gas and
the fluid and solid phase statistics drastically change. In this regime, the collision rate
is high and dominated by the normal relative velocity among particles.
C
2016 AIP
Publishing LLC. [http://dx.doi.org/10.1063/1.4942518]
I. INTRODUCTION
The transport of particles in flows is relevant to many industrial applications and environmental
processes. Examples include sediment transport in rivers, avalanches and pyroclastic flows, as well
as many oil industry and pharmaceutical processes. Often the flow regime encountered in such
applications is turbulent due to the high flow rates and it can be substantially affected by the pres-
ence of the solid phase. Depending on the features of both fluid and solid phases, many different
scenarios can be observed and the understanding of such flows is still incomplete.
The rheological properties of these suspensions have mainly been studied in the viscous Stoke-
sian regime and in the low speed laminar regime. Even limiting our attention to monodisperse rigid
neutrally buoyant spheres suspended in Newtonian liquids, we find interesting rheological behaviors
such as shear thinning or thickening, jamming at high volume fractions, and the generation of
high effective viscosities and normal stress differences.
1–3
It is known that the effective viscosity
of a suspension µ
e
changes with respect to that of the pure fluid µ due to the modification of the
response of the complex fluid to the local deformation rate.
4
In the dilute regime, an expression for
the effective viscosity µ
e
with the solid volume fraction φ has first been proposed by Einstein
5,6
and then corrected by Batchelor
7
and Batchelor and Green.
8
As the volume fraction increases, the
mutual interactions among particles become more important and the effective viscosity increases
until the system jams.
9
At high volume fractions, the variation of the effective viscosity µ
e
is
described exclusively by semi-empirical laws such as those by Eiler and Krigher and Dougherty
1
that also capture the observed divergence at the maximum packing limit,
10
φ
m
= 0.58–0.62. In
laminar flows, shear-thickening or normal stress differences occur due to inertial effects at the
particle scale. Indeed, when the particle Reynolds number Re
a
is non-negligible, the symmetry
of the particle pair trajectories is broken and the microstructure becomes anisotropic, leading to
macroscopical behaviors such as shear-thickening.
11–13
Finally, in the highly inertial regime, the
effective viscosity µ
e
increases linearly with shear rate due to augmented particle collisions.
14
1070-6631/2016/28(3)/033301/19/$30.00 28, 033301-1 © 2016 AIP Publishing LLC
Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 185.5.62.153
On: Tue, 01 Mar 2016 17:58:53
033301-2 Fornari et al. Phys. Fluids 28, 033301 (2016)
Another important feature observed in viscous flows is shear-induced migration. When consid-
ering a pressure-driven Poiseuille flow, either in a tube or in a channel, the particles irreversibly
migrate toward the centerline, i.e., from high to low shear rate regions.
4,15
Interestingly when iner-
tial effects become important, a different kind of migration occurs as the particles tend to move
radially away from both the centerline and the walls, toward an intermediate equilibrium position.
This type of migration was first observed in a tube
4,16
and was named tubular pinch. It is mechanisti-
cally unrelated to the rheological properties of the flow and results from the fluid-particle interaction
within the conduit. The case of the laminar square duct flow has also been studied to identify the
particle equilibrium positions.
17,18
It was found that finite-size particles migrate toward the corners
or to the center of edges depending on the bulk Reynolds number. At high Reynolds numbers (but
still in the laminar regime), some particles were also found in an inner region near the center of the
duct.
Typically, as the Reynolds number is increased, inertial effects become important and the un-
laden flow undergoes a transition from laminar to turbulent conditions. The presence of the solid
phase may alter this process by either increasing or reducing the critical Reynolds number above
which the transition to the turbulent regime occurs. The case of a dense suspension of particles in a
pipe flow has been studied experimentally
19
and numerically.
20
It has been suggested that transition
depends upon the pipe to particle diameter ratios and the volume fraction. For larger particles, tran-
sition shows a non-monotonic behavior that cannot be solely explained in terms of an increase of
the effective viscosity. For smaller neutrally buoyant particles instead, the critical Reynolds number
increases monotonically with the solid volume fraction due to the raise in effective viscosity.
The transition in dilute suspensions of finite-size particles in plane channels has been studied
by Lashgari et al.
21
and Loisel et al.
22
It has been shown that the critical Reynolds number above
which turbulence is sustained is reduced. At fixed Reynolds number and solid volume fraction, the
initial arrangement of particles is important to trigger the transition. Lashgari et al.
23
also investi-
gated numerically a channel flow laden with solid spherical particles at higher volume fractions and
for a wide range of Reynolds numbers. These authors identified three different regimes for different
values of the solid volume fraction φ and the Reynolds number Re. In each regime (laminar, turbu-
lent, and inertial shear-thickening), the flow is dominated by different components of the total stress
(viscous, turbulent, or particle stresses, respectively).
Regarding the fully turbulent regime, most of the previous studies have focused on dilute or
very dilute suspensions of particles smaller than the hydrodynamic scales and heavier than the fluid.
In the one-way coupling regime
24
(i.e., when the solid phase has a negligible effect on the fluid
phase) and limiting our attention to wall-bounded flows, it has been shown that particles migrate
from regions of high to low turbulence intensities.
25
This phenomenon is known as turbophoresis
and it has been shown to be stronger when the turbulent near-wall characteristic time and the
particle inertial time scale are similar.
26
Small-scale clustering has also been observed in this kind
of inhomogeneous flows,
27
leading together with turbophoresis to the formation of streaky particle
patterns.
28
In the two-way coupling regime (i.e., when the mass density ratios are high and the
back-reaction of the dispersed phase on the fluid cannot be neglected), the solid phase has been
shown to reduce the turbulent near-wall fluctuations increasing their anisotropy
29
and eventually
reducing the total drag.
30
When the suspensions are dense, it is of fundamental importance to consider particle-particle
interactions and collisions. Indeed, the chaotic dynamics of the fluid phase affects the rheological
properties of the suspension, especially at high Reynolds numbers. This is known as a four-way
coupling regime. Increasing the particle size directly affects the turbulent structures at smaller
and comparable scales
31
thereby modulating the turbulent field. In a turbulent channel flow, it has
been reported that finite-size particles larger than the dissipative length scale increase the turbulent
intensities and the Reynolds stresses.
32
Particles are also found to preferentially accumulate in
the near-wall low-speed streaks.
32
This has also been observed in open channel flows laden with
heavy finite-size particles. In this case, the flow structures are found to have a smaller streamwise
velocity.
33,34
Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 185.5.62.153
On: Tue, 01 Mar 2016 17:58:53
033301-3 Fornari et al. Phys. Fluids 28, 033301 (2016)
Concerning turbulent channel flows of neutrally buoyant particles, recent studies with φ ≃ 7%
report that due to the attenuation of the large-scale streamwise vortices, the fluid streamwise ve-
locity fluctuation is reduced. When the particles are heavier than the carrier fluid and therefore
sediment, the bottom wall acts as a rough boundary which makes the particles resuspend.
35
Recent
simulations from our group have shown that the overall drag increases as the volume fraction is
increased from φ = 0% up to 20%. This trend cannot be solely explained in terms of the increase of
the suspension effective viscosity. It is instead found that as particle volume fraction increases, the
velocity fluctuation intensities and the Reynolds shear stresses decrease while there is a significant
increase of the particle induced stresses. The latter, in turn, lead to a higher overall drag.
36
As noted by Prosperetti,
37
however, results obtained for solid to fluid density ratios R =
ρ
p
/ρ
f
= 1 cannot be easily extrapolated to other cases (e.g., when R > 1). In the present study, we
therefore investigate numerically the effects of varying the density ratio R of the suspended phase
and consequently the mass fraction χ for different volume fractions. The aim is to understand sepa-
rately the effects of excluded volume and (particle and fluid) inertia on the statistical observables
of both phases. To isolate the effects of different density ratios R on the macroscopical behavior of
the suspension, we consider an ideal situation where the effect of gravity is neglected, leaving its
analysis to future studies.
We consider a turbulent channel flow laden with rigid spheres of radius a = h/18, where h is
the half-channel height (see Picano et al.
36
). Direct numerical simulations (DNS) fully describing
the solid phase dynamics via an immersed boundary method (IBM) are performed as in Lucci
et al.
38
and Kidanemariam et al.
33
among others. First, cases at fixed mass fractions χ = 0.2 are
examined and compared to cases with constant volume fraction φ = 5% and density ratios R rang-
ing from 1 to 10. It is observed that the influence of the density ratio R on the statistics of both
phases is less important than that of an increasing volume fraction φ. The main effects at density ra-
tio R ∼ 10 are shear-induced migration towards the centerline of the channel and slight reduction of
the fluid velocity fluctuations in the log-layer. The results drastically change when further increasing
R (up to ∼1000). It is found that for sufficiently high R (&100), the solid phase behaves as a dense
gas uncorrelated to the details of the carrier fluid flow.
II. METHODOLOGY
A. Numerical method
Different methods have been proposed in the last years to perform direct numerical simulations
of multiphase flows. In the present study, simulations have been performed using the algorithm orig-
inally developed by Breugem
39
that fully describes the coupling between the solid and fluid phases.
The Eulerian fluid phase is evolved according to the incompressible Navier-Stokes equations,
∇ · u
f
= 0, (1)
∂u
f
∂t
+ u
f
· ∇u
f
= −
1
ρ
f
∇p + ν∇
2
u
f
+ f, (2)
where u
f
, ρ
f
, and ν = µ/ρ
f
are the fluid velocity, density, and kinematic viscosity, respectively (µ
is the dynamic viscosity), while p and f are the pressure and a generic force field (used to model
the presence of particles). The particles centroid linear and angular velocities, u
p
and ω
p
are instead
governed by the Newton-Euler Lagrangian equations,
ρ
p
V
p
du
p
dt
= ρ
f
∂V
p
τ · n dS, (3)
I
p
dω
p
dt
= ρ
f
∂V
p
r × τ · n dS, (4)
where V
p
= 4πa
3
/3 and I
p
= 2ρ
p
V
p
a
2
/5 are the particle volume and moment of inertia; τ = −pI +
2µE is the fluid stress, with E =
∇u
f
+ ∇u
T
f
/2 the deformation tensor; r is the distance vec-
tor from the center of the sphere while n is the unity vector normal to the particle surface
Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 185.5.62.153
On: Tue, 01 Mar 2016 17:58:53
033301-4 Fornari et al. Phys. Fluids 28, 033301 (2016)
∂V
p
. Dirichlet boundary conditions for the fluid phase are enforced on the particle surfaces as
u
f
|
∂V
p
= u
p
+ ω
p
× r.
In the numerical code, an immersed boundary method is used to couple the fluid and solid
phases. The boundary condition at the moving particle surface (i.e., u
f
|
∂V
p
= u
p
+ ω
p
× r) is
modeled by adding a force field on the right-hand side of the Navier-Stokes equations. The fluid
phase is therefore evolved in the whole computational domain using a second order finite difference
scheme on a staggered mesh while the time integration is performed by a third order Runge-Kutta
scheme combined with a pressure-correction method at each sub-step. The same integration scheme
is also used for the Lagrangian evolution of Eqs. (3) and (4). Each particle surface is described by
uniformly distributed N
L
Lagrangian points. The force exchanged by the fluid on the particles is
imposed on each lth Lagrangian point and is related to the Eulerian force field f by the expression
f(x) =
N
L
l=1
F
l
δ
d
(x − X
l
)∆V
l
. In the latter, ∆V
l
represents the volume of the cell containing the lth
Lagrangian point while δ
d
is the Dirac delta. This force field is calculated through an iterative algo-
rithm that ensures a second order global accuracy in space. In order to maintain accuracy, Eqs. (3)
and (4) are rearranged in terms of the IBM force field,
ρ
p
V
p
du
p
dt
= −ρ
f
N
l
l=1
F
l
∆V
l
+ ρ
f
d
dt
V
p
u
f
dV, (5)
I
p
dω
p
dt
= −ρ
f
N
l
l=1
r
l
× F
l
∆V
l
+ ρ
f
d
dt
V
p
r × u
f
dV, (6)
where r
l
is the distance from the center of a particle while the second terms on the right-hand sides
are corrections to account for the inertia of the fictitious fluid contained within the particle volume.
Particle-particle interactions are also considered. When the gap distance between two particles is
smaller than twice the mesh size, lubrication models based on Brenner’s asymptotic solution
40
are
used to correctly reproduce the interaction between the particles. A soft-sphere collision model is
used to account for collisions between particles and between particles and walls. An almost elastic
rebound is ensured with a restitution coefficient set at 0.97. These lubrication and collision forces
are added to the right-hand side of Eq. (5). For more details and validations of the numerical code,
the reader is referred to previous publications.
36,39,41,42
B. Flow configuration
We consider a turbulent channel flow between two infinite flat walls located at y = 0 and
y = 2h, where y is the wall-normal direction while x and z are the streamwise and spanwise
directions. The domain has size L
x
= 6h , L
y
= 2h , and L
z
= 3h and periodic boundary conditions
are imposed in the streamwise and spanwise directions. A fixed value of the bulk velocity U
0
is
achieved by imposing a mean pressure gradient in the streamwise direction. The imposed bulk
Reynolds number is equal to Re
b
= U
0
2h/ν = 5600 (where ν represents the kinematic viscosity of
the fluid) and corresponds to a Reynolds number based on the friction velocity Re
τ
= U
∗
h/ν = 180
for the single phase case. The friction velocity is defined as U
∗
=
τ
w
/ρ
f
, where τ
w
is the stress at
the wall. A cubic staggered mesh of 864 × 288 × 432 grid points is used to discretize the domain.
All results will be reported either in non-dimensional outer units (scaled by U
0
and h) or in inner
units (with the superscript “+,” scaled by U
∗
and δ
∗
= ν/U
∗
).
The solid phase consists of non-Brownian rigid spheres with a radius to channel half-width
ratio fixed to a/h = 1/18. For a volume fraction φ = 5%, this radius corresponds to about 10 plus
units. In Figure 1, we display the instantaneous streamwise velocity on four orthogonal planes
together with the finite-size particles dispersed in the domain. Each particle is discretized with
N
l
= 746 Lagrangian control points while their radii are 8 Eulerian grid points long. Using an
Eulerian mesh consisting of 8 grid points per particle radius (∆x = 1/16) is a good compromise in
terms of computational cost and accuracy. We have performed a simulation with a finer mesh (12
points per particle radius, ∆x = 1/24), R = 1 and φ = 5%. We find indeed that the friction Reynolds
number Re
τ
changes by 1%, and the velocity fluctuations change locally at most by 4%.
Reuse of AIP Publishing content is subject to the terms at: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 185.5.62.153
On: Tue, 01 Mar 2016 17:58:53
