Simulating dense, rate-independent suspension rheology using LAMMPS

Dense suspensions are widespread in nature, manufacturing and process engineering. Particle-based simulations have proven to be an invaluable complement to experimental rheological characterisation, serving as a virtual rheometer that enables rapid exploration of parameter space and detailed scrutiny of microscopic dynamics. To maximise the utility of such simulations, it can be advantageous to exploit pre-existing, well-optimised, well-documented codes. Here we provide a simple description of how to use LAMMPS to study the rheology of dense, granular suspensions.


Introduction
Suspensions of micron-sized particles with solid and liquid mixed in roughly equal proportion present intriguing flow properties that challenge physicists and engineers of all kinds [1].A useful starting point for characterising their rheology is to understand the rate-independent behaviour (relevant for solid particles of ≈ 10-1000 µm size), demonstrated experimentally by Boyer et al. [2] and later reviewed by many others (e.g.[3]).A number of particle-based simulations reproduce the rate-independent rheology (e.g.[4,5]), providing (i) corroboration of the experimental result; (ii) a source of particle-resolved data inaccessible experimentally; and (iii) a platform for examining systematically more complex microphysics, for instance particle inertia [6] (accepting that the forces as written below apply strictly only in the viscous limit), stress-induced friction [7] and adhesion [8].Many simulation techniques and codes are available, and we do not review these here.In what follows, we describe how to generate numerical rheology data consistent with the rateindependent result using LAMMPS [9,10].
The simulation technique shares many details with traditional molecular dynamics, and the reader is expected to be familiar with basic concepts including contact detection, neighbour listing, timestepping and so on.Our specific model is more commonly labelled as a 'discrete element method' B Christopher Ness chris.ness@ed.ac.uk 1 School of Engineering, University of Edinburgh, Edinburgh EH9 3JL, UK owing to its similarity to approaches used in granular physics, specifically the absence of thermal forces and the inclusion of particle-particle friction (both appropriate for the size range mentioned above).The basic approach for obtaining a numerical rheology measurement is to (i) initialise the system with a packing of non-overlapping spherical particles of desired size distribution at a desired volume fraction φ; (ii) evaluate the trajectory of each particle i by numerically solving "F = ma" in the presence of a prescribed background fluid velocity gradient ∇u ∞ and a set of pairwise hydrodynamic and contact interactions.(We assume φ 0.4 throughout, otherwise a more detailed account of the hydrodynamics is required.)When desired, a bulk stress tensor ˚is calculated from the interaction forces and particle positions, thus generating rheology data, viz. the stress ˚as a function of deformation rate E (with E ≡ 1 2 (∇u ∞ +∇u ∞T )) and volume fraction φ.

Dimensional analysis for rate-independent suspensions
It is instructive to consider first a dimensional analysis.This introduces the quantities (and their units) that we will define when setting up a simulation.We consider a dense, granular (so we omit k b T from the following) suspension of stiff spheres under a flow with imposed deformation rate.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 η f [mass/(length×time)] and the particleparticle friction coefficient μ [dimensionless], relevant for micron-sized (and larger) particles.The relevant macroscopic quantities are the size of the simulation box L [length] and the volume fraction φ [dimensionless] therein.The background fluid flow is characterised by a velocity field u ∞ [length/time] and its gradient (a tensor) ∇u ∞ [1/time] (that we specify, and take to be spatially uniform), the time t for which the flow was applied, and a stress tensor [mass/(time 2 ×length)] (that we measure).We write a scalar velocity gradient as γ (≡ ∂u x /∂ y) and a scalar stress as xy (the x y component of ˚).A list (others are possible) of non-dimensional parameters necessary to fully define a given suspension under given flow conditions is then: Setting (i)-(iii) to be 1 ensures, respectively, bulk conditions, stiff particles and no particle inertia.Under these conditions, and assuming μ is constant (i.e.particle friction is Coulombic) and we shear to steady state ( γ t → ∞), we have simply that η r = η r (φ), hence the label 'rate-independent'.This is consistent with the result of Boyer et al. [2] and will be the focus of our example below.

