Simulating the rheology of dense suspensions using pairwise formulation of contact, lubrication and Brownian forces

Abstract Dense suspensions of solid particles in viscous liquid are ubiquitous in both industry and nature, and there is a clear need for efficient numerical routines to simulate their rheology and microstructure. Particles of micron size present a particular challenge: at low shear rates, colloidal interactions control their dynamics while at high rates, granular-like contacts dominate. While there are established particle-based simulation schemes for large-scale non-Brownian suspensions using only pairwise lubrication and contact forces, common schemes for colloidal suspensions generally are more computationally costly and thus restricted to relatively small system sizes. Here, we present a minimal particle-based numerical model for dense colloidal suspensions that incorporates Brownian forces in pairwise form alongside contact and lubrication forces. We show that this scheme reproduces key features of dense suspension rheology near the colloidal-to-granular transition, including both shear thinning due to entropic forces at low rates and shear thickening at high rates due to contact formation. This scheme is implemented in LAMMPS, a widely used open source code for parallelised particle-based simulations, with a runtime that scales linearly with the number of particles, making it amenable for large-scale simulations.


I. INTRODUCTION
Dense suspensions of Brownian and non-Brownian solid particles in viscous liquid present intriguing flow properties, and understanding their rheology is a subject of both fundamental and technological relevance [1,2].Of particular interest are suspensions comprising particles with radius a ≈ 1 µm, or more broadly in the range 0.1−10 µm.These are present in numerous applications, in all areas of food science [3] and consumer products, as well as across manufacturing and construction sectors [4] and indeed in many geophysical contexts [5].Often their physics are challenging because their Brownian diffusion time may be comparable to the processing or macroscopic timescales involved in their use, so that they sit at the boundary of colloidal and granular systems [6].
Particle-based simulation offers a promising route to better understand the physics of these materials, providing highly-resolved information complementary to what can be obtained by experiment.With simultaneous access to particle trajectories and bulk rheology, one might devise new micromechanical constitutive equations [7] or develop microstructural insight that could guide the future analysis of experimental data.Numerical models might also be useful for exploring parameter space and systematically linking aspects of particle-level physics (friction [8], adhesion [9], roughness [10]) to bulk flow behaviour.As a result, one might aim to optimise industrial processes such as mixing and extrusion or indeed to optimise the design of materials themselves through additives, using insight gained through particle-based simulation.
Stokesian Dynamics (SD) [11,12] is a computational method used to simulate the rheological behavior of colloidal and granular particles suspended in a viscous fluid, addressing the special case of inertia-free flow where the Stokes number is zero [13].The method involves balancing all of the forces on each particle by evaluating their velocities via a grand mobility matrix containing information on the relative positions of every particle in the system.Despite accurately capturing the long and short range hydrodynamic interactions between particles, SD has not been adopted widely as a predictive tool in applied and industrial settings in the same way as other particle-based simulation methods have, due to the complexity of its implementation and its computational expense (notwithstanding recent developments that significantly speed it up [14,15]).
The discrete element method (DEM) [16], on the other hand, is a particle-based computational method (a variant of molecular dynamics) that is used widely to simulate the behavior of granular materials including powders, particles and grains taking into account their pairwise interactions.In contrast to SD, DEM does not balance forces on each particle.Instead inertia is present, and one simply sums the forces and the resultant leads to an acceleration that can be realised through a conventional time-stepping algorithm such as Velocity-Verlet.This approach has proven to be useful for studying overdamped suspensions under shear flow nonetheless [17,18], where one introduces short ranged lubrication forces and sets the Stokes number to be O(10 −2 ) or smaller.This approach is pragmatic in the sense that the physics associated with flowing dense suspensions can be implemented in existing, widely-used codes with large user bases, so that they have a clear path to adoption in engineering and other applied contexts.To date there is not, to our knowledge, a discrete element method simulation that includes the relevant physics of dense suspensions at the colloidal-to-granular interface, accounting for short-ranged hydrodynamics, Brownian forces, and (frictional) particle-particle contacts.
Here we present a minimal particle-based simulation model for predicting the rheology of dense Brownian and non-Brownian suspensions.Our model comprises hydrodynamic lubrication, particle-particle contacts and Brownian forces.After first describing the model in detail, we present some aspects of the effective interactions and diffusion that arise, before giving a detailed account of the rheological predictions of the model.The model reproduces well the main features of the experimentally observed rheology of dense suspensions, namely a low shear rate plateau that gives way to shear thinning and later shear thickening as the shear rate is increased, with the relative viscosity of the suspension increasing sharply with solids volume fraction and particle-particle friction coefficient.

