1 Introduction
Suspensions of solid particles in fluids can be found in many environmental and industrial applications. Sediment transport in estuaries (Mehta Reference Mehta2014), blood flow in the human body, pyroclastic flows from volcanoes and pulp fibres in paper-making (Lundell, Söderberg & Alfredsson Reference Lundell, Söderberg and Alfredsson2011) are among the examples of flows that deserve further investigations. The presence of solid rigid particles alters the global transport and rheological properties of the mixture in complex (and often unpredictable) ways. Many efforts have therefore been devoted to quantify the effects of the particles in these flows, starting from the simpler case of monodisperse rigid neutrally buoyant spherical particles. The first studies of suspensions under laminar conditions can be traced back to Einstein (Reference Einstein1906, Reference Einstein1911), who analytically derived an expression for the effective viscosity
$\unicode[STIX]{x1D708}_{e}$
of a suspension of rigid spheres in the dilute and viscous limit,
$\unicode[STIX]{x1D708}_{e}/\unicode[STIX]{x1D708}=1+(5/2)\unicode[STIX]{x1D719}$
, where
$\unicode[STIX]{x1D719}$
is the volume fraction and
$\unicode[STIX]{x1D708}$
is the kinematic viscosity of the suspending fluid. A quadratic correction, accounting for particle–particle interactions, was later proposed for higher volume fractions (Batchelor Reference Batchelor1970; Batchelor & Green Reference Batchelor and Green1972). The rheology of dense suspensions is usually characterized by semi-empirical formulae for the effective viscosity (Stickel & Powell Reference Stickel and Powell2005; Guazzelli & Morris Reference Guazzelli and Morris2011).
Inertial effects, yet in laminar flows, are shown to induce significant modifications of the suspension microstructure and to create a local anisotropy responsible for shear thickening (Kulkarni & Morris Reference Kulkarni and Morris2008; Picano et al.
Reference Picano, Breugem, Mitra and Brandt2013), and thus lead to a change of the macroscopic suspension dynamics. Shear thickening and particle migration towards regions of low shear has been observed in several previous studies for dense suspensions at low Reynolds number (Hampton et al.
Reference Hampton, Mammoli, Graham, Tetlow and Altobelli1997; Brown & Jaeger Reference Brown and Jaeger2009; Yeo & Maxey Reference Yeo and Maxey2011). The highly inertial regime was considered in the pioneering work of Bagnold (Reference Bagnold1954), who showed that shearing closely spaced particles induces an effective viscosity that increases linearly with the shear rate, resulting in a normal or dispersive stress in addition to the shear stress (Hunt et al.
Reference Hunt, Zenit, Campbell and Brennen2002). Recently, Lashgari et al. (Reference Lashgari, Picano, Breugem and Brandt2014, Reference Lashgari, Picano, Breugem and Brandt2016) documented the existence of three different regimes when changing the volume fraction
$\unicode[STIX]{x1D719}$
of neutrally buoyant spherical particles and the Reynolds number
$Re$
: a laminar-like regime at low
$Re$
and low to intermediate
$\unicode[STIX]{x1D719}$
, where the viscous stress dominates dissipation; a turbulent-like regime at high Reynolds number and low to intermediate
$\unicode[STIX]{x1D719}$
, where the turbulent Reynolds stress plays the main role in the momentum transfer across the channel; and a third regime at higher
$\unicode[STIX]{x1D719}$
, denoted as inertial shear thickening, characterized by a significant enhancement of the wall shear stress due to the particle-induced stresses.
When the Reynolds number is sufficiently high, the flow becomes turbulent, exhibiting chaotic and multiscale dynamics. The presence of the finite-size particles (particles comparable in size to or larger than the smallest hydrodynamic scales of the flow) can change the turbulent structures at or below the particle size (Naso & Prosperetti Reference Naso and Prosperetti2010; Homann, Bec & Grauer Reference Homann, Bec and Grauer2013). These interactions modulate the whole process by inducing non-trivial effects on the turbulence; see, for example, the studies in homogeneous isotropic turbulence by Lucci, Ferrante & Elghobashi (Reference Lucci, Ferrante and Elghobashi2010) and Fornari, Picano & Brandt (Reference Fornari, Picano and Brandt2016b
), the latter including sedimentation. The first simulations of finite-size particles in a turbulent channel flow, like those discussed here, were performed by Pan & Banerjee (Reference Pan and Banerjee1996), revealing that turbulent fluctuations and stresses increase in the presence of the solid phase. Matas, Morris & Guazzelli (Reference Matas, Morris and Guazzelli2003), Loisel et al. (Reference Loisel, Abbas, Masbernat and Climent2013), Yu et al. (Reference Yu, Wu, Shao and Lin2013) and Lashgari, Picano & Brandt (Reference Lashgari, Picano and Brandt2015) considered the turbulence onset in suspensions of neutrally buoyant spherical particles and reported a decrease of the critical Reynolds number for transition to turbulence in the semi-dilute regime. The simulations by Shao, Wu & Yu (Reference Shao, Wu and Yu2012) revealed a decrease of the fluid streamwise velocity fluctuations due to an attenuation of the large-scale streamwise vortices in a turbulent channel flow. Kidanemariam et al. (Reference Kidanemariam, Chan-Braun, Doychev and Uhlmann2013) considered heavy finite-size particles and showed accumulation in the near-wall low-speed streaks at low
$\unicode[STIX]{x1D719}$
. Picano, Breugem & Brandt (Reference Picano, Breugem and Brandt2015) studied dense suspensions of neutrally buoyant particles in turbulent channel flow up to volume fraction
$\unicode[STIX]{x1D719}=20\,\%$
. These authors showed that the velocity fluctuation intensities and the Reynolds shear stress gently increase with
$\unicode[STIX]{x1D719}$
and then sharply decrease at
$\unicode[STIX]{x1D719}=20\,\%$
, even though the overall drag still increases. They attributed the drag increase to the enhancement of turbulence activity for
$\unicode[STIX]{x1D719}\leqslant 10\,\%$
and then to the particle-induced stresses that govern the dynamics at high
$\unicode[STIX]{x1D719}$
. Costa et al. (Reference Costa, Picano, Brandt and Breugem2016) showed that the turbulent drag of a suspension of spherical particles is always higher than that predicted by accounting only for the effective suspension viscosity. This is attributed to the formation of a particle–wall layer, a layer of particles forming near the wall. Fornari et al. (Reference Fornari, Formenti, Picano and Brandt2016a
) investigated the role of fluid and particle inertia, also in the semi-dilute regime, and show that the excluded volume is responsible for the turbulence modulations, while the particle inertia is negligible for solid to fluid density ratios below 10.
The dynamics of suspension in the presence of non-spherical particles is less well understood (Prosperetti Reference Prosperetti2015), the majority of previous studies in the turbulent regime having dealt with point-like spheroidal particles. Most of these investigations further assume dilute conditions and neglect the feedback on the flow. In this so-called one-way coupling regime, spherical particles display accumulation near the wall (turbophoresis) and preferential sampling of low-speed regions (see e.g. Sardina et al. Reference Sardina, Picano, Schlatter, Brandt and Casciola2011, Reference Sardina, Schlatter, Brandt, Picano and Casciola2012). Turbulent channel flow of non-spherical particles has been investigated by several authors (Zhang et al. Reference Zhang, Ahmadi, Fan and McLaughlin2001; Mortensen et al. Reference Mortensen, Andersson, Gillissen and Boersma2008; Marchioli, Fantoni & Soldati Reference Marchioli, Fantoni and Soldati2010; Challabotla, Zhao & Andersson Reference Challabotla, Zhao and Andersson2015a ), with focus on the particle dynamics. Challabotla, Zhao & Andersson (Reference Challabotla, Zhao and Andersson2015b ) investigated the rotational motion of inertia-free spheroids in turbulent channel flow using equations given by Jeffery (Reference Jeffery1922) for particle rotation. These authors showed that oblate spheroids preferentially align their symmetry axes normal to the wall, whereas prolate ones are preferentially parallel to the wall. The mean particle rotation was also reported to reduce when increasing the particle aspect ratio. Far from the wall, where the mean shear vanishes, this preferential alignment disappears and the behaviour is similar to that observed in homogeneous and isotropic turbulence (Voth Reference Voth2015). Kulick, Fessler & Eaton (Reference Kulick, Fessler and Eaton1994) studied small particles at higher concentrations, modelling feedback on the flow (two-way coupling), and reported that the fluid turbulence is attenuated by the addition of particles, while the turbulence anisotropy increases. This effect was reported to increase with the particle Stokes number, particle mass loading and distance to the wall. Paschkewitz et al. (Reference Paschkewitz, Dubief, Dimitropoulos, Shaqfeh and Moin2004) and Gillissen et al. (Reference Gillissen, Boersma, Mortensen and Andersson2008) showed drag reduction in suspensions of rigid fibres, similarly to what was reported for dilute polymer solutions (Ptasinski et al. Reference Ptasinski, Boersma, Nieuwstadt, Hulsen, Van den Brule and Hunt2003; Dubief et al. Reference Dubief, White, Terrapon, Shaqfeh, Moin and Lele2004).
Despite these previous efforts, the turbulent flow of finite-size non-spherical particles is still unexplored. This is therefore the object of the present study. In particular, we consider turbulent channel flow of finite-size oblate spheroids at volume fractions up to
$\unicode[STIX]{x1D719}=15\,\%$
. Aspect ratio (ratio of polar over equatorial radius)
${\mathcal{A}}{\mathcal{R}}=1/3$
is chosen for the particles to depart adequately from sphericity, where the effect of shape is more noticeable. We show that, unlike spherical particles, oblate particles cause drag reduction as the volume fraction
$\unicode[STIX]{x1D719}$
increases within the investigated range. We attribute the drag reduction to the absence of a particle–wall layer and to an attenuation of the near-wall turbulence, explained by the preferential orientation of particles and reduced rotation near the wall.
The paper is organized as follows. The governing equations and the flow geometry are introduced in § 2, followed by the results of the numerical simulations in § 3. The main conclusions and final remarks are drawn in § 4. Results for laminar flow are reported in appendix A for comparison.
2 Methodology
Several approaches for performing interface-resolved direct numerical simulations (DNS) of particle-laden flows have been proposed in recent years. Among these methods are force coupling (Lomholt & Maxey Reference Lomholt and Maxey2003), front tracking (Unverdi & Tryggvason Reference Unverdi and Tryggvason1992), Physalis (Zhang & Prosperetti Reference Zhang and Prosperetti2005; Sierakowski & Prosperetti Reference Sierakowski and Prosperetti2016), algorithms based on the lattice Boltzmann method for resolving the fluid phase (Ladd Reference Ladd1994a ,Reference Ladd b ) and the immersed boundary method (IBM) (Peskin Reference Peskin1972) as used here. Several algorithms for IBM have been proposed (Mittal & Iaccarino Reference Mittal and Iaccarino2005; Uhlmann Reference Uhlmann2005; Breugem Reference Breugem2012; Kempe & Fröhlich Reference Kempe and Fröhlich2012) since the original work by Peskin (Reference Peskin1972). The possibility of exploiting efficient computational algorithms for solving the Navier–Stokes equations on a Cartesian grid has made IBM a popular tool to investigate particle suspensions. The IBM algorithm proposed by Breugem (Reference Breugem2012) has been recently extended to ellipsoidal particles by Ardekani et al. (Reference Ardekani, Costa, Breugem and Brandt2016b ), using lubrication, friction and collision models for the short-range particle interactions. In this work we use the same numerical model to simulate dense suspensions of oblate spheroidal particles in turbulent plane channel flow.
2.1 Governing equations
The incompressible Navier–Stokes equations describe the flow field in the Eulerian phase:
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{f}\left(\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}\right)=-\unicode[STIX]{x1D735}p-\unicode[STIX]{x1D735}p_{e}+\unicode[STIX]{x1D707}_{f}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+\unicode[STIX]{x1D70C}_{f}\boldsymbol{f}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0. & \displaystyle\end{eqnarray}$$
Here
$\boldsymbol{u}$
is the fluid velocity,
$p_{e}$
is the contribution to the total pressure from a constant pressure gradient that drives the flow,
$p$
is the modified pressure (the total pressure minus
$p_{e}$
and the contribution from the hydrostatic pressure), and
$\unicode[STIX]{x1D70C}_{f}$
and
$\unicode[STIX]{x1D707}_{f}$
are the density and dynamic viscosity of the fluid. The extra term
$\boldsymbol{f}$
on the right-hand side of (2.1) is the IBM force field, active in the immediate vicinity of a particle surface to enforce no-slip and no-penetration boundary conditions.
The motion of rigid spheroidal particles is described by Newton–Euler Lagrangian equations:
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{p}V_{p}\frac{\text{d}\boldsymbol{U}_{p}}{\text{d}t}=\oint _{\unicode[STIX]{x2202}S_{p}}\unicode[STIX]{x1D749}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}A-V_{p}\unicode[STIX]{x1D735}p_{e}+(\unicode[STIX]{x1D70C}_{p}-\unicode[STIX]{x1D70C}_{f})V_{p}\boldsymbol{g}+\boldsymbol{F}_{c}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}(\unicode[STIX]{x1D644}_{p}\unicode[STIX]{x1D74E}_{p})}{\text{d}t}=\oint _{\unicode[STIX]{x2202}S_{p}}\boldsymbol{r}\times (\unicode[STIX]{x1D749}\boldsymbol{\cdot }\boldsymbol{n})\,\text{d}A+\boldsymbol{T}_{c}. & \displaystyle\end{eqnarray}$$
Here
$\boldsymbol{U}_{p}$
and
$\unicode[STIX]{x1D74E}_{p}$
are the particle translational and the angular velocity,
$\unicode[STIX]{x1D70C}_{p}$
,
$V_{p}$
and
$\unicode[STIX]{x1D644}_{p}$
are the mass density, volume and moment-of-inertia tensor of a spheroidal particle,
$\boldsymbol{r}$
indicates the position vector with respect to the centre of the spheroid and
$\boldsymbol{n}$
is the outward unit normal vector at the particle surface
$\unicode[STIX]{x2202}S_{p}$
where the stress tensor
$\unicode[STIX]{x1D749}=-p\unicode[STIX]{x1D644}+\unicode[STIX]{x1D707}_{f}(\unicode[STIX]{x1D735}\boldsymbol{u}+\unicode[STIX]{x1D735}\boldsymbol{u}^{\text{T}})$
acting on the surface of the particle is integrated. The force terms
$-\unicode[STIX]{x1D70C}_{f}V_{p}\boldsymbol{g}$
and
$V_{p}\unicode[STIX]{x1D735}p_{e}$
account for the hydrostatic pressure and a constant pressure gradient
$\unicode[STIX]{x1D735}p_{e}$
, with
$\boldsymbol{g}$
the gravitational acceleration. Finally,
$\boldsymbol{F}_{c}$
and
$\boldsymbol{T}_{c}$
are the force and torque resulting from particle–particle (particle–wall) collisions.
2.2 Numerical algorithm
The flow field is resolved on a uniform (
$\unicode[STIX]{x0394}x=\unicode[STIX]{x0394}y=\unicode[STIX]{x0394}z$
) staggered Cartesian grid, while particles are represented by a set of Lagrangian points, uniformly distributed on the surface of each particle. The number of Lagrangian grid points
$N_{L}$
on the surface of each particle is defined such that the Lagrangian grid volume
$\unicode[STIX]{x0394}V_{l}$
becomes equal to the volume of the Eulerian mesh
$\unicode[STIX]{x0394}x^{3}$
.
Taking into account the inertia of the fictitious fluid phase inside the particle volumes, Breugem (Reference Breugem2012) showed that (2.3) and (2.4) can be rewritten as
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}_{p}V_{p}\frac{\text{d}\boldsymbol{U}_{p}}{\text{d}t}\approx -\unicode[STIX]{x1D70C}_{f}\mathop{\sum }_{l=1}^{N_{L}}\boldsymbol{F}_{l}\unicode[STIX]{x0394}V_{l}+\unicode[STIX]{x1D70C}_{f}\frac{\text{d}}{\text{d}t}\left(\int _{V_{p}}\boldsymbol{u}\,\text{d}V\right)+(\unicode[STIX]{x1D70C}_{p}-\unicode[STIX]{x1D70C}_{f})V_{p}\boldsymbol{g}+\boldsymbol{F}_{c}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}(\unicode[STIX]{x1D644}_{p}\unicode[STIX]{x1D74E}_{p})}{\text{d}t}\approx -\unicode[STIX]{x1D70C}_{f}\mathop{\sum }_{l=1}^{N_{L}}\boldsymbol{r}_{l}\times \boldsymbol{F}_{l}\unicode[STIX]{x0394}V_{l}+\unicode[STIX]{x1D70C}_{f}\frac{\text{d}}{\text{d}t}\left(\int _{V_{p}}\boldsymbol{r}\times \boldsymbol{u}\,\text{d}V\right)+\boldsymbol{T}_{c}. & \displaystyle\end{eqnarray}$$
The point force
$F_{l}$
is calculated at each Lagrangian point using the difference between the particle surface velocity (
$\boldsymbol{U}_{p}+\unicode[STIX]{x1D74E}_{p}\times \boldsymbol{r}$
) and the interpolated first prediction velocity at the same point. The first prediction velocity is obtained by advancing (2.1) in time without considering the force field
$\boldsymbol{f}$
.
The forces
$\boldsymbol{F}_{l}$
integrate to the force field
$\boldsymbol{f}$
using the regularized Dirac delta function
$\unicode[STIX]{x1D6FF}_{d}$
of Roma, Peskin & Berger (Reference Roma, Peskin and Berger1999):
$$\begin{eqnarray}\boldsymbol{f}_{ijk}=\mathop{\sum }_{l=1}^{N_{L}}\boldsymbol{F}_{l}\unicode[STIX]{x1D6FF}_{d}(\boldsymbol{x}_{ijk}-\boldsymbol{X}_{l})\unicode[STIX]{x0394}V_{l},\end{eqnarray}$$
with
$\boldsymbol{x}_{ijk}$
and
$\boldsymbol{X}_{l}$
referring to an Eulerian and a Lagrangian grid cell. This smooth delta function essentially replaces the sharp interface with a thin porous shell of width
$3\unicode[STIX]{x0394}x$
; it preserves the total force and torque on the particle provided that the Eulerian grid is uniform. An iterative algorithm is employed to calculate the force field
$\boldsymbol{f}$
, allowing for a better estimate of no-slip and no-penetration boundary conditions (Breugem Reference Breugem2012). Equations (2.5) and (2.6) for the particle motion and (2.1) and (2.2) for the flow are integrated in time using an explicit low-storage Runge–Kutta method with the pressure correction scheme used in Breugem (Reference Breugem2012) to project the velocity field in the divergence-free space.
When the distance between particles (or a particle and a wall) is smaller than one Eulerian grid size, the lubrication force is underpredicted by the IBM. To compensate for this inaccuracy and to avoid computationally expensive grid refinements, a lubrication model based on the asymptotic analytical expression for the normal lubrication force between unequal spheres (Jeffrey Reference Jeffrey1982) is used. Here we approximate the two spheroidal particles with two spheres with the same mass and radius corresponding to the local curvature at the points of contact. Using these approximating spheres, a soft-sphere collision model with Coulomb friction takes over the interaction when the particles touch. The restitution coefficients used for normal and tangential collisions are 0.97 and 0.1, with Coulomb friction coefficient set to 0.15. More details about the models and validations can be found in Costa et al. (Reference Costa, Boersma, Westerweel and Breugem2015) and Ardekani et al. (Reference Ardekani, Costa, Breugem and Brandt2016b ).
2.3 Flow geometry
We study a pressure-driven plane channel flow in a computational domain of size
$L_{x}=6h$
,
$L_{y}=2h$
and
$L_{z}=3h$
in the streamwise, wall-normal and spanwise directions, where
$h$
is half the channel height. The bulk velocity
$U_{b}$
is fixed to guarantee a constant bulk Reynolds number
$Re_{b}=2hU_{b}/\unicode[STIX]{x1D708}=5600$
corresponding to a friction Reynolds number
$Re_{\unicode[STIX]{x1D70F}}=U_{\ast }h/\unicode[STIX]{x1D708}=180$
for the single-phase case with
$\unicode[STIX]{x1D708}$
, the kinematic viscosity of the fluid phase, and
$U_{\ast }=\sqrt{\unicode[STIX]{x1D70F}_{w}/\unicode[STIX]{x1D70C}_{f}}$
, the friction velocity, calculated with the shear stress
$\unicode[STIX]{x1D70F}_{w}$
at the wall. Periodic boundary conditions are imposed for both fluid and particles in the streamwise,
$x$
, and spanwise,
$z$
, directions, while the no-slip and no-penetration boundary conditions are employed at the walls. We note here that simulations in larger domains would be quite expensive computationally and have not been performed for the case of oblate particles. Nevertheless, the domain used in this study is larger than the minimal-unit channels adopted to identify the physical mechanisms underlying self-sustaining turbulence in Newtonian fluids (Hamilton, Kim & Waleffe Reference Hamilton, Kim and Waleffe1995) as well as polymer suspensions where drag reduction is also observed (Xi & Graham Reference Xi and Graham2010). In addition, as also reported in Picano et al. (Reference Picano, Breugem and Brandt2015), the presence of particles tends to disrupt long flow structures so that we believe the main conclusions from the simulations, focusing on drag reduction and particle dynamics, would not be different in longer domains.
We consider non-Brownian neutrally buoyant rigid spheroidal particles with aspect ratio
${\mathcal{A}}{\mathcal{R}}=1/3$
(ratio of polar over equatorial radius). The particle equivalent diameter
$D_{eq}$
, i.e. the diameter of a sphere with the same volume, is
$h/D_{eq}=9$
to compare with the results of Picano et al. (Reference Picano, Breugem and Brandt2015) for spheres with diameter
$D=D_{eq}$
. The corresponding maximum and minimum diameters of the oblate particle used here are
$D_{1}=h/6.24$
and
$D_{2}=h/18.72$
. The particle Reynolds number based on the local shear,
$Re_{p}\equiv \dot{\unicode[STIX]{x1D6FE}}D_{eq}^{2}/4\unicode[STIX]{x1D708}$
, ranges between approximately 95 close to the walls and 0 at the centreline of the channel, with
$\dot{\unicode[STIX]{x1D6FE}}$
approximated by the wall-normal gradient of the mean velocity.
We perform simulations at four different volume fractions
$\unicode[STIX]{x1D719}=5$
, 7.9, 10 and 15 %, corresponding to 7500 particles at
$\unicode[STIX]{x1D719}=15\,\%$
; we also consider the case of spheres at
$\unicode[STIX]{x1D719}=10\,\%$
and the unladen case for a direct comparison. We reproduced the results at
$\unicode[STIX]{x1D719}=10\,\%$
in Picano et al. (Reference Picano, Breugem and Brandt2015) with the collision model discussed in Ardekani et al. (Reference Ardekani, Costa, Breugem and Brandt2016b
), accounting also for friction between the particles (and the wall), and at higher resolution, the same used for the oblate particles. The results show a 4 % difference in the friction Reynolds number
$Re_{\unicode[STIX]{x1D70F}}$
, at twice the grid resolution.
The simulations are performed over a uniform Cartesian grid with the resolution of 32 grids per equivalent diameter
$D_{eq}$
, and
$N_{L}=3720$
and
$N_{L}=3219$
Lagrangian points on the surface of oblate and spherical particles, respectively. A summary of the simulated cases is given in table 1. The simulations start from the laminar Poiseuille flow with random distribution of the particles’ positions and orientations. The noise introduced by the presence of the particles rapidly triggers the transition to a fully turbulent state, after which the statistics are collected for approximately 16 large-eddy turnover times
$h/U_{\ast }$
.
Instantaneous snapshots of streamwise velocity
$u$
in the presence of particles on orthogonal planes
$xz$
,
$yz$
and
$xy$
for (a) oblates at
$\unicode[STIX]{x1D719}=5\,\%$
, (b) oblates at
$\unicode[STIX]{x1D719}=10\,\%$
, (c) oblates at
$\unicode[STIX]{x1D719}=15\,\%$
and (d) spheres at
$\unicode[STIX]{x1D719}=10\,\%$
. For clarity, just a fraction of particles lying within the
$xz$
and
$yz$
planes are displayed.
Summary of the different simulation cases:
$N_{p}$
indicates the number of particles with equivalent diameter
$D_{eq}=h/9$
;
$N_{x}$
,
$N_{y}$
and
$N_{z}$
are the number of grid cells in each direction; and
$\unicode[STIX]{x0394}x^{+}$
and
${D_{eq}}^{+}$
are the Eulerian grid size and the particle equivalent diameter, given in viscous wall units for each simulation.
3 Results
We first display snapshots of the fluid flow and particles; see figure 1, where the instantaneous streamwise velocity
$u$
is depicted on different horizontal and vertical planes for the cases with oblates at
$\unicode[STIX]{x1D719}=5$
,
$10$
and
$15\,\%$
and for spheres at
$\unicode[STIX]{x1D719}=10\,\%$
. For clarity, just a fraction of the particles (those lying within the visualized
$xz$
and
$yz$
planes) are displayed. Streamwise low-speed streaks, characteristic of wall-bounded turbulence, can be observed near the wall for all cases, being however more noisy in the flow laden with spheres.
3.1 Fluid-phase statistics: drag reduction
For the statistics presented hereafter, the fluid velocity is considered only when outside the particles, whereas values pertaining to rigid-body motion of the solid phase are used for the particle statistics. To distinguish between phases, we use a phase indicator function that determines whether a grid point is located inside the fluid or solid phase. The mean fluid velocity profiles, scaled in outer units
$U_{f}/U_{b}$
and inner units
$U_{f}^{+}\equiv U_{f}/U_{\ast }$
, are depicted in figure 2 versus
$y/h$
and
$y^{+}\equiv yU_{\ast }/\unicode[STIX]{x1D708}$
for the different cases. The results in figure 2(a) show that the mean velocity at the core of the channel increases when increasing the solid volume fraction (see inset), whereas it decreases around
$y/h\approx 0.2$
, that is at the beginning of the log layer (
$y^{+}\approx 35$
) for the studied Reynolds number. As the simulations assume constant bulk velocity
$U_{b}$
, this is a first indication of reduced drag. The mean fluid velocity profile for spheres at
$\unicode[STIX]{x1D719}=10\,\%$
deviates from the unladen case more than for oblates at the same volume fraction. Displaying data in inner scaling (see figure 2
b), we identify a region (
$50<y^{+}<150$
) where the mean profiles follow a log law of the type
$$\begin{eqnarray}U_{f}^{+}=\frac{1}{\unicode[STIX]{x1D705}}\ln (y^{+})+\unicode[STIX]{x1D6FD}.\end{eqnarray}$$
The values of the effective von Kármán and additive constants,
$\unicode[STIX]{x1D705}$
and
$\unicode[STIX]{x1D6FD}$
, that best fit our results are reported in table 2 together with the friction Reynolds number
$Re_{\unicode[STIX]{x1D70F}}$
. Note that a reduced
$\unicode[STIX]{x1D705}$
is a sign of drag reduction, while a reduced
$\unicode[STIX]{x1D6FD}$
indicates the opposite (Virk Reference Virk1975). The data in table 2 show that
$\unicode[STIX]{x1D705}$
and
$\unicode[STIX]{x1D6FD}$
decrease with increasing volume fraction
$\unicode[STIX]{x1D719}$
for oblate particles. The largest reduction of
$\unicode[STIX]{x1D6FD}$
is however found for spheres, while the value of
$\unicode[STIX]{x1D705}$
is close to that for oblates at
$\unicode[STIX]{x1D719}=10\,\%$
. The combination of the two effects results in drag reduction for oblate particles and a drag enhancement for the spheres as confirmed by the friction Reynolds number
$Re_{\unicode[STIX]{x1D70F}}$
, shown in figure 3(a) versus
$\unicode[STIX]{x1D719}$
for oblate and spherical particles. For oblate particles, the effect of the solid phase is small for
$\unicode[STIX]{x1D719}\lesssim 5\,\%$
; the total drag clearly decreases below the value of single-phase flow for larger volume fractions:
$Re_{\unicode[STIX]{x1D70F}}=172.5$
for
$\unicode[STIX]{x1D719}=15\,\%$
, corresponding to a drag reduction of 8.2 % with respect to the single-phase case (
$(\unicode[STIX]{x1D70F}_{w}|_{\unicode[STIX]{x1D719}=15\,\%}-\unicode[STIX]{x1D70F}_{w}|_{\unicode[STIX]{x1D719}=0\,\%})/\unicode[STIX]{x1D70F}_{w}|_{\unicode[STIX]{x1D719}=0\,\%}$
). The drag reduction percentage is depicted versus the volume fraction in figure 3(b). The data in Picano et al. (Reference Picano, Breugem and Brandt2015) for perfect spheres, conversely, shows a clear increase of the friction Reynolds number when increasing the volume fraction.
Mean fluid velocity profiles
$U_{f}$
in the streamwise direction for the different cases under investigation. The data are scaled with (a) outer units,
$U_{f}/U_{b}$
versus
$y/h$
, and (b) inner units,
$U_{f}^{+}\equiv U_{f}/U_{\ast }$
versus
$y^{+}\equiv yU_{\ast }/\unicode[STIX]{x1D708}$
, depicted in semi-logarithmic scale. The inset in (a) displays the increase of the mean velocity at the core of the channel.
(a) Friction Reynolds number
$Re_{\unicode[STIX]{x1D70F}}$
for oblate particles and spherical particles from Picano et al. (Reference Picano, Breugem and Brandt2015) and (b) the drag reduction percentage for oblate particles versus volume fraction
$\unicode[STIX]{x1D719}$
. The case denoted ‘Oblates (NR)’ (no rotation) will be addressed later in this section.
The friction Reynolds number
$Re_{\unicode[STIX]{x1D70F}}$
from all simulations at fixed bulk Reynolds number
$Re_{b}=5600$
and the von Kármán constant
$\unicode[STIX]{x1D705}$
and the additive constant
$\unicode[STIX]{x1D6FD}$
of the logarithmic law, calculated in the range of
$50<y^{+}<150$
.
Drag reduction with small (smaller than the Kolmogorov length scale) rigid fibre additives has been reported in the literature (e.g. Paschkewitz et al.
Reference Paschkewitz, Dubief, Dimitropoulos, Shaqfeh and Moin2004; Gillissen et al.
Reference Gillissen, Boersma, Mortensen and Andersson2008). These authors associate the drag reduction with the attenuation of the turbulence and the increase of its anisotropy, which results in higher streamwise velocity fluctuations and lower spanwise and wall-normal velocity fluctuations (when scaled in inner units) with respect to the single-phase flow, similarly to what was observed for polymer additives (De Angelis, Casciola & Piva Reference De Angelis, Casciola and Piva2002; Nowbahar et al.
Reference Nowbahar, Sardina, Picano and Brandt2013). De Angelis et al. (Reference De Angelis, Casciola and Piva2002) explained this by revealing that, in most of the field, polymers are extracting energy from the turbulence, resulting in a reduction of spanwise and wall-normal velocity fluctuations. In the streamwise direction, however, an increase in the velocity fluctuations is caused by larger streamwise vortices located further from the wall. The root-mean-square (r.m.s.) of the fluid velocity fluctuations and the Reynolds shear stresses for suspensions of finite-size oblates are depicted in figure 4 for the different cases studied here. The data reveal that the peak of the turbulent velocity fluctuations is reduced in the presence of oblate particles with respect to the single-phase flow for all three components. In the cross-stream directions, the fluctuations decrease with the volume fraction of oblate particles, while an increase is observed for spheres at
$\unicode[STIX]{x1D719}=10\,\%$
. The peak in the streamwise velocity fluctuations displays a reduction for both spherical and oblate particles with respect to the single-phase flow, while a slight increase of the fluctuations is observable in the regions far from the wall. This slight increase, present for both spherical and oblate particles, can be related to the movement of particles in the wall-normal direction, where they accelerate or decelerate the flow based on their streamwise velocity. A clear turbulence attenuation is shown by the Reynolds shear stress profiles of oblate particles. The turbulence attenuation is weak at
$\unicode[STIX]{x1D719}=5\,\%$
and becomes more pronounced at higher volume fractions; a considerable reduction in the Reynolds shear stress is observed for the case at
$\unicode[STIX]{x1D719}=15\,\%$
. On the contrary, a significant increase in the turbulence activity is observed for spheres at
$\unicode[STIX]{x1D719}=10\,\%$
. Interestingly,
$v_{f}^{\prime }$
and
$w_{f}^{\prime }$
change significantly with respect to the single-phase flow in the close vicinity of the wall, which is peculiar to spherical particles and missing in the case of oblates. This indicates that the effect of particles on the turbulence very close to the wall is considerably higher for spherical particles than for oblates due to the absence of a particle layer close to the wall, as documented later. As discussed in Costa et al. (Reference Costa, Picano, Brandt and Breugem2016), accounting for the particle dynamics in this region is critical for an accurate prediction of the overall drag.
Root-mean-square velocity fluctuations and Reynolds shear stress for the fluid phase, scaled in outer units: (a) streamwise
${u_{f}^{\prime }}_{\!rms}$
; (b) wall-normal
${v_{f}^{\prime }}_{\!rms}$
; (c) spanwise
${w_{f}^{\prime }}_{\!rms}$
; and (d) Reynolds shear stress
$\langle u_{f}^{\prime }v_{f}^{\prime }\rangle$
.
Correlations of the velocity fluctuations versus the spanwise separation, normalized by the major radius of the oblate particles
$z/R_{max}$
, for the different cases under consideration: (a) streamwise–streamwise component
$R_{uu}$
at
$y^{+}=20$
; (b)
$R_{uu}$
at
$y^{+}=40$
; (c) wall-normal–wall-normal component
$R_{vv}$
at
$y^{+}=20$
; and (d)
$R_{vv}$
at
$y^{+}=40$
.
Finally, we examine how the presence of the oblate particles affects the turbulent structures near the wall by computing the two-point spatial correlation of the velocity field. The autocorrelations of the streamwise and wall-normal velocity,
$R_{uu}$
and
$R_{vv}$
, are depicted in figure 5 versus the spanwise separation, normalized by the major radius of the oblate particles. Figure 5(a,b) shows
$R_{uu}$
at two different distances from the wall,
$y^{+}=20$
and
$y^{+}=40$
, while figure 5(c,d) reports
$R_{vv}$
at the same wall distances. The data reveal the characteristics of typical turbulence structures near the wall, i.e. quasi-streamwise vortices and low-speed streaks (Kim, Moin & Moser Reference Kim, Moin and Moser1987; Pope Reference Pope2000; Brandt Reference Brandt2014). The results at
$y^{+}=20$
show an increase in the minimum location of
$R_{uu}$
, i.e. widening of the near-wall streaks, whose size becomes independent of the volume fraction, and appears to be more correlated to the particle size. This is in contrast with the results for spherical particles reported in Picano et al. (Reference Picano, Breugem and Brandt2015), where the spacing increases monotonically with the volume fraction. We also observe a decrease of the minimum of the correlation, i.e. less pronounced streaks, as the volume fraction increases, while further away from the wall (
$y^{+}=40$
) this effect disappears; the peak is still distinct, but the streamwise velocity streaks are considerably wider in the particle-laden flows. The results for
$R_{uu}$
are consistent with the turbulence attenuation reported above for oblate particles as the width and the spacing of the streamwise velocity streaks increase in the drag-reducing turbulent flows (Stone, Waleffe & Graham Reference Stone, Waleffe and Graham2002), corresponding to an increase in the extent of the buffer layer. The wall-normal autocorrelations
$R_{vv}$
indicate a progressive smoothing of the local minima with the volume fraction in the cases of oblate particles at both wall-normal locations under consideration, indicating a more random flow in terms of coherent turbulence structures.
3.2 Effective viscosity and stress budget
To better understand the effect of particles on the fluid turbulence, one may consider the idealized case of a suspension with effective viscosity
$\unicode[STIX]{x1D708}_{e}$
. The effective viscosity of a particle suspension is always higher than the viscosity for the single-phase flow: the ratio between the effective suspension viscosity and the fluid viscosity
$\unicode[STIX]{x1D708}_{r}=\unicode[STIX]{x1D708}_{e}/\unicode[STIX]{x1D708}$
is typically estimated via empirical relations at volume fractions larger than a few per cent, an example being the Eilers fit (Stickel & Powell Reference Stickel and Powell2005),
$\unicode[STIX]{x1D708}_{r}=[1+1.25\unicode[STIX]{x1D719}/(1-\unicode[STIX]{x1D719}/0.6)]^{2}$
. The effective viscosity estimated from simulations of spherical and oblate particles in laminar flow at
$Re_{b}=1000$
and
$\unicode[STIX]{x1D719}=10\,\%$
(see appendix A) is in good agreement with the mentioned fit (2 % underpredicted for oblates and 4 % overpredicted for spheres), and this relation is therefore used to estimate the effective viscosity used for the analysis of suspensions in turbulent flow.
The results in Picano et al. (Reference Picano, Breugem and Brandt2015) show that spherical particles increase the turbulence activity in the flow up to
$\unicode[STIX]{x1D719}=10\,\%$
, an effect that overcomes the increase of the effective viscosity of the suspension. In fact, the increase in the effective viscosity can result in turbulence attenuation; however, the cases in Picano et al. (Reference Picano, Breugem and Brandt2015) show an enhancement in the presence of spherical particles. For oblate particles, instead, we observe a clear turbulence attenuation and an associated overall drag reduction. This raises the question of whether the increased viscosity of the suspension is the only reason for the turbulence attenuation or whether the specific particle shape also contributes to dampen the turbulence. To answer this question, we perform simulations of single-phase flow at bulk Reynolds number equal to the effective bulk Reynolds number
$Re_{b}^{e}=2hU_{b}/\unicode[STIX]{x1D708}_{e}$
, calculated for the corresponding particle volume fraction.
(a) Reynolds stress of the combined phase, scaled in outer units for oblate particles at
$\unicode[STIX]{x1D719}=10\,\%$
, compared to a single-phase flow simulation with
$\unicode[STIX]{x1D708}=\unicode[STIX]{x1D708}_{e}|_{\unicode[STIX]{x1D719}=10\,\%}$
. The solid lines show the linear fitting of the slope of the profiles at the centreline (
$y/h=1$
). (b) The effective turbulent friction Reynolds number
$Re_{T}^{e}$
and the effective friction Reynolds number
$Re_{\unicode[STIX]{x1D70F}}^{e}$
versus the volume fraction
$\unicode[STIX]{x1D719}$
. These values are compared to the effective turbulent friction Reynolds number
$Re_{T}^{e_{0}}$
and the effective friction Reynolds number
$Re_{\unicode[STIX]{x1D70F}}^{e_{0}}$
obtained from single-phase simulations of a fluid with viscosity equal to the effective viscosity of the suspension.
The ratio between effective suspension viscosity and the fluid viscosity
$\unicode[STIX]{x1D708}_{r}$
, the effective friction Reynolds number
$Re_{\unicode[STIX]{x1D70F}}^{e}\equiv hU_{\ast }/\unicode[STIX]{x1D708}_{e}$
, the effective turbulent friction Reynolds number
$Re_{T}^{e}\equiv hU_{\ast }^{\text{T}}/\unicode[STIX]{x1D708}_{e}$
, the effective suspension bulk Reynolds number
$Re_{b}^{e}\equiv 2hU_{b}/\unicode[STIX]{x1D708}_{e}$
, and the effective turbulent friction Reynolds number and effective friction Reynolds number
$Re_{T}^{e_{0}}$
and
$Re_{\unicode[STIX]{x1D70F}}^{e_{0}}$
obtained from single-phase simulations by accounting only for the effective suspension viscosity, for all cases with oblate particles.
Figure 6(a) shows the Reynolds stress of the combined phase
$\langle u_{c}^{\prime }v_{c}^{\prime }\rangle =\unicode[STIX]{x1D6F7}\langle u_{p}^{\prime }v_{p}^{\prime }\rangle +(1-\unicode[STIX]{x1D6F7})\langle u_{f}^{\prime }v_{f}^{\prime }\rangle$
for the case with oblate particles at
$\unicode[STIX]{x1D719}=10\,\%$
, compared to a single-phase flow simulation with
$\unicode[STIX]{x1D708}=\unicode[STIX]{x1D708}_{e}|_{\unicode[STIX]{x1D719}=10\,\%}$
. The data reveal that, in the presence of oblate particles, the turbulence activity is even lower than that predicted by the simulation at
$Re_{b}^{e}$
, indicating that not only the effective viscosity is responsible for the turbulence attenuation shown above.
As in Picano et al. (Reference Picano, Breugem and Brandt2015), to quantify the level of turbulence activity, we define the turbulent friction velocity
$U_{\ast }^{\text{T}}=\sqrt{\text{d}\langle u_{c}^{\prime }v_{c}^{\prime }\rangle /\text{d}(y/h)}|_{y/h=1}$
, calculated as the square root of the wall-normal derivative of the Reynolds stress profile at the centreline of the channel. The slope of the Reynolds stress profile at the centreline is shown in figure 6(a). Note that the turbulent friction velocity approximates the wall friction velocity for single-phase flows at high bulk Reynolds number, the error scaling with
$1/Re_{b}$
(Pope Reference Pope2000). This estimate accounts only for the effect of turbulence and not for the viscous or the particle-induced stress and therefore helps to quantify the role of the turbulence activity on the overall drag.
The effective turbulent friction Reynolds number, defined as
$Re_{T}^{e}=U_{\ast }^{\text{T}}h/\unicode[STIX]{x1D708}_{e}$
, is reported in table 3 for all cases under consideration together with the turbulent friction Reynolds number obtained from the single-phase simulations
$Re_{T}^{e_{0}}$
(with
$\unicode[STIX]{x1D708}=\unicode[STIX]{x1D708}_{e}$
). In addition to these two, we also depict in figure 6(b) the friction Reynolds number based on the effective viscosity,
$Re_{\unicode[STIX]{x1D70F}}^{e}=Re_{\unicode[STIX]{x1D70F}}/\unicode[STIX]{x1D708}_{r}$
, and
$Re_{\unicode[STIX]{x1D70F}}^{e_{0}}$
, the friction Reynolds number obtained from the single-phase simulations. These are shown for increasing volume fraction
$\unicode[STIX]{x1D719}$
and oblate particles. We conclude from figure 6(b) that the oblate particles reduce the turbulence activity to lower values than those obtained by accounting only for the effective suspension viscosity, resulting in an overall drag reduction. The observed drag reduction is therefore related to the specific dynamics of the oblate particles and their interactions with the turbulent velocity field.
Momentum budget, normalized with
$\unicode[STIX]{x1D70C}_{f}U_{\ast }^{2}$
, for oblate particles at (a)
$\unicode[STIX]{x1D719}=5\,\%$
, (b)
$\unicode[STIX]{x1D719}=10\,\%$
and (c)
$\unicode[STIX]{x1D719}=15\,\%$
, and for spheres at (d)
$\unicode[STIX]{x1D719}=10\,\%$
. Here,
$\unicode[STIX]{x1D70F}$
is the total stress,
$\unicode[STIX]{x1D70F}_{V}$
denotes viscous stress,
$\unicode[STIX]{x1D70F}_{T}=-\langle u_{c}^{\prime }v_{c}^{\prime }\rangle$
the turbulent Reynolds shear stress of the combined phase and
$\unicode[STIX]{x1D70F}_{P}$
the particle-induced stress. The particle Reynolds shear stress
$\unicode[STIX]{x1D70F}_{T_{p}}=-\langle u_{p}^{\prime }v_{p}^{\prime }\rangle$
is also depicted with dots.
The results presented so far indicate that the turbulence attenuation and the absence of a particle layer close to the wall (documented in the next section) are responsible for the drag reduction in the flow laden with oblate particles. Picano et al. (Reference Picano, Breugem and Brandt2015) showed that, for spherical particles at
$\unicode[STIX]{x1D719}=20\,\%$
, the turbulence activity reduces while the total drag still increases. This is attributed to the increase of the particle-induced stress at high volume fractions. To better understand the effect of particle-induced stress for the cases with oblate particles, we perform here a stress budget analysis similar to that in Lashgari et al. (Reference Lashgari, Picano, Breugem and Brandt2014) and Picano et al. (Reference Picano, Breugem and Brandt2015).
Based on the formulation proposed in Zhang & Prosperetti (Reference Zhang and Prosperetti2010), we can write the mean momentum balance in the channel as
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D70F}}{\unicode[STIX]{x1D70C}_{f}}=U_{\ast }^{2}\left(1-\frac{y}{h}\right)=\unicode[STIX]{x1D708}(1-\unicode[STIX]{x1D6F7})\frac{\text{d}U_{f}}{\text{d}y}-[\unicode[STIX]{x1D6F7}\langle u_{p}^{\prime }v_{p}^{\prime }\rangle +(1-\unicode[STIX]{x1D6F7})\langle u_{f}^{\prime }v_{f}^{\prime }\rangle ]+\frac{\unicode[STIX]{x1D6F7}}{\unicode[STIX]{x1D70C}}\langle \unicode[STIX]{x1D70E}_{xy}^{p}\rangle ,\end{eqnarray}$$
where
$\unicode[STIX]{x1D70F}$
is the total stress; see the Appendix in Picano et al. (Reference Picano, Breugem and Brandt2015) for a derivation. The first term on the right-hand side of the budget above is the viscous shear stress
$\unicode[STIX]{x1D70F}_{V}$
, the second and third term are together the turbulent Reynolds shear stress of the combined phase
$\unicode[STIX]{x1D70F}_{T}=-\langle u_{c}^{\prime }v_{c}^{\prime }\rangle$
, and the fourth term indicates the particle-induced stress
$\unicode[STIX]{x1D70F}_{P}$
. The momentum transfer pertaining to each term, normalized by
$\unicode[STIX]{x1D70C}_{f}U_{\ast }^{2}$
, is depicted versus
$y/h$
in figure 7 for oblate particles at
$\unicode[STIX]{x1D719}=5$
, 10 and 15 % and for spheres at
$\unicode[STIX]{x1D719}=10\,\%$
. The contribution due to the particle-induced stress is always lower than that from the Reynolds stress for oblate particles. As the volume fraction increases, the relative momentum transfer due to the Reynolds stress decreases, yet the particle-induced stress does not increase enough to compensate for the reduction in turbulence activity. The particle-induced stress is considerably higher near the wall for spheres than for oblate particles, as expected by the absence of the particle layer close to the wall in the cases of oblates. The momentum budget analysis performed here shows that, unlike the case of spheres, the effect of the particle-induced stress on the total drag is small for oblate particles and the Reynolds shear stress is the main factor determining the overall drag. This confirms that the specific dynamics of the oblate particles is closely related to the turbulence attenuation and therefore drag reduction.
3.3 Particle dynamics
The mean local volume fraction
$\unicode[STIX]{x1D6F7}(y)$
, normalized by the bulk volume fraction
$\unicode[STIX]{x1D719}$
, is depicted in figure 8(a). Spherical particles display a local maximum at a distance slightly larger than one particle radius from the wall. Picano et al. (Reference Picano, Breugem and Brandt2015) attribute this local maximum to the formation of a particle layer at the wall, due to the wall–particle interactions that stabilize the particle position. Costa et al. (Reference Costa, Picano, Brandt and Breugem2016) explain that the presence of this particle layer always results in drag increase, which is therefore higher than what can be predicted by accounting only for the effective suspension viscosity. Interestingly, the particle layer at the wall is not present in the flow with oblate particles. The local volume fraction,
$\unicode[STIX]{x1D6F7}$
, is considerably lower in the region close to the wall (
$y/h<0.1$
) for oblates than for spherical particles. A minor migration towards the channel centre is observed for oblate particles in turbulent flow and this effect is more pronounced when increasing the volume fraction. Interestingly, the oblate distribution is uniform in laminar flow, while spherical particles tend to migrate towards the channel centre; see figure 15 in appendix A.
Figure 8(b) shows the mean particle velocity profiles. The difference between the particle and the fluid mean velocity, the local slip, is depicted in the inset. As the velocity of the particles is not zero at the wall, significant slip velocity is found close to the wall. Interestingly, a local minimum of the velocity difference is observed for spheres at the location of the particle layer close to the wall. Spherical particles, trapped in the particle layer, experience a smaller slip velocity due to the higher local volume fraction at this location, thus creating the mentioned local minimum in the mean velocity difference. This particle layer disappears for the case with oblates and the difference in the mean velocity decays monotonically with decreasing wall-normal distance.
Profiles of solid-phase averaged data versus
$y/h$
: (a) mean local volume fraction
$\unicode[STIX]{x1D6F7}(y)$
, normalized by total volume fraction
$\unicode[STIX]{x1D719}$
, and (b) mean particle velocity profiles, normalized by the fluid bulk velocity
$U_{b}$
. The inset in (b) shows the difference between the particle and the fluid mean velocity.
Contours of the mean forces and torques acting on the particles as a function of distance to the wall,
$y/h$
, and the streamwise component of the particle orientation vector,
$\hat{\boldsymbol{o}}_{p}^{x}$
. (a) Drag force
$F_{D}$
and (b) lift force
$F_{L}$
in units of
$\unicode[STIX]{x1D70C}_{f}U_{b}^{2}D_{eq}^{2}$
. (c) Spanwise torque
$T_{z}$
in units of
$\unicode[STIX]{x1D70C}_{f}U_{b}^{2}D_{eq}^{3}$
. (d) Schematic of the action of the forces and torques for the cases of positive and negative
$\hat{\boldsymbol{o}}_{p}^{x}$
.
Intensity of the solid-phase velocity fluctuation components, scaled in outer units (bulk velocity
$U_{b}$
), for the different cases under consideration: (a) streamwise
${u_{f}^{\prime }}_{\!rms}$
, (b) wall-normal
${v_{f}^{\prime }}_{\!rms}$
and (c) spanwise
${w_{f}^{\prime }}_{\!rms}$
components. (d) The ratio between the turbulent kinetic energy of the fluid and that of the solid phase
$K_{f}/K_{p}$
.
The reasons behind the absence of the particle layer for oblates are investigated here using a force and torque analysis in the case at
$\unicode[STIX]{x1D719}=15\,\%$
. We compute the mean drag force
$F_{D}$
, lift force
$F_{L}$
and spanwise torque
$T_{z}$
acting on the particles as a function of distance to the wall,
$y/h$
, and of the streamwise component of the particle orientation vector
$\hat{\boldsymbol{o}}_{p}^{x}$
. The orientation vector
$\hat{\boldsymbol{o}}_{p}$
is defined as the unit vector parallel to the particle symmetry axis and pointing towards the channel centre. Contours of the normalized
$F_{D}$
,
$F_{L}$
and
$T_{z}$
are depicted in figure 9(a–c). The results show that, when particles are sufficiently close to the wall (high mean velocity gradients), the forces and the torques acting on the particle change according to their orientation as sketched in figure 9(d). Particles, on average, have a negative spanwise angular velocity due to the mean flow gradient, and the analysis presented here, interestingly, indicates that, when these have a positive
$\hat{\boldsymbol{o}}_{p}^{x}$
, they are lifted by the flow towards the channel centre whereas the opposite is true when they have a negative
$\hat{\boldsymbol{o}}_{p}^{x}$
. Note, however, that the magnitude of the upward lift is significantly higher than the downward force. This difference in the magnitude can explain the absence of a particle layer for oblates. It can also be concluded from figure 9(d) that particles near the wall with a negative
$\hat{\boldsymbol{o}}_{p}^{x}$
tend to align parallel to the wall and accelerate towards it; particles with positive
$\hat{\boldsymbol{o}}_{p}^{x}$
, conversely, accelerate towards the channel centre. This indicates that oblate particles are most likely parallel to the wall in its vicinity. This stable configuration and its importance for drag reduction are addressed later in this section.
To gain further insight into the turbulence attenuation in the presence of oblates, we investigate the particle collective behaviour. Figure 10 shows the r.m.s. of the solid-phase velocity fluctuations. The data clearly reveal that all three components of the particle velocity fluctuations are considerably smaller than in the case of spherical particles, except for a very small region close to the wall, where the rare collisions with the wall have a large impact on the statistics. It should also be noted that, contrary to the case of spheres, oblate particles cannot roll on the wall and any collision with the wall can create a large particle velocity fluctuation. A maximum of the wall-normal particle velocity fluctuations
$v_{p}^{\prime }$
very close to the wall and a local minimum at
$y/h\approx 0.1$
are observed for oblate particles, whereas a local maximum and a weak local minimum at
$y/h\approx 0.15$
can be seen for spheres. The ratio between the turbulent kinetic energy of the fluid and of the solid phase
$K_{f}/K_{p}=({u_{f}^{\prime }}^{2}+{v_{f}^{\prime }}^{2}+{w_{f}^{\prime }}^{2})/({u_{p}^{\prime }}^{2}+{v_{p}^{\prime }}^{2}+{w_{p}^{\prime }}^{2})$
, depicted in figure 10(d), shows that particle velocities tend to fluctuate less than the fluid one at the same position except for the region close to the wall where velocities fluctuate more due to the absence of a no-slip condition. This tendency is considerably higher for oblate particles, as shown by the peak of the velocity fluctuations close to the wall,
$y\approx 0.1$
.
(a) Cosine of the mean particle inclination angle
$\unicode[STIX]{x1D703}$
, measured with respect to the wall, versus
$y/h$
. Mean particle absolute value of angular velocity, given in outer units for: (b) spanwise direction
$\overline{|\unicode[STIX]{x1D6FA}_{z}|}$
, (c) streamwise direction
$\overline{|\unicode[STIX]{x1D6FA}_{x}|}$
and (d) wall-normal direction
$\overline{|\unicode[STIX]{x1D6FA}_{y}|}$
.
The cosine of the mean particle inclination angle,
$\unicode[STIX]{x1D703}$
, measured with respect to the wall, is reported in figure 11(a) versus the distance from the wall. Note that values of
$\cos \unicode[STIX]{x1D703}$
close to 1 indicate that the particles tend to be aligned with their semi-minor (symmetry) axis normal to the wall, whereas the oblates sit with the major axis normal to the wall if
$\unicode[STIX]{x1D703}\approx 0$
. From figure 11(a), we note the clear tendency of oblate particles to be, on average, parallel to the wall; far from the wall they tend to be more isotropic, even though they still show some preferential orientation. This preferential orientation far from the wall is more pronounced in the laminar regime when the fluid velocity fluctuations are weak (see appendix A). For spheres, the poles are fixed arbitrarily to show their uniform orientation.
The mean absolute value of the particle angular velocities is reported in figure 11. As seen in figure 11(b), displaying the rotation rate in the spanwise direction,
$\overline{|\unicode[STIX]{x1D6FA}_{z}|}$
, oblate particles have significantly lower rotation rates than spheres close to the wall. The mean particle angular velocity in the streamwise and wall-normal directions (see figure 11
c,d) are lower in the case of oblate particles, except for the near wall-region where the wall-normal angular velocity is larger than that of spheres.
As depicted in figure 11(a), oblate particles are, on average, parallel to the wall; to fully characterize the particle relative orientation we therefore consider the properties of their order-parameter tensor (Prost Reference Prost1995),
$$\begin{eqnarray}\unicode[STIX]{x1D64C}=\langle \hat{\boldsymbol{o}}\otimes \hat{\boldsymbol{o}}-\unicode[STIX]{x1D644}/3\rangle ,\end{eqnarray}$$
where
$\hat{\boldsymbol{o}}$
is the unit vector associated with the particle symmetry axis and
$\unicode[STIX]{x1D644}$
the identity tensor. The eigenvalues (
$\unicode[STIX]{x1D706}_{i}$
) of the tensor
$\unicode[STIX]{x1D64C}$
are equal to zero in the case of a fully isotropic particle orientation, while the tendency to align in a certain direction (nematic order) is reflected in non-vanishing real eigenvalues. The eigenvalues (
$\unicode[STIX]{x1D706}_{1}>\unicode[STIX]{x1D706}_{2}>\unicode[STIX]{x1D706}_{3}$
) of the symmetric tensor
$\unicode[STIX]{x1D64C}$
can be expressed in the following form:
$$\begin{eqnarray}\unicode[STIX]{x1D706}_{1}={\textstyle \frac{2}{3}}\unicode[STIX]{x1D706},\quad \unicode[STIX]{x1D706}_{2}=-{\textstyle \frac{1}{3}}(\unicode[STIX]{x1D706}+\unicode[STIX]{x1D701}),\quad \unicode[STIX]{x1D706}_{3}=-{\textstyle \frac{1}{3}}(\unicode[STIX]{x1D706}-\unicode[STIX]{x1D701}),\end{eqnarray}$$
with
$\unicode[STIX]{x1D706}$
the scalar nematic order parameter and
$\unicode[STIX]{x1D701}$
the biaxial index. It can be shown (Prost Reference Prost1995) that
$\unicode[STIX]{x1D706}$
varies between
$-0.5<\unicode[STIX]{x1D706}<1$
, where a positive
$\unicode[STIX]{x1D706}$
indicates that particles tend to be oriented in one direction, while a negative
$\unicode[STIX]{x1D706}$
shows that the orientation vector preferentially lies in one plane. The parameter
$\unicode[STIX]{x1D701}$
shows the biaxiality of the order-parameter tensor:
$\unicode[STIX]{x1D701}\approx 0$
indicates that there is only one preferential direction associated with the eigenvector of largest eigenvalue (
$\unicode[STIX]{x1D706}_{1}$
), and a non-zero value of
$\unicode[STIX]{x1D701}$
creates a nematic phase that is considered biaxial, meaning that there exist two preferential directions. Here we divide the half-channel height
$h$
into 20 regions and compute the order-parameter tensor for each of these regions. The values of
$\unicode[STIX]{x1D706}$
and
$\unicode[STIX]{x1D701}$
obtained in each slab are displayed in figure 12(a,b) versus the channel height for all the investigated cases. The results reveal, as expected from figure 11(a), a large value of
$\unicode[STIX]{x1D706}$
close to the wall (
$\unicode[STIX]{x1D706}\approx 1$
), which decays at the centre of the channel. Thus, particles are preferentially aligned parallel to the wall in its vicinity. From the values of
$\unicode[STIX]{x1D701}$
throughout the channel, we find a region with larger
$\unicode[STIX]{x1D701}$
(
$0.2<y/h<0.4$
), where the tensor
$\unicode[STIX]{x1D64C}$
can be considered biaxial. This means that particles tend to be preferentially oriented also in a second direction, which is found to be in the spanwise direction (particles are rolling with their symmetry axis, oriented in the spanwise direction). From a phenomenological point of view, it appears that particles close to the wall leave the area, rolling in the streamwise direction; therefore particles preferentially align with their symmetry axis in the spanwise direction. This information can be obtained from the eigenvector associated with
$\unicode[STIX]{x1D706}_{2}$
.
(a) The nematic order parameter
$\unicode[STIX]{x1D706}$
and (b) the biaxial parameter
$\unicode[STIX]{x1D701}$
versus
$y/h$
for the simulated cases. Results for spheres are also depicted to indicate a suspension with fully isotropic particle orientation.
Orientational correlation function
$OCF(r)$
, versus centre separation
$r/D_{eq}$
for the simulated cases in three regions, where regions I to III refer to
$y/h<1/3$
,
$1/3<y/h<2/3$
and
$2/3<y/h<1$
, respectively.
Next, we analyse the particle relative orientation by calculating the orientation correlation function (
$OCF$
),
$$\begin{eqnarray}OCF(r)=\langle 2|\hat{\boldsymbol{o}}_{p}\cdot \hat{\boldsymbol{o}}_{q}|-1\rangle ,\end{eqnarray}$$
where
$\hat{\boldsymbol{o}}_{p}$
and
$\hat{o} _{q}$
denote the orientations of particles
$p$
and
$q$
at distance
$r$
between their centres. Value
$OCF(r)=1$
indicates particle pairs at distance
$r$
perfectly aligned, while
$OCF(r)=-1$
indicates particles with the symmetry axes orthogonal to each other. This observable is zero for a suspension with random particle orientation. Figure 13 depicts
$OCF$
versus the separation
$r/D_{eq}$
for the simulated cases in three regions, where regions I to III refer to the near-wall region,
$y/h<1/3$
, a region between the wall and the core,
$1/3<y/h<2/3$
, and the region around the channel centreline,
$2/3<y/h<1$
. The results show that particles become more aligned with respect to each other as the volume fraction
$\unicode[STIX]{x1D719}$
increases; the alignment increases as the distance from the wall reduces. A small peak is observed in figure 13 around
$1.2<r/D_{eq}<1.5\approx 2r_{max}$
, which indicates particles sitting on top of each other, parallel to the wall and with their centres shifted by less than the major oblate diameter. The particles therefore seem to form structures that are more reluctant to rotate.
The results so far indicate that oblate particles experience considerably smaller angular velocities close to the wall (except for
$\overline{|\unicode[STIX]{x1D6FA}_{y}|}$
) and stay prevalently with their major axes parallel to the wall. These two facts may explain the turbulence attenuation observed in the presence of oblate particles. As the volume fraction
$\unicode[STIX]{x1D719}$
increases, oblate particles create a sort of strong shield that dampens the turbulent activity close to the wall, preventing the outer layer turbulence from directly interacting with the wall. Indeed, in the simulations reported in Ardekani et al. (Reference Ardekani, Costa, Breugem and Brandt2016a
), we examined the effect of the particle rotation on the turbulence dynamics. To this end, we have performed a simulation at
$\unicode[STIX]{x1D719}=10\,\%$
in which the oblate particles can translate but not rotate and are kept always parallel to the wall (semi-minor axis normal to the wall).
Reynolds stress profile of the combined phase for the non-rotating (NR) oblate case, compared to the previous cases with free oblates.
The results of this simulation display an even larger turbulence attenuation; see figure 14 where the Reynolds stresses of the combined phase of the non-rotating (NR) and freely rotating oblate particles are compared. As a consequence of the reduced Reynolds stresses, the friction Reynolds number considerably decreases,
$Re_{\unicode[STIX]{x1D70F}}=163$
, corresponding to 18 % drag reduction, as also reported in figure 3. It should be noted that it might be feasible to operate with a suspension of particles kept at a preferential orientation using an external magnetic field (Rosensweig Reference Rosensweig2013).
4 Final remarks
We have reported results from simulations of turbulent channel flow of suspensions of finite-size oblate spheroidal particles with
${\mathcal{A}}{\mathcal{R}}=1/3$
at different volume fractions, 5, 7.9, 10 and 15 %. The numerical approach proposed by Ardekani et al. (Reference Ardekani, Costa, Breugem and Brandt2016b
) is used for the simulations of oblate spheroidal particles presented here. This is based on the IBM method for the fluid–solid interactions with lubrication and contact models for the short-range particle–particle (particle–wall) interactions.
Drag reduction with respect to single-phase turbulence is found in the presence of oblate particles, contrary to the results for spherical particles shown in Picano et al. (Reference Picano, Breugem and Brandt2015), for which the drag increases with the particle volume fraction. The two important factors determining the overall drag are the turbulence activity in the suspension and the particle-induced stresses, which are due not only to the increase of the suspension effective viscosity in the presence of particles but also to their distribution across the channel. In fact, Picano et al. (Reference Picano, Breugem and Brandt2015) have shown that for suspensions of spheres a particle layer forms close to the wall, causing high particle-induced stresses in this region. The particle-induced stresses compensate for the reduced Reynolds shear stress observed at high volume fractions (
$\unicode[STIX]{x1D719}>10\,\%$
), resulting in an overall drag increase. Costa et al. (Reference Costa, Picano, Brandt and Breugem2016) indeed showed that a combination of higher effective viscosity and the presence of the mentioned particle layer is at the origin of the overall drag increase in suspensions of spheres. Interestingly, such a particle layer is not present in the flow laden with oblate particles and this makes the turbulence activity the most important factor determining the overall drag.
Attenuation of the turbulence activity is observed in the presence of oblate particles, an effect more pronounced as the volume fraction increases. It is well known that the effective viscosity of a particle suspension is always higher than that of the single-phase flow, which can cause turbulence attenuation. In this study, however, we show that the presence of oblate particles reduces the turbulence activity to lower values than those obtained by accounting only for the effective suspension viscosity. To show this, we perform simulations of single-phase flow at bulk Reynolds numbers calculated with the effective viscosity for the corresponding particle volume fraction, concluding that the specific dynamics of the oblate particles and their interactions with the turbulent velocity field considerably reduce the turbulence activity of the suspension. In fact, turbulence attenuation is sufficiently high that, despite the increase in effective viscosity, overall drag reduction is achieved.
We explain the turbulence attenuation observed in the presence of oblate particles by noting that oblate particles experience considerably lower angular velocities close to the wall and stay prevalently aligned parallel to the wall. At high volume fractions, the oblate particles create a kind of strong shield that dampens the turbulent fluctuations close to the wall, separating the outer layer turbulence from direct interaction with the wall. To show the isolated effect of the particle rotation and the orientation on the turbulent dynamics, we perform a simulation at
$\unicode[STIX]{x1D719}=10\,\%$
in the same flow geometry, in which the oblate particles are free to translate but cannot rotate and are kept always parallel to the wall (semi-minor axis normal to the wall). This simulation reveals an even larger turbulence attenuation and therefore increased drag reduction when compared to the cases where oblate particles are free to rotate.
Examining the particle relative orientation, e.g. by calculating the nematic order and the orientation correlation function of the centre separation
$r$
between the particles, we show that the particles are preferentially parallel to the wall in its vicinity and are also increasingly parallel to each other as the volume fraction
$\unicode[STIX]{x1D719}$
increases, with the highest relative alignment in the region close to the wall. The relative alignment in this region is of utmost importance in reducing the local turbulence fluctuations. The analysis of the forces acting on the particles further shows that particles aligning with the wall are attracted towards it, while they are lifted towards the centreline when inclined.
We have shown here that finite-size oblates reduce drag in turbulent channel flow as already observed for small fibres and polymers. The physical mechanisms are, however, different, as the latter two seem to act on the small scales in the flow. In the future, it will therefore be interesting to study the interactions among finite-size fibres or prolate particles and turbulence as done in Do-Quang et al. (Reference Do-Quang, Amberg, Brethouwer and Johansson2014) at low volume fractions.
Acknowledgements
This work was supported by the European Research Council grant no. ERC-2013-CoG-616186, TRITOS. The authors acknowledge computer time provided by SNIC (Swedish National Infrastructure for Computing) and support from the COST Action MP1305: Flowing matter.
Appendix. Laminar flow of suspensions
In this appendix, we report results for the laminar flow of suspensions of spheres and oblate particles. The data are compared here and used as reference when discussing the turbulent cases above. The simulations are performed at
$Re_{b}=1000$
and
$\unicode[STIX]{x1D719}=10\,\%$
in the same pressure-driven plane channel.
Mean local volume fraction
$\unicode[STIX]{x1D6F7}(y)$
, normalized by total volume fraction
$\unicode[STIX]{x1D719}$
, versus
$y/h$
for the laminar flow of spherical and oblate particles in channel.
The suspension effective viscosity, quantified here by the friction Reynolds number,
$Re_{\unicode[STIX]{x1D70F}}$
, in analogy to the turbulent cases, is slightly larger for the oblate particles,
$Re_{\unicode[STIX]{x1D70F}}=44.6$
, against
$Re_{\unicode[STIX]{x1D70F}}=43.3$
for spheres (i.e. 15 % increase for oblates and 12 % for spheres with respect to the laminar single-phase flow). Figure 15 displays the wall-normal profiles of the local volume fraction,
$\unicode[STIX]{x1D6F7}$
, for the two types of particles under consideration. The data clearly show that spheres migrate towards the channel centre, displaying a local maximum close to the wall, as also reported by Lashgari et al. (Reference Lashgari, Picano, Breugem and Brandt2016); interestingly, this migration disappears for disc-like particles whose distribution is more uniform throughout the channel. Two local maxima are observed for oblates: one close to the wall (
$y/h=0.1$
) and the other one close to the channel centre (
$y/h=0.8$
). The higher local volume fraction close to the wall can explain the higher drag for the laminar flow of oblate particles. The analysis of the forces and torques presented in the text has been repeated for laminar flow to observe that the difference between the positive and negative lift forces is much less than in the turbulent regime, probably due to the lower velocity gradients close to the wall, which may explain why the particles have a more or less uniform distribution.
(a) The cosine of the mean particle inclination angle
$\unicode[STIX]{x1D703}$
, measured with respect to the wall, versus
$y/h$
. Mean particle absolute value of angular velocity, normalized with
$U_{b}/h$
in (b) spanwise direction
$\overline{|\unicode[STIX]{x1D6FA}_{z}|}$
, (c) streamwise direction
$\overline{|\unicode[STIX]{x1D6FA}_{x}|}$
and (d) wall-normal direction
$\overline{|\unicode[STIX]{x1D6FA}_{y}|}$
.
(a) The nematic order parameter
$\unicode[STIX]{x1D706}$
and (b) the biaxial parameter
$\unicode[STIX]{x1D701}$
versus
$y/h$
for the laminar cases studied here. Results for spheres are also depicted to confirm their fully isotropic particle orientation.
Orientational correlation function
$OCF(r)$
versus centre separation
$r/D_{eq}$
for the laminar cases in three regions, where regions I to III refer to
$y/h<1/3$
,
$1/3<y/h<2/3$
and
$2/3<y/h<1$
, respectively.
The mean particle orientation and rotation in laminar flow are compared in figure 16. The tendency of oblate particles to align with the gradient direction in the regions far from the wall is observed to be more pronounced in the laminar regime than in turbulent flow. The mean absolute value of the particle angular velocities, depicted in figure 16(b–d), shows higher values for oblate particles in the streamwise and the wall-normal directions close to the wall and smaller for rotation rates in the spanwise direction.
The nematic order parameter
$\unicode[STIX]{x1D706}$
and the biaxial parameter
$\unicode[STIX]{x1D701}$
are depicted in figure 17 versus the channel height. The profile of
$\unicode[STIX]{x1D706}$
in figure 17(a) is similar to the mean orientation profile of oblate particles in figure 16(a); values close to 1 in both figures show the tendency of the particles to be aligned with the symmetry axis normal to the wall. Figure 17(b) displays two local maxima of the biaxial parameter
$\unicode[STIX]{x1D701}$
: one close to the wall and the other one close to the channel centre. Interestingly, these local maxima appear to correspond to the local maxima of the local volume fraction profiles in figure 15, revealing that the tendency to be oriented also in the spanwise direction increases in those locations with higher local volume fractions.
Finally, figure 18 reports the orientational correlation function,
$OCF$
, versus the particle centre separation
$r/D_{eq}$
in three regions, where regions I to III refer to
$y/h<1/3$
,
$1/3<y/h<2/3$
and
$2/3<y/h<1$
. The results for laminar flow show that the particle orientations are more correlated with respect to the turbulent cases – in other words that the order is disrupted by turbulent mixing.
