The Dynamic Lattice Method - Vibrations and Wave Propagation in Discontinuous and Heterogeneous Media


 The development of a new dynamic lattice element method (dynamicLEM) as well as its application in the simulation of wave propagation in discontinuous and heterogeneous media is the focus of this research paper. The conventional static lattice models are efficient numerical methods to simulate crack initiation and propagation in cemented geomaterials. The advantage of the LEM and the developed dynamic solution, such as simulation of arbitrary crack initiation and propagation, illustration and simulation of existing inherent material heterogeneity as well as stress redistribution upon crack opening, opens a new engineering field and tool for material analysis. To realize the time dependency of the dynamic LEM, the governing Newton's second law is solved while using the Newmark-β method and implementing the non-linear Newton-Raphson Jacobian. The method validation is done according to the results of a boundary element method (BEM) in the plane P-SV-wave propagation within a plane strain domain. Further validation tests comparing the generated wave types, simulation and study of crack discontinuities as well as inherent heterogeneities in the geomaterials are conducted to illustrate the accurate applicability of the new dynamic lattice method. The results indicate that with increasing heterogeneity within the material, the wave field becomes significantly scattered and further analysis of wave fields according to the wavelength/heterogeneity ratio become indispensable. Therefore, in a heterogeneous medium, the application of continuum methods in relation to structural health monitoring should be precisely investigated and improved. The developed dynamic lattice element method is an ideal simulation tool to consider particle scale irregularities, crack distributions and inherent material heterogeneities and can be easily implemented in various engineering applications.