II. METHODOLOGY
We consider a model system of nearly-monodisperse solid spheres, dispersed at high solids volume fraction ϕ in a density-matched Newtonian liquid.The microscopic physics included in our model represent a minimal set of ingredients necessary to make useful predictions of the rheology of suspensions comprising particles with radius in the range 10 −7 to 10 −4 m.The trajectories of individual particles with translational and rotational motion are governed by Langevin equations that comprise three force (F ) and torque (T ) contributions: direct particle contacts (F C , T C ), hydrodynamics (F H , T H ), and Brownian noise (F B , T B ).The equations of motion for translation and rotation of the particles are written, respectively, as where x i represents the position of particle i, Ω i represents its rotational velocity, and a i and m i are its radius and mass respectively.The subscript i represents single-body forces and torques acting on particle i, while the subscript i, j represents pairwise forces and torques acting between particles labelled i and j.The superscripts C, H, B, D, L refer to the force and torque components arising due to contacts (C), hydrodynamics (H) and Brownian (B) effects, with the latter two acting both through drag (D) and lubrication (L).Each of these force and torque terms is described in detail below.These equations of motion can be understood as Langevin equations in which the ⟨•⟩ C terms represent particle-particle interactions; the ⟨•⟩ H terms represent configuration-dependent viscous friction (i.e.dissipative forces linear in the particle velocities); and the ⟨•⟩ B terms represent configuration-dependent (multiplicative) noise.Although particle inertia is present in the model we omit fluid inertia [19], arguing that for the regimes of interest the principle contributions to the overall bulk rheology will come from particle-particle contact and hydrodynamic lubrication interactions.Particles are subjected to a liquid flow field given by U ∞ (acting through the body force F H,D i as described below), leading to a rate of strain tensor E = 1 2 ∇U ∞ + (∇U ∞ ) T .Pairwise forces and torques are summed over the neighbours j of each particle i, and the positions, velocities and acceleration are updated in a stepwise manner following the Velocity-Verlet algorithm [20].Below we describe each of the force and torque contributions in detail; shown in Figure 1 are illustrative schematics of each of the forces.

A. Contact forces and torques
The particle-particle contact force F C follows a conventional granular-type interaction [16], and is activated for any two particles i and j for which the centre-to-centre distance |r i,j | is smaller than the sum of the radii a i + a j .Contact forces include a repulsive part acting normal to the pairwise centre-to-centre vector r i,j (we define a unit vector n i,j = r i,j /|r i,j |), and a tangential part.For simplicity we model contacts as linear springs, so that particle pairs experience repulsive contact forces proportional to their scalar overlap, defined once in contact as δ i,j = (a i +a j )−|r i,j |.The implementation of our model within LAMMPS [21] nonetheless allows straightforward implementation of more complex δ i,j dependence.Tangential forces are linear in ξ i,j , a vector describing the accumulated displacement of the particle pair perpendicular to n i,j since the initiation of the contact.Contact force and torque magnitudes are controlled by normal and tangential stiffness constants k n and k t that set the hardness of the particles.The force and torque are given respectively by: We additionally introduce a static friction coefficient µ that constrains the tangential force to |k t ξ i,j | ≤ µk n δ i,j .For larger values of k t ξ i,j the tangential part of the force and the torque are truncated, and particle contacts transition from a rolling to a sliding regime.We present data for µ = 0 throughout, except in Figure 6(b)-(c) where we explore the role of contact friction.Each pairwise contact between particles i and j contributes to the overall contact stress of the system with a tensorial stresslet given by the outer product −F C i,j ⊗ r i,j .The contact stress Σ C i,j is obtained by summing this quantity over all contacting particle pairs and dividing by the system volume and dimension.
Contact forces of this kind have successfully been deployed in numerical models for rate-independent granular suspension rheology [17,22,23] and for models of shear thickening suspensions [8] (in the latter case rate dependence arises from a 'critical load' that the contact force must exceed before static friction is activated).

B. Hydrodynamic forces and torques
In general, hydrodynamic interactions in suspensions appear as single particle drag forces F H,D , pairwise nearcontact lubrication forces F H,L , and many-body long range forces.In high volume fraction dense suspensions, however, it is argued by many authors that the hydrodynamic interactions are dominated by near-contact lubrication interactions [24] (which diverge on close approach) and that long range interactions are effectively screened by intervening particles [8,25].We follow this reasoning and therefore omit long range hydrodynamics from our model.Below we describe in detail the drag and lubrication forces deployed in the model.Single particle drag forces and torques are given by T H,D where we use the isolated-particle Stokes terms and, for simplicity, do not introduce volume fraction dependent hindrance functions.Here η is the liquid viscosity, U ∞ (x i ) is the value of the liquid streaming velocity at the position of the centre of mass of particle i, and Ω ∞ = 1 2 (∇ × U ∞ ) (spatially uniform assing U ∞ is uniform in space).The drag forces lead to a per particle stress given by Σ H,D i = 20  3 πηa 3 i E. For pairwise lubrication forces and torques acting between interacting particles i and j we start from the conventional representation given by Kim and Karrila [26] as where R is the resistance matrix containing tensorial operations that linearly couple particle forces (torques) to velocities (rotational velocities), taking into account relative particle positions.After some algebra and omitting terms that vanish with the size of the interparticle gap (see Radhakrishnan [27] for details) one can obtain the forces in a simplified pairwise form as where F H,L i,j is the force acting on particle i by particle j; N = n i,j ⊗n i,j is a tensorial normal operator; T = I−n i,j ⊗n i,j is a tensorial projection operator; n i,j is the unit vector pointing from particle i to particle j; U i is the velocity of particle i; Ω i is the rotational velocity of particle i; and I is the identity tensor in three dimensions.The scalar prefactors X and Y encode the geometry of the interacting pair, namely the size of the interparticle gap and the size ratio of the interacting particles.Their superscripts A, B and subscripts 11, 22 are more appropriate to the labelling convention used by Kim and Karrila [26] but nonetheless we retain them here for ease of referencing to that work.The particle size ratio is written as β = a j /a i and the dimensionless interparticle gap is ξ = 2 (|r i,j | − (a i + a j )) /(a i + a j ).The scalar prefactors are given by Meanwhile the torques on particle i and j as a result of their interaction with particles j and i respectively are written as with scalar prefactors given by Similar expressions may be obtained for the elements of the hydrodynamic lubrication stress tensor, though these can be shown to be equivalent (up to an order ξ term in the normal stresses) to the form used for the contact forces.The contribution to the hydrodynamic stress coming from each pairwise interaction is thus given by Σ H,L i,j = −F H,L i,j ⊗ r i,j .To mitigate against divergence in the scalar prefactors at particle contacts (that is, where ξ → 0) we use ξ eff = 10 −3 in the calculation whenever ξ < 10 −3 .We do not calculate pairwise lubrication forces when particles are separated by gaps ξ > 0.05, having verified that this choice does not affect our conclusions.

