Unifying homogeneous and inhomogeneous rheology of dense suspensions Open Access

Dense suspensions consist of Brownian or non-Brownian solid particles suspended in viscous fluids in roughly equal proportions. Such materials, for example, cornstarch in water, slurries, blood, etc., have widespread applications in both daily life and industry, and their flow is observed in many natural phenomena. Therefore, understanding their rheology is pivotal [1–4]. However, the rheology of these materials is extremely complex and lacks a universal description due to its dependence on various parameters such as the size, asphericity, and surface roughness of the solid particles, the nature of the suspending fluid, as well as the complexity involved in the external driving. Various constitutive models have been introduced based on defining appropriate dimensionless numbers, but their validity is often limited to specific scenarios [5–10]. Therefore, identifying suitable quantities to establish constitutive laws to describe the rheology of different types of dense suspensions under various driving conditions has been a topic of intense research [6,11–17].

Under the homogeneous scenario, where the strain rate is spatially uniform, the rheology of dense suspensions is well described by three dimensionless quantities: the solid volume fraction $ϕ$⁠; the ratio of the viscous time scale $η/P$ to shear time scale $1/ γ ˙$ known as viscous number $J=η γ ˙/P$⁠; and the ratio of shear stress $σ x y$ to normal stress $P$ known as effective or macroscopic friction coefficient $μ= σ x y/P$ [6]. Here, $η$ and $γ ˙$ are suspending the fluid viscosity and strain rate, respectively. The interdependence of these dimensionless numbers forms the constitutive laws for the homogeneous rheology of dense suspensions. However, the validity or structure of such constitutive laws might alter depending on the various properties of the constituent particles. Nevertheless, in real systems, the nature of the rheology is often inhomogeneous due to the spatial variation of the strain rate, so the aforementioned dimensionless numbers are not sufficient [11]. Under homogeneous straining, $J$ decreases with decreasing $μ$ and increasing $ϕ$ and, eventually, vanishes when $μ$ and $ϕ$ approach their limiting value $μ J$ and homogeneous shear jamming volume fraction $ϕ J H$⁠, respectively. This physically signifies the cessation of the flow, however, in the case of inhomogeneous flow, one can observe $ϕ> ϕ J H$ and $μ< μ J$ even for finite $J$ due to shear-induced particle migration [18–28] and nonlocal effects [12,18,29–31].

Nonlocal phenomena in the rheology of various soft materials are studied extensively and pictured as a process where flow in regions with $μ> μ J$ facilitates flow in regions with $μ< μ J$ via diffusion of local fluidity of the system [32]. Such diffusion of local fluidity can be described by an inhomogeneous Helmholtz-like equation, which suggests a cooperative motion controlled by an inherent length scale [33,34]. In the context of dry granular systems, in [35], granular fluidity is macroscopically defined as $g= γ ˙/μ$⁠. Later, in [36], Zhang and Kamrin established the microscopic definition of fluidity, denoted as $g ~$⁠, in terms of fluctuations of particle velocity, $δu$⁠. $g ~$ is uniquely determined by the local $ϕ$ and can be expressed as $g ~=ga/δu=F ( ϕ )$⁠, where $a$ is the particle diameter. Such fluidity is found to be independent of $ϕ$ at low volume fraction [37,38] but decreases with $ϕ$ at sufficiently large volume fraction, vanishing as $ϕ$ approaches random close packing, $ϕ RCP$ [39], irrespective of whether the flow is homogeneous or inhomogeneous. Moreover, in [14], a new dimensionless number, granular temperature $Θ=ρδ u 2/DP$⁠, is introduced. Using power-law scaling, it is demonstrated that $μ$⁠, properly scaled by $Θ$⁠, unifies homogeneous and inhomogeneous rheology with the inertial number $I$ being the scaling variable, thus replacing the conventional $μ(I)$ rheology by a $μ(I,Θ)$ rheology for inhomogeneous flows in dry granular systems. Here, $I$ is the counterpart of $J$ for dry granular systems and $D$ is the spatial dimension.

