Assessment of the Length and Depth of Delamination-Type Defects Using Ultrasonic Guided Waves

: Fiber-reinforced composite laminates are being increasingly used in various engineering components in the sectors of aerospace and green energy. Due to impacts throughout the service life of the structure, matrix breakage and delaminations signiﬁcantly altering the structural integrity of the laminate can occur. Hence, robust guided wave structural health monitoring systems are required to ensure continuous safety of engineering structures. In this paper, the ultrasonic method based on the analysis of A 0 mode reﬂecting within the defected area has been proposed to extract the length and the depth of the delamination-type defect. The technique proposed in this study extracts the depth of the damage by analyzing the magnitude variations of direct A 0 mode which are caused by the di ﬀ erence of wave velocities in the upper and lower sub-laminates. This results in an altering and frequency-dependent forward-scattered amplitude of direct A 0 mode. Furthermore, the proposed approach uses previously obtained information about the depth of the defect, which allows for the determination of the phase velocities of A 0 and S 0 modes in the upper and lower sub-laminates. As a result, the accuracy of the damage length estimation is increased. The performance of the proposed method was proven with 2D and 3D numerical simulations and experiments on samples with artiﬁcial defects. The method validation results showed that the proposed method with some limitations is capable of extracting the length of delamination with an approximate error below 6%.


Introduction
Delamination is one of the most common defects found in fiber-reinforced composite laminates due to their weak transverse tensile and inter-laminar shear strengths [1]. Such defects may be caused by an impact and tend to develop internally between the neighboring fiber layers of the laminate. This results in loss of the strength of the structure without any visual evidence. As a result, sophisticated non-destructive testing and signal processing techniques are of great significance for the structural integrity assessment of fiber reinforced laminates. Ultrasonic guided wave structural health monitoring (SHM) emerged as one of the most powerful tools which uses a network of embedded sensors to detect and continuously track the development of such defects in various industrial fields [2]. SHM systems usually employ ultrasonic guided waves which propagate in elongated thin-walled structures, possess low attenuation and are sensitive to the integral elastic properties of the material; thus, they can be used to infer large structures with only a few coupled transducers. Numerous researches have targeted non-linear, multimodal behavior, scattering and energy leakage of guided waves for the detection and sizing of delamination-type defects [3][4][5][6][7].
By interaction of guided waves with delamination, waves are scattered and converted to other modes [8]. The mechanism of wave scattering depends on the type of the guided wave mode as well as on the depth and size of the defect. Various investigations have demonstrated that A 0 and S 0 modes of guided waves propagate in the two regions divided by the delamination individually and then interact with each other after exiting the defect [9,10]. The anti-symmetric modes are found to be more suitable for the assessment of delaminations as they show higher sensitivity to delamination located at different through-thickness positions compared to the symmetric mode [11]. A few studies on A 0 mode interaction with asymmetrical and symmetrical delaminations in laminated composites were presented by Ramadas et al. [12,13]. It was found that, in the case of symmetrical delamination, incident A 0 mode and converted A 0 (A 0 ›S 0 ›A 0 ) mode propagate after the exit of delamination. In the case of an asymmetrical defect, an additional S 0 mode (A 0 ›S 0 ›S 0 ) appears upon the interaction with the defect. The authors also proposed the method for the assessment of the delamination size based on the time-of-flight measurement of A 0 mode [14]. However, such a method appeared to be valid at low dispersion regions of A 0 mode as the velocities in the upper and lower sub-laminate layer of the defect were assumed to be equal to the velocity in the entire structure. In case of sub-surface delaminations, the velocities in sub-laminates can be essentially different, which leads to high errors.
Lamb wave scattering and conversion at the leading and trailing edge of delamination is another topic which has been extensively studied by numerous research groups. For example, Schaal et al. [15] estimated frequency-dependent scattering coefficients for fundamental Lamb wave modes separately at the leading and trailing edges of delamination. Hu et al. [16] developed a technique to locate the delamination-type defects by using two embedded transducers and an S 0 mode by measuring signals reflected from the tip of the defect. A study on S 0 mode interaction with the delamination-type defect was presented, which showed that, at each end of delamination, the transmitted and reflected A 0 /S 0 modes are being generated with each mode having different reflection and transmission coefficients depending on the local bending stiffness of the structure. An extensive study on Lamb wave mode conversion in a 2D plate containing a finite length delamination parallel with the surface was presented by Shkerdin and Glorieux [17]. Mode conversions between Lamb wave modes and between Rayleigh wave and Lamb modes were analyzed. It was found that the transmission coefficients of Lamb waves are both defect length-and depth-dependent, and they vary for different transmitted modes. Birt et al. [18] studied the dependencies between the magnitude of the reflected S 0 mode and the width of delamination.
As a result of such complex scattering and mode conversion, various linear acoustic tools for the detection and localization of delamination-type defects can be developed. For example, by measuring the time-of-flight (ToF) of the reflected and transmitted waves, the existence and location of delamination can be estimated [8,19]. Additionally, by analyzing the ToF and velocities of the modes propagating at the layers divided by delamination, the length of the defect can be assessed [20,21]. However, the success of the above-mentioned methods is strongly influenced by the dispersion and superposition of the modes of guided waves. The defects are mostly weak reflectors, while the monitoring system tends to detect and locate them in a large structure by using a small number of sensors. This leads to a low signal-to-noise ratio (SNR), which further complicates the signal analysis and the estimation of the time of flight (ToF).
Most of the researches in the field of linear acoustics are currently targeting the detection, localization and sizing of delamination-type defects. However, the assessment of the depth of a defect is no less important for the evaluation of defect severity. There are only a few works available which discuss the depth estimation of a delamination-type defect. Recently, Eremin et al. [22] analyzed the relationship between the length and depth of delamination with the scattering resonance frequencies of delamination. Munian et al. [23] analyzed the reflection of A 0 mode as a function of the depth-wise position of delamination. It was concluded that mid-plane delamination scatters less energy of the asymmetric mode, meanwhile, off mid-plane delaminations additionally introduce a reflection of converted S 0 mode. Gupta et al. [24] exploited S 0 mode to analyze its reflection from delamination situated at different depths and ply layups. Reflection coefficients as a function of the through-thickness location of delamination were estimated for various ply layups, which showed the importance of the Appl. Sci. 2020, 10, 5236 3 of 19 defect position and the composition of the material to the detectability of delamination. Recently, Murat et al. [25,26] analyzed the scattering of A 0 mode from delaminations of different sizes which were situated at different depths. It was observed that the delamination length exerts strong influence on the scattering directivity, while its depth influences the magnitude of the scattered wave. Despite some promising results, the abovementioned works in general have discussed the fundamental influence of the defect depth on the scattering of guided wave modes. There is still a lack of methods to extract the depth of the defect without any initial knowledge about its through-thickness location when dispersion is involved.
This paper proposes the SHM method to extract the absolute depth and size of asymmetrical delamination-type defects in plate-like structures with parallel surfaces. The technique proposed in this study extracts the depth of the damage by analyzing the amplitude variation patterns of direct A 0 mode at different excitation frequencies. These magnitude variations are caused by the difference of wave velocities in the upper and lower sub-laminates due to their uneven thicknesses. This leads to the varying phase differences of the wave-packets at the trailing edge of the defect and results in an altering forward scattered amplitude of the direct A 0 mode. Hence, in comparison with the predefined reference dependencies, the absolute values of the delamination depth can be extracted. Such magnitude variation is related to the excitation frequency as the phase velocity of A 0 mode is frequency and thickness dependent. As a result, different magnitude variation patterns at different excitation frequencies can be observed. Having such patterns at different frequencies allow to increase the accuracy of damage depth estimation. This is very important in the presence of strong dispersion, as the rapid changes in velocity of guided wave modes can introduce large errors. Furthermore, the information on the damage depth can be exploited for the assessment of the defect length as the actual depth of the defect determines the phase velocities of A 0 and S 0 modes in the upper and lower sub-laminates. Then, the delamination size can be evaluated by measuring the ToF between the wave packets of the converted and direct A 0 modes. Hence, the approach proposed in this research can provide information about the absolute depth and the length of the damage employing a few time traces measured at a series of different excitation frequencies. The technique proposed in this study is essentially different from the currently existing ones as it uses multiple excitation frequencies, which leads to different patterns of A 0 mode magnitude variation and, subsequently, to better resolution. Most of the available research studies use constant a-priori wave velocities for damage size estimation which is not valid in the presence of dispersion. The technique proposed in this paper can be used even in the presence of the dispersion of A 0 mode as it exploits the dependency of the phase velocity on the thickness of the layers above and below the defect. This paper focuses on the proof of the concept of the proposed delamination depth and length estimation method by analyzing a single delamination-type defect. However, in further studies, it is expected to extend the presented method to extract features of multiple parallel delaminations which are a common outcome of low velocity impact damage.