C. Brownian forces and torques
To satisfy fluctuation-dissipation theorem, we must produce Brownian forces that follow where F B is a list of the Brownian forces and torques, R is the overall resistance operator for the system (taking into account both one body and pairwise hydrodynamic dissipation that we describe separately below), k b T is the thermal energy and ∆t is the computational timestep (discussed in more detail below).
For one-body Brownian forces we need 6 random numbers (i.e. two vectors in three-dimensional space ψ i and φ i ) to satisfy the translational and rotational degrees of freedom of each particle i.The elements of the random vectors ψ i , φ i are drawn from a Gaussian distribution and satisfy ⟨φ α φ β ⟩ = ⟨ψ α ψ β ⟩ = δ αβ and they are uncorrelated with each other so that ⟨φ α ψ β ⟩ = 0.The following forces and torques satisfy fluctuation-dissipation theorem (we label them as Brownian drag 'B,D' to align with the hydrodynamic drag forces and torques defined above).The one-body Brownian force and torque on particle i are given by: ⟩ over many realisations of the vectors ψ i and φ i leads, respectively, to 2k b T ∆t 6πηa i I and 2k b T ∆t 8πηa 3 i I as required (with I the identity matrix in three dimensions).Pairwise Brownian forces and torques similarly require two random vectors θ i,j and χ i,j (independent of ψ i and φ i but with the same properties) to satisfy the relative translational and rotational motion of two interacting particles.The pairwise forces and torques must be constructed in such a way that, for particles i and j, averaging ⟨F B,L i,j ⊗F B,L i,j ⟩ and ⟨T B,L i,j ⊗ T B,L i,j ⟩ over many realisations of θ i,j and χ i,j recovers the form of the pairwise hydrodynamic lubrication forces and torques described above.Doing so, which involves exploiting that the normal and projection operators present in the definition of the lubrication forces and torques are idempotent (i.e.⟨(N i,j θ i,j ) ⊗ (N i,j θ i,j )⟩ = N i,j ) and orthogonal (i.e.⟨(N i,j θ i,j ) ⊗ (T i,j θ i,j )⟩ = 0), one obtains the following expressions for the pairwise Brownian force and torque Our model thus involves computing Equations 1 and 2 to evaluate the trajectory of each particle, subject to imposed forces given by Equations 3, 5, 8, 20, and 22, and torques given by Equations 4, 6, 13, 14, 21, 23 and 24.

D. Brownian stress calculation
One can similarly obtain from fluctuation-dissipation theorem an expression for the Brownian stress resulting from the pairwise interaction between particles i and j that averages over many realisations so that ⟨Σ B,L i,j ⊗ Σ B,L i,j ⟩ recovers the form of the hydrodynamic lubrication stress, but as described above this can similarly be shown to be equivalent to Σ B,L i,j = −F B,L i,j ⊗ r i,j .Since the pairwise Brownian force term contains the normal operator N i,j acting on the random vector θ i,j , one obtains a prefactor in the stress containing the dot product n i,j • θ i,j .This quantity will always approach zero when averaged over many realisations of θ i,j , so that the Brownian stress computed in this way averages to zero.Nonetheless, particle pairs do experience non-zero Brownian forces acting at all timesteps that will influence their trajectories so that the resulting contact and lubrication stresses will be altered by the presence of the Brownian forces.Below we describe a method that allows us to estimate the contribution of Brownian motion to the overall stress.
It is important to note here that our method, in which particle inertia is accounted for, is fundamentally different to other computational approaches, notably Stokesian Dynamics (SD) [11,13,28] in which the trajectories are evolved with a timestep longer than the inertial one.In the latter methods (see in particular Banchio and Brady [12]) the Brownian stress for the overall system is obtained as ), in practice using a midpoint scheme in which the positions and velocities of every particle are sampled at some increment of the overall timestep.Here R SU and R F U represent parts of the overall resistance matrix that couple, respectively, stresses to velocities and forces to velocities.Our method described above is based on the Langevin equation so that particle inertia is small but present, and force balance is not strictly achieved at each timestep.In order to obtain an estimate of the Brownian contribution to the stress, we deploy a structural method that exploits the anisotropy of the radial distribution function, using the approach described by Brady [29].The Brownian stress attributable to the pair i, j can be written as where p 1/1 (x j |x i ) is the probability density for finding a particle at x j given that there is a particle at x i , and n = N/V is the number density of particles in the suspension (where V and N are the system volume and particle number respectively).The integral is over the surface of contact S 2 of two touching particles.
To compute this function we sum for each particle the diadic product of its unit vector with each of its neighbours within a thin shell ∆ = 0.05a i , so that for a given configuration the Brownian contribution to the stress is [30]: The stress measured by this approach is not added to the hydrodynamic and contact stresses computed in our model, but is available to provide insight into the role of Brownian motion in setting the overall material response.

