Linear and Nonlinear Guided Wave Imaging of Impact Damage in CFRP Using a Probabilistic Approach

The amount and variety of composite structures that need to be inspected for the presence of impact damage has grown significantly in the last few decades. In this paper, an application of a probabilistic ultrasonic guided wave imaging technique for impact damage detection in carbon fiber-reinforced polymers (CFRP) is presented. On the one hand, a linear, baseline-dependent, technique utilizing the well-known correlation-based RAPID method and an array of piezoelectric transducers is applied to detect impact-induced damage in plate-like composite structures. Furthermore, a baseline-independent nonlinear extension of the standard RAPID method is proposed, and its performance is demonstrated both numerically and experimentally. Compared to the conventional RAPID, the baseline-free version suffers from a somewhat lower imaging quality. However, this drawback is compensated by the fact that no damage-free (intact) baseline is necessary for successful imaging of damage.


Introduction
Conventional nondestructive testing (NDT) methods, such as visual inspection, ultrasonic C-scan, thermography or shearography, are well established and frequently used in the field of quality control and defect inspection. Unfortunately, as one can imagine from an example, such as an ultrasonic C-scan of a composite aircraft wing, these methods are quite bulky and time consuming. Unlike point-by-point inspection techniques, guided wave (GW) techniques used for material characterization and evaluation are driven by their ability to inspect large areas using only a limited amount of ultrasonic transducers. Therefore, the use of GW-based NDT could significantly decrease the inspection time and costs compared to conventional NDT methods. Moreover, the transducers within the transmitting-receiving network can be kept in place during the operational lifetime of the component, because they are small enough not to influence the mechanical performance of the component. Thanks to the improvement in technology, the durability of such systems becomes more and more recognized. Therefore, GW-based methods represent a perfect choice for structural health monitoring (SHM) and defect imaging.
Lamb waves are one type of guided waves usually associated with elastic wave propagation in thin plates. First described in 1917 by British mathematician Horace Lamb in [1], they have drawn attention ever since. Worlton [2,3] was perhaps the first author who recognized the potential of Lamb waves for NDT. Shortly after, Grigsby [4] described and highlighted the most important properties of Lamb waves with respect to the inspection of thin plates. However, the modern era of Lamb wave research and practical applications dates back to the influential book written by Viktorov in 1967 [5].

Guided Wave Imaging Using RAPID
Conventional RAPID, as described by Gao et al. [27], utilizes data from an ultrasonic sparse array consisting of n e permanently-attached PZT transducers. The direct line-of-sight coverage of such an array is visualized in Figure 1. In the conventional RAPID, signals between each transmitter-receiver pair are acquired in two different states: baseline (intact, without damage) and damaged (after the damage has been introduced) [28]. Let the signal transmitted from the sparse array element i to element j be denoted B ij and D ij for the baseline and damage state, respectively. If the sparse array consists of n e elements, then the total number of acquired signals is n e (n e − 1). The total number of signals can be reduced down to n e 2 (n e − 1), if the reciprocity of the system (B ij = B ji ) is assumed [17,30]. Fundamentally, the RAPID algorithm can be broken down into two major parts: the SDC calculation and the imaging. The SDC values represent a measure of the dissimilarity between each two signals obtained in two different states, i.e., between B ij and D ij . Even though being the second step in the algorithm, the imaging part will be described first, since it is general and common for both the baseline-dependent (linear) and baseline-free (nonlinear) version of the algorithm. Let us first thus assume that the SDC ij coefficients are known for all T-R pairs. We start by dividing the inspected area of the sample into a rectangular equidistant mesh. Next, we define a priori probability distribution s ij (see Figure 2b) for each T-R pair and every point [x, y] of the mesh as follows: where β stands for a free-to-choose shape factor that defines the area influenced by one T-R pair [33]. R ij (x, y) is a geometrical function defined as: which is representing the ratio of the distances |TP| + |PR| and the |TR|. As a result, the numerator β − R ij (x, y) from (1) actually describes an ellipse with focal points located in T, R and a major axis length β |TR| 2 , as illustrated in Figure 2. This ellipse forms the borderline of the region influenced by a particular T-R pair. If a selected point P lies within the area defined by the borderline, its s ij decreases with the distance from the straight line that connects points T and R. In other words, the further the point P from the direct path between the transmitter and the receiver, the lower its s ij value will be. The s ij is evaluated for all T-R pairs and at all points of the rectangular grid. Once combined with the individual SDC ij values, the final damage index heat map is calculated as: where contributions from all T-R pairs are summed up at a particular grid point with the SDC ij values serving as geometrical weights. A cluster of points with a high damage index P ij then indicates the potential location of a flaw in the sample.

