This paper presents the local interaction simulation approach (LISA) for efficient modeling of linear and nonlinear ultrasonic guided wave active sensing of complex structures. Three major modeling challenges are considered: material anisotropy with damping effects, nonlinear interactions between guided waves and structural damage, as well as geometric complexity of waveguides. To demonstrate LISA's prowess in addressing such challenges, carefully designed numerical case studies are presented. First, guided wave propagation and attenuation in carbon fiber composite panels are simulated. The numerical results are compared with experimental measurements obtained from scanning laser Doppler vibrometry (SLDV) to illustrate LISA's capability in modeling damped wave propagation in anisotropic medium. Second, nonlinear interactions between guided waves and structural damage are modeled by integrating contact dynamics into the LISA formulations. Comparison with commercial finite element software reveals that LISA can accurately simulate nonlinear ultrasonics but with much higher efficiency. Finally, guided wave propagation in geometrically complex waveguides is studied. The numerical example of multimodal guided wave propagation in a rail track structure with a fatigue crack is presented, demonstrating LISA's versatility to model complex waveguides and arbitrary damage profiles. This article serves as a comprehensive, systematic showcase of LISA's superb capability for efficient modeling of transient dynamic guided wave phenomena in structural health monitoring (SHM).

## Introduction

The success of guided wave-based structural health monitoring (SHM) systems substantially relies on the effective design of sensing array as well as the insightful interpretation of active sensing signals. However, identifying the optimum design parameters, such as transducer dimensions, interrogating frequency, preferable wave mode, and sensing locations, is a challenging task. In addition, the multimodal, dispersive nature of guided waves as well as their complex interaction dynamics with structural features and damage further give rise to the difficulty of properly analyzing the sensing signals. Thus, highly efficient computational models of guided wave-based SHM procedures are of critical importance for both SHM system design and signal interpretation. Several demanding requirements for the computational models can be identified: accuracy for high frequency, short wavelength, long propagating distance waves, efficiency in terms of computational time and computer resources, and versatility to explore a wide range of design parameters such as different material and geometric scenarios. However, commercially available numerical modeling tools that are based on a finite element method (FEM) cannot satisfy all these requirements, as ultrasonic wave propagation in large-scale, complex structures usually results in a computationally prohibitive problem.

In order to satisfy the desperate needs for guided wave-based SHM modeling tasks, research efforts in developing accurate, efficient, versatile modeling schemes have been widely pursued. Shen and Giurgiutiu formulated the analytical expressions with transfer functions for guided wave active sensing in simple isotropic plate structures, which was further consolidated as the analytical framework called WaveFormRevealer [1]. Rahani and Kundu [2] and Glushkov et al. [3] developed the Green's function-based distributed point source method and successfully modeled wave generation, propagation, and scattering in a highly efficient manner. Hybrid modeling techniques have also been proposed to develop efficient simulation schemes. The semi-analytical finite element method has been used and combined with local FEM to simulate wave interaction with damage in one-dimensional wave propagation problems [4–6]. Moreau and Castaings have used orthogonally analytical relation to reduce the size of FEM to obtain three-dimensional (3D) guided wave scattering features [7]. Gresil and Giurgiutiu investigated the hybrid analytical/FEM modeling concept in the time domain and achieved promising results [8,9]. Shen and Giurgiutiu further realized this hybrid concept in the frequency domain and conducted experimental verifications [10]. A recent contribution from Hafezi et al. put forward a peri-ultrasound technique for modeling ultrasonic response, which considered both linear and nonlinear wave damage interaction scenarios [11].

Nevertheless, the aforementioned modeling schemes mainly focused on relatively simple SHM configurations. Challenges would generally arise when material anisotropy, damping, structural geometric complexity, and nonlinear wave–damage interactions need to be considered. A powerful modeling technique based on the finite difference formulation and sharp interface model, known as the local interaction simulation approach (LISA), has been developed, acquiring more and more modeling capabilities to handle complex simulation tasks. Delsanto et al. first derived the one-dimensional, two-dimensional, and 3D LISA formulations for isotropic, heterogeneous materials executed on Connection Machines [12–14]. LISA underwent considerable progress during the past decade, with its application in metallic structures [15,16], extension to general anisotropic materials [17–19], coupled field capabilities [20], hybridization with other numerical methods [21,22], and execution on powerful graphics processing units (GPUs) with compute unified device architecture (CUDA) technology [23–25]. In the authors' recent study, a penalty method was deployed to introduce contact dynamics into the LISA formulation to simulate contact acoustic nonlinearity [26]. The GPU implementation has substantially facilitated low-cost supercomputing for ultrasonic elastic waves in large-scale structures.