Interaction of A 0 Mode with Delamination-Type Defect
Let us analyze a composite plate with length l and thickness h which is shown in Figure 1. It will be investigated with the 2D approach, which means that the plate is infinite along the y axis. Let us consider that the delamination-type defect is present in the structure, which is l 1 long and situated at distances of l 0 and l 2 from the source (excitation) and monitoring (reception) points, respectively, at a depth of h 1 with respect to the top surface. The delamination is parallel to the surface of the sample (Figure 1).
such frequencies, the wavelengths of A0 mode varies from 17.2 mm (at 50 kHz), to 5.6 mm (at 200 kHz) at 4mm GFRP plate (Ex = 10 GPa, Ez = 35.7 GPa; υxz = 0.325, υzx = 0.091, υyx = 0.35; Gxz = 2.8 GPa; ρ = 1800 kg/m 3 ). This shows that both A0, S0 and higher order A1 and S1 modes can exist at such a constant wavelength (see Figure 2a). However, using a Plexiglas wedge (2700 m/s) as a coupling medium and normal incidence angles, the fundamental A0 mode is being generated most effectively (see Figure  2b), and sophisticated mode isolation techniques were not applied. The signal of A0 mode excited at the source point and measured at the monitoring point after propagation through the delaminated area will be the subject of all the upcoming investigations. Upon the interaction of A0 mode with the delamination-type defect, the reflection, transmission and mode conversion occur at each end of the damage. In such a way, at the leading edge of the damage, the incident A0 mode reflects back, and in the forward direction it splits into wave packets which accordingly propagate on the sub-laminates above and below the defect. Moreover, mode conversion occurs at the leading edge; therefore, part of the energy transforms into S0 mode as well (see Figure  3a). Similarly, at the trailing edge of the damage, both A0 and S0 modes are reflecting back, propagating forward and converting to each other (see Figure 3b). This representation does not take into account the scattering coefficients of each reflected/converted mode. The fundamental A 0 mode has been selected in this paper. As multiple guided wave modes are being generated simultaneously, different strategies exist to excite isolated modes, such as angle beam excitation [27], element width selection [28], inter-element spacing or interdigital transducers [29,30], in-phase and out-of-phase excitation [31], etc. In this research the fundamental dispersion zone has been selected and the excitation frequencies were limited from 50 kHz up to 200 kHz. At such frequencies, the wavelengths of A 0 mode varies from 17.2 mm (at 50 kHz), to 5.6 mm (at 200 kHz) at 4 mm GFRP plate (E x = 10 GPa, E z = 35.7 GPa; υ xz = 0.325, υ zx = 0.091, υ yx = 0.35; G xz = 2.8 GPa; ρ = 1800 kg/m 3 ). This shows that both A 0 , S 0 and higher order A 1 and S 1 modes can exist at such a constant wavelength (see Figure 2a). However, using a Plexiglas wedge (2700 m/s) as a coupling medium and normal incidence angles, the fundamental A 0 mode is being generated most effectively (see Figure 2b), and sophisticated mode isolation techniques were not applied.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 4 of 19 [29,30], in-phase and out-of-phase excitation [31], etc. In this research the fundamental dispersion zone has been selected and the excitation frequencies were limited from 50 kHz up to 200 kHz. At such frequencies, the wavelengths of A0 mode varies from 17.2 mm (at 50 kHz), to 5.6 mm (at 200 kHz) at 4mm GFRP plate (Ex = 10 GPa, Ez = 35.7 GPa; υxz = 0.325, υzx = 0.091, υyx = 0.35; Gxz = 2.8 GPa; ρ = 1800 kg/m 3 ). This shows that both A0, S0 and higher order A1 and S1 modes can exist at such a constant wavelength (see Figure 2a). However, using a Plexiglas wedge (2700 m/s) as a coupling medium and normal incidence angles, the fundamental A0 mode is being generated most effectively (see Figure  2b), and sophisticated mode isolation techniques were not applied. The signal of A0 mode excited at the source point and measured at the monitoring point after propagation through the delaminated area will be the subject of all the upcoming investigations. Upon the interaction of A0 mode with the delamination-type defect, the reflection, transmission and mode conversion occur at each end of the damage. In such a way, at the leading edge of the damage, the incident A0 mode reflects back, and in the forward direction it splits into wave packets which accordingly propagate on the sub-laminates above and below the defect. Moreover, mode conversion occurs at the leading edge; therefore, part of the energy transforms into S0 mode as well (see Figure  3a). Similarly, at the trailing edge of the damage, both A0 and S0 modes are reflecting back, propagating forward and converting to each other (see Figure 3b). This representation does not take into account the scattering coefficients of each reflected/converted mode. The signal of A 0 mode excited at the source point and measured at the monitoring point after propagation through the delaminated area will be the subject of all the upcoming investigations. Upon the interaction of A 0 mode with the delamination-type defect, the reflection, transmission and mode conversion occur at each end of the damage. In such a way, at the leading edge of the damage, the incident A 0 mode reflects back, and in the forward direction it splits into wave packets which accordingly propagate on the sub-laminates above and below the defect. Moreover, mode conversion occurs at the leading edge; therefore, part of the energy transforms into S 0 mode as well (see Figure 3a). Similarly, at the trailing edge of the damage, both A 0 and S 0 modes are reflecting back, propagating forward and converting to each other (see Figure 3b). This representation does not take into account the scattering coefficients of each reflected/converted mode. Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 19 (a) (b) In this research, only the forward scattered wave-packets of A0 mode will be considered; they are squared red in Figure 3b. Mathematically, the spectrum of the signal of A0 mode UIA 0 (f) which arrives first to the monitoring (reception) point can be expressed as: where UsA 0 (f) is the frequency spectrum of the excitation signal of A0 mode; H(f,cp,l) is the transfer function of the waves propagating in the corresponding layer; cpA 0 (f), cpS 0 (f) are the phase velocities of A0 and S0 modes, respectively; l is the propagation distance (l0, l2-propagation distances in the defectfree area before and after the defect, respectively, l1-length of the defect); L is the separation between the source and the monitoring points (L = l0 + l1 + l2); and α(f) is the frequency dependent attenuation function; kT01A, kT02A, kT10S and kT20S are the transmission coefficients for A0 and S0 modes at the leading and trailing edges of the delamination, respectively, for the layer above and below the defect. The estimation of such reflection and transmission coefficients was out of the scope of this study. They essentially depend on the guided wave mode, the delamination depth, the material properties and the total thickness of the plate, and they should be determined for each particular case separately. The wave packet described by Equation (1) is A0 mode, which is a superposition of the signals traveling at upper and lower sub-laminates. This wave-packet converts from A0 mode to S0 mode at the leading edge and back to A0 at the trailing edge of the defect. Due to this reason, this wave packet UIA 0 (f) arrives earlier compared to the direct transmission of A0 mode which does not convert within the defected area. This happens due to the reason that S0 mode possesses a greater group velocity at low frequencies compared to A0 mode. In this research, such a wave-packet shall be referred to as "first A0 arrival". Consequently, the frequency spectrum of the direct transmission of A0 mode can be written as: The "direct A0 transmission" described by Equation (3) is the wave-packet which does not convert to other modes at the leading and trailing edges of delamination. In this research, only the forward scattered wave-packets of A 0 mode will be considered; they are squared red in Figure 3b. Mathematically, the spectrum of the signal of A 0 mode U IA0 (f ) which arrives first to the monitoring (reception) point can be expressed as: where U sA0 (f ) is the frequency spectrum of the excitation signal of A 0 mode; H(f,c p ,l) is the transfer function of the waves propagating in the corresponding layer; c pA0 (f ), c pS0 (f ) are the phase velocities of A 0 and S 0 modes, respectively; l is the propagation distance (l 0 , l 2 -propagation distances in the defect-free area before and after the defect, respectively, l 1 -length of the defect); L is the separation between the source and the monitoring points (L = l 0 + l 1 + l 2 ); and α(f ) is the frequency dependent attenuation function; k T01A , k T02A , k T10S and k T20S are the transmission coefficients for A 0 and S 0 modes at the leading and trailing edges of the delamination, respectively, for the layer above and below the defect. The estimation of such reflection and transmission coefficients was out of the scope of this study. They essentially depend on the guided wave mode, the delamination depth, the material properties and the total thickness of the plate, and they should be determined for each particular case separately. The wave packet described by Equation (1) is A 0 mode, which is a superposition of the signals traveling at upper and lower sub-laminates. This wave-packet converts from A 0 mode to S 0 mode at the leading edge and back to A 0 at the trailing edge of the defect. Due to this reason, this wave packet U IA0 (f ) arrives earlier compared to the direct transmission of A 0 mode which does not convert within the defected area. This happens due to the reason that S 0 mode possesses a greater group velocity at low frequencies compared to A 0 mode. In this research, such a wave-packet shall be referred to as "first A 0 arrival". Consequently, the frequency spectrum of the direct transmission of A 0 mode can be written as: Appl. Sci. 2020, 10, 5236 6 of 19 The "direct A 0 transmission" described by Equation (3) is the wave-packet which does not convert to other modes at the leading and trailing edges of delamination.
The scattering of the incident A 0 mode at the delamination-type defect is further illustrated with B-scans in Figure 4 which were simulated by using the Finite Element plane-strain model on a 2000 mm × 4 mm 2D GFRP plate. These B-scans represent a centered delamination-type defect located at a depth of 1.5 mm with respect to the top surface of the sample. A 0 mode was excited with a Gaussian envelope tone burst of three cycles and a central frequency of 80 kHz. The signals which are analyzed in this paper are denoted as u IIA0 (t) ("direct A 0 transmission") and u IA0 (t) ("first A 0 arrival") in Figure 4.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 19 The scattering of the incident A0 mode at the delamination-type defect is further illustrated with B-scans in Figure 4 which were simulated by using the Finite Element plane-strain model on a 2000 mm×4 mm 2D GFRP plate. These B-scans represent a centered delamination-type defect located at a depth of 1.5 mm with respect to the top surface of the sample. A0 mode was excited with a Gaussian envelope tone burst of three cycles and a central frequency of 80 kHz. The signals which are analyzed in this paper are denoted as uIIA 0 (t) ("direct A0 transmission") and uIA 0 (t) ("first A0 arrival") in Figure  4.

