Micro–Macro Relationships in the Simulation of Wave Propagation Phenomenon Using the Discrete Element Method

The present work is aimed to investigate the capability of the discrete element method (DEM) to model properly wave propagation in solid materials, with special focus on the determination of elastic properties through wave velocities. Reference micro–macro relationships for elastic constitutive parameters have been based on the kinematic hypothesis as well as obtained numerically by simulation of a quasistatic uniaxial compression test. The validity of these relationships in the dynamic analysis of the wave propagation has been checked. Propagation of the longitudinal and shear wave pulse in rectangular sample discretized with discs has been analysed. Wave propagation velocities obtained in the analysis have been used to determine elastic properties. Elastic properties obtained in the dynamic analysis have been compared with those determined by simulation of the quasistatic compression test.


Introduction
Wave propagation is a fundamental phenomenon which is encountered frequently in different natural processes, for instance, earthquakes, and engineering problems such as non-destructive testing or structures subjected to impact loading. A range of numerical schemes are used in literature to discretize the wave equation in space and investigate the wave propagation phenomenon, e.g., finite different schemes as elastodynamic finite integration technique [1,2] and local interaction simulation approach [3], finite element method [4], spectral element method [5], discontinuous Galerkin method [6], or boundary element method [7]. The discrete element method (DEM) is often used for the analysis of different problems of geomechanics involving wave propagation [8][9][10]. Owing to its simplistic mathematical framework based on Newtonian equations of motion, the DEM has emerged as a numerical tool of frequent choice for investigating the systems with inhomogeneities, or discontinuities existing in the material or appearing under loading. Although the main reason to use the DEM in such problems is not its capability to represent elastic waves, but this feature is important since wave propagation is an inherent phenomenon in any dynamic problem.
In the DEM, a material is represented by a large assembly of discrete elements interacting with each other by contact. Various contact laws representing physical effects such as elasticity, viscosity, damage and friction [11][12][13][14] can be used to define the inter-particle contact interactions. It is the cumulative behaviour of these micro level interactions that determine the bulk properties of the material. Perhaps, the most important advantage of DEM with respect to the previously mentioned methods is the ability to investigate the physical phenomena occurring at the micro level and its relationship with the macroscopic behaviour. It gives an advantage to DEM over other methods in modelling problems of practical nature because the expected behaviour of a system can be obtained numerically, simply by choosing appropriate microscopic parameters. The influence of the particle level interaction in the DEM on macroscopic characteristics of wave propagation has been studied by many authors, for instance, Sadd et al. in [8], studied the effects of contact laws on wave attenuation and dispersion behaviour of granular material, whereas in [9] Sadd et al. mainly focused on studying the influence of material microstructure on wave propagation behaviour. Mouraille and Luding [10] investigated dispersion and frequency dependence of wave propagation properties of a regular granular media by exploiting micro-macro transition [15] between particle level interactions and global behaviour. In [16], O'Donovan and O'Sullivan presented a detailed study of the wave velocities and inter-particle contact stiffness using an ideal and relatively simple hexagonal assembly of uniform sized particles. O'Donovan et al. compared experimental results on a model cubical cell of soil with DEM and continuum analysis in [17] and hence was limited only for particular material properties.
In the above mentioned studies, wave propagation in cohesionless granular media has been considered. Propagation of seismic waves in a cohesive material (an intact rock) has been simulated by Toomey and Bean [18] using a regular hexagonal configuration of circular discs. The results have been compared with a high-order finite difference solution of the wave equation. Propagation of elastic waves in rocks and material dynamic properties have been investigated using discrete element particle models by Resende et al. [19]. As a matter of fact, modelling of wave propagation is an inherent part of all discrete element simulations of any problem involving dynamic loading, such as blasting [20], rock cutting [21,22], or rock fragmentation [23].
Micro-macro relationships are essential in material modelling with the DEM. A standard procedure consists in the calibration of the microscopic (contact) model simulating a laboratory quasi-static strength test, such as unconfined compression, shear, Brazilian, or triaxial compression test and apply the calibrated model to simulation of the investigated problem, cf. [21]. It is not always possible to verify if the effective macroscopic properties are represented properly in the discrete element model. The wave propagation problem gives such a possibility. Propagation of ultrasound waves is commonly used as a measurement technique to determine elastic properties. We will exploit the analogous possibility in the numerical analysis to verify if the micro-macro relationships for elastic constants obtained in the simulation of wave propagation correspond to the relationships determined by other methods.
This work presents 2D discrete element modelling of elastic wave propagation in a rock-type cohesive material. Performance of an irregular disc configuration in the simulation of longitudinal and shear elastic waves will be investigated for different values of the shear to normal contact stiffness ratio. Two separate cases will be investigated for the longitudinal mode-propagation in a bulk solid and propagation under uniaxial stress conditions e.g., in a bar. Numerical results will be compared to theoretical predictions for equivalent macroscopic elastic properties determined using the micro-macro relationships obtained by simulation of the uniaxial compression test.
The outline of the paper is as follows. Basic notions of elastic wave propagation are given in Section 2. The formulation of the discrete element method is briefly described in Section 3. Micro-macro constitutive relationships in the DEM have been introduced in Section 4. The micro-macro relationships for elastic properties, Young's modulus, shear modulus, and Poisson's ratio, have been determined theoretically and numerically by simulation of the uniaxial compression tests in Section 5. Finally, the results of DEM simulation of the elastic wave propagation have been presented in Section 6. The macroscopic elastic properties in the discrete model have been evaluated based on the longitudinal and shear wave velocities. The corresponding micro-macro relationships have been compared with those obtained by the simulation of the uniaxial compression test (UCT).