This article presents the recent progress of the LISA framework for efficient modeling of guided wave phenomena considering material, structural, and wave–damage interaction complexities. The LISA fundamentals and its GPU implementation will be introduced. Three numerical case study examples will follow to demonstrate LISA's prowess in handling various aspects of the modeling challenges. First, guided wave propagation in anisotropic composite panels will be modeled and compared with experiential measurements, which aims to show LISA's capability in considering material anisotropy and the guided wave attenuation due to damping. Second, guided wave interaction with a fatigue crack will be modeled and compared with the results obtained from the commercial FEM package ANSYS, which would present LISA's ability to simulate nonlinear interactions between interrogating waves and the damage. Finally, guided wave propagation in a rail track with a fatigue crack will be presented, which demonstrates LISA's versatility in handling wave propagation in geometrically complex waveguides and interaction with arbitrary damage profiles.

## Local Interaction Simulation Approach Framework and Graphics Processing Units Implementation

Figure 1 presents the LISA framework, showing its derivation procedure as well as the new features. LISA approximates the partial differential elastodynamic wave equations with finite difference quotient expressions. The coefficients in LISA iterative equations depend only on the local physical material properties. A sharp interface model was used to enforce the stress and displacement continuity between the neighboring cells and nodes. Therefore, changes of material properties in the cells surrounding a computational node can be captured through these coefficients. The details of formulation derivation can be found in Ref. [18]. Guided wave generation can be achieved with an efficient frequency domain local FEM. Details of this hybrid approach can be found in Ref. [22]. Damping effects are considered based on the 3D Kelvin–Voigt viscoelasticity model. A viscosity matrix is introduced for a generic lamina with arbitrary stacking angle to capture the directional and coupling damping effects. It should be noted that the iterative equations with damping effects require the results from *previous three time steps* to determine the displacement result at current time step. A commercial preprocessor from ANSYS was integrated seamlessly into the framework for computational grid generation and material properties allocation, which enables easily modeling of structures with complex geometric features and material distribution. As mentioned above, a penalty function method was implemented in the LISA formulation to model the contact acoustic nonlinearity during wave damage interactions. A Coulomb model was used to capture the stick–slip motion at the crack interface. A nonreflective and absorbing boundary condition was enabled to minimize the model size and avoid undesired boundary reflections when required.

There are two major characteristics of the current LISA formulation that enables the computation to be expedited on GPUs. First, LISA formulations are massively parallel. This is because the computation of a general node or a contact node only depends on the solutions of its 18 neighboring nodes at the previous three time steps. Thus, the behavior of each node is independent from the others at the target time step, i.e., the computation of each node can be carried out individually in parallel. Second, the wave propagation simulation tasks usually require dense discretization of the structure, resulting in a computationally intensive problem. GPUs, with their massive concurrent threading feature, are suitable to handle such large size problems by distributing the workloads among a large number of functional units and carry out highly efficient parallel computing. In order to take advantage of the nice parallelizable feature of LISA and the superb computational capability of the powerful GPU device, the LISA procedure was implemented using CUDA. All the parameters are first established in the host memory (RAM). Then, a copy of these parameters is sent to the device memory (GPU global memory) for it to be processed. The computation of each node is assigned to a functional thread, i.e., each thread will gather the displacements of its 18 neighboring nodes at the previous three time steps, process the material properties in the eight surrounding cells, and execute the kernel to compute the displacement of this node at the current time step. Since one of the bottlenecks of a CUDA program is the data transfer between the device memory and host memory, only the required step results are gathered (every 20–30 steps depending on the frequency of the propagating waves) from the GPU to the CPU to minimize such data transfer.

