Modulation of turbulence by finite-size particles in statistically steady-state homogeneous shear turbulence

Abstract We perform interface-resolved simulations to study the modulation of statistically steady-state homogeneous shear turbulence by neutrally buoyant finite-size particles. We consider two shapes, spheres and oblates, and various solid volume fractions, up to 20%. The results show that a statistically steady state is not exclusive to single-phase homogeneous shear turbulence as the production and dissipation rates of the turbulent kinetic energy are also statistically in balance in particle-laden cases. The turbulent kinetic energy shows a non-monotonic behaviour with increasing solid volume fraction: increasing turbulence attenuation up to a certain concentration of solid particles and then enhancement of the turbulent kinetic energy at higher concentrations. This behaviour is observed at lower volume fractions for oblate particles than for spheres. The attenuation of the turbulence activity at lower volume fractions is explained through the enhancement of the dissipation rate close to the surface of particles. At higher volume fractions, however, particle pair interactions induce regions of high Reynolds shear stress, resulting in the enhancement of the turbulence activity. We show that the oblate particles of the considered size have larger rotational rates than spheres with no preferential orientation. This is in contrast to previous studies in wall-bounded flows where preferential orientation close to the wall and reduced rotation rates result in turbulence attenuation and thus drag reduction. Our results shed some light on the effect of rigid particles, smaller than the near-wall turbulent structures but still comparable to the viscous length scale, on the dynamics of the equilibrium logarithmic layer in wall-bounded flows.


Introduction
The occurrence of dispersed two-phase flows is widespread in many environmental and industrial applications. Sediment transport in estuaries (Mehta 2013), blood flow in the human body, pyroclastic flows from volcanoes and pulp fibers in paper making (Lundell, Söderberg & Alfredsson 2011) are among the examples of flows that deserve further particle kinetic energy increase with the particle aspect ratio and are substantially larger than those for spherical particles.

