From discrete particles to continuum fields near a boundary

An expression for the stress tensor near an external boundary of a discrete mechanical system is derived explicitly in terms of the constituents' degrees of freedom and interaction forces. Starting point is the exact and general coarse graining formulation presented by Goldhirsch in [I. Goldhirsch, Gran. Mat., 12(3):239-252, 2010], which is consistent with the continuum equations everywhere but does not account for boundaries. Our extension accounts for the boundary interaction forces in a self-consistent way and thus allows the construction of continuous stress fields that obey the macroscopic conservation laws even within one coarse-graining width of the boundary. The resolution and shape of the coarse-graining function used in the formulation can be chosen freely, such that both microscopic and macroscopic effects can be studied. The method does not require temporal averaging and thus can be used to investigate time-dependent flows as well as static and steady situations. Finally, the fore-mentioned continuous field can be used to define 'fuzzy' (highly rough) boundaries. Two discrete particle method simulations are presented in which the novel boundary treatment is exemplified, including a chute flow over a base with roughness greater than a particle diameter.


Introduction
The main topic of this paper is the issue of coarse-graining, near a boundary. We consider the bulk method described by Isaac Goldhirsch [1], and extend it to account for boundary forces due to the presence of a wall or base. The Goldhirsch special edition of Granular Matter is an appropriate place to present some of the ideas that we have developed in this area.
Continuum fields often need to be constructed from from discrete particle data. In molecular dynamics [2] and granular systems [3,4], these discrete data are the positions, velocities and forces of each atom or particle. In contrast, in the case of smooth particle hydrodynamics [5], the continuum system itself is approximated by a discrete set of fluid parcels. In all these methods, a crucially important issue is how to compute the continuum fields in the most appropriate way. Several techniques have been developed to calculate the continuum fields, see [10] and references therein. Particularly the stress tensor is of interest: the techniques include the Irvin-Kirkwood's approach [6] or the method of planes [7]. Here, we use the coarse-graining approach (CG) as first described in [8].
The CG method [8,1] has several advantages over other methods, including: i) the fields automatically satisfy the conservation equations of continuum mechanics; ii) it is not assumed that the particles are rigid or spherical, and iii) the results are valid for single particles (no averaging over ensembles of particles is required). The only assumptions are: each particle pair has a single point of contact, the contact area can be replaced by a contact point, and collisions are not instantaneous.
In Sect. 2, we use the derivation of [1] to extend the CG method to account for the presence of a boundary. Explicit expressions for the resulting continuum fields are derived. In Sect. 2.5, an alternative stress definition is proposed extend-ing the stress field into the boundary region. In Sect. 3, the approach is tested with two DPM simulations, and in Sect. 4 we draw conclusions.

Assumptions and notation
We are interested in deriving macroscopic fields, such as density, velocity and the stress tensor from averages of microscopic variables such as the positions, velocities and forces of the constituents. Averaging will be done such that the continuum fields, by construction, satisfy conservation laws. Vectorial and tensorial components are denoted by Greek letters in order to distinguish them from the Latin particleindices i, j. Bold vector notation will be used when appropriate. We will follow the derivation of [1], but extend it by introducing two types of particle: N flowing particles {1, 2, . . . , N} and K boundary particles {N + 1, . . ., N + K}.
Each particle i has mass m i , center of mass position r iα and velocity v iα . The force f iα acting on particle i is a combination of the sum of the interaction force f i jα with another particle j, the interaction force f ikα with a boundary particle k, and a body force b iα (e.g., gravity), The interaction forces are binary and anti-symmetric such that action equals reaction, f i jα = − f jiα , i, j ≤ N. We assume that each particle pair (i, j), i ≤ N, j ≤ N + K has, at most, a single contact point, c i jα , at which the contact forces act. The positions of the boundary particles are fixed, as if they had infinite mass. The trajectories of the flowing particles are governed by Newton's second law and if tangential forces and torques are present, rotations follow from the angular form of Newton's law.
In the following sections, we commence from Ref. [1] to derive definitions of the continuum fields. To be precise, a body force density is introduced to account for body forces, and to incorporate boundary effects an interaction force density (IFD) is introduced. While the idea of an IFD is more generally applicable (e.g., for mixtures), it is employed here to account for the presence of a boundary.