E. Additional simulation details
We simulate O(10 3 ) spherical particles of radius a and 1.4a (mixed approximately equally by volume) in a cubic periodic simulation box of length L. For each set of flow conditions we carried out between 10 and 800 realisations in order to obtain satisfactory ensemble averages.The principle particle properties (these set the length, mass and time scales) are the characteristic particle radius a [length], the particle density ρ [mass/length 3 ] (taken throughout to be equal to the fluid density so that the particles are neutrally buoyant), and the particle normal stiffness k n [mass/time 2 ] (this has a tangential counterpart k t ).With respect to these quantities, 1 time unit corresponds to the inverse frequency of a mass ρa 3 = 1 on a linear spring with stiffness k n = 1.The remaining material properties to be defined are the fluid viscosity η [mass/(length×time)] and the particle-particle friction coefficient µ [dimensionless], relevant for micron sized (and larger) particles.The thermal energy scale in the system is set by k b T .
The simulation box is deformed according to a specified ∇U ∞ .For instance, when the only nonzero element of ∇U ∞ is an off-diagonal (say γ), shearing is applied by tilting the triclinic box (at fixed volume) according to L xy (t) = L xy (t 0 ) + L γt.When the strain (γ = γt, with t the time for which the simulation has run) reaches 0.5 in this example, the system is remapped to a strain of -0.5.This has no effect on the particle-particle forces or on the stress, and is simply a numerical tool to permit unbounded shear deformation while preventing the domain from becoming elongated in one axis [31].Reported in the following is the relative viscosity of the suspension η r = Σ xy /η γ, with Σ xy the shear component of the stress tensor, γ the shear rate and η the fluid viscosity.

F. The timescales that appear in the simulation
The full list of dimensional parameters taken as inputs to the model is then a, L, t, ρ, k n , k b T , η and γ.Taking a/L ≪ 1 and γt ≫ 1, dimensional analysis dictates that we require three non-dimensional groups to fully characterise this system.In other words, a measured non-dimensional quantity e.g. the reduced viscosity η r = Σ xy /η γ, can be a function of at most three non-dimensional control parameters.This is in addition to non-dimensional inputs viz. the volume fraction ϕ and the friction coefficient µ.Central to our work will be the study of viscosities as a function of Peclet number, since this latter quantity will control the colloidal to granular crossover.It is desirous to choose the remaining two non-dimensional control parameters such that particles are effectively hard and non-inertial.To obtain an appropriate set of non-dimensional control parameters, we consider the following list of timescales present in the model (in which we only include dimensional elements for simplicity): The contact time τ C is a characteristic time spent by two particles in contact (assuming contacts are describable as linear springs), in the absence of other forces playing a role.It is obtained by solving the following equation of motion for the overlap δ between contacting particles: ρa 3 (d 2 δ/dt 2 ) = k n δ.The inertial relaxation time τ I is the characteristic time taken for the velocity of a particle to reach that of the background fluid in the absence of other forces.It is obtained by solving the following equation of motion for the velocity v of a particle: ρa 3 (dv/dt) = ηav.The Brownian time τ B is the characteristic time take for a particle to diffuse by a distance equal to its own radius under thermal motion in the absence of other forces.The convective timescale τ S is simply the inverse of the shear rate.To resolve each of these timescales accurately within the simulation we chose the numerical timestep to be substantially smaller than the smallest of the timescales listed above.The Peclet number (P e) described above is given by 6πτ B /τ S = 6πηa 3 γ/k b T , and we vary this quantity across a broad range from 0.01 to 100000, aiming to explore the colloidal to granular transition.
The contact time τ C should be chosen to be sufficiently small that overlaps between particles are orders of magnitude smaller than the particle radii, such that particles be considered hard spheres.To do this we ensure throughout that τ c is at least an order of magnitude smaller than the next smallest timescale.The role of particle inertia can be expressed via (i) a particle Reynolds number τ I /τ S = ρa 2 γ/η, and (ii) an inertia-diffusion ratio τ I /τ B = ρk b T /η 2 a. Below we explore how small each of these quantities need to be set in order to ensure inertia plays no significant role in the measured results.

III. RESULTS: INTERACTIONS AND DIFFUSION A. Two-particle simulations measuring the effective potential
To evaluate the net pairwise potential resulting from the particle-level forces described above, we carried out O(10 3 ) simulations of two particles with radii a in a cubic periodic box of length 4a (see snapshot in Figure 2(a) Inset) subject to all of the forces described above, and with U ∞ = 0. We calculate the radial distribution function g(r) with r = |r i,j | and averaged this across timesteps in the steady state and across all realisations (Figure 2(a)), then obtained the potential of mean force as U (r)/k b T = − ln(g(r)), Figure 2(b).The result confirms that there is no net potential acting between particles when they are not in contact (i.e. when r > 2a), so the lubrication and Brownian forces do not introduce an overall repulsion or attraction.When particles are in contact (r/(a i + a j ) < 1) there is a steep repulsive potential that, as expected, is related to the stiffness of our contacts defined above as U (r)/k b T = 0.5k n δ 2 i,j .The model thus approximates a suspension of colloidal hard spheres, in which the particle-particle interaction is zero and infinite for non-contacts and contacts respectively.

