Explainable, Flexible, FRF-Based Parametric Surrogate for Guided Wave-Based Evaluation in Multiple Defect Scenarios

1. Introduction

Within the context Non Destructive Evaluation (NDE) and SHM applications of thin plate and shell structures, Lamb waves [1] offer several favorable properties. Firstly, their short wavelength enables interaction with small defects and, therefore, facilitates defect detection in early damage stages [2]. Further, compared to bulk waves, they exhibit reduced geometrical damping allowing to cover large areas with only few sensors. Additionally, the availability of commonly used, low cost PZT-sensors, renders Lamb wave-based SHM schemes accessible for a wide range of applications [3]. Their flexibility allows for monitoring of various defects, such as fatigue cracks [4,5], notches [6,7], delaminations [8,9] or corrosion [10].

Due to their multi-modality and their dispersive nature, Lamb waves form complex patterns. To simplify these patterns, many Lamb wave-based SHM-schemes rely on baseline subtraction to obtain only the response influenced by scattering from defects [11]. In order to keep the scattering response simpler [12], the excitation of single modes is often attempted, which can be achieved either through specially designed actuator setups [13] or through specific frequency tuning [14]. However, the localization of defects still calls for advanced analysis methods, which can be largely classified into the data driven and model-based classes.

Popular data driven methods include Delay-And-Sum (DAS) [15,16,17] approaches and the Sparse Reconstruction (SR) [18,19], while further data driven methods rely on the Time-Reversal Method (TRM) [20] and the Time of Flight (ToF) [21,22]. More recently, Deep Learning (DL) approaches have gained ground with respect to more essential data-driven schemes. A comparison of deep learning models developed for the specific task of GW-based defect detection, including Dense Neural Networks (DNNs), Convolutional Neural Networks (CNNs), Recurrent Neural Networks (RNNs) and Long Short-Term Memory (LSTM) networks is given by Rautela and Gopalakrishnan [23]. In a relevant application which considers a single delamination Azadi et al. [24] use a total of ten Carbon Fiber Reinforced Polymer (CFRP) laminates – nine damaged and one healthy to train their 1D-CNN. A laser sequentially creates a localized thermal expansion, allowing for broad band excitation in a grid of points. This excitation, in combination with piezoelectric sensors, positioned in three locations, allows to record an experimental dataset with more than 2 million measurements. Using this dataset and hyperparameter optimization, the 1D-CNN proves to offer accurate delamination classification from raw signal data. The aspect of training is particularly important in establishing robust DL-based schemes. In this vein, Rautela et al., [25] introduce a two-step detection and localization approach. They specifically train their 1D-CNNs on data collected from CFRP plates with attached aluminum disks, which are employed to simulate damage. From this dataset, physical features, such as the time of flight are extracted. In further work, in an effort to reduce the training data, Postorino et al.[26] apply transfer learning with optimal transport on data generated with a numerical model, which was validated with experiments.

In the present work, data driven methods exploiting the principle of FRFs or transfer functions are of special interest, as the FRF forms a central concept in the development of the proposed surrogate. Sun et al. [27] use the concept of electrical transfer functions to verify the rationality of their measurements. Electrical transfer functions are further used by Wang and Shen [28], in the context of the time-reversal method. The transfer function between the excitation and the measured signal is modified in order to compensate for the influence of the transducer. A drawback of their method lies in the necessity for mode separation, which is often not feasible, especially if the distance between the actuator and sensor is small. The surrogate proposed in the present paper shows the ability to track individual modes and, therefore, it allows to detect specific modes, from which reflections at different parts of the model stem from. In the Decomposition of the Time-Reversal Operator and MUltiple SIgnal Classification (DORT-MUSIC) technique, the transfer function is evaluated for a set of sensor-transducer pairs and, based on individual transfer functions, a transfer matrix is constructed. This matrix is decomposed by means of a Singular Value Decomposition (SVD), which allows to extract reflectivity characteristic features from signals originating from different sources, e.g. scatterers. To obtain the transfer matrix, He and Yuan [29] use a sensor setup, where the actuators are perpendicular to the sensors, while Fan et al. [30] rely on a single piezoelectric array. To improve the localization of debonding in CFRP-reinforced steel structures, An et al.[31] derive a Waveform Feature Index (WFI) based on changes in the Transfer function. A method to experimentally derive FRF is given by Yang et al. [32], who use the step function as an excitation. They show that this approach is efficient for FRFs up to 900 $\mathrm{k}$$\mathrm{Hz}$.A further application of the FRF, is in the combination with the wave focusing approach. In this approach, the excitation is modified in order to account for the dispersive nature of Lamb waves [33]. Based on use of the FRF, Yang et al. [34] reduced the number of required physical recordings needed for wave focusing. Further, the utilization of the FRF allows to seamlessly predict the response for varying excitations. This flexibility of the FRF in dealing with various excitations, is a key property, which is exploited in the creation of the surrogate proposed in the present paper.

Despite their attractive trait in using minimal prior knowledge, data-driven methods often remain limited to the setups they were trained on. To increase the flexibility of SHM-schemes, mechanical knowledge in the form of the governing equations of the wave propagation problem can be introduced [35]. Such model-based approaches come with the bottleneck of a high computational time, particularly for Lamb-wave problems, as the high frequencies and short wave lengths involved require fine discretizations in space and time for mesh-based solvers [3]. As a consequence of the high associated computational cost, inverse algorithms need to be efficient with respect to the number of required evaluations of the model. Such an approach is introduced by Bürchner et al. [36], who perform Full Waveform Inversion (FWI) for the 1D wave equation. They introduce a scaling factor for the density of the structure, with the scaling factor tending to 0 for void regions and 1 for solid regions. This allows to calculate gradients efficiently during the optimization process. In turn, exploiting these gradients in a gradient-based optimization scheme, leads to a fast convergence of the inverse problem. This approach was extended to solve elastodynamic problems in geophysics applications [37], where scaling is applied on the shear modulus, rather than the density. The robustness of the approach is further demonstrated for a numerical example with artificial noise. Based on the use of density scaling, Bürchner et al. [38] introduced a model-based approach, where the inverse problem is solved in two steps. In the first step, a coarser resolution is used to identify the area of the defects (in the form of voids) and, in a second step, in the area around the defect, the material properties are modeled in a refined manner, to obtain a sharp resolution for the representation of the void. Similarly, the Finite Cell Method (FCM) allows to model features in the plate without the use of extensive mesh refinement. This concept is combined with Particle Swarm Optimization (PSO) by Zakian et al. [39] for damage detection and localization. An approach, where Reduced Order Models (ROMs) are exploited in the inverse setting is given by Rautela and Gopalakrishnan[40], who use a ROM for the Spectral Finite Element Method (SFEM) in frequency domain in combination with CNNs for damage detection and RNNs for damage localization. This approach allows for damage detection, localization and severity estimation.