Method to Estimate the Depth of Defect
As it was shown by Equation (3), the "direct A0 transmission" is a superposition of the waves traveling at the sub-laminates below and above the damage which does not convert to other modes, neither at the leading, nor at the trailing edge of the defect. If the delamination is asymmetrical across the thickness of the laminate (h1 ≠ h2 in Figure 1), the phase velocities at sub-laminates are different (cp 1 ≠ cp 2 ), and they depend on the excitation frequency and thickness of each sub-laminate. Beyond the defect, the modes from the upper and lower sub-laminates interfere with each other, which results in the single forward scattered A0 mode uIIA 0 (t). The magnitude of this mode is proportional to the result of interference which can be either in phase, out of phase or intermediate. If the transducer is driven by a set of different frequencies, a collection of magnitude variation patterns can be obtained and used for damage depth extraction. The use of multiple excitation frequencies enables one to exploit the frequency dependent phase velocity of A0 mode, especially in the case of the presence of dispersion. Thus, the idea of the proposed method to estimate the depth of damage relies on the magnitude measurements of "direct A0 transmission" uIIA 0 (t) at burst signals with different central frequencies (f1, f2, … fn). We should note that this technique is only valid if the delamination is asymmetrical h1 ≠ h/2 ≠ h2 and cp 1 ≠ cp 2 . For ideally symmetrical delamination, there will be no out-ofphase interference of A0 mode as the velocities at the layer above and below the defect will be very close.
To extract the absolute values of the defect depth, a database of reference dependencies is required which would represent the collection of magnitude variation patterns at different excitation frequencies, the depth and the length of delaminations. Then, the estimated reference dependency can be compared to the experimental measurement in order to extract the depth of delamination. It is presumed that the reference dependency which best matches the experiments is the one that gives the closest definition of the damage depth in the structure. The step-by-step procedure of the damage depth estimation can be outlined as follows:

Method to Estimate the Depth of Defect
As it was shown by Equation (3), the "direct A 0 transmission" is a superposition of the waves traveling at the sub-laminates below and above the damage which does not convert to other modes, neither at the leading, nor at the trailing edge of the defect. If the delamination is asymmetrical across the thickness of the laminate (h 1 h 2 in Figure 1), the phase velocities at sub-laminates are different (c p1 c p2 ), and they depend on the excitation frequency and thickness of each sub-laminate. Beyond the defect, the modes from the upper and lower sub-laminates interfere with each other, which results in the single forward scattered A 0 mode u IIA0 (t). The magnitude of this mode is proportional to the result of interference which can be either in phase, out of phase or intermediate. If the transducer is driven by a set of different frequencies, a collection of magnitude variation patterns can be obtained and used for damage depth extraction. The use of multiple excitation frequencies enables one to exploit the frequency dependent phase velocity of A 0 mode, especially in the case of the presence of dispersion. Thus, the idea of the proposed method to estimate the depth of damage relies on the magnitude measurements of "direct A 0 transmission" u IIA0 (t) at burst signals with different central frequencies (f 1 , f 2 , . . . f n ). We should note that this technique is only valid if the delamination is asymmetrical h 1 h/2 h 2 and c p1 c p2 . For ideally symmetrical delamination, there will be no out-of-phase interference of A 0 mode as the velocities at the layer above and below the defect will be very close.
To extract the absolute values of the defect depth, a database of reference dependencies is required which would represent the collection of magnitude variation patterns at different excitation frequencies, the depth and the length of delaminations. Then, the estimated reference dependency can be compared to the experimental measurement in order to extract the depth of delamination. It is presumed that the reference dependency which best matches the experiments is the one that gives the closest definition of the damage depth in the structure. The step-by-step procedure of the damage depth estimation can be outlined as follows: 1.
The source of guided waves E is driven by burst with Gaussian envelope of central frequency f 1 to introduce A 0 mode in the investigated structure.

2.
Time trace u A0 (t) is received with sensor R ref , which represents the structure without the damage (the reference signal). Meanwhile, the other receiver R 1 captures time history u IIA0 (t) which represents the response of the structure at the monitoring point beyond the defect. 3.
Ratio u A0expf 1 (t) of peak-to-peak amplitudes of waveforms u A0 (t) and u IIA0 (t) is estimated at the burst central frequency f 1 . 4.
Steps 1-3 are repeated over for other excitation frequencies where f N denotes the burst central frequencies used for the excitation of A 0 mode. 5.
The experimental dataset U A0exp (f ) is compared to the prescribed database of reference dependencies U A0ref (f,h,l) (where h, l are the depth and the length of damage, respectively). Based upon this concept, the goal is to select the single reference dependency that is the closest to the experimental data. The criteria to characterize the similarity of two datasets include standard deviation as follows: where where h opt is the value which corresponds to the delamination depth where the reference dependency best matches the experimental dataset; U A0ref (f,h,l) is the calculated reference dependency; U A0exp (f ) is the experimental dataset; h and l are the depth and the length of delamination, respectively; and "std" and "mean" denote the standard deviation and the arithmetic average, respectively.
If the structure is defect-free, function U A0exp (f ) is close to monotonic and does not feature essential amplitude variations. Meanwhile, if the variation in the magnitude ratio is observed, it can be related to the presence of a defect at a particular depth. The proposed depth estimation technique relies under the assumption, that velocities of direct A 0 mode in the upper and lower sub laminates are similar and the single wave packet is captured after the exit of delamination. In case of very large and asymmetric defects, two wave packets of A 0 mode can appear after the exit of the defect. In such a case, the inspection frequency must be decreased, or the number of signal cycles increased in order to maintain a single wave packet and to observe magnitude variation patterns.

