Effects of the finite particle size in turbulent wall-bounded flows of dense suspensions

TL;DR

This study uses detailed simulations to analyze how particle size influences turbulence and particle behavior in dense suspensions within wall-bounded flows, revealing size-dependent effects on flow statistics and particle dynamics.

Contribution

It provides new insights into finite-size effects on turbulence and particle dynamics in dense suspensions, especially near walls, extending existing flow scaling laws.

Findings

Particle-wall layer influences suspension micro-structure over a distance proportional to particle size.

Significant slip velocity and extreme shear stress events occur in the particle-wall layer.

Finite-size effects are more pronounced for larger particles in the bulk flow.

Abstract

We use interface-resolved simulations to study finite-size effects in turbulent channel flow of neutrally-buoyant spheres. Two cases with particle sizes differing by a factor of 2, at the same solid volume fraction of 20% and bulk Reynolds number are considered. These are complemented with two reference single-phase flows: the unladen case, and the flow of a Newtonian fluid with the effective suspension viscosity of the same mixture in the laminar regime. As recently highlighted in Costa et al. (PRL 117, 134501), a particle-wall layer is responsible for deviations of the statistics from what is observed in the continuum limit where the suspension is modeled as a Newtonian fluid with an effective viscosity. Here we investigate the fluid and particle dynamics in this layer and in the bulk. In the particle-wall layer, the near wall inhomogeneity has an influence on the suspension micro-structure over a distance proportional to the particle size. In this layer, particles have a significant (apparent) slip velocity that is reflected in the distribution of wall shear stresses. This is characterized by extreme events (both much higher and much lower than the mean). Based on these observations we provide a scaling for the particle-to-fluid apparent slip velocity as a function of the flow parameters. We also extend the flow scaling laws in to second-order Eulerian statistics in the homogeneous suspension region away from the wall. Finite-size effects in the bulk of the channel become important for larger particles, while negligible for lower-order statistics and smaller particles. Finally, we study the particle dynamics along the wall-normal direction. Our results suggest that 1-point dispersion is dominated by particle-turbulence (and not particle-particle) interactions, while differences in 2-point dispersion and collisional dynamics are consistent with a picture of shear-driven interactions.

Figures

Tables

Table 1. Table 1: Parameters for lubrication and contact models for short-range particle-particle and particle-wall interactions.

interaction	$e_{n,d}$	${\epsilon }_{\Delta x}$	${\epsilon }_{\sigma }$	$N$
particle-wall	$0.97$	$0.075$	$0.001$	$10$
particle-particle	$0.97$	$0.025$	$0.001$	$10$

Table 2. Table 2: Physical and computational parameters of the DNS data sets. δ v s ​ p ​ h superscript subscript 𝛿 𝑣 𝑠 𝑝 ℎ \delta_{v}^{sph} ( ≳ δ v greater-than-or-equivalent-to absent subscript 𝛿 𝑣 \gtrsim\delta_{v} ) denotes the viscous wall unit for the corresponding single-phase flow at the same \Rey b subscript \Rey 𝑏 \Rey_{b} , estimated from \Rey τ s ​ p ​ h = 0.09 ​ \Rey b 0.88 superscript subscript \Rey 𝜏 𝑠 𝑝 ℎ 0.09 superscript subscript \Rey 𝑏 0.88 \Rey_{\tau}^{sph}=0.09\Rey_{b}^{0.88} (Pope, 2001 ) . N p subscript 𝑁 𝑝 N_{p} denotes the number of particles, N x / y / z subscript 𝑁 𝑥 𝑦 𝑧 N_{x/y/z} denotes the number of grid points in the x / y / z 𝑥 𝑦 𝑧 x/y/z directions.

Case	${\text{Rey}}_b$	$h/D_p$	$D_p/{\delta }_v^{sph}$	$\Phi (\%)$	$N_x\times N_y\times N_z$	$N_p$
D10	$\mathrm{12\hspace{0.17em}000}$	$36$	$9.7$	$20$	$3456\times 1152\times 1728$	$\mathrm{640\hspace{0.17em}000}$
D20	$\mathrm{12\hspace{0.17em}000}$	$18$	$19.4$	$20$	$1728\times 576\times 864$	$\mathrm{80\hspace{0.17em}000}$
SPR	$\mathrm{12\hspace{0.17em}000}$	–	–	–	$1728\times 576\times 864$	$0$
CLR	$\mathrm{6\hspace{0.17em}400}$	–	–	–	$1152\times 384\times 576$	$0$

Table 3. Table 3: Friction Reynolds numbers for the different cases studied. \Rey τ = u τ ​ h / ν subscript \Rey 𝜏 subscript 𝑢 𝜏 ℎ 𝜈 \Rey_{\tau}=u_{\tau}h/\nu is the typical friction Reynolds number defined from the wall friction velocity and fluid viscosity, \Rey τ e = \Rey τ ​ ν / ν e superscript subscript \Rey 𝜏 𝑒 subscript \Rey 𝜏 𝜈 superscript 𝜈 𝑒 \Rey_{\tau}^{e}=\Rey_{\tau}\nu/\nu^{e} is instead defined from the suspension effective viscosity and \Rey τ e ⁣ ∗ = u τ ​ ( 1 − δ p ​ w / h ) 1 / 2 ​ ( h − δ p ​ w ) / ν e = \Rey τ e ​ ( 1 − δ p ​ w / h ) 3 / 2 superscript subscript \Rey 𝜏 𝑒 subscript 𝑢 𝜏 superscript 1 subscript 𝛿 𝑝 𝑤 ℎ 1 2 ℎ subscript 𝛿 𝑝 𝑤 superscript 𝜈 𝑒 superscript subscript \Rey 𝜏 𝑒 superscript 1 subscript 𝛿 𝑝 𝑤 ℎ 3 2 \Rey_{\tau}^{e*}=u_{\tau}(1-\delta_{pw}/h)^{1/2}(h-\delta_{pw})/\nu^{e}=\Rey_{\tau}^{e}(1-\delta_{pw}/h)^{3/2} is defined from a wall friction velocity and channel height corrected for finite-size effects and the effective suspension viscosity.

Case	${\text{Rey}}_{\tau }$	${\text{Rey}}_{\tau }^e$	${\text{Rey}}_{\tau }^{e\ast }$
SPR (unladen)	$350$	$350$	$350$
D10	$395$	$209$	$203$
D20	$406$	$215$	$203$
CLR	$201$	$201$	$201$

Equations

$$\nabla \cdot \mathbf{u}$$

$$\frac{\mathrm{D}\mathbf{u}}{\mathrm{D}t}$$

$$m_p\frac{\mathrm{d}\mathbf{U}}{\mathrm{d}t}$$

$$I_p\frac{\mathrm{d}\mathbf{\Omega }}{\mathrm{d}t}$$

$$\mathbf{U}+\mathbf{\Omega }\times \mathbf{r}=\mathbf{u}\forall \mathbf{x}\in \partial V.$$

$${\mathbf{F}}_{ij,n}=-k_n{\mathbf{\delta }}_{ij,n}-{\eta }_n{\mathbf{u}}_{ij,n}$$

$$k_n=\frac{m_e\left({\pi }^2+{\ln }^2{e}_{n,d}\right)}{(N\Delta t{)}^2},{\eta }_n=-\frac{2m_e\ln e_{n,d}}{(N\Delta t)},$$

$${⟨o⟩}^{f/s}(y)=\frac{{\sum }_{ik}{o}_{ijk}(x,y,z){\varphi }_{ijk}^{f/s}(x,y,z)}{{\sum }_{ik}{\varphi }_{ijk}^{f/s}(x,y,z)}$$

$${⟨o⟩}^c(y)=\frac{{\sum }_{ik}\left(o_{ijk}(x,y,z){\varphi }_{ijk}^f(x,y,z)+o_{ijk}(x,y,z){\varphi }_{ijk}^s(x,y,z)\right)}{{\sum }_{ik}\left({\varphi }_{ijk}^f(x,y,z)+{\varphi }_{ijk}^s(x,y,z)\right)}.$$

$$\frac{u}{u_{\tau }^{\ast }}=$$

$$\frac{U_c-u}{u_{\tau }^{\ast }}=-$$

