Non-Dissipative Structural Evolutions in Granular Materials

The structure of the contact network in granular assemblies can evolve due to either dissipative mechanisms such as sliding at contact points, or non-dissipative mechanisms through the phenomenon of contact gain and loss. Being associated with negligible deformations, non-dissipative mechanisms is actually active even in the small strain range of ∼ 10−3, especially in the case of densely packed assemblies. Hence, from a constitutive modelling point of view, it is crucial to be able to estimate such non-dissipative evolutions since both elastic and plastic properties of granular assemblies highly depend on contact network characteristics. The current study proposes an analytical scheme that allows us to estimate the non-dissipative contact gain/loss regime in terms of directional changes in the average contact force. The probability distribution of contact forces is used to compute the number of lost contact for each direction. Similarly, the number of newly formed contacts is estimated by considering the probability distribution of the gap between neighbouring particles. Based on the directional contact gain/loss computed, the changes in coordination number and fabric anisotropy can be found which, together with statistical treatments of Love-Weber stress expression, form a complete system of equations describing the evolution of other controlling microvariables. Finally, the results of the calculations have been compared with DEM simulations which verify the accuracy of the proposed scheme.


Introduction
The internal structure of jammed granular materials can change under external loading according to various mechanisms that can be classified into three broad categories: • Non-dissipative-reversible (e): which refers to the socalled elastic deformation at interparticle contact points.Such processes conserve energy, contribute to overall strain, and are recoverable upon load reversal.
• Dissipative-irreversible (p): which includes the plastic deformation associated with sliding at contact points.Such processes dissipate energy, contribute to overall strain, and are, by nature, irrecoverable upon load reversal.
• Non-dissipative-irreversible (δ): which includes the changes in contact network topology due to contact gain and loss events.These events do not contribute to overall strain, nor do they dissipate energy.However, such structural alterations are irreversible and the previous structure cannot be recovered upon load reversal.
Unlike the common elastic/plastic dichotomy of mechanisms in continuum mechanics, the presence of nondissipative/irreversible events in granular media has given rise to many peculiarities in granular media.Being associated with no deformation, such events are capable of altering the contact structure within minute strain ranges which has prompted the rapid adaptability of the internal structure with respect to the external loading observed in initially dense granular assemblies [1,2].Moreover, being non-dissipative by nature, such contact loss/gain mechanisms can theoretically coexist with purely elastic deformations of the domain which leads to a non-reversible elasticity and the concept of locked energy in closed stress loading paths [3].This is in stark contrast with common beliefs that elastic response of granular materials should, by definition, preclude any evolution of the contact network [4,5].
From a constitutive modelling point of view, it is crucial to be able to calculate the evolution of internal structure mainly because almost all micromechanical models (including models for elasticity of granular materials) heavily depends on the structure of contacts often expressed in terms of coordination number and fabric anisotropy.However, there have been very few studies in the literature on the analytical calculation of the contact gain/loss regimes and their contribution to ensuing structural evolutions [6].
As such, the current study proposes an analytical scheme to compute the change in fabric anisotropy and coordination number due to non-dissipative/irreversible mechanisms.Starting from a generic contact structure, the directional change in average contacts forces is calculated which, together with the directional distributions of contact forces and interparticle gaps, are used to estimate contact loss and gain regimes through a probabilistic approach.

Decomposing Fabric Evolution
The second-order fabric tensor describing the directional distribution of contact normals is classically defined as follows: where N c is the total number of contacts, n i is the contact normal vector, and a c is fabric anisotropy.The directional probability density function of normal contacts can also be expressed in terms of fabric anisotropy as: where the angle θ is measured from the major principal direction of F i j .The incremental change in fabric tensor can be decomposed into contributions from lost (a l c ), gained (a g c ), and persisting contacts (a p c ), with a δ c = a g c − a l c being the contributions of δ mechanism, and a p c of the p mechanism.The outcome of such decomposition in a typical 2D, biaxial DEM test on an initially dense granular assembly is given in Fig. 1.Similar results have been reported in [2].It can be seen that δ mechanism dominates the overall structural evolution during the initial stages of deviatoric loading.The dissipative mechanism (p) becomes significant only after a threshold for coordination number is reached, as argued in [1].

Analytical Modelling
This section deals with calculating the directional distribution of lost and gained contacts through probabilistic treatments.Interested readers are referred to [7] for more details.We start by calculating the change in average normal contact force for each direction in terms of contact and force anisotropies.Next, assuming the probability distribution function (PDF) for magnitude of normal contact forces to be known for each direction, the number of lost contacts is calculated by integrating over contact forces with magnitudes smaller than the change in average normal force.A similar procedure is followed for gained contact with the integration being over the interparticle gaps smaller than the average deformation for each direction.It should be mentioned that, for the sake of simplicity, all the calculations are presented for 2D condition with no principal stress rotation.The calculation procedure, however, can be easily extended to be applicable to more general cases.

