From discrete to continuum description of weakly inertial bedload transport

Granular bed motion induced by a liquid shear flow is considered as a model of bedload. This flow is characterized by a localization of the granular flow close to the upper surface of the bed. In such a framework, the paper aims at discussing and proposing continuum effective models describing the physical properties, as the rheology, of the two-phase granular suspension at a so-called mesoscale larger than the grain size. The main questions addressed here concern the relevance and the extension of effective rheological models for granular suspensions often extracted from simple-shear flow configurations in a viscous limit, to situations of both granular shear localization as in bedload and weakly inertial flow regime. For this purpose, we consider the steady transport of a granular bed by a laminar Couette fluid flow above it and close to the onset of motion, for a grain-to-fluid density ratio of 2.5, as silica in water, and a range of particle Reynolds numbers ${\text{Re}}_p\in [0.1,10]$ and Shields numbers $\theta \in [0.1,0.7]$. To provide accurate continuum models, the dynamics into the granular shear layer has to be first known down to a scale smaller than each grain, a so-called microscale. Numerical simulations are thus performed at the microscale at which individual grain dynamics is resolved, using an immersed boundary method (IBM) coupled to a discrete element method (DEM) granular solver. Upscaling is then performed to obtain the equivalent momentum balance at the mesoscale, characterized by continuum phases, using a spatial averaging method on volume element larger than the grain diameter. This approach allows us obtaining stresses, strains, and their relationships for the fluid phase, the granular phase, and the equivalent mixture, independently. The main contribution of this work is threefold: (i) we highlight the relevance of mesoscopic rheological continuum law for localized granular shear flow; (ii) we extract rheological models from direct numerical simulations (IBM/DEM) in a weakly inertial regime, going beyond purely viscous situations; and (iii) we extend Coulomb-like model $\mu \left(I\right)$ of a granular suspension to incorporate fluid/particle inertial effects showing a different dependence of fluid phase and granular phase contributions with dimensionless numbers.

Figure 1

(a) Scheme of the considered problem at the microscale. (b) Domain of simulation. (c) Scheme of the considered problem at the mesoscale. (d) Vertical profile of the streamwise velocity in the fluid phase (blue, ${\langle u\rangle }^f$), the granular phase (red, ${\langle u\rangle }^p$), and solid volume fraction (black, $\varphi$) for ${\text{Re}}_p=1$ and $\theta =0.67$. Note that the velocity profile of the granular phase is represented with a solid line for $\varphi ⩾0.025$ and with a dashed line otherwise. Here, we use a characteristic length of the weighting function $h_g=d$. $h_{\text{bed}}$ is the mean height of the granular bed defined as the barycentre of the local granular flow rate $h_{\text{bed}}=\int ydq/\int dq=\int y\varphi u_{p}dy/\int \varphi u_{p}dy$ as suggested by Durán et al. (2012) [32].

Figure 2

Time evolution of the instantaneous granular flow rate per unit depth scaled by $V_s$ with $V_s=({\rho }_p-{\rho }_f)g{d}^2/\left(18{\eta }_f\right)$ the Stokes velocity (${\text{Re}}_p=1, \theta =0.67, \text{Ga}=1.22$). Black line: IBM-DEM results. Red line: mean value at steady state. Blue line: $q_{\text{conv}}$ as defined using Eq. (8) with $t_{\text{max}}=1300{\dot{\mathrm{\Gamma }}}^{-1}$).

Figure 3

Vertical profiles of the vertical projection of the equivalent-fluid momentum balance at the mesoscale (42) and (43), for ${\text{Re}}_p=1, \theta =0.67, \text{Ga}=1.22$. (a) Vertical profiles of all the terms in Eq. (42), labeled from (1) to (4). The black dotted line corresponds to the “balance error” profile ${\zeta }_b^{f,y}$ (see text for details). (b) Vertical profiles of all the terms in Eq. (43), namely ${\mathrm{\Sigma }}_{yy}^f$ (blue), $\epsilon {\langle {\sigma }_{yy}\rangle }^f$ (cyan), $n{\langle s_{yy}^f\rangle }^p$ (magenta), and $-{\rho }_f\epsilon {\langle u_y^{\prime 2}\rangle }^f$ (yellow). The black dotted line is ${\mathrm{\Sigma }}_{yy}^f$ obtained by vertically integrating Eq. (42), indicating the “integral cumulative” error ${\zeta }_{ic}^{f,y}$.

Figure 4

