In the literature, two different frameworks exist for describing the rheology of solid/liquid suspensions: (1) the “viscous” framework in terms of the relative suspension viscosity, $ \eta r$, as a function of the reduced solid volume fraction, $\varphi / \varphi m$, with $ \varphi m$ the maximum flowable packing fraction, and (2) the “frictional” framework in terms of a macroscopic friction coefficient, $\mu $, as a function of the viscous number, $ I v$, defined as the ratio of the viscous shear to the wall-normal particle stress. Our goal is to compare the two different frameworks, focusing on the effect of friction between particles. We have conducted a particle-resolved direct numerical simulation study of a dense non-Brownian suspension of neutrally buoyant spheres in slow plane Couette flow. We varied the bulk solid volume fraction from $ \varphi b=0.1$ to 0.6 and considered three different Coulomb friction coefficients: $ \mu c=0$, 0.2, and 0.39. We find that $ \eta r$ scales well with $\varphi / \varphi m$, with $ \varphi m$ obtained from fitting the Maron–Pierce correlation. We also find that $\mu $ scales well with $ I v$. Furthermore, we find a monotonic relation between $\varphi / \varphi m$ and $ I v$, which depends only weakly on $ \mu c$. Since $ \eta r=\mu / I v$, we thus find that the two frameworks are largely equivalent and that both account implicitly for Coulomb friction. However, we find that the normal particle stress differences, $ N 1$ and $ N 2$, when normalized with the total shear stress and plotted against either $\varphi / \varphi m$ or $ I v$, remain explicitly dependent on $ \mu c$ in a manner that is not yet fully understood.

## I. INTRODUCTION

In many industrial processes, solid–liquid suspensions play an important role, such as in food processing, 3D printing of solid rocket fuels, and solid–liquid fluidized-bed reactors. Next to that, suspensions are commonly encountered in daily life as, for example, in cosmetic products, paints, and fresh concrete. To be able to produce and transport these suspensions in an efficient manner, a thorough understanding of the suspension rheology is required, i.e., the behavior of the suspension under imposed shear.

For spherical particles, which we consider in this work, the relative viscosity has been extensively studied both experimentally [17–20] and numerically [21–27] over the years. This has led to a large amount of data for various kinds of particles with different size, roughness, and friction coefficients. The literature has shown that the data collapse fairly well to a master curve when the volume fraction is normalized by the maximum flowable packing fraction $ \varphi m$, which is dependent on the properties of the particles. The rheology described as a function of the normalized volume fraction $\varphi / \varphi m$ will be referred to as the *viscous rheology* in the remainder of this work.

In experiments, it is difficult to independently vary the solid friction of the particles; therefore, numerical simulations provide useful insight into the effect of solid friction on the relative viscosity. Sierou and Brady [21] used accelerated Stokesian dynamics to investigate the flow of hard spheres in shearing flow, reaching volume fractions of up to 60 vol. %. They concluded that the microstructure and rheology are extremely sensitive to particle contacts. Additionally, they observed ordering at volume fractions over $\varphi =0.5$, which affected the relative viscosity. The layering of particles decreased the relative viscosity at these volume fractions with respect to a disordered suspension. Yeo and Maxey [22] performed fully three-dimensional direct numerical simulations (DNSs), employing a force-coupling method to study the rheology of suspensions without solid friction. For the core region of the planar Couette geometry, the relative viscosity showed good agreement with the Eilers correlation [14]; however, for the wall region, layering of particles caused disparate results. Gallier *et al.* [23] studied the effect of the solid friction coefficient on the rheology using DNS, employing a fictitious domain method, finding that the viscosity increases with increasing Coulomb friction coefficient. By separating the hydrodynamic and contact contributions, they showed that the strong increase in relative viscosity near jamming is caused by the contact between particles. More recently, Chèvremont *et al.* [27] used discrete element method (DEM) simulations to study the effect of friction by solving for the particle contact and lubrication forces based on the assumption that the long-range hydrodynamic interactions do not play a role at high volume fractions [25]. They found a similar dependence of the rheology on the friction coefficient and good agreement with previous simulations. Additionally, the experimental data of Boyer *et al.* [19] showed good agreement with their simulations when a solid friction coefficient of $ \mu c\u22480.5$ was used.

The relative viscosity is not sufficient to describe the complete rheology of suspensions and is, therefore, not sufficient for constitutive modeling. The microstructure induces anisotropy in the particle stresses under shear, which leads to non-Newtonian behavior [4]. Therefore, the normal stress differences are important for the description of the suspension rheology, which are defined as $ N 1= \Sigma 11\u2212 \Sigma 22$ and $ N 2= \Sigma 22\u2212 \Sigma 33$ for the first and the second normal stress difference, respectively. Here, the subscripts 1, 2, and 3, refer to the direction of the mean flow, velocity gradient, and vorticity, respectively. The normal stress differences are difficult to measure, explaining the sparse experimental literature data and the large discrepancies among different datasets. Experiments [17,18,28,29] and numerical simulations [21–23,26,30] are largely in agreement for the second normal stress difference $ N 2$, which is negative and increases in magnitude with increasing volume fraction. There is, however, less consensus on the first normal stress difference $ N 1$. Large variations between datasets and experimental uncertainty give rise to inconclusive results on the sign of $ N 1$, but for low solid fractions, the magnitude is much smaller than that of $ N 2$. Some results show a strong increase in magnitude at higher volume fractions, either to positive or negative values. Gallier *et al.* [24] showed that the large positive values at high volume fractions are caused by wall effects in confined suspensions. Recently, Badia *et al.* [31] made an effort to develop a continuum model for suspension flow. In that work, they deduced material functions from simulation data in the literature, which include polynomial descriptions of the contact normal stresses and normal stress differences. These polynomials are a function of the normalized volume fraction $\varphi / \varphi m$ and provide a good fit to a set of data compiled from multiple literature sources [8,23,32], suggesting that the normal stresses scale with the normalized volume fraction.

*et al.*[19] using insights from the granular flow community. This description, which we will refer to as

*frictional rheology*, describes the suspension rheology in terms of the macroscopic friction coefficient $\mu $ and the viscous number $ I v$. Here, $\mu $ is defined as the ratio between the total shear stress and the wall-normal particle stress, $\mu = \sigma 12/ \sigma p , 22$, and $ I v$ is defined as the ratio of the viscous shear stress over the wall-normal particle stress, $ I v= \eta 0 \gamma \u02d9/ \sigma p , 22$. Using a custom pressure-imposed shear cell on a rheometer, Boyer

*et al.*[19] performed experiments measuring the contribution of the particles to the wall-normal stress. The shear cell was designed such that the fluid was able to pass through the top plate while keeping this top plate freely moving to impose a normal force to the suspension. Based on the experimental results, Boyer

*et al.*[19] proposed the following relation for the macroscopic friction coefficient:

*et al.*[23] took a different approach to obtain the hydrodynamic contribution to the macroscopic friction coefficient, based on the rheological model of Krieger–Dougherty [16], giving

*et al.*[27] proposed changes to the hydrodynamic component of the model by Gallier

*et al.*[23],

*et al.*[27] found very good agreement with existing data for an $ I v$ range of 5 decades. The literature shows that frictional rheology is suitable to describe the relation between the shear stress and the wall-normal particle stress for different suspension characteristics and particle properties, with slightly altered model coefficients. Nevertheless, the frictional rheology framework is, to our knowledge, not thoroughly tested against the normal stresses and stress differences that are needed for a full constitutive set of equations for suspension rheology.

In this paper, particle-resolved DNS is used to study the rheology of dense non-Brownian suspensions of neutrally buoyant spheres in a Newtonian carrier fluid in plane Couette flow. The Navier–Stokes equations are solved for the fluid phase and the Newton–Euler equations for the solid phase, where the coupling between the phases is numerically accomplished through a computationally efficient immersed boundary method (IBM) [34]. The aim of this study is to compare the viscous to the frictional rheology, i.e., the scaling of the rheology with $\varphi / \varphi m$ and the viscous number, respectively.

One of the main novelties of our work lies in the advanced numerical method used to investigate suspension rheology. We use particle-resolved DNS combined with a frictional soft-sphere collision model for particle contacts, which has been extensively validated against experimental data for oblique particle-wall and particle-particle collisions in viscous fluids [35]. Compared to the previous particle-resolved DNS study of Yeo and Maxey [22], our contact model includes tangential friction and stick/slip transition behavior. Gallier *et al.* [23] used a similar contact model compared to ours in their particle-resolved DNS study, but with a subtle difference in the definition of particle contact related to particle roughness, and their contact force expressions are different from our model except for the Coulomb friction force in the slip regime. Furthermore, we use a larger domain size in our simulations (domain height of 13.5 vs 10 particle diameters) and consider a wider range of bulk solid volume fraction ( $ \varphi b=0.1\u22120.6$).

In our study, we consider three different values for the Coulomb coefficient of sliding friction, $ \mu c=0$, $0.2$, and $0.39$, to investigate to what degree $\varphi / \varphi m$ and $ I v$ account for variations in interparticle friction. The results are used to perform a comprehensive comparison with previous work for the microstructure, stresses, shear rheology, and frictional rheology. To our knowledge, such comprehensive comparison has not been performed since the work of Gallier *et al.* [23], after which new models were introduced for the frictional rheology [27,33]. We consider the separate contributions from hydrodynamic and contact interactions to the rheology modeling and investigate in what way the solid friction coefficient affects these contributions. This shows the nuances of modeling the rheology for both the dilute regime, where hydrodynamics dominates, and the dense regime, where the contact interactions are most important. The presentation of our results as a function of both $\varphi / \varphi m$ and $ I v$ provides insight into the relation between the two rheology descriptions.

## II. METHODS

### A. Governing equations and contact model

*et al.*[35]. This contact model is based on a linear spring-dashpot model in which rigid particles are computationally allowed to slightly overlap. The normal contact force is determined using [35]

The contact model coefficients $ k n$, $ k t$, $ \eta n$, and $ \eta t$ depend on the dry normal ( $ e n$) and tangential ( $ e t$) coefficients of restitution and the contact time ( $ T c o n$) [35]. In the present DNS study, we have fixed the latter at $ e n=0.97$, $ e t=0.1$, and $ T c o n=4\Delta t$ with $\Delta t$ the computational time step. To improve the temporal resolution at which the lubrication and contact forces are resolved, substepping is employed [35] with the number of substeps within a computational time step set to 40 in the current simulations.

We have extensively validated our contact model against experimental data for oblique particle-wall and particle-particle collisions in viscous fluids in a previous study [35]. Our model is more advanced compared to the work of Yeo and Maxey [22] as it includes also tangential friction. Gallier *et al.* [23] used a similar contact model in the sense that their model does include tangential forces and stick–slip transition behavior. However, there is a subtle difference in the definition of particle contact related to particle roughness, and their contact force expressions are different from our model except for the Coulomb friction force in the slip regime. They define the start of particle contact by the moment when the distance between the centroids of two spheres drops below $d+ h r$, with $ h r$ the particle roughness height, while in our case, we assume particle contact to happen when the distance drops below $d$. In our model, the effect of particle roughness is included explicitly in the lubrication force corrections and implicitly in the parameters of the collision model, notably the Coulomb friction coefficient.

### B. Numerical method and computational setup

The governing equations for the fluid phase are solved on a fully staggered, Cartesian grid using a finite-volume/central-differencing approach combined with a predictor/correcter and third-order Runge–Kutta 3 (RK3) scheme for integration in time. This method is explained in detail by Breugem [34]. For the present simulations in the Stokes regime, the semi-implicit Crank–Nicolson scheme was implemented for integrating the diffusion term in Eq. (12b). An efficient IBM is used for the fluid/solid coupling, by means of locally adding forces to the Navier–Stokes equations near the boundary of a particle, to enforce Eq. (14) by good approximation. In the present simulations, this coupling is semi-implicit based on the approach of Tschisgale *et al.* [36].

The integration of the Newton–Euler equations for the solid phase is performed using the same RK3 scheme. For particles in close proximity, lubrication force corrections are used to account for subgrid hydrodynamic interactions, for which a contribution is added to the contact force. This lubrication force is capped for interparticle distances smaller than the particle roughness, which is assumed to be $0.001d$. The grid resolution used for the simulations is $d/\Delta x=16$, meaning there are 16 grid cells for the fluid phase over one particle diameter. This resolution was successfully used in the rheological study by Picano *et al.* [37] in the weakly inertial regime.

In this work, we consider neutrally buoyant spherical particles, $ \rho r=1$, that are monodisperse in size. The main parameters that we vary are the bulk volume fraction, which ranges from $ \varphi b=0.1$ to $ \varphi b=0.6$, and the Coulomb coefficient of sliding friction, which is set to $ \mu c=0$, $ \mu c=0.2$, and $ \mu c=0.39$. The value of $ \mu c=0.39$ is based on the value measured for polystyrene spheres in recent tilted-flume experiments by Shajahan *et al.* [38]. This gives a total of 24 cases, as specified in Table I. The same friction coefficient is used for interparticle and particle-wall contacts. The computational domain is a rectangular box with dimensions: $ L x= L z=13.5d$ and $ L y=27d$. A schematic of the computational domain is shown in Fig. 1. The domain is periodic in the x and y directions, while a no-slip/no-penetration boundary condition is imposed on the walls. The shear Reynolds number is set to $R e \gamma \u02d9=0.1$ for this study. This value is chosen to be in the Stokes regime while still being high enough to allow particles to move over sufficiently large distances during the course of a simulation such that the flow and spatial particle distribution become independent from the initial condition. The particles are initialized by a random positioning in the domain, while making sure there is no overlap between particles and of particles with the wall. For concentrations above 30%, this random initialization is not feasible anymore. Therefore, a body centered cubic (BCC) structure, with random perturbations in the particle positions, is used for the higher concentrations. The spacing of the BCC structure is adapted to fit the desired number of particles in the domain. The flow is initialized by imposing a linear velocity (plane Couette) profile in the streamwise direction over the entire domain. The initial translational and angular particle velocities are set to the local fluid velocity and half of the local fluid vorticity at the position of the particle centroid. Parameters for the computational setup, as well as derived quantities, are listed in Table I.

