2.1 Eulerian–Lagrangian model

The Eulerian–Lagrangian model is based on the explicit solution of the dynamic equation for each individual particle using the 3D DEM open-source code YADE (Šmilauer et al. Reference Šmilauer, Catalano, Chareyre, Dorofeenko, Duriez, Dyck, Elias, Er, Eulitz and Gladky2015). The DEM solution is spatially averaged and explicitly coupled with the fluid phase momentum balance equations (2.1)–(2.2) through the drag term $n\langle f_{x}^{D}\rangle ^{p}$ and the solid volume fraction $\unicode[STIX]{x1D719}$ .

For each particle $p$ at position $\boldsymbol{x}^{p}$ , the Newton equations are solved considering nearest-neighbours interactions:

where the sum of the contact forces $\boldsymbol{f}_{c}^{pk}$ is made over the ensemble of nearest neighbours $\mathscr{N}$ , $\boldsymbol{f}_{g}^{p}$ is the gravity force, $\boldsymbol{f}_{b}^{p}$ is the buoyancy force and $\boldsymbol{f}_{D}^{p}$ is the 3D drag force applied by the fluid on particle $p$ . For each contact, the contact force is computed explicitly using the classical spring–dashpot contact law (Schwager & Pöschel Reference Schwager and Pöschel2007):

where $F_{n}$ and $F_{t}$ are the normal and tangential contact forces between particle $p$ and $k$ , $\unicode[STIX]{x1D6FF}_{n}$ and $\unicode[STIX]{x1D6FF}_{t}$ are the normal and tangential overlaps, $k_{n}$ and $k_{s}$ are the normal and tangential contact stiffnesses, $c_{n}$ is the normal viscous damping and $\unicode[STIX]{x1D707}^{p}$ is the tangential friction coefficient. This contact law is well suited for granular flow analysis and allows one to define a unique restitution coefficient $e_{n}$ (Schwager & Pöschel Reference Schwager and Pöschel2007). Similarly, the rotation of the particles is solved from the Newton equations of motion.

The interaction with the fluid phase is restricted to buoyancy and drag forces:

where the average fluid velocity and the fluid pressure are taken at the centre of particle $p$ . The drag coefficient $C_{D}$ depends on the particle Reynolds number $Re_{p}=\Vert \langle \boldsymbol{u}\rangle _{\boldsymbol{x}^{\boldsymbol{p}}}^{f}-\boldsymbol{v}^{\boldsymbol{p}}\Vert d/\unicode[STIX]{x1D708}^{f}$ and takes into account hindrance effects (DallaValle Reference Dallavalle1948; Richardson & Zaki Reference Richardson and Zaki1954): $C_{D}=(0.4+24.4/Re_{p})(1-\unicode[STIX]{x1D719})^{-3.1}$ .

Considering steady uniform turbulent flows, the buoyancy force reduces to its wall-normal component and (2.2) leads to a hydrostatic fluid pressure distribution, the wall-normal average drag force being negligible. The solution of the streamwise fluid phase momentum balance (2.1) requires the evaluation of both the spatially averaged solid phase volume fraction and the momentum transmitted from the fluid to the particles through drag forces:

In the DEM processing, the solid phase averaging $\langle \bullet \rangle ^{s}$ is used instead of the particle phase averaging $\langle \bullet \rangle ^{p}$ . While the two are identical provided that there is scale separation (Anderson & Jackson Reference Anderson and Jackson1967; Jackson Reference Jackson2000), the importance of the wall-normal gradients in bedload transport requires a weighting function length scale lower than the particle diameter in order to define an independent averaging (Maurin Reference Maurin2015; Maurin et al. Reference Maurin, Chauchat, Chareyre and Frey2015). Therefore, a small wall-normal weighting function length scale has been adopted (typically $d/30$ ), and this choice has been validated through an experimental comparison (Maurin et al. Reference Maurin, Chauchat, Chareyre and Frey2015).

The 3D DEM and the fluid model are solved as transient problems applying a fixed bottom boundary condition for both the fluid ( $\langle \boldsymbol{u}\rangle ^{f}(z=0)=0$ ) and the particle phase (fixed random particles) and imposing the position of the water free surface ( $\text{d}\langle u_{x}\rangle ^{f}/\text{d}z(z=h)=0$ ). In order to achieve a stable integration, the DEM time step is bounded by the propagation time of the fastest wave over a particle diameter (Maurin Reference Maurin2015; Maurin et al. Reference Maurin, Chauchat, Chareyre and Frey2015). The fluid resolution time step corresponds to a typical characteristic evolution time scale of the granular medium and is taken much larger than the DEM one (Maurin et al. Reference Maurin, Chauchat, Chareyre and Frey2015): $\unicode[STIX]{x0394}t_{f}=10^{-2}~\text{s}$ with respect to $\unicode[STIX]{x0394}t_{p}\sim O(10^{-4}{-}10^{-5})~\text{s}$ . Therefore, the coupling between the fluid and granular phases ensures momentum conservation of the system on average. The model has been compared with experiments and has shown its ability to describe accurately the granular depth structure in turbulent bedload transport (Maurin et al. Reference Maurin, Chauchat, Chareyre and Frey2015).

