Modelling inelastic granular media using Dynamical Density Functional Theory

We construct a new mesoscopic model for granular media using Dynamical Density Functional Theory (DDFT). The model includes both a collision operator to incorporate inelasticity and the Helmholtz free energy functional to account for external potentials, interparticle interactions and volume exclusion. We use statistical data from event-driven microscopic simulations to determine the parameters not given analytically by the closure relations used to derive the DDFT. We numerically demonstrate the crucial effects of each term in the DDFT, and the importance of including an accurately parametrised pair correlation function.


Introduction
Granular media is ubiquitous in industrial and natural processes [4,41], but very difficult to accurately model in large systems. Microscopic models [13,7] generally produce accurate results, but are usually restricted to simple systems, a relatively small number of particles, or short simulation times. This is due to computational cost scaling poorly with the number of particles in the system. Models which approximate the media as a macroscopic continuum [30,46] are not inhibited by the number of particles in the system. However, continuum models need to be supplied with constitutive equations for bulk properties, such as the particle stress tensor. Such constitutive equations, in the dilute flow regime, are usually obtained by invoking the kinetic-collisional theory [30]. The disadvantages of such constitutive equations lie mainly in what is believed to be its inability to treat the meso-scale: issues such as cluster formation (and breakage), for instance, have been discussed at large [27], and attempts to solve those issues have been proposed (e.g., by adjusting the Navier-Stokes equations [27,35]).
Recently, approaches using Dynamical Density Functional Theory (DDFT) have produced promising results in the field of modelling complex fluids [32,3,31] at a mesoscopic level. Derived from a particlebased model, DDFTs are continuum models that utilise the well-studied Helmholtz free energy functional [22] and can include interparticle and external potentials, volume exclusion [44], hydrodynamic interactions [40,20] and multiple species [18]. However, current models do not account for inelastic (or indeed elastic) collisional dynamics, which have crucial dissipative effects in granular media.
We introduce a DDFT adapted for granular media in Section 2 which incorporates particle collisions at the mesoscopic level, and present some numerical results in Section 3, which show the potential of DDFT for modelling granular media. There are many areas for further investigation, which we discuss at the end of this paper.
2 Derivation of the model 2

.1 Microscopic Dynamics
A set of N particles, each of mass m, in d dimensions can be modelled via Langevin [15] or Newton equations of motion: for positions r N = (r 1 , r 2 , ..., r N ) and momenta p N = (p 1 , p 2 , ..., p N ), the dynamics are given by where Here V ext (r i , t) represents any external potential, for example gravity. Pairwise interactions are modelled by an interparticle potential V 2 (r i , r j , t) for i, j = 1, ..., N , and we analogously include higher order interparticle potentials such as V 3 (r i , r j , r k ) for i, j, k = 1, ..., N . The second term on the right hand side of the equation for momentum in eq. (1) represents external frictional effects, for example if the particles are moving in a bath, where γ > 0 determines the strength of the effect. The final term a(t) is a Brownian motion term [10] from thermal fluctuations in the bath, with strength (mk B T γ) 1/2 determined by a fluctuation-dissipation theorem [16], where k B is Boltzmann's constant and T is the bath temperature. Granular media particles are usually assumed to be unaffected by thermal fluctuations, however for modelling purposes this term is sometimes included as a thermostat [33]. For point-like particles, eq. (1) fully describes the dynamics of the system. In this paper, we assume that particles are spherical with diameter σ. We include the effects of collisions in the dynamics of particles by restricting their movement to the hard sphere domain: When two particles that are moving toward one another come into contact, we must then instantaneously change their momenta to avoid particle overlap. We assume that collisions are binary and instantaneous; i.e. a collision between the i th and j th particles occurs at time t if r i (t) − r j (t) = σ, and (p i (t) − p j (t)) · (r i (t) − r j (t)) < 0. To resolve the collision we apply a collision rule. A standard collision rule which maps pre-collisional velocities p in i , p in j post-collisional velocities p out i , p out j is given by [6]: where ν i,j = (r i − r j )/ r i − r j and α ∈ (0, 1] is the coefficient of restitution. If α = 1, the collisions are perfectly elastic and no energy is dissipated, and the component of velocity in the direction of the collision is reflected. If α < 1, energy is lost via a reduction of the velocity component in the direction of the collision. We note that linear and angular momentum are conserved by this collision rule: However, kinetic energy is not conserved when α < 1: For simplicity, in this derivation we do not investigate the effect of angular momentum on the DDFT model [14] when collisional effects are included. In principle, eq. (1), on the domain eq. (3) with collision rule eq. (4) can accurately model a system of N particles. Particle based methods can produce very accurate results, but for large N , or when the system is very dense, reaching the desired simulation time is generally infeasible, or the simulation becomes too memory-intensive. Furthermore, when α < 1 a system of particles obeying the above microscopic dynamics can experience inelastic collapse, where an infinite number of collisions can occur in finite time, effectively jamming simulations [34].