## Damped Wave Propagation in Anisotropic Composite Plates

Composite materials are increasingly used in aerospace and automotive industries due to their strong, light-weight, fatigue-tolerant properties. However, such material advancement brought considerable challenges for guided-wave SHM techniques, because the wave characteristics (amplitude, velocity, and attenuation) are heavily direction dependent, which also imposes great difficulty in accurately capturing such highly anisotropic wave phenomena in computational models.

This section presents the numerical case study of damped guided wave propagation in anisotropic composite panels using LISA as well as the comparison with experimental measurements using scanning laser Doppler vibrometry (SLDV). The multilayered composite specimens used in this study were a 12-layer unidirectional [0]_{12T} carbon fiber reinforced polymer (CFRP) composite plate and a 12-layer cross-ply [0/90]_{3S} CFRP composite plate. They are manufactured from 0.125 mm thick pre-impregnated composite tape made from IM7 fibers and Cy-com 977-3 resin. Thus, the thickness of these composite panels is 1.5 mm. A 12.8-mm diameter and 0.23-mm thick circular piezowafer was bonded on the center of the plate surface for wave generation. It should be noted that the fiber volume fraction does not enter the model directly. It is considered via the parameters such as stiffness matrix, damping matrix, and density. Details of the model and its parameters can be found in Refs. [18,22]. The material mechanical properties associated with our specimen are given in Table 1, and the piezoelectric properties are given in Table 2. The method of obtaining viscosity coefficients of solids itself is a challenging branch of study. Available data are also limited for CFRP composite materials in current literature. Thus, in this study, the viscosity coefficients were obtained by updating our LISA model toward the experimental measurements. The chosen viscosity matrix coefficient for the unidirectional CFRP lamina is given in Table 3.

To demonstrate the effectiveness of the anisotropic damping formulation for an arbitrary lamination angle, the unidirectional composite panel case was first considered. Figure 2 presents the LISA simulation results of guided wave propagation and attenuation in the unidirectional carbon fiber composite plate at 75 kHz. The damping effect was realized using the Kelvin–Voigt viscoelastic formulation, details of which can be found in Ref. [22]. The first row of Fig. 2 compares the simulation results without damping effects and with damping effects. It was found that when material damping is considered, guided waves undergo considerable attenuation. Then, the unidirectional plate was rotated with various orientation angles, and the simulation results are shown in the second row. The consistence of guided wave propagation and attenuation profiles demonstrates LISA's capability in modeling wave propagation and attenuation in composite panels with an arbitrary lamination orientation.

Figure 3 shows the comparison between LISA solutions and SLDV measurements of guided wave propagation in composite panels. The first row corresponds to the results of guided wave propagation in the unidirectional [0]_{12T} CFRP composite plate at 75 kHz. The second row corresponds to the results of guided wave propagation in the cross-ply [0/90]_{3S} CFRP composite plate at 75 kHz. The first column presents the group velocity directivity curves of the fundamental wave modes S0, A0, and SH0 in these two composite plates. The second column shows the wavefield pattern from SLDV measurements, while the third column presents the LISA simulation results. It can be observed that the LISA simulation results agree well with experimental measurements, demonstrating LISA's capability and accuracy in capturing wave propagation in arbitrary lamination scenarios.

To further demonstrate the importance of considering wave attenuation due to material damping, sensing signals at the coordinate location of (0 mm, 100 mm) in Fig. 3 for the unidirectional composite panel case are presented in Fig. 4, comparing three situations: experimental measurement, LISA result with damping, and LISA result without damping. The signal amplitudes are normalized to the maximum oscillation at the transducer. It can be seen that the conventional LISA without damping effects would overestimate the A0 mode amplitude. On the other hand, the viscoelastic LISA formulation can capture the anisotropic damping behavior for guided wave attenuation. It should also be noted that the damping effects have more influence on the A0 mode than the S0 mode. This is related to the motion characteristics of their mode shapes and wavelengths.

## Nonlinear Ultrasonics Modeling Capability