2.2 Eulerian–Eulerian model

The Eulerian–Eulerian model is based on the numerical solution of the fluid phase momentum balance equations (2.1) and (2.2) coupled with the spatially averaged granular phase momentum balance equations in the streamwise and wall-normal directions. For steady uniform flow conditions the granular phase momentum equations read (Chauchat et al. Reference Chauchat, Revil-Baudard and Hurther2015)

where $\unicode[STIX]{x1D70F}_{xz}^{p}$ is the spatially averaged granular phase shear stress. Like in the Eulerian–Lagrangian model, the buoyancy is applied along the wall-normal direction and leads to hydrostatic pressure distribution for the granular phase. In the continuous formalism the spatially averaged streamwise drag force is given by

where the extra term $\boldsymbol{u}_{\boldsymbol{d}}$ is the drift velocity which represents the dispersion effect due to the averaged fluid–particle velocity fluctuations. This term is responsible for sediment transport suspension and allows one to reproduce the Rouse profile provided that the wall-normal component is taken as $(u_{d})_{z}=-\unicode[STIX]{x1D708}^{t}/\unicode[STIX]{x1D719}~\text{d}\unicode[STIX]{x1D719}/\text{d}z$ (Chauchat et al. Reference Chauchat, Revil-Baudard and Hurther2015). The streamwise component is taken as $(u_{d})_{x}=0$ in this paper.

The granular shear stress tensor is modelled using the $\unicode[STIX]{x1D707}(I)$ rheology presented in the introduction. Combining equations (1.1) and (1.2), the shear and normal components of the granular stress tensor can be expressed as

Inverting equation (1.3), the inertial number is obtained as a function of the solid volume fraction, and the two equations can be directly expressed as a function of the solid volume fraction and the average solid velocity gradient.

The numerical algorithm is based on a transient solution of the governing equations (2.1), (2.2), (2.11) and (2.12), together with the continuity equation for the solid phase. Fixed boundary conditions are imposed at the bottom of the sample at $z=0$ ( $\langle \boldsymbol{u}\rangle ^{f}=0$ ) where the granular phase is considered to be at rest ( $\langle \boldsymbol{v}^{p}\rangle ^{s}=0$ , $\text{d}\unicode[STIX]{x1D719}/\text{d}z=0$ ), while the positions of the water and granular free surfaces are imposed at $z=h$ ( $\text{d}\langle u_{x}\rangle ^{f}/\text{d}z=0$ , $\langle u_{z}\rangle ^{f}=0$ , $\langle v_{z}^{p}\rangle ^{s}=0$ , $\text{d}\unicode[STIX]{x1D719}/\text{d}z=0$ ). The governing equations are discretised using a finite volume technique for the mass conservation equation and a finite difference technique for the momentum balance. A staggered grid is used with the velocities located at the cell face and the scalar quantities (e.g. volume fractions or viscosities) located at the cell centre. For the pressure–velocity coupling a projection method is used. The numerical schemes are Euler implicit for the time derivative, upwind for the advection terms and central difference for the diffusion terms. For additional informations on the numerical model, the reader is referred to Chauchat et al. (Reference Chauchat, Guillou, Pham Van Bang and Dan Nguyen2013) and Chauchat et al. (Reference Chauchat, Revil-Baudard and Hurther2015).

Figure 1.

Scheme of the numerical set-up and its equivalent average unidirectional picture with typical fluid velocity $\langle \boldsymbol{u}^{f}\rangle =\langle u_{x}^{f}\rangle (z)\boldsymbol{e}_{\boldsymbol{x}}$ , solid volume fraction $\unicode[STIX]{x1D719}$ and solid velocity $\langle \boldsymbol{v}\rangle ^{s}=\langle v_{x}\rangle ^{s}(z)\boldsymbol{e}_{\boldsymbol{x}}$ depth profiles. The inclined 3D biperiodic granular description is coupled with a unidirectional fluid momentum balance using imposed fixed random bottom boundary condition and water free-surface position. The particle colour is representative of the velocity intensity.

4.1 Parametrisation of the $\unicode[STIX]{x1D707}(I)$ rheology