$${\text{Rey}}_{\tau }=\frac{{\text{Rey}}_b}{2{\xi }_{pw}^{1/2}}{\left(\frac{1}{\kappa }\left[\ln \left({\text{Rey}}_{\tau }{\chi }^e{\xi }_{pw}^{3/2}\right)-1\right]+B+B_d\right)}^{-1},$$

$${\text{Rey}}_{\tau }=\frac{0.09{\left({\text{Rey}}_b{\chi }^e{\xi }_{pw}\right)}^{0.88}}{{\xi }_{pw}^{3/2}{\chi }^e},$$

$$g(\theta ,\phi )=\frac{1}{r^2\sin (\theta )}\frac{{\mathrm{d}}^2{N}_{\theta ,\phi }}{\mathrm{d}\theta \mathrm{d}\phi }\frac{1}{n_0};n_0=\frac{N(N-1)}{2V}$$

$$N_{pw}{D}_p{L}_w\mu \frac{U_{pw}}{D_p}=\rho u_{\tau }^2{L}_x{L}_z.$$

$$\Phi \frac{L_w}{D_p}\frac{U_{pw}}{D_p}\sim \frac{u_{\tau }^2}{\nu },$$

$$\frac{U_{pw}}{u_{\tau }}=C_{pw}R{e}_{\tau }^p{\Phi }^{-1/2},$$

$$⟨\tau ⟩=u_{\tau }^2\left(1-\frac{y}{h}\right)=\underset{{\tau }_{T_f}}{\underbrace{(1-\varphi )⟨-u_f^{\prime }{v}_f^{\prime }⟩}}+\underset{{\tau }_v}{\underbrace{(1-\varphi )\nu \frac{\mathrm{d}{U}_f}{\mathrm{d}y}}}+\underset{{\tau }_{T_p}}{\underbrace{\varphi ⟨-u_p^{\prime }{v}_p^{\prime }⟩}}+{\tau }_p,$$

$$⟨\Delta z_1^2⟩(\Delta t,y)=⟨{\left(z(t+\Delta t)-z(t)\right)}^2⟩$$

$$⟨\Delta z_2^2⟩(\Delta t,y)=⟨{\left(\delta z(t_0+\Delta t)-\delta z(t_0)\right)}^2⟩,$$

$$g(r,y)=\frac{1}{4\pi r^2}\frac{\mathrm{d}{N}_r}{\mathrm{d}r}\frac{1}{n_0},$$

$$\Delta v^{n,-}(r,y)=\max \left\{0,-({\mathbf{u}}_j-{\mathbf{u}}_i)\cdot \frac{{\mathbf{x}}_j-{\mathbf{x}}_i}{|{\mathbf{x}}_j-{\mathbf{x}}_i|}\right\};$$

$$\Delta v^{n,-}\propto D_p\frac{\partial u}{\partial y}\sim D_p\frac{u_{\tau }}{h},$$

Full text

1 Introduction

Turbulent transport of solid particles in a suspending Newtonian fluid is often found in natural and industrial contexts. Few of many well-known examples are sediment transport in rivers, sand storms, slurries, and the flocculation and sedimentation processes in the treatment of drinking water. In these examples and many others the flow can be dense, i.e., the solid volume fraction sufficiently high that short-range particle-particle interactions are dynamically significant. Moreover, the particles often have a finite size, i.e. comparable to or larger than the smallest relevant length scale of the flow.

The suspended particles may vary in size, shape, density, mechanical properties, etc. Moreover, the fluid can be non-Newtonian and the flow dynamics influenced by external fields (e.g. gravitational or magnetic), in addition to other specific characteristics of the different applications that can make the flow even more complex. Taking all these effect into account at once in a single configuration makes it difficult to disentangle the fundamental role of each parameter on the flow characteristics. Even the simple canonical case of the turbulent flow of mono-dispersed, rigid non-Brownian, neutrally-buoyant and spherical particles suspended in an incompressible Newtonian fluid is still subject of active research, with numerous fundamental questions remaining unanswered (Prosperetti, 2015). The present work addresses this flow. In particular, we aim at describing finite-size effects on a dense turbulent suspension. Finite-size effects cause deviations of certain observables from the case where the dispersed phase can be modeled by an effective viscosity or by a more complex constitutive law for the stresses (i.e. for additives much smaller than the smallest relevant flow scale as colloidal particles, polymers and fibers).

Modeling the effects of the suspending particles in the flow through a localized source/sink of momentum (Crowe et al., 1977) is unrealistic when finite-size effects are important. In these cases the Navier-Stokes equations governing the suspending fluid should be coupled to the Newton-Euler equations governing the particle motion directly through the imposition of no-slip and no-penetration boundary conditions at the particle surface. This results in a relatively complex problem, difficult to tackle without three-dimensional and time-resolved flow data.

Throughout the years many studies on laminar shear flows laden with non-Brownian solid particles have been reported. In these cases one can take advantage of the linearity of the Stokes equations and achieve a rich variety of results by e.g. superposition of several canonical solutions. An iconic example is the effective viscosity of a suspension of non-Brownian spheres in the dilute and Stokes limit, derived by Einstein (1905) to be $\mu_{e}/\mu=1+(5/2)\Phi$, with $\mu$ being the viscosity of the suspending fluid, and $\Phi$ the bulk solid volume fraction. Many studies followed throughout the years (see e.g. Stickel & Powell, 2005; Guazzelli & Morris, 2011; Brown & Jaeger, 2014, for more rheological studies).

Most of the experimental works on particle suspensions have been limited to integral quantities, in particular the total wall shear. Often, the torque measured in a Taylor-Couette system required to keep a certain shear rate is used to measure effective viscosities and obtaining insight in the suspension rheology (see e.g. Bagnold, 1954; Hunt et al., 2002; Stickel & Powell, 2005). Although many important results have been extracted from this approach, more detailed measurements of important features, such as the micro-scale organisation or the particle dynamics are challenging to obtain.

Lack of direct measurements at the micro-scale level gave room for important analytical studies that relate the bulk suspension behavior to the particle dynamics. A well-known example is the work of Leighton & Acrivos (1987)\pc, who introduced the concept of shear-induced migration due to irreversible inter-particle interactions\pc, to explain the migration of particles to the fluid reservoir in their Taylor-Couette experiments.

Great progress on the understanding of the flow dynamics at the particle scale has been achieved through numerical simulations. Here one can measure directly the suspension micro-structure and particle dynamics. In particular for flows in the Stokes regime, one can use very accurate and relatively inexpensive particle-based methods such as stokesian dynamics (Brady & Bossis, 1988) to understand in detail the bulk suspension behavior from a micro-scale perspective. Recent advances in experimental techniques made also possible direct measurements of a suspension micro-structure (see e.g. Blanc et al., 2013).

When inertial effects become significant, the governing equations for the fluid phase are non-linear and therefore analytical descriptions, numerical simulations and even experiments become more challenging. Moreover, if the Reynolds number is sufficiently high, the flow becomes turbulent, exhibiting chaotic and multi-scale dynamics. This makes the understanding of these flows even more difficult.

Despite this, significant progress has been made with regard to wall-bounded particle-laden flows in the turbulent regime, in the point-particle limit. Several studies have been carried out analytically (e.g. Maxey & Riley, 1983; Reeks, 1977), numerically (Soldati & Marchioli, 2009; Eaton & Fessler, 1994; Sardina et al., 2012) and experimentally (Fessler et al., 1994). As a result, a deep understanding of the mechanisms e.g. leading to preferential accumulation of particles smaller than the Kolmogorov scale depending on their inertia and the local flow characteristics has been achieved. Similar maturity for the cases where the feedback of particles in the flow or finite size effects are relevant is still far from being accomplished.

The previously mentioned challenges and limitations of experimental and theoretical approaches makes the use of advanced numerical tools a necessity for obtaining detailed information; despite the well-known limitations in terms of Reynolds numbers that can be reached in simulations (Prosperetti, 2015). Lately, several groups have been successfully using numerical algorithms for interface-resolved direct numerical simulations (DNS) of different turbulent flows laden with finite size particles: examples are suspensions in isotropic turbulence (Ten Cate et al., 2004; Lucci et al., 2010), vertical channel flow (Uhlmann, 2008), sedimentation (Chouippe & Uhlmann, 2015; Fornari et al., 2016b), bed load transport (Kidanemariam & Uhlmann, 2014; Vowinckel et al., 2014), channel transport of mono-dispersed particles (Wang et al., 2016; Yu et al., 2016; Wang et al., 2017), and recently of poly-disperse \pc(Lashgari et al., 2017; Fornari et al., 2018) and non-spherical particles \pc(Ardekani et al., 2017; Eshghinejadfard et al., 2017). Likewise, the present work uses such simulations to study turbulent channel transport of neutrally-buoyant finite size spheres.