Despite the efficiency of the demonstrated model-based approaches, the computational toll remains high, as long as a high-fidelity model needs to be solved, which prohibits their application in real time or near real time monitoring. ROMs in contrast allow for expedited evaluations of the forward model. One class of ROMs is projection-based Model Order Reduction (MOR), where the matrices of the Full Order Model (FOM) are projected into a lower dimensional space [41]. In this space, the governing equations are solved efficiently, while the solution of the system is subsequently transformed into the original space to recover the full solution. For wave propagation problems, however, the Kolmogorov n-width, which is a measure of the accuracy with which low-dimensional spaces can represent the solution of a high-dimensional problem, decays slowly. Thus, the size of the reduced bases needed to achieve the acceptable levels of accuracy might not lead to substantial computational gains [42]. An efficient use of projection based ROMs for wave propagation, can be realized through local projection bases [43]. Thereby, the initially collected snapshots, based on which projection bases are constructed, are clustered, with a separate basis being constructed for each cluster. This leads to bases of smaller size, however it also requires the selection of an optimal basis at each time-step of the online solution. While projection-based ROMs aim to reduce the number of Degrees of Freedom (DOFs), an alternative approach, proposed by Bigoni [44] aims at reducing the number of linear system solutions required, which are also the most computationally demanding part of such a simulation. To this end, the ROM is evaluated in the frequency domain for the principal frequencies, which are significantly less than the number of time steps that would be required to solve the problem in time domain. Such a frequency domain based ROM was used by Bigoni et al. [45] for damage localization. Building on a similar concept, the surrogate demonstrated in the present paper is evaluated in the frequency domain - only for a discrete subset of frequencies of interest. This concept was further applied to develop a ROM based on FRFs for detection of single defects in the authoring team’s previous work [46]. Through geometrical transformations, defects of equivalent FRF values were found and then the FRFs are interpolated exploiting the Matching Pursuit algorithm. An alternate approach lies in the idea of decoupling the time and space domains of the underlying differential equations. This approach, named Proper Generalized Decomposition (PGD), is known to avoid the the curse of dimensionality [47]. Using this scheme, a global-local model for beams is presented by Yan et al. [48], who model the main parts of the beams using the wave finite-element method, while the section around the defect is modeled via conventional Finite Element Method (FEM). Another surrogate, which relies on the combination of Bayesian optimization with a machine learning approach, enables uncertainty quantification with respect to the material parameters in composite plates [49]. Further, in the framework of SHM, Postorino et al. [50] create a metamodel with radial basis functions. As this metamodel can rapidly evaluate the response for different defect locations, it can be combined with PSO to localize damage in a single defect scenario.

Surrogate models based on ray tracing track individual wave packets through the plate, allowing to gather further information about the wave packets at the sensor location. The sensor response is then found by summing the single wave packets together. Further, the computational time of such surrogates is independent of the size of the plate. Ray tracing can be performed using the pencil method, where the wave is tracked with an axial ray and a paraxial ray is introduced to track its amplitude [51]. The method is able to model reflections and mode conversion precisely, even at curved edges and allows to model wave sources of finite size. This method was extended to model laminate plates [52]. Mardanshani et al. [53] combine the pencil method with the wave finite element method to model wave propagation in plates with joints. In this concept, the wave finite element method models the transition at the joints, while the pencil method tracks the energy of the wave. Based on an analytical scattering model, Xu et al. [54] provide a ROM that uses a dictionary approach for the scattering, which is combined with the modeling of edge reflections through mirroring. Localizing a single defect with this model, outperforms the DAS method regarding the necessary receivers to accurately localize the defect. Similarly to this approach, we explore a surrogate with ray tracing, relying on a numerical model to create a dictionary for the scattering behavior. Allouko et al. [55] model the scatterer with finite elements and the area of the pristine plate through semi-analytical solutions for half-spaces. The decoupling of the space around a defect (near field) and the far field is also the base of the model introduced by Lozano et al. [56]. They model both fields with different versions of the Scaled Boundary Finite Element Method (SBFEM) and couple the models with the mortar method. Another analytical scattering model is introduced by Briand et al. [57], who explore a damage index based on the amplitude of the scattering model for an estimation of the damage size. A fast analytical predictor is also given by Poddar and Giurgiutiu [58], who split the area around notches and cracks in several different domains and adjust the boundaries of the domains according to the neighboring domains. This model is called Complex Modes Expansion with vector Projection (CMEP).

Generally, as seen in the aforementioned ray-tracing models, the scattering behavior of defects and flaws needs to be investigated carefully, as it significantly influences the sensor response. Such investigations have been performed numerically [6,59,60,61,62,63], analytically [60,63] and experimentally [61,62,63] for different types of defects. The results of these analyzes indicate that the reflection is dependent on the mode and frequency of the incident wave [62]. Sedaghati et al. [60] derive an analytical formulation for the scattering of the $S_0$ wave mode by coupling the Poisson plate theory with the Mindlin plate theory. They analyze symmetric and asymmetric double-sided blind holes and conclude that the thickness and the type of the hole are relevant to the scattering behavior. This, in turn, reveals that conclusions on the size and type of defects can be drawn on the basis of characteristics of the scattering signals. In their paper, they provide a summary of further analytical solutions of the wave scattering phenomena for different types of defects. The investigation of the $A_0$-mode shows that the Scattering Directivity Pattern (SDP) is asymmetric for angled cracks [6] and that for some angles of incident waves there are no reflections [61]. The latter seems to be also apparent for the circular defects investigated in the present paper. Fromme and Rouge [61] conclude that angles with no reflections need to be considered, when designing and mounting an SHM-system. Further, Soleimanpour and Soleimani [59] investigate the clapping phenomenon for cracks in an aluminum plate for the $S_0$-mode. The clapping leads to a nonlinearity in the scattering behavior, which creates higher harmonics. The paper focuses on the first two harmonics. It shows that the observed nonlinear effects are well suited for damage detection. If a setup where linear waves are investigated is used they recommend to use a pitch-catch configuration, where the actuator and the sensor are in different sides of the defect. In the present contribution, we will focus on linear waves that are indeed examined in such a pitch-catch setup. In their study on $S_0$ waves with cylindrical holes, Diligent et al. [63] investigate relations between the size of the defect and the scattering coefficient, as well as between the distance from the defect and the scattering coefficient. The latter relation is described through a simple formula, reflecting geometrical damping, which is used herein to find the scattering waves for sensors at any distance to a scatterer.