Similarly, motivated by the concept of $Θ$ in [14], a recent study [13] defined a new dimensionless quantity, the suspension temperature $Θ s=ηδu/aP$⁠. Higher values of $Θ s$ stand for higher mobility of the particles, thus enabling a higher degree of softness of the system from the rheological point of view. Additionally, the suspension temperature measures the competition between convective flow, driven by external forcing, and diffusive flow driven by collisions with other particles. When particle trajectories are governed mostly by collisions, their trajectories will deviate more from the affine flow, so their suspension temperature will be large. Using this novel quantity along with the other dimensionless quantities $ϕ,J,$ and $μ$⁠, a set of constitutive laws is identified that unifies the homogeneous and inhomogeneous rheology of dense suspensions of frictionless spherical particles. However, such an approach remains unexplored for frictional and aspherical particles. Interestingly, the constitutive laws for homogeneous rheology differ for frictional and aspherical particles compared to frictionless spherical particles. For frictional particles, as the particle friction coefficient $μ p$ increases, the sliding between particle pairs becomes increasingly restricted, imposing constraints on both rotational and translational degrees of freedom. This, in turn, leads to a reduction in $ϕ J H$ [40]. With asphericity, the effect on $ϕ J H$ is nonmonotonic. Specifically, for rod-shaped particles, $ϕ J H$ decreases when the aspect ratio (⁠ $AR=$ length/diameter) exceeds an intermediate value of approximately 1.5. This decrease is due to increased entanglement, or a higher number of contacts per particle, which results in inefficient random packing. However, when the aspect ratio is below this threshold, $ϕ J H$ increases because the reduced entanglement allows for more efficient packing [41–45]. Therefore, the validity of $Θ s$ for systems of frictional and aspherical particles remains unclear and requires further investigation.

In this work, using particle-based simulations, we unify the homogeneous and inhomogeneous rheology of dense suspensions of frictional spherical particles and frictionless rod-shaped particles [46–50]. In combination with the suspension temperature defined in [13] for frictionless spherical particles, along with other dimensionless numbers used to describe the homogeneous rheology, we identify new scaling relations that collapse the data of the homogeneous and inhomogeneous rheology. We find that some of the scaling relations identified here retain the same mathematical form but exhibit different exponent values across systems involving frictionless and frictional spherical particles, as well as frictionless rod-shaped particles. However, other scaling relations apply only to specific systems, highlighting intrinsic differences in their rheology. We further validate these scaling relations by demonstrating their ability to predict various dimensionless numbers in previously unexamined simulation results.

In Figs. 1(b)–1(f), steady-state profiles in $y$ of the coarse-grained particle velocity $u x$ (flow direction), strain rate $γ ˙=∂ u x/∂y$⁠, velocity fluctuations $δu$⁠, pressure $P$ (⁠ $=(−1/3)Tr( Σ)$⁠) and the normal stresses, and the shear stress $σ x y$ are shown for frictional spherical particles with $ϕ ¯=0.60$ and $μ p=0.1$⁠, with each plotted point representing a block. The particle velocity profile and applied streaming velocity follow a similar trend, as expected, but the former is flattened at the regions of largest $ϕ$ leading to significant deviations between $γ ˙$ and $γ ˙ ∞$ (⁠ $=∂ U x ∞/∂y$⁠). The pressure becomes spatially uniform in the steady state, and the normal stresses exhibit weak anisotropy at the regions where $μ> μ J$⁠. The shear stress follows a similar spatial variation to the shear rate.

In Fig. 2, the spatial variation of the dimensionless numbers in the steady state is presented for two different $ϕ ¯$⁠, close to and far from $ϕ J H$⁠. In Fig. 2(a), the viscous number $J$ is presented. Since the pressure is uniform in the steady state, $J$ looks similar to $γ ˙$ shown in Fig. 1(e). Although we start our simulation from a uniform volume fraction, with straining, particles move toward the region with smaller strain rate due to normal stress $σ y y$ imbalance and accumulate there [21,22,60]. In the steady state $ϕ$ attains a maximum at $J=0$ and decreases as $J$ increases. $μ$ and $Θ s$ have a similar variation of $σ x y$ and $δu$⁠, respectively.

