Lamb-Wave-Based Multistage Damage Detection Method Using an Active PZT Sensor Network for Large Structures

A multistage damage detection method is introduced in this work that uses piezoelectric lead zirconate titanate (PZT) transducers to excite/sense the Lamb wave signals. A continuous wavelet transformation (CWT), based on the Gabor wavelet, is applied to accurately process the complicated wave signals caused by the damage. For a network of transducers, the damage can be detected in one detection cell based on the signals scattered by the damage, and then it can be quantitatively estimated by three detection stages using the outer tangent circle and least-squares methods. First, a single-stage damage detection method is carried out by exciting a transducer at the center of the detection cell to locate the damaged subcell. Then, the corner transducers are excited in the second and third stages of detection to improve the damage detection, especially the size estimation. The method does not require any baseline signal, and it only utilizes the same arrangement of transducers and the same data processing technique in all stages. The results from previous detection stages contribute to the improvement of damage detection in the subsequent stages. Both numerical simulation and experimental evaluation were used to verify that the method can accurately quantify the damage location and size. It was also found that the size of the detection cell plays a vital role in the accuracy of the results in this Lamb-wave-based multistage damage detection method.


Introduction
Among structural health monitoring (SHM) techniques, Lamb-wave-based methods have shown great potential for fast damage detection in plate-like structures [1]. Lamb waves can propagate a distance with low attenuation and are highly sensitive to small imperfections [2,3]. These waves have been significantly utilized for SHM and determination of surface defects in large metallic structures [4]. Accurately estimating the damage not only requires the construction of damage images for its localization but also an accurate quantification of both the location and size of the damage. The presence of damage in a structure results in the forward and backward scattering of waves. The change in Lamb wave propagation due to this reflection and transmission can be detected and analyzed to evaluate the damage [5]. There are a variety of methods to excite and receive Lamb wave signals. These methods can be grouped into five broad categories based on the transducers [6]: ultrasonic probes, lead zirconate titanate (PZT) transducers and piezocomposite transducers, laser-based ultrasonic transducers, interdigital transducers (IDTs), and fiber-optic transducers. Comparing these methods, PZT transducers can be directly mounted on the surface of the target structure to excite/receive is introduced in this research to accurately quantify not only the location but also the size of the damage. The method introduces a strategy to arrange PZT transducers in the form of a network of detection cells on a large plate-like structure. The CWT method using the Gabor wavelet is applied to analyze the dynamic response caused by the damage. The damage location and size are evaluated in three stages through an algorithm based on the outer tangent circle method. A single-stage method provides an initial estimate of the damage size and locates the damaged cell. It further locates the damaged subcell inside the damaged cell. The estimation of damage location and size inside the damaged subcell significantly improves in the second and third stages of detection. The arrangement of transducers and the signal-processing method remain the same in all stages of damage detection. The method does not require the baseline signal from the undamaged plate, and the improvements in the subsequent stages can be achieved without the need for additional transducers. These features make the method less complicated, low-cost, and less time consuming.
The organization of this article is as follows: Section 2 introduces the theoretical background of the outer tangent circle method. The detection strategy of the single-stage method, including the dispersion curves, FE simulation model, and influence of the size of the detection cell, are presented in Section 3. The multistage stage damage detection in Section 4 includes two-and three-stage detection strategies. Section 5 covers the application of the method for multiple damage cases. The verification of the method through experiments is presented in Section 6. At the end, the article is summed up with a discussion of the results in Section 7, and conclusive remarks in Section 8.

Outer Tangent Circle Method
Consider a plate with the central PZT transducer (P0) acting as an actuator, as shown in Figure 1, to induce waves that propagate radially outward towards the free boundaries of the plate. Imagine circular damage located at coordinates (x, y) from the excitation point. Four other PZT transducers (P1-P4 in Figure 1) act as sensors and are arranged in a square configuration beyond the damage location. The four corner transducers receive the direct signal from P0 and also the signals scattered by the damage. For example, the direct signal from P0 propagates a distance d 0 and arrives at P4 with a time of flight t 040 . However, there is another signal received at P4 by the scattering from the damage that travels a longer distance d 1 + d 2 and requires a longer arrival time t 041 + t 042 . The difference in the time of flight between the direct and damage-scattered signals is the delay time ∆t 04 , as shown in Figure 1, which can be calculated from where v is the group velocity. The distances d 0 , d 1 , and d 2 can be calculated using the geometric relations in Figure 1 by where d is the distance between the corner transducers. Substitution into Equation (1) results in the relation For a known value of ∆t 04 , Equation (3) represents an ellipse with focal points at P0 and P4, as shown with the blue dashed line in Figure 1. Similarly, if the delay times for P1, P2, and P3 are ∆t 01 , ∆t 02 , and ∆t 03 , respectively, then there will be three more ellipses with focal points at P0-P1, P0-P2, and P0-P3, respectively. Therefore, a total of four ellipses can be drawn by x 2 + y 2 + (x 1 − x) 2 + (y 1 − y) 2 = v · ∆t 01 + √ 2 2 d x 2 + y 2 + (x 2 − x) 2 + (y 2 − y) 2  Since the Lamb wave is scattered from the outer periphery of the damage, these ellipses will remain tangent to the damage periphery. A circle externally tangent to these ellipses will exactly cover the damage area. An algorithm based on the least-squares method was designed to draw the four ellipses using Equation (4), and to find the center and size of the externally tangent circle. The center of the externally tangent circle gives the location (x, y) of damage, and the diameter of this circle provides the size.

