An analytical approach for frequency modulated thermal wave imaging for testing and evaluation of glass fiber reinforced polymers

Thermal non-destructive testing and evaluation (NDT&E) is a reliable, non-contact, whole field testing technique widely used for the assessment of defects such as voids and delaminations. Thermal NDT enables to observe surface temperature gradient on the sample surface using infrared camera to detect presence of any subsurface anomalies. Active infrared thermography (AIRT) is an attractive NDT tool used for detecting flaws caused during various stages from manufacturing to in-service applications [1–3]. The popular optically stimulated active thermal NDT (TNDT) techniques are pulsed thermography (PT) [2–5], pulse phased thermography (PPT) [6–11] and lock-in thermography (LT) [12–14].

PT employs short duration input pulse as heat flux over test sample to measure thermal gradient. Although, PT is a simple and fast technique, its temperature response gets affected due to variations in surface emissivity and non-uniform heating [4, 5] over the test specimen. PPT works on same principle (experimental setup) as PT except the data analysis scheme which is carried out using Fourier transform to extract phase related information from reflected thermal wave. The results obtained from the phase analysis approach are being less sensitive to surface features such as emissivity and non-uniform heating [7–11]. The modulated thermographic technique such as LT utilizes a mono-frequency input thermal excitation of considerably low peak power to compute magnitude and phase information of reflected thermal wave. But LT demands repetition of experiments at different input frequencies to inspect defects located at different depths inside the test sample which make it a time consuming technique [12–14]. Therefore, to overcome the limitations of conventional TNDT techniques such as peak power requirement, penetration depth and resolution, a non-periodic linear frequency modulated (LFM) thermal excitation scheme was introduced.

The fundamental theory of linear FMTWI (LFMTWI) allows heating the test sample with a low peak power heat sources having suitable band of frequencies chosen on the basis of thermal properties and its geometry [15–19]. The thermal waves diffuse into the test sample due to the absorbed heat over the sample. Therefore, a change in surface temperature gradient can be observed in the presence of any subsurface discontinuity. The present paper focuses on the quantitative analysis of the temperature gradient at different depth locations in a glass fiber reinforced polymer (GFRP) sample using a novel analytical tool by analyzing the temperature T(v,t) as a function of position and time observed at the surface v = 0 of test sample. The present paper also emphasizes the time domain Hilbert transform data processing approach using LFMTWI which provides better depth resolution.

The proposed work in this paper has been arranged as follows: section 2 presents principle of LFMTWI technique and analytical solution for heat transfer model under adiabatic boundary conditions. The GFRP sample has been modeled and simulated using finite element modeling software is presented in section 3. Further, time domain data analysis scheme has been explained in this section. Section 4 deals with results obtained from the proposed approach and its comparison with simulated results. Final conclusions are reported in section 5.

In LFMTWI, frequency modulated thermal waves diffuses into the test sample within a pre-defined band of frequencies having same energy. This applied heat flux diffuse into the test sample to generate thermal waves and produces a time-varying temperature distribution over the test sample. The existence of sub-surface anomalies may alter heat distribution inside the solid leading to a thermal contrast over the test sample.

The thermal waves generated in a homogeneous finite thickness GFRP sample are studied analytically with the help of a theoretical model. The one-dimensional heat conduction equation describing the temperature distribution without any external source or sink in the test sample assuming constant thermal properties is given by [16]:

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{\rm{T}}\left({\rm{v}},{\rm{t}}\right)}{\partial {{\rm{v}}}^{2}}=\displaystyle \frac{1}{\alpha }\displaystyle \frac{\partial {\rm{T}}\left({\rm{v}},{\rm{t}}\right)}{\partial {\rm{t}}};0\leqslant {\rm{v}}\leqslant {\rm{l}},{\rm{t}}\gt 0\end{eqnarray} \tag{ 1 }$$

where T(v, t) depicts instantaneous thermal profile along spatial coordinate v (m) at time t (s). α denotes thermal diffusivity (m² s⁻¹) of GFRP sample given by α = k_p/ρC_p. k_p- thermal conductivity (W m⁻¹-K⁻¹), ρ-density (kg m⁻³) and C_p-specific heat (kJ kg⁻¹-K⁻¹).