In the current paper, a novel approach is proposed, which combines several of the aforementioned aspects, to achieve flexible surrogate models, able to predict the sensor response in plates with several flaws. An online-offline approach is adopted, where the scattering behavior for different flaws is firstly learned in the form of transmissibility functions and next stored in a dictionary during an offline phase. The created dictionary of transmissibility functions is integrated in a ray tracing approach, allowing to consequently evaluate the sensor response.

The remainder of the paper is structured as follows. In Section 2, the theory underling the surrogate is presented, followed by the description of the proposed approach in Section 3, pertaining to both the online and offline phase of the model, including algorithmic implementation details. In Section 4, the method is demonstrated on a numerical example and results for two different kinds of defects are investigated. The paper is concluded by a discussion on the advantages and limitations of the model and an outlook on further research.

2. Theory

2.1. Lamb Waves

Lamb waves are multi-modal waves propagating in thin plates, which means that, generally, more than one mode is propagating (Figure 1), with an infinite number of these modes existing in theory. Generally, these modes are classified, based on the resulting stress distribution along the cross section, into symmetric and anti-symmetric modes [64]. However, most modes comprise a cut-off frequency below which they do not propagate. The cut-off frequencies depend on the material and thickness of the plate. Two Lamb wave modes, the fundamental symmetric, $S_0$, and the fundamental anti-symmetric, $A_0$, mode exist, which propagate across all frequencies [65]. For most SHM-applications the excitation is chosen so that they operate at the frequency range where only these two fundamental modes propagate [8,21,66]. In the presence of perturbations within the plate, such as defects [67] or in the presence of boundaries [68], mode conversion can occur. Further, Lamb waves are characterized by their dispersive nature, which implies that their propagation speed is frequency dependent (Figure 1). Since wave packets comprise of multiple frequencies, rather than a single one, different propagation speeds are present within the wave packet. As a consequence, a wave packet will change its form over time, while its amplitude will drop.

2.2. Frequency Domain Representation

In this section, an arbitrary time signal is considered, as illustrated in Figure 2. The Fourier transform $\mathcal{F}$ of a real function shifted in time $f(t-\Delta t)$ results in a complex, periodic function:

$$\mathcal{F}\left(f(t-\Delta t)\right)=F\left(\omega \right)·e^{-i\omega \Delta t}$$

(1)

where $F\left(\omega \right)=\mathcal{F}\left(f\right(t\left)\right)$, e denotes Euler’s number, $\omega$, and $i^2=-1$. This relationship can be used to describe a wave measured by a sensor placed at a varying distance from the source [70] and will be used in what follows to describe the scattering behavior of waves. As a comment, it should be mentioned that the Fourier transformation requires stationary signals.

2.3. Frequency Response Function (FRF) and Transmissibility Function

FRFs and transmissibility functions, which are well established in SHM [71], allow to evaluate a system’s response for various time history inputs. For the case of Lamb-wave propagation investigated in this paper, both functions can be derived starting from the discretized form of the equilibrium equations of a Linear Time Invariant (LTI) Multi Degree of Freedom (MDOF) system:

$$\mathit{M}\ddot{\mathit{u}}\left(t\right)+\mathit{C}\dot{\mathit{u}}\left(t\right)+\mathit{K}\mathit{u}\left(t\right)={\mathit{f}}_{\mathit{s}}·f\left(t\right)$$

(2)

where $\mathit{M}$,$\mathit{C}$, and $\mathit{K}$$\in {\mathbb{R}}^{n\times n}$ describe the discretized mass, damping, and stiffness matrices, respectively, vectors $\mathit{u}$, $\dot{\mathit{u}}$, and $\ddot{\mathit{u}}$$\in {\mathbb{R}}^{n\times 1}$ denote the deformation vector and its first and second derivative with respect to time t. The right hand side of Equation 2 describes the excitation of the system, where ${\mathit{f}}_{\mathit{s}}$ denotes the spatial distribution of the excitation and $f\left(t\right)$ its time history. This split of the excitation is possible, when the spatial distribution of the excitation is time invariant, which is usually the case for Lamb wave excitation [72]. Applying the Fourier transform to Equation 2 one obtains:

$$(\mathit{K}+i\omega \mathit{C}-{\omega }^2\mathit{M})\mathit{u}\left(\omega \right)={\mathit{f}}_{\mathit{s}}·f\left(\omega \right)$$

(3)

In the above, the same symbols are used for displacement $\mathit{u}$ and excitation f, regardless of whether they refer to frequency or time domain, since the corresponding domain is evident from their dependence on time or frequency respectively.

Reorganizing Equation 3 yields the FRF:

$$FRF\left(\omega \right)={\displaystyle \frac{\mathit{u}\left(\omega \right)}{f\left(\omega \right)}}={(\mathit{K}+i\omega \mathit{C}-{\omega }^2\mathit{M})}^{-1}{\mathit{f}}_{\mathit{s}}$$

(4)

For discrete frequencies ${\omega }_i$ the FRF can be evaluated numerically:

$$FRF\left({\omega }_i\right)={\displaystyle \frac{\mathit{u}\left({\omega }_i\right)}{f\left({\omega }_i\right)}}$$

(5)

similarly, if the FRF is available, the response for any given excitation can be computed as:

$$\mathit{u}\left({\omega }_i\right)=FRF\left({\omega }_i\right)f\left({\omega }_i\right)$$

(6)

Next, the transmissibility function is derived for discrete frequencies from Equation 3. First, the equation is simplified by substituting $\widehat{\mathit{K}}\left(\omega \right)=\mathit{K}+i\omega \mathit{C}-{\omega }^2\mathit{M}$:

$$\widehat{\mathit{K}}\left(\omega \right)\mathit{u}\left(\omega \right)={\mathit{f}}_{\mathit{s}}·f\left(\omega \right)$$

(7)

Next, it is assumed that the response is known at a single DOF and should be calculated for all remaining locations. To this end, the system is partitioned in two sets of equations:

$$\left[\begin{array}{cc}{\widehat{\mathit{K}}}_{uu}\left(\omega \right)& {\widehat{\mathit{K}}}_{uk}\left(\omega \right)\\ {\widehat{\mathit{K}}}_{ku}\left(\omega \right)& {\widehat{K}}_{kk}\left(\omega \right)\end{array}\right]\left[\begin{array}{c}{\mathit{u}}_u\left(\omega \right)\\ u_k\left(\omega \right)\end{array}\right]=\left[\begin{array}{c}{\mathit{f}}_{s,u}\\ f_{s,k}\end{array}\right]·f\left(\omega \right)$$

(8)