N body equations for rigid particles
If collisions are neglected, associated with eq. (1) is the Kramer's equation (or in absence of thermal fluctuations, the Liouville equation) [42], a partial differential equation (PDE) which models the dynamics of the N -particle distribution function f (N ) (r N , p N , t), the probability of finding N particles with positions r N and momenta p N at time t: When constructing eq. (8) for deterministic dynamics, the microscopic dynamics are assumed to be smooth. However, particles which undergo instantaneous collisions have discontinuities in their velocity profile, and so their trajectories are not smooth. In [47], the weak formulation of the Liouville equation is derived for a system of elastic spherical particles that obey linear dynamics, using distribution theory. Without any additional assumptions, careful dissection of the phase space in the weak formulation leads to an additional collisional term: where for k ≤ N , andn is the outward unit normal of ∂D N (r 1 , t) with H(r N ) the Hausdorff measure on ∂D N (r 1 , t). When considering the associated Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy [5] (in weak form), the additional term integrates to the well-known Boltzmann collision operator for α = 1. We note that this is an alternative method to derivations where the collision term is constructed by additional assumptions on an interaction force at the level of the BBGKY hierarchy [24]. It is a topic of future work to see how long range potentials, friction and inelasticity affect this derivation. In this work we make an assumption that is popular in the literature; interactions between particles can be split into a term which describes 'soft' interactions, and a term for 'hard' interactions, e.g. a collision operator. This assumption is heuristically validated by the derivation in [47] where the collision operator occurs due to geometric properties of the phase space. The collision operator derived in [47] does not include soft interparticle interactions, so the potential term V 2 must be incorporated separately.
Equation (8) (or an analogous weak formulation including eq. (9)) replaces a system of 2N differential equations with a single PDE. However, the (spatial-momentum) dimension of f N (r N , p N , t) is 2dN . If we were to simulate this system on a discretized domain with M points in each direction, we would require M 2dN points for simulation, which quickly becomes computationally intractable. However, it is known [11] that the N -particle distribution function is a functional of the one-body position density: We include the characteristic function to stress that the phase space of the system does not allow particles to overlap, ensuring that we integrate over the hard-sphere domain while keeping the first position variable free. In contrast, in derivations where particles interact purely via soft potentials, the integral in position is over R d(N −1) . Heuristically, when integrating eq. (8) the inclusion of this characteristic function leads to collisional terms in the BBGKY hierarchy. Rigorously, the weak formulation admits the Boltzmann collision operator [47].
To arrive at an equation to model ρ(r, t) we first define the n-reduced phase space particle distribution function by where r N −n = (r n+1 , ..., r N ), r n = (r 1 , r 2 , ..., r n ), and similar for p N −n , p n . To ease notation in the derivation we write r = r 1 , p = p 1 . By integrating eq. (8) with respect to r N −1 and p N −1 , and appealing to the symmetry of arguments in f (N ) , we arrive at the one-body Kramer's equation, the first equation in the BBGKY hierarchy: where L coll (f (2) ) incorporates binary collisions via a collision operator. We consider the following inelastic collision operator [8]: where p i is the pre-collisional velocity associated to p i , determined using eq. (4). Long range interactions are included in where, for example, v 2 (r, r 2 ) relates to the two-body potential V 2 (r i , r j ), where the prefactor of 1 2 has been absorbed by a symmetry argument.

