Scaling description of frictionless dense suspensions under inhomogeneous flow

Predicting the rheology of dense suspensions under inhomogeneous flow is crucial in many industrial and geophysical applications, yet the conventional `$\mu(J)$' framework is limited to homogeneous conditions in which the shear rate and solids fraction are spatially invariant. To address this shortcoming, we use particle-based simulations of frictionless dense suspensions to derive new constitutive laws that unify the rheological response under both homogeneous and inhomogeneous conditions. By defining a new dimensionless number associated with particle velocity fluctuations and combining it with the viscous number, the macroscopic friction and the solids fraction, we obtain scaling relations that collapse data from homogeneous and inhomogeneous simulations. The relations allow prediction of the steady state velocity, stress and volume fraction fields using only knowledge of the applied driving force.

Introduction.Dense suspensions are an important class of soft matter system comprising Brownian or non-Brownian particles mixed roughly equally by volume with viscous fluid [1].Their rheology attracts sustained interest from physicists due to the manifold complex phenomena that arise with apparently simple constituents [2,3].These include non-equilibrium absorbing state transitions [4], shear thickening [5], thinning [6], and yield stress behaviour [7].As well as being of fundamental interest, characterising this complexity is key to the extensive use of dense suspensions in various formulation and processing industries.
A useful model with which to build rheological understanding is the non-Brownian suspension [8], an especially appealing system when one considers the case of inertialess hard spheres.By analogy to dry granular systems [9], a recent study successfully obtained constitutive laws for this system [10], confirming their rateindependence and finding one-to-one relations between the volume fraction ϕ and each of two dimensionless rheological quantities, the viscous number J = η γ/P and the macroscopic friction coefficient µ = σ xy /P .Here η is the suspending liquid viscosity, γ is the shear rate, P is a measure of the particle contribution to the normal stress, and σ xy is the shear stress.This important result, the so-called µ(J)-rheology, forms the basis of subsequent models that introduce rate-dependence through additional stress scales [11,12].
The applicability of µ(J) becomes limited when considering inhomogeneous flows in which γ varies spatially [13][14][15].In particular, the lower limit of µ (which we denote µ J ) is non-zero in all homogeneously flowing systems irrespective of the particle-particle friction coefficient µ p [16][17][18] but can by construction vanish when mechanical balance dictates sign changes in σ xy such as along pipe centrelines.In such scenarios regions that would otherwise be jammed (i.e. with µ < µ J and J = 0) can have non-zero γ thanks to facilitation by nearby flowing regions [19,20].This non-local effect has been extensively studied in amorphous solids [21] and dry granular systems [22], often by formulating a fluidity field with diffusive behaviour characterised by an inhomogeneous Helmholtz equation.Microscopically it is conceptualized that the fluidity originates from an activated process that diffuses through the system in a cooperative way controlled by an inherent length scale [19,[21][22][23][24]. Recent works in dry granular matter [25][26][27] interpret the fluidity in terms of particle velocity fluctuations δu and density ρ, defining a fourth dimensionless quantity Θ = ρδu 2 /P and seeking constitutive relations linking it to ϕ, µ and I [9] (the dry counterpart to J).This successfully collapses data from homogeneous and inhomogeneous simulations onto a master curve, but is limited in that the Θ fields required to make predictions thereafter must be obtained by simulation.Naturally such findings raise the question of whether similar constitutive equations exist to unify homogeneous and inhomogeneous dense suspension rheology.
Here we use particle-based simulation [28] to model dense suspensions under homogeneous and inhomogeneous conditions, achieving the latter through an imposed Kolmogorov flow following the approach of Saitoh and Tighe [19].We seek to unify the rheology under both sets of conditions by first defining a dimensionless suspension temperature based on particle velocity fluctuations, as Θ = ηδu/aP , analogous to the granular temperature [26], and then obtaining relations among the four dimensionless numbers ϕ, J, µ and Θ.Although the µ(J) framework was devised based on frictional millimetric grains, recent experiments demonstrate it is nonetheless applicable to frictionless ones [29], and we focus here on the latter.Doing so we find scalings that can collapse homogeneous and inhomogeneous rheology data onto a set of master curves that can then be used to predict the rheology of other flow types.
Simulations details.We simulate a mixture of frictionless, non-Brownian spheres of radius a and 1.4a mixed in equal number in a periodic box of dimensions L x , L y , L z , using LAMMPS [30,31] (see Fig. 1(a)).Particles are suspended in a density (ρ) matched viscous liquid, and we impose pairwise contact and hydrodynamic forces as described by Ref. [18].Briefly, the hy-  Inhomogeneous flow of a frictionless dense suspension.Shown are (a) a typical configuration of the system for φ = 0.60, with the red region highlighting a coarse-graining box; and the steady-state profiles in y of (b) the x-components of the externally applied liquid velocity field u ∞ x (green line) and the coarse-grained velocity field of the particles ux (red points).Velocity is presented here in units of κ; (c) the expected shear rate for a Newtonian fluid γ∞ = ∂u ∞ x /∂y (green line) and the measured shear rate γ (red points), both in units of κ/a; (d) the velocity fluctuations δu in units of κ; (e) the local volume fraction ϕ, noting that the higher values at low γ demonstrate particle migration has taken place; (f) the pressure P expressed in units of ηκ/a; (g) the shear stress σxy computed from the particle interactions (red points) and by integrating over the left hand side of Eq. 4 (green points), in the same units as P .
drodynamic lubrication force for particles of radius a i and a j , with center-to-center vector r i,j , is given by , where u i,j is the relative velocity of the particles and h = (a i + a j ) − |r i,j |.F h i,j is not computed for h > 0.05a, and it saturates to ∼ (1/h c )u i,j for h ≤ h c (with h c = 0.001a), allowing particles to come into contact.Contact forces arise only when |r i,j | < (a i +a j ) and are given by where k is a spring constant and n i,j = r i,j /|r i,j |.Particles additionally experience dissipative drag due to motion relative to the fluid, given by , with u i the velocity of particle i and u ∞ (y i ) the liquid streaming velocity at the position of particle i.
Flow is generated by specifying u ∞ to induce particle motion through drag.We obtain homogeneous rheology data for fixed-volume systems of ϕ = 0.48 to 0.65 by generating simple shear via u ∞ (y) = γyδ x , with y the direction of the velocity gradient and δ x the unit vector along x.We chose our parameters such that ρ γa 2 /η ≪ 1 and γ ρa 3 /k ≪ 1, recovering rate-independence [10].To obtain inhomogeneous flow we specify a spatially dependent liquid velocity as u ∞ (y) = κ sin (2πy/L y ) δ x (see Fig. 1(b), and the gradient γ∞ in Fig. 1(c)), and later test the model with u ∞ (y) = κ sin 3 (2πy/L y )δ x .We run simulations with L y = 50a, 100a and 200a (with L x , L z = 20a) and systems containing O(10 4 ) particles (we verified that larger systems produce equivalent rheology results).We simulated systems with mean volume fraction φ = 0.5 to 0.63 (achieved by varying the particle number), and κ is a constant with dimensions of velocity, chosen so that the measured ρ γa 2 /η remains < 0.01 throughout and particle inertia is negligible.The stress (a tensor) is computed on a per-particle basis as Σ i = j (F * i,j ⊗ r i,j ), counting both contact and hydrodynamic forces.
We aim to compare the spatially-variant values of J, µ, ϕ and Θ obtained via inhomogeneous flow with the spatially-invariant ones obtained via homogeneous flow (the latter follow closely our previous results [18]).Doing so requires computing the variation in y of the stress and velocity fields under inhomogeneous flow, which we do by binning particle data in blocks of width a and volume V b = L x aL z , with the per-block value of a quantity being simply the mean of the per-particle quantities of the particles with centers lying therein.We compute the velocity fluctuation (necessary for calculating the Θ field) of each particle as δu i = |u i,x − u † i,x | where u i,x is the xcomponent of u i and u † i,x is the average x velocity of all particles with centers lying in a narrow window ±ϵ (taking ϵ = O(0.1a)) of y, and we then bin δu i per block.As all three components of the velocity fluctuations are statistically equivalent we have used only the x values to compute Θ.In what follows we report steady state data only [32], averaging across 6 realizations and at least 500 configurations per realization.
Results.Shown in Fig. 1 (e) ( f Relations between the dimensionless control parameters.Shown in the top row are the relations between the dimensionless viscous number J and (a) the volume fraction ϕ for a range of homogeneous ϕ (black data) and inhomogeneous φ; (b) the effective friction coefficient µ and; (c) the suspension temperature Θ.In the bottom row are the collapses using the scaling Eqns. 1 (d), 2 (e) and 3 (f), for different φ and L. In (d) we show data for L/a = 50 to highlight its deviation from the scaling relation.Black triangles represent homogeneous data (simple shear) and all other points are for inhomogeneous flow at different φ.
state ϕ exhibits spatial variation set up by particle migration to balance the normal stress [13,14,33].The velocity profile follows a similar trend to the applied force, as expected, but is flattened at the regions of largest ϕ leading to significant deviations between γ and γ∞ .The pressure becomes spatially uniform, and the shear stress follows the shear rate in sign.Since P is spatially invariant in the steady state, one can deduce that the variation of the quantities η γ/P , σ xy /P and ηδu/aP follow γ, σ xy and δu respectively.
We analyse inhomogeneous data by computing the dimensionless control parameters in each block, defining the scalar shear rate and stress components on the basis of invariants of the respective tensor quantities so that J, µ > 0. This is done for a range of φ, with parametric plots of J(y), ϕ(y), µ(y) and Θ(y) shown in Figs.2(a)-(c).Each plotted point represents a y-coordinate, and colors represent different φ.Shown also (in black) are homogeneous data.Reading across the data points of a single color from right-to-left represents moves from regions of high-to-low γ in the inhomogeneous domain.
The homogeneous ϕ(J) and µ(J) relations follow qualitatively the result of Boyer et al. [10], though our frictionless particles render ϕ J and µ J dissimilar.Θ(J) follows a power-law relation, as in dry granular matter [26] though with a different exponent.In general large-J inhomogeneous data approximately match homogeneous data, though they deviate with decreasing J demonstrating the shortcomings of the existing constitutive laws.
With the help of scaling theory, we next attempt to find constitutive laws that simultaneously describe the rheology under homogeneous and inhomogeneous flow.We focus first on how the inverse viscosity J/µ = η γ/σ xy vanishes as ϕ approaches the jamming point ϕ J .This trend is followed by all the homogeneous and inhomogeneous simulations, leading to our first scaling relation plotted in Fig. 2(d) with α = 4.1 and ϕ J = 0.6555.The next scaling relation is motivated by Kim and Kamrin [26].In homogeneous flow, within the range of our data we find µ 2.5 ∼ J (Fig. 2(b)) and Θ 1.44 ∼ J (Fig. 2(c)).Since for the range of φ explored here inhomogeneous data follow homogeneous laws at large J, we expect a scaling of the form µ 2.5 Θ 1.44 ∼ F 1 (J).Indeed this results in a good collapse as shown in Fig. 2(e), in which data are described by the relation with β = 3 and ϑ = 0.06.The final scaling relation is motivated by the relation between granular fluidity and ϕ reported for dry granular matter.Zhang and Kamrin [25] write a non-dimensional granular fluidity g = gd/δu, where g = γ/µ, and d is the spatial dimension.We define an equivalent quantity in terms of the previously discussed dimensionless numbers, namely J/µΘ, though we find a better collapse is achieved through a change to the exponents as with indicating limits to the range of applicability.An issue in the former case may be that our simplified hydrodynamics, accounting only for lubrication, becomes nonphysical at lower ϕ and that a more highly resolved fluid field is required.
Given a profile of one of the dimensionless numbers, one could therefore fully characterise the rheology of the system.In our simulations, however, the only known input is the externally applied force, which we recall is defined through u ∞ .To use the scaling relations we need to establish another relation that can provide us one of these dimensionless numbers from the knowledge of the applied force profile.Considering the inertia-free momentum balance ∇ • Σ = −f per unit volume, we can write the following equation for the k th block of the simulation cell (which we verified in Fig. 1(g)): Here N k , u ∞ x,k , u x,k and σ xy,k are the particle number in the block, the liquid streaming velocity at the centre of the block, and the particle velocity and stress averaged over the block, which has volume V b .A is an order unity quantity necessary to account for small variations in u ∞ x,k across the block.The first term of Eq. 4 represents the net applied force and the second represents the net viscous force exerted by the fluid due to drag.The resultant of these is balanced by the net stress gradient inside the block.Using the definition of our dimensionless numbers, Eq. 4 can be rewritten for the streaming velocity at y as with u ′∞ x (y) = u ∞ x (y)ηA/aP and asterisks representing multiplication by sgn( γ∞ (y)), noting that P is uniform at steady state and using ϕ(y) = (4/3)πa 3 N (y)/V b , acknowledging our earlier comment about phase separation [32].Equation 5 thus relates the externally applied liquid flow field to the profiles of J, µ and ϕ.
For a known u ∞ we solve Eqs. 1, 2, 3 and 5 numerically in the following way.We first guess a ϕ (y) profile by assuming accumulation at points where the spatial derivative of the imposed force vanishes, starting with a simple form as ϕ(y) = np j=1 a j /[(y − y 0 j ) 2 + b 2 j ] + ϕ 0 , with mass conserved through φ = 1 Ly Ly 0 ϕ (y) dy.Here y 0 j are the coordinates of the point where the first derivative of the applied force vanishes, n p is the number of such points and b j is the width of the Lorentzian function peaked at y 0 j .We then compute directly J, µ and Θ using Eqs. 1, 2 and 3, before attempting to balance Eq. 5.The imbalance of Eq. 5 reflects the accuracy of our guess.We refine ϕ(y) by tuning ϕ 0 , a j and b j until Eq. 5 is satisfied (up to some tolerance).Shown in Fig. 3 are predicted results compared against 'unseen' simulation data (i.e.data not used to obtain the scaling exponents) with φ = 0.55, 0.57 and u ∞ (y) = κ sin 3 (2πy/L y )δ x , demonstrating the degree of success of the scaling relations for predicting y-profiles of ϕ, J, µ and Θ. Considering the highly nonlinear nature of the scaling relations, the quality of the predictions is reasonably good.Conclusions.Using particle-based simulation we seek universality amongst flows of dense, frictionless suspen-sions.Along with canonical suspension rheology control parameters ϕ, J and µ, we introduce a fourth quantity Θ characterising velocity fluctuations, inspired by recent studies in dry granular physics [26].We find a trio of scaling relations among these quantities that collapse data for homogeneous and inhomogeneous flow.Utilising a momentum balance we show that from knowledge of the externally applied force, one can use the relations to predict the features of a general inhomogeneous flow.Our work raises manifold avenues for future work.In particular, the microscopic origin of the exponents is not understood, nor is their generalisation to the broader class of suspensions that includes polydisperse particles (for which colloidal forces may become relevant [34]), nonspheres and other complexities.Meanwhile the question of a diverging lengthscale -apparently a staple of nonlocal rheology in dry granular matter [22,24,35]-remains open.Computing a granular fluidity field from our data, we find, similar to [19] (b)-(g) are, respectively, steady-state profiles in y of the coarse-grained velocity (in x) u x , shear rate γ = ∂u x /∂y, velocity fluctuations δu, volume fraction ϕ, pressure P (= (1/3)Tr(Σ)), and shear stress σ xy , for φ = 0.60, with each plotted point representing a block.Although at initialisation the particle density is homogeneous (i.e.ϕ ̸ = ϕ(y)), in the steady

Fig. 2 (
Fig.2(f)) and ϵ = −10.98,ϕ m = 0.618, ζ = 0.0004 and λ = 1.533.We thus have three scaling relations, Eqs. 1, 2 and 3, that relate ϕ, J, µ and Θ.The collapse appears poorer for φ = 0.5 (Fig.2(f)) and L/a = 50 (Fig.2(d)), indicating limits to the range of applicability.An issue in the former case may be that our simplified hydrodynamics, accounting only for lubrication, becomes nonphysical at lower ϕ and that a more highly resolved fluid field is required.Given a profile of one of the dimensionless numbers, one could therefore fully characterise the rheology of the system.In our simulations, however, the only known input is the externally applied force, which we recall is defined through u ∞ .To use the scaling relations we need to establish another relation that can provide us one of these dimensionless numbers from the knowledge of the applied force profile.Considering the inertia-free momentum balance ∇ • Σ = −f per unit volume, we can write the following equation for the k th block of the simulation cell (which we verified in Fig.1(g)):

FIG. 3 .
FIG.3.Predictions of the scaling relations against simulation data not used for obtaining the scaling exponents, with u ∞ (y) = κ sin 3 (2πy/Ly)δx.Shown are (a) the volume fraction ϕ; (b) the viscous number J; (c) the effective friction coefficient µ; and (d) the suspension temperature Θ, with predictions given by solid lines and simulation data in points, for φ = 0.55 (red) and 0.57 (green).
, no divergence in the characteristic lengthscale, which remains O(a) everywhere.This raises an important open question regarding what are the minimal conditions required for a diverging lengthscale in inhomogeneous particulate flows.B.P.B. acknowledges support from the Leverhulme Trust under Research Project Grant RPG-2022-095; C.N. acknowledges support from the Royal Academy of Engineering under the Research Fellowship scheme.We thank Ken Kamrin, Martin Trulsson, Mehdi Bouzid, Romain Mari and Jeff Morris for useful discussions.