Peri-Elastodynamic Simulations of Guided Ultrasonic Waves in Plate-Like Structure with Surface Mounted PZT

Peridynamic based elastodynamic computation tool named Peri-elastodynamics is proposed herein to simulate the three-dimensional (3D) Lamb wave modes in materials for the first time. Peri-elastodynamics is a nonlocal meshless approach which is a scale-independent generalized technique to visualize the acoustic and ultrasonic waves in plate-like structure, micro-electro-mechanical systems (MEMS) and nanodevices for their respective characterization. In this article, the characteristics of the fundamental Lamb wave modes are simulated in a sample plate-like structure. Lamb wave modes are generated using a surface mounted piezoelectric (PZT) transducer which is actuated from the top surface. The proposed generalized Peri-elastodynamics method is not only capable of simulating two dimensional (2D) in plane wave under plane strain condition formulated previously but also capable of accurately simulating the out of plane Symmetric and Antisymmetric Lamb wave modes in plate like structures in 3D. For structural health monitoring (SHM) of plate-like structures and nondestructive evaluation (NDE) of MEMS devices, it is necessary to simulate the 3D wave-damage interaction scenarios and visualize the different wave features due to damages. Hence, in addition, to simulating the guided ultrasonic wave modes in pristine material, Lamb waves were also simulated in a damaged plate. The accuracy of the proposed technique is verified by comparing the modes generated in the plate and the mode shapes across the thickness of the plate with theoretical wave analysis.


Introduction
In Structural Health Monitoring (SHM) research, Lamb waves are widely used for damage detection in the metallic plate-like structures [1,2]. High-frequency ultrasonic actuators and sensors are strategically mounted on the plate-like structure to detect, localize and characterize the damages [3]. Symmetric (S 0 ) and Antisymmetric (A 0 ) Rayleigh-Lamb wave modes while travel through the plate, interacts with the boundaries and the discontinuities [4] and are subjected to mode conversion. Efficient diagnostic and prognostic algorithms are then employed to estimate the severity of the damage and the damage growth.
Sensor signals play a critical role in quantifying the extent of damage within the structure. In most practical cases with SHM, the damage state of the material is unknown and the sensor signals are the observables. There could be infinite possibilities of damage states in the material and it is impossible to experimentally obtain the understanding of the sensor signals due to the varying was used to displace the material points located at one end of the plate to generate the surface skimming Rayleigh wave. Only recently, Hafezi et al. [28,33,34] employed PD theory to model the in plane longitudinal ultrasonic wave in an aluminum plate [17,29]. In these works, the elastic wave propagation was simulated by considering only one layer of the material points without out of plane deformation. Hence, the method is inefficient and insufficient to simulate the Lamb wave modes. Additionally, comparing the peridynamic and continuum formulation of strain energy density, the mathematical equation for the bond constant used in these articles [28,33,34] is incorrect, which is correctly presented in this article. Studies reported were neither validated with any analytical solution nor compared with any experimental results. In this article, it is argued with corrected bond constant that one layer of the material point is not only insufficient but also inaccurate to simulate the Lamb wave modes. Monolayer simulation cannot accommodate the out of plane deformation. Therefore, in this article, it is shown that with specific discretization schemes, at least three layers of spatial material points are required to correctly capture the fundamental Lamb wave modes (S 0 and A 0 ) in a plate.
In addition to the SHM of plate structure, NDE is a crucial step for the ensured functionality and fail-safe design of specific plate-like MEMS devices. Such virtual NDE experimental scenarios for MEMS devices could be simulated using Peri-elastodynamic as proposed in this article. From here onwards in this article, the plate-like structure for SHM and the plate-like structure in MEMS devices for their characterization are commonly termed as plate-like structure or material.

Basics of Bond-Based Peridynamic Formulation
Peridynamic is a meshless simulation method where a material body is discretized into a series of material points. To illustrate the kinetics, an undeformed and a deformed state of two particles are shown in Figure 1a,b, respectively. The deformation between two material points produces a pairwise interaction along the bond (Figure 1b). The equation of equilibrium at the material point x at time t can be written as follows [26], where, u, , b, f and V are the displacement, density, body force per unit volume at the material point, the pairwise force acting along the bond between x and x and the volume of the material point at x , respectively. The integral of the forces acting at the parent material point x is performed over a finite region H, called Horizon in the peridynamic theory. The material points inside a Horizon are called the family members of the parent material point. The parent material point interacts with the other material points within its family while the interactions with the points outside the Horizon are negligible. Figure 1c,d represent a family of a material point in three-dimension and in two-dimension, respectively. Relative distance and the displacement between the two material points in the reference configuration ( Figure 1a) can be viz., Relative displacement between the two material points in the deformed configuration can be expressed as (Figure 1), (ξ + η) = (u(x , t) + x ) − (u(x, t) + x) (4) Constitutive law in the peridynamic approach is expressed as follows, where, c and s are the bond constant and the stretch of the bond, respectively. The bond constant for a two-dimensional material body is calculated by balancing the strain energy density from the continuum mechanics and peridynamic formulation viz. [32], where E, h and δ are the Young's Modulus, the half thickness of the plate and the radius of the Horizon H, respectively. In peridynamics, the stretch is the ratio of the change in the length of the peridynamic bond due to the deformation of the initial bond length, expressed as follows, Sensors 2018 , 18, 274  4 of 17 where, c and s are the bond constant and the stretch of the bond, respectively. The bond constant for a two-dimensional material body is calculated by balancing the strain energy density from the continuum mechanics and peridynamic formulation viz. [32], where E, h and δ are the Young's Modulus, the half thickness of the plate and the radius of the Horizon H, respectively. In peridynamics, the stretch is the ratio of the change in the length of the peridynamic bond due to the deformation of the initial bond length, expressed as follows,

