%20-%20https://sketch.com%20--%3e%3ctitle%3e3ECBBF50-FD08-412A-943D-2CFC5E025A16%3c/title%3e%3cdesc%3eCreated%20with%20sketchtool.%3c/desc%3e%3cg%20id='Stylesheet'%20stroke='none'%20stroke-width='1'%20fill='none'%20fill-rule='evenodd'%3e%3cg%20transform='translate(-926.000000,%20-668.000000)'%20id='BUTTONS'%3e%3cg%20transform='translate(100.000000,%20597.000000)'%3e%3cg%20id='BTN-/-RIGHTS'%20transform='translate(816.000000,%2065.000000)'%3e%3cg%20id='ICON-/-Rights'%20transform='translate(10.000000,%206.000000)'%3e%3cpath%20d='M13.5752841,0%20L1.34943182,0%20C0.59375,0%200,0.586666667%200,1.33333333%20L0,22.6666667%20C0,23.4133333%200.620738636,24%201.34943182,24%20L17.6505682,24%20C18.40625,24%2019,23.3866667%2019,22.6666667%20L19,5.36%20L13.5752841,0%20Z%20M13.5752841,1.89333333%20L17.0838068,5.36%20L13.5752841,5.36%20L13.5752841,1.89333333%20Z%20M1.34943182,22.6666667%20L1.34943182,1.33333333%20L12.1988636,1.33333333%20L12.1988636,5.36%20C12.1988636,6.10666667%2012.7926136,6.69333333%2013.5482955,6.69333333%20L17.6235795,6.69333333%20L17.6505682,22.6666667%20L1.34943182,22.6666667%20Z'%20id='Shape'%20fill='%23595959'%20fill-rule='nonzero'%3e%3c/path%3e%3ccircle%20id='Oval'%20stroke='%23595959'%20stroke-width='1.5'%20fill='%23595959'%20cx='9.5'%20cy='13.5355339'%20r='5.5'%3e%3c/circle%3e%3cpolyline%20id='Path-2'%20stroke='%23FFFFFF'%20stroke-width='1.5'%20transform='translate(9.500000,%2013.035534)%20rotate(45.000000)%20translate(-9.500000,%20-13.035534)%20'%20points='8%2016.0355339%2011%2016.0355339%2011%2010.0355339'%3e%3c/polyline%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/svg%3e)
1. Introduction
Dispersed two-phase flows are present in many industrial and geophysical applications such as fluidized beds, slurry flows or sediment transport. Our ability to predict the dynamics of the system as a whole relies on our understanding of the fine-scale physical processes such as particle–particle interactions or fluid–particle interactions. One of the key challenges is the coupling between the particles and the carrier phase turbulence, the so-called turbulence-particle interactions (Balachandar & Eaton Reference Balachandar and Eaton2010). The modelling methodology has to be carefully chosen depending on the available computational resources, flow regime and turbulence-particle interaction regime.
For particles having a response time 
$t_s$ smaller than the Kolmogorov time scale 
$t_\eta$ associated with the smallest turbulent scales 
$\eta$ (typically 
$St<0.2$ with 
$St=t_s/t_\eta$ the Stokes number), the particles will follow almost exactly the carrier phase turbulence at all scales. For this regime, the equilibrium-Eulerian approach is a good approximation to model the particles dynamics and only mass and momentum conservation equations for the carrier phase are solved together with a relaxation equation for the particle phase velocity and the particle phase mass conservation equation (Ferry & Balachandar Reference Ferry and Balachandar2005). In many geophysical or industrial flows, the Stokes number may exceed 0.2 (
$St>0.2$) and the particles no longer follow the carrier phase turbulence exactly (Balachandar & Eaton Reference Balachandar and Eaton2010). In this situation, more sophisticated methodologies such as fully resolved direct numerical simulation (DNS), Eulerian–Lagrangian point-particle models or Eulerian–Eulerian two-fluid models are required to take into account the couplings between the particles and the carrier phase turbulence (two-way coupling) and the particle-particle interactions (four-way coupling).
The most accurate method to account for turbulence-particle interactions is the fully resolved DNS (e.g. Kidanemariam et al. Reference Kidanemariam, Chan-Braun, Doychev and Uhlmann2013; Vowinckel, Kempe & Fröhlich Reference Vowinckel, Kempe and Fröhlich2014; Vowinckel et al. Reference Vowinckel, Nikora, Kempe and Fröhlich2017). The interface between the carrier phase and the particles and, by extension, the fluid–particle interactions are explicitly resolved. In order to use this method, two constraints on the grid need to be satisfied: (i) the grid size needs to be everywhere of the order of the Kolmogorov length scale 
$\eta$, and (ii) the grid size should not be larger than one tenth of the particle size 
$d_p$ (
$\varDelta \sim d_p/10$, with 
$d_p$ the particle diameter). Putting together these constraints, fully resolved particle-laden boundary layer flow DNS is only achievable for bulk Reynolds numbers up to approximately 
$10^{3}$, with at most a few million particles and a billion grid points. In order to achieve Reynolds numbers relevant to realistic sediment transport conditions 
$O(10^{5})$ for medium to very coarse sand, simulations would require of the order of 
$10^{12}$ to 
$10^{14}$ grid points. Such simulations are not possible with today's computational resources, and a compromise has to be found in terms of modelling strategy.
Concerning the Lagrangian point-particle approach, particles are considered punctual, their interactions with the carrier phase are modelled and Newton's second law is used to predict their trajectories (Maxey & Riley Reference Maxey and Riley1983). The limitations are twofold, on the one hand, the computational grid size 
$\varDelta$ has to be much greater than the particle size and, on the other hand, the domain size is limited by the maximum number of particles achievable in the simulation. As an example, for the two-fluid simulation of scour around cylinders at the laboratory scale by Mathieu et al. (Reference Mathieu, Chauchat, Bonamy and Nagel2019) and Nagel et al. (Reference Nagel, Chauchat, Bonamy, Liu, Cheng and Hsu2020), the number of particles involved would be of the order of 2 billion which is beyond the current computational power capacity. Furthermore, a separation of scale has to be satisfied and particles should be smaller than the Kolmogorov length scale (
$d_p/\eta <1$). For finite-size particles (
$d_p/\eta >1$), sub-particle scale processes need to be modelled in order to accurately predict the particle dynamics (Finn & Li Reference Finn and Li2016).
Contrary to the fully resolved DNS and the Lagrangian point-particle methodology, the Eulerian–Eulerian two-fluid model has no limitations in term of number of particles. According to Finn & Li (Reference Finn and Li2016), for particle-laden boundary layer simulations, the two-fluid approach is only suited for a narrow range of Stokes numbers 
$0.2 < St < 1$. Indeed for 
$St > 1$, the uniqueness of the Eulerian particle phase velocity field is not guaranteed (Ferry & Balachandar Reference Ferry and Balachandar2001). In other words, for a given fluid phase velocity field, the particles can follow different paths (i.e. the particles velocity field is not unique) depending on the initial condition. Nevertheless, uniqueness of the particle phase velocity is not crucial considering time-averaged particle phase quantities (e.g. concentration, velocity) and assuming ergodicity. More precisely, time-averaged variables are issued from multiple realization of the flow and, therefore, multiple particles trajectories. However, similarly to the point-particle approach, for finite-sized particles, additional sub-particle scale correction models are required (Finn & Li Reference Finn and Li2016).
Over the last three decades, turbulence-resolving two-fluid models have been developed to simulate fluidized beds (O'Brien & Syamlal Reference O'Brien and Syamlal1993; Agrawal et al. Reference Agrawal, Loezos, Syamlal and Sundaresan2001; Heynderickx et al. Reference Heynderickx, Das, De Wilde and Marin2004; Igci et al. Reference Igci, Andrews, Arthur, Sundaresan, Pannala and O'Brien2008; Wang, van der Hoef & Kuipers Reference Wang, van der Hoef and Kuipers2009; Ozel, Fede & Simonin Reference Ozel, Fede and Simonin2013). In this context, the particles are usually inertial (
$St>1$) and smaller than the Kolmogorov length scale (
$d_p/\eta <1$). The clear separation of scale between the fluid flow and the particles allows us to perform two-fluid DNS to fully resolve the turbulent spectrum without approximation. In fluidized beds particles show preferential concentration behaviour resulting in the formation of mesoscale structures such as clusters or streamers that can be captured by the two-fluid model (Agrawal et al. Reference Agrawal, Loezos, Syamlal and Sundaresan2001). Such structures have length scales of the order of 10 to 100 particle diameters and significantly impact the flow dynamics at large scale (Agrawal et al. Reference Agrawal, Loezos, Syamlal and Sundaresan2001). When performing large eddy simulation (LES) in the framework of the two-fluid model, the effect of the unresolved mesoscale structures needs to be incorporated through sub-grid scale closures to accurately predict the two-phase flow dynamics (Agrawal et al. Reference Agrawal, Loezos, Syamlal and Sundaresan2001). Several sub-grid models have been tested by Ozel et al. (Reference Ozel, Fede and Simonin2013) in this context and the functional model for the sub-grid drag force has been shown to perform better. Furthermore, Ozel et al. (Reference Ozel, Fede and Simonin2013) showed that the effect of unresolved mesoscale structures vanishes for a filter width 
$\varDelta$ of the order of the particle diameter.
Recently, Cheng, Hsu & Chauchat (Reference Cheng, Hsu and Chauchat2018) applied the two-fluid LES approach with the functional sub-grid drag model from Ozel et al. (Reference Ozel, Fede and Simonin2013) to reproduce the unidirectional sheet flow experiment from Revil-Baudard et al. (Reference Revil-Baudard, Chauchat, Hurther and Barraud2015). The major difference between fluidized bed configurations mentioned above and the sheet flow configuration comes from the fact that, in the latter, particles are finite sized (
$d_p/\eta >1$). In order to obtain accurate predictions of the flow and the particles dynamics, Cheng et al. (Reference Cheng, Hsu and Chauchat2018) had to use a grid size slightly smaller than the particle diameter (
$d_p/\varDelta \geq 1$). This simulation allowed us to explain, among other things, the physical origin of the modulation of the carrier phase turbulence induced by the presence of particles as being due to the turbulent drag work. However, Cheng et al. (Reference Cheng, Hsu and Chauchat2018) observed an under-estimation of the time-averaged sediment concentration in suspension and a strong sensitivity of the simulation results to the grid resolution. The sub-grid drag model from Ozel et al. (Reference Ozel, Fede and Simonin2013) was originally designed to take into account the effect of unresolved particle clusters and streamers of the order of 10 to 100 particle diameters for coarse-grid simulations (
$d_p/\varDelta <0.1$) and not the effect of mesoscale structures for over-resolved simulations (
$d_p/\varDelta \geq 1$). Therefore, the sub-grid closure used by Cheng et al. (Reference Cheng, Hsu and Chauchat2018) was probably not ideal for this situation, thus explaining the under-prediction of the sediment concentration and the strong sensitivity to the grid resolution. As mentioned earlier, particles are bigger than the Kolmogorov length scale in this configuration. Finite-size effects probably play an important role and shall be modelled.
The modelling of interactions between the carrier phase and finite-size particles has been extensively studied in the literature. Experimental studies (Voth et al. Reference Voth, La Porta, Crawford, Alexander and Bodenschatz2002; Qureshi et al. Reference Qureshi, Bourgoin, Baudet, Cartellier and Gagne2007; Xu & Bodenschatz Reference Xu and Bodenschatz2008) provided evidence that finite-size particle dynamics is substantially affected by turbulent flow scales smaller than the particles compared with particles smaller than the Kolmogorov length scale. All the studies agreed on the facts that the variance of the acceleration probability density functions decreases for increasing particle size. These experimental observations have been further confirmed by numerical studies (Voth et al. Reference Voth, La Porta, Crawford, Alexander and Bodenschatz2002; Calzavarini et al. Reference Calzavarini, Volk, Bourgoin, Lévêque, Pinton and Toschi2009; Homann & Bec Reference Homann and Bec2010; Gorokhovski & Zamansky Reference Gorokhovski and Zamansky2018). One way to recover some of the features of experimental and numerical finite-size particle acceleration probability density functions with the point-particle methodology is to include the Faxén correction term in the fluid–particle interaction model (Calzavarini et al. Reference Calzavarini, Volk, Bourgoin, Lévêque, Pinton and Toschi2009). The Faxén correction term takes into account the non-uniformity of the flow at the particle scale. While this method is suitable for Lagrangian simulations, the methodology developed by Gorokhovski & Zamansky (Reference Gorokhovski and Zamansky2018), taking into account finite-size effect through an effective viscosity at the particle scale included in the expression of the drag force, is more suitable for volume-averaged two-phase flow models.
In this paper the two-fluid LES approach is applied to dilute suspension of finite-sized particles transported in a turbulent boundary layer flow. A finite-size correction model for the two-fluid approach inspired from the model proposed by Gorokhovski & Zamansky (Reference Gorokhovski and Zamansky2018) will be developed and tested against experimental data for particle-laden boundary layer flow configurations having 
$d_p/\eta >1$. In § 2 the two-fluid LES model formulation is presented. In § 3 the numerical results for one clear water configuration and three particle-laden flow configurations are presented with and without the finite-size correction model. In § 4 the sensitivity of the model results to the finite-size correction model components, grid resolution and second filter size are discussed. Finally, conclusions are drawn in § 5.
2. Model formulation
2.1. Filtered two-phase flow equations
To perform LES with a two-phase flow model, a separation between the large turbulent flow scales (low frequency) and the small ones (high frequency) is operated by a filter (operator 
$\bar \cdot$). In analogy with compressible flows, a change of variable called Favre filtering is used to obtain filtered two-phase flow equations. Any variable 
$\psi (x_i, t)$ with 
$x_i=(x, y, z)^\textrm {T}$ the position vector and 
$i=1, 2, 3$ representing three spatial components can be decomposed into the sum 
$\psi (x_i, t) = \tilde \psi (x_i, t) + \psi ^{''}(x_i, t)$, with 
$\tilde \psi (x_i, t)$ the resolved Favre-filtered part and 
$\psi ^{''}(x_i, t)$ the unresolved sub-grid part. Favre-filtered fluid and solid velocities, 
$\tilde u^{f}_i = (\tilde u^{f}, \tilde v^{f}, \tilde w^{f})^\textrm {T}$ and 
$\tilde u^{s}_i = (\tilde u^{s}, \tilde v^{s}, \tilde w^{s})^\textrm {T}$, are defined as
\begin{equation} \tilde u^{f}_i = \frac{\overline{(1-\phi)u^{f}_i}}{(1-\bar\phi)},\quad \tilde u^{s}_i = \frac{\overline{\phi u^{s}_i}}{\bar\phi}, \end{equation}
with 
$\phi$ the solid phase volumetric concentration and 
$u^{f''}_i = u^{f}_i - \tilde u^{f}_i$ and 
$u^{s''}_i = u^{s}_i - \tilde u^{s}_i$ are the sub-grid scale velocity fluctuations.
The filtered two-phase flow equations are composed of the filtered fluid and solid phase continuity equations (2.2) and (2.3) and the filtered fluid and solid phase momentum equations (2.4) and (2.5),
\begin{gather} \frac{\partial (1-\bar\phi)}{\partial t} + \frac{\partial (1- \bar\phi)\tilde u_i^{f}}{\partial x_i} = 0, \end{gather}
\begin{gather}\frac{\partial \bar\phi}{\partial t} + \frac{\partial\bar\phi\tilde u_i^{s}}{\partial x_i} = 0, \end{gather}
\begin{gather} \frac{\partial\rho^{f}(1-\bar\phi)\tilde u^{f}_i}{\partial t} + \frac{\partial\rho^{f}(1-\bar\phi)\tilde u^{f}_i \tilde u^{f}_j}{\partial x_j} ={-}(1-\bar\phi) \frac{\partial \bar P^{f}}{\partial x_i} + \frac{\partial}{\partial x_j} \left(\tilde T^{f}_{ij} + \sigma^{f, sgs}_{ij}\right) + \bar I_i \nonumber\\ \quad + \rho^{f}(1-\bar\phi)g_i + \varPhi^{f,sgs}_i, \end{gather}
\begin{gather} \frac{\partial\rho^{s}\bar\phi\tilde u^{s}_i}{\partial t} + \frac{\partial\rho^{s}\bar\phi\tilde u^{s}_i \tilde u^{s}_j}{\partial x_j} ={-}\bar\phi\frac{\partial \bar P^{f}}{\partial x_i} - \frac{\partial \bar P^{s}}{\partial x_i} + \frac{\partial}{\partial x_j}\left(\tilde T^{s}_{ij} + \sigma^{s, sgs}_{ij}\right) - \bar I_i + \rho^{s}\bar\phi g_i \nonumber\\ \quad +\varPhi^{s,sgs}_i, \end{gather}
with 
$\rho ^{f}$ and 
$\rho ^{s}$ the fluid and solid densities, 
$g_i$ the acceleration of gravity, 
$\bar P^{f}$ and 
$\bar P^{s}$ the filtered fluid and solid pressures, 
$\tilde T^{f}_{ij}$ and 
$\tilde T^{s}_{ij}$ the filtered fluid and solid phase shear stress tensors, 
$\sigma ^{f,sgs}_{ij}$, 
$\sigma ^{s,sgs}_{ij}$, 
$\varPhi ^{f,sgs}_i$ and 
$\varPhi ^{s,sgs}_i$ the fluid and solid sub-grid scale stress tensors and other sub-grid scale contributions respectively presented in § 2.2 and 
$\bar I_i$ the filtered momentum exchange term between the two phases.
The filtered fluid phase shear stress tensor is defined as:
\begin{equation} \tilde T^{f}_{ij} = \rho^{f}(1-\bar\phi)\nu^{f}\left(\frac{\partial\tilde u^{f}_i}{\partial x_j} + \frac{\partial\tilde u^{f}_j}{\partial x_i} - \frac{2}{3} \frac{\partial\tilde u^{f}_k}{\partial x_k}\delta_{ij}\right), \end{equation}
with 
$\nu ^{f}$ the fluid viscosity and 
$\delta _{ij}$ the Kronecker symbol. The filtered solid phase pressure 
$\bar P^{s}$ and shear stress tensor 
$\tilde T^{s}_{ij}$ are calculated using the kinetic theory for granular flows as detailed in § 2.3.
The filtered momentum exchange term 
$\bar I_i$ between the two phases is composed of the drag, lift and added mass forces 
$\bar D_i$, 
$\bar L_i$ and 
$\bar A_i$, respectively, following the expression
\begin{align} \bar I_i = \bar D_i + \bar L_i + \bar A_i \quad \mbox{with}\ \left\{\begin{array}{lll} \bar D_i = \dfrac{\rho^{s}\bar\phi}{\widetilde{t_s}}\left(\tilde u^{f}_i - \tilde u^{s}_i\right), \\ \bar L_i = \bar\phi(1-\bar\phi)C_l\rho^{m}\Vert\tilde u^{f}_i - \tilde u^{s}_i\Vert\epsilon_{ijk} \dfrac{\partial\tilde u^{m}_k}{\partial x_j}, \\ \bar A_i = \bar\phi(1-\bar\phi)C_a\rho^{f}\left[\dfrac{\partial\tilde u^{f}_i}{\partial t} + \dfrac{\partial\tilde u^{f}_i \tilde u^{f}_j}{\partial x_j} - \dfrac{\partial\tilde u^{s}_i}{\partial t} + \dfrac{\partial\tilde u^{s}_i \tilde u^{s}_j}{\partial x_j}\right], \end{array}\right.\end{align}
where 
$C_l=0.5$ and 
$C_a=0.5$ are the lift and added mass coefficients, 
$\rho ^{m}=\bar \phi \rho ^{s} + (1-\bar \phi )\rho ^{f}$ is the volume-averaged mixture density, 
$\tilde u^{m}_i = \bar \phi \tilde u^{s}_i + (1-\bar \phi )\tilde u^{f}_i$ is the mixture velocity and 
$\tilde t_s$ is the particle response time following the drag law proposed by Gidaspow (Reference Gidaspow1986), i.e.
\begin{equation} \tilde t_s = \frac{4}{3}\frac{\rho^{s}}{\rho^{f}} \frac{d_p}{C_D\Vert\tilde u^{f}_i - \tilde u^{s}_i\Vert}(1-\bar\phi)^{2.65} \quad \mbox{with}\ \left\{\begin{array}{ll} C_D = \dfrac{24}{{{Re}}_p}\left(1+0.15{{Re}}_p^{0.687}\right), \\ {{Re}}_p = \dfrac{d_p\Vert\tilde u^{f}_i - \tilde u^{s}_i\Vert}{\nu^{f}}, \end{array}\right. \end{equation}
with 
$C_D$ the drag coefficient from Schiller & Naumann (Reference Schiller and Naumann1933).
2.2. Sub-grid scale modelling
As a direct result of the filtering of the two-phase flow equations, additional sub-grid terms appear in the momentum equations. The fluid and solid phase sub-grid stress tensors 
$\sigma ^{f,sgs}_{ij}=\rho ^{f}(1-\bar \phi )(\widetilde {u^{f}_iu^{f}_j} - \tilde u^{f}_i \tilde u^{f}_j)$ and 
$\sigma ^{s,sgs}_{ij}=\rho ^{s}\bar \phi (\widetilde {u^{s}_iu^{s}_j} - \tilde u^{s}_i \tilde u^{s}_j)$ come from the filtering of the nonlinear advection terms in the momentum equations. Whereas Cheng et al. (Reference Cheng, Hsu and Chauchat2018) modelled the sub-grid stress tensors using the dynamic procedure proposed by Germano et al. (Reference Germano, Piomelli, Moin and Cabot1991) and Lilly (Reference Lilly1992), in the present contribution they are modelled using the dynamic Lagrangian procedure proposed by Meneveau, Lund & Cabot (Reference Meneveau, Lund and Cabot1996). The main difference is that model coefficients are averaged over streamlines (details can be found in appendix A) and not plane-averaged over homogeneous flow directions. The dynamic Lagrangian procedure has the advantage of getting rid of the necessity to have homogeneous directions and of preserving a certain locality in space making it applicable to more complex geometries and inhomogeneous flows in future research (Meneveau et al. Reference Meneveau, Lund and Cabot1996). The sub-grid stress tensors are written as
\begin{equation} \sigma^{f,sgs}_{ij} = 2\rho^{f}(1-\bar\phi)\varDelta^{2} \vert \tilde{\boldsymbol S}^{f} \vert \left(C_1^{f}\tilde S_{ij}^{f} - \tfrac{1}{3}C_2^{f}\tilde S_{kk}^{f}\right), \end{equation}and
\begin{equation} \sigma^{s,sgs}_{ij} = 2\rho^{s}\bar\phi\varDelta^{2} \vert \tilde{\boldsymbol S}^{s} \vert\left(C_1^{s}\tilde S_{ij}^{s} - \tfrac{1}{3}C_2^{s} \tilde S_{kk}^{s}\right), \end{equation}
with 
$\tilde S^{f}_{ij}$ and 
$\tilde S^{s}_{ij}$ the fluid and solid resolved strain rate tensor, respectively, and 
$C_1^{f}$, 
$C_2^{f}$, 
$C_1^{s}$, 
$C_2^{s}$ the dynamically computed model coefficients (details in appendix A).
Other Eulerian–Eulerian sub-grid contributions resulting from the filtering of the pressure, stress and momentum exchange terms are represented by 
$\varPhi ^{f,sgs}_i$ and 
$\varPhi ^{s,sgs}_i$. These sub-grid terms are taking into account the effect of unresolved particle clusters and streamers having length scales smaller than the filter width 
$\varDelta$. Cheng et al. (Reference Cheng, Hsu and Chauchat2018) modelled the sub-grid momentum exchange term using a drift velocity model proposed by Ozel et al. (Reference Ozel, Fede and Simonin2013), but since the typical size of the smallest mesoscale structures is of the order of 10 to 100 particle diameters (Agrawal et al. Reference Agrawal, Loezos, Syamlal and Sundaresan2001), the sub-grid terms taking into account these effects should vanish for filter sizes of the order of the particle size. This has been confirmed by Ozel et al. (Reference Ozel, Fede and Simonin2013) whom quantitatively reported the relative importance of sub-grid terms by explicitly filtering two-phase Eulerian–Eulerian DNS results for different filter sizes. In all the simulations presented in this paper, 
$\varDelta$ is always of the order of the particle size or smaller and, therefore, the sub-grid contributions 
$\varPhi ^{f,sgs}_i$ and 
$\varPhi ^{s,sgs}_i$ can be considered as negligible.
2.3. Particle stress modelling
For solid particles in a boundary layer flow, the solid phase volume fraction changes by several orders of magnitudes from the outer part of the flow to the bottom boundary. Therefore, the modelling methodology used to describe the disperse phase hydrodynamic has to be valid for a wide range of volume fractions from the dilute limit where the interaction with the carrier phase is dominant to higher volume fractions and collision-dominated regimes.
The filtered solid phase pressure 
$\bar P^{s}$ and shear stress tensor 
$\tilde T^{s}_{ij}$ are given by
\begin{equation} \bar P^{s} = \rho^{s}\bar\phi\left[1 + 2(1+e_c)\bar\phi g_{s0}\right]\bar\varTheta - \rho^{s}\lambda\frac{\partial\tilde u^{s}_k}{\partial x_k}\delta_{ij} \end{equation}and
\begin{equation}\tilde T^{s}_{ij} = \rho^{s}\bar\phi\nu^{s}\left(\frac{\partial\tilde u^{s}_i}{\partial x_j} + \frac{\partial\tilde u^{s}_j}{\partial x_i} - \frac{2}{3} \frac{\partial\tilde u^{s}_k}{\partial x_k}\delta_{ij}\right), \end{equation}
respectively, from Gidaspow (Reference Gidaspow1994) with the particle phase shear viscosity 
$\nu ^{s}$ and bulk viscosity 
$\lambda$ given by
\begin{equation}\nu^{s} = d_p\sqrt{\bar\varTheta}\left[\frac{4\bar\phi^{2}g_{s0}(1+e_c)}{\sqrt{5{\rm \pi}}} + \frac{\sqrt{\rm \pi}g_{s0}(1+e_c)^{2}(2e_c-1)\bar\phi^{2}}{15(3-e_c)} + \frac{\sqrt{\rm \pi}\bar\phi}{6(3-e_c)}\right] \end{equation}and
\begin{equation}\lambda = \frac{4}{3}\bar\phi^{2}\rho^{s}d_pg_{s0}(1+e_c)\sqrt{\frac{\bar\varTheta} {\rm \pi}},
\end{equation}
respectively, where 
$g_{s0} = (2-\bar \phi )/2(1-\bar \phi )^{3}$ is the radial distribution function for dense rigid spherical particle gases from Carnahan & Starling (Reference Carnahan and Starling1969), 
$e_c=0.8$ is the restitution coefficient for binary collisions and 
$\bar \varTheta$ is the filtered granular temperature. According to Février, Simonin & Squires (Reference Février, Simonin and Squires2005), 
$\bar \varTheta$ represents the pseudo-thermal kinetic energy associated with the uncorrelated random motions of the particles and should not be confused with the turbulent sub-grid scale turbulent kinetic energy (TKE) 
$k^{s}$ associated with correlated motions of the particles.
The filtered granular temperature 
$\bar \varTheta$ is obtained by solving the following transport equation:
\begin{equation} \frac{3}{2}\left[\frac{\partial\bar\phi\rho^{s}\bar\varTheta}{\partial t} + \frac{\partial\bar\phi\rho^{s}\tilde u^{s}_i\bar\varTheta}{\partial x_j}\right] = \varPi_R + \varPi_q + J_{int} - \gamma + \varPhi^{sgs}_\varTheta. \end{equation}
Here 
$\varPi _R$ is the production of granular temperature by resolved flow scales given by
\begin{equation} \varPi_R = \left(-\bar P^{s}\delta{ij} + \tilde T^{s}_{ij}\right) \frac{\partial \tilde u^{s}_i}{\partial x_j}, \end{equation}
$\varPi _q$ is the divergence of the granular temperature flux analogous to Fourier's law of conduction given by
\begin{gather}\varPi_q = \frac{\partial}{\partial x_i} \left[{-}D_\varTheta\frac{\partial\bar\varTheta}{\partial x_j}\right], \end{gather} where 
$D_\varTheta$ is the conductivity of the granular temperature calculated following
\begin{gather}D_\varTheta = \rho^{s}d_p\sqrt{\bar\varTheta}\left[\frac{2\bar\phi^{2}g_{s0}(1+e_c)}{\sqrt{\rm \pi}} + \frac{9\sqrt{\rm \pi}g_{s0}(1+e_c)^{2}(2e_c-1)\bar\phi^{2}}{2(49-33e_c)} + \frac{5\sqrt{\rm \pi}\bar\phi}{2(49-33e_c)}\right], \end{gather}
$\gamma$ is the dissipation rate of granular temperature given by
\begin{gather}\gamma = 3(1-e_c^{2})\phi^{2}\rho^{s}g_{s0}\bar\varTheta\left[\frac{4}{d_p} \sqrt{\frac{\varTheta}{\rm \pi}}-\frac{\partial u^{s}_j}{\partial x_j}\right], \end{gather}
$J_{int}$ is the fluid–particle interaction term and 
$\varPhi ^{sgs}_\varTheta$ is the sub-grid term.
The fluid–particle interaction term given by
\begin{equation} J_{int} = 3\frac{\rho^{s}\bar\phi}{\tilde t_s}(\bar\varTheta^{f} - \bar\varTheta)\end{equation}
represents the balance between the production of granular temperature due to the fluid pseudo-thermal kinetic energy 
$\bar \varTheta ^{f}$ and the dissipation due to drag. The formulation is similar to the interaction term presented in Fox (Reference Fox2014) but transcripted in the LES formalism with the difference that, for finite-size particles, the fluid pseudo-thermal kinetic energy is confounded with the sub-grid TKE. It is given by 
$\varTheta ^{f} = \frac {2}{3}\alpha k^{f}$, with the sub-grid fluid TKE 
$k^{f}=(C_1^{f}/C_\varepsilon )^{2/3}\vert \tilde {\boldsymbol S}^{f} \vert$, where 
$C_\varepsilon$ is an empirical constant of the order of unity coming from the Kolmogorov theory (Yoshizawa & Horiuti Reference Yoshizawa and Horiuti1985) taken to be 
$C_\varepsilon =1.048$ in recent studies (Arshad et al. Reference Arshad, Gonzalez-Juez, Dasgupta, Menon and Oevermann2019; Chatzimichailidis et al. Reference Chatzimichailidis, Argyropoulos, Assael and Kakosimos2019; Ries et al. Reference Ries, Li, Nishad, Dressler, Ziefuss, Mehdizadeh, Hasse and Sadiki2021) and 
$\alpha =\textrm {e}^{-BSt}$ is a coefficient characterizing the degree of correlation between particles and fluid velocity fluctuations (Hsu, Jenkins & Liu Reference Hsu, Jenkins and Liu2004). The empirical parameter B is a tuning coefficient for Reynolds average models set to 1 in the present simulations.
The sub-grid term 
$\varPhi ^{sgs}_\varTheta =\rho ^{s}\phi \varepsilon ^{s}$, where 
$\varepsilon ^{s}$ is the dissipation rate of solid sub-grid TKE, represents the production of granular temperature from the energy transfer between the correlated solid sub-grid TKE and the granular temperature. Similarly to the sub-grid terms in the momentum equation, 
$\varPhi ^{sgs}_\varTheta$ should vanish for grid sizes of the order of the particle diameter (Agrawal et al. Reference Agrawal, Loezos, Syamlal and Sundaresan2001; Ozel et al. Reference Ozel, Fede and Simonin2013) and is therefore neglected. A summary of the energy transfers between fluid and solid resolved, unresolved, correlated and uncorrelated scales of the flow is presented in figure 1.
Schematic representation of the energy transfers between fluid and solid resolved, unresolved, correlated and uncorrelated scales of the flow with 
$K^{f}$ and 
$K^{s}$ the resolved fluid and solid TKE, 
$k^{f}$ and 
$k^{s}$ the fluid and solid sub-grid TKE, 
$\varTheta ^{f}$ the fluid pseudo-thermal kinetic energy and 
$\varTheta$ the granular temperature. Terms in red are neglected because the grid size 
$\varDelta$ is of the order of the particles diameter 
$d_p$ (Agrawal et al. Reference Agrawal, Loezos, Syamlal and Sundaresan2001; Ozel et al. Reference Ozel, Fede and Simonin2013).
2.4. Finite-size correction model
Compared with particles smaller than the Kolmogorov length scale, finite-size particles do not only act as a temporal filter of the turbulent flow scales through the drag force but also as a spatial filter (Qureshi et al. Reference Qureshi, Bourgoin, Baudet, Cartellier and Gagne2007; Calzavarini et al. Reference Calzavarini, Volk, Bourgoin, Lévêque, Pinton and Toschi2009). In order to take into account the finite-size effect of the particles in the Eulerian–Eulerian two-phase flow model, a distinction is made between turbulent eddies having larger or smaller length scales than the particle diameter 
$d_p$ (blank and hatched zones, respectively, of the idealized turbulent spectrum represented in figure 2). Following observations made by Qureshi et al. (Reference Qureshi, Bourgoin, Baudet, Cartellier and Gagne2007) and Calzavarini et al. (Reference Calzavarini, Volk, Bourgoin, Lévêque, Pinton and Toschi2009), turbulent eddies larger than the particle diameter will contribute to the relative velocity between the two phases in the drag force as fluid velocity ‘seen’ by the particles whereas smaller eddies are assumed to (1) modify the particle response time by increasing the viscosity ‘seen’ by the particles by defining an effective turbulent viscosity at the particle scale following Gorokhovski & Zamansky (Reference Gorokhovski and Zamansky2018), and (2) contribute to particle agitation by increasing the production of granular temperature to be consistent with the energy transfers between correlated and uncorrelated solid phase velocity fluctuations (Février et al. Reference Février, Simonin and Squires2005; Fox Reference Fox2014). The idealized turbulent spectrum presented in figure 2 summarizes the different contributions of the turbulent fluid flow scales to the solid phase dynamics.
Schematic representation of an idealized turbulent spectrum including the different flow scales and their contributions to the particles dynamics (
$\eta$: Kolmogorov scale, 
$\varDelta$: filter width, 
$d_p$: particle diameter, 
$\breve \varDelta$: second filter width, 
$L$: integral scale of turbulence).
To take into account finite-size effects, the filtered effective drag force is rewritten as
\begin{equation} \bar D_i = \frac{\rho^{s} \bar\phi}{\breve t_s}(\breve u^{f}_i - \tilde u^{s}_i), \end{equation}
with 
$\breve u^{f}_i$ the fluid velocity ‘seen’ by the particles corresponding to the resolved fluid phase velocity 
$\tilde u^{f}_i$ filtered at a scale 
$\breve \varDelta \sim O(d_p)$. According to Kidanemariam et al. (Reference Kidanemariam, Chan-Braun, Doychev and Uhlmann2013) the value of 
$\breve \varDelta$ should not be too large to still be relevant to predict the particles motion but not too small to be sufficiently free from the local flow disturbances generated by the presence of the particles. To be able to determine the filter length 
$\breve \varDelta$, they reported the ratio between the averaged magnitude of the flow velocity around spheres and the undisturbed flow field as a function of the distance from the centre of the sphere for different particle Reynolds numbers. From their analysis, around 
$80\,\%$ of the undisturbed mean flow velocity is recovered with a filter width taken as twice the diameter of the particle. Therefore, to compute the fluid velocity ‘seen’ by the particles, the filter size is first chosen to be 
$\breve \varDelta = 2 d_p$. A sensitivity analysis to the filter size is presented in § 4.2.
Whereas the turbulent scales smaller than the particle diameter are usually unresolved, due to the mesh refinement close to the wall, these turbulent scales are composed of both resolved and unresolved eddies in this region. In the present configuration, 
$\breve \varDelta =\varDelta _{x,z}$ in the streamwise and spanwise directions but 
$\breve \varDelta >\varDelta _y$ in the wall-normal direction. To calculate 
$\breve u^{f}_i$, a weighted average of the resolved fluid velocity in the wall-normal direction is performed using a Gaussian distribution with standard deviation 
$\breve \varDelta$ to compute the weighting coefficients.
The new particle response time 
$\breve t_s$ still follows the drag law given by (2.8) but the relative velocity between the two phases is calculated using the filtered fluid velocity 
$\breve u^{f}_i$, and the expression for the particle Reynolds number is modified according to Gorokhovski & Zamansky (Reference Gorokhovski and Zamansky2018) to take into account the effect of turbulent scales smaller than the particles by the mean of a turbulent eddy viscosity 
$\nu ^{t}_p$ at the scale of the particles following
\begin{equation} {{Re}}_p = \frac{d_p\Vert\breve u^{f}_i - \tilde u^{s}_i\Vert}{\nu^{f}+\nu^{t}_p}. \end{equation} The turbulent viscosity at the particle scale can be calculated using Kolmogorov scaling and Prandtl's mixing length hypothesis following 
$\nu ^{t}_p \sim \varepsilon _p^{1/3}d_p^{4/3}$, with 
$\varepsilon _p$ the dissipation of TKE at the particle scale (Gorokhovski & Zamansky Reference Gorokhovski and Zamansky2018). By assuming that the turbulent scales between 
$\breve \varDelta$ and 
$d_p$ are in the inertial range of the turbulent spectrum, the approximation 
$\varepsilon _{\breve \varDelta } = \varepsilon _{p}$ can be made with 
$\varepsilon _{\breve \varDelta }$ the dissipation rate at the filter scale.
The expression of the dissipation rate at the filter scale 
$\varepsilon _{\breve \varDelta }$ is estimated following the expression from Yoshizawa & Horiuti (Reference Yoshizawa and Horiuti1985) defined as a function of the filter width 
$\breve \varDelta$ and the total TKE below 
$\breve \varDelta$ defined as the sum of 
$\breve k=\frac {1}{2}\tilde u^{f''}_i\tilde u^{f''}_i$ the resolved TKE (from 
$\breve \varDelta$ to 
$\varDelta$), with 
$\tilde u^{f''}_i=\tilde u^{f}_i - \breve u^{f}_i$ and 
$k^{f}$ the sub-grid TKE (from 
$\varDelta$ to 
$\eta$), i.e.
\begin{equation} \varepsilon_{\breve\varDelta} = C_\varepsilon\frac{(\breve k + k^{f})^{3/2}}{\breve\varDelta}. \end{equation}Eventually, the particle response time with finite-size correction is written as
\begin{equation} \breve t_s = \frac{4}{3}\frac{\rho^{s}}{\rho^{f}} \frac{d_p}{C_D\Vert\breve u^{f}_i - \tilde u^{s}_i\Vert}(1-\bar\phi)^{2.65} \quad \mbox{with}\ \left\{\begin{array}{ll} C_D = \dfrac{24}{{{Re}}_p}\left(1+0.15{{Re}}_p^{0.687}\right), \\ {{Re}}_p = \dfrac{d_p\Vert\breve u^{f}_i - \tilde u^{s}_i\Vert} {\nu^{f}+\varepsilon_{\breve\varDelta}^{1/3}d_p^{4/3}}. \end{array}\right. \end{equation} Furthermore, the turbulent flow scales below 
$\breve \varDelta$ contribute to increase the production of granular temperature isotropically. The fluid–particleinteraction term 
$J_{int}$ in (2.15) includes the resolved sub-particle TKE, i.e.
\begin{equation} J_{int} = \frac{\rho^{s}}{\breve t_s}\frac{\phi}{1-\phi}\left[2\alpha (\breve k+k^{f}) - 3\bar\varTheta\right]. \end{equation} It shall be mentioned that the proposed model tends to the two-fluid model in its traditional formulation for particles smaller than the Kolmogorov length scale (
$d_p/\eta <1$). Indeed, if 
$\breve \varDelta \leq \eta$ then 
$\breve u^{f}_i=\tilde u^{f}_i$ and, therefore, 
$\breve k = k^{f} =0$. The turbulent viscosity at the particle scale vanishes 
$\nu ^{t}_p=0$ and eventually 
$\breve t_s = \tilde t_s$. Furthermore, the proposed finite-size correction model allows us to approach the exact solution for vanishingly small mesh sizes. Indeed, the correction model allows us to decorrelate the filter width associated with the particle size in the drag law and the filter size imposed by the mesh for the LES making the solution mesh-independent.
2.5. Numerical implementation
The present model is adapted from the turbulence averaged two-phase flow solver sedFoam (https://github.com/sedFoam/sedFoam) (Chauchat et al. Reference Chauchat, Cheng, Nagel, Bonamy and Hsu2017; Cheng, Hsu & Calantoni Reference Cheng, Hsu and Calantoni2017). It is implemented in the open-source computational fluid dynamics toolbox OpenFoam (Jasak & Uroić Reference Jasak and Uroić2020) and solves the Eulerian–Eulerian two-phase flow mass and momentum equations using a finite volume method and a pressure-implicit with splitting of operators (PISO) algorithm for velocity-pressure coupling (Rusche Reference Rusche2002). In the PISO algorithm, at each time step, intermediate velocities are first computed by solving the momentum equations without the pressure gradient term. Then, the Poisson equation for the pressure is solved in order to calculate the corrected pressure field and ensure mass conservation. Eventually, the velocity fields are corrected based on the new pressure field. Several steps can be applied to the velocity prediction-correction to increase convergence (nCorrectors in OpenFoam). In the present simulations, nCorrectors = 2 is sufficient for convergence. More information about consistency, algorithm and numerical implementation can be found in Chauchat et al. (Reference Chauchat, Cheng, Nagel, Bonamy and Hsu2017).
In the present simulations the same numerical schemes as in Cheng et al. (Reference Cheng, Hsu and Chauchat2018) are used to provide a second-order accuracy in both space and time. A second-order implicit backward scheme is used for temporal derivatives (denoted as backward in OpenFoam) and a second-order total variation diminish scheme is used for the mass conservation equation and the granular temperature transport equation (denoted as limitedLinear in OpenFoam). For the advection terms in the momentum equations, a second-order centred scheme is used for which high frequency filtering of the oscillations induced by second-order discretization is performed by introducing a small amount of upwind scheme (denoted as filteredLinear in OpenFoam). The gradient are computed using a second-order centred scheme (denoted as linear in OpenFoam).
3. Results
In this section numerical simulations performed on different flow configurations are presented in order to assess the two-fluid LES model presented in § 2. First, a clear water configuration, i.e. without particles (
$\bar \phi =0$), is presented and compared with existing experimental and numerical DNS data to validate the model, the choice of the grid resolution and the numerical schemes. Second, three particle-laden flow configurations involving finite-sized particles are reproduced numerically to evaluate the capability of the two-fluid LES model, including the finite-size correction model, to predict the turbulent suspension of particles.
3.1. Clear water configuration
The clear water configuration from Kiger & Pan (Reference Kiger and Pan2002) consists of a closed unidirectional channel flow with Reynolds number 
${ {Re}}_\tau =u_\tau h/\nu ^{f}=560$ based on the wall-friction velocity 
$u_\tau = 2.8\times 10^{-2}\ \textrm {m}\ \textrm {s}^{-1}$ and channel half-height 
$h=0.02\ \textrm {m}$.
The numerical domain is a bi-periodic rectangular box (figure 3) with cyclic boundary conditions in 
$x$ and 
$z$ directions and no slip boundary condition at the top and bottom boundaries for the velocities. The gradient of any other quantities is set to zero at the walls. The flow is driven by a pressure gradient along the 
$x$-axis dynamically adjusted at each time step in order to match the experimental bulk velocity 
$U_b=0.51\ \textrm {m}\ \textrm {s}^{-1}$. The mesh is composed of 
$314\times 220\times 160$ elements corresponding to a total of 
$11,105,280$ cells. The spanwise and streamwise resolution is constant with 
$\varDelta _x^{+}\approx \varDelta _z^{+} \approx 11$ wall units (
$+$ symbol with 
$\psi ^{+} = \psi u_\tau /\nu ^{f}$). The mesh is stretched along the 
$y$-axis with 
$\varDelta _y^{+}\approx 1$ at the wall and 
$\varDelta _y^{+} \approx 6$ at the centreline. The time step is fixed to 
${\rm \Delta} t=10^{-4}\ \textrm {s}$ to ensure a maximum Courant-Friedrichs-Lewy number (CFL) lower than 
$0.3$ for stability reasons. In a recent publication, Montecchia et al. (Reference Montecchia, Brethouwer, Wallin, Johansson and Knacke2019) performed a sensitivity analysis to the CFL number (
$\textrm {CFL}=0.1$, 0.2 and 0.3) and the results did not show strong sensitivities.
Sketch of the geometry and boundary conditions of the numerical domain for the simulation of the clear water configuration and configuration GB from Kiger & Pan (Reference Kiger and Pan2002).
All the simulations presented in this paper are initialized with a fully developed turbulent boundary layer flow obtained from a preliminary simulation. A first run is conducted to let the turbulence develop until the wall-friction velocity and the integral of the total flow kinetic energy have reached a steady state. This corresponds to approximately 
$200 T_b$, with 
$T_b=h/U_b$ the non-dimensional bulk time scale of the flow. Then, a second run is performed to compute turbulence statistics and Favre-averaged quantities over a duration of 
$200 T_b$. The Favre-averaging procedure is represented by the operator 
$\langle \cdot \rangle _F$ (details can be found in appendix B). In clear water flow conditions, Favre averaging is equivalent to ensemble averaging denoted as 
$\langle \cdot \rangle$.
Similarly to what has been done by Kiger & Pan (Reference Kiger and Pan2002), the average profiles obtained experimentally and numerically are compared with the profiles from the DNS of Moser, Kim & Mansour (Reference Moser, Kim and Mansour1999) with 
${ {Re}}_\tau =590$. Since the Reynolds number in the configuration from Kiger & Pan (Reference Kiger and Pan2002) is close to the DNS, it is reasonable to compare the profiles between the two configurations.
The averaged velocity profiles, Reynolds stress and root-mean-square (r.m.s.) of the streamwise velocity fluctuations 
$\tilde u^{f'}_{rms}$ and wall-normal velocity fluctuations 
$\tilde v^{f'}_{rms}$ are presented in figure 4 in wall units. In the simulations, the friction velocity is calculated based on the time-averaged streamwise pressure gradient as
\begin{equation} u_\tau = \sqrt{h \left\langle \frac{\partial \bar P^{f}}{\partial x} \right\rangle}. \end{equation}
The computed wall-friction Reynolds number is equal to 
${ {Re}}_\tau = 544$ (
$u_\tau =2.72\times 10^{-2}\ \textrm {m}\ \textrm {s}^{-1}$) which corresponds to an error below 
$3\,\%$ compared with the experiments.
Average profiles of velocity in (a), Reynolds stress in (b), r.m.s. of streamwise velocity fluctuations in (c) and r.m.s. of wall-normal velocity fluctuations in (d) from the two-fluid model (T.F. model) compared with the numerical results from Moser et al. (Reference Moser, Kim and Mansour1999) (DNS) and experimental data from Kiger & Pan (Reference Kiger and Pan2002) (Exp.).
The present clear water simulation produces profiles of averaged velocity and turbulence statistics that agree very well with the DNS and experimental data. However, especially for the Reynolds stress, some discrepancies between experimental measurements and the simulations appear near the wall. Kiger & Pan (Reference Kiger and Pan2002) stated that their measurements can be considered highly reliable in the outer log layer with less than 
$5\,\%$ error for 
$y^{+}> 50$ and up to 
$25\,\%$ variability for 
$y^{+} < 50$.
The agreement between numerical and experimental data confirms that without solid particles, the two-phase flow model behaves exactly as a single-phase flow model. The accurate prediction of the flow hydrodynamics and turbulent statistics allows us to validate the model implementation, the choices of numerical parameters and gives confidence to perform particle-laden simulations in the next sections.
3.2. Particle-laden configurations
In this section particle-laden configurations involving spherical glass beads (GB) from Kiger & Pan (Reference Kiger and Pan2002), natural sediment (NS) particles from Muste et al. (Reference Muste, Yu, Fujita and Ettema2005) and almost neutrally buoyant sediment (NBS) particles from Muste et al. (Reference Muste, Yu, Fujita and Ettema2005) are reproduced numerically. The flow and particles parameters are presented in table 1.
Flow and particles parameters for configurations GB, NS, NBS and NS*.
The targeted configurations correspond to the turbulent dilute suspended sediment transport boundary layer flows. In this situation, particles are entrained into suspension by the turbulent coherent flow structures and under steady-state flow conditions, an equilibrium concentration profile across the water depth establishes as the result of an equilibrium between the gravity driven settling flux 
$v_s\langle \phi \rangle$, with 
$v_s$ the settling velocity of the particles, and the turbulent Reynolds sediment flux 
$\langle v^{s'}\phi '\rangle$ (Rouse Reference Rouse1938), with 
$v^{s'}$ the solid phase vertical velocity fluctuations and 
$\phi '$ the sediment concentration fluctuations. By analogy with Fickian diffusion, this Reynolds sediment flux is modelled using a gradient diffusion model. Introducing this model in the Reynolds-averaged sediment mass balance leads to the equation
\begin{equation} v_s\langle\phi\rangle - \frac{\nu^{f}_t}{S_c}\frac{{\textrm{d}}\langle\phi\rangle}{\textrm{d}y} = 0, \end{equation}
with 
$\nu ^{f}_t$ the turbulent eddy viscosity (or turbulent momentum diffusivity) and 
$S_c$ the turbulent Schmidt number representing the efficiency of the sediment diffusion relative to 
$\nu _t^{f}$. For 
$S_c <1$, sediment particles are dispersed more efficiently by turbulence than fluid parcels.
Given that 
$\nu ^{t}_f=l_m^{2} {\textrm {d}}\langle u^{f} \rangle /{\textrm {d}y}$ using Prandtl's mixing length 
$l_m=\kappa y$, with 
$\kappa =0.41$ the von Kármán constant, and using the log-law-of-the-wall to describe the mean velocity profile, (3.2) can be integrated analytically to give the following expression for the Reynolds-averaged particle concentration profile:
\begin{equation} \frac{\langle\phi\rangle}{\phi_0} = \left(\frac{y_0}{y}\right)^{Ro}. \end{equation}
Here 
$\phi _0$ is a reference concentration at a given reference elevation 
$y_0$ and 
$R_o=S_c v_s/u_\tau \kappa$ is the Rouse number. For open-channel flows, a free surface correction to the Prandtl's mixing length is introduced 
$l_m=\kappa y \sqrt {1-y/h}$, with 
$h$ representing the water depth, and the Reynolds-averaged particle concentration profile reads as
\begin{equation} \frac{\langle\phi\rangle}{\phi_0} =\left[\frac{y}{h-y}\frac{h-y_0}{y_0} \right]^{{-}Ro}. \end{equation} These two analytical solutions provide a reference with which the two-fluid LES model results can be compared. The value of 
$S_c$ is still debated in the sediment transport community (Lyn Reference Lyn2008), the most widely accepted model is the one from van Rijn (Reference van Rijn1984) relating the turbulent Schmidt number to the suspension number, 
$S_c=(1+2(v_s/u_\tau )^{2})^{-1}$. Nevertheless, a lot of scatter is observed on existing experimental data and no satisfactory explanation exists to support van Rijn's empirical formula (Lyn Reference Lyn2008).
The hydrodynamic configuration, numerical domain and parameters for configuration GB are the same as the clear water case presented in § 3.1. The only difference comes from the addition of a given amount of particles in the flow corresponding to a mean volumetric concentration of particles in the channel 
$\phi _{tot}=2.31\times 10^{-4}$. The particles are spherical and mono-dispersed with diameter 
$d_p=195\ \mathrm {\mu } \textrm {m}$ (
$d_p^{+}\approx 5.5$) and density 
$\rho ^{s} = 2600\ \textrm {kg}\ \textrm {m}^{-3}$. For such particles, the computed fall velocity in still water using the drag law from (2.8) is 
$v_s=2.4\times 10^{-2} \ \textrm {m}\ \textrm {s}^{-1}$ (
$v_s/u_\tau =0.85$).
Configurations NS and NBS from Muste et al. (Reference Muste, Yu, Fujita and Ettema2005) consist of turbulent particle-laden open-channel flows with water depth 
$h=0.021\ \textrm {m}$ in which finite-sized particles with density 
$\rho ^{s}=2650$ and 
$\rho ^{s}=1025\ \textrm {kg}\ \textrm {m}^{-3}$, respectively, are seeded. The NS and NBS hydrodynamic conditions are the same with a bulk velocity 
$U_b=0.84\ \textrm {m}\ \textrm {s}^{-1}$ and a targeted friction velocity 
$u_\tau =4.2\times 10^{-2}\ \textrm {m}\ \textrm {s}^{-1}$ corresponding to a Reynolds number based on the wall-friction velocity 
${ {Re}}_\tau =882$. Both types of particles have the same diameter 
$d_p=230 \ \mathrm {\mu } \textrm {m}$ (
$d_p^{+}\approx 9.7$) resulting in a larger fall velocity (i.e. larger suspension number) for NS 
$v_s=2.4\times 10^{-2}\ \textrm {m}\ \textrm {s}^{-1}$ (
$v_s/u_\tau =0.57$) compared with NBS 
$v_s=6\times 10^{-4}\ \textrm {m}\ \textrm {s}^{-1}$ (
$v_s/u_\tau =0.01$). For both configurations, the mean volumetric concentration of sediment is equal to 
$\phi _{tot}=4.6\times 10^{-4}$. The computational domain is a rectangular box with bi-periodic boundary conditions along the 
$x$- and 
$z$-axis and no slip boundary conditions at the bottom boundary (figure 5).
Sketch of the geometry and boundary conditions of the numerical domain for configurations NS, NBS and NS* from Muste et al. (Reference Muste, Yu, Fujita and Ettema2005).
The mesh is composed of 8 323 000 cells with uniform streamwise and spanwise grid resolutions 
$\varDelta _x^{+} = \varDelta _z^{+} = 19$. The mesh resolution is stretched along the 
$y$-axis with the first grid point located at 
$\varDelta _y^{+} \approx 1$ and 
$\varDelta _y^{+} \approx 3$ at the top. Here 
${ {Re}}_p$, 
$St$ and 
$d_p/\eta$ are calculated based on the scaling analysis from Finn & Li (Reference Finn and Li2016).
A first set of simulations for each configuration is performed in order to evaluate the predictive capability of the two-fluid model without finite-size correction. The visualization of the instantaneous turbulent coherent structures using an isocontour of Q-criterion (figure 6a) and volume rendering of the concentration (figure 6b) from the GB configuration shows the imprint of turbulence on the sediment concentration field. The differences of sediment concentration relative to the coherent structures highlights the importance of turbulence-particle interactions. The averaged solid phase concentration profiles obtained experimentally and numerically are compared in figure 7. For GB (figure 7a), experimental and numerical concentration profiles are normalized by the reference concentration 
$\phi _0$ taken at 
$y_0=0.06h$. For both configurations GB and NS (figures 7a and 7b), the volume fraction of particles in suspension is significantly under-estimated compared with the experimental data. However, for the NBS configuration (figure 7c), the average concentration profile predicted by the two-fluid model fits perfectly well the experimental results. For 
$v_s/u_\tau \ll 1$, the weight of the particles is entirely supported by turbulence (Berzi & Fraccarollo Reference Berzi and Fraccarollo2016). The two-phase flow model in its original formulation correctly reproduces the vertical balance between settling and Reynolds fluxes. For this flow and these particle parameters, finite-size effects can be considered as negligible and the two-fluid model shows very good predictive capabilities without the finite-size correction model. In the following, configurations GB and NS for which the suspension number is higher are further investigated to understand the physical origin of the observed discrepancies.
Visualization of the instantaneous turbulent coherent structures using an isocontour of Q-criterion coloured by the velocity (a) and volume rendering of the concentration (b) from the GB configuration.
Solid phase volumetric concentration profiles from the experiments (Exp.) and two-phase flow simulations (T.F. model) from configurations GB (a), NS (b) and NBS (c). In panel (a) experimental and numerical concentration profiles are normalized by the reference concentration 
$\phi _0$ taken at 
$y_0=0.06h$.
The research hypothesis developed in this work is that the discrepancies observed in figure 7 are due to finite-size effects. One could also argue that these discrepancies are due to missing fluid–particle forces such as added mass and lift forces. A simulation including these two forces have been performed for the configuration GB and the averaged concentration profiles are compared with the experiments and the analytical concentration profile from (3.3) in figure 8. Since this expression is derived for an infinite boundary layer, one has to keep in mind that for closed channel flows, this expression could become less accurate near the centreline of the channel. The comparison between the simulation including only the drag force and the simulation including drag, lift and added mass forces indicates that the drag force is the dominant interaction force for this configuration. The contributions from lift and added mass forces are almost negligible in this problem. The concentration profiles from both simulations show a power law that fits with (3.3) for 
$Ro=2.85$, whereas 
$Ro=2.04$ in the experiments (figure 8b). The Basset history force does not appear in the momentum exchange term between the two phases since it is defined from a purely Lagrangian point of view. It would therefore be very difficult to obtain a volume average expression of the history force in the Eulerian formalism. To the best of the authors’ knowledge, there is no references in the literature showing the Eulerian expression of the Basset history force. The authors believe that the Basset history force would be significant very near the bottom boundary where wall-particle collisions occur but should not affect too much the vertical distribution of particles in the upper part of the channel where particle acceleration is weaker. It is only through a detailed comparison with Lagrangian point-particle simulations including the Basset history force that the role of this force could be investigated which is beyond the scope of the present work.
Solid phase volumetric concentration profiles from the experiment (Exp.), two-phase flow simulation including only the drag force (T.F. model (
$\bar D_i$)) and two-phase flow simulation including drag, lift and added mass forces (T.F. model (
$\bar D_i + \bar L_i + \bar A_i$)) from configuration GB in semi-log scale in panel (a) and in log-log compared with analytical profiles from (3.3) with 
$Ro=2.04$ and 
$Ro=2.85$ in (b).
The two-phase flow model in its initial formulation, using a standard drag law, added mass and lift forces, can not reproduce the turbulent suspension of particles in this configuration. In the following, the role of unresolved turbulent length scales smaller than the particle size is investigated.
3.3. Evaluation of the finite-size correction model
From the scaling analysis proposed by Balachandar (Reference Balachandar2009), the relative velocity between the fluid phase and inertial particles is mainly influenced by an eddy having the same time scale as the particles with the corresponding length scale 
$l^{*}=\varepsilon ^{f1/2}t_s^{3/2}$, with 
$\varepsilon ^{f}$ the dissipation rate of fluid TKE. According to Finn & Li (Reference Finn and Li2016), if 
$l^{*}>\varDelta >d_p$, all the relevant flow scales are resolved and the particle dynamics can be accurately predicted. The average value of 
$l^{*}$ for configuration GB is calculated in the simulation and plotted in figure 9. It can be seen that 
$l^{*}/d_p<1$ and that 
$l^{*}$ decreases by one order of magnitude from the wall to the centreline of the channel. This result shows that, for this configuration, turbulent scales smaller than the particles can have a significant effect on the particle dynamics and may be responsible for the observed discrepancies.
Average profile of the length scale 
$l^{*}$ associated with the turbulent eddy having the same time scale as the particles as a fraction of 
$d_p$ from the two-phase simulation of configuration GB.
Given the broad range of length and time scales involved in a particle-laden horizontal boundary layer, multiple types of turbulence-particle interactions occur at different locations of the boundary layer. It is therefore crucial to develop a model applicable over a wide range of turbulence-particle interaction regimes. In the following, the finite-size correction model presented in § 2.4 is tested for the three configurations GB, NS and NBS. The results of the simulations for the averaged concentration profile for configuration GB with and without the finite-size correction model are compared in figure 10. The prediction of the concentration profile by the two-phase flow model is significantly improved by the finite-size correction model without any tuning coefficient. The Rouse number predicted with the finite-size correction model is 
$Ro=1.86$ which is much closer to the experimental value compared with the prediction without correction (figure 10b). However, in the experiment the concentration profile is well described by the power law across the water depth, whereas in the simulation the concentration decreases more rapidly toward the centreline of the channel.
Solid phase volumetric concentration profiles from the experiment (Exp.), two-phase flow simulation with finite-size correction model (T.F. model (FS)) and two-phase flow simulation without finite-size correction model (T.F. model) from configuration GB in semi-log scale in (a) and in log-log compared with analytical profiles from (3.3) with 
$Ro=2.04$, 
$Ro=2.90$ and 
$Ro=1.86$ in (b).
In order to further evaluate the finite-size correction model, the configurations NS and NBS are reproduced numerically using the two-phase flow model with finite-size correction. Analytical, experimental and numerical averaged concentration profiles are compared in figure 11 for both configurations.
Solid phase volumetric concentration profiles from the experiment (Exp.), two-phase flow simulation with finite-size correction model (T.F. model (FS)), two-phase flow simulation without finite-size correction model (T.F. model) and analytical profiles from (3.4) with 
$Ro=0.83$, 
$Ro=1.00$, 
$Ro=2.29$ and 
$Ro=0.07$ from configuration NS in (a) and configuration NBS in (b).
For configuration NS (figure 11a), the same conclusions as for configuration GB can be drawn. The finite-size correction model significantly improves the prediction of the turbulent suspension of particles without any tuning coefficient. The predicted Rouse number (
$Ro=0.83$) is closer to the experimental value (
$Ro=1.00$). The modelled concentration profile obtained using finite-size correction is in very good agreement with the experimental data compared with the simulation without the finite-size correction model.
For configuration NBS (figure 11b), the finite-size correction model does not alter the results predicted without finite-size correction. Almost no differences can be observed between the concentration profiles obtained with and without finite-size correction confirming that finite-size effects are negligible for configurations with a low suspension number.
As a partial conclusion, it has been demonstrated that finite-size effects are important to predict turbulent suspension of inertial particles in a boundary layer flow when the suspension number is of the order of unity. The finite-size correction model proposed in this work significantly improves the model prediction for the average sediment concentration profile without the use of tuning parameter to fit the experimental data.
3.4. Lag velocity
Another interesting feature of turbulent suspension of inertial particles is the existence of a velocity lag between the average streamwise velocity of the fluid and of the particles (Kaftori, Hetsroni & Banerjee Reference Kaftori, Hetsroni and Banerjee1995; Niño & Garcia Reference Niño and Garcia1996; Kiger & Pan Reference Kiger and Pan2002; Righetti & Romano Reference Righetti and Romano2004; Muste et al. Reference Muste, Yu, Fujita and Ettema2005; Kidanemariam et al. Reference Kidanemariam, Chan-Braun, Doychev and Uhlmann2013). Kidanemariam et al. (Reference Kidanemariam, Chan-Braun, Doychev and Uhlmann2013), based on fully resolved DNS, have been able to clearly identify the physical origin of this velocity lag as being due to the preferential concentration of suspended particles in low speed regions of the fluid flow which can be identified with ejection events. This velocity lag is not observed for particle-laden flows with a low suspension number (Muste et al. Reference Muste, Yu, Fujita and Ettema2005) such as NBS but can be as high as 20 % of the bulk fluid velocity (Kidanemariam et al. Reference Kidanemariam, Chan-Braun, Doychev and Uhlmann2013).
The averaged fluid and solid velocity profiles obtained numerically with and without the finite-size correction for configurations GB, NS and NBS are shown in the top panels of figure 12. The velocity profiles are in very good agreement with the experiments and they do not show much sensitivity to the finite-size correction model. The velocity difference is too small to be visible on these graphs, the lag velocity 
$u_{lag} = \langle \tilde u^{f}_i \rangle - \langle \tilde u^{s}_i \rangle$ is shown in the bottom panels of figure 12. For configurations GB and NS, the lag velocity is positive and of the order of 5–10 % of the bulk fluid velocity. The two-fluid model predicts the correct sign and order of magnitude for the lag velocity. The major discrepancy is observed in the near-wall region 
$y/h<0.2$ where the two-fluid model predicts a peak that is not observed in Kiger & Pan (Reference Kiger and Pan2002) experiments. For the NS configuration, the lag velocity decreases linearly with the distance to the free surface. This is probably a free surface effect that is not fully captured by the symmetry plane boundary condition used in the present simulation, nevertheless, the model predictions are very satisfactory. Given the role played by the coherent structures of the flow in the velocity lag, the fact that the flow is over-resolved near the bottom boundary could explain the observed discrepancies for the GB case. However, since this difference is not observed for the other configurations, it might also be due to the high variability of the velocity measurements of Kiger & Pan (Reference Kiger and Pan2002) near the bottom (up to 25 % for 
$y^{+}<50$). In both GB and NS configurations the finite-size correction model has a small influence on the lag velocity. In the NBS configuration the experimental data reveals a negligible lag velocity that even becomes negative. The two-fluid model with and without the finite-size correction model predicts a zero lag velocity except very near the bottom wall. From these three configurations, one can conclude that the existence of a lag velocity is not due to finite-size effects. More importantly, the fact that the model is able to recover the absence of lag velocity for NBS means that the two-fluid LES captures the physical mechanism correctly and can be used as a predictive tool to study this mechanism.
Averaged fluid and solid velocity profiles in (a,b,c) and lag velocity in (d,e,f) from the two-fluid model with finite-size correction (T.F. model (FS)) and the two-fluid model without finite-size correction (T.F. model) for configurations GB (a,d), NS (b,e) and NBS (c,f) compared with experimental data (Exp.).
3.5. Turbulent statistics
Among the three configurations, the most accurate measurements of turbulent statistics have been obtained for configuration GB. In the following, this configuration is analysed in detail for the fluid and particle phase flow statistics.
The wall-friction velocity for configuration GB predicted with and without finite-size correction is 
$u_\tau =2.70\times 10^{-2}$ and 
$u_\tau =2.72\times 10^{-2}\ \textrm {m}\ \textrm {s}^{-1}$, respectively. Whereas the numerical wall-friction velocity is similar between the clear water and particle-laden configurations, the experiments suggest an increase of the friction velocity up to 
$u_\tau =2.99\times 10^{-2}\ \textrm {m}\ \textrm {s}^{-1}$. Averaged fluid and solid Reynolds stress and the r.m.s. of streamwise and wall-normal velocity fluctuation profiles from configuration GB with and without finite-size correction are compared with experimental data in figure 13. From figure 13(a), the two-fluid model slightly under-estimates the fluid Reynolds shear stress compared with the experiments explaining the lower friction velocity in the simulations. However, experimental and numerical results are similar: the solid phase Reynolds shear stress is slightly greater than the fluid Reynolds shear stress away from the bottom wall. The maximum value for the solid Reynolds shear stress predicted by the two-fluid LES model is the same as in the experiments but the location is different. The r.m.s. of streamwise and wall-normal velocity fluctuations are in very good agreement with experimental results (figures 13b and 13c). As for the fluid Reynolds shear stress, the fluid phase velocity fluctuations are slightly under-estimated by the two-fluid model for 
$y/h>0.1$. For both experimental and numerical profiles, the r.m.s. of streamwise solid phase velocity fluctuations are equal near the centreline of the channel and becomes smaller in the near bottom wall region. Similarly to the observation of Kidanemariam et al. (Reference Kidanemariam, Chan-Braun, Doychev and Uhlmann2013) in their fully resolved DNS, the two-fluid model predicts stronger wall-normal solid velocity fluctuations compared with the fluid away from the wall, whereas experimental solid and fluid profiles are similar close to the channel centreline. The r.m.s. of the solid velocity fluctuations decreases more rapidly than the fluid ones towards the wall. Overall, the turbulent statistics are not significantly affected by the finite-size correction model.
Average profiles of fluid and solid Reynolds stress in (a), r.m.s. of streamwise velocity fluctuations in (b) and r.m.s. of wall-normal velocity fluctuations in (c) from the two-phase model with finite-size correction (T.F. model (FS)) and the two-phase model without finite-size correction (T.F. model) compared with experimental data (Exp.) from configuration GB.
The slight differences observed between the experimental and numerical fluid phase turbulent statistics come from the modulation of the turbulence by the particles. Again, from the scaling analysis by Finn & Li (Reference Finn and Li2016) and given the parameters of the configuration from Kiger & Pan (Reference Kiger and Pan2002), the presence of the particles is expected to dampen fluid turbulence, whereas in the experiments a slight increase of the Reynolds stress and velocity fluctuations are observed compared with the clear water configuration. Indeed, according to Balachandar (Reference Balachandar2009), the turbulence enhancement due to the presence of the particles comes from the combined action of the oscillating wakes behind particles having a high particle Reynolds number. The conjugate action of all the wakes of the particles participates to increase the overall fluid turbulence.
To be able to predict the turbulence enhancement, the two-fluid model should have the capacity to capture the vortex shedding behind the particles by fully resolving the fluid/solid interface which is not the case for the Eulerian–Eulerian two-phase flow model. Nevertheless, the turbulence enhancement due to the particles is not a dominant mechanism in this configuration. According to Finn & Li (Reference Finn and Li2016), the net production of turbulence by the particles is dominant for particle Reynolds numbers higher than Reynolds number 
${ {Re}}_p=400$ even if oscillatory wakes behind particles can be observed for lower 
${ {Re}}_p$ depending on flow properties, particle shape or distance from the wall for example. In the present configuration, the particle Reynolds number based on the scaling relations from Finn & Li (Reference Finn and Li2016) is equal to 
${ {Re}}_{p} = 4.8$ and the maximum particle Reynolds number predicted in the simulation is 
${ {Re}}_{p,max}\approx 20$, which is significantly below the threshold Reynolds number of 400.
Flow hydrodynamics and turbulent statistics are in good agreement with experimental data and the overall relative behaviour between the fluid and solid phase is correctly captured by the two-phase flow model. The fact that the two-fluid flow model does not resolve the particle–fluid interface implies that the turbulence enhancement induced by the presence of the particles is not resolved. However, for such flow and particle parameters, according to the scaling analyses from Balachandar (Reference Balachandar2009) and Finn & Li (Reference Finn and Li2016), this mechanism is not dominant. The lower fluid velocity fluctuations predicted by the two-fluid model near the channel centreline only results in a slight under-estimation of the sediment concentration in the same region compared with the experiments.
From the resolved and sub-grid TKE profiles for the fluid and solid phases presented in figure 14, it appears that most of the TKE is resolved (
$k^{f}/K^{f} < 5\,\%$ and 
$k^{s}/K^{s}< 5\,\%$). The fact that the solid phase sub-grid TKE is equal to zero through the channel height shows that the resolution is very close to DNS and validates the hypothesis made in § 2 to neglect the sub-grid terms. Furthermore, 
$K^{s} \gg \varTheta$ in the channel except very near the solid boundary (
$y/h<0.05$) showing that kinetic and collisional dispersive forces should not be dominant compared with the drag force to suspend the solid particles for such dilute configurations. This hypothesis is further investigated in § 4.
Resolved and sub-grid fluid phase TKE 
$K^{f}$ and 
$k^{f}$ in (a), resolved and sub-grid solid phase TKE 
$K^{s}$ and 
$k^{s}$ and granular temperature 
$\varTheta$ in (b) made dimensionless by the friction velocity 
$u_\tau$ from configuration GB.
3.6. Volume fraction sensitivity
An additional simulation (configuration NS*) is performed to evaluate the robustness of the proposed model to an increase of the mass loading. Experimental data from Muste et al. (Reference Muste, Yu, Fujita and Ettema2005) using NS having higher volume fractions compared with the previous NS configuration is reproduced numerically using the two-fluid model. The hydrodynamic and particle parameters are the same as for the NS configuration but the total solid phase volume fraction is multiplied by a factor 3.5 (
$\phi _{tot}=16.2\times 10^{-4}$). The concentration can still be considered dilute and particles do not form a settled bed at the bottom of the channel.
The averaged concentration profile from the simulation of configuration NS with higher volume fractions is compared with experimental data in figure 15. As observed in the previous sections, the agreement with experimental data is significantly improved by the finite-size correction for natural sediments. This result is even more spectacular considering that without the correction model, the particles settle almost completely at the bottom of the channel resulting in an even larger under-estimation of the suspension of particles at higher volume fractions by the original two-fluid model.
Solid phase volumetric concentration profiles from the experiment (Exp.), two-phase flow simulation with finite-size correction model (T.F. model (FS)), two-phase flow simulation without finite-size correction model (T.F. model) from configuration NS*.
4. Discussion
In this section the sensitivity of the model to the different components of the finite-size correction model are discussed as well as the sensitivity to the grid/second filter resolution is presented.
4.1. Relative influence of the different terms of the finite-size correction model
In order to evaluate the relative influence of the modified drag law and the modified production of granular temperature, a new simulation is performed for which finite-size effects are taken into account only in the production term of the granular temperature transport equation. In other words, the simulation is performed using the drag law from (2.8) and the production of granular temperature from (2.25). The average concentration profile obtained from the simulation including finite-size effects only in the production of granular temperature equation is compared with the concentration profile from the simulations with and without finite-size correction in figure 16.
Solid phase volumetric concentration profiles from the two-phase flow simulation with finite-size correction model (T.F. model (FS)), two-phase flow simulation with finite-size correction only in the production term of granular temperature (T.F. model (FS-
$J_{int}$ only) and two-phase flow simulation without finite-size correction model (T.F.model) from configuration GB in semi-log scale.
The concentration profile obtained from the two-fluid simulation including finite-size effects only in the production term of the granular temperature is similar to the profile without the finite-size correction model. Indeed, for dilute particles, fluid–particle interactions are dominant compared with particle–particle interactions. The slope of the concentration profile in dilute regions of the flow is shaped by the drag force and the modification of the production term of granular temperature has almost no effect. However, the effect of the modification of the granular temperature transport equation could become dominant for higher concentrations. It should be noted that the modification of the granular temperature transport equation is necessary because a simulation including finite-size effects only in the drag law and not in the production term of granular temperature was shown to be highly unstable. The observed instability is not the result of a numerical issue but rather the consequence of a physical inconsistency in the energy transfers between the two phases. Indeed, including finite-size effects only in the drag law and not in the production term of granular temperature is not physically accurate. Some of the turbulent flow scales are not taken into account in the energy budget between the phases making the simulation unstable.
The averaged wall-normal momentum balance for the solid phase can be written as
\begin{equation}
\underbrace{\left\langle\bar\phi \rho^{s}
g\right\rangle}_{{\rm Gravity}} =
\left\langle\underbrace{\frac{\partial R^{s}_{xy}}{\partial
x} + \frac{\partial R^{s}_{yy}}{\partial y} +
\frac{\partial R^{s}_{zy}}{\partial
z}}_{\substack{{\rm Reynolds} \\ {\rm stresses}}}
\right\rangle +
\underbrace{\left\langle\bar\phi\frac{\partial \bar
P^{f}}{\partial y}\right\rangle }_{\substack{{\rm Fluid}\\ {\rm pressure} \\ {\rm (Buoyancy)}}} +
\underbrace{\left\langle\frac{\partial \bar P^{s}}{\partial
y}\right\rangle }_{\substack{{\rm Granular} \\
{\rm pressure}}} + \underbrace{\left\langle\bar
I_y\right\rangle}_{{\rm Drag}},
\end{equation}
with 
$R_{ij}=\partial \langle \rho ^{s} \bar \phi \tilde u_i^{s'}\tilde u_j^{s'}\rangle /\partial x_j$ the Reynolds shear stress coming from averaging of the nonlinear advection terms with 
$\tilde u_i^{s'}= \tilde u_i^{s} - \langle \tilde u_i^{s} \rangle$ the solid phase resolved velocity fluctuations. In figure 17 the four terms on the right-hand side of (4.1) are plotted in dimensionless form, normalized by 
$\langle \rho ^{s}\bar \phi g\rangle$, in semi-log scale for 
$y/h$. This figure shows that the predicted suspended particle concentration profile results from a balance between gravity, buoyancy and drag forces in the upper part of the channel. In such dilute systems the effect of dispersive kinetic and collisional forces are not significant except very near the wall 
$y/h<0.05$ (
$y/d_p<5$). This supports the hypothesis that the discrepancies observed in the original model can not be due to a flaw in the kinetic theory formulation but are due to fluid–particle interaction forces.
Averaged contributions to the wall-normal momentum budget 
$\langle C_y \rangle$ on the right-hand side (RHS) of the momentum balance (4.1) and their sum as a fraction of the gravity force 
$\langle \bar \phi \rho ^{s}g\rangle$ from configuration GB.
4.2. Second filter size sensitivity
As mentioned in § 2.4, the width of the second filter 
$\breve \varDelta$ should not be too small to be free from disturbances generated by the presence of the particles. On the other hand, the second filter size should not be too large in order to provide an accurate representation of the velocity ‘seen’ by the particles. The minimum filter width 
$\breve \varDelta _{min}=2d_p$ has been determined from Kidanemariam et al. (Reference Kidanemariam, Chan-Braun, Doychev and Uhlmann2013), but there is no clear criteria for the maximum filter width. However, for computational efficiency, since the second filter width depends on the spatial discretization in the streamwise and spanwise directions for this configuration, it can be crucial to determine the maximum acceptable filter width to accurately predict the average concentration profile with coarser grid resolutions.
Additional simulations for configuration GB with different mesh resolutions have been performed to measure the influence of the spatial discretization and the second filter width 
$\breve \varDelta$ on the sediment concentration profile prediction. The mesh characteristics for the different simulations are presented in table 2. The comparison between Reynolds fluxes and concentration profiles with and without finite-size correction obtained with mesh M1, M2, M3 and M4 are presented in figure 18.
Averaged Reynolds flux (a) and solid phase volumetric concentration (b) profiles from two-phase flow simulations with finite-size correction model (T.F. model (FS)) from configuration GB using mesh M1, M2, M3 and M4 and without finite-size correction model (T.F. model) using mesh M1 and M4.
Mesh characteristics for the second filter size sensitivity test.
The turbulent dispersion of the particles increases for a coarser resolution (figure 18a). As a consequence, the amount of suspended particles in the water column predicted by the two-fluid model increases with increasing filter width (figure 18b). The difference between the concentration profiles from simulations using the finite-size correction model with 
$\breve \varDelta =2 d_p$ (mesh M1) and 
$\breve \varDelta =3 d_p$ (mesh M2) is negligible and the agreement can still be considered as acceptable for a filter width of 
$\breve \varDelta =4 d_p$ (mesh M3). However, discrepancies become important for a larger filter width (
$\breve \varDelta =8 d_p$ with mesh M4). More quantitatively, the concentration profile converges at first order with the reference simulation results from mesh M1 for increasing vertical resolution.
Even without the finite-size correction model, the Reynolds flux is increased between simulations using mesh M1 and mesh M4 (figure 18a) suggesting that the over-prediction of the concentration does come from the finite-size correction model only but also from the modification of the flow hydrodynamic for coarser grid resolutions.
As a conclusion, in order to accurately predict the concentration profile, the sensitivity analysis suggests that the grid resolution at the wall should not exceed 4 wall units 
$\varDelta _y^{+} < 4$. It is mandatory to have at least one grid point in the laminar sub-layer in order to resolve the turbulent coherent flow structures in the near-wall region and to use a second filter smaller than 
$8 d_p$ (
$\breve \varDelta < 8d_p$) to accurately resolve the fluid velocity ‘seen’ by the particles.
5. Conclusion
Turbulence-particle interactions may play a key role in particle-laden flows by modifying the turbulent dispersion of particles by turbulent eddies and by the feedback of particles on the turbulent eddies. From a modelling point of view, a specific challenge is the huge range of cascading turbulent eddy sizes 
$O(10^{-1}\text {--}10^{-4})$ m and their interactions with different grain sizes 
$O(10^{-3}\text {--}10^{-5})$ m. The very wide range of length scales involved does not allow us to systematically use the turbulence-resolving approach at the particle scale to address this problem due to its prohibitive computational cost and turbulence-resolving continuum approaches, such as the two-fluid LES approach, are needed.
In this paper the two-fluid LES method has been tested against experimental data and a finite-size correction model has been developed. The new model has been validated against available experimental data for the dilute turbulent suspension of finite-sized particles transported by a boundary layer flow. The improved model has been shown to accurately predict the suspended particle concentration profile as well as the existence of a streamwise lag velocity for heavier-than-fluid particles without the use of any tuning parameter to fit the experimental data. In the proposed correction model, a distinction is made between turbulent flow scales larger or smaller than the particle diameter. The velocity field ‘seen’ by the particles in the drag law is filtered at a scale 
$\breve \varDelta \geq 2d_p$ and smaller turbulent scales contribute to reduce the particle response time by the addition of a sub-particle scale eddy viscosity to the molecular viscosity in the particle Reynolds number definition. The second effect of the correction is to increase the production of granular temperature by a modification of the source term in the granular temperature equation. While modification of the drag law is more important for the accurate prediction of the suspended particle concentration profile in the dilute configuration investigated herein, the modification of the granular temperature equation is mandatory for the physical consistency and the numerical stability of the model. At last, the sensitivity analysis of the model results for the second filter size 
$\breve \varDelta$ has shown that the grid resolution could be as high as 4 particle diameters without loss of accuracy as long as one grid point is located in the laminar sub-layer.
The work presented herein is an important step towards two-fluid LES of more complex applications such as scour around hydraulic structures, wave-driven sediment transport or turbidity currents to cite a few geophysical flow examples. In all the aforementioned applications, additional complex interaction mechanisms such as interactions with a sediment bed occur. The two-fluid approach should therefore be validated against configurations involving higher sediment volume fractions compared with the present configurations.
Furthermore, detailed measurements of turbulence-particle interactions are really challenging and using turbulence-resolving simulations in addition to measurements is probably the only way to improve our understanding of the role of these mechanisms on particle transport dynamics.
Acknowledgements
J.C. and A.M. would like to thank R. Zamansky for the fruitful discussions about the finite-size correction model. Most of the computations presented in this paper were performed using the GENCI infrastructure under allocations A0060107567 and A0080107567 and the GRICAD infrastructure. LEGI is member of Labex TEC21 (Investissements dAvenir Grant Agreement ANR-11-LABX-0030) and Labex Osug@2020 (Investissements dAvenir Grant Agreement ANR-10-LABX-0056).
Funding
The authors would like to acknowledge the financial support from Agence de l'Innovation de Défense (AID), Shom through project MEPELS and Agence Nationale de la Recherche (ANR) through project SheetFlow (ANR-18-CE01-0003). T.-J.H. and J.C. also like to acknowledge support from the Munitions Response Program of the Strategic Environmental Research and Development Program under Project MR20-S1-1478.
Declaration of interest
The authors report no conflict of interest.
Appendix A. Dynamic Lagrangian model from Meneveau et al. (Reference Meneveau, Lund and Cabot1996)
As a result of the nonlinear advection terms filtering, additional sub-grid terms need to be modelled 
$\sigma ^{a,sgs}_{ij}=\rho ^{a}\bar \phi ^{a}(\widetilde {u^{a}_iu^{a}_j} - \tilde u^{a}_i \tilde u^{a}_j)$, with subscript 
$a=\{f,s\}$ denoting the fluid or solid phase, 
$\rho ^{a}$ and 
$\bar \phi ^{a}$ the density and filtered volumetric concentration of phase 
$a$, respectively. The most common way to model the sub-grid stress is to use the Smagorinsky model,
\begin{equation} \sigma^{a,sgs}_{ij} = 2\rho^{a}\bar\phi^{a}\varDelta^{2} \vert \tilde{\boldsymbol S}^{a} \vert \left(C_1^{a}\tilde S_{ij}^{a} - \tfrac{1}{3}C_2^{a}\tilde S_{kk}^{a}\delta_{ij}\right), \end{equation}
with 
$\varDelta$ the filtered width, 
$\tilde {\boldsymbol S}^{a}$ resolved strain rate tensor of phase 
$a$ and 
$C_1^{a}$ and 
$C_2^{a}$ the model coefficients. To adjust the model coefficients, a dynamic procedure samples the turbulent stress from the smallest resolved scales and extrapolate to determine the turbulent stress associated with unresolved turbulent scales below 
$\varDelta$.
The starting point to determine the first coefficient 
$C_a^{a}$ is the algebraic identity
\begin{equation} \mathcal{L}_{ij}^{a} = \mathcal{T}_{ij}^{a} - \mathcal{\tau}_{ij}^{a}, \end{equation}
relating the turbulent stress associated to two different filter widths 
$\varDelta$ and 
$\hat \varDelta =2\varDelta$ with
\begin{equation} \mathcal{L}_{ij}^{a}=\widehat{\tilde u^{a}_i \tilde u^{a}_j} - \hat{\tilde u}^{a}_i \hat{\tilde u}^{a}_j,\quad \mathcal{T}_{ij}^{a}=\widehat{\widetilde{u^{a}_iu^{a}_j}} - \hat{\tilde u}^{a}_i \hat{\tilde u}^{a}_j \quad \mbox{and}\quad \mathcal{\tau}_{ij}^{a}=\widehat{\widetilde{u^{a}_iu^{a}_j}} - \widehat{\tilde u^{a}_i \tilde u^{a}_j}. \end{equation} The Smagorinsky model is used to model the turbulent stress 
$\mathcal {\tau }_{ij}^{a}$ at scale 
$\varDelta$ and 
$\mathcal {T}_{ij}^{a}$ at scale 
$2\varDelta$:
\begin{gather} \mathcal{\tau}_{ij}^{a} ={-}2C_1^{a}\varDelta^{2}\vert \tilde{\boldsymbol S}^{a} \vert{\tilde S_{ij}^{a}}, \end{gather}
\begin{gather}\mathcal{T}_{ij}^{a} ={-}2C_1^{a}(2\varDelta)^{2}\vert \widehat{\tilde{\boldsymbol S}^{a}} \vert\widehat{\tilde S_{ij}^{a}}. \end{gather} Replacing expressions (A4) and (A5) in the identity (A2) and minimizing the mean square error between the resolved identity and the Smagorinsky model leads to the expression for the coefficient 
$C_1^{a}$,
\begin{equation} C_1^{a}=\frac{\mathcal{F}_{LM}^{a}}{\mathcal{F}_{MM}^{a}}, \end{equation}
with 
$\mathcal {F}_{LM}^{a}$ and 
$\mathcal {F}_{MM}^{a}$ the products 
$\mathcal {L}_{ij}^{a}\mathcal {M}_{ij}^{a}$ and 
$\mathcal {M}_{ij}^{a}\mathcal {M}_{ij}^{a}$ with 
$\mathcal {M}_{ij}^{a}=2\varDelta ^{2}[\widehat {\vert \tilde S^{a}\vert \tilde S_{ij}^{a}}-4\vert \hat {\tilde S}^{a}\vert \hat {\tilde S}_{ij}^{a}]$ averaged over streamlines following the expressions:
\begin{align} \mathcal{F}_{LM}^{a}=\int^{t}_{-\infty}\mathcal{L}_{ij}^{a}\mathcal{M}_{ij}^{a}(z(t'), t')W(t-t')\,{\textrm{d}}t', \end{align}
\begin{align} \mathcal{F}_{MM}^{a}=\int^{t}_{-\infty}\mathcal{M}_{ij}^{a}\mathcal{M}_{ij}^{a} (z(t'), t')W(t-t')\,{\textrm{d}}t'. \end{align}
The weighting function 
$W(t-t')$ is used to control the relative importance of the events near time 
$t$ with those of earlier times. The exponential form 
$W(t-t')=T^{-1}\,\textrm {e}^{-(t-t')/T}$, with 
$T=({3\varDelta }/{2})(\mathcal {F}_{LM}^{a}\mathcal {F}_{MM}^{a})^{-({1}/{8})}$, makes the integrals (A7) and (A8) solutions of the transport equations
\begin{equation} \frac{\partial \bar\phi^{a} \mathcal{F}_{LM}^{a}}{\partial t} + \frac{\partial\bar\phi^{a} \tilde u_j^{a}\mathcal{F}_{LM}^{a}}{\partial x_i} = \frac{2}{3\varDelta}(\mathcal{F}_{LM}^{a}\mathcal{F}_{MM}^{a})^{{1}/{8}} (\mathcal{L}_{ij}^{a}\mathcal{M}_{ij}^{a} - \mathcal{F}_{LM}^{a}) \end{equation}and
\begin{equation} \frac{\partial \rho^{a} \bar\phi^{a} \mathcal{F}_{MM}^{a}}{\partial t} + \frac{\partial \rho^{a} \bar\phi^{a} \tilde u_j^{a}\mathcal{F}_{MM}^{a}}{\partial x_i} = \frac{2}{3\varDelta}(\mathcal{F}_{LM}^{a}\mathcal{F}_{MM}^{a})^{{1}/{8}} (\mathcal{M}_{ij}^{a}\mathcal{M}_{ij}^{a} - \mathcal{F}_{MM}^{a}). \end{equation} To compute the second model coefficient 
$C_2^{a}$, a similar procedure is used. Using the following identity for the spherical part of the sub-grid scale shear stress tensor:
\begin{equation} \mathcal{L}^{*,a} =\mathcal{T}^{*,a} - \mathcal{\tau}^{*,a}. \end{equation}Here
\begin{equation} \mathcal{L}^{*,a}=\tfrac{1}{3}\textrm{tr}(\widehat{\tilde u^{a}_i \tilde u^{a}_j} - \hat{\tilde u}^{a}_i \hat{\tilde u}^{a}_j),\quad \mathcal{T}^{*,a}=\tfrac{1}{3}\textrm{tr}(\widehat{\widetilde{u^{a}_iu^{a}_j}} - \hat{\tilde u}^{a}_i \hat{\tilde u}^{a}_j) \quad \mbox{and}\quad \mathcal{\tau}^{*,a}=\tfrac{1}{3}\textrm{tr}(\widehat{\widetilde{u^{a}_iu^{a}_j}} - \widehat{\tilde u^{a}_i \tilde u^{a}_j}), \end{equation}modelled following
\begin{gather} \mathcal{\tau}^{*,a} ={-}\tfrac{2}{3}C_2^{a}\varDelta^{2}\vert \tilde{\boldsymbol S}^{a} \vert{\tilde S_{kk}^{a}\delta_{ij}}, \end{gather}
\begin{gather}\mathcal{T}^{*,a} ={-}\tfrac{2}{3}C_2^{a}(2\varDelta)^{2}\vert \widehat{\tilde{\boldsymbol S}^{a}} \vert\widehat{\tilde S_{kk}^{a}}\delta_{ij}. \end{gather} Minimizing the mean square error between the resolved identity and the Smagorinsky model leads to the expression for the coefficient 
$C_2^{a}$,
\begin{equation} C_2^{a}=\frac{\langle\mathcal{L}^{*,a}\rangle_C}{\langle\mathcal{M}^{*,a}\mathcal{M}^{*,a}\rangle_C}, \end{equation}
with 
$\mathcal {M}^{*,a}=-\frac {2}{3}\varDelta ^{2}[\widehat {\vert \tilde S^{a}\vert \tilde S_{kk}^{a}\delta _{ij}}-4\vert \hat {\tilde S}^{a}\vert \hat {\tilde S}_{kk}^{a}\delta _{ij}]$ and operator 
$\langle \cdot \rangle _C$ representing average over the cell faces.
Appendix B. Averaging procedure
The given variable 
$\psi$ can be decomposed into the sum of the Favre-averaged variable 
$\langle \psi \rangle _F$ and the associated fluctuation 
$\psi '$. Favre-averaging operations on the variables 
$\psi ^{f}$ or 
$\psi ^{s}$ correspond to perform an ensemble average (operator 
$\langle \cdot \rangle$) of the variables weighted by the ensemble-averaged phase concentration 
$\langle 1-\phi \rangle$ or 
$\langle \phi \rangle$ following
\begin{equation} \langle \psi^{f} \rangle_F = \frac{\langle(1-\phi)\psi^{f}\rangle}{\langle 1-\bar\phi\rangle}, \quad \langle \psi^{s} \rangle_F = \frac{\langle\phi\psi^{f}\rangle}{\langle\phi\rangle}. \end{equation} Numerically, averaged variables are calculated by performing a spatial averaging operation in the streamwise and spanwise directions of a temporally averaged variable 
$\langle \psi \rangle _t$ following
\begin{equation} \langle \psi\rangle = \frac{1}{L_xL_z}\int_0^{L_x}\int_0^{L_z} \langle \psi \rangle_t {\textrm{d}x}\,{\textrm{d}}z, \end{equation}
with 
$L_x$ and 
$L_z$ the lengths of the numerical domain in the streamwise and spanwise directions, respectively.
The temporal averaging operation is performed using an iterative procedure at each time step with the temporal average value of the variable 
$\psi$ at time 
$t_{n+1}$ given by
\begin{equation} \langle\psi(t_{n+1})\rangle_{t} = \frac{\psi(t_{n+1}) + n\langle\psi(t_n)\rangle_t}{n+1}. \end{equation} Second-order statistical moments such as r.m.s. of the velocity fluctuations or Reynolds stresses are obtained by calculating the fluid or solid Favre-averaged covariance tensor 
$\langle \psi '_i\psi '_j\rangle _F$ following
\begin{equation} \langle\psi'_i\psi'_j\rangle_F = \langle\psi_i\psi_j\rangle_F - \langle\psi_i\rangle_F\langle\psi_j\rangle_F. \end{equation}One can notice that in clear water conditions (without solid phase), fluid phase Favre averaging is equivalent to ensemble averaging.