Dispersion Curves
An infinite number of guided wave modes can be excited in a thin plate with two surfaces. The phase and group velocities of Lamb waves vary as the product of signal frequency and thickness of the plate varies. The dispersion curves for Lamb waves in a thin aluminum plate were drawn using the reverberation-ray matrix method [38], which shows both symmetric (S) and antisymmetric (A) modes propagating simultaneously, as shown in Figure 2. Since the Lamb wave is scattered from the outer periphery of the damage, these ellipses will remain tangent to the damage periphery. A circle externally tangent to these ellipses will exactly cover the damage area. An algorithm based on the least-squares method was designed to draw the four ellipses using Equation (4), and to find the center and size of the externally tangent circle. The center of the externally tangent circle gives the location (x, y) of damage, and the diameter of this circle provides the size.

Dispersion Curves
An infinite number of guided wave modes can be excited in a thin plate with two surfaces. The phase and group velocities of Lamb waves vary as the product of signal frequency and thickness of the plate varies. The dispersion curves for Lamb waves in a thin aluminum plate were drawn using the reverberation-ray matrix method [38], which shows both symmetric (S) and antisymmetric (A) modes propagating simultaneously, as shown in Figure 2. In the low-frequency range, only two modes of Lamb waves (S0, A0) propagate, with the group velocity of the S0 mode being significantly higher than the A0 mode. If one uses a modified signal with a narrow frequency band, the S0 and A0 modes can be identified separately with the peak of the S0 mode arriving first and the A0 mode arriving second. In the mode conversion, after reflection from the damage, S0 will reflect a faster symmetric mode S0 and a slower antisymmetric mode A0 [39]. In the majority of studies, the S0 mode is selected for damage detection because it has lower attenuation, fast propagation velocity, and lower dispersion in comparison with the A0 mode. Therefore, the scattered S0 mode is also easy to identify as the one propagating faster [40,41], and it was exclusively considered for the damage detection method proposed here.

FE Simulation Model
The numerical simulations were performed on an 800 × 800 × 1.293 mm square aluminum plate using the commercial FE software ABAQUS explicit. The material properties used were E = 71 GPa, υ = 0.3, and ρ = 2700 kg/m 3 . Each PZT transducer had a diameter of 5.4 mm and could excite or sense the Lamb wave signals in radially outward directions. In order to accurately characterize the scattering of Lamb waves after interacting with the damage, an FE mesh having at least 8-12 nodes per Lamb wavelength was required [42]. The plate was meshed with 3D shell elements S4R with the element size of 0.5 mm in all simulation models in this work, which gave at least 28 elements per wavelength. The time step for dynamic calculation must be less than the ratio of the minimum distance of any two adjoining nodes to the maximum wave velocity [6], which, for this case, was the velocity of the S0 mode. Therefore, a time step of 0.02 μs was used to accurately simulate the Lamb wave propagation for a total time of 150 μs. The simulation model was first validated by calculating the S0 mode group velocity in Equation (4) by using the transducer configuration shown in Figure  3. The central transducer at point P0 provided the excitation, and the transducers at P1 and P2 were used to acquire the response signals. In the low-frequency range, only two modes of Lamb waves (S 0 , A 0 ) propagate, with the group velocity of the S 0 mode being significantly higher than the A 0 mode. If one uses a modified signal with a narrow frequency band, the S 0 and A 0 modes can be identified separately with the peak of the S 0 mode arriving first and the A 0 mode arriving second. In the mode conversion, after reflection from the damage, S 0 will reflect a faster symmetric mode S 0 and a slower antisymmetric mode A 0 [39]. In the majority of studies, the S 0 mode is selected for damage detection because it has lower attenuation, fast propagation velocity, and lower dispersion in comparison with the A 0 mode. Therefore, the scattered S 0 mode is also easy to identify as the one propagating faster [40,41], and it was exclusively considered for the damage detection method proposed here.

