2.1 Governing equations

The Navier–Stokes equations for a dilute polymer solution are,

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{f}=0, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}_{f}}{\unicode[STIX]{x2202}t}}+\boldsymbol{u}_{f}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}_{f}=-\unicode[STIX]{x1D735}p+{\displaystyle \frac{\unicode[STIX]{x1D6FD}}{\mathit{Re}}}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{f}+{\displaystyle \frac{(1-\unicode[STIX]{x1D6FD})}{\mathit{Re}}}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}+\boldsymbol{F}, & \displaystyle\end{eqnarray}$$

where subscript ‘ $f$ ’ identifies the fluid phase. The quantity $\unicode[STIX]{x1D6FD}\equiv \unicode[STIX]{x1D707}_{s}^{\ast }/(\unicode[STIX]{x1D707}_{s}^{\ast }+\unicode[STIX]{x1D707}_{p}^{\ast })$ is the ratio of the solvent to the total viscosities, the latter defined as the sum of the solvent $\unicode[STIX]{x1D707}_{s}^{\ast }$ and polymer $\unicode[STIX]{x1D707}_{p}^{\ast }$ contributions. In addition to the viscoelastic stress tensor $\unicode[STIX]{x1D64F}$ , a generic force field $\boldsymbol{F}$ is also included in (2.2). The viscoelastic stress is expressed in terms of the polymer conformation tensor $\unicode[STIX]{x1D658}$ , and using the FENE-P model,

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D658}}{\unicode[STIX]{x2202}t}}+\boldsymbol{u}_{f}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D658}=\unicode[STIX]{x1D658}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}_{f}+(\unicode[STIX]{x1D658}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}_{f})^{\text{T}}-\unicode[STIX]{x1D64F}, & \displaystyle\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D64F}={\displaystyle \frac{1}{Wi}}\left({\displaystyle \frac{\unicode[STIX]{x1D658}}{\unicode[STIX]{x1D713}}}-{\displaystyle \frac{\unicode[STIX]{x1D644}}{a}}\right);\quad \unicode[STIX]{x1D713}=1-{\displaystyle \frac{\text{tr}(\unicode[STIX]{x1D658})}{L_{max}^{2}}};\quad a=1-{\displaystyle \frac{3}{L_{max}^{2}}}, & \displaystyle\end{eqnarray}$$

where $L_{max}$ is the maximum extensibility of the polymer chains and $\unicode[STIX]{x1D644}$ refers to the unit tensor.

Particles are assumed to undergo rigid-body motion, and therefore no-slip boundary conditions are imposed at their surfaces $\unicode[STIX]{x2202}{\mathcal{V}}_{p}$ ,

(2.5) $$\begin{eqnarray}\displaystyle \boldsymbol{u}_{f}|_{\unicode[STIX]{x2202}{\mathcal{V}}_{p}}=\boldsymbol{u}_{p}+\unicode[STIX]{x1D74E}_{p}\times \boldsymbol{r}, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}_{p}$ and $\unicode[STIX]{x1D74E}_{p}$ are the particle translational and angular velocities, and $\boldsymbol{r}$ is the distance vector emanating from the centre of the sphere. The motion of the particles is governed by the Newton–Euler Lagrangian equations,

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}{\mathcal{V}}_{p}{\displaystyle \frac{\text{d}\boldsymbol{u}_{p}}{\text{d}t}}=\oint _{\unicode[STIX]{x2202}{\mathcal{V}}_{p}}\unicode[STIX]{x1D748}_{f}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}{\mathcal{A}}+\boldsymbol{F}_{c}, & \displaystyle\end{eqnarray}$$

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}{\mathcal{I}}_{p}{\displaystyle \frac{\text{d}\unicode[STIX]{x1D74E}_{\boldsymbol{p}}}{\text{d}t}}=\oint _{\unicode[STIX]{x2202}{\mathcal{V}}_{p}}\boldsymbol{r}\times \unicode[STIX]{x1D748}_{f}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}{\mathcal{A}}, & \displaystyle\end{eqnarray}$$

(2.8) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D748}_{f}={\displaystyle \frac{\unicode[STIX]{x1D6FD}}{\mathit{Re}}}[-p\unicode[STIX]{x1D644}+2\unicode[STIX]{x1D640}]+{\displaystyle \frac{1-\unicode[STIX]{x1D6FD}}{\mathit{Re}}}\unicode[STIX]{x1D64F}, & \displaystyle\end{eqnarray}$$

where ${\mathcal{I}}_{p}=\unicode[STIX]{x03C0}d_{p}^{5}/60$ and ${\mathcal{V}}_{p}=\unicode[STIX]{x03C0}d_{p}^{3}/6$ are the moment of inertia and the volume of a spherical particle, while $\unicode[STIX]{x1D748}_{f}$ is the fluid phase stress tensor. The force $\boldsymbol{F}_{c}$ results from particle–particle or particle–wall collisions and $\boldsymbol{n}$ is the outward unit vector at the particle surface. Lastly, $\unicode[STIX]{x1D640}$ is the strain-rate tensor of the fluid phase.

2.2 Numerical methodology

The solution of (2.1) and (2.2) is obtained numerically using a control-volume formulation and a fractional-step algorithm (Rosenfeld, Kwak & Vinokur Reference Rosenfeld, Kwak and Vinokur1991). The advection terms are treated explicitly using an Adams–Bashforth scheme, while an implicit Crank–Nicolson scheme is adopted for the treatment of the diffusion term and polymer stress.

The algorithm for the numerical solution of the conformation-tensor equations (2.3) was adopted in a number of high-fidelity simulations of instability waves and transition (Lee & Zaki Reference Lee and Zaki2017; Hameduddin, Gayme & Zaki Reference Hameduddin, Gayme and Zaki2019) as well as fully turbulent flows (Hameduddin et al. Reference Hameduddin, Meneveau, Zaki and Gayme2018; Hameduddin & Zaki Reference Hameduddin and Zaki2019). The equations are marched in time using a third-order accurate Runge–Kutta method where the advection and polymer stretch terms are treated using an Adams–Bashforth discretization within each sub-step. A semi-implicit method is adopted for the relaxation term, as suggested by Dubief et al. (Reference Dubief, Terrapon, White, Shaqfeh, Moin and Lele2005). A variant of slope-limiting approach for spatial derivatives proposed by Vaithianathan et al. (Reference Vaithianathan, Robert, Brasseur and Collins2006) is used to ensure numerical stability without introducing any artificial diffusion. In brief, the advection term is evaluated using second-order accurate central difference scheme as long as the conformation tensor remains positive definite; when its smallest eigenvalue approaches zero, the discretization is replaced by the one-sided difference which maximizes the minimum eigenvalue (for details, see appendix B by Hameduddin et al. Reference Hameduddin, Meneveau, Zaki and Gayme2018).

The no-slip boundary conditions at the particle surfaces are enforced using a variant of the sharp-interface, robust immersed boundary (IB) method described by Nicolaou, Jung & Zaki (Reference Nicolaou, Jung and Zaki2015). The method was adapted for freely moving particles in viscoelastic fluids and the algorithm was optimized to take advantage of the geometric simplicity of the particles for computational efficiency. The IB force is included on the right-hand side of (2.2), $\boldsymbol{F}=\boldsymbol{F}^{IB}$ , in order to satisfy the no-slip conditions (2.5). By integrating (2.2) over the volume of a particle and using the divergence theorem, the surface integral in (2.6) can be replaced by the volume integral (appendix B in Uhlmann Reference Uhlmann2005),