μ_{c}
. | ϕ_{b}
. | N_{p}
. | $ t t \gamma \u02d9$ . | $ t s \gamma \u02d9$ . | L_{core}/d
. | L_{wall}/d
. | ϕ_{c}
. | $ \gamma \u02d9 c/ \gamma \u02d9$ . |
---|---|---|---|---|---|---|---|---|

0 | 0.1 | 940 | 94 | 102 | 2.5 | 5.5 | 0.104 | 0.980 |

0 | 0.2 | 1880 | 94 | 102 | 2.5 | 5.5 | 0.208 | 0.980 |

0 | 0.3 | 2819 | 63 | 53 | 2.5 | 5.5 | 0.309 | 0.969 |

0 | 0.4 | 3759 | 94 | 61 | 2.5 | 5.5 | 0.410 | 0.901 |

0 | 0.45 | 4229 | 109 | 105 | 2.5 | 5.5 | 0.454 | 0.876 |

0 | 0.5 | 4699 | 125 | 89 | 2 | 5.75 | 0.486 | 0.831 |

0 | 0.55 | 5169 | 156 | 83 | 1.6 | 5.95 | 0.517 | 0.974 |

0 | 0.6 | 5639 | 125 | 114 | 1.6 | 5.95 | 0.565 | 0.939 |

0.2 | 0.1 | 940 | 94 | 94 | 2.5 | 5.5 | 0.109 | 0.963 |

0.2 | 0.2 | 1880 | 94 | 94 | 2.5 | 5.5 | 0.210 | 0.971 |

0.2 | 0.3 | 2819 | 63 | 61 | 2.5 | 5.5 | 0.312 | 0.954 |

0.2 | 0.4 | 3759 | 94 | 52 | 2.5 | 5.5 | 0.410 | 0.885 |

0.2 | 0.45 | 4229 | 109 | 97 | 2.5 | 5.5 | 0.458 | 0.817 |

0.2 | 0.5 | 4699 | 94 | 120 | 2 | 5.75 | 0.503 | 0.719 |

0.2 | 0.55 | 5169 | 125 | 84 | 1.6 | 5.95 | 0.504 | 0.834 |

0.2 | 0.6 | 5639 | 125 | 55 | 1.6 | 5.95 | 0.561 | 0.879 |

0.39 | 0.1 | 940 | 94 | 102 | 2.5 | 5.5 | 0.105 | 0.975 |

0.39 | 0.2 | 1880 | 94 | 102 | 2.5 | 5.5 | 0.210 | 0.971 |

0.39 | 0.3 | 2819 | 55 | 44 | 2.5 | 5.5 | 0.312 | 0.948 |

0.39 | 0.4 | 3759 | 63 | 55 | 2.5 | 5.5 | 0.414 | 0.860 |

0.39 | 0.45 | 4229 | 109 | 105 | 2.5 | 5.5 | 0.460 | 0.788 |

0.39 | 0.5 | 4699 | 94 | 66 | 1.8 | 5.85 | 0.510 | 0.644 |

0.39 | 0.55 | 5169 | 156 | 83 | 1.6 | 5.95 | 0.523 | 0.473 |

0.39 | 0.6 | 5639 | 188 | 45 | 1.6 | 5.95 | 0.580 | 0.771 |

μ_{c}
. | ϕ_{b}
. | N_{p}
. | $ t t \gamma \u02d9$ . | $ t s \gamma \u02d9$ . | L_{core}/d
. | L_{wall}/d
. | ϕ_{c}
. | $ \gamma \u02d9 c/ \gamma \u02d9$ . |
---|---|---|---|---|---|---|---|---|

0 | 0.1 | 940 | 94 | 102 | 2.5 | 5.5 | 0.104 | 0.980 |

0 | 0.2 | 1880 | 94 | 102 | 2.5 | 5.5 | 0.208 | 0.980 |

0 | 0.3 | 2819 | 63 | 53 | 2.5 | 5.5 | 0.309 | 0.969 |

0 | 0.4 | 3759 | 94 | 61 | 2.5 | 5.5 | 0.410 | 0.901 |

0 | 0.45 | 4229 | 109 | 105 | 2.5 | 5.5 | 0.454 | 0.876 |

0 | 0.5 | 4699 | 125 | 89 | 2 | 5.75 | 0.486 | 0.831 |

0 | 0.55 | 5169 | 156 | 83 | 1.6 | 5.95 | 0.517 | 0.974 |

0 | 0.6 | 5639 | 125 | 114 | 1.6 | 5.95 | 0.565 | 0.939 |

0.2 | 0.1 | 940 | 94 | 94 | 2.5 | 5.5 | 0.109 | 0.963 |

0.2 | 0.2 | 1880 | 94 | 94 | 2.5 | 5.5 | 0.210 | 0.971 |

0.2 | 0.3 | 2819 | 63 | 61 | 2.5 | 5.5 | 0.312 | 0.954 |

0.2 | 0.4 | 3759 | 94 | 52 | 2.5 | 5.5 | 0.410 | 0.885 |

0.2 | 0.45 | 4229 | 109 | 97 | 2.5 | 5.5 | 0.458 | 0.817 |

0.2 | 0.5 | 4699 | 94 | 120 | 2 | 5.75 | 0.503 | 0.719 |

0.2 | 0.55 | 5169 | 125 | 84 | 1.6 | 5.95 | 0.504 | 0.834 |

0.2 | 0.6 | 5639 | 125 | 55 | 1.6 | 5.95 | 0.561 | 0.879 |

0.39 | 0.1 | 940 | 94 | 102 | 2.5 | 5.5 | 0.105 | 0.975 |

0.39 | 0.2 | 1880 | 94 | 102 | 2.5 | 5.5 | 0.210 | 0.971 |

0.39 | 0.3 | 2819 | 55 | 44 | 2.5 | 5.5 | 0.312 | 0.948 |

0.39 | 0.4 | 3759 | 63 | 55 | 2.5 | 5.5 | 0.414 | 0.860 |

0.39 | 0.45 | 4229 | 109 | 105 | 2.5 | 5.5 | 0.460 | 0.788 |

