A one parameter damageable contact law for DEM, with application to frictional-cohesive granular materials

A novel type of damageable cohesive law is presented for paired particle interactions in a DEM granular arrangement. It is designed in the spirit of a mixed breakage criterion for solid cohesion interaction, which can be implemented in parallel to a granular-frictional contact law. The evolution of damage at the contact level can be be easily modulated to enable a progressive transition from an initially linear elastic response to a loss of cohesion, by using a single parameter χ. In a straightforward numerical implementation, the effect of this contact model is presented for a 3D periodic boundary condition DEM code. The results from a series of simulations show that, for a constant peak resistance of the cohesion, a more progressive damage result in an increase of the peak stress in a particle assembly, as well as a continuous transition in the stiffness of the stress-strain response around the peak stress.


Introduction
The discrete element methods (DEM) are well suited to the study of granular packing of various materials, where the kinematics of individual particles in a granular assembly can be simultaneously resolved by the equations of motion. The sum of forces acting on the particles is dependant on their morphological properties, the contact laws governing their interactions, as well as the spacial organisation of the particles in a complex connectivity network. Even in the simplest implementations of individual particle geometries and contact laws, such as dry frictional contacts for spherical particles, the spacial organisation and evolution of granular systems during loading have proven to successfully reproduce the rich phenomenological behavior of various types of granular materials. The DEM has therefore found numerous applications in the study of geological processes, agronomic transformations and powder technologies.
For applications involving different types of granular media, different enriched DEM contact models have been developed to take into account complex physico-chemical inter-particle interactions and coupling effects with the interstitial solid and fluid phases in the porous media [1]. The contact models are usually formulated to capture the physical behavior of a particular granular media while retaining a low computational cost. This last requirement is of great importance for the DEM framework to remain effective in its application to large modeling domains. In this case long simulation times are expected due to the large number of particle interactions being computed and updated at every time step.
In this context, we present a damageable cohesive contact model which accounts for the progressive loss of cohesion between particles using a single additional damage parameter. It is here generally introduced within the limiting cohesion criterion for mixed interactions presented in [2], with a parallel implementation of the granularfrictional and cohesive parts of the contact. Hence, the model intends to describe spherical particles bonded together with solid-matter bridges, which can become damaged by erosion and micro-cracking due to deformations in the contact. We then present the results for a simple implementation of this modified contact law for DEM simulations performed using elementary volumes (REV).

Damageable contact law
Within a general three-dimensional DEM framework, based on the approach introduced by [3] and sometimes referred to as "molecular dynamics", the particles are modeled as soft spheres which are, in the event of a detected contact, allowed to overlap and generate repulsive forces. These forces are usually decomposed between the normal direction to the contact orientation and the frictional forces acting in the tangential plane perpendicular to the contact orientation (figure 1a). Here, a linear repulsion in the normal direction, with the classical Coulomb friction in the tangential direction, is employed; in principle, there is nothing to prohibit the use of more sophisticated models if appropriate. For example, a Hertz-Mindlin model could be used for particle repulsion, and a Bowden-Tabor model [4] for adhesive friction if surface roughness are to be considered. Additionally, a state of solid cohesion may exist between two particles, for which the forces are also decomposed in the same normal direction and tangential plane. The components of the resulting force vector for each interaction is therefore decomposed both in the spacial directions (figure 1a) and in a parallel action of two independent mechanisms, namely the granular-frictional and the cohesive parts of the contact law (figure 1b). According to this decomposition, and introducing a damage coefficient D to the cohesive part, the resulting force between two particles, "q" and "r", can be written as where the denomination "gran" indicates the part of the interaction generating a granular-frictional force, i.e. the normal repulsive and tangential frictional forces; while the denomination "coh" indicates the part of the interaction contributed by the (damageable) cohesive forces. The damage variable D evolves depending on the state of the contact between 0 to 1, according to the cohesion criteria detailed in Sect. 2. The contact law governing the mechanical response between each paired interaction relates the forces to the relative motion of the particles in a paired interaction. In the present model, a simplified approach to the Hertzian contact model is considered, using an independent linear elastic relation for each of the force component. The calculation of the 3D displacements δ n and incremental displacement δ t are based on an objective formulation [5]. Figure 2 illustrates the behavior of the contact laws. As illustrated in this schematic representation, the granularfrictional contact state is completely determined by the relative position of the two overlapping particles. For this part of the contact, the forces in the normal and tangential directions vanish if no contact is detected. Contrarily, a cohesion criterion is necessary to obtain the state of the solid cohesive part of the interaction and, if considered, the value of the associated damage variable (D).

