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, η r, as a function of the reduced solid volume fraction, ϕ / ϕ m, with ϕ m the maximum flowable packing fraction, and (2) the “frictional” framework in terms of a macroscopic friction coefficient, μ, 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 ϕ b = 0.1 to 0.6 and considered three different Coulomb friction coefficients: μ c = 0, 0.2, and 0.39. We find that η r scales well with ϕ / ϕ m, with ϕ m obtained from fitting the Maron–Pierce correlation. We also find that μ scales well with I v. Furthermore, we find a monotonic relation between ϕ / ϕ m and I v, which depends only weakly on μ c. Since η r = μ / 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 ϕ / ϕ m or I v, remain explicitly dependent on μ c in a manner that is not yet fully understood.

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.

The rheology of suspensions has been studied extensively for many years since pioneering work by Albert Einstein [1,2]. A good overview of the subject is given in several review articles [3–5]. Suspensions show an increase in mixture viscosity, with respect to the viscosity of the carrier fluid, for increasing solid volume fraction. Dense suspensions can exhibit shear rate dependent behavior such as yield stresses [6], shear thinning [7,8], and (dis)continuous shear thickening [9–12]. The focus of the present study is on the rheology of monodisperse suspensions of spheres in a Newtonian carrier fluid subject to simple shear. We consider high concentrations but well below the jamming limit. In this case, the suspension viscosity, η, depends on a number of physical parameters, such as the densities of the fluid ( ρ f) and solid phase ( ρ p), the dynamic fluid viscosity ( η 0), the sphere diameter ( d), the number density of particles per unit volume ( n), parameters related to frictional particle contacts (especially the Coulomb coefficient of sliding friction, μ c), the thermal energy k T (with k the Boltzmann constant and T the temperature), the imposed shear rate ( γ ˙), and the time passed since the moment this shear rate was imposed ( t) [3],
η = f ( ρ f , ρ p , η 0 , d , n , μ c , k T , γ ˙ , t ) .
(1)
Based on dimensional analysis, these physical parameters can be combined to form the following dimensionless groups:
η r = η η 0 , ϕ = n 1 6 π d 3 , ρ r = ρ p ρ f , P e γ ˙ = 3 π η 0 d 3 γ ˙ k T , R e γ ˙ = ρ f d 2 γ ˙ η 0 , t r = t k T η 0 d 3 , μ c ,
(2)
where ϕ is the particle volume fraction, P e γ ˙ is the Péclet number, and R e γ ˙ is the particle Reynolds number based on the applied shear rate. In this work, we will consider neutrally buoyant particles for which ρ r = 1, non-Brownian particles for which P e γ ˙ 1, simple shear in the Stokes regime at R e γ ˙ = 0.1, and steady-state conditions t r 1. Given these assumptions, dimensional analysis thus suggests that the relative viscosity depends only on ϕ and μ c,
η r = f ( ϕ , μ c ) .
(3)
In his seminal work in the early 1900s, Einstein [1,2] showed theoretically that for very dilute suspensions up to ϕ 0.05, when hydrodynamic and contact particle interactions can be ignored, the relative viscosity increases linearly with the particle volume fraction. It was many years later that Batchelor and Green [13] extended his work to moderate concentrations up to ϕ 0.1 by including hydrodynamic interactions, resulting in a second-order correction term in the expression for the relative viscosity. For simple shear flows, the result is given by
η r = 1 + 5 2 ϕ + 5.2 ϕ 2 .
(4)
For concentrations higher than ϕ 0.1, contact particle-particle interactions cannot be neglected anymore and become increasingly important with increasing concentration. For very high concentrations, the suspension viscosity diverges toward a state in which the suspension stops to flow, known as the jamming transition. The divergence of the viscosity is described by several empirical correlations [14–16] of which the Maron–Pierce correlation works particularly well for high volume fractions [15],
η r = ( 1 ϕ ϕ m ) 2 .
(5)
Here, ϕ m is the maximum flowable packing fraction or the volume fraction at which the suspension ceases to flow as the relative viscosity diverges. The maximum flowable packing fraction typically has a value of ϕ m = O ( 0.58 0.65 ), while the precise value is usually adjusted to provide a best fit with experimental data. Note that the Maron–Pierce correlation does not explicitly depend on μ c as one would expect based on Eq. (3) for the dense regime where particle contacts are frequent. However, this dependency is generally assumed to be implicitly contained in ϕ m, whose value is typically lower for higher values of μ c [10,11]. Furthermore, we remark that the Maron–Pierce correlation does not recover Einstein’s analytical expression for the dilute regime, and it suggests, though counter-intuitive, that there is still a dependency on ϕ m in this limit. On the other hand, the Krieger–Dougherty [16] relation for the suspension viscosity
η r = ( 1 ϕ ϕ m ) 2.5 ϕ m
(6)
does recover the Einstein viscosity in the dilute limit but appears less accurate in the very dense regime. Equations (5) and (6) are equivalent if and only if the maximum flowable packing fraction ϕ m = 0.8, although actual values are closer to 0.6.

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 ϕ m, which is dependent on the properties of the particles. The rheology described as a function of the normalized volume fraction ϕ / ϕ 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 ϕ = 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 μ c 0.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 = Σ 11 Σ 22 and N 2 = Σ 22 Σ 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 ϕ / ϕ 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.

An alternative description of the rheology, not hampered by the divergence of the viscosity near jamming, was introduced by Boyer 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 μ and the viscous number I v. Here, μ is defined as the ratio between the total shear stress and the wall-normal particle stress, μ = σ 12 / σ p , 22, and I v is defined as the ratio of the viscous shear stress over the wall-normal particle stress, I v = η 0 γ ˙ / σ 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:
μ ( I v ) = μ 1 + μ 2 μ 1 1 + I 0 / I v + I v + 5 2 ϕ m I v 1 / 2 ,
(7)
where μ 1 = 0.32, μ 2 = 0.7, and I 0 = 0.005 are empirically determined constants. In Eq. (7), the first two terms correspond to the contact contribution, whereas the last two terms correspond to the hydrodynamic contribution. Furthermore, they proposed a relation between the solid volume fraction and the viscous number,
ϕ ( I v ) = ϕ m 1 + I v 1 / 2 .
(8)
The frictional rheology is related to the previously discussed viscous rheology through η r = μ / I v, and from Eqs. (7) and (8), Einstein’s relation can be recovered in the limit of very large I v or very small ϕ. Note that this frictional description of the rheology does not explicitly account for interparticle friction. However, both μ and I v are strongly influenced by interparticle friction in the dense regime, and it is thus assumed to be implicitly accounted for.
Following the introduction of this frictional rheology framework, attempts have been made to improve the model description by adapting the contributions of the hydrodynamic and contact interactions [23,27,33]. Gallier 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
μ h = I v + [ η ] ϕ m I v 1 n and ϕ = ϕ m 1 + I v n ,
(9)
with [ η ] 2.4 and n = 0.4. Lecampion and Garagash [33] adapted both contributions to the macroscopic friction coefficient without significantly changing the total, leading to
μ ( ϕ ) = μ 1 + ϕ m β ( 1 ϕ ϕ m ) + ( I v + ( 5 2 ϕ m + 2 ) I v 1 / 2 ) ( 1 ϕ ϕ m ) 2 ,
(10)
with β = 0.158, μ 1 = 0.3, and ϕ m = 0.585. Chèvremont et al. [27] proposed changes to the hydrodynamic component of the model by Gallier et al. [23],
μ = μ 1 + μ 2 μ 1 1 + I 0 / I v + I v ( 1 ϕ m 1 1 + I v 1 2 ) [ η ] ,
(11)
with μ 1 = 0.36, μ 2 = 0.7, I 0 = 0.0133, and [ η ] = 5 / 2. Chèvremont 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 ϕ / ϕ 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 ( ϕ b = 0.1 0.6).