Actuator Modelling with in Plane Excitation
In SHM, PZT actuator is attached to the host structure as shown in Figure 2a. In this article, a 300 mm × 200 mm plate is used for the wave propagation simulation (Figure 2a). A plate with a hole is also considered to study the wave damage interactions (Figure 2b). Lamb waves are generated by applying a standard tone-burst voltage signal to the PZT actuator [3]. Except the interface between the PZT actuator and the plate, all other boundaries of the plate are considered stress-free in the simulation. The voltage signal actuates the PZT and transformed the energy into in plane mechanical strain. The in-plane strain causes the rapid localized displacement in the host structure, which results in the Lamb wave propagation in the plate. The Lamb wave modes create out of plane displacements and hence, to accommodate such deformation three layers of material points are used in the modeling (Figure 2c). Figure 2d shows the discretization used in this study which is further described in Section 2.4. Displacement in the structure varies linearly along the length of the PZT and attains a maximum value at the boundaries [37] as shown in Figure 2e. In this study, a square PZT with a dimension of 2.4 × 2.4 mm is modeled by applying a maximum of 1 μ m in plane radial displacement to the circumferential material points, shown in the Figure 2e. The displacements of the material points at the center of the PZT were enforced to zero. The equation for variation of displacement due to the application of a tone burst signal is expressed as.

Actuator Modelling with in Plane Excitation
In SHM, PZT actuator is attached to the host structure as shown in Figure 2a. In this article, a 300 mm × 200 mm plate is used for the wave propagation simulation (Figure 2a). A plate with a hole is also considered to study the wave damage interactions (Figure 2b). Lamb waves are generated by applying a standard tone-burst voltage signal to the PZT actuator [3]. Except the interface between the PZT actuator and the plate, all other boundaries of the plate are considered stress-free in the simulation. The voltage signal actuates the PZT and transformed the energy into in plane mechanical strain. The in-plane strain causes the rapid localized displacement in the host structure, which results in the Lamb wave propagation in the plate. The Lamb wave modes create out of plane displacements and hence, to accommodate such deformation three layers of material points are used in the modeling (Figure 2c). Figure 2d shows the discretization used in this study which is further described in Section 2.4. Displacement in the structure varies linearly along the length of the PZT and attains a maximum value at the boundaries [37] as shown in Figure 2e. In this study, a square PZT with a dimension of 2.4 × 2.4 mm is modeled by applying a maximum of 1 µm in plane radial displacement to the circumferential material points, shown in the Figure 2e. The displacements of the material points at the center of the PZT were enforced to zero. The equation for variation of displacement due to the application of a tone burst signal is expressed as.
u(x, t) = U(x)e −pt 2 /2 sin(ω c t) (8)  where, ω c and U(x) are the central frequency and the maximum displacement amplitude of the excitation signal given to different material points, respectively. The parameter p in the Equation (8) is expressed as p = (2khω c /N c ) 2 (9) where, k, h and N c are the signal shape factor, the half thickness of the specimen and the number of cycles of the actuation signal, respectively. To select the desired excitation frequency, the tuning curves for S 0 and A 0 modes in the Aluminum 6061-T6 were obtained from an open source software "Waveform Revealer" [38], developed by the LAMSS laboratory at the University of South Carolina (USC). As shown in Figure 3a,b, the central frequency of 150 kHz is chosen to make sure that only the S 0 and A 0 modes are excited when the modal amplitudes are comparable but nonequal. In the present study, a 3.5 count tone-burst signal with the central frequency (ω c ) of 150 kHz is used. Figure 3c,d show the time domain signal and its frequency content, respectively. where, ω c and U(x) are the central frequency and the maximum displacement amplitude of the excitation signal given to different material points, respectively. The parameter p in the Equation where, k , h and N c are the signal shape factor, the half thickness of the specimen and the number of cycles of the actuation signal, respectively. To select the desired excitation frequency, the tuning curves for S0 and A0 modes in the Aluminum 6061-T6 were obtained from an open source software "Waveform Revealer" [38], developed by the LAMSS laboratory at the University of South Carolina (USC). As shown in Figure 3a,b, the central frequency of 150 kHz is chosen to make sure that only the S0 and A0 modes are excited when the modal amplitudes are comparable but nonequal. In the present study, a 3.5 count tone-burst signal with the central frequency ( c ω ) of 150 kHz is used.    where, ω c and U(x) are the central frequency and the maximum displacement amplitude of the excitation signal given to different material points, respectively. The parameter p in the Equation where, k , h and N c are the signal shape factor, the half thickness of the specimen and the number of cycles of the actuation signal, respectively. To select the desired excitation frequency, the tuning curves for S0 and A0 modes in the Aluminum 6061-T6 were obtained from an open source software "Waveform Revealer" [38], developed by the LAMSS laboratory at the University of South Carolina (USC). As shown in Figure 3a,b, the central frequency of 150 kHz is chosen to make sure that only the S0 and A0 modes are excited when the modal amplitudes are comparable but nonequal. In the present study, a 3.5 count tone-burst signal with the central frequency ( c ω ) of 150 kHz is used.