Introduction
Current dynamic simulation methods in civil engineering and material science are mainly continuum based methods, such as boundary element, finite element, finite difference and volume element methods. These methods have been developed while considering a wide range of boundary conditions from the 70's until nowadays. However, the simulation of structural vibrations or wave propagation in solids under the consideration of damaged or heterogeneous structures can not be easily performed using the continuum based methods. In particular, the consideration of local damages or a stochastically distributed heterogeneity is difficult to realize. The structural health monitoring methods are widely considered and implemented in engineering applications to ensure the usability as well as the state analysis of structures, and finally the failure of structures [1]. In particular, the right identification of the structural state needs a consideration of all disturbances within the structure, like small damages or material heterogeneity. In general, the dynamical methods differ in structural dynamics based on the analysis of waves to determine existing damages and wave based methods to identify the location of the discontinuities [2], often with high frequency waves to detect the damages [3,4]. So far during the past decades, mainly continuum based approaches were developed to detect the damages in the structures [5]. Meshram et al. (2015) [6] utilized finite element analysis to analyse and detect cracks in a cantilever beam. Advanced wavelet-based neural networks are considered in the detection of damages in beam-like structures [7]. Recently, a numerical scheme is proposed by Agathos et al. (2018) [8] for the detection of multiple cracks in three dimensional (3D) structures.
Besides the continuum based methods, the development of discrete methods have become very popular to model micro and meso-scaled problems, mostly to understand the phenomenon of material behaviour. Discrete methods, such as cohesive Discrete Element Method (CDEM), are used, e.g. to simulate the fracking in cemented geomaterials [9,10]. The standard Discrete Element Method (DEM) on the other hand, as an inherent Lagrange based method, is considered for simulation of large particle movements in non-cohesive granular materials [11]. Application of CDEM is also expanded in different scientific fields, such as simulation of crack propagation in a vitreous biopolymer material [12]. Coupled FED-DEM methods are developed to tackle the computational costs that are associated with discrete methods and are able to simulate fracking while considering the irregularities that exist in geomaterials [13]. However, mainly due to the great numerical costs, the discrete simulation methods were barley considered until now to simulate a full dynamic or wave field problem in engineering applications. To decrease the computational costs and lower the problems of complexity that exist for analysing discrete structures, the so-called lattice element method (LEM) has been developed. The initial idea was to analyse the crack propagation in solid and heterogeneous materials, like concrete [14].
The LEM can be considered as a hybrid method, where the generated domain can be either seen as a continuum or packed discrete particles [15]. The lattice elements are able to simulate the frack initiation and propagation in solids and cemented geomaterials, where small particle movements or displacements are expected [16,17,18,19]. One of the main advantages of the LEM related to other numerical methods is its ability to simulate the stress redistribution and concentration upon the cracking process. The domain can be represented in the most easy way by a series of springs [16], more complex than Euler Bernoulli or Timoshenko beams [20]. The generalized beam lattice model, where mechanical bond-aggregate interfere properties are implemented in a lattice element, is also used to simulate the crack growth in concrete [21]. The Rigid-body-spring networks are also similar to lattice elements, representing the domain with a series of spring elements. In these methods, between each unit cell, three spring elements representing a resistance against axial, transverse and rotation displacements are considered [22,23,24]. The initial step in application of lattice model is to grant the mesh-independence outcome of the simulations.   [20] introduced the embedded strong discontinuity into lattice elements, which results in meshindependent computations of failure response. In another approach, the correlation between single element and continuum properties is achieved using the strain energy theory, where a strain energy stored in a unit cell is equal to the strain energy stored in a continuum [25,26]. The crack simulation in lattice elements is performed with removing an element reaching its strength threshold (brittle failure) or decreasing the stiffness according to a bi-linear softening scheme. The failure of elements are defined based on the Mohr-Coulomb failure model with tension cut off [21,22] or based on mode I and II fracture mechanisms, where a critical stress intensity factor is used to calculate the released strain energy rate [27,16,28]. In all the named developments the focus was on the study of fracture propagation in materials.
During recent years, the application of lattice elements has been extended to simulate the coupled thermo-hydro-mechanical processes. The heat transfer and change of effective thermal conductivities in cemented geomaterials [29], in unsaturated granular geocomposites [30] or modified geomaterials [31] have all been investigated. Grassl (2009) [32,33] developed hydro-mechanical dual lattice models to simulate the flow and permeability changes in geomaterials. The dual lattice model was also implemented to investigate the fluid driven percolation and developed flow paths in salt and clay stone samples under anisotropic confining stresses [15]. The first initial approach to extend the common lattice element method into the field of dynamics for use of wave field simulation in heterogeneous materials was done by   [34]. Here, the dynamic LEM is used to simulate wave fields considering a progressive fracking process in brittle and quasi-brittle material under dynamic loading. Also, Rizvi et al. (2020) applied the dynamic LEM to solve the problem of mechanical waves in rock mass or cemented granular material under dynamic loadings [35]. The simulation of crack propagation under dynamical forces by an embedded strong discontinuity approach is also implemented by   [36]. The simulation results of the dynamicLEM are used to train the artificial neural networks (ANN) model in order to detect the location of the discontinuities [37], which reduces the computational costs.
This research paper extends the application of the lattice element method in dynamic structural analysis to tackle the effect of discontinuities and inherent material heterogeneities on wave fields. The inherent mesh irregularity of the lattice element and stochastic fracking process makes the model suitable for simulation of discontinuities under expected small deformations. Initially, the implemented mathematical methodology is explained, where Newton's second law within the equation of motion is solved in time domain to determine the accelerations, velocities and displacements of generated nodes. Afterwards, validation benchmarks are presented and the results of the dynamic lattice for a 2D domain is compared to the boundary element method (BEM) solution. Within that benchmark, the dynamic structural behaviours are compared in the frequency domain to demonstrate the accuracy of the method. To show the influence of discontinuities, series of full wave field simulations are performed and analysed. The visualization of the wave types within propagating wave field, the wave field scattering and dispersion at heterogeneities, as well as the wave field shadows around discontinuities illustrate the excellent ability of the new dynamic method and its future applications in structural dynamics and structural health monitoring.
2 Methodology and mathematical framework