Method to Assess the Length of Defect
If the depth of the defect is determined, its length can be estimated from the delay between "direct A 0 transmission" u IIA0 (t) and "first A 0 arrival" u IA0 (t) which can be denoted as ∆t A0 . Delay ∆t A0 is proportional to the excitation frequency and the propagation path of S 0 mode, which itself is actually determined by the length of the defect ( Figure 5). The time delay between two neighboring wave packets can be estimated according to: where u IA0 (t) and u IIA0 (t) are the wave packets of "first A 0 arrival" and "direct A 0 transmission", and HT denotes the Hilbert transform. Then, the length of the defect at a particular depth h can be calculated as: Appl. Sci. 2020, 10, 5236 The Equation (7) calculates mean defect length from a set of defect lengths obtained at different excitation frequencies. If multiple burst central frequencies are used to excite the guided waves, the group velocity values of A 0 and S 0 modes should be taken at different central frequencies, which allow one to consider the dispersion.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 19 Figure 5. Bscan image of the vertical component of particle velocity, which shows delay between converted and direct A0 modes beyond the defect.
As the "direct A0 transmission" uIIA 0 (t) is the result of the interference between the wave packets propagating above and below the defect, the maximum of the wave after the exit of delamination (uIIA 0 (t)) may be shifted and no longer correspond to the wave energy velocity. Hence the delay ΔtA 0 in Equation (6) will be estimated with some error. The best way here would be to calculate the delay according to the mean energy (mass center) of the signal [32]. However, in such case the time window covering entire signal should be selected. Such a task is not trivial due to the dispersion of A0 mode signals, presence of other modes and reflections from object boundaries. Additionally, the different frequency components of the dispersed A0 mode are redistributed in the envelope of the signal. The lower frequency components are concentrated at the tail of the signal. So, the energy of wave packet uIIA 0 (t) should be related to a particular frequency and velocity which are used for estimation of the time of flight. By using a maximum or Hilbert envelope of the signal it can be assumed that it corresponds to the dominant frequency of the signal and gives sufficient ToF accuracy.
It is noteworthy that, for the reliable operation of the method, direct uIIA 0 (t) and converted modes uIA 0 (t) must be separated in the time domain to an adequate extent in order to be able to measure the time delay between the two wave packets. Thus, there exists the minimum size of a defect that can be evaluated by using this approach. Such a minimum detectable defect is frequency dependent, so higher frequency signals can be used to improve sizing accuracy. It is important that the depth of the defect strongly influences the accuracy of the delamination length estimation, especially in the presence of significant dispersion. If the depth of the defect is unknown, such an approach becomes valid in non-dispersive regions only. However, by using the proposed delamination depth estimation technique, the thicknesses and velocities in the upper and lower sub-laminate can be determined, which leads to the improved accuracy of the delamination length evaluation.

Estimation of Reference Dependencies
As mentioned previously, the proposed technique to estimate the depth of delamination relies on reference calculations, which means that the experimental signals must be compared to the modelbased predictions. In order to extract the depth of the damage, the dependencies between the magnitude of "direct A0 transmission" UA 0 ref(f,h,l), excitation frequency (f) and depth (h) at fixed defect lengths (l) will be calculated. In this study, for the sake of simplicity, the estimation of reference dependencies shall be demonstrated on the 2D model of a GFRP plate.

Theoretical Formulation of Reference Dataset
In the case of depth assessment, our goal is to predict the waveform of A0 mode which As the "direct A 0 transmission" u IIA0 (t) is the result of the interference between the wave packets propagating above and below the defect, the maximum of the wave after the exit of delamination (u IIA0 (t)) may be shifted and no longer correspond to the wave energy velocity. Hence the delay ∆t A0 in Equation (6) will be estimated with some error. The best way here would be to calculate the delay according to the mean energy (mass center) of the signal [32]. However, in such case the time window covering entire signal should be selected. Such a task is not trivial due to the dispersion of A 0 mode signals, presence of other modes and reflections from object boundaries. Additionally, the different frequency components of the dispersed A 0 mode are redistributed in the envelope of the signal. The lower frequency components are concentrated at the tail of the signal. So, the energy of wave packet u IIA0 (t) should be related to a particular frequency and velocity which are used for estimation of the time of flight. By using a maximum or Hilbert envelope of the signal it can be assumed that it corresponds to the dominant frequency of the signal and gives sufficient ToF accuracy.
It is noteworthy that, for the reliable operation of the method, direct u IIA0 (t) and converted modes u IA0 (t) must be separated in the time domain to an adequate extent in order to be able to measure the time delay between the two wave packets. Thus, there exists the minimum size of a defect that can be evaluated by using this approach. Such a minimum detectable defect is frequency dependent, so higher frequency signals can be used to improve sizing accuracy. It is important that the depth of the defect strongly influences the accuracy of the delamination length estimation, especially in the presence of significant dispersion. If the depth of the defect is unknown, such an approach becomes valid in non-dispersive regions only. However, by using the proposed delamination depth estimation technique, the thicknesses and velocities in the upper and lower sub-laminate can be determined, which leads to the improved accuracy of the delamination length evaluation.

Estimation of Reference Dependencies
As mentioned previously, the proposed technique to estimate the depth of delamination relies on reference calculations, which means that the experimental signals must be compared to the model-based predictions. In order to extract the depth of the damage, the dependencies between the magnitude of "direct A 0 transmission" U A0ref (f,h,l), excitation frequency (f ) and depth (h) at fixed defect lengths (l) will be calculated. In this study, for the sake of simplicity, the estimation of reference dependencies shall be demonstrated on the 2D model of a GFRP plate.

Theoretical Formulation of Reference Dataset
In the case of depth assessment, our goal is to predict the waveform of A 0 mode which propagates within the defect. Equation (3) can be rewritten in a short form as: where U A0ad (f ) and U A0bd (f ) represent the propagation of A 0 mode above and below the delamination, respectively. Waveforms u A0ad (t) and u A0bd (t) can be predicted according to the theory of linear acoustics which states that the output waveform of the system is the convolution of the input signal and the system impulse response. Consequently, each of the signals u A0ad (t) and u A0bd (t) can be expressed by the equations: where u ref (t, f l ) is the theoretical incident signal at a particular excitation frequency (l = 1 ÷ M, M is the total number of frequencies used to drive the emitter); α(f ) is the attenuation function; l k is the length of the defect (k = 1 ÷ N, where N is the total number of length incrementations); c p (f,h 1m ) and c p (f,h 2m ) are the phase velocities for the layer above (h 1m ) and below the defect (h 2m ) which depend on the level of defect asymmetry; h 1m and h 2m are the thicknesses of the layer above and below the defect; and FT denotes the Fourier transform. The in phase or out of phase interference of these signals defines the peak to peak amplitude of "direct A 0 transmission" u A0ref (l k ,h 1m ,h 2m ,f l ) at excitation frequency f l : To illustrate the estimation of reference dependency, let us consider a GFRP plate of thickness h = 4 mm with delamination-type defects situated at various depths as follows: h 1m = mh/8, h 2m = h − h 1m , where m = 1, 2, . . . ,3. Let us assume that the excitation burst central frequency of the incident A 0 mode varies starting from 50 kHz to 200 kHz with increments of 1 kHz. The graphic illustration of the magnitude of A 0 mode u A0ref (l k ,h 1m , h 2m ,f l ) as the function of excitation frequency f l at discrete depths h 1m = mh/8, h 2m = −h 1m is illustrated in Figures 6 and 7. The results are presented at fixed defect lengths of l 1 = 50 mm, l 2 = 70 mm, l 3 = 90 mm and l 4 = 110 mm. Each graph in Figures 6 and 7 represent magnitude variation of forward scattered A 0 mode at different excitation frequencies. It can be observed that the magnitude variation of "direct A 0 transmission" depends not only on the depth of the defect (which is represented by solid, dashed and dash-dot lines) but on its length as well. As one might expect, the size of the defect could be estimated from the same magnitude variations. However, these magnitude variations due to the length of the defect repeat. For example, if the defect is 0.5 mm below the surface (solid line), the magnitude variation versus frequency is almost identical at defect lengths of l 3 = 90 mm (Figure 7a) and l 4 = 110 mm (Figure 7b). It means that multiple defect lengths correspond to the same magnitude variation of A 0 mode. On the other hand, the observed magnitude variations due to the depth are always different. Therefore, in this paper it is suggested to use magnitude variation for damage depth estimation only and then apply time of flight measurement techniques to estimate the size of defect. By knowing the depth of the damage, the accuracy of damage sizing can be increased, as the actual depth of defect determines the phase velocities of A 0 and S 0 modes in the upper and lower sub-laminates.

