Contact Dynamics modeling of viscoelastic granular materials using irregular polyhedral particles

Viscoelastic granular materials are present in several disciplines. One example is asphalt mixture employed in road construction. In the last three decades, discrete element modeling has been positioned as a valid tool for the analysis of this multiphase material at the grain-scale. All this despite the simplification of the shape of the particles used in these studies. In this work, it is proposed a simplified procedure for the generation of viscoelastic granular samples composed of irregular polyhedra. The numerical aggregates were generated by a Poisson-Voronoi tessellation based on the particle size distribution (PSD) and statistic data of aggregates, without using complex imaging technics. This procedure set the porosity of the packing, while controlling the PSD. Using this procedure implies a significant computational-time reduction by skipping several preparation stages for polyhedral samples, such as deposition by gravity and compaction. This approach can be used for the study granular materials as inclusion in a solid matrix as concrete or asphalt mixtures, particle breaking, and fatigue damage of viscoelastic materials.


Introduction
Granular materials are present in different disciplines such as agriculture, the pharmaceutical industry and civil engineering, among many others. In construction, crushed rock aggregates are used as filling materials, and for manufacturing concrete and asphalt mixtures. For the latter, solid grains are surrounded by a viscoelastic matrix composed of mastic of bitumen and filler. The mechanical behavior of this multiphase material is highly dependent on the properties of its individual components in interaction.
Laboratory experiments are able to quantify the viscoelastic properties of these materials, but struggle to get a micromechanical insight. To study the properties of these materials at the particle-scale, the discrete element approach was employed over the last 30 years. In most of these studies, the particles were modeled as spheres or clumps composed of spheres [1][2][3][4], neglecting the crucial role of particle shape on granular fabric properties such as porosity, contact anisotropy, force chains network, etc. [5][6][7][8][9]. On the other hand, the preparation of numerical samples composed of polyhedral particles is a particularly time-consuming process. Complex algorithms are used in contact detection and repulsive force computation, during the preparation stages such as gravity deposition, compaction and stabilization of numerical samples.
In this work, a simple methodology is proposed for the preparation of viscoelastic granular samples composed of irregular polyhedral particles, generated as standard Voronoi tessellations. This approach reduces the computational time spent during the preparation process, while * e-mail: juan-carlos.quezada@insa-strasbourg.fr controlling the particle size distribution (PSD) and the porosity of the packing.

Numerical protocol
To identify the mechanical properties of viscoelastic materials such as asphalt mixtures, a common practice is to perform a complex modulus test in a two point bending (2PB) configuration [10]. To validate the proposed numerical procedure, in this study we propose to confront the experimental and numerical data of complex modulus tests obtained for frequencies ranging from 3 up to 40 Hz, and temperatures varying between -10 and 30 • C. The viscoelastic properties to analyze will be the dynamic modulus |E * | and phase angle Φ obtained from these trials.

Particles generation
The numerical samples are composed of rigid irregular polyhedral particles generated by a standard Poisson-Voronoi tessellation of the trapezoidal prism for a 2PB test, using the NEPER software [11], which generates randomly polycrystals as tessellations. About 3,800 tessellations are created following the experimental PSD (figure 1) cut at 2 mm, in order to reduce the total quantity of elements in the sample, where the fines are included in the mortar phase. Figure 2a displays the generated tessellations using the Poisson-Voronoi method. Then, the vertices of each tessellation are recuperated to build convex polyhedrons, generating triangular faces placed on the convex hull. The created particles have an average number of vertices equals to 14.81 and an average face number of 25.62. These quantities of vertices and faces are quite similar to those identified from digitalized crushed aggregates [9,12]. This protocol creates convex polyhedrons, based on simplified shapes of actual aggregates, without using complex imaging technics. A snapshot of the generated numerical sample with this approach is displayed in figure 2b. After the numerical generation stage, each sample is fixed to a bottom and a top plate. The latter simulates the mobile plate in the experimental set-up for a 2PB test.
To take in account the void content of 4.7%, the particle volumes are decreased, to generate voids between them. Then, the density of particles is set to 3556 kg.m −3 , to reach the total sample mass of 0.6 kg, concentrating the mass of the aggregates and the mortar phase within the particles. The particle-particle and particle-wall coefficient of friction is set to 0.7, which is a typical value for crushed aggregates.