0.39 | 0.5 | 4699 | 94 | 66 | 1.8 | 5.85 | 0.510 | 0.644 |

0.39 | 0.55 | 5169 | 156 | 83 | 1.6 | 5.95 | 0.523 | 0.473 |

0.39 | 0.6 | 5639 | 188 | 45 | 1.6 | 5.95 | 0.580 | 0.771 |

### C. Statistical averaging

*average*stress within a particle and may, therefore, deviate from the

*actual*local stress within a particle when the particle stress distribution is strongly

*heterogeneous*. As a result, the particle stress profiles may show artificial fluctuations near the wall where the gradients in the solid volume concentration are large due to particle layering. Nonetheless, in the present study, we focus on the rheology in the

*homogeneous*core region of the channel for which Eq. (20) is a good approximation.

## III. RESULTS

The best way to get a first insight into the suspension behavior is to look at the particle positions and velocities in the domain. A snapshot of the simulations with $ \varphi b=0.5$ is shown in Fig. 2 for a time instance at which the flow is fully developed. The figure shows the spatial particle distribution and particle velocities in a 3D-slab and two perpendicular cross-sections of the flow domain, for a friction coefficient of $ \mu c=0$ (a) and $ \mu c=0.39$ (b). As is visible from the colors, the particles at the top of the domain move in the positive y-direction with the wall velocity and the particles at the bottom move the opposite way with the same velocity. The particles in the interior show a smooth velocity gradient from the bottom to the top of the domain with a velocity close to zero at the centerline. Figure 2(a), corresponding to a friction coefficient of $ \mu c=0$, shows that there is a clear particle layering along the walls that extends into the domain over a distance of approximately 4 particle diameters, after which the layering is less distinct. Furthermore, outside this region of particle layering, the particles align in trains that are oriented approximately from the top left toward the bottom right of the domain. This alignment is along the compressional axis of the flow which is the direction in which the most particle contacts are expected for a random particle distribution under shear. The addition of solid friction, shown in Fig. 2(b), results in less distinct wall-induced layering that extends into the domain over a distance of approximately 2 particle diameters. The alignment of particles along the compressional axis is, therefore, also more clearly visible, as the region of the domain where this occurs is larger in this case.

Figure 3 shows a snapshot of the particles in the flow for the case with $ \varphi b=0.6$ and $ \mu c=0.39$. This case shows a pronounced wall-induced layering that extends throughout the entire domain, where the layers of particles at the wall are very neatly aligned in the flow directions. Around the centerline of the channel the layers are somewhat distorted but remain clearly visible. The alignment of the particles in the streamwise direction differs strongly from the case with $ \varphi b=0.5$ where the alignment of the particles is along the compressional axis of the flow, indicating a change in dynamics when the solid volume fraction is high enough for the wall-induced ordering to span the entire domain. Moreover, looking at the xz-plane shown in the right pane of Fig. 3, additional structuring of the particles is visible. In this plane, the particles are arranged in a hexagonal packing, which is clearly visible in the entire domain.

As explained in Sec. II C on statistical averaging, the results of the simulations are plane-averaged in planes parallel to the channels walls. This leads to statistics as a function of the wall-normal coordinate in the channel, which are then time-averaged over the simulated time period in which the flow is fully developed. The resulting concentration profiles are a good way to get an overview of the effect of volume fraction and solid friction on the particle distribution in the suspensions. Figure 4 shows the concentration profiles for all cases. These profiles affirm the observation we made in Fig. 2 that there is particle layering near the wall. The layering is more severe with higher bulk volume fraction and extends into the domain with multiple clear layers for $ \varphi b>0.4$. For the cases with a concentration of $ \varphi b\u22650.55$, the layering even extends throughout the entire domain, as was also observed in Fig. 3 for the highest concentration. In addition, the concentration profiles show that the wall-induced layering changes from $ \varphi b=0.55$ to $ \varphi b=0.6$, for $ \mu c=0.39$. This can be seen from the number of peaks in the profiles, equal to 15 and 16 peaks (versus 23 and 24 at the start of the simulation) for $ \varphi b=0.55$ and $ \varphi b=0.6$, respectively. This is in line with the hexagonal packing that was observed in Fig. 3. The change in layering occurs at a lower volume fraction, between $ \varphi b=0.5$ and $ \varphi b=0.55$, for the suspensions with a lower friction coefficient, showing that the layering is partially mitigated by the increase of solid friction. The increase in friction coefficient also leads to less distinct peaks near the wall, which was previously observed by Gallier *et al.* [24].

Next to the concentration profiles, we also look at the velocity profiles. The mean intrinsic fluid velocity normalized by the wall velocity is shown in Fig. 5 for all three friction coefficients. The fluid velocity for the cases without friction shows a slight dependence on the solid volume fraction. For increasing volume fraction, the velocity gradient or shear rate decreases for the region around the center of the channel. This effect is stronger with increasing friction coefficient and clearly visible at $ \mu c=0.39$. There is a minimum in the shear rate as a function of the bulk volume fraction after which the shear rate increases again. This increase corresponds to the cases with layering throughout the entire domain, showing that this also affects the fluid velocity profile.

In Fig. 6, we compare the mean intrinsic fluid and particle velocity profiles for the cases with $ \varphi b=0.4$ and for friction coefficients of $ \mu c=0$ and $ \mu c=0.39$. While in the core region, the particle and fluid velocity are almost exactly equal, macroscopic particle/fluid slip is observed in the wall region, which is more clearly illustrated in the inset. At the wall ( $z/d=0$) the fluid moves with the velocity of the wall due to the no-slip condition that is enforced there. The particle phase does not have to obey a no-slip condition and shows a slightly lower velocity at the wall. Moving away from the wall the profiles fluctuate somewhat until a position of $z\u22481.25d$ where the velocity profiles coincide again, showing that the slip only occurs in a thin layer near the wall with a thickness on the order of a particle diameter. Even though the velocity profiles change with friction coefficient, the slip velocity is very similar as well as the width of the layer where the slip occurs for both friction coefficients shown in Fig. 6.