B. Mean square displacement
We next verify that our simulated particles follow statistically the anticipated trajectories by computing their mean squared displacement (MSD) under various conditions.An isolated particle with motion governed by the single body drag and Brownian forces described above is expected to follow a trajectory with a short-time ballistic part and a long-time diffusive part that leads to an overall MSD given by [32,33]: with m = (4/3)πρa 3 and γ = 6πηa.This expression gives ⟨x 2 ⟩ ∼ t 2 and ⟨x 2 ⟩ ∼ t at small and large times respectively.It can equivalently be written in terms of our characteristic timescales defined above as: Shown in Figure 2(c) are MSDs for a dilute sample with ϕ = 0.001 in which pairwise particle-particle interactions are absent.In terms of our model timescales, we set τ S = ∞ (i.e.no shear); τ C = 10 −3 ; τ I = 10 −1 ; and we vary τ B to explore the behaviour at different temperatures.We measure the elapsed time in units of τ I , so that the crossover from ballistic to diffusive behaviour begins in each case at t/τ I ∼ 1.As expected based on the expression above, increasing temperature (which decreases τ B ) while keeping all other variables constant simply shifts the MSD result vertically with ⟨x 2 ⟩ ∼ k b T .
We next calculate the MSD for a series of larger ϕ, with results shown in Figure 2(d)-(e).In all cases the particles follow a ballistic trajectory at short times that is roughly independent of ϕ.The longer time behaviour shows a decreasing diffusion coefficient (D = d/dt(⟨x 2 ⟩)) with increasing ϕ, a consequence of pairwise hydrodynamic and contact interactions resisting particle motion.For all volume fractions below jamming D approaches a constant at long time scales, confirming the presence of a diffusive regime.
In order for inertia to play a negligible role in our model, it is important for the diffusive timescale to be longer than the inertial relaxation one.In other words, the time taken for a particle velocity to relax to that of the background fluid should be much shorter than the time taken for the particle to diffuse by its own radius.To understand quantitatively how to achieve this, we measured D for varying τ I /τ B across a broad range of ϕ.The normalized long time diffusion coefficient (D(ϕ)/D 0 ) is shown in Figure 2(f), with D 0 = k b T /πηa(= a 2 /τ B ).Our result shows that when τ I /τ B is smaller than 0.17, D(ϕ)/D 0 becomes independent of temperature and follows a linearly decreasing trend.This suggests a criteria for the maximum value of τ I /τ B , which we check under shearing conditions in the following.

IV. RESULTS: RHEOLOGY
In the following we first describe the need for substantial ensemble averaging, especially when Brownian motion dominates, and we demonstrate the convergence of the measured rheology with the size of the sampling window.We next go on to expose the role of particle inertia in our model under shear, and establish the parameter range in which it can be assumed negligible.We then present rheology data showing η r as a function of P e, highlighting the breakdown of the individual contributions (hydrodynamic, contact and Brownian) and their variation with volume fraction.We finally demonstrate the role of particle contact friction and a short-ranged repulsive potential.

A. Averaging method
All of the rheology simulations described in the following were carried out with O( 103 ) particles, comprising an approximately equi-volume mixture of those with radius a and 1.4a.Given the comparatively small number of particles (compared to a real experimental system, for instance) and the random nature of the Brownian forces added to the system, the stress signals output by a single simulation are extremely noisy, especially at low P e. (The same is true for inertia-free simulations [34], though the error bars are rarely reported.)Thus the number of realisations that must be averaged over to obtain smooth data and reliable estimates of the true rheology increases as P e is reduced.
Shown in Figure 3(a) is the range of measured η r as a function of the number of steady state snapshots averaged over, for 3 different P e.At low P e one must sample the system ≈ 10 8 times to obtain a measurement of η r with standard deviation less than 10%, whereas for large P e 10 5 samples are sufficient.Importantly, the time taken to reach steady state also differs drastically with P e.For systems dominated by thermal fluctuation (i.e.low P e) the approach to steady state is set by the passage of Brownian time as opposed to the accumulated strain, with systems at ϕ = 0.54 and below taking 3-4 Brownian times to reach steady state at P e = 0.01.For larger ϕ this timescale is stretched rapidly, likely due to the proximity of glassy physics.At very large P e, meanwhile, steady states are reached for strains γt of 1-2 [35].

B. Brownian stress at zero shear rate
To obtain the Brownian contribution to the viscosity in the limit of zero shear rate, we apply the Green-Kubo method [36] by calculating the time autocorrelation function of the shear stress, taking as input the Brownian stress computed as described above, for unsheared simulations.The Brownian viscosity is written as: where Σ B xy is the shear component of the Brownian stress tensor.The stress correlation decreases exponentially with increasing ∆t so that the Brownian viscosity can be modelled as η B,GK (t) = η ϕ (1−e −∆t/τ ϕ ).As shown in Figure 3(b), the correlation time τ ϕ is short and weakly varying for ϕ < 0.5, so that η ϕ can be measured using readily accessible data for which ∆t/τ ϕ is large.For ϕ > 0.5, however, τ ϕ grows quickly and we estimate η ϕ by extrapolation.The rapid growth of the correlation time τ ϕ is likely indicative of a nearby glass transition, though we defer detailed analysis of this behaviour to future work.By this approach we obtain an estimate of the Brownian contribution to the viscosity at zero shear rate, which we discuss further in the following.
Given the large quantity of data required for obtaining smooth results, it is worth considering the scaling of the simulation run time with the system size.To estimate the scaling of the run time we carried out simulations with N = 10 1 , 10 2 , 10 3 , 10 4 particles with ϕ = 0.5, running a serial build of LAMMPS on one core for 10 7 timesteps.The result shown in Figure 3(c) confirms that our simulation has complexity O(N).