The spatial variation of $∂ γ ˙/∂y$⁠, as shown in Fig. 1(e), highlights the inhomogeneous nature of the flow in our setup, as for homogeneous flow, $∂ γ ˙/∂y=0$ throughout. Given that the velocity profile $U ∞$ follows a sinusoidal pattern [as shown in Fig. 1(b)], regions with larger $J$ have smaller $∂ γ ˙/∂y$⁠, suggesting that the data from these regions might align with homogeneous (simple shear) flow data. However, regions where both $J$ and $∂ γ ˙/∂y$ are small are less likely to correspond to homogeneous data. These regions exhibit reduced inhomogeneity, as the entire region moves together (creeping motion), as evidenced by the flatness of the velocity profile in Fig. 1(b). The $ϕ$ and $μ$ here go above and below their homogeneous limit, $ϕ J H$ and $μ J$⁠, respectively, and the local flow is primarily controlled by $Θ s$⁠, as we will see later. Moreover, in [13], the inhomogeneity in the flow is quantified for a setup similar to the one studied here by demonstrating a growing length scale associated with the cooperative diffusion of fluidity. This underpins the true inhomogeneous nature of the flow in our setup.

Here, $ν≈2$⁠. For homogeneous rheology $ϕ J= ϕ J H$ monotonically decreases from $ϕ RCP$ with increasing $μ p$ [40]. For inhomogeneous flow, $ϕ J= ϕ J I H$ independent of $μ p$ and is found to be close to $ϕ RCP$⁠. However, this power law divergence of viscosity for inhomogeneous flow is only valid for $ϕ< ϕ J H( μ p)$⁠, although the viscosity diverges at $ϕ R C P$⁠. $η 0$ is a $μ p$ dependent coefficient with different values for homogeneous $( η 0 H )$ and inhomogeneous $( η 0 I H )$ rheology. We find $( η 0 I H/ η 0 H)≈2$ and $3$ for $μ p=0.1$ and 0.5, respectively. The reason to have such dependence is the following. Since some part of our inhomogeneous simulation box is strained homogeneously (at large $J$⁠) data from this region fall on top of the homogeneous flow data. Thus, we have two different power laws for homogeneous and inhomogeneous flow, which start from the same point and diverge with the same exponent but at different volume fractions. Therefore, in our scaling relation, we have a prefactor that not only depends on the $μ p$ but also on the nature of the rheology.

Thus, we effectively have three scaling relations. The second and third scaling relations, given by Eqs. (15) and (17), are valid across a wide range of volume fractions, both above and below $ϕ J H$⁠. In contrast, the first scaling relation is divided into two parts. The first part, Eq. (11), applies for $ϕ< ϕ J H$⁠, while the second part, Eq. (13), is valid above $ϕ J H$⁠. Later, we will demonstrate how these scaling relations, combined with the momentum balance equation, allow us to predict all the relevant dimensionless numbers based solely on the applied fluid flow. In addition to the scaling relations discussed above, using Eqs. (13) and (17), it is possible to deduce another relation $μ 1.1∼ F 1 b(ϕ)/ F 3(ϕ)$⁠. The data collapse for $μ p=0.1$ using this scaling relation is presented in the inset of (d).

For the system of frictional spherical particles, we find decoupling of homogeneous and inhomogeneous shear jamming volume fraction (i.e., $ϕ J H< ϕ RCP≈ ϕ J I H$⁠) due to the frictional constraints. Similar decoupling is also expected for rod-shaped particles due to the constraints imposed by asphericity quantified by the aspect ratio. The relations between different dimensionless numbers for this system, shown in Fig. 5, are similar to those for frictional spherical particles. The homogeneous data follow the same functional form given by Eqs. (7)–(10) and are, therefore, not discussed here. Both $ϕ J H$ and $μ J$ exhibit nonmonotonic dependence on $AR$⁠, with an maximum value of $AR≈1.5$⁠, shown in the inset of Fig. 5(b).