FE Simulation Model
The numerical simulations were performed on an 800 × 800 × 1.293 mm square aluminum plate using the commercial FE software ABAQUS explicit. The material properties used were E = 71 GPa, υ = 0.3, and ρ = 2700 kg/m 3 .
Each PZT transducer had a diameter of 5.4 mm and could excite or sense the Lamb wave signals in radially outward directions. In order to accurately characterize the scattering of Lamb waves after interacting with the damage, an FE mesh having at least 8-12 nodes per Lamb wavelength was required [42]. The plate was meshed with 3D shell elements S4R with the element size of 0.5 mm in all simulation models in this work, which gave at least 28 elements per wavelength. The time step for dynamic calculation must be less than the ratio of the minimum distance of any two adjoining nodes to the maximum wave velocity [6], which, for this case, was the velocity of the S 0 mode. Therefore, a time step of 0.02 µs was used to accurately simulate the Lamb wave propagation for a total time of 150 µs. The simulation model was first validated by calculating the S 0 mode group velocity v in Equation (4) by using the transducer configuration shown in Figure 3. The central transducer at point P0 provided the excitation, and the transducers at P1 and P2 were used to acquire the response signals.  The excitation signal was an in-plane force with a 5.5-cycle sinusoidal tone burst modified by a Hanning window at a central frequency of 383 kHz (as shown in Figure 4), which was generated at the center of transducer P0 to excite the Lamb waves. For the excitation frequency of 383 kHz and the plate thickness of 1.293 mm, the frequency-thickness product . was 514 Hz-m, shown as a red dashed line in Figure 2. This point lay in the low-frequency range of the dispersion curve, and both the S0 and A0 wave modes existed. The S0 mode was clearly the one traveling with higher group velocity, and it could be easily separated from the A0 mode at this frequency. To extract the proper wave information from the complicated wave signals, a CWT, Gabor wavelet [43], was chosen for signal processing. The signals received at P1 and P2 are plotted in Figure 5a after converting them to the nondimensional form by normalizing with the maximum amplitude, and their CWTs are presented in Figure 5b.  The first arrival peaks were the direct excitation signal S0 at P1 and P2, received at the time = 25.69 μs and = 71.96 μs , respectively, and the remaining two peaks were the boundary reflections. Therefore, the traveling time of Lamb waves from P1 to P2 was the difference between and , namely, 46.26 μs. Because the distance between these points was 250 mm, the velocity was calculated as 5.404 × 10 3 m s ⁄ , which was almost the same as that obtained from the dispersion curve in Figure 2 (red dash line; i.e., 5.391 × 10 3 m s ⁄ ). This means that the numerical model had good accuracy in calculating Lamb wave propagation, and CWT was suitable for accurately calculating the delay time which was essential for elliptical reconstruction.  The excitation signal was an in-plane force with a 5.5-cycle sinusoidal tone burst modified by a Hanning window at a central frequency of 383 kHz (as shown in Figure 4), which was generated at the center of transducer P0 to excite the Lamb waves. For the excitation frequency of 383 kHz and the plate thickness of 1.293 mm, the frequency-thickness product f .d was 514 Hz-m, shown as a red dashed line in Figure 2. This point lay in the low-frequency range of the dispersion curve, and both the S 0 and A 0 wave modes existed. The S 0 mode was clearly the one traveling with higher group velocity, and it could be easily separated from the A 0 mode at this frequency. To extract the proper wave information from the complicated wave signals, a CWT, Gabor wavelet [43], was chosen for signal processing. The signals received at P1 and P2 are plotted in Figure 5a after converting them to the nondimensional form by normalizing with the maximum amplitude, and their CWTs are presented in Figure 5b.  The excitation signal was an in-plane force with a 5.5-cycle sinusoidal tone burst modified by a Hanning window at a central frequency of 383 kHz (as shown in Figure 4), which was generated at the center of transducer P0 to excite the Lamb waves. For the excitation frequency of 383 kHz and the plate thickness of 1.293 mm, the frequency-thickness product . was 514 Hz-m, shown as a red dashed line in Figure 2. This point lay in the low-frequency range of the dispersion curve, and both the S0 and A0 wave modes existed. The S0 mode was clearly the one traveling with higher group velocity, and it could be easily separated from the A0 mode at this frequency. To extract the proper wave information from the complicated wave signals, a CWT, Gabor wavelet [43], was chosen for signal processing. The signals received at P1 and P2 are plotted in Figure 5a after converting them to the nondimensional form by normalizing with the maximum amplitude, and their CWTs are presented in Figure 5b.  The first arrival peaks were the direct excitation signal S0 at P1 and P2, received at the time = 25.69 μs and = 71.96 μs , respectively, and the remaining two peaks were the boundary reflections. Therefore, the traveling time of Lamb waves from P1 to P2 was the difference between and , namely, 46.26 μs. Because the distance between these points was 250 mm, the velocity was calculated as 5.404 × 10 3 m s ⁄ , which was almost the same as that obtained from the dispersion curve in Figure 2 (red dash line; i.e., 5.391 × 10 3 m s ⁄ ). This means that the numerical model had good accuracy in calculating Lamb wave propagation, and CWT was suitable for accurately calculating the delay time which was essential for elliptical reconstruction.  The first arrival peaks were the direct excitation signal S 0 at P1 and P2, received at the time t 1 = 25.69 µs and t 2 = 71.96 µs, respectively, and the remaining two peaks were the boundary reflections. Therefore, the traveling time of Lamb waves from P1 to P2 was the difference between t 1 and t 2 , namely, 46.26 µs. Because the distance between these points was 250 mm, the velocity was calculated as 5.404 × 10 3 m/s, which was almost the same as that obtained from the dispersion curve in Figure 2 (red dash line; i.e., 5.391 × 10 3 m/s). This means that the numerical model had good accuracy in calculating Lamb wave propagation, and CWT was suitable for accurately calculating the delay time which was essential for elliptical reconstruction.

Detection Strategy
The transducers can be arranged as a network of measuring points on a large structure, such as the wing of an aircraft or a wind turbine blade. In fact, many modern structures need to be monitored continuously during their service life, and transducers have already been embedded for in situ structural health detection [44,45]. The current method proposes a simple network of transducer detection cells where each square cell is comprised of five transducers: four at the corners of the square and one at the center, as shown in Figure 6a. At the first stage of damage detection, only the central transducer of each cell is excited to locate the damaged cell. Therefore, we only considered the single damaged cell with four PZT transducers located at points P1-P4, as shown in Figure 6b.

Detection Strategy
The transducers can be arranged as a network of measuring points on a large structure, such as the wing of an aircraft or a wind turbine blade. In fact, many modern structures need to be monitored continuously during their service life, and transducers have already been embedded for in situ structural health detection [44,45]. The current method proposes a simple network of transducer detection cells where each square cell is comprised of five transducers: four at the corners of the square and one at the center, as shown in Figure 6a. At the first stage of damage detection, only the central transducer of each cell is excited to locate the damaged cell. Therefore, we only considered the single damaged cell with four PZT transducers located at points P1-P4, as shown in Figure 6b.