Coarse graining
From statistical mechanics, the microscopic mass density of the flow at a point r α at time t is defined by where δ (r r r) is the Dirac delta function. We use the following definition of the macroscopic density, i.e., we have replaced the Dirac delta function by an integrable 'coarse-graining' function W whose integral over the domain is unity.

Mass balance
The coarse-grained momentum density is defined by Hence, the macroscopic velocity field V α (r r r,t) is defined as the ratio of momentum and density fields, V α (r r r,t) = p α (r r r,t)/ρ(r r r,t). It is straightforward to confirm that ρ α and p α satisfy the continuity equation (c.f. [1,8]

Momentum balance
Subsequently, we will consider the momentum conservation equation with the aim of establishing the macroscopic stress field, σ αβ . As we want to describe boundary stresses as well as internal stresses, the boundary interaction force density (IFD), t α , has been included, as well as the body force density, b α , which are not present in the original derivation, [1]. The desired momentum balance equations take the form, To determine the stress it is required to compute the temporal derivative of (4), where f iα = m i dv iα /dt is the total force on particle i. Using (1), the first term in (7) can be expanded as with the abbreviation W i = W (r r r − r r r i ). The first term, which represents the bulk particle interactions, satisfies since f i jα = − f jiα and because the dummy summation indices can be interchanged. It follows from (9) that Substituting (10) into (8) yields Next, A α is rewritten using Leibnitz's rule to obtain a formula for the stress tensor. The following identity holds for any continuously differentiable coarse-graining function W where r i jα = r iα − r jα is the vector from r jα to r iα . Substituting identities (12) into (11) yields In Ref. [1], it is shown that the second term in (7) can be expressed as where v ′ iα is the fluctuation velocity of particle i, given by Substituting (13) and (14) into momentum balance (6) yields where the kinetic and bulk contact contributions to the stress tensor are defined as Here, the underlined terms in (16) are not in the original derivation presented in Ref. [1] and account for the presence of the boundary.
Expression (16) can be satisfied by defining the last term on the right hand side as the IFD. This however has the disadvantage that the boundary IFD is located around the center of mass of the flowing particles. The more natural physical location of the boundary IFD would be at the interface between the flowing and boundary particles.
Therefore, we move the IFD to the contact points, c ikα , between flowing and boundary particles: similar to (12), where W ik = W (r r r − c c c ik ) and a ikα = r iα − c ikα . Substituting (18) into the last term in (16) we obtain Thus, substituting (19) into (16), we define the stress by where the contribution to the stress from the contacts between flow and boundary particles is and the IFD is Equations (20) differ from the standard result of [1] by an extra term, σ w αβ , that accounts for the additional stress created by the interaction of the boundary with the flow. The definition (21) gives the boundary IFD applied by the flowing particles; i.e., it has been constructed such that in the limit w → 0, the IFD acts at the contact points between boundary and flow.
Note, that this framework is general and can be used to compute more than boundary IFDs. For example, one can obtain the drag between two different species of interacting particles by replacing the flowing and boundary particles with the particles of the two species in the definition of the continuum fields. By placing an IFD at the contact points, the IFDs of both species are exactly antisymmetric and thus disappear in the momentum continuity equation of the combined system. In mixture theory, e.g. [11], such interaction terms appear in the governing equations for the individual constituents and are called interaction body forces. These interaction body forces are an exact analog to the IFDs. Therefore, our approach can interpreted as treating the system as a mixture of boundary and flow particles and the IFD is the interaction body force between different species of particle.
Further, we note that the integral of the stress in (17) and (20) over the domain Ω satisfies the virial definition of mechanical stress,

Extending the stress profile into a base or wall
In contrast to the previous subsection, an alternative stress definition is presented here, where the IFD and the stress are combined into a single tensor. Similar to (12) and (18), the following identity holds, since the coarse-graining function W satisfies W (|r r r| → ∞) = 0. Substituting (23) into (16) we can obtain an alternative solution with zero IFD, t ′ α = 0, where the stress is given by σ αβ This stress definition is not identical to the one in (20) and (17). It eliminates the IFD term entirely and provides a natural extension of the stress into the boundary. However, the extended stress does contain contributions from both internal and external forces, and the spatial integral of the stress components has to be extended to infinity. In singular special cases this can lead to artificial results. Another disadvantage of Eq. (24) is its difficult interpretation due to the long-ranging integral. One could see it as the stress inside a 'virtual/fake' wall-material on which the body-force is not acting (equivalent to foam with zero mass-density). However, this is far fetched and not realistic, so that we rather stick to the formulation in Sect. 2.4. It is also possible to extend the stress tensor to the boundary by other means, such as mirroring the stress at the boundary, or using a one-sided coarse-graining function. This is not discussed further since the first method requires a definition of the exact location of the boundary, while the second method can introduce a spatial shift in the stress field due to spatial inhomogeneities.

Contact model
For illustrational purposes, we simulate a granular system with contact interaction forces. The statistical method, however, is based only on the assumptions in §2 and therefore can be applied more generally. We use a viscoelastic force model with sliding friction as described in detail in [9,4]. The parameters of the system are nondimensionalised such that the flow particle diameters are d = 1, their mass m = 1, and the magnitude of gravity g = 1. The normal spring and damping constants are k n = 2 ·10 5 and γ n = 50, respectively; thus, the collision time is t c = 0.005 d/g and the coefficient of restitution is ε = 0.88. The tangential spring and damping constants are k t = 2/7k n and γ t = γ n , such that the frequency of normal and tangential contact oscillation, and the normal and tangential dissipation are equal. The microscopic friction coefficient is set to µ p = 0.5. We integrate the resulting force and torque relations in time using the Velocity-Verlet algorithm with a time step ∆t = t c /50.
We take the coarse-graining function to be a Gaussian of width, or variance, w. Other coarse-graining functions are allowed, but the Gaussian has the advantage that it produces smooth fields and the required integrals can be performed exactly.

Quasi-static example in two dimensions
In order to visualise definitions (20) and (21), we firstly consider a two-dimensional configuration consisting of five fixed boundary particles and five flowing bulk particles, with gravity in the z-direction, see Fig. 1. The flow is relaxed until the flowing particles are static; hence, the only contribution to the stress is due to the enduring contacts. To visualise the spatial distribution of the IFD and stress, the norms |t t t| = t 2 α and |σ σ σ | = max |x x x|=1 |σ σ σ x x x| (the maximum absolute eigenvalue), are displayed in Fig. 1. A very small coarse graining width, w = d/8, is chosen to make the spatial averaging visible: the IFD, Eq. (21), centers around the contact points between flowing and static particles, r ikα , while the stress, Eqs. (20) and (17), is distributed along the contact lines, r iα r jα and r iα c ikα .

Three-dimensional steady chute flow
Secondly, we consider a three-dimensional simulation of a steady uniform granular chute flow, see Fig. 2 and Ref. [12].   chute contains 1000 flowing particles, which are initially randomly distributed. The simulation is computed until a steady state is reached. A screen shot of the steady-state system is given in Fig. 2 After averaging in x, y and t directions, this yields Integrating over (z, ∞), we obtain Thus, setting α = z in (27), the extended stress component σ ′ zz can be obtained without computing the semi-infinite line integral. The depth profile for the downward normal stress σ zz is shown in Fig. 3. Since the rough boundary is not at a fixed height, the stress gradually decreases at the bottom due to the decreasing number of bulk particles near the base. Due to the coarse graining, the stress tensor has a gradient even in the case of a flat wall, but the gradient disappears as w → 0. Using the extended stress definition, the bed and surface locations can be defined as the line where the downwards normal stress σ ′ zz vanishes and where it reaches its maximum value (to within 2%), see Fig. 3. Additionally, since the flow is steady and uniform, (6) yields the so called lithostatic balance, σ ′ zz = − ∞ z g z ρ dr z , which is satisfied with good accuracy.

Conclusions
We have derived explicit expressions for the stress tensor and the interaction force density (IFD) near an external boun-dary of a discrete mechanical system. These expressions were obtained by coarse-graining the microscopic equations and therefore exactly satisfy the governing balance laws of mass (5) and momentum (6). A boundary IFD was computed using the contact points between the flow and the basal particles. Our results can be extended to other IFDs, for example, the drag between two different species of particles. The power of our extension to Goldhirsch's method has been demonstrated by computing stress profiles for a chute flow over a fuzzy boundary. It avoids the problems inherent in other methods and gives the expected linear lithostatic profile all the way to the base.
The present formulation for boundary interaction forces allows us to draw the analogy to electrostatics, where the divergence of the electric field (analogous to the divergence of stress) is compensated by a charge-density source like our interaction force density (21). The analogy can also be made to mixture theory where, by treating the system as a mixture of boundary and flow particles, the IFD is then interpreted as the interaction body force between the two species.