Wall-bounded turbulent flows
In wall-bounded flows, the turbulence modulation by finite-size particles is more complex and less predictive as the confinement effect of the wall creates additional implications. The first simulations of finite-size particles in a turbulent channel flow were performed by Pan & Banerjee (1996). Those authors revealed that turbulent fluctuations and stresses increase in the presence of the solid phase. Matas, Morris & Guazzelli (2003), Loisel et al. (2013) and Yu et al. (2013) reported a decrease of the critical Reynolds number for transition to turbulence in the semi-dilute regime with neutrally buoyant spherical particles. Picano, Breugem & Brandt (2015) investigated dense suspensions in turbulent channel flow up to a volume fraction of 20 %. Their study revealed that the overall drag increase is due to the enhancement of the turbulence activity up to a certain volume fraction (φ ≤ 10 %) and to significant particle-induced stresses at higher concentrations. Costa et al. (2016) explained that the turbulent drag of sphere suspensions is always higher than that predicted by only accounting for the effective suspension viscosity for particle sizes of the order of 20 viscous units. They attributed this increase to the formation of a particle-wall layer, a layer of spheres forming near the wall in turbulent suspensions. Based on the thickness of the particle-wall layer, they proposed a relation able to predict the friction Reynolds number as a function of the bulk Reynolds number (Costa et al. , 2018. Ardekani, Rosti & Brandt (2019) showed in a numerical experiment that removing the particle-wall layer results in turbulence attenuation with respect to single-phase flow, while the presence of this layer contributes to larger velocity fluctuations close to the wall for lower particle volume fractions.
Studies of finite-size non-spherical particles in a turbulent channel flow are more scarce in the literature (Do-Quang et al. 2014;Ardekani et al. 2017;Eshghinejadfard, Hosseini & Thévenin 2017;Eshghinejadfard, Zhao & Thévenin 2018;. Those studies showed that prolate and oblate spheroids preferentially align with the wall in its vicinity, experiencing considerably smaller rotational rates with respect to spheres. This dampens the wall-normal velocity fluctuations and thus attenuates the turbulence. For prolate particles, this effect is less pronounced since their larger angular velocities create additional counteracting stresses ).

Homogeneous shear turbulence
In the limit of very large Reynolds numbers, particle-resolved simulations of wall-bounded multiphase flows are no longer feasible. However, investigating particle suspension in homogeneous shear turbulence (HST) (Tavoularis & Corrsin 1981a,b;Pumir 1996;Mashayek 1998;Sekimoto, Dong & Jiménez 2016;Rosti et al. 2019) can provide us with clues on multiphase flow behaviour in this regime. In HST, the flow remains statistically homogeneous in all spatial directions while the turbulence is sustained through a natural energy production mechanism, owing to the presence of a mean velocity gradient. Given a linear mean velocity (constant shear rate), HST simulations reproduce the dynamics of the equilibrium logarithmic layer in wall-bounded turbulence (Sekimoto et al. 2016). Even though the ideal HST is self-similar with unbounded energy growth (Sukheswalla, Vaithianathan & Collins 2013), considering a finite computational domain bounds the large eddies and affects the flow similarly to the confinement effect enforced by a wall. Indeed, Pumir (1996) showed that a statistically stationary state can be reached over long periods of time, denoted SS-HST. Most of the simulations of Pumir (1996) were performed in a cubic box, while Sekimoto et al. (2016) revealed that an appropriate box aspect ratio is essential to reproduce one-and two-point statistics that agree with those in the logarithmic layers in turbulent channel flows. Dong et al. (2017) further compared the dynamics between SS-HST and turbulent channel flow.
Most recently, suspensions of finite-size solid spherical particles in transient HST have been investigated in the dilute regime (Tanaka & Teramoto 2015;Tanaka 2017). It has been shown in those numerical studies that the turbulence kinetic energy is attenuated in the presence of spherical particles. Motivated by this, we study here the modulation of SS-HST at high particle concentration and by particles of different shapes, for the first time. We use interface-resolved simulations to investigate neutrally buoyant spherical and oblate particles (with aspect ratio 1/3) in SS-HST up to volume fractions of 20 % and 10 %, respectively. The particle size is ≈ 20 Kolmogorov length scales, of the order of the Taylor microscale. The results are compared with those for wall-bounded flows . The paper is organised as follows. The governing equations, numerical method and the flow geometry are introduced in § 2, followed by the results of the numerical simulations in § 3. The main conclusions and final remarks are presented in § 4.

Governing equations
The evolution of the fluid phase is described by the incompressible Navier-Stokes equations for a Newtonian fluid: where u is the fluid velocity vector, ρ and ν the fluid density and kinematic viscosity and p the pressure. The last term on the right-hand side of (2.2) accounts for the presence of the particles through immersed boundary method (IBM) forcing, active close to their surface (see § 2.2.3). By decomposing the velocity field into u = U + u , where U = (Sy, 0, 0) is the mean shear flow, with S the shear rate, and u = (u , v , w ) the velocity fluctuations in the streamwise, normal and spanwise directions, the equations for the fluctuations are Here, the third term on the left-hand side of (2.4) denotes the advection of the mean shear by the velocity fluctuations and is responsible for the turbulent energy production in HST. The dynamics of the rigid particles is governed by Newton-Euler equations for the conservation of linear and angular momentum: where u p and ω p are the particle linear and angular velocity vectors, m p and I p denote the particle mass and moment of inertia, r is the position vector with respect to the particle centre and n is the outward-pointing normal to the particle surface ∂Ω p . The fluid stress tensor is given by τ = −pI + νρ(∇u + ∇u T ). Finally, F col and T col denote the force and torque resulting from short-range particle-particle interactions, such as lubrication and collisions.
The sets of equations governing each phase are coupled through the no-slip and no-penetration condition at the particle surface, i.e. (2.7) 2.2. Numerical method The governing equations are solved numerically, using the method proposed by Gerz, Schumann & Elghobashi (1989), which employs the shear-periodic boundary condition. This was later modified by Tanaka (2017) to also handle the dispersed phase, using IBM, and is explained briefly here.

Fluid phase
To solve the equation for the fluid velocity fluctuations (equation (2.4)), we treat the advection by the mean shear flow Sy(∂u /∂ x) separately, by means of discrete Fourier interpolation. All other terms are evolved using a three-step Runge-Kutta method, except the pressure gradient term, for which the Crank-Nicolson scheme is used. The equations are discretised in space on a uniform, staggered Cartesian grid with the finite-volume method in which spatial derivatives are estimated with the central-differencing scheme; see Breugem (2012) for details of convergence proof, order of accuracy and validation of the flow solver. In particular, the first prediction velocity u * is obtained as where RHS ≡ −Sv x − ∇ · (u u ) + ν∇ 2 u , q = 1, 2, 3 denotes Runge-Kutta substeps and α q = (8/15, 5/12, 3/4) and β q = (0, −17/60, −5/12) are the Runge-Kutta coefficients. Note that the term RHS q−2 is advected by the mean shear. This modification, introduced by Tanaka (2017), enhances the stability and accuracy of the scheme and is performed by discrete Fourier interpolation: (2.9) In the next step, the first prediction velocity and pressure are advected by the mean shear flow, again using discrete Fourier interpolation: The second prediction velocity u * * is then obtained as where the operator ∇ q−1/2 evaluates the pressure gradient at the substep q − 1/2 and is defined as in Tanaka (2017) as (2.13) Given the source term f , which accounts for the interaction between the dispersed phase and the carrier fluid using IBM (Breugem 2012;§ 2.2.3), the third prediction velocity u * * * is obtained as (2.14) Finally, after solving the Poisson equation for the correction pressurep, the velocity and the pressure are corrected as

Shear-periodic boundary condition
The advection by the mean shear flow makes it impossible to seek for periodic solutions in the normal direction. In this direction, the so-called shear-periodic boundary condition holds (Gerz et al. 1989), which, for an arbitrary quantity h, reads as where L x and L y are the size of the computational domain in the streamwise and normal directions.

Dispersed phase and IBM forcing
The fluid and solid phases interact through the direct forcing IBM; see, for example, Kajishima & Takiguchi (2002) where the authors used a volume of fluid approach to account for the presence of the particles and Uhlmann (2005) and Breugem (2012) where the surface of the particles is tracked via a set of Lagrangian points. The method has been validated and used extensively for unbounded (see e.g. Ardekani et al. 2016;Fornari, Ardekani & Brandt 2018) and wall-bounded (see e.g. Lashgari et al. 2017; flows laden with finite-sized spheroidal particles. In this approach, the particles are discretised by a set of Lagrangian points, uniformly distributed along their surface. The procedure to integrate the Newton-Euler equations is similar to that for the fluid phase. First, we decompose the particle velocity u p into the mean fluid velocity at the centroid of the particleū(x p ) and a fluctuation velocity u p : is the position vector of the particle centroid. Next, the particles and the Lagrangian grid points attached to them are advected by the mean shear flow: where X l denotes the position vector of the lth Lagrangian grid point. The next steps to calculate the IBM force are similar to Breugem (2012): first, interpolation of the first fluid prediction velocity u * from the Eulerian to the Lagrangian grid; then, calculation of the IBM force, using the slip velocity between the fluid velocity fluctuations and that of the particles at the location of each Lagrangian grid point; finally, spreading the resulting IBM force from the Lagrangian to the Eulerian grid. The interpolation and spreading operations are done using the regularised Dirac delta function of Roma, Peskin & Berger (1999), which acts over three grid points in all coordinate directions. When the gap between two particles is smaller than the grid spacing, the IBM fails to resolve the hydrodynamic interactions. Therefore, we use a lubrication correction model, based on the asymptotic analytical expression for the normal lubrication force between spheres of different sizes (Jeffrey 1982), for subgrid hydrodynamic interactions. Spheroidal particles are approximated as spheres with radius equal to the local radius of curvature of the spheroidal particle (Ardekani et al. 2016). This lubrication force is kept constant below a second threshold for the distance between particles, to account for the surface roughness of the particles. When particles are in collision, the lubrication force is turned off, and a collision force based on the soft sphere model is activated. The collision model works based on a mass-spring-damper system in the directions normal and tangential to the contact line between the overlapping particles, and calculates the collision force based on the particle relative velocity and overlap. Details of the collision model are provided in Costa et al. (2015), later adapted by Ardekani et al. (2017) to model the close-range interactions between spheroidal particles.

Computational set-up
In this study, we simulate a suspension of rigid neutrally buoyant spherical/oblate particles, subjected to a uniform mean shear flow. The spherical and oblate particles have the same volume V and their characteristic length is denoted by D eq = (6 V/π) 1/3 , i.e. the diameter of a sphere with the same volume. The dimensions of the computational box are L x × L y × L z = 40D eq × 20D eq × 19D eq , with N x × N y × N z = 1280 × 640 × 608 Eulerian grid points in the streamwise, normal and spanwise directions. The aspect ratios of the computational box are chosen as L x /L z ≈ 2.1 and L y /L z ≈ 1.05; the steady-state simulations of HST are considered as minimal in the spanwise direction, i.e. containing on average only a few large-scale structures along the spanwise direction (Rogers & Moin 1987;Sekimoto et al. 2016). The flow is periodic in the streamwise and spanwise directions, with the shear-periodic boundary condition imposed at the top and bottom boundaries.
The non-dimensional numbers that characterise the fluid phase are the Taylor microscale Reynolds number Re λ and the shear-rate parameter S * , which manifests the ratio of the dk is the integral length scale, with E(k) the energy spectrum at each wavenumber k; Sk = τ p /τ η denotes the Stokes number, with τ p = (ρ f + 2ρ p )D 2 eq /(36ρ f ν) the particle response time and τ η = (ν/ ) 1/2 the Kolmogorov time scale. The reported quantities are statistically averaged when the flow has reached the stationary state.
'eddy turnover' time 2K/3 to the time scale of the mean deformation 1/S (Lee, Kim & Moin 1990): where the Taylor microscale is defined as λ ≡ 10νK/ , K = 1/2 u 2 1/2 denotes the turbulent kinetic energy per unit mass, = ν ω 2 is the energy dissipation rate, ω = ∇ × u is the fluctuating part of the vorticity vector and · indicates statistical average. To obtain the initial field, we start a single-phase flow case from a homogeneous isotropic turbulence velocity field with a prescribed energy spectrum and a random phase at the non-dimensional time St = 0 (Tanaka 2017). The initial microscale Reynolds number Re λ (St = 0) and shear-rate parameter S * (St = 0) are set to 113 and 2.9 and change as the turbulent field develops. When the statistically stationary state (SS-HST) is reached, the velocity field is saved and used for the particle-laden cases (see § 2.4). The ratio between the grid spacing and the Kolmogorov length scale -defined as η = (ν 3 / ) 1/4 -is equal to 0.16 at the beginning of the simulation and reaches to 0.78 at the steady state, which guarantees that all scales are well resolved. The particles are neutrally buoyant, with a relative size of D eq /η ≈ 20 at St = 0. We consider particles with aspect ratio (ratio of polar over equatorial radius) AR = 1 (spheres) and AR = 1/3 (oblates). The surface of the particles is tracked, using 3219 Lagrangian grid points for spheres and 3720 in the case of oblates. The particles are introduced randomly into the computational domain, with initial velocity equal to the local mean flow velocity. The physical and computational parameters of the main simulation cases are summarised in table 1. The cases pertaining to spherical particles are denoted as Spx, whereas Obx is used for oblate particles. The number x defines the solid volume fraction. We perform simulations at four different volume fractions φ = 1 %, 5 %, 10 % and 20 % for the spheres and at two volume fractions φ = 5 % and 10 % for the oblates, with a single-phase case for direct comparison. Due to the higher computational cost of the interface-resolved simulations of oblate particles, we do not simulate a case with φ = 20 % and φ = 1 %; in the latter case, however, it has been shown in previous studies that shape effects, like excluded volume effect, are not significant for volume fractions φ < 5 % (see e.g. Fornari et al. 2016a;Ardekani et al. 2017). Hence, the statistics would be relatively close to those pertaining to the case of spheres at φ = 1 % and to the single-phase flow.

Validation and characteristics of statistically stationary state
To test the numerical code for the specific case of HST, the single-phase HST case Si is assessed against the results of Pumir (1996). The initial homogeneous isotropic turbulence field, described in the previous section, is subjected to the shear rate S at St = 0 and statistics are collected afterwards. Figure 1(a) shows the time history of the box-averaged turbulent kinetic energy (black line) and enstrophy Ω = ω i ω i (green line), normalised by their time-averaged values. Two distinct states are distinguishable through the time evolution of the flow. First, until St ≈ 30, the turbulent kinetic energy grows faster than enstrophy, which indicates excess production over dissipation. After the initial transient state, the flow reaches a statistically stationary state when the production and the dissipation rates of the turbulent energy are almost in balance. This is characterised by a sequence of spikes of the turbulent energy, followed by spikes of enstrophy with a delay of approximately 5St, which is evident in our results and in those of Pumir (1996). Figure 1(b) shows the normalised probability density function of the velocity components. The figure shows that the distributions of the normal and spanwise components have more extended tails than that of the streamwise component. This is in agreement with the results of Pumir (1996), suggesting that the strong anistropy may arise from the anisotropic forcing of HST. In particular, we simulated Run No. 2 in Pumir (1996) in a cubic box with 256 grid points in each direction. The velocity anisotropy tensor components b ij = u i u j /u k u k − δ ij /3 in our simulation are b 11 = 0.231, b 22 = 0.129 and b 12 = 0.147, and have a maximum difference below 5 % from those reported in the cited reference.
When the single-phase case Si reaches the statistically stationary state, in which the production and the dissipation rates of the turbulent kinetic energy are statistically in balance P ≈ (see the black line in figure 1c), the flow field is used to initialise the multiphase cases. At this point, the turbulence flow parameters evolve from the initial values to Re λ = 103 and S * = 2.4. The rigid particles are introduced randomly with volume fraction ranging from 1 % to 20 %. The time history of the ratio between the production and the dissipation rates of the turbulent kinetic energy for case Sp5 (the red line in figure 1c) confirms that the stationary state is not exclusive to the single-phase HST; after the introduction of the dispersed phase, the flow goes through a relatively short transient state and reaches a second steady state, in which the production and the dissipation rates of the turbulent energy are statistically in balance. In their recent study, Rosti et al. (2019) documented the existence of a steady state in the presence of deformable droplets for a similar initial field and geometry.

Results
We start by showing snapshots of the flow field and particles. In figure 2 we show the two cases Sp10 and Ob10, with the instantaneous streamwise velocity u contours depicted on the vertical and horizontal planes. Only particles with streamwise centre coordinate larger than 0.75L x are shown in the figure for better clarity. We observe that the level of the streamwise fluctuations is higher for the case of oblates, suggesting higher turbulent activity and hence Re λ . Also, more distinct patches of low-speed velocity fluctuations are present for the case Ob10 compared to Sp10.

Turbulence modulation
In this first section, we quantify the turbulence modulation by the presence of the particles by discussing the statistically averaged flow parameters for all the cases. The averaged Eulerian fluid statistics reported here correspond to mean intrinsic averages. The intrinsic average of a quantity ξ is computed as where Ψ ijk,t is the fluid volume fraction at the grid cell ijk and instant t. Figure 3(a) depicts the Taylor microscale Reynolds number Re λ = (2K/3) 1/2 λ/ν as a function of the particle volume fraction φ. The results show that increasing the volume  fraction of spheres up to φ = 10 % decreases Re λ with respect to the single-phase flow. However, at φ = 20 % a considerable jump (≈ 17 %) in Re λ is observed, indicating a non-monotonic effect of the volume fraction. The same behaviour can also be observed for the oblate particles, with the increase in Re λ happening at a lower volume fraction, i.e. Re λ is lower for Ob5 compared to Sp5, but attains a higher value than single phase (Si) and Sp10 at a concentration of 10 % (Ob10). Figure 3(b,c) displays the change of the turbulent kinetic energy K = u i u i /2 and of the Taylor microscale λ = 10νK/ , the two parameters defining Re λ , as a function of the solid volume fraction. The data reveal that the variation of Re λ can be mainly attributed  to the modulation of the turbulent energy by the particles, as K shows the same trend as Re λ (cf. figure 3a). The Taylor microscale λ, conversely, decreases monotonically with increasing volume fraction for all the cases; this variation, however, is only ∼3 %. Note that a reduction of the turbulent energy in a dilute suspension of rigid spheres (φ = 0.5 %) was also observed in the study by Tanaka & Teramoto (2015) in transient HST. Those authors show that the increase of the viscous dissipation is responsible for the decrease of the turbulent energy. Finally, figure 3(d) shows the variation of the shear-rate parameter S * = 2SK/3 as a function of the volume fraction for all cases under consideration. The trend is similar to that of the Taylor microscale that can be seen in figure 3(b). The shear-rate parameter is defined using the dissipation length l d = (2K) 3/2 / , the length scale associated with energy-containing eddies (Lee et al. 1990). Therefore, we infer from figure 3(c,d) that the length and time scales of the energy-containing eddies are decreasing with increasing solid volume fraction.
The averaged spectra of the turbulent kinetic energy are reported in figure 4. As expected, we observe a range at intermediate wavenumbers where E ∝ k −5/3 , similar to that found in the logarithmic layer of wall-bounded flows (Sekimoto et al. 2016). Nonetheless, the presence of the solid particles has a significant effect on the energy spectrum; they increase the energy level at large wavenumbers (small scales), while they slightly reduce the energy of the low-wavenumber (large) structures. This effect is amplified by increasing the volume fraction of the solid phase, while the shape effects (spherical versus oblate particles) are observed to be rather negligible.   Pumir 1996), i.e. the streamwise component has a higher root-mean-square (r.m.s.) value than the normal and spanwise ones. Interestingly, the difference between the spanwise and the normal components of the fluid velocity increases in the presence of the solid particles (cf. figure 1b).
To understand the role of particle volume fraction, we first compare the cases laden with spherical particles. Up to φ = 10 %, the normal velocity has a similar distribution to the single-phase HST, whereas the streamwise and spanwise components have stronger tails. Increasing the volume fraction to 20 % leads to an increase in the probability of strong fluctuations in all velocity components, which consequently results in a higher magnitude of the turbulent kinetic energy and Re λ . On the other hand, the main difference for oblate particles is in the distribution of the normal velocity for φ = 10 %, where the probability of extreme fluctuations has increased significantly. The higher value of the p.d.f. in the range |v |/(SD eq ) > 5 again results in the higher values of the turbulent energy and Re λ documented above for oblates when increasing the particle volume fraction.
This can also be observed in figure 5(d), where the r.m.s. velocity fluctuations are averaged and depicted for all the cases under investigation. Interestingly, the effect of the volume fraction is more pronounced on u than on the normal and spanwise components of the velocity fluctuations for spheres up to φ = 10 % and for oblates at φ = 5 %, i.e. the reduction of the turbulent energy for these cases is associated with the reduction of u rather than of the two cross-stream components. However, for cases Sp20 and Ob10 all the components experience a significant increase with respect to the single-phase case, owing to the long tails of the p.d.f.s. The relation between these results and the particle dynamics will be explained in § 3.2.
Thus far, we have discussed the modification of the fluctuating velocities and related parameters to quantify the turbulence modulation by the presence of the particles. To have a better picture of the events responsible for the modification of the turbulence, we present the weighted Reynolds shear stress contours, computed by multiplying the absolute value of the Reynolds shear stress by the joint probability density of its occurrence in the u -v plane (Zhou et al. 1999). This method considers the contribution of the different events to the total stress by locating them on the four quadrants of the u -v plane, denoted Q 1-4 . The events on the Q 2 (u < 0, v > 0) and Q 4 (u > 0, v < 0) quadrants result in production of turbulence, whereas Q 1 (u > 0, v > 0) and Q 3 (u < 0, v < 0) are responsible for attenuation.
The contours of the weighted Reynolds shear stress, statistically averaged in the fluid phase for all the cases, are displayed in figure 6. For all cases, the diagrams are symmetric with respect to the origin, which fulfils the symmetry of the HST under the transformation (x, y, z) → (−x, −y, z), and, expectedly, confirms that the inhomogeneity observed in the quadrant analysis of the wall-bounded flows (see e.g. Wallace 2016) vanishes in unbounded ones. Figure 6(a-d) shows the results for the cases laden with spheres. For cases Sp1 (figure 6a), Sp5 (figure 6b) and Sp10 (figure 6c), the contribution of events Q 2 and Q 4 to the production and the quenching events Q 1 and Q 3 have not changed considerably (see iso-line of 0.6 in figure 6a-c). The top contributor events of the Q 2 and Q 4 quadrants have higher magnitude of the streamwise component u compared to the normal one v . Conversely, in case Sp20 (see figure 6d), the production associated with Q 2 and Q 4 events becomes more predominant as the magnitude of the velocity fluctuations for the top contributor events increases for both the streamwise and the normal components; at the same time, the attenuation induced by Q 1 and Q 3 events is much weaker than for the cases with lower volume fraction of spherical particles. Note also that the importance of the streamwise and normal components of the velocity vector in the production and damping mechanisms is more balanced in case Sp20.
To conclude this analysis, figure 6(e, f ) displays the contours of the weighted Reynolds shear stress for cases Ob5 and Ob10. Here, we see that the area of the contours contributing to both production and attenuation increases slightly when increasing the volume fraction, indicating the importance of the contribution from rare but highly energetic events. When comparing the flow laden with oblate particles with the case of spheres at the same volume fraction, we note that the contours are more skewed towards higher values of the velocity, which confirms the increased role of the higher-intensity events in the case of oblate particles.

Turbulent kinetic energy budget
To quantify the effect of the dispersed phase on the modulation of the turbulence, we look at the turbulent kinetic energy budget for the fluid phase. The governing equation for the evolution of the turbulent kinetic energy K reads

4)
Sp10 Ob5 Ob10 FIGURE 7. Contribution of turbulent production P, dissipation rate and the interphase interaction term I to the turbulent kinetic energy budget for the different cases under consideration. Each term is normalised by the product of the shear rate S and the averaged turbulent kinetic energy of the single-phase case Si.
where I denotes the energy transfer through the interphase interaction (Tanaka & Teramoto 2015), which can act as a source or a sink of turbulent energy (Ferrante & Elghobashi 2003). At the steady state, the rate of change of K is obviously zero and the remaining terms are in balance. For a more clear comparison, we normalise each term on the right-hand side of (3.2) by the product of the shear rate S and the averaged turbulent kinetic energy of the singe-phase case Si: (3.8) The relative contribution of each term σ 1-3 to the turbulent kinetic energy budget for all the cases under consideration is displayed in figure 7. The data reveal that the production and the dissipation rates are almost in perfect balance for all cases, and that the contribution of the interphase interaction term is less than 3.5 % of the total. Note that despite the very small contribution of σ 3 , its presence is necessary to have a correct energy balance. Also, the presence of the solid phase affects the production and dissipation rates of the turbulent kinetic energy, i.e. the variations of σ 1 and σ 2 due to the presence of particles.
Comparing to the single-phase case Si, the production and the dissipation rates decrease and the interphase interaction increases monotonically with the sphere volume fraction for φ ≤ 10 %, whereas at φ ≤ 20 % the production and the dissipation rates are greater than in the single-phase flow and the interphase interaction contributes to the production of the kinetic energy, instead of being a sink of energy as at lower volume fractions.
In the case of the oblate particles, the same trend is observed at lower volume fractions, i.e. case Ob10 has greater production and dissipation rates than cases Ob5 and Si, and the contribution of the interphase interaction to the dissipation of the kinetic energy is lower.