In addition to the modeling challenges from material anisotropy and damping effects, the need for accurately capturing nonlinear interactions between guided waves and structural damage would further add to the complexity of the model. This is especially true when the modeling endeavors aim at simulating nonlinear ultrasonic SHM techniques [27–29].

To demonstrate LISA's capability in modeling contact acoustic nonlinearity problems, a benchmark simulation described in Fig. 5 was conducted, and the solution was compared with the result from ANSYS 14. The model consists of a 500 mm long, 20 mm wide, and 5 mm thick aluminum strip. A 10 mm long through-thickness breathing crack is located in the center of the strip. A pair of in-plane, in-phase line prescribed displacements (1 *μ*m peak to peak value) at both top and bottom of the plate surface were used to generate 100 kHz 10-cycle tone burst S0 guided waves into the structure. It should be noted that whether a certain wave mode can be generated effectively or not depends on the agreement between the excitation profile and the target wave mode shape. S0 mode mainly undergoes in-plane motion which is symmetric across the thickness. Thus, the in-plane, in-phase, top and bottom surface excitation pair can generate S0 mode. On the other hand, out-of-plane, in-phase, top and bottom surface excitation pair will generate A0 mode. In this study, the generated S0 waves will propagate along the structure, interact with the breathing crack, bring the crack information with them, and are finally picked up at the sensing point. The out-of-plane displacement at 10 mm right after the crack was recorded.

The ANSYS contact model used SOLID45 eight node structural element to discretize the geometry and CONTA173/TARGE170 contact elements to model the crack surfaces. In the LISA model, the cell size was 1 mm in the in-plane direction and 0.5 mm across the thickness. The mesh size $le$ was checked against the shortest wavelength $\lambda min$ possible at the third harmonics, i.e., $le\u2264\lambda min/20$. A time step $\Delta t$ of 0.125 *μ*s was utilized for the time marching Newmark-beta integration to ensure the accuracy up to the third harmonic frequency, i.e., $\Delta t\u22641/20fmax$. The time marching step was 0.06 *μ*s obtained from the Courant–Friedrichs–Lewy (CFL) number requirement [18]. It should be noted that a denser mesh is used in the thickness direction as an effort to accurately depict the mode shapes of guided waves. The computational grid on the mode shape dimensions imposes much influence on the guided wave modeling accuracy. A systematic study on this topic has been conducted by Nadella [30].

Figure 6 shows the verification of the contact LISA model against the ANSYS solution. It can be noticed that the signals agree with each other in time domain with merely slight differences. The frequency spectra shown in Fig. 6(b) also demonstrated that the frequency components of fundamental excitation, low-frequency DC component, second, third, even fourth higher harmonics compare very well with each other. Differences only appeared at very high-frequency range, where the time step is not sufficiently small enough to generate accurate results for meaningful comparison between the two methods. This numerical verification against commercial finite element code attested the capability and validity of the new contact LISA model.

It should be noted that the new contact LISA model achieved much higher computational efficiency over the conventional nonlinear FEM simulation. Both computational tasks were conducted on an Asus ESC2000 G2 workstation with a 2.00 GHz Intel Xeon E2-2650 processor, 32 GB of 1.60 GHz memory, and an Nvidia GeForce GTX Titan graphics processor with 2688 CUDA cores. The FEM simulation with 279,900 degrees-of-freedom (DOFs) took around 19 s for each time step, resulting in a total computational time of 8 h for 1500 time steps. On the other hand, LISA simulation with 648,120 degrees-of-freedom merely consumes around 0.043 s for each time step, resulting in a total computational time of 2.15 min for 3000 time steps. Thus, it is apparent that the new contact LISA model is much more efficient than the conventional nonlinear FEM simulation, while achieving comparable results.