From the concentration and velocity profiles, we conclude that there is a distinct difference in suspension dynamics between the wall region and the core of the channel. The first notable difference is the wall-induced layering, which originates from the constraint that particles cannot overlap with the flat solid walls and which is amplified by the perfect monodispersity of the spheres. Second, there is significant macroscopic particle-fluid slip in the wall region. This is caused by the no-slip condition for the fluid, while the particles are allowed to slip. The choice of monodisperse spheres and flat walls in the present study was motivated by the simplicity of the system and of the related implementation in our in-house particle-resolved DNS code. In the analysis of the suspension rheology, we carefully distinguished between the wall-affected layers on the one hand and the core region with a homogeneous mean concentration and linear mean velocity profile on the other hand. Furthermore, as wall-induced layering becomes progressively more important for increasing bulk concentration, we ignored the cases with the highest bulk concentrations where a homogeneous concentration in the core is absent. In the literature, different methods have been used to counteract the tendency of particle layering along the walls, such as (1) using a bidisperse or polydisperse particle size distribution as these are naturally less prone to ordering [39,40], (2) exchanging the flat walls for walls with a roughness on the order of the particle size, which is widely used in suspension rheometry [41–43] and also applied in numerical investigations [44–46], and (3) using Lees–Edwards (periodic and sliding) boundary conditions [47–49] for the particles, which removes the constraint that particles cannot overlap with the walls. The use of such methods is recommended for future research, in particular, to extend the currently explored range of bulk concentrations in our rheology analysis toward the jamming limit.

The width of the core region, $ L c o r e/d$, is defined around the centerline with a value set to $ L c o r e/d=2.5$ for the cases up to $ \varphi b=0.45$, while set somewhat smaller for the cases with $ \varphi b\u22650.5$, see Table I. We would like to stress the fact that splitting the domain in a wall and a core region makes it possible to obtain an accurate description of the suspension rheology in the statistically homogeneous core with a uniform concentration profile. The cases with $ \varphi b\u22650.55$ for which the wall-layering extends throughout the entire domain are excluded from the analysis and model fits. Figures 4 and 5 show that the concentration and shear rate in the core region differ slightly from the bulk quantities. For the core, we, therefore, determined the average shear rate ( $ \gamma \u02d9 c$) and average concentration ( $ \varphi c$), of which the values are also listed in Table I. The cases with $ \varphi b=0.4$ will be used as an example throughout this work because these cases nicely illustrate the suspension in the dense regime.

*et al.*[51], obtained using a Couette rheometer. Our results are also in good agreement with the Stokesian Dynamics simulations of Foss and Brady [52], for their high Péclet cases, and of Sierou and Brady [21]. From the figures, it becomes clear that fewer particles reside near the reference particle in the extensional flow direction, in the so-called depletion zones. The size of these zones becomes smaller with increasing volume fraction and shifts to a position closer to 45 $ \xb0$ from the positive y-axis. Furthermore, the spots of high probability at $r\u22481d$ become more prominent with higher volume fraction and shift toward the y-axis, indicating more frequent contact of particles in this direction. This corresponds to the snapshots of the spatial particle distribution and the concentration profiles shown before. For the dense cases, secondary rings of higher than average probability appear at a position of $r\u22482d$ corresponding to multiparticle interactions. These secondary rings are more clearly visible in the compressional quadrants of the figure, corresponding to the particle trains that we observed in Fig. 2. Additionally, in the cases without friction horizontal streaks of higher probability are visible for $ \varphi b>0.3$ showing the preference of particles aligning parallel to the wall. It is important to note that this effect is much stronger in the cases without friction, showing that solid friction decreases the tendency for wall-induced layering in dense suspensions.

A convenient way to compare the PPDFs of different cases is to radially average over a thin shell between $1<r/d<1.05$ around the particle as a function of the position on the particle. This gives a good indication of the position on the particle where close interactions with neighboring particles occur. Figure 8(a) shows the PPDFs in this thin shell for the same cases that are shown in Fig. 7. The angle $\theta $ is defined clockwise from the negative y-axis as depicted in the figure. The first feature that stands out is that the lower theta-bound of the depletion zones ( $g<1$) coincides for all the volume fractions investigated, with and without friction, and is located at approximately $3\pi /4$ and $7\pi /4$, respectively. The width of this depletion zone decreases with increasing volume fraction but seems insensitive to the friction coefficient. Following this depletion zone, a peak develops at $\theta =\pi $ with increasing volume fraction. This peak indicates the alignment of particles parallel to the walls. This alignment is stronger for the cases without friction, which is evident from the larger peak amplitude. The plateaus that are visible for all cases correspond to the compressional quadrants of the shear flow. In this plateau, the probability increases with increasing volume fraction and is mostly insensitive to the friction coefficient. Figure 8(b) shows the PPDF in the shell for the cases with $ \varphi b=0.4$ over a range of $0<\theta <\pi $. The numerical results of Gallier *et al.* [23] and Yeo and Maxey [22] are shown for comparison. Our data show very good agreement with the data from the literature.

*et al.*[23], as well as the experimental data from Zarraga

*et al.*[17] and Dbouk

*et al.*[18]. We observe good agreement between our simulations and the results from the literature up to volume fractions of $ \varphi b=0.5$; however, the experimental results show consistently higher relative viscosity over the entire range. The dashed lines show the relative viscosity derived from the frictional rheology model and will be discussed later. The analytical result of Batchelor and Green [13], valid for concentrations up to $\varphi \u22480.1$ and defined in Eq. (4), is shown for reference as the dashed-dotted line. We observe that our results at $ \varphi b=0.2$ are only slightly higher than predicted from this relation, suggesting a still minor role of contacts at this concentration. From Fig. 11(a) we observe that the data does not collapse onto a single curve. A clear influence is present of the friction coefficient, in line with Eq. (1) based on dimensional analysis. It is only for the lower concentration range, where contacts play a minor role, that the data seem to collapse. For the cases with $ \varphi b\u22650.55$, we observe a deviation from the trend with a decrease in relative viscosity. This deviation can be explained by the layering of particles that we observe for these cases. The particle layering makes it easier for particles to slide past each other, decreasing the resistance to shearing motion and, therefore, decreasing the relative viscosity. The deviation from the trend at high volume fraction is even more visible in Fig. 11(b). This figure shows the relative viscosity as a function the volume fraction normalized by the maximum flowable packing fraction $ \varphi m$. The maximum flowable packing fraction is determined by fitting the Maron–Pierce correlation [15], given by Eq. (5), to the data, excluding the cases with layering throughout the domain. All data collapse to a single curve, showing that the effect of solid friction can be accounted for using the maximum flowable packing fraction. The values of the maximum flowable packing fraction that we find are $ \varphi m=0.70$, $ \varphi m=0.66$, and $ \varphi m=0.64$ for the cases with $ \mu c=0$, $ \mu c=0.2$, and $ \mu c=0.39$, respectively. These values are somewhat higher than those of the experimental datasets from the literature but are in line with the values reported in the numerical studies of Yeo & Maxey [22] and Gallier

*et al.*[23]. It is worthwhile to note that the values for $ \varphi m$ do not change significantly when the cases with $ \varphi b=0.5$ are excluded from the fit, indicating that the core region of these cases is not affected by wall-induced layering.