Numerical and Experimental Verification of Reference Dependencies
In the following chapter, the reference dependency estimation technique presented in the previous section is verified with the appropriate numerical simulations and experiments on GFRP and aluminum samples. To achieve the aim of this study, the numerically and experimentally estimated amplitude variations in the case of different lengths of a defect are compared with the predictions calculated by Equation (11). It is important to distinguish that in this verification, magnitude variations of forward scattered A0 mode are being estimated at different lengths of defect as it was impossible to iteratively change the depth of defect in experiments. The proposed damage detection technique itself uses magnitude variation for damage depth assessment, while the length is extracted from time of flight (ToF) measurements.

Numerical Verification
At first, the reference dependency estimation technique is verified numerically. A total number of 50 2D linear structural mechanics FE models were employed for an anisotropic GFRP plate with the dimensions of x = 2000 mm, y = 4 mm. In this case, the plain strain model was implemented by assuming that the object is infinite in the z direction. In each of the 50 models, delamination-type

Numerical and Experimental Verification of Reference Dependencies
In the following chapter, the reference dependency estimation technique presented in the previous section is verified with the appropriate numerical simulations and experiments on GFRP and aluminum samples. To achieve the aim of this study, the numerically and experimentally estimated amplitude variations in the case of different lengths of a defect are compared with the predictions calculated by Equation (11). It is important to distinguish that in this verification, magnitude variations of forward scattered A0 mode are being estimated at different lengths of defect as it was impossible to iteratively change the depth of defect in experiments. The proposed damage detection technique itself uses magnitude variation for damage depth assessment, while the length is extracted from time of flight (ToF) measurements.

Numerical Verification
At first, the reference dependency estimation technique is verified numerically. A total number of 50 2D linear structural mechanics FE models were employed for an anisotropic GFRP plate with the dimensions of x = 2000 mm, y = 4 mm. In this case, the plain strain model was implemented by assuming that the object is infinite in the z direction. In each of the 50 models, delamination-type

Numerical and Experimental Verification of Reference Dependencies
In the following chapter, the reference dependency estimation technique presented in the previous section is verified with the appropriate numerical simulations and experiments on GFRP and aluminum samples. To achieve the aim of this study, the numerically and experimentally estimated amplitude variations in the case of different lengths of a defect are compared with the predictions calculated by Equation (11). It is important to distinguish that in this verification, magnitude variations of forward scattered A 0 mode are being estimated at different lengths of defect as it was impossible to iteratively change the depth of defect in experiments. The proposed damage detection technique itself uses magnitude variation for damage depth assessment, while the length is extracted from time of flight (ToF) measurements.

Numerical Verification
At first, the reference dependency estimation technique is verified numerically. A total number of 50 2D linear structural mechanics FE models were employed for an anisotropic GFRP plate with the dimensions of x = 2000 mm, y = 4 mm. In this case, the plain strain model was implemented by assuming that the object is infinite in the z direction. In each of the 50 models, delamination-type defects were introduced 1.5 mm below the surface (h 1 = 1.5 mm, h 2 = 2.5 mm) by separating the appropriate nodes of the mesh. The length of delamination l varied from 10 to 500 mm with an increment of 10 mm, thus yielding a total number of 50 numerical models. The delamination was centered horizontally along the x axis of the investigated sample for all the defect lengths. Throughout the simulations, 2D structural solid plane-42 finite elements were used. The spatial size of the element was equal to 0.5 mm, which corresponds to 26 nodes per wavelength if the slowest A 0 mode at 80 kHz central frequency is considered. The excitation signal in all cases was a Gaussian envelope tone burst of three cycles and a central frequency of 80 kHz. Such an excitation regime was deliberately selected to converge with the experimental verification, which will be described further in this section. A 0 mode was introduced by applying out-of-plane displacements to the first node column of the model located on the sample edge (see Figure 8). Such a force excites both fundamental A 0 and S 0 modes at a given frequency band and acts similarly to a simple thickness mode transducer. The integration step in the time domain was 0.625 µs, which is 1/20 of the period at 80 kHz central frequency. Throughout the simulation, a generalized stiffness matrix was used to describe a symmetric GFRP laminate (biaxial: 0 • and 90 • /bias: ±45 • /biaxial: 0 • and 90 • ) S , while employing the research conducted by Ryan et al. [33] Hence, the stiffness matrix was obtained by summing up the contribution of each ply in terms of their respective thicknesses. The following material properties of GFRP were used: Young's modulus (E x = 10 GPa, E z = 35.7 GPa); Poisson's ratio (υ xz = 0.325, υ zx = 0.091, υ yx = 0.35); shear modulus G xz = 2.8 GPa; density ρ = 1800 kg/m 3 . The sketch of the numerical model used in this study is presented in Figure 8. defects were introduced 1.5 mm below the surface (h1 = 1.5 mm, h2 = 2.5 mm) by separating the appropriate nodes of the mesh. The length of delamination l varied from 10 to 500 mm with an increment of 10 mm, thus yielding a total number of 50 numerical models. The delamination was centered horizontally along the x axis of the investigated sample for all the defect lengths. Throughout the simulations, 2D structural solid plane-42 finite elements were used. The spatial size of the element was equal to 0.5 mm, which corresponds to 26 nodes per wavelength if the slowest A0 mode at 80 kHz central frequency is considered. The excitation signal in all cases was a Gaussian envelope tone burst of three cycles and a central frequency of 80 kHz. Such an excitation regime was deliberately selected to converge with the experimental verification, which will be described further in this section. A0 mode was introduced by applying out-of-plane displacements to the first node column of the model located on the sample edge (see Figure 8). Such a force excites both fundamental A0 and S0 modes at a given frequency band and acts similarly to a simple thickness mode transducer. The integration step in the time domain was 0.625 μs, which is 1/20 of the period at 80 kHz central frequency. Throughout the simulation, a generalized stiffness matrix was used to describe a symmetric GFRP laminate (biaxial: 0° and 90°/bias: ±45°/biaxial: 0° and 90°)S, while employing the research conducted by Ryan et al. [33] Hence, the stiffness matrix was obtained by summing up the contribution of each ply in terms of their respective thicknesses. The following material properties of GFRP were used: Young's modulus (Ex = 10 GPa, Ez = 35.7 GPa); Poisson's ratio (υxz = 0.325, υzx = 0.091, υyx = 0.35); shear modulus Gxz = 2.8 GPa; density ρ = 1800 kg/m 3 . The sketch of the numerical model used in this study is presented in Figure 8. The variable monitored in this study was the vertical component of particle velocity (y) at the nodal locations before (x1 = 500 mm, y1 = 4 mm) and beyond (x2 = 1500 mm, y2 = 4 mm) the defect (see Figure 8). The obtained peak to peak amplitude variations versus the length of delamination are presented in Figure 9. In the same figure, the predicted curve obtained by using Equation (11) with identical parameters is presented as the solid line.  The variable monitored in this study was the vertical component of particle velocity (y) at the nodal locations before (x 1 = 500 mm, y 1 = 4 mm) and beyond (x 2 = 1500 mm, y 2 = 4 mm) the defect (see Figure 8). The obtained peak to peak amplitude variations versus the length of delamination are presented in Figure 9. In the same figure, the predicted curve obtained by using Equation (11) with identical parameters is presented as the solid line.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 19 defects were introduced 1.5 mm below the surface (h1 = 1.5 mm, h2 = 2.5 mm) by separating the appropriate nodes of the mesh. The length of delamination l varied from 10 to 500 mm with an increment of 10 mm, thus yielding a total number of 50 numerical models. The delamination was centered horizontally along the x axis of the investigated sample for all the defect lengths. Throughout the simulations, 2D structural solid plane-42 finite elements were used. The spatial size of the element was equal to 0.5 mm, which corresponds to 26 nodes per wavelength if the slowest A0 mode at 80 kHz central frequency is considered. The excitation signal in all cases was a Gaussian envelope tone burst of three cycles and a central frequency of 80 kHz. Such an excitation regime was deliberately selected to converge with the experimental verification, which will be described further in this section. A0 mode was introduced by applying out-of-plane displacements to the first node column of the model located on the sample edge (see Figure 8). Such a force excites both fundamental A0 and S0 modes at a given frequency band and acts similarly to a simple thickness mode transducer. The integration step in the time domain was 0.625 μs, which is 1/20 of the period at 80 kHz central frequency. Throughout the simulation, a generalized stiffness matrix was used to describe a symmetric GFRP laminate (biaxial: 0° and 90°/bias: ±45°/biaxial: 0° and 90°)S, while employing the research conducted by Ryan et al. [33] Hence, the stiffness matrix was obtained by summing up the contribution of each ply in terms of their respective thicknesses. The following material properties of GFRP were used: Young's modulus (Ex = 10 GPa, Ez = 35.7 GPa); Poisson's ratio (υxz = 0.325, υzx = 0.091, υyx = 0.35); shear modulus Gxz = 2.8 GPa; density ρ = 1800 kg/m 3 . The sketch of the numerical model used in this study is presented in Figure 8. The variable monitored in this study was the vertical component of particle velocity (y) at the nodal locations before (x1 = 500 mm, y1 = 4 mm) and beyond (x2 = 1500 mm, y2 = 4 mm) the defect (see Figure 8). The obtained peak to peak amplitude variations versus the length of delamination are presented in Figure 9. In the same figure, the predicted curve obtained by using Equation (11) with identical parameters is presented as the solid line.