Formulation of the Discrete Element Method
Numerical simulations have been performed in this work using the discrete element code DEMPack [26], which has been validated in various applications [21,22,[27][28][29]. The underlying formulation of discrete element method implementation following main assumptions of Cundall and Strack [30] is presented below.

Equations of Motion
The DEM considers the dynamics of a particulate system. In this work, 2D DEM models employing initially bonded cylindrical particles (discs) have been used. The translational and rotational motion of the discrete elements is described by means of the Newton-Euler equations of rigid body dynamics. For the i-th element, where u i is the element centroid displacement in a fixed (inertial) coordinate frame x, ω ω ω i -the angular velocity, m i -the element mass, J i -the moment of inertia. Vectors F i and T i are respectively composed of the forces and moments due to the external load F ext i , due to the contact interaction with adjacent particles, f c ij and t c ij , and due to the external damping, F : where n c i is the number of elements in contact with the i-th discrete element, and s c ij is the vector connecting the centre of mass of the i-th element with the contact point with the j-th element ( Figure 1). In the present work, only the force-type contact interaction will be considered, resulting in zero values for the interaction moments t c ij of Equation (12). in Equations (11) and (12) in the present work are of non-viscous type and are given by: where α t and α r , are respective damping factors for translational and rotational motion.

Time Integration Scheme
A second order explicit central difference scheme [31] is employed to integrate Equations (9) and (10). For translation motion, the time integration operator at the n-th time step is given as: Similarly, the time integration scheme for the rotational motion can be formulated as: If required, the rotational configuration can be determined, however, for the disc elements used in this work, the evaluation of rotational configuration is not essential. The explicit time integration scheme used in DEM imposes a limitation on the time step due to the conditional numerical stability. The time step ∆t must not be larger than the critical time step ∆t cr , which is determined by the highest natural frequency of the system, ξ max as, The highest frequency ξ max can be evaluated by solving the eigenvalue problem defined for the entire system of connected particles, however, this would be computationally expensive. Analogously to the standard simplification proposed for the explicit FEM [32], the maximum frequency of the full system in the DEM can be estimated by natural frequencies of subsets of connected particles surrounding each particle, cf. [33], where ξ (i) max is the maximum natural frequency of the system of connected particles surrounding the i-th particle. The problem of the critical time evaluation can be simplified further by considering equivalent single degree mass-spring systems with the natural frequency where k (i) eff is the effective stiffness governing the motion of the i-th particle. Hence, the limit on the time step can be given by The effective stiffness k (i) eff depends on the normal and tangential contact stiffnesses, the number of particles connected to the i-th particle as well as contact directions, cf. [33]. In practice, the time step can be estimated approximately taking, cf. [34] ∆t ≤ α m min k max (25) where m min is the minimum mass and k max is the largest normal or tangential contact stiffness and α is the user specified parameter accounting for multiple contacts for each mass. For regular packings of equal particles with the same stiffness for all the contacts the parameter α can be determined analytically [33]. For irregular packing a safe value of the parameter α can be based on the results of numerical simulations [34].

Contact Model
The formulation of the contact model employs the decomposition of the contact force between two elements f c into the normal and tangential components, f n and f t , respectively: where n is the unit normal vector at the contact point ( Figure 1). The normal and tangential contact forces can be evaluated assuming different models [12,35]. Granular materials are usually modelled assuming cohesionless frictional contact [36], while rock-like materials, as well as various other materials, require cohesive contact models [13]. The present work has been focused on the elastic behaviour of the materials modelled with bonded particles, therefore a cohesive contact model is presented. Nevertheless, the formulation presented is valid for a cohesionless contact model, as well. An initial bonding between the neighbouring particles is assumed. The bonds are established between the neighbouring particles i and j, if the gap between the particles g satisfies the condition: (27) where g 0 max is a tolerance in the contact verification and the gap g is given by In the soft contact approach used here, the impenetrability condition is satisfied approximately and a certain overlap between the contact particles h is allowed such that, The normal and tangential particle interactions are often modelled by linear springs connected in parallel with dashpots, providing an additional mechanism to dissipate contact oscillations. However, in this work the main focus is to investigate the elastic wave propagation, therefore the elastic bonds are assumed here and no additional dissipative mechanisms in the contact model are taken into account ( Figure 2). The normal and tangential contact force components are evaluated assuming the linear relationships: where k n is the interface stiffness in the normal direction, k t -interface stiffness in the tangential direction, and u t -the relative displacement at the contact point in the tangential direction. The relative tangential displacement u t must be evaluated incrementally, cf. [33]: where u n t is the vector of the relative tangential displacement from the previous time step rotated to the present contact plane and ∆u t is the incremental relative tangential displacement with v n+1/2 t being the relative tangential velocity at the contact point determined as where (v c r ) n+1/2 is the relative velocity at the contact point and v n+1/2 rn its projection on the Cohesive bonds are broken instantaneously when the interface strength is exceeded in the tangential direction by the tangential contact force or in the normal direction by the tensile contact force where Φ n in the interface strength in the normal direction, and Φ t the interface strength in the tangential direction. After decohesion the contact is treated assuming a standard contact model with Coulomb friction. The normal contact force can be compressive only (F n ≤ 0) and the tangential contact force is limited by µ| f n | f t ≤ µ| f n | (38) where µ is the Coulomb friction coefficient.

Micro-Macro Relationships
Obtaining a required macroscopic behaviour by using suitable interparticle contact model and relevant parameters is one of the main difficulties in the DEM. The effect of microscopic parameters in the DEM on macroscopic response have been investigated considerably in literature [37][38][39]. The present work is aimed to investigate the micro-macro relationships between elastic macroscopic properties defined in terms of effective elastic moduli and microscopic DEM parameters. It has been widely indicated in the literature that the macroscopic stiffness parameters depend upon the microscopic parameters such as normal and tangential contact stiffness, k n and k t , respectively, particle sizeR, porosity, e, coordination number [40]. The micro-macro relationships in the DEM can be obtained using various theoretical [41] or numerical [42] methods.

Analytical Micro-Macro Relationships
Effective macroscopic properties of granular materials can be estimated in terms of micromechanical parameters using the kinematic Voigt or static Reuss hypotheses [43]. The kinematic assumption of uniform strain gives the following analytical formulae for the Young's modulus and Poisson's ratio for the isotropic packing of equal sized cylindrical discs [41]: where N c is the total number of inter particle contacts in the volume V. For the non-uniform particle sized assembly such as used in this work, the arithmetic mean of the radii is generally used. The specimen volume V can further be expressed in terms of the particle volume, V p and specimen porosity e as follows: where N p is the number of particles in the specimen and for the disc particle with the thickness, t. Thus, by taking into account Equations (41) and (42), the Equation (39) can be written as: where n c is the coordination number, a parameter defined as an average number of contacts per particle, By rearranging Equation (43), we get the micro-macro relationship for the Young's modulus in the form,

Numerical Micro-Macro Constitutive Relationships
Micro-macro constitutive relationships can be obtained numerically by performing the DEM simulations of laboratory tests, such as the unconfined compression test [44], triaxial test [45], Brazilian test [40], or shear test [46]. Dimensional analysis with the Buckingham π-theorem provides a suitable framework to establish the numerical micro-macro relationships from the DEM simulations in the dimensionless form [44,47].
According to [42] the dimensionless micro-macro constitutive relationships for the Young's modulus and Poisson's ratio can be proposed in the following form: where Ψ E and Ψ ν are certain scaling functions of dimensionless parameters k t /k n and porosity e. Following [48] the effect of the specimen porosity and its characteristics can be taken into account using the results of micromechanical considerations expressed by Equations (40) and (45). Thus, the dimensionless (46) and (47) can be rewritten as follows whereΨ E andΨ ν are scaling functions of the dimensionless parameter k t /k n .

Determination of the Micro-Macro Relationships
The dimensionless micro-macro relationships (48) and (49) will be determined numerically by simulation of the uniaxial compression test (UCT) and compared with the analytical ones given by Equations (40) and (45), which have been derived on the basis of Voigt's kinematic hypothesis. Figure 3 shows a square sample 50 mm by 50 mm discretized with 4979 disc shaped elements of nonuniform size with an average radius of 0.370 mm, the maximum and minimum radii being 0.652 mm and 0.218 mm, respectively. The particle packing is characterized by the coordination number, n c = 5.8 and porosity e = 0.096. The loading has been introduced by the flat plates moving with a constant velocity 5 mm/s and compressing the specimen through the contact with its top and bottom sides. The microscopic parameters used in these simulations have been the following: particle density ρ = 2000 kg/m 3 , normal contact stiffness k n = 7 × 10 10 N/m. Cohesive bonds strengths φ n = φ t = 2.9 × 10 4 N and Coulomb friction coefficient µ = 0.83 have been assumed. Additionally, a non-viscous background damping has been used in this example assuming the damping factors α t = α r = 0.2. Simulations have been conducted for values of the tangential to normal contact stiffness ratio k t /k n ranging between 0.0 to 1.0 in steps of 0.1, obtained by changing the tangential contact stiffness k t for each simulation run.   The stress-strain curve has been used to determine the macroscopic Young's modulus E: by taking the stress range from 0 to 50% of the maximum stress level in the considered simulation. The macroscopic Poisson's ratio has been calculated as where the increments of the strain components correspond to the range used in the determination of the Young's modulus. The macroscopic stress, σ has been evaluated in this work by using an averaging procedure involving the internal contact forces in the RVE [43]: where R RVE is the volume of RVE which in the present case is assumed to be the volume, V of the sample, and N c is the number of contacts in the RVE. f c denotes the contact force vector at the contact c, whereas L c is the vector connecting centroids of two contacting particles and known as the branch vector. The macroscopic strain tensor for the discrete element assembly has been calculated using the procedure proposed in [40]. Averaging is performed over a triangular mesh generated over the centres of the particles forming the specimen. This is a two level averaging procedure. First, a constant strain ε k in all the triangles are determined using the formula derived from the averaging equation: where S k is the area of an elementary cell. Applying the divergence theorem, the surface integral in Equation (53) can be transformed into the line integral where L k is the closed boundary of the triangular element, u-the displacement field and n-the unit normal vector outward to the element. The line integral in Equation (54) is evaluated in terms of nodal displacements and geometric parameters characterizing the triangular element. Having determined the strain ε k ij in each element the average strain tensor in the whole specimen is obtained by the weighted averagingε With known porosity e, normal contact stiffness k n , coordination number n c , unit particle thickness t, and the Young's modulus determined from Equation (50), the scaling functionΨ E for the given ratio k t /k n can be calculated according to Equation (48). The values of the scaling functionΨ ν defined by Equation (49)  The numerical values have been compared with the analytical predictions according to Equations (40) and (45), which have been derived on the basis of Voigt's kinematic hypothesis. The close agreement between analytical and numerical values of dimensionless elastic parameters observed here verifies the correctness of dimensionless micro-macro relationships obtained in a static unconfined compression test simulation.

Discrete Element Simulation Setup
The propagation of the plane elastic longitudinal and shear waves has been simulated using a rectangular sample (cf. Figure 7) discretized with 682 bonded disc elements with parameters shown in Table 1. Simulations have been performed for the ratio k t /k n ranging between 0.0 to 1.0 in steps of 0.1. Particles forming the left edge of the sample are free, whereas the right edge of the sample is fixed. Longitudinal waves have been simulated for propagation in a bulk medium and in a bar. For the bulk model, the elements comprising the top and bottom edges are allowed to move only in the direction of wave propagation, whereas free boundary conditions have been applied at the top and bottom edges to obtain uniaxial stress condition for simulating the longitudinal bar wave propagation. The periodic boundary conditions have been imposed for these edges for simulation of the shear wave. The wave impulse has been induced by assigning initial displacements to the particles on the free edge using the following function: where the position of the particles in x-direction is bounded within, 0 ≤ x ≤ L/2 in reference to the left edge of the sample. An amplitude A = 0.01 mm and wavelength L = 8 mm are assumed, resulting in number of elements per wavelength for the maximum particle radius, R max = 0.145 mm approximately equal to 28 (L/(2 × R max )). Thus it was ensured that the recommended 20 elements per wavelength are used analogously to the FEM [49]. The initial displacements in the x direction have generated the longitudinal wave, while the initial displacements in the y direction have been set to induce the shear wave. Wave propagation has been simulated assuming zero damping conditions. Breakage of cohesive bonds has been impeded by setting very high values of bond strength.   Figures 8 and 9 show the propagation of the longitudinal wave pulse in a bulk solid and in a bar, respectively, in terms of particle displacement vectors at different time steps for the ratio k t /k n = 0 . The shear wave pulse propagation in the discrete sample for the ratio k t /k n = 0 in terms of particle displacement vectors at different time steps is shown in Figure 10.      Figure 7) at shear wave propagation.
The peak-to-peak method has been used to evaluate numerical wave velocities from the evolution of particle displacements with respect to time. Displacement time graphs for selected particle pairs for the longitudinal bulk and bar wave have been plotted in Figures 11 and 12, respectively. Displacement time graph for selected particle pairs for the shear wave has been plotted in Figure 13. The time ∆t taken by the wave to travel between chosen nodes has been shown in these plots. An average of velocities for five pairs of discrete elements (highlighted with green color in Figure 7) has been used as the velocity for a particular k t /k n ratio. From the DEM simulations results shown in Figures 8-13, it can be observed that the longitudinal bulk wave pulse travels with the fastest velocity and the shear wave pulse propagates with the lowest velocity amongst all the three cases. This verifies that in fact by using the DEM the expected wave propagation behaviour can be reproduced. The computational cost (CPU time) for simulating wave propagation cases investigated is quite small. For instance, the CPU time required to simulate the propagation of shear wave impulse through the length of the discrete sample shown in Figure 7, is equal to 7 s with k t /k n = 0.0, critical time step ∆t cr = 1.941 × 10 −8 s and required 503 time steps. Dependence of the wave velocities on the stiffness ratio k t /k n has been shown in Figure 14 for the longitudinal bar wave and in Figure 15 for the shear wave. Numerical results presented in Figurs 14 and 15 have been compared with analytical ones. Analytical wave velocities have been determined as functions of the ratio k t /k n using Equations (2) and (6) with macroscopic Young's modulus E and Poisson's ratio ν of the discrete sample deduced from Equations (39) and (40), respectively for a given stiffness ratio, k t /k n with no. of contacts, N c = 1852 and average coordination number, n c = 5.43 for the DEM sample shown in Figure 7. A very good agreement between numerical results and theoretical predictions can be observed in Figures 14 and 15. Similarly, in Figure 16 where longitudinal bulk to shear wave velocity ratio, c bulk l /c s as a function of the ratio k t /k n has been plotted, overall numerical and analytical results agree well, except a few minor deviations at lower values of ratio k t /k n . The ratio of longitudinal bar to shear wave velocity, c bar l /c s as a function of the ratio k t /k n has been plotted in Figure 17, where a good agreement between numerical and analytical results can be observed as well.  The wave velocities obtained numerically for different k t /k n ratios have been used to evaluate macroscopic elastic parameters for the DEM sample. Equations (3) and (4) have been used to determine effective macroscopic Young's modulus and Poisson's ratio from the simulations of the longitudinal bulk wave propagation. Similarly, macroscopic Young's modulus and Poisson's ratio for longitudinal wave propagation in a bar have been analysed using Equations (7) and (8). Shear modulus for the transverse wave propagation has been determined from Equation (5). The elastic moduli determined by DEM simulations of the elastic wave propagation have been plotted as functions of the k t /k n ratio in Figures 18 and 19 in comparison with the results obtained by simulations of the uniaxial compression test. The Young's modulus and Poisson's ratio corresponding to the simulation of compression test have been evaluated using the dimensionless scale functions given in Figures 5 and 6, respectively, and these parameters have been used subsequently in the relationship, G = E/2(1 + ν) for the evaluation of shear modulus.
A very good agreement between the Young's moduli obtained through dynamic and quasistatic numerical simulations in the full range of ratio k t /k n between 0.0 and 1.0 can be seen in Figure 18. Similarly, a close concurrence between Shear moduli obtained from the dynamic and quasistatic numerical simulations for the entire range of the ratio k t /k n can be noticed in Figure 19. A certain deviation of dynamic results from the quasistatic data can be observed in Figure 20 for some values of the Poisson's ratio. This can be explained by analysing the form of Equation

Conclusions
Numerical simulations have confirmed the capability of the discrete element method to represent properly phenomenon of the elastic wave propagation in a solid material both in longitudinal and shear modes. The wave propagation velocities agree well with theoretical predictions. This allows to use wave velocities to determine elastic macroscopic properties of solid media discretized with discrete elements. Comparisons with micro-macro relationships obtained numerically by simulations of the quasistatic compression test show that the values of the Young's and shear moduli based on the wave velocities agree very well with those determined by simulations of the quasistatic compression test, however, there is a certain difference for the Poisson's ratio in the range of its low values (lower than 0.1). Calculation of the Poisson's ratio in this range is very sensitive to inaccuracies in evaluation of wave velocities.