Particle-level forces and shearing
Rate-independent rheology is obtained by subjecting particles to three types of force and torque: Stokes drag, pairwise lubrication and pairwise contact.The full form of these is reported by several authors [5,[11][12][13] and need not be repeated here.Instead, we describe the forces in simplified terms.
The Stokes drag (Fig. 1a) on particle i (radius a i ) is proportional to the difference between its velocity u i and the fluid streaming velocity at its centre u ∞ (x i ): This force is essentially what induces flow in the simulation, causing particles to conform to the streaming velocity set by u ∞ .Similarly, a torque acts to cause the particles to rotate with angular velocity ( i ), set by Neighbouring particles i and j with centre-to-centre vector r i, j (Fig. 1(b)[i]) experi-ence lubrication forces (see [14,15]) dependent on the gap h between them and their relative velocity (Fig. 1(b)[ii]).The leading term of the force on particle i (assuming it has equal radius to particle j) scales with 1/h and the normal component of the pairwise velocity difference: while the leading term of the torque is: The lubrication forces oppose relative motion between particle pairs.They are prevented from diverging at contact by setting a lower limit on the allowed value of h (typically O(10 −3 a i )) and are cut off at an outer distance of 0.05a i .At the volume fractions of interest, all particles have multiple neighbours within this range, so that setting a larger cut-off does not significantly affect the result.Overlapping particle pairs i and j (Fig. with the torque given by c i, j = a i k t (r i, j /|r i, j | × ξ ). ( The friction coefficient μ sets an upper bound on ξ through |ξ | ≤ μk n δ/k t .The stress contribution from drag forces is proportional to E. The α, β component of the stress due to lubrication and contact is found, respectively, by summing (F l,α i, j r )/2 and F c,α i, j r β i, j over all pairs.The forces are summed on each particle, and the trajectories are then updated according to Newtonian dynamics, using a numerical scheme with timestep chosen to be small compared to ρa 3 /k n and ρa 2 /η f .In LAMMPS, the simulation box deforms according to the specified ∇u ∞ .For instance, when the only nonzero element of ∇u ∞ is an off-diagonal, shearing is applied by tilting the triclinic box (at fixed volume) according to L xy (t) = L xy (t 0 ) + L y γ t, Fig. 1(d).When the strain (γ = γ t) 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.

LAMMPS inputs and outputs
The above physics are implemented in LAMMPS (Largescale Atomic/Molecular Massively Parallel Simulator) [9], a classical molecular dynamics code written in C++.The LAMMPS documentation should be referred to at all times [10].A skeletal set of instructions for downloading and compiling LAMMPS and running the scripts below is provided at Ref [16].Our strategy is to separate the generation and shearing of suspensions into two distinct simulations.The first defines dimensionless numbers (i) and (v) above, producing non-overlapping particle packings in a cuboidal, periodic domain of fixed φ; the second defines (ii), (iii), (iv) and (vii) then applies the deformation allowing us to measure (vi).This decomposition allows one to build a library of configurations at different φ that can be reused many times for different deformation protocols.LAMMPS takes as its input a text file containing a list of commands and their arguments.In the following, we don't describe each of these in detail, but instead provide minimal scripts, indicating where the physics above enters.Each of the commands is fully described in the LAMMPS documentation [10] Creating particle packings.-Forthe time being, it isn't necessary to specify the 'full' physics described above, nor do we need to output the stresses.Rather, we need just enough detail to create assemblies of non-overlapping particles.We therefore omit lubrication forces at this stage for simplicity.We first generate particles (of two types, each having a different radius) with random coordinates in a box of set dimensions.Their overlaps generate contact forces (following Eq. 4) that lead to motion; damping against a stationary background fluid (i.e.u ∞ = 0 in Eq. 1) extracts energy until the system comes to rest.The properties of the configuration of particles produced by this script don't really matter: we are not trying to sample an 'equilibrium' configuration (this is not relevant for granular systems) but simply create a packing with no (or minimal) overlaps.Shown in the in.createpanel is an example input script ('#' indicates comments) to generate a suspension with φ = 0.5.
Highlighted in red from top to bottom are (with units as stated earlier): (i) the size of the cubic simulation box L = 14.6381; (ii) the numbers and radii of particles of types 1 and 2: N 1 = N 2 = 100, 2a 1 = 2 and 2a 2 = 2.8, thus setting the volume fraction as φ = 4/3π(N 1 a 3  1 + N 2 a 3 2 )/L 3 = 0.5; (iii) the particle density ρ; (iv) the particle stiffness k n = 10000 (with k t = 7k n /10); (v) the particle-particle friction coefficient μ = 0.1; (vi) the timestep; (vii) the fluid viscosity η f = 0.1 (the number we put in the script is 6πη f a); (viii) the number of timesteps to run.Highlighted in blue are the seeds used to generate the initial particle positions.New realisations can be generated by rerunning the simulation with different numbers here.This script produces data.file,containing a snapshot of the system after the final timestep, to be read in by future scripts.The file contains a list of the particle ID, diameter, density, coordinates (x, y, z) and velocity components.It also produces a dump file (create.dump) that can be visualised using, for example, Ovito [17].
Shearing particle packings.-Thesecond script (in.run, see panel below) takes data.fileas an input and applies a deformation to the sample.The key inputs to this script are the parameters related to the particle-particle interaction (k n , k t , μ and η f ) and the shear rate, set by specifying the components of ∇u ∞ .Through these, we set the values of dimensionless control parameters (ii), (iii), (iv) and (vii) listed above.# SPECIFY THE OUTPUTS thermo_style custom time c_str [1] c_str [2] c_str [3] c_str [4] c_str [5] c_str [6] thermo 10000 dump id all custom 10000 run.dumpid x y z radius log run.log# SPECIFY THE TIMESTEP, THE INTEGRATION SCHEME AND RUN timestep 0.0001 fix 1 all nve/sphere fix 2 all deform 1 xy erate 0.001 remap v run 30000000 The remaining content of the input script is concerned with specifying the bulk stress calculation and requesting it as an output, necessary for obtaining (vi).
Highlighted from top to bottom (units as before) are the fluid viscosity η f = 0.1, the particle stiffness k n = 10000 (k t = 7000) and friction coefficient μ = 0.1 and the shear rate γ = 0.001.Using the command fix deform we have specified just the x y component of ∇u ∞ , γ as defined above.The other components are 0 by default, so this leads to a simple shear with flow in x and gradient in y i.e. u ∞ = ( γ y, 0, 0).The script runs at this γ for 30000000 timesteps, each of duration 0.0001.Thus, the total shear strain is γ t = 0.001 × 30000000 × 0.0001 = 3.Note that the fix viscous command is not required here because the drag force is applied within the lubrication pair style.
Outputs from shearing simulation.-There are two different types of output produced by LAMMPS during a simulation run, log files and dump files.Log files are typically described as containing thermodynamic data, but for our purposes we can interpret this as bulk suspension or derived properties, usually the components of the stress tensor ˚, but also e.g. the average particle contact number.As specified in in.run above, the log file (run.log)contains the accumulated simulation time t, followed by the 6 unique components of the stress tensor in order (x x, yy, zz, x y, xz, yz), output every 10000 timesteps (specified by the thermo command): run.logTime c_str [1] c_str [2] c_str [3] c_str [4] c_str [5] c_str [6] 0 0.00073391 -0.0018116 -0.00028190 -0.001589 0.00018953 -0.00058983 1 0.00063973 0.0014631 0.00086064 -8.3279e-05 -5.7537e-05 -0.00020223 2 0.00040622 0.00045020 0.00066278 -0.00015159 9.5433e-05 -0.00012521 3 0.00029818 0.00032181 9.1296e-05 -0.00044916 6.2575e-05 -0.00012202 . . .These stress components have units [mass/(time 2 ×length)].In order to express in dimensionless units (the "reduced viscosity" as it is conventionally represented) for comparison to experimental data, one must divide by η f γ (and off-diagonals must be multiplied by −1).Shown in Fig. 1(e) is a plot of the reduced viscosity as a function of the accumulated strain.The reduced viscosity is in this instance taken simply as the x y component of the stress ( xy , the 4th column of the stress outputs in the log file) divided by η f γ .whereas the strain is the accumulated time t multiplied by γ .From the content of the log file, one might also compute, e.g. the viscous number η f γ /P (with P the mean of the diagonal components of ( x x , yy , zz ), the normal stress differences x x − yy , yy − zz and so on.Dump files contain particle-level information (positions, velocities, radii) or contact-level information (forces, relative positions), usually output at fixed intervals.In the script above this is specified (by the dump command) to be every 10000 timesteps.This example script produces run.dump, which lists the particle IDs, positions and radii at specified time intervals.This file might be used for post-processing, for instance, to compute structural properties, to follow particle trajectories or to be read directly into various visualisation packages (Fig. 1(d) was generated using Ovito [17], for instance).By repeating this pair of simulations at a range of φ (achieved by changing the value of L in in.create) and taking time (and ensemble) averages of the components of ˚, one can reproduce the rate-independent rheology of Boyer et al. [2], as demonstrated in Fig. 1(f)-(g).Moreover, one may relax the conditions specified in the dimensional analysis to explore additional physics: (i) introducing particle inertia by increasing the value of ρ γ a 2 /η f ; (ii) introducing particle softness by increasing the value of γ ρa 3 /k n , and so on.

Closing remarks
We have provided a brief description of how to use the molecular dynamics code LAMMPS to generate dense, granular suspension rheology data.Examples of the use of this approach to study the physics of dense suspensions can be found in Refs [18][19][20][21][22][23].Moreover, the flexibility of the code allows one to simulate more complex geometries (such as extrusion [24]), both within LAMMPS [25,26] and within derivative codes such as LIGGGHTS [27].Suspensions of broad polydispersity might be simulated using more complex pairwise hydrodynamics such as those described by [28], while rate-dependence can be introduced via the introduc- 1(c)[i]) experience repulsive contact forces dependent upon the scalar overlap δ (Fig. 1(c)[ii]) and the tangential displacement accumulated over the duration of the contact ξ (Fig. 1(c)[iii]):

Fig. 1
Fig.1Particle-level physics in a simplified dense suspension.a A particle (radius a i , position x i , velocity u i ) in a fluid with streaming velocity u ∞ ; b pairwise lubrication interaction showing [i] particle velocities u i , u j , centre-to-centre vector r i j and surface separation h; [ii] relative velocity (u j − u i ) and its components normal (u j − u i ) n and tangential (u j − u i ) t to r i j ; c contact force showing [i] particles approaching; [ii] particles entering contact with overlap δ and initiation of sliding with ξ = 0; [iii] contact with overlap δ and accumulated, nonzero, tangential PARTICLE CONFIGURATION read_data data.file# SPECIFY THE PARTICLE-PARTICLE INTERACTION pair_style hybrid/overlay granular lubricate/bmpoly 0.1 1 1 0.001 0.05 1 0 pair_coeff * * granular hooke 10000 0 tangential linear_history 7000 0 0.1 pair_coeff * * lubricate/bmpoly # DO THE STRESS CALC compute str all pressure NULL pair .