Linear RAPID
In the case of the conventional implementation of RAPID, the SDC ij is calculated using the standard correlation coefficient defined as: where σ(B ij ) and σ(D ij ) are the signal variances and: defines the covariance of the baseline and the damaged signals. Values µ ij denote the mean value of the corresponding signal between transmitter i and receiver j, and the index [k] simply indicates that all signals are discretely sampled at a sampling rate of f s = 1 ∆ t (∆t is the sampling interval). The SDC ij value that quantifies the difference in signals can then be calculated using ρ ij as: The lower the correlation between the baseline and damaged signals, the higher the SDC coefficient will be and the more profound the dissimilarity is between the two signals.
The basic idea behind the use of SDC values for imaging purposes can be illustrated with the following example. Imagine that a damage has been introduced to a sample after the acquisition of baseline signals. The presence of the damage triggers changes in the wave propagation through the sample, and if the damage lies on or close to a direct path between a selected T-R pair of the network, the wave propagation will be altered significantly, resulting in a rise of the corresponding SDC ij (correlation decreases). However, if the damage is sufficiently far from the sensor pair, its influence on the wave propagation will be negligible, and the corresponding SDC ij value is low. Thus, by combining the input from all T-R pairs, the damage index P ij map can be constructed, and the damage can be visualized.

Nonlinear RAPID
The main advantages and disadvantages of the conventional RAPID methodology, outlined in the previous section, were already presented in the Introduction. It is clear that the weakest spot of this technique is the baseline signal selection process. Changes in the environmental conditions, at which the B ij and D ij were measured, can result in severe deterioration of the imaging quality. Hence, in this section, we propose a method that transforms the conventional RAPID into a baseline-free method by introducing a SDC parameter based on the scaling subtraction method (SSM) [34,35]. This will eliminate the requirement for an intact baseline signal, which should improve the real-world applicability of RAPID.
First, some major assumptions have to be presented before the method itself is described. It is assumed that: • the defect behaves nonlinearly with increasing excitation amplitude, meaning that it responds differently to a low amplitude excitation than to a high amplitude excitation due to the existence of a mesoscopically nonlinear stress-strain behavior at the defect location. • the nonlinearity of the equipment and transducers is low.
These assumptions have to be made, otherwise the new localization concept would make no sense. If the nonlinearity of the transmitter and receiver is higher than the nonlinearity of the defect, the imaging would simply fail. Since the nonlinearity of the transducer response depends on the excitation amplitude, this has to be kept sufficiently low not to provoke a distortion of the transmitted signal. This acceptable level can be found experimentally by direct monitoring of the output waveform distortion as a function of input amplitude.
The details of the baseline-free RAPID method are as follows. First, a low amplitude excitation signal is applied to the sparse array transducers, and the corresponding response signals B ij are obtained. The response signals B ij will act as reference (defect-free) signals. Next, a high amplitude excitation is applied, and the response signals D ij are acquired under the same measurement conditions. Both sets of measurements can be obtained at the same time (a matter of seconds difference). Hence, the influence of the environmental conditions, such as temperature variation or mechanical loading, is negligible. The only difference is that the excitation amplitude has been up-scaled by a scaling factor: where A B , A D are the amplitudes of the excitation signals for low and high amplitude, respectively. The new, nonlinearity-based, SDC coefficient is given by the mean square difference between the high amplitude response signals and the up-scaled low amplitude response: If the inspected system is purely linear, the response signals scale up perfectly and are simply equal to k s · B ij . Hence, the value of the SDC will be equal to zero in this case. The imaging will result in a blurred unfocused image without any defect indication. However, if a nonlinear defect is present in the interrogated sample, the SDC value attains a non-zero value for some of the T-R paths.
Based on this principle, the SDC in (3) can be simply replaced by Equation (8), and the new RAPID method basically becomes baseline-independent, because it does not need an input from a state measured before the damage took place.