Numerical Time Integration
In Peri-elastodynamics approach, the plate is spatially discretized into a finite number of material points. Each material point has finite volume in the reference configuration. Material volume for 3D uniform discretization grid is calculated as dl 3 , where dl is the element length [32,39]. By replacing the integration with a finite summation over all the material points inside the Horizon, the equation of motion at material point i after the time step n can be expressed as, The net force (f) acting on a material point is calculated by summing the peridynamic forces on the parent material point due to all the neighboring points inside its Horizon. N f represents the number of material points within the Horizon enclosing the parent material point i. The Velocity-Verlet integration [31] scheme is employed in this study to calculate the displacements in the time domain for given boundary and initial conditions as follows, Stability of the numerical solution can be obtained for a small-time step ∆t and a spatial discretization step ∆s. To have a convergence of the displacements, the detailed procedure to select the time step (∆t) and spatial discretization (∆s) is discussed in Section 2.4.

Peri-Elastodynamic Spatial and Temporal Discretization
Proper spatial and temporal discretization are the critical parameters for the convergence of the solution in wave propagation simulation. Maximum spatial discretization (∆s) must meet the criterion below [40], λ min = c min f ; ∆s = λ min /10 (12) where, λ min is the minimum wavelength of the Lamb wave modes and c min is the minimum phase velocity of the simulated modes at the excitation frequency f. Phase velocity at the excitation frequency can be easily obtained from the theoretical dispersion curves. In this work, the phase velocity of the A 0 mode at 150 kHz is used as c min in Equation (12). The Courant-Friedrichs-Levy condition is used to obtain a numerically stable time step (∆t) [40].
where, c max is the maximum phase wave velocity of the propagating modes. In this work, the phase velocity of the S 0 mode at 150 kHz is the c max in Equation (13). Thus, to obtain a converging solution, the spatial and the temporal step sizes are chosen to be 1.2 mm and 0.01 µS, respectively, satisfying the Equations (12) and (13). Material points are chosen in a grid fashion with a spacing of ∆s to model the plate with a layer spacing of 1 mm between each layer L1, L2 and L3. 41,750 material points were used in each layer in the pristine plate. A total of 125,250 material points was used in the simulation including all the three layers L1, L2 and L3. In case of the damaged plate (with hole), there was 41,610 material points in each layer and a total of 124,830 material points were used for the Peri-elastodynamic simulation. Each parent point is assigned with a family based on its Horizon and bonds were established between each pair of material points within the family. A 3.015∆S was used as the Horizon size in the simulation. To model a plate with a through-thickness hole, a pristine discretization is performed and then the material points are removed from the geometry to produce the hole.

Lamb Wave Dispersion Relation
Dispersion curves of various Lamb wave modes are used to predict the existence of various modes at a particular excitation frequency [3]. Two types of Lamb wave modes exist in a plate based on the particle motion, named Symmetric modes (i.e., S 0 , S 1 , S 2 . . . ) and Antisymmetric modes (i.e., A 0 , A 1 , A 2 . . . ). Generation of the Lamb wave modes in a plate depends on the frequency of excitation, the thickness of the plate and the material properties (Density, Young's modulus or Shear Modulus and Poisson's ratio) of the material. Dispersion of various Lamb wave modes is obtained by solving Rayleigh-Lamb wave equations. Rayleigh-Lamb wave equation for symmetric Lamb wave modes is expressed by, For antisymmetric modes, equation is written as follows, Parameters in the above equations are expressed as, where ω, k, c s , c L , ν, µ, and H are the angular frequency, wavenumber, the shear wave velocity, the longitudinal wave velocity, the poison's ratio, the shear modulus, the density and the half thickness of the plate, respectively. In this work, the dispersion curves for various Lamb wave modes in an Aluminum 6061-T6 plate were calculated using the commercially available 'Disperse' software [41], designed by the Imperial College, London, UK, as shown in Figure 3a. The plate thickness was set to 2 mm and the material properties were set to the values listed in Table 1.