Figure 7 shows the granular shear to normal stress ratio and the solid volume fraction as a function of the inertial number, for simulations with Shields numbers from $\unicode[STIX]{x1D703}^{\ast }=0.04$ to $\unicode[STIX]{x1D703}^{\ast }=0.6$ and for particle diameters $d\in \{3;6;12\}$ mm. In these new configurations closer to realistic bedload transport, most of the data are seen to collapse for both the solid volume fraction and the granular stress ratio up to high inertial numbers. Above $I\sim 3$ and under $I\sim 10^{-2}$ , the curves split apart similarly to the general case. No dependence upon the particle diameter is observed over the whole range of inertial numbers, pointing out the negligible effect of the $Re_{p}$ and $St$ numbers in the problem. The lowest-Shields-number cases (lightest points) exhibit some specific behaviour with oscillations, lower shear to normal stress ratio and solid volume fraction at given inertial numbers. These cases correspond to the sharpest transitions between a quasistatic granular bed and a few particles in motion on top of it. Therefore, in these configurations the continuous assumption is not strictly valid and the limits of the rheological description are reached.

Figure 7.

Shear to normal stress ratio and solid volume fraction as a function of the dry inertial number, for realistic bedload transport simulations with specific density $\unicode[STIX]{x1D70C}^{p}/\unicode[STIX]{x1D70C}^{f}-1=1.5$ , Shields number from incipient motion to $0.7$ and particle diameter of $d=3$  mm,  $d=6$  mm and $d=12$  mm. The black lines represent the best fit obtained with equations (1.2) and (4.1) ( $\unicode[STIX]{x1D707}_{1}=0.38$ , $\unicode[STIX]{x1D707}_{2}=0.64$ , $I_{0}=0.279$ , $\unicode[STIX]{x1D719}^{max}=0.61$ and $b=0.31$ ). The symbol colours are functions of the simulation Shields number: the darker the colour, the higher the Shields number.

Overriding these limitations close to the critical Shields number and assuming that the high-inertial-number granular behaviour is not dominant in bedload transport, the $\unicode[STIX]{x1D707}(I)$ constitutive laws are pragmatically best-fitted on the Eulerian–Lagrangian simulation results over the range of inertial numbers $I_{dry}\in [10^{-2},3]$ . The classical solid volume fraction expression (1.3) predicts negative values at high inertial numbers and is unable to describe the present results. Consequently, the expression is modified and is chosen following Aussillous et al. (Reference Aussillous, Chauchat, Pailha, Médale and Guazzelli2013) and Revil-Baudard & Chauchat (Reference Revil-Baudard and Chauchat2013) as

The best fits of equations (1.2) and (4.1) on the data give the following parameters: $b=0.31$ , $\unicode[STIX]{x1D707}_{s}=0.35$ , $\unicode[STIX]{x1D707}_{2}=0.97$ and $I_{0}=0.69$ , with $\unicode[STIX]{x1D719}^{max}$ imposed to take the value measured in the Eulerian–Lagrangian simulations: $\unicode[STIX]{x1D719}^{max}=0.61$ . Developing the solid volume fraction expression around $I=0$ , the value of $b$ obtained is consistent with both the parametrisation of Da Cruz et al. (Reference Da Cruz, Emam, Prochnow, Roux and Chevoir2005) in the dry regime and the one of Boyer et al. (Reference Boyer, Guazzelli and Pouliquen2011) and Trulsson et al. (Reference Trulsson, Andreotti and Claudin2012) in the viscous regime. The values of $\unicode[STIX]{x1D707}_{2}$ and $I_{0}$ are logically higher than in the original $\unicode[STIX]{x1D707}(I)$ rheology, as no real saturation of the stress ratio is observed in our data. The best fits describe very well the results obtained within the collapsing range of $I_{dry}$ , as can be seen in figure 7. Outside this interval, they do not reproduce the non-local effects at low inertial numbers ( $I_{dry}<10^{-2}$ ) and the drop in both the shear to normal stress ratio and the solid volume fraction for $I_{dry}>3$ .

Assuming the dense granular behaviour to be dominant, the two relationships (1.2)–(4.1) together with the values determined from the best fits are proposed as a $\unicode[STIX]{x1D707}(I)$ -based granular rheology for turbulent bedload transport.

4.2 Eulerian–Eulerian simulations

In order to evaluate the relevance of the proposed granular rheology it has been implemented in the Eulerian–Eulerian model presented in § 2. Simulations have been performed using exactly the same parameters for both the Eulerian–Lagrangian and the Eulerian–Eulerian models ( $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{f}-1=1.5$ ; $d=6$  mm, $\unicode[STIX]{x1D6FC}=0.05$ , $h$ ). The particle diameter has been varied ( $d\in 3,6,12$  mm) but does not influence the results, so that only the case $d=6$  mm is presented for clarity.

Figure 8.