Correction for Direct Path Propagation
The original RAPID algorithm brings about a second problem. For the imaging to be successful, the SDC value for each T-R pair should be calculated only from the part of the recorded signal that corresponds to the direct propagation path (extended by parameter β) between the transducers. This can be easily done in isotropic materials, but it is considerably more difficult to fulfil this restriction in orthotropic samples.
In order to restrict the useful signal and include only the direct propagation, a simple manipulation can be suggested for the received signal. Assuming an anisotropic plate, let us define the distance between the transmitter and receiver by x and the angle of propagation with respect to the coordinate system of the anisotropic plate by θ. Furthermore, we assume that the entire recorded signal contains n s samples. To assure that only the direct propagation is taken into account, a reduced amount of samples n r can be considered. Here, n r (n r < n s ) can be easily calculated from an estimate of the TOF between transmitter i and receiver j as follows: where f s is the sampling frequency and n c is the number of cycles in the excitation waveform at a frequency f . TOF ij is the time-of-flight defined as: where c ph (θ) is the phase velocity of the selected GW mode in the direction given by the angle of propagation θ in the anisotropic plate. As can be seen from (9), a number of cycles n c that is equal to the arrival time of the center of the wave packet is taken into account. The phase velocity can be determined from numerical dispersion data, provided most of the elastic constants of the sample are known. This approach is used to ensure a proper imaging strategy in anisotropic materials.

Damage Index Image Processing
The resulting damage index P(x, y) heat map can be improved by thresholding the SDC values prior to being used in the imaging step. In order to do this, the SDC values are first rescaled to a closed interval [0,1]. Next, the values below a specified threshold are set to zero, and the remaining ones are left intact. Clearly, the number of contributing T-R pairs that form the image will drop, and only the most influenced ones remain. In a more advanced manner, the value of the threshold t SDC can be determined empirically or linked to the mean value and standard deviation of the initial SDC values. The general description of this thresholding operation is: Finally, in order to facilitate an easier visual perception, the resulting image can be simply thresholded at a fixed level t P : thereby creating a binary image that immediately pinpoints the areas with the highest damage index P ij . If a fixed threshold is applied to the result, the size and location of the suspected defect(s) can be easily estimated.

Numerical Simulations
Numerical simulations were carried out in order to verify the proposed baseline-free RAPID methodology. The simulations were performed using COMSOL Multiphysics c (COMSOL AB, Stockholm, Sweden).

Model
The test sample was simulated as a simple bulk square CFRP T300/924C plate with orthotropic symmetry and dimensions of 288 mm × 296 mm × 2.7 mm (see Figure 3a). The density of the material is ρ = 1548 kg m −3 , and its elastic properties are summarized in Table 1. The simulation model includes damping that was implemented using a Rayleigh damping model, which is a standard included in COMSOL Multiphysics c . In our case, the mass damping parameter α dM of the model was set to zero, and the stiffness damping parameter: was defined based on the quality factor Q. The quality factor was chosen to be constant (Q = 20) for all simulations.   Figure 4 shows the dependency of the phase velocity on the angle of propagation for the fundamental GW modes at a fixed frequency of 50 kHz (calculated using the Legendre polynomial approach [36]). These data will be used to carry out the direct propagation correction.