In our study, we consider three different values for the Coulomb coefficient of sliding friction, μ c = 0, 0.2, and 0.39, to investigate to what degree ϕ / ϕ 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 ϕ / ϕ m and I v provides insight into the relation between the two rheology descriptions.

The fluid phase is governed by the Navier–Stokes equations that are made dimensionless with the particle diameter as a length scale, l ref = d, the reciprocal of the shear rate as a time scale, t ref = 1 / γ ˙, and the fluid density, ρ f. These primary scales can be use to derive a velocity scale, u ref = γ ˙ d, and a stress scale, σ r e f = ρ f d 2 γ ˙ 2. This leads to the following nondimensional equations:
u = 0 ,
(12a)
u t + u u = p + 1 R e γ ˙ 2 u ,
(12b)
in which u is the velocity, and p is the pressure. The motions of the spherical particles are governed by the Newton–Euler equations
ρ r d u c d t = 6 π ( A p ( σ f n ) d A + F c ) ,
(13a)
ρ r d ω c d t = 60 π ( A p r × ( σ f n ) d A + T c ) ,
(13b)
in which u c and ω c are the translational and rotational velocity of a particle, respectively; the fluid stress tensor is defined as σ f = p I + 1 R e γ ˙ ( ( u ) + ( u ) T ); n is the outward pointing normal vector; r is the position vector with respect to the particle centroid; and the final terms, F c and T c, are the force and torque working on the particle as a result of contacts with other particles or the channel walls. The fluid and solid phase are coupled through a no-slip/no-penetration condition on the surface of the particles
u = u c + ω c × r on A p .
(14)
Contact between particles is handled using a soft-sphere contact model, including lubrication corrections for short-range particle-particle and particle-wall interactions, which is described in detail by Costa 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]
F i j , n = k n δ i j , n η n u i j , n ,
(15)
where δ i j , n is the overlap distance between two particles, u i j , n is the relative normal velocity of the particles, k n is the normal spring coefficient, and η n is the normal damping coefficient. The effect of solid friction is incorporated in the tangential contact force, which accounts for stick/slip transition according to [35]
F i j , t = min ( k t δ i j , t η t u i j , t , μ c F i j , n ) t i j .
(16)
Here, δ i j , t is the tangential displacement, u i j , t is the relative tangential velocity, k t is the tangential spring coefficient, η t is the tangential damping coefficient, μ c is the Coulomb coefficient of sliding friction, and t i j is the surface tangent unit vector.

The contact model coefficients k n, k t, η n, and η 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 Δ t with Δ 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.

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.001 d. The grid resolution used for the simulations is d / Δ 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, ρ r = 1, that are monodisperse in size. The main parameters that we vary are the bulk volume fraction, which ranges from ϕ b = 0.1 to ϕ b = 0.6, and the Coulomb coefficient of sliding friction, which is set to μ c = 0, μ c = 0.2, and μ c = 0.39. The value of μ 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.5 d and L y = 27 d. 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 γ ˙ = 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.

FIG. 1.

Schematic representation of the computational domain indicating the velocities of the top and bottom wall. The numbers along the coordinate axes denote the flow (1), velocity gradient (2), and vorticity (3) directions.

FIG. 1.

Schematic representation of the computational domain indicating the velocities of the top and bottom wall. The numbers along the coordinate axes denote the flow (1), velocity gradient (2), and vorticity (3) directions.

Close modal
TABLE I.

Parameters of the simulated cases: μc is the Coulomb coefficient of sliding friction, ϕb is the bulk solid volume fraction, Np is the number of particles, t t γ ˙ is the transient time of the simulation until the flow is fully developed, t s γ ˙ is the duration of the simulated steady-state, Lcore is the width of the core region, Lwall is the width of the wall region, ϕc is the average solid volume fraction in the core region, and γ ˙ c is the average shear rate in the core region.

μcϕbNp t t γ ˙ t s γ ˙Lcore/dLwall/dϕc γ ˙ c / γ ˙
0.1 940 94 102 2.5 5.5 0.104 0.980 
0.2 1880 94 102 2.5 5.5 0.208 0.980 
0.3 2819 63 53 2.5 5.5 0.309 0.969 
0.4 3759 94 61 2.5 5.5 0.410 0.901 
0.45 4229 109 105 2.5 5.5 0.454 0.876 
0.5 4699 125 89 5.75 0.486 0.831 
0.55 5169 156 83 1.6 5.95 0.517 0.974 
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 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ϕbNp t t γ ˙ t s γ ˙Lcore/dLwall/dϕc γ ˙ c / γ ˙
0.1 940 94 102 2.5 5.5 0.104 0.980 
0.2 1880 94 102 2.5 5.5 0.208 0.980 
0.3 2819 63 53 2.5 5.5 0.309 0.969 
0.4 3759 94 61 2.5 5.5 0.410 0.901 
0.45 4229 109 105 2.5 5.5 0.454 0.876 
0.5 4699 125 89 5.75 0.486 0.831 
0.55 5169 156 83 1.6 5.95 0.517 0.974 
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 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 
The data from the DNS consist of detailed three-dimensional velocity and concentration fields. These fields need to be averaged over space and time to gain useful insights into the flow dynamics. For plane Couette flow, the mean flow only varies in the wall-normal direction. Therefore, we apply plane-averaging over planes parallel to the walls, followed by averaging over time. For the so-called superficial average of the velocity, this can be expressed as
u f ¯ ( z ) = 1 L x L y x y u f ( x , y , z , t ) ( 1 α ( x , y , z , t ) ) Δ x Δ y ¯ ,
(17a)
u p ¯ ( z ) = 1 L x L y x y u p ( x , y , z , t ) α ( x , y , z , t ) Δ x Δ y ¯ ,
(17b)
for the fluid and solid phase, respectively. The brackets and bar denote the spatial averaging and ensemble averaging of samples collected over time, respectively. Here, α is the solid phase indicator function (0 in the fluid and 1 in the particle phase), which is also used to determine the mean local concentration as a function of the wall-normal coordinate using
ϕ ¯ ( z ) = 1 L x L y x y α ( x , y , z , t ) Δ x Δ y ¯ .
(18)
Using the superficial averages and the mean local concentration, the intrinsic (phase-averaged) velocities can be computed as
u f f ¯ = u f ¯ 1 ϕ ¯ ,
(19a)
u p p ¯ = u p ¯ ϕ ¯ .
(19b)
In addition to the averaging, the superficial (macroscopic) particle stresses need to be computed to describe the rheology of the suspensions. The microscopic particle stress distribution within the particles is not explicitly solved in the DNS as this is not required for rigid particles. However, it can be approximated from the average microscopic stress within a particle, based on the stress distribution on the particle’s surface and the contact forces acting on the particle according to
σ p , i j 6 π ( A p r j σ f , i l n l d A + r j F c , i ) ,
(20)
where the surface integral is equal to the hydrodynamic particle stresslet for freely rotating particles with zero net torque and the summation accounts for the contribution from the particle contacts. Similar to the superficial velocity averages, the superficial particle stress can be determined as a function of the wall-normal coordinate using
σ p , i j ¯ ( z ) = 1 L x L y x , y σ p , i j ( x , y , z , t ) α ( x , y , z , t ) Δ x Δ y ¯ .
(21)
Note that σ p , i j at the right-hand side of Eq. (21) is approximated by the 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.

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 ϕ 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 μ c = 0 (a) and μ 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 μ 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.