DDFT Derivations
As the BBGKY equations are hierarchical, they do not constitute a small enough closed set of equations for efficient simulation. We therefore must truncate the hierarchy, and introduce additional assumptions to close the remaining set of equations.
Equation (13) involves integrals with higher order distribution functions which must be approximated. We consider moments of eq. (13), and close this system by approximating higher order moments and distributions using lower order counterparts. We first approximate the many-body interactions (eq. (15)) in the nonequilibrium fluid by those of an equilibrium fluid with the same one body density profile [3]: where δF ex [ρ(r)]/δρ(r) is the functional derivative of the excess part of the Helmholtz free energy functional. We also assume that higher order distributions are uncorrelated in velocity: for 2 ≤ k ≤ N . We note that eq. (4) implies there is correlation in velocity, however we expect these correlations to be short range and therefore dominated by correlations in position. We note that, upon truncation of the hierarchy and approximation of higher order distributions in terms of low order counterparts, information from the Liouville equation (the N th equation in the BBGKY hierarchy), in particular volume exclusion effects, are lost. It is therefore necessary to include a term which approximates volume exclusion from the Liouville equation. It is popular to include a pairwise interaction 'potential' which forbids overlap: The Helmholtz free energy functional can be generalised to include volume exclusion due to hard particles via a suitable modification of F ex . In one dimension we use the exact functional for volume exclusion derived by Percus [45], while Fundamental Measure Theory (FMT) is used to accurately approximate volume exclusion for spheres [43] for d > 1. However, neither the 'potential' in (18) nor FMT directly include the effect of collisions in the system, which must be included using a collision operator.
The collision operator is also an alternative way of including volume exclusion effects. In [33], in one dimension, volume exclusion effects are incorporated by approximating the correlation function g (2) in the RET collision operator using an analytic form: where η(x) is the local packing fraction. As the packing fraction approaches 1, the value of g given by eq. (19) blows up and, in an analogous way to the Percus free energy in eq. (29), causes volume exclusion in the model. However, the numerics in the present work show that this approximation is not accurate for dynamics with inelastic collisions, where, in particular, the value of the correlation function at contact increases for low local densities, rather than decreasing as in the elastic case. In Section 3 we provide an example which shows that when the correlation function is approximated by experimental data and a volume exclusion free energy term is absent, the local density can exceed physical limits (see fig. 5).
Returning to the derivation of a continuum model, we define ρ, v and E as the number density, the local average velocity and the granular temperature of the system respectively: We finally assume that the one particle distribution function can be approximated by a local equilibrium Maxwell-Boltzmann distribution [21]: It is well known that the local equilibrium of a granular fluid is in fact not Maxwellian [9], and other approximations are a topic of current research [17]. However, the assumption eq. (21) allows us to write the second moment of f (1) (r, p, t) as a product of the density ρ(r, t) and local average velocity v(r, t). We expect other functional forms of the local equilibrium can be implemented in the same manner, but to introduce the model we use eq. (21). The continuity equation is then derived by integrating eq. (13) with respect to p, under the assumptions stated: By multiplying eq. (13) by p, then integrating with respect to p, standard calculus results lead to the momentum equation [3], which now includes the granular temperature, and the first moment of the collision operator: and also includes the Helmholtz free energy functional: where Λ is the (irrelevant) thermal de Broglie wavelength. Finally, when considering the third moment by multiplying by p ⊗ p then integrating with respect to momentum, by using eq. (23) and eq. (22), we arrive at an equation describing the evolution of the granular temperature: Equations (23) and (25) include centred moments of the collision operator: where the argument of L coll has changed to account for the assumption eq. (17). For eq. (14), by applying eq. (21), we can write eqs. (26) and (27) in terms of Gaussians and error functions. The exact forms used in simulations are given in appendix A. Given the centred moments of L coll (f (1) , g (2) ), and the correlations g (k) , the set of equations eqs. (22), (23) and (25) then constitute a closed model for granular media, incorporating volume exclusion due to hard particles, external and inter-particle potentials, and (in)elastic collisions.
There are some important differences between the DDFT model here and existing DDFTs in the literature. Firstly, we include moments of the collision operator eqs. (26) and (27), which must be included to incorporate dissipative effects due to inelastic collisions. Our numerical experiments (see fig. 3a) show that the collision terms do affect the dynamics, and that the effect is pivotal when α < 1. We also include an additional moment eq. (25) of eq. (13), as the effects of the collision term are evident in the granular temperature of the system; eq. (4) reduces the variance of particle velocities when α < 1, so we expect it to have a dissipative effect on E(r, t). In particular under our assumptions, in one dimension, Thus inelasticity has a small effect on eq. (23), but can be incorporated by including an additional moment. Although DDFT derivations involving collision terms [33] and temperature gradients [48] have been studied, it is clear that for granular media both terms play important roles. When comparing to results in kinetic theory, the addition of the free energy term allows us to include effects both from interparticle interactions and volume exclusion, by considering the interactions at the particle level.
We also note the importance in the choice of initial condition; the initial density, velocity and granular temperature must satisfy the physical restrictions of the PDE; in one dimension this corresponds to not exceeding the packing fraction limit ρ v = 1, and that the density must be non-negative ρ(r) ≥ 0 for all r.
In one dimension, we can incorporate volume exclusion exactly using the Percus free energy [38]: We note that in one dimension for constant ρ, as ρ → 1 σ , we approach the maximum density and so δFex[ρ] δρ → ∞, i.e. the chemical potential of the system blows up. This can be seen as a constraint on adding particles to the system, and implicitly stops the volume of the system increasing beyond the physical limit; eq. (29) is defined so that it is limited by the close packing value. Thes effects are not present when solely including a collision operator with arbitrary g, such as those obtained from microscopic simulations.