Detection Strategy
The transducers can be arranged as a network of measuring points on a large structure, such as the wing of an aircraft or a wind turbine blade. In fact, many modern structures need to be monitored continuously during their service life, and transducers have already been embedded for in situ structural health detection [44,45]. The current method proposes a simple network of transducer detection cells where each square cell is comprised of five transducers: four at the corners of the square and one at the center, as shown in Figure 6a. At the first stage of damage detection, only the central transducer of each cell is excited to locate the damaged cell. Therefore, we only considered the single damaged cell with four PZT transducers located at points P1-P4, as shown in Figure 6b. The minimum size of the estimated damage should be equal to or greater than half of the wavelength of the Lamb wave. For 383-kHz frequency excitation and 5.404 × 10 3 m/s group velocity, the wavelength was 14 mm. So, 7-mm-diameter circular damage, as a through-hole in the plate, was located at (90,−30) mm, and the distance d between the corner transducers was 240 mm, as shown in Figure 6b. The central transducer at P0 excited the Lamb wave signals, and the four corner transducers recorded the signals in the damaged plate. The signal received at P2 in the damaged plate is presented in Figure 7a, which included two fluctuations. The first one was the direct excitation signal S 0 , while second one was the signal scattered by the damage. The result after CWT is shown in Figure 7b, and the difference between the first and second peaks was the delay time ∆t 02 = 32.76 µs. The minimum size of the estimated damage should be equal to or greater than half of the wavelength of the Lamb wave. For 383-kHz frequency excitation and 5.404 × 10 3 m s ⁄ group velocity, the wavelength was 14 mm. So, 7-mm-diameter circular damage, as a through-hole in the plate, was located at (90,−30) mm, and the distance d between the corner transducers was 240 mm, as shown in Figure 6b. The central transducer at P0 excited the Lamb wave signals, and the four corner transducers recorded the signals in the damaged plate. The signal received at P2 in the damaged plate is presented in Figure 7a, which included two fluctuations. The first one was the direct excitation signal S0, while second one was the signal scattered by the damage. The result after CWT is shown in Figure 7b, and the difference between the first and second peaks was the delay time ∆ = 32.76 μs. In order to verify the delay times, the signal received in the damaged plate was compared with the baseline signal from the undamaged plate. The undamaged signal, shown in Figure 8a, was received at P2 in the intact plate, while the damaged signal at the same point is shown in Figure 8b. The difference between undamaged and damaged signals (shown in Figure 8c) clearly indicated the presence of damage in the plate. The CWTs of undamaged and difference signals are shown in Figure  8d, and the delay time was calculated using the first peaks of both signals. The delay time ∆ was 32.76 μs, which was same as that from Figure 7b. It indicated that the strategy adopted in Figure 7b to calculate the delay time without the baseline signal was accurate. Similarly, the delay time ∆ ( = 1, … ,4) could be calculated and is listed in Table 1.  In order to verify the delay times, the signal received in the damaged plate was compared with the baseline signal from the undamaged plate. The undamaged signal, shown in Figure 8a, was received at P2 in the intact plate, while the damaged signal at the same point is shown in Figure 8b. The difference between undamaged and damaged signals (shown in Figure 8c) clearly indicated the presence of damage in the plate. The CWTs of undamaged and difference signals are shown in Figure 8d, and the delay time was calculated using the first peaks of both signals. The delay time ∆t 02 was 32.76 µs, which was same as that from Figure 7b. It indicated that the strategy adopted in Figure 7b to calculate the delay time without the baseline signal was accurate. Similarly, the delay time ∆t 0i (i = 1, . . . , 4) could be calculated and is listed in Table 1.  The four ellipses described by Equation (4) are drawn in Figure 9a using the delay times given in Table 1, and the known values of group velocity and transducer spacing. The damage was assumed to be a circle tangent to these four ellipses. The center of this circle is the damage location, and the diameter of the circle represents its size. The blue circle in Figure 9b represents the damage estimation using this method, while the red circle is the actual damage. The damage location estimated through the single-stage damage detection method was (87,−29) mm, while the size was 2.18 mm. The relative error in damage location and size were 3.33% and 68.92%, respectively. It was proposed to excite the central transducer of each cell separately to find the damaged cell in the first stage of the damage detection. In the subsequent stages, the damage was evaluated locally inside the damaged cell. The four ellipses described by Equation (4) are drawn in Figure 9a using the delay times given in Table 1, and the known values of group velocity and transducer spacing. The damage was assumed to be a circle tangent to these four ellipses. The center of this circle is the damage location, and the diameter of the circle represents its size. The blue circle in Figure 9b represents the damage estimation using this method, while the red circle is the actual damage. The damage location estimated through the single-stage damage detection method was (87,−29) mm, while the size was 2.18 mm. The relative error in damage location and size were 3.33% and 68.92%, respectively. It was proposed to excite the central transducer of each cell separately to find the damaged cell in the first stage of the damage detection. In the subsequent stages, the damage was evaluated locally inside the damaged cell. The four ellipses described by Equation (4) are drawn in Figure 9a using the delay times given in Table 1, and the known values of group velocity and transducer spacing. The damage was assumed to be a circle tangent to these four ellipses. The center of this circle is the damage location, and the diameter of the circle represents its size. The blue circle in Figure 9b represents the damage estimation using this method, while the red circle is the actual damage. The damage location estimated through the single-stage damage detection method was (87,−29) mm, while the size was 2.18 mm. The relative error in damage location and size were 3.33% and 68.92%, respectively. It was proposed to excite the central transducer of each cell separately to find the damaged cell in the first stage of the damage detection. In the subsequent stages, the damage was evaluated locally inside the damaged cell.

Influence of the Detection Cell Size
The effect of the detection cell size on the measuring accuracy was studied by estimating the same damage, 7-mm diameter at (90,−30) mm, for different transducer spacings d. The spacings considered were d = 220, 240, 260, 300 340, 380, and 420 mm. The relative error in evaluating the damage is plotted in Figure 10. The results indicated that there was a significant improvement in the estimation of damage size as the size of the detection cell increased beyond 300 mm. The reason is that the longer distance reduced the measurement error of delay time over the spacing ∆t/d, as both sides of the nondimensional Equation (4) are divided by the distance d to calculate the damage location (x, y). The other reason may be that the S 0 mode was better separated from the A 0 mode beyond d = 300 mm and reduced the inaccuracy in calculating the delay time. There was also slight improvement in the accuracy of location estimation. A recent study has also presented a similar trend of a decrease in error as the distance from the actuator increases [46]. Therefore, the damage location and size can be better estimated using a bigger cell if the attenuation of the Lamb wave is lower.