Experimental Verification
In addition to numerical calculations, experiments on a custom-made aluminum alloy plate with the dimensions of 1250 mm × 600 mm × 5 mm were carried out as well. The aluminum plate was selected for the verification of the method due to a few reasons. At first, the properties of aluminum and the dispersion relations of guided waves in it are well known. The second reason is that it was the most reliable way to create thin delamination with a strictly defined length. Afterwards, the method can be applied on any other types of materials with known properties. In the experiments, a special aluminum sample with a machined groove was produced in order to obtain the response adequate to that obtained with the delamination-type defect. In the course of the production of the sample, two aluminum sheets with thicknesses of 2 mm and 3 mm were bonded adhesively to produce a total thickness of 5 mm. Prior to bonding, the 3 mm sheet was machined to produce a wedge-shaped groove with a depth of 500 µm. The lateral dimensions of the groove ranged from 35 to 335 mm, which enabled the delamination-type defect of a variable size to be simulated. The sketch of the experimental setup used for validation purposes is presented in Figure 10a.
The experimental investigation was conducted by using custom-made low frequency ultrasonic transducers (130 kHz) operating in a thickness mode with an active diameter of 10 mm. The transducers possess the bandwidth of 100 kHz at −6 dB level and were driven by a 50 V tone-burst with a Gaussian envelope of 3 periods and a central frequency of 130 kHz. The number of cycles were selected to maintain a good balance between the amplitude measurements for depth estimation which requires narrowband excitation and ToF measurements for defect length assessment where a wideband excitation is preferred for better resolution. The best compromise between the propagation distance, the number of cycles, the wavelength and the resolution was selected according to the concept of the minimum resolvable distance (MRD), which was described by Wilcox et al. [34] In the experiment, the transducers were arranged in a pitch-catch configuration and attached perpendicularly to the surface of the sample to introduce A 0 mode. The transmitter and the receiver were linked together and scanned along the short edge of the artificial defect as it is shown in Figure 10a. To achieve a reliable and uniform acoustic contact between the transducer and the specimen, custom-made spring type adjusters were used (Figure 10b). As a result of the experiment, a set of waveforms at different defect lengths were collected. The peak-to-peak amplitude of each time trace was compared to the reference time history captured at the same basis on the defect-free region.

Experimental Verification
In addition to numerical calculations, experiments on a custom-made aluminum alloy plate with the dimensions of 1250 mm × 600 mm × 5 mm were carried out as well. The aluminum plate was selected for the verification of the method due to a few reasons. At first, the properties of aluminum and the dispersion relations of guided waves in it are well known. The second reason is that it was the most reliable way to create thin delamination with a strictly defined length. Afterwards, the method can be applied on any other types of materials with known properties. In the experiments, a special aluminum sample with a machined groove was produced in order to obtain the response adequate to that obtained with the delamination-type defect. In the course of the production of the sample, two aluminum sheets with thicknesses of 2 mm and 3 mm were bonded adhesively to produce a total thickness of 5 mm. Prior to bonding, the 3 mm sheet was machined to produce a wedge-shaped groove with a depth of 500 μm. The lateral dimensions of the groove ranged from 35 to 335 mm, which enabled the delamination-type defect of a variable size to be simulated. The sketch of the experimental setup used for validation purposes is presented in Figure 10a.
The experimental investigation was conducted by using custom-made low frequency ultrasonic transducers (130 kHz) operating in a thickness mode with an active diameter of 10 mm. The transducers possess the bandwidth of 100 kHz at −6 dB level and were driven by a 50 V tone-burst with a Gaussian envelope of 3 periods and a central frequency of 130 kHz. The number of cycles were selected to maintain a good balance between the amplitude measurements for depth estimation which requires narrowband excitation and ToF measurements for defect length assessment where a wideband excitation is preferred for better resolution. The best compromise between the propagation distance, the number of cycles, the wavelength and the resolution was selected according to the concept of the minimum resolvable distance (MRD), which was described by Wilcox et al. [34] In the experiment, the transducers were arranged in a pitch-catch configuration and attached perpendicularly to the surface of the sample to introduce A0 mode. The transmitter and the receiver were linked together and scanned along the short edge of the artificial defect as it is shown in Figure  10a. To achieve a reliable and uniform acoustic contact between the transducer and the specimen, custom-made spring type adjusters were used (Figure 10b). As a result of the experiment, a set of waveforms at different defect lengths were collected. The peak-to-peak amplitude of each time trace was compared to the reference time history captured at the same basis on the defect-free region. The comparison between mathematical predictions and experimental measurements can be seen in Figure 11. The results demonstrate a fairly good agreement between the experimental values and the mathematical predictions, especially below the defect lengths of 200 mm. The mismatch in the amplitude ratio between the mathematical and experimental estimations may appear due to the The comparison between mathematical predictions and experimental measurements can be seen in Figure 11. The results demonstrate a fairly good agreement between the experimental values and the mathematical predictions, especially below the defect lengths of 200 mm. The mismatch in the amplitude ratio between the mathematical and experimental estimations may appear due to the mismatch of dispersion curves, ignored transmission and reflection coefficients in Equations (9) and (10), and by diffraction losses. The discrepancies between magnitude variations and length of defect that are visible in Figures 9 and 11 have a minor impact for damage depth estimation. The actual size of the defect is obtained from ToF measurement, using the damage depth information from magnitude variation of forward scattered direct A 0 mode.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 19 mismatch of dispersion curves, ignored transmission and reflection coefficients in Equations (9) and (10), and by diffraction losses. The discrepancies between magnitude variations and length of defect that are visible in Figures 9 and 11 have a minor impact for damage depth estimation. The actual size of the defect is obtained from ToF measurement, using the damage depth information from magnitude variation of forward scattered direct A0 mode.

Verification with 2D Model
To demonstrate the overall performance of the proposed delamination depth and length estimation method, two separate 2D GFRP finite element models, both with dimensions of 2000 mm × 4 mm shall be considered. For the sake of better understanding, let us denote them as sample No. 1 and sample No. 2. Each of the considered models contains delaminations of different sizes and through-thickness locations. In sample No. 1, the delamination with a length of x1 = 70 mm is introduced 1 mm below the top surface (h11 = 1 mm, h21 = 3 mm). Meanwhile, sample No. 2 has a defect of x2 = 90 mm located at a depth of 1.5 mm (h12 = 1.5 mm, h22 = 2.5 mm). A sketch of the considered cases is presented in Figure 12.
Throughout the simulations, 2D structural solid plane-42 finite elements were used. The same properties of GFRP laminate were maintained as described in the previous section. In each of the investigated cases, a Gaussian envelope tone burst of three cycles was used for excitation. The excitation frequencies were changed within the range from f1 = 50 kHz to fn = 200 kHz with increments of 10 kHz. Hence, 16 different excitation frequencies were used for each of the two considered models. The spatial size of the element in all cases was constant and equal to 0.5 mm, which corresponds to 12 up to 37 nodes per wavelength for the slowest A0 in a frequency range of f1, f2, …, fn. A0 mode was introduced into the structure by applying the out-of-plane nodal displacements to the edge of the sample in the same way as it was demonstrated in Figure 8. The integration step in the time domain varied from 1 μs for the lowest excitation frequency f1 = 50 kHz to 0.25 μs for fn = 200 kHz. The variable monitored in this study was a vertical component of the particle velocity (y) at nodal locations r1 (x1 = 750 mm, y1 = 4 mm) and r2 (x2 = 1250 mm, y2 = 4 mm).
(a) (b) Figure 11. Magnitude variation of A 0 mode due to the size of defect at fixed depth: comparison between analytical calculations and experiment.