Analytical solution for equation (1) was obtained under adiabatic boundary conditions (BC's) i.e. negligible heat flux is exchanged at the other end 'l' of the sample thickness. The heat conduction problem for finite thickness was formulated as per the following time domain boundary conditions:

$$\begin{eqnarray}&&{\left.-{{\rm{k}}}_{{\rm{p}}}\displaystyle \frac{\partial {\rm{T}}}{\partial {\rm{v}}}\right|}_{{\rm{v}}=0}={\rm{W}}({\rm{t}}),t\gt 0\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\left.-{{\rm{k}}}_{{\rm{p}}}\displaystyle \frac{\partial {\rm{T}}}{\partial {\rm{v}}}\right|}_{{\rm{v}}={\rm{l}}}=0,t\gt 0\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&T(v,t)=To,t=0\end{eqnarray} \tag{ 4 }$$

where T_o is the test sample initial temperature and W(t) defines LFM input heat flux shown induced at the surface of the GFRP sample expressed as [16, 18]:

$$\begin{eqnarray}&&-{{\rm{k}}}_{{\rm{p}}}{\left.\displaystyle \frac{\partial {\rm{T}}({\rm{v}},{\rm{t}})}{\partial {\rm{v}}}\right|}_{{\rm{v}}=0}={{\rm{W}}}_{{\rm{o}}}\,* \,\left[\sin \left\{2\pi \left({{\rm{f}}}_{{\rm{i}}}{\rm{t}}+\displaystyle \frac{{{\rm{Bt}}}^{2}}{2}\right)\right\}\right]\end{eqnarray} \tag{ 5 }$$

where W_o depicts applied heat flux peak value, B is given by B = (f_i − f_o)/T_d is the ratio of difference between initial and the final frequency (bandwidth) to the total duration. The frequency sweep is chosen in order to have large diffusion lengths allowing deep penetration into the entire test sample and T_d is the total duration for the applied heat flux.

The fundamental solution of the homogeneous heat conduction equation for finite thickness solid sample by imposing time dependent adiabatic boundary conditions (2–4) is given in the form of:

$$\begin{eqnarray}&&{\rm{T}}\left({\rm{v}},{\rm{t}}\right)={{\rm{T}}}_{{\rm{a}}}\left({\rm{v}},{\rm{t}}\right)+{{\rm{T}}}_{{\rm{b}}}\left({\rm{v}},{\rm{t}}\right)\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&{{\rm{T}}}_{{\rm{a}}}\left({\rm{v}},{\rm{t}}\right)=\displaystyle \sum _{{\rm{n}}=0}^{\infty }{{\rm{R}}}_{{\rm{n}}}\,\cos \left[\lambda {\rm{v}}\right]{{\rm{e}}}^{-t{\rm{\alpha }}{\lambda }^{2}};\end{eqnarray} \tag{ 7 }$$

and

$$\begin{eqnarray}&&{{\rm{T}}}_{{\rm{b}}}\left({\rm{v}},{\rm{t}}\right)=\displaystyle \sum _{{\rm{n}}=0}^{\infty }{\int }_{0}^{{\rm{t}}}{{\rm{S}}}_{{\rm{n}}}\left({\rm{\Pi }}\right)\cos \left[\lambda {\rm{v}}\right]{{\rm{e}}}^{-\left[\left({\rm{t}}-{\rm{\Pi }}\right)\alpha {\lambda }^{2}\right]}{\rm{d}}{\rm{\Pi }}\end{eqnarray} \tag{ 8 }$$

where λ are the Eigen values given by; λ = nπ/l with n = 0,1,2,...∞. Each value of λ corresponds to an independent solution satisfying the given heat conduction equation and boundary conditions stated above. П is a new introduced parameter unrelated to t in order to obtain the auxiliary solution. R_n and S_n are the integration constant values obtained by satisfying the initial condition as follows:

$$\begin{eqnarray}&&{{\rm{R}}}_{{\rm{n}}}=\displaystyle \frac{{\rm{a}}}{{\rm{l}}}\displaystyle {\int }_{0}^{{\rm{l}}}\left({\rm{T}}({\rm{v}},0)\cos [\lambda {\rm{v}}]\right){\rm{dv}};\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{{\rm{S}}}_{{\rm{n}}}({\rm{\Pi }})=\displaystyle \frac{{\rm{b}}}{{\rm{l}}}\displaystyle {\int }_{0}^{{\rm{l}}}\left({\rm{P}}({\rm{v}},{\rm{t}})\cos [\lambda {\rm{v}}]\right){\rm{dv}}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\rm{P}}({\rm{v}},{\rm{t}})=\displaystyle \frac{\alpha }{{{\rm{k}}}_{{\rm{p}}}}\displaystyle \frac{{\rm{W}}({\rm{t}})}{{\rm{l}}}+{\rm{W}}^{\prime} ({\rm{t}})\left(\displaystyle \frac{{{\rm{v}}}^{2}}{2{\rm{l}}}-{\rm{v}}\right)\end{eqnarray} \tag{ 11 }$$