Numerical Computation and Results
The Peri-elastodynamic simulations were performed on a workstation with two Intel Xeon (R) CPU E5-2650 V3 2.30 GHz processors with total 128 GB RAM, in a single core. One simulation was completed within a reasonable time of~46.5 h. Note that this problem is highly parallelizable and can be implemented with distributed clusters, GPUs, multiple threads, or a combination of these methods. Preliminary translation of this MATLAB program into multithreaded C++ resulted in a 20 times speedup.

Lamb Wave Propagation in the Pristine Plate
In this article, fundamental Lamb wave modes (S 0 and A 0 ) are simulated which are widely used in the damage detection with ultrasonic SHM. Figure 4 shows the out of plane displacements u z (x, y, t) and, in plane displacements u x (x, y, t) and u y (x, y, t) for the Lamb wave propagation plotted at t = 20,  40 and 60 µS, respectively. Each stack of the figure is plotted with the respective displacement pattern in the layers L1, L2 and L3. It can be seen that the S 0 mode travels faster than the A 0 mode and confirms the dispersion curves. In wave propagation simulation, most inaccuracy comes from the boundary reflections if the results are not converged. The best approach to judge if the simulation of the wave propagation is converged is to evaluate the boundary reflections. In the present simulation, the reflected modes from the plate boundaries are clearly visible in the figures. Usually, there is a possibility of divergence of the solution at the boundaries due to the nonconvergence of the solution. However, with the specific steps with the Peri-elastodynamic process described in this paper, the results will be converged and bounded boundary reflections will be achieved.

Lamb Wave Propagation in the Pristine Plate
In this article, fundamental Lamb wave modes (S0 and A0) are simulated which are widely used in the damage detection with ultrasonic SHM. Figure 4 shows the out of plane displacements u (x,y,t) z and, in plane displacements u (x,y,t) x and u (x,y,t) y for the Lamb wave propagation plotted at t = 20, 40 and 60 μS, respectively. Each stack of the figure is plotted with the respective displacement pattern in the layers L1, L2 and L3. It can be seen that the S0 mode travels faster than the A0 mode and confirms the dispersion curves. In wave propagation simulation, most inaccuracy comes from the boundary reflections if the results are not converged. The best approach to judge if the simulation of the wave propagation is converged is to evaluate the boundary reflections. In the present simulation, the reflected modes from the plate boundaries are clearly visible in the figures. Usually, there is a possibility of divergence of the solution at the boundaries due to the nonconvergence of the solution. However, with the specific steps with the Peri-elastodynamic process described in this paper, the results will be converged and bounded boundary reflections will be achieved.  In the top (L1) and the bottom (L3) layers in Figure 4(a-1-a-3) and Figure 4(b-1-b-3), both Symmetric and Antisymmetric modes (S 0 and A 0 ) are visible in the wave fields composed of in plane displacements. The contribution of the A 0 mode in the in-plane motion at the middle layer (L2) of the plate is negligible but the out of plane motion is dominant. As shown in Figure 4(c-1-c-3), the out of plane displacement of the A 0 mode is visible (i.e., contribution from u z (x, y, t)) while the S 0 mode is barely noticeable. This is because the displacements of the in-plane particles of the S 0 mode dominates over the out of plane motion of the particles. Also, the contribution of the S 0 mode in the Time-space representation of the in plane (u x (x, t)) and the out of plane displacement (u z (x, t)) are presented in Figure 5. Displacements at the top, middle and the bottom surfaces (L1, L2 and L3), are presented to investigate the existence of the different Lamb wave modes and their contribution to the displacement in each layer. In Figure 5(a-1-a-3), it is observed that the S 0 mode contributes to the in-plane displacement in all layers and the A 0 mode contributes only to the top and the bottom layers. Time-space representation of the out of plane displacement is shown in Figure 5(b-1-b-3), which show that the A 0 mode had a higher amplitude than the S 0 mode. A 0 mode contributed to the displacement of all the layers (L1, L2 and L3), whereas, S 0 mode contributed only to the top (L1) and the bottom (L3) layers. Video S1 in the supplementary material shows the guided wave propagation through the u z (x, y, t) displacement of the layer L1 of the pristine state of the plate used in this article.
Symmetric and Antisymmetric modes (S0 and A0) are visible in the wave fields composed of in plane displacements. The contribution of the A0 mode in the in-plane motion at the middle layer (L2) of the plate is negligible but the out of plane motion is dominant. As shown in Figure 4(c-1-c-3), the out of plane displacement of the A0 mode is visible (i.e., contribution from u (x,y,t) z ) while the S0 mode is barely noticeable. This is because the displacements of the in-plane particles of the S0 mode dominates over the out of plane motion of the particles. Also, the contribution of the S0 mode in the displacement of the middle layer (L2) of the plate is negligible, because, in the S0 mode, the middle layer (L2) remains undisturbed.
Time-space representation of the in plane ( u (x,t) x ) and the out of plane displacement ( u (x,t) z ) are presented in Figure 5. Displacements at the top, middle and the bottom surfaces (L1, L2 and L3), are presented to investigate the existence of the different Lamb wave modes and their contribution to the displacement in each layer. In Figure 5(a-1-a-3), it is observed that the S0 mode contributes to the inplane displacement in all layers and the A0 mode contributes only to the top and the bottom layers. Time-space representation of the out of plane displacement is shown in Figure 5(b-1-b-3), which show that the A0 mode had a higher amplitude than the S0 mode. A0 mode contributed to the displacement of all the layers (L1, L2 and L3), whereas, S0 mode contributed only to the top (L1) and the bottom (L3) layers. Video S1 in the supplementary material shows the guided wave propagation through the u (x,y,t) z displacement of the layer L1 of the pristine state of the plate used in this article.