Suspensions of neutrally-buoyant particles close to the onset of turbulence have also been explored both experimentally (Matas et al., 2003) and numerically (Loisel et al., 2013). A common feature of the results from these studies is the fact that solely an increase in effective viscosity of the suspension does not explain the observed phenomenology, in particular for large particles.

In particle suspension flows a new mechanism for momentum transport emerges in the form of a particle stress (Batchelor, 1970). It therefore makes sense to follow the idea of Lashgari et al. (2014), who distinguished three different flow regimes, depending on the relative importance of viscous, Reynolds and particle stresses to the total stress of the suspension: a laminar (low Reynolds number and low volume fraction), turbulent (high Reynolds number and moderate to low volume fractions) and inertial shear-thickening (high volume fractions), the latest regime being characterized by a significant increase in wall shear that is not accompanied by an increase in the magnitude of the turbulent stresses.

Recently, Picano et al. (2015) presented detailed direct numerical simulations of turbulent channel flow laden with neutrally-buoyant finite size spheres. They showed that, for fixed Reynolds number and particle size, particle stresses at a volume fraction of about $20\%$ are responsible for a non-monotonic behavior of the near-wall peak in Reynolds shear stresses as a function of the volume fraction. The associated turbulence attenuation is higher than what predicted only by accounting for an effective viscosity. The decrease in the magnitude of the Reynolds stresses is accompanied by a more dominant increase in particle-induced stresses, which ultimately results in an overall drag increase.

This observation is consistent with flow regime map of Lashgari et al. (2014) and supported by the recent work of Costa et al. (2016). The authors built up upon the work of Picano et al. (2015) by extending the data set to higher Reynolds number and smaller particle sizes. They showed that a layer of near-wall particles causes the suspension to deviate from the continuum limit, where its dynamics is well represented by an effective suspension viscosity. Away from this layer, the suspension mean flow is shown to be well described by an effective suspension viscosity. Based on this idea Costa et al. (2016) successfully derived the scaling laws for the mean flow in the overlap region, and accurate correlations capable of predicting the overall drag that the suspension experiences as a function of the three governing parameters: $\Rey_{b}$, $\Phi$, and $D_{p}/h$, respectively the bulk Reynolds number, solid volume fraction and particle size ratio.

In this work we investigate how finite size effects change the flow and particle dynamics near the wall, and up to which point it is actually important to consider them in the bulk of the channel. We related the observed mesoscale behavior to the local micro-scale dynamics. This allowed us to present a clear picture of the sources of finite size effects in a dense turbulent suspension, needed for future modeling efforts.

2 Methods and computational setup

We performed interface-resolved DNS of a particle suspension flowing through a plane channel in the turbulent regime. The numerical algorithm solves the continuity and Navier-Stokes equations for an incompressible Newtonian fluid with density $\rho$ and kinematic viscosity $\nu$,

[TABLE]

where $\mathbf{u}$ is the fluid velocity vector, $p+p_{e}$ the fluid pressure with respect to an arbitrary constant reference value and $\boldsymbol{\nabla}p_{e}$ corresponds to a constant pressure gradient that may serve as driving force for the flow.

This set of equations is solved together with the Newton-Euler equations governing the motion of a spherical particle with mass $m_{p}$, and moment of inertia $I_{p}$,

[TABLE]

where $\mathbf{U}$ and $\boldsymbol{\Omega}$ denote respectively the particle linear and angular velocity vectors, $\boldsymbol{\sigma}\equiv-(p+p_{e})\mathbf{I}+\rho\nu\left(\nabla{\mathbf{u}}+\nabla{\mathbf{u}}^{T}\right)$, $\mathbf{r}$ a position vector with respect to the particle centroid, $\mathbf{n}$ the outward-pointing unit vector normal to the particle surface, $A$ the surface area of the particle, and $\mathbf{F}_{c}$ and $\mathbf{T}_{c}$ correspond to external forces and torques associated with short-range inter-particle or particle-wall interactions (such as solid-solid contact).

Equations (1)-(2) and (3)-(4) are coupled through the imposition of no-slip and no-penetration boundary conditions at the particle surface:

[TABLE]

A straightforward way of satisfying equation (5) is by discretizing the equations in a grid that conforms to the particles’ surfaces – \pca so-called body-fitted method. Despite the recent progress towards efficient implementation of such methods (see Vreman, 2016), such implementation remains of prohibitive computational demand for the dense cases considered here. Instead, we use the direct forcing immersed boundary method (IBM) developed in Breugem (2012). The idea of an IBM is that, instead of directly applying the boundary conditions expressed in equation (5) by conforming the computational grid to the \pcsurface of all the particles, one applies a smooth force distribution in strategic regions of the domain such that the boundary condition is satisfied with sufficient accuracy. This way one can benefit from efficient solvers for the Navier-Stokes equations on a regular Cartesian grid.

The fluid flow is solved with a finite-volume pressure correction scheme, coupled with an IBM. Our IBM requires a quasi-2D Lagrangian grid that discretizes the particle surface. The fluid prediction velocity is interpolated from the Eulerian to a Lagrangian grid. There the force required for satisfying no-slip and no-penetration is computed. Subsequently, the force is spread back to the Eulerian grid. A regularized Dirac delta function with support of 3 grid cells is used to perform interpolation/spreading operations (Roma et al., 1999; Uhlmann, 2005). Regularization of the particle-fluid interface results in first-order spacial accuracy. Slight inward retraction of the Lagrangian grid circumvents this issue, allowing for second-order accuracy in space. Errors due to spreading operations from different neighboring Lagrangian grid points to the same Eulerian grid points are reduced via a multi-direct forcing scheme. We refer to Breugem (2012) for more details on this IBM and Navier-Stokes solver.

Short-range hydrodynamic particle-particle and -wall interactions are reproduced by the IBM when the gap sizes are sufficiently well resolved. However, since the method uses a fixed grid, a lubrication model is inevitably needed once the inter-particle gap width becomes smaller than the grid spacing. The closures for short-range particle-particle and particle-wall interactions used in the present study \pchave been validated in Costa et al. (2015), \pcand are summarized below.

\pc

Normal lubrication interactions unresolved by the IBM, and roughness effects are accounted for by a two-parameter model. When a sphere approaches a wall or another sphere, reaching a normalized gap distance $\varepsilon_{\Delta x}$ (with $\varepsilon\equiv\delta_{ij,n}/R_{p}$, with $\delta_{ij,n}$ being the gap distance between particles $i$ and $j$), a lubrication correction force is computed. The parameter $\varepsilon_{\Delta x}$ corresponds to the non-dimensional threshold distance below which the IBM cannot resolve canonical lubrication interactions in the Stokes regime. The closure uses a correction force $\Delta F_{lub}=-6\pi\mu R_{p}u_{ij,n}(\lambda(\varepsilon)-\lambda({\varepsilon_{\Delta x}}))$, added to the integration of the particle momentum equation (3), where $\lambda$ is the Stokes amplification factor. To account for rougnhess effects, this correction is made independent of the gap-width for $0<\varepsilon<\varepsilon_{\sigma}$: $\Delta F_{lub}\propto u_{ij,n}(\lambda(\varepsilon_{\sigma})-\lambda{\varepsilon_{\Delta x}})$. When the two surfaces start to overlap, the lubrication correction is set to zero and a soft-sphere collision model takes over. These corrections are based on asymptotic expansions of analytical solutions for short-range particle-particle and particle-wall interactions in the Stokes regime, see e.g. Dance & Maxey (2003) for a review.

\pc

The soft-sphere model computes the normal collision force from the following linear spring-dashpot model:

[TABLE]

\pc

where $\boldsymbol{\delta}_{ij}$ and $\mathbf{u}_{ij,n}$ are the gap-width and relative particle velocity projected in the line-of-centers $\mathbf{n}_{ij}=(\mathbf{x}_{j}-\mathbf{x}_{i})/(||\mathbf{x}_{j}-\mathbf{x}_{i}||)$; $\mathbf{x}_{i/j}$ corresponds to the centreline position of particles $i/j$. The spring and dashpot coefficients are given by: \pc