C. The role of inertia
The particle-particle contact timescale τ C is set sufficiently small that it does not compare to any other timescale in the system under any conditions, so that particles can always be considered to be hard.We verify this in Figure 5(b) by showing that the relationship between η r and P e measured under different values of τ C does not vary.It is, however, crucial that in varying P e one maintains acceptable values of τ I /τ S and τ I /τ B .To determine sufficiently small values of these two ratios so that inertia may be neglected, we carried out two sets of simulations.In the first we simulate shear flow with ϕ = 0.55 and k b T = 0 (so we don't need to consider the Brownian timescale τ B ), while varying the dimensionless shear rate τ I /τ S from 5 × 10 −3 to 10.To be in the limit in which inertia is negligible, we require a linear relation between the shear stress and the shear rate i.e. a Stokes flow.In other words, we are correctly simulating an inertia-free flow if η r is independent of τ I /τ S .From Figure 4(a) we can observe that this holds for τ I /τ S ⪅ 10 −1 .In what follows, we therefore ensure that this inequality holds for all parameter sets.Our result here is qualitatively consistent with prior simulations [17] and experiments [37,38], though the value of the Stokes number at the crossover is apparently highly sensitive to system details.
In the second we simulate shear flow with τ I /τ S < 0.01, ϕ = 0.45 and at a range of P e, exploring the relative importance of inertia by varying the timescale ratio τ I /τ B .This control parameter essentially sets the distance a particle will typically cover under ballistic motion.In order for inertia to be negligible in the model, we expect that this distance should be at least an order of magnitude smaller than the particle size, so that a typical Brownian kick to a particle does not lead it to collide with a distant neighbour.From our result in Figure 2(b) we find that a ballistic to diffusive crossover occurs at ⟨x 2 ⟩/a 2 = 0.01 for τ I /τ B = O(10 −2 ).Our shear simulations (Figure 4(b)) similarly show that η r is a function of τ I /τ B only when the latter quantity is > 0.01.Therefore, in what follows we carry out simulations with τ I /τ B < 0.01 and τ I /τ S < 0.1.

D. Flow curves
Our main rheology results are presented in Figure 5.We simulated a broad range of P e (10 −2 − 10 4 ), focussing on three different volume fractions ϕ (Figures 5(a)-(c)) and adhering to the constraints on τ I obtained above.To achieve this range of P e it was necessary to vary both the shear rate γ and the thermal energy k b T .We present in Table I a full list of the parameters used to generate the result in Figure 5(a).In each rheology figure we break the overall viscosity down into its contributions from hydrodynamic, contact and Brownian stresses.The stresses obtained by taking the outer product of the pairwise vectors and forces evaluated during the simulation run are the hydrodynamic one, the contact one, and the 'instantaneous' Brownian stress.The latter (not shown in Figure 5), as described earlier, averages to zero so does not lead to a viscosity contribution.The total stress (shown in black in Figure 5(a)-(c)) is therefore just the sum of the hydrodynamic and contact parts.
As a post-processing step we make an estimate of the effective Brownian stress (approximating the one that would be measured in a Stokesian Dynamics simulation), following the calculation based on structural anisotropy described earlier.This gives us the red lines in Figure 5(a)-(c).Interestingly the Brownian stress maps quite closely to the contact stress for low P e, indicating that the surge in contact stress observed in this range may be due to short-lived contacts induced by the Brownian kicks.Indeed the formulation of the Brownian stress is similar to that of the   contact stress, differing only in the presence of the contact overlap δ i,j appearing in the latter.
Overall we find that the predicted rheology corresponds well with canonical results, both in the experimental literature [39,40] and those obtained by Stokesian Dynamics simulation [41] and similar numerical methods [34].At all volume fractions there is a shear thinning region for P e < 1 that gives way to shear thickening beyond P e > 1, with the η r values at large P e tending towards those reported for non-Brownian suspensions under a very similar numerical framework [23].For ϕ = 0.45 and ϕ = 0.5 we observe a low P e plateau, whereas at ϕ = 0.55, η r apparently continues to increase with decreasing P e.The latter effect is perhaps an artefact of proximity to a glass transition, though we defer a more detailed study of this effect to future work due to the diverging timescales involved.The hydrodynamic stress increases weakly with increasing P e, whereas the contact stress qualitatively follows the overall stress in its shape.The increase in contact viscosity at high P e may be attributed to the onset of contact force chains as the system approaches the non-Browian limit and can be considered granular [42], while at low P e it is related to the Brownian forces as described above.
Shown in Figure 5(d)-(f) are slices through the three dimensional radial distribution function g(r i,j ), showing the flow-gradient (xy) plane at P e = 0.01, P e = 1 and P e = 10 4 .The general shape of the pairwise distributions is consistent with literature data [43], showing increased anisotropy with increasing P e and a sharpening of the peaks at a + a, a + 1.4a and 1.4a + 1.4a.

E. Viscosity variation with volume fraction
In order to understand better the limiting behaviour at small P e, we determine the behaviour of η r as a function of ϕ.To do so we first evaluate the Brownian contribution to the zero shear viscosity using the Green-Kubo method described above.To obtain an estimate of the full viscosity, we take the value of the hydrodynamic viscosity at the smallest (non-zero) measured P e, and add this to the Green-Kubo prediction of the Brownian stress (assuming the latter to be a good proxy for the contact stress, as was assumed by Brady [29] and is supported by our simulation data in Figure 5).Doing so at a range of ϕ, and comparing the result to the minimum η r measured at P e = 4.6 for each ϕ as well as the large P e limit, we obtain Figure 6(a).
In both the low and high P e limits, we find that η r , particularly at large ϕ, can be fit relatively well with a simple relation as η r ≈ (1 − ϕ/ϕ J ) −λ , with ϕ J (P e → 0) = 0.587 and ϕ J (P e = 10 4 ) = 0.642 (and λ ≈ 1.5, similar to Mari et al. [44]).At intermediate P e, η r is reduced relative to its value in the the non-Brownian limit, and the value of ϕ J is marginally increased.The large P e value of ϕ J will be highly sensitive to details of the particle-particle contact interaction, especially the presence of a static friction coefficient as we have reported elsewhere [23,45].In particular, for large friction coefficients the large P e value of ϕ J (usually denoted ϕ m ) will likely drop below the low P e value.In this scenario one expects flow curves near jamming to be diverging at both low and high P e, with finite η r at intermediate P e.We leave this complexity to be explored in future work, and in the following we examine the role of friction for a small range of ϕ.

F. Role of particle-particle friction and short-ranged repulsion
In the context of experimental work by Guy et al. [6], it is important to consider the role of particle friction at the colloidal-to-granular transition.Since granular particles are large, micron size objects they will likely have a static friction coefficient, which may constitute both sliding and rolling components [45,46].So far we have only considered a model system of frictionless particles.It is well-established that the presence of static sliding friction means that each particle-particle contact will constrain more than one degree of freedom of each particle, so that for large friction coefficients (in practice µ ⪆ 0.5) a rigid packing can be obtained with a per particle contact number of ≈ 4 (as opposed to 6 for frictionless spheres), with limiting volume fraction ϕ m ≈ 0.57.In Figure 6(b) we report rheology predictions from simulations of suspensions with a range of particle-particle friction coefficients µ (black data), demonstrating that the presence of friction leads to a dramatic increase in η r at large P e.This behaviour, and its sensitivity to ϕ demonstrated in Figure 6(c) (black data), is qualitatively consistent with the large literature on friction-driven shear thickening e.g.Mari et al. [44].Notably, η r at lower P e is unaffected by friction, suggesting that Brownian forces suppress the mobilisation of static friction for all µ, at least at ϕ = 0.5.In this respect the Brownian forces act analogously to a weak repulsive potential, inhibiting the formation of sustained particle contacts and rendering the suspension effectively frictionless even when µ > 0. This leads to a bulk viscosity with rate dependence qualitatively similar to that of shear thickening suspensions with load-activated friction describable by the canonical model of Wyart and Cates [47].Importantly, though, is it not clear that the shear thickening transition, when controlled by Brownian motion, is governed by a single stress scale.In particular, the range of P e over which the transition happens in Figure 6(b) (black data) is rather broad (occurring over 4-5 orders of magnitude in P e), especially when compared to Mari et al. [44] in which the transition takes at most 2 orders of magnitude in shear rate.To explore this we introduce a short ranged repulsive force defined by with κ = 0.01(a i + a j ).We show results of this model for A/k b T = 0, 1, 10 2 , 10 3 in Figure 6(b) and for several ϕ at µ = 0.5 in Figure 6(c).Introducing a sufficiently large repulsive force scale (in practice we required A/k b T ≈ 100) narrows the range of P e over which shear thickening occurs, and shifts the transition to larger P e.This result suggests not only an additive effect of Brownian and repulsive forces as reported by Mari et al. [34], but rather a qualitative change in the functionality of η r with P e when the onset of contacts is set by the magnitude of Brownian or repulsive forces.Examining the subtly in more detail is a promising area in which our model might be deployed.Thus with the introduction of particle-particle friction and a short ranged repulsive force we can control in our model the position and extent of shear thickening, providing a flexible starting point from which to make predictions of the rheology in more specific contexts.

V. CONCLUDING REMARKS
In conclusion, we have implemented a minimal numerical model for the rheology of dense suspensions that incorporates sufficient microscopic physics to predict the colloidal to granular crossover as a function of P e.The model is implemented in LAMMPS [21] so that its run time scales linearly with the number of particles.The Brownian component of our model differs from that in Stokesian Dynamics in that we resolve the fluctations at a much shorter, inertial timescale.The naively-calculated Brownian stress therefore averages to zero over realisations and instead we compute an estimation of the Brownian contribution to the stress based on the structural statistics measured from the simulation.This stress follows closely the contact stress that we measure directly from the pairwise forces and relative positions.The model predicts shear thinning at low P e, with a low P e plateau (in some cases) that increases with volume fraction.At larger P e a Brownian regime gives way to a contact dominated regime in which particle-particle interactions proliferate and friction (if present) becomes important.In this latter regime shear thickening is observed even for zero particle friction, though its extent increases with increasing friction coefficient.We finally introduced into our model a short range repulsive force, a crucial prerequisite for shear thickening in the paradigmatic model of non-Brownian suspensions [44].This keeps particles separated and inhibits the contact contribution to the stress, thus broadening the intermediate P e viscosity plateau (as observed by Cwalina and Wagner [48]) or equivalently shifting the value of P e at which particle contacts become important.
We have focussed in this article on steady, simple shear rheology.Broadening the work to inhomogeneous conditions (such as those described by Gillissen and Ness [49]) and to dynamic simple shear to measure the frequency-dependent (and indeed amplitude dependent [50]) response are promising lines of future research that will provide additional scope for constitutive model development and for validation against experimental data.
In future we anticipate deploying our code in mixed systems, in which one population of particles are Brownian and another are non-Brownian [48].This is motivated by numerous real world examples such as geophysical flows and many scenarios in chemical engineering and manufacturing.In such systems the small, Brownian particles (in some industries these are referred to as superplasticizers) will simultaneously contribute a Brownian stress but improve the efficiency of packing, so that their overall effect on the rheology is non-trivial and likely to be non-monotonic and P e-dependent.Mapping out this complexity as functions of the small and large particle sizes and their relative numbers requires a tractable numerical model, and will likely rely on the implementation of more advanced neighbour listing algorithms such as those by Shire et al. [51].Extending this further to systems with continuous, broad size distributions remains on open challenge [52].

FIG. 1 .
FIG. 1. Snapshot of simulation box, and schematics of the leading pairwise force terms present in the model.In all cases the leading component of the pairwise force acts along the (positive or negative) direction of the unit vector ni,j pointing from the centre of particle i to the centre of particle j.(a) Snapshot of the simulation box showing particles of radius a (red) and 1.4a (blue).Particles are mixed in approximately equal volume.(b) Contact force, with force acting along ni,j; shown from left to right are [i] sketch of contacting particles i and j, with the overlap δi,j shown in green; [ii] the unit vector ni,j pointing from the centre of particle i to j; [iii] repulsive contact F C forces acting along the positive and negative directions of ni,j.(c) Pairwise hydrodynamic lubrication force, with force set by the component of the relative particle velocity acting along ni,j; Shown from left to right are [i] sketch showing the velocity vectors U of neighbouring particles i and j; [ii] the relative velocity Uj −Ui breaks down into tangential (yellow) and normal (red) components, with the latter pointing along ni,j; [iii] lubrication forces act along the positive and negative directions of ni,j, proportional to the normal part of the relative velocity.(d) Pairwise Brownian lubrication force with a random pairwise vector θi,j projected onto ni,j.Shown from left to right are [i] neighbouring particles with centre-to-centre unit vector ni,j; [ii] random vectors θi,j drawn from a Gaussian distribution (green) are projected onto ni,j by the tangential (yellow) and normal (red) operators; [iii] Brownian forces act along the positive and negative directions of ni,j, proportional to the normal part of the random vector.

2 ⟩/a 2 t
a i + a j ) U(r)/k b T r − (a i + a j ) −ln(g(r)) 0.5k n δ 2 i, j /k b T ⟨x