FIG. 2.

Visualization of the particles in a 3D-slab and two perpendicular cross-sections of the flow domain for the cases with ϕ b = 0.5 and a friction coefficient of μ c = 0 (a) and μ c = 0.39 (b). The particles are colored by their centroid velocity in the flow direction normalized by the wall velocity.

FIG. 2.

Visualization of the particles in a 3D-slab and two perpendicular cross-sections of the flow domain for the cases with ϕ b = 0.5 and a friction coefficient of μ c = 0 (a) and μ c = 0.39 (b). The particles are colored by their centroid velocity in the flow direction normalized by the wall velocity.

Close modal

Figure 3 shows a snapshot of the particles in the flow for the case with ϕ b = 0.6 and μ 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 ϕ 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.

FIG. 3.

Visualization of the particles in two perpendicular cross-sections of the flow domain for the case with ϕ b = 0.6 and a friction coefficient of μ c = 0.39. The particles are colored by their centroid velocity in the flow direction normalized by the wall velocity.

FIG. 3.

Visualization of the particles in two perpendicular cross-sections of the flow domain for the case with ϕ b = 0.6 and a friction coefficient of μ c = 0.39. The particles are colored by their centroid velocity in the flow direction normalized by the wall velocity.

Close modal

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 ϕ b > 0.4. For the cases with a concentration of ϕ b 0.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 ϕ b = 0.55 to ϕ b = 0.6, for μ 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 ϕ b = 0.55 and ϕ 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 ϕ b = 0.5 and ϕ 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].

FIG. 4.

Mean solid volume fraction for cases with μ c = 0 (a), μ c = 0.2 (b), and μ c = 0.39 (c) as a function of the wall-normal coordinate.

FIG. 4.

Mean solid volume fraction for cases with μ c = 0 (a), μ c = 0.2 (b), and μ c = 0.39 (c) as a function of the wall-normal coordinate.

Close modal

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 μ 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.

FIG. 5.

Mean intrinsic fluid velocity normalized by the wall velocity for cases with μ c = 0 (a), μ c = 0.2 (b), and μ c = 0.39 (c) as a function of the wall-normal coordinate.

FIG. 5.

Mean intrinsic fluid velocity normalized by the wall velocity for cases with μ c = 0 (a), μ c = 0.2 (b), and μ c = 0.39 (c) as a function of the wall-normal coordinate.

Close modal

In Fig. 6, we compare the mean intrinsic fluid and particle velocity profiles for the cases with ϕ b = 0.4 and for friction coefficients of μ c = 0 and μ 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 1.25 d 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.

FIG. 6.

Mean intrinsic velocity normalized by the wall velocity for the fluid and solid phase as a function of the wall-normal coordinate, shown for the cases with ϕ b = 0.4 and for friction coefficients of μ c = 0 and μ c = 0.39. The black dashed line indicates the single phase reference, and the inset shows the slip between phases close to the wall.

FIG. 6.

Mean intrinsic velocity normalized by the wall velocity for the fluid and solid phase as a function of the wall-normal coordinate, shown for the cases with ϕ b = 0.4 and for friction coefficients of μ c = 0 and μ c = 0.39. The black dashed line indicates the single phase reference, and the inset shows the slip between phases close to the wall.

Close modal

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 ϕ b = 0.45, while set somewhat smaller for the cases with ϕ b 0.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 ϕ b 0.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 ( γ ˙ c) and average concentration ( ϕ c), of which the values are also listed in Table I. The cases with ϕ b = 0.4 will be used as an example throughout this work because these cases nicely illustrate the suspension in the dense regime.

In order to investigate the suspension microstructure in the core region, we computed the particle pair distribution function (PPDF). This expresses the likelihood of finding a neighboring particle at a certain distance and orientation, r, with respect to an arbitrary reference particle in the core region
g ( r ) = V / V v o x N ( N 1 ) p p d δ ( r ( x p d x p ) ) ,
(22)
where V is the total volume of the core region, V v o x is the voxel volume of the regular voxel grid used for the calculation, N is the total number of particles in the core region, x p is the centroid position of the reference particle, and x p d is the centroid position of the neighboring particle. The function δ gives a value of 1 if the particle centroid of the neighboring particle is in the voxel of interest and 0 otherwise. In Eq. (22), g ( r ) is normalized such that g ( r ) = 1 for a statistically homogeneous particle distribution in space. Figure 7 shows the PPDFs for a selection of cases with a bulk volume fraction up to ϕ b = 0.5 and for friction coefficients of μ c = 0 and μ c = 0.39, calculated for the core region of the domain. The color shows the probability of finding another particle centroid at that location, with blue being a lower than average probability and red being higher. Note that the scale of the colorbar is split in two parts at a value of 1, with both red and blue scaling linearly but with different factors. The top row shows the cases with μ c = 0, while the bottom row shows the cases with μ c = 0.39. The bulk volume fraction increases from left to right. The PPDFs obtained from our simulations show good qualitative agreement with the experiments of Parsi and Gadala-Maria [50] and Blanc 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 ° from the positive y-axis. Furthermore, the spots of high probability at r 1 d 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 2 d 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 ϕ 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.
FIG. 7.

PPDF shown for the core region and shear plane. The different panels correspond to different concentrations ( ϕ b = 0.2 0.5, increasing from left to right) and different friction coefficients ( μ c = 0 for top and μ c = 0.39 for bottom panels). The black disc denotes the reference particle. The shell of 0.5 < r / d < 1 is blanked out as it is an excluded volume inaccessible to neighboring particles. Note that the scale of the colorbar is linear, but split in two parts, with an eleven times steeper increase for values above 1 than for values below 1.

FIG. 7.

PPDF shown for the core region and shear plane. The different panels correspond to different concentrations ( ϕ b = 0.2 0.5, increasing from left to right) and different friction coefficients ( μ c = 0 for top and μ c = 0.39 for bottom panels). The black disc denotes the reference particle. The shell of 0.5 < r / d < 1 is blanked out as it is an excluded volume inaccessible to neighboring particles. Note that the scale of the colorbar is linear, but split in two parts, with an eleven times steeper increase for values above 1 than for values below 1.

Close modal

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 θ 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 π / 4 and 7 π / 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 θ = π 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 ϕ b = 0.4 over a range of 0 < θ < π. 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.