Vector Field Representation of the Lamb Wave Modes
To prove the accuracy of the Peri-elastodynamic simulation, the characteristics of the Lamb wave modes (S0 and A0), in plane and out of plane particle motion across-the-thickness are plotted in

Vector Field Representation of the Lamb Wave Modes
To prove the accuracy of the Peri-elastodynamic simulation, the characteristics of the Lamb wave modes (S 0 and A 0 ), in plane and out of plane particle motion across-the-thickness are plotted in Figure 6, after 41 µS. Out of plane (u z (x, y, t)) and in plane (u x (x, y, t)) displacement distribution of A 0 mode are extracted along the cross-sections C 1 − C 2 and C 1 − C 2 , of the plate. Similarly, the out of plane and the in-plane displacement distribution of S 0 mode is plotted along the cross-section lines D 1 − D 2 and D 1 − D 2 , respectively. Vector fields and displacement distributions are shown in Figure 6a-d. It is observed in Figure 6a that all particles moved either upwards (+Z) or downwards (−Z) with variable amplitude (like bending motion) due to the generation of the A 0 mode. In Figure 6c, the top and the bottom layers are symmetrically displaced with respect to the mid-plane and the displacement of the mid-plane is almost zero due to the generation of the S 0 mode. In plane particle motion in A 0 and S 0 modes are also shown in Figure 6b,d. In Figure 6b, the particles at the top and the bottom layers are moved in the opposite directions along the in-plane direction and the displacement of the mid-plane is zero due to the generation of the A 0 mode. In Figure 6c, the particle displacements are constant across the thickness due to the S 0 mode. Vector fields and the mode shapes in Figure 6 indicate that the Peri-elastodynamics simulated the Lamb waves accurately. displacement of the mid-plane is almost zero due to the generation of the S0 mode. In plane particle motion in A0 and S0 modes are also shown in Figure 6b,d. In Figure 6b, the particles at the top and the bottom layers are moved in the opposite directions along the in-plane direction and the displacement of the mid-plane is zero due to the generation of the A0 mode. In Figure 6c, the particle displacements are constant across the thickness due to the S0 mode. Vector fields and the mode shapes in Figure 6 indicate that the Peri-elastodynamics simulated the Lamb waves accurately. Figure 6. Peri-elastodynamics simulation, vector field and displacement distribution of the S0 and A0 modes across the thickness of the plate: (a) Vector field of the A0 mode for out of plane motion, (b) Vector field of the A0 mode for in plane motion, (c) Vector field of the S0 mode for out of plane motion, (d) Vector field of the S0 mode for in plane motion.

Frequency-Wavenumber Analysis: Verification of the Simulation Results
Multidimensional Fourier transform is widely used to separate the different Lamb wave modes [42]. Two-dimensional and three-dimensional Fourier transforms (2D or 3D FFT) are performed on space and time domain data. The equation that transforms the time-space wavefield data into the frequency-wavenumber representation of the wave field and can be expressed as follows [42], where, take the values, 1 and 2, take values 1, 2 and 3, ( , ) designates the -th displacement after time at the point located on the 2D x-y plane. , designates the -th displacement at frequency in reciprocal space of wavenumbers at the point located on the 2D reciprocal − plane. Here index, 1, 2 and 3 stands for the coordinate x, y and z. In multi-modal wave propagation analysis, distinguishing the different modes from a time domain signal is difficult, especially on a small plate where the wave modes tend to overlap. In this work, frequency-wavenumber plots are presented to visualize the different modes separately. This is also verified by comparing the simulated dispersion results with the theoretical dispersion curves.