(2.9) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}{\mathcal{V}}_{p}{\displaystyle \frac{\text{d}\boldsymbol{u}_{p}}{\text{d}t}}=-\int _{{\mathcal{V}}_{p}}\boldsymbol{F}^{IB}\,\text{d}{\mathcal{V}}_{p}+{\displaystyle \frac{\text{d}}{\text{d}t}}\left(\int _{{\mathcal{V}}_{p}}\boldsymbol{u}_{f}\,\text{d}{\mathcal{V}}_{p}\right)+\boldsymbol{F}_{c}; & & \displaystyle\end{eqnarray}$$

similarly, (2.7) can be rewritten as,

(2.10) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70C}{\mathcal{I}}_{p}{\displaystyle \frac{\text{d}\unicode[STIX]{x1D74E}_{\boldsymbol{p}}}{\text{d}t}}=-\int _{{\mathcal{V}}_{p}}\boldsymbol{r}\times \boldsymbol{F}^{IB}\,\text{d}{\mathcal{V}}_{p}+{\displaystyle \frac{\text{d}}{\text{d}t}}\left(\int _{{\mathcal{V}}_{p}}\boldsymbol{r}\times \boldsymbol{u}_{f}\,\text{d}{\mathcal{V}}_{p}\right). & & \displaystyle\end{eqnarray}$$

Following Breugem (Reference Breugem2012), the second terms on the right-hand sides of equations (2.9) and (2.10) are evaluated over the fluid field occupying the solid domain. A short-range repulsive force proposed by Glowinski et al. (Reference Glowinski, Pan, Hesla, Joseph and Périaux2001) is used to treat the particle–particle and particle–wall collisions and the conformation tensor is set to unity inside the solid domain.

Our numerical tool has been extensively validated in Newtonian turbulence (Lee et al. Reference Lee, Jung, Sung and Zaki2013), moving-body problems in Newtonian fluids (Nicolaou et al. Reference Nicolaou, Jung and Zaki2015), and single-phase viscoelastic turbulence (Lee & Zaki Reference Lee and Zaki2017). In addition, we have verified the accuracy of our numerical tools in three reference configurations: a particle rotating in Newtonian and viscoelastic simple shear; lateral migration of an isolated particle in laminar Poiseuille flow; and Newtonian particle-laden turbulent flow in a channel. Details of these simulations are provided in appendix B.

Figure 1. Schematic of the particle-laden simulations. No-slip boundary conditions $\boldsymbol{u}_{f}=0$ are imposed at $y=\{0,2\}$ , and periodicity is enforced in the $x$ and $z$ directions.

Table 1. Physical and computational parameters of the simulations. Reynolds number $\mathit{Re}=h^{\ast }U_{b}^{\ast }/\unicode[STIX]{x1D708}^{\ast }$ and Weissenberg number $Wi=\unicode[STIX]{x1D706}^{\ast }U_{b}^{\ast }/h^{\ast }$ are based on the channel half-height $h^{\ast }$ , the bulk velocity $U_{b}^{\ast }$ , total kinematic viscosity $\unicode[STIX]{x1D708}^{\ast }$ and viscoelastic fluid relaxation time $\unicode[STIX]{x1D706}^{\ast }$ . The domain sizes are $L_{\{x,y,z\}}$ in the $\{x,y,z\}$ directions, and the number of grid cells is $N_{\{x,y,z\}}$ . The noted line types and symbols are adopted in figures, unless otherwise stated.

2.3 Flow configuration

The computational domain is a channel geometry with the streamwise, spanwise and wall-normal directions aligned with $x$ , $z$ and $y$ coordinates, as illustrated in figure 1. The bulk velocity $U_{b}^{\ast }$ and channel half-height $h^{\ast }$ are selected as the characteristic velocity $({\mathcal{U}}^{\ast }=U_{b}^{\ast })$ and length $({\mathcal{L}}^{\ast }=h^{\ast })$ for outer scaling. Where inner scaling is adopted, it is denoted by ‘ $+$ ’ symbol and, unless otherwise stated, is defined based on the friction velocity $u_{\unicode[STIX]{x1D70F}}^{\ast }\equiv \sqrt{\langle \unicode[STIX]{x1D70F}_{wall}^{\ast }/\unicode[STIX]{x1D70C}_{f}^{\ast }\rangle }$ and the length $l_{\unicode[STIX]{x1D70F}}^{\ast }\equiv \unicode[STIX]{x1D708}^{\ast }/u_{\unicode[STIX]{x1D70F}}^{\ast }$ , where $\unicode[STIX]{x1D70F}_{wall}^{\ast }$ is the total wall shear stress. The flow is maintained at constant mass flux, no-slip boundary conditions are imposed at the bottom and top walls $y=\{0,2\}$ , and periodicity is enforced in the two horizontal directions.

Physical and computational parameters of the simulations are provided in table 1. One Newtonian and two viscoelastic conditions are simulated, each with and without the presence of neutrally buoyant spherical particles. A designation is introduced in table 1 to identify each simulation: letter W is followed by the Weissenberg number $Wi$ and P by the bulk solid volume fraction of the system $\unicode[STIX]{x1D711}(\%)$ in the associated case. In all simulations the grid cells are uniformly distributed in the streamwise and spanwise directions. In the wall-normal direction, the grid spacing is uniform in the particle-laden conditions and a hyperbolic tangent grid stretching is adopted in the single-phase simulations. In the particle-laden flows, the grid size is dictated by the required number of grid cells to correctly predict the flow field in the vicinity of the particles. Following previous studies (Goyal & Derksen Reference Goyal and Derksen2012; Picano et al. Reference Picano, Breugem and Brandt2015), the ratio of particle diameter to the grid cell size is set to $d_{p}/\unicode[STIX]{x0394}x=16$ .

The single-phase Newtonian simulation is initialized using a superposition of laminar Poiseuille flow and small-amplitude random fluctuations which trigger breakdown to turbulence, while turbulence is triggered in the particle-laden case naturally due to the presence of the randomly seeded particles. The viscoelastic conditions are initialized with the velocity and pressure fields from their Newtonian counterparts in the statistically stationary state. Beyond the initial transients, ensemble-averaged statistics of any mixture quantity $\langle o\rangle$ were evaluated, conditioned on the fluid and particle phases, and denoted $\langle o_{\{\,f,p\}}\rangle _{\{\,f,p\}}$ . The unconditional average can be reconstructed,

(2.11) $$\begin{eqnarray}\displaystyle \langle o\rangle =\langle (1-\unicode[STIX]{x1D712})o_{f}\rangle +\langle \unicode[STIX]{x1D712}o_{p}\rangle =(1-\unicode[STIX]{x1D719})\langle o_{f}\rangle _{f}+\unicode[STIX]{x1D719}\langle o_{p}\rangle _{p}, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D712}=\{0,1\}$ is a phase indicator of the fluid and solid, and $\unicode[STIX]{x1D719}=\langle \unicode[STIX]{x1D712}\rangle$ is the solid volume fraction. Averaging is performed in the homogeneous streamwise and spanwise directions, as well as in time. Note that the subscripts $\langle \bullet \rangle _{\{f,p\}}$ are henceforth dropped and conditional averaging is implied. Fluctuations in any fluid or particle quantities are defined relative to their respective means, $o_{\{f,p\}}^{\prime }=o_{\{f,p\}}-\langle o_{\{f,p\}}\rangle$ . Statistics were collected for sufficiently long periods (e.g. 1800 convective times for case W15P5) in order to ensure convergence which was verified by comparing results from half and the total number of samples.

3.1 Drag modulation and stress budget