FIG. 8.

PPDFs in the core region, averaged over a thin shell around the reference particle ( 1 < r / d < 1.05) and shown for a selection of cases with bulk volume fractions up to ϕ b = 0.5 (a) and comparison of the cases with ϕ b = 0.4 with the literature (b).

FIG. 8.

PPDFs in the core region, averaged over a thin shell around the reference particle ( 1 < r / d < 1.05) and shown for a selection of cases with bulk volume fractions up to ϕ b = 0.5 (a) and comparison of the cases with ϕ b = 0.4 with the literature (b).

Close modal
The applied shear results in stresses inside the microstructure of the suspensions, which can be used to determine the rheology of the suspension. The total shear stress is defined as [cf. Eq. (4.1) in Batchelor [53] and Eq. (7.3) in Guazzelli and Morris [54]]
σ 12 ¯ = 1 R e γ ˙ d v m ¯ d z + σ p , 12 ¯ v f w f ¯ ρ r v p w p ¯ ,
(23)
where the first term on the right-hand side is the viscous shear stress with v m ¯ = v p ¯ + v f ¯ the mixture velocity, the second term on the right hand side is the particle shear stress, and the last two terms on the right hand side correspond to the Reynolds stresses from fluctuating fluid and particle motions, respectively. The Reynolds stresses are very close to zero for the present simulations, as expected in the Stokes regime, and are hence omitted in our analysis of the stresses. Figure 9 shows the shear stresses for the cases with ϕ b = 0.4. Comparing the total stress of the three friction coefficients shows that the shear stress increases with increasing friction coefficient. For all three cases, the total stress is constant in the core region, as expected for a fully developed flow, but fluctuations are observed in the wall regions. The fluctuations are caused by the approximation of the microscopic particle stress within the particles by the particle-averaged value, resulting in a bias when the stress distribution is strongly inhomogeneous as is the case in the wall regions, see the discussion in Sec. II C. Irrespective of the friction coefficient, the largest part of the shear stress is the particle stress. This particle stress also accounts for the increase in total stress with increasing friction coefficient, as the viscous stress is similar for all three cases. The particle stress consists of two separate contributions originating from the hydrodynamic interactions σ p , 12 h y d r o ¯ and the contacts σ p , 12 c o n t a c t ¯ between particles. These are shown in Fig. 9 as the dashed and dotted line, respectively. The hydrodynamic contribution is larger than the contact contribution for the cases with μ c = 0 and μ c = 0.2, whereas the two contributions are equal for the case with μ c = 0.39. By increasing the friction coefficient, only the contact contribution to the particle stress increases. This shows that the increase in total shear stress is completely caused by the stress originating from the particle contacts.
FIG. 9.

Mean shear stresses σ i , 12 ¯ for the cases with ϕ b = 0.4 for a Coulomb friction coefficient of μ c = 0 (a), μ c = 0.2 (b), and μ c = 0.39 (c), plotted as a function of the wall-normal coordinate. Dashed lines show the decomposition of the particle stress in a hydrodynamic and contact contribution.

FIG. 9.

Mean shear stresses σ i , 12 ¯ for the cases with ϕ b = 0.4 for a Coulomb friction coefficient of μ c = 0 (a), μ c = 0.2 (b), and μ c = 0.39 (c), plotted as a function of the wall-normal coordinate. Dashed lines show the decomposition of the particle stress in a hydrodynamic and contact contribution.

Close modal
Next to the shear stress components, the normal stress components are important for the description of the suspension rheology. Figure 10 shows the wall-normal stresses for the cases with ϕ b = 0.4. The total wall-normal stress is defined as
σ 22 ¯ = p f ¯ + σ p , 22 ¯ ,
(24)
where the first term is the fluid pressure and the second term is the wall-normal particle stress. We remark that in the DNS a zero reference pressure at the walls is used. Similar to the shear stress, the total stress magnitude increases with increasing friction coefficient. The particle stress provides the largest contribution to the wall-normal stress and moreover is the term that increases with increasing friction coefficient. Decomposing the particle stress into the hydrodynamic contribution σ p , 22 h y d r o ¯ and the contact contribution σ p , 22 c o n t a c t ¯, shown as the dashed and dotted line, respectively, shows that the contacts give the largest contribution to the particle stress. Moreover, the hydrodynamic contribution has the opposite sign to the total particle stress and is close to zero.
FIG. 10.

Mean wall-normal stresses σ i , 22 ¯ for the cases with ϕ b = 0.4 for a Coulomb friction coefficient of μ c = 0 (a), μ c = 0.2 (b), and μ c = 0.39 (c), plotted as a function of the wall-normal coordinate. Dashed lines show the decomposition of the particle stress in a hydrodynamic and a contact contribution.

FIG. 10.

Mean wall-normal stresses σ i , 22 ¯ for the cases with ϕ b = 0.4 for a Coulomb friction coefficient of μ c = 0 (a), μ c = 0.2 (b), and μ c = 0.39 (c), plotted as a function of the wall-normal coordinate. Dashed lines show the decomposition of the particle stress in a hydrodynamic and a contact contribution.

Close modal
From the stress budget in the shear plane we can determine the relative viscosity using the total shear stress divided by the viscous shear stress
η η 0 = 1 + σ p , 12 ¯ 1 R e γ ˙ v m ¯ z .
(25)
Doing this, for all the cases, we can determine the relative viscosity as a function of the volume fraction, which is shown in Fig. 11(a). Our data are shown together with the numerical results of Yeo and Maxey [22] and Gallier 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 ϕ 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 ϕ 0.1 and defined in Eq. (4), is shown for reference as the dashed-dotted line. We observe that our results at ϕ 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 ϕ b 0.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 ϕ 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 ϕ m = 0.70, ϕ m = 0.66, and ϕ m = 0.64 for the cases with μ c = 0, μ c = 0.2, and μ 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 ϕ m do not change significantly when the cases with ϕ b = 0.5 are excluded from the fit, indicating that the core region of these cases is not affected by wall-induced layering.
FIG. 11.

Relative viscosity as a function of the volume fraction (a) and the volume fraction normalized by the maximum flowable packing fraction (b), where ϕ m is determined by fitting the Maron–Pierce correlation [15], Eq. (5), to the data indicated by the solid line. Colored markers show our data, where the open markers indicate the cases with ϕ b 0.55, exhibiting pronounced particle layering, that were excluded from the fit. Black markers show data from the literature [17,18,22,23]. The dashed lines show the relative viscosity derived from the frictional rheology model defined in Eq. (30) using η r = μ / I v.

FIG. 11.

Relative viscosity as a function of the volume fraction (a) and the volume fraction normalized by the maximum flowable packing fraction (b), where ϕ m is determined by fitting the Maron–Pierce correlation [15], Eq. (5), to the data indicated by the solid line. Colored markers show our data, where the open markers indicate the cases with ϕ b 0.55, exhibiting pronounced particle layering, that were excluded from the fit. Black markers show data from the literature [17,18,22,23]. The dashed lines show the relative viscosity derived from the frictional rheology model defined in Eq. (30) using η r = μ / I v.