Young's Modulus [GPa] Shear Modulus [GPa] Poisson's Ratio
For the simulations, the sample contained either one or two (depending on the configuration) 20 mm × 20 mm nonlinear delaminations located at 1/4th of the plate's thickness (see Figure 3a). Each defect was simulated as a clapping system (kissing bond) based on the nonlinear stress-strain model developed by Delrue and Van Den Abeele [37]. The graphical representation of this model is illustrated in Figure 3b. The nonlinear dynamic behavior of the delamination is controlled by a set of virtual spring and damper forces at both sides of the delamination, as illustrated in Figure 3b. These forces can be expressed as functions of the gap distance ∆Z. Above a certain separation threshold, the two sides are completely separated (stress free surfaces); below the threshold, particular formulations of the van der Waals forces are implemented. When the surfaces are close to each other, they will be attracted to each other. However, when they tend to be too close, the attraction force will turn into a repulsive force, trying to separate the surfaces again. The piecewise behavior of these additionally introduced elastic contact forces for a kissing bond is illustrated in Figure 5. Apart from the elastic contact force, damping forces were also implemented, which are acting against the velocity of separation. These forces make sure that the surfaces of the delamination are not opening too abruptly, avoiding a destruction of the material. At the same time, they assure that the surfaces are not closing to fast so that the two surfaces cannot overlap. These forces are shown in Figure 3b. The nonlinear viscoelastic behavior at the delamination level mimics the clapping behavior of a delamination. Depending on the displacement amplitude of the wave passing by the defect, this nonlinearity will be activated or not, creating a distortion in the wave propagation, which can be measured in the received signal after appropriate signal processing. The implementation of the nonlinear spring-damper forces is performed by introducing dynamic boundary conditions in COMSOL Multiphysics c , at those nodes that correspond to the delamination surface. At these positions, the nodes are split in pairs, and the following analytical formula's are implemented for the spring and damper forces [37]: where Z 0 is a small characteristic distance between the faces of the delamination, a, b are free parameters defining the separation conditions (1 < a < b), γ is a damping coefficient and v z t and v z b stand for the normal velocities of the top and bottom interfaces. k 1 , k 2 and k 3 are the virtual spring constants that are connected by: Details on the parameter values for the nonlinear delamination model that was used in the numerical simulation are summarized in Table 2.

Results
In order to test the viability of the proposed baseline-free modification, four numerical evaluate cases were defined. In the test cases, we test the performance of the imaging algorithm using a variable number and locations of delaminations and different forms of excitation. The test cases are summarized in Table 3. The imaging was carried out using a sparse network of eight virtual transducers placed in a rectangular array. The locations of the transducers are depicted in Figure 3a. The excitation waveform consisted of a 20-cycle (n c = 20) Hanning windowed sine burst at 50 kHz. The excitation was implemented as a prescribed normal displacement at the area of the transmitting virtual element. The received signal was sampled with a sampling frequency of f s = 10 MHz. The adaptive signal length correction was applied using the dispersion data from Figure 4. Figures 6-9 show the results of the baseline-free RAPID detection for a CFRP plate containing a single or two nonlinear delaminations with the properties described in the previous subsection. Except for Test Case 3, all delaminations were successfully detected, as can be clearly seen from the binarized images. In Test Case 3, due to the mutual position of the sparse array elements and the defects, it was necessary to lower the threshold by 5% in order to successfully detect both delaminations (see Figure 8c). However, it can be easily observed from Figure 8a that both defect zones are clearly highlighted.
The fourth and last test case demonstrates the effect of signal duration (number of sine cycles) on the imaging performance. The test case is similar to the first test case with one exception. The number of cycles in the excitation signal was decreased from n c = 20 to n c = 3. The effect of this change can be observed in Figure 9. It is obvious from a comparison between Figures 6 and 9 that the imaging output is different for different values of n c . A higher number of cycles clearly improves the performance of the baseline-free imaging algorithm.   [30,−70]. The binarization threshold was set to t P = 0.8, the shape factor to β = 1.015 and the SDC threshold to t SDC = 0.25. The actual location of the delamination is marked with a white square in (a) and a red square in (b).
The main reason for the inferior imaging performance at low n c is the fact that a shorter signal is inherently more broadband, and the more broadband the signal is, the higher the dispersion will be. As a natural consequence, we obtain a less precise estimation of the TOF and a poorer adaptive signal length correction. On the other hand, a larger n c instigates unavoidable problems with reflections coming back from the edges of the plate. Therefore, an optimal n c that avoids reflections and provides narrowband excitation has to be found for every application, based on the sample geometry and the sample's dispersion characteristics.  The binarization threshold was set to t P = 0.8, the shape factor to β = 1.015 and the SDC threshold to t SDC = 0.25. The actual location of the delamination is marked with a white square in (a) and a red square in (b).
Another reason may be linked to the behavior of the nonlinear delamination and the related SDC value. At low n c , only a very small portion of the signal is distorted due to the nonlinear delamination, because the faces of the delamination undergo just three cycles instead of 20. Therefore, the nonlinear effects influence the result only slightly at n c = 3, whereas they have more time to build up for n c = 20. According to Delrue and Van Den Abeele [37], the combination of SSM and the nonlinear delamination model performs very well for larger n c ≈ 40. The rigorous study of the mode behavior for the small n c has not been published yet.

