Terahertz Spectroscopic Characterization and Thickness Evaluation of Internal Delamination Defects in GFRP Composites

The use of terahertz time-domain spectroscopy (THz-TDS) for the nondestructive testing and evaluation (NDT&E) of materials and structural systems has attracted significant attention over the past two decades due to its superior spatial resolution and capabilities of detecting and characterizing defects and structural damage in non-conducting materials. In this study, the THz-TDS system is used to detect, localize and evaluate hidden multi-delamination defects (i.e., a three-level multi-delamination system) in multilayered GFRP composite laminates. To obtain accurate results, a wavelet shrinkage de-noising algorithm is used to remove the noise from the measured time-of-flight (TOF) signals. The thickness and location of each delamination defect in the z-direction (i.e., through-the-thickness direction) are calculated from the de-noised TOF signals considering the interaction between the pulsed THz waves and the different interfaces in the GFRP composite laminates. A comparison between the actual and the measured thickness values of the delamination defects before and after the wavelet shrinkage denoising process indicates that the latter provides better results with less than 3.712% relative error, while the relative error of the non-de-noised signals reaches 16.388%. Also, the power and absorbance levels of the THz waves at every interface with different refractive indices in the GFRP composite laminates are evaluated based on analytical and experimental approaches. The present study provides an adequate theoretical analysis that could help NDT&E specialists to estimate the maximum thickness of GFRP composite materials and/or structures with different interfaces that can be evaluated by the THz-TDS. Also, the accuracy of the obtained results highlights the capabilities of the THz-TDS for the NDT&E of multilayered GFRP composite laminates.


Introduction
Glass fiber reinforced polymer-matrix (GFRP) composites have the advantages of combining several types of properties that are not usually found together in a single material (e.g. lightweight, high strength, low thermal conductivity, high resistance to corrosion, good thermal and electrical insulation performance, etc.), making them some of the most widely used materials for the manufacturing of structural parts in the aerospace, missile engine housing, electrical insulators, and construction industries [1][2][3]. However, both the design process and the Page 2 of 21 Nsengiyumva et al. Chinese Journal of Mechanical Engineering (2023) 36:6 material composition make GFRP composites vulnerable to several types of defects and structural damage such as cracks, fiber breakage, voids, and delamination which unavoidably affect their design performance during their in-service stage. These types of flaws are frequently induced by several factors including internal stress release, poor material handling during the manufacturing process, water vaporization, burns from high temperatures, impact damage, etc. In this context, structural components made of GFRP composites such as wind turbine aerofoils, and hydrofoils for marine propulsors, as well as helicopter rotor blades, and aircraft auxiliary fuel tanks [4,5], must undergo extensive quality evaluation tests after the manufacturing process to ensure all defects are eliminated prior to being introduced into service and must continuously be monitored during operation to extend their service lives [2,6]. To achieve this goal while equally preserving the structural integrity of the test structure, nondestructive testing and evaluation (NDT&E) techniques must be used. To achieve the aforementioned goals, various NDT&E techniques capable of characterizing damage and defects in fiber-reinforced composites have been developed. In our recent study [2], we indicated that ultrasonic testing is one of the most widely used methods to evaluate defects and structural damage in fiber-reinforced composite materials, and it performs well in detecting both surface and deeply buried damage and defects. However, this method does have several limitations as it requires a coupling medium (except for airborne ultrasonic), lacks the sensitivity to shallow surface fracture defects, and is subjected to large attenuations of the sound waves when propagating in multilayered composite materials such as GFRP composites [5,7]. Additional NDT&E techniques have also been used for the testing of GFRP composite materials, including X-ray radiography, infrared thermography (IRT), optical coherence tomography (OCT), optical interferometric techniques (e.g., moiré interferometry, holographic interferometry, speckle interferometry, electronic speckle pattern interferometry (ESPI), digital speckle correlation, and shearography and their subsequent variations), acoustic emission (AE), vibration testing, strain monitoring, and microwave, etc. [2,3]. Each NDT&E technique has its characteristics, but most methods have some limitations in large-scale sensing, imaging, and corresponding safety considerations [2,6]. To solve this challenge and provide nondestructive testing (NDT) engineers and practitioners with the most up-to-date NDT systems, some NDT&E techniques are currently being combined to improve the quality and accuracy of the test results [6] and guarantee the adequate performance of the GFRP composite materials.
In the field of materials and structural systems testing, terahertz (THz) technology has recently become one of the most promising NDT techniques, thanks to the adaptive physical modeling of the transmission process of the THz waves and their high-penetration capabilities in non-conducting materials [8][9][10]. The THz part of the electromagnetic spectrum extends from 0.1 to 10 THz (i.e., corresponding to the wavelength ranging from 3 mm down to 30 μm), and lies between the microwave and infrared regions of the electromagnetic spectrum. Figure 1 illustrates the region of the electromagnetic spectrum covered by the THz waves. Among the many advantages of the THz waves include the fact that they pose little to no harmful effects to biological tissues (i.e., they can also be used to inspect tissues in biomedical engineering) [11,12], and they can provide high spatial resolution due to their relatively shorter wavelengths. All these remarkable properties confer THz waves the capabilities to detect and evaluate various types of defects in GFRP composites [13][14][15], measure their optical and dielectric properties [8,16], and the fiber orientation in unidirectional GFRP composites [17] and the stacking sequence of angled GFRP composites [18], etc. THzbased NDT applications have been delayed for many years due to some technical issues that limited the development of efficient THz emission and detection devices, and the THz region of the electromagnetic spectrum was formerly known as the "THz gap". However, this gap has already been lifted following the fast development of highly performing semiconductors, ultrafast electronics, and laser systems, as well as ultra-micro machining technologies. Indeed, these fast developments have equally contributed to the development of highly efficient THz time-domain spectroscopy (THz-TDS) that is now widely used in NDT applications in many different fields [2,9,19]. In all these applications, however, the configurations and sample scanning processes are often presented differently and Table 1 summarizes the most common configurations and sample scanning processes available in the literature today.
There are currently many different studies outlining the use of THz-TDS for the identification, visualization, and characterization of defects and structural damage in fiber-reinforced composite structures, and the fact that THz-based inspections do not require any couplant between the testing device and the material system under test (i.e., non-contact method) is another added advantage. In Ref. [20], for example, the authors used the THz-TDS to detect and characterize simulated delamination defects in a 12 layers carbon fiber-reinforced polymermatrix (CFRP) composite plate and reported having achieved successful test and accurate results. Although this study was successful and the test results were accurate, several authors have indicated that the conductivity of carbon fibers generally hinders the penetration of THz waves in CFRP composites, and THz waves are only capable of detecting the impact-induced matrix cracking on the surface and subsurface of CFRP composite materials [21][22][23]. In Ref. [13], the authors utilized a PCA-based THz-TDS system configured in reflection and transmission modes to visualize a delamination defect simulated in a GFRP composite plate consisting of eight glass fiber layers and reported successful results. In a similar study [24], a fiber-coupled THz-TDS based on PCA was used to visualize multi-delamination defects and their thickness in a GFRP composite plate. However, the study used averaging method to remove the noise from the THz signals which drastically reduced the amplitude of the measured THz signals making it difficult to determine the absorbance of the THz wave at each interface within the GFRP composite sample.
Additional studies featuring the application of the THz-TDS system for the NDT of GFRP composites include but are not limited to Refs. [25,26]. In Ref. [25], for example, the authors used a fiber-coupled pulsed THz-TDS system featuring PCA-based THz waves emitter and receiver and they were able to extract occluded textual content from a packed stack of paper pages down to nine pages without human supervision. In a recent study [26], the authors used a PCA-based THz-TDS system to detect and evaluate simulated delamination defects in a GFRP plate with the THz wave operating at an incident angle of 25° in the reflection mode and reported successful results. However, all these studies used the averaging method to remove the noise from the measured THz signals which, in some cases, resulted in the amplitude of the measured THz signals being reduced, which subsequently limited the accuracy of the obtained results. In particular, the use of averaging method to remove the noise from the measured THz signals may not provide accurate results when the THz-TDS is used to measure the thickness of internal delamination defect/damage in GFRP composites because the accuracy of these types of measurements highly depends on accurate temporal localization of the different peaks in the measured signals.
In this study, we use a fiber-coupled PCA-based THz-TDS system to detect and evaluate the thicknesses of hidden multi-delamination defects in GFRP composite laminates. The optical parameters of the different materials involved in the sample makeup (i.e., compact GFRP composite and polytetrafluoroethylene, shorten to PTFE, for simulated delamination defects) are first determined. Then, the thickness and location of each delamination defect in the z-direction are successfully determined using analytical and experimental approaches. As opposed to previous studies using the THz-TDS system to detect GFRP composite materials' delamination defects and the averaging signal processing method to remove the noise from the measured THz signals [13,20], the present study applies the wavelet shrinkage de-noising algorithm to remove the noise and improve the signal-to-noise ratio (SNR) and preserve the amplitude of the measured TOF signals which enables accurate determination of the thickness and location of each delamination defect. The measured and calculated results are compared to determine the accuracy of the proposed method and the effect of the wavelet shrinkage de-noising algorithm on the evaluation of GFRP composite materials' delamination defects. Also, the attenuation of the THz waves at every interface in the through-the-thickness direction is calculated using Fresnel equations and compared with the experimental values. This information is important as it can help NDT engineers and practitioners to determine the effective thickness of the GFRP composite parts that can be evaluated by the THz-TDS system when considering different options or selecting adequate NDT tools for their NDT tests. The rest of the paper is organized as follows: Section 2 provides the theoretical formulation and THz signals analysis framework to determine the materials' optical parameters as well as relevant signal processing techniques used to process all the measured THz signals. Section 3 outlines the experimental methods and system description. Section 4 provides the results of the study and the discussion to provide adequate explanations of all the results obtained and Section 5 outlines the conclusions of the study.