Particle dynamics
In this section, we investigate the link between the dynamics of the particles and the modulation of the turbulence by examining the Lagrangian statistics of the solid particles. Figure 8(a) displays the normalised translational kinetic energy of the solid particles, defined as K p ≡ 0.5(u 2 p + v 2 p + w 2 p ), as a function of the volume fraction φ, whereas figure 8(b) depicts K p normalised by the turbulent kinetic energy of the carrier fluid. The amplitude of the particle velocity fluctuations follows the same trend as that of the carrier fluid, displayed in figure 3(b). The kinetic energy of the spherical particles first decreases when increasing the solid volume fraction until φ = 10 %, and then increases significantly when the volume fraction is increased to φ = 20 %. In the case of oblate particles, the increase of the particle kinetic energy occurs at a lower volume fraction, and in fact we only observe an increase of K p for the two volume fractions considered here, φ = 5 % and 10 %. Scaling the particle kinetic energy by the turbulent kinetic energy of the carrier fluid (cf. figure 8b), we see that rigid particles tend to fluctuate less than the fluid. This is similar to the results of Tanaka & Teramoto (2015), who show that the ratio between the magnitude of the fluctuating particle velocity and that of the fluid for dilute suspensions of spherical particles in transient HST is around 0.85, and to the observations of Picano et al.  (2015) for a channel flow laden with neutrally buoyant spherical particles for the region which lies in the near-wall range 30 < y + < 70 in inner units. By comparing figures 8(a) and 8(b), we infer that increasing the volume fraction of the solid particles, for spheres from 10 % to 20 % and for oblates from 5 % to 10 %, has a larger impact on the increase of the fluid velocity fluctuations than on that of the particle velocity fluctuations. The spanwise component of the particle mean angular velocity, normalised by the rotation associated with the mean shear flow −S/2, is displayed in figure 8(c). First, note that the average spanwise rotation rates are around 8 % more for the spheres than for the oblates. Moreover, for both spheres and oblates, the changes of the averaged spanwise angular velocity by the increase of the volume fraction φ do not show a clear trend. The value obtained by Tanaka & Teramoto (2015) for spherical particles in the dilute regime in transient HST is around 0.9, larger than those obtained here at higher φ. Finally, the magnitude of the particle angular velocity vector |Ω p | = ω 2 x,p + ω 2 y,p + ω 2 z,p is shown in figure 8(d). The rotation rate |Ω p | is around 20 % larger for the oblates than for the spheres and an increase of the volume fraction of the solid particles only slightly increases its value. A comparison between figures 8(c) and 8(d) reveals that the oblate particles have larger angular momentum in the shear-plane directions than spherical particles and also have a clearer tendency to rotate opposite to the mean shear flow vorticity vector (see figure 10c).
Next, we wish to give further insight into the behaviour of the solid phase by examining the velocity probability functions. To this end, we display in figure 9 the p.d.f.s of the three particle velocities (figure 9a-c) together with that of the particle translational kinetic energy K p (figure 9d); see table 2 in the appendix for the values of the different statistical moments. Note that since particles have a considerable size with respect to the flow structures (D eq ≈ 25η) and follow a rigid-body motion, they are able to couple the fluctuations of the fluid velocity in the cross-flow directions, which therefore exhibit a more isotropic behaviour. This is evident in table 2 in the appendix, where the r.m.s. of the particle velocity fluctuations in the cross-flow directions are shown to assume values approximately equal to the average of the r.m.s. of the normal and spanwise fluid velocity fluctuations. Moreover, we observe no significant change in the distributions of the normal and spanwise components of the particle velocity when increasing the volume fraction of the spherical particles up to φ = 10 %, so that the decrease of the translational kinetic energy of the particles K p (figure 9d) can be associated with the decrease of the streamwise component ( figure 9a). At the highest volume fraction considered, φ = 20 % (case Sp20), however, we observe an increase in all the p.d.f.s of the velocity components in the range 5 < |u p,i /(SD eq )| < 10. This isotropic increase in all components can be an indication that at this solid-phase concentration, particle-particle interactions become relevant. The effect of the shape of the particles is evident comparing the distributions for the cases Sp10 and Ob10. The data indicate a higher probability of larger values of the particle translational kinetic energy in the range 100 < K p /(SD eq ) 2 < 200 in the flow laden with oblate particles, which relates to the higher probability of extreme fluctuations of the streamwise velocity component u p .
The p.d.f.s of the components of the particle angular velocity are displayed in figure 10(a-c), whereas that of the magnitude of the angular velocity vector |Ω p | is shown in figure 10(d); the values of the statistical moments can be found in table 2 in the appendix. The data in the figure reveal that an increase of the volume fraction alone only extends the tails of the p.d.f.s slightly, while the shape of the particles has a more noticeable effect, i.e. the cases laden with oblate particles have significantly longer tails compared to the cases with spheres. The spanwise angular velocity of the particles is of particular interest, as this is driven by the imposed shear. The p.d.f.s are clearly skewed towards negative values, which indicates that particles, expectedly, rotate in the direction imposed by the mean shear for most of the instances; the mean values are reported in table 2 in the appendix. One can also note that oblates happen to rotate opposite to the direction of the angular velocity of the mean shear more frequently and at higher rates; this explains why they have lower averaged values of ω z,p than spheres, as shown in figure 8(c).
Next, we investigate possible particle clustering by computing the pair distribution function P(r), which denotes the conditional probability of finding a particle at distance r given one at the origin (see e.g. Kulkarni & Morris 2008): where θ is the polar angle measured from the positive z axis, ψ is the azimuthal angle measured counterclockwise from the positive x axis, n is the average particle number density, t s is the total number of sampling points, ΔV = r 2 sin(θ )Δθ ΔψΔr is the volume of the sampling bin and H(r, θ, ψ) is the histogram of particle pairs. More specifically, the particle-pair histogram H(r, θ, ψ) is progressively built by discretising the pair space in r, θ and ψ and putting each separation vector r into the corresponding bin of size ΔV at each sampling time.
The pair distribution function has been calculated in a spherical shell with diameter of 5D eq around the mass centre of the particles and displayed in figure 11, where the data are shown in the x-y plane (θ = π/2) for the cases Sp10, Sp20 and Ob10. Note that in this figure and subsequent sections, the oblate particles are visualised by a sphere with a diameter equal to their minor axis. All panels show a relative accumulation of particles in the compressional quadrants and a depletion in the extensional quadrants, similar to results from previous experiments and simulations of suspensions in simple shear flow at low Reynolds number (Morris 2009). On increasing the volume fraction, the narrow strips with high values of P become more uniform around the surface of the particle (see figure 11b). This can be another indication that at φ = 20 %, particle-particle interactions significantly affect the particle dynamics. In the case of oblates, the contour levels are below unity around the reference particle and no clustering is observed (see figure 11c). Still, the contours of P have lower values in the extensional quadrants compared to compressional ones. The absence of clustering together with the higher probability of extreme events observed in the p.d.f.s of the linear and angular velocities of the oblates can be indicative of strong collisions between them, as they depart rapidly after contact. x/D eq FIGURE 11. Pair distribution function P(r, ψ) in the plane normal to the spanwise direction (θ = π/2) for the cases (a) Sp10, (b) Sp20 and (c) Ob10.
The dynamics of the particle pairs is not only determined by their relative position, i.e. by the pair distribution function, but also by their relative velocity (Lashgari et al. 2016). Following the study by Sundaram & Collins (1997), we therefore compute the normal relative velocity of the particle pairs as a function of the separation distance r (the distance between the centres). Considering particles p and q, the normal relative velocity of the particle pair is obtained as the inner product of their relative velocity and relative distance: (3.10) The normal relative velocity can be either negative, Δv − n , for approaching particles or positive, Δv + n , when the two particles depart from each other. Figure 12 shows the magnitude of the negative relative velocity |Δv − n | as a function of the separation distance r for some different cases. The relative velocity increases almost monotonically with r as the pairs are more likely to approach with higher velocity when farther away. The approaching velocities are observed to increase slightly with increasing volume fraction, meaning that particles experience stronger pair interactions at higher volume fractions. Note that the abrupt decrease of the approaching velocity for separation distances close to r/D eq = 1 shows the effect of the lubrication force on the dynamics of the spherical particles. In addition, oblate particles can get closer to each other (r/D eq < 1) due to their shape, while maintaining high relative velocities. This creates stronger pair interactions for oblate particles, resulting in higher turbulent kinetic energy when the number of pair interactions increases (high volume fractions).
To better characterise the effect of the shape of the particles, we examine the orientation with respect to the underlying flow of the oblate particles. Figure 13(a) shows the p.d.f.s of the magnitude of the components of the unit vector associated with the particle symmetry axisô = (ô x ,ô y ,ô z ) for the cases Ob5 and Ob10. The data reveal that oblate particles have an almost random distribution of orientations, with a weak tendency to be aligned with the normal and the spanwise directions. Figure 13(b) shows the orientation correlation function OCF(r), which quantifies the relative orientation of close particles: whereÔ p andÔ q denote the orientation of particles p and q at distance r between their centres. The value OCF(r) = 0 indicates a suspension with random particle orientation, while OCF(r) = 1 means particle pairs at distance r perfectly aligned. The results show that oblate particles tend to be more aligned with each other for separation distances in the r/D eq FIGURE 12. Magnitude of relative (approaching) normal velocity |Δv − n | as a function of the separation distance between particle pairs. range 0.5 < r/D eq < 1.5, i.e. as long as they are in contact with each other they move as a small cluster with similar alignment -something that, however, the data in figure 11(c) show to happen rarely. The correlation function vanishes very quickly when increasing the separation distance, showing decorrelation in the alignment of the particles for distances above 1.5D eq . Ardekani et al. (2017) performed the same analysis for oblate particles in a turbulent channel flow. Those authors showed that the correlation distance is significantly smaller in the core of the channel, compared to the near-wall particles experiencing strong shear. The data in figure 13(b) are comparable to the results of the wall-bounded flow for particles in the core region, indicating the negligible effect of the mean shear on the particle orientation.