Figure 12 shows the decomposition of the relative viscosity into the hydrodynamic and contact components, where the hydrodynamic part is the sum of the contribution from the viscous stress and the hydrodynamic contribution to the particle stress. The lines show the relative viscosity derived from the frictional rheology model and will be discussed later. In Fig. 12(a), it can be seen that for low volume fractions the hydrodynamic component accounts for the largest part of the viscosity. With increasing volume fraction the contact part becomes more important and accounts for the strong diverging nature close to the maximum flowable packing fraction. The effect of friction is not visible in the hydrodynamic part, for which the data overlap in the range of volume fractions considered here. For the contact part, it becomes clear that friction causes a divergence of the relative viscosity at a lower volume fraction. For the frictionless case, it is interesting to note that the hydrodynamic contribution is the largest over the entire range of volume fractions investigated. On the contrary, for the cases with friction the contact contribution exceeds the hydrodynamic for $ \varphi c\u22730.45$. This figure also shows that the breakdown of the diverging trend for the high volume fractions originates primarily from the contact contribution, which decreases for the cases with pronounced layering in the core region. We again included the analytical expression of Batchelor and Green [13], Eq. (4), for reference as the dashed-dotted line. The close agreement between the our data and this relation at $ \varphi b\u22640.2$ indicates that at low volume fractions the hydrodynamic contribution can be described using this relation. Additionally, this close agreement suggests that the hydrodynamic contribution scales with $ \varphi c$ rather than $ \varphi c/ \varphi m$. Figure 12(b) shows the decomposition of the relative viscosity as a function of the normalized volume fraction. The normalization causes a collapse of the contact part up to $ \varphi c/ \varphi m\u22480.7$, while the effect of the friction coefficient now becomes apparent in the hydrodynamic part of the relative viscosity. This shows that although the relative viscosity collapses with $ \varphi c/ \varphi m$, the separate contributions do not completely collapse.

*et al.*[23] is the closest match. The deviation of our data from the models might be caused by the choice for the particle pressure when calculating the viscous number and macroscopic friction coefficient. In the experiments of Boyer

*et al.*[19], a constant particle pressure is imposed; however, it is not entirely clear whether this corresponds to the full normal particle stress or only the contact contribution. Gallier

*et al.*[23] noticed that the agreement between their simulations and the model by Boyer

*et al.*[19] is better when using only the contact contribution to the normal stress. To check the effect on the frictional rheology modeling, we tried using the complete wall-normal stress, the wall-normal particle stress, and the contact contribution to the wall-normal particle stress. We decided to use the wall-normal particle stress, based on the results shown in the Appendix. Finally, we note that the deviation seen for high viscous number from the model of Gallier

*et al.*[23] might also be related to subtle differences in the contact modeling; both $\mu $ and $ I v$ are normalized with the wall-normal particle stress, which becomes small for high $ I v$ and, hence, is sensitive to how particle contact is dealt with.

*et al.*[23] matches the closest to our data; however, there is still a slight deviation that increases with increasing friction coefficient and is more visible toward the high and low end of the viscous number range. The relations used by Lecampion and Garagash [33], and Chèvremont

*et al.*[27] are identical to the original model proposed by Boyer

*et al.*[19], defined in Eq. (8). To model our results, we start with this relation between the viscous number and the solid volume fraction, which we determine by fitting the relation used by Gallier

*et al.*[23] to our data

. | μ_{c} = 0
. | μ_{c} = 0.2
. | μ_{c} = 0.39
. | Boyer et al. [19]
. | Gallier et al. [23]
. | Chèvremont et al. [27]
. |
---|---|---|---|---|---|---|

ϕ_{m} | 0.72 | 0.69 | 0.67 | 0.585 | 0.64 | f (μ_{c}) |

n | 0.3664 | 0.3421 | 0.3309 | 0.5 | 0.4 | 0.5 |

[η] | 3.53 | 3.27 | 3.07 | 2.5 | 2.4 | 2.5 |

μ_{1} | 0.2 | 0.3 | 0.36 | 0.32 | 0.32 | 0.36 |

μ_{2} | 0.71 | 0.88 | 1.04 | 0.7 | 0.7 | 0.7 |

I_{0} | 0.0194 | 0.0368 | 0.0584 | 0.005 | 0.005 | 0.0133 |

. | μ_{c} = 0
. | μ_{c} = 0.2
. | μ_{c} = 0.39
. | Boyer et al. [19]
. | Gallier et al. [23]
. | Chèvremont et al. [27]
. |
---|---|---|---|---|---|---|

ϕ_{m} | 0.72 | 0.69 | 0.67 | 0.585 | 0.64 | f (μ_{c}) |

n | 0.3664 | 0.3421 | 0.3309 | 0.5 | 0.4 | 0.5 |

[η] | 3.53 | 3.27 | 3.07 | 2.5 | 2.4 | 2.5 |

μ_{1} | 0.2 | 0.3 | 0.36 | 0.32 | 0.32 | 0.36 |

μ_{2} | 0.71 | 0.88 | 1.04 | 0.7 | 0.7 | 0.7 |

I_{0} | 0.0194 | 0.0368 | 0.0584 | 0.005 | 0.005 | 0.0133 |

*et al.*[23] shows the best agreement to our data, especially for low $ I v$. With increasing viscous number, we see that the hydrodynamic contribution is underestimated by all models. Using the model of Gallier

*et al.*[23], as defined in Eq. (9), we can fit the hydrodynamic contribution to the macroscopic friction coefficient $ \mu h$. Here, $ \varphi m$ and $n$ are already determined, leaving $[\eta ]$ as the only fitting coefficient. This fit gives the values of $[\eta ]$ for all three solid friction coefficients, which are also listed in Table II. The colored solid lines in Fig. 13(c) show that these fits provide a very good description of our data. It is interesting to note that the effect of friction plays a role in the hydrodynamic contribution at low viscous numbers, showing that the hydrodynamic part of the particle shear stress decreases with increasing friction coefficient. At high viscous numbers, the fits overlap and show that solid friction does not affect the hydrodynamic contribution in the dilute limit, as expected. The contact part $ \mu c$ is shown in Fig. 13(d), which only shows weak dependence on the viscous number over the range investigated in this work. However, the plateau that exists at high viscous number depends on the Coulomb coefficient of sliding friction $ \mu c$. This contradicts the results of Chèvremont

*et al.*[27], who found a plateau value of 0.7 irrespective of the solid friction. They did observe a deviation from the master curve for cases with a low $ \mu c$ at lower viscous numbers. There is a limited agreement between our data and the models, which seem to have the right shape but not the right coefficients to describe our data. Our data on the contact contribution show a dependence on the solid friction coefficient, which was also found by Chèvremont

*et al.*[27] but dismissed in their modeling. To be able to fit the contact contribution defined as