Average Normal Contact Force
The average normal contact force for direction θ is defined which can change as either the sum of the forces or the number of contacts along θ changes.However, in order for the variables to be comparable at the initial state, here we are only interested in the change in average normal force due to the former effect.As such, the change in average normal force in the absence of contact gain/loss, d fn * (θ), can be calculated as: Adopting the notation and anisotropy parameters introduced in [8], the directional distribution of average normal contact force can be approximated by fn (θ) = fn (1 + a n cos 2θ) with fn and a n being the average normal force and normal force anisotropy, respectively.The change in average normal force can be written as: Moreover, recalling that N c (θ) = N ct 2π (1 + a c cos 2θ) from fabric anisotropy relation, the relative change in number of contact along a given direction can be calculated as follows: with dN ct /N ct = dz/z, z being coordination number defined as average number of contacts per particle.By inserting Eqs. ( 4) and ( 5) into (3) and truncating second and higher orders of anisotropy parameters, the change in average normal force for a given direction can be expressed in terms of the change in prominent micro-variables, i.e.
On the other hand, based on expressions in [8], the stress increments can be written in terms of microvariables: As a further simplification, it is temporarily assumed that the contribution of da t (anisotropy of tangential forces) is negligible, and as such, Eq. ( 7) can be substituted into (6) to replace the microvariables with boundary stresses, i.e.The value of Δθ can be calculated by setting Eq. ( 8) to zero.Thus,

Lost and Gained Contacts
The PDF of the magnitude of contact forces along a direction θ is assumed to be given by P n ( f n / fn (θ)) whose form does not change with direction.A typical form of such PDF is given in Fig. 2 where the PDF is shown to differ for all the contacts and contacts along minor stress direction.However, the lost contacts are mainly oriented close to minor stress direction and the dependency of PDF on the direction can be arguably ignored.The number of lost contacts along a direction θ can now be calculated by integrating over the normal forces whose magnitude is smaller than d fn * (θ): It should be noticed that the integration is performed only over the minute part of force PDF, within which the distribution can be accurately replaced by an straight line.The total number of lost contacts can be calculated by integrating Eq.(10) over θ: The contribution of such a contact loss regime to fabric anisotropy change can also be calculated as follows: Turning to gained contacts, the number of newly formed contacts along θ is calculated by integration over particle pairs with interparticle gaps smaller than the average elastic deformation of contacts along direction θ.Such an average deformation is approximated by d fn * (θ)/k n with k n being normal contact stiffness.Compared to the procedure for calculating lost contacts, the normal force PDF is replaced here by PDF of interparticle gaps, P ξ (δ p /2r), with δ p being face-to-face distance between neighbouring particles, and r the average radius of particles.Figure .3 shows the form of P ξ obtained from the 2D biaxial DEM simulation.Similar to the lost contacts calculation, the directional variation of P ξ has been ignored.The number of gained contact along direction θ can now be calculated as follows: The parameter β nc determines the ratio of non-contacting neighbours 1 to the total number of contacts, which can be expressed in terms of coordination number assuming that the total number of neighbours (contacting and non-contacting) remains constant [7].The contribution of gained contacts to coordination number (dz g ) and fabric anisotropy (da g c ) changes can also be calculated similar to lost contacts, with the directional integration performed over Ω + .

System of Equations
The above analyses involve five microvariables describing the statistics of contact and force network; z, a c , fn , a n , and a t .Based on the decomposition given in Fig. 1, we can safely assume that microstructure of granular material changes solely due to non-dissipative δ mechanisms, and hence, the net changes in fabric anisotropy and coordination number calculated from contact gain and loss can be assumed to represent the total changes: dz = dz g − dz l and da c = da g c − da l c .Moreover, given the stress increments, Eqs.(7) provide us with two more equations.The last equation for closing the system of equations comes from the observation that the parameter a n is strictly proportional to the deviatoric stress ratio q/p at the initial stages, with the proportionality coefficient remaining the same irrespective of the loading path [1,7].The five above-mentioned relationships can now be solved in terms of stress increments to find the evolution of the five microvariables.

Model Verification
The model's predictions of the microvariables evolution have been compared with DEM simulation of 2D biaxial tests on an initially dense assembly.The simulation parameters are: number of particles 20, 000, particle size distribution D max /D min = 1.25, interparticle friction μ = 0.5, initial void ratio e 0 = 0.182, initial coordination number z 0 = 4.06, and contact stiffness k/p = 10 4 which has been chosen to be the same for normal and tangential directions.The simulation results have been used to find the PDF's of force and interparticle gaps which enter the calculation procedure as an input.Fig. 4 shows the comparison between model predictions and DEM results where a satisfactory agreement is observed especially for the initial part of the loading.As anticipated, the simulation results diverge from the predictions as the deviatoric load becomes large enough, which indicates the onset for prevalence of dissipative p mechanisms.

Conclusion
The role of non-dissipative/irreversible mechanisms such as contact gain and loss events in altering the internal structure of granular media has been investigated.It has been shown, through decomposition of fabric anisotropy, that such mechanisms are the main cause behind the rapid adaptation of the granular microstructure with respect to external stresses, especially at the initial stages of deviatoric loading.The irrecoverable fabric realignment associated with non-dissipative/irreversible mechanisms is found to coexist with elastic mechanisms at strain levels as small as 10 −2 , which compromises the possibility of deriving elastic deformational characteristics of granular materials from a energy potential.The study also puts forward an analytical procedure that relates the directional distribution of contact gain and loss regimes to the probability distribution functions of contact forces and interparticle gaps.The model has been validated by comparing the prediction of prominent microvariables to 2D DEM simulations.

Figure 1 .
Figure 1.Decomposition of fabric anisotropy into contribution from different mechanisms.

Figure 2 .
Figure 2. PDF of normal contact force magnitude for all the contacts, and contacts along the minor principal stress direction.

Figure 3 .
Figure 3. PDF of interparticle gaps for a 2D assembly.