Parametrisation of g (2) (σ)
To use eqs. (22), (23) and (25) derived in Section 2 we need accurate approximations for the functions g (k) , in particular the pair correlation function g (2) , if we assume that particle interactions are pairwise. Analytical approaches to finding g (2) are generally restricted to simple systems [26]. One can construct an additional coupled 2d-dimensional PDE for g (2) via moment closure schemes [25], but this increases the dimensionality of the problem and so introduces a prohibitively larger computational cost when d > 1.
It is known [40] that higher order correlations equilibriate much faster than the density. This result validates the adiabatic approximation; correlations can be approximated by their local equilibrium values. Many empirical forms for g (2) at equilibrium (or quasi-equilibrium) are constructed and used in the literature, although it is reasonable to expect that g (2) depends on properties of the particles and dynamics. In particular, inelasticity leads to the effect of particle streaming, which should be visible in the correlation function (see fig. 2).
We therefore empirically construct g (2) for the system of interest. By applying statistical methods on synthetic data generated by extensive particle simulations with small N , we parametrise g (2) without simulating the entire system, avoiding excessive computational cost. We can then incorporate it in the model eqs. (22), (23) and (25). Analogously, parametrisations could be performed with experimental data.
We present an example of this methodology using a system of 100 deterministic hard rods (d = 1) on a 2π-periodic domain in the absence of external or interparticle potentials, with γ = 2 and σ = 2π/25, 000. In the absence of a thermostat, the trajectories of individual particles can then be solved analytically. Therefore, instead of a numerical method involving a timestep for microscopic simulation, we predict collision times of particles, then sort and schedule and process these events before advancing simulation. This methodology is at the centre of Event Driven Particle Dynamics (EDPD), which was first developed as early as the 1950s [2], but is still a modern topic of research.
A naive EDPD algorithm will predict future collisions between all particles after each collision has been processed, producing an algorithm with computational complexity O(N 2 ) per collision. Our simulations use the cell method [1] to reduce the computational cost of event prediction to O(N ). We also update the position of each particle asynchronously [28], which reduces the cost of advancing simulation to O(1). In addition, data structures such as binary search trees [39] and bounded increasing priority queues [37] can be implemented to decrease the computational cost of event sorting and scheduling to O(1) per collision. Combining methods for prediction and scheduling gives dynamics with a cost of O(1) per collision. Software packages such as DynamO [7] implement all of these methods, but are currently limited to systems with friction coefficient γ = 0. Finally, methods to efficiently parallelise EDPD algorithms have also been constructed [23].
In addition, several numerical methods are available to avoid inelastic collapse in EDPD. A review of these methods is available in [29], and in our simulations, we implement the TC model, which renders collisions elastic when a particle undergoes a collision a small time t c after a previous collision. Together with polydispersity, the TC method stops sharp peaks from forming in the radial correlation function, but also accurately approximates the dynamics of the system. We give two examples of situations that undergo inelastic collapse in one dimension in fig. 1. When α < 1, an infinite number of collisions occurs in finite time, so the dynamics are 'jammed'. The TC model allows the particle to 'vibrate' when the collisions become very frequent. We note that, by design, volume exclusion is incorporated in the algorithm; collisions are predicted and processed so that particles do not overlap. Therefore, unlike in continuum modelling, we do not require any additional potential in these dynamics to include volume exclusion effects.
Using EDPD to construct 5000 samples of each system with a range of values of α and solid volume fractions ρ v , we construct a parametrisation of g (2) (σ). Examples of systems with α = 0.9 and α = 0.5 are given in fig. 2. For d = 1 in the absence of friction, external and interparticle potentials, the radial correlation function blows up when particles are monodisperse and inelastic; at equilibrium all particles will be moving in contact. We therefore include a variance in the diameter of particles σ v = 0.1σ. Equations (22), (23) and (25) can be adapted to take into account poly-dispersity, but we expect the effect to be negligible in this case and so we ignore it in the DDFT. In each sample, the initial velocities of the particles are normally distributed with mean 0 and variance 1, and positions are drawn from a uniform distribution in the domain. We evolve the system until 99.9% of the energy of the particles has dissipated due to inelasticity and friction, then construct a near equilibrium parametrisation of g (2) ( r 1 − r 2 ). We note that when ρ v is small and α is close to 1, few collisions happen before the effect of friction causes particles to lose all their energy, and when ρ v = 1, the system is fully dense, so the value of the correlation function is independent of α.
In fig. 3a, we display g (2) ( r 1 − r 2 ) for different densities, as well as g at time t = 0. We note that for low densities and α = 0.5 the radial correlation function has several peaks. This is evidence of particle streaming, where inelastic collisions cause particles to move at the same velocity, near to one another.   Figure 2: Samples that can be used to construct g (2) (σ) at different points in time. The red lines are trajectories of the centres of mass of individual particles. We note that the domain is periodic, so trajectories can disappear and reappear at the top and bottom of the y axis. Left: α = 0.9, right: α = 0.5. In these examples the initial conditions are the same, but display characteristic differences in their dynamics, in particular when α = 0.5 particle streaming is more evident. Under the assumption that g (2) = g (2) ( r 1 − r 2 ), only the value of g (2) at r 1 − r 2 = σ is included in the collision terms. In fig. 3b we use the results of particle simulations to parametrise g (2) (σ) for different α and ρ v , using cubic smoothing spline interpolation [12]. We omit the point when σ v = 1 from the interpolation to improve the curve fit for values used in the DDFT simulation.