The friction Reynolds number $\mathit{Re}_{\unicode[STIX]{x1D70F}}\equiv \mathit{Re}\,u_{\unicode[STIX]{x1D70F}}^{\ast }/U_{b}^{\ast }$ is reported in figure 2. Particle-laden cases systematically experience higher $Re_{\unicode[STIX]{x1D70F}}$ relative to their single-phase counterparts. A similar trend was observed in previous studies of Newtonian flows laden with dilute and dense suspensions of spherical particles (Picano et al. Reference Picano, Breugem and Brandt2015; Lashgari et al. Reference Lashgari, Picano, Breugem and Brandt2016; Costa et al. Reference Costa, Picano, Brandt and Breugem2018). Those studies generally concur that the particle stress is responsible for the drag enhancement. A curious observation from figure 2 is that the particle drag enhancement is higher in the low Weissenberg case ( $15.8\,\%$ for W15P5), in comparison to both the Newtonian ( $5.5\,\%$ for W0P5) and high Weissenberg ( $6.5\,\%$ for W25P5) cases. In addition, the effectiveness of the polymer to reduce drag in the low Weissenberg particle-laden case is appreciably diminished in the particle laden case, $\mathit{Re}_{\unicode[STIX]{x1D70F}}^{\text{W15}}/\mathit{Re}_{\unicode[STIX]{x1D70F}}^{\text{W0}}<\mathit{Re}_{\unicode[STIX]{x1D70F}}^{\text{W15P5}}/\mathit{Re}_{\unicode[STIX]{x1D70F}}^{\text{W0P5}}$ . At the high Weissenberg condition, however, the propensity of viscoelasticity to reduce drag appears largely unaffected by the presence of the particle phase.

The above observations will be examined in detail and explained, starting by an assessment of the different contributions to the total stress $\unicode[STIX]{x1D70F}_{tot}$ at every wall-normal position,

(3.1) $$\begin{eqnarray}\displaystyle & & \displaystyle \overbrace{-[(1-\unicode[STIX]{x1D719})\langle u_{f}^{\prime }v_{f}^{\prime }\rangle +\unicode[STIX]{x1D719}\langle u_{p}^{\prime }v_{p}^{\prime }\rangle ]}^{\unicode[STIX]{x1D70F}_{Re}}+\overbrace{{\displaystyle \frac{\unicode[STIX]{x1D6FD}}{Re}}(1-\unicode[STIX]{x1D719}){\displaystyle \frac{\text{d}\langle u_{f}\rangle }{\text{d}y}}}^{\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D707}}}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\underbrace{{\displaystyle \frac{(1-\unicode[STIX]{x1D6FD})}{Re}}(1-\unicode[STIX]{x1D719})\langle \unicode[STIX]{x1D61B}_{xy}\rangle }_{\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D6FD}}}+\underbrace{\unicode[STIX]{x1D719}\langle \unicode[STIX]{x1D70E}_{p,xy}\rangle }_{\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D719}}}=\underbrace{\langle \unicode[STIX]{x1D70F}_{wall}\rangle (1-y)}_{\unicode[STIX]{x1D70F}_{tot}}.\end{eqnarray}$$

From left to right, these contributions are the turbulent Reynolds stress $\unicode[STIX]{x1D70F}_{Re}$ , the viscous stress $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D707}}$ , the polymer stress $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D6FD}}$ and the particle stress $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D719}}$ . Figure 3 compares the contribution of each term for different cases in outer (left axis) and inner (right axis) scaling. As anticipated for the single-phase viscoelastic flows, $\unicode[STIX]{x1D70F}_{Re}$ is significantly attenuated and the incurred viscoelastic stress $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D6FD}}$ is relatively small (figures 3(a), 3(b) and 3(c)); the net effect is a decrease of the total stress. The addition of particles to the Newtonian fluid leads to an induced stress which plays a major role in drag enhancement, while the share of turbulent Reynolds stress is almost unchanged relative to W0. In contrast, the addition of particles appreciably reduces the contribution of $\unicode[STIX]{x1D70F}_{Re}$ to the total stress in the viscoelastic cases – an effect that may seem counter-intuitive at first given that the total drag has increased in W15P5 and W25P5 relative to W15 and W25. The drastic increase in the polymer stresses explains the overall increase in drag in the presence of the particles. In cases W15P5 and W25P5, $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D6FD}}$ is comparable to the viscous stress contribution near the wall and takes a significant share of the total stress throughout the channel. On balance, the attenuation of the turbulence shear stress and the polymer stress enhancement in the vicinity of the particles result in a less effective drag reduction in case W15P5 relative to W15, and an almost equal level of drag reduction in cases W25P5 and W25.

The pronounced increase of polymer stress in the particle-laden viscoelastic configurations arises due to the polymer stretching mechanism which differs from the single-phase flow. In the latter, polymer chains are stretched by the turbulent fluctuations. In the two-phase configuration, particles alter the velocity field in their vicinity and cause local stretching of polymer chains and enhancement of the associated stresses (figure 4); this effect is examined in detail in § 4.1.

Figure 4. Instantaneous snapshot of particle positions, contours of $u_{f}$ (vertical plane) and $u_{f}-\langle u_{f}\rangle$ (horizontal plane) for case W15P5. For clarity, only 20 % of particles are shown. Marked area is magnified and reproduced with contours of $\unicode[STIX]{x1D61B}_{xy}$ .

The turbulent stresses $\unicode[STIX]{x1D70F}_{Re}$ in the reference cases W15 and W25 are within $25\,\%$ of one another. And while the addition of the particle phase reduces $\unicode[STIX]{x1D70F}_{Re}$ at both elasticities, it is much more influential at the higher Weissenberg number: W15P5 experiences a twofold reduction compared to fourfold in case W25P5. This prominent dissimilarity is a sign of intrinsic differences in particle–turbulence interactions in W15P5 and W25P5, and calls for further scrutiny of the time-dependent effects that contribute to this average stress (§ 3.3). Lastly, in all particle-laden cases, a peak in the particle stress is recorded near $y=0.1$ , approximately one diameter from the wall, where the particles are expected to experience high relative velocity with respect to the fluid. However, the contribution of the particle stress is more evenly distributed throughout the channel in cases W15P5 and W25P5, relative to the Newtonian W0P5 where it is most relevant in the near-wall region. The persistent contribution in the outer flow in the viscoelastic cases reinforces a potential change in the role of particles in the wall-normal momentum transport.

Figure 5. Mean streamwise fluid velocity in (a) outer and (b) inner units. Line types as listed in table 1: (

) W0; (

) W15; (

) W25; (

) W0P5; (

) W15P5; (

) W25P5. (a) Marked area is magnified in the inset. (b) Dashed line (

) is the Virk asymptote (Virk et al. Reference Virk, Mickley and Smith1970); dotted line (

) Newtonian law of the wall $\langle u_{f}^{+}\rangle =2.5\log (y^{+})+5.0$ ; dashed dotted line (

) laminar Poiseuille profile $\langle u_{f}^{+}\rangle =y^{+}$ .

3.2 Mean fluid velocity and Reynolds stress profiles