Verification with 2D Model
To demonstrate the overall performance of the proposed delamination depth and length estimation method, two separate 2D GFRP finite element models, both with dimensions of 2000 mm × 4 mm shall be considered. For the sake of better understanding, let us denote them as sample No. 1 and sample No. 2. Each of the considered models contains delaminations of different sizes and through-thickness locations. In sample No. 1, the delamination with a length of x 1 = 70 mm is introduced 1 mm below the top surface (h 11 = 1 mm, h 21 = 3 mm). Meanwhile, sample No. 2 has a defect of x 2 = 90 mm located at a depth of 1.5 mm (h 12 = 1.5 mm, h 22 = 2.5 mm). A sketch of the considered cases is presented in Figure 12.
magnitude variation of forward scattered direct A0 mode.

Verification with 2D Model
To demonstrate the overall performance of the proposed delamination depth and length estimation method, two separate 2D GFRP finite element models, both with dimensions of 2000 mm × 4 mm shall be considered. For the sake of better understanding, let us denote them as sample No. 1 and sample No. 2. Each of the considered models contains delaminations of different sizes and through-thickness locations. In sample No. 1, the delamination with a length of x1 = 70 mm is introduced 1 mm below the top surface (h11 = 1 mm, h21 = 3 mm). Meanwhile, sample No. 2 has a defect of x2 = 90 mm located at a depth of 1.5 mm (h12 = 1.5 mm, h22 = 2.5 mm). A sketch of the considered cases is presented in Figure 12.
Throughout the simulations, 2D structural solid plane-42 finite elements were used. The same properties of GFRP laminate were maintained as described in the previous section. In each of the investigated cases, a Gaussian envelope tone burst of three cycles was used for excitation. The excitation frequencies were changed within the range from f1 = 50 kHz to fn = 200 kHz with increments of 10 kHz. Hence, 16 different excitation frequencies were used for each of the two considered models. The spatial size of the element in all cases was constant and equal to 0.5 mm, which corresponds to 12 up to 37 nodes per wavelength for the slowest A0 in a frequency range of f1, f2, …, fn. A0 mode was introduced into the structure by applying the out-of-plane nodal displacements to the edge of the sample in the same way as it was demonstrated in Figure 8. The integration step in the time domain varied from 1 μs for the lowest excitation frequency f1 = 50 kHz to 0.25 μs for fn = 200 kHz. The variable monitored in this study was a vertical component of the particle velocity (y) at nodal locations r1 (x1 = 750 mm, y1 = 4 mm) and r2 (x2 = 1250 mm, y2 = 4 mm). Throughout the simulations, 2D structural solid plane-42 finite elements were used. The same properties of GFRP laminate were maintained as described in the previous section. In each of the investigated cases, a Gaussian envelope tone burst of three cycles was used for excitation. The excitation frequencies were changed within the range from f 1 = 50 kHz to f n = 200 kHz with increments of 10 kHz. Hence, 16 different excitation frequencies were used for each of the two considered models. The spatial size of the element in all cases was constant and equal to 0.5 mm, which corresponds to 12 up to 37 nodes per wavelength for the slowest A 0 in a frequency range of f 1 , f 2 , . . . , f n . A 0 mode was introduced into the structure by applying the out-of-plane nodal displacements to the edge of the sample in the same way as it was demonstrated in Figure 8. The integration step in the time domain varied from 1 µs for the lowest excitation frequency f 1 = 50 kHz to 0.25 µs for f n = 200 kHz. The variable monitored in this study was a vertical component of the particle velocity (y) at nodal locations r 1 (x 1 = 750 mm, y 1 = 4 mm) and r 2 (x 2 = 1250 mm, y 2 = 4 mm).