Close modal

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 ϕ c 0.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 ϕ b 0.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 ϕ c rather than ϕ c / ϕ 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 ϕ c / ϕ m 0.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 ϕ c / ϕ m, the separate contributions do not completely collapse.

FIG. 12.

Decomposition of the relative viscosity into the hydrodynamic and contact contribution as a function of the volume fraction (a) and the volume fraction normalized by the maximum flowable packing fraction (b). The markers show the DNS results, where the open markers indicate the cases with ϕ b 0.55, exhibiting pronounced particle layering. The lines show the relative viscosity derived from the frictional rheology model defined in Eq. (30) using η r = μ / I v.

FIG. 12.

Decomposition of the relative viscosity into the hydrodynamic and contact contribution as a function of the volume fraction (a) and the volume fraction normalized by the maximum flowable packing fraction (b). The markers show the DNS results, where the open markers indicate the cases with ϕ b 0.55, exhibiting pronounced particle layering. The lines show the relative viscosity derived from the frictional rheology model defined in Eq. (30) using η r = μ / I v.

Close modal
As discussed in the Introduction section, an alternative way to describe the shear rheology of suspensions is using the frictional rheology description, where the macroscopic friction coefficient μ is expressed as a function of the viscous number I v. These can be determined from our simulations using the viscous shear stress over the wall-normal particle stress
I v = 1 R e γ ˙ v m ¯ z σ p , 22 ¯
(26)
for the viscous number and the total shear stress over the wall-normal particle stress
μ = 1 R e γ ˙ v m ¯ z + σ p , 12 ¯ σ p , 22 ¯
(27)
for the macroscopic friction coefficient. Figure 13(a) shows the macroscopic friction coefficient as a function of the viscous number. Our data are denoted by the colored markers, the colored lines present a model description of our data, and the black lines are models from the literature [19,23,27,33]. Comparing our data to the models from the literature shows good agreement at low viscous number (high solid volume fraction); however, at high I v, there is only a qualitative agreement between our data and the models. At these high viscous numbers, our simulations give a consistently higher macroscopic friction than the models predict, for which the model of Gallier 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 μ 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.
FIG. 13.

Frictional rheology with (a) the macroscopic friction coefficient, (b) the normalized volume fraction, (c) the hydrodynamic part of the macroscopic friction coefficient, and (d) the contact part of the macroscopic friction coefficient as a function of the viscous number. Colored lines show fits to the data, and black lines show models from the literature [19,23,27,33]. Open markers indicate the cases exhibiting pronounced particle layering in the core and were excluded from fits.

FIG. 13.

Frictional rheology with (a) the macroscopic friction coefficient, (b) the normalized volume fraction, (c) the hydrodynamic part of the macroscopic friction coefficient, and (d) the contact part of the macroscopic friction coefficient as a function of the viscous number. Colored lines show fits to the data, and black lines show models from the literature [19,23,27,33]. Open markers indicate the cases exhibiting pronounced particle layering in the core and were excluded from fits.

Close modal
A slight deviation from literature models is also present for the relation between the solid volume fraction and the viscous number, as shown in Fig. 13(b). Here, the model of Gallier 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
ϕ = ϕ m 1 + I v n .
(28)
Here, ϕ m and n are the fitting parameters. To see how well the viscous number collapses the data, we fit Eq. (28) to each friction coefficient separately. The fits are shown in Fig. 13(b) as the colored solid lines, and the fitting coefficients are listed in Table II. The results for all three friction coefficients can be modeled very well separately using this relation, giving values for ϕ m (see Table II) that are slightly higher than those determined using the Maron–Pierce correlation [15] for the relative viscosity [see legend of Fig. 11(b)]. However, the similarity in fitting coefficient n indicates that our data might also be modeled using a single value for n.
TABLE II.

Fitting coefficients for the frictional rheology model compared to the coefficients found in the literature.

μc = 0μc = 0.2μc = 0.39Boyer 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 
I0 0.0194 0.0368 0.0584 0.005 0.005 0.0133 
μc = 0μc = 0.2μc = 0.39Boyer 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 
I0 0.0194 0.0368 0.0584 0.005 0.005 0.0133 
To model the relation between the macroscopic friction coefficient and the viscous number, the macroscopic friction coefficient is decomposed into its separate contributions. The hydrodynamic part μ h is shown in Fig. 13(c), together with the corresponding contribution of the literature models. Here, we can see that the model of Gallier 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 μ h. Here, ϕ m and n are already determined, leaving [ η ] as the only fitting coefficient. This fit gives the values of [ η ] 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 μ 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 μ 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 μ 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
μ c = μ 1 + μ 2 μ 1 1 + I 0 / I v ,
(29)
we need to choose appropriate values for μ 1 as we have no data at very low viscous numbers. To do this, we used the result of Chèvremont et al. [27] by looking at the low I v plateau in the contact contribution. The values we chose are μ 1 = 0.2, μ 1 = 0.3, and μ 1 = 0.36 for the cases with μ c = 0, μ c = 0.2, and μ c = 0.39, respectively. The resulting fits are shown in Fig. 13(d), and the fitting coefficients μ 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 μ 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.
Using Eqs. (9), (28), and (29), we can derive a relation for the relative viscosity based on the frictional rheology model
η r = μ c + μ h I v = μ 1 ( ϕ m / ϕ 1 ) 1 / n + μ 2 μ 1 ( ϕ m / ϕ 1 ) 1 / n + I 0 + 1 + [ η ] ϕ ( 1 ϕ ϕ m ) 1 .
(30)
This results in a relative viscosity that is not only a function of ϕ / ϕ m but also depends on the model coefficients that are, in turn, a function of the solid friction. Equation (30) does, however, provide a very good description of the relative viscosity as a function of the solid volume fraction, as can be seen from the dashed lines in Fig. 11(a). At low volume fractions, Eq. (30) nearly coincides with the Batchelor and Green relation [13], shown by the dashed-dotted line, although in the limit of ϕ 0 our model does overpredict the Einstein viscosity since [ η ] > 2.5. Normalizing the volume fraction by ϕ m found by the Maron–Pierce fit to the data also shows a collapse for Eq. (30); however, this does not result in a single curve for different solid friction coefficients, as is shown by the dashed lines in Fig. 11(b), although the difference is small.

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 η r = μ / I v and ϕ / ϕ m is a monotonically decreasing function in I v, which is only weakly affected by μ 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.