Dimensionless sediment transport rate as a function of the Shields number for two-phase continuous simulations (Eulerian–Eulerian) and coupled fluid–discrete-element simulations (Eulerian–Lagrangian), with realistic parameters for turbulent bedload transport: $\unicode[STIX]{x1D70C}^{p}/\unicode[STIX]{x1D70C}^{f}-1=1.5$ and $d=6$  mm. The black crosses and the red line represent the data and the fit of Meyer-Peter & Müller (Reference Meyer-Peter and Müller1948) respectively.

Figure 9.

Solid depth profile comparison between the two-phase continuous model (Euler–Euler) and the coupled fluid–discrete-element model (Euler–Lagrange) for three different Shields numbers $\unicode[STIX]{x1D703}^{\ast }$ with $d=6$  mm and $\unicode[STIX]{x1D70C}^{p}/\unicode[STIX]{x1D70C}^{f}-1=1.5$ . The solid velocity $\langle v_{x}^{p}\rangle ^{s}$ and transport rate density $\langle q_{s}\rangle ^{s}$ are expressed in $\text{m}~\text{s}^{-1}$ while the solid volume fraction is dimensionless. The fluid mechanics convention is used.

Figure 8 presents the dimensionless sediment transport rate versus the Shields number for the Eulerian–Eulerian model, the Eulerian–Lagrangian model and the experimental data and empirical law of Meyer-Peter & Müller (Reference Meyer-Peter and Müller1948). A good agreement with the experimental data is obtained using both models. The proposed granular rheology enables the Eulerian–Eulerian model to accurately reproduce the critical Shields number ( $\unicode[STIX]{x1D703}_{c}^{\ast }\sim 0.04$ ) and the global trend over a wide range of Shields numbers ( $\unicode[STIX]{x1D703}^{\ast }\in [\unicode[STIX]{x1D703}_{c}^{\ast },0.7]$ ). Quite surprisingly, the Eulerian–Eulerian model better predicts the smooth transition around the critical Shields number, in a region where the continuous assumption breaks down. Conversely, at high Shields numbers the increased importance of the dilute region, which is not well described by the proposed rheology (§ 4.1), induces a slight discrepancy between the Eulerian–Eulerian and Eulerian–Lagrangian models.

To go further in the analysis, the comparison is extended to the granular depth profiles by decomposing the sediment transport rate into the integral of the sediment transport rate density $\langle q_{s}\rangle (z)=\langle v_{x}\rangle ^{s}(z)\unicode[STIX]{x1D719}(z)$ (Bagnold Reference Bagnold1956):

Figure 9 shows the comparison of the solid volume fraction, the solid velocity and the transport rate density profiles obtained with the Eulerian–Eulerian and Eulerian–Lagrangian models for three different Shields numbers $\unicode[STIX]{x1D703}^{\ast }=\{0.05;0.1;0.2\}$ . The Eulerian–Eulerian model reproduces accurately the granular depth profiles obtained with the Eulerian–Lagrangian model for solid volume fractions approximately higher than 0.3 ( $\unicode[STIX]{x1D719}\gtrsim 0.3$ ). This value roughly corresponds to the domain for which the proposed granular rheology describes the Eulerian–Lagrangian results very well (see figure 7). Only the lower quasistatic velocity profile – negligible in terms of sediment transport rate – is not reproduced, as non-local effects have not been taken into account in the granular rheology. In the dilute region ( $\unicode[STIX]{x1D719}\lesssim 0.3$ ), the particle velocity is overestimated by the Eulerian–Eulerian model for the three different Shields numbers, showing an increasing discrepancy towards the lower solid volume fraction region. This is associated with an underestimation of the solid volume fraction close to incipient motion and an overestimation at high Shield numbers. These behaviours point out both the ill-posedness of the $\unicode[STIX]{x1D707}(I)$ rheology as the solid volume fraction goes to zero and the inability of the proposed rheology to describe the dilute region accurately. Despite these limitations, the good lower boundary condition given by the accurate description of the dense granular flow constrains the dilute behaviour and enables a good description of the sediment transport rate over a wide range of Shields numbers with the Eulerian–Eulerian model.

The present results show that a $\unicode[STIX]{x1D707}(I)$ rheology is able to reproduce the main features of turbulent bedload transport in terms of sediment transport rate and dense granular depth profiles, in steady uniform configurations. Considering the much lower computational cost of Eulerian–Eulerian simulations compared with Eulerian–Lagrangian ones, it will allow an important step to be made in the upscaling process towards applications. In this perspective, the fluid phase description might have to be upgraded, using a two-equation turbulence model such as $k$ – $\unicode[STIX]{x1D716}$ or $k$ – $\unicode[STIX]{x1D714}$ (Hsu et al. Reference Hsu, Jenkins and Liu2004; Amoudry Reference Amoudry2014; Lee, Low & Chiew Reference Lee, Low and Chiew2016) or eddy resolving simulations.