Domain Discretization and Discontinuity
The general discretization of a domain is performed using the vectorizable random lattice [38]. In this method, the generated irregularity of a mesh is controlled using the defined randomness factor (α f ), where α f = 1 is the maximum irregularity and α f = 0 generates a regular mesh. After the generation of the nodes, which can all act generally as crack nuclei, the Voronoi Tessellation and Delaunay Triangulation are considered to mesh the domain by generating polygonal cells and their lattice connectivities ( Fig.1a and Fig.1b). Within the lattice model, these elements represent the domain and carry the mechanical loads [29]. The location, orientation and length of the generated discontinuities are in a stochastic manner as shown in Fig.1c. One of the advantages of considering lattice model for crack simulations is the simple generation of a heterogeneous granular domain. Therefore, this numerical method can easily be used in crack simulations in e.g. concrete material, where granular particles are cemented with each other [39,40]. In these models, granular particles, cement (bond material) and interfaces are defined with different mechanical properties. Fig.2 illustrates two different granular packing layouts with uniform and nonuniform granular distributions. The arbitrarily distributed heterogeneity with four different minerals is shown in Fig.2c. To use the dynamic LEM in engineering applications, such as structural health monitoring, the visualisation and quantification of wave fields in heterogeneous materials has a great importance.

Composition of the Static Lattice Method
Similar to the Rigid-Body-Spring Network methods, the resistance to axial, transverse and rotational displacements are modeled using three spring elements [41]. Implementing the spring elements eliminates the simulation difficulties that arose from short beam elements, such as small beam aspect ratio. Each generated nuclei has three degrees of freedom and Fig.3 illustrates the generated Voronoi cells and the defined spring elements among them.
where A is the cross-section area (A = t × l), t is the thickness of the domain, l is the length of elements or Euclidean distance between nuclei, k s is the transverse stiffness, k n is the axial stiffness, k φ is the rotational stiffness, G is the shear modulus and E is the Young's modulus. While assuming the linear elasticity, according to Hooke's law, where T is the transformation matrix, K l is the stiffness matrix in local coordinates, f g is the vector of forces, K g is the stiffness matrix and u g is the vector of displacements in the global Cartesian coordinate.

Figure 3: The generated Voronoi cells and the defined spring elements
The regularisation of a lattice model is carried out according to the strain energy theory and energy balance, where the stored energy in a continuum (U R ) is considered to be equal to the stored strain energy of a unit cell (U cell ) [25]. For a spring element, the strain energy stored in a unit cell is determined according to the total number of bond elements (N b ), the elements response forces (f b ) and the corresponding displacements (u b ). In a continuum, the stored energy is calculated with the applied stresses (σ R ) and strains (ε R ) throughout its volume (V R ).
For a unit cell, the stored strain energy is calculated based on the length of the elements (l), first stiffness coefficient (R ′ ) and orientation of the unit normal vector (n i,j,k,m ).
For a continuum, the stored strain energy is calculated by where C R is a continuum stiffness matrix. Finally, a correlation between continuum properties such as Young's Modulus (E R ), Shear Modulus (G R ) and Bulk Modulus (K R ), with single lattice properties (E, G, K), is driven based on the continuum stress condition and with the linear elastic theory. According to this correlation, the desired Poisson's ratio (ν) of a continuum will also be achieved. The irregularity of the generated mesh (α f ) has an influence on the correlated ratios (α 1 , α 2 , α 3 ).
Under the assumption of linear elastic fracture mechanics (LEFM) and while considering the Mode I and II fracture mechanism within the model, the crack initiation and propagation can be simulated. The implemented failure criteria here is based on the Mohr-Coulomb failure model with tension cut off [42,43]. Fig. 4a depicts the failure criteria and defined cohesion (c) and tensile strength (f t ). The elements failure under compression forces is permitted to grant crack propagation. The compression strength (f c ) of the elements is assumed to be 10 times the tensile strength, where σ n is the axial stress, α ′ is a scaling parameter equal to 0.005 [44], M i,j are the moments in node i and j, h is the cross-section length of the elements, φ ′ is the friction angle and f s is the shear strength. Modeling a quasi-brittle material, such as concrete or rock, requires implementation of a softening scheme [44]. Fig.4b illustrates the implemented bi-linear softening scheme, where based on the softening ratio (SR), the stiffness of the element is reduced.
where σ p is the peak stress, ε p is the peak strain, ε f is the failure strain and ε is the elements strain. When SR = 1, then a brittle failure is simulated and the element stiffness is deducted to zero (removed).