FIG. 2 .
FIG. 2. Evaluating the potential of mean force and the diffusion properties that arise from the particle-level forces described above, in the absence of shear flow.(a) The radial distribution function g(r) (with r = |ri,j|) computed from a two particle simulation [Inset: snapshot of simulation]; (b) Potential of mean force U (r), showing measured result (points) and the input particle stiffness (solid line); (c)-(d) Mean squared displacement as a function of elapsed time for (c) three values of the timescale ratio τI /τB at ϕ = 0.001; (d) three values of ϕ at τI /τB = 1.7;(e) Diffusion coefficient as a function of elapsed time for a range of ϕ at τI /τB = 1.7.The solid line in (c)-(e) represents the predictions of Equation 31; (f) Long time diffusion coefficient at a broad range of ϕ and τI /τB.

FIG. 3 .PeFIG. 4 .
FIG. 3. Computing the suspension viscosity ηr under sheared and non-sheared conditions, and the scaling of computational run time with system size.Shown in (a) is the convergence of the measured ηr as a function of the number of snapshots averaged over, for P e = 0.01 (red), 1 (green) and 100.(blue) The noisy stress signal when Brownian motion dominates necessitates large numbers of realisations.In (b) is ηB,GK measured via the Green-Kubo relation taking the autocorrelation of the Brownian shear stress as input, plotted as a function of the correlation time; (c) Simulation run time versus number of particles for ϕ = 0.5, P e = 1, when running a serial compilation of LAMMPS on a single processor.We show data for a short simulation comprising 10 7 timesteps.