The efforts in capturing the nonlinear ultrasonic phenomena are important, because the nonlinear ultrasonic signals may contain rich diagnostic information to qualify and quantify the damage. This aspect is further illustrated with Fig. 7, which presents the simulation results and compares the sensing signals among the pristine, notch, and the fatigue crack cases. The notch is modeled by replacing a thin layer of material elements with air cells. The fatigue crack is introduced using a discrete meshing strategy which creates contact pairs at the crack interface. Detailed steps to generate the fatigue crack can be found in Ref. [26]. It can be noticed that the time trace of the notch case resembles the waveform of the pristine signal with an amplitude drop and balanced positive and negative oscillatory amplitudes. In contrast, the time trace of the fatigue crack case is heavily distorted. The contact–impact between the crack surfaces induced high-amplitude spikes. The unbalanced waveform results in a low-frequency modulation on the entire wave packet. The waveform also shows zigzags unlike the smooth pattern present in the notch crack case. The frequency spectra of these two signals show great difference too. It can be observed that the pristine and notch signal spectra show only the excitation frequency component at 100 kHz. No higher harmonics nor DC component exists. On the other hand, the fatigue crack signal presents distinctive nonlinear higher harmonics and DC component in its spectrum. The generation of higher harmonics and DC component is a classical nonlinear phenomenon, and it serves as the basis for many nonlinear ultrasonic inspection methodologies. This numerical example shows that LISA formulations with contact dynamics are capable of simulating such nonlinear ultrasonic signals.

## Wave Propagation in Geometrically Complex Waveguides

The most commonly encountered modeling challenge comes from the large, complex structural geometries. It may not cause a concern in commercial FEM packages, but it remains a challenge for many of the newly developed computationally efficient schemes. LISA's seamless integration with commercial preprocessors enables the handling of geometrically complex waveguides. To substantiate this aspect, this section presents the numerical case study of utilizing LISA to model guided wave propagation in a rail track.

### LISA Model for a Rail Track.

Figure 8 presents the LISA numerical model for the investigation of guided wave propagation and interaction with a fatigue crack in a rail track. Two lines of surface traction forces are utilized to simulate a pair of transducers generating guided waves into the rail structure. A 100 kHz 10-cycle tone burst is used for wave generation. The ultrasonic waves propagate along the rail track, interact with a fatigue crack, and are finally picked up at the sensing location. Absorbing layers with increasing damping is used on both ends of the model to eliminate boundary reflections, enabling the simulation of wave propagation in an infinite long rail track with a finite dimensional model. A 1 mm mesh size is adopted for the cross-sectional plane, while a 2 mm mesh sized is deployed for the track direction. The time step according to the CFL condition is 110.33 ns, which corresponds to a CFL number of 0.99. The 3D LISA grid shows that the adopted discretization can capture the complex geometric details of the rail track, i.e., LISA can handle very complex structural geometries.

In order to depict the fatigue damage, a parametric fatigue crack zone is used in the LISA model. Figure 9(a) shows the profile of a typical fatigue zone in a rail track. Figure 9(b) presents the parametric fatigue zone characterized by the size of *a*. In this study, fatigue cracks of progressive sizes are investigated, i.e., *a* = 0 mm, 5 mm, 10 mm, 15 mm, and 20 mm, corresponding to pristine case, incipient damage case, median damage case, and sever damage case.

### Guided Wave Propagation and Interaction With Damage in the Rail Track.

Figure 10 presents the LISA simulation result of guided wave generation, propagation, and interaction with a fatigue crack in the rail track. Figure 10(a) shows the snapshots of guided waves generation by the line traction forces at 50 *μ*s, guided wave propagation, mode development, and interaction with the crack at 150 *μ*s, as well as guided wave absorption by the absorbing layers with increasing damping boundaries at 300 *μ*s. It can be seen that the guided wave modes are very complex in the rail track. Figure 10(b) shows the crack open and close phenomenon during the wave damage interaction, which introduces the contact acoustic nonlinearity into the sensing signal.