Implementation of the Dynamic Lattice Method
The dynamic lattice method is based on the general equation of motion, known as Newton's second law, where the accelerations (ü), velocities (u) and displacements (u) of nodes in each time step (t) are calculated. The mass matrix (M ) is assembled according to the consistent mass matrix (CMM) by lumping the mass at the nodes. The global stiffness matrix (K g ) is assembled as previously shown in Eq.2 and Eq.1. The contact forces in the dynamic lattice method in contrast to discrete element models is equal to zero, as no contact search between the elements is required. This reduces the computational costs of the dynamic lattice models. The damping matrix (C) is also specified, although the contact damping is neglected here.
where f ext is a vector of external dynamic loads and ∆t is the time step. Different explicit and implicit methods can be used to solve the equation of motion, such as the central difference method, implicit linear acceleration method and Newmark-β method with incremental formulation. In this paper, Newton's second law is solved using Newmark-β method with incremental formulation as given below.
where γ, β are the Newmark-β parameters and δu, δu, δü are the increments of displacement, velocity and acceleration, respectively. The incremental formulation of Newton's second law is given as, The incremental formulation of the equation of motion can be solved under the assumption of different γ and β parameters. When β = 1 4 and γ = 1 2 , the Newmark method is unconditionally stable and implicit. If β = 1 6 and γ = 1 2 , the Newmark method is similar to the linear acceleration method. Under the assumption of β = 0 and γ = 1 2 the Newmark method is identical to the central difference method [45]. The Newton-Raphson Jacobian is implemented here to solve the system of nonlinear equations with multiple variables.
where n is the numerical iteration number and J −1 is the inverse of the Jacobian matrix, which is equal to the derivative of the F (δu n ) with respect to each unknown node displacements. The convergence of the dynamic lattice method in time domain depends on the wave length and the length of the lattice elements. The higher frequencies lead to smaller wave lengths, which requires smaller elements sizes and increases the computational costs. The developed dynamic lattice method also considers the Sommerfeld radiation conditions for infinite domains by use of Perfectly Matched Layer (PML) [46] to absorb the waves' energy at the boundaries and avoid boundary reflections, if it is required.