where ${\rm{a}}=\left\{\begin{array}{l}1\,{\rm{for}}\,{\rm{n}}=0;\\ \,2\,{\rm{elsewhere}}\end{array}\right.,{\rm{b}}=\left\{\begin{array}{l}1\,{\rm{for}}\,{\rm{n}}=0;\\ \,2\,{\rm{elsewhere}}\end{array}\right.$

The temperature gradient obtained over the sample by neglecting the radiation and convection losses as a function of time at a given depth by applying LFM heat flux is expressed as:

$$\begin{eqnarray}&&\begin{array}{l}{\rm{T}}\left(0,{\rm{t}}\right)={\rm{T}}\left({\rm{v}},0\right)-{{\rm{T}}}_{{\rm{v}}}\left(0,{\rm{t}}\right)\displaystyle \frac{{\rm{l}}}{3}\\ +\,\displaystyle \frac{{{\rm{W}}}_{{\rm{o}}}}{2\sqrt{{\rm{B}}}{{\rm{k}}}_{{\rm{p}}}}\left(\begin{array}{l}\cos \left[{\chi }^{2}\pi /2\right]\left(-\,{\rm{FresnelS}}\left[\chi \right]+{\rm{FresnelS}}\left[\varepsilon \right]\right)+\\ \sin \left[{\chi }^{2}\pi /2\right]\left({\rm{FresnelC}}\left[\chi \right]-{\rm{FresnelC}}\left[\varepsilon \right]\right)\end{array}\right)\\ +\,\displaystyle \frac{{\left(-1\right)}^{n}{{\rm{W}}}_{{\rm{o}}}{\rm{l}}}{16\sqrt{{\rm{B}}}{{\rm{k}}}_{{\rm{p}}}{{\rm{n}}}^{3}{\pi }^{3}}\left(\begin{array}{l}{\eta }_{1}\left(\left({e}^{\displaystyle \frac{{\lambda }^{2}\alpha {{\rm{f}}}_{{\rm{i}}}}{2{\rm{B}}}}\left({{\varsigma }_{{\rm{a}}}}^{-1}-{{\varsigma }_{{\rm{b}}}}^{-1}\right)\right.\right)+\left(\left.{e}^{i2\pi t\left({{\rm{f}}}_{{\rm{i}}}+{\rm{Bt}}\right)}\left({{\varsigma }_{{\rm{c}}}}^{-1}-{{\varsigma }_{{\rm{b}}}}^{-1}\right)\right)\right)+\\ {\left(-1\right)}^{1/4}{\eta }_{2}\left(\left({{\rm{e}}}^{{\psi }_{1}+{\psi }_{2}}\left({\rm{Erf}}\left[{\gamma }_{1}\right]-{\rm{Erf}}\left[{\phi }_{1}\right]\right)\right.\right)+\\ \left.\left.{\left({{\rm{e}}}^{{\psi }_{1}+{\psi }_{2}}\left({\varsigma }_{{\rm{a}}}{\varsigma }_{{\rm{b}}}\right.\right)}^{-1}\left({\varsigma }_{{\rm{a}}}{\rm{Erf}}\left[{\rm{i}}{\gamma }_{2}\right]-{\rm{i}}{\varsigma }_{{\rm{b}}}{\rm{Erf}}\left[{\rm{i}}{\phi }_{2}\right]\right)\right)\right)\end{array}\right)\end{array}\end{eqnarray} \tag{ 12 }$$

where S represents FresnalS, C represents FresnalC respectively. While acquiring solution for surface temperature T(0,t), the following set of equations has been obtained:

$$\begin{eqnarray}&&\begin{array}{l}\chi =\displaystyle \frac{{{\rm{f}}}_{{\rm{i}}}}{\sqrt{{\rm{B}}}};\\ \varepsilon =\displaystyle \frac{2{\rm{Bt}}+{{\rm{f}}}_{i}}{\sqrt{{\rm{B}}}};\\ {\phi }_{1}=\left(\displaystyle \frac{1+{\rm{i}}}{4\sqrt{{\rm{B}}\pi }}\right)\left(\alpha {\lambda }^{2}+2{\rm{i}}\pi {{\rm{f}}}_{i}\right);\\ {\phi }_{2}=\left(\displaystyle \frac{1+{\rm{i}}}{4\sqrt{{\rm{B}}\pi }}\right)\left(\alpha {\lambda }^{2}-2{\rm{i}}\pi {{\rm{f}}}_{{\rm{i}}}\right);\end{array}\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray*}&&\begin{array}{l}{\gamma }_{1}=\left(\displaystyle \frac{1+{\rm{i}}}{4\sqrt{{\rm{B}}\pi }}\right)\left(\alpha {\lambda }^{2}+{\rm{i2}}\pi ({{\rm{f}}}_{{\rm{i}}}+2{\rm{Bt}})\right);\\ {\gamma }_{2}=\left(\displaystyle \frac{1+{\rm{i}}}{4\sqrt{{\rm{B}}\pi }}\right)\left(\alpha {\lambda }^{2}-{\rm{i2}}\pi ({f}_{i}+2{\rm{Bt}})\right);\\ {\eta }_{1}={(-1)}^{n}{\rm{in}}\pi 16\sqrt{{\rm{B}}};\\ {\eta }_{2}={(-1)}^{n}{\rm{in}}\pi \alpha {\lambda }^{2}\sqrt{32}\end{array}\end{eqnarray*}$$

This analytical solution is useful in the analysis of temperature variations at various depth locations as a function of time for identification of subsurface defects presents at any given location v > 0 as shown in figure 1.

Zoom In Zoom Out Reset image size

The present work explores a finite element method (FEM) based modeling and simulation approach to validate the proposed analytical method. A three dimensional model based on FEM has been designed for GFRP sample having four flat bottom holes each having diameter of 6 mm having thermal properties shown in table 1. The layout of the modeled GFRP test sample is shown in figure 2. Initially, the test sample of finite thickness 0 ≤ v ≤ l is placed at a particular temperature T_o and while the top surface was subjected to 2 kW LFM input heat flux with a sweep of 0.01–0.1 Hz for a duration of 100 s. The mathematical formulation of the time-dependent heat conduction equation has been considered assuming no internal heat generation within the solid sample. The temperature measurements are recorded at various depths ranging from 0.5–1.2 mm illustrated in figure 2.

Zoom In Zoom Out Reset image size

The temperature gradient was recorded for 100 s duration to validate the proposed analytical method. Additionally, zero mean thermal response for the above mentioned defect depths was formulated by using an appropriate polynomial fitted curve in order to carry out any post-processing technique.

Thus obtained mean removed temperature response is then processed by computing the time domain based cross-correlation (CC) approach between the chosen reference signals (deepest depth) with respect to the other depth locations to determine the potency of relation between the chosen variables. The cross-correlation for the time domain is given as follows:

$$\begin{eqnarray}&&{\rm{CC}}(\tau )={\rm{T}}\left({\rm{t}}\right)\otimes {\rm{R}}\left({\rm{t}}+\tau \right)\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{\rm{CC}}(\tau )=\displaystyle {\int }_{-\infty }^{+\infty }{\rm{T}}\left({\rm{t}}\right)\,* \,{\rm{R}}\left({\rm{t}}+\tau \right){\rm{dt}}\end{eqnarray} \tag{ 15 }$$

where R(t) is the reference signal, T(t) is the temperature response for various defect depths, ⨂ denotes correlation, t is the time and τ is the time delay parameter.

These mean removed temperature distribution curves are then used to construct the compressed pulses with the help of cross-correlation technique [19]. The Fourier domain CC was computed between the reference (sound) signal and the temperature recorded for different defect depths given as [20–22]:

$$\begin{eqnarray}&&{\rm{CC}}(\tau )={\Im }^{-1}\left[{\rm{R}}{(\omega )}^{\ast }{\rm{T}}(\omega )\right]\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{\rm{R}}(\omega )=\Im \left[{\rm{R}}({\rm{t}})\right];{\rm{T}}(\omega )=\Im \left[{\rm{T}}({\rm{t}})\right]\end{eqnarray} \tag{ 17 }$$

where * denotes complex conjugate, ℑ and ℑ⁻¹ denotes fast Fourier transform (FFT) and inverse fast FT (IFFT) between the reference signal R(t) and the defect depths temperature response T(t).

The CC output facilitates the entire energy of the FM signal into a narrow sinc shaped signals of high amplitude. This obtained cross-correlated image enables high depth resolution of the surface and sub-surface defects. The lag time for main lobe and corresponding sample defect depths can be calculated in the form of:

$$\begin{eqnarray}&&{{\rm{Lag}}}_{{\rm{cc}}}={\left.\tau \right|}_{{\rm{cc}}(\tau )=\max [{\rm{cc}}(\tau )]}\end{eqnarray} \tag{ 18 }$$

Further, CC-phase for the captured thermal response was obtained using Hilbert transform (HT) where all the frequency components are shifted by an angle of 90 degrees [17, 18].

$$\begin{eqnarray}&&{\rm{HT}}({\rm{t}})=\displaystyle \frac{1}{\pi }\displaystyle {\int }_{-\infty }^{+\infty }\left[\displaystyle \frac{{\rm{R}}({\rm{v}},\tau )}{{\rm{t}}-\tau }\right]{\rm{d}}\tau \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{{\rm{CC}}}_{{\rm{HT}}}(\tau )={\Im }^{-1}\left[{{\rm{R}}}_{{\rm{HT}}}{(\omega )}^{\ast }{\rm{T}}(\omega )\right];\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{{\rm{R}}}_{{\rm{HT}}}(\omega )=\Im \left[{\rm{HT}}\left\{{\rm{R}}({\rm{t}})\right\}\right]\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{{\rm{R}}}_{{\rm{HT}}}(\omega )=-{\rm{i}}\mathrm{sgn}(\omega ){\rm{R}}(\omega )\end{eqnarray} \tag{ 22 }$$

where R_HT is the Hilbert transform of the reference sound signal and sgn denotes the signum function.

The time domain CC-phase can be constructed between the chosen sound signal and captured defect depths given as follows:

$$\begin{eqnarray}&&{\rm{\Phi }}={\left.{\tan }^{-1}\left[\displaystyle \frac{{\rm{CC}}(\tau )}{{{\rm{CC}}}_{{\rm{HT}}}(\tau )}\right]\right|}_{\tau =0}\end{eqnarray} \tag{ 23 }$$

The block diagram for constructing the time domain CC-coefficient and phase information using Hilbert transform is shown in figure 3.

Zoom In Zoom Out Reset image size

An analytical study for the LFM incident heat flux over the surface of the GFRP test sample with frequency sweep 0.01–0.1 Hz for 100 s duration has been carried out. The GFRP sample has four flat bottom holes in the range of 0.5 mm-1.2 mm each having 6 mm diameter is shown in figure 2. The study was carried out under adiabatic boundary conditions and temperature profile for various defect depth locations has been recorded using equation (12) with the thermal properties as given in table 1. Figure 4(a) depicts the analytically obtained thermal responses at different depth locations. Since the obtained thermal profiles have mean (rise) shifted, therefore in order to carry out the proposed post-processing, the zero mean thermal responses have been obtained using the first-order polynomial fit as shown in figure 4(b).

Zoom In Zoom Out Reset image size

The time domain based CC coefficient approach was adopted on the zero mean temporal thermal profiles to obtain the pulse compression for highly depth resolved imaging of subsurface defects shown in figures 5(a), (b). The CC was calculated for different defect depths of the range 0.5 mm- 1.2 mm in a GFRP sample by considering the 6.95 mm as the reference thermal response.

Zoom In Zoom Out Reset image size

Figure 6 depicts the normalized correlation coefficients calculated for different defect depths for the data obtained from the FEM simulator and proposed analytical solution. It is clear from the trend obtained from the CC-coefficient thermal profiles; there is a significant reduction in the coefficient values as a function of depth.

Zoom In Zoom Out Reset image size

Further, time domain based phase information is retrieved by the Hilbert transform based post processing to reconstruct the imaginary component of the reference temperature distribution for analyzing the phase response for LFMTWI as shown in figure 3. Finally, depth resolvability for LFMTWI was demonstrated in figures 7(a), (b) in terms of (a) CC-lag and (b) CC-phase while keeping reference as 6.95 mm thickness of GFRP test sample.

Zoom In Zoom Out Reset image size

The results depicted in figure 7(a) indicate high depth resolvability for the computed analytical solution in comparison to the FEM based simulator. The CC-lag time (ms) induced to the presence of defects in the GFRP test sample calculated from equation (18) decreases with the increase in depth and converge to zero for the deepest defect.

Furthermore, the CC-phase contrast for LFMTWI model is depicted in figure 7(b) for FEM simulator and proposed solution. A lower phase value was observed initially for the shallowest defect and then leads to higher phase value with the increase in defect depths. The presented results suggest its capability for better or comparable defect depth resolvability for LFMTWI input heat flux.

The present paper highlights a novel analytical treatment for LFMTWI and validated with the commercially available FEM software with adiabatic BC's which can be treated for practical experimental scenario. The results obtained for a GFRP test sample having four flat bottom defect holes are considered and the proposed post-processing schemes have been analyzed to prove the defect detection capabilities of the FMTWI. Depth resolvability for the above solution in terms of cross-correlation coefficient and its resultant time lag and phase are discussed. Furthermore, the provided analytical representation has been validated with the numerical results obtained from the three dimensional finite element models to prove the applicability of the proposed analytical model.

No new data were created or analysed in this study.