2.1. The FDM: a summary

The FDM has been detailed in the past (Gallier et al. Reference Gallier, Lemaire, Lobry and Peters2014a), including the specific version that is used in the present article (Orsi et al. Reference Orsi, Lobry and Peters2023). The main features of the method are recalled here. In the FDM, the flow of a Newtonian liquid, with dynamic viscosity $\eta$ and density $\rho$, is computed in the whole meshed simulation domain $\mathcal {D}$ including the particles. The particles are rigid spheres with the same density as the liquid. The rigidity of the particles is enforced using a force density $\rho \boldsymbol{\lambda}$ that acts only in the volume of each particle $\mathcal {D}_p$. In the present work, inertia is neglected for both fluid and particle dynamics. As a consequence, fluid velocity $\boldsymbol {u}$ and pressure $P$ obey modified Stokes equations

(2.1)$$\begin{gather} \left.\begin{array}{c} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\sigma} + \rho \boldsymbol{\lambda}= \boldsymbol{0},\\ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol u = 0, \end{array}\right\} \end{gather}$$

(2.2)$$\begin{gather}\boldsymbol{\sigma} ={-} P\boldsymbol{\delta} + \eta ( \boldsymbol{\nabla} \boldsymbol u + \boldsymbol{\nabla} \boldsymbol u^T), \end{gather}$$

constraining the fluid inside each particle ($\kern 1.5pt p$) to undergo rigid body motion,

(2.3)\begin{equation} \boldsymbol u(\boldsymbol x) = \boldsymbol{U}_p + \boldsymbol{\varOmega}_p\times (\boldsymbol x - \boldsymbol x_p ) \quad \text{for} \ \boldsymbol x \in \mathcal{D}_p, \end{equation}

where $\boldsymbol {U_p}$ and $\boldsymbol{\varOmega} _p$ are the particle translational and rotational velocities, respectively. In addition, various boundary conditions for the pressure and the velocity may be applied, depending on the considered specific flow. They are given in § 2.3 for the channel flow of interest in the present article.

In the present article, the particles are only subjected to hydrodynamic and contact forces. Then, Newton's equations write for each particle $(p)$,

(2.4)\begin{equation} \left.\begin{gathered} \boldsymbol{F}^h_p + \boldsymbol{F}^c_p = \boldsymbol{0},\\ \boldsymbol{T}^h_p + \boldsymbol{T}^c_p = \boldsymbol{0}, \end{gathered}\right\} \end{equation}

where $\boldsymbol{F}^h_p$ and $\boldsymbol{F}^c_p$, respectively, stand for the hydrodynamic and contact forces exerted on the particle $(p)$, and $\boldsymbol{T}^h_p$ and $\boldsymbol{T}^c_p$ are the corresponding torques. Contact forces and torques will be briefly described in the next section. Hydrodynamic forces and torques are split into two contributions. The first part corresponds to the moments that originate in the fluid flow computed by the solver, and are hereafter denoted by $\boldsymbol{F}^{FDM}_p$ and $\boldsymbol{T}^{FDM}_p$. These are computed as follows:

(2.5)\begin{equation} \left.\begin{gathered} \boldsymbol{F}^{FDM}_p={-} \rho \int_{\mathcal{D}_p} \boldsymbol{\lambda}\, \mathrm{d}\mathcal{V},\\ \boldsymbol{T}^{FDM}_p ={-} \rho \int_{\mathcal{D}_p} (\boldsymbol x-\boldsymbol x_p )\times \boldsymbol{\lambda}\, \mathrm{d}\mathcal{V}. \end{gathered}\right\} \end{equation}

However, the solver cannot properly compute the flow between particles when their surfaces are closer than the mesh size. As a consequence, the stress moments applied to the particles must be corrected using theoretical expressions. The method we follow in the present article has been fully explained and validated by Orsi et al. (Reference Orsi, Lobry and Peters2023). The generalized particle velocity and subgrid force vectors are, respectively, denoted by $\boldsymbol {\mathcal {U}} = ( \boldsymbol{U}_1, \boldsymbol{U}_2, \ldots, \boldsymbol{\varOmega} _1, \boldsymbol{\varOmega} _2, \ldots )$ and $\boldsymbol {\mathcal {F}}^{SG} = ( \boldsymbol{F}_1^{SG}, \boldsymbol{F}_2^{SG}, \ldots, \boldsymbol{T}_1^{SG}, \boldsymbol{T}_2^{SG}, \ldots )$. The latter may be computed from the former using the subgrid resistance matrix $\boldsymbol {\mathcal {R}}^{SG}_{\mathcal {FU}}$,

(2.6)\begin{equation} \boldsymbol{\mathcal{F}}^{SG} ={-} \boldsymbol{\mathcal{R}}^{SG}_{\mathcal{FU}} \boldsymbol{\cdot} \boldsymbol{\mathcal{U}}, \end{equation}

where the subgrid resistance matrix depends on the fluid viscosity, particle positions and radii, and it is computed pairwise as the difference between the theoretical two-particle resistance matrix and its counterpart as measured by the uncorrected FDM solver. Finally, the total hydrodynamic force and torque read

(2.7)\begin{equation} \left.\begin{gathered} \boldsymbol{F}^h_p = \boldsymbol{F}^{FDM}_p + \boldsymbol{F}^{SG}_p, \\ \boldsymbol{T}^h_p = \boldsymbol{T}^{FDM}_p + \boldsymbol{T}^{SG}_p. \end{gathered}\right\} \end{equation}

Equation (2.6) deserves several comments. First, unlike the standard subgrid corrections, no mention of the ambient velocity field is made, meaning that only the particle velocities are needed to compute the subgrid forces. Such an expression has been fully justified by Orsi et al. (Reference Orsi, Lobry and Peters2023). It stems from the fact that only lubrication flows that involve relative velocities of particle surfaces are missed by the fluid solver, and consequently have to be considered for corrections. This formulation is especially well suited to the present flow, where the velocity profile depends on the solid volume fraction profile, which evolves in time and is not given in advance. It should also be noted that the same type of corrections has to be applied for particles close to bounding walls when particle–wall interactions have to be considered (Orsi et al. Reference Orsi, Lobry and Peters2023).

Equations (2.1)–(2.6) must be solved at each time step. The unknown quantities are the pressure and velocity fields, ($P, \boldsymbol {u}$), the force density, $\rho \boldsymbol{\lambda}$, and the particle velocities, $\boldsymbol{U}_p$ and $\boldsymbol{\varOmega} _p$. The main issue lies in the determination of the force density that allows proper accounting for the constraints imposed on the particles, namely particle rigidity – equation (2.3) – and particle momentum balance – (2.4) together with (2.5)–(2.7). The numerical method that is followed to solve this problem has been explained in detail by Orsi et al. (Reference Orsi, Lobry and Peters2023). In particular, it has been implemented in the frame of the OpenFOAM toolbox, which easily allows the solving of (2.1) and (2.2) using high-performance computing facilities taking advantage of the parallelization.

2.1.1. Suspension rheology: dipole moment

To compute the effective stress in a suspension of solid particles, the first moment of the hydrodynamic surface stress acting on each particle $\boldsymbol {D}_p$ (Batchelor Reference Batchelor1970; Lhuillier Reference Lhuillier2009; Nott et al. Reference Nott, Guazzelli and Pouliquen2011) is needed,

(2.8)\begin{align} \boldsymbol{\mathsf{D}}^{FDM}_p &= \int_{\partial \mathcal{D}_p} \boldsymbol{\sigma} \boldsymbol{\cdot} \boldsymbol n \otimes (\boldsymbol x -\boldsymbol x_p)\, \mathrm{d}\mathcal{S}\nonumber\\ &={-}\rho \int_{\mathcal{D}_p} \boldsymbol{\lambda} \otimes (\boldsymbol x - \boldsymbol x_p)\, \mathrm{d}\mathcal{V} -\left(\int_{\mathcal{D}_p} P\, \mathrm{d}\mathcal{V}\right) \boldsymbol{\delta}. \end{align}

The second equality is obtained using (2.1) and (2.2) (Orsi et al. Reference Orsi, Lobry and Peters2023).

The antisymmetric part of the force dipole $\boldsymbol {D}^{FDM}_p$ is a tensor whose components are the components of the hydrodynamic torque $\boldsymbol {T}_p^{FDM}$ exerted on particle $(p)$, while the symmetric part is usually called stresslet. The latter is split into the deviatoric stresslet $\boldsymbol {S}^{FDM}_p$ and the isotropic part $s^{FDM}_p/3\ \boldsymbol{\delta}$, where $s^{FDM}_p$ stands for the trace of the force dipole. As a consequence, the FDM traceless stresslet, stresslet trace and force dipole write as

(2.9) \begin{equation} \left.\begin{gathered} \boldsymbol{\mathsf{S}}_p^{FDM}={-}\rho \int_{\mathcal{D}_p} \left\lbrace \tfrac{1}{2}[ \boldsymbol{\lambda} \otimes(\boldsymbol x -\boldsymbol x_p) +(\boldsymbol x -\boldsymbol x_p) \otimes \boldsymbol{\lambda} ] -\tfrac{1}{3}\boldsymbol{\lambda} \boldsymbol{\cdot}(\boldsymbol x -\boldsymbol x_p)\boldsymbol{\delta} \right\} \mathrm{d}\mathcal{V},\\ s_p^{FDM} ={-}\int_{\mathcal{D}_p}[ \rho \boldsymbol{\lambda} \boldsymbol{\cdot}(\boldsymbol x -\boldsymbol x_p) +3 P ]\, \mathrm{d}\mathcal{V},\\ \boldsymbol{\mathsf{D}}^{FDM}_p = \boldsymbol{\mathsf{S}}_p^{FDM} +\dfrac{s_p^{FDM}}{3}\boldsymbol{\delta} - \dfrac{1}{2}\boldsymbol{\epsilon} \boldsymbol{\cdot} \boldsymbol{T}_p^{FDM}. \end{gathered}\right\} \end{equation}