Validation of the Dynamic Lattice Element method
The validation of the dynamic lattice method is performed by comparison to different benchmark tests. Firstly, the results obtained from the presented dynamic lattice method are compared with the solution of the Boundary Element Method (BEM) for an elastodynamic problem. The benchmark considered here is a plane P-SV-wave propagation within a plane strain domain. The BEM is solved according to the boundary integral equation method (BIEM) and with use of collocation procedure, whereas it can be written as [47]: Here, ω is the circular frequency; x, ξ are the coordinates of the source and receiver; U * , T * are the elastodynamic fundamental solutions for displacement and traction, respectively; t is the traction acting on the domain boundary Γ; and c is the jump-term which depends on the geometry of the source point.
The fundamental solutions for time-harmonic elasticity are as follows: where Ψ =K 0 (k 2 r) + 1 The notations C 1 and C 2 are the longitudinal and shear wave velocities, respectively; K n is the Bessel function of the second kind and n th order; k 1 = i(ω/C 1 ), k 2 = i(ω/C 2 ) with i = √ −1; and the subscripts l and k are the direction of the source load and the receiver response, respectively. Further description of the classical formulation can be found in [47]. After discretization, Eq.15 is rewritten in matrix notation as: where H and G are the influence matrices; and u, t are the vectors of nodal displacement and traction, respectively. Eq.18 is solved by: (1) assembling the influence matrices components corresponding to the unknown nodal values as matrix A in the left hand side; and (2) multiplying the prescribed boundary condition values with their corresponding influence matrix components and summing them in vector f in the right-hand side, which is given as: Here, x is a vector containing the unknown values. The geometry of the benchmark problem is a square with dimension of 6x6 m (Fig.5). The material properties are as follows: shear modulus G = 10 6 N/m 2 , density ρ=100 kg/m 2 , damping ratio 5 %, and Poisson's ratio ν=0.25. The prescribed boundary conditions are defined as: (1) zero perpendicular displacement on bottom boundary, (2) zero transverse displacement on the side boundaries, and (3) applied uniform traction of 100 N/m on top boundary. The BEM solution is obtained using quadratic line elements. An analytical solution of the resonance frequencies is given for discrete frequencies as: where L is the travel distance, C 1 is the P wave velocity and n = 0, 1, 2 ..., which yields ω n ≈ 45.40, 136.20, 227.00, 317.81 ... rad/s [47].  With increasing the mesh size of the lattice model, the obtained results as shown in Fig. 6 match the BIEM solution. The major factor that affects the final output of the solutions is the considered wavelength to elements length ratio [48]. When the ratio is kept higher than 10, the obtained results are accurate.

Simulation of wave fields in cracked and heterogeneous materials
For the analysis of wave fields in arbitrarily damaged and heterogeneous materials it is essential that the wave field is simulated while considering all the wave field effecting perturbations. The following essential factors can be categorized for the modeling: • Crack dependent criteria: length, thickness, orientation and location • Domain dependent criteria: characteristics and properties, such as stiffness, anisotropy and heterogeneity factors • Excitation dependent criteria: source, wavelength, frequency and magnitude • Model dependent criteria: element length size and mesh irregularity The effect of single crack orientation, multiple discontinuities, particle heterogeneity similar to concrete material and randomly distributed mineral heterogeneity as can be found in rock geomaterials on disturbance of displacement wave fields is investigated here. The Fig.7 depicts the generated cantilever beam element, where the receiver sensors are located on the outer surfaces (boundaries). The dynamic excitation is also carried out through these predefined receiver (reference) sensors. The arrival of first and second displacement wavefronts can be measured in each receiver sensor, which can be accessible to detect the discontinuities [49]. In all of the simulated results, the considered frequency of a single rectangular pulse is 0,2 MHz. The Young's modulus of a homogeneous medium is assumed to be equal to 3 GPa, where the Poisson's ratio is equal to 0.2. The randomness factor of the generated mesh is 0.5. The mesh size is kept constant and equal to 800x160, where the total number of generated nodes and lattice elements are equal to 128000 and 382081, respectively. Therefore, the minimum ratio of wavelength to element length is kept as low as 22. The maximum equivalent diameter of polygonal Voronoi cells is equal to 0.125 mm and the damping ratio is assumed to be equal to zero. The applied magnitude of the excitation pulse is equal to 10 N, which is low enough to avoid any frack occurrence and propagation. The non-linear dynamic frack simulation and propagation is not in the scope of this research work. The natural frequency of a simulated domain is dependent on the continuum effective stiffness and assigned masses on the Voronoi Cells. In each model setup, the same mesh with similar Voronoi masses are simulated. Therefore, the natural frequency of a system (cantilever beam) then only varies with the assigned different effective stiffness values.