Frequency-Wavenumber Analysis: Verification of the Simulation Results
Multidimensional Fourier transform is widely used to separate the different Lamb wave modes [42]. Two-dimensional and three-dimensional Fourier transforms (2D or 3D FFT) are performed on space and time domain data. The equation that transforms the time-space wavefield data into the frequency-wavenumber representation of the wave field and can be expressed as follows [42], where, j take the values, 1 and 2, p take values 1, 2 and 3, u p x j, t designates the p-th displacement after time t at the point x j located on the 2D x-y plane. u p k j , ω designates the p-th displacement at frequency ω in reciprocal space of wavenumbers at the point k j located on the 2D reciprocal k x − k y plane. Here index, 1, 2 and 3 stands for the coordinate x, y and z. In multi-modal wave propagation analysis, distinguishing the different modes from a time domain signal is difficult, especially on a small plate where the wave modes tend to overlap. In this work, frequency-wavenumber plots are presented to visualize the different modes separately. This is also verified by comparing the simulated dispersion results with the theoretical dispersion curves. For this purpose, 2D and 3D Fast Fourier Transforms (FFT) were performed on the simulated displacement wave field to obtain the frequency-wavenumber representations.
To perform the 2D-FFT, out of plane (u z (x, y, t)) and in plane (u x (x, y, t)), displacement data are obtained across-the-thickness of the plate along the selected red dotted line shown in Figure 4(a-1). 163 spatial points with a resolution of 1.2 mm along the red line shown in Figure 4(a-1) were used in the analysis. Matrix size used to store the displacements wave field was 163 × 8000. Note that, the 2D FFT was performed on the displacement vectors obtained from all the three material layers (L1, L2 and L3).
Frequency-wavenumber domain representation of the in-plane displacement (u x (x, y, t)) is depicted in Figure 7(a-1-a-3), respectively. Both the S 0 and A 0 modes are identified at the top and the bottom material layers. This is because they both significantly contributed to the energy of the in-plane wave motion. The amplitudes of the A 0 mode are slightly greater than that of the S 0 mode. A similar phenomenon is predicted from the tuning curve of the plate at 150 kHz. The contribution of the in-plane motion of the A 0 mode to the energy of the middle layer is almost zero. Similarly, Figure 7(b-1-b-3) are obtained from the wavenumber-frequency domain representation of the out of plane (u z (x, y, t)) displacements at the top (L1), middle (L2) and the bottom layer (L3), respectively. The energy distribution of the A 0 mode is higher than the S 0 mode at all the material layers. The S 0 mode is visible only at the top and the bottom layers of very low amplitude. This is because the out of plane motion of the particles in S 0 mode is very low and the displacement at the midplane is almost zero.
To perform the 2D-FFT, out of plane ( u (x,y,t) z ) and in plane ( u (x,y,t) x ), displacement data are obtained across-the-thickness of the plate along the selected red dotted line shown in Figure 4(a-1). 163 spatial points with a resolution of 1.2 mm along the red line shown in Figure 4(a-1) were used in the analysis. Matrix size used to store the displacements wave field was 163 × 8000. Note that, the 2D FFT was performed on the displacement vectors obtained from all the three material layers (L1, L2 and L3).
Frequency-wavenumber domain representation of the in-plane displacement ( u (x,y,t) x ) is depicted in Figure 7(a-1-a-3), respectively. Both the S0 and A0 modes are identified at the top and the bottom material layers. This is because they both significantly contributed to the energy of the in-plane wave motion. The amplitudes of the A0 mode are slightly greater than that of the S0 mode. A similar phenomenon is predicted from the tuning curve of the plate at 150 kHz. The contribution of the in-plane motion of the A0 mode to the energy of the middle layer is almost zero. Similarly, Figure 7(b-1-b-3) are obtained from the wavenumber-frequency domain representation of the out of plane ( u (x,y,t) z ) displacements at the top (L1), middle (L2) and the bottom layer (L3), respectively. The energy distribution of the A0 mode is higher than the S0 mode at all the material layers. The S0 mode is visible only at the top and the bottom layers of very low amplitude. This is because the out of plane motion of the particles in S0 mode is very low and the displacement at the midplane is almost zero. Next, the 3D-FFT is employed to transform the 3D displacement data ( u (x,y,t) x , u (x,y,t) y and u (x,y,t) z ) into the frequency-wavenumber domain ( u (k ,k ,ω) x x y , u (k ,k ,ω) y x y and u (k ,k ,ω) z x y ) and are shown in Figure 8a-c, respectively. In this work, frequency transformation is Next, the 3D-FFT is employed to transform the 3D displacement data (u x (x, y, t), u y (x, y, t) and u z (x, y, t)) into the frequency-wavenumber domain (u x (k x , k y , ω), u y (k x , k y , ω) and u z (k x , k y , ω)) and are shown in Figure 8a-c, respectively. In this work, frequency transformation is performed only on the data obtained from the top surface (L1). The size of the matrix used to store the 3D displacement data was 163 × 250 × 8000. Wavenumber plots at the frequencies 110 kHz, 150 kHz, 185 kHz and 225 kHz, are presented in the Figures. It is seen that both the S 0 and the A 0 modes appeared in the form of two concentric circular rings. The radius of the circles corresponds to the wave numbers at the respective frequencies. Wavenumbers of the S 0 and A 0 modes at the 150 kHz are obtained from the Peri-elastodynamic simulation and are 0.56 rad/mm and 0.187 rad/mm, respectively. A smaller circle corresponds to the S 0 while the larger corresponds to the A 0 mode. It is also observed that the energy of the modes at the frequencies 110 kHz, 185 kHz and 225 kHz are lower compared to that of at the 150 kHz. This is because most of the energy of the modes is concentrated around the excitation frequency (150 kHz).
the wave numbers at the respective frequencies. Wavenumbers of the S0 and A0 modes at the 150 kHz are obtained from the Peri-elastodynamic simulation and are 0.56 rad/mm and 0.187 rad/mm, respectively. A smaller circle corresponds to the S0 while the larger corresponds to the A0 mode. It is also observed that the energy of the modes at the frequencies 110 kHz, 185 kHz and 225 kHz are lower compared to that of at the 150 kHz. This is because most of the energy of the modes is concentrated around the excitation frequency (150 kHz). To verify the directional dependency of the Lamb wave propagation, 2D wavenumber plots (u k , k , u k , k and u k , k ) at ω = 150 kHz, are obtained from the 3D FFT and were compared with those obtained from the theoretical predictions using Disperse software at 150 kHz. Theoretical wavenumber plot is superimposed on the numerically obtained wavenumber plots in Figure 9a-c. Good agreements between the numerical and analytical results are obtained. It shows that the Peri-elastodynamic can predict the dispersion relation of fundamental Lamb wave modes accurately in all directions. Therefore, the strong evidences discussed above demonstrate that the Peri-elastodynamics would be a potential tool to effectively simulate the Lamb wave modes for NDE and SHM applications.  To verify the directional dependency of the Lamb wave propagation, 2D wavenumber plots (u x k x , k y , u y k x , k y and u z k x , k y ) at ω c = 150 kHz, are obtained from the 3D FFT and were compared with those obtained from the theoretical predictions using Disperse software at 150 kHz. Theoretical wavenumber plot is superimposed on the numerically obtained wavenumber plots in Figure 9a-c. Good agreements between the numerical and analytical results are obtained. It shows that the Peri-elastodynamic can predict the dispersion relation of fundamental Lamb wave modes accurately in all directions. Therefore, the strong evidences discussed above demonstrate that the Peri-elastodynamics would be a potential tool to effectively simulate the Lamb wave modes for NDE and SHM applications.
are obtained from the Peri-elastodynamic simulation and are 0.56 rad/mm and 0.187 rad/mm, respectively. A smaller circle corresponds to the S0 while the larger corresponds to the A0 mode. It is also observed that the energy of the modes at the frequencies 110 kHz, 185 kHz and 225 kHz are lower compared to that of at the 150 kHz. This is because most of the energy of the modes is concentrated around the excitation frequency (150 kHz). To verify the directional dependency of the Lamb wave propagation, 2D wavenumber plots (u k , k , u k , k and u k , k ) at ω = 150 kHz, are obtained from the 3D FFT and were compared with those obtained from the theoretical predictions using Disperse software at 150 kHz. Theoretical wavenumber plot is superimposed on the numerically obtained wavenumber plots in Figure 9a-c. Good agreements between the numerical and analytical results are obtained. It shows that the Peri-elastodynamic can predict the dispersion relation of fundamental Lamb wave modes accurately in all directions. Therefore, the strong evidences discussed above demonstrate that the Peri-elastodynamics would be a potential tool to effectively simulate the Lamb wave modes for NDE and SHM applications.