The mean streamwise fluid velocity profiles are reported in figure 5, in outer and inner scaling. The characteristic velocity for the latter scaling is based on the viscous component of the wall shear stress, i.e.  $u_{\unicode[STIX]{x1D707}}^{\ast }\equiv \sqrt{\langle \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D707},wall}^{\ast }/\unicode[STIX]{x1D70C}_{f}^{\ast }\rangle }$ . In the viscoelastic cases, $\langle u_{f}\rangle$ deviates from the Newtonian turbulent channel-flow profile and approaches the laminar solution. The deviation from the Newtonian turbulent profile is a signature of suppressed wall-ward momentum mixing. As described by White et al. (Reference White, Dubief and Klewicki2018), polymer stress takes over the role of Reynolds stress with respect to inertial mixing, but is less effective; a reduced level of momentum mixing in the wall-normal direction eventually results in the modification of the law of the wall. The addition of particles intensifies this effect to the point where the profile from W25P5 tends to Virk’s maximum drag reduction asymptote (Virk et al. Reference Virk, Mickley and Smith1970). This shift can be an indication that the particles effectively suppress the turbulent fluctuations and, in turn, the associated momentum mixing. Nevertheless, the impact of particles on the mean streamwise velocity profile is almost negligible in W0P5. Previous studies have shown that in semi-dilute particle-laden Newtonian flows, except in the near-wall region, particles render the mean streamwise velocity profile similar to the laminar case (Picano et al. Reference Picano, Breugem and Brandt2015; Costa et al. Reference Costa, Picano, Brandt and Breugem2018). In the near-wall region, however, particles are known to accelerate the fluid velocity due to their finite slip (Shao et al. Reference Shao, Wu and Yu2012). This effect is weakly apparent in case W0P5 and overshadowed by the strong modification of the mean streamwise velocity profile due to the turbulence attenuation in cases W15P5 and W25P5 (inset of figure 5(a)).

Figure 6. Components of the Reynolds stress tensor, where $u_{f}^{\prime }$ , $v_{f}^{\prime }$ and $w_{f}^{\prime }$ are the fluctuations in the streamwise, wall-normal and spanwise fluid velocities. Line types as listed in table 1: (

) W0; (

) W15; (

) W25; (

) W0P5; (

) W15P5; (

) W25P5.

The impact of the particle phase on the turbulence fluctuations is captured by the Reynolds stresses (figure 6). First recall that viscoelasticity alone is known to swell the streamwise streaks thus resulting in an upward shift of the peak of streamwise fluctuations; in addition, relative to Newtonian flows, the peak value of the streamwise stress increases at low/moderate elasticity and reduces at high elasticity (Dallas et al. Reference Dallas, Vassilicos and Hewitt2010; Xi & Graham Reference Xi and Graham2012). Figure 6(a) shows that the peak streamwise stresses for W15 and W25 shift away from the wall, without any significant change in its value. Particles weaken the streamwise vortices and decrease the peak streamwise velocity fluctuations in the Newtonian and viscoelastic cases. This effect, associated with the particle slip velocity acting as a sink of kinetic energy in the coherent streaks, is substantially intensified in the viscoelastic flows. For the wall-normal and spanwise stresses, however, the picture is quite divergent. In the near-wall region, where the particles are the main source of cross-flow fluctuations due to their finite slip velocity, $\langle v_{f}^{\prime }v_{f}^{\prime }\rangle$ and $\langle w_{f}^{\prime }w_{f}^{\prime }\rangle$ are higher in all particle-laden cases relative to the single-phase counterparts. However, qualitative differences in the impact of particles on the Newtonian and viscoelastic flows arise away from the wall: while in W0P5, particles slightly enhance the cross-stream fluctuations, in W15P5 and W25P5 they attenuate the fluctuations in most of the channel. In the Newtonian case, the particles redistribute momentum from the streaks to the cross-stream directions by inducing particle-scale vortices (Shao et al. Reference Shao, Wu and Yu2012). Since the streaks are appreciably weakened in viscoelastic flows and, as a result, this effect is muted in W15P5 and W25P5, the cross-stream fluctuations are also reduced. Moreover, in contrast to the Newtonian case the presence of particles in viscoelastic flows reduces, rather than enhances, the turbulent shear stress – a matter that we examine in detail in § 3.3.

Figure 7. Time series of $\mathit{Re}_{\unicode[STIX]{x1D70F},xz}$ . Line types as listed in table 1: (

) W0; (

) W15; (

) W25; (

) W0P5; (

) W15P5; (

) W25P5.

3.3 Transient dynamics and conditional statistics

The regeneration cycle of wall turbulence involves (i) formation of streamwise-elongated streaks, (ii) breakdown of the streaks and (iii) generation of streamwise vortices which enhance the three-dimensionality of the flow and are followed by a surge in the wall shear rates. In Newtonian fluids, particles reinforce the generation of streamwise vortices and promote the wall shear stress (Wang, Abbas & Climent Reference Wang, Abbas and Climent2018). On the other hand, viscoelasticity alters the cycle by promoting extended periods of low shear rates known as ‘hibernating’ states, which sporadically give way to bursts of high shear-rate ‘active’ states of turbulence. Xi & Graham (Reference Xi and Graham2012) showed that once viscoelasticity surpasses a critical value, the stretching of polymer chains weakens the turbulent fluctuations and in turn contributes to the dominance of the hibernating state. Due to the lower shear stress associated with hibernating turbulence, polymer solutions experience a lower turbulent drag in an average sense. In an effort to determine the influence of the particle phase on this process, we examine the transient dynamics in our simulations within a similar framework.

Figure 7 contrasts the time evolutions of the wall-averaged instantaneous friction Reynolds numbers, $\mathit{Re}_{\unicode[STIX]{x1D70F},xz}\equiv \mathit{Re}\sqrt{(1/\unicode[STIX]{x1D70C}_{f}^{\ast })\langle \unicode[STIX]{x1D70F}_{wall}^{\ast }\rangle _{xz}}/U_{b}^{\ast }$ , from the different simulations. In both Newtonian flows, W0 and W0P5, very minimal fluctuations are observed which is indicative of an active turbulence state. While the hibernating state also exists in Newtonian flows (Xi & Graham Reference Xi and Graham2010), it is relatively rare and its influence is overshadowed by the strong dominance of active turbulence. In the low Weissenberg cases, W15 and W15P5, the fluctuating behaviour of $\mathit{Re}_{\unicode[STIX]{x1D70F},xz}$ signals a competition between active and hibernating states. Infrequent peaks of $\mathit{Re}_{\unicode[STIX]{x1D70F},xz}$ in W15P5 indicate that the hibernating state prevails. At high Weissenberg, W25 demonstrates a cyclic behaviour akin to that of W15, but W25P5 appears dominated by the hibernating state. In both W15P5 and W25P5, the particles complement the influence of viscoelasticity by promoting hibernation, or suppressing active turbulence. The drag itself is, nonetheless, increased in presence of the particles due to the increase in the polymer stress, as noted in § 3.1.

Figure 8. Wall-normal profiles of the overall and conditional averages of streamwise velocity for (a) W15 and (b) W15P5. Averages are computed over the homogeneous streamwise and spanwise directions and (black lines,

) the full sampling time, (blue lines,

) hibernating states and (red lines,

) active states. Each curve is normalized by its respective $u_{\unicode[STIX]{x1D707}}^{\ast }$ .

Figure 9. Wall-normal profiles of the overall and conditional averages of contributions to the shear stress (3.1). Averages are computed over the homogeneous streamwise and spanwise directions and (black lines,

) the full sampling time, (blue lines, 

) hibernating states and (red lines,

) active states. (a–c) case W15; (d–f) case W15P5.

Figure 10. Particle migration during the hibernating and active turbulence states in case W15P5. Time evolution of contours of streamwise- and spanwise-averaged (top) particle wall-normal velocity, (middle) phase indicator and (bottom) $\mathit{Re}_{\unicode[STIX]{x1D70F},xz}$ .