where subindex k refers to the DOF, where the response is known, and the subindex u to the unknown DOFs. This means $u_k\left(\omega \right)\in {\mathbb{R}}^{1\times 1}$ and ${\mathit{u}}_u\left(\omega \right)\in {\mathbb{R}}^{m\times 1}$, where $m=n-1$. Solving the second row of Equation 8 for $f\left(\omega \right)$ and substituting $f\left(\omega \right)$ in the first row yields:

$${\widehat{\mathit{K}}}_{uu}\left(\omega \right){\mathit{u}}_u\left({\omega }_i\right)+{\widehat{\mathit{K}}}_{uk}\left(\omega \right)u_k\left(\omega \right)={\displaystyle \frac{{\mathit{f}}_{s,u}}{f_{s,k}}}{\widehat{\mathit{K}}}_{ku}\left(\omega \right){\mathit{u}}_u\left(\omega \right)+{\widehat{K}}_{kk}\left(\omega \right)u_k\left(\omega \right)$$

(9)

Reorganizing Equation 9 and solving for the ratio of ${\displaystyle \frac{{\mathit{u}}_u\left(\omega \right)}{u_k\left(\omega \right)}}$ results in the definition of transmissibility function:

$$\mathit{T}\left(\omega \right)={\displaystyle \frac{{\mathit{u}}_u\left(\omega \right)}{u_k\left(\omega \right)}}={\left({\widehat{\mathit{K}}}_{uu}\left(\omega \right)-{\displaystyle \frac{{\mathit{f}}_{s,u}}{f_{s,k}}}{\widehat{\mathit{K}}}_{ku}\left(\omega \right)\right)}^{-1}\left({\displaystyle \frac{{\mathit{f}}_{s,u}}{f_{s,k}}}{\widehat{K}}_{kk}\left(\omega \right)-{\widehat{\mathit{K}}}_{uk}\left(\omega \right)\right)$$

(10)

For discrete frequencies ${\omega }_i$ the transmissibility function $\mathit{T}\left({\omega }_i\right)$ can be evaluated numerically:

$$\mathit{T}\left({\omega }_i\right)={\displaystyle \frac{{\mathit{u}}_u\left({\omega }_i\right)}{u_k\left({\omega }_i\right)}}$$

(11)

while, if the transmissibility function is available, the response for any given response at DOF k can be obtained as:

$${\mathit{u}}_u\left({\omega }_i\right)=\mathit{T}\left({\omega }_i\right)u_k\left({\omega }_i\right)$$

(12)

3. Method

The solver presented in this paper is based on the principle of ray tracing. Following this approach, each ray is represented by either discrete FRFs or discrete transmissibility functions – one function for each discrete frequency of interest. For each feature (e.g. defects or rivet holes) within the plate, several parameters need to be evaluated in an offline phase, which is described in subsection 3.1, followed by a section on the use of these parameters in the online stage (subsection 3.2). At this point, the following assumptions are put in place:

• The material is linear-elastic and isotropic.
• The features are assumed to be of circular shape, and therefore the response is not affected by the direction of the incident wave.
• The size of the features is small compared to the size of the distances between actuator, sensors and features.

As the assumption of circular features introduces a strong limitation, it should be stated here that it is assumed that this assumption can be relaxed in future work. At this point, this assumption allows to train the ROM with very few simulations. A higher number of simulations is expected to allow the ROM to be trained for small features of any form.

3.1. Offline Phase

As a first step, the scattering behavior of an arbitrary feature is numerically investigated to extract the relevant information needed for the online phase. To extract this information, an idealized setup of an infinite plate with a single actuator, one feature and one sensor (Figure 3) is chosen, where only a single mode is excited. In this setup, a direct wave packet from the actuator and one or several wave packets stemming from the feature can be measured at the sensor location. Several wave packets with their origin at the feature might be present, as mode conversion might occur at the feature – depending on the frequency and the feature’s shape.

In what follows we investigate the different wave packets, which reach a sensor, separately. Firstly, a baseline simulation is performed using the pristine setup, i.e., one without defects. Through baseline subtraction, the waves stemming from the feature can be extracted, given that the direct wave is known from the baseline. To differentiate between the different scattering waves from the feature, the deformations or internal forces on the cross section need to be investigated. Although this task is challenging in an experimental setting, it can be easily accomplished in simulation. Using this approach, all possible cases of wave packets propagating from the actuator to the feature and then to the sensor should be investigated. This means that in a setup where the two fundamental modes ($S_0$ and $A_0$) are present and mode conversion occurs at the feature, four cases must be investigated. These cases are described by mode pairs, where the first mode is the one excited at the actuator and the second mode is the one investigated at the sensor: $S_0-S_0$, $S_0-A_0$, $A_0-S_0$, and $A_0-A_0$. Due to the linearity assumption, these cases can simply be superposed in a later step to obtain the full sensor response.

Next, several sensors are used in a radial line from the defect, as shown in Figure 4. The first sensor needs to have a certain distance from the feature to prevent evanescent waves from affecting the sensor response. A single mode wave packet, originating from the feature (due to either reflection of the same mode or mode conversion), changes its shape in between two sensors due to two effects, namely dispersion and geometric damping. Further, the wave packets arrive later in time as the distance between sensor and feature grows. These three influences (dispersion, geometrical damping, time delay) on the sensor response are investigated next for discrete signals in frequency domain. Firstly, the geometric damping is known to be inversely proportional to $\sqrt{d}$[63,73], where d is the distance to the origin of the wave, and independent of the frequency. In contrast, dispersion is frequency dependent and does not occur if the wave is composed of a single frequency. In the described method, the signal is separately investigated for each discrete frequency, implying that the dispersion is inherently accounted for, as it is a purely frequency dependent phenomenon. As for the differentiation in terms of the time of arrival, it is shown in subsection 2.2 that the response will be periodic. The time shift $\Delta t$ in Equation 1 thereby reflects the difference in distance $\Delta d$. The relation between $\Delta t$ and $\Delta d$ is expressed as:

$$\Delta t={\displaystyle \frac{\Delta d}{c\left(\omega \right)}}$$

(13)

where $c\left(\omega \right)$ is the frequency dependent wave speed for the investigated mode and frequency ($\omega$).

This means that the response at a sensor location for a discrete frequency and a certain mode can be described as:

$$u\left(\omega \right)={\displaystyle \frac{A}{\sqrt{d}}}·e^{i\left({\displaystyle \frac{d}{c\left(\omega \right)}}\omega +\theta \right)}$$

(14)