FIG. 5 .
FIG. 5. Rheology and microstructure of dense suspensions at the transition from Brownian to non-Brownian flow.Shown in (a)-(c) are the suspension viscosity ηr as functions of P e, for frictionless particles with (a) ϕ = 0.45; (b) ϕ = 0.5 (shown also are results for two additional values of kn); (c) ϕ = 0.55, showing the contributions from contacts, hydrodynamics and Brownian forces.In (d) is the total viscosity as a function of P e and ϕ.In (e)-(g) are slices through the three-dimensional radial distribution function g(ri,j) showing the flow-gradient (xy) plane under steady state simple shearing conditions for ϕ = 0.5 and (e) P e = 0.01; (f) P e = 1; (g) P e = 10000.

FIG. 6 .
FIG. 6.The viscosity variation with volume fraction and particle-particle contact friction.(a) Variation of ηr with volume fraction ϕ at three P e, showing fits to ηr = (1 − ϕ/ϕJ ) −λ ; (b) ηr as a function of P e for several particle-particle friction coefficients µ and repulsive force magnitudes A/k b T at a volume fraction of ϕ = 0.5; (c) ηr as a function of P e for several ϕ and A/k b T , with µ = 0.5.The colour legend in (b) refers also to (c).

TABLE I .
The parameters used to generate each data point in Figure5(a).