The mean streamwise velocity profiles were evaluated conditioned on the state of the turbulence, and the hibernating and active states are contrasted to the conventional mean in figure 8, for W15 and W15P5. In line with the observations by Wang, Shekar & Graham (Reference Wang, Shekar and Graham2017), at both Weissenberg numbers the profile of the hibernating state is closer to Virk’s asymptote, and that of active turbulence is closer to Newtonian turbulence. In W15P5, the overall average is biased towards the hibernating state-conditioned profile, in agreement with the scarcity of the active turbulence instances in figure 7. Figure 9 compares different components of the stress budget during the hibernating and active turbulence states, for W15 and W15P5. Without and with particles, during the hibernating state the Reynolds stress $\unicode[STIX]{x1D70F}_{Re}$ is reduced and the viscous stress $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D707}}$ increased. Similar to previous studies of single-phase viscoelastic turbulence (Xi & Graham Reference Xi and Graham2010; Tamano, Graham & Morinishi Reference Tamano, Graham and Morinishi2011; Wang et al. Reference Wang, Shekar and Graham2017), the share of the polymer stress in W15 increases during the active state. The effect has been interpreted in terms of the polymer chains being stretched and relaxed during the active and hibernating states, respectively, giving rise to the similar trend in polymer and turbulent stresses (Wang et al. Reference Wang, Shekar and Graham2017). Interestingly, however, $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D6FD}}$ exhibits the opposite behaviour in W15P5: in most of the channel the contribution of the polymer stress is slightly higher during the hibernating state relative to the active one. This change in character will be ascribed to the presence of particles, which will be demonstrated in § 4.1 to be the primary cause of polymer stretch. In other words, $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D6FD}}$ in W15P5 is inherently different from its single-phase counterpart.

The importance of the particle phase becomes evident when their behaviour is examined as the turbulence shifts from the hibernating to active state. Figure 10 shows the time dependence of the spatially averaged particle velocity $\langle v_{p}\rangle _{x,z}$ and phase indicator $\langle \unicode[STIX]{x1D712}\rangle _{x,z}$ during this shift. Particles gradually migrate towards the channel centre during the hibernating state. Due to the subdued turbulence activity, momentum mixing is weak and the mean streamwise velocity tends towards the parabolic profile. As a result, the particles’ migration towards the channel centre is similar in nature to their lateral migration in laminar viscoelastic Poiseuille flow (see e.g. Huang et al. Reference Huang, Feng, Hu and Joseph1997; Li, McKinley & Ardekani Reference Li, McKinley and Ardekani2015, and appendix B). As they migrate to the centre, their streamwise velocity will increase due to the higher local momentum; however their fluctuating motion is attenuated during the hibernating state due to the weaker turbulence activity. Migration towards the channel centre continues until the near-wall region is almost devoid of particles, $\langle \unicode[STIX]{x1D712}\rangle _{xz}<2\,\%$ . At this stage, active turbulence is reinitiated and the particles are pushed back towards the wall by the fluctuating motions. As a result, the streamwise fluid velocity profile is flattened by turbulent mixing, and particle fluctuations are once again revived.

Particle migration is key to understanding the increased contribution of the polymer stress during hibernation (figure 9(f)). Upon closer examination of figures 9(f) and 10, we observe that this increase takes place only where the particles are more concentrated, i.e. in the outer region ( $y>0.2$ ) during hibernation. In fact, near the wall, the contribution of the polymer stress is larger during the active state due to the coexistence of turbulence-induced polymer stretch and the return of the particles towards the wall as the turbulence is revived.

4.1 Particle dynamics

Figure 11. Wall-normal profiles of (a) solid volume fraction and (b) mean streamwise particle and fluid velocities. Line types as listed in table 1: (

) W0P5; (

) W15P5; (

) W25P5. (b) Dashed lines correspond to the fluid velocities; Inset shows the mean relative particle velocities.

The focus in now placed on the particle dynamics, and its connection to the Newtonian and viscoelastic flow statistics. The solid volume fraction profiles are reported in figure 11(a). In the Newtonian case, a flat profile in the outer region signals an even distribution of particles, due to strong turbulent mixing. In viscoelastic conditions, the solid volume fraction increases away from the wall, and at the channel centre of case W25P5 exceeds twice its bulk value ( $\unicode[STIX]{x1D711}=5\,\%$ ). This finding is in line with the particle migration mechanism explained in § 3.3. Note that this migration behaviour is due to fluid elasticity. It is therefore inherently different from the phenomenologically similar migration towards the channel centre that can take place in Newtonian turbulence, at high solid volume fractions $\unicode[STIX]{x1D719}\geqslant 30\,\%$ or density ratios $\unicode[STIX]{x1D70C}_{p}^{\ast }/\unicode[STIX]{x1D70C}_{f}^{\ast }\approx 10$ (Fornari et al. Reference Fornari, Formenti, Picano and Brandt2016; Lashgari et al. Reference Lashgari, Picano, Breugem and Brandt2016); that effect is due to a particle normal stress imbalance in the wall-normal direction. A stabilized layer of particles is formed near the wall in cases W15P5 and W25P5, and is due to the particle–wall lubrication effect and asymmetric interactions in that region. The same mechanism is present in Newtonian particle-laden flows at higher solid volume fractions (Yeo & Maxey Reference Yeo and Maxey2011; Lashgari et al. Reference Lashgari, Picano, Breugem and Brandt2016), and is observed here in the viscoelastic cases despite a highly dilute concentration of particles near the wall.

Figure 12. Probability density function (PDF) of normalized (a,b) perturbation shell velocity and (c,d) slip velocity. The median of each PDF is marked by a vertical dashed line. Line types as listed in table 1: (

) W0P5; (

) W15P5; (

) W25P5. The shell and slip velocities are defined in (4.1) and (4.2). Data are conditionally sampled for particles at (a,c) $y=1/9$ and (b,d) $y=1/3$ .

The mean streamwise particle and fluid velocity profiles are plotted in figure 11(b). In each case, the mean fluid and particle velocity profiles are in good agreement near the core of the channel; however, the particle streamwise velocity is larger than that of the fluid in close vicinity to the wall ( $y\leqslant 0.5d_{p},~y^{+}\leqslant 10$ ). A positive particle slip, loosely defined as the particle velocity relative to the local fluid velocity, is evident in the viscous sublayer and is due to the particles’ finite rigid-body motion while the fluid velocity vanishes. Surprisingly, a negative slip velocity is observed near $0.75d_{p}\leqslant y\leqslant 1.5d_{p}$ in W15P5 and W25P5 (inset of figure 11(b)). Its interpretation must, however, be carefully considered in light of the heterogeneity of the particle distribution which can introduce bias if the simple definition of slip as $\langle u_{p}\rangle -\langle u_{f}\rangle$ is adopted.

A more precise definition of the velocity experienced, or seen, by the particles is required in order to estimate the slip velocity. We adopt the definition proposed by Kidanemariam et al. (Reference Kidanemariam, Chan-Braun, Doychev and Uhlmann2013): the instantaneous fluid velocity is integrated over a spherical shell, ${\mathcal{S}}$ , with diameter equal to $3d_{p}$ and centred around the particle. In order to avoid sampling bias due to heterogeneity in $y$ , the shell surface is trimmed by two wall-parallel planes located at one radius below and above the particle position. Formally, the shell and particle slip velocities are defined as,

(4.1) $$\begin{eqnarray}\displaystyle & \displaystyle u_{{\mathcal{S}}}={\displaystyle \frac{\displaystyle \oint _{{\mathcal{S}}}(1-\unicode[STIX]{x1D712})u_{f}\,\text{d}{\mathcal{A}}}{\displaystyle \oint _{{\mathcal{S}}}(1-\unicode[STIX]{x1D712})\,\text{d}{\mathcal{A}}}}, & \displaystyle\end{eqnarray}$$

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle u_{slip}=u_{p}-u_{{\mathcal{S}}}. & \displaystyle\end{eqnarray}$$