with the unknown parameters being the amplitude A and phase shift $\theta$. The frequency of the periodic term with respect to the distance ${\omega }_{dist}={\displaystyle \frac{\omega }{c\left(\omega \right)}}$ can be evaluated, using the dispersion characteristics, as displayed in Figure 1. To avoid confusion, it should be emphasized at this point that ${\omega }_{dist}$ is not a frequency with respect to time, but a spatial frequency with respect to the distance from the origin of the wave. In what follows, it is assumed that the sensor response along a radial line from a source can be adequately described by Equation 14.

Next, the amplitude A and phase shift $\theta$ need to be determined. With the amplitude A and the distance to the source d being strictly positive, the magnitude of the sensor response in Equation 14 becomes:

$$|u\left(\omega \right)|={\displaystyle \frac{A}{\sqrt{d}}}·|e^{i\left({\displaystyle \frac{d}{c\left(\omega \right)}}\omega +\theta \right)}|={\displaystyle \frac{A}{\sqrt{d}}}$$

(15)

And, therefore, as $|e^{i\left({\displaystyle \frac{d}{c\left(\omega \right)}}\omega +\theta \right)}|=1$:

$$A=|u\left(\omega \right)|·\sqrt{d}$$

(16)

The former implies that the amplitude could be extracted from a single sensor along the investigated radial line. However, within a practical setting, it is preferable to evaluate this for several virtual sensors and evaluate the mean so as to account for the fact that the size of the feature is small, but not 0. The only parameter that remains unknown is the phase $\theta \in [0,2\pi ]$, which can be inversely evaluated through use of an optimization scheme with only few sensors needed. For the results section of this paper a PSO employing 8 virtual sensors showed satisfying results. As Equation 14 is expected to approximate the response at a sensor location for a discrete frequency well, the PSO determines $\theta$ with good approximation. At this point, it is possible to to describe the response at various sensor locations along a line by only three parameters: the spatial frequency (${\omega }_{dist}$), the amplitude (A), and the phase shift ($\theta$). This allows to, firstly, significantly reduce the amount of stored training data and, secondly, to accurately extrapolate beyond the distance of the last sensor used for training.

As the response along a radial line due to a specific excitation can be accurately described by Equation 14, it is of interest to obtain the response due an arbitrary excitation. To achieve this, the concepts of Frequency Response Functions (FRFs) and transmissibility functions show promising results. These concepts, which are introduced in subsection 2.3 and described by Equation 6 and Equation 12, respectively, allow to evaluate the response at a sensor for inputs with different time–dependence. As stated above, the response for the sensors in Figure 4 due to the wave packets stemming from the feature, is found through baseline subtraction. Following Equation 11, the transmissibility function for each sensor with respect to the baseline can be numerically evaluated as:

$${\mathit{T}}^{bl}\left({\omega }_i\right)={\displaystyle \frac{{\mathit{u}}_u^{bl}\left({\omega }_i\right)}{u_k^{bl}\left({\omega }_i\right)}}$$

(17)