Influence of the Detection Cell Size
The effect of the detection cell size on the measuring accuracy was studied by estimating the same damage, 7-mm diameter at (90,−30) mm, for different transducer spacings d. The spacings considered were d = 220, 240, 260, 300 340, 380, and 420 mm. The relative error in evaluating the damage is plotted in Figure 10. The results indicated that there was a significant improvement in the estimation of damage size as the size of the detection cell increased beyond 300 mm. The reason is that the longer distance reduced the measurement error of delay time over the spacing ∆t/d, as both sides of the nondimensional Equation (4) are divided by the distance d to calculate the damage location (x, y). The other reason may be that the S0 mode was better separated from the A0 mode beyond d = 300 mm and reduced the inaccuracy in calculating the delay time. There was also slight improvement in the accuracy of location estimation. A recent study has also presented a similar trend of a decrease in error as the distance from the actuator increases [46]. Therefore, the damage location and size can be better estimated using a bigger cell if the attenuation of the Lamb wave is lower.

Multistage Damage Detection Method
The central transducer at P0 was excited during the single-stage damage detection, and the four corner transducers received the signals. However, each PZT transducer could, in fact, both excite and receive the Lamb wave signals. If each corner PZT transducer was also excited, and the other corner transducers were used to detect the signal, then it provided six additional delay times, listed in Table  2, which provided six additional ellipses around the damage, as shown in Figure 11a. To analyze the results, the detection cell can be divided into four triangular subcells formed by the transducer locations. These subcells are represented by Rijk in Figure 11b, where the subscripts "ijk" are the identifiers of the three transducers forming the triangular subcell. Even though the relative error of the damage size estimation was unacceptably large in the single-stage detection, the method could accurately locate the damage location and, thus, the damaged subcell. For the simulated case, the results obtained in single-stage damage detection indicate that the damage was contained in subcell R041, as shown in Figure 11b

Multistage Damage Detection Method
The central transducer at P0 was excited during the single-stage damage detection, and the four corner transducers received the signals. However, each PZT transducer could, in fact, both excite and receive the Lamb wave signals. If each corner PZT transducer was also excited, and the other corner transducers were used to detect the signal, then it provided six additional delay times, listed in Table 2, which provided six additional ellipses around the damage, as shown in Figure 11a. To analyze the results, the detection cell can be divided into four triangular subcells formed by the transducer locations. These subcells are represented by R ijk in Figure 11b, where the subscripts "ijk" are the identifiers of the three transducers forming the triangular subcell. Even though the relative error of the damage size estimation was unacceptably large in the single-stage detection, the method could accurately locate the damage location and, thus, the damaged subcell. For the simulated case, the results obtained in single-stage damage detection indicate that the damage was contained in subcell R 041 , as shown in Figure 11b.

The Two-Stage Damage Detection Method
The damaged cell and subcell R041 were located in the first stage of detection. In the second stage, the damage was estimated through three ellipses formed in the damaged subcell to improve the size estimation. These ellipses were formed using the delay time values Δ , Δ , and Δ between the points P0-P1, P0-P4, and P1-P4, respectively. The first two ellipses were provided by the single-stage detection method, while the additional third ellipse was formed between the corner points P1 and P4 using the delay time Δ . The delay time Δ was calculated by exciting any of the corner transducers (P1 or P4) and sensing the signal at the other corner one. The ellipses used for the secondstage are shown in Figure 12a. The location and size of the damage in the second-stage were calculated as (92,−28) mm and 5.41 mm, respectively. The estimated (blue) and actual (red) damage is shown in Figure 12b, with relative errors in damage location and size of, respectively, 2.98% and 22.67%. It is evident that in the second-stage, both the damage location and size were estimated more accurately than in the single-stage detection method.

The Three-Stage Damage Detection Method
The single-stage damage detection method located the damaged cell and subcell, while the damage location and size estimation were improved in the second-stage by one additional corner ellipse. In order to further improve the damage size estimation in the third-stage of damage detection, besides the additional ellipse in the second stage, two more ellipses were constructed between the

The Two-Stage Damage Detection Method
The damaged cell and subcell R 041 were located in the first stage of detection. In the second stage, the damage was estimated through three ellipses formed in the damaged subcell to improve the size estimation. These ellipses were formed using the delay time values ∆t 01 , ∆t 04 , and ∆t 14 between the points P0-P1, P0-P4, and P1-P4, respectively. The first two ellipses were provided by the single-stage detection method, while the additional third ellipse was formed between the corner points P1 and P4 using the delay time ∆t 14 . The delay time ∆t 14 was calculated by exciting any of the corner transducers (P1 or P4) and sensing the signal at the other corner one. The ellipses used for the second-stage are shown in Figure 12a. The location and size of the damage in the second-stage were calculated as (92,−28) mm and 5.41 mm, respectively. The estimated (blue) and actual (red) damage is shown in Figure 12b, with relative errors in damage location and size of, respectively, 2.98% and 22.67%. It is evident that in the second-stage, both the damage location and size were estimated more accurately than in the single-stage detection method.

The Two-Stage Damage Detection Method
The damaged cell and subcell R041 were located in the first stage of detection. In the second stage, the damage was estimated through three ellipses formed in the damaged subcell to improve the size estimation. These ellipses were formed using the delay time values Δ , Δ , and Δ between the points P0-P1, P0-P4, and P1-P4, respectively. The first two ellipses were provided by the single-stage detection method, while the additional third ellipse was formed between the corner points P1 and P4 using the delay time Δ . The delay time Δ was calculated by exciting any of the corner transducers (P1 or P4) and sensing the signal at the other corner one. The ellipses used for the secondstage are shown in Figure 12a. The location and size of the damage in the second-stage were calculated as (92,−28) mm and 5.41 mm, respectively. The estimated (blue) and actual (red) damage is shown in Figure 12b, with relative errors in damage location and size of, respectively, 2.98% and 22.67%. It is evident that in the second-stage, both the damage location and size were estimated more accurately than in the single-stage detection method.