Wave Fields in the Fractured Materials
In order to investigate the disturbance of wave fields in fractured materials, a cantilever beam body as shown in Fig.7 is considered. Initially, three homogeneous domains with different discontinuity conditions are simulated: • NC -homogeneous beam element with no damage • SC1 -homogeneous beam body with predefined single crack, where the crack orientation is perpendicular to the loading direction • SC2 -homogeneous beam body with predefined single crack, where the crack orientation builds an angle of 45 • with the horizon A single rectangular pulse is excited through receiver sensor R5, EX5. The simulated results according to dynamicLEM are provided in Fig.8 and Fig.9. Fig.8 illustrates the displacement wave field through the simulated beam body (qualitative results). In these results, the P-wave, SV-wave and Rayleigh surface waves are clearly visible and detectable. The time history of the displacements in the reference sensors (R 4 , R 5 , R 6 , R 13 , R 14 , and R 15 ) are plotted in Fig.9 (quantitative results). The disturbance of the wave fields due to the predefined discontinuities in SC1 and SC2 models are visible when compared to the NC model. Additionally, the arrival of the first and second wavefronts can be used to detect the location, orientation and length of the discontinuities in the domain. The wave field covers direct, diffracted and reflected waves around the crack within the domain. In these simulations, the time step (∆t) is equal to 0.0000001. The convergence of the dynamic solution depends on the size of the time steps.  In the next series of the results, multiple predefined discontinuities are generated inside the homogeneous beam element as shown in Fig.10. The ability of the proposed dynamic element method to image even complex wave field patterns is investigated here. In the presented setup, five cracks (MC1:MC5) of different locations, lengths and orientations are randomly generated (Fig.10). Three excitation positions, Ex2, Ex5 and Ex8 are used to excite the domain independently whereas the wave fields at the receivers can be analyzed. Similar to the previous results, qualitative and quantitative results of the simulations are presented in Fig.11 and Fig.12. Due to higher discontinuity of the simulated domain, a greater disturbance of wavefronts and the recorded time histories of the displacements are observed. For the identification of multiple discontinuities, analysis of the recorded time histories under multiple excitation sources are required. In the provided results, besides the excitation source, the recorded time history at R 11 , R 13 , R 14 , R 15 and R 17 sensors are plotted.