Contact Dynamics method
The numerical simulations were performed using the LMGC90 software, which handles the interaction of rigid polyhedral particles [9,[13][14][15][16][17] . This software is based on the Contact Dynamics (CD), which is a discrete element approach for the simulation of non-smooth granular dynamics [18][19][20]. In this method, the equations of motion for each particle are formulated as differential inclusions in which velocity jumps replace accelerations [21,22]. The principal difference between this method and classic DEM approaches lies in the expression of the contact laws as complementary relations between the contact forces and velocities. This approach ensures basically the mutual exclusion between particles without introducing regularized laws often used in explicit methods such as the distinct element method [23][24][25] or molecular dynamics [26][27][28][29]. The unilateral contact interactions and Coulomb friction law are treated as complementarity relations or set-valued contact laws.
The determination of the contact set for polyhedral particles is obtained by a "bounding box" method used to sort a list of neighboring particle pairs. Then, for each pair, the overlaps are calculated using the "common-plane method" [25]. In the case of an overlap, the contact plane is determined by means of the intersection between the two particles. This detection procedure is fairly quick and allows simulating large samples composed of polyhedral particles [12].

Viscoelastic contact model
To simulate a viscoelastic granular material, two conditions are employed for particles at contact. For rigid particles, to respect the mutual exclusion criterion, a Coulomb friction law together with the Signorini condition are applied. On the other hand, to model the mastic phase surrounding the particles, a viscoelastic model is used for distant contacts, creating a viscoelastic rod between two particles at contact. This viscoelastic model is based on the Burgers model, which is composed by a Maxwell model putted in series with a Kelvin-Voigt model (figure 3). At the macroscale, the stiffness and viscosities for the Maxwell and the Kelvin-Voigt parts correspond to K m , C m , K k , and C k respectively. At the particle-scale, the normal contact components are calculated as the product of the macroscale parameters by the contact cross-section over the initial gap value. Here, the cross-section is assessed taking the minimum radius of two particles at contact. More details about the numerical implementation for this contact model can be found in [30]. The best fit parameters for this contact model are displayed in table 1. These values are obtained from the experimental master curves for the dynamic modulus and the phase angle.  Figure 4 displays the comparison between the experimental and numerical master curves for the dynamic modulus |E * | and phase angle Φ respectively. Experimental and numerical tests were performed using the same sample overall properties, such as geometry, mass and density, where the employed PSD is displayed in figure 2. The dynamic modulus is obtained as the absolute value of the ratio between the peak stress to peak strain for a complex modulus test. The phase angle Φ can be described as the time lag between the sinusoidal stress and strain signals.

Validation of the numerical approach
To build each master curve, a reference temperature T re f equals to 15 o C was chosen. Then the translation of all isotherm values is performed by calculating the reduced frequency as: a T f = f ×10 a T , where f corresponds to each frequency value in each curve and a T is the shift factor, which translates the data applying the time-temperature superposition principle. These values are obtained from the fit with the temperatures ones using the Williams-Landel-Ferry (WLF) equation (Eq. 1): where C 1 and C 2 are 28 o C and 206.8 o C respectively, and T is the temperature. As expected, experimental and numerical data of |E * | collapse in the same curve. On the other hand, numerical Φ values follow the trend of experimental ones, despite some fluctuations around the average values.

Concluding remarks
In this work contact dynamic simulations of a viscoelastic granular material composed of irregular polyhedrons were presented. The particle-generation procedure was based on a Poisson-Voronoi tessellation using statistic criteria to reproduce the number of vertices and faces identified from digitalized aggregates, without using a 3D scan or other complex imaging technics. In these numerical samples, the particle size distribution and the porosity were controlled.
This procedure reduces significantly the computational time during the preparation stages. To validate this numerical approach, the dynamics modulus and the phase angle were calculated to analyze the viscoelastic behavior during a 2PB test. Regarding the master curves, the numerical results of these simulations were in a good agreement with experimental data.
In conclusion, the generation procedure presented in this work can be seen as an interesting alternative to the classic procedures for generating polyhedral samples. Future works will focus on the improvement of the tessellation process for particle generation, to consider the properties of granular fabric, and on the effect of void content on the mechanical response.

Acknowledgments
The work presented in this paper was financially supported by the French National Research Agency (ANR -Sol-DuGri project ANR-14-CE22-0019). The authors want also to acknowledge CST Colas France for providing the experimental data.