Having concluded that the frictional rheology and viscous rheology are equivalent with regard to the relative shear viscosity, we turn our attention to the normal stresses. Looking at all three normal stress components, we observe non-Newtonian phenomena that occur in suspensions, which are the anisotropy in the normal stresses and the related anisotropy in the microstructure. Figure 14 shows the mean normal particle stresses in the core region of the channel normalized by the mean total shear stress and plotted as a function of the normalized solid volume fraction. The chosen normalizations are consistent with the viscous rheology framework, noting that the total shear stress scales with the viscous shear stress times a function of ϕ c / ϕ m as found above from the shear rheology analysis. Looking at the three different components of the normal stress, σ p , 11 ¯, σ p , 22 ¯, and σ p , 33 ¯ in Figs. 14(a)14(c), respectively, it becomes clear that the normal particle stresses are anisotropic. The effect of friction does not seem to be accounted for in the normalization of the solid volume fraction, ϕ c / ϕ m, as the data of the three friction coefficients does not collapse to a single curve, except for the wall-normal particle stress σ p , 22 ¯. Plotting the normal stresses as a function of the viscous number would also not result in a collapse of the data with different solid friction coefficients, due to the monotonic relation between ϕ c / ϕ m and I v that is only weakly dependent on μ c as observed from Fig. 13(b). We remark that the normalization of the normal stresses with the total shear stress is also consistent with the frictional rheology framework, as the total shear stress scales with the wall-normal particle stress times a function of I v as earlier found from the shear rheology analysis. The effect of layering in the domain is best visible in Fig. 14(a) for the streamwise particle stress σ p , 11 ¯, where the cases with ϕ b 0.55, indicated by the open markers, show a clear deviation from the trend. The lines in the figure represent polynomial fits to the data based on the work of Badia 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
σ p , 11 ¯ σ 12 ¯ = a 1 ( ϕ c ϕ m ) a 2 ,
(31a)
σ p , 22 ¯ σ 12 ¯ = σ p , 11 ¯ σ 12 ¯ ( b 1 + b 2 ϕ c ϕ m + b 3 ( ϕ c ϕ m ) 2 ) ,
(31b)
σ p , 33 ¯ σ 12 ¯ = σ p , 11 ¯ σ 12 ¯ ( c 1 + c 2 ϕ c ϕ m + c 3 ( ϕ c ϕ m ) 4 ) .
(31c)
The fitting parameters used to describe our data are shown in Table III for each of the friction coefficients together with the parameters used by Badia 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 ϕ b = 0.4, the particle stress is almost completely governed by the contact contribution. Figure 14(d) shows the normal stress ratios λ 2 = σ p , 22 ¯ / σ p , 11 ¯ and λ 3 = σ p , 33 ¯ / σ p , 11 ¯ 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, λ 2 shows an increase with increasing volume fraction for all three friction coefficients, whereas λ 3 is more or less constant around a value of 0.5 with a slight increase at high volume fractions.
FIG. 14.

Mean normal stresses σ p , i i ¯ normalized by the mean total shear stress and plotted as a function of the normalized volume fraction. Lines show fits of Eqs. (31a)(31c) to the data. Open markers indicate the cases with where the open markers indicate the cases with ϕ b 0.55, exhibiting strong particle layering in the core and were excluded from the fits.

FIG. 14.

Mean normal stresses σ p , i i ¯ normalized by the mean total shear stress and plotted as a function of the normalized volume fraction. Lines show fits of Eqs. (31a)(31c) to the data. Open markers indicate the cases with where the open markers indicate the cases with ϕ b 0.55, exhibiting strong particle layering in the core and were excluded from the fits.

Close modal
TABLE III.

Fitted coefficients of the correlations for the mean normal stresses given by Eqs. (31a)(31c). The columns show the results for μc = 0, μc = 0.2, and μc = 0.39 and reference values from Badia et al. [31].

μc = 0μc = 0.2μc = 0.39Badia et al. [31]
a1 −2.3022 −2.3305 −2.3648 −2.4247 
a2 2.7915 2.9724 3.0688 4.1280 
b1 0.9717 0.5938 0.3792 0.3750 
b2 −2.0377 −0.2760 0.4415 0.0366 
b3 2.5687 0.9074 0.3524 0.4846 
c1 0.5895 0.4977 0.3361 2.1446 
c2 −0.5932 −0.3755 −0.0983 −2.7234 
c3 1.8373 1.2661 1.0189 1.5759 
μc = 0μc = 0.2μc = 0.39Badia et al. [31]
a1 −2.3022 −2.3305 −2.3648 −2.4247 
a2 2.7915 2.9724 3.0688 4.1280 
b1 0.9717 0.5938 0.3792 0.3750 
b2 −2.0377 −0.2760 0.4415 0.0366 
b3 2.5687 0.9074 0.3524 0.4846 
c1 0.5895 0.4977 0.3361 2.1446 
c2 −0.5932 −0.3755 −0.0983 −2.7234 
c3 1.8373 1.2661 1.0189 1.5759 
A common way to describe the anisotropy in the normal particle stresses is using the mean normal stress difference
N 1 = σ p , 11 ¯ σ p , 22 ¯ ,
(32a)
N 2 = σ p , 22 ¯ σ p , 33 ¯ .
(32b)
Figure 15 shows the mean normal stress differences, normalized by the mean total shear stress and plotted as a function of the normalized solid volume fraction, together with data and models from the literature. The dashed lines show the mean normal stress differences computed from the fits of Eqs. (31a)(31c) to our data. The first normal stress difference, shown in Fig. 15(a), is slightly negative for volume fractions up to ϕ c / ϕ m 0.8 after which there is a strong increase to high positive values. For low volume fractions, our data show good agreement with numerical results from the literature [23,26] and qualitative agreement to the polynomial description of Badia 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 σ p , 11 ¯. 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 μ 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.
FIG. 15.

Mean normal stress differences N 1 (a) and N 2 (b), both normalized by the mean total shear stress and plotted as a function of the normalized solid volume fraction ϕ c / ϕ m. Our data are shown by the colored markers, and the dashed lines are the normal stress differences calculated from the fits of Eq. (31). Black markers indicate data from the literature [23,26,31], and the solid line gives a polynomial fit used by Badia et al. [31].

FIG. 15.

Mean normal stress differences N 1 (a) and N 2 (b), both normalized by the mean total shear stress and plotted as a function of the normalized solid volume fraction ϕ c / ϕ m. Our data are shown by the colored markers, and the dashed lines are the normal stress differences calculated from the fits of Eq. (31). Black markers indicate data from the literature [23,26,31], and the solid line gives a polynomial fit used by Badia et al. [31].

Close modal

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 ϕ / ϕ m and I v depends only weakly on the friction coefficient. Considering the normalized wall-normal particle stress component σ p , 22 ¯ / σ 12 ¯, we do observe from Fig. 14 a well-nigh collapse as a function of normalized volume fraction, ϕ c / ϕ m, for the three friction coefficients. In passing, we remark that this is consistent with the collapse for μ as a function of I v in the frictional rheology framework, as μ = σ 12 ¯ / σ p , 22 ¯ 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 σ p , 33 ¯ 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.

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 ϕ b 0.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 ϕ c / ϕ 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 ϕ c / ϕ m 0.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 ϕ c / ϕ 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.

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).

The authors have no conflicts to disclose.

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].

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, ϕ 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.

FIG. 16.

Core volume fraction (a) and macroscopic friction coefficient (b) as a function of the viscous number, based on different three choices for the normalization of the wall-normal stress and shown for cases with μ c = 0 and μ c = 0.39. The different colors denote normalization with the full wall-normal stress, the full wall-normal particle stress, and the contact contribution to the wall-normal particle stress, respectively. Dots and diamonds denote μ c = 0 and μ c = 0.39, respectively. Only data for cases with ϕ b = 0.1 0.5 is shown due to the pronounced wall-induced layering observed in the cases with ϕ b 0.55.