The Three-Stage Damage Detection Method
The single-stage damage detection method located the damaged cell and subcell, while the damage location and size estimation were improved in the second-stage by one additional corner ellipse. In order to further improve the damage size estimation in the third-stage of damage detection, besides the additional ellipse in the second stage, two more ellipses were constructed between the

The Three-Stage Damage Detection Method
The single-stage damage detection method located the damaged cell and subcell, while the damage location and size estimation were improved in the second-stage by one additional corner ellipse. In order to further improve the damage size estimation in the third-stage of damage detection, besides the additional ellipse in the second stage, two more ellipses were constructed between the corresponding corner transducers of subcell R 041 . These ellipses were formed using the delay times ∆t 13 and ∆t 24 between the transducers, respectively, at points P1-P3 and P2-P4. These three ellipses between the corner transducers of subcell R 041 are shown in Figure 13a and were used to calculate the damage location and size of, respectively, (90,−30) mm and 7.19 mm. The actual and estimated damage in the third-stage damage detection are shown in Figure 13b, indicating almost vanishing errors, respectively, of 0% and 2.66% for damage location and size. corresponding corner transducers of subcell R041. These ellipses were formed using the delay times Δ and Δ between the transducers, respectively, at points P1-P3 and P2-P4. These three ellipses between the corner transducers of subcell R041 are shown in Figure 13a and were used to calculate the damage location and size of, respectively, (90,−30) mm and 7.19 mm. The actual and estimated damage in the third-stage damage detection are shown in Figure 13b, indicating almost vanishing errors, respectively, of 0% and 2.66% for damage location and size.

Application for Multiple Damage Detection
In the case of multiple damage in a plate-like structure, the arrangement of transducers in the form of a network of cells, as proposed in Figure 6a, can help to estimate each damage in its cell. For the case when multiple damages exist in one cell, the further division of a cell into four triangular subcells will probably separate them, as the size of the subcell is very small compared with the whole structure. However, in an extreme case, there can be multiple damages so close to each other that they may exist in one subcell. For such a case, there will be as many damage-scattered peaks in Figure 7 as the number of damages in the plate, and each peak will provide the delay time caused by the corresponding damage. Therefore, based on the delay times calculated from the damage-scattered wave signals, all the ellipses would converge around their corresponding damage. An example multiple damage existing in one subcell is presented in Figure 14, where two damages, D1 and D2, of 10 and 8 mm

Application for Multiple Damage Detection
In the case of multiple damage in a plate-like structure, the arrangement of transducers in the form of a network of cells, as proposed in Figure 6a, can help to estimate each damage in its cell. For the case when multiple damages exist in one cell, the further division of a cell into four triangular subcells will probably separate them, as the size of the subcell is very small compared with the whole structure. However, in an extreme case, there can be multiple damages so close to each other that they may exist in one subcell.
For such a case, there will be as many damage-scattered peaks in Figure 7 as the number of damages in the plate, and each peak will provide the delay time caused by the corresponding damage. Therefore, based on the delay times calculated from the damage-scattered wave signals, all the ellipses would converge around their corresponding damage. An example multiple damage existing in one subcell is presented in Figure 14, where two damages, D 1 and D 2 , of 10 and 8 mm diameters, existed at (20,−40) and (−60,−100) mm, respectively. The signal received at P2 when the central transducer at P0 was excited in the first stage of detection is presented in Figure 15a, which shows two damage-scattered peaks. The CWT of this signal is presented in Figure 15b, which helped to calculate the delay time for each damage using the strategy discussed in Section 3.2.
Similarly, there were two delay time values for each actuator-receiver pair, which formed a total of eight ellipses in the first stage using four actuator-receiver pairs. Each ellipse from each pair converged about one damage, which made four ellipses converge about D 1 and the other four about D 2 . It is not important to know that which delay time or ellipse belongs to which damage; the designed algorithm estimates each damage using the four ellipses converging at one location. Figure 16a shows the multiple damage detection in the single stage, and the comparison of estimated and actual damage is shown in Figure 16b for D 1 and Figure 16c for D 2 . Figure 16d shows the ellipses formed in second stage, and the comparison of estimated and actual damage for D 1 and D 2 is given in Figure 16e,f, respectively. Similarly, Figure 16g shows the ellipses formed in the third stage, and the comparison of estimated and actual damage for D 1 and D 2 is in Figure 16h,i, respectively.
In the case of multiple damage in a plate-like structure, the arrangement of transducers in the form of a network of cells, as proposed in Figure 6a, can help to estimate each damage in its cell. For the case when multiple damages exist in one cell, the further division of a cell into four triangular subcells will probably separate them, as the size of the subcell is very small compared with the whole structure. However, in an extreme case, there can be multiple damages so close to each other that they may exist in one subcell. For such a case, there will be as many damage-scattered peaks in Figure 7 as the number of damages in the plate, and each peak will provide the delay time caused by the corresponding damage. Therefore, based on the delay times calculated from the damage-scattered wave signals, all the ellipses would converge around their corresponding damage. An example multiple damage existing in one subcell is presented in Figure 14, where two damages, D1 and D2, of 10 and 8 mm  Figure 15a, which shows two damage-scattered peaks. The CWT of this signal is presented in Figure 15b, which helped to calculate the delay time for each damage using the strategy discussed in Section 3.2. Similarly, there were two delay time values for each actuator-receiver pair, which formed a total of eight ellipses in the first stage using four actuator-receiver pairs. Each ellipse from each pair converged about one damage, which made four ellipses converge about D1 and the other four about D2. It is not important to know that which delay time or ellipse belongs to which damage; the designed algorithm estimates each damage using the four ellipses converging at one location. Figure  16a shows the multiple damage detection in the single stage, and the comparison of estimated and actual damage is shown in Figure 16b for D1 and Figure 16c for D2. Figure 16d shows the ellipses formed in second stage, and the comparison of estimated and actual damage for D1 and D2 is given in Figure 16e,f, respectively. Similarly, Figure 16g shows the ellipses formed in the third stage, and the comparison of estimated and actual damage for D1 and D2 is in Figure 16h and (i), respectively. diameters, existed at (20,−40) and (−60,−100) mm, respectively. The signal received at P2 when the central transducer at P0 was excited in the first stage of detection is presented in Figure 15a, which shows two damage-scattered peaks. The CWT of this signal is presented in Figure 15b, which helped to calculate the delay time for each damage using the strategy discussed in Section 3.2. Similarly, there were two delay time values for each actuator-receiver pair, which formed a total of eight ellipses in the first stage using four actuator-receiver pairs. Each ellipse from each pair converged about one damage, which made four ellipses converge about D1 and the other four about D2. It is not important to know that which delay time or ellipse belongs to which damage; the designed algorithm estimates each damage using the four ellipses converging at one location. Figure  16a shows the multiple damage detection in the single stage, and the comparison of estimated and actual damage is shown in Figure 16b for D1 and Figure 16c for D2. Figure 16d shows the ellipses formed in second stage, and the comparison of estimated and actual damage for D1 and D2 is given in Figure 16e,f, respectively. Similarly, Figure 16g shows the ellipses formed in the third stage, and the comparison of estimated and actual damage for D1 and D2 is in Figure 16h and (i), respectively. The relative error in the location and size estimation for each damage is presented in Table 3. The results indicate that there was improvement in damage estimation in the second and third stages of the detection for multiple damages in one subcell. However, the relative errors in the case of multiple damages were slightly higher than the single damage case.