[TABLE]

where $e_{n,d}$ is the dry coefficient of restitution, and $m_{e}=\left(m_{i}^{-1}+m_{j}^{-1}\right)^{-1}$ the reduced mass of the particles. $N\Delta t$ is the collision time, set as a multiple $N$ of the time step of the Navier-Stokes solver, $\Delta t$. This allows the flow solution to gradually adapt to the sudden changes in particle velocity (Costa et al., 2015; Biegert et al., 2017). The model for the tangential component of the collision force in Costa et al. (2015) is not described here, as in the present study we consider the particles to be frictionless; as discussed later, frictional interactions can be neglected for the flow governing parameters at stake. The values of $\varepsilon_{\Delta x}$, $\varepsilon_{\sigma}$ and $e_{n,d}$ used are given in table 1. \pc

\pc

The parameters in table 1 are taken from Costa et al. (2015), where a thorough validation is performed against several canonical cases, like the trajectory of a sphere colliding onto a planar surface in a viscous liquid, and wet coefficients of restitution for head-on particle-wall and particle-particle collisions.

The flow dynamics is governed by three parameters: the bulk Reynolds number $\Rey_{b}\equiv U_{b}(2h)/\nu$, the particle size ratio $D_{p}/h$, and the bulk volume fraction of solid particles $\Phi=N_{p}V_{p}/V_{t}$ (note that the particles are neutrally-buoyant), where $U_{b}$ is the flow bulk velocity (forced to be constant in the numerical algorithm), $h$ the half channel height, $D_{p}$ the particle diameter, $N_{p}$ the total number of particles, and $V_{p}$ and $V_{t}$ the volumes of a particle and of the computational domain. In the present work the bulk volume fraction of solid particles is fixed to $\Phi=20\%$, and the Reynolds number to $\Rey_{b}=12\,000$, in order to ensure sufficient inner-to-outer scale separation Costa et al. (2016). Two different particle sizes are considered as reported in table 2, together with other relevant physical and computational parameters.

\pc

We explore the dataset of dense suspension flow reported in (Costa et al., 2016). The simulations were performed in a domain with dimensions $L_{x}/h\times L_{y}/h\times L_{z}/h=6\times 2\times 3$, using a grid spacing dictated by the number of grid cells required to resolve the flow conforming the spheres: $\Delta x/D_{p}=1/16$\pc, which results in a distribution of $746$ Lagrangian grid points on the surface of each particle. The interface-resolved simulations were complemented with two single-phase reference cases: the unladen case at the same bulk Reynolds number, denoted SPR (single-phase reference), and the continuum limit where the flow dynamics can be reproduced by a single-phase fluid with the effective viscosity of the suspension at $\Phi=20\%$, denoted as CLR (continuum limit reference). This simulation corresponds thus to a single-phase flow with bulk Reynolds number $\Rey_{b}\nu/\nu^{e}\approx 6\,400$ and effective viscosity $\nu^{e}=\nu(1+(5/4)\Phi/(1-\Phi/\Phi_{max}))^{2}$ obtained from Eilers fit Stickel & Powell (2005), with $\Phi_{max}=0.6$.

\pc

Some remarks on the choice of the collision model parameters should be made at this point. First, the choice of $e_{n,d}=0.97$ is motivated by values measured experimentally for dry collisions of rigid spheres, (Foerster et al., 1994; Joseph et al., 2001). However, the results should not be sensitive to small changes in $e_{n,d}$. The reason is that, as we will see, the particle inertia is sufficiently small that the resulting impact Stokes number $\mathrm{St}\equiv(1/9)\rho_{p}u_{ij,n}D_{p}/\mu<10$. In this scenario, hydrodynamic viscous effects dominate the collision dynamics and the particles typically do not rebound, despite the value $e_{n,d}\approx 1$ (Legendre et al., 2006; Yang & Hunt, 2006). Second, the effect of particle-particle and particle-wall solid friction in the suspension dynamics is small. For particle-particle interactions, the mean separation distance ($\sim\Phi^{-1/3}D_{p}$) is sufficiently large that sustained frictional contacts are unlikely. In regard to particle-wall interactions, there is no mechanism causing a large net force for neutrally-buoyant particles, pushing them towards the wall. Since the frictional force is proportional to this normal force, solid-solid friction at the wall is negligible. In conclusion, the main role of the collision model in this flow is ensuring that particles occupy their own volume; i.e. once particles are in the imminence of contact, they cannot overlap. Then the solid-fluid coupling dictate the flow dynamics. Finally, we do not implement lubrication corrections in the tangential direction. These choices are supported by the success of this approach in reproducing the effective viscosity of a suspension up to $\Phi=30\%$ in Picano et al. (2013), where friction and tangential lubrication corrections were neglected.

\pc

Unless otherwise stated, spatial averages of an observable $o$ pertaining to a certain phase are obtained from the following intrinsic volume average:

[TABLE]

where $\phi(x,y,z)$ is a phase-indicator function ($\phi_{ijk}^{f}+\phi_{ijk}^{s}=1$), and the superscripts $f$ and $s$ denote fluid and solid, respectively. For simplicity of notation we will drop the brackets $\left<\right>$. Profiles pertaining to the combined phase are obtained similarly,

[TABLE]

When computing the profiles of particle velocity, the Eulerian grid points in the solid domain are populated by the velocity corresponding to rigid-body motion: $\mathbf{U}_{p}(x,y,z)=\mathbf{U}+\boldsymbol{\Omega}\times(\mathbf{x}(x,y,z)-\mathbf{x}_{c})$, with $\mathbf{x_{c}}$ the particle centroid position.

Due to the large size of the simulations, all spatial averages were computed during run time and averaged in time a posteriori, with samples obtained from a time interval of about $1500h/U_{b}$, spaced by $7h/U_{b}$. A total of about 7 million core hours in the supercomputer CURIE (thin nodes; B510 bullx) at CEA, France, and Beskow (Cray XC40) at KTH, Sweden were used.

3 Results

Near-wall dynamics – the main cause of finite size effects

We will here refer to finite-size effects as the causes for the meso or macroscale-averaged flow statistics to be different from those obtained from the continuum limit approximation in which the suspension is modeled as a Newtonian fluid with an effective viscosity due to the presence of the particles.

Costa et al. (2016) showed that a layer of particles, flowing near the wall with significant (apparent) slip velocity, is responsible for an additional increase in drag which cannot be modeled by an effective suspension viscosity. By accounting for this effect, the authors were able to scale the profiles of mean velocity and velocity defect ($U_{c}-u$, where $U_{c}$ is the mean centerline velocity) for a wide range of combinations of $\Rey_{b}$, $\Phi$ and $D_{p}/h$. The theory in Costa et al. (2016) assumes that the domain can be split into two regions, bounded at a wall-normal distance $y=\delta_{pw}$: a region away from the wall \pcwith uniform particle concentration, denoted homogeneous suspension region (HSR) ($y>\delta_{pw}$) where the mean flow of the suspension is well represented by the continuum limit of a Newtonian fluid with effective viscosity $\nu^{e}$, and a region close to the wall denoted particle-wall layer (PWL) ($y<\delta_{pw}$) where the difference between the mean flow of the two phases, at the mesoscale level, makes this assumption invalid. By exploiting the stress budget, the authors could derive scaling laws for the mean velocity and velocity defect in the HSR, provided that, like in a single-phase turbulent channel, inner-to-outer scale separation is satisfied. In other words, the friction Reynolds number based on the scaling parameters of the homogeneous suspension region needs to be sufficiently high. These are a channel height corrected for a virtual-wall origin, $h-\delta_{pw}$, the friction velocity taken from the profile of the total stresses at $y=\delta_{pw}$, $u_{\tau}^{*}=\sqrt{u_{\tau}^{2}(1-\delta_{pw}/h)}$ and the effective suspension viscosity $\nu^{e}$, yielding $\Rey_{\tau}^{e*}=u_{\tau}^{*}(h-\delta_{pw})/\nu^{e}$. The relations derived in Costa et al. (2016) read,

[TABLE]

with $\delta_{v}^{e*}=\nu^{e}/u_{\tau}^{*}$; the coefficients $\kappa$, $B$ and $B_{d}$ retain the values of single-phase flows, whereas the effective viscosity is obtained from Eilers fit $\nu^{e}/\nu=(1+(5/4)\Phi/(1-\Phi/\Phi_{max}))^{2}$ (Stickel & Powell, 2005). The virtual wall origin of the HSR is given by $y=\delta_{pw}=C(\Phi/\Phi_{max})^{1/3}D_{p}$. These scaling laws can be further used to derive a master equation that accurately predicts the overall drag of the suspension, here expressed in terms of a friction Reynolds number:

[TABLE]

or, in an explicit form based on a similar correlation for single-phase flow from Pope (2001),

[TABLE]

where $\xi_{pw}=(1-\delta_{pw}/h)$ and $\chi^{e}=\nu/\nu^{e}$.

Table 3 presents the mean wall shear, expressed in terms of the friction Reynolds number $\Rey_{\tau}=u_{\tau}h/\nu$ from the different simulations considered here. As expected from eqs. (12) and (13), the addition of finite size neutrally-buoyant particles results in an increase in drag with respect to the value of 350 of the single phase flow, which (when fixing the other governing parameters) increases with increasing particle size. We first note that the smaller the particles, the more the suspension friction Reynolds number $\Rey_{\tau}^{e}=\Rey_{\tau}\nu/\nu^{e}$ approaches the value obtained in the continuum limit, CLR; in other words using an effective viscosity provides a better prediction of the total drag. This is a consequence of the reduced finite-size effects that typically occur near the wall (Costa et al., 2016).

Finally, the third column in the table shows the suspension Reynolds number of the HSR, $\Rey_{\tau}^{e*}=\Rey_{\tau}^{e}(1-\delta_{pw}/h)^{3/2}$, obtained considering both an effective suspension viscosity and a virtual wall origin as explained above (with $\delta_{pw}$ computed with $C=1.5$ and $\Phi_{\max}=0.6$). The friction Reynolds numbers from the interface-resolved simulations are, in this case, equal for both particle sizes and very close to the value of $\Rey_{\tau}\nu/\nu^{e}$ for the CLR case. This strongly supports the proposed correction for finite-size effects.

Figure 1 displays snapshots of the flow for the two laden cases, D10 and D20. The top panels show planar sections colored by the streamwise flow velocity, and particle positions only shown for streamwise locations smaller than $0.4L_{x}$ and wall-normal distance $y<h$. The bottom plane (also shown in the two planar views, in the bottom panels) corresponds to a wall-normal position $y=1.5D_{p}$, corresponding to about $y/\delta_{v}=15$ and $30$ for cases D10 and D20. The figures show the typical near-wall low-speed streaks found in single-phase wall-bounded turbulence, also present in these suspensions (Picano et al., 2015). Differences in velocity contrast between the two cases can be attributed to the choice of a plane at wall-normal distance that scales with the particle diameter, meaning that, in viscous units, the plane corresponding to case D10 is closer to the wall. However, in both cases the maximum streak amplitude is achieved at relatively larger wall-normal distances than what is expected for single-phase flows, as illustrated in figure 2. This figure depicts contours of autocorrelations of the streamwise suspension velocity (i.e. computed from the fluid flow field with rigid body motion inside the particles) $R_{uu}^{z}(y,\Delta z)\equiv\left<u^{\prime}(y,z)u^{\prime}(y,z+\Delta z)\right>/\left<u^{\prime 2}(y)\right>$, averaged in time and in the streamwise direction for different spanwise separations. Panels (b) and (c) show the results for the interface-resolved simulations D10 and D20. Clearly the near-wall minima – footprint of the low- and high-speed streaks (Kim et al., 1987), and near-wall property that should scale in inner units – are shifted upwards, and the larger the particles are, the larger this shift is.

Figure 2 also shows that the smaller the particles are, the closer the suspension dynamics resemble the continuum limit, CLR. \pcThis can be seen when comparing panel (d), corresponding to the continuum limit CLR to panels (b) – D10 and (c) – D20. There is a clear qualitative discrepancy between the CLR and the bigger particles (D20), contrasting with a significant quantitative agreement between CLR and the case with small particles (D10). For case D20, the wall-normal minimum is located at such a large distance to the wall ($\approx 0.3h$) that inner-to-outer scale separation is compromised. Also, the autocorrelations for the interface-resolved cases show a non-monotonic trend at wall-normal distances $y\lesssim D_{p}$, where the autocorrelation (for fixed $y$) reaches a local maximum close to the wall. This can be even better observed from the profile of the integral length scale $L_{xz}=\int_{-L_{z}/2}^{+L_{z}/2}R_{uu}^{z}\,\mathrm{d}(\Delta z)$, indicated by the dashed lines in the figure. This confirms that the flow dynamics in this region is qualitatively different from what a simple continuum rheological description would predict.

Flow dynamics near the particle-wall layer (PWL)

Next we investigate in more detail the flow dynamics near the particle-wall layer. Panel (a) of figure 3 depicts the mean streamwise fluid velocity in inner units. Note that the viscous wall unit used for inner-scaling is defined, for consistency, with the effective suspension viscosity at the same volume fraction, as we are ultimately interested in the deviations of the flow dynamics from the continuum limit. For wall-normal distances $y\gtrsim 10\delta_{v}^{e}$, the profile for the case D10 (smaller particles) clearly shows a logarithmic scaling, with a von Kármán constant $\kappa=0.36$. On the other hand, case D20 does not show a clear logarithmic region because of the larger extent of the PWL (Costa et al., 2016).

Panel(b) of figure 3 presents the difference between the fluid- and solid-phase velocity profiles, expressed in terms of a particle Reynolds number $Re_{p}=(u_{p}-u_{f})D_{p}/\nu$. Profiles pertaining to the solid phase are obtained by averaging over the rigid body motion of the particles. Two regions can be clearly distinguished in the figure, roughly separated by the line marking the wall-normal distance $y=D_{p}$. For $y\gtrsim D_{p}$, the difference between the velocity of the two phases is small, whereas the profiles clearly deviate for wall-normal distances $y\lesssim D_{p}$, reaching the highest (apparent) particle-to-fluid slip velocity at the wall. This is a signature of particle layering due to the kinematic constraint that the wall imposes on the particles. Particles flowing at the wall acquire most of their linear momentum at wall-normal distances higher than their radius and this is transported uniformly throughout their volume. Conversely, the fluid momentum must vanish at the wall. Notably, Figure 3 (b) shows (see also its inset) that the particle Reynolds number in the bulk is virtually zero for case D10, while it still assumes values of $O(1)$ for case D20. This evidence of finite particle inertia in the bulk can result in inertial shear-thickening effects (see Picano et al., 2013). We will come back to this when the dynamics of the homogeneous suspension region are discussed. \pcNote that in this figure both the effect of particle translation and rotation is accounted for. To discern the effect of rotation, panel (c) of figure  3 compares the profile of particle and fluid angular velocity for the two cases. A similar trend is observed: larger deviations with respect to the fluid profile is found for larger particles, and the highest discrepancy between fluid and particle statistics occurs close to the wall.

The particle near-wall layer impacts the scale separation on which the scaling laws for mean velocity are based. It is therefore interesting to understand how the near-wall inhomogeneity affects the particle structure, and how far this inhomogeneity extends in the wall-normal direction. We can quantify this through the angular distribution function at contact (adf). The adf measures the probability of finding a particle pair at a fixed distance $r$, as a function of the polar ($\theta$) and azimuthal ($\varphi$) angles, normalized by the values of a random particle distribution. Formally,

[TABLE]

where $N_{\theta,\varphi}$ denotes the number of particles on a segment of a spherical surface of radius $r=D_{p}/2$, and polar and azimuthal angles within the ranges $\theta^{\prime}\in[0,\theta]$ and $\varphi^{\prime}\in[0,\varphi]$. The adf is computed from the simulation data using bins with wall-normal extent $D_{p}$, centered at $y_{\pc{bin}}/D_{p}=1$, $2$, and $3$, i.e. at the same wall-normal distance if scaled with the particle size \pc(see the bottom-left corner of figure 4).

The projection of the adf onto the $x-y$ plane is displayed in figure 4 for the two cases under investigation. For a suspension in homogeneous shear, the particles tend to be preferentially attracted towards each other when located in a compression region (negative local strain) and reciprocally repelled from each other when located in a stretching region (Morris, 2009). This results in a distribution with two planes of symmetry: two regions of higher clustering in the $xy>0$ quadrants and two of lower accumulation in the $xy<0$ quadrants. The inhomogeneity introduced by the wall, bounding the particle trajectories, induces a deviation from this picture. The wall-effect vanishes for wall-normal distances $y/D_{p}\gtrsim 3$ for both simulations. Interestingly, both adf agree quantitatively. As the distance from the wall is normalized by the particle diameter, for fixed volume fraction, the wall inhomogeneity affects the particle dynamics at distances proportional to $D_{p}$, despite the different local behavior (different wall-normal location in viscous units) of the suspending fluid.