Figure 11(a) presents the comparison of time domain sensing signals between the pristine case (*a* = 0 mm) and the severely damaged case (*a* = 20 mm). The signals are very complex with multiple wave packets, arisen from the multimodal and dispersive nature of guided waves. According to the conventional linear ultrasonic techniques based on phenomena such as amplitude change, phase shift, or scattering, it is very difficult to make an appropriate assessment of the damage existence and severity. On the other hand, from the perspective of nonlinear ultrasonic features, the damage signature is much more evident. Figures 11(b) and 11(c) show the time–frequency domain sensing signals of the pristine case and the severely damaged case, respectively. For the pristine case, only the excitation frequency *f*_{0} exists. However, for the damaged case, in addition to *f*_{0}, superharmonics are present at 2*f*_{0}, 3*f*_{0}, 4*f*_{0}, and so on. Such a unique nonlinear feature allows the early detection of incipient fatigue damage in structures. Superharmonics is a classical nonlinear ultrasonic phenomenon which has been widely reported using experiments in the literatures [27–29]. Again, the focus of this research is on developing an efficient numerical approach to model such a nonlinear effect. Experimental verifications on this rail model should be conducted in a future work.

This numerical case study demonstrated LISA's capability in handling guided wave propagation in geometrically complex waveguides and nonlinear interactions with structural damage with arbitrary profiles. It should be emphasized again that for this case study, the model size reached 10,192,107 degrees-of-freedom, which is computationally prohibitive for most conventional FEM packages. It is LISA's parallel GPU implementation that facilitates its highly efficient modeling capability.

### LISA Model Data for Damage Quantification Analysis.

To further quantify the damage severity, a damage index (DI) can be developed based on the nonlinear energy proportion in the LISA model data. Figure 12(a) shows the 3D view of the time–frequency sensing signal of the severely damaged case. The magnitude represents the energy participation of each frequency component along the time history. In order to consider all the wave modes participating in the wave propagation, the time domain integration is used at the fundamental and superharmonic frequencies to construct the DI, which is given by

where $P(t,nf0)$ denotes the magnitude at the *n*th superharmonic location at time $t$. It should be noted that the magnitude should be elevated above zero by subtracting the minimum base value in the spectral graph. The DI represents the nonlinear energy ratio, designating the degree of nonlinearity in the sensing signal. The first term consists of the integration over the “peaks” of the spectrum, while the second term is the integration over the “valleys” of the spectrum serving as the inherent baseline. By subtracting the inherent baseline, this DI becomes a indicative measurement of signal nonlinearity. Figure 12(b) presents the damage index results for various crack sizes. It can be noticed that the DI is of high sensitivity; even at incipient damage level (*a* = 5 mm) the DI increased noticeably. Thus, such high sensitivity will allow the detection of small cracks at the early stage. In addition, the DI also grows monotonically with the damage severity, which allows crack growth monitoring and the quantification of the damage severity.

## Concluding Remarks

This article presented the local interaction simulation approach for efficient modeling of linear and nonlinear ultrasonic guided waves for active sensing of complex structures. Three major modeling challenges were addressed: material anisotropy with damping effects, nonlinear interactions between guided waves and structural damage, as well as the geometric complexity of waveguides. First, guided wave propagation in unidirectional and cross-ply carbon fiber composite panels was modeled and compared with experimental measurements. The LISA solutions agreed well with the experimental results obtained by scanning laser Doppler vibrometry, which demonstrated LISA's capability in modeling the directivity of guided wave propagation and attenuation in composite materials. Second, the nonlinear interactions between guided waves and structural damage were considered. The verification against the commercial FEM package ANSYS not only showed LISA's accuracy in capturing the contact acoustic nonlinearity but also highlighted LISA's superb computational efficiency over the conventional FEM approach. Third, the case study of guided wave propagation and interaction with a fatigue crack in a rail track was conducted. The simulation results showed that LISA was capable of handling very complex structure and damage geometries. The further analysis of LISA generated sensing signals also allowed the development of a nonlinear damage index for the quantification of the fatigue damage. This article systematically demonstrated LISA's prowess for efficient modeling of transient dynamic guided wave phenomena in SHM procedures.

## Acknowledgment

This work was also partially sponsored by the National Rotorcraft Technology Center (NRTC) Vertical Lift Rotorcraft Center of Excellence (VLRCOE) at the University of Michigan, with Mahendra J. Bhagwat as the program manager. Opinions, interpretations, conclusions, and recommendations are those of the authors and are not necessarily endorsed by the U.S. Government.

## Funding Data

The National Natural Science Foundation of China (Contract No. 51605284).