Vertical profiles of the horizontal projection of the equivalent-fluid momentum balance at the mesoscale (44) and (45), for ${\text{Re}}_p=1, \theta =0.67, \text{Ga}=1.22$. (a) Vertical profiles of all the terms in Eq. (44), labeled from (1) to (3). The black dotted line corresponds to the “balance error” profile ${\zeta }_b^{f,x}$ (see text for details). (b) Vertical profiles of all the terms in Eq. (45), namely ${\mathrm{\Sigma }}_{xy}^f$ (blue), $\epsilon {\langle {\sigma }_{xy}\rangle }^f$ (cyan), $n{\langle s_{xy}^f\rangle }^p$ (magenta), and $-{\rho }_f\epsilon {\langle u_x^{\prime }{u}_y^{\prime }\rangle }^f$ (yellow). The black dotted line is ${\mathrm{\Sigma }}_{xy}^f$ obtained by vertically integrating Eq. (44), indicating the “integral cumulative” error ${\zeta }_{ic}^{f,x}$.

Figure 5

Vertical profiles of the vertical projection of the equivalent-granular momentum balance at the mesoscale (46) and (47), for ${\text{Re}}_p=1, \theta =0.67, \text{Ga}=1.22$. (a) Vertical profiles of all the terms in Eq. (46), labeled from (1) to (5). The black dotted line corresponds to the “balance error” profile ${\zeta }_b^{p,y}$ (see text for details). (b) Vertical profiles of all the terms in Eq. (47), namely ${\mathrm{\Sigma }}_{yy}^p$ (red), $n{\langle s_{yy}^p\rangle }^p$ (magenta), and $-{\rho }_p\varphi {\langle u_{p|y}^{\prime 2}\rangle }^f$ (yellow). Note that here ${\mathrm{\Sigma }}_{yy}^p\approx n{\langle s_{yy}^p\rangle }^p$. The black dotted line is ${\mathrm{\Sigma }}_{yy}^p$ obtained by vertically integrating Eq. (46), indicating the “integral cumulative” error ${\zeta }_{ic}^{p,y}$.

Figure 6

Vertical profiles of the horizontal projection of the equivalent-granular momentum balance at the mesoscale (48) and (49), for ${\text{Re}}_p=1, \theta =0.67, \text{Ga}=1.22$. (a) Vertical profiles of all the terms in Eq. (48), labeled from (1) to (4). The black dotted line corresponds to the “balance error” profile ${\zeta }_b^{p,x}$ (see text for details). (b) Vertical profiles of all the terms in Eq. (49), namely ${\mathrm{\Sigma }}_{xy}^p$ (red), $n{\langle s_{yy}^p\rangle }^p$ (magenta), and $-{\rho }_p\varphi {\langle u_{p|y}^{\prime 2}\rangle }^f$ (yellow). Note that here ${\mathrm{\Sigma }}_{xy}^p\approx n{\langle s_{xy}^p\rangle }^p$. The black dotted line is ${\mathrm{\Sigma }}_{xy}^p$ obtained by vertically integrating (48), indicating the “integral cumulative” error ${\zeta }_{ic}^{p,x}$.

Figure 7

Apparent viscosity ${\eta }_{\text{eq}}^f/{\eta }_f$ of the fluid phase vs solid volume fraction $\varphi$: ($•$), present IBM-DEM simulation for ${\text{Re}}_p=1, \theta =0.67, \text{Ga}=1.22$ [Eq. (50)]; (dashed black) Einstein's model [8, 9] ${\eta }_{\text{eq}}^f/{\eta }_f=1+\frac{5}{2}\varphi$; (solid black) Gibilaro et al. (2007)'s correlation [11] ${\eta }_{\text{eq}}^f/{\eta }_f={(1-\varphi )}^{-\beta }$ with $\beta =2.8$; and (red) adjustment with ${\eta }_{\text{eq}}^f/{\eta }_f={(1-\varphi )}^{-1.71}$. Confidence intervals of 95% are computed using the geometric standard deviation of the apparent viscosity for each uncorrelated time sample.

Figure 8

(a) Apparent viscosity ${\eta }_{\text{eq}}^f/{\eta }_f$ of the fluid phase vs solid volume fraction $\varphi$ for all the cases (${\text{Re}}_p, \theta$, Ga) considered in the present work: see Table 1 for detail of the various symbols. Green solid lines show the obtained fit of the constitutive law ${\eta }_{\text{eq}}^f/{\eta }_f={(1-\varphi )}^{-\beta }$ for $({\text{Re}}_p,\theta ,\text{Ga})=(0.1,0.67,0.22)$ (black crosses) and $({\text{Re}}_p,\theta ,\text{Ga})=(10,0.67,3.9)$ (red diamonds), respectively. Green dashed lines in the insert are the small-$\varphi$ expansion (52). (b) Corresponding exponent $\beta$ as a function of $\theta$. (c) $\beta -{\beta }_0$ as a function of $r\theta \sqrt{{\text{Re}}_p}$ with ${\beta }_0=0.565$. Lines correspond to $2/3$ power-law functions with scale factor 1 (dotted line) and 0.85 (dashed line).

Figure 9

Vertical distribution of the normal components of the equivalent granular phase (${\text{Re}}_p=1, \theta =0.67, \text{Ga}=1.22$): (solid black) ${\mathrm{\Sigma }}_{yy}^p$; (blue) ${\mathrm{\Sigma }}_{xx}^p$; (red) ${\mathrm{\Sigma }}_{zz}^p$; and (dashed black) $P^p=({\mathrm{\Sigma }}_{xx}^p+{\mathrm{\Sigma }}_{yy}^p+{\mathrm{\Sigma }}_{zz}^p)/3$.

Figure 10

Constitutive law of the equivalent granular phase (${\text{Re}}_p=1, \theta =0.67, \text{Ga}=1.22$). (a) Solid volume fraction $\varphi$ versus $\mathcal{K}$, (b) apparent friction coefficient ${\mu }_{\text{eq}}^p$ versus $\mathcal{K}$, (c) apparent friction coefficient ${\mu }_{\text{eq}}^p$ versus $\varphi$: ($•$) present IBM-DEM simulation; [red solid lines in panels (a)–(c)] $\varphi \left(\mathcal{K}\right), {\mu }_{\text{eq}}^p\left(\mathcal{K}\right)$, and ${\mu }_{\text{eq}}^p\left(\varphi \right)$ given by Eqs. (37), (38), and (39), respectively, with ${\varphi }_c=0.615, b=0.7, {\mu }_1=0.23, {\mu }_2=0.8$, and $\sqrt{{\mathcal{K}}_0}=0.15$; (black thin dash line) in panel (a), same as red line but for $b=1$; (green dashed lines) in panels (b), (c), same as red line but for fitting parameters ${\mu }_2$ and $\sqrt{{\mathcal{K}}_0}$ as in Fig. 11. Confidence intervals of $95\%$ are computed using the geometric standard deviation of the solid fraction (a) or the apparent friction coefficient (b) for each uncorrelated time sample.

Figure 11

Constitutive laws $\varphi \left(\mathcal{K}\right)$ (a), ${\mu }_{\text{eq}}^p\left(\mathcal{K}\right)$ (b), and ${\mu }_{\text{eq}}^p\left(\varphi \right)$ (c) of the equivalent granular phase for all the cases (${\text{Re}}_p, \theta$, Ga) considered in the present work: see Table 1 for detail of the various symbols and Fig. 10 for legend. Green solid lines are Eqs. (37), (38), and (39), respectively, with ${\varphi }_c=0.615, b=0.7, {\mu }_1=0.23, {\mu }_2=0.75$, and $\sqrt{{\mathcal{K}}_0}=0.18$. Note that the cases corresponding to the smaller Shield numbers considered here, namely $\theta \le 0.22$, are not shown because they are too close to the threshold of motion of the granular phase.

Figure 12

Apparent viscosity ${\eta }_{\text{eq}}^m/{\eta }_f$ of the equivalent mixture phase vs solid volume fraction $\varphi$: ($•$) present IBM-DEM simulation for ${\text{Re}}_p=1, \theta =0.67, \text{Ga}=1.22$ [Eq. (54)]; (dotted black) Einstein's model [8, 9] ${\eta }_{\text{eq}}^m/{\eta }_f=1+\frac{5}{2}\varphi$; (dashed black) Batchelor and Green's model [10] ${\eta }_{\text{eq}}^m/{\eta }_f=1+\frac{5}{2}\varphi +7.6{\varphi }^2$; (dashed red) Krieger and Dougherty's model [49] ${\eta }_{\text{eq}}^m/{\eta }_f={(1-\varphi /{\varphi }_m)}^{-\frac{5}{2}{\varphi }_m}$ with ${\varphi }_m=0.615$; (dashed blue) Boyer et al.'s model [12] ${\eta }_{\text{eq}}^m/{\eta }_f=1+\frac{5}{2}\varphi {(1-\varphi /{\varphi }_c)}^{-1}+{\mu }_{\text{eq}}^p\left(\varphi \right){[\varphi /({\varphi }_c-\varphi )]}^2$ using fittings parameters as obtained for the granular phase (see legend of Fig. 11). Insert: linear representation. Confidence intervals of 95% are computed using the geometric standard deviation of the apparent viscosity for each uncorrelated time sample. Insert: comparison with the resolved simulations of Vowinckel et al. (2021) [26] (colored circles, see the legend of Fig. 13 for key).

Figure 13

(a) Apparent friction coefficient for the mixture phase ${\mu }_{\text{eq}}^m$ as a function of the viscous number $\mathcal{J}$: a comparison between the experiments of Houssais et al. (2016) [16] here denoted H16 (blue dots), the Euler-Lagrange simulations of Pahtz and Duran (2018) [20] here denoted PD18 (gray triangles), the resolved simulations of Vowinckel et al. (2021) [26] here denoted V21 (gray dots) and the present simulations (black diamonds). The dot line corresponds to the apparent friction coefficient of a suspension flow proposed by Ref. [12], ${\mu }^{\text{BGP}}=\mathcal{J}+\frac{5}{2}{\varphi }_m\sqrt{\mathcal{J}}$, with ${\varphi }_m=0.585$. Red line corresponds to the combination of effective friction models for granular and fluid phases, ${\mathcal{F}}^{-1}{\eta }_{\text{eq}}^f/{\eta }_f+{\mu }_{\text{eq}}^p$, as developed in the present paper, using Eqs. (32) and (51) to model the equivalent fluid phase, and Eqs. (37) and (38) and fitting parameters from Fig. 10 to model the equivalent granular phase. ${\mathcal{F}}^{-1}$ is computed using Eq. (57). (b) Constitutive law ${\mu }_{\text{eq}}^m\left(\mathcal{K}\right)$ of the equivalent mixture phase for all the cases (${\text{Re}}_p, \theta$, Ga) considered in the present work: see Table 1 for detail of the various symbols and Fig. 10 for legend.

Figure 14

(a) Apparent viscosity ${\eta }_{\text{eq}}^m/{\eta }_f$ of the mixture phase vs solid volume fraction $\varphi$ for all the cases considered in the present work: see Table 1 for detail of the various symbols (solid gray line correspond to Einstein's model (28) at small $\varphi$). (b) Values of ${\varphi }_m$ from a best fit of Krieger and Dougherty's model [49] ${\eta }_{\text{eq}}^m/{\eta }_f={(1-\varphi /{\varphi }_m)}^{-\frac{5}{2}{\varphi }_m}$; ($▾$), ${\text{Re}}_p=0.1$; ($▸$), ${\text{Re}}_p=1$; ($▴$), ${\text{Re}}_p=10$; for all $\theta$ with the same markers as in Fig. 8. The dashed line indicates the corresponding mean value: ${\varphi }_m\simeq 0.622$.

Figure 15

$({\eta }_{\text{eq}}^p/{\eta }_f)/{\mu }_{\text{eq}}^p\equiv \mathcal{F}$ as a function of $\varphi$ for all the cases considered in the present work: blue solid line corresponds to model (56) for ${\text{Re}}_p=0$ and green solid line is Eq. (57) (see text for details).

Figure 16

Comparison of the predicted apparent viscosity of the mixture versus that obtained from present simulations ${\eta }_{\text{eq}}^m/{\eta }_f$: (a) Krieger and Dougherty (1959) ${\eta }_{\text{eq}}^{KD}/{\eta }_f$; (b) Boyer et al. (2011) ${\eta }_{\text{eq}}^{\text{BGP}}/{\eta }_f$; (c) model (55) and (56) ${\eta }_{\text{eq}}^F/{\eta }_f$; (c) model (55, 56, 57) ${\eta }_{\text{eq}}^{FLBC}/{\eta }_f$: see Table 1 for detail of symbols. The two dashed lines represents the discrepancy ratio of 2 and 0.5, respectively.

Figure 17

Vertical profiles of streamwise granular velocity (filled symbols), fluid velocity (open symbols), and solid volume fraction (solid lines). ($•$) case A; ($\blacksquare$) case B; ($▴$) case C; ($⧫$) case D.

Figure 18

Effect of the spatial resolution and size of the domain on the constitutive laws of the equivalent (a) fluid phase, (b), (c) granular phase, and (d), (e) mixture phase: (black) case A; (magenta) case B; (red) case C; (blue) case D.

Figure 19

Comparison of vertical distribution of all the terms of the momentum balance Eqs. (42) and (44) of the equivalent granular phase in the (a) vertical direction, (b) streamwise direction, for the cases A (pale) and B (bright). Same legend as Figs. 5 and 6. See Eqs. (46) and (48).

Figure 20

(a) Velocity profile of the fluid phase (dashed) and the granular phase (solid) for various of $h_g=\{d/10,d/4,d/2,d\}$. Velocity is scaled by $V_s$ with $V_s=({\rho }_p-{\rho }_f)g{d}^2/\left(18{\eta }_f\right)$ the Stokes velocity. (b) Solid volume fraction for the same values of $h_g$. Here, ${\text{Re}}_p=1$ and $\theta =0.67$ ($\text{Ga}=1.22$) and we use the setup of case A of Table 3.

Figure 21

Apparent rheology of the mesoscopic fluid phase (top row), granular phase (second and third rows) and mixture phase (bottom row): (a, blue) $h_g=d/2$, (b, magenta) $h_g=d/4$; (c, red) $h_g=d/10$. Results for $h_g=d$ are also plotted as reference (black). For the apparent viscosity of the fluid phase, the model ${\eta }_{\text{eq}}^f/{\eta }_f={(1-\varphi )}^{-1.71}$ is also plotted (solid black line). Here, ${\text{Re}}_p=1$ and $\theta =0.67$ ($\text{Ga}=1.22$) and we use the setup of case A of Table 3.

Figure 22

$\mathcal{F}\equiv P^p/\left({\eta }_f\dot{\gamma }\right)$ as a function of $\varphi$ for $h_g=d$ (black symbols), $h_g=d/2$ (blue symbols), $h_g=d/4$ (magenta symbols), $h_g=d/10$ (red symbols). As $d/4$ and $d/10$ are shown to be dispersed (see Fig. 21), $h_g=d$ and $h_g=d/2$ symbols are slightly larger for clarity. Here, ${\text{Re}}_p=1$ and $\theta =0.67$ ($\text{Ga}=1.22$) and we use the setup of case A of Table 3. Dash line and solid are reported from Fig. 15, i.e., from Ref. [12] and the extended model (57).

Figure 23

Vertical profile of the drag force on the grains compared with Richardson and Zaki's drag law [58] given in Eq. (C2) (${\text{Re}}_p=1, \theta =0.67, \text{Ga}=1.22$). Left: $\mathrm{\Delta }{U}^{*}$ vs $F^{*}/{\mathrm{\Phi }}^{*}$. Right: $F^{*}/\mathrm{\Delta }{U}^{*}$ vs ${\mathrm{\Phi }}^{*}$ with $F^{*}=\frac{{\mathbf{F}}_D}{{\eta }_f\dot{\mathrm{\Gamma }}/d}, \mathrm{\Delta }{U}^{*}=\frac{{\langle \mathbf{u}\rangle }^f-{\langle \mathbf{u}\rangle }^p}{\dot{\mathrm{\Gamma }}d}$ and ${\mathrm{\Phi }}^{*}=\frac{18\varphi }{{(1-\varphi )}^{2.65}}$. Each row corresponds to one value of the characteristic length of the weighting function $h_g$, namely (from top to bottom) $h_g=\{d,d/2,d/4,d/10\}$.

Figure 24

Effect of including a lubrication force on the phases rheologies: (a) Apparent viscosity ${\eta }_{\text{eq}}^f/{\eta }_f$ of the equivalent fluid phase vs solid volume fraction $\varphi$. (b) Apparent friction coefficient ${\mu }_{\text{eq}}^p$ of the equivalent granular phase vs $\mathcal{J}$. (c) Apparent viscosity ${\eta }_{\text{eq}}^m/{\eta }_f$ of the equivalent mixture phase vs solid volume fraction $\varphi$. (d) Apparent friction coefficient ${\mu }_{\text{eq}}^m$ of the equivalent mixture phase vs $\mathcal{J}$: present IBM-DEM simulations for ($◊$) ${\text{Re}}_p=1, \theta =0.67, \text{Ga}=1.2$ and ($*$) ${\text{Re}}_p=1, \theta =0.33, \text{Ga}=1.7$, with a lubrication force (small symbols) and without lubrication force (large symbols).