\pc

The local volume fraction profile also presented in figure 4 illustrates that the near-wall inhomogeneity vanishes above $y\sim 4D_{p}$. For higher wall-normal distances this profile is uniform (Costa et al., 2016). It is interesting to note that a denser near-wall layer was existing at lower Reynolds numbers for the same volume fraction (see e.g. Picano et al., 2015; Wang et al., 2016; Yu et al., 2016). This decrease of layering is attributed to the stronger mixing in this flow due to a higher Reynolds number.

Let us now investigate the influence of the particle-wall layer on the distribution of wall shear-stress. To this end, we report in figure 5 contours of the instantaneous wall-shear for the case D20, together with the particle positions close to the wall (shown with transparency). The data reveal a strong correlation between regions of high shear and the particle location. This can explain the increase in drag with respect to predictions from the continuum limit: high particle-to-fluid apparent slip velocity close to the wall corresponds locally to events of high wall shear stress, which are not present in the ‘equivalent’ single-phase flow with a modified suspension viscosity. Panel (b) of the same figure shows the mean wall shear conditioned to the particle positions. Results for small and large particles show a similar qualitative behavior: the shear is higher than the mean around the particles, except close to the surface of the particle, in the spanwise direction. This is a consequence of the reduction of streamwise velocity as the fluid moves around the particle. The contours of the shear stress for the smaller particles show a higher local increase in shear stress, which is, on the other hand, more localized. In fact, the contribution of these hot spots of shear stresses to the total shear \pc(computed from a conditional average of the wall shear stresses in a $2D_{p}\times 2D_{p}$ underneath the particle) is the same for both cases, D10 and D20, and quite high: $69\%$.

The probability density function (pdf) of the shear stress distribution for cases D20, D10 and SPR, presented in figure 6, shows a clear difference between the typical distribution of shear stresses for a single-phase flow, and for the particle-laden cases. In the presence of particles, the tails of the pdf are greatly widened, with very high probability of finding values of the shear stress up to $4$ times the mean value, and likewise for values lower than the mean, including negative events. We note there that it is reported for single-phase flow that $\tau_{w}^{+,rms}=\tau_{w}^{rms}/\left<\tau_{w}\right>\approx 0.4$ and that the pdf typically follows a log-normal distribution (Örlü & Schlatter, 2011). Indeed, we obtained a value of $\tau_{w}^{+,rms}=0.39$ for case SPR, with a pdf well-fitted by a log-normal distribution (not shown). The hot-spots of high shear cause a striking difference for the laden cases: $\tau_{w}^{+,rms}=1.1$ and $1.0$ for cases D10 and D20. Moreover, the pdf no longer fits a log-normal distribution. The hot-spots of high wall shear induce a kink in the pdf for $\tau\gtrsim\left<\tau\right>$, resulting in a clear exponential tail, fitted by $\lambda\exp(-\lambda\tau/\left<\tau\right>)$ with $\lambda\approx 1$. Further, case D10 shows higher probability of extreme events, (which, again, are more localized). The largest values attained by the shear stress that in single-phase flows do not contribute significantly to the mean shear, account for about $10\%$ of the total in particle-laden channel flows, as shown by the cumulative probability distribution function (cdf) displayed in panel (b).

The results presented so far shed light on the relation between the wall slip velocity $U_{pw}=u_{p}(y=0)$ and the wall friction $u_{\tau}$ (recall figure 3 b). Assuming that (1) the average shear in the hot-spots matches the average wall shear and (2) the particles impose a mean shear that scales with $U_{pw}/D_{p}$ over an area $D_{p}L_{w}$ ($L_{w}$ being the streamwise extent of the hot-spot, which as seen in figure 5 is significantly larger than $D_{p}$), we can write:

[TABLE]

with \pc$N_{pw}$ the number of particles in this layer. Further noting that $N_{pw}D_{p}^{3}/(D_{p}L_{x}L_{z})\sim\Phi$, we have

[TABLE]

and finally, assuming that $L_{w}/D_{p}$ scales as the mean-free-path in a two-dimensional layer of particles, i.e., $\propto\Phi^{-1/2}$, we obtain

[TABLE]

with $Re_{\tau}^{p}=u_{\tau}D_{p}/\nu$ a particle friction Reynolds number. The data in figure 7, including previous cases in literature, confirm the validity of this relation with $C_{pw}\approx 0.2$. The data points that do not follow the predicted scaling correspond to very large particles ($h/D_{p}=5$), for which the influence of the outer scales on the near-wall particle dynamics compromises the first assumption above, and to low volume fractions.

Homogeneous suspension region (HSR)

The contribution of the different stresses to the suspension streamwise momentum transport (i.e. the stress budget) can be derived from volume- and ensemble-averaging the mixture streamwise momentum equation (see Marchioro et al., 1999; Picano et al., 2015). It reads for plane channel flow

[TABLE]

where the terms on the right-hand-side denote the fluid Reynolds stresses $\tau_{{T}_{f}}$, the viscous stresses $\tau_{v}$, the particle Reynolds stresses $\tau_{{T}_{p}}$, and the particle stresses $\tau_{p}$, the latter related to the stresslet, moment of the material acceleration acting on a particle and inter-particle collisions (Guazzelli & Morris, 2011); in the expression above\pc, $\phi$ denotes the mean local solid volume fraction. The wall-normal profiles of the different contributions to the total stress (obtained dividing equation (18) by $u_{\tau}^{2}(1-y/h)$) are reported in figure 8.

As also observed in Lashgari et al. (2016, 2014); Picano et al. (2015) the profiles of particle stresses reach a maximum at $y=D_{p}/2$, at the location of the local maximum of the profiles of the local volume fraction (see figure 4), followed by a minimum at $y\sim 2D_{p}$. The Reynolds stresses are relatively small in this region, and therefore the local maximum/minimum of the particle stresses is compensated by a local minimum/maximum of the viscous shear stresses. This is caused by the strong shear that the first layer of particles, flowing with significant slip velocity, imposes on the fluid above it. Further away from the wall, the profiles follow a trend that resembles the one of the shear-stress profile, suggesting a linear (Newtonian) scaling of the particle stresses with the shear rate.

The deviation from the continuum limit (CLR) can be examined by comparing the sum of the particle and viscous stresses to the profile of viscous stresses in panel (d) of figure 8 (note that an effective viscosity incorporates the effects of both viscous and particle stresses). As expected, the case with smaller particles is much closer to this continuum limit than the case with larger particles (see also Costa et al., 2016). For the latter, while the profile pertaining to the viscous contribution in the bulk is still close to the one of case D10, the particle stress contribution is much larger. This difference can be attributed to the inertial-shear-thickening mechanism due to excluded volume effects proposed in Picano et al. (2013), as the particle Reynolds number based on the particle slip velocity for case D20 is much larger (recall Figure 3 (b)). We should stress here the importance of defining the viscous wall unit from the effective suspension viscosity: if the profiles of cases D10 and D20 would be compared to CLR using the classical inner scaling (where the fluid viscosity and wall friction velocity are used to scale the wall-normal distance) these would not match, as the abscissa on the upper axis $y/\delta_{v}^{e}$ would have to be multiplied by $\nu^{e}/\nu=1.89$ in panels (a) and (b).

As for the other quantities presented so far, we expect the second-order Eulerian statistics for the laden cases to approach those of CLR with decreasing particle size. These are shown in figure 9 where we depict the profiles of the root-mean-square (rms) of the fluctuating fluid velocity and Reynolds stresses from the different simulations. When comparing to the single-phase flow statistics we also observe, as reported in Picano et al. (2015), that the fluctuating fluid velocity close to the wall is enhanced by the presence of the particles.