Flow statistics around particles
In this section, we focus on the conditionally averaged flow statistics around the surface of the particles. The ensemble averaging is carried out in a cubic box of side 3D eq around the  mass centre of the particles with a frequency of 0.5St (one sample every 0.5St) and a total duration of ∼ 150St, after the statistically steady state has been reached. For normalisation, a space averaging has also been carried out in the cubic box of side 3D eq around the particles in addition to the ensemble averaging, and the averaged quantities are denoted by · sur . Figure 14 shows contours of the turbulent kinetic energy K (figure 14a-c), the turbulent production P (figure 14d-f ) and the dissipation rate (figure 14g-i) normalised by their conditional averaged values K sur , P sur and sur , for the cases Sp5, Sp20 and Ob5. Contours are plotted in the x-y plane passing through the centre of the reference particle. Figure 14(a-c) shows therefore the modulation of the turbulent kinetic energy around the surface of the particles. The iso-contours of K/ K sur form almost spherical shells around the surface of the particles and are slightly elongated towards the compressional quadrants. The value of K/ K sur is around 0.88 close to the surface of the particles and goes asymptotically to 1 over a distance of ≈ 0.45D eq . Figure 14(d-f ) shows the contours of the production P and figure 14(g-i) those of the dissipation rate . The production is stronger in the extensional quadrants with a value of P/ P sur = 1.1 for the spheres and 1.3 for the oblates close to the surface of the particles. On the other hand, in the compressional quadrants the production is lower than the mean value and is reduced to P/ P sur = 0.85 for the spheres and 0.7 for the oblates.
As regards the dissipation, previous experimental and numerical studies have shown that the energy dissipation is locally enhanced near the surface of the particles (see e.g. Lucci et al. 2010;Tanaka & Eaton 2010). Our results (see figure 14g-i) show the same trend. The dissipation rate is enhanced near the surface of the particles and reaches values five times larger than the mean value in the compressional quadrants in the case of spherical particles. The enhancement of the dissipation rate is more uniform around the surface of the oblates with a value of / sur ≈ 2.3 uniformly in all directions. Tanaka (2017) reported similar behaviour for the modulation of the production P and the dissipation rate in transient HST and dilute suspensions of spherical particles. More details of the origins of flow modulation by the solid phase can be obtained by the flow structures conditionally averaged around the surface of the particles. Figure 15 shows the fluctuating velocity vectors in the x-y plane, passing through the centre of the reference particle. In the plot, the velocity vectors are superposed onto the contours of the magnitude of the fluctuating velocity vector |u | = √ u 2 + v 2 + w 2 ; contours are normalised by the velocity vector magnitude, conditionally averaged in the cubic box of side 3D eq around the particles, |u | sur . Figure 15(a,b) reports the results for cases Sp5 and Ob5. For both particle shapes, four vortical structures are visible around the reference particle. To discuss these flow structures in relation to the compressional/extensional quadrants around the surface of the particle more easily, a schematic is presented in figure 15(c), where the compressional/extensional quadrants are denoted by χ 1-4 , while Γ 1-4 indicate the vortical structures. Figure 15 reveals that the averaged flow field around the particle surface consists of two pairs of counter-rotating vortices, with the induced velocity field reminiscent of the one induced by the force dipole (stresslet) in the vanishing inertia regime (see e.g. Graham 2018). As shown previously in figure 14(d-f ), the production is highly modulated in all quadrants χ 1-4 . Hence, any pair of adjacent counter-rotating spanwise vortices formed around the particle surface modulate the tangential Reynolds shear stress and the turbulence production in the area between them in all quadrants. In previous studies, Kida & Tanaka (1994) and Ahmed & Elghobashi (2000) showed that in HST the regions of high energy production are mainly enclosed by quasi-streamwise counter-rotating vortices, whereas here the vortices are mainly parallel to the spanwise direction.
Further, Ahmed & Elghobashi (2000) showed that particles increase the dissipation rate by creating local velocity gradients which leads to an increase of the local strain and dissipation rates. In particular, the increase of the extensional strain has a prominent effect on the enhancement of the dissipation rate , although compressive strain is also increased by the presence of the solid particles. Our results show an increase of the dissipation rate in all quadrants (see figure 14g-i), but mainly in χ 2 and χ 4 where the vortex pairs Γ 1-2 and Γ 3-4 are inducing extensional strain rate.
Note that for the oblate particles (figure 15b), the velocity magnitude in quadrants χ 1 and χ 3 is slightly larger than that in the other quadrants. This difference appears because the dynamics of the oblate particles and the flow field around them reflect their anistropic shape; i.e. the particles sample different flow patterns based on their orientation. To demonstrate this, we display in figure 16 the joint probability density of the particle orientation and the spanwise angular velocity. Here, the horizontal axis indicates the streamwise component of the particle orientation vector multiplied by the sign function of the normal component, so that positive values correspond to instances when the oblate particles are perpendicular to the mean shear direction and negative values to instances when they are aligned with mean shear. The data in the figure reveal the higher probability of sampling particles that are inclined with the shear direction when they experience larger (c) Schematic showing the resulting vortical structures formed from the interaction between a rotating particle and a linear shear flow; the compressional and extensional quadrants around the particle are marked by χ 1-4 and the vortical structures by Γ 1-4 . angular velocities. This creates a small asymmetry in the flow pattern that is superposed on top of the symmetric flow pattern that appears around the spherical particles (see figure 15a). This is consistent with the findings of Ardekani et al. (2017) for the force and torque experienced by oblate particles in the vicinity of a wall where the shear rate is strong.