Given $U x ∞$⁠, we solve Eq. (26) and the scaling relations [Eqs. (11), (13), (15), and (17)] for frictional spherical particles, and Eqs. (19), (21), and (23) for frictionless rod-shaped particles] numerically by the following method. Initially, we assume a $ϕ ( y )$ profile, hypothesizing accumulation of particles at points where the spatial derivative of $U ∞$ is zero, starting with a simple Lorentzian form $ϕ(y)= ∑ k = 1 n p a k/$ $[ ( y − y k 0 ) 2+ b k 2]+ ϕ 0$, with mass conserved through $ϕ ¯=(1/ L y) ∫ 0 L yϕ ( y )dy$⁠. Here, $n p$ is the number of points where the first derivative of $U ∞$ is zero, $y k 0$ is the coordinate of such a point, and $b k$ is the width of the Lorentzian function centered at $y k 0$⁠. Next, we compute $J$⁠, $μ$⁠, and $Θ s$ directly using the scaling relations, before attempting to balance Eq. (26). The imbalance in Eq. (26) indicates the accuracy of our initial guess. We refine $ϕ(y)$ by adjusting $ϕ 0$⁠, $a k$⁠, and $b k$ until Eq. (26) is satisfied within an acceptable tolerance.

The results, shown in Fig. 7, compare predicted outcomes against previously unseen simulation data (i.e., data not used for obtaining the scaling exponents) with $ϕ ¯=0.595$ for frictional spherical particles, $0.615$ for frictionless rod-shaped particles, and $U ∞(y)=κsin(2πy/ L y) x ^$⁠, demonstrating the success scaling relations in predicting $y$ profiles of $ϕ$⁠, $J$⁠, $μ$⁠, and $Θ s$⁠. Despite the highly nonlinear nature of the scaling relations and many orders of magnitude spread of $J$ and $Θ s$⁠, the predictions are reasonably accurate.

Through particle-based simulations, we establish a universal description of the flow behavior of dense suspensions, of frictional spherical and frictionless rod-shaped particles. In addition to the standard control parameters solid volume fraction $ϕ$⁠, viscous number $J$⁠, and macroscopic friction coefficient $μ$⁠, we introduce a novel parameter, suspension temperature $Θ s$⁠, representing velocity fluctuations, inspired by concepts from dry granular materials. Our findings reveal scaling relations among these parameters that successfully collapse data for both homogeneous and inhomogeneous flows. Using the momentum balance, we demonstrate that the characteristics of general homogeneous and inhomogeneous flow can be predicted based on the applied external force. It is important to note that in this work, we primarily focused on finding the existing scaling relations from the simulation data rather than a thorough investigation of the physical origin of them. Other choices of dimensionless quantities, for instance, based on the rotational degrees of freedom may warrant further investigation. The exponents reported here are found by using an ad hoc method to obtain the best data collapse. The physical origin of these exponents and the validity of the identified scaling relations in flows with more complex geometries such as hopper flow [38] remain unexplored here but represent natural and important avenues for future work, alongside extending the proposed scaling relations to situations where the gradient of the flow rate extremely large [13].

B.P.B. acknowledges support from the Leverhulme Trust under Research Project Grant No. RPG-2022-095; C.N. acknowledges support from the Royal Academy of Engineering under the Research Fellowship scheme. We are grateful to Anoop Mutneja, Eric Breard, and Ken Kamrin for discussions.

The authors have no conflicts to disclose.

All the scripts used to generate the data reported in this study are available from the corresponding author upon reasonable request.

Shown in Tables I and II are the numerical values of fitting parameters used in the scaling relations reported above.