Where the superscript ${\left(·\right)}^{bl}$ denotes the baseline, ${\mathit{u}}_u^{bl}\left({\omega }_i\right)$ describes the frequency response at the sensors and $u_k^{bl}\left({\omega }_i\right)$ the frequency response in the feature location. Similarly for the case where the feature is present (denoted through the superscript ${\left(·\right)}^{fe}$:

$${\mathit{T}}^{fe}\left({\omega }_i\right)={\displaystyle \frac{{\mathit{u}}_u^{fe}\left({\omega }_i\right)}{u_k^{fe}\left({\omega }_i\right)}}$$

(18)

For a feature of size 0, the frequency response at the feature is the same for the baseline and the case where the feature is present ($u_k^{fe}\left({\omega }_i\right)=u_k^{bl}\left({\omega }_i\right)$). Therefore, the baseline subtracted transmissibility function (denoted through the superscript ${\left(·\right)}^{bs}$ is found as:

$${\mathit{T}}^{bs}\left({\omega }_i\right)={\displaystyle \frac{{\mathit{u}}_u^{fe}\left({\omega }_i\right)-{\mathit{u}}_u^{bl}\left({\omega }_i\right)}{u_k^{bl}\left({\omega }_i\right)}}$$

(19)

This enables evaluating the response at each sensor location, if the excitation at the actuator location is different or if the distance between the actuator and the feature varies.

In a similar way, the response corresponding to a line of sensors that are radially placed with respect to the actuator can be compressed and simplified. The only difference to a wave packet stemming from a feature lies in the application of the principle of FRFs rather than the principle of transmissibility functions. In this case, $\mathit{u}\left({\omega }_i\right)$ in Equation 5 reflects the response at the sensors and $f\left({\omega }_i\right)$ the excitation at the actuator.

At this point, it is possible to evaluate the response for any sensor along the radial line from the feature by performing the following steps:

• evaluate the response due to the direct wave

  1.1

  evaluate the FRF for the distance $d_{AS}$ from the actuator to the sensor

  1.2

  multiply the obtained FRF with the excitation

• evaluate the virtual response at the feature;

  2.1

  evaluate the FRF for the distance $d_{AF}$ from the actuator to the feature

  2.2

  multiply the obtained FRF with the excitation

• evaluate the response for reflected and mode-converted wave packets

  3.1

  evaluate the transmissibility functions for all modes of interest for the distance $d_{FS}$ from the feature to the sensor

  3.2

  multiply the obtained transmissibility functions with the response obtained in step 2.

• sum up the response obtained in step 1.2 and the response obtained in step 3.2

In what is described so far, each computed response corresponds to a sensor placed along a specific radial line. In order to more generally devise an approach to densely recover the response in the full field of an infinite, isotropic plate with a single feature, actuators along multiple radial lines can be employed covering a range of actuator-feature-sensor angles. The compression of each radial line of sensors into an amplitude (A) and a phase ($\theta$) proves advantageous in this respect. Using the same simulations that were previously used to infer the parameters (A and $\theta$) for a single line, an arbitrary number of lines can be compressed. Afterwards, the amplitude and the phase can be interpolated for each angle of interest. The limiting factors are the accuracy of the original simulation, the amount of data stored for the sensors, and the time needed for the optimization of the phase. These limitations do not hinter the evaluation and interpolation of a high number of parameters. Owing to the assumption of isotropic behavior, the evaluation for angles in the interval [0, π] is sufficient.

The offline process is summarized in Algorithm 1.

Algorithm 1: Offline phase

3.2. Online Phase

The compressed FRFs and transmissibility functions can now be used to rapidly obtain the response at any plausible sensor location within a plate. For simplicity, in a first step, a scenario will be investigated where only a single wave mode is propagating, e.g. the $S_0$–mode, in an infinite plate with two features $F_1$ and $F_2$ and a single sensor (Figure 5). Therefore, three points of interest exist to which a wave packet, excited at an actuator, will travel: the two feature and the sensor locations. We first investigate the wave from the actuator to the sensor. The response of the ray, representing this propagating wave packet, can be evaluated through the FRF determined in the offline part. Further, the two rays, pointing from the actuator to the features need to be investigated, as these initiate reflections at the features. When further tracking the ray, which points to feature $F_1$, reflections to two points of interest occur (dashed rays in Figure 5), where the first ray points to the sensor. To compute the sensor response of this path, three steps are necessary: (1) Compute the FRF response of the ray pointing from the actuator to the feature $F_1$, (2) Evaluate the transmissibility function for the ray pointing from the feature $F_1$ to the sensor and (3) Multiply this evaluated transmissibility function with evaluate FRF from step (1). The further reflection from the feature $F_1$ points to the feature $F_2$, where reflections again occur in two points of interest: the sensor location and once again back in the feature $F_1$ position (dotted in Figure 5). The response at the sensor location can be evaluated by multiplying the FRF and the respective transmissibility function belonging to the respective ray on this path. This ray tracing can be understood as a tree, which is displayed in 5(b). In theory, the branches of this tree can be tracked to infinity, however, the amplitude reduces with every reflection and the cumulative length of relevant rays. Moreover, propagation across longer distances requires more time for the wave to arrive at the sensor, which is eventually longer than the investigated time period. Therefore, from a practical viewpoint, it is sufficient to investigate a finite number of reflections.

Next, this concept should be generalized to a plate, where several modes occur. Therefore, a plate of infinite size comprising one actuator, one feature, and one sensor will be investigated, where the frequency of the excitation is chosen so that only the two fundamental modes ($S_0$ and $A_0$) propagate, as displayed in Figure 6. As two modes might be excited, for each point of interest, two rays need to be investigated, one for each mode. The two possible modes entail the possibility of mode conversion at the feature. This leads to a total of four wave packets propagating from the feature to the sensor. These wave packets pertain to the reflected $S_0$ and mode converted $A_0$ modes, stemming from the arriving $S_0$ mode, as well as the reflected $A_0$ and mode converted $S_0$ modes, stemming from the arriving $A_0$ mode. At the sensor location, this leads to a total of six incoming wave packets, which need to be summed. These wave packets can be determined through the evaluated FRFs for the direct waves and the multiplication of the FRFs for the waves propagating to the feature, with the corresponding transmissibility functions for the reflected and mode-converted waves.

3.3. Algorithmic Implementation

As shown in subsection 3.2, the ray tracing process can be represented as a tree structure. Within this tree structure, each feature can be seen as the start of a new tree, with nearly the same structure. This allows to implement the algorithm in a recursive fashion, as overviewed in algorithm 2. The complete algorithm for the online part, including the first call of the recursive function, is displayed in algorithm 3.

Figure 6. An infinite plate with one feature and one sensor. Possible travel paths for propagating $S_0$ and $A_0$ modes. Thin lines: $S_0$ mode; Fat lines: $A_0$ mode; Dashed: $S_0$ and $A_0$ modes, stemming from mode conversion at the feature.

Figure 6. An infinite plate with one feature and one sensor. Possible travel paths for propagating $S_0$ and $A_0$ modes. Thin lines: $S_0$ mode; Fat lines: $A_0$ mode; Dashed: $S_0$ and $A_0$ modes, stemming from mode conversion at the feature.

Algorithm 2: Recursive function

Algorithm 3: Online Phase

4. Results

The application of the described model to a synthetic example is described next. The performance of the surrogate is evaluated in terms of precision and speedup compared to the FOMs.

4.1. Full Order Model (FOM)

To train the described surrogate and to evaluate its efficiency, high fidelity simulations are necessary. Several spatial discretization methods are available for offering high-order Ansatz spaces that can effectively represent wave propagation according to Equation 2. See for instance the FEM with high polynomial degree (p-FEM) [74,75], the SBFEM [76], or IsoGeometric Analysis (IGA) [77]. The Spectral Element Method (SEM) [78], which is also adopted in this work, is often one of the preferred alternatives due to its optimal mass lumping [79], rendering the use of fast explicit solvers (i.e., integrators in the time domain) straightforward.

Upon selecting an appropriate discretization method, an external meshing software, such as gmsh [80], can be used to model the domain. If a complex structure and/or various damage configurations are present, this can become a tedious and time-intensive task, which can render automatized damage-detection algorithms infeasible. With this, as well as other meshing-related issues, in mind, several enhancements of the aforementioned methods have been proposed, which enable to supplement the mesh with additional geometrical definitions, thus easing the reliance on the performance of meshing software.

In the present contribution, the SEM was supplemented with a mesh-independent representation of voids, whose influence on Equation 2 is accounted for via customized space integration rules. This approach, known as the Spectral Cell Method [81,82,83,84,85], allows to maintain a level of efficiency that is on par with the SEM, however it also leads to a loss of the diagonal mass matrix property for elements that are intersected by voids. To mitigate this issue, several alternatives have been proposed, including diagonal scaling [86], or moment fitting [87] schemes, as well as use of implicit-explicit solvers [88,89]. In this contribution, we have adopted the latter approach, implemented in the C++ language. The corresponding code has been validated by comparison against a numerical implementation adopting shell elements [90], against benchmark analytical solutions for simple setups [87], as well as corroborated using Laser Doppler Vibrometer (LDV) experimental measurements [85].

4.2. Plate Configuration

For all given results, a 2 $\mathrm{m}$$\mathrm{m}$ aluminum plate with its material properties displayed in Table 1 is examined. As previously described, a regular mesh is used for all simulations, with features modeled as voids via customized integration. Each element of the mesh is 3 $\mathrm{m}$$\mathrm{m}$× 3 $\mathrm{m}$$\mathrm{m}$× 2 $\mathrm{m}$$\mathrm{m}$ with slightly smaller elements for sections whose dimensions are not a perfect multiple of 3 $\mathrm{m}$$\mathrm{m}$.

To avoid boundary reflections, Absorbing Layers using Increasing Damping (ALID) are applied. The ALID scheme refers to several layers of elements at the boundary of the domain, with each layer comprising a higher damping coefficient as the distance from the plate increases, which allows for absorption of the incoming waves. Here, we apply ALID through 30 layers, each with a width of one element, and an increasing mass proportional damping. The damping coefficients are determined through a quadratic function, with the maximum damping coefficient assigned as 10⁷ in the 30^th layer, following the recommendations in [91].

4.3. Excitation

The excitation in the given examples stems from two point loads applied at the top and the bottom of the plate in different locations. This introduces some practical questions. If both loads comprise the same magnitude and are oriented to the same or opposite direction, the plate is in theory excited purely anti-symmetrically or symmetrically, respectively. However, it should be noted that even in a frequency range, where only the fundamental Lamb wave modes ($S_0$ and $A_0$) are propagating, shear horizontal modes are activated. In the shown example this introduces negligible error and is, therefore, ignored. If the shear horizontal modes are of interest, it should be possible to targetedly excite these and consequently treat them similarly to Lamb Wave modes. We acknowledge that point loads are not thoroughly depicting realistic excitation sources within an actual monitoring setup. In this initial example, we follow the assumption that the excitation source, the features and the sensors have quasi zero size. While this could be modeled for the excitation and sensors, this is not possible for the features, which need to be modelled with finite size. The features in the presented example are therefore modeled with non-zero size, which introduces errors. However, these errors are small, if the ratio between the size of the feature and the distance to the feature is small. For a non point-wise excitation, the spatial dimension of the excitation will still be small compared to the size of the wave propagation distances and, therefore, also non point-wise excitations are expected to be possible. When seeking to model arbitrary excitation sources, i.e., not symmetric or anti-symmetric, we can always express these via superposition of a symmetric and anti-symmetric counterpart, as the behavior of the plate is expected to be linear-elastic.

Two different excitation scenarios are applied. Firstly, in a training phase a sudden impulse is applied as excitation, where the loading at the first time step is a non zero, while for all further time steps the loading is zero. This is possible, as an impulse induces a broadband excitation across all discrete frequencies, thus allowing for the derivation of the FRFs and transmissibility functions for all frequencies of interest. For the testing phase, a Hann-window pulse is applied, as commonly used in SHM-applications [13,28]. The Hann-window pulse is described as:

$$f(t)=\left\{\begin{array}{cc}(1-cos(2\pi t{\displaystyle \frac{f_c}{c}}))\left(sin\left(2\pi t{f}_c\right)\right)\hfill & if0<t<{\displaystyle \frac{c}{f_c}}\hfill \\ 0\hfill & else\hfill \end{array}\right.$$

(20)

with $f_c$ denoting the carrier frequency of the pulse and c the number of cycles. A Hann-window excitation in frequency domain shows a narrow main lobe with flat side lobes. The frequency range covered by the main lobe can be evaluated through [92]:

$$f_c\pm 2{\displaystyle \frac{f_c}{c}}$$

(21)

A narrow and distinct frequency range allows to deal with limited dispersion of the wave packet. For the following examples, the used Hann-windows will have a carrier frequency $f_c=300\mathrm{k}\mathrm{Hz}$ and $c=5$ cycles, resulting in a frequency range of interest between 180 $\mathrm{k}$$\mathrm{Hz}$ and 420 $\mathrm{k}$$\mathrm{Hz}$. In this frequency range the $S_0$–mode is nearly non-dispersive, while the $A_0$–mode is significantly dispersive (Figure 1).

Table 1. Material properties.

Table 1. Material properties.

4.4. Offline Phase

In the offline phase, the plate should be configured in order to learn the scattering behavior efficiently. The setup of the plate is shown in Figure 7, while the mesh details in the area of the features are illustrated in Figure 8. As two different features (a through hole and a notch) are investigated in what follows, a total of six simulations are needed. These simulations cover the three cases of a pristine plate, a plate with a through hole, and a plate with a notch, where for each case, a simulation with symmetric and anti-symmetric excitation, respectively, is needed. In total, a number of 808 sensors is assumed as positions where the response is to be reconstructed by the delivered ROM. More specifically, we assume the sensor placed in an arrangement of eight circles around the feature, with 101 sensors per circle.

For these simulations, the time step size is set at $2.5\times {10}^{-8}$ $\mathrm{s}$, slightly less than the critical time step for the regular mesh. The simulations are run for 7999 time steps or a total of $2\times {10}^{-4}-2.5\times {10}^{-8}\mathrm{s}$. The total simulation time is chosen so that the slowest wave in the investigated frequency range reaches the outermost sensors. With these settings, a frequency step of 5 $\mathrm{k}$$\mathrm{Hz}$ is achieved, where the first 160 discrete frequencies (up to 795 $\mathrm{k}$$\mathrm{Hz}$) are investigated.

Based on these six simulations, a dictionary is formed. For each of the 101 radial directions and each of the 160 frequencies a sine is fitted through the eight radial points, with the phase found through a PSO. In a second step, the obtained amplitudes and phases are linearly interpolated through the circumference.

In Figure 9 Figure 10, and Figure 11 the scattering characteristics and accuracy of the created dictionary are assessed. Firstly, Figure 9 shows the scattering amplitude of the notch for different angles for the frequency of 400 $\mathrm{k}$$\mathrm{Hz}$. It is visible that for some angles the scattering is very low. Figure 10 shows an instance of the FRF for an angle of 126° between the incident and the scattered wave. The decaying periodic function, described in Equation 14, is clearly visible. The surrogate accurately represents the behavior of the FOM, despite being trained only at 8 points along a line. Only for some points close to the notch larger deviations can be observed. The error for individual training points can be seen in Figure 11. In comparison with Figure 9, it should be noticed that lines of high errors often correspond to angles of low amplitude and, therefore, a small denominator, when computing relative errors.

4.5. Online Phase

In Figure 12, a sketch of the plate configuration mesh used to test the online phase of the proposed approach is illustrated. This setup reflects a square plate with a single actuator, two features and seven sensors. The origin of the coordinate system is assumed to coincide with the position of the actuator with a distance of $0.05$ $\mathrm{m}$ to the left, $0.15$ $\mathrm{m}$ to the right, and $0.1$ $\mathrm{m}$ to the top and bottom boundary of the plate. The features are located at ($0.05$ $\mathrm{m}$| $0.04$ $\mathrm{m}$) and ($0.07$ $\mathrm{m}$| $-0.02$ $\mathrm{m}$). Around the plate there are 30 layers of ALID, to prevent boundary reflections.

We first examine the scenario where Feature 1 corresponds to a hole (as shown in Figure 8(a)) and Feature 2 is inactive. This means that a pure reflection without mode conversion is expected. Further, at most sensor locations two wave packets are expected - the direct and the reflected wave packet - which may or may not overlap. For the sensor at location (0.9|0.9), the wave packet stemming from the feature is expected to result in a direction that is opposing the direct wave, as the sensor lies in the shadow of the feature. The plots for the sensors at locations (0.05|-0.08) and (0.09|0.09) are shown in Figure 13 and Figure 14, respectively.

It can be seen in Figure 13, that there is an additional wave packet for both symmetric and antisymmetric excitation, and that the additional wave packets are captured well by the surrogate, both in amplitude and in phase. For the sensor at location (0.09|0.09), which is placed in the shadow of the hole, it can be seen that the response without a hole would have a higher amplitude than for the case with a hole; an effect which is well captured by the surrogate. However, in the case of anti-symmetric excitation (Figure 14b), traces of a periodic signal are discernible in the response of the surrogate. These are artifacts, which likely result from a poorly fitted frequency for a certain angle. This seems tolerable, as the overall response is fitted well.

Next, in Figure 15 and Figure 16 the response is displayed for the case where Feature 1 is a notch rather than a hole. Figure 15a allows to clearly distinguish three arriving wave packets, compared to the two, which were visible in the case of the hole (Figure 13a). Thereby, the additional, third wave packet stems from mode conversion from the incident symmetric wave packet into an anti-symmetric wave packet. Further, it can be seen that the amplitude of the reflected symmetric wave packet is smaller in this case, since the notch only reduces the cross section, while the hole goes through the whole cross section. On the other hand, the reduction of the cross section results to a minor drop in the amplitude of the direct wave for the case of the sensor being in the shadow of the feature (Figure 15) compared to the case with the hole. All these effects are captured well by the surrogate. For the sensor located in location (0.09|0.09), the wave packets originating from the notch are hardly distinguishable from the direct wave (Figure 16a). However, the surrogate still matches the behavior predicted by the FOM well. Although the effect of the notch is clearly visible for the case of symmetric excitation, there is hardly any effect for the case of the anti-symmetric excitation. This is expected due to the high difference of the amplitude for symmetric and anti-symmetric modes for the out-of-plane displacement (here approximately two orders of magnitude).

For the plots in Figure 17, both features are modeled as holes. Focusing on the case of anti-symmetric excitation in Figure 17b, four wave packets are now visible. The first wave packet stems from the direct wave, while the second and third wave packets correspond to the reflections from the two holes. The fourth wave packet stems from a wave path, which includes reflections at both holes (i.e. Actuator - Feature 1 - Feature 2 - Sensor or Actuator - Feature 2 - Feature 1 - Sensor), while the wave packet from the other path that includes both wave packets is not clearly visible. This may be attributed to an overlap with one of the other wave packets or to the amplitude for this path dropping significantly, to the point that it is not visibly discernible. Finally, the results for the case where feature 1 corresponds to a hole and feature 2 is a notch can be seen in Figure 18, where again the main wave packets are captured in amplitude and phase.

4.6. Speedup

All surrogate simulations were performed on a single Intel(R) Xeon(R) Gold 5318Y CPU @ $2.10$ $\mathrm{G}$$\mathrm{Hz}$ core using Matlab 2024a [93], while the FOM simulations were run on the ETH Euler cluster with eight core AMD EPYC 7742 processors (2.25 GHz nominal, peak of 3.4 GHz). The surrogate simulations with a single feature take approximately $0.5$ $\mathrm{s}$ and with two features approximately $4.4$ $\mathrm{s}$ for the seven sensors displayed in Figure 12. On the other hand, the computational time for the FOM is nearly independent of the features, corresponding to $1.6$ $\mathrm{h}$ for each simulation. For the cases considered and when using the described hardware, this results in a speedup on the order of ${10}^3$ to ${10}^4$ . In addition, tests including three and four features were performed, where the simulations were run for approximately 16 $\mathrm{s}$ and 49 $\mathrm{s}$, respectively, resulting in a speed-up in the order of ${10}^2$ . Generally, the computational time for the surrogate depends heavily on the number of features that are assumed present, but only indirectly depends on the size of the plate (through the attenuation and the propagation time). Conversely, computing the FOM is heavily dependent on the size of the plate, while the effect of the features on the numerical toll is only minor albeit implies a required mesh modification.

5. Discussion

We introduce an approach to efficiently model Guided-Wave based sensing scenarios by delivering a surrogate for rapid evaluation of Lamb Wave propagation. The method can flexibly account for diverse geometrical features and is easily extendable to more complex setups. In this study, which serves more for introduction of the approach, we place specific assumptions, namely:

• The plate is assumed to comprise a linear-elastic and isotropic material.
• The area where the excitation is applied, the features and the sensors are considered small compared to their interconnecting distance and can, therefore, be assumed as point-wise.
• The explored frequency range is chosen so that only fundamental modes are propagating and only the out-of-plane displacement is investigated.
• No viscous damping is assumed.

Most of these assumptions can be relaxed, while others are allowed for the usual cases explored within the NDE and SHM contexts, given the small amplitude of actuated Lamb waves. The method introduces a new approach to modeling GW-based assessment, featuring the following advantages:

• The computational cost of the surrogate is significantly lower to the FOMs and therefore suited for inverse problem settings, as mandated in condition monitoring and damage diagnosis tasks. In the provided examples the speedup is in the order of ${10}^3$ to ${10}^4$ .
• The method can be extended to account for further intricacies, including geometries that accommodate diversified plate setups. The authors expect that an extension to anisotropic plates is possible with more training simulations.
• The full time history at the sensor location is reconstructed, and not only the arrival time or the amplitude, as furnished by commonly adopted alternatives.
• The mechanics behind the surrogate is well understood.

The rooting of the method in mechanics principles reduces the training effort, with only few simulations required for training, while ensuring that the model estimates can be extrapolated. Further, the process allows to track single wave packets through various possible propagation paths, which allows to gaining information on the individual wave packets present in the sensor response. This allows a better understanding of the wave propagation phenomena, which is superior compared to the classical FOMs.

Building on the proposed method, several promising avenues for further research can be explored. While the current approach focuses on circular defects, a modified training strategy could enable its application to more complex geometries, significantly expanding its versatility. To further enhance efficiency, investigating the relationship between scattering wave characteristics and defect size could reduce the need for extensive training simulations. Notably, Diligent et al. [63] observed that for holes, the amplitude of scattering waves is nearly proportional to their diameter—an insight that could be leveraged to refine the model. For real-world applications, incorporating boundary reflections is essential. A key step in this direction is analyzing the interaction between different wave modes at the boundary, such as shear horizontal and anti-symmetric modes [51]. Additionally, extending the method to anisotropic materials presents an exciting opportunity to account for direction-dependent wave propagation speeds, further broadening its applicability. Optimizing sensor placement also holds great potential. While the current offline approach employs a regular grid, a strategically optimized placement could minimize errors, as observed in Figure 14b. Finally, enhancing computational efficiency through parallelization and adaptive stopping strategies would further accelerate the surrogate model, making it even more effective for practical use.