Estimation of Defect Depth
In order to estimate the depth of the defect, the procedure described in Section 2 is employed. The reference dependency is calculated according to Equations (9)-(11) for defect lengths l 1 = 10 mm, l 2 = 11 mm, . . . , l N = 100 mm and defect depths h 11 = 0.5 mm, h 12 = 1 mm and h 13 = 1.5 mm for the same set of excitation frequencies (f 1 = 50 kHz, f 2 = 60 kHz, . . . , f n = 200 kHz). Then, h opt is calculated by using Equation (4) in order to estimate the depth of the damage. The graphic illustration of u std (h,l) at different lengths and depths is presented in Figure 13. The criteria to estimate the actual depth of the defect are based on the selection of the minimum of u std (h,l It can be observed that, for sample No. 1 (Figure 13a), hopt is close to zero for depth h12 = 1 mm (the dash-dot line in Figure 13a). Hence, the estimated depth of delamination for sample No. 1 is 1 mm, which corresponds to the actual depth of the damage. Similar observations may be made for sample No. 2. hopt possesses the lowest values at depth h13 = 1.5 mm (the dashed line in Figure 12b), which is in agreement with the actual depth of the damage. In the presented case, the actual depth of the damage was included in the database. In order to evaluate the performance of the method in case of absence of actual depth in database, another set of magnitude patterns has been generated with depth increments of 0.4 mm (h11 = 0.4 mm, h12 = 0.8 mm, h13 = 1.2 mm, h14 = 1.6 mm, h15 = 2 mm). In the case for sample No. 1 the minimum values of the hopt were close for both h12 = 0.8 mm and h13 = 1.  It can be observed that, for sample No. 1 (Figure 13a), h opt is close to zero for depth h 12 = 1 mm (the dash-dot line in Figure 13a). Hence, the estimated depth of delamination for sample No. 1 is 1 mm, which corresponds to the actual depth of the damage. Similar observations may be made for sample No. 2. h opt possesses the lowest values at depth h 13 = 1.5 mm (the dashed line in Figure 12b), which is in agreement with the actual depth of the damage. In the presented case, the actual depth of the damage was included in the database. In order to evaluate the performance of the method in case of absence of actual depth in database, another set of magnitude patterns has been generated with depth increments of 0.4 mm (h 11 = 0.4 mm, h 12 = 0.8 mm, h 13 = 1.2 mm, h 14 = 1.6 mm, h 15 = 2 mm). In the case for sample No. 1 the minimum values of the h opt were close for both h 12 = 0.8 mm and h 13 = 1.2 mm, however the absolute minimum was obtained for depth h 12 = 0.8 mm. For sample No. 2 h opt possesses the lowest values at depth h 14 = 1.6 mm. The amount of error of course depends on the depth increments in the database. Such database in general is computationally inexpensive (in this case U A0ref (11 × 5 × 90), hence can be calculated small depth increments according to the desired depth estimation accuracy.

Evaluation of Defect Length
As the depth of the defect for both samples has been determined, the length is estimated by measuring the delay between the "first A 0 arrival" u IA0 (t) and "direct A 0 transmission" u IIA0 (t) of A 0 mode. The delays are measured according to the maximum envelope values of each wave packet by using Equation (6). Then, the length of the defect is estimated according to Equation (7).
In this particular case, the delay values are estimated at a single excitation frequency f 9 = 130 kHz. The vertical component of the particle velocity at nodal point r 2 for sample No. 1 and sample No. 2 is presented in Figure 14a,b. The picture shows the wave packets of "first A 0 arrival" u IA0 (t) and "direct A 0 transmission" u IIA0 (t), along with the estimated delay values according to the Hilbert envelope.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 15 of 19 Figure 13. Deviation values of difference between the experimental and reference dependency: the results for sample No. 1 with a defect located 1 mm below the surface (a), the results for sample No. 2 with delamination situated at a depth of 1.5 mm (b).

Evaluation of Defect Length
As the depth of the defect for both samples has been determined, the length is estimated by measuring the delay between the "first A0 arrival" uIA 0 (t) and "direct A0 transmission" uIIA 0 (t) of A0 mode. The delays are measured according to the maximum envelope values of each wave packet by using Equation (6). Then, the length of the defect is estimated according to Equation (7).
In this particular case, the delay values are estimated at a single excitation frequency f9 = 130 kHz. The vertical component of the particle velocity at nodal point r2 for sample No. 1 and sample No. 2 is presented in Figure 14a,b. The picture shows the wave packets of "first A0 arrival" uIA 0 (t) and "direct A0 transmission" uIIA 0 (t), along with the estimated delay values according to the Hilbert envelope.
According to Equation (7), the delay time values correspond to the defect length of 67.6 mm for sample No. 1, meanwhile, for sample No. 2, the delay indicates the length of 83.4 mm. The results demonstrate that the average error in the length estimation is approximately 6% relative to the actual size. Such error is valid for frequency range between 50 and 200 kHz for defect sizes of 70 mm and 90 mm. The amount of error depends both on the accuracy of damage depth estimation and the delay time measurement. For different defect sizes and excitation frequencies, the errors may change. However, the accuracy of defect size and depth estimation can still be improved by adding scattering coefficients to the reference model and measuring delay time according to the mass center of the signal instead of using Hilbert envelope.

Verification with 3D Model
In order to verify the applicability of the method for the real measurement, the 3D case was investigated. The reason for that is that in the case of 2D modeling of guided waves, the diffraction does not affect the signals, in contrary to the 3D case when diffraction takes place, and signals can be overlapped with boundary reflections. There are many signal processing methods which work perfectly in the case of 1D or 2D models; however, they still fail when they are applied for real structures or 3D models. So, the verification of the method for a 3D case is very important for the assessment of the potential of the method for in situ applications.
For this purpose, a 3D GFRP numerical model with the dimensions of 600 mm × 300 mm × 4 mm was employed with a square (90 mm × 90 mm) delamination-type defect, located at 1.6 mm below the top surface. The defect was centered both according to the x and z axes ( Figure 15). The structure was According to Equation (7), the delay time values correspond to the defect length of 67.6 mm for sample No. 1, meanwhile, for sample No. 2, the delay indicates the length of 83.4 mm. The results demonstrate that the average error in the length estimation is approximately 6% relative to the actual size. Such error is valid for frequency range between 50 and 200 kHz for defect sizes of 70 mm and 90 mm. The amount of error depends both on the accuracy of damage depth estimation and the delay time measurement. For different defect sizes and excitation frequencies, the errors may change. However, the accuracy of defect size and depth estimation can still be improved by adding scattering coefficients to the reference model and measuring delay time according to the mass center of the signal instead of using Hilbert envelope.

Verification with 3D Model
In order to verify the applicability of the method for the real measurement, the 3D case was investigated. The reason for that is that in the case of 2D modeling of guided waves, the diffraction does not affect the signals, in contrary to the 3D case when diffraction takes place, and signals can be overlapped with boundary reflections. There are many signal processing methods which work perfectly in the case of 1D or 2D models; however, they still fail when they are applied for real structures or 3D models. So, the verification of the method for a 3D case is very important for the assessment of the potential of the method for in situ applications.
For this purpose, a 3D GFRP numerical model with the dimensions of 600 mm × 300 mm × 4 mm was employed with a square (90 mm × 90 mm) delamination-type defect, located at 1.6 mm below the top surface. The defect was centered both according to the x and z axes ( Figure 15). The structure was meshed by using SOLID64 elements which had eight nodes each with three degrees of freedom. The properties of the GFRP structure were the same as it was described in the previous sections. The spatial size of the element was set to 0.8 mm, which corresponds to 16 nodes per wavelength for A 0 mode at 90 kHz central frequency. In total, four 3D models were created, each with a different excitation frequency which was varying from 90 kHz to 150 kHz with an increment of 20 kHz. The excitation zone was a square with dimensions of 10 mm by 10 mm. It was centered according to the x axis and had an offset of 200 mm from the origin of the z axis. The excitation signal in all the cases was a Gaussian envelope tone burst of three cycles in order to correspond to the previous experiments. A 0 mode was introduced by applying the normal out-of-plane nodal displacements to the top surface of the sample at a predefined area of 10 mm 2 . The integration step in the time domain was 0.625 µs, which is 1/18 of the period at 90 kHz central frequency. The sketch of the numerical model used in this study is presented in Figure 15.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 16 of 19 meshed by using SOLID64 elements which had eight nodes each with three degrees of freedom. The properties of the GFRP structure were the same as it was described in the previous sections. The spatial size of the element was set to 0.8 mm, which corresponds to 16 nodes per wavelength for A0 mode at 90 kHz central frequency. In total, four 3D models were created, each with a different excitation frequency which was varying from 90 kHz to 150 kHz with an increment of 20 kHz. The excitation zone was a square with dimensions of 10 mm by 10 mm. It was centered according to the x axis and had an offset of 200 mm from the origin of the z axis. The excitation signal in all the cases was a Gaussian envelope tone burst of three cycles in order to correspond to the previous experiments. A0 mode was introduced by applying the normal out-of-plane nodal displacements to the top surface of the sample at a predefined area of 10 mm 2 . The integration step in the time domain was 0.625 μs, which is 1/18 of the period at 90 kHz central frequency. The sketch of the numerical model used in this study is presented in Figure 15.
The magnitude variation of A0 mode was estimated at various excitation frequencies (f1 = 90 kHz, f2 = 110 kHz, f3 = 130 kHz, f4 = 150 kHz) and compared to the database of predictions at different defect depths (h11 = 0.4 mm, h12 = 0.8 mm, h13 = 1.2 mm, h14 = 1.6 mm) and lengths (l1, l2, … ln). The comparison between the measured magnitude variation from the 3D model and the prediction at the defect size of 90 mm is presented in Figure 16. As it can be observed, the approximated curve which goes through the measured magnitude values (the solid line) is close to the prediction at depth h14 (the dash double-dot line). This depth of the damage corresponds to the actual depth used in the model. To measure the match between the measurement and the reference dependency, the parameter hopt described by Equation (4) was used. The results are summarized in Table 2.   The magnitude variation of A 0 mode was estimated at various excitation frequencies (f 1 = 90 kHz, f 2 = 110 kHz, f 3 = 130 kHz, f 4 = 150 kHz) and compared to the database of predictions at different defect depths (h 11 = 0.4 mm, h 12 = 0.8 mm, h 13 = 1.2 mm, h 14 = 1.6 mm) and lengths (l 1 , l 2 , . . . l n ). The comparison between the measured magnitude variation from the 3D model and the prediction at the defect size of 90 mm is presented in Figure 16. As it can be observed, the approximated curve which goes through the measured magnitude values (the solid line) is close to the prediction at depth h 14 (the dash double-dot line). This depth of the damage corresponds to the actual depth used in the model. To measure the match between the measurement and the reference dependency, the parameter h opt described by Equation (4) was used. The results are summarized in Table 2.    The results presented in Table 2 show that the correct damage depth estimation was observed by using the parameter h opt . For the damage length estimation, the delay between direct and converted A 0 modes was estimated at a frequency of f 4 = 150 kHz. The measured delay was equal to ∆t A0 = 38.1 µs. It was compared to the estimated delay at a defect depth of h 14 = 1.4 mm. It can be observed that the estimated delay values correspond to the defect length of 84.1 mm. In this case, the length estimation error is approximately 6.5% relative to the actual size of the defect. We should note that, in the 3D case, the size of the defect is evaluated in the direction of wave propagation only.

Conclusions
A model-based method to estimate the absolute depth and length of delamination was proposed which is based on the magnitude variation pattern and the time-of-flight measurement of A 0 mode at different excitation frequencies.
A reference dependency estimation technique was developed which allows predicting the magnitude variation patterns of A 0 mode at different depths and lengths of a delamination-type defect. It was shown that, in comparison with the experimental measurements, the reference dependencies can be used to extract the absolute depth of the delamination. The proposed delamination depth estimation technique exploits the phase velocity dependency of A 0 mode on the thickness of the sub-laminate, thus, different magnitude variation patterns at different excitation frequencies are observed. As a result, the proposed method provides improved depth estimation accuracy even in the presence of significant dispersion. Furthermore, finite element verification demonstrated that the defect depth can be estimated with a high accuracy even in those cases when the thickness of the sub-laminate that is created by the presence of the defect is relatively thin in comparison with the wavelength of A 0 mode. The performance of the proposed method was demonstrated and verified in detecting and describing delamination-type defects of various lengths and depths.
It was demonstrated that the length of delamination at a fixed depth can be reliably estimated from the delay time between the converted and direct A 0 modes by exploiting the conversion of the incident A 0 mode to S 0 mode within the defect. The proposed approach uses previously obtained information about the depth of the defect, which allows for the determination of the phase velocities of A 0 and S 0 modes in the upper and lower sub-laminates. The results of our numerical validation demonstrated that the size of the defect can be reconstructed with an error of approximately 6%.
The main value of the proposed technique lies in multi-frequency excitation, which leads to different magnitude patterns of A 0 mode for delamination depth assessment and more accurate velocity values for defect length estimation. Hence, such a technique can be used even in the presence of the dispersion of A 0 mode and provide information about the absolute depth and length of the damage by employing a few time traces measured at slightly different frequencies.