Let us take one step further and test the scaling arguments of Costa et al. (2016) for the second-order statistics. These can be compared directly to the case CLR, as it is meant to be the same flow, finite size effects aside. We therefore correct the wall-normal distance for the presence of the particle-wall layer through a virtual wall origin, $\delta_{pw}$, and correct velocity scale for inner-scaling, $u_{\tau}^{*}=u_{\tau}(1-\delta_{pw}/h)^{1/2}$ (panels a,c,e and g in the figure). The magnitude of near-wall peaks in streamwise velocity rms, typically located at $y=12\delta_{v}\lesssim D_{p}$, are not recovered in the interface-resolved simulations. Despite this, away from the wall (where the dynamics is only weakly affected by the particle-wall layer due to sufficient scale separation) the profiles show good agreement. Particularly, the profiles pertaining case D10 show a remarkable collapse over the entire outer region for all three components of the velocity rms and Reynolds stresses, even without correction. The same cannot be safely stated for the case with larger particles (D20) where finite size effects, also present away from the particle-wall layer, have an influence on the statistics in the bulk of the flow. Despite these small deviations for case D20, it is clear that, away from the wall, the second-order statistics are described with good approximation by those of a Newtonian fluid with an effective suspension viscosity.

Particle Dynamics

\pc

Analysing the mean flow solely in an Eulerian framework does not give direct insights on the dynamics of individual particles. To \pcbetter understand the particle dynamics, we further explore the DNS dataset by computing Lagrangian statistics. We first focus on the evolution of the single-point mean-square displacement of particles in the spanwise direction to prevent biases due to statistical inhomogeneity in the wall-normal, and non-zero mean flow in the streamwise direction (see also Lashgari et al., 2016). This is defined as

[TABLE]

where $\left<\right>$ denotes the time average over all particles located in a bin centered at a wall-normal position $y$ at time $t$ and wall-normal extent $D_{p}$. This observable is displayed in figure 10 versus time in units of a characteristic integral scale of the turbulent fluid motion, $h/u_{\tau}$.

One can \pcdistinguish two well-known regimes: the ballistic regime where the high temporal correlation results in a mean square displacement $\left<\Delta z_{1}^{2}\right>\propto\Delta t^{2}$ and the diffusive regime where the motion decorrelates from the initial sampling instant and $\left<\Delta z_{1}^{2}\right>\propto\Delta t$ (see e.g. Sierou & Brady, 2004). In the ballistic regime, particles reach the fastest dispersion at the distance from the wall corresponding to the peak in particle velocity rms (see figure 9). The diffusion coefficient (slope of the profile in the diffusive regime) is nearly the same for the both cases, away from the wall. Two main mechanisms are responsible for the spanwise particle displacements and the subsequent self-diffusion: (1) short-range inter-particle interactions and (2) interactions with turbulent structures of dimension larger than the particle size. One can picture a sequence of these events as successive random-walk steps during the particle motion. Indeed, Lashgari et al. (2016) observed in the turbulent regime an increase of the diffusion coefficient with increasing Reynolds number for fixed volume fraction, and a milder increase with increasing volume fraction at fixed Reynolds number. The smaller diffusion coefficients observed near the wall can therefore be explained by the constraint that the wall confinement imposes on the particle motion, and by the reduction of the characteristic turbulent integral scale, estimated as $l_{m}(y)\sim\min(\kappa y,0.1h)$ \pcnote: I changed to $l_{m}(y)$, with $l_{m}$ being a mixing length Pope (2001). The turbulent length scales are hindered up to a point where $l_{m}\sim D_{p}$, the value below which turbulent fluctuations have a weak influence on the particle kinematics. Despite the relatively small particle inertia, the temporal filtering due to the higher response time of the larger particles (Hinze, 1975) can also cause a decrease in diffusion coefficient. We expect however the latter effect to be milder, as increasing particle inertia through an increased mass density (while fixing the other governing parameters) has a small influence on spanwise particle dispersion, as shown in Fornari et al. (2016a).

The diffusion coefficient is larger for the flow with smaller particles (D10), even though the collision probability is larger for the largest particles (as suggested by the asymptotic limit of vanishing particle size at fixed volume fraction, and will be confirmed later). We can therefore conclude that, at the volume fraction under consideration here, the turbulent fluid motion is the main source of self-diffusion and, consequently, larger particles result in slower dispersion. Note also that the diffusive regime is reached away from the wall at times $\Delta t=O(h/u_{\tau})$ of the same order as the turnover time of large eddies in the bulk, consistently with this observation.

Further insight into the particle dynamics can be gained by examining the pair-dispersion statistics. The two-point mean spanwise square displacement is displayed in figure 11 for two particles at contact at $t=t_{0}$,

[TABLE]

where $\delta z$ denotes the spanwise interparticle distance and $t_{0}$ the instant at which the particles are in contact. For short time scales, lubrication dominates and the particles display a highly correlated motion. The duration of this regime increases with the distance to the wall. This slower relative dispersion is linked to the decreasing relative inter-particle velocity towards contact with increasing wall-normal distance, to be quantified later. The large particles, case D20, show faster pair dispersion, also due to higher relative inter-particle velocity towards contact: it turns out, as elaborated later when the collisional dynamics is addressed, that this quantity scales with the particle Reynolds number. At later times, the particle dispersion $\left<\Delta z_{2}^{2}\right>\propto\Delta t^{\alpha}$, with $2<\alpha<3$. Note that for tracer particles in a homogeneous isotropic turbulent flow $\alpha=3$ for separation distances within the inertial subrange (Salazar & Collins, 2009). For larger separation distances the diffusive limit is recovered.

The mean square displacement of tracer particles in turbulence with sub-Kolmogorov separation distance (i.e. in the dissipation subrange) grows exponentially in time as the relative particle velocity and the separation are proportional (Batchelor, 1952; Salazar & Collins, 2009): i.e. $\left<\delta z_{2}^{2}\right>=\left<\delta z_{0}^{2}\right>\mathrm{e}^{2\Delta t/t_{c}}$, where $t_{c}$ is a characteristic inter-particle response time. Although in the present study $D_{p}$ is $O(10)$ times larger than the smallest turbulence scale, we observe a clear exponential growth at short times, just after the first highly correlated regime when $\left<\Delta z_{2}^{2}\right>$ is approximately constant. The time at which this exponential regime sets in corresponds to a mean square separation of about ${2D_{p}^{2}/3}$. This exponential regime is highlighted in the bottom panels of figure 11 where the time on the horizontal axis is divided by $\Delta t^{l}=\Delta t|_{<\Delta z_{2}^{2}>=2D_{p}^{2}/3}$. This growth is due to a different mechanism than the one of tracer particles, possibly the combination of a uniform shear force (note that under uniform shear $\delta z$ is proportional to the relative velocity) following a short-range lubrication interaction.

Finally, we investigate the particle collisional dynamics. A collision event takes place when (1) the particles are at sufficiently close distance and (2) their relative velocity drives them towards contact. These two factors are investigated separately. The first by the radial distribution function (rdf), which quantifies the probability of finding a second particle at distance $r$ normalized by the probability of a random distribution of particles:

[TABLE]

where $N_{r}$ denotes the number of particle pairs in a spherical volume of radius $r$. Thus, if $g(r,y)$ assumes values larger than $1$, particles are preferentially sampled. The second observable is the distribution of the relative particle velocity projected in the direction of the line-of-centers $\Delta v^{n,-}$, given for two particles $i$, $j$ with velocities $\mathbf{u}_{i/j}$ and positions $\mathbf{x}_{i/j}$ by:

[TABLE]

where the $\max\{\}$ operator samples the relative velocity towards contact, as the superscript ‘$-$’ suggests. The product of these quantities measures the rate at which particles approach each other. Its value for $r=D_{p}$ is the so-called collision kernel, $\kappa_{c}$, the probability of a collision event.

As for the dispersion statistics, we take into account the inhomogeneity in the wall-normal direction by averaging in wall-parallel bins with wall-normal extent $2D_{p}$. Figure 12 presents the radial distribution function, negative particle relative velocity and the approach rate at contact $r=D_{p}$.

The profile of $\Delta v^{n,-}$ shows that the dominant mechanism driving particles towards each other is shear from particles at different wall-normal locations. Note that the shear rate away from the PWL is similar for both simulations considered (see figure 8 and table 3). Assuming that particles are driven towards each other by shearing a layered arrangement of particles:

[TABLE]