Experimental Evaluation
The accuracy of the proposed multistage damage detection method was further verified through experimental investigation. The plate dimensions, material properties, and transducer configuration were exactly the same as those in the numerical simulation shown in Figure 6b. The damage was a punched hole of 7-mm diameter at the location of (90,−30) mm. The PZT wafers used in the experiments were made of PSN-33 with a density of 7.70 × 10 3 kg/m 3 and a diameter of 5.4 mm. The PZT wafers had a resonance frequency of 383 kHz and had the ability to excite/sense the Lamb wave signals in radially outward directions. The experiments were performed for two PZT configurations with d = 240 mm and d = 300 mm, respectively. The schematic diagram of the experimental setup for multistage damage detection is shown in Figure 17a, and the real-time experimental setup is shown in Figure 17b. A signal generator (RIGOL, DG1022) was programmed to generate the signal shown in Figure 4a at the excitation transducer P0, and the signal was amplified to 45 V using a power amplifier (KROHN-HITE, 7602M). A four-channel oscilloscope (LeCroy, LC574AL) was used to The relative error in the location and size estimation for each damage is presented in Table 3.
The results indicate that there was improvement in damage estimation in the second and third stages of the detection for multiple damages in one subcell. However, the relative errors in the case of multiple damages were slightly higher than the single damage case.

Experimental Evaluation
The accuracy of the proposed multistage damage detection method was further verified through experimental investigation. The plate dimensions, material properties, and transducer configuration were exactly the same as those in the numerical simulation shown in Figure 6b. The damage was a punched hole of 7-mm diameter at the location of (90,−30) mm. The PZT wafers used in the experiments were made of PSN-33 with a density of 7.70 × 10 3 kg/m 3 and a diameter of 5.4 mm. The PZT wafers had a resonance frequency of 383 kHz and had the ability to excite/sense the Lamb wave signals in radially outward directions. The experiments were performed for two PZT configurations with d = 240 mm and d = 300 mm, respectively. The schematic diagram of the experimental setup for multistage damage detection is shown in Figure 17a, and the real-time experimental setup is shown in Figure 17b. A signal generator (RIGOL, DG1022) was programmed to generate the signal shown in Figure 4a at the excitation transducer P0, and the signal was amplified to 45 V using a power amplifier (KROHN-HITE, 7602M). A four-channel oscilloscope (LeCroy, LC574AL) was used to record the response of each transducer. The oscilloscope captured the response signals over a sampling time of 200 µs at a sampling rate of 1 GHz. Every signal was acquired for an average of 500 times to reduce the noise.
sampling time of 200 μs at a sampling rate of 1 GHz. Every signal was acquired for an average of 500 times to reduce the noise.
During the single-stage of damage detection, the central transducer at point P0 was excited to induce the Lamb wave, and transducers P1-P4 (d = 240 mm) and P5-P8 (d = 300 mm) were used to measure the received signals. In multistage damage detection, the corner PZT transducers at P1, P2, P5, and P6 were also used to excite Lamb waves for the second and third stages of damage detection. The signal received at P2 was processed in dimensionless form and is shown in Figure 18a, which showed fluctuations caused by the scattering from the damage and the plate's edge. The acquired signal shown in Figure 18a was then processed by Sym8 wavelet function in MATLAB to reduce the noise. The denoised signal is presented in Figure 18b, which looks similar to the signal received at the same point during the simulation (as shown in Figure 7b). The CWTs of recorded and denoised signals are presented in Figure 18c,d, respectively, and the delay time was calculated as ∆ = 31.27 μs, which was close to the delay time ∆ = 32.76 μs calculated through FE simulations. During the single-stage of damage detection, the central transducer at point P0 was excited to induce the Lamb wave, and transducers P1-P4 (d = 240 mm) and P5-P8 (d = 300 mm) were used to measure the received signals. In multistage damage detection, the corner PZT transducers at P1, P2, P5, and P6 were also used to excite Lamb waves for the second and third stages of damage detection.
The signal received at P2 was processed in dimensionless form and is shown in Figure 18a, which showed fluctuations caused by the scattering from the damage and the plate's edge. The acquired signal shown in Figure 18a was then processed by Sym8 wavelet function in MATLAB to reduce the noise. The denoised signal is presented in Figure 18b, which looks similar to the signal received at the same point during the simulation (as shown in Figure 7b). The CWTs of recorded and denoised signals are presented in Figure 18c,d, respectively, and the delay time was calculated as ∆t 02 = 31.27 µs, which was close to the delay time ∆t 02 = 32.76 µs calculated through FE simulations. The delay times, similarly calculated through experiments for signals received at all transducers, are shown in Table 4.