Simulation of Wave Damage Interaction
The study of wave-damage interaction is an important part of wave propagation simulation. Thus, to demonstrate the feasibility, in this article, a simulation is performed on a plate with a circular hole as shown in Figure 2. The in-plane displacements in the simulated wave field demonstrate the reflection and transmission of the respective wave modes in three different layers across the depth. Figure 10a-c, show the displacement wave fields in the plate with circular hole after 48 µS. Only the vertical displacement pattern has the reflection from the hole in the wave fields at all the three layers. However, the reflections from the hole are observed in the wave fields with in plane displacements at the top and the bottom layers only. Figure 11a,b shows the 2D-FFT analysis of the u x (x, y, t) and u z (x, y, t) displacements at the top surface of the plate, respectively. Reflected and transmitted S 0 and A 0 wave modes are clearly visible but much stronger for u x (x, y, t) displacement in Figure 11a and weaker for u z (x, y, t) displacement in Figure 11b. Boundary reflection of the same modes are also clearly visible. Based on the simulation results, a PZT sensor at a distance of 45 mm from the center of the circular hole along any direction will be suitable to detect the damage in the PZT mounted plate-like structure. If the PZT sensor is placed towards the positive X axis from the hole, the transmitted wave signals will be received by the sensor but if the PZT sensor is placed towards the negative X axis from the hole, the reflected wave signal will be received. Based on the simulation results presented in Figure 10, it can be said that if a PZT with 1-3 and 1-1 polarization is used, any location is suitable to detect the hole with both transmitted and reflected wavefronts. Because the superposition of the transmitted and reflected wavefront will contain both the in plane u x (x, y, t) and u y (x, y, t) and the out of plane deformation u z (x, y, t) of the PZT attached to the plate.
However, the reflections from the hole are observed in the wave fields with in plane displacements at the top and the bottom layers only. Figure 11a,b shows the 2D-FFT analysis of the u (x,y,t) x and u (x,y,t) z displacements at the top surface of the plate, respectively. Reflected and transmitted S0 and A0 wave modes are clearly visible but much stronger for u (x,y,t) x displacement in Figure 11a and weaker for u (x,y,t) z displacement in Figure 11b. Boundary reflection of the same modes are also clearly visible. Based on the simulation results, a PZT sensor at a distance of 45 mm from the center of the circular hole along any direction will be suitable to detect the damage in the PZT mounted plate-like structure. If the PZT sensor is placed towards the positive X axis from the hole, the transmitted wave signals will be received by the sensor but if the PZT sensor is placed towards the negative X axis from the hole, the reflected wave signal will be received. Based on the simulation results presented in Figure 10, it can be said that if a PZT with 1-3 and 1-1 polarization is used, any location is suitable to detect the hole with both transmitted and reflected wavefronts. Because the superposition of the transmitted and reflected wavefront will contain both the in plane u (x,y,t) x and u (x,y,t) y and the out of plane deformation u (x,y,t) z of the PZT attached to the plate.  Video S2 in the supplementary material shows the guided wave interaction with a hole in the plate through the u z (x, y, t) displacement of the layer L1 used in this article.