We note here that a reference has to be defined for the pressure. It has been chosen as the mean pressure over the whole simulation volume, meaning that $\int _D P\, \mathrm {d}\mathcal {V}=0$ at each time step.

As in the case of force and torque, a subgrid correction must be added to the FDM stresslet on the particles, which is computed using specific resistance subgrid matrices (Orsi et al. Reference Orsi, Lobry and Peters2023). The total hydrodynamic traceless stresslet, stresslet trace and force dipole are then computed as

(2.10)\begin{equation} \left.\begin{gathered} \boldsymbol{\mathsf{S}}_p^{h} = \boldsymbol{S}_p^{FDM} + \boldsymbol{S}_p^{SG},\\ s_p^{h} = s_p^{FDM} + s_p^{SG},\\ \boldsymbol{\mathsf{D}}_p^{h} = \boldsymbol{D}_p^{FDM} + \boldsymbol{D}_p^{SG}. \end{gathered}\right\} \end{equation}

2.2. Contact forces

In the present article, contact forces reflect both particle roughness and friction. They are described in detail in Gallier et al. (Reference Gallier, Lemaire, Peters and Lobry2014b), Peters et al. (Reference Peters, Ghigliotti, Gallier, Blanc, Lemaire and Lobry2016) and Orsi et al. (Reference Orsi, Lobry and Peters2023), and we give here only a short account of them. Contact between particles $(p)$ and $(q)$ occurs as soon as the distance between their surfaces is smaller than the roughness height $h_r$, resulting in elastic forces exerted by particle $(q)$ on particle $(p)$, both normal $(\boldsymbol {F}^{c,n}_{pq})$ and tangential $(\boldsymbol {F}^{c,t}_{pq})$. Sliding is triggered as soon as the tangential force reaches the value $\Vert \boldsymbol {F}^{c,t}_{pq} \Vert = \mu _s \Vert \boldsymbol {F}^{c,n}_{pq} \Vert$, where $\mu _s$ stands for the static friction coefficient. In the sliding regime, the latter relation between normal and tangential forces holds, meaning that no difference is made here between the static and dynamic friction coefficients. The contact torques on particles $(p)$ and $(q)$ read

(2.11)\begin{equation} \boldsymbol{T}^{c}_{pq} = \dfrac{a_{p}}{a_{p}+a_{q}} \boldsymbol{x}_{pq} \times \boldsymbol{F}^{c,t}_{pq} = \dfrac{a_p}{a_q}\boldsymbol{T}^{c}_{qp}, \end{equation}

where $a_p$ and $a_q$ stand for the particle radii, and $\boldsymbol x_{pq} = \boldsymbol x_q - \boldsymbol x_p$ is the distance between the particle centres. Contact forces also induce an additional stresslet which is given for the contacting particles $(p)$ and $(q)$ by

(2.12)\begin{equation} \boldsymbol{\mathsf{S}}^{c}_{pq} = \dfrac{1}{2} \dfrac{a_{p}}{a_{p}+a_{q}} ( \boldsymbol{F}^{c}_{pq} \otimes \boldsymbol{x}_{pq} + \boldsymbol{x}_{pq} \otimes \boldsymbol{F}^{c}_{pq} ) = \dfrac{a_p}{a_q}\boldsymbol{S}^{c}_{qp}. \end{equation}

The induced contact force dipole is written for particles $(p)$ and $(q)$ as follows:

(2.13)\begin{equation} \boldsymbol{\mathsf{D}}^{c}_{pq} = \boldsymbol{\mathsf{S}}^{c}_{pq} - \dfrac{1}{2}\boldsymbol{\epsilon} \boldsymbol{\cdot} \boldsymbol{T}^c_{pq} =\dfrac{a_p}{a_q}\boldsymbol{D}^{c}_{qp}, \end{equation}

where $\boldsymbol {\epsilon }$ stands for the Levi-Civita symbol. The stresslet in (2.12) is not traceless, and its trace determines the contribution of contact forces to particle pressure. Finally, it should be noted that the same type of frictional contact force is exerted by the bounding walls on contacting particles and contributes to the particle force dipole.

2.3. Numerical set-up