Final remarks
We have used interface-resolved simulations to study the turbulence modulation of particle-laden HST at statistically steady state. Neutrally buoyant spherical and oblate particles (AR = 1/3) are investigated up to volume fractions of 20 % and 10 %, respectively. A particle size of ≈ 20 Kolmogorov length scales in the single-phase flow is considered, which is equal to λ/D eq ≈ 0.9 in terms of the Taylor microscale.
We show that a statistically steady state is not exclusive to single-phase HST as the production and dissipation rates of the turbulent kinetic energy are observed to be statistically in balance in the studied particle-laden cases. The simulations show a non-monotonic trend of the turbulent kinetic energy with increasing volume fraction. Both particle shapes decrease the turbulent kinetic energy with respect to the single phase up to a FIGURE 16. Joint p.d.f. of the particle orientation (streamwise component of the particle orientation vector times the sign function of the normal one) and the particle spanwise angular velocity for the case Ob5.
certain volume fraction before the turbulence activity enhances at higher volume fractions. This increase of the turbulent kinetic energy is observed for both spherical and oblate particles, starting, however, at a lower volume fraction for the latter (20 % for spheres and 10 % for oblate particles). Analysis of the turbulent kinetic energy budget shows a very similar attenuation and then an increase of the production and dissipation rates of the turbulent energy, while the contribution from the interphase interaction remains below 3 %. We have documented two main mechanisms for turbulence modulation in the presence of particles. First we show that the dissipation rate is significantly larger in the vicinity of solid particles (see figure 14) due to the enlargement of extensional strain rate in these regions (see also Ahmed & Elghobashi 2000). This leads to turbulence attenuation as observed in the results of this study for lower volume fractions. This is in line with the numerical study of Tanaka & Teramoto (2015) of transient HST, where the turbulence kinetic energy is attenuated in the dilute regime. In agreement, the additional reduction of the turbulence activity for oblate particles at lower volume fractions reported here can be related to the larger surface area of the oblate particles with respect to spheres at fixed volume. On the other hand, we also observe that two pairs of counter-rotating vortices form around the particle surface, thus increasing the production of turbulence kinetic energy, locally. Theses vortices, which have also been observed in the study by Tanaka & Teramoto (2015), increase the Reynolds shear stress in the region between them (Tanaka 2017). When the volume fraction is large enough, this mechanism (see figure 15) and the particle-particle interactions (see figure 11b) generate strong velocity fluctuations (see figure 9) that overall enhance the turbulence activity. Oblate particles, on the other hand, have higher rotation rates compared to spheres and also a tendency to rotate more in the opposite direction to the mean shear flow, which increases mixing and the level of turbulent fluctuations more than in the case of spheres. Even though clustering is not observed for oblate particles, it should be interpreted as a sign of strong collisions between them as they depart from each other rapidly after a collision takes place. Also, due to their shape, oblate particle pairs can have interactions at smaller separation distance between their centres (see figure 12). The associated strong fluctuations can be observed in the tails of the p.d.f.s of particle linear and angular velocities (see figures 9 and 10).
As discussed in previous works on SS-HST, this configuration reproduces the dynamics of the equilibrium logarithmic layer in wall-bounded turbulence (Sekimoto et al. 2016). As a consequence, we can interpret the results presented here as indicative of the behaviour of small particles, smaller than the logarithmic layer in wall-bounded high-Reynolds turbulence, but still finite size in viscous units. We recall that the simulations of finite-size oblate particles, of size of the order of the streak spacing and logarithmic layer width, in turbulent channel flow show that oblate particles with small enough aspect ratio experience small rotational rates while aligning with the wall in its vicinity . This shields the walls from the bulk region and thus induces turbulence attenuation by lowering the wall-normal velocity fluctuations. This decrease in the turbulence activity manifests itself as drag reduction in wall-bounded flows. We show in this study that oblates in HST do not show any noticeable clustering or preferential orientation and they instead increase the turbulence activity. Therefore we can speculate that as the size of oblate particles reduces with respect to the logarithmic-layer extension, the observed drag reduction would disappear. This is expected to occur when the Reynolds number is relatively large, so in a range not yet reachable by interface-resolved numerical simulations; nevertheless, this prediction is consistent with the recent experimental measurements of Zade (2019) in high-Reynolds pipe flow laden with oblate particles. TABLE 2 (cntd). First four central moments of the density probability functions of the fluid and particle velocity components, together with the Reynolds shear stress u v , the particle kinetic energy and magnitude of the angular velocity vector for different cases. Values smaller than 10 −3 in magnitude are set to zero.