Perspectives and Peri-Elastodynamic Application to the MEMS Devices
Similar ultrasonic characterization of the plate and membrane-like structures are used in microelectro-mechanical systems (MEMS). Characterization of MEMS requires an understanding of wave propagation and wave damage interactions with the fundamental Symmetric and Antisymmetric wave modes at their respective scales. For MEMS characterization, ultrasonic waves are generated by high frequency focused transducers or by Laser ultrasonic devices. It is a daring task to understand the ultrasonic signal generated from such inspections of MEMS devices and requires a virtual analysis of wave propagation in MEMS for their enhanced or more accurate characterization. NDE is a crucial step for their ensured functionality and fail-safe design of few specific plate-like MEMS. Such virtual NDE experimental scenarios for MEMS could be simulated using Peri-elastodynamic as Figure 11. Space-time wavefield representations for top surface of the plate: (a) u x (x, y, t) for plate with a through-thickness hole (b) u z (x, y, t) for plate with a through-thickness hole.

Perspectives and Peri-Elastodynamic Application to the MEMS Devices
Similar ultrasonic characterization of the plate and membrane-like structures are used in micro-electro-mechanical systems (MEMS). Characterization of MEMS requires an understanding of wave propagation and wave damage interactions with the fundamental Symmetric and Antisymmetric wave modes at their respective scales. For MEMS characterization, ultrasonic waves are generated by high frequency focused transducers or by Laser ultrasonic devices. It is a daring task to understand the ultrasonic signal generated from such inspections of MEMS devices and requires a virtual analysis of wave propagation in MEMS for their enhanced or more accurate characterization. NDE is a crucial step for their ensured functionality and fail-safe design of few specific plate-like MEMS. Such virtual NDE experimental scenarios for MEMS could be simulated using Peri-elastodynamic as proposed in this article.

Conclusions
In this article, a numerical wave field computational tool called Peri-elastodynamics is proposed to simulate the guided waves in a plate-like structure with surface mounted PZT. Feasibility of the method is proved by simulating an SHM problem with PZT induced Lamb wave propagation in an isotropic aluminum plate. Fundamental symmetric (S 0 ) and antisymmetric (A 0 ) Lamb wave modes were generated. Further, their characteristics were investigated and compared with the theoretical predictions. Particle displacements due to S 0 and A 0 mode propagation were visualized through the vector-field plots across-the-thickness of the plate. Lamb wave modes simulated by the numerical technique were presented in the frequency-wavenumber domain and compared with those obtained from the analytical predictions. It can be concluded that if the process described in the paper is adopted meticulously, the Peri-elastodynamics can simulate the Lamb wave propagation accurately. The computational time for the SHM problem presented in this paper is approximately~48 h but can be easily accelerated by implementing the parallel computing. In this article, no parallel computing facility was used. Based on the reported computation time in the literature, it is anticipated that with the proposed formulation the wave propagation simulation using Peri-elastodynamic will be time efficient. But it is difficult to comment however, if there has been a gain in the computation time compared to the existing methods, as a direct comparison of a similar problem solved using different method is necessary and it is left for the future research. Such numerical computation is valuable for the computational NDE of structures and devices, where experimentally studying the wave damage interactions are expensive tasks. In future, virtual wave propagation in structure and MEMS devices using the complementary tools will help avoid expensive experiments but extract the right wave feature to first characterize and then certify the materials and devices.
Supplementary Materials: The following are available online at www.mdpi.com/1424-8220/18/1/274/s1, Video S1: Peri-elastodynamic simulation showing Symmetric and Antisymmetric wave in an aluminum plate. Video S2: Peri-elastodynamic simulation showing Symmetric and Antisymmetric wave interaction with a circular hole in an aluminum plate.