A shell perturbation velocity can be defined relative to the local mean fluid velocity, $u_{{\mathcal{S}}}-\langle u_{f}\rangle$ . Probability density functions (PDFs) of this quantity are reported in figure 12 at two heights from the wall, normalized by $\langle u_{f}\rangle$ . According to figure 12(a), a wide range of positive and negative shell perturbation velocities are experienced by the particles in both the Newtonian and viscoelastic cases, and a slight bias towards negative values is present for the latter. This trend signals a preferential concentration of particles in low-speed streaks in the buffer layer ( $15\leqslant y^{+}\leqslant 21$ ), which is precipitated by sweeping high-speed vortices that displace the particles into low-speed regions. This phenomenon has been reported for Newtonian wall turbulence, in the case of point particles (Kaftori et al. Reference Kaftori, Hetsroni and Banerjee1994; Pan & Banerjee Reference Pan and Banerjee1997) and small particles with size $d_{p}^{+}\approx 5{-}10$ (Kidanemariam et al. Reference Kidanemariam, Chan-Braun, Doychev and Uhlmann2013; Wang et al. Reference Wang, Abbas and Climent2018), but has not been observed for large particles $(d^{+}\approx 30)$ (Hetsroni & Rozenblit Reference Hetsroni and Rozenblit1994). Since turbulence structures, in particular streaks, are inflated by viscoelasticity and thickened in the enlarged buffer layer, the relative particle size is smaller than in the Newtonian configuration. As a result, the particles further adapt to the fluid motion induced by the near-wall streamwise vortices, preferentially residing in the low-speed streaks. In the log layer ( $35\leqslant y^{+}\leqslant 63$ ), the distribution of normalized shell velocity perturbation is centred at zero and predictably narrowed (figure 12(b)), in both the Newtonian and viscoelastic cases.

Figure 12 also reports PDFs of the normalized slip velocity, which are also broad in the buffer layer and narrower in the log layer. In the former region (figure 12(c)), viscoelasticity enhances the probability of particles experiencing negative slip. Together figures 12(a) and 12(c) clarify that the negative relative mean velocities reported in figure 11(b) in the viscoelastic cases stem from the combination of two effects: (i) preferential location of particles in the low-speed streamwise velocity zones and (ii) mean negative slip velocity of particles with respect to their shell velocity. In the log layer (figure 12(d)), the PDFs of normalized slip velocity are centred at zero. Although narrow relative to their counterparts in the buffer region, the widths of all distributions are still non-negligible particularly when compared to the normalized turbulence fluctuations. Specifically, in all particle-laden cases the variance of $u_{slip}/\langle u_{f}\rangle$ is comparable in magnitude to $\sqrt{\langle u_{f}^{\prime }u_{f}^{\prime }\rangle }/\langle u_{f}\rangle$ . This observation is notable for two reasons: (i) it shows that despite the zero average slip velocity, particles still act as effective sinks of kinetic energy in the log layer, which explains their significant role in the reduction of the fluid streamwise velocity fluctuations (figure 6); (ii) a non-zero slip velocity gives rise to the stretching of polymer chains in the neighbourhood of the particles.

Figure 13. PDF of the normalized trace of the conformation tensor, sampled for fluid points at $y=0.15$ . Solid grey line (

): W15; black line (

): W15P5; black dashed line (

): sampled points are also constrained to be within a spherical region concentric with particles and twice their diameter in case W15P5.

In § 3.1, the pronounced increase of polymer stress in the presence of particles was ascribed to stretching of the polymer in their vicinity (figure 3). This assertion is supported by comparing the PDFs of the normalized trace of the conformation tensor, $\text{tr}(\unicode[STIX]{x1D658})/L_{max}^{2}$ , in cases W15 and W15P5 (figure 13). The comparison is made at $y=0.15$ , where the polymer stress contribution differs dramatically in the two cases (compare figures 3(b) and 3(e)). The evident dissimilarity between the shapes of the distributions in figures 13 highlights the qualitative difference in polymer stretching in the single-phase and particle-laden cases. In the absence of particles, the stretching of polymer chains, measured by the trace of the conformation tensor, has nearly a Gaussian distribution. In contrast, in the particle-laden case the distribution is highly skewed towards more stretched polymer chains. By exclusively considering fluid points located in the vicinity of the particles, specifically within spheres concentric with the particles and twice their diameter, the distribution is further skewed towards higher values, which points to the role of particles in the stretching of polymer chains.

Figure 14. Comparison of particle and fluid velocity fluctuations in the (a) streamwise, (b) wall-normal and (c) spanwise directions. (d) Comparison of particle mean rotation in the spanwise direction and fluid vorticity; inset depicts the ratio of the two quantities in the near-wall region. Line types as listed in table 1: (

) W0P5; (

) W15P5; (

) W25P5. Quantities in the fluid phase are represented by dashed lines.

Statistics of fluctuations in the particle and fluid velocities are plotted in figure 14; the fluid curves are reproduced from figure 6. In the streamwise direction, a non-monotonic behaviour is observed in the outer region: low viscoelasticity enhances the particle fluctuations beyond the Newtonian results but higher Weissenberg number reverses this effect and yields lower $\langle u_{p}^{\prime }u_{p}^{\prime }\rangle$ than Newtonian. This trend parallels the effect on fluid fluctuations, which was discussed in connection with figure 6: elasticity expands the turbulent structures thus shifting the peak $\langle u_{f}^{\prime }u_{f}^{\prime }\rangle$ towards the core of the channel, and ultimately attenuates the magnitude of turbulent fluctuations.

In the Newtonian fluid, all three components of the particle fluctuations are lower than the fluid ones throughout most of the channel, and the trend reverses near the wall. For the streamwise and spanwise components the trend reverses in a very small near-wall layer, and for the wall-normal component the reversal appears at wall distances below one diameter. In contrast to the Newtonian case, in the viscoelastic fluids all three components of the particle velocity fluctuations are commensurate with or exceed the fluid ones. The streamwise component exceeds the fluid values over a thick near-wall region of similar size to the particle diameter; the wall-normal particle fluctuations are notably higher than the fluid ones over a large part of the channel; and the spanwise fluctuations are essentially equal to the fluid ones except in the close vicinity of the wall where the particle slip velocity is finite.

Figure 15. Particle fluctuations in (a) streamwise and (b) wall-normal directions in case W15P5. Particle averages were performed using (

) full rigid-body motion and (

) motion of particle centres only; (

) fluid fluctuations.

Since the mean vorticity vector is aligned with the spanwise coordinate, the particles sustain a mean angular velocity in that direction. The two quantities, namely the spanwise vorticity and particle angular velocity, are compared in figure 14(d). The impact of viscoelasticity on the vorticity distribution is also observed in the particle angular velocities, which are increased in the bulk of the channel relative to the Newtonian case. When normalized, $\langle \unicode[STIX]{x1D714}_{p,z}\rangle /\langle \unicode[STIX]{x1D714}_{f,z}\rangle$ is similar for all three cases within the region $y\leqslant 0.2$ (inset of figure 14(d)).

