1. Introduction
Suspensions of solid particles in liquid flows are widely encountered in industrial application and environmental problems. Sediment transport, avalanches, slurries, pyroclastic flows, oil industry and pharmaceutical processes represent typical examples where a step forward in the understanding and modelling of these complex fluids is essential. Given the high flow rates typically encountered in these applications, inertia strongly influences the flow regime that may be chaotic and turbulent. The main aim of the present work is therefore to investigate the interactions between the phases of a suspension in the turbulent regime.
Suspensions are often constituted by a Newtonian liquid laden with solid particles that may differ for size, shape, density and stiffness. Even restricting our analysis to monodisperse rigid neutrally buoyant spheres, the laminar flow of these suspension shows peculiar rheological properties, such as high effective viscosities, normal stress differences, shear thinning or thickening and jamming at high volume fractions, see e.g. Stickel & Powell (Reference Stickel and Powell2005), Morris (Reference Morris2009) and Wagner & Brady (Reference Wagner and Brady2009) for recent reviews on the topic. In particular, still dealing with simple laminar flows, the suspended phase alters the response of the complex fluid to the local deformation rate leading, for example, to an increase of the effective viscosity of the suspension
${\it\mu}_{e}$
with respect to that of the pure fluid
${\it\mu}$
(Guazzelli & Morris Reference Guazzelli and Morris2011). A first attempt to characterise this effect can be traced back to Einstein (Reference Einstein1906, Reference Einstein1911) who provided a linear estimate of the effective viscosity
${\it\mu}_{e}={\it\mu}\,(1+2.5\,{\it\Phi})$
, with
${\it\Phi}$
the volume fraction, valid in the dilute regime. Few decades later, Batchelor (Reference Batchelor1970) and Batchelor & Green (Reference Batchelor and Green1972) derived and proposed a quadratic correction that partially accounts for the mutual interactions among particles, which become increasingly critical when increasing the volume fraction. Indeed, the suspension viscosity increases by more than one order of magnitude in the dense regime, until the system jams behaving as a glass or a crystal (Sierou & Brady Reference Sierou and Brady2002). For dense cases only semi-empirical laws exist for the effective viscosity; the mixture viscosity has been observed to diverge when the system approaches the maximum packing limit
${\it\Phi}_{m}=0.58{-}0.62$
(Boyer, Guazzelli & Pouliquen Reference Boyer, Guazzelli and Pouliquen2011), as reproduced by empirical fits such as those by Eilers and Kriegher & Dougherty (Stickel & Powell Reference Stickel and Powell2005).
The rheological properties of suspensions have been often studied in the viscous Stokesian regime where inertial effects are negligible and can be safely neglected. Nonetheless in several applications the flow Reynolds number is high enough that the inertia is significant at the particle scale. The seminal work of Bagnold (Reference Bagnold1954) on the highly inertial regime revealed how the increase of the particle collisions induces an effective viscosity that increases linearly with the shear rate. Even if the macroscopic flow is viscous and laminar, inertial effects at the particle scale may induce shear-thickening (Kulkarni & Morris Reference Kulkarni and Morris2008b
; Picano et al.
Reference Picano, Breugem, Mitra and Brandt2013) or normal stress differences (Zarraga, Hill & Leighton Reference Zarraga, Hill and Leighton2000). This change of the macroscopic behaviour is due to a strong modification of the particle microstructure, i.e. the relative position and velocity of the suspended particles (Morris Reference Morris2009; Picano et al.
Reference Picano, Breugem, Mitra and Brandt2013). A finite particle-scale Reynolds number,
$\mathit{Re}_{a}>0$
, breaks the symmetry of the particle pair trajectories (Kulkarni & Morris Reference Kulkarni and Morris2008a
; Picano et al.
Reference Picano, Breugem, Mitra and Brandt2013) and induces an anisotropic microstructure, in turns responsible of shear-thickening.
It is well established that the macroscopic flow behaviour changes dramatically from the laminar conditions to the typical chaotic dynamics of transitional and turbulent flows when increasing the Reynolds number, still for single-phase fluids. The effect of a dense suspended phase on the transition to turbulence in pipe flows has been investigated experimentally by Matas, Morris & Guazzelli (Reference Matas, Morris and Guazzelli2003). These authors report a non-trivial behaviour of the critical Reynolds number at which transition is observed. The critical Reynolds number for relatively large particles is found to first decrease and then increases, with a minimum in the range
${\it\Phi}\sim 0.05{-}0.1$
. This non-monotonic behaviour cannot be explained only in terms of the increase of the suspension effective viscosity. These experiments have been numerically reproduced in Yu et al. (Reference Yu, Wu, Shao and Lin2013). Recently, Lashgari et al. (Reference Lashgari, Picano, Breugem and Brandt2014) showed that the flow behaviour is more complex than that pertaining to unladen flows: three different regimes coexist with different probability when changing the volume fraction
${\it\Phi}$
and the Reynolds number
$\mathit{Re}$
. In each regime the flow is dominated by viscous, turbulent and particle stresses, respectively.
As far as the turbulent regime is concerned, most part of the previous studies pertains to the dilute or very dilute regimes. In the very dilute regime, the particle concentration is so small that the solid phase has a negligible effect on the flow. In this, so-called, one-way coupling regime, the main object of most investigations is the particle transport properties. In particular, inertia affects the particle turbulent dispersion leading to preferential migration. Small-scale clustering has been observed both in isotropic (see e.g. Toschi & Bodenschatz Reference Toschi and Bodenschatz2009) and inhomogeneous flows (see e.g. Sardina et al. Reference Sardina, Schlatter, Brandt, Picano and Casciola2012). It amounts to a segregation of the particles in fractal sets (Bec et al. Reference Bec, Biferale, Cencini, Lanotte, Musacchio and Toschi2007; Toschi & Bodenschatz Reference Toschi and Bodenschatz2009) induced by the coupling of the turbulent flow dynamics (dissipative) and the particle inertia when the time scales of the two phenomena are similar. In wall-bounded flow, particle inertia induces a mean particle drift towards the wall, so-called turbophoresis (Reeks Reference Reeks1983). This effect is most pronounced when the particle inertial time scale almost matches the turbulent near-wall characteristic time (Soldati & Marchioli Reference Soldati and Marchioli2009). Clustering and turbophoresis interact leading to the formation of streaky particle patterns (e.g. Sardina et al. Reference Sardina, Picano, Schlatter, Brandt and Casciola2011).
Increasing the solid phase concentration, while still keeping small the volume fraction and the particle diameter with respect to the flow length scales, the flow satisfies the so-called two-way coupling approximation, see among others Ferrante & Elghobashi (Reference Ferrante and Elghobashi2003) and Balachandar & Eaton (Reference Balachandar and Eaton2010). This regime is characterised by high mass density ratios, i.e. the ratio between the mass of the solid phase and the fluid one, and low volume fractions (Balachandar & Eaton Reference Balachandar and Eaton2010) in the limit of high mass fractions; this occurs typically for solid particles or droplets dispersed in a gas phase when the density ratio between particles and fluid is high (approximately 1000). In this regime the dispersed phase back-reacts on the carrier fluid exchanging momentum, with inter-particle interactions and excluded volume effects being negligible given the small volume fractions. In homogeneous and isotropic flows, Squires & Eaton (Reference Squires and Eaton1991) and Elghobashi & Truesdell (Reference Elghobashi and Truesdell1993) observe an attenuation of the turbulent kinetic energy at large scales accompanied by an energy increase at small scales. Sundaram & Collins (Reference Sundaram and Collins1999) and Ferrante & Elghobashi (Reference Ferrante and Elghobashi2003) also performed systematic studies to understand the effect of the particle inertia and of the mass fraction on the flow. Gualtieri et al. (Reference Gualtieri, Picano, Sardina and Casciola2013) report that the particle segregation in anisotropic fractal sets induces an alternative mechanism to directly transfer energy from large to small scales. Similar results have been reported for wall-bounded turbulent flows. Kulick, Fessler & Eaton (Reference Kulick, Fessler and Eaton1994) showed that the solid phase reduces the turbulent near-wall fluctuations increasing their anisotropy, see also Li et al. (Reference Li, McLaughlin, Kontomaris and Portela2001). Zhao, Andersson & Gillissen (Reference Zhao, Andersson and Gillissen2010) showed how these interactions may lead to drag reduction.
If the dispersed phase is not constituted by elements smaller than the hydrodynamic scales, the suspended phase directly affects the turbulent structures at scales similar or below the particle size (Naso & Prosperetti Reference Naso and Prosperetti2010; Bellani et al. Reference Bellani, Byron, Collignon, Meyer and Variano2012; Homann, Bec & Grauer Reference Homann, Bec and Grauer2013). As the system is nonlinear and chaotic, these large-scale interactions modulate the whole process inducing non-trivial effects on the turbulence cascade (Lucci, Ferrante & Elghobashi Reference Lucci, Ferrante and Elghobashi2010; Yeo et al. Reference Yeo, Dong, Climent and Maxey2010) where increase or decrease of the spectral energy distribution depends on the particle size and mass fraction. Pan & Banerjee (Reference Pan and Banerjee1996) were the first to simulate the effect of finite-size particles in a turbulent channel flow showing that when these are larger than the dissipative scale turbulent fluctuations and stresses become larger. The open-channel flow laden with heavy finite-size particles has been investigated in the dilute regime by Kidanemariam et al. (Reference Kidanemariam, Chan-Braun, Doychev and Uhlmann2013) and Kidanemariam & Uhlmann (Reference Kidanemariam and Uhlmann2014) showing that the solid phase preferentially accumulates in near-wall low-speed streaks, the flow structures characterised by smaller streamwise velocity.
Increasing the volume fraction, the coupling among the phases becomes richer and particle–particle hydrodynamic interactions and collisions cannot be neglected. In this dense regime, so-called four-way coupling, the rheological properties of the suspension interact with the chaotic dynamics of the fluid phase when the flow inertia is sufficiently large, i.e. at high Reynolds number. Few studies investigate dense suspensions in the highly inertial regime: Matas et al. (Reference Matas, Morris and Guazzelli2003), Loisel et al. (Reference Loisel, Abbas, Masbernat and Climent2013) and Yu et al. (Reference Yu, Wu, Shao and Lin2013) show the effect on transition in wall-bounded flows showing a decrease of the critical Reynolds number in the semidilute regime. Concerning the turbulent regime of relatively dense suspensions of wall-bounded flows, Shao, Wu & Yu (Reference Shao, Wu and Yu2012) report results for channel flow up to 7 % volume fraction both considering neutrally buoyant and heavy particles. These authors document a decrease of the fluid streamwise velocity fluctuation due to an attenuation of the large-scale streamwise vortices. In the case of heavy, sedimenting, particles, the bottom wall behaves as a rough boundary with particles free to resuspend. Different regimes have been observed when the importance of the particle buoyancy is varied in the recent study of Vowinckel, Kempe & Fröhlich (Reference Vowinckel, Kempe and Fröhlich2014).
In this context, even restricting to the case of neutrally buoyant particles, little is know on the effect of a dense suspended phase on the fully turbulent regime. The main reason can be ascribed to the well known difficulties to tackle this case either experimentally or numerically. As noticed previously, the dense regime is characterised by a complex particle microstructure that induces non-trivial macroscopic features. When the large-scale inertia is high enough, the interaction between the suspension microstructure, i.e. rheology, and turbulence dynamics is expected to significantly alter the macroscopic flow dynamics. This is the object of the present study.
To this end, we consider turbulent channel flows laden with finite-size particles (radius
$a=h/18$
with
$h$
the half-channel height) up to a volume fraction
${\it\Phi}=0.2$
. We use data from a direct numerical simulation (DNS) that fully describe the solid phase dynamics via an immersed boundary method (IBM). We show that the classical laws for the turbulent mean velocity profiles are modified in the presence of the particles and the overall drag increases. At the highest volume fraction investigated,
${\it\Phi}=0.2$
, the velocity fluctuation intensities and the Reynolds shear stresses are found to suddenly decrease. We consider the mean momentum budget to show that the particle-induced stress is responsible of the overall drag increase at high
${\it\Phi}$
, while the turbulent drag decreases.
2. Methodology
2.1. Numerical algorithm
During recent years, different methods have been proposed to perform accurate DNS of dense multiphase flows. Fully Eulerian methods have been adopted to deal with two-fluid flows, such as front-tracking, sharp or diffuse interface methods (see e.g. Tryggvason et al. Reference Tryggvason, Bunner, Esmaeeli, Juric, Al-Rawahi, Tauber, Han, Nas and Jan2001; Bray Reference Bray2002; Benzi et al. Reference Benzi, Sbragaglia, Succi, Bernaschi and Chibbaro2009; Celani et al. (Reference Celani, Mazzino, Muratore-Ginanneschi and Vozella2009); Magaletti et al. Reference Magaletti, Picano, Chinappi, Marino and Casciola2013), whereas mixed Lagrangian–Eulerian techniques are found to be the most appropriate for solid–liquid suspensions (Ladd & Verberg Reference Ladd and Verberg2001; Takagi et al. Reference Takagi, Oguz, Zhang and Prosperetti2003; Lucci et al. Reference Lucci, Ferrante and Elghobashi2010; Kidanemariam & Uhlmann Reference Kidanemariam and Uhlmann2014; Vowinckel et al. Reference Vowinckel, Kempe and Fröhlich2014). In this framework, the present simulations have been performed with a numerical code that fully describes the coupling between the solid and fluid phases (Breugem Reference Breugem2012). The Eulerian fluid phase evolves according to the incompressible Navier–Stokes equations,
$$\begin{eqnarray}\displaystyle & \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}_{f}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial \boldsymbol{u}_{f}}{\partial t}+\boldsymbol{u}_{f}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{f}=-\frac{1}{{\it\rho}}\boldsymbol{{\rm\nabla}}p+{\it\nu}{\rm\nabla}^{2}\boldsymbol{u}_{f}+\boldsymbol{f}, & \displaystyle\end{eqnarray}$$
$\boldsymbol{u}_{f}$
is the fluid velocity,
$\boldsymbol{f}$
a generic force field,
$p$
the pressure,
${\it\nu}={\it\mu}/{\it\rho}$
the kinematic viscosity of the pure fluid with
${\it\mu}$
the dynamic viscosity and
${\it\rho}$
the fluid density (same as the particle density in this study). The solid phase consists of neutrally buoyant rigid spheres whose centroid linear and angular velocities,
$\boldsymbol{u}_{p}$
and
${\bf\omega}_{p}$
, are governed by the Newton–Euler Lagrangian equations, 
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}V_{p}\frac{\text{d}\boldsymbol{u}_{p}}{\text{d}t}={\it\rho}\oint _{\partial \mathscr{V}_{p}}{\bf\tau}\boldsymbol{\cdot }\boldsymbol{n}\text{d}S, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle I_{p}\frac{\text{d}{\bf\omega}_{p}}{\text{d}t}={\it\rho}\oint _{\partial \mathscr{V}_{p}}\boldsymbol{r}\times {\bf\tau}\boldsymbol{\cdot }\boldsymbol{n}\text{d}S, & \displaystyle\end{eqnarray}$$
$a$
is the particle radius and
$V_{p}=4{\rm\pi}a^{3}/3$
the particle volume; the fluid stress is
${\bf\tau}=-p\boldsymbol{I}+2{\it\mu}\unicode[STIX]{x1D640}$
with
$\unicode[STIX]{x1D640}=\left(\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{f}+\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{f}^{\text{T}}\right)\!/2$
the deformation tensor. In (2.4),
$I_{p}=(2/5){\it\rho}\mathscr{V}_{p}a^{2}$
represents the moment of inertia,
$\boldsymbol{r}$
the distance vector from the centroid of the sphere and
$\boldsymbol{n}$
the unity vector normal to the particle surface
$\partial \mathscr{V}_{p}$
. On the particle surfaces, Dirichlet boundary conditions for the fluid phase are enforced as
$\boldsymbol{u}_{f}|_{\partial \mathscr{V}_{p}}=\boldsymbol{u}_{p}+{\bf\omega}_{p}\times \boldsymbol{r}$
.
In the simulations reported in this paper, the coupling between the two phases is obtained by using an IBM: this amounts to adding a force field
$\boldsymbol{f}$
on the right-hand side of (2.2) to mimic the actual boundary condition at the moving particle surface, i.e.
$\boldsymbol{u}_{f}|_{\partial \mathscr{V}_{p}}=\boldsymbol{u}_{p}+{\bf\omega}_{p}\times \boldsymbol{r}$
. The fluid phase is evolved solving (2.1) and (2.2) in a domain containing all of the particles, without the need to adapt the mesh to the current particle position, using a second-order finite difference scheme on a staggered mesh. The time integration is performed by a third-order Runge–Kutta scheme combined with a pressure-correction method on each substep. The Lagrangian evolution of (2.3) and (2.4) is performed using the same Runge–Kutta scheme of the Eulerian solver. The particle surface is tracked using
$N_{L}$
Lagrangian points uniformly distributed on the surface of the spheres on which the forces exchanged with the fluid phase are imposed. To maintain accuracy, the right-hand side of (2.3) and (2.4) are rearranged in terms of the IBM force field and take into account the mass of the fictitious fluid phase occupied by the particle volumes
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}V_{p}\frac{\text{d}\boldsymbol{u}_{p}}{\text{d}t}=-{\it\rho}\mathop{\sum }_{l=1}^{N_{L}}\boldsymbol{F}_{l}\,{\rm\Delta}V_{l}+{\it\rho}\frac{\text{d}}{\text{d}t}\left(\int _{V_{p}}\boldsymbol{u}_{f}\text{d}V\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle I_{p}\frac{\text{d}{\bf\omega}_{p}}{\text{d}t}=-{\it\rho}\mathop{\sum }_{l=1}^{N_{L}}\boldsymbol{r}_{\boldsymbol{l}}\times \boldsymbol{F}_{l}\,{\rm\Delta}V_{l}+{\it\rho}\frac{\text{d}}{\text{d}t}\left(\int _{V_{p}}\boldsymbol{r}\times \boldsymbol{u}_{f}\text{d}V\right), & \displaystyle\end{eqnarray}$$
${\rm\Delta}V_{l}$
is the volume of the cell around the
$l$
Lagrangian point and
$\boldsymbol{r}_{l}$
the distance from the particle centre. Here
$\boldsymbol{F}_{l}$
is the force acting on the
$l$
Lagrangian point on the particle and is related to the Eulerian force field
$\boldsymbol{f}\!:\boldsymbol{f}(\boldsymbol{x})=\sum _{l=1}^{N_{L}}F_{l}\,{\it\delta}_{d}(\boldsymbol{x}-X_{l}){\rm\Delta}V_{l}$
. The procedure to determine the force field from the boundary conditions at the particle surface follows an iterative algorithm that allows the code to achieve second order global accuracy in space. All of the details of this implementation are presented in Breugem (Reference Breugem2012).
The numerical method models the interaction among the particles also when their gap distance is of the order of or below the grid size. In particular, lubrication models based on the Brenner’s asymptotic solution (Brenner Reference Brenner1961) are used to correctly reproduce the interaction between particles when their gap distance is smaller than twice the mesh size. When particles collide with the wall or among themselves a soft-collision model ensures an almost elastic rebound with a restitution coefficient set at 0.97. A complete discussion of these models can be found in Breugem (Reference Breugem2012) and Lambert et al. (Reference Lambert, Picano, Breugem and Brandt2013) where several test cases are presented as validation.
To avoid duplication of published material, we provide here only evidence for the ability of the present numerical tool to accurately simulate dense suspensions. Figure 1 displays the relative viscosity, the ratio between the effective viscosity of the suspension and the viscosity of the fluid phase
${\it\nu}_{r}={\it\nu}_{e}/{\it\nu}$
, in laminar flows for two different volume fractions,
${\it\Phi}=0.2$
and
${\it\Phi}=0.3$
. The configuration where this is measured is the Couette flow at vanishing Reynolds number where the wall-to-wall distance is 10 times the particle radius. A cubic mesh is used to discretise the computational domain with eight points per particle radius,
$a$
. The streamwise and spanwise length of the computational domain are 1.6 times the wall-normal width, i.e.
$16a$
. The relative viscosity extracted after the initial transient phase is measured by the friction at the wall and perfectly matches previous numerical investigations (Yeo & Maxey Reference Yeo and Maxey2010a
,Reference Yeo and Maxey
b
) and empirical fits of experimental data, such as the Eilers fit (Stickel & Powell Reference Stickel and Powell2005).
Relative viscosity
${\it\nu}_{r}$
versus the volume fraction
${\it\Phi}$
in a Couette flow in absence of inertia,
$\mathit{Re}=0$
. The present data are compared with the results by Yeo & Maxey (Reference Yeo and Maxey2010a
,Reference Yeo and Maxey
b
) and the Eilers fit
$(1+1.25{\it\Phi}/(1-{\it\Phi}/0.65))$
.
2.2. Flow configuration
In this work we study a pressure-driven channel flow between two infinite flat walls located at
$y=0$
and
$y=2h$
with
$y$
the wall-normal direction. Periodic boundary conditions are imposed in the streamwise,
$x$
, and spanwise,
$z$
, directions for a domain size of
$L_{x}=6h$
,
$L_{y}=2h$
and
$L_{z}=3h$
. A mean pressure gradient acting in the streamwise direction imposes a fixed value of the bulk velocity
$U_{0}$
across the channel corresponding to a constant bulk Reynolds number
$\mathit{Re}_{b}=U_{0}2h/{\it\nu}=5600$
, with
${\it\nu}$
the kinematic viscosity of the fluid phase; this value corresponds to a Reynolds number based on the friction velocity
$\mathit{Re}_{{\it\tau}}=U_{\ast }h/{\it\nu}=180$
for the single-phase case where
$U_{\ast }=\sqrt{{\it\tau}_{w}/{\it\rho}}$
with
${\it\tau}_{w}$
the stress at the wall. As reported in table 1, the bulk Reynolds number based on the suspension effective viscosity
$\mathit{Re}_{e}$
varies with the volume fraction following the increase of the effective viscosity
${\it\nu}_{e}={\it\nu}_{r}\,{\it\nu}$
where
${\it\nu}_{r}$
is the relative viscosity estimated by the Eilers fit (Stickel & Powell Reference Stickel and Powell2005). The domain is discretised by a cubic mesh of
$864\times 288\times 432$
points in the streamwise, wall-normal and spanwise directions. Hereafter all of the variables have been made dimensionless with
$U_{0}$
and
$h$
, except those with the superscript ‘
$\text{}^{+}$
’ that are scaled with
$U_{\ast }$
and
${\it\delta}_{\ast }={\it\nu}/U_{\ast }$
(inner scaling).
Summary of the DNS reported here. They pertain to suspensions of
$N_{p}$
particles of radius
$a/h=1/18$
at different volume fractions
${\it\Phi}$
. Here
$N_{x},N_{y},N_{z}$
indicate the number of grid points in each direction and the bulk Reynolds number is defined as
$\mathit{Re}_{b}=U_{0}\ast 2h/{\it\nu}$
. The relative viscosity, i.e. the ratio between effective suspension viscosity and the fluid viscosity
${\it\nu}_{r}={\it\nu}_{e}/{\it\nu}=[1+1.25\ast {\it\Phi}/(1-{\it\Phi}/0.6)]^{2}$
has been estimated via the Eilers fit (Stickel & Powell Reference Stickel and Powell2005). The effective bulk Reynolds number is defined as
$\mathit{Re}_{e}=U_{0}2h/{\it\nu}_{e}=U_{0}2h/({\it\nu}\,{\it\nu}_{r})=\mathit{Re}_{b}/{\it\nu}_{r}$
.
Non-Brownian spherical neutrally buoyant rigid particles are considered. The ratio between the particle radius and the channel half-width is fixed to
$a/h=1/18$
, corresponding to 10 plus units for the lowest volume fraction considered and 12 for the largest. Three different volume fractions,
${\it\Phi}=0.05;0.1;0.2$
, have been examined in addition to the single phase case for a direct comparison. The highest volume fraction here addressed requires 10 000 finite-size particles in the computational domain with
$N_{l}=746$
Lagrangian control points on the surface of each sphere and eight Eulerian grid points per particle radius, see table 1. The simulations were run on a Cray XE6 system using 2048 cores for a total of approximately
$10^{6}$
CPU hours for each case.
The simulation starts from the laminar Poiseuille flow for the fluid phase and a random positioning of the particles. Transition naturally occurs at the fixed Reynolds number because of the noise added by the presence of the particles. Statistics are collected after the initial transient phase.
3. Results
3.1. Single-point flow and particle velocity statistics
Snapshots of the suspension flow are shown in figure 2 for the different nominal volume fractions
${\it\Phi}$
under investigation. The instantaneous streamwise velocity is represented on different orthogonal planes with the bottom plane located in the viscous sublayer to highlight the low- and high-speed streaks characteristic of near-wall turbulence. Finite-size particles are displayed only on one half of the domain to give a visual feeling on how dense the solid phase is for the different
${\it\Phi}$
. Indeed, at the highest volume fraction,
${\it\Phi}=0.2$
, the particles are so dense that completely hide the bottom wall. The cases with
${\it\Phi}=0.05$
and
${\it\Phi}=0.1$
show velocity contours similar to those of the unladen case where it is possible to recognise the typical near-wall streamwise velocity streaks; these are however more noisy and characterised by significant small-scale modulations (of particle size). At
${\it\Phi}=0.2$
the small-scale noise is stronger and the streaks become wider.
Instantaneous snapshots of the streamwise velocity on different orthogonal planes together with the corresponding particle position represented only on one half of the domain. The four panels represent the different values of the volume fraction under investigation, (a)
${\it\Phi}=0$
, (b)
${\it\Phi}=0.05$
, (c)
${\it\Phi}=0.1$
and (d)
${\it\Phi}=0.2$
.
The mean fluid velocity profiles are shown in figure 3. The statistics conditioned to the fluid phase have been calculated considering only the points located out of the volume occupied by the particles in each field (phase-ensemble average). Figure 3(a) reports the velocity in outer units
$U_{f}$
, indicating that the maximum velocity at the mid-plane grows with
${\it\Phi}$
(note that the flow rate is constant in these simulations). In general, the mean velocity more closely resembles the laminar parabolic profile when increasing the volume fraction: the velocity increases in the centre of the channel at higher
${\it\Phi}$
, whereas it decreases near the wall, up to
$y\sim 1/2$
. The higher the volume fraction the more intense this effect is. Figure 3(b) displays the mean fluid velocity profiles scaled in inner units in the log-linear scale,
$U_{f}^{+}=U_{f}/U_{\ast }$
versus
$y^{+}=y/{\it\delta}_{\ast }$
where the friction velocity and viscous length are
$U_{\ast }=\sqrt{{\it\tau}_{w}/{\it\rho}}$
and
${\it\delta}_{\ast }={\it\nu}/U_{\ast }$
with
${\it\tau}_{w}$
the wall stress. The progressive decrease of the profiles with the volume fraction
${\it\Phi}$
indicates that the overall drag increases. Analysing the flow in terms of the canonical classification of wall turbulence, we can still recognise for all cases a region (
$y^{+}>40{-}50$
) where the mean profile follows a log-law:
$$\begin{eqnarray}U^{+}=(1/k)\log (y^{+})+B\end{eqnarray}$$
with
$k$
the von Kármán constant and
$B$
the additive coefficient. Fits of these constants and the corresponding Friction Reynolds number
$\mathit{Re}_{{\it\tau}}=U_{\ast }h/{\it\nu}$
are reported in table 2 for all of the volume fractions investigated.
Mean fluid velocity profiles for the different volume fractions under investigations in (a) outer units (hereafter wall normal distances without the superscript
$\text{}^{+}$
are assumed to be rescaled by
$h$
) and (b) inner units:
$U_{f}^{+}=U_{f}/U_{\ast }$
versus
$y^{+}=y/{\it\delta}_{\ast }$
, with
$U_{\ast }$
and
${\it\delta}_{\ast }$
the friction velocity and viscous length scale, see the definition in the text.
The von Kármán constant
$k$
and additive constant
$B$
of the log-law estimated from the present simulations for the different volume fractions
${\it\Phi}$
examined. Here
$B$
and
$k$
have been fitted in the range
$y^{+}\in [50,150]$
. The friction Reynolds number
$\mathit{Re}_{{\it\tau}}=U_{\ast }h/{\it\nu}$
, and the effective friction Reynolds number, defined as
$\mathit{Re}_{{\it\tau}}^{e}=U_{\ast }h/{\it\nu}_{e}=\mathit{Re}_{{\it\tau}}/{\it\nu}_{r}$
, are also reported together with an estimate of the effective friction Reynolds number based on the correlation
$\mathit{Re}_{{\it\tau}}^{\prime e}\simeq 0.09\mathit{Re}_{e}^{0.88}$
, see e.g. Pope (Reference Pope2000).
The friction Reynolds number computed from the simulation data differs from what can be estimated using the rheological properties of the suspension, that is using the relative viscosity
${\it\nu}_{r}$
, see table 1. The values of
$\mathit{Re}_{{\it\tau}}^{e}=U_{\ast }h/{\it\nu}_{e}=\mathit{Re}_{{\it\tau}}/{\it\nu}_{r}$
in table 2 are computed using the measured wall friction and the effective viscosity of the suspension. Considering the bulk effective Reynolds number
$\mathit{Re}_{e}=\mathit{Re}_{b}/{\it\nu}_{r}$
, computed in a similar way, it is also possible to estimate an expected value of the friction Reynolds number using the correlation valid in Newtonian flows,
$\mathit{Re}_{{\it\tau}}^{\prime e}\simeq 0.09\mathit{Re}_{e}^{0.88}$
(see Pope Reference Pope2000). The data in table 2 clearly indicate that the effective friction Reynolds number
$\mathit{Re}_{{\it\tau}}^{e}=\mathit{Re}_{{\it\tau}}/{\it\nu}_{r}$
is always higher than what expected considering only the effective viscosity of the suspension, i.e.
$\mathit{Re}_{{\it\tau}}^{\prime e}$
. This fact (
$\mathit{Re}_{{\it\tau}}^{\prime e}<\mathit{Re}_{{\it\tau}}^{e}$
) implies that the particles alter the turbulence and induce an additional dissipation mechanism, as shown by the higher measured wall friction. As shown later, the increased friction can be explained by an increase of the turbulent activity for
${\it\Phi}\leqslant 0.1$
, whereas this is no more the case for the highest volume fraction considered. The increased dissipation at this higher
${\it\Phi}$
may be explained by an increased particle induced stress, i.e. inertial shear thickening (Morris Reference Morris2009; Picano et al.
Reference Picano, Breugem, Mitra and Brandt2013). Inertial shear thickening occurs in a dense suspension when inertial effects are present at the particle scale (finite-particle Reynolds number) and amounts to an increase of the effective viscosity with respect the value obtained by rheological experiments at vanishing inertia (Reynolds number) and same volume fraction. The relatively high Reynolds number of the present turbulent cases triggers inertial effects in the transported particles.
The slope of the log-layer increases, i.e. the von Kármán constant
$k$
decreases, while the additive constant
$B$
decreases. At
${\it\Phi}=0.2$
the differences with respect to the unladen case become critical with
$B$
strongly negative and
$k$
about half of the value for the single phase flow. These two behaviours act in an opposite way: a reduced von Kármán constant
$k$
usually denotes drag reduction (Virk Reference Virk1975), while a small or negative additive constant
$B$
an increase of the drag. The combination of these two counteracting effects lead to an increase of the overall drag for the present cases as demonstrated by the increase of the friction Reynolds number
$\mathit{Re}_{{\it\tau}}$
. The decrease of the additive constant
$B$
appears to be linked to particle–fluid interactions occurring near the wall. In particular, focusing on the case at
${\it\Phi}=0.2$
, we note a sudden change in the mean velocity profile after the first layer of particles, i.e.
$y^{+}\sim 20\sim d_{p}^{+}$
. The near wall dynamics is therefore influenced by the particle layering induced by the wall. A similar behaviour has been observed in turbulent flows over porous media (Breugem, Boersma & Uittenbogaard Reference Breugem, Boersma and Uittenbogaard2006) suggesting that the near-wall layers of particles may act as a porous media for the fluid phase.
It is worth commenting at this point that increasing the bulk Reynolds number usually leads to a widening of the log-law region and, consequently, to a stronger impact of the slope of the log-law on the overall mean velocity profile. Assuming that the constant
$B$
does not change significantly upon increasing the bulk Reynolds number (at fixed
$d^{+}$
), the overall mass flux may increase leading to drag reduction if the log region is long enough for the mean velocity at
${\it\Phi}=0.2$
to become larger than the corresponding values for the single-phase fluid near the channel centreline. This is just a speculation and its proof is out of the scope of the present investigation where we consider only a fixed bulk Reynolds number. Simulations at higher Reynolds number and fixed particle size (in plus units) are currently computationally too expensive and out of our reach.
The root-mean-square (r.m.s.) of the fluid velocity fluctuations and the Reynolds shear stress in outer units are reported in figure 4. We note that despite the increase of the friction Reynolds number the peak of the streamwise velocity r.m.s.,
${u_{f}^{\prime }}_{rms}$
, decreases with
${\it\Phi}$
, while a non-monotonic behaviour is apparent in the bulk of the flow. For values of
${\it\Phi}\leqslant 0.1$
the intensity of the cross-stream velocity fluctuations increases with respect to the single-phase cases, displaying also higher peak values. This indicates that the particle presence redistributes energy towards a more isotropic state. Interestingly, at the highest volume fraction considered,
${\it\Phi}=0.2$
, we note a decrease of the level of fluctuations with respect to all the other cases, with the exception of a thin region close to wall, which will be discussed more in detail in the following. At this high volume fraction we therefore note a reduced turbulence activity, as confirmed by considering the variations of the Reynolds stress in the presence of particles in figure 4(d). Note that the Reynolds stresses represent the main engine for the production of turbulent fluctuations. While these stresses increase for
${\it\Phi}=0.05$
and
${\it\Phi}=0.1$
, they decrease at
${\it\Phi}=0.2$
despite the increase of the friction Reynolds number. At first sight, this aspect may appear controversial, however, as we will discuss in detail in § 3.2, the reduction of the turbulent activity at
${\it\Phi}=0.2$
is associated with an increase of the stresses induced by the solid phase which results in enhanced drag.
Intensity of the fluctuation velocity components and the Reynolds shear stress for the fluid phase in outer units for different volume fractions
${\it\Phi}$
: (a) streamwise
${u_{f}^{\prime }}_{rms}$
; (b) wall-normal
${v_{f}^{\prime }}_{rms}$
; (c) spanwise
${w_{f}^{\prime }}_{rms}$
velocity fluctuations; (d) shear stress
$\langle u_{f}^{\prime }v_{f}^{\prime }\rangle$
.
Further insight into the near wall dynamics can be gained by displaying the same quantities scaled in inner units, see figure 5. The peak of the fluctuation intensity reduces for all of the velocity components when divided by the friction velocity with the only exception of the spanwise component. More importantly, we observe that the fluctuation level monotonically increases with
${\it\Phi}$
in the viscous sublayer. This enhancement of the near-wall fluctuation can be explained by considering the squeezing motions occurring between the wall and an incoming or outgoing particle. We also note that the peak of the Reynolds stresses decreases monotonically (when scaled by the friction velocity squared) becoming about half of the expected value for the highest volume fraction considered here. The reduction of the Reynolds stress in inner units indicates that the increase of the drag is not due to an enhancement of the turbulence activity, rather that it is linked to the solid phase dynamics.
Intensity of the fluctuation velocity components and Reynolds shear stress for the fluid phase in inner units for different volume fractions
${\it\Phi}$
: (a) streamwise
${u_{f}^{\prime }}_{rms}^{+}$
; (b) wall-normal
${v_{f}^{\prime }}_{rms}^{+}$
; and (c) spanwise
${w_{f}^{\prime }}_{rms}^{+}$
velocity component; (d) shear stress
$\langle u_{f}^{\prime }v_{f}^{\prime }\rangle ^{+}$
.
To analyse the solid phase behaviour, we report the mean local volume fraction
${\it\phi}(y)$
and the mean particle velocity
$U_{p}$
in figure 6. The mean local volume fraction (figure 6
a) shows a first local maximum around
$y=0.06{-}0.1$
, a value slightly larger than one particle radius (
$y=1/18$
). Increasing the bulk volume fraction
${\it\Phi}$
the intensity of the peak grows, while a local minimum appears at
$y\sim d_{p}=h/9$
. As also observed in dense laminar regimes (Yeo & Maxey Reference Yeo and Maxey2010a
), a particle layer forms at the wall and becomes more intense when increasing the bulk volume fraction
${\it\Phi}$
. It should be noted however that these near-wall maxima are smaller or similar to the bulk concentration, hence they are not related to the turbophoretic drift typically observed in dilute suspensions when particles are heavier than the fluid (Reeks Reference Reeks1983). Instead, these near-wall layers are induced by the planar symmetry of the wall and the excluded finite volume of the solid spheres. We believe that the formation of this particle layer follows a mechanics similar to that usually observed in laminar Poiseuille and Couette flows (Yeo & Maxey Reference Yeo and Maxey2010a
, Reference Yeo and Maxey2011; Picano et al.
Reference Picano, Breugem, Mitra and Brandt2013). Once a particle reaches the wall the strong wall–particle lubrication interaction stabilises the particle wall-normal position that is therefore mainly affected by the collisions with other particles. Hence, it becomes difficult for the particles belonging to the first layer to escape from it. Figure 6(b) depicts the mean particle velocity
$U_{p}^{+}$
in inner units (solid lines) where the fluid velocity is also reported with symbols for a close comparison. As shown in the figure, solid and fluid phases flow with the same mean velocity in the whole channel with the exception of the first particle layer near the wall,
$y^{+}\leqslant 20$
, where particles have a mean velocity larger than the surrounding fluid. It should be considered here that while the velocity at the wall is zero for the fluid, this is not the case for the solid phase as particles can have a relative tangential motion.
Particle average data for different nominal volume fractions
${\it\Phi}$
: (a) mean local volume fraction
${\it\phi}$
versus the wall-normal coordinate
$y$
(maximum statistical error
$\pm 0.01$
); (b) mean particle velocity profile,
$U_{p}^{+}=U_{p}/U_{\ast }$
in viscous units
$y^{+}=y/{\it\delta}_{\ast }$
(lines) (maximum statistical error
$\pm 0.25\,U_{\ast }$
). The mean fluid velocity
$U_{f}^{+}$
is also reported for comparison (symbols).
The fluctuation intensities, r.m.s., of the particle velocities are shown in figure 7(a–c), in inner units. The streamwise component
${u_{p}^{\prime }}_{rms}$
shows similar fluctuation levels for both phases and all
${\it\Phi}$
with some small differences close to the wall where the solid phase fluctuations do not vanish. Considering the three velocity components we generally observe that particles tend to fluctuate less than the fluid at the same position except for the region close to the wall. This behaviour is summarised in figure 7(d) where we display 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})$
. Apart from a thin region close to the wall, the fluid turbulent kinetic energy is higher than the energy of the solid phase by about
$10{-}20\,\%$
. The higher particle fluctuation level in the near-wall region, due to the absence of a no-slip condition at the wall, suggests that this is the cause of the near-wall enhancement of the fluid fluctuation level (compared with the single-phase flow) discussed above. One last remark concerns the local peak of the wall-normal particle velocity fluctuation close to wall. This maximum originates from particles that reach and leave the first layer at the wall. In this region the fluid velocity fluctuations increase with
${\it\Phi}$
, although the maximum for the solid phase decreases. This is not contradictory, as it just indicates that at small volume fractions the incoming/leaving particles are fewer, but faster; with increasing
${\it\Phi}$
, more particles enter and leave the first layer although at smaller velocity as it is more crowded.
Intensity of the fluctuation velocity components for the solid phase in inner units for the different volume fraction
${\it\Phi}$
studied (maximum statistical error
$\pm 0.06\,U_{\ast }$
). (a) streamwise
${u_{p}^{\prime }}_{rms}^{+}$
; (b) wall-normal
${v_{p}^{\prime }}_{rms}^{+}$
; and (c) spanwise
${w_{p}^{\prime }}_{rms}^{+}$
component. Symbols represent the fluctuation levels of the fluid phase. (d) Displays the wall-normal profile of the ratio between the turbulent kinetic energy of the fluid and of the solid phase.
Figure 8 reports the mean particle angular velocity
${\it\Omega}_{z}$
, panel (a), and the particle angular velocity fluctuation r.m.s. in the spanwise
${\it\omega}_{z\,rms}^{\prime }$
, panel (b), streamwise
${\it\omega}_{x\,rms}^{\prime }$
, panel (c) and wallnormal
${\it\omega}_{y\,rms}^{\prime }$
, panel (d), directions. The mean particle angular velocity
${\it\Omega}_{z}$
is maximum close to the wall and vanishes in the centerline for symmetry. This behavior indicates that the particle belonging to the layer close to the wall tend to roll on the wall minimising their local slip velocity, which as previously discussed is in principle not vanishing. The slight reduction of the maximum rotation observed when increasing the volume fraction
${\it\Phi}$
is induced by the more intense particle–particle interactions occurring in the first layer. Interestingly, at
${\it\Phi}=0.2$
, in the bulk of the flow, the mean angular velocity is higher than in the other cases. This can be explained by the higher fluid velocity gradient exhibited in this region at
${\it\Phi}=0.2$
, see for instance figure 3. Concerning the fluctuation levels of the particle angular velocity, we note that the maximum of each component occurs near the wall showing values that are about
$15{-}25\,\%$
of the mean value. Near the peak, the spanwise fluctuating component shows higher intensity
${\it\omega}_{z\,rms}^{\prime }$
driven by the inhomogeneity of the mean angular velocity, while the three components become of similar magnitude near the centreline (isotropy). The densest case shows slightly smaller fluctuations whereas the flows at
${\it\Phi}=0.05$
and
$0.1$
exhibit almost the same values.
Particle angular velocity statistics in outer units for the different volume fractions
${\it\Phi}$
studied: (a) mean angular velocity (spanwise component)
${\it\Omega}_{z}$
; (b) fluctuation r.m.s. of the spanwise fluctuation component
${\it\omega}_{z\,rms}^{\prime }$
; (c) fluctuation r.m.s. of the streamwise component
${\it\omega}_{x\,rms}^{\prime }$
; (d) fluctuation r.m.s. of the wall-normal component
${\it\omega}_{y\,rms}^{\prime }$
.
Momentum budget for the different bulk volume fractions
${\it\Phi}$
under investigation. The wall is at
$y=0$
, whereas
$y=1$
is the channel centreline. Here
${\it\tau}_{V}$
,
${\it\tau}_{T}$
and
${\it\tau}_{P}$
represent the viscous, turbulent and particle induced stresses. Here
${\it\tau}_{T_{p}}$
is the particle Reynolds stress and
${\it\tau}=U_{\ast }^{2}\,(1-y)$
the total stress, (a)
${\rm\Phi}=0$
, (b)
${\rm\Phi}=0.05$
, (c)
${\rm\Phi}=0.1$
and (d)
${\rm\Phi}=0.2$
.
(a) Wall-normal profiles of the shear stress of the combined phase
$\langle u_{C}^{\prime }v_{C}^{\prime }\rangle /U_{0}^{2}$
(symbols) and linear fitting of the slope of the profile at the centreline,
$y=1$
, depicted with solid lines. (b) Friction Reynolds number
$\mathit{Re}_{{\it\tau}}=u_{\ast }h/{\it\nu}$
and turbulent friction Reynolds number
$\mathit{Re}_{T}=U_{\ast }^{T}h/{\it\nu}$
versus the bulk volume fraction
${\it\Phi}$
from the simulations presented here.
Correlations of the velocity fluctuations versus the spanwise separation
${\rm\Delta}z^{+}$
for different bulk volume fractions
${\it\Phi}$
. Streamwise–streamwise component
$R_{uu}$
at (a)
$y=d=h/18\simeq 20{\it\delta}_{\ast }$
and (b) at
$y=2d=h/9\simeq 40{\it\delta}_{\ast }$
. Wall-normal component
$R_{vv}$
at (c)
$y=d=h/18\simeq 20{\it\delta}_{\ast }$
and (d) at
$y=2d=h/9\simeq 40{\it\delta}_{\ast }$
.
3.2. Total stress balance
The understanding of the momentum exchange between the two phases in dense particle-laden turbulent channel flows is conveniently addressed by examining the streamwise momentum budget, i.e. the average stress budget. Following the rationale on the mean momentum balance given in appendix A (see also Marchioro, Tankslay & Prosperetti (Reference Marchioro, Tankslay and Prosperetti1999) and Zhang & Prosperetti (Reference Zhang and Prosperetti2010) for more details), we can write the whole budget as the sum of three terms:
$$\begin{eqnarray}{\it\tau}={\it\tau}_{V}+{\it\tau}_{T}+{\it\tau}_{P},\end{eqnarray}$$
where
${\it\tau}=U_{\ast }^{2}\,(1-y)$
is the total stress,
${\it\tau}_{V}={\it\nu}(1-{\it\phi})(\text{d}U_{f})/(\text{d}y)$
is the viscous stress,
${\it\tau}_{T}=-\langle u_{c}^{\prime }v_{c}^{\prime }\rangle$
is the turbulent Reynolds shear stress of the combined phase
$\langle u_{c}^{\prime }v_{c}^{\prime }\rangle ={\it\phi}\langle u_{p}^{\prime }v_{p}^{\prime }\rangle +(1-{\it\phi})\langle u_{f}^{\prime }v_{f}^{\prime }\rangle$
(with the particle Reynolds stress
${\it\phi}\langle u_{p}^{\prime }v_{p}^{\prime }\rangle ={\it\tau}_{T_{p}}$
) and
${\it\tau}_{P}=({\it\phi}/{\it\rho})(\langle {\it\sigma}_{p\,xy}\rangle )$
the particle induced stress.
Figure 9 reports the stress balance given in (3.2) from the simulations for the four bulk volume fractions
${\it\Phi}$
presented here and normalised by the corresponding friction velocity squared,
$U_{\ast }^{2}$
(the particle-induced stress has been indirectly calculated from the balance). As already known for the single-phase flow (Pope Reference Pope2000), the total stress
${\it\tau}$
is mainly given by the turbulent Reynolds stress term for
$y\geqslant 0.2$
. The relevance of the viscous stress increases approaching the wall, becoming the leading term as the Reynolds stress is zero at the wall. At
${\it\Phi}=0.05$
, see figure 9(b), the basic picture remains unaltered with the particle-induced stress
${\it\tau}_{P}$
showing a non-negligible contribution only near the wall; note that the particle turbulent Reynolds stress is still negligible in this configuration. Increasing the volume fraction to
${\it\Phi}=0.1$
(figure 9
c), the particle-induced stress becomes of the same order of magnitude as the other terms in the near wall region,
$y\approx 0.05$
, which roughly corresponds to a particle radius. The contribution from the particle stress, although still subleading with respect to the turbulent stress
${\it\tau}_{T}$
, is important throughout the whole channel. Note also that the turbulent stress associated to the solid phase alone,
${\it\tau}_{T_{p}}$
, amounts to
${\sim}10\,\%$
of the total
${\it\tau}_{T}$
, scaling almost linearly with the volume fraction. For the highest volume fraction considered,
${\it\Phi}=0.2$
(see figure 9
d), the near-wall dynamics is dominated by the particle-induced stresses. This is now the leading term around
$y=0.05$
. Moreover, the total stress in the bulk of the flow,
$y>0.2$
, is transmitted by the turbulent shear stress
${\it\tau}_{T}$
and the particle induced stress
${\it\tau}_{P}$
in similar shares. In other words, the turbulent shear stress amounts to about half of the total stress in the bulk of the flow. This indicates that the turbulent dynamics is strongly altered by the dense particle concentration: although the system is still turbulent, the particle-induced stress becomes crucial in transferring the mean stress through the channel.
This behaviour is consistent with the decrease of the turbulence activity discussed previously for the flow with the highest particle number,
${\it\Phi}=0.2$
. As mentioned above, although the turbulence intensities and the Reynolds shear stress are attenuated, the total drag, i.e. the friction Reynolds number increases. One can therefore conclude that this increase of the total drag is not associated to a turbulence enhancement, but to an increase of the particle-induced stress, or borrowing rheological terms, to an increase of the effective viscosity of the flowing medium.
In order to quantify the level of turbulence activity, we can define the turbulent friction velocity as
$$\begin{eqnarray}U_{\ast }^{T}=\left.\sqrt{\frac{\text{d}\langle u_{c}^{\prime }v_{c}^{\prime }\rangle }{\text{d}y}}\right|_{y=1},\end{eqnarray}$$
that is the square root of the wall-normal derivative of the Reynolds stress profile at the centreline (
$y=1$
). This quantity has been chosen because it can be shown that the turbulent friction velocity well approximates the wall friction velocity for unladen cases at high bulk Reynolds number,
$U_{\ast }^{T}=U_{\ast }+O(1/\mathit{Re})$
, see e.g. Pope (Reference Pope2000).
Figure 10(a) reports the turbulent Reynolds stress of the combined phase
$\langle u_{c}^{\prime }v_{c}^{\prime }\rangle$
in outer units together with a straight line indicating the slope at
$y=1$
. The intercept of this line originating at
$(y,\langle u_{c}^{\prime }v_{c}^{\prime }\rangle )=(1,0)$
with the vertical axis provides the value of the turbulent friction velocity,
$U_{\ast }^{T}$
as defined above. As clear from the figure,
$U_{\ast }^{T}$
increases when adding the solid phase until
${\it\Phi}=0.1$
and then decreases at
${\it\Phi}=0.2$
. Using these values, we can then define a turbulent friction Reynolds number:
$\mathit{Re}_{T}=U_{\ast }^{T}\,h/{\it\nu}$
. Since
$\mathit{Re}_{T}$
is proportional to
$U_{\ast }^{T}$
, it follows that
$\mathit{Re}_{T}=\mathit{Re}_{{\it\tau}}+O(1/\mathit{Re})$
for high-Reynolds-number single-phase turbulent channel flows. Figure 10(b) depicts the friction Reynolds number
$\mathit{Re}_{{\it\tau}}$
and the turbulent Reynolds number
$\mathit{Re}_{T}$
just introduced versus the bulk volume fraction
${\it\Phi}$
. The values of the two Reynolds numbers in the unladen case,
${\it\Phi}=0$
, are close, as expected. Increasing the volume fraction, both
$\mathit{Re}_{{\it\tau}}$
and
$\mathit{Re}_{T}$
increase up to
${\it\Phi}=0.1$
,
$\mathit{Re}_{T}$
at a slower rate. Interestingly, at
${\it\Phi}=0.2$
, the turbulent friction Reynolds number
$\mathit{Re}_{T}$
suddenly decreases, whereas the friction Reynolds number based on the actual wall-shear still increases.
The friction velocity and Reynolds number are a measure of the overall drag as they are proportional to the imposed pressure gradient, while the turbulent friction velocity and corresponding Reynolds number introduced here indicate only the portion of the drag directly induced by the turbulent activity. We therefore conclude that in dense cases, i.e.
${\it\Phi}=0.2$
, a turbulent drag reduction indeed occurs and this is related to a reduced turbulence activity. Nonetheless, this turbulent drag reduction does not reflect in a decrease of the total drag at the Reynolds number investigated here because the particle-induced stress (increased viscosity of the suspension) more than counteracts the positive effect due to the reduced turbulent mixing. The observations emerging from our analysis of the momentum budget explain and are consistent with the large reduction of the von Kármán constant
$k$
found for this dense case and reported in table 2.
3.3. Velocity correlations
Further understanding of the effect of the solid phase on the turbulent channel flow is obtained by examining the two-point spatial correlation of the velocity field. It is well known that the autocorrelations of the streamwise and wall-normal velocity along the spanwise direction,
$$\begin{eqnarray}\displaystyle & \displaystyle R_{uu}(y,{\rm\Delta}z)=\frac{\langle u^{\prime }(x,y,z,t)u^{\prime }(x,y,z+{\rm\Delta}z,t)\rangle }{{u^{\prime }}_{rms}^{2}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle R_{vv}(y,{\rm\Delta}z)=\frac{\langle v^{\prime }(x,y,z,t)v^{\prime }(x,y,z+{\rm\Delta}z,t)\rangle }{{v^{\prime }}_{rms}^{2}} & \displaystyle\end{eqnarray}$$
${\rm\Delta}z^{+}=60{-}80$
and
${\rm\Delta}z^{+}=30{-}40$
, respectively, for a single-phase turbulent flow.
These values reflect the typical structures of wall-bounded turbulence, i.e. quasi-streamwise vortices and low-speed streaks that sustain the turbulence process (Kim, Moin & Moser Reference Kim, Moin and Moser1987; Waleffe Reference Waleffe1997; Pope Reference Pope2000; Brandt Reference Brandt2014). It has also been observed that in drag-reducing turbulent flows the width and the spacing of these characteristic structures increases (De Angelis, Casciola & Piva Reference De Angelis, Casciola and Piva2002; Stone, Waleffe & Graham Reference Stone, Waleffe and Graham2002), leading to an increase of the spanwise separation of these minima.
The streamwise autocorrelation
$R_{uu}$
is shown in figure 11(a,b), where it is evaluated at two wall-normal distances,
$y=d\simeq 20{\it\delta}_{\ast }$
and
$y=2d\simeq 40{\it\delta}_{\ast }$
. The correlations are here calculated for the combined phase, but they do not differ appreciably if calculated only for the fluid phase. At
$y=d$
we note a progressive increase of the separation distance with the particle volume fraction together with a smoothening of the minimum, indicating a less-evident width of the near-wall flow structures. Further away from the wall,
$y=2d$
, we observe the formation of wider streamwise velocity streaks for the flow with
${\it\Phi}=0.2$
, with a separation of the minimum of the autocorrelation,
${\rm\Delta}z^{+}$
, that is almost twice that pertaining to single-phase near-wall turbulence. The system tends therefore to form streaks twice as large as those in single-phase turbulent channel flows. These larger structures are also seen by the shift of the lowest minimum of
$R_{uu}$
to
$y=2d$
instead of
$y=d\sim 20{\it\delta}_{\ast }$
where the single-phase channel flow shows the sharper minimum in the autocorrelation functions. The wall-normal autocorrelations
$R_{vv}$
are shown in figure 11(c,d) for the same two wall-parallel planes. Increasing the volume fraction
${\it\Phi}$
we observe less sharp minima that completely disappear at
${\it\Phi}=0.2$
. This suggests a significant alteration of the structure of the wall turbulence at high volume fractions with a flow much less organised in coherent structures. Similar observations are reported in Loisel et al. (Reference Loisel, Abbas, Masbernat and Climent2013) for transitional flows at lower
${\it\Phi}$
. The behaviour of the velocity autocorrelations is consistent with what found in turbulent drag-reducing flows (growth of the buffer region). Hence, it appears once more that despite the total drag increase, the turbulent induced drag reduces at least at high volume fraction.
4. Final remarks
We report data from the numerical simulations of turbulent channel flow laden with finite-size particles at high volume fractions. The simulations have been performed using an efficient implementation of the IBM that enable us to fully resolve the fluid–structure interactions. We provide a statistical analysis to assess the effect of an increasing solid volume fraction (up to
${\it\Phi}=0.2$
) on a turbulent channel flow at fixed bulk Reynolds number, i.e.
$\mathit{Re}_{b}=U_{0}\,2h/{\it\nu}$
.
The finite-size particles interact with the turbulent motions altering the near-wall turbulence regeneration process. For the two lowest volume fractions considered,
${\it\Phi}\leqslant 0.1$
, we still observe the classic behaviour of near-wall turbulence, modulated however by the particle presence. At
${\it\Phi}=0.2$
the solid phase is so dense that several aspects of turbulent wall flows are lost: the mean velocity profile is strongly altered, the turbulent fluctuations decrease, the velocity autocorrelations show streamwise elongated structures twice as wide as in single-phase channel flows and the absence of a negative correlation of the wall-normal velocity, in addition to a more isotropic distribution of the velocity fluctuations.
The law of the wall is modified by the presence of a solid phase but can still be recognised at the Reynolds number of our simulations for all of the volume fractions investigated. The von Kármán and additive constants,
$k$
and
$B$
, therefore assume different values. In particular, increasing the volume fraction we report a reduction of
$k$
, increase of the slope, and a strong decrease of
$B$
, increased near-wall dissipation. The reduction of
$k$
usually denotes turbulent drag reduction. However, in the present cases we always observe an increase of the overall drag due to the decrease of the additive constant
$B$
. This is also confirmed by the increase of the friction Reynolds number,
$\mathit{Re}_{{\it\tau}}$
, when increasing the volume fraction at constant mass flux.
We evaluate the streamwise momentum balance for the flows under investigation and show that the additional stress due to the presence of the particles becomes increasingly relevant when increasing the particle volume fraction. As expected the Reynolds transport term dominates at zero and low
${\it\Phi}$
, while at
${\it\Phi}=0.2$
the particle stress becomes of the same order of magnitude.
Examining the turbulent shear stress and the streamwise momentum balance, we thus note that the turbulence activity and the related stress reduce at the highest volume fraction considered here, i.e.
${\it\Phi}=0.2$
. In order to characterise the turbulent drag, we define a turbulent friction Reynolds number
$\mathit{Re}_{T}$
whose friction velocity is based on the slope of the Reynolds shear stress profile at the centreline. This parameter approximates the usual
$\mathit{Re}_{{\it\tau}}$
in unladen turbulent channel flow. Using this turbulent friction Reynolds number, we quantitatively show that the turbulent drag (measured by
$\mathit{Re}_{T}$
) first gently increases with
${\it\Phi}$
and then sharply decreases at
${\it\Phi}=0.2$
, even though the overall drag still increases.
These results suggest that further increasing the Reynolds number while keeping constant the particle size in inner units,
$d^{+}$
, may lead to an overall drag reduction in dense cases as shown here. The main assumption behind this conjecture is that the near-wall turbulence–particle dynamics remain similar when the bulk Reynolds number is increased, as might occur when the particle size in inner units remains constant (i.e. the friction particle Reynolds number). Indeed, we show here that increasing the bulk Reynolds number the turbulent induced drag increases its weight in the stress balance. Hence, the reduced turbulence activity and the consequent reduced turbulent drag should induce a decrease of the total drag at high enough Reynolds. This should appear as an extension of the log-layer with almost the same reduced
$k$
and
$B$
as reported here. New and even larger simulations would be needed in the future to test this hypothesis. In the meanwhile, we hope to stimulate new experimental investigations towards this direction.
This study reports detailed statistics of particle-laden channel flow at high volume fractions, accessible only recently (Kidanemariam & Uhlmann Reference Kidanemariam and Uhlmann2014; Vowinckel et al. Reference Vowinckel, Kempe and Fröhlich2014), and it could therefore be extended in many non-trivial directions. Two-body particle statistics, such as collisions rates and clustering, have not yet been considered because it is out of the scope of the present work. In addition, the effect of the particle shape (Bellani et al. Reference Bellani, Byron, Collignon, Meyer and Variano2012) and deformability (e.g. Clausen, Reasor & Aidun Reference Clausen, Reasor and Aidun2011) surely deserves attention as it will add new interesting physics to our current understanding.
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 the support from the COST Action MP1305: Flowing matter.
Appendix A. Total stress of the suspension mixture
In this work we use the framework developed by Prosperetti and coworkers to examine the stresses in suspension mixtures, see e.g. Marchioro et al. (Reference Marchioro, Tankslay and Prosperetti1999) and Zhang & Prosperetti (Reference Zhang and Prosperetti2010) for more details.
We assume the same density
${\it\rho}$
for the fluid and the particles and consider dimensional variables for all of the calculations presented in this appendix. Following Zhang & Prosperetti (Reference Zhang and Prosperetti2010), we define the phase indicator as
${\it\xi}=0$
in the fluid phase and
${\it\xi}=1$
in the solid phase. Defining the phase-ensemble average, ‘
$\langle \,\rangle$
’, as the ensemble average (implicitly) conditioned to the phase considered (particulate, fluid and combined), we can calculate the local volume fraction in a point as
$$\begin{eqnarray}{\it\phi}=\langle {\it\xi}\rangle .\end{eqnarray}$$
Considering a generic observable of the combined phase
$o_{c}={\it\xi}o_{p}+(1-{\it\xi})o_{f}$
, constructed in terms of
$o_{p/f}$
, the same observable in the particulate and fluid phases, it holds that
$$\begin{eqnarray}\langle o_{c}\rangle =\langle {\it\xi}o_{p}\rangle +\langle (1-{\it\xi})o_{f}\rangle ={\it\phi}\langle o_{p}\rangle +(1-{\it\phi})\langle o_{f}\rangle .\end{eqnarray}$$
Note that we are not using different symbols for the different phase ensemble averages, but implicitly assume that the phase conditioning is indicated by the subscript inside the brackets.
The force balance for the volume
$\mathscr{V}$
delimited by the surface
$\mathscr{S}(\mathscr{V})$
is
$$\begin{eqnarray}{\it\rho}\int _{\mathscr{V}}{\it\xi}\boldsymbol{a}_{p}+(1-{\it\xi})\boldsymbol{a}_{f}\text{d}V=\oint _{\mathscr{S}(\mathscr{V})}[{\it\xi}{\bf\sigma}_{p}+(1-{\it\xi}){\bf\sigma}_{f}]\boldsymbol{\cdot }\boldsymbol{n}\text{d}S,\end{eqnarray}$$
with
$\boldsymbol{n}$
the outer unity vector normal to the surface
$\mathscr{S}(\mathscr{V})$
, the subscripts ‘
$f$
’ and ‘
$p$
’ denoting fluid and particle phases,
$\boldsymbol{a}_{i}$
and
${\bf\sigma}_{i}$
the acceleration and the general stress in the phase
$i$
. Applying the phase ensemble average to (A 3), we obtain
$$\begin{eqnarray}{\it\rho}\int _{\mathscr{V}}\langle {\it\xi}\boldsymbol{a}_{p}\rangle +\langle (1-{\it\xi})\boldsymbol{a}_{f}\rangle \text{d}V=\int _{\mathscr{V}}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }[\langle {\it\xi}{\bf\sigma}_{p}\rangle +\langle (1-{\it\xi}){\bf\sigma}_{f}\rangle ]\text{d}V,\end{eqnarray}$$
where we used the divergence theorem to the differentiable integrand on the right-hand side. Since the last equation holds for any mesoscale volume
$\mathscr{V}$
, we can use the corresponding differential form of the equation,
$$\begin{eqnarray}{\it\rho}\langle {\it\xi}\boldsymbol{a}_{p}\rangle +{\it\rho}\langle (1-{\it\xi})\boldsymbol{a}_{f}\rangle =\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }[\langle {\it\xi}{\bf\sigma}_{p}\rangle +\langle (1-{\it\xi}){\bf\sigma}_{f}\rangle ].\end{eqnarray}$$
Considering the identities (A 1) and (A 2), we can further simplify the expression above
$$\begin{eqnarray}{\it\rho}{\it\phi}\langle \boldsymbol{a}_{p}\rangle +{\it\rho}(1-{\it\phi})\langle \boldsymbol{a}_{f}\rangle =\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\phi}\langle {\bf\sigma}_{p}\rangle +(1-{\it\phi})\langle {\bf\sigma}_{f}\rangle ).\end{eqnarray}$$
Assuming the constitutive law of a Newtonian fluid
${\bf\sigma}_{f}=-p\boldsymbol{I}+2{\it\mu}\unicode[STIX]{x1D640}$
with
$p$
the pressure and
$\unicode[STIX]{x1D640}=(\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{f}+\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{f}^{\text{T}})/2$
the symmetric part of the fluid velocity gradient tensor and considering that both the fluid and particle velocity fields are divergence-free, equation (A 6) can be rewritten as
$$\begin{eqnarray}\displaystyle & & \displaystyle {\it\phi}\frac{\langle \boldsymbol{u}_{p}\rangle }{\partial t}+{\it\phi}\langle \boldsymbol{u}_{p}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{p}\rangle +(1-{\it\phi})\frac{\langle \boldsymbol{u}_{f}\rangle }{\partial t}+(1-{\it\phi})\langle \boldsymbol{u}_{f}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{f}\rangle \nonumber\\ \displaystyle & & \displaystyle \quad =\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\phi}\langle {\bf\sigma}_{p}/{\it\rho}\rangle )-\boldsymbol{{\rm\nabla}}[(1-{\it\phi})\langle p/{\it\rho}\rangle ]+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }[(1-{\it\phi})2{\it\nu}\langle \unicode[STIX]{x1D640}\rangle ].\end{eqnarray}$$
We next denote the statistically stationary mean fluid and particle velocities as
$\boldsymbol{U}_{f/p}=\langle \boldsymbol{u}_{f/p}\rangle$
and the fluctuations around these mean values as
$\boldsymbol{u}_{f/p}^{\prime }=\boldsymbol{U}_{f/p}-\langle \boldsymbol{u}_{f/p}\rangle$
, so that the average momentum equation becomes
$$\begin{eqnarray}\displaystyle & & \displaystyle {\it\phi}\boldsymbol{U}_{p}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{U}_{p}+{\it\phi}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\langle \boldsymbol{u}_{p}^{\prime }\boldsymbol{u}_{p}^{\prime }\rangle +(1-{\it\phi})\boldsymbol{U}_{f}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{U}_{f}+(1-{\it\phi})\langle \boldsymbol{u}_{f}^{\prime }\boldsymbol{u}_{f}^{\prime }\rangle \nonumber\\ \displaystyle & & \displaystyle \quad =\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }({\it\phi}\langle {\bf\sigma}_{p}/{\it\rho}\rangle )-\boldsymbol{{\rm\nabla}}[(1-{\it\phi})P/{\it\rho}]+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }[(1-{\it\phi})2{\it\nu}\langle \unicode[STIX]{x1D640}\rangle ],\end{eqnarray}$$
with
$P$
the mean pressure.
Exploiting the symmetries of a fully developed parallel channel flow, characterised by two homogeneous directions, the streamwise
$x$
and spanwise
$z$
, we project (A 8) in the inhomogeneous wall-normal direction
$y$
,
$$\begin{eqnarray}\frac{\text{d}}{\text{d}y}\left[(1-{\it\phi})\langle v_{f}^{\prime 2}\rangle +{\it\phi}\langle v_{p}^{\prime 2}\rangle +(1-{\it\phi})\frac{P}{{\it\rho}}-\frac{{\it\phi}}{{\it\rho}}\langle {\it\sigma}_{p\,yy}\rangle \right]=0.\end{eqnarray}$$
Integrating (A 9) in the
$y$
direction and denoting the wall pressure by
$P_{w}(x)$
, we obtain
$$\begin{eqnarray}(1-{\it\phi})\langle v_{f}^{\prime 2}\rangle +{\it\phi}\langle v_{p}^{\prime 2}\rangle +\frac{P_{T}}{{\it\rho}}=\frac{P_{w}}{{\it\rho}},\end{eqnarray}$$
where we also introduced the mean total pressure
$P_{T}=(1-{\it\phi})(P/{\it\rho})-{\it\phi}\langle {\it\sigma}_{P\,yy}\rangle /{\it\rho}$
. It should be noted that
$P_{T}$
coincides with
$P_{w}$
at the wall and that
$$\begin{eqnarray}\frac{\partial P_{T}}{\partial x}=\frac{\partial P_{w}}{\partial x}.\end{eqnarray}$$
Projecting (A 8) in the streamwise direction
$x$
, we have
$$\begin{eqnarray}\frac{\text{d}}{\text{d}y}\left[(1-{\it\phi})\langle u_{f}^{\prime }v_{f}^{\prime }\rangle +{\it\phi}\langle u_{p}^{\prime }v_{p}^{\prime }\rangle -{\it\nu}(1-{\it\phi})\frac{\text{d}U_{f}}{\text{d}y}-\frac{{\it\phi}}{{\it\rho}}\langle {\it\sigma}_{p\,xy}\rangle \right]=-\frac{\text{d}}{\text{d}x}\left(\frac{P_{w}}{{\it\rho}}\right),\end{eqnarray}$$
where we neglect the terms
$(\partial /\partial x)[({\it\phi}/{\it\rho})(\langle {\it\sigma}_{p\,xx}-{\it\sigma}_{p\,yy}\rangle )]$
because of the streamwise homogeneity.
Integrating (A 12) in the wall normal direction and denoting the Reynolds shear stress of the combined phase
$\langle u_{C}^{\prime }v_{C}^{\prime }\rangle =(1-{\it\phi})\langle u_{f}^{\prime }v_{f}^{\prime }\rangle +{\it\phi}\langle u_{p}^{\prime }v_{p}^{\prime }\rangle$
, we obtain the equation for the total stress
${\it\tau}(y)$
,
$$\begin{eqnarray}{\it\tau}(y)=-\langle u_{C}^{\prime }v_{C}^{\prime }\rangle +{\it\nu}(1-{\it\phi})\frac{\text{d}U_{f}}{\text{d}y}+\frac{{\it\phi}}{{\it\rho}}\langle {\it\sigma}_{p\,xy}\rangle =\left.{\it\nu}\frac{\text{d}U_{f}}{\text{d}y}\right|_{w}\left(1-\frac{y}{h}\right),\end{eqnarray}$$
where we considered the boundary condition at the wall,
${\it\tau}_{w}={\it\tau}(0)={\it\nu}(\text{d}U_{f}/\text{d}y)|_{y=0}$
. Equation (A 13) shows that the total stress of a turbulent suspension in a channel geometry is given by three contributions: the viscous part,
${\it\tau}_{V}={\it\nu}(1-{\it\phi})(\text{d}U_{f}/\text{d}y)$
, the turbulent part
${\it\tau}_{T}=-\langle u_{C}^{\prime }v_{C}^{\prime }\rangle =-(1-{\it\phi})\langle u_{f}^{\prime }v_{f}^{\prime }\rangle -{\it\phi}\langle u_{p}^{\prime }v_{p}^{\prime }\rangle$
and the particle-induced stress,
${\it\tau}_{P}={\it\phi}\langle {\it\sigma}_{p\,xy}\rangle /{\it\rho}$
. It should be noted that the turbulent stress accounts for the coherent streamwise and wall-normal motion of both fluid and solid phases. The particle-induced stress is originated by the total stress exerted by the solid phase, see (A 3), and takes into account hydrodynamic interactions and collisions. In the absence of particles,
${\it\phi}\rightarrow 0$
, (A 13) reduces to the classic momentum balance for single-phase turbulence (see Pope Reference Pope2000).
