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.

1. 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).

2. Basic Notions of Elastic Wave Propagation

An elastic wave in a solid body can propagate in two modes, in the form of the longitudinal (also called compression) and transverse (also called shear) waves. In the longitudinal wave, the motion of the material points is in the direction of propagation whereas, in the shear wave, the motion of the material points is in a plane perpendicular to the direction of propagation. Wave propagation velocity, being one of the main parameters characterizing waves, in solid media depends on elastic constants viz. Young’s modulus E, Poisson’s ratio , and shear modulus G, with the following relationships [24]: for longitudinal wave velocity, and for the shear wave velocity, with denoting the bulk density. Knowing the wave velocities, the Poisson’s ratio and the Young’s modulus can be deduced by combining Equation (1) with Equation (2) in the following form: and from the Equation (2), the shear modulus is obtained as:

The effective macroscopic elastic properties will be determined from the simulations of the bulk longitudinal and shear wave propagation using Equations (3)–(5).

Alternatively, the elastic properties can be determined combining the simulation of the shear wave and the simulation of the longitudinal wave propagation under the uniaxial stress conditions (like in a bar). The relationship for shear wave velocity remains the same as given by Equation (2), and the longitudinal bar wave velocity is given by [25]: The Young’s modulus can be obtained directly from Equation (6): and combining relationships (6) and (2) the Poisson’s ratio can be obtained as: The shear modulus remains the same as defined by Equation (5).

3. 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.

3.1. 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 is the element centroid displacement in a fixed (inertial) coordinate frame , —the angular velocity, —the element mass, — the moment of inertia. Vectors and are respectively composed of the forces and moments due to the external load , due to the contact interaction with adjacent particles, and , and due to the external damping, and : where is the number of elements in contact with the i-th discrete element, and 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 of Equation (12).

Figure 1

Definition of the inter-particle interaction.

The damping terms and in Equations (11) and (12) in the present work are of non-viscous type and are given by: where and , are respective damping factors for translational and rotational motion.

3.2. 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 must not be larger than the critical time step , which is determined by the highest natural frequency of the system, as, The highest frequency 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 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 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 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] where is the minimum mass and 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].

3.3. Contact Model

The formulation of the contact model employs the decomposition of the contact force between two elements into the normal and tangential components, and , respectively: where 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: where 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 is the interface stiffness in the normal direction, —interface stiffness in the tangential direction, and —the relative displacement at the contact point in the tangential direction. The relative tangential displacement must be evaluated incrementally, cf. [33]: where is the vector of the relative tangential displacement from the previous time step rotated to the present contact plane and is the incremental relative tangential displacement with being the relative tangential velocity at the contact point determined as where is the relative velocity at the contact point and its projection on the normal direction:

Figure 2

Rheological scheme of the bonded particle interaction model.

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 in the interface strength in the normal direction, and 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 () and the tangential contact force is limited by where is the Coulomb friction coefficient.

4. 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, and , respectively, particle size , porosity, e, coordination number [40]. The micro–macro relationships in the DEM can be obtained using various theoretical [41] or numerical [42] methods.

4.1. 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 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, and specimen porosity e as follows: where 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 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,

4.2. 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 and are certain scaling functions of dimensionless parameters 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 and are scaling functions of the dimensionless parameter .

5. 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, = 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 kg/m, normal contact stiffness N/m. Cohesive bonds strengths N and Coulomb friction coefficient have been assumed. Additionally, a non-viscous background damping has been used in this example assuming the damping factors . Simulations have been conducted for values of the tangential to normal contact stiffness ratio ranging between to in steps of , obtained by changing the tangential contact stiffness for each simulation run.

Figure 3

Uniaxial compression test (UCT) of an irregular configuration of nonuniform size particles used for establishing dimensionless constitutive relationships.

Figure 4 presents the macroscopic stress–strain curve for one of the analyses.

Figure 4

The stress-strain curve for the uniaxial compression numerical test for .

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 is the volume of RVE which in the present case is assumed to be the volume, V of the sample, and is the number of contacts in the RVE. denotes the contact force vector at the contact c, whereas 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 in all the triangles are determined using the formula derived from the averaging equation: where 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 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 in each element the average strain tensor in the whole specimen is obtained by the weighted averaging

With known porosity e, normal contact stiffness , coordination number , unit particle thickness t, and the Young’s modulus determined from Equation (50), the scaling function for the given ratio can be calculated according to Equation (48). The values of the scaling function defined by Equation (49) are obtained directly from Equation (51). The values of the dimensionless parameter involving the Young’s modulus and Poisson’s ratio for all the simulated cases in the full range of the ratio have been plotted in Figure 5 and Figure 6, respectively.

Figure 5

Dimensionless micro–macro relationship for the Young’s modulus as a function of ratio .

Figure 6

Dimensionless micro–macro relationship for the Poisson’s ratio as a function of ratio .

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.