The particles’ mean spanwise angular velocity strongly affects their velocity fluctuations in the $x$ and $y$ directions. This effect is isolated by comparing the velocity fluctuations of the particle centres to those of the particles in full rigid-body motion (figure 15); only case W15P5 is considered for brevity since the results are qualitatively similar to the other conditions. The effect of particles’ rotation on their streamwise velocity fluctuations appears minimal (figure 15(a)). In the wall-normal component, however, the impact is much more pronounced: the fluctuations of the particle centres are at the same level as those in the fluid phase in most of the channel, and the large peak observed in the full rigid-body results near $y\approx 0.08$ nearly vanishes (figure 15(b)). Particle rotation thus plays an important role in establishing this peak when the full rigid-body motion is considered. In both the streamwise and wall-normal directions, the fluctuations of the particle centres have small local maxima near $y\approx 0.17$ , which are due to the bouncing motion when particles with high inertia collide with the stabilized layer of particles close to the wall. This assertion is supported by a large peak in the profile of fluctuations in the collision force at the same $y$ position (not shown). Moreover, particles leaving the near-wall stabilized layer experience a rapid increase in their streamwise velocity, due to the large mean gradient in the region $0.1\leqslant y\leqslant 0.2$ (see figure 11(b)), which accounts for the maximum in the fluctuations of the particle centres.

4.2 Particle clustering and microstructure analysis

Even in the simple case of two settling spheres in a quiescent flow, the attraction/repulsion of a particle pair in a viscoelastic fluid depends on various factors including elasticity, inertia, shear thinning/thickening, initial distance and configuration. Despite these complex features, both experiments and simulations confirm that particle aggregation is enhanced by normal stresses induced by elasticity (cf. Zenit & Feng Reference Zenit and Feng2018, for detailed review). Since particle clustering and wake dynamics, even at low concentrations, significantly affect the turbulence (Uhlmann & Doychev Reference Uhlmann and Doychev2014; Capecelatro, Desjardins & Fox Reference Capecelatro, Desjardins and Fox2018), we dedicate this section to the study of particle microstructures in our simulations.

The radial distribution function (RDF) is a common first-order measure of particle clustering in suspensions (Shaw Reference Shaw2003). It is defined as the average number of particles $\langle \text{d}N\rangle$ located at a spherical shell of radius $r$ and thickness $\text{d}r$ , and which is centred at the reference particle position $\boldsymbol{P}_{ref}$ . In a homogeneous system of point particles, the RDF is commonly normalized by the same quantity for a random distribution of points in the space (Reade & Collins Reference Reade and Collins2000). In our system, however, due to the presence of a strong heterogeneity in the $y$ direction, as well as the non-overlapping condition for finite-sized particles, an unbiased analysis calls for a tailored normalization. We define $g(r)$ in the limit $\text{d}r\rightarrow 0$ by,

(4.3) $$\begin{eqnarray}\displaystyle & \displaystyle g(r)=\langle \text{d}N(r)\rangle /\langle \,\text{d}N(r)^{R}\rangle , & \displaystyle\end{eqnarray}$$

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle \text{d}N(r)=\mathop{\sum }_{m=1}^{N_{p}}\unicode[STIX]{x1D6FF}(r-|r_{m}|), & \displaystyle\end{eqnarray}$$

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{r}_{m}=\boldsymbol{P}_{m}-\boldsymbol{P}_{ref}, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{P}$ and $\unicode[STIX]{x1D6FF}(r)$ are the position vector and Dirac delta function. The averaging operation is performed over all particles at a given $y$ location at different instances of time. The superscript ‘ $R$ ’ indicates that the quantity is measured for a random distribution of points seeded in the same domain as our DNS dataset, with the following two additional constraints: (i) the minimum distance between two points is always larger than a particle diameter (non-overlapping condition); (ii) the distribution of the points in the $y$ direction matches that of the DNS dataset within 0.5 % error. The latter condition ensures that $g(r)$ is not biased by the heterogeneity in the $y$ direction and that any deviation of $g(r)$ from unity is a signature of a local aggregation/segregation of particles.

Figure 16 shows the RDF in the near-wall and channel centre regions. The value of $g(r)$ is expectedly high at a separation approximately equal to one particle diameter in the Newtonian and viscoelastic cases, because the collision force in that region stabilizes particle pairs. The peak is, however, more than twofold higher in the viscoelastic cases, which signals an increase in the formation of particle pairs. At the channel centre, the propensity of elasticity is slightly reduced. The RDF rapidly tends to the unity at both wall-normal heights, be that for the Newtonian or viscoelastic flows, which indicates absence of any long-range particle structure.

Figure 16. Radial distribution function (4.3) with reference position in the region (a) $0.8<y<1.2$ and (b) $y<0.4$ . Line types as listed in table 1: (

) W0P5; (

) W15P5; (

) W25P5.

In laminar viscoelastic flows, particle pairs initially aligned in a cross-stream configuration experience a strong attraction due to the normal stresses (Goyal & Derksen Reference Goyal and Derksen2012; Zenit & Feng Reference Zenit and Feng2018). If the initial configuration is alignment in the streamwise direction, repulsion or attraction of the particle pairs is determined by their initial distance: a distance larger than a critical value, which itself is inversely proportional to viscoelasticity, results in repulsion due to the negative wake of the leading particle. Therefore, provided elasticity is sufficiently high, repulsion takes place regardless of initial separation (D’Avino, Hulsen & Maffettone Reference D’Avino, Hulsen and Maffettone2013). In suspensions, the co-existence of these effects may give rise to more complex structures such as particle chains (Zenit & Feng Reference Zenit and Feng2018). The microstructure of particles, and specifically preferential directionality in particle-pair alignment, is examined in our simulations.

We define the particle-pair distribution function, $q(r,\unicode[STIX]{x1D713},\unicode[STIX]{x1D703})$ , which depends on the pair separation $r$ , the polar angle $\unicode[STIX]{x1D713}$ relative to the positive $z$ axis and the azimuthal angle $\unicode[STIX]{x1D703}$ measured counterclockwise from the positive $x$ axis,

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle q(r,\unicode[STIX]{x1D713},\unicode[STIX]{x1D703})=\langle \text{d}N(r,\unicode[STIX]{x1D713},\unicode[STIX]{x1D703})\rangle /\langle \,\text{d}N(r,\unicode[STIX]{x1D713},\unicode[STIX]{x1D703})^{R}\rangle , & \displaystyle\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\displaystyle & \displaystyle \text{d}N(r,\unicode[STIX]{x1D703})=\mathop{\sum }_{m=1}^{N_{p}}\unicode[STIX]{x1D6FF}(r-|\boldsymbol{r}_{m}|)\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}_{m}),\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D713}-\unicode[STIX]{x1D713}_{m}), & \displaystyle\end{eqnarray}$$

(4.8) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{r}_{m}=\boldsymbol{P}_{ref}-\boldsymbol{P}_{m}, & \displaystyle\end{eqnarray}$$

(4.9a,b ) $$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D713}_{m}=\cos ^{-1}\left({\displaystyle \frac{\boldsymbol{r}_{m}\boldsymbol{\cdot }\boldsymbol{e}_{z}}{|\boldsymbol{r}_{m}|}}\right);\quad \unicode[STIX]{x1D719}_{m}=\arctan \left({\displaystyle \frac{\boldsymbol{r}_{m}\boldsymbol{\cdot }\boldsymbol{e}_{y}}{\boldsymbol{r}_{m}\boldsymbol{\cdot }\boldsymbol{e}_{x}}}\right), & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{e}_{\{x,y,z\}}$ are the unit vectors along the Cartesian coordinates. Similar to the RDF, averaging is performed over the ensemble of reference particles located at a given $y$ location at different instances of time, and superscript ‘ $R$ ’ indicates a random distribution conditioned in a similar manner as in the computation of the RDF. Figure 17(a–c) compares $q$ for the different cases at the channel centre. Due to symmetry at this location, the azimuthal angle is only plotted in the range $0<\unicode[STIX]{x1D703}<\unicode[STIX]{x03C0}$ . The narrow stripe with high values of $q$ immediately surrounding the particle, near $r=1$ , highlights the strong effect of particle collisions in stabilizing particle pairs. This effect is intensified in the viscoelastic cases, in line with the RDF results. It should be noted that, due to the absence of shear and slip velocity, we do not observe any significant preferential alignment of particle pairs near the channel centre in the Newtonian or viscoelastic cases.

