1. Introduction
The turbulent flow of suspensions over porous media is crucial in various environmental phenomena and industrial technologies. These include sedimentation in rivers, pyroclastic flows, flow through unconsolidated sand (Zhai et al. Reference Zhai2005), and hydraulic fracturing (Song, Jinzhou & Yongming Reference Song, Jinzhou and Yongming2014). In these applications, both inertia and the presence of a porous surface play a critical role in the suspension dynamics and fluid–solid interactions. This study therefore focuses on understanding how the permeability of a porous medium affects the dynamics of suspensions and the behaviour of particles in a turbulent channel flow with fixed porosity and porous layer thickness.
1.1. Turbulent flows over porous media
The flow over and through a porous medium with a complex structure can be examined by modelling the porous medium as a continuum with given physical properties (i.e. porosity and permeability). The porosity 
$\varepsilon$ of a porous medium is defined as the ratio of its void volume to its bulk volume. Fractions and percentages can be used to express this dimensionless quantity. A porous material's permeability 
$K$, on the other hand, determines how well fluids can pass through it. A porous material with high permeability allows fluids to flow more easily than one with low permeability. The permeability of a porous medium is related to the porosity, shape and level of connections of the pores (Bear Reference Bear1988).
One of the first studies on flow through porous media was published by Darcy (Reference Darcy1856). In his experiments, he demonstrated that the average flow rate flowing through sand (porous media) is proportional to the pressure gradient across it. Darcy's experimental data demonstrated clearly a linear relationship between average discharge rate and pressure gradient. Darcy's law consists of two terms: one term can be interpreted as an average drag felt by the fluid as it passes through the porous medium and is balanced by the second term, which is the pressure gradient term. For the low porous Reynolds number for which Darcy's law holds, i.e. 
$Re_p \ll 1$, the flow is governed by Stokes’ equation, in which the viscous term is responsible for the Darcy drag. A deviation from linearity is observed as 
$Re$ increases. A nonlinear relationship was proposed by Forchheimer, with a quadratic drag term that better fits experimental data; the equation is called Forchheimer equation (Joseph, Nield & Papanicolaou Reference Joseph, Nield and Papanicolaou1982; Whitaker Reference Whitaker1996; Giorgi Reference Giorgi1997; Lage Reference Lage1998). As a result of inertial effects in the pores, the quadratic term appears in this equation. When the flow is no longer governed by Stokes’ equation, the quadratic drag term is due to the convective term, 
$(u \boldsymbol {\cdot } \boldsymbol {\nabla }) u$, in the Navier–Stokes equations.
Subsequently, Brinkman (Reference Brinkman1947) suggested extending Darcy's law to high-permeability materials by adding a sort of viscous term in Stokes’ equation. The Brinkman model has been studied over the years by different researchers interested in many different applications (see, among others, Neale & Nader Reference Neale and Nader1974; Vafai & Thiyagaraja Reference Vafai and Thiyagaraja1987; Vafai & Kim Reference Vafai and Kim1990; Hahn, Je & Choi Reference Hahn, Je and Choi2002; Mirbod, Andreopoulos & Weinbaum Reference Mirbod, Andreopoulos and Weinbaum2009; Mirbod, Wu & Ahmadi Reference Mirbod, Wu and Ahmadi2017; Wu & Mirbod Reference Wu and Mirbod2018, Reference Wu and Mirbod2019; Haffner & Mirbod Reference Haffner and Mirbod2020; Kang & Mirbod Reference Kang and Mirbod2021; Hooshyar, Yoshikawa & Mirbod Reference Hooshyar, Yoshikawa and Mirbod2022). Importantly, Beavers & Joseph (Reference Beavers and Joseph1967) improved the Darcy equation to introduce interface jump conditions between liquid and porous media. In 1996, Whitaker derived the so-called volume-averaged Navier–Stokes (VANS) equations from the Navier–Stokes equations in a porous media (Whitaker Reference Whitaker1996). The VANS equations contain new terms, including the Darcy term, an average adhesion resistance, and the Forchheimer term, a second-order tensor due to the inertial effects of the resistance function of the geometric structure of the porous medium. It should be noted, however, that defining a model for the fluid-porous interface remains a challenge. In fact, the structure of a porous medium modifies rapidly within the small layer next to the fluid-porous interface, i.e. the so-called Brinkman layer (Neale & Nader Reference Neale and Nader1974); therefore, the Navier–Stokes averaging over a volume containing the interface is not consistent. To properly connect the VANS and Navier–Stokes equations, Ochoa-Tapia & Whitaker (Reference Ochoa-Tapia and Whitaker1995) introduced momentum-transfer conditions, which improved the computation of the transferred stress to both the fluid phase and the porous media. Here, in particular, we employ the Navier–Stokes equation for the free fluid flow and the VANS equations inside the porous media. These equations are connected using the interface jump condition between liquid and porous media as proposed by Ochoa-Tapia & Whitaker (Reference Ochoa-Tapia and Whitaker1995). It is worth noting that over the past 20 years, several groups have also tried to model the interface layer and control the momentum transfer (see e.g. Alazmi & Vafai Reference Alazmi and Vafai2001; Goyeau et al. Reference Goyeau, Lhuillier, Gobin and Velarde2003; Chandesris & Jamet Reference Chandesris and Jamet2006; Valdés-Parada, Goyeau & Ochoa-Tapia Reference Valdés-Parada, Goyeau and Ochoa-Tapia2007).
Both numerical and experimental works have been conducted to study the effect of isotropic and anisotropic porous materials on turbulent flows. Using a continuum approach for porous media, laminar and turbulent flows were investigated by Breugem, Boersma & Uittenbogaard (Reference Breugem, Boersma and Uittenbogaard2004). The authors showed that the permeability of the porous wall depends not only on the porosity but also on the gradient of the volume-averaged velocity inside the porous media. The same group later used direct numerical simulations (DNS) to investigate the turbulent channel flow over three-dimensional (3-D) regular cubes mimicking porous structures (Breugem & Boersma Reference Breugem and Boersma2005). They showed quantitatively that the VANS equations are able to simulate accurately the turbulent flows over and through permeable walls. They also analysed the impact of a packed-bed porous wall on the turbulent channel flow (Breugem, Boersma & Uittenbogaard Reference Breugem, Boersma and Uittenbogaard2006) and found the near-wall structures to be significantly different: streamwise-correlated streaks almost disappear, while the spanwise vorticity near the wall increases for a highly permeable surface. This is consistent with the decrease in the peak value of the streamwise root-mean-square (r.m.s.) velocity normalized by the friction velocity at the permeable wall. Furthermore, they showed that the turbulence near a highly permeable wall is dominated by relatively large vertical structures generated by a Kelvin–Helmholtz instability. Turbulent flows over porous media considering the porous layer as a continuum, without a need for knowledge of the porous media microstructure, have also been examined by Rosti, Cortelezzi & Quadrio (Reference Rosti, Cortelezzi and Quadrio2015) and Rosti, Brandt & Pinelli (Reference Rosti, Brandt and Pinelli2018). These authors also found that permeability plays a critical role and showed that porous materials with high wall-normal and low wall-parallel permeabilities are characterized by increased turbulence isotropy (Rosti et al. Reference Rosti, Brandt and Pinelli2018). In a more recent work, Rosti et al. (Reference Rosti, Pramanik, Brandt and Mitra2020) studied the flow inside a deformable porous medium for the case of a viscous fluid and an incompressible viscoelastic material. It was found that in a poroelastic medium, the flow permeability is a function of the elastic shear modulus of the deformable surface and of the imposed pressure difference, suggesting that the dynamics of the flow over such a complex medium cannot be deduced directly from studies over rigid porous walls.
While most of the previous work discussed above considered isotropic porous media, anisotropic porous layers have been found to be effective for passive control of turbulence (see e.g. Gomez-de Segura, Sharma & García-Mayoral Reference Gomez-de Segura, Sharma and García-Mayoral2017; Rosti et al. Reference Rosti, Brandt and Pinelli2018). Kuwata & Suga (Reference Kuwata and Suga2017) used multiple-relaxation-time lattice Boltzmann DNS to simulate the turbulent flow over a porous layer. Considering anisotropic permeable walls and squared pore arrays aligned with the Cartesian axes, they found that streamwise and spanwise permeabilities considerably enhance the Reynolds stresses over the porous walls when compared to the vertical permeability. The turbulent channel flow over porous layers has also been studied experimentally by Suga et al. (Reference Suga, Matsumura, Ashitaka, Tominaga and Kaneda2010) and Suga, Nakagawa & Kaneda (Reference Suga, Nakagawa and Kaneda2017); the authors showed that the spanwise scales of the measured structures can be reasonably correlated with the wall-normal distance plus the zero-plane displacement, estimated from the mean velocity profile.
1.2. Turbulent suspension flows over impermeable walls
Particle-laden turbulent flows over impermeable walls have been investigated extensively. Many earlier studies considered small heavy particles in the very dilute regime. It has been shown that inertia induces particle preferential sampling and migrations that lead to clustering in both homogeneous and inhomogeneous flows (see, among many, Bec et al. Reference Bec, Biferale, Cencini, Lanotte, Musacchio and Toschi2007; Toschi & Bodenschatz Reference Toschi and Bodenschatz2009; Sardina et al. Reference Sardina, Schlatter, Brandt, Picano and Casciola2012).
Pan & Banerjee (Reference Pan and Banerjee1996) were the first to study the effect of finite-size particles in a turbulent channel flow at very low concentrations and showed that smaller particles decrease turbulence forces and Reynolds stresses, while larger particles increase turbulence intensities and stresses. These findings were in line with the other works that considered the effect of the particles on the spectral energy distribution and the turbulent cascade (Francesco Lucci & Elghobashi Reference Francesco Lucci and Elghobashi2010; Yeo et al. Reference Yeo, Dong, Climent and Maxey2010). The appearance of different numerical algorithms has facilitated an increasing number of studies on finite-size particles in turbulence. Among these, Kidanemariam et al. (Reference Kidanemariam, Chan-Braun, Doychev and Uhlmann2013) examined open channel flow with heavy particles of finite size in the dilute regime, and showed that the solid phase is preferably assembled near the wall in low-speed regions. Later, the same approach was used by these authors to study pattern formation in turbulent channel flow over a mobile sediment bed (Kidanemariam & Uhlmann Reference Kidanemariam and Uhlmann2014). Concerning neutrally buoyant particles in channel flow, Picano, Breugem & Brandt (Reference Picano, Breugem and Brandt2015) studied neutrally buoyant spheres, while Ardekani & Brandt (Reference Ardekani and Brandt2019) considered non-spherical particles, showing a reduced drag. The effect of particles on the critical Reynolds number for turbulence onset in semi-dilute suspension flows was studied in Matas, Morris & Guazzelli (Reference Matas, Morris and Guazzelli2004), Loisel et al. (Reference Loisel, Abbas, Masbernat and Climent2013) and Yu et al. (Reference Yu, Wu, Shao and Lin2013), and the chaotic flow of relatively dense suspensions of neutrally buoyant and heavy particles in channel flows, for a volume fraction of 7 %, was examined in Shao, Wu & Yu (Reference Shao, Wu and Yu2012). Vowinckel, Kempe & Fröhlich (Reference Vowinckel, Kempe and Fröhlich2014) reported on various regimes that differ depending on the particle buoyancy. Recent studies conducted by us have examined both suspensions and solvent, and their related instabilities over and through various porous media models, in Bagheri et al. (Reference Bagheri, Moradi Bilondi, Taheri, Anand, Schatz and Mirbod2022), Haffner et al. (Reference Haffner, Kang, Shapley and Mirbod2019), Mirbod & Shapley (Reference Mirbod and Shapley2022), Bagheri, Kang & Mirbod (Reference Bagheri, Kang and Mirbod2019), Moradi Bilondi et al. (Reference Moradi Bilondi, Mastroianni, Brandt and Mirbod2022) and Kang & Mirbod (Reference Kang and Mirbod2021). However, the impact of different porous media on non-Brownian, non-colloidal, neutrally buoyant, rigid, spherical particles flowing in a channel remains unknown.
Thanks to recent developments in algorithms and computational power, suspension flows with inertia can be studied by interface-resolved DNS, with ten thousand finite-size particles. Using fully resolved DNS and the immersed boundary method (IBM) for neutrally buoyant spheres as conducted here, Lashgari et al. (Reference Lashgari, Picano, Breugem and Brandt2014) reported a transition akin to that found in the pure fluid case at low volume fractions, i.e. from the laminar base state directly to turbulence as the Reynolds number was increased. However, a ‘shear-thickening’ regime (distinct from laminar flow and turbulence) was identified when the volume fraction was sufficiently large, dominated by the particle-induced stresses. Later, Picano et al. (Reference Picano, Breugem and Brandt2015) and Costa et al. (Reference Costa, Picano, Brandt and Breugem2016) observed that under fully turbulent conditions, significant particle-induced stress first appears at lower solid volume fraction 
${{\varPhi _b}}$ in a layer close to the wall of size similar to the particle diameter; with increasing 
${{\varPhi _b}}$, the wall layer is characterized by large particle-induced stress and found to be responsible for the increased drag. These authors also reported that because of these non-homogeneities in concentration distribution, the increase in the overall drag cannot be explained simply by using the notion of an increased effective suspension viscosity. The study by Ardekani et al. (Reference Ardekani, Costa, Breugem, Picano and Brandt2017) considered non-spherical oblate particles and showed that the drag is reduced, as these particles do not accumulate in the particle layer but remain away from the wall, confirming the importance of the particle dynamics in the near-wall region for the global flow behaviour (Peng, Ayala & Wang Reference Peng, Ayala and Wang2019, for the role of particle rotation, see also). Numerous DNS have also been performed to examine the sediment transport of particle-laden flows over mobile and immobile beds (Herwig, Kempe & Fröhlich Reference Herwig, Kempe and Fröhlich2011; Kempe, Vowinckel & Fröhlich Reference Kempe, Vowinckel and Fröhlich2014; Vowinckel et al. Reference Vowinckel, Kempe and Fröhlich2014, Reference Vowinckel, Nikora, Kempe and Fröhlich2017, Reference Vowinckel, Biegert, Luzzatto-Fegiz and Meiburg2019a,Reference Vowinckel, Withers, Luzzatto-Fegiz and Meiburgb, Reference Vowinckel, Biegert, Meiburg, Aussillous and Guazzelli2021; Biegert, Vowinckel & Meiburg Reference Biegert, Vowinckel and Meiburg2017; Jain, Tschisgale & Fröhlich Reference Jain, Tschisgale and Fröhlich2020, Reference Jain, Tschisgale and Fröhlich2021; Köllner et al. Reference Köllner, Meredith, Nokes and Meiburg2020; Papadopoulos et al. Reference Papadopoulos, Nikora, Vowinckel, Cameron, Jain, Stewart, Gibbins and Fröhlich2020). Few studies, on the other hand, have investigated the particle flow over complex walls. Recently, Rosti, Mirbod & Brandt (Reference Rosti, Mirbod and Brandt2021) have investigated the rheology of particle suspensions flowing over porous walls in a plane Couette flow and in a laminar regime; they found that the porous walls induce a progressive decrease in the suspension effective viscosity as the wall permeability increases. This behaviour was explained by the weakening of the wall-blocking effect and by the appearance of a slip velocity at the interface of the porous medium.
1.3. Outline
This study employs DNS to simulate turbulent suspension flows of monodisperse, rigid, spherical particles over various porous media where particles interact with the porous layer but do not enter into it. The solid volume fraction 
${{\varPhi _b}}$ considered varies in the range 0–0.2. We fully describe the solid phase dynamics via an IBM for the finite-size particles, with a ratio between the radius 
$a$ and the channel half-width 
$h$ equal to 
$1/10$, while solving the VANS equations for flow inside the porous media. The porous layer porosity and thickness are kept constant, while their permeability varies. The organization of the paper is as follows. Section 2 shows mathematical formulations, numerical procedures, system configuration and related parameters. The results and discussions are presented in § 3. Section 4 addresses the summary and conclusion of the study.
2. Methodology
2.1. Governing equations
We study the fully developed flow of monodispersed suspensions of rigid, spherical particles of radius 
$a$ suspended in a Newtonian fluid in a channel bounded by identical, rigid, homogeneous and isotropic porous media. The incompressible Navier–Stokes equations are used to describe the Eulerian fluid phase:
\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol u}_f=0, \end{gather}
\begin{gather} \frac{\partial {\boldsymbol u}_f}{\partial t}+ {\boldsymbol u}_f \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol u}_f={-}\frac{1}{\rho}\,\boldsymbol{\nabla} p+\frac{\mu}{\rho}\, {\nabla}^2 {\boldsymbol u}_f+ {\boldsymbol f}. \end{gather}
Here, 
${\boldsymbol u}_f$ and 
$p$ represent the velocity and pressure fields, respectively, and the term 
${\boldsymbol f}$ is a generic force field. In addition, the parameters 
$\mu$ and 
$\rho$ are the dynamic viscosity and the density of the fluid phase, respectively, and their ratio 
$\nu = {\mu }/{\rho }$ indicates the kinematic viscosity of the fluid. The centroid linear and angular velocities of neutrally buoyant particles, 
${\boldsymbol u}_p$ and 
${\boldsymbol \omega }_p$, are governed by the Newton–Euler Lagrangian equations, defined as
\begin{gather} \rho_{p} V_{p}\,\frac{{\rm d}{\boldsymbol u}_p}{{\rm d}t}={\oint_{\partial V_p} ({-}p{\boldsymbol I} +\mu (\boldsymbol{\nabla} {\boldsymbol u}_f+\boldsymbol{\nabla} {\boldsymbol u}_f^{\rm T}))\boldsymbol{\cdot} {\boldsymbol n}\,{\rm d}A}+{\boldsymbol F}_p, \end{gather}
\begin{gather} I_{p}\,\frac{{\rm d}{\boldsymbol \omega}_p}{{\rm d}t}={\oint_{\partial V_p} {\boldsymbol r}\times ({-}p{\boldsymbol I} +\mu (\boldsymbol{\nabla} {\boldsymbol u}_f+\boldsymbol{\nabla} {\boldsymbol u}_f^{\rm T}))\boldsymbol{\cdot} {\boldsymbol n}\,{\rm d}A}+{\boldsymbol T}_p, \end{gather}
where 
$\rho _p$ indicates the particle density, and 
$V_p$ is the particle volume, defined as 
$4{\rm \pi} a^3/3$. In addition, 
${\boldsymbol I}$ is the unit tensor, 
$I_p=\tfrac {2}{5}\rho _p V_p a^2$ is the moment of inertia, 
${\boldsymbol r}$ is the distance vector from the centre of the particle, and 
${\boldsymbol n}$ is the unit vector normal to the particle surface 
$\partial V_p$. The terms 
${\boldsymbol F}_p$ and 
${\boldsymbol T}_p$ represent the force and torque acting on the particles as a result of the particle–particle and particle–wall interactions. To couple the solid and fluid phases, an IBM is used. This accounts for a body force 
${\boldsymbol f}$ added to the right-hand side of (2.2) to mimic the actual no-slip and no-penetration boundary conditions at the moving particle surface (Breugem Reference Breugem2012; Picano et al. Reference Picano, Breugem and Brandt2015; Ardekani et al. Reference Ardekani, Costa, Breugem and Brandt2016):
\begin{equation} {\boldsymbol u}_f|_{\partial V_p} = {\boldsymbol u}_p + {\boldsymbol \omega}_p \times {\boldsymbol r}. \end{equation} We denote the porous layer with the porosity 
$\varepsilon$ and the permeability with the tensor 
$K_{ij}$. When the porous medium is isotropic, the permeability reduces to a single scalar value, i.e. 
$K$. In this study, we employ the VANS equations to describe the fluid motion in the porous medium, as proposed by Whitaker (Reference Whitaker1969, Reference Whitaker1986). For isotropic porous media with negligible fluid inertia and large-scale separation, the VANS equations can be stated as (Kang & Mirbod Reference Kang and Mirbod2019; Rosti et al. Reference Rosti, Mirbod and Brandt2021)
\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \langle {\boldsymbol u}_f \rangle_s =0, \end{gather}
\begin{gather} \frac{\partial \langle {\boldsymbol u}_f \rangle_s}{\partial t}={-}\frac{\varepsilon}{\rho}\, \boldsymbol{\nabla} \langle p \rangle_f +\nu\,{\nabla}^2 \langle {\boldsymbol u}_f \rangle_s-\frac{\nu \varepsilon}{K}\,\langle {\boldsymbol u}_f \rangle_s , \end{gather}
Notably, the VANS equations are based on two different averaging methods (see Ochoa-Tapia & Whitaker Reference Ochoa-Tapia and Whitaker1995), namely the superficial volume average 
$\langle \psi \rangle _s = ({1}/{V})\int _{V_f}\psi \,{\rm d}V_f$, and the intrinsic volume average 
$\langle \psi \rangle _f = ({1}/{V_f})\int _{V_f}\psi \,{\rm d}V_f$, where the superficial and intrinsic volume averages are linearly related to each other using the condition 
$\varepsilon$, 
$\langle \psi \rangle _s = ({V_f}/{V})\langle \psi \rangle _f = \varepsilon \langle \psi \rangle _f$. Here, 
$V_f$ is the volume occupied by the fluid within the averaging volume 
$V$. Note that superficial and intrinsic volume averages are employed for the velocity and pressure fields as reported by Quintard & Whitaker (Reference Quintard and Whitaker1994) and Whitaker (Reference Whitaker1996).
In our analysis, we assume pressure and velocity continuity at the suspension-porous interface, while the shear stress displays a jump that can be controlled by the parameter 
$\tau$ (Ochoa-Tapia & Whitaker Reference Ochoa-Tapia and Whitaker1995; Goyeau et al. Reference Goyeau, Lhuillier, Gobin and Velarde2003; Rosti et al. Reference Rosti, Mirbod and Brandt2021). The value of 
$\tau$ determines the stress transfer between the suspending fluid and the porous media. Consequently, the boundary conditions at the fluid-porous interface can be written as
\begin{gather} {\boldsymbol u}_f = \langle {\boldsymbol u}_f \rangle_s, \quad p = \langle p \rangle_f , \end{gather}
\begin{gather} \sigma\left(\frac{\partial u_f}{\partial y}-\frac{1}{\varepsilon}\,\frac{\partial \langle u_f \rangle_s}{\partial y}\right) ={\pm} \tau u_f, \end{gather}
\begin{gather} \sigma\left(\frac{\partial w_f}{\partial y}-\frac{1}{\varepsilon}\,\frac{\partial \langle w_f \rangle_s}{\partial y}\right) ={\pm} \tau w_f . \end{gather}
In (2.8a,b), 
${\boldsymbol u}_f$ is the fluid flow velocity vector, while 
${u_f}$ and 
${w_f}$ in (2.9) and (2.10) are the streamwise and spanwise components of the fluid velocity, respectively. In (2.9) and (2.10), on the right-hand side, a positive sign indicates that the pure fluid region is above the porous media, while a negative sign indicates the opposite. The dimensionless parameter 
$\sigma = {\sqrt {K}}/{h}$ is the so-called Darcy number or the non-dimensional permeability parameter, where 
$h$ is the half-channel height. The coefficient 
$\tau$ represents the stress transfer between the suspending fluid and the porous media, and takes varying values (zero, positive or negative) depending on how the porous material's structure varies in the heterogeneous transition layer, and how the interface is fabricated. A detailed comparison of the different boundary conditions can be found in (Alazmi & Vafai Reference Alazmi and Vafai2001). There were also several attempts to calculate 
$\tau$ theoretically (Goyeau et al. Reference Goyeau, Lhuillier, Gobin and Velarde2003; Deng & Martinez Reference Deng and Martinez2005; Min & Kim Reference Min and Kim2005; Chandesris & Jamet Reference Chandesris and Jamet2006, Reference Chandesris and Jamet2007; Valdés-Parada et al. Reference Valdés-Parada, Goyeau and Ochoa-Tapia2007, Reference Valdés-Parada, Alvarez-Ramírez, Goyeau and Ochoa-Tapia2009; Carotenuto & Minale Reference Carotenuto and Minale2011). Recently, also attempts to measure the jump coefficient were carried out by Carotenuto & Minale (Reference Carotenuto and Minale2011) and by Bagheri & Mirbod (Reference Bagheri and Mirbod2022). In the present study, we assume 
$\tau =+1$, 
$0$ and 
$-1$. Note that 
$\tau = 0$ implies that the stress carried by the suspending fluid is fully transferred to the fluid flow inside the porous layer. More details can be found in § 3.3.
2.2. Numerical implementation
Our numerical simulations are based on a 3-D solver that adopts an IBM for coupling the fluid and particle phases in the pure Newtonian fluid region, 
$0 < y < 2h$. The IBM was first developed by Peskin (Reference Peskin1972), who simulated blood flow patterns around heart valves. Since then, it has been widely modified and improved, as reported in Mittal & Iaccarino (Reference Mittal and Iaccarino2005). Among others, Uhlmann (Reference Uhlmann2005) proposed a direct forcing IBM and showed the ability of the method to simulate large-scale configurations with many particles. Here, we use the IBM as first proposed by Breugem (Reference Breugem2012), while the various modifications introduced by this author enabled us to also reach second-order accuracy in space. It is worth noting that the IBM has first-order accuracy in the local velocity field, and may have apparent second-order accuracy in force/torque as reported by Peng & Wang (Reference Peng and Wang2020). Moreover, the particle–particle and particle–wall interaction force and torque in (2.3) and (2.4) include a lubrication correction and a soft collision model as reported in Costa et al. (Reference Costa, Boersma, Westerweel and Breugem2015). In particular, we use Brenner's asymptotic solution (Brenner Reference Brenner1961) to correct the lubrication force when the distance between solid objects is less than a certain threshold and cannot be resolved accurately by the numerical mesh. At such small distances, surface roughness is accounted for by saturating this force. When spherical particles are in contact, both the normal and tangential contact force components are obtained from the overlap and the relative velocity. It is worth noting that in this analysis, we considered that the particles do not enter the porous layer. In other words, the particles approach an impermeable wall at the suspension-porous interface where the above lubrication and collision models can be applied.
The governing differential equations are solved on a staggered grid using a second-order central finite-difference scheme, while the dispersed phase is represented by a set of Lagrangian points that are distributed uniformly on the surface of each particle. The forces were computed at each Lagrangian grid point on the particle surface based on the difference between the computed velocity at the particle surface and the first prediction velocity of the fluid at the same point in space; the first prediction velocity was obtained by time integration of the free-flow equations before imposing the continuity constraint. The forces are extrapolated on the Eulerian grid used for the momentum equations and added to the first prediction velocity, followed by the correction pressure scheme (Breugem Reference Breugem2012). For the time integration, we use the explicit fractional-step method, while all the other terms are advanced with the third-order Runge–Kutta scheme. More details in terms of the numerical procedures for both the IBM and VANS equations and the code validations can be found in Lashgari et al. (Reference Lashgari, Picano, Breugem and Brandt2014), Picano et al. (Reference Picano, Breugem and Brandt2015), Ardekani et al. (Reference Ardekani, Costa, Breugem and Brandt2016), Rosti & Brandt (Reference Rosti and Brandt2017), Izbassarov et al. (Reference Izbassarov, Rosti, Ardekani, Sarabian, Hormozi, Brandt and Tammisola2018) and Rosti et al. (Reference Rosti, Mirbod and Brandt2021).
2.3. System configuration
Herein, we investigate turbulent suspensions flowing in a channel between two infinite walls located at 
$y = -h_p$ and 
$y = 2h+h_p$. The two porous layers have thickness 
$h_p$ placed at the top and bottom of the channel, while the suspension-porous interface is located at 
$y = 0$ and 
$y = 2h$; see figure 1. Periodic boundary conditions are applied in the wall-parallel directions over a domain of size 
$L_x = 6h$, 
$L_y = 2(h + h_p)$ and 
$L_z = 3h$ in the streamwise, wall-normal and spanwise directions, respectively. A time-dependent pressure gradient acting in the streamwise direction imposes a constant bulk velocity 
$U_0$ corresponding to a constant bulk Reynolds number 
$Re = {2hU_0}/{\nu }=5600$. For the simulations presented here, we use 960 grid points in the streamwise direction 
$x$, 
$320+32$ grid points in the wall-normal direction 
$y$ (where 32 is the number of points used for the porous layers), and 480 points in the spanwise direction 
$z$. The non-dimensional parameters and their values are as follows. The non-dimensional permeability parameter is 
$\sigma \in [0.63, 6.3] \times 10^{-3}$, while the ratio between the porous and the channel height is 
$h_p/h=1/10$, and the porosity of the porous layers is constant at 0.6. The ratio of particle diameter to channel height is 
$a/h=1/10$, and the particle volume fraction ranges from 
$0$ to 
$20\,\%$. Table 1 reports all different conditions considered in this study. For the sake of comparison, the simulation parameters are chosen to be similar to those in previous works on single-phase turbulent channel flows over porous walls, and turbulent suspension flows in smooth channels, as reported in Rosti et al. (Reference Rosti, Cortelezzi and Quadrio2015), Picano et al. (Reference Picano, Breugem and Brandt2015) and Rosti & Brandt (Reference Rosti and Brandt2017). All the simulations are started from a fully developed turbulent channel flow. After the flow has reached statistical steady state, the calculations are continued for an interval 
$480h/U_0$ time units, during which 
$96$ full flow fields are stored for further statistical analysis. To verify the convergence of the statistics, we have computed them using different numbers of samples, and verified that the differences are negligible.
Figure 1. Sketch of the 3-D channel flow with 50 particles overlying porous media followed by the coordinate system adopted in this study.
Table 1. List of simulations performed. The porosity is fixed at 
$\varepsilon =0.6$, and the bulk Reynolds number is 
$Re=5600$, whereas the non-dimensional permeability parameter and the solid volume fraction vary as indicated. The parameter 
$N_p$ indicates the number of particles at each particle volume fraction.
3. Results
We study the turbulent suspension flow of monodispersed, neutrally buoyant, rigid, spherical particles over porous media in a channel where the particles do not enter inside the porous layers. Figure 2 compares the contours of the normalized streamwise velocity component from instantaneous snapshots for the three different volume fractions under investigation. From the qualitative results shown in the figure, it can already be seen that the suspension flow streamwise velocity is sensitive to the particle volume fraction 
${{\varPhi _b}}$. In particular, the turbulence activities are greatly enhanced throughout the channel due to the disturbances from the particles inside the channel. As the particle volume fraction increases, the streamwise velocity increases at the centre of the channel; however, the velocity intensity decreases near the interface, suggesting that the slip velocity decreases when increasing the number of particles in the flow. Visualizations of the normalized streamwise velocity are shown in figure 3 for three different wall non-dimensional permeability parameters 
$\sigma$. The wall permeability also has a critical impact on the turbulence suspension flows. By increasing the wall permeability, the turbulence activity is greatly enhanced throughout the channel. Clearly, as the wall permeability increases, the streamwise velocity intensity increases near the interface, which suggests that the slip velocity increases by increasing the wall-normal velocity at the interface.
Figure 2. Contours of the instantaneous turbulent streamwise velocity 
$u$ in the 
$x$–
$y$ plane for a fixed wall non-dimensional permeability parameter 
$\sigma = 0.63 \times 10^{-3}$ and three different particle volume fractions: (a) 
${{\varPhi _b}} = 0.05$, (b) 
${{\varPhi _b}}=0.1$, and (c) 
${{\varPhi _b}}= 0.2$.
Figure 3. Contours of the instantaneous turbulent streamwise velocity 
$u$ in the 
$x$–
$y$ plane for a fixed particle volume fraction 
${{\varPhi _b}} = 0.10$ and three different wall non-dimensional permeability parameters: (a) 
$\sigma = 0.63 \times 10^{-3}$, (b) 
$\sigma = 2 \times 10^{-3}$, and (c) 
$\sigma = 6.3 \times 10^{-3}$.
3.1. Flow statistics
We start our quantitative analysis of the different suspension flows with the normalized pressure drop needed to drive the suspensions at a constant suspension flow rate. Figure 4 shows the friction Reynolds number 
$Re_{\tau }$ for all cases considered in this study as a function of both the wall non-dimensional permeability parameter 
$\sigma$ and the volume fraction 
${{\varPhi _b}}$. The friction Reynolds number is based on the fluid viscosity that can be defined as 
$Re_\tau = u_\tau h / \nu$, where 
$u_\tau = \sqrt {\tau _{interface} / \rho }$ is the friction velocity, and 
$\tau _{interface}$ is the total mean shear stress at the interface. The total mean shear stress at the interface is the sum of the viscous shear stress, the Reynolds shear stress and the particle contribution (see § 3.2 for more details). In the absence of particles, i.e. 
${{\varPhi _b}} = 0$, 
$Re_{\tau }$ increases with the wall permeability, as was also observed in previous studies, e.g. Rosti et al. (Reference Rosti, Cortelezzi and Quadrio2015). Although some works (Hahn et al. Reference Hahn, Je and Choi2002; Itoh et al. Reference Itoh, Tamano, Iguchi, Yokota, Akino, Hino and Kubo2006) showed that a surface with preferential streamwise permeability might have a drag-reducing effect, permeable surfaces have generally been reported to increase the turbulent drag (Jimenez et al. Reference Jimenez, Uhlmann, Pinelli and Kawahara2001; Breugem et al. Reference Breugem, Boersma and Uittenbogaard2006; Rosti et al. Reference Rosti, Cortelezzi and Quadrio2015; Kuwata & Suga Reference Kuwata and Suga2016). Our results are then consistent with previous works, observing a rise in turbulent drag compared to a smooth wall. Note that the skin-friction coefficient of a pure fluid flow on porous walls was first studied by Ruff & Gelhar (Reference Ruff and Gelhar1972), who showed experimentally that the skin-friction coefficient is higher than that of impermeable walls. Later, several researchers (Kong & Schetz Reference Kong and Schetz1982; Zippe & Graf Reference Zippe and Graf1983) measured the increase in the skin-friction coefficient over the porous layers and compared the results with the case of impermeable walls. These authors showed that the presence of a permeable wall reduces the turbulence intensity, Reynolds shear stress, pressure and vorticity fluctuations throughout the channel except very near the wall. The results in figure 4 also reveal that adding particles to the flow leads to a further increase in the friction Reynolds number 
$Re_{\tau }$. It is worth noting that the increase in drag has been associated with the appearance of Kelvin–Helmholtz-like rollers over the surface that increase both momentum transfer and the Reynolds stresses near the interface, as shown in figure 5 and reported in Jimenez et al. (Reference Jimenez, Uhlmann, Pinelli and Kawahara2001), Breugem et al. (Reference Breugem, Boersma and Uittenbogaard2006), Rosti et al. (Reference Rosti, Cortelezzi and Quadrio2015) and Kuwata & Suga (Reference Kuwata and Suga2016). This figure shows snapshots of the fluctuating streamwise velocity at the free-flow porous interface for the cases of a pure solvent and flow with 
${{\varPhi _b}}=0.1$ at fixed 
$\sigma = 0.63 \times 10^{-3}$. The high- and low-velocity regions represent the footprints of the streaky pattern of the flow over the porous layer, shown in figure 5(a), that are similar to those observed over smooth walls. A Kelvin–Helmholtz-like pattern is indicated by these streaky patterns, which are wider and have coherence structures when pure solvent is present. When particles are added to the flow, the coherent structures appear shredded over the porous layer, while the longitudinal streamwise vortices become more twisted. The Kelvin–Helmholtz-like rollers have also been analysed for the solvent flowing over both isotropic and anisotropic porous media in Rosti et al. (Reference Rosti, Brandt and Pinelli2018). These structures can be further detected and analysed from the two-point correlation functions, as has been studied by others (Jimenez et al. Reference Jimenez, Uhlmann, Pinelli and Kawahara2001; Kuwata & Suga Reference Kuwata and Suga2016), which is the subject of our future investigations.
Figure 4. Friction Reynolds number 
$Re_{\tau }$ as a function of (a) wall non-dimensional permeability parameter 
$\sigma$, for different values of the particle volume fraction 
${{\varPhi _b}}$, and (b) the particle volume fraction 
${{\varPhi _b}}$ for different values of wall non-dimensional permeability parameter 
$\sigma$.
Figure 5. Instantaneous streamwise velocity fluctuation contours in the 
$x$–
$z$ plane at the interface 
$y = 0$ (flow going from bottom to top): (a) for solvent 
${{\varPhi _b}}=0.0$, and (b) for suspension with the particle volume fraction 
${{\varPhi _b}}=0.1$ at a fixed wall non-dimensional permeability parameter 
$\sigma = 0.63\times 10^{-3}$.
Figure 6 shows the mean streamwise fluid velocity profiles for different particle volume fractions 
${{\varPhi _b}}$ and the non-dimensional permeability parameter 
$\sigma = 0.63\times 10^{-3}$, the lowest considered in this study. Here, we observe that while the mean fluid velocity is blunted in the pure turbulent flow, it increases at the channel centreline when increasing the volume fraction, leading to a more laminar-like flow. In contrast, near the wall, the fluid velocity decreases by increasing the particle volume fraction 
${{\varPhi _b}}$; this effect is more pronounced for highly viscous suspension flows, as shown in Rosti et al. (Reference Rosti, Mirbod and Brandt2021). Increasing the volume fraction does not modify the velocity inside the porous layer, i.e. for 
$y < 0$. This is because the wall permeability is very low and constant in this case, and the particles are not moving inside the flow either; therefore, any disturbances produced by the particles vanish rapidly inside the porous layer. On the other hand, as expected and can be observed in figure 7 for a fixed particle volume fraction 
${{\varPhi _b}} = 0.10$, by increasing the non-dimensional permeability parameter, the velocity inside the porous layer increases. In this case, for high wall permeability, the fluid moves easily inside the porous layer. Larger values of 
$\sigma$ also result in higher slip velocity at the suspension-porous interface (the velocity at 
$y = 0$). This increase in the slip velocity when increasing the wall non-dimensional permeability parameter from 
$\sigma = 0.63\times 10^{-3}$ to 
$\sigma = 6.3\times 10^{-3}$ (see figure 7) implies that the local shear rate decreases with 
$\sigma$. Nevertheless, it appears that the parameter 
$\sigma$ does not significantly affect the maximum velocity in the centre of the channel and, in general, the velocity profile in the fluid region, 
$0< y/h<1$. Similar observations can be made by examining the data for the other volume fractions.
Figure 6. Mean fluid streamwise velocity component 
$u_f$ normalized by the bulk velocity 
$U_0$ as a function of the wall-normal distance 
$y$ for various particle volume fractions 
${{\varPhi _b}}$ at a fixed wall non-dimensional permeability parameter 
$\sigma = 0.63\times 10^{-3}$.
Figure 7. Mean fluid streamwise velocity component 
$u_f$ normalized by the bulk velocity 
$U_0$ as a function of the wall-normal distance 
$y$ for different values of the wall non-dimensional permeability parameter 
$\sigma$ at a fixed particle volume fraction 
${{\varPhi _b}} = 0.10$.
Both the presence of particles in the flow and the existence of the porous layers alter the velocity fluctuations in the bulk as well as at the wall. The r.m.s. values of the three fluid velocity fluctuations, normalized by the bulk velocity, are displayed in figure 8 for the wall non-dimensional permeability parameter 
$\sigma = 0.63\times 10^{-3}$ and the different values of the particle volume fraction 
${{\varPhi _b}}$. We note that a reduction in the fluid streamwise velocity fluctuations occurs with increasing volume fraction to 
${{\varPhi _b}} = 0.20$. On the other hand, the intensity of the cross-stream velocity fluctuations increases by adding particles in suspensions up to 
${{\varPhi _b}}=0.1$, indicating that the presence of particles creates a more isotropic flow, and then decreases at the largest 
${{\varPhi _b}}$ considered, in agreement with the results for the flow over a rigid wall reported by Lashgari et al. (Reference Lashgari, Picano, Breugem and Brandt2014) and Picano et al. (Reference Picano, Breugem and Brandt2015) for smaller particles. Moreover, in the presence of particles, the Reynolds shear stress increases consistently up to 
${{\varPhi _b}}=0.10$, as shown in figure 8(d). Upon further increasing the particle volume fraction to 
${{\varPhi _b}} = 0.20$, the Reynolds shear stress decreases, confirming a reduction in turbulent activity as shown in Picano et al. (Reference Picano, Breugem and Brandt2015) for the case of a solid wall. This reduction, however, is not associated with a decreased drag because the larger particle stresses at higher 
${{\varPhi _b}}$ more than compensate for the lower turbulence levels. The results pertaining to the flow over a porous wall confirm this picture, with some differences, as will be discussed further below.
Figure 8. Intensity of the fluid fluctuation velocity components normalized by 
$U_0$ and the Reynolds shear stress for different particle volume fractions and non-dimensional permeability parameter 
$\sigma = 0.63\times 10^{-3}$: (a) streamwise 
$u_{rms}$ velocity fluctuations, (b) wall-normal 
$v_{rms}$ velocity fluctuations, (c) spanwise 
$w_{rms}$ velocity fluctuations, and (d) shear stress 
$(uv)'$.
The variation of the r.m.s. values of the fluid velocity fluctuations with the wall non-dimensional permeability parameter 
$\sigma$ is depicted in figure 9 for the flows at particle volume fraction 
${{\varPhi _b}} = 0.10$. The magnitude of the r.m.s. values of all the fluid velocity fluctuations increases with the wall permeability, consistent with previous studies for single-phase flows (Breugem et al. Reference Breugem, Boersma and Uittenbogaard2006; Rosti et al. Reference Rosti, Cortelezzi and Quadrio2015) where the increase in the streamwise velocity component is larger than that of the wall-normal and spanwise velocity fluctuations. In the boundary layer above the permeable wall, a sharp increase in the r.m.s. velocities is observed when compared to the layer below the permeable wall, as also seen for the spanwise and wall-normal wall fluctuations. This can be attributed to the strengthening of the viscous effects near a permeable wall (Perot & Moin Reference Perot and Moin1995). Note that due to the fact that the r.m.s. values of true local velocities inside the permeable walls are negligible, we can assume that the turbulent flow in the channel is not affected significantly either by the thickness of the porous layer or by the distance to the solid wall.
Figure 9. Intensity of the fluid fluctuation velocity components normalized by 
$U_0$ and the Reynolds shear stress for 
${{\varPhi _b}} = 0.10$ and different values of the wall non-dimensional permeability parameter: (a) streamwise 
$u_{rms}$ velocity fluctuations, (b) wall-normal 
$v_{rms}$ velocity fluctuations, (c) spanwise 
$w_{rms}$ velocity fluctuations, and (d) shear stress 
$(uv)'$.
The results suggest that the particle dynamics is indeed modified by the wall permeability and, in turn, the particle motion modifies the flow. Figure 9(d) shows the profiles of the Reynolds stresses for different values of the wall non-dimensional permeability parameter 
$\sigma$, normalized by the bulk velocity, 
$U_0^{2}$. The maximum shear stress occurs near the interface with the permeable wall. The data also show that the Reynolds shear stress increases with the wall permeability. This confirms that for the higher wall permeability, the weakening of the wall-blocking effect causes a strong increase in the Reynolds shear stress, as reported by Breugem et al. (Reference Breugem, Boersma and Uittenbogaard2006). In other words, fluid with high streamwise momentum may penetrate inside the porous surface, while to satisfy the conservation of mass, the fluid with relatively low streamwise momentum is transported from the porous region into the channel.
The slip velocity 
$u_s$ at the interface 
$y=0$, normalized by the bulk fluid velocity, is shown in figure 10(a) for different values of the wall non-dimensional permeability parameter and the particle volume fraction. As expected, by increasing the non-dimensional permeability parameter 
$\sigma$, the slip velocity increases, while with increasing particle volume fraction, it decreases for all cases studied here. Interestingly, for 
${{\varPhi _b}} \le 0.10$, the slip velocity is weakly dependent on the volume fraction, and it is determined mainly by the permeability under dilute conditions. Indeed, a more significant variation in the 
$u_s$ values is observed when further increasing the dispersed phase concentration to 0.2. These results confirm that the suspension slip velocity in turbulent flows depends on the porous permeability, and when the porous non-dimensional permeability parameter is low, i.e. of the order of 
$\sigma ={O}(10^{-4})$, the slip velocity depends mainly on the suspension concentration.
Figure 10. (a) Normalized interface slip velocity as a function of the wall non-dimensional permeability parameter 
$\sigma$ for the different particle volume fractions under investigation. (b) Normalized interface wall-normal fluctuation velocity as a function of the wall non-dimensional permeability parameter 
$\sigma$ for different particle volume fractions.
We then characterized the normalized wall-normal velocity fluctuations at the interface for different values of the wall non-dimensional permeability parameter 
$\sigma$ and the particle volume fraction 
${{\varPhi _b}}$, as shown in figure 10(b). The wall-normal velocity fluctuations increase with both the particle volume fraction and the wall permeability, from zero for small wall permeability to the largest value for the highest 
${{\varPhi _b}}$. Note that the velocity fluctuations at the suspension-porous interface do not change significantly with the values of 
$\sigma$ investigated; instead, the effect of the particle volume fraction is more pronounced with a more significant increase for 
${{\varPhi _b}}=0.2$. Comparing figures 10(a,b), one can conclude that the mean slip at the wall is determined mainly by the wall permeability, while the velocity fluctuations are significantly influenced by the particle dynamics in the free-flow region.
To understand further how turbulent structures are modified by the presence of both porous media and particles, we examine some basic features of turbulent coherent structures. As shown in Rosti & Brandt (Reference Rosti and Brandt2020), the turbulent fluctuations are generally reduced and the flow is very anisotropic in the case of drag reduction. On the other hand, the presence of particles causes the opposite effect, i.e. more isotropic structures. Here, in particular, we study the effect of wall permeability on the distribution of the wall-normal and streamwise fluctuation velocity. Figure 11 shows the joint probability distribution functions (p.d.f.s) of the wall-normal and streamwise fluctuation velocity at the interface of the porous layer and suspensions between the centre of the channel and the interface, and also at the centreline of the channel for the flows at particle volume fraction 
${{\varPhi _b}} = 0.10$ and the three wall non-dimensional permeability parameters under consideration. The ejection and sweep events are disrupted by increasing the wall permeability, leading to a decrease in the streamwise velocity fluctuations 
$u'$ and an increase of the wall-normal component 
$v'$. These results show that permeability creates a more isotropic flow; this has been associated with drag increase, as reported in several previous studies, e.g. Breugem et al. (Reference Breugem, Boersma and Uittenbogaard2006). Note also that the differences are minor in the channel centreline. Similarly, figure 12 shows the joint p.d.f.s of the wall-normal and streamwise fluctuation velocity at the suspension flow-porous interface, between the centre of the channel and the interface, and at the centreline of the channel, for the flows at wall non-dimensional permeability parameter 
$\sigma = 0.63 \times 10^{-3}$ for different particle volume fractions. As observed, the presence of particles induces a disruption of the structures in the near-wall region, resulting in a decrease in 
$u'$ and an increase in 
$v'$, leading to an overall more isotropic state. In other words, both wall permeability and the presence of dispersed particles induce a more isotropic flow. However, as discussed above, the increase in wall permeability also leads to an increase in the turbulent Reynolds stresses, while the more isotropic turbulence due to the presence of the particles is associated with less correlated motions, i.e. lower Reynolds stresses (cf. figures 8d and 9d).
Figure 11. Normalized probability distribution functions (p.d.f.s) of the streamwise fluctuation velocity 
$u'$ and wall-normal fluctuation velocity 
$v'$ for a fixed particle volume fraction 
${{\varPhi _b}} = 0.10$ and various non-dimensional permeability parameters: (a–c) 
$\sigma = 0.63\times 10^{-3}$, (d–f) 
$\sigma = 2.0\times 10^{-3}$, and (g–i) 
$\sigma = 6.3\times 10^{-3}$. (a,d,g) At the flow-porous interface, (b,e,h) between the centre of the channel and the interface, and (c, f,i) at the centreline of the channel.
Figure 12. Normalized probability distribution functions (p.d.f.s) of the streamwise fluctuation velocity 
$u'$ and wall-normal fluctuation velocity 
$v'$ at the interface for a fixed wall non-dimensional permeability parameter 
$\sigma = 0.63\times 10^{-3}$ and various particle volume fractions: (a–c) 
${{\varPhi _b}} = 0.05$, (d–f) 
${{\varPhi _b}} = 0.10$, and (g–i) 
${{\varPhi _b}} = 0.20$. (a,d,g) At the flow-porous interface, (b,e,h) between the centre of the channel and the interface, and (c, f,i) at the centreline of the channel.
3.2. Particle dynamics and stress budget
In this subsection, we investigate the particle dynamics in turbulent channel flow in the presence of porous walls. In the case of spherical particles, the formation of a particle wall layer and the rotation rate close to the wall are shown to have a significant impact on the turbulent global drag (Costa et al. Reference Costa, Picano, Brandt and Breugem2016, Reference Costa, Picano, Brandt and Breugem2018; Peng et al. Reference Peng, Ayala and Wang2019). Indeed, in the presence of an elastic wall, particles are shown to migrate away from the wall, which significantly reduces the overall pressure drop (Ardekani, Rosti & Brandt Reference Ardekani, Rosti and Brandt2019). Figure 13(a) shows the mean local volume fraction 
$\phi$ for the case of wall non-dimensional permeability parameter 
$\sigma = 0.63\times 10^{-3}$ and the three nominal volume fractions 
${{\varPhi _b}}$ investigated, while figures 13(b,c) display how the particle distribution changes with the permeability for two different values of the particle volume fraction. Note that the oscillations in figures 13(b,c) can be attributed to particle layering close to the interface, and as the particle diameter is 0.2, these oscillations occur at one, two or three times the particle diameter.
Figure 13. Average particle concentration 
$\phi$ as a function of the wall-normal distance 
$y$ for: (a) different particle volume fractions at fixed non-dimensional permeability parameter 
$\sigma = 0.63\times 10^{-3}$; (b,c) various wall non-dimensional permeability parameters 
$\sigma$ at fixed particle volume fractions 
${{\varPhi _b}} = 0.10$ and 
$0.20$.
Overall, we note the formation of a first local maximum close to the interface, associated with the particle-wall layers in the free-flow region. As also observed in highly dense suspensions in laminar flows (Yeo & Maxey Reference Yeo and Maxey2010), the intensity of the particle wall layer increases as the particle volume fraction increases (see figure 13a). Its formation is explained by the solid lubrication interactions when a particle reaches the wall, which stabilizes the near-wall particle position; once close to the wall, the motion of finite-size particles is influenced mainly by collisions with other particles. Consequently, it is difficult for particles belonging to the first layer to escape.
In the presence of a permeable wall, we note weak migrations of particles away from the wall and a less pronounced wall layer for 
${{\varPhi _b}}=0.05$ and 0.10, with the peak value close to the wall found to decrease with the permeability, consistent with the observation of increased wall-normal fluctuations (Ardekani et al. Reference Ardekani, Rosti and Brandt2019, for the case of suspension flows over elastic walls, see also). Note that these near-wall maxima are of the order of the bulk concentrations and are not related to the turbophoretic drift observed in dilute suspensions of heavy particles (Reeks Reference Reeks1983). The same trend observed here was reported in Rosti et al. (Reference Rosti, Mirbod and Brandt2021) for laminar Couette flow of suspensions over a porous substrate. Finally, we also see an increase in the local volume fraction 
$\phi$ towards the channel centreline for the largest volume fraction examined, 
${{\varPhi _b}}=0.20$. This migration toward the channel core was observed experimentally in turbulent flows over solid walls and related to inertial effects in Zade et al. (Reference Zade, Costa, Fornari, Lundell and Brandt2018). The accumulation towards the centreline was also associated with quenching of the turbulence fluctuations and reduced mixing (Ardekani et al. Reference Ardekani, Al Asmar, Picano and Brandt2018; Yousefi et al. Reference Yousefi, Ardekani, Picano and Brandt2021).
To examine further the global suspension behaviour, we look at the streamwise momentum budget. The total momentum budget in the free-flow region (i.e. 
$0 < y < 2h$) can be written as the sum of the viscous stress, the Reynolds stress, and the particle contribution (Picano et al. Reference Picano, Breugem and Brandt2015, for more details and derivations, see) as
\begin{equation} \tau = \tau^{visc} + \tau^{reyn} + \tau^{part}. \end{equation}
Here, 
$\tau$ is the mean total shear stress balancing the imposed pressure gradient, 
$\tau ^{visc} = \mu ({{\rm d}u}/{{\rm d} y})$ is the viscous shear stress, 
$\tau ^{reyn} = -\rho (uv)'$ is the Reynolds shear stress associated with the fluid and particle velocity correlation, and 
$\tau ^{part}$ is the particle contribution, including particle collisions, which was found by subtracting the other two components from the total shear stress.
The different contributions to the total shear stress in the turbulent suspension flows are displayed for the different particle volume fractions and wall non-dimensional permeability parameter 
$\sigma = 0.63\times 10^{-3}$ in figure 14(a), and for 
$\sigma = 6.3\times 10^{-3}$ in figure 14(b). Overall, we note that the relative viscous stress varies weakly with both the volume fraction and the wall permeability, whereas the percentage contribution of turbulence decreases with 
${{\varPhi _b}}$. Although the level of fluctuations increases up to 
${{\varPhi _b}}=0.1$, the turbulent stress decreases because these chaotic motions are less correlated. The contribution associated with the particle stresses increases with 
${{\varPhi _b}}$, differently for different permeabilities. In particular, we see a significant increase when increasing from 0.05 to 0.1 at the lowest permeability, and the most significant growth at 
${{\varPhi _b}}=0.2$ at the highest permeability studied. This suggests that the weakening of the particle wall layer is responsible for the reduced turbulent transport at low values of 
${{\varPhi _b}}$, which is more significant at lower values of the wall permeability. On the other hand, the particle accumulation towards the channel core by inertial collisions is responsible for the reduction of the turbulence activity at 
${{\varPhi _b}}=0.2$, an effect more weakly dependent on the wall permeability; cf. figures 14(a,b). To summarize, the main contribution to the total shear stress is the particle stress at high concentrations, while the Reynolds stress contribution is higher for lower concentrations. The viscous shear stress remains constant for all the volume fractions and the different wall permeabilities considered here. The values of viscous stress are small compared to the two other components because of the turbulent flow considered here, despite the moderate Reynolds number. Focusing on the effect of permeability, the particle-stress contribution decreases with increasing wall non-dimensional permeability parameter 
$\sigma$, whereas the Reynolds stress shows the opposite behaviour, i.e. it decreases with the volume fraction 
${{\varPhi _b}}$ and increases with the wall non-dimensional permeability parameter 
$\sigma$ because of the increase in fluctuations at the interface.
Figure 14. Histograms showing the different components of the mean shear stress balance as a function of the particle volume fraction 
${{\varPhi _b}}$ for (a) 
$\sigma = 0.63\times 10^{-3}$, and (b) 
$\sigma = 6.3\times 10^{-3}$. Blue, red and yellow show the viscous stress, Reynolds shear stress, and particle contributions, respectively.
Figure 15 displays the spanwise component of the mean particle angular velocity for different particle volume fractions as a function of the wall permeability, where the values are averages over all the particles close to the interface, i.e. 
$y/h=0.05$. It can be revealed that the particle-averaged angular velocity decreases by increasing the wall non-dimensional permeability parameter 
$\sigma$. As expected, this is because the near-wall shear decreases with permeability, as discussed above. In addition, we note that the spanwise component of the mean particle angular velocity decreases with the particle volume fraction 
${{\varPhi _b}}$, as observed in previous studies (Ardekani & Brandt Reference Ardekani and Brandt2019). However, this effect is small in comparison with the variations induced by the wall permeability. In brief, with increasing wall permeability, the particle rotation decays; consequently, this causes a reduction in particle-induced interactions. The inset of figure 15 plots the ratio of particle angular velocity with respect to the fluid shear rate as a function of bulk particle volume fraction 
${{\varPhi _b}}$ for various permeability parameter 
$\sigma$. Interestingly, the data in the main plot of figure 15 for both 
${{\varPhi _b}} = 0.05$ and 
$0.1$ are now collapsed onto a curve 
$\omega _z/[{({\rm d}u/{{\rm d} y})}_{interface}] = 0.03346 - 4.419\sigma + 313.3{\sigma }^2$ that has been determined from a polynomial data fit. There is, however, a higher ratio for 
${{\varPhi _b}}=0.2$. At 
${{\varPhi _b}}=0.2$, the fluid shear rate is lower, as shown in figure 6, and the induced particle velocity and angular velocity are both affected by particles gathered at the interface.
Figure 15. Spanwise component of the near-wall particle angular velocity as a function of the wall non-dimensional permeability parameter 
$\sigma$ for different volume fractions 
${{\varPhi _b}} = 0.05$, 
$0.10$ and 
$0.20$.
Figure 16(a) shows the mean streamwise particle velocity profiles for different volume fractions 
${{\varPhi _b}}$ at a fixed wall non-dimensional permeability parameter 
$\sigma = 0.63\times 10^{-3}$, while figure 16(b) displays the particle velocity profile for different wall non-dimensional permeability parameters and a fixed particle volume fraction 
${{\varPhi _b}} = 0.10$. Comparing the mean particle velocity profiles (lines) with the mean fluid velocity values (symbols) indicates that the mean particle velocity and mean fluid velocity profiles are comparable inside the free-flow region, except near the porous wall where the mean particle velocity is larger than the mean local fluid velocity. Although the fluid velocity at the interface is almost zero, the particle velocity differs due to the lack of the no-slip condition, resulting in a greater particle velocity than the fluid velocity.
Figure 16. Mean particle streamwise velocity component 
$u_p$ normalized by the bulk velocity 
$U_0$ as a function of the wall-normal coordinate 
$y$ (lines): (a) for various particle volume fractions 
${{\varPhi _b}}$ at a fixed wall non-dimensional permeability parameter 
$\sigma = 0.63\times 10^{-3}$; and (b) for various wall non-dimensional permeability parameters 
$\sigma$ at a fixed particle volume fraction 
${{\varPhi _b}} = 0.10$. The symbols represent the mean fluid velocity profiles.
3.3. Effect of the momentum transfer coefficient 
$\tau$
We now investigate the modification induced on the turbulence statistics by varying the coefficient 
$\tau$, which is used for the momentum transfer conditions at the flow-porous interface. This coefficient can be related to the unknown structure of the suspending fluid-porous interface. As mentioned above, the value of 
$\tau$ determines the stress transfer between the suspending fluid and the porous media. This was first introduced by Ochoa-Tapia & Whitaker (Reference Ochoa-Tapia and Whitaker1995) as an additional boundary condition at the interface; it is order one and can be positive or negative. It was characterized further by Ochoa-Tapia & Whitaker (Reference Ochoa-Tapia and Whitaker1995), Valdés-Parada et al. (Reference Valdés-Parada, Goyeau and Ochoa-Tapia2007, Reference Valdés-Parada, Aguilar-Madera, Ochoa-Tapia and Goyeau2013), Minale (Reference Minale2014a,Reference Minaleb) and Chen, Wang & Yang (Reference Chen, Wang and Yang2016). In particular, Minale (Reference Minale2014a,Reference Minaleb) showed that a negative 
$\tau$ quantifies the amount of stress transferred from the free fluid to the porous matrix, while a positive 
$\tau$ quantifies the amount of stress transferred from the porous matrix to the free fluid, whereas 
$\tau = 0$ represents the case in which the stress carried by the free fluid is fully transferred to the fluid saturating the porous matrix. Here, we consider the flow of a pure fluid, 
${{\varPhi _b}} = 0$, and a dense suspension, 
${{\varPhi _b}} = 0.20$, for comparison; we expect that the turbulence statistics for other volume fractions fall within these two cases. The porous layer non-dimensional permeability parameter is fixed at 
$\sigma = 0.63\times 10^{-3}$, and the porosity is 
$\varepsilon = 0.6$, again with Reynolds number 
$Re = 5600$. We compare data from three values of 
$\tau$, i.e. 
$-1$, 
$0$ and 
$1$, chosen because the momentum transfer coefficient 
$\tau$ falls in this range for most of the porous media. Note that the case 
$\tau = 0$, as reported in the previous subsections, is used as the reference case.
Figure 17(a) shows the friction Reynolds number 
$Re_{\tau }$, where we observe a significant drag increase for a non-zero momentum transfer coefficient 
$\tau$. Note that 
$Re_{\tau }$ is calculated using the mean pressure gradient needed to keep 
$U_0$ constant. Decreasing the slip velocity and particle rotation causes a decrease in the particle concentration near the interface, as shown in figure 17(b), where we note that the particle wall layer is significantly less pronounced when there is a net stress transferred from the porous matrix to the free fluid (
$\tau \ne 0$). The data in the figure also show that for non-zero values of 
$\tau$, the local maximum of the particle distribution near the channel centreline increases considerably.
Figure 17. (a) Friction Reynolds number 
$Re_{\tau }$ as a function of 
$\tau$, for different particle volume fractions 
${{\varPhi _b}}=0$ and 0.2. (b) Average particle concentration 
$\phi$ as a function of the wall-normal distance 
$y$ for different 
$\tau$ at a fixed non-dimensional permeability parameter 
$\sigma = 0.63 \times 10^{-3}$.
Figure 18 shows the mean velocity profiles for the three momentum transfer coefficients considered. Figure 18(a) shows the velocity profile pertaining to the single-phase flow over porous media, 
${{\varPhi _b}} = 0$, whereas figure 18(b) reports the data for 
${{\varPhi _b}} = 0.20$. The dashed line represents the case with positive 
$\tau$, the dash-dot-dotted line shows the case with negative 
$\tau$, and the solid line shows our reference case where the momentum transfer coefficient is 
$\tau = 0$. The effects of the momentum transfer coefficient 
$\tau$ on the mean velocity profile are clear, especially at the suspension-porous interface: a non-zero value of 
$\tau$ induces a significant decrease in the velocity at the suspension-porous interface for both the highly dense case (
${{\varPhi _b}} =0.20$) and the single-phase case (
${{\varPhi _b}} =0$), while a different behaviour is observed at the centreline for the two cases. For the pure fluid 
${{\varPhi _b}} =0$, the non-zero 
$\tau$ induces an increase in the velocity at the centreline. However, for the highly-dense case, a velocity decrease can be observed at the centreline for the non-zero values of the momentum transfer coefficient. The results for the pure fluid flow over the porous layer are in agreement with previous works (Tilton & Cortelezzi Reference Tilton and Cortelezzi2008; Rosti et al. Reference Rosti, Cortelezzi and Quadrio2015).
Figure 18. Mean fluid streamwise velocity component 
$u_f$ normalized by the bulk velocity 
$U_0$ as a function of the wall-normal distance 
$y$ for different values of momentum transfer coefficient 
$\tau$ and wall non-dimensional permeability parameter 
$\sigma = 0.63 \times 10^{-3}$: (a) pure fluid (
${{\varPhi _b}} = 0$), (b) 
${{\varPhi _b}} = 0.20$.
The normalized slip velocities for different cases studied here are shown in figure 19(a). This figure shows that the slip velocities for non-zero momentum transfer coefficient 
$\tau$ are much smaller than for the case 
$\tau =0$ in both single-phase flow, i.e. 
${{\varPhi _b}} = 0$, and 
${{\varPhi _b}} = 0.20$. Note also that even for larger permeabilities, previous works observed the same trend for the slip velocity for both positive and negative 
$\tau$; see, for example, the results for a single-phase flow over porous media in Rosti et al. (Reference Rosti, Cortelezzi and Quadrio2015). Figure 19(b) displays the normalized spanwise component of the mean particle angular velocity for the different momentum transfer coefficients 
$\tau$ under consideration. Note that the values have been averaged over all particles close to the interface. These data reveal that the particle average angular velocity decreases significantly for values of the momentum transfer coefficient 
$\tau$ different from zero, consistent with the trend exhibited by the slip velocity.
Figure 19. (a) Normalized interface slip velocity as a function of 
$\tau$ for different particle volume fractions under investigation at wall non-dimensional permeability parameter 
$\sigma = 0.63 \times 10^{-3}$. (b) Normalized spanwise component of the particle angular velocity as a function of 
$\tau$ for particle volume fractions 
${{\varPhi _b}} = 0.20$.
The variations in the r.m.s. values of the fluid velocity fluctuations with the momentum transfer coefficient 
$\tau$ are depicted in figure 20 for the flows at particle volume fraction 
${{\varPhi _b}} = 0.20$. The magnitude of the r.m.s. values of the streamwise fluid velocity fluctuations decreases for non-zero 
$\tau$. The same can be seen for the wall-normal and spanwise velocity fluctuations, although the reduction appears to be less pronounced and more evident outside the viscous wall layer. The results therefore confirm the reduction in turbulent transport observed from the momentum budget shown in figure 21.
Figure 20. Intensity of the fluid fluctuation velocity components normalized by 
$U_0$ and the Reynolds shear stress for different 
$\tau$ and a fixed particle volume fraction 
${{\varPhi _b}} = 0.20$ and non-dimensional permeability parameter 
$\sigma = 0.63\times 10^{-3}$: (a) streamwise 
$u_{rms}$ velocity fluctuation; (b) wall-normal 
$v_{rms}$ velocity fluctuation; (c) spanwise 
$w_{rms}$ velocity fluctuation; and (d) shear stress 
$(uv)'$.
Figure 21. Histograms showing the different components of the mean shear stress balance as a function of 
$\tau$ for (a) 
${{\varPhi _b}} = 0.0$ and (b) 
${{\varPhi _b}} = 0.20$. Blue, red and yellow show the viscous stress, Reynolds shear stress and particle contributions, respectively.
To better understand the modifications on the dynamics and on the overall drag induced by a different stress condition at the wall, we also performed the momentum budget analysis; see (3.1) and relative discussion. The results of this analysis are displayed in figure 21. For the pure fluid, 
${{\varPhi _b}} = 0$, by changing the momentum transfer coefficient from zero to 
${-1}$ and 
$1$, 
$\tau ^{visc}$ increases and 
$\tau ^{reyn}$ decreases. As a result, despite the attenuation of the turbulence, the overall drag increases due to the increase in the viscous stresses. For the case of a dense suspension, 
${{\varPhi _b}} = 0.20$, we see that the attenuation of the turbulence is relatively more important, which can be explained by the increased particle migration towards the channel core region (see figure 17b). However, the formation of a relatively densely packed region causes an increase in the particle stresses that more than balances the decrease in turbulent transport. As a consequence, the drag increase at 
${{\varPhi _b}} = 0.20$ is more pronounced for 
$\tau = \pm 1$. These results indicate that the condition at the suspension-porous interface has a significant impact on the flow and particle dynamics for the permeability of the porous medium considered here.
4. Conclusion
We employed direct numerical simulations to study turbulent suspension flows in a channel covered with isotropic porous media on the top and bottom surfaces. The VANS equations are used to describe the volume-averaged flow inside the porous layer, with the permeable wall characterized by permeability and porosity. We then coupled VANS with the IBM to resolve particle–particle and particle–fluid interactions. Here, we have kept the Reynolds number 
$Re = {U_0h}/{\nu } = 5600$ and investigated the effect of particle volume fraction 
${{\varPhi _b}}$, non-dimensional permeability parameter 
$\sigma$, and the momentum transfer coefficient 
$\tau$ of the porous layer for spherical and neutrally buoyant rigid particles.
Using the Reynolds turbulent friction number 
$Re_{\tau }$, i.e. the response to determining the turbulent drag based on the friction velocity, we show quantitatively that turbulent stress increases with both particle volume fraction 
${{\varPhi _b}}$ and wall non-dimensional permeability parameter 
$\sigma$. Overall, we observed that wall permeability is associated with higher drag due to increased turbulence activity compared to pure solvent. Turbulent suspension flows over hyperelastic materials have also been reported in Rosti & Brandt (Reference Rosti and Brandt2017). Although some experiments for pure fluid flowing over porous surfaces show that wall permeability causes an increase in skin friction in the turbulent-flow regime (Zagni & Smith Reference Zagni and Smith1976; Zippe & Graf Reference Zippe and Graf1983), others reveal a decrease in the skin friction in the laminar-flow regime (e.g. Beavers & Joseph Reference Beavers and Joseph1967).
Because the channel top and bottom walls have an isotropic porous medium, the velocity profile in the free-flow region remains the same as if the walls were smooth, as reported in Picano et al. (Reference Picano, Breugem and Brandt2015), but a slip velocity 
$u_s$ is generated at the suspension-porous interface. In our study, we found that the slip velocity 
$u_s$ decreases with increasing particle volume fraction 
${{\varPhi _b}}$. By contrast, when the permeability of the wall increases, the slip velocity increases. In general, the decrease in particle–particle interaction results from increasing the slip velocity 
$u_s$ with the wall permeability parameter 
$\sigma$.
Our analysis of the flow streamwise momentum balance showed that as particle volume fraction increases, the additional stress due to particles increases. For instance, at 
${{\varPhi _b}}=0.20$, the particle stress dominates the flow, while Reynolds transport dominates at 
${{\varPhi _b}}=0$ and 
${{\varPhi _b}} < 0.2$. It is also shown that adding more particles increases particle-induced stresses, which causes drag to increase. While turbulent stresses decay rapidly with increasing particle volume fractions, viscous stresses remain nearly constant.
A final analysis was performed on the statistics of turbulence caused by momentum transfer coefficients 
$\tau$. Slip velocity decreases noticeably when momentum transfer coefficient 
$\tau$ is changed. Based on our studies, we found that for a non-zero momentum transfer coefficient, the average particle angular velocity near the wall is lower than 
$\tau = 0$. The other noticeable effect of the momentum transfer coefficient 
$\tau$ on the concentration profile is where non-zero 
$\tau$ leads to the migration of particles to the centre of the channel, which means that fewer particles gather near the free flow and porous interface compared to 
$\tau = 0$. Increasing particle contribution also decreases Reynolds shear stress, thereby reducing total shear stress. The results also show that 
$\tau$ has a considerable effect on the behaviour of turbulent suspension flow over porous media; therefore, we need to consider the momentum transfer coefficient as an essential factor in future studies.
To improve this study further, other important factors can be considered. We should certainly pay attention to the effect of particle size (Lashgari et al. Reference Lashgari, Picano, Costa, Breugem and Brandt2017), shape (Bellani et al. Reference Bellani, Byron, Collignon, Meyer and Variano2012) and deformability (Alghalibi, Rosti & Brandt Reference Alghalibi, Rosti and Brandt2019) as new insights to our current study. In addition, experiments are needed to verify our results. We assumed that particles in this study would not be able to enter the porous layer due to their size. Understanding whether the same behaviour occurs when particles move inside the porous media, and the impact of anisotropic porous media, will be the subjects of our future investigations.
Funding
P.M. was supported by the National Science Foundation (award no. 1854376) and by the Army Research Office (award no. W911NF-18-1-0356). The work of S.M.A. and A.M.B. was funded by the National Science Foundation (award no. 1854376). L.B. acknowledges financial support from the Swedish Research Council (VR) through the INTERFACE research environment (grant no. VR 2016-06119). MER was supported by the Okinawa Institute of Science and Technology Graduate University (OIST) with subsidy funding from the Cabinet Office, Government of Japan.
Declaration of interests
The authors report no conflict of interest.