*et al.*[27] by looking at the low $ I v$ plateau in the contact contribution. The values we chose are $ \mu 1=0.2$, $ \mu 1=0.3$, and $ \mu 1=0.36$ for the cases with $ \mu c=0$, $ \mu c=0.2$, and $ \mu c=0.39$, respectively. The resulting fits are shown in Fig. 13(d), and the fitting coefficients $ \mu 2$ and $ I 0$ are listed in Table II. The figure shows that Eq. (29) is able to describe our data very well with fitting coefficients that are of the same order as those found in the literature. Adding the two contributions together provides a description of the complete macroscopic friction coefficient as shown in Fig. 13(a). The colored solid lines are the result of fitting the separate contributions to the contact and hydrodynamic parts, which shows very good agreement to the data. The model, which is essentially the one proposed by Boyer

*et al.*[19] with adjusted coefficients, only shows dependence on the solid friction coefficient in the limit of low viscous numbers. The decomposition of $\mu $ shows that the effect of friction is present in both the hydrodynamic and contact contributions, whereas the total does not show this dependence due to compensating effects of friction that cancel each other.

The relative viscosity derived from the frictional rheology model can also be decomposed into two separate contributions, with the first two terms in Eq. (30) accounting for the contacts while the final two terms account for the hydrodynamic contribution. The result of that decomposition is shown as the lines in Fig. 12. The good agreement between the model and the data show that the frictional rheology and viscous rheology are equivalent with regard to the relative shear viscosity; since $ \eta r=\mu / I v$ and $\varphi / \varphi m$ is a monotonically decreasing function in $ I v$, which is only weakly affected by $ \mu c$, the two frameworks are actually interchangeable. Moreover, the decomposition shows that the hydrodynamic contribution can be described using Eq. (4) for low volume fractions. Nevertheless, a more detailed investigation on the effect of the Coulomb coefficient of sliding friction on the frictional rheology model is needed. Specifically, investigating how the fitting coefficients change for different Coulomb coefficient of sliding friction.

*et al.*[31] who used these correlations to describe the normal contact stresses in a continuum modeling approach for the suspension rheology. The correlations are defined as

*et al.*[31] to fit data from the literature. There is some deviation from the coefficients found by Badia

*et al.*[31], which might be caused by the fact that we use the total particle stresses, whereas the original authors used the correlations for the contact part of the stresses. This should only amount to a slight difference as we have seen that, for $ \varphi b=0.4$, the particle stress is almost completely governed by the contact contribution. Figure 14(d) shows the normal stress ratios $ \lambda 2= \u27e8 \sigma p , 22 \u27e9 \xaf/ \u27e8 \sigma p , 11 \u27e9 \xaf$ and $ \lambda 3= \u27e8 \sigma p , 33 \u27e9 \xaf/ \u27e8 \sigma p , 11 \u27e9 \xaf$ as a function of the normalized volume fraction. This figure shows that the normal stress ratios are similar for all friction coefficients, but do not completely collapse on a single curve. Additionally, $ \lambda 2$ shows an increase with increasing volume fraction for all three friction coefficients, whereas $ \lambda 3$ is more or less constant around a value of 0.5 with a slight increase at high volume fractions.

. | μ_{c} = 0
. | μ_{c} = 0.2
. | μ_{c} = 0.39
. | Badia et al. [31]
. |
---|---|---|---|---|

a_{1} | −2.3022 | −2.3305 | −2.3648 | −2.4247 |

a_{2} | 2.7915 | 2.9724 | 3.0688 | 4.1280 |

b_{1} | 0.9717 | 0.5938 | 0.3792 | 0.3750 |

b_{2} | −2.0377 | −0.2760 | 0.4415 | 0.0366 |

b_{3} | 2.5687 | 0.9074 | 0.3524 | 0.4846 |

c_{1} | 0.5895 | 0.4977 | 0.3361 | 2.1446 |

c_{2} | −0.5932 | −0.3755 | −0.0983 | −2.7234 |

c_{3} | 1.8373 | 1.2661 | 1.0189 | 1.5759 |

. | μ_{c} = 0
. | μ_{c} = 0.2
. | μ_{c} = 0.39
. | Badia et al. [31]
. |
---|---|---|---|---|

a_{1} | −2.3022 | −2.3305 | −2.3648 | −2.4247 |

a_{2} | 2.7915 | 2.9724 | 3.0688 | 4.1280 |

b_{1} | 0.9717 | 0.5938 | 0.3792 | 0.3750 |

b_{2} | −2.0377 | −0.2760 | 0.4415 | 0.0366 |

b_{3} | 2.5687 | 0.9074 | 0.3524 | 0.4846 |

c_{1} | 0.5895 | 0.4977 | 0.3361 | 2.1446 |

c_{2} | −0.5932 | −0.3755 | −0.0983 | −2.7234 |

c_{3} | 1.8373 | 1.2661 | 1.0189 | 1.5759 |

*et al.*[31] represented by the solid black line. The shift toward positive values is also present for the literature data and the model; however, the strong increase in our data is caused by the particle layering in the domain at high volume fractions affecting $ \u27e8 \sigma p , 11 \u27e9 \xaf$. The first normal stress difference computed from the normal stress fits, shown by the dashed lines, shows good agreement with our data and clearly highlights the effect of friction. With increasing Coulomb friction coefficient, $ N 1$ becomes less negative and the upswing toward positive values is less strong. The second normal stress difference, shown in Fig. 15(b), shows a clear trend toward negative values for high volume fractions. The effect of friction is clearly visible, resulting in a more negative $ N 2$ with increasing friction coefficient, which is stronger at higher volume fractions. This is in good agreement with the data from the literature [23]. Furthermore, the data for the cases with $ \mu c=0.39$ coincide with the polynomial description of Badia

*et al.*[31], shown by the solid black line, up to the point where the model shows an upswing. The volume fraction at which this happens corresponds to our cases exhibiting a strong layering effect, preventing a good comparison to the model this close to the maximum packing fraction. The second normal stress difference, computed from the normal stress fits and shown by the dashed lines, also shows good agreement with our data. The upswing in the model at high volume fractions is outside the range of our data; however, this upswing is similar to that of the fit of Badia

*et al.*[31] and shifts to lower volume fraction for lower friction coefficient.

The viscous and frictional rheology frameworks fail to collapse the data to a master curve for the normal stresses and, therefore, the normal stress differences, as observed in Fig. 15 and noting from Fig. 13(b) that the monotonic relation between $\varphi / \varphi m$ and $ I v$ depends only weakly on the friction coefficient. Considering the normalized wall-normal particle stress component $ \u27e8 \sigma p , 22 \u27e9 \xaf/ \u27e8 \sigma 12 \u27e9 \xaf$, we do observe from Fig. 14 a well-nigh collapse as a function of normalized volume fraction, $ \varphi c/ \varphi m$, for the three friction coefficients. In passing, we remark that this is consistent with the collapse for $\mu $ as a function of $ I v$ in the frictional rheology framework, as $\mu = \u27e8 \sigma 12 \u27e9 \xaf/ \u27e8 \sigma p , 22 \u27e9 \xaf$ is actually the inverse of the normalized wall-normal particle stress. The other two normal stresses do not collapse as a function of normalized volume fraction, with $ \u27e8 \sigma p , 33 \u27e9 \xaf$ showing the largest effect of the Coulomb friction coefficient. The effect of friction on the normal stress differences we observe matches with data from the literature [23,26], where there is also no collapse with maximum flowable packing fraction. To complete the description of the suspension rheology, the effect of friction on the normal stress differences needs to be studied in more detail.