with $u_{\tau}/h$ an estimate of the shear rate in the bulk of the flow. This scaling is tested in the inset of panel (b) and (c), yielding a better agreement of the profiles of $\Delta v^{n,-}$ and $\kappa_{c}$, despite the differences in $g(D_{p})$. Note that $u_{\tau}/h$ is approximately the same for both cases (table 3). The higher values of $g$ at close distances for case D20 are also compatible with this picture, as larger particles are more prone to be driven towards each other by a shear-induced relative motion.

\pc

The mechanism governing energy dissipation during particle collisions can be predicted from the impact Stokes number, which can be computed directly from the mean relative approach velocity close to contact: $\mathrm{St}=(1/9)\rho_{p}\Delta v^{n,-}_{max}D_{p}/\mu$. The Stokes number corresponding to the maximum value of $\Delta v^{n,-}$, shown in figure 12(b) (case D20), is $\mathrm{St}\approx 0.28$. This average Stokes number is therefore much lower than the critical value of $10$ everywhere in the domain, suggesting that the collisional dynamics are governed by short-range hydrodynamic interactions, with strong viscous dissipation and no elastic rebound; see Legendre et al. (2006) and Yang & Hunt (2006).

Figure 13 shows contour plots of the same quantities as in figure 12, now as a function of both the separation distance $r$ and the wall-normal coordinate $y$. The contours of $g$ show that this quantity is weakly dependent on the wall-normal coordinate, with local maxima at $r/D_{p}\approx 1$ and $2$, consistent with the presence of statistically significant particle pairs and triplets. It is also interesting to notice that the global maximum of $g$ occurs at a finite distance to the wall. Comparing the two simulations, the maximum is located at the same wall-normal distance when scaled with the particle diameter $y\sim 3D_{p}$, in agreement with the observation above that wall-confinement effects are noticeable at distances proportional to the particle diameter (see figure 4). The maxima of $g(r)$ correspond therefore to the optimal trade-off between high shear (driving particles at different wall-normal locations towards each other) and low confinement; indeed further away from the wall, where the mean shear is relatively low, $g$ becomes almost independent of $y$. The maxima of $g$ close to the wall also suggest that larger particles have higher probability of forming particle pairs, and lower probability of forming triplets.

Panels (c) and (d) of the same figure display contours of $\Delta v^{n,-}$. For fixed separation distance $r$, the average inter-particle approaching velocity decreases with the wall-normal distance \pc$y$, which can be explained by a decreasing local shear rate: the lower the shear, the smaller shear-induced differences of the particle relative velocities. \pcConsistently, the variations of $\Delta v^{n,-}$ with $y$ are larger close to the wall. The negative relative particle velocity is larger for the largest particles, case D20, over the entire flow, \pcas the scaling suggested in equation (23) predicts. Finally, we report the rate-of-approach $\Delta v^{n,-}g$ in panels (e) and (f). Clearly, the differences in $g$ discussed above induce significant differences in this quantity only close to contact. At larger separation, the behavior of $\Delta v^{n,-}g$ is dictated by $\Delta v^{n,-}$.

4 Conclusions and outlook

We performed interface-resolved direct numerical simulations of turbulent channel flow of suspensions of neutrally-buoyant spherical particles. Two flow cases are considered, with same Reynolds number and volume fraction, and particle size varied by a factor of $2$. The simulations are compared to two single-phase reference simulations: (1) the unladen case and (2) the continuum limit of a Newtonian fluid with a viscosity corresponding to the effective viscosity of the suspension under investigation in a laminar shear flow.

As observed in Costa et al. (2016), the main finite size effect in the zero- and first-order Eulerian statistics originate from the near-wall dynamics, in the near-wall particle layer. We show here that the larger the particles are, the further away from the wall the effects of the particle layer are significant. The near-wall inhomogeneity due to the geometrical constraints is felt at distances that scale with the particle size at fixed volume fraction, seemingly independent of the local fluid flow dynamics. This results in a smaller inner-to-outer scale separation for the flow with larger particles, clearly featured in several flow statistics.

Particles at the wall with significant particle-to-fluid slip velocity create hot-spots of high local wall shear (on average about $4$ times larger than the mean) which contribute significantly to the mean wall-stress. These hot-spots show higher magnitude of the wall stress for small particles, which are, consistently, more localized. Their contribution to the mean wall shear is nearly the same for both particle sizes considered, about $70\%$. Also interesting to note is the highly increased probability of local shear stresses lower than the mean, and even negative (i.e., instantaneous flow reversal). These hot-spots change considerably the distribution of shear stresses, resulting in a pdf with a wide exponential tail and rms values, $\tau_{w}^{rms}\approx\tau_{w}$, in contrast to what found in canonical single-phase wall-bounded turbulent flows, where the pdf is well fitted by a log-normal distribution and $\tau_{w}^{rms}\approx 0.4\tau_{w}$. To quantify, the probability of finding a value of the shear stress of about $2$ times the total mean is negligible in the reference single-phase flow, whereas it becomes of the order of $10\%$ in the particle-laden flows. These findings are used to derive a scaling law for the wall particle-to-fluid apparent slip velocity as a function of the flow governing parameters.

Profiles of particle Reynolds number based on the particle slip velocity present values in the bulk of $O(1)$ for the case with large particles, and very close to zero for the suspension with smaller particles. This finite particle Reynolds number for case D20 can explain the significantly larger contribution of particle stresses to the total stresses with respect to the continuum limit reference. This can be attributed to the inertial-shear-thickening mechanism described in Picano et al. (2013), where finite inertia effectively increases the particles’ excluded volume.

Finite-size effects in the bulk of the flow are apparent in the second-order statistics, as shown by the profiles of the fluid velocity rms and Reynolds stresses. For these, we have tested the scaling arguments of Costa et al. (2016) and found that accounting for the particle-wall layer is a sufficient correction for the case with smaller particles, $D_{p}/\delta_{v}\approx 10$ (in fact, the finite size effects are relatively small in the bulk flow and using an effective viscosity only can still provide reasonable estimates). The flow with larger particles ($D_{p}/\delta_{v}\approx 20$), however, shows clear deviations from the continuum limit in the bulk of the channel also with the proposed rescaling, despite the fact that the correction for the particle-wall layer is sufficient for the lower-order statistics.

When investigating the particle dynamics, we compute statistics in wall-parallel bins of small wall-normal extent to take into account the effects of an inhomogeneous mean flow. Single-point dispersion statistics show that the spanwise particle dispersion coefficient is fairly independent of the wall-normal location, except very close to the wall. The spanwise particle dispersion is attributed to particle-turbulence interactions because of, first, the smaller diffusion coefficient near the wall and second, the larger diffusion coefficients for the smaller particles. The first effect is a consequence of the smaller spectrum of turbulence scales capable of disrupting the particle motion, whereas the second to the wider range of scales able to significantly displace the particles.

Conversely, the two-point dispersion statistics strongly depend on the wall-normal location and thus on the local shear rate. Higher shear rates induce larger relative velocities which result in faster dispersion. Even though the particles have a finite size, their dispersion statistics at short times show an exponential growth of the absolute displacement with time – as observed for point-particles with sub-Kolmogorov separation distances. The mechanism is, however, different. In this case short-range inter-particle interactions are likely the cause. Faster pair-dispersion for the larger particles is linked to the larger inter-particle interaction velocity.

Finally we investigated the particle collisional dynamics. Larger particles show higher mean values of relative velocity towards contact, consistent with the picture of \pcviscous-dominated, shear-induced contacts. This also explains the higher probability of finding pair of larger particles at close distance. Wall-normal variations in collision probability are therefore a consequence of the variations in local shear, and thus the larger particles collide more frequently than the small particles.

In this work, we have explored the advantages of massively parallel simulations, which allow for a multi-scale, three-dimensional and time-resolved picture of a system with well-defined physical parameters. This and similar studies show that the community have reached a point when simulations of interface-resolved particle-laden flows are possible. Such simulations, yet computationally expensive, can serve as valuable tool for validation of simpler two-way coupling algorithms and perhaps Eulerian models; we believe this type of simulations will spawn several investigations aiming at better models.

5 Acknowledgements

This work was supported by the Portuguese Foundation for Science and Technology under the grant No. SFRH/BD/85501/2012, by the European Research Council Grant No. ERC-2013-CoG-616186, TRITOS, by the Swedish Research Council (VR) and by COST Action MP1305: Flowing Matter. We acknowledge computer time provided by SNIC (Swedish National Infrastructure for Computing) and PRACE project 2014112543 for awarding us access to resource CURIE based in France at Genci/CEA.