DDFT simulation
For a continuum approach we simulate eqs. (22), (23) and (25). We use pseudospectral code provided in [36], which is available at [19]. We consider a periodic domain [0, 100] with 100 computational points. Using more than 100 points has little effect on the result of the simulation. To match microscopic simulations, we set γ = 2 and σ = 2π/25000, so that there are 120, 000 particles in the system (a fully packed domain would hold 400, 000 particles in this case). We include the Percus free energy functional [38] to incorporate volume exclusion in the dynamics.
Before continuing, we note that the system can be made dimensionless by considering the following scaling: where L is a length scale (in our simulations we use the domain length, but the particle diameter could also be considered), T is a time scale and M is a mass scale (the mass of a particle). Furthermore, the granular temperature is dimensionless. The scaling of variables and parameters are constructed by considering the the microscopic dynamics and the distributions f (N ) as in [33]. The initial conditions are given by where N c is a normalisation constant, and ρ v = 0.3 is the total solid volume fraction. The initial conditions are chosen such that areas of higher density will move toward one another and 'collide'. Figure 4 displays results of model eqs. (22), (23) and (25) at different times t when α = 0.5. The results show that every term is necessary for accurate dynamics: If the Percus term is neglected the density is sometimes overestimated as particle volume exclusion of hard particles is not incorporated. If the collision term is included but g (2) (σ) = 1 (i.e. the uncorrelated case) the inelastic effects are not noticeable, and Figure 4: Results from the DDFT simulation for α = 0.5. Each simulation has the same initial condition (black, dashed) at time 0. The black line neglects the collision operator and the free energy term. The blue line includes the free energy term but not the collision operator. The green line includes both terms, with g (2) (σ) = 1, and the red line includes both terms and uses g (2) (σ) determined by particle simulations, shown in fig. 3b. the high density areas are reflected upon 'collision'. When all terms are present with g (2) constructed from particle simulations we see that the two higher density areas coalesce, an intuitive result for inelastic dynamics. To show the importance of including a volume exclusion free energy term in the DDFT, we consider an example with modified initial conditions: we set N = 175, 000, and v 0 (r) = 26 sin 2πr 100 .
The results in fig. 5 show that the system reached an unphysical density in finite time if the Percus free energy is not included. Finally, we perform the same dynamics for different α with initial conditions eq. (31). The results in fig. 6 show that particles coalesce more for smaller α. In fig. 7 we provide the density near equilibrium for different coefficients of restitution. We note that in this example, the long time behaviour of the density is similar for all coefficients of restitution. This is because the effect of the collision operator is small when the local average velocity is small, so in this example where energy is not added into the system using any other external potentials, the effect of friction determines the dynamics for long times.