Brittle cohesion criterion
The state of solid cohesion, which may exist between any two particles in the REV, is tested against a cohesion breakage criterion at every time-step. This criterion can be formulated as a surface in a hyperspace accounting for the Figure 2. Normal and tangential action of the granular-frictional and cohesive contact laws.
combined contribution of the different deformation mechanisms acting at the contact point of the paired interaction [2]. In a formulation accounting for the normal and tangential effects, it is expressed as where δ n and δ t are respectively the normal and tangential displacements. The subscript "0" in equation (2) denotes the intercepts of the function with the principal axis and are prescribed parameters influencing the level of resistance of the cohesive forces to the different deformation mechanisms. In a 2D parametric space, equation (2) is represented by an open ended quadratic surface around the δ n axis (see figure 3). For this breakage surface, the criterion cannot be satisfied from the single action of normal compression at the contact point. The integer parameter α determine the shape of this surface and thus the resistance of the bond to the combined effect of the different types of motions.
In the instance of an instantaneous loss of cohesion, i.e. no progressive damage, the cohesive part of the contact vanishes as the breakage surface is attained within any given time step. Therefore, the damage variable holds a binary value, where the cohesive contact is either linear elastic (D = 0) or broken (D = 1). However, in such an implementation, the abrupt change in the cohesive state between a pair of particles induces a discontinuity in the local contact force.
In fact this type of event instantaneously releases stored strain energy from individual contacts into the particulate arrangement, and thus may result in kinematic instabilities as well as large fluctuations in the neighboring contact forces with the propagation of elastic waves. Such event may result from small relative displacement in the contact or numerical noise as the breakage surface is approached.

Damage cohesion criterion
To mitigate the non-physical effect of an abrupt discontinuity in the contact force, a damageable cohesive law is introduced at the contact level. This extended formulation of the initial cohesive induces a progressive transition between the peak resistance of the cohesive contact, as the breakage criterion is first satisfied, and a complete loss of cohesion. The modified damage criterion from equation (2) takes the form where χ ∈ [1, χ ] introduces a homothetic transformation of the initial criterion. The damage threshold parameter χ , at which the cohesive part of the contact is lost (i.e. D = 0), is a user-defined constant for the model. The evolution of the damage variable D, associated with each cohesive contact, is therefore relative to the current surface (χ) and the maximum surface (χ ). For each time step, where "i" denotes the value of the variables from the previous time step and "f" the value of the variables calculated in the current time step, if G D (δ f n , δ f t , χ i ) > 0, then the value of χ f is updated to restore the equality in equation (4). Otherwise, it remains unchanged, i.e. χ f = χ i . Additionally, as the cohesive contact is incrementally damaged, the cohesive state is not restored in the case of unloading, and retains the maximum damaged state experienced by the contact, i.e. χ f ≥ χ i . The damage variable D in the cohesive part of the contact law is thus calculated as with the different cases (1-4) illustrated in Figure 3, for an idealized loading path in the δ t − δ n plane. The damage variable evolves from: (1)   The damage variable D of a contact is therefore completely determined by the maximum value of χ experienced in the contact, relative to the evolution of the normal and tangential displacements, as well as the single additional model parameter χ . With a value of χ = 1, the damageable model is seen to reduce to the binary brittle model. As higher values of χ are selected, the reduction of the cohesive force becomes more progressive before breakage.
The presented damageable contact law can be easily extended to account for additional deformation mechanisms, such as relative rotation of the particles. Additional terms can simply by added to equation (4). Other types of hypersurfaces may also be used to form breakage and damage criteria, which could account for the effect of compressive forces.