Figure 17. Particle-pair distribution function $q(r,\unicode[STIX]{x1D713},\unicode[STIX]{x1D703})$ , averaged in the range $-\unicode[STIX]{x03C0}/8<\unicode[STIX]{x1D713}-\unicode[STIX]{x03C0}/2<\unicode[STIX]{x03C0}/8$ . Reference particle is located in the region (a,b,c) $0.8<y<1.2$ and (d,e,f) $y<0.4$ .

The particle-pair distribution function near the wall is reported in figure 17(d–f). In all three panels, two narrow stripes with high values of $q$ are observed immediately around the particle in the two quadrants $\unicode[STIX]{x03C0}/2<\unicode[STIX]{x1D703}<\unicode[STIX]{x03C0}$ and $3\unicode[STIX]{x03C0}/2<\unicode[STIX]{x1D703}<2\unicode[STIX]{x03C0}$ . They arise, respectively, due to the deficit and excess particle phase velocity relative to the local fluid in the presence of mean shear – an effect that has been reported for suspensions with a Newtonian carrier fluid (Morris Reference Morris2009). A relative depletion of particles occurs in the straining quadrants $0<\unicode[STIX]{x1D703}<\unicode[STIX]{x03C0}/2$ and $\unicode[STIX]{x03C0}<\unicode[STIX]{x1D703}<3\unicode[STIX]{x03C0}/2$ , evidenced by a region of lower probability of finding a particle pair. Viscoelasticity results in two unique features: (i) the particle-free regions are markedly extended in the streamwise direction and (ii) $q$ increases at the outer edges of these two regions which indicates that an inclined alignment with respect to the streamwise coordinate is promoted. In order to explain these effects, we evaluated the particle-conditioned mean streamwise velocity $\langle u_{f}\rangle _{pc}$ relative to $\langle u_{f}\rangle$ . A positive and a negative isosurface of $\langle u_{f}\rangle _{pc}-\langle u_{f}\rangle$ are shown in figure 18 for case W15P5. The isosurfaces extend downstream and upstream of the reference particle, respectively, and their effect is to expel leading and trailing particles, consistent with the streamwise extended regions of low $q$ in figure 17(b).

4.3 Flow visualization

Figure 18. Particle-conditioned mean streamwise fluid velocity $\langle u_{f}\rangle _{pc}$ relative to $\langle u_{f}\rangle$ , for case W15P5. The conditional average is performed for particles located at $y=0.2$ .

Figure 19. Particle-conditioned mean flow field relative to the unconditional fluid velocity, $\langle u_{f}\rangle _{pc}-\langle u_{f}\rangle$ . The conditional average is performed for particles located at $y=0.1$ and which experience $u_{slip}=0.05$ . Contours show the streamwise component, and the vectors correspond to the in-plane velocities. (a) W0P5; (b) W15P5.

In order to obtain a precise interpretation of the flow in the vicinity of the particle, we evaluate the average flow field conditioned on the particle position and also the particle slip velocity. Figure 19 shows results for particles whose centre is at $y=0.1$ and whose slip velocity is $u_{slip}=0.05$ , from the Newtonian and W15P5 cases. The conditionally averaged in-plane velocity vectors and contours of their streamwise components are reported after subtracting the mean fluid velocity $\langle u_{f}\rangle$ . The flow fields are dramatically different in the wake regions: a classic backwash behind the particle is observed in case W0P5; in contrast, a negative zone in the trail of the particle in case W15P5 signals the presence of a reverse wake which is enfolded by a conical entraining zone of accelerated fluid. It should be noted that the negative zone in the contours of $\langle u_{f}\rangle _{pc}-\langle u_{f}\rangle$ is a signature of decelerated streamwise velocity only, while $u_{f}$ is almost never negative in an absolute sense. Similar wake patterns have been previously observed in experiments (Kemiha et al. Reference Kemiha, Frank, Poncin and Li2006) and simulations (Frank & Li Reference Frank and Li2006; Goyal & Derksen Reference Goyal and Derksen2012) of an isolated rising bubble or falling solid sphere in extensional flows. Compared to those studies, however, the negative wake and the positive cone are spatially truncated due to the presence of neighbouring particles in our cases. According to the 2-D axisymmetric numerical simulations by Frank & Li (Reference Frank and Li2006), where the Reynolds number based on the particle diameter is more than unity, the opening angle of the conical zone is independent of the Reynolds number, while it is weakly and inversely correlated with the Weissenberg number. The wake characteristics in figure 19(b) explain the peculiar particle microstructure observed in the viscoelastic cases (cf. figure 17): while the negative wake repels the trailing particles, the entrainment zone stabilizes particle pairs with an inclined alignment with respect to the streamwise direction.

Figure 20. Instantaneous contours of the streamwise component of the polymer stress tensor $\unicode[STIX]{x1D61B}_{xx}$ at $y=0.1$ . (a): W15; (b) W15P5. Lines are isocontours of (solid lines

) positive and (dash-dotted lines

) negative streamwise fluid velocity perturbations, $u_{f}-\langle u_{f}\rangle$ . Pink circles mark the positions of particles that are cut by the $xz$ plane.

Figure 20 shows instantaneous plan views of the streamwise component of the polymer stress tensor at $y=0.1$ , from cases W15 and W15P5; two isolevels of streamwise velocity perturbation are also marked. In case W15, the high- and low-stress regions have similar shapes and length scale as the isolines of velocity perturbations; the co-location of the stress and velocity perturbations is qualitatively indicative of the connection between turbulent structures and polymer stretch. It should be noted that, for the sake of clarity of the visualization, the isolines are plotted for one positive and one negative value of the velocity perturbation only, and thus do not represent all the turbulent structures. In contrast, there is no conspicuous similarity between region of high/low polymer stress regions and the isolines of velocity perturbations in case W15P5. As a matter of fact, the polymer stress is mostly concentrated in the vicinity of the particles, which highlights the change in the dominant stretch mechanism in particle-laden conditions (cf. § 4.1).

Figure 21. Instantaneous contours of $u_{f}-\langle u_{f}\rangle$ in the (a,d) Newtonian flows, and in the lower Weissenberg viscoelastic cases during (b,e) active and (c,f) hibernating states. (a–c) Single-phase simulations; (d–f) two-phase simulations. Particle positions are visualized in two sub-regions of the domain, near the wall and at channel centre. The horizontal planes are visualized at $y^{+}\approx 8$ in all cases.

Instantaneous visualizations of streamwise flow structures and particle positions are provided in figure 21. The figure contrasts realizations from the single-phase and particle-laden simulations, for the Newtonian fluids and at $Wi=15$ ; for the latter case, the active and hibernating states are shown. These instantaneous fields encapsulate some important conclusions drawn in the previous sections: (i) particles weaken the streamwise streaks in Newtonian and viscoelastic flows; (ii) the modulation of turbulent structures is relatively minimal in the Newtonian case and pronounced in the viscoelastic flows; (iii) particles migrate towards the channel centre during hibernating states, and back towards the wall during active intervals. Additionally, the wake-induced fluid fluctuations are discernible in the near-wall region of the viscoelastic condition, particularly during the active turbulence period.