Conclusions and future work
We have constructed a new model for granular media, which can incorporate inelastic collisions using classical collision operators, and interparticle interactions using DDFT methods. We have presented a simple example which displays the importance of each term in the model, but the model can also be used for systems in 2D or 3D with more complicated dynamics such as adhesion between particles. Our results show that our methodology is successful; small scale, inexpensive particle dynamics can be used to fine-tune parameters in the mesoscopic model, such as the radial correlation function. Figure 5: Results from the DDFT simulation, using g (2) (σ) determined by particle simulations for α = 1, and the same initial condition (black, dashed). The blue line gives the result when the free energy term is neglected, and the green is the same simulation with the Percus free energy term included. Any values of ρ which are above the physical limit of ρ = 1 are coloured red. In this example we used 600 computational points, to ensure that the volume exclusion effects are numerically stable. Figure 6: Results from the DDFT simulation, using g (2) (σ) determined by particle simulations for different α, using the same initial condition (black, dashed). Much of the current research on DDFT for complex fluids can be adapted to the system of equations in this work, including extension to poly-disperse or multi-species systems [18], and inclusion of more complicated drag forces due to interactions between the particles and fluid in the system using a hydrodynamic interaction tensor [20]. Further work on fundamental derivations in the style of [47] will be beneficial to construct collision operators from more complicated dynamics or particles. Furthermore, inclusion of a more accurate local (non-Maxwellian) equilibrium approximation for granular media could improve the model.
The synthetic data presented is an example of how modern computational and data-scientific methods can be applied to fine-tune parameters in continuum models; we parametrise using statistics from particle simulations. For systems with more complicated interactions between particles we will need to use state of the art particle simulation methods, but our modelling approach avoids the computational bottleneck caused by simulating large numbers of particles. Figure 8: Results from the DDFT simulation, using g (2) (σ) given by eq. (19), using the same initial condition (black, dashed), for different coefficients of restitution α.
B DDFT using an analytic approximation of g In [33], an analytic approximation of the radial correlation function is considered, which is independent of α. Although our numeric investigation of the correlation function for different α disagrees with the analytic approximation, for comparison in fig. 8, we provide results using this approximation of g, with the same initial configurations considered in Section 3, in fig. 6. The results show that, under this approximation, the introduction of inelasticity plays a much smaller role in the dynamics. This is in contrast to the large effects seen in the microscopic simulations.