DEM implementation
The contact model presented in Sect. 2 is implemented in a DEM code named PBC3D [6]. This code was specifically developed as a compact and efficient DEM algorithm to be used either as a standalone program or as part of a multiscale implementation, e.g. in a FEM-DEM architecture. The domain geometry consists of a three-dimensional elementary volume (REV) of spherical particles, bounded by periodic condition in all directions.
To evaluate the effect of the damageable cohesive law coupled with the granular-frictional law, a series of simulations is performed in biaxial compression loading conditions, where˙ 2 = 0 andσ 3 = 0, while the vertical axial strain 1 is monotonically increased. The subscripts 1,2,3 denote the principal directions of stress and strain, aligned with the geometry of the REV. The constant parameters of the model are summarized in table 1.  At the REV level, this results in larger fluctuations in the stress-strain response as the level of cohesion increases. This behavior is noticeable around the peak stress and during the post-peak regime. In contrast, in the second series of simulations, where χ increases, the maximum cohesive force in each paired interaction remains unchanged, but the cohesive forces are decreasing progressively. The higher peak stress in this case is attributable to the stabilization of forces in damaging interactions, which occurs by transferring forces to less damaged interactions in the assembly. As χ increases, it can also be noted that the fluctuations in the stress-strain response tend to decrease. This result in a smoother transition around the peak stress, where a stable hardening phase can be identify during the late pre-peak section of the curves. The parameter χ also influences REV response softening rate.  The influence of the damage parameter on the stress level and shape of the stress-strain curve should be considered when calibrating the cohesion parameters of the contact law. The network effect of individual interactions in a granular assembly of spherical particles is shown here to be important in the redistribution of forces during the progression of damage and loss of cohesion. Therefore, the addition of a damage parameter in the cohesive contact law is seen as a more physical representation of a state and evolution of solid cohesion between the particles.
The increased stability of the REV in the presence of a damage parameter also has advantageous implications for hierarchical numerical modeling using DEM to model the behavior of cohesive material [6,7]. In the context of a FEM-DEM implementation, the stress-strain relation of individual REVs at every FEM time step serves as a constitutive model at the integration points of the FEM mesh [8,9]. As such, the solution of the implicit scheme, solved at the macroscopic FEM level, relies on the stability of the REV to obtain a continuous spacial distribution of the nodal forces. Therefore, by curbing large fluctuations in the stress state, a better global convergence can be achieved in the transition towards the post-peak regime.

Conclusion
In this communication, we presented a one parameter damageable cohesive law, which extends the class of mixed interaction breakage criteria for cohesive bonding of discrete particles. This damage parameter can be modulated to, on one hand, reduce to a brittle loss of cohesion, or on the other hand, provide an increasingly progressive loss of cohesion at the contact level. The usefulness of this enriched criterion is demonstrated in its application to an open-ended superquadratic breakage surface for a normal and tangential relative motion of the particles. It can however be extended, without technical difficulties, to other types of surfaces and additional deformation mechanisms at the contact level.
In the presented application to DEM, for a 3D periodic boundary condition elementary volume, the damage criterion is shown to result in an increased contribution of the cohesion in the macroscopic response. It results in an increase of the peak stress, and a smoother transition into the post-peak regime. A stable evolution of the response, encouraged by the effect of a continuous damage of the cohesive interaction, is desirable for applications which require a stable convergence of the resulting stresses, such as in the case of FEM-DEM coupling.