Cell Size
Delay Times (µs) The delay times, similarly calculated through experiments for signals received at all transducers, are shown in Table 4. The relative error in damage estimation at different stages of detection through simulations and experiments are compared with the actual location and size of the damage in Table 5. The experimental results of damage location and size for single-stage damage detection in the smaller detection cell had relative errors, respectively, of 15.28% and 128.56%. The relative errors in the same cell were reduced in the second stage to, respectively, 5.13% and 56.64%, and in the third stage to, respectively, 3.16% and 24.20%. A similar trend in the improvement of damage detection was found in the larger detection cell. The damage assessment accuracy of the multistage method was significantly higher than the single-stage detection in both cases; in particular, there was a significant improvement in the damage size detection.

Discussion of Results
During the numerical simulation, the excitation of the central PZT transducer in the single-stage damage detection method caused Lamb waves to propagate radially outward in all directions, and be scattered by any damage in the plate. The single-stage method efficiently located the damage (3.33% of error for d = 240 mm). The relative error in damage size estimation in the single-stage method was unacceptable (68.92% for d = 240 mm). However, the accuracy of damage location can help to identify the damaged cell, and the triangular subcell containing the damage. So, the single-stage damage detection method can directly obtain the size and location information of the damage, and it can locate the damaged subcell for the subsequent detection stages, especially for improving the estimation accuracy of damage size.
Based on the same transducer configuration, the second stage of damage detection could be conducted to improve the damage estimation by two ellipses from the single-stage detection, and one additional corner ellipse around the triangular subcell. The additional ellipse was obtained by exciting a corner transducer, and receiving at the other corner transducer separately in this subcell. With the three ellipses, the second stage of damage detection could reduce the relative errors in damage location and size, respectively, from 3.33% to 2.98% and from 68.92% to 22.67%.
Similarly, using the same transducer configuration, the third stage of damage detection was carried out by using the additional ellipse from the second stage and two extra ellipses, which were obtained by exciting the two corner transducers of the damage subcell and receiving at the other two corner transducers in the diagonally opposite direction. Using these three ellipses could reduce the relative errors from the second to the third stage in damage location and size, respectively, from 2.98% to 0% and 22.67% to 2.66%.
From the trend shown in Figure 10, it can be concluded that the bigger ellipses obtained by larger detection cells could provide a better estimation of damage. Therefore, compared with the fast detection of the single-stage method with only a one-time excitation, the second stage could obtain better results by including one additional large ellipse from one extra excitation, and the third stage had the best damage estimation from two additional large ellipses by two extra excitations. This is because the single-stage method obtained the ellipses with transducer spacing √ 2 2 d, the additional ellipse in the second stage was obtained by the corner transducers with large spacing d, and the third stage included two extra ellipses with larger transducer spacing √ 2d. It can be noticed that the bigger ellipse touching the damage, as shown in Figure 13b, was almost a straight line and more accurately defined the tangent to the damage periphery than the smaller ellipse, as shown in Figure 12b. As a result, depending on the location of the damage, the accuracy of the damage size was limited in the single-stage damage detection method by the size of the detection cell. In addition, the longer distance reduced the error in calculating the delay time by Equation (4). However, this was restricted by the attenuation of Lamb wave propagation.
The same conclusion can be confirmed by the experimental results in Table 5; that is, the bigger ellipses formed in the large detection cell (d = 300 mm) could obtain better results than the small ones (d = 240 mm) at the same stage, and the detection accuracy could be improved at the higher stage in each detection cell. Even if the three-stage detection could provide accurate damage information, the excitation of the central transducer in single-stage damage detection was still important to locate the damaged cell and subcell used in the second and third damage detection stages. Actually, the utilization of third-stage detection depends on the nature of the problem and the accuracy requirements in practice.
Comparing Figure 18c,d, the difference is hard to identify. This means the CWT had a better signal-to-noise ratio, and could extract the time-frequency feature exactly from the noisy signal. This is also the reason why CWT could process the Lamb wave signals of the damaged plate directly without a baseline. Generally, a greater value of error is expected in the larger cell due to the dispersion of Lamb waves, which leads to the inaccuracy in determining the delay time. However, the accuracy of this damage detection method improves as the size of the detection cell increases. This is because the CWT is a better time-frequency representative and can depress the dispersive characteristics effectively.

Conclusions
The Lamb-wave-based multistage damage detection method introduced in this research can accurately quantify both the location and size of damage. According to the PZT transducer network in practice, each square transducer cell can be used to identify the damage locally. The single stage can be performed at each cell simultaneously in a large PZT network to estimate the damage location and size and, further, to identify the damaged cell and subcell. In order to further improve the damage size detection, the second stage can be carried out only in the damaged subcell. If further improvement in size estimation is required, the third stage can be utilized for higher accuracy achievement. Compared with existing ellipse-based methods, the amount of data collected in this method is very small due to a small number of transducers used to form the damaged cell. Four excitations are required from the single to the third stage: one excitation each in the single and second stages, and two excitations in the third stage of detection. The data collected in one stage is also reused in subsequent stages without any other additional tools or techniques. Therefore, the proposed method can save both computational time and evaluation cost during the NDT/SHM process. Further, the proposed method does not require a baseline signal or any additional transducers in higher stages of detection, which makes it suitable for on-line monitoring. The transducer configuration, signal-processing, and outer tangent circle method introduced in the single-stage damage detection method remain exactly the same in the second and third stages of damage detection. The method is relatively simple and straightforward to implement for detection of circle-like damage in thin plate structures used in civil, automobile, aerospace, ship building, and wind turbine industries.

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