Wave Fields in Heterogeneous Materials
One of the main advantages of the lattice element method over conventional continuum methods is its ability to simulate the discontinuities while accounting for the inherent heterogeneity and irregularities in the particle scale. The irregularities such as the shape factors are already implemented while considering the irregular meshes and the defined randomness factor. In solid geomaterials, such as concrete or rock, the domain is composed of granular particles, cement (bond material) and the bond-particle interfaces as shown in Fig.2. Two different heterogeneous domains, one similar to concrete and the other one similar to rock geomaterial, are simulated here to analyse and study the effect of the heterogeneity ratio on wave disturbance. As a result, not only the effect of the heterogeneity ratio on the wave disturbance can be investigated, but also the mineral cluster effect can be studied. Similar to concrete composition, a particle packing procedure is implemented here to generate a heterogeneous domain composed of two main components: aggregates and cement matrix. A rectangular beam element with non-uniform packing and the particles diameter varying from 0.5∼4 mm is generated as shown in Fig.13. The same boundary condition as previous setups is considered here and the rectangular pulse with an amplitude of 10 N is excited from the R 5 receiver. The stiffness of aggregates (k p ), bond cement (k b ) and aggregate-bond interface (k i ) are then assigned for each corresponding lattice element. Five different heterogeneity (stiffness) ratios are simulated: Figure 13: The generated concrete beam structure composed of aggregates and cement The qualitative and quantitative results of DynamicLEM in a heterogeneous concrete body are presented in Fig.14 and Fig.15. Under the assumption of different heterogeneity ratios, the natural frequency of the simulated domains are also affected. According to the results, the wave fields exhibit a larger diffusion and noise with a heterogeneity ratio increment in the domain. Due to the increased wave velocity in CH2 and CH3, the simulation time step is reduced to ∆t = 0.00000001, however the frequency is kept constant and equal to 0.2 MHz. The higher the difference between the stiffness of the mediums, the greater the magnitude of the reflected wave fronts. Visualization and investigation of particle scale heterogeneity behaviour by ordinary continuum based methods is not possible, which opens a new research field to study these effects in detail by means of dynamic LEM. Next, for a medium with a randomly distributed heterogeneity, similar to what can be found in rock materials, a qualitative and quantitative analysis of wave field disturbance under multiple heterogeneity ratios are analysed. In this setup, a stochastic distribution of heterogeneity in each generated Voronoi cell is carried out. Therefore, instead of a cluster of cells representing a similar component as shown in the previous example, in this example, each Voronoi cell individually represents a unique mineral. The beam body without any predefined discontinuity is composed of four different minerals as shown with dark blue (DB), light blue (LB), red (R) and orange (O) cells in Fig.16. The interface values between two minerals (i and j) are determined based on the equivalent value: Here, E int is the Young's modulus of interface element between minerals i and j. This also highlights the importance of micro and nano-scale THM properties in heterogeneous material, which requires further investigations. Figure 16: The randomly distributed heterogeneity inside the beam body to mimic a rock geomaterial In order to investigate the effect of a heterogeneity ratio on wave fields and determine the critical heterogeneity ratio threshold, the following ratios are considered and simulated:  Fig.17 and Fig.18. In all of the simulations a similar setup with same mesh size (same masses), randomness factor and stochastic heterogeneity is considered. The adapted effective stiffnesses of domains results in different natural frequencies in each setup. The simulated damping ratio is equal to zero, which grants continuous vibration of the beam structure. Besides the simulation of wave propagation under isotropic conditions, the dynamic LEM is able to simulate an anisotropic media as well. Based on the presented results for both concrete and rock type geomaterials, it can be concluded that: • Increasing the stiffness heterogeneity between particle, bond matrix and aggregate-cement interface induces excessive disruption on the wave fronts. The magnitude of the reflected wave fronts is increased when the the ratio of the stiffness is enlarged. The recorded and plotted data are dependent on the natural frequency of a beam structure. In all of the setups presented for concrete body, the assigned masses (generated mesh) on Voronoi cells are constant. With increasing the stiffness of a domain, the wave velocity also becomes greater. Although the increment of the heterogeneity induced wave dispersion in the concrete body, the evaluation of the outcomes based on the arrival of P and SV waves is still possible.
• In a rock domain, the effect of heterogeneity ratios as well as number of minerals on wave disturbance is analysed. Similar to the concrete body, increasing the stiffness heterogeneity induces excessive disruption on the wave fronts. In contrast to the concrete body, after a certain stiffness heterogeneity ratio in rock domain, the evaluation of the outcomes based on the arrival of P and SV waves is not possible. This is clearly visible from the plotted time histories, starting from RH5, where at around 0.9 × 10 −4 the accumulated noises disrupt the wave fronts.
• The total number of minerals in a rock domain has an effect on the disturbance of the wave fronts. This is clearly visible when comparing the RH8 and RH9 results, where the maximum stiffness heterogeneity ratio is equal to 200. In RH8, the disturbance of wave fronts starts from 0.45 × 10 −4 , where in RH9, the wave dispersion begins at 0.4 × 10 −4 . It should be noted that the natural frequencies of these domains are not equal, which can also induce different vibration responses. The results clearly depict that increasing the number of minerals induces larger disturbance on the wavefronts.
• While comparing the rock and concrete bodies it is clear that in heterogeneous rock body in a same stiffness heterogeneity ratio, the wave dispersion and scattering is higher. In contrast to CH4, in RH9, where in both setups the maximum stiffness heterogeneity ratio is equal to 200, the evaluation of the arrival of P and SV waves is not possible. This again emphasises the importance of considering the inherent heterogeneity in numerical simulations. Identification of discontinuities based on arrival of first and second P and S waves will be then a challenging task, which requires further investigations. For a similar max stiffness ratios, in cluster mineral distribution (concrete) the heterogeneity of a domain is lower than individual mineral distribution (Rock).
• A small stiffness transition zone (mineral to mineral) produces a larger variation of P and S wave velocities, which leads to larger reflections and dispersion of waves. Increasing the thickness of the transition zone (cluster of minerals) results in gradual variation of seismic velocities. This also reduces the energy of reflected and diffracted waves, which eventually has less impact on the disturbance of wave fronts.