Theoretical Formulation and Analysis of THz Signals
In this section of the study, the theoretical background and the mathematical formulations used to evaluate the material's optical properties and the thickness of the delamination defects are presented. The section also details the mathematical formulations pertinent to the development of the stationary wavelet transform (SWT) based algorithm that is used to enhance the measured TOF signals and obtain accurate test results.
To be able to fully characterize the sample and evaluate the position and sizes of the different delamination defects, parameters such as the thickness of each delamination, their location in the sample, as well as the actual attenuation level of the THz radiation/waves in the GFRP composite sample are determined.

Extraction of the Sample's Optical Parameters
The optical parameters of the GFRP composite samples and the Teflon film used to simulate the delamination defects in the GFRP composite samples are determined using the reference and sample signals measured in reflection and transmission modes. The reference and sample time-domain signals are captured by placing the gold-painted mirror and the GFRP composite sample at the sample holder platform on the THz-TDS system, respectively. These time-domain signals are captured in both the transmission and reflection modes and their corresponding frequency-domain signals are calculated by applying the Fast Fourier transformation (FFT) using the following equation: where E(t) and E(ω) denote the THz time-domain signal and its corresponding frequency domain signal, respectively. After the transformation of the THz time-domain signals for both the samples E s (t) and the reference E r (t) into their frequency-domain counterparts, the optical properties of the GFRP composite samples, and the Teflon film used in the study were calculated. The optical parameters of the GFRP composite materials represent their macroscopic characteristics at THz frequencies.
The representation of these optical parameters is generally subsumed in the following equation of the material's complex refractive index: where n(ω) denotes the real part of the refractive index and κ(ω) the extinction coefficient of the material's complex refractive index. The material's extinction coefficient κ(ω) is directly proportional to the absorption coefficient expressed as α(ω) = 2ωκ(ω)/c with ω being the angular frequency, and c the in-vacuo propagation speed of light (i.e., THz wave in this case) [13]. In the reflection mode, a difference in the optical path length occurs based on the different propagation velocities of the THz waves in the sample and air. Considering a sample of thickness d , the time delay spent by the THz wave passing through the GFRP composite sample from the time it impinges the bottom surface to the time it reaches the rear surface of the same sample can be expressed as: where the time delay t 1 and t 2 correspond to the peak time delays of the reflections at the bottom and rear surfaces of the GFRP composite sample, respectively. Therefore, by transforming Eq. (3) it is now easy to obtain the refractive index of the GFRP composite sample.
At this stage, it is observed that the time-resolved detection scheme of THz-TDS can also be directly applied to measure the depth characteristics or features to obtain the information describing the multilayer structure and determine the in-depth localization of the multi-delamination defects of our GFRP composite samples using the following expression: n(ω) = �tc 2d . where d represents the thickness of the sample or the inclusion in the sample, t the time elapsed between two consecutive peaks or two different reflections at the times t i and t i+1 , θ denotes the incident angle of the THz waves impinging the surface of the test GFRP composite, and n the refractive index of the sample or the inclusion in the sample [13]. It is noted that the value of the thickness d may also be the localization of the defect or the damage in the through-the-thickness direction (i.e., delamination, inclusions, etc.) or the depth location of an interface of a different refractive index. After measuring THz timedomain signals for the sample E s (t) and reference E r (t) and calculating their corresponding frequency-domain signals E s (ω) and E s (ω) , respectively, the absorption coefficients for the GFRP composite samples and the Teflon inclusions were calculated as follows: Subsequently, the transmission function H (ω) is also calculated using the following expression: where β = ω ñ(ω) − 1 c d is a complex expression that represents the propagation coefficients of the electromagnetic wave propagating through the GFRP composite sample at an angular frequency ω over the propagation distance d that is equal to the thickness of the GFRP composite sample. The parameters E r (ω) and E s (ω) are the reference signal and the sample signals, and the reference signals, respectively, ϕ(ω) is the phase difference between the reference and the sample, and ρ(ω) is the magnitude of the ratio of the sample to the reference THz signals. Using the expression of the transmission function given in Eq. (7), the equations for the optical parameters of the GFRP composite sample are obtained as follows [16]:

The Determination of the Thickness of the Sample and Defect
The time-resolved detection scheme of THz-TDS is directly applicable to measuring the depth information of a multi-layer sample. When pulsed THz waves are incident on an object, the reflected THz waveform consists of a series of pulses reflected from the interfaces [21]. The thickness of the sample or defect can be calculated in the THz reflection scan mode via the following equation: where T GFRP represents the thickness of the sample; t is the time between successive reflections; c is the speed of light in air; n GFRP is the refractive index of the sample, and θ is the incident angle of the THz waves. The factor of one-half arises since the THz waves are measured in the reflection mode. In this study, the angle θ is assumed to be 0 • because the incident THz wave is quasinormal and the expression in Eq. (11) will now change to T GFRRP = �tc/2n GFRP without the factor cos (θ) since cos (0) = 1.

The Fraction of Incident THz Power
The fraction of the incident THz power that is reflected from the interface is given by the reflectance. To evaluate the path of the THz waves, the reflection and transmission power of THz waves at the interface of two media (i.e., materials 1 and 2) can be calculated using the Fresnel equations as follows [24]: where R s represents the reflectance of s-polarized THz waves; R p is the reflectance of p-polarized THz waves; z 1 and z 2 are the wave impedances of the media involved (i.e., materials 1 and 2); µ 1 and µ 2 are the magnetic permeability of materials 1 and 2, respectively; ε 1 and ε 2 are the electric permittivities of the materials 1 and 2, respectively, at the frequency of the THz waves, θ i is the incident THz-wave angle, θ t is the reflection THz-wave angle. In the case of non-magnetic media (i.e., materials for which µ 2 = µ 1 = µ 0 (where µ 0 denotes the (11) permeability of free space), impedances z 1 and z 2 can be expressed as z 1 = z 0 /n 1 and z 2 = z 0 n 2 , respectively. Replacing the expressions of these two parameters into Eqs. (12) and (13) and simplifying the resulting expressions, the following equations are obtained for the s and p polarized reflectances of the THz waves as follows [24]: The total reflectance R of the THz waves in the GFRP composite samples can now be calculated as: In this case, it follows that the total transmittance ( T ) of the THz waves can be expressed as: The THz-wave power can be calculated using Eq. (17), which corresponds to the transmission power of the THz-wave at the surface of the GFRP composites.

The Absorbance of THz Waves in Refraction Material
In refraction materials such as GFRP composites, the absorbance of the THz waves denotes the attenuation of the transmitted THz wave power in the material at the different interfaces. In this context, the THz waves impinging the top surface of the GFRP composite sample propagate through the sample until they reach the interface between the GFRP and the upper side of the top delamination defect. At this level, the incident THz beam is partly reflected at the upper side of the delamination defect while the remaining portion of the THz beam is transmitted into the remaining sections of the GFRP composite sample. The composition of our GFRP composite samples is such that the next section of the sample adjacent to the interface of the first reflection is the delamination defect. In this case, the remaining THz waves after that first reflection are transmitted into the delamination defect and the aforementioned process is repeated at every interface until the THz waves reach the rear surface of the sample. However, when the THz waves pass through the specimen, transmissions and reflections are repeated at every interface between two different media for every reflected portion of the THz, and only the first reflection and/or transmission at every interface reaches the detector to be measured. The transmitted/ reflected THz waves that cannot reach the detector or that cannot be easily visualized are generally referred to as absorbed THz power. The absorbed power of the THz waves can be calculated using Eq. (18) as follows [26]: where A denotes the absorbance of the THz waves in the GFRP composite material; I denotes the incident THzwave power; α denotes the absorption coefficient of the GFRP composite material or simply the absorption coefficient of our GFRP composite samples, and d is the thickness of the GFRP composite sample. This equation is only applicable to a single set of composite materials and cannot be used when we have several materials superposed one on top of the other (i.e., the indicated parameters are for one type of material, meaning that the Teflon inclusions are not included in the equation and should be calculated separately). As indicated in the previous paragraph, this scenario will repeat itself at every interface between the delamination defect and the compact section of the GFRP composite sample (i.e., media of different refractive indices) where only the effects of two parts of E in,i , namely E in,i+1 and R i+1 are visualized and quantified and the absorbed THz power at each of these interfaces can be calculated using Eq. (18).

The Stationary Wavelet Denoising Algorithm
To obtain good quality signals and accurate measurement results, all measured signals are subjected to a stationary wavelet transform (SWT) based signal denoising algorithm for denoising the THz signals before using them to calculate the samples' features of interest. In its detailed implementation, the SWT is the translation-invariance Unlike the general DWT, the SWT is always up-sampled (i.e., never sub-sampled) at each level of decomposition. Inspired by the close similarity between the THz pulse and the common wavelet basic functions, its application in THz signal processing practices was first presented in Ref. [28] whereby the authors used the SWT to extract the tomographic results (i.e., the spectroscopic information about each reflecting layer of a sample). In most signal processing practices where the use of SWT is involved, the latter decomposes a single-dimensional signal x(n) into approximation coefficients vector cA k,l and detail coefficients cD k,l by convolving with a lowpass filter ψ i (n) and a high-pass filter φ i (n) (i.e., where i = 1, 2, 3...p and p denotes the maximum level of decomposition of the signal by the SWT) along the temporal axis [29]. It is noted that the low-frequency information is concentrated in the approximation coefficients while the high-frequency information is concentrated in the detail coefficients. Figure 2 illustrates the decomposition process of the SWT whereby the single-dimensional signal passes through a filter bank for the convolution process.
It is noted that the wavelet coefficients with small absolute values are generally considered the noise while those with large absolute values are generally regarded as the main featured information of the signal. Applying the thresholding process will remove the small absolute value coefficients, and depending on the decomposition level and the type of wavelet transform used, the reconstruction of the signal from these coefficients is projected to generate a signal in which the contribution of noise has been greatly reduced [30]. In general, SWT de-noising with soft thresholding is performed through 4 main steps [29,31]. The wavelet coefficients are determined by taking the SWT of the input signal: Assuming the noise level in the THz time-domain signal is equal to the parameter σ , and n the number of sampling points in the THz time-domain signals, the value of the threshold parameter T is calculated using the following equation: In this particular study, the values of the wavelet coefficients cD k,l are thresholded by using the soft-thresholding process (i.e., which is one of the thresholding processes for the SWT operations) and the following results are obtained: At this stage, the inverse stationary wavelet transform (ISWT) is performed to obtain the denoised THz timedomain signal x(n) as follows: In the present study, we only considered the detailed coefficients to ensure all the information in the signal is captured and obtain highly accurate results.

Experimental Methods and System Description
In this study, the THz-TDS is used to evaluate and characterize a series of simulated delamination defects in layered GFRP composite samples. The THz-TDS system uses THz radiation for the spectroscopic characterization and imaging of the GFRP composite samples which are generally generated and detected using methods such as optical rectification, photoconductive antenna (PCA), and the surface electric field of a semiconductor [32], which may be done coherently or incoherently depending on whether both the amplitude and phase of the field are measured or only the intensity of the THz radiation is measured [24]. In general, coherent detection is closely associated with the generation technique for the THz waves as they share the underlying mechanisms and key components (i.e., utilize the same light source). In this study, PCA is used to coherently generate and detect the THz waves and the following section provides a detailed description of the operation of our THz-TDS system.

THz-TDS and Imaging System
The THz-TDS (TeraView TPS 4000) is used to detect and characterize hidden delamination defects in GFRP composite samples. To measure the optical parameters and characterize the delamination defects of the GFRP samples, the system is configured to perform tests in reflection and transmission modes. These two modes are not simultaneously configured but rather one is configured after the other depending on which mode is needed at the different stages of the study. This THz-TDS system features a scanning range of up to 1200 ps and its resolution can reach 0.1 ps at a rapid scanning frequency of 50 Hz. In its effective range of frequencies (0.01-4.5 THz), the system's SNR is generally more than 90 dB at the peak and this remains at more than 70 dB during the entire sample measurement process. The system's main components include the femtosecond laser, the optical delay line, and the X-Y imaging stage as well as the two photoconductive antennas (PCA), i.e., the emitter and receiver, (21) and several optical components (Figure 3). The sample housing unit is filled with Nitrogen (N 2 ) during the measurement to eliminate the effect of moisture/humidity on the test results and the power of the THz radiation is kept below 1 mW to avoid any potential thermal strain on the sample.
The operation of our THz-TDS is such that the PCAbased THz emitter gets excited by an ultrafast femtosecond laser and produces single-cycle THz pulses that are coherently detected by the PCA-based THz detector in reflection or transmission mode. The sample is raster-scanned by a set of motorized stages moving in an X-Y direction and the generated THz pulses intercept the sample along the vertical plane and get reflected by the same sample (reflection mode) or transmit through the sample (transmission mode) before it is collected and sent to the laser-gated PCA-based detector. A current proportional to the THz electric field of the THz radiation is produced and measured (i.e., amplitude and phase) by gating the photoconductive gap to the femtosecond pulse that is synchronized with the emission of the THz radiation. After amplification and the processing of the obtained THz signal at each pixel, a final THz image is obtained. The existence of distinct refractive indexes in the material (i.e., compact GFRP composite sample and the TPFE discs/films) produces multiple reflections at the different interfaces and these reflections are presented in form of several peaks at the corresponding time delays in the transmitted and/or reflected THz signals. The sequence and thicknesses of the various laminates of the test sample can be reconstructed easily by analyzing these peaks with knowledge of the refractive indices of the material constituents. The absorbance of the THz waves by the sample is also calculated by analyzing the difference between the incident power of the THz waves and the power at every interface in the GFRP composite sample. Finally, the THz time-domain signal can be extended to its frequency-domain version using the FFT for additional analysis.

Sample Preparation Procedure and Materials
The GFRP composite samples used in this study were manufactured using unidirectional (UD) glass fiber prepreg (G15000/6509/33%, GF contents: 67 wt%, Weihai GuangWei Composites Co. Ltd, Shandong-China). All the samples were laminated by the hand layup process following a unidirectional pattern. A series of Teflon film disks of different diameter sizes (i.e., 5 mm, 10 mm, and 15 mm) were embedded in each sample to simulate the location of three different levels of delamination in between the stacked layers of glass fiber prepregs during the manufacturing process ( Figure 4). Each disk had a thickness of 130 µm and their centers were thoroughly aligned. Prior to deciding the thickness of the different delamination defects, we first conducted an overview of the literature to know the usual thickness of simulated delamination in GFRP composite materials. Our literature survey indicates that most of the simulated delamination defects present an estimated thickness between 20 µm to 250 µm. In this context, we decided to fabricate the simulated embedded delamination defects using Teflon film with a thickness of 130 µm. The thicknesses of the different simulated delamination defects found in the literature are also listed in Table 1 and readers are directed to this part of the work for more information.
A total of 13 layers of unidirectional glass fiber prepreg with an estimated thickness of 120 µm each were used to manufacture a 220 mm × 220 mm GFRP composite plate. After the laminating process by hand-layup, the GFRP composite plate was cured in a vacuum oven for 2 hours at 125 °C. The pressure in the vacuum oven was maintained at 0.6 MPa during the entire time of the curing process. After the curing process, the sample was left to cool down and a micrometer (DL321025S, Deli Group Co. Ltd.-China, range: 0-25 ± 0.001 mm) was used to measure its thickness.
The thickness of the test sample was measured at 8 different locations and the standard error was calculated. The sample thickness was found to be equal to 1.532 ± 0.045 mm. It is noted that no delamination was visible at the edge of the sample nor was the delamination visually apparent or detectible by thickness variation measured by the micrometer. After measuring the thickness, the 220 mm × 220 mm × 1.532 mm GFRP composite plate was cut into 25 different pieces of 40 mm × 40 mm × 1.532 mm each to obtain the samples used in our experiments. It is noted that the PTFE or Teflon disks were randomly placed in between the different layers of the composite plate but each time, the disks were carefully arranged to form a three-level  multi-delamination system. To obtain the best samples, a preliminary scanning of all 25 samples obtained by cutting the manufactured GFRP composite plate was conducted using the THz-TDS system. This preliminary scan revealed that only 5 samples had 3 layers of delamination defects of different sizes at different depth locations each as we intended to use them (Figure 4(b)). As such, only these 5 samples were considered for further tests, and the rest of them were discarded. The reference GFRP composite sample (i.e., the GFRP composite sample manufactured following the aforementioned process but does not have any inclusions) and a piece of Teflon film are used to determine the optical properties of these two material systems prior to their utilization in our study.

Optical Parameters of the Samples
The optical parameters of the samples are first determined as the starting point by using the time-domain THz signals of the reference and sample measured in reflection and transmission modes. As such, both the reference and sample THz time-domain signals, i.e., E r (t) and E s (t) , respectively, are first captured by placing the gold-painted mirror and the GFRP composite sample at the sample holder unit (also known as the X-Y stage of the THz-TDS), respectively. The same procedure is used for both the GFRP composite reference sample and Teflon film whereby the reference sample is replaced by a piece of PTFE (i.e., results not shown here for simplicity and brevity of the article). To remove the influence of the noise and ensure the high-quality timedomain signals are used in our calculations, we first applied the wavelet shrinkage de-noising to process the THz time-domain signals before the transformation in the frequency-domain. To do this, we choose the Symlet (Sym4) type of wavelets, which are a modified form of Daubechies wavelets. A maximum level of 9 is used for the wavelet decomposition because there were no noticeable changes observed when using higher levels of wavelet decomposition to justify the extra computational expenses. This level of decomposition is maintained for all our wavelet shrinkage denoising operations used in this study for consistency purposes. The corresponding THz signals in the frequency domain, i.e., E r (ω) and E s (ω) , are calculated from the aforementioned THz time-domain signals by using Eq. (1). In the time domain, the noise is predominantly observed in the region with greater fluctuations just after the main THz pulse of the measured signals before the denoising process. A typical example of the time-domain THz sample signal and its frequency-domain spectrum before and after the wavelet de-noising operation are shown in Figure 5. It is observed that the region with larger fluctuations/noise is located in the temporal region between 6.5 and 33 ps. Figure 5 also indicates that the use of SWT helps remove fluctuations in the signal without losing much energy in the main THz pulse in the time domain while major noise dips are successfully suppressed in the frequency domain.
The aforementioned process was also used to denoise the signals and determine the FFT for the different signals used to evaluate the different features of the GFRP composite samples. Representative signals with a significant reduction of noise for different regions of the sample with and without delamination defects were obtained after the wavelet shrinkage denoising process. After the calculation of the results in Eqs. (8)- (10), the values of refractive indices for both the reference GFRP composite laminate and Teflon film were obtained as presented in The average value of the calculated real part of the refractive index n of the GFRP composite sample and the PTFE film is found to be 2.157 and 1.462, respectively, at the frequency of 0.45 THz. Figure 6 indicates that the refractive indices for both materials remain almost steady (i.e., they vary very little in the considered range of frequencies) while the absorption coefficients of these two materials (i.e., the GFRP composite sample and the Teflon film) continuously increase in the same range of frequencies. Interestingly, the values of the refractive indices obtained for both materials are also consistent with those reported in Refs. [33][34][35] for the GFRP composites and the Teflon films, respectively. Figure 7 indicates that when the THz beam ( E in ) goes through a multilayered GFRP composite sample with multi-delamination defects (i.e., materials with different layers where the individual layers are composed of different material systems), it will be partly reflected ( R i ) and partly transmitted ( E in,i ) at each interface encountered. This is because the different interfaces present different refractive indices which cause the refraction angles of the THz wave propagation to change. In this context, the THz wave will also propagate through the GFRP composite sample passing through the top delamination defect, the middle delamination defect, and the bottom delamination defect (i.e., these delamination defects are considered the different layers) until it reaches the back surface or rear of the GFRP composite sample as illustrated in Figure 7. The signals reflected at each interface ( R i with i = 1, 2, . . . , 8 ) will travel back through the multilayer structure to the PCA detector, while the transmitted signal will continue through successive interfaces ( E in,i with i = 1, 2, . . . , 8 ) until it reaches the interface between the back surface of the composite and the air (Figure 7). The reflected signals will appear in the THz time-domain waveform at different times of flight according to the sequence of mediums penetrated by the THz wave. This timedomain waveform will be used to detect the different delamination defects and to estimate their relative indepth location based on the TOF between the different interfaces.

Thickness Characterization of the Internal Delamination Defects
To visualize the position of the delamination defects in the GFRP composite sample, the wavelet shrinkage denoising was also applied to process the THz time-domain signals (i.e., the sample signals) in both the regions with and without delamination. To do this, we chose the Symlet (Sym4) type of wavelets, which are a modified form of the Daubechies wavelets. Likewise, a maximum level of 9 was used for the wavelet decomposition because there were still no noticeable changes that could be observed when using higher levels to justify the extra computational expenses. After the wavelet shrinkage denoising process, we obtained representative signals with a significant reduction of noise for different regions of the sample with and without delamination defects. Figure 8 presents a comparative representation of the different THz timedomain waveforms collected from the regions of the GFRP composite samples with and without delamination defects (i.e., see the individual graphs for the different labels). In Figure 8(a), the data obtained from the delamination-free or pristine region of the GFRP composite sample is presented for both clarity and comparison purposes. This signal has two peaks, the first peak appears at t 1 = 6.271892 ps and corresponds to the first interface which is also the top surface of the GFRP composite sample, while the second peak appears at t 2 = 28.33535 ps and corresponds to the bottom surface or the rear of the GFRP composite sample. The TOF between these two different peaks was found to be 21.663458 ps. This value of TOF between these two interfaces of the pristine region of the GFRP composite sample was particularly important because it was used later in our calculations to determine the thickness of the GFRP composite sample which was then compared with the real thickness of the GFRP composite sample for the validation of the method.
In Figure 8(b), the 21.44864-23.5741 ps interval in the THz signal highlighted in cyan represents the TOF covered by the THz pulse while passing through the bottom delamination defect (i.e., the delamination defect closest to the rear of the GFRP composite sample). Figure 8(b) also shows that the bottom delamination defect is the biggest in diameter as it covers the whole area of the other two delamination defects (middle and top delamination defects). However, the fraction of the THz beam reaching this point is small, and that is why the amplitude of the THz signal coming out of this delamination defect is also smaller than the other two. In Figure 8(c), the 12.68007-14.40954 ps interval in the THz signal highlighted in green equally represents the TOF covered by the THz pulse while passing through the middle delamination defect. At this position, the middle delamination overlaps with the bottom delamination defect because the latter is smaller in diameter than the former and it is not possible to select any region within the middle delamination without selecting the top delamination at the same time. In Figure 8(d), the 7.976565-9.609055 ps interval in the section of the waveform highlighted in blue represents the characteristics of the top delamination defect. The selection of the region with the top delamination entails that we also unintentionally select the regions of the other two delamination defects because the top delamination defect is smaller in diameter (9 mm) than the other two (12 and 15 mm for the middle and bottom delamination defects). In this context, Figure 8(d) represents the data from the delaminated region (i.e., the region on the GFRP composite sample with multiple delamination defects). In summary, the various subsidiary peaks observed at the optical delay times greater than 10 ps correspond to the different pulses that have bounced back off various interfaces between the compact GFRP composite and its delaminated regions. That is, Figure 8 and the reason we can see that is because of the relatively large refractive index mismatch between the resin matrix and the PTFE at that level. In the images indicated in Figure 8(b)-(d), the PTFE film inserted into the GFRP composite to create artificial delamination defects caused the composite's laminas at the three different levels of the delamination defects to deform and the corresponding peaks in the time-domain THz signals are shifted to larger optical delays. In these two signals, the features outside the delamination defects and in between the individual delamination defects are similar. However, the plies in the regions between the top and bottom delamination defects are substantially compressed due to the fabrication process that ensured a uniform overall thickness of the entire GFRP composite laminate even after adding the 3 levels of the PTFE films. At this stage, a detailed investigation into the measured data was conducted to obtain quantitative information describing the internal features of the GFRP composite sample and delamination thicknesses. The aforementioned distinguishable reflected signals characterized by their optical delays in conjunction with the refractive indices of the constituent materials as calculated in the previous section of the present paper enabled us to determine the exact location and thickness of the individual delamination defects within GFRP composite samples. In particular, a signal denoising process based on wavelet shrinkage was performed which helped us to eliminate all the cases of overlapping echoes and unclear peaks in the THz time-domain signals before calculating the quantitative information describing the internal features of the composite and delamination thicknesses. Using the calculated refractive indices of the GFRP composite sample and the PTFE film (Figure 7) as well as the information extracted from the wavelet denoised signals, we were able to calculate the thickness of the portions of the GFRP composite sample situated between every two consecutive reflections be it for the delamination or the GFRP composite sample itself. In this context, it was also possible to determine the exact location of the individual delamination defects within the GFRP composite sample. As indicated earlier, the thickness of the entire GFRP composite sample is determined by first identifying the peaks of the time delays between the front and the rear surfaces of the GFRP composite and substituting the required values of the different parameters in Eq. (11) for subsequent calculations. However, since our THz system uses quasi-normal incident THz waves, the angle of incidence θ is neglected and assumed to be equal to 0° in all the calculations in this study, suggesting that cos θ = 1 . As such, the overall thickness of the composite is calculated as follows: where c = 2.998 × 10 8 m/s is the in-vacuo speed of light, ∆t=22.604314 ps is the time between the reflections at the bottom surface and the rear of the GFRP composite, and the factor of one-half arises since the system is used in reflection mode (i.e., the THz pulse passes through the composite twice). The calculated thickness of the GFRP composite is equal to 1.571 mm, which is close to the measured value of 1.532 ± 0.045 mm. The difference between the thicknesses of the GFRP composite sample measured by the micrometer and the THz-TDS system originates from the difference in refractive indices in the GFRP composite. In this calculation, we did not consider the refractive index of the delamination defects as being different from that of the compact GFRP composite. Instead, we considered the entire composite to be a compact unit whilst it has some delamination defects that have a different refractive index. Indeed, the best practice would have been to consider the different TOF and the corresponding refractive index and add up all the different thicknesses to obtain the correct thickness of the GFRP composite (i.e., individual thicknesses of the delamination defects and the inter-delamination regions should be considered separately). Considering the different refractive indices, the newly obtained thickness of our GFRP composite is 1.536 mm, which is the same value as the thickness we measured using the micrometer with the standard error considered. Similarly, using the calculated refractive index for the PTFE film and the time delays between the reflections associated with the interfaces between the individual PTFE films and the adjoining composite layers, we can easily determine the thickness of the individual layers of PTFE films in the delaminated GFRP composite samples. In the case of the top delamination defect, for example, the thickness is calculated as follows: where ∆t = 1.272837 ps and n Teflon = 1.462 corresponds to the time delay of the reflections at the interfaces between the front and the rear sides of the PTFE film and the refractive index of the Teflon film, respectively. Again, the calculated value of the Teflon film is quite close to the measured value of 130 µm. At this stage, we use the same formula to calculate the thickness of the individual delamination defects in the GFRP composite sample. Table 2 presents the calculated results of the individual delamination defects. As indicated in Section 2, all the PTFE discs inserted in the composite sample had the same thickness of 130 µm each. The only difference is their diameter sizes which were 9, 12, and 15 mm for the top, middle, and bottom delamination defects, respectively. It is observed that the experimental values of the thicknesses of the individual PTFE discs strongly agreed with the measured values, further confirming the accuracy of our THz-TDS measurements.
To demonstrate the benefits of using the SWT to process the THz time-domain signals, we used the THz signal as measured directly from our THz-TDS system and   selected several regions on a single GFRP composite sample, and calculated the average thicknesses of the individual delamination defects before and after the denoising process by the SWT. Table 3 presents the results of the individual delamination defects calculated from the THz time-domain signals before and after the signal denoising by the SWT-based algorithm. The values of the thicknesses of the individual delamination defects calculated from the measured THz signals indicate that the measured and calculated values are very close and in some cases equal to the real values, which further confirms the accuracy of our measurements. Using the same logic, we can also calculate the exact position of the different delamination defects in the through-the-thickness direction of the GFRP composite samples. In our earlier discussion, we indicated that the different peaks in the time-domain signals indicate the different interfaces in the GFRP composite samples. To calculate the exact in-depth location of the top delamination defect, we evaluated the thickness of the composite between the first and second peaks in the time-domain signal (i.e., the thickness of the composite before the top delamination). The THz pulse considered in this case indicates that the first peak is located at the time t 1 = 6.279424 ps while the second peak is located at t 2 = 8.978687 ps.
where c and T loc−td denote the in-vacuo speed of light and the in-depth location of the top delamination defect, respectively. The parameter t denotes the TOF between the front surface of the composite and the first interface of the top delamination defect and is equal to 2.699263 ps . Using the same equation, the middle and the bottom delamination defects are found to be situated at exactly 530.74 µm and 1.081 mm from the front surface of the GFRP composite, respectively. It is noted that the locations of the middle and the bottom delamination defects are calculated in different steps following the passage of the THz pulses as they go through the layers of different refractive indices in the GFRP composite, viz. the refractive index of the inter-delamination layers or compact GFRP composite ( n GFRP ) and the refractive index of the delamination defects ( n Teflon ). In a similar context, we can now determine the location of all the delamination defects in 5 different GFRP composite samples. Table 4 reports the values of the distances or thicknesses indicating the locations of the individual delamination in the through-thickness direction of the GFRP composite sample considered. It is indicated that these delamination defects were created at the same depth location when

Absorbance and Fraction of the THz Waves
To investigate the optical path of the THz wave passing through our multilayered GFRP composite sample with multiple delamination defects, the power of the reflected THz wave was calculated using the Fresnel equations as presented in Eqs. (14)- (17). It is reminded that the THz-TDS was used in both the transmission and reflection modes to obtain the THz signals used to measure the optical parameters of the materials and the reflection mode to obtain the signals used for the nondestructive evaluation of the delamination defects. In all these cases, the incident THz wave impinging the surface of the GFRP composite sample was considered quasi-normal (i.e., the incident angle is equal to 0°) so the expression for sin (θ) = 0 . Using the aforementioned Fresnel equations, we found that the THz wave was transmitted into the GFRP composite sample with 72.61% of the waves transmitted through the surface, while 27.39% of the power of the THz wave was reflected. This reflected THz component was captured by the PCA receiver as the power of the reflected THz waves as represented in Figure 9. The intensity of the first peak of the THz pulse from the PCA emitter is represented by the signal from the gold-painted mirror, while the first peak of the THz wave reflected from the surface of the GFRP composite sample is represented by the GFRP composite plate signal.
In analyzing the absorbance of the THz wave by the GFRP composite plates and the different interfaces we will use the information outlined in Figure 7. The latter indicates that there are 8 reflections in total in line with the number of delamination defects present in the GFRP composite sample. Out of these eight reflections, seven (i.e., from E in,1 to E in,7 ) are subjected to the THz  [36]. The power of the THz wave measured at this level is estimated to be 3.824% for the third peak, which is about a − 6.54% maximum error. This error might have been caused by the multiple calculations performed to reach this stage.
In the case of the middle delamination defect, several transmissions and reflections must be taken into consideration. In primis, the first, second, and third peaks were obtained from the surface and the first delamination, respectively, and the calculation details are provided in the previous paragraph. In the same context, the fourth and fifth peaks originated from the middle delamination defect. The second and fourth peaks share the same phase because they are both reflected from the upper side of the top and middle delamination defects, respectively. However, the fifth peak presents a reverse phase with the fourth peak because the former is reflected from the fixed end while the latter is reflected from the free end, corresponding to the bottom and upper sides of the second delamination, respectively. One of the main advantages of the THz waves is that they can successfully penetrate overlapping delamination defects to allow users to easily detect deeply buried defects beneath other defects/damage. Using the same calculation patterns, we can also calculate the powers of the THz wave at the interfaces associated with the middle delamination defects. In this context, the power of the THz wave penetrating the top delamination and reaching its bottom side is estimated to be 28.003% (corresponding to 0.7261 3 × 0.7315). Considering the attenuation of the THz wave at the upper side of the middle delamination defect, the power of the transmitted THz wave penetrating the GFRP was determined to be 20.048% (corresponding to 0.7261 3 × 0.7315 2 ). In this case, the power of the transmitted THz wave when being reflected at the upper side of the middle delamination was estimated to be 5.611% (corresponding to 0.7261 3 × 0.7315 2 × 0.2739). The reflected THz wave was transmitted again into the compact GFRP, and its power was determined to be In this context, the intensity of the THz power in the fourth peak was found to be 1.149% (corresponding to 0.7266 6 × 0.7315 4 × 0.2739) with three times the transmittance and one time the attenuation. Similarly, the power of the THz wave in the fifth peak was estimated to be 0.606% (0.7261 8 × 0.7315 4 × 0. 2739), with two times the transmittance when compared to that of the fourth peak. The power of the THz wave measured at this level is estimated to be 0.998 and 0.532%, for the fourth and fifth peaks, respectively. As indicated in Figure 7, the centers of all three delamination defects overlapped. The reflections from R 1 to R 5 have already been calculated and the power of the THz wave at these interfaces has already been calculated. At this stage, we can now calculate the power of the THz wave for the sixth and seventh reflections. These two reflections/peaks originated from the bottom delamination defect. Since the THz wave can easily penetrate the compact GFRP composite sample and the three overlapping delamination defects, it is easy to evaluate the sizes of all three delamination defects. Similarly, the sixth and the second peaks had the same phase because they are both reflected from the upper side of the bottom delamination defect. However, the seventh peak presents a reverse phase because it was reflected from the bottom side of the bottom delamination defect. The THz wave went through the compact GFRP composite sample, in the same way, it did for the middle delamination defect. At the bottom side of the middle delamination defect, the power of the transmitted THz wave was found to be 10.800% (corresponding to 0.7261 5 × 0.7315 2 ). Subsequently, the THz wave penetrated the compact GFRP composite with a power of 8% (corresponding to 0.7261 5 × 0.7315 3 ) following the attenuation in the compact section of the GFRP composite sample. Similarly, the transmitted THz wave with attenuation was reflected at the upper side of the bottom delamination defect, and the power at this level was estimated to be 2.164% (corresponding to 0.7261 5 × 0.7315 3 × 0.2739). The power of the sixth peak was found to be 0.171% (corresponding to 0.7261 10 × 0.7315 6 × 0.2739). This power is calculated using the same process as above considering the transmittance and attenuation, five times and three times, respectively. The power of the seventh peak was found to be 0.090% (corresponding to 0.7261 12 × 0.7315 6 × 0.2739) considering the two times transmittance and compared with that of the sixth peak above.

Conclusions
In this work, the THz-TDS system is used to characterize multilayered GFRP composite laminates and evaluate the thickness and location of simulated internal delamination defects. The study also presents a method to calculate the power of the THz waves at every interface in the through-the-thickness direction to help NDT&E engineers and practitioners estimate the maximum thickness of the multilayered GFRP composites that can be evaluated by the THz-TDS system. The specific results of the research and conclusions are summarized as follows: (1) The multi-delamination defects were successfully simulated by inserting layers of PTFE in the GFRP composite samples and the optical parameters for the different materials in the samples (i.e., compact GFRP composite and PTFE) were calculated. (2) The SWT-based algorithm was used to denoise the measured THz signals and both the thickness and in-depth location of each delamination defect were successfully calculated using analytical and experimental approaches. (3) The wavelet shrinkage de-noising algorithm was proven to remove the fluctuations in the THz signals effectively while equally preserving the features associated with the multi-delamination defects. The denoised THz signals are clearer than the measured THz signals and both the in-depth locations and thicknesses of the different delamination defects were accurately calculated. (4) A comparison between the actual and the measured thickness values of the delamination defects before and after the wavelet shrinkage denoising process indicates that the latter is highly accurate with less than 3.712% relative error, while the same error of non-de-noised signals reaches 16.388%. (5) This method could potentially be of great interest to the NDT&E engineers and practitioners using THz-TDS to remove undesired features associated with atmospheric material resonances obscuring them from correctly visualizing the features of the samples to be detected. That is, the difficulty and expense of purging the sample housing unit with dry nitrogen to physically remove the humidity features from the THz signals can be replaced by adequate signal processing such as this in many cases involving laboratory and field-based THz-TDS NDT&E applications. (6) Finally, the theoretical analysis provided in the present study would help NDT&E engineers and practitioners to estimate the maximum thickness of the GFRP composite materials and structural systems that can be evaluated by the THz-TDS.