## IV. CONCLUSIONS

The particle-resolved simulations of dense suspensions in plane Couette flow performed in this study give a detailed insight into the microstructure and stresses that occur in these suspensions. The effects of solid volume fraction and Coulomb coefficient of sliding friction were shown while comparing the data with experimental and numerical data from the literature. Moreover, we commented on the use of the maximum flowable packing fraction to account for particle characteristics and compared this approach to the frictional rheology framework.

From our results, we can conclude that wall-induced particle layering extends into the entire domain for the cases with volume fractions of $ \varphi b\u22650.55$. This layering affects the rheology by changing the microstructure in a way that particles can slide over each other more easily, lowering the relative viscosity and breaking the diverging trend of the suspension viscosity with increasing concentration. We observed that increasing the solid friction coefficient partially mitigates the layering, which we can see from both the concentration profiles and the PPDFs.

The total shear stress increases with increasing solid volume fraction and increases with increasing friction coefficient. Considering the separate contributions to the total shear stress, we can see that the particle stress accounts for the largest part. This particle stress consists of two parts of which the hydrodynamic part is larger for the low volume fractions, while the contact part dominates for higher volume fractions. The volume fraction at which the contact part becomes larger than the hydrodynamic part decreases with increasing friction coefficient. This can also be seen in the decomposition of the relative viscosity, where the divergence at the lower volume fraction is caused by the increase in the contact particle stress. The relative viscosity, which is in good agreement with the literature, can be collapsed using the maximum flowable packing fraction determined by fitting the Maron–Pierce correlation [15].

The normal stresses show anisotropy leading to nonzero normal stress differences. The normal particle stresses can be described using polynomial fits used by Badia *et al.* [31] and only show a collapse in $ \varphi c/ \varphi m$ for the wall-normal direction. The normal stress differences show good agreement with the data from the literature, where $ N 1$ is slightly negative for volume fractions up to $ \varphi c/ \varphi m\u22480.8$ after which it becomes positive with a large magnitude due to pronounced wall-induced layering. Increasing the friction coefficient decreases the magnitude of $ N 1$ as was also observed in the literature. $ N 2$ is negative over the entire range investigated here and increases in magnitude with volume fraction. The effect of friction is opposite for $ N 2$, increasing the magnitude with increasing friction coefficient.

Our data can be described using the frictional rheology model from Gallier *et al.* [23] with adapted fitting coefficients. Additionally, the solid friction coefficient should be taken into account when using the frictional rheology model, and further investigation is needed to incorporate this into the model. The effect of friction is visible when the macroscopic friction coefficient is decomposed into the hydrodynamic and contact parts. These are modeled individually and show that friction slightly lowers the hydrodynamic contribution while it increases the contact contribution. Therefore, the total macroscopic friction coefficient shows the collapse of the data, due to compensating effects of friction on the two individual contributions, with a slight deviation at very low viscous numbers. The relative viscosity derived from the frictional rheology model gives a very good agreement with our DNS data. This indicates that the viscous and frictional rheology frameworks are largely equivalent, primarily because the relative volume fraction $ \varphi c/ \varphi m$ is a monotonic function of $ I v$ and depends only weakly on the friction coefficient. Both frameworks capture the effect of solid friction on the shear rheology. This is, however, not the case for the streamwise and spanwise normal particle stresses as well as the normal particle stress differences, which still show a significant dependency on the friction coefficient that is not yet fully understood.

## ACKNOWLEDGMENTS

The authors gratefully acknowledge Delft University of Technology, the funding of this project under a TU Delft/TNO agreement, and the support of SURF Cooperative for the use of the Dutch national e-infrastructure (Snellius) under Computing Grant No. 2020.036 from the Dutch Research Council (NWO).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are openly available in 4TU. ResearchData repository at https://doi.org/10.4121/bbb1d553-7230-4a60-80a2-23d5697d89a5 [55].

### APPENDIX: ON THE NORMALISATION WITH THE NORMAL STRESS IN THE FRICTIONAL RHEOLOGY FRAMEWORK

The deviation of our DNS results for the frictional rheology from models proposed in the literature might be partially explained by the different choices made in the literature for the normalization with the wall-normal particle stress. In the experiments of Boyer *et al.* [19], a constant particle pressure is imposed; however, it is not entirely clear whether this corresponds to the full normal particle stress or only the contact contribution. Gallier *et al.* [23] noticed that the agreement between their simulations and the model by Boyer *et al.* [19] is better when using only the contact contribution to the normal stress. To check the effect on the frictional rheology modeling we tried using the full wall-normal stress [given by Eq. (24)], the full wall-normal particle stress, and the contact contribution to the wall-normal particle stress. Figure 16(a) shows the results for the volume fraction in the core region, $ \varphi c$, as a function of viscous number based on the three different choices for the wall-normal stress. The colors represent the different stress contributions and different symbols are used to distinguish between frictionless and frictional particles. The figure shows that all stress contributions give comparable results. The normalization of the viscous number with the full wall-normal stress shows the closest match with the model of Boyer *et al.* [19]; however, toward low viscous number the deviation increases. Using the wall-normal particle stress results in a relatively good agreement with the model of Gallier *et al.* [23] over the entire range of viscous numbers investigated. The normalization of the viscous number with the contact contribution to the wall-normal particle stress shows the largest deviation from the models at both high and low viscous numbers. At the low viscous number, the deviation might be caused by the pronounced wall-induced layering.

Figure 16(b) shows the macroscopic friction coefficient as a function of the viscous number, where both $\mu $ and $ I v$ are defined based on either the full wall-normal stress, the full wall-normal particle stress, or the contact contribution to the wall-normal particle stress. Here, the contributions again give similar results. The particle stress and the contact contribution follow the same curve except at very low $ I v$. At the high viscous number, the best match with the models is found when using the full wall-normal stress.

Our comparison study shows that the wall-normal stress contribution that is used to determine the frictional rheology has an effect on the results for the frictional rheology, although this effect is limited. Since the flow dynamics in the DNS depends on the pressure gradient instead of pressure itself, and the mean value of the pressure at the walls is arbitrarily set to zero in our simulations for reference, we have chosen not to use the full wall-normal stress but the wall-normal particle stress in the definition of $\mu $ and $ I v$.

## REFERENCES

*A Physical Introduction to Suspension Dynamics*