Discussion and Conclusion
The application of a 2D dynamic lattice in the simulation of displacement wave fields in the fracked and heterogeneous medium is studied. The provided static dynamic lattice result for a compact tension specimen (CTS) emphasises the importance of consideration of irregularities in the particle scale, as the fracking path and post failure behavior are affected by the inter-particle stress distribution. The proposed dynamic lattice model based on the general equation of motion is not only used to simulate the displacement wave fields in homogeneous domain but also to illustrate qualitatively and quantitatively the disturbance of wave fields in the discontinuous and heterogeneous medium. It is also possible to detect the defects using the dynamic lattice model, in which a high-frequency wave can be generated and recorded in the reference sensors. The dynamic lattice is also capable of detecting the location of the initiated crack in a stressed structure. In order to validate the in-house developed 2D dynamic lattice, the displacement vs. angular frequency results of a 2D plain strain domain are compared to the boundary element method (BEM) solution. It is shown that with larger wavelength to element length ratios, the lattice results provide great accuracy.
The advantage of dynamic lattice in the displacement wave field simulation lies in its embedded irregularity in domain discretization, its ability to define discontinuities and particle heterogeneity, simulate frack initiation and propagation as well as the stress redistribution and concentration upon crack propagation. Similar to discrete methods, the lattice method also considers the inherent heterogeneity in particle scale. To better illustrate this advantage, the developed dynamicLEM is used to simulate the displacement wave fields in heterogeneous concrete and rock geomaterials.
The results indicate that with increasing the maximum heterogeneity ratio (stiffness ratio), the disturbance of wave fields becomes greater. After a certain heterogeneity threshold, the quantitative evaluation of dynamic results with conventional methods is not possible. It is also shown that increasing the number of mineral compositions in a rock body has a considerable effect on the wave field disturbance. Comparing the rock and concrete beam bodies indicates the importance of mineral distribution scheme and particle sizes on dispersion of wavefronts. It is shown that in a same stiffness heterogeneity ratio, the reflection and diffraction of waves in rock domain (individual mineral distribution) are stronger than concrete domain, where a cluster of Voronoi cells represent an individual aggregate. The results coheres with the theory, where a small stiffness transition zone (mineral to mineral) produces a larger variation of P and SV wave velocities compared to a larger transition zone (cluster of minerals). The simulations are performed on a Desktop-Pc with a Xeon processor (2.10 GHz) with a total number of 16 cores and the computational time for a single simulation is approximately 24 hours.
The dynamic lattice element is capable of modeling the wave field's disturbance and dispersion when a discontinuity and heterogeneity in the domain is present. The identification of discontinuities in heterogeneous materials is a challenging task, where particle scale models similar to dynamicLEM can provide accurate outcomes. However, the simulation of large domains with high frequencies requires a substantial number of elements, not only to model the particle size effect but also to reach greater wavelength to elements length ratio and avoid numerical difficulties. One way to overcome this obstacle is to consider the artificial neural networks (ANN) models to predict the location of the discontinuities in the homogeneous and heterogeneous bodies [37]. The recorded displacement wave fields in each receiver sensor is considered as training data to develop an ANN method to not only decrease the computational costs but also increase the accuracy of the crack predictions. Eventually, for any given wave spectrum that can be obtained from the field applications, the ANN method is able to predict the location, length and orientation of the existing discontinuities.