FIG. 16.

Core volume fraction (a) and macroscopic friction coefficient (b) as a function of the viscous number, based on different three choices for the normalization of the wall-normal stress and shown for cases with μ c = 0 and μ c = 0.39. The different colors denote normalization with the full wall-normal stress, the full wall-normal particle stress, and the contact contribution to the wall-normal particle stress, respectively. Dots and diamonds denote μ c = 0 and μ c = 0.39, respectively. Only data for cases with ϕ b = 0.1 0.5 is shown due to the pronounced wall-induced layering observed in the cases with ϕ b 0.55.

Close modal

Figure 16(b) shows the macroscopic friction coefficient as a function of the viscous number, where both μ 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 μ and I v.

1.
Einstein
,
A.
, “Eine neue bestimmung der moleküldimensionen,” Ph.D. thesis (Universität Zürich, Zürich, 1905).
2.
Einstein
,
A.
, “
Berichtigung zu meiner arbeit: Eine neue bestimmung der moleküldimensionen
,”
Ann. Phys.
339
(
3
),
591
592
(
1911
).
3.
Stickel
,
J. J.
, and
R. L.
Powell
, “
Fluid mechanics and rheology of dense suspensions
,”
Annu. Rev. Fluid Mech.
37
,
129
149
(
2005
).
4.
Guazzelli
,
É.
, and
O.
Pouliquen
, “
Rheology of dense granular suspensions
,”
J. Fluid Mech.
852
,
P1
(
2018
).
5.
Ness
,
C.
,
R.
Seto
, and
R.
Mari
, “
The physics of dense suspensions
,”
Annu. Rev. Condens. Matter Phys.
13
,
97
117
(
2022
).
6.
Singh
,
A.
,
S.
Pednekar
,
J.
Chun
,
M. M.
Denn
, and
J. F.
Morris
, “
From yielding to shear jamming in a cohesive frictional suspension
,”
Phys. Rev. Lett.
122
(
9
),
098004
(
2019
).
7.
Guy
,
B. M.
,
M.
Hermes
, and
W. C. K.
Poon
, “
Towards a unified description of the rheology of hard-particle suspensions
,”
Phys. Rev. Lett.
115
(
8
),
088304
(
2015
).
8.
Lobry
,
L.
,
E.
Lemaire
,
F.
Blanc
,
S.
Gallier
, and
F.
Peters
, “
Shear thinning in non-Brownian suspensions explained by variable friction between particles
,”
J. Fluid Mech.
860
,
682
710
(
2019
).
9.
Seto
,
R.
,
R.
Mari
,
J. F.
Morris
, and
M. M.
Denn
, “
Discontinuous shear thickening of frictional hard-sphere suspensions
,”
Phys. Rev. Lett.
111
(
21
),
218301
(
2013
).
10.
Mari
,
R.
,
R.
Seto
,
J. F.
Morris
, and
M. M.
Denn
, “
Shear thickening, frictionless and frictional rheologies in non-Brownian suspensions
,”
J. Rheol.
58
(
6
),
1693
1724
(
2014
).
11.
Singh
,
A.
,
R.
Mari
,
M. M.
Denn
, and
J. F.
Morris
, “
A constitutive model for simple shear of dense frictional suspensions
,”
J. Rheol.
62
(
2
),
457
468
(
2018
).
12.
Morris
,
J. F.
, “
Shear thickening of concentrated suspensions: Recent developments and relation to other phenomena
,”
Annu. Rev. Fluid Mech.
52
,
121
144
(
2020
).
13.
Batchelor
,
G. K.
, and
J. T.
Green
, “
The determination of the bulk stress in a suspension of spherical particles to order c 2
,”
J. Fluid Mech.
56
(
3
),
401
427
(
1972
).
14.
Eilers
,
H.
, “
Die viskosität von emulsionen hochviskoser stoffe als funktion der konzentration
,”
Kolloid Z.
97
(
3
),
313
321
(
1941
).
15.
Maron
,
S. H.
, and
P. E.
Pierce
, “
Application of ree-eyring generalized flow theory to suspensions of spherical particles
,”
J. Colloid Sci.
11
(
1
),
80
95
(
1956
).
16.
Krieger
,
I. M.
, and
T. J.
Dougherty
, “
A mechanism for non-newtonian flow in suspensions of rigid spheres
,”
Trans. Soc. Rheol.
3
(
1
),
137
152
(
1959
).
17.
Zarraga
,
I. E.
,
D. A.
Hill
, and
D. T.
Leighton
, Jr.
, “
The characterization of the total stress of concentrated suspensions of noncolloidal spheres in newtonian fluids
,”
J. Rheol.
44
(
2
),
185
220
(
2000
).
18.
Dbouk
,
T.
,
L.
Lobry
, and
E.
Lemaire
, “
Normal stresses in concentrated non-Brownian suspensions
,”
J. Fluid Mech.
715
,
239
272
(
2013
).
19.
Boyer
,
F.
,
É.
Guazzelli
, and
O.
Pouliquen
, “
Unifying suspension and granular rheology
,”
Phys. Rev. Lett.
107
(
18
),
188301
(
2011
).
20.
Dagois-Bohy
,
S.
,
S.
Hormozi
,
É.
Guazzelli
, and
O.
Pouliquen
, “
Rheology of dense suspensions of non-colloidal spheres in yield-stress fluids
,”
J. Fluid Mech.
776
,
R2
(
2015
).
21.
Sierou
,
A.
, and
J. F.
Brady
, “
Rheology and microstructure in concentrated noncolloidal suspensions
,”
J. Rheol.
46
(
5
),
1031
1056
(
2002
).
22.
Yeo
,
K.
, and
M. R.
Maxey
, “
Dynamics of concentrated suspensions of non-colloidal particles in Couette flow
,”
J. Fluid Mech.
649
,
205
231
(
2010
).
23.
Gallier
,
S.
,
E.
Lemaire
,
F.
Peters
, and
L.
Lobry
, “
Rheology of sheared suspensions of rough frictional particles
,”
J. Fluid Mech.
757
,
514
549
(
2014
).
24.
Gallier
,
S.
,
E.
Lemaire
,
L.
Lobry
, and
F.
Peters
, “
Effect of confinement in wall-bounded non-colloidal suspensions
,”
J. Fluid Mech.
799
,
100
127
(
2016
).
25.
Gallier
,
S.
,
F.
Peters
, and
L.
Lobry
, “
Simulations of sheared dense noncolloidal suspensions: Evaluation of the role of long-range hydrodynamics
,”
Phys. Rev. Fluids
3
(
4
),
042301
(
2018
).
26.
Seto
,
R.
, and
G. G.
Giusteri
, “
Normal stress differences in dense suspensions
,”
J. Fluid Mech.
857
,
200
215
(
2018
).
27.
Chèvremont
,
W.
,
B.
Chareyre
, and
H.
Bodiguel
, “
Quantitative study of the rheology of frictional suspensions: Influence of friction coefficient in a large range of viscous numbers
,”
Phys. Rev. Fluids
4
(
6
),
064302
(
2019
).
28.
Singh
,
A.
, and
P. R.
Nott
, “
Experimental measurements of the normal stresses in sheared Stokesian suspensions
,”
J. Fluid Mech.
490
,
293
320
(
2003
).
29.
Boyer
,
F.
,
O.
Pouliquen
, and
É.
Guazzelli
, “
Dense suspensions in rotating-rod flows: Normal stresses and particle migration
,”
J. Fluid Mech.
686
,
5
25
(
2011
).
30.
Singh
,
A.
, and
P. R.
Nott
, “
Normal stresses and microstructure in bounded sheared suspensions via Stokesian dynamics simulations
,”
J. Fluid Mech.
412
,
279
301
(
2000
).
31.
Badia
,
A.
,
Y.
D’Angelo
,
F.
Peters
, and
L.
Lobry
, “
Frame-invariant modeling for non-Brownian suspension flows
,”
J. Non-Newtonian Fluid Mech.
309
,
104904
(
2022
).
32.
Peters
,
F.
,
G.
Ghigliotti
,
S.
Gallier
,
F.
Blanc
,
E.
Lemaire
, and
L.
Lobry
, “
Rheology of non-Brownian suspensions of rough frictional particles under shear reversal: A numerical study
,”
J. Rheol.
60
(
4
),
715
732
(
2016
).
33.
Lecampion
,
B.
, and
D. I.
Garagash
, “
Confined flow of suspensions modelled by a frictional rheology
,”
J. Fluid Mech.
759
,
197
235
(
2014
).
34.
Breugem
,
W.-P.
, “
A second-order accurate immersed boundary method for fully resolved simulations of particle-laden flows
,”
J. Comput. Phys.
231
(
13
),
4469
4498
(
2012
).
35.
Costa
,
P.
,
B.
Jan Boersma
,
J.
Westerweel
, and
W.-P.
Breugem
, “
Collision model for fully resolved simulations of flows laden with finite-size particles
,”
Phys. Rev. E
92
(
5
),
053012
(
2015
).
36.
Tschisgale
,
S.
,
T.
Kempe
, and
J.
Fröhlich
, “
A non-iterative immersed boundary method for spherical particles of arbitrary density ratio
,”
J. Comput. Phys.
339
,
432
452
(
2017
).
37.
Picano
,
F.
,
W.-P.
Breugem
,
D.
Mitra
, and
L.
Brandt
, “
Shear thickening in non-Brownian suspensions: An excluded volume effect
,”
Phys. Rev. Lett.
111
(
9
),
098302
(
2013
).
38.
Shajahan
,
T.
,
T.
Schouten
,
S. K. R.
Raaghav
,
C.
van Rhee
,
G. H.
Keetels
, and
W.-P.
Breugem
, “Characteristics of slurry transport regimes: Insights from experiments and interface-resolved direct numerical simulations,” available at SSRN: https://ssrn.com/abstract=4556143 or http://dx.doi.org/10.2139/ssrn.4556143.
39.
Rastogi
,
S. R.
,
N. J.
Wagner
, and
S. R.
Lustig
, “
Microstructure and rheology of polydisperse, charged suspensions
,”
J. Chem. Phys.
104
(
22
),
9249
9258
(
1996
).
40.
Pednekar
,
S.
,
J.
Chun
, and
J. F.
Morris
, “
Bidisperse and polydisperse suspension rheology at large solid fraction
,”
J. Rheol.
62
(
2
),
513
526
(
2018
).
41.
Barnes
,
H. A.
, “
A review of the slip (wall depletion) of polymer solutions, emulsions and particle suspensions in viscometers: Its cause, character, and cure
,”
J. Non-Newtonian Fluid Mech.
56
(
3
),
221
251
(
1995
).
42.
Carotenuto
,
C.
, and
M.
Minale
, “
On the use of rough geometries in rheometry
,”
J. Non-Newtonian Fluid Mech.
198
,
39
47
(
2013
).
43.
Pawelczyk
,
S.
,
M.
Kniepkamp
,
S.
Jesinghausen
, and
H.-J.
Schmid
, “
Absolute rheological measurements of model suspensions: Influence and correction of wall slip prevention measures
,”
Materials
13
(
2
),
467
(
2020
).
44.
Chun
,
B.
,
I.
Kwon
,
H. W.
Jung
, and
J. C.
Hyun
, “
Lattice boltzmann simulation of shear-induced particle migration in plane Couette-Poiseuille flow: Local ordering of suspension
,”
Phys. Fluids
29
(
12
),
121605
(
2017
).
45.
Chun
,
B.
,
J. S.
Park
,
H. W.
Jung
, and
Y.-Y.
Won
, “
Shear-induced particle migration and segregation in non-Brownian bidisperse suspensions under planar poiseuille flow
,”
J. Rheol.
63
(
3
),
437
453
(
2019
).
46.
Di Vaira
,
N. J.
,
Ĺ.
Ĺaniewski-Wołłk
,
R. L.
Johnson
,
S. M.
Aminossadati
, and
C. R.
Leonardi
, “
Influence of particle polydispersity on bulk migration and size segregation in channel flows
,”
J. Fluid Mech.
939
,
A30
(
2022
).
47.
Lees
,
A. W.
, and
S. F.
Edwards
, “
The computer study of transport processes under extreme conditions
,”
J. Phys. C: Solid State Phys.
5
(
15
),
1921
1928
(
1972
).
48.
Lorenz
,
E.
,
A. G.
Hoekstra
, and
A.
Caiazzo
, “
Lees-Edwards boundary conditions for lattice boltzmann suspension simulations
,”
Phys. Rev. E
79
(
3
),
036706
(
2009
).
49.
Khalili
,
B.
, and
M.
Ashrafizaadeh
, “
Implementation of Lees–Edwards periodic boundary conditions for three-dimensional lattice boltzmann simulation of particle dispersions under shear flow
,”
J. Comput. Sci.
68
,
101982
(
2023
).
50.
Parsi
,
F.
, and
F.
Gadala-Maria
, “
Fore-and-Aft asymmetry in a concentrated suspension of solid spheres
,”
J. Rheol.
31
(
8
),
725
732
(
1987
).
51.
Blanc
,
F.
,
E.
Lemaire
,
A.
Meunier
, and
F.
Peters
, “
Microstructure in sheared non-Brownian concentrated suspensions
,”
J. Rheol.
57
(
1
),
273
292
(
2013
).
52.
Foss
,
D. R.
, and
J. F.
Brady
, “
Structure, diffusion and rheology of Brownian suspensions by Stokesian dynamics simulation
,”
J. Fluid Mech.
407
,
167
200
(
2000
).
53.
Batchelor
,
G. K.
, “
The stress system in a suspension of force-free particles
,”
J. Fluid Mech.
41
(
3
),
545
570
(
1970
).
54.
Guazzelli
,
E.
, and
J. F.
Morris
,
A Physical Introduction to Suspension Dynamics
(
Cambridge University
,
Cambridge
,
2011
), Vol. 45.
55.
Peerbooms
,
W.
,
T.
Nadorp
,
A. E. D. M.
van der Heijden
, and
W.-P.
Breugem
, “
Interparticle friction in sheared dense suspensions: Comparison of the viscous and frictional rheology descriptions
,”
4TU.ResearchData
, dataset (
2024
).