6. Simulation of the Wave Propagation

6.1. 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 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, in reference to the left edge of the sample. An amplitude mm and wavelength mm are assumed, resulting in number of elements per wavelength for the maximum particle radius, = 0.145 mm approximately equal to 28 (). 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.

Figure 7

Discrete element sample used in numerical simulations of wave propagation. Disc elements in green color indicate the nodes used to measure the wave velocity in the sample.

Table 1

Parameters of discrete element method (DEM) sample used in numerical studies of wave propagation.

Symbol	Parameter	Value	Units
N_p	no. of particles	682	-
R_max.	max. radius	0.145	mm
R_min.	min. radius	0.1	mm
l	sample length	16.54	mm
w	sample width	2.3	mm
e	porosity	0.11	-
ρ_p	particle density	2000.0	kg/m^3
ρ_avg.	sample average density	1784.26	kg/m^3
k_n	normal contact stiffness	1 × 10^10	N/m

6.2. Numerical Results

Figure 8 and Figure 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 = 0. The shear wave pulse propagation in the discrete sample for the ratio = 0 in terms of particle displacement vectors at different time steps is shown in Figure 10.

Figure 8

Longitudinal bulk wave propagation through DEM sample. Snapshots of particle displacement vectors captured at time—(a) s, (b) s, (c) s, (d) s.

Figure 9

Longitudinal bar wave propagation through DEM sample. Snapshots of particle displacement vectors captured at time—(a) s, (b) s, (c) s, (d) s.

Figure 10

Shear wave propagation through DEM sample. Snapshots of particle displacement vectors captured at time—(a) s, (b) s, (c) s, (d) s.

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 Figure 11 and Figure 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 ratio. From the DEM simulations results shown in Figure 8, Figure 9, Figure 10, Figure 11, Figure 12 and Figure 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 = 0.0, critical time step = 1.941 × s and required 503 time steps. Dependence of the wave velocities on the stiffness ratio has been shown in Figure 14 for the longitudinal bar wave and in Figure 15 for the shear wave. Numerical results presented in Figure 14 and Figure 15 have been compared with analytical ones. Analytical wave velocities have been determined as functions of the ratio 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, with no. of contacts, = 1852 and average coordination number, for the DEM sample shown in Figure 7. A very good agreement between numerical results and theoretical predictions can be observed in Figure 14 and Figure 15. Similarly, in Figure 16 where longitudinal bulk to shear wave velocity ratio, as a function of the ratio has been plotted, overall numerical and analytical results agree well, except a few minor deviations at lower values of ratio . The ratio of longitudinal bar to shear wave velocity, as a function of the ratio has been plotted in Figure 17, where a good agreement between numerical and analytical results can be observed as well.

Figure 11

Evolution of x-displacement with time for nodes A and B cf. Figure 7 at longitudinal bulk wave propagation.

Figure 12

Evolution of x-displacement with time for nodes A and B cf. Figure 7 at longitudinal bar wave propagation.

Figure 13

Evolution of y-displacement with time for nodes A and B (cf. Figure 7) at shear wave propagation.

Figure 14

Comparison of longitudinal bar wave velocity, in DEM sample with analytical results as a function of ratio .

Figure 15

Comparison of shear wave velocity in bulk solid, in DEM sample with analytical results as a function of ratio .

Figure 16

Comparison of longitudinal bulk to shear wave velocity ratio, in DEM sample with analytical results as a function of ratio .

Figure 17

Comparison of longitudinal bar to shear wave velocity ratio, in DEM sample with analytical results as a function of ratio .

The wave velocities obtained numerically for different 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 ratio in Figure 18 and Figure 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 Figure 5 and Figure 6, respectively, and these parameters have been used subsequently in the relationship, for the evaluation of shear modulus.

Figure 18

Comparison of Young’s modulus obtained using longitudinal bulk wave velocity with shear wave velocity and longitudinal bar wave velocity with shear wave velocity in DEM sample, with UCT numerical results as a function of ratio .

Figure 19

Comparison of shear modulus obtained using shear wave velocity in DEM sample with UCT numerical results as a function of ratio .

A very good agreement between the Young’s moduli obtained through dynamic and quasistatic numerical simulations in the full range of ratio 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 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 (3) for the wave propagation in a bulk solid and Equation (8) for the wave propagation in a bar. It can be deduced that the numerical evaluation of Poisson’s ratio is highly sensitive to the numerical evaluation of the wave velocity ratios and . Even a small numerical error in the evaluation of wave velocities leads to the square of the error in Poisson’s ratio calculation. For small values of the Poisson’s ratio (close to zero) the numerical error can be of the order of the evaluated parameter value.

Figure 20

Comparison of Poisson’s ratio obtained using longitudinal bulk wave velocity with shear wave velocity and longitudinal bar wave velocity with shear wave velocity in DEM sample, with UCT numerical results as a function of ratio .

7. 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.