Experimental Results
The above discussed linear and nonlinear versions of RAPID were further experimentally verified on an orthotropic CFRP plate.

Test Sample
The test sample to be investigated is a 400 mm × 400 mm CFRP plate consisting of 11 separate layers/plies (see Figure 10a). The stacking sequence of the laminate, as well as the ply type are given in Table A1. The laminated composite can be concisely described using the standard stacking sequence notation as [(0, 45) 2 , 90, −45, 0, −45, 90, 0] T [39]. The elastic properties of the plies are collected in Table A2. The total thickness of the plate is approximately h tot = 4.3 mm. Due to the (non-symmetric) stacking sequence, the global (homogenized) in-plane elastic properties of the plate cannot be assumed simply quasi-isotropic. Isotropy can only be assumed for the lowest antisymmetric mode A 0 , as can be seen in Figure 10b. The angular profile of c ph for the A 0 mode is nearly circular. The phase velocities of the other modes, such as S 0 and SH 0 , are strongly dependent on the direction of propagation. This is a fact that has to be kept in mind during the damage location or beam-forming imaging of laminated plates with non-symmetric layups.

Experimental Setup
For the experimental validation of RAPID, an experimental setup for GW imaging was developed, consisting of a single channel arbitrary waveform generator (AWG) NI PXI-5412 (100 MHz, 14-bit, single channel, National Instruments Corporation, Austin, TX, USA) that was routed to a multiplexing unit via an AR150A100B amplifier (Amplifier Research, Souderton, PA, USA). The signal from the switch board was connected to the individual PZTs of the sparse array network using a 50 Ω shielded coaxial cable. The input from the amplifier is connected to the transmitting element, while the other elements are connected to the receiving digitizers: NI PXI-5122 (100 MHz, 14-bit, two channels, National Instruments Corporation, Austin, TX, USA). All transmitting and receiving cards were hosted in a single chassis NI PXI-1022 (National Instruments Corporation, Austin, TX, USA). The transmitting element is disconnected from the receiving stage during the excitation to avoid leakage and crosstalk of this signal to other channels. The PZT transducers used as sparse array elements were DuraAct c P-876 flexible PZTs (PI Ceramic GmbH, Lederhose, Germany), in either rectangular or circular shape, covered in a protective polymeric layer. The entire measurement was controlled using a virtual instrument (VI) program written in LabVIEW c .