We investigate three values of the mean volume fraction $\bar {\phi } = [0.3, 0.4,0.5]$. The suspensions are bidisperse ($a_2 = 1.4a_1$, $\overline {\phi _1} = \overline {\phi _2} = \bar {\phi }/2$) to avoid crystallization. Particles and walls are frictional, with a friction coefficient $\mu _s= 0.5$ and roughness height $h_r = 0.005a_1$. The latter is a value commonly found in experiments with model spherical particles (Blanc et al. Reference Blanc, Peters and Lemaire2011; Tanner & Dai Reference Tanner and Dai2016). The particles are initially placed at random positions in the simulation cell using the same procedure as in Howard et al. (Reference Howard, Dong, Patel, D'Elia, Maxey and Stinis2023). The size of the simulation cell is $(L_x,\, L_y,\, L_z) = (40 a_1,\, 40 a_1,\, 20 a_1)$ (figure 1). The fluid equations are discretized on a fixed Cartesian grid with mesh size $0.2 a_1$. The trade-off between precision and computation speed has been discussed in Gallier et al. (Reference Gallier, Lemaire, Lobry and Peters2014a). The flow is driven by a uniform pressure gradient in the negative $x$ direction so that (2.1) together with (2.2) reads

(2.14)\begin{equation} \eta \Delta \boldsymbol{u} - \boldsymbol{\nabla} p - \dfrac{\mathrm{d}P_0}{\mathrm{d}\kern0.07em x}\boldsymbol{e}_x + \rho \boldsymbol{\lambda} = \boldsymbol{0}. \end{equation}

The plane boundaries normal to the $y$ direction are rigid walls at rest, while the other boundaries are given periodic conditions, for both the pressure perturbation $p$ and the velocity $\boldsymbol {u}$. The particle distribution is made periodic along the same directions. Length, shear rate, time, velocity and stresses are made dimensionless as follows:

(2.15a–e)\begin{equation} y^{*} = \dfrac{y}{L_{y}}, \quad \dot{\gamma}^{*} = \dfrac{\dot{\gamma}}{\overline{\dot{\gamma}}_{0}}, \quad t^{*} = \overline{\dot{\gamma}}_{0}t, \quad \boldsymbol{u}^{*} = \dfrac{\boldsymbol{u}}{\overline{\dot{\gamma}}_{0}L_{y}}, \quad \varSigma^{*} = \dfrac{\varSigma}{\eta\overline{\dot{\gamma}}_{0}}, \end{equation}

where $\bar {\dot {\gamma }}_0$ is the mean shear rate in a pressure-driven flow of pure liquid,

(2.16)\begin{equation} \bar{\dot{\gamma}}_{0} = \dfrac{L_{y}}{4\eta}\left\vert\dfrac{\mathrm{d}P_{0}}{\mathrm{d}\kern0.07em x}\right\vert. \end{equation}

Figure 1. Simulation cell.

In the following, the required equations are given in dimensional form, for the sake of physical clarity. The dimensionless equations may be recovered from the dimensional ones by letting $\eta =1$, $L_y=1$, $\bar {\dot {\gamma }}_{0}=1$, meaning for instance that $\mathrm {d}P_0/\mathrm {d}\kern0.07em x = -4$. However, figures display dimensionless quantities, where the superscript (*) has been dropped, though. All numerical values given in the text also correspond to dimensionless quantities.

As we will see, the actual shear rate in the channel is always lower than $\bar {\dot {\gamma }}_{0}$ due to the enhanced viscosity induced by the particles. The time step $\Delta t$ – meaning $\bar {\dot {\gamma }}_{0}\Delta t$ – is chosen as $\Delta t = 2 \times 10^{-3}$ for all simulations. An estimation of the mean accumulated strain as felt by the fluid is classically computed from the mean distance covered by the fluid (Lyon & Leal Reference Lyon and Leal1998; Howard et al. Reference Howard, Maxey and Gallier2022) as

(2.17)\begin{equation} \bar \gamma(t)= \dfrac{2}{Ly}\int_0^t \left[ \dfrac{1}{V} \int_\mathcal{D} u_x\, \mathrm{d}\mathcal{V} \right] \mathrm{d}t= \dfrac{2 L_t}{Ly}. \end{equation}

The simulation parameters are gathered in table 1.

Table 1. List of simulation parameters: $\bar {\phi }$, mean volume fraction; $N_p$, number of particles; $\Delta t$, time step; $\bar {\gamma }_{total}$, total accumulated strain; $\bar {\gamma }_{steady}$, accumulated strain at steady flow start; $\Delta \bar {\gamma }_{steady}$, accumulated strain interval for steady averaging; $\Delta t_{steady}$, time interval for steady state averaging.

4.1. Local averaging

The procedure yielding the various locally averaged quantities is explained here. To this purpose, we follow the general method proposed by Yeo & Maxey (Reference Yeo and Maxey2011). It is formally different from the method developed by Jackson (Reference Jackson1997), which is based on the convolution of the point quantities with some weighting function. However, we checked that both methods yield very similar quantities, provided that the width of the weighting function is conveniently chosen.

Due to the periodicity of the system in the $x$ and $z$ directions, the only relevant coordinate for the averaged quantities is $y$, and the variation with the former coordinates is averaged out. In addition, we are not interested in variation with $y$ at a scale significantly smaller than the radius $a_1$. As a consequence, the simulation volume is split into slices, $\Delta y=a_1$ in width, and the various quantities are averaged over each slice, which results in a sampled version of their averaged counterpart (figure 2).

Figure 2. The various quantities are averaged over each position bin, either over the whole slice volume, only over the particle volume or over the liquid volume.

The various quantities may be grouped into different classes according to the way they are computed. First, the suspension velocity is computed as

(4.1)\begin{equation} \boldsymbol{u}_{s}(y_n) = \dfrac{1}{L_{x}L_{z}\Delta y}\int_{y_n- \Delta y / 2}^{y_n + \Delta y /2}\iint \boldsymbol{u}(\boldsymbol{x})\, \mathrm{d}\kern0.07em x\, \mathrm{d}z\, \mathrm{d}\kern0.05em y, \end{equation}

where $\boldsymbol {u}(\boldsymbol {x})$ is the velocity field from the solver and the surface integral is computed over a plane at constant $y$. Then, some solid phase quantities result from averaging only over the particle volume. This is the case of the averaged solid volume fraction, which is computed from the particle positions and radii according to

(4.2)\begin{equation} \phi(y_n) = \dfrac{1}{L_{x}L_{z}\Delta y}\int_{y_n- \Delta y / 2}^{y_n + \Delta y /2}\iint\sum_{p} \chi_{p}(\boldsymbol{x})\, \mathrm{d}\kern0.05em x\, \mathrm{d}z\, \mathrm{d}\kern0.05em y, \end{equation}

where $\chi _p(\boldsymbol {x})$ is the indicator function of particle $p$, equal to 1 inside the particle and 0 outside. Similarly, the particle phase velocity is computed as

(4.3)\begin{equation} \phi(y_n) \boldsymbol{u}_{p}(y_n) = \dfrac{1}{L_{x}L_{z}\Delta y}\int_{y_n- \Delta y / 2}^{y_n + \Delta y /2}\iint \sum_{p} \boldsymbol{U}_{p} \chi_{p}(\boldsymbol{x})\, \mathrm{d}\kern0.07em x\, \mathrm{d}z\, \mathrm{d}\kern0.05em y. \end{equation}

We note here that each point inside a particle is assumed to move at the particle centre velocity $\boldsymbol {U}_p$, meaning that the particle angular velocity is not accounted for in the definition of $\boldsymbol {u}_{p}$. We also considered another definition for the particle velocity field, where $\boldsymbol {U}_p$ is replaced by the true velocity $\boldsymbol {U}_p + \boldsymbol{\varOmega} _p \times (\boldsymbol {x} - \boldsymbol {x}_p)$. It turned out that the computed quantities were nearly identical, except for a small discrepancy in the vicinity of the bounding walls, where particle layering occurs (Orsi et al. Reference Orsi, Lobry and Peters2023). Eventually, we kept the expression in (4.3) since its implementation is easier and its computation faster. It should be finally noted from the latter equation that $\phi (y_n) \boldsymbol {u}^{p}(y_n)$ is the solid volume flux density, averaged over a slice, consistent with the solid volume balance equation (3.1).

The fluid phase pressure is computed similarly, except that the volume integration is performed over the volume outside the particles,

(4.4)\begin{equation} [1-\phi(y_{n})]\langle P \rangle^{f}(y_{n}) = \dfrac{1}{L_{x} L_{z}\Delta y} \int_{y_{n}-\Delta y/2}^{y_{n}+\Delta y/2}\iint P(\boldsymbol{x}) [1-\chi_{p}(\boldsymbol{x})]\, \mathrm{d}\kern0.07em x\, \mathrm{d}z\, \mathrm{d}\kern0.05em y. \end{equation}

The general particle phase averaged quantity, $n\langle Y \rangle _p$, has the meaning of the volume density of a quantity $Y_p$ defined for each particle $(p)$, and $n$ is the number density of particles. For instance, $n \langle \boldsymbol {f} ^h \rangle _p$ is the hydrodynamic force density exerted on the particles and is computed from the hydrodynamic force on each particle $\boldsymbol {F}^h_p$. In principle, such a quantity should be computed as the volume average of superimposed Dirac-type functions $\sum _p Y_{p} \delta (\boldsymbol x - \boldsymbol x_p)$. To get a smoother averaged quantity, and since we are not interested in variations at a very small scale, the density $Y_{p} \delta (\boldsymbol x - \boldsymbol x_p)$ is smeared over the particle volume, so that

(4.5)\begin{equation} n \langle Y \rangle _p(y_{n}) = \dfrac{1}{L_{x}L_{z}\Delta y} \int_{y_{n}-\Delta y/2}^{y_{n}+\Delta y/2}\iint \sum_{p} \dfrac{Y_{p}}{\frac{4}{3}\pi a_{p}^{3}} \chi_{p}(\boldsymbol{x})\, \mathrm{d}\kern0.05em y\, \mathrm{d}\kern0.07em x\, \mathrm{d}z. \end{equation}

Equation (4.5) holds for the hydrodynamic and contact forces, and force dipole densities, $n \langle \boldsymbol {f} ^i \rangle _p$ and $n \langle \boldsymbol{\mathsf{D}}^i \rangle _p$ of interest in (3.5)–(3.10).

Finally, the contribution of the external pressure $P_0 = ({\mathrm {d}P_0}/{\mathrm {d}\kern0.07em x})x$ to the stress needs to be tackled separately, since it does not depend on the $y$ coordinate, but instead, its gradient is directed along the $x$ direction. The contribution of $P_0$ cannot be directly computed using the integrals above, since its gradient would be averaged out in the process – see for instance (4.4). However, it cannot be overlooked, since it drives the flow, and has to be accounted for in the stress balance. In more detail, the pressure is only involved in two terms in the expression of the fluid stress tensor, (3.10), namely the fluid pressure $(1 - \phi ) \langle P \rangle ^f \boldsymbol{\delta}$ and the hydrodynamic force dipole trace $n \langle s^{FDM} \rangle _p/3 \boldsymbol{\delta}$ (2.9). Inspection of (2.9), (4.4) and (4.5) shows that the sum of the two contributions involves the local average of $P_0$ over the suspension and therefore yields approximately $-P_0 \boldsymbol{\delta}$, which is quite intuitive. The latter term is thus simply added to the fluid phase stress, while the fluid pressure and the hydrodynamic force dipole are computed from the perturbation pressure $p$ only. It should also be noted that the hydrodynamic force exerted on the particles originating in the applied pressure $P_0$ is inherently accounted for in (2.5).

4.2. Material functions

The standard material functions in shear flow are locally measured at each position in the channel. First, the suspension stress is computed from (3.6), (3.10) and (3.11):

(4.6)\begin{align} \boldsymbol{\varSigma} = \boldsymbol{\varSigma}^ f+\boldsymbol{\varSigma}^ c &={-} \dfrac{\mathrm{d}P_0}{\mathrm{d}\kern0.07em x} x \boldsymbol{\delta} - (1- \phi)\langle p \rangle ^f \boldsymbol{\delta} \nonumber\\ &\quad + \eta \partial_y u^s_x (\boldsymbol{e}_x \otimes \boldsymbol{e}_y + \boldsymbol{e}_y \otimes \boldsymbol{e}_x) + n\langle \boldsymbol{D}^h \rangle _p + n\langle \boldsymbol{D}^c \rangle _p, \end{align}

where the last term is the contact contribution to the stress $\boldsymbol{\varSigma} ^c$, while the remaining terms compose the fluid (or hydrodynamic) stress. We note here that, when considering the tangential ($xy$) component of the total stress, it is only needed to take into account the symmetric part of the force dipole densities since the antisymmetric part of the hydrodynamic and contact dipoles cancel each other out exactly due to the lack of angular inertia, which is written in (2.4). In addition, we checked that the antisymmetric part of each locally averaged dipole density, either fluid or contact in nature, is very small compared with the corresponding symmetric part.

As shown in § 5.2, the shear rate may be computed indifferently from the particle or suspension velocity profiles,

(4.7)\begin{equation} \dot{ \gamma} = \vert \partial_y u^s_x \vert \approx \vert \partial_y u^p_x \vert. \end{equation}

In rate-independent suspensions in HSF, all stresses are proportional to each other and the shear rate $\dot \gamma$, so that the material functions are defined as the ratio of the relevant stresses, and only depend on the volume fraction. The reduced suspension viscosity and normal stress differences are defined as usual,

(4.8)$$\begin{gather} \eta_S = \dfrac{ \varSigma_{xy}}{\eta \partial _y u^s_x} = 1 + \dfrac{n \langle {\mathsf{S}}_{xy}^{h}\rangle _p + n \langle {\mathsf{S}}_{xy}^{c}\rangle _p }{ \eta \partial _y u^s_x}, \end{gather}$$

(4.9a,b)$$\begin{gather}\hat N_1 = \dfrac{\varSigma_{xx} - \varSigma_{yy}}{\vert \varSigma_{xy}\vert},\quad \hat N_2 = \dfrac{\varSigma_{yy} - \varSigma_{zz}}{\vert \varSigma_{xy}\vert}, \end{gather}$$

and the reduced contact tangential and normal stresses write as

(4.10a–d)\begin{equation} \hat \varSigma_{12}^c = \dfrac{\varSigma_{xy}^c}{\vert \varSigma_{xy}\vert},\quad \hat \varSigma_{11}^c = \dfrac{\varSigma_{xx}^c}{\vert \varSigma_{xy}\vert},\quad \hat \varSigma_{22}^c = \dfrac{\varSigma_{yy}^c}{\vert \varSigma_{xy}\vert},\quad \hat \varSigma_{33}^c = \dfrac{\varSigma_{zz}^c}{\vert \varSigma_{xy}\vert}, \end{equation}

where $\varSigma _{ij}^c =n \langle S_{ij}^{c}\rangle _p$.

An alternative description, namely the $\mu (J)$ rheology, is now standard (Boyer et al. Reference Boyer, Guazzelli and Pouliquen2011; Guazzelli & Pouliquen Reference Guazzelli and Pouliquen2018), where the material functions depend on the viscous number $J$:

(4.11a,b)\begin{equation} J = \dfrac{\eta \dot \gamma}{\varSigma^c_{yy}},\quad \phi = \phi(J),\quad \mu = \dfrac{\varSigma_{xy}}{\varSigma^c_ {yy}}=\dfrac{1}{\hat{\varSigma}_{22}^c}= \mu(J). \end{equation}

5.1. Transient particle migration

We focus first on the transient flow during which the volume fraction and velocity profiles develop. Figure 3 displays the time evolution of the volume fraction at various positions in the channel for the mean volume fraction $\bar {\phi }=0.4$. The particle migration towards the channel centre is evidenced. Steady-state is reached approximately at time $t \approx 1000$, corresponding to the mean accumulated strain $\bar \gamma \approx 134$, except at the centre of the channel, where long-time migration is slower and a steady state seems to be reached at $t\approx 2000$ ($\bar \gamma \approx 291)$. The same data are used to plot the evolution of the volume fraction profile over time in figure 4(a). Each profile results from an averaging over the time interval $\Delta t=100$, which is enough to get rid of part of the statistical fluctuations. This interval is extended up to $\Delta t=1000$ for the steady state profile. The resulting profiles are reasonably smooth, except in the vicinity of the bounding walls, where layering occurs as expected, even though the suspension is bidisperse (Howard et al. Reference Howard, Maxey and Gallier2022, Reference Howard, Dong, Patel, D'Elia, Maxey and Stinis2023). We note that layering is not completely apparent due to our loose sampling. However, we clearly evidenced it in a previous article (Orsi et al. Reference Orsi, Lobry and Peters2023), where it was shown to extend over around $3a_1$ at $\bar {\phi }=0.45$. The corresponding particle velocity profiles are displayed in figure 4(b). A parabolic profile typical of a uniform viscosity field classically turns into a more blunted profile due to the viscosity variations that reflect the enhanced volume fraction at the channel centre. A significant wall slip is also evidenced, which originates in the mentioned layering. Such a wall slip has often been observed in experiments (Jana, Kapoor & Acrivos Reference Jana, Kapoor and Acrivos1995; Sarabian et al. Reference Sarabian, Firouznia, Metzger and Hormozi2019) or in simulations of monodisperse (Yeo & Maxey Reference Yeo and Maxey2011; Gallier et al. Reference Gallier, Lemaire, Peters and Lobry2014b, Reference Gallier, Lemaire, Lobry and Peters2016) and bidisperse (Howard et al. Reference Howard, Maxey and Gallier2022; Orsi et al. Reference Orsi, Lobry and Peters2023; Peerbooms et al. Reference Peerbooms, Nadorp, van der Heijden and Breugem2024) suspensions. We checked that the particle and suspension instantaneous velocity profiles are very close to each other, with a somewhat larger discrepancy in the layering area, as also shown at steady flow in § 5.2.

Figure 3. Time evolution of the volume fraction at different $y$-positions in the channel for $\bar {\phi } = 0.4$. Averaging for steady flow starts from the vertical red line.

Figure 4. Volume fraction and velocity profiles during particle migration for $\bar {\phi } = 0.4$. The vertical dashed lines indicate the volume boundaries.

In addition to the variation of the overall particle volume fraction, the spatial distribution of the small and large particles is computed. Figure 5 displays the time evolution of the volume fraction of small and large particles. A running average has been performed with a kernel of width $\Delta T = 20$ to decrease the amplitude of the fast fluctuations. The small particle volume fraction approximately follows the evolution of the mean volume fraction, although at a slower rate. The variation of the large particle volume fraction is somewhat more complicated, with a faster variation at a small time, and long-time reorganization. On the whole, it seems that the latter is not complete at the end of the simulation.

Figure 5. Small ($\phi _1$) and large ($\phi _2$) particles volume fraction as a function of time. Same simulation and legend as in figure 3. Upper scale is typical accumulated strain (2.17). High-frequency fluctuations have been averaged out using a simple running average ($\Delta T = 20$) to enhance readability. Same positions in the channel as in figure 3.

To get a clearer understanding of the evolution of both populations, the transient volume fraction profiles of the small and large particles are displayed in figure 6, both at the beginning (figure 6a,b) and at the end (figure 6c,d) of the simulation. At the beginning, a faster migration of the large particles takes place towards the channel centre, compared with the small particles, resulting in a peak of the large particle volume fraction at the centre, together with dips at the base of the peak. At larger time, reorganization of the particles occurs, with a deepening of the dips in the large particle volume fractions, together with the building of secondary peaks at the same position in the small particle volume fraction. It should be noted that during the time interval $[900, 3000]$, the overall volume fraction profile only weakly evolves, mostly at the centre of the channel (figure 4).

Figure 6. Evolution with time of the small (a,c) and large (b,d) particle volume fraction profiles at the beginning (a,b) and at the end (c,d) of the simulations. $\bar {\phi }=0.4$. Vertical dashed lines are boundary positions.

This is in qualitative agreement with the simulations of bidisperse suspensions, $0.3$ in volume fraction, by Chun et al. (Reference Chun, Park, Jung and Won2019) who showed that particle segregation is much slower than the overall particle migration, except for a peak of the large particle volume fraction at the channel centre. We stress again that long-time particle reorganization seems to affect only weakly the total solid volume fraction and flow velocity profiles (figure 4). The time and accumulated strain for steady averaging are given in table 1 for each simulation.

5.2. Steady profiles

In the present section, we focus on the steady volume fraction and velocity profiles. Figure 7 displays the steady volume fraction profiles for all three values of the mean volume fraction $\bar {\phi }$, together with the specific volume fraction profile for the large and small particles. The jamming volume fraction $\phi _J = 0.592$ deduced from simulations of homogeneous bidisperse suspensions in simple shear flow (§ 4) is also displayed. In the less concentrated case $\bar {\phi }=0.3$, the volume fraction at the channel centre does not reach $\phi _J$, in agreement with several experimental (Lyon & Leal Reference Lyon and Leal1998; Rashedi et al. Reference Rashedi, Sarabian, Firouznia, Roberts, Ovarlez and Hormozi2020) and numerical (Yeo & Maxey Reference Yeo and Maxey2011; Chun et al. Reference Chun, Park, Jung and Won2019; Howard et al. Reference Howard, Maxey and Gallier2022) studies. This is, however, in contrast to SBM predictions, according to which the volume fraction $\phi _J$ is reached at the channel centre whatever the mean volume fraction $\bar {\phi }$. Turning to the distribution of the large and small particles ($\phi _2$ and $\phi _1$, respectively), both migrate towards the centre, with approximately uniform composition $\phi _{2}/\phi \approx 0.5$ except at the bounding walls where only small particles are found. An explanation for the latter behaviour has been proposed by Howard et al. (Reference Howard, Maxey and Gallier2022) based on interactions between large and small particles close to a wall.

Figure 7. (a–c) Total volume fraction ($\phi$), large ($\phi _2$) and small ($\phi _1$) particles volume fraction profiles at a steady flow. (d–f) Ratio of the large particles to overall solid volume fraction. Here (a,d) $\bar \phi = 0.3$; (b,e) $\bar \phi = 0.4$; (c,f) $\bar \phi = 0.5$.

The simulations at larger mean volume fraction $\bar {\phi }$ are somewhat different from the previous case. The overall volume fraction exceeds $\phi _J$ at the channel centre (the ‘plug’), over a width significantly larger than the particle size: at $\bar {\phi } = 0.5$, the volume fraction reaches $0.65$ and exceeds $\phi _J = 0.592$ over $10a_1$. Such a feature has been evidenced by Oh et al. (Reference Oh, Song, Garagash, Lecampion and Desroches2015) in their experiments on migration in non-Brownian suspensions of monodisperse particles in a pipe flow. They show that, for a large enough mean volume fraction, the volume fraction at the centreline of the tube levels off to a value close to the random close packing fraction $\phi _{rcp} \approx 0.64$. We note here that this value of the volume fraction is also known from simulation (Gallier et al. Reference Gallier, Lemaire, Peters and Lobry2014b; Mari et al. Reference Mari, Seto, Morris and Denn2014; Singh et al. Reference Singh, Mari, Denn and Morris2018; Chèvremont et al. Reference Chèvremont, Chareyre and Bodiguel2019) to be approximately the suspension jamming volume fraction when frictionless particles are considered. Regarding the local composition of the suspension, the large particle volume fraction presents a peak at the channel centre, as expected, but also two dips at a distance from the centre around $5a_1$. The small particles migrate towards the centre too, although with a peak at the positions of the previously mentioned dips. Leaving aside the layering area at the bounding wall, we note that the proportion of large particles $\phi _2/\phi$ varies over a quite wide range $[0.35, 0.74]$. Despite these variations, the overall volume fraction smoothly varies over the channel width.

The relative volume fraction profile at $\bar {\phi }=0.4$ is in qualitative agreement with the data computed by Howard et al. (Reference Howard, Maxey and Gallier2022) for frictionless particles using the FCM, with slightly different ratio $a_2/a_1$ though, which parameter they show to have a clear influence on segregation. We note that the dips in the large particle volume fraction profile are less pronounced in their data. To check the influence of friction, we also performed a single simulation for frictionless particles, and we compare in figure 8 the volume fraction profile with data by Howard (Reference Howard2018), for the same ratio $a_2/a_1 = 1.4$ and mean volume fraction $\bar {\phi }=0.4$

Figure 8. Nearly steady total volume fraction profile ($\phi$), large ($\phi _2$) and small ($\phi _1$) particles volume fraction profiles. Frictionless particles $\mu _s=0$. Symbols are short time average from FDM simulations. Lines are FCM simulations by Howard (Reference Howard2018). Vertical dashed lines are boundary positions.

We checked that the total volume fraction has approximately reached steady state at the observed total strain $\bar {\gamma } = 325$, while the small and large particle concentration is still somewhat evolving. The agreement is very good. In particular, the marked dips in the large volume fraction profile observed close to the centre plane for frictional particles are significantly weaker for frictionless particles. Slight discrepancies may be observed, though, mainly in the total volume fraction profile. This may be explained by the larger system size that Howard (Reference Howard2018) and Howard et al. (Reference Howard, Maxey and Gallier2022) consider, namely $L_y = 56 a_1$. We note that in both simulations, the volume fraction remains lower than the jamming volume fraction for frictionless particles, i.e. $\phi _J^{\mu _s=0}\approx 0.64\hbox{--} 0.65$ (Mari et al. Reference Mari, Seto, Morris and Denn2014; Gallier, Peters & Lobry Reference Gallier, Peters and Lobry2018), even at the channel centre. This latter feature highlights the influence of friction, which increases the contact contribution to the stress (Gallier et al. Reference Gallier, Lemaire, Peters and Lobry2014b; Peerbooms et al. Reference Peerbooms, Nadorp, van der Heijden and Breugem2024), thereby enhancing particle migration, as implied by (3.15).

Back to frictional particles, figure 9 displays the particle and suspension axial velocity profiles, together with their relative difference. As already mentioned, the velocity profiles are blunted compared with the parabolic profile observed in a homogeneous Newtonian liquid. This is all the more the case as a larger value of $\bar {\phi }$ is considered. It is also clear that the suspension and particle velocities are very close to each other, except in the layering area at the boundaries, where the difference grows while oscillating as the boundary is approached (the corresponding points are out of the scale in figure 9d–f).

Figure 9. (a–c) Particle ($u_x^p$) and suspension ($u_x^s$) velocity profiles at a steady flow. (d–f) Profile of the relative difference between suspension and particle velocity. Same $\bar {\phi }$ as in figure 7.

In particular, there is no significant difference between the shear rate as computed either from one or the other velocity profiles, except in the layering area. This is clearly illustrated in figure 10, where the shear rate profiles from the time-averaged velocity gradients over the steady flow duration, are displayed. The shear rate varies over two or three orders of magnitude, quite smoothly over the whole channel, except in the vicinity of the bounding walls, due to particle layering, and also at the channel centre, the latter only in the case of the largest mean volume fraction $\bar {\phi }=0.5$. It should be stressed that the reliability of the shear-rate measurement is of importance since it is involved in the local determination of the constitutive laws. In particular, the time-averaged velocity gradient is expected to be exactly zero at the very centre of the channel and to be very small in the vicinity of the centreline. The relevance of such a small shear rate will now be addressed.

Figure 10. Shear rate profiles from the particle and suspension velocity profiles. Same values of $\bar {\phi }$ as in figure 7.

To gain some qualitative insight into this matter, the standard deviation of the velocity gradient is computed and compared with the average velocity gradient. In more detail, the velocity gradient profile is computed at each time step during the steady flow according to (4.1), yielding the time-averaged velocity gradient $\overline {\partial _{y}u_x^{s}}$ and the standard deviation $\sigma (\partial _{y}u_x^{s})$ profiles. The ratio of the latter to the former is displayed in figure 11 for all three values of the mean volume fraction $\bar {\phi }$. In the region of the flow where this ratio is larger than typically unity, the fluctuations of the shear rate profile are sufficiently intense for the instantaneous shear rate to change sign frequently. In this case, the flow behaviour may differ significantly from a steady simple shear flow, yielding quite a different microstructure. While such large fluctuations are limited to a narrow region close to the centreline, typically $3-4a_1$ in width, at $\bar {\phi }=0.4$, they occur over a $10a_1$-wide region at $\bar {\phi }=0.5$. In both cases, the high fluctuations area also corresponds to the plug region where the volume fraction exceeds the jamming volume fraction as computed in a HSF, suggesting that the suspension is jammed. In the same line, the computed accumulated strain during steady state in that region is smaller than unity. At $\bar {\phi } = 0.3$, the relative standard deviation is always lower than unity, consistent with the largest observed volume fraction ($\approx$0.5) lower than $\phi _J$.

Figure 11. Relative standard deviation of the suspension velocity gradient: (a) $\bar {\phi }=0.3$; (b) $\bar {\phi }=0.4$; (c) $\bar {\phi }=0.5$.

It is worth noting that the velocity gradient fluctuations computed from the particle velocity field (not displayed) do not differ from the suspension velocity gradient fluctuations.

6.1. Momentum balance in steady flow

Before the local constitutive laws are addressed, the mechanical consistency of the numerical model will be checked. In the following section, the momentum balance for the suspension, i.e. (3.12) together with (4.6), is examined, in both the velocity and velocity gradient directions.

6.1.1. Tangential stress $\varSigma _{xy}$

The $x$-component of (3.12) is written as

(6.1)\begin{equation} {-}\dfrac{\mathrm{d}P_0}{\mathrm{d}\kern0.07em x} + \partial_y ( \eta \partial_y u_x^s +n\langle {\mathsf{S}}_{xy}^h \rangle _p + n\langle {\mathsf{S}}_{xy}^c \rangle _p )=0. \end{equation}

The first term is the imposed pressure gradient. The second and third terms constitute, respectively, the hydrodynamic contribution of the pure fluid and of the particles to the stress, while the last term is the contact force contribution. Equation (3.12) classically reflects that the tangential stress is driven by the pressure gradient. Its validity will be checked at steady flow in the following.

The total steady dimensionless shear stress profile is displayed in figure 12(a), together with the contact and hydrodynamic contributions, for the mean volume fraction $\bar {\phi } = 0.4$. The theoretical shear stress $({\mathrm {d}P_0}/{\mathrm {d}\kern0.07em x})y$ (i.e. $4y$ in dimensionless form) computed from the imposed pressure gradient using (6.1) is also displayed. Except in the layering area, a very good agreement is found between the total shear stress profile and its theoretical expression, including symmetry and cancelling at the channel centre, reflecting the symmetry of the averaged stress boundary conditions. The momentum balance equation may also be examined in its direct, (6.1), and not integrated form in figure 12(b), where the force densities (i.e. the stress gradients) have been computed using numerical differentiation. The uniform tangential stress gradient again reflects the good local behaviour of the FDM numerical model.

Figure 12. (a) Steady tangential stress profile, together with contributions detailed in (6.1). (b) Corresponding force densities. $\bar {\phi }= 0.4$.

The quantitative relation between material functions and volume fraction will be addressed in detail in § 6.2. However, figure 12 allows qualitative inspection of the relative importance of the different contributions. The ratio of the contact contribution to the total shear stress increases with the volume fraction, and the contact contribution exceeds the hydrodynamic contribution when the volume fraction exceeds approximately $0.42$, as expected from the usual constitutive laws in HSF (§ 6.2). We also note that the force density from the hydrodynamic shear stress exceeds the contact force density for $\vert y^* \vert \geq 0.11$, i.e. for $\phi \leq 0.5$.

Finally, the same curves as in figure 12 may be drawn for the other values of $\bar {\phi }$. They are displayed in Appendix A (figures 28 and 29), with the same features as described above. It should be noted that the same general behaviour is observed during the transient migration, whether instantaneous quantities are considered, or short-time averaging is performed, with more intense fluctuations, though.

6.1.2. Normal stress $\varSigma _{yy}$

The $y$-component of (3.12) is written as

(6.2)\begin{equation} \partial_{y}\varSigma_{yy} = \partial_{y}\varSigma_{yy}^{f} + \partial_{y}\varSigma_{yy}^{c} = 0. \end{equation}

The total $yy$ component of the total stress is expected to be uniform, meaning that the variation of the fluid stress should compensate for the variation of the contact contribution. This behaviour is observed in figure 13(a), except in the layering area close to the walls, both at the beginning of the simulation and during steady flow, when migration is complete. Uniformity of the total stress is especially striking when the solid volume fraction is homogeneous as the simulation starts, since both the fluid and solid contributions are strongly non-uniform, in agreement with the SBM according to which the particles are pushed towards the channel centre due to the divergence of the contact stress, (3.15).

Figure 13. Total stress profile (a) together with fluid (b) and contact (c) contributions, at the beginning of the simulation and during steady flow: $\bar {\phi } = 0.4$.

When particle migration is complete, both the fluid and contact contributions are uniform, in agreement with the SBM, except at the channel centre, where a weak stress gradient persists. The latter is not predicted by the SBM. However, the fluid and contact contributions still counterbalance each other.

Again, the corresponding curves for $\bar {\phi } = 0.3$ and $0.5$ are displayed in Appendix A.2 (figures 30 and 31). The total stress is uniform as well for both volume fraction values. However, in the case of the lower mean volume fraction, $\bar {\phi } = 0.3$, the contact and fluid stresses are slightly non-uniform after the migration is finished. This will be discussed in § 7.

6.2. Local constitutive laws

The goal of the present section is to determine to what extent the main constitutive laws evidenced in HSF simulations still locally hold in the pressure-driven flow. To keep the noise level low, all needed quantities are averaged over the steady flow duration. The main results displayed here will be discussed in § 6.3. Also displayed are the correlation laws for the normal stress differences and for the reduced contact stresses by Badia et al. (Reference Badia, D'Angelo, Peters and Lobry2022) which were fitted onto particle resolved simulations of bidisperse suspensions for the same friction coefficient $\mu _s=0.5$. We note that such correlation laws, with slightly different numerical coefficients, were recently found to satisfactorily describe simulated quantities for various values of the friction coefficients (Peerbooms et al. Reference Peerbooms, Nadorp, van der Heijden and Breugem2024). Before turning to the material functions, the measurement resolution is briefly discussed.

6.2.2. Shear viscosity

The steady relative viscosity is computed from the time-averaged shear stress and the time-averaged suspension velocity gradient $\eta _S = \overline {\varSigma _{xy}}/(\eta \overline {\partial _y u^s_x})$ at each position in the channel. It is displayed as a function of the solid volume fraction $\phi$ for $\bar {\phi }=0.4$ in figure 14(a), together with the corresponding data from HSF simulations (Orsi et al. Reference Orsi, Lobry and Peters2023). The agreement is very good in an intermediate region, not too close to the wall or the centreline. Discrepancies are evidenced in the vicinity of the bounding walls (low volume fraction), where layering occurs, and in a central area (high volume fraction), approximately $9a_1$ in width. The data for which reasonable agreement is empirically found are highlighted in figure 14(a). The corresponding region is denoted by ‘the intermediate area’ in the following. The same procedure is repeated for the other simulations. The local viscosity data over the whole channel are displayed for all values of $\bar {\phi }$ in Appendix B.1 (figure 32), and the data from the intermediate region are gathered in figure 14(b), showing again very good agreement over a wide volume fraction range $0.27 \leq \phi \leq 0.55$. Also displayed is the Maron–Pierce correlation law, which has been fitted onto the data from HSF simulations,

(6.3)\begin{equation} \eta_S = \dfrac{\alpha}{(1 - \phi / \phi_J )^2}, \end{equation}

yielding the jamming volume fraction $\phi _J = 0.592$.

Figure 14. Viscosity as a function of local volume fraction. Here (a) $\bar {\phi }=0.4$. (b) All simulations. Here ($\bullet$, blue) $\bar {\phi }=0.3$, ($\blacklozenge$, orange) $\bar {\phi }=0.4$, ($\blacktriangledown$, green) $\bar {\phi }=0.5$. Black symbols are HSF simulations; plain line is (6.3) with $\phi _J = 0.592$ and $\alpha = 0.8$.

We note (Appendix B.1, figure 32) that in the central region, for all simulations, the computed reduced viscosity is below its counterpart in HSF. In addition, no divergence is observed at the channel centre, which is allowed by the simultaneous vanishing of the shear stress and the shear rate. Finally, $\eta _S$ seems to be defined in the ‘jammed’ phase where $\phi >\phi _J$. However, as mentioned previously, the computed shear rate suffers from very high fluctuations in that region, making the reliability of its averaged value questionable.

For the sake of completeness, the various relative contributions to the viscosity, or equivalently the relative contributions to the shear stress, are displayed in figure 15 for the intermediate region. The relative hydrodynamic contribution to the viscosity is split into the pure fluid contribution equal to unity, the particle solver and subgrid contributions, denoted, respectively, by $\eta _{FDM}$ and $\eta _{SG}$:

(6.4)\begin{equation} \eta_{\boldsymbol{\mathsf{S}}} = 1 + \dfrac{n \langle {\mathsf{S}}_{xy}^{FDM}\rangle _p + n \langle S_{xy}^{SG}\rangle _p + n \langle S_{xy}^{c}\rangle _p } { \eta\ \langle \partial _y u^s_x \rangle} = 1+\eta_{FDM}+\eta_{SG}+\eta_{c}. \end{equation}

The data are in excellent agreement with those from HSF simulations. The relative contribution of the contact forces to the viscosity, i.e. $\hat \varSigma _{12}^c$ (4.10), is displayed for the three values of $\bar {\phi }$ in figure 35 (Appendix B.3). It turns out that in the central area, where the usual expression for the viscosity does not hold, the relative contact contribution is smaller than expected. This feature, together with the low value of the computed relative viscosity in the same region, suggests that the force microstructure is significantly different compared with what prevails in HSF. It should be noted that this behaviour is not limited to the region where $\phi > \phi _J$.

Figure 15. Relative contributions to the shear viscosity in the intermediate region: (a) contact and hydrodynamic contributions; (b) FDM contribution; (c) subgrid contribution; (d) pure fluid contribution. Same colour code as in figure 14(b).

6.2.3. Normal stress differences

The reduced first ($\hat N_1$) and second ($\hat N_2$) normal stress differences in the intermediate region are displayed in figure 16. The data from the pressure-driven flow simulations are in reasonable agreement with the HSF data, although over a narrower position range compared with the shear viscosity. We note that, for $\bar {\phi } = 0.3$ and $0.4$ the discrepancy starts at the same position as the volume fraction gradient significantly increases, which corresponds, for $\bar {\phi } = 0.4$, to the position where the small particle volume fraction exceeds its counterpart for the large particles. However, a significant discrepancy is observed over the whole channel at $\phi =0.5$ for the first normal stress difference. Regarding the central region (Appendix B.2, figure 33), we also note that the values of the first normal stress difference largely differ from the HSF correlation, including its sign, which again suggests a different microstructure in this area.

Figure 16. Reduced normal stress differences in the intermediate region. Same colour code as in figure 14(b). Plain line is correlation law by Badia et al. (Reference Badia, D'Angelo, Peters and Lobry2022).

6.2.4. Contact normal stresses

The present section is devoted to the constitutive laws that define contact normal stresses. Since the latter, especially $\hat \varSigma _{22}^c$ in the present case, are involved in the particle transverse velocity, (3.15), and consequently in the solid volume transport equation, (3.18), it is worth checking that the laws deduced from HSF simulations still hold in the present pressure-driven flow. All three contact stresses are displayed in figure 17 for all simulations, in the volume fraction range that corresponds to the intermediate region of each run. Again, a good agreement with HSF simulations is observed, although the reduced normal stresses as computed here are slightly larger than the HSF data.

Figure 17. Reduced contact normal stresses in the intermediate region. Same colour code as in figure 14(b). Plain line is correlation law by Badia et al. (Reference Badia, D'Angelo, Peters and Lobry2022).

In the central region, the HSF constitutive laws for the contact normal stresses hold better compared with the shear viscosity (Appendix B.4, figures 36–38). Even though at the centre of the channel, where the tangential stress vanishes, the reduced contact stresses diverge, consistent with the non-vanishing contact normal stresses at this position, it is, however, striking that at $\bar {\phi } = 0.3$, the discrepancy is only noticeable for $\vert y \vert \leq a_1$ or only for $\phi \geq \phi _J$ in the two most concentrated cases.

6.2.5. The $\mu (J)$ rheology

Finally, the $\mu$-rheology (Boyer et al. Reference Boyer, Guazzelli and Pouliquen2011) is computed. The macroscopic friction coefficient $\mu = \varSigma _{xy}/\varSigma _{yy}^c$ and volume fraction are displayed as a function of the viscous number $J=\eta \dot {\gamma }/\varSigma _{yy}^c$ in figure 18(a,b), respectively, together with the HSF results. A good agreement is found in the intermediate region for both quantities. Several correlation laws are available in the literature. In the present paper, we compare our simulation data with the experimental laws by Boyer et al. (Reference Boyer, Guazzelli and Pouliquen2011) and to the correlation laws fitted on our own earlier simulations for a microscopic friction coefficient $\mu _s=0.5$ (Badia et al. Reference Badia, D'Angelo, Peters and Lobry2022). We note that recent numerical studies extensively investigated the influence of the friction coefficient $\mu _s$ on these correlation laws (Chèvremont et al. Reference Chèvremont, Chareyre and Bodiguel2019; Peerbooms et al. Reference Peerbooms, Nadorp, van der Heijden and Breugem2024). The displayed correlation laws deserve a few words. The correlation laws by Badia et al. (Reference Badia, D'Angelo, Peters and Lobry2022) have been fitted onto data concerning moderately concentrated to concentrated suspensions ($\phi \le 0.56$) at $\mu _s=0.5$, which explains why they are close to the HSF data, in particular the macroscopic friction coefficient $\mu$, for the whole range of viscous number $J\in [2.10^{-3}, 2]$. However, the limit of $\mu$ for vanishing $J$ may not be considered, since it was not confirmed by simulation data. On the other hand, the experimental correlation laws by Boyer et al. (Reference Boyer, Guazzelli and Pouliquen2011) were tuned for small values of $J$ ($\lesssim 0.2$) and thus may not perfectly fit the high $J$ range (d'Ambrosio, Blanc & Lemaire Reference d'Ambrosio, Blanc and Lemaire2023).

Figure 18. Macroscopic friction coefficient $\mu$ (a) and volume fraction (b) as a function of viscous number $J$ in the intermediate region. Same colour code as in figure 14(b). Plain line is correlation by Badia et al. (Reference Badia, D'Angelo, Peters and Lobry2022), dashed line is correlation by Boyer et al. (Reference Boyer, Guazzelli and Pouliquen2011)

In the central region, the macroscopic friction coefficient $\mu$ is lower than its counterpart in HSF, as displayed in Appendix B.5 (figures 39 and 40). This may be qualitatively checked in figures 12 and 13, where the shear stress $\varSigma _{xy}$ and the shear rate $\dot \gamma$ vanish at the channel centre, while the contact stress $\varSigma _{yy}^c$ does not, meaning that $\mu$ goes to zero at vanishing $J$, while it is expected to level off to a value close to $0.3$ in HSF (Boyer et al. Reference Boyer, Guazzelli and Pouliquen2011; Gallier et al. Reference Gallier, Lemaire, Peters and Lobry2014b). In the same region, the local volume fraction is larger than predicted by the HSF law.

This is to be compared with the numerical study of highly concentrated frictionless non-Brownian suspensions by Bhowmik & Ness (Reference Bhowmik and Ness2024), who use Kolmogorov flow to induce non-homogeneous shear rate and particle migration. They also observe discrepancies between HSF constitutive laws at positions where the shear rate vanishes or is small. More specifically, the macroscopic friction coefficient $\mu$ decreases towards zero as the viscous number decreases, as in the present work.

7.1. Transverse velocity

The transverse particle velocity relative to the suspension, which is related to the particle volume flux density (4.3), is displayed in figure 19 at various times during the transient flow. It is directed towards the channel centre, with a higher intensity at the start of the flow, and as expected vanishes as the steady state is reached. Its magnitude is very low compared with the amplitude of the main flow velocity, at most $0.1\,\%$. We checked that $u_y^s$, which is computed from the solver velocity, is nearly zero, e.g. at $\bar {\phi }=0.4$ around $2 \times 10^{-6}$, in agreement with the incompressibility of the flow that is stated in (3.4). We checked that the fluctuations of the transverse particle velocity averaged over the whole steady flow duration are significantly larger than this residual transverse suspension velocity.

Figure 19. Transverse particle relative velocity, $\bar {\phi } =0.4$. Averaging time is $\Delta t = 200$.

7.2. Various force densities exerted on the particle phase

As recalled by Lhuillier (Reference Lhuillier2009), the momentum balance equation for the particle phase, (3.5), is equivalent to (3.8). The latter does not need checking, since it directly originates from the single particle momentum balance equation (2.4), which was double-checked during the validation of the simulation code. However, other relations, such as the relation between the contact force density and the divergence of the contact stress (3.7), deserve a closer look. In the following, the various force densities at play are displayed, and their relations to the relevant stress gradients are evidenced.

7.2.1. Contact force density

According to the SBM, particle migration is driven by the contact force density $n \langle f_y^c \rangle _p$, which is counterbalanced by the hydrodynamic drag force, yielding the particle volume flux. The former should also be equal to the divergence of the contact stress, $\partial _y \varSigma _{yy}^c$ (3.7), according to a theoretical analysis for pairwise interaction forces between particles (Irving & Kirkwood Reference Irving and Kirkwood1950; Jackson Reference Jackson1997). Equation (3.7) is expected to hold provided that the length scale defining the variation of averaged quantities is sufficiently large compared with the particle radius. It will be now examined. The contact force density along the $y$ direction is displayed in figure 20 for all simulations at the beginning of the flow and during steady flow after the migration is complete, together with the divergence of the contact stress along the same direction. The latter is computed from the contact stress by numerical differentiation (central difference). A very good agreement is observed, except in the layering region though, which discrepancy may be easily understood, given that all averaged quantities strongly vary at the particle scale in that region (Howard et al. Reference Howard, Maxey and Gallier2022; Orsi et al. Reference Orsi, Lobry and Peters2023). We stress that the force density and contact stress are independently computed, the former from the total contact force exerted on each particle, and the latter from the total contact stresslet. In addition, only two instants are considered here for the sake of figure clarity, but we checked that (3.7) holds at any time.

Figure 20. Transverse contact force density and its counterpart from the contact stress: solid symbol, $t\in [0,\, 200]$; open symbol, steady flow (averaging time in table 1); (a) $\bar {\phi }=0.3$; (b) $\bar {\phi }=0.4$; (c) $\bar {\phi }=0.5$.

As expected, the contact force density pushes the particles towards the channel centre at the beginning of the simulation, and approximately vanishes at the end of transient migration, at least in the intermediate region. In the central region, two force peaks remain after migration is complete at $\bar {\phi }=0.3$ and $\bar {\phi }=0.4$, which correspond to the dip observed in the contact stress profile, as already mentioned in § 6.1.2. At $\bar {\phi }=0.3$, a small force density also remains in the intermediate region after migration completion, which does not induce any net particle volume flux. This will be discussed in the following.

7.2.2. Hydrodynamic force density

The hydrodynamic force density $n \langle f_y^h \rangle _p$ is exactly the opposite of the contact force density, due to the absence of inertia in the motion of each particle, (2.4). In addition, (3.9), which connects the hydrodynamic force density exerted on the particles to the divergence of the fluid stress in the $y$ direction, may be understood as the consequence of (3.7), (3.8), (3.11) and (3.12), which have already been shown to apply. However, inspection of the various contributions of the normal stress $\varSigma _{yy}^f$ is instructive. Equations (2.10), (2.9) and (4.6) yield

(7.2)\begin{align} \varSigma_{yy}^f &={-}(1-\phi)\langle p \rangle ^f + n\left\langle S_{yy}^{FDM} + \dfrac{s^{FDM}}{3}\right\rangle _p + n\left\langle S_{yy}^{SG} + \dfrac{s^{SG}}{3}\right\rangle _p \nonumber\\ &={-}(1-\phi)\langle p \rangle ^f +\varSigma_{yy}^{FDM} +\varSigma_{yy}^{SG}, \end{align}

where $\varSigma _{yy}^{FDM}$ is computed from the force density $\boldsymbol{\lambda}$ and pressure perturbation $p$ inside the particles (2.9). The external pressure $P_0$ has been omitted, since it does not depend on $y$. Here $\varSigma _{yy}^{FDM}$ is displayed in figure 21 together with the contributions appearing in (7.2) for $\bar {\phi }=0.4$.

Figure 21. Contributions to the fluid normal stress $\varSigma _{{yy}}^f$: (a) time-averaged stress $t \in [0,\, 200]$; (b) steady profiles $t\in [2000,\, 3000]$; $\bar {\phi }=0.4$.

As already mentioned, at the beginning of the flow, a strong fluid stress gradient is observed, which reflects the hydrodynamic hindering force applied to the particles against their motion. Most of the stress gradient originates in the fluid solver contribution $-(1-\phi )\langle p \rangle ^f + \varSigma _{yy}^{FDM}$, meaning that the subgrid corrections have very weak effects on this stress component, unlike the tangential stress (figure 15). It should also be noted that the fluid pressure contribution $-(1-\phi )\langle p \rangle ^f$ results in approximately half the total stress gradient, which highlights the important effect of such a pore pressure, and more generally of the fluid flow, in the migration dynamics (see Guazzelli & Pouliquen (Reference Guazzelli and Pouliquen2018) and references therein). After the migration is complete, the remaining total stress gradient is very weak, and the intensities of the various contributions are similar, except at the channel centre, where the residual stress gradient is the strongest, and mostly originates in the solver contribution.

The same trends may be observed for the other values of the mean volume fraction $\bar {\phi }$ (Appendix C, figures 41 and 42), with slight variations. At $\bar {\phi }=0.3$, the subgrid contributions are relatively higher, while the opposite occurs at $\bar {\phi }=0.5$. We also note that in the latter case, the fluid stress is uniform in the central region as well, consistent with the same behaviour exhibited by the contact stress.

7.3. Relation between transverse particle flux and hydrodynamic force density

The main hypothesis of the SBM is the relation between the hydrodynamic force density applied to the particle phase and the transverse velocity of the particles relative to the suspension. In more detail, as already mentioned, it is usually assumed that the hydrodynamic force density applied on the particle phase is a drag force, meaning that the transverse relative velocity $u_y^p - u_y^s$ is proportional to the hydrodynamic force density (7.1). The usually assumed lack of a non-drag contribution in the hydrodynamic force deserves to be checked, which we intend to do now. For the sake of clarity, the transient migration and the steady flow will be separately examined in the following sections.

7.3.1. Transient flow

The transverse relative velocity $u_y^p - u_y^s$ is displayed in figure 22 as a function of the modified hydrodynamic force $2 \bar {a}^2/9 \eta \times f(\phi )/\phi \times n\langle f_y^h \rangle _p$. All data, e.g. $u_y^p$, $u_y^s$, $\phi$, etc., result from time-averaging at a particular $y$-position in the channel over the time $\Delta t_{{averaging}} = 200$. In addition, only the data from the intermediate region are displayed. Each colour corresponds to an instant during the transient flow. The mean radius has been arbitrarily chosen as $\bar {a} = (a_1 +a_2 )/2 = 1.2a_1$, and the exponent in the expression of the hindered settling function, (3.14), has been chosen as $n = 4.8$ since this particular value seems to allow best overall agreement for the whole set of simulations.

Figure 22. Transverse relative particle velocity as a function of hydrodynamic hindering force during transient flow: $\bar {\phi }=0.3$ (a); $\bar {\phi }=0.4$ (b); $\bar {\phi }=0.5$ (c). See symbol colours online. Data from the intermediate region only. Plain line is theoretical relation (3.13).

For all simulations, the transverse velocity obeys (3.13) with reasonable accuracy, especially for $\bar {\phi }=0.4$ and $\bar {\phi }=0.5$. In the case $\bar {\phi }=0.3$, and to a lesser extent $\bar {\phi }=0.4$, a significant discrepancy is observed. It turns out that, at large times, the velocity goes faster towards zero than predicted from the driving force using (3.13). In other words, at the end of the transient flow, a force remains with no particle velocity, equivalent to a small non-drag force. This will be confirmed in the next section, which is dedicated to the steady flow.

Since some stress constitutive laws were shown not to hold in an extended region at the centre of the channel, this particular area deserves specific attention. Figure 23 displays the same quantities as in figure 22, computed from the data in the central region of the channel. Both the velocity and the modified driving force are lower than in the intermediate region, and the data suffer from more noise. However, the same general behaviour as in the intermediate region is observed, which deserves several comments. First, at $\bar {\phi }=0.5$, the drag law is valid during the whole transient flow, while the volume fraction exceeds $\phi _J$ from $t = 600$. This suggests that the relation between flux and force is very robust, and hardly depends on the actual suspension microstructure, which is not too surprising, given that Richardson and Zaki correlation primarily applies in particle sedimentation, and is now shown to be valid in shear flow as well. This may be understood given that no long-range order is expected in a sheared suspension, but only preferential orientation of contacting pairs (Blanc et al. Reference Blanc, Lemaire, Meunier and Peters2013), at least for volume fractions not too close to the jamming point. As a consequence, the permeability of such a sheared suspension is expected to be close to that of a random suspension. Second, it should be noted that, at a large time, the force peaks at $\bar {\phi }=0.3$ and $\bar {\phi }=0.4$, which constitute most of the central region, contributing to the non-drag force. This will be confirmed in the next section.

Figure 23. Transverse relative particle velocity as a function of hydrodynamic hindering force during transient flow: (a) $\bar {\phi }=0.3$; (b) $\bar {\phi }=0.4$; (c) $\bar {\phi }=0.5$. See symbol colours online. Data from the central region only. Plain line is theoretical relation (3.13).

7.3.2. Steady flow

It has been noted in § 6.1.2 that at steady flow the contact and fluid stresses were both uniform across the channel at $\bar {\phi } = 0.5$, but a gradient of these quantities was observed in the central region in the other cases, and even in the intermediate region at $\bar {\phi } = 0.3$. In other words, a finite hydrodynamic force density may occur at some position in the channel while no particle transverse velocity is observed, at odds with the prediction of the standard SBM.

This is well illustrated in figure 24, where the same quantities as in the previous section are displayed for $\bar {\phi } = 0.3$ and $\bar {\phi } = 0.4$. The data at each position in the central and intermediate regions are averaged over the time interval at the end of the simulation, during which the variations of the volume fraction and velocity profiles seem to be insignificant. The specific time interval is quoted in table 1 for each run. The error bars are computed from the standard deviation of each time-averaged quantity, e.g. the particle velocity $u^p_y(y)$ at the $y$-position in the channel, which is estimated using standard methods (Allen & Tildesley Reference Allen and Tildesley2017; Orsi et al. Reference Orsi, Lobry and Peters2023). In the case of the lowest mean volume fraction $\bar {\phi }=0.3$, a small but significant force density persists, while the relative particle velocity vanishes. The data from positions in the central region are highlighted in figure 24, and contribute to this force density. Such a residual force density may be interpreted as a non-drag contribution to the hydrodynamic force. However, it should be noted that in our case, this contribution is quite small compared with the overall force intensity at the beginning of the simulation. At $\bar {\phi } = 0.4$, the same trend may be observed, less intense and slightly shadowed by noise. Regarding the velocity that does not appear to vanish exactly in figure 24, it should be kept in mind that the spurious suspension $y$-velocity from the solver (§ 7.1) is of the order $2\times 10^{-6}$, comparable in magnitude to the velocity displayed in figure 19, and smaller than its standard deviation. Thus, the small remaining velocity seems irrelevant.

Figure 24. Transverse relative particle velocity as a function of hydrodynamic hindering force during steady flow. Data from the intermediate and central regions. The latter is highlighted ($\square$). Plain line is theoretical relation from (3.13).

Nott et al. (Reference Nott, Guazzelli and Pouliquen2011) argued that the non-drag contribution to the hydrodynamic force could originate in pairwise interaction forces, which would be written as the divergence of stress, the latter built from the particle pair interaction forces in the same way as the contact stress is built from the contact forces ((2.12), (2.13) and (3.6)). The latter relation holds in the case of the subgrid forces, since, from the particular expression of the subgrid resistance matrix (Orsi et al. Reference Orsi, Lobry and Peters2023), it is possible to show that the subgrid stresslet applied on a particle pair is simply the first moment of the subgrid forces on this particle pair, so that

(7.3)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\varSigma} ^{SG} = n \langle\,\boldsymbol{f}^{SG} \rangle _p, \end{equation}

which we have checked in the same way as (3.7) (§ 7.2.1). It is therefore tempting to identify the non-drag force to this pairwise subgrid force, meaning that the drag force density would be the remaining force from the solver $n \langle \boldsymbol {f}^h \rangle _p^{drag} \approx n \langle \boldsymbol {f}^{FDM} \rangle _p$ while the ‘driving’ force would originate in the short scale pairwise forces $n \langle \boldsymbol {f}^{c} + \boldsymbol {f}^{SG} \rangle _p$. This identification cannot have more than a qualitative scope, since the separation between subgrid and solver forces (and stresslets) depends on the computational mesh size. To test this conjecture, the relation between the steady velocity and the steady solver force density is displayed in figure 25. At $\bar {\phi } = 0.3$ most of the points are gathered significantly closer to the frame origin than in figure 24, meaning that the solver force density is approximately zero as expected, which seems to confirm the identification of the pairwise subgrid lubrication forces as the source of a non-drag hydrodynamic force. The data from the central region lie further from the theoretical curve, showing that the proposed explanation for the non-drag force seems to fail in this area. At $\bar {\phi } = 0.4$, the same argument approximately holds, but the difference between solver and total hydrodynamic forces is less important, due to the weaker relative contribution of the subgrid forces to the normal stress for the considered volume fraction range, as the comparison between figures 21 and 41 (Appendix C) shows. We also note that at $\bar {\phi }=0.5$, where the probed volume fraction range is higher, the subgrid forces are relatively even weaker (Appendix C, figure 42), meaning that the non-drag force is approximately zero. Surprisingly, this is also true in the central region, which may be related to the high volume fraction $\phi \geq \phi _J$.

Figure 25. Transverse relative particle velocity as a function of solver hydrodynamic hindering force during steady flow. Data from the intermediate and central regions. The latter is highlighted ($\square$). Plain line is theoretical relation from (3.13).

Finally, as mentioned in § 5.1, while the overall migration is complete during what we call the steady flow, segregation is not, and particle reorganization is still occurring. However, we do not present extensive curves for the segregation dynamics since the data suffers a larger amount of noise, possibly due to the smaller number of particles that constitute each population.

7.4. The SBM migration predictions

For the sake of completeness, the predictions of the SBM are now shortly compared with the simulation data. The governing equations are recalled in Appendix D. They are based on (3.16)–(3.18), without any non-drag force. In addition, to prevent the volume fraction from exceeding the jamming volume fraction, the particle flux is forced to vanish as $\phi _J$ is reached (Snook et al. Reference Snook, Butler and Guazzelli2016; Badia et al. Reference Badia, D'Angelo, Peters and Lobry2022). This procedure allows carrying out the computations, but obviously cannot account for the large volume fraction observed at the channel centre in the simulations. A specific modelling would be required for this purpose, which we keep for a future study.

The predictions of the SBM are compared with the simulation data in figures 26 and 27 for the mean volume fraction $\bar {\phi } = 0.4$ and their counterparts for $\bar {\phi } = 0.3$ and $\bar {\phi } = 0.5$ are displayed in Appendix D.2. The profiles from the SBM have been obtained at a high spatial sampling rate (i.e. $\Delta y =0.2 a_1$) and then resampled down to the simulation profiles rate (i.e. $\Delta y =a_1$) in the same spirit as the averaging procedure in (4.2). This results in the spreading of the central peak, whose initial growth rate is decreased with respect to the raw profiles. Thanks to this procedure, no peak narrower than a particle radius is observed, in agreement with physical intuition.

Figure 26. Evolution of the (a,b) volume fraction and (c,d) suspension velocity profiles as a function of the mean accumulated strain $\bar {\gamma }$; $\bar {\phi } = 0.4$. The profiles are averaged over a strain interval equal to half the strain sampling period: (a,c) FDM data; (b,d) SBM predictions.

Figure 27. Evolution of (a) volume fraction and (b) suspension velocity as a function of mean accumulated strain at various positions in the channel; $\bar {\phi } = 0.4$. A running average has been applied to the data to get rid of the high-frequency fluctuations. Plain line is FDM data and the dashed line is SBM predictions.

A fair agreement is found between the simulation data and the SBM predictions, despite significant fluctuations in the former, due to the finite size of the simulation volume. This reflects the good agreement of the simulation data and the main constitutive laws involved in the SBM, at least in the intermediate region. A point-to-point comparison is not easy due to the quite significant volume fraction dispersion at $\bar {\gamma }=0$ that is sometimes observed, reflecting the mentioned fluctuations. For instance, in the case $\bar \phi = 0.4$ (figure 27), $\phi (\bar {\gamma }=0)$ extends over the range [0.39, 0.45]. In addition, in the layering area at the boundaries, the volume fraction and the velocity profiles do not match to the SBM prediction, since the main stress constitutive laws are not fulfilled, and the simulated coarse-grained volume fraction is smaller than predicted, as usually observed. In the central region, again, the volume fraction from the simulations either exceeds the jamming volume fraction ($\bar {\phi }=0.4$ and $0.5$), or does not reach it ($\bar {\phi }=0.3$), which has implications on the volume fraction – and velocity – profiles over the whole channel, making the detailed comparison difficult. On the other hand, the SBM predicts a very fine peak at the channel centre that develops very quickly. Such a peak, when it is narrower than the particle size, may not have a clear signification. The resampling procedure that was adopted aims at smoothing out this feature, but does not manage to do so in the most concentrated case (Appendix D.2, figure 45a,b), where the volume fraction at the channel centre as predicted by the SBM evolves faster than measured in the simulation.

Despite these issues, the overall migration dynamics are satisfactorily predicted by the SBM in the intermediate region between the layering area and the channel centre. However, we stress that such comparison yields far less precise information than the local computation of constitutive laws presented before, since the volume fraction and the velocity at a given position are affected by their counterpart at each position in the channel and for all earlier times.