Linear Conventional RAPID
In order to illustrate the conventional baseline dependent RAPID, a test was conducted on a CFRP plate. It was equipped with 16 rectangular DuraAct c transducers (PI Ceramic GmbH, Lederhose, Germany) that were coupled using Hysol c adhesive (Henkel AG, Düsseldorf, Germany). The baseline signal was recorded at standard room temperature in laboratory conditions using the experimental setup described in Section 4.2. It was subsequently damaged by a 21 J impact in the upper right corner (see Figure 11c for the result of a traditional C-scan analysis after impact). Figure 11. Baseline-dependent RAPID imaging on a CFRP plate, excitation waveform: five-cycle sine burst f = 50 kHz, β = 1.015 and t SDC = 0.5 (a) damage index P(x, y), (b) thresholded binary RAPID image using a threshold level at t P = 0.9 max{P(x, y)} (c) ultrasonic C-scan of the plate.
Before and after the impact, with a difference of more than a month, the sample was excited using a five-cycle sine burst at 50 kHz with a Hanning window. According to Figure 10b, only the A 0 and S 0 modes are excited at this frequency range. Moreover, using the experimental dispersion analysis, it was found that the A 0 mode is the most dominant mode at this excitation frequency. The received signals were sampled at 10 MHz, and n s = 16,384 samples were acquired. The input amplitude to the amplifier was 0.1 V. The "current" set of signals D ij was acquired using the same test stand at ambient room temperature without any precise temperature adjustments.
For the subsequent analysis, the length of the signal between every T-R pair was restricted to include only the direct propagation of the A 0 mode [32]. This was easily done based on the angular dispersion data and the phase and group velocities from Figure 10b. Hence, knowing the propagation distance, direction and velocity, the signal can be properly time windowed.
The resulting damage index image is depicted in Figure 11a. The impact zone is well indicated in the binarized image (Figure 11b).

Nonlinear Baseline-Free RAPID
The experimental setup for the nonlinearity-based baseline-free RAPID validation was the same as described in Section 4.2. The lower input amplitude to the amplifier was set to A lin = 0.01 V, and the SSM scaling coefficient was k s = 10. The actual voltage input to the PZT was approximately 8 V and 80 V for the lower and higher amplitude excitation, respectively. The excitation frequency was f = 50 kHz, and the waveform was a three-cycle Hanning windowed tone burst. Such a low cycle count was selected due to the position of the defect and due to the shape of the sparse array. A longer duration of excitation would result in a mixing of directly propagating signals and signals that are reflected from the edges of the sample. Hence, the short excitation is unavoidable in this case. The signals were sampled with a sampling frequency f s = 10 MHz, and each signal was averaged 512 times to get a good SNR for the low amplitude excitation. Figure 12 shows the result of the nonlinear RAPID imaging of the impact damage, applying the following imaging parameters: t SDC = 0.15, β = 1.015, t P = 0.8. The defect was satisfactorily highlighted in the damage index map and subsequently segmented using the fixed threshold binarization.

Discussion and Conclusions
Based on the linear and nonlinear measurements and subsequent analysis, we come to the following conclusions: First, the conventional linear RAPID, employing correlation-based SDC coefficients, proved to be a considerably reliable imaging technique when used in the laboratory conditions. The predicted location of the damage corresponds very well with the ultrasonic C-scan data depicted in Figure 11c. On the other hand, the size of the defect is somewhat overestimated. However, it is a very satisfactory result especially considering the extended period between which the baseline and damaged signals were taken.
Second, the nonlinear baseline-free RAPID was able to detect the same defect using SSM-based SDC coefficients. Important to note here is that the coefficient and the baseline signals were acquired after the impact damage was introduced in the plate and that no reference to an intact state is required. Figure 12 suggests that most of the nonlinear effects are generated when waves propagate through the damaged region in the y-direction. This can be attributed to the distribution of the partial damage caused by the impact along the dominant material symmetry axes (0 • , 90 • ). However, the shape of the damage distribution inferred by the nonlinear RAPID method does not entirely conform with the result from the ultrasonic C-scan. Therefore, another technique with higher spatial resolution, for instance uCT, should be used to confirm whether the micro-cracking is indeed present in a wider area than predicted by the C-scan.
Compared to the conventional RAPID, the baseline-free version suffers from lower imaging quality. However, this drawback is compensated by the fact that no damage-free (intact) baseline is necessary for successful imaging of damage. This promising result can be further investigated and developed in view of achieving more robust and more environment interference-resistant SHM methods in the future.
Author Contributions: J. Hettler and M. Tabatabaeipour conceived of the idea of baseline-free RAPID and designed the experiments. S. Delrue performed the numerical simulations, and K. Van Den Abeele contributed with a critical analysis of the data.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: