Numerical simulation of subsurface defect identification by pulsed thermography and improvement of this technique for noisy data

Pulsed thermography is an active non-destructive technique which uses optical excitation source to stimulate heating of the object under investigation. This work is devoted to the simulation of the pulsed thermography method in a steel plate with a ceramic coating containing artificial defects of various depths and sizes. The simulation has been carried out on the base of a model which takes into account complex heat exchange of the sample with the surrounding by convection, conduction and radiation. Comparison of the temperature contrast with the experimental data has shown that the results are in a good qualitative and quantitative agreement in all stages of the cooling process. Due to the fact that the temperature contrast is often susceptible to surface noise of various nature the Kalman-based signal processing technique was developed. The comparative analysis has shown that the proposed filtration technique provides a better signal-to-noise ratio in comparison to the considered well-known techniques of signal reconstruction when proper calibration of the filtration parameters is carried out.


INTRODUCTION
eramic coatings are widely applied in various fields of industry as thermal barriers which protect structures from high-temperature influence. They provide corrosion defense, improve the reliability of metallic details and have high resistance to fuels, lubricants and other chemical substances. The technological process of the coating consists of the adding of the one or several layers of ceramic particles on a solid surface. Apparently, adhesion of the coating with the main structure is not uniform. Initiation of structural defects such as delamination, voids, inclusions, cracks and other heterogeneities can substantially decrease the efficiency of the protection. Accurate identification of coating defects will improve the operational reliability of the structure and reduce the risk of disastrous consequences induced by its failure. Pulsed thermography is an active non-destructive technique which uses an optical excitation source (usually, xenon lamp or laser) to stimulate heating of the object under investigation. The registration of the temperature during cooling is carried out by an infrared camera which can be positioned on the same side with the excitation source (reflection mode) or on the opposite side (transmission mode). The duration of the heating is significantly less than the observation time. Thermal properties of the subsurface defect (air) are different from the properties of the main material. Therefore, the diffusion rate of the structure in case of the presence of the defect is reduced which reflects in higher values of the surface temperature above the discontinuity. The main parameter which is applied to the analysis of the data obtained during pulsed thermography is a temperature contrast. This value is estimated as the difference between the mean surface temperature of the area above the discontinuity and the mean temperature of the sound area (temperature in the reference region without defects). Temperature contrast depends on the size of the specific defect and its depth and usually susceptible to the surface noise induced by non-uniform heating, variation of emissivity and reflections from the environment [1][2]. Therefore, accurate prediction of the defect location requires the development of relatively simple and effective techniques of signal processing. Pulsed thermography can be applied to the identification of subsurface defects, their depths and sizes in metals, ceramics and composite materials [3][4][5][6] while mathematical models can provide a deeper understanding of this process and obtain optimal parameters for the experimental investigations. In [7] numerical simulations were performed to obtain optimal values of excitation parameters with respect to the evaluation of the paint coating thickness. A numerical model of the pulsed thermography for heterogeneous materials has been developed in [8] on the base of the finite volume method. The significant influence of such parameters as irradiation power density, non-uniform heating and geometric characteristics of the defects on the value of the temperature contrast was established with the use of the proposed model. M. Grosso et al. [9] applied a pulsed thermography model to the investigation on the limits of this method with regard to the adhesive composite joint. The sensitivity of this method to the thickness of the material, heating time and excitation energy was studied and the limiting values of the defect depths were obtained. R. C. Waugh et al. [10] developed a finite-element model of pulsed thermography procedure and applied it to several materials (aluminum, carbon fibre-reinforced plastic and adhesively bonded joints). The obtained results have shown that the main limiting factor of the pulsed thermography is the uncertainty in thermal material properties. In [11] an algorithm for an accurate simulation of the pulsed thermography is proposed. This algorithm updates several parameters (material properties, pulse duration, heat power and ambient temperature) using smoothed and refined experimental data. The proposed algorithm demonstrated a large optimization of the calculation time. However, not all of the material properties were defined properly and the accuracy of the prediction was affected by signal noise. More recent approaches for identification of subsurface defects are based on adaptive algorithms and artificial neuron networks. In [12] an artificial neuron network is applied to the estimation of the defect's depth. Results of the numerical simulation obtained according to the mathematical model of the pulsed thermography procedure were used for generation of ideal (without noise) sets of training data. Verification of the proposed algorithm was based on its application to the experimentally obtained data. It has been shown that the accuracy of the prediction in the worst case is equal to 90%. In this work, a three-dimensional numerical model of the pulsed thermography method is proposed for the evaluation of the temperature contrast in steel plate with a ceramic coating containing artificial defects of various depths and sizes which represent delamination of the coating from the main structure. The developed model takes into account complex heat exchange which includes heat transfer by convection, conduction and radiation. The verification of the model was based on the experimental results presented in [13]. The first part of the paper is devoted to the investigation of influence of various forms and durations of the heating pulse on the dynamics of the temperature evolution. In the second part of the paper, the processing technique for thermal contrast data has been developed. The approach is based on the Kalman filter adapted to the problem of subsurface defect detection by pulsed thermography. The proposed methodology has C been compared to the main well-known and widely used techniques of signal reconstruction (simple moving average method, Gaussian filter, Savitzky-Golay filter and median filter). Values of the temperature contrast evaluated by the numerical simulations in the first part of the paper were used as ideal reference signals and thus, allow us to estimate the efficiency of the proposed approach.

Mathematical formulation
he evolution of the temperature field during heating of the solid is described by the heat transfer equation: is the initial temperature of the solid. The following boundary conditions are applied to the heated surface F of the solid: where n is the unit normal to The fundamental base of non-destructive techniques is the temperature contrast. The most commonly used definition of this parameter refers to the absolute temperature contrast and states that it is equal to the difference between the mean temperature in the defective area and the mean temperature in the non-defective (sound) area [2]: where C is the value of absolute temperature contrast, d T is the mean temperature in the defective area, s T is the mean temperature in the sound area. This parameter characterizes the visibility of the specific defect but it is susceptible to surface noise [2]. Therefore, an accurate prediction of this value using simple and effective signal processing techniques should be carried out. In this work, model (1)-(5) was applied to the estimation of the "ideal" (without noise) temperature contrast evolution in sample with subsurface defects of various depths and sizes.

Validation of the model
The verification of the model was based on the simulation of the experiment published in [13]. Pulsed thermography was applied to evaluate of the temperature contrast in rectangular AISI 316 L steel plate containing artificial square defects of various sizes located at various depths from the surface. The plate has a length of 150 mm and a width of 100 mm. The thickness of the sample is 3.5 mm. A schematic representation of the specimen is given in Fig. 1 (a). T The numerical simulation was carried out with the finite-element software COMSOL Multiphysics®. The plate was heated by a square pulse with the amplitude o Q equal to 4·10 5 W/m 2 . The descending part of the pulse was smoothed to simulate gradual decrease of the incident heat flux. The temperature was registered for a time interval of 2 seconds with a step size of 0.01 s. The following boundary conditions were applied to the specimen: (3)-(5) to the front surface and (4)- (5) to the lateral and rear surfaces. To simulate non-uniform heating of the front surface induced by the convexity of the heating source, the incident heat flux was multiplied by the piecewise function which scales the amplitude o Q from 0.88 o Q to 1.03 o Q along the width of the specimen. Subsurface defects were represented as parts of the specimen filled with air. Continuity conditions were applied at the interface between steel and air. Physical properties as well as values of parameters used in boundary and initial conditions are given in Tab. 1. The modelled specimen was discretized by tetrahedral finite elements with various sizes. A more refined mesh was used in defective zones of the sample. The convergence analysis has shown that the optimal size of the elements in areas adjacent to the defects is equal to 0.6 mm. The maximum size of the elements in the finite-element model was 8 mm. The complete mesh consisted of 1250000 elements. The mesh is presented in Fig. 1 Table 1: Thermophysical properties of air, steel and parameters of heating.
. It can be seen that maximum values of the temperature contrast have defects located at the depth of 0.4 mm (the closest to the surface). The peak contrast time for the first row of the subsurface defects is reached at the beginning of the cooling process. A further decrease in the temperature contrast is due to diffusion. The obtained results illustrate that the peak contrast time depends on the depth of the defect that confirms the experimental results [4]. Defects located at depths of 3.36 mm and 3.17 mm from the surface could not been observed at the given time of 2 s. It should be noted that defects with the smallest size (2 mm) are not visible even at the depth of 0.4 mm from the surface.   3 shows the temperature contrast in defects of various sizes located at the depth of 1.13 mm obtained as a result of numerical simulation and experimental study conducted in [13]. It can be seen that the results are in a good qualitative and quantitative agreement in all stages of the cooling process.

Effect of excitation parameters on temperature contrast of the bi-layered sample with subsurface defects
Results provided by Fig. 3 shows that model (1)-(5) can be used for the prediction of the temperature contrast and peak contrast time. Therefore, it will be used for the investigation of heating time and shape of the heating pulse on values of these parameters for the same specimen coated with ceramic layer. The physical properties of the ceramic layer (Yttria stabilized zirconia) are given in Tab. 2 [14].

Property Value Unit
Density 5600 kg/m 3 Thermal conductivity 0.9 W/(m·K) Heat capacity 505 J/(kg·K) The structure was considered as an object consisting of three materials: steel (the main detail), ceramics (top coating) and air (subsurface defects of various sizes representing delamination of the coating from the main detail). The sizes of the main detail and the defects correspond to the sizes presented in Fig. 1 (a). Fig. 4 shows a schematic representation of the cross-sectional area of the considered layered specimen. Subsurface defects have different thickness which varies from 0.14 mm to 3.5 mm (corresponds to the through defect). The following boundary conditions were applied: (3)- (5) to the coating and (4)-(5) to the lateral and rear surfaces of the object. Continuity conditions were used at the interface between the considered materials. The amplitude of the heating pulse o Q was equal to 4·10 5 W/m 2 in all cases. The finite-element model of the specimen with coating of 0.6 mm is shown in Fig. 5. As in the previous example, the mesh was considerably refined in the coating and in areas adjacent to the defects and sound area. The minimum element length in the coating was 0.07 mm. The maximum element length of the main detail was 10 mm. The complete mesh consisted of approximately 3·10 6 elements.    6 shows the effect of the heating duration on the temperature increase ΔT=T-T0 obtained after 2 seconds of cooling. The impulse had a square shape with a smoothed descending part to simulate gradual decrease in the heating power. The smoothing interval was equal to 6.6 ms for all considered cases. The thickness of the coating was equal to 0.6 mm. It can be seen that all subsurface defects are visible. Defects of the smallest thickness (the first column) have minor values of temperature increment. Fig. 7 demonstrates the temperature contrast for two utmost groups of defects located in the first column (defects with the smallest thickness) and the last column (through defects). It is interesting to note that defects with sizes of 10 mm, 8mm and 6 mm have nearly the same curves of temperature contrast evolution. Sensitivity of the temperature contrast to the size of the defect starts with the size equal to 4 mm. As it was expected, the smallest values of C in both columns have defects with the size of 2 mm. The difference in maximum contrast between them is approximately 7%. It should also be noted that there is a discrepancy in the peak contrast time: defects with smaller thickness reach the peak at earlier stages.  However, even with the relatively high value of the heating time (4 ms) the temperature contrast is less than 1 degree for the smallest size of the defect. In real experiments, the identification of this defect will be complicated due to the noise which can level off the peak.
The smallest defect (size of 2 mm) in the first column has been chosen as the representative of the worst case for the observations. Thus, further investigations will be carried out for it. The heating duration was assumed to be 4 ms due to the acceptable value of the obtained temperature contrast. Fig. 9 (a) illustrates the effect of the coating thickness on the temperature contrast for this defect. It can be seen that the maximum value of the temperature contrast is a non-linear function of the coating thickness while the peak contrast time is nearly linearly dependent on the thickness of the ceramic layer ( Fig. 9 (b)). It is also worth noting that the temperature contrast becomes almost insensitive to the coating thickness of more than 1 mm. In general, the obtained results show that pulsed thermography can be successfully applied to the defect identification in ceramic coatings because typical values of the thickness lies in the range from 0.1 mm to 0.6 mm [14]. For these values of the thickness, the maximum temperature contrast is reasonably high and is reached for the short period of time. For values of more than 1 mm the contrast is poor and the defect can be lost in noise.   10 illustrates the effect of the shape of the heating pulse on the value of C. The reason for the difference arises from the fact that the heating source can gain the maximum power gradually and cool down in the same manner. In the limiting case, this shape tends to be triangular. A comparative study of three shapes of heating pulses (square pulse without smoothing, square pulse with the smoothed descending branch and completely smoothed pulse close to the triangle) with the same duration of the heating by the maximum power (which is equal to 4 s) demonstrates that the difference between the obtained results is considerable. The temperature contrast obtained by the smoothed pulse is three times higher than the contrast achieved by the heating with the non-smoothed pulse. Hence, accurate description of the temperature contrast requires specific information on the inertia of the heating source. The conducted research has shown that pulsed thermography is sensitive to the size of the defect, its depth and the heating time. The smaller the defect is and the deeper it lies the higher power of heating is required to identify it. On the other hand, the noise that affects the signal will increase too [15]. Thus, it is necessary to develop simple and effective approaches for its reduction. The results obtained in this section will be used as a "pure" (reference) signals for the investigation of the efficiency of the filtration technique proposed in the next section.

Kalman-based filtration technique
he developed methodology is a continuation of the work [16] on Kalman-based filtration technique for signals obtained by thermal non-destructive testing. The main idea of Kalman filter is to estimate the current state of the dynamic system using information about its previous state and the value of its current measurement. The algorithm can be divided into two steps. In the first step (prediction) filter extrapolates state variables and their uncertainties. In the second step (correction) the result is refined. The algorithm can be used for the in-situ estimation of the object's state using only current measurements, data on its previous state and its uncertainty. The detailed description of Kalman filter can be found, for example, in [17]. In our case, the state variable is the temperature contrast C. The value of the temperature contrast k C on the k th time step can be estimated according to the following equation [18]: where   A t is the state transition function, 1 k C is the state of the temperature contrast on the previous step k-1, k w is the process noise. In general, A is the mxm matrix. A measurement k z of the k C is presented in the form: where H is a variable that relates the state to the measurement k z , k v is the measurement noise. In general, H is the lxm matrix which can change its value with each time step. Effective application of Kalman algorithm requires formulation of the model for the considered process. The model is based on two analytical solutions which are widely applied to the problems of pulsed thermal non-destructive testing. In case of Dirac pulse heating of the plate without subsurface defects, the temperature on the front surface can be evaluated according to equation [2]: T where e Q is the absorbed energy,  e c  is the effusivity, t is the time.
To find the temperature response of the surface in case of the presence of subsurface defects, the following relation can be applied [13]: where r R is the reflection coefficient, L is the depth of the subsurface defect,  / ( c)    is the thermal diffusivity.
Therefore, the temperature contrast can be evaluated as: (11) Hence, evolution of the temperature contrast can be written in the form: Using (12) and finite-difference form of the time derivative we can present Eqn. (7) in the form: where t is the step size. Thus, our system is described by Eqns. (13) and (8) which are the base for the Kalman filtration technique.
According to the Kalman algorithm [19], the initial values of the optimal state 0 a C and the error covariance 0 a P should be provided: where 0 C is the initial state, 0 P is the initial error covariance.
The next step is the prediction where values f k C and f k P are calculated: is the optimal estimation of k C on the k-1 time step, Q is the process noise covariance parameter.
Correction step of the Kalman filter includes calculation of the Kalman gain k K : where R is the measurement noise covariance. Updated values of the optimal estimation a k C and error covariance k P are calculated as: In order to improve values of the temperature contrast obtained by Kalman filter the Rauch-Tung-Steilbel smoother algorithm [20] was applied as a postprocessing treatment of the obtained data. The algorithm can be summarized as: where  1,1 k n , ˆk K is the smoothed gain, ˆk C is the smoothed value of the temperature contrast, ˆk P is the error covariance. Let's illustrate the efficiency of the proposed filtration algorithm. In the previous section, we have obtained reference (noise-free) signals. However, experimental values of temperature contrast are distorted by surface noise induced by nonuniform heating, variation of emissivity, reflections from the environment and so on. In order to simulate experimentally obtained signals we add Gaussian white noise to the reference signal. In this case, Eqn. (8) has the form: where k z represents noisy signal, k C is the reference signal obtained by numerical simulation, k v denotes white Gaussian noise. It is easy to notice that  1 H . Reliability of non-destructive techniques can be characterized by the signal to noise ratio. The higher this value is the smaller defects can be detected. Therefore, the main goal of the developed technique is to provide a relatively high value of the signal to noise ratio which can help to detect the presence of the defect. Signal to noise ratio (SNR) is defined as [21]: where    2 1 n signal k k P X is the power of the signal, is the power of the noise, k X is the k-th point of the discrete reference signal without noise, ' k X is the k-th point of the discrete signal obtained after the processing, n is the number of points in the discrete signal. Fig. 11 illustrates application of the proposed filtration technique to the signal with random additive Gaussian noise. Evolution of the temperature contrast obtained for the defect of 8 mm located at the depth of 0.6 mm (a green curve in Fig. 7 (a)) was used as a reference signal. Values of R=5 and Q=0.001 allows us to achieve a high SNR value which is equal to 929 and reconstruct the shape of the reference signal ( Fig. 11 (a)). Fig. 11 (b) displays one more application of the proposed technique to the filtration of the same reference signal corrupted with another noise. Results of the filtration also demonstrate a high value of SNR which is equal to 1150. Therefore, the proposed technique is insensitive to the initial noise and can provide accurate reconstruction of the signal for large sizes of defects. Fig. 12 presents results of filtration of the reference signal obtained for the defect with the minimal size of 2 mm (a purple curve in Fig. 7 (a)) using values of R=7 and Q=0.003. In this case, we have a much smaller value of the optimal SNR which is equal to 241. The reason for this decline is that the used analytical model for the temperature evaluation (10) does not take into account the size of the defect and depends only on the depth. To achieve more accurate results of filtration it is necessary to employ more precise analytical models for evolution of the temperature contrast. The analytical model (12) contains the depth of the defect as one of the governing parameter defining the behavior of the temperature contrast evolution. However, in real applications of pulsed thermography for defect detection this value is often unknown. Fig. 13 displays the effect of the false depth on the reconstructed signal. The reference signal corresponds to the numerical results of the temperature contrast obtained for 2 mm subsurface defect located at the depth of 2 mm. Fig. 13 (a) shows results of the reconstruction attained with the use of the precise depth of the defect and values of R=9, Q=0.03. It can be seen that the reconstructed signal repeats the shape of the reference signal with high accuracy). Fig. 13 (b) provides results obtained for the false depth of the defect L which was equal to 0.6 mm. Calibration of the R parameter allows us to reconstruct the shape of the signal close to the original. The applied value of R was equal to 6, the value of Q was the same. Thus, in case when the applied value of the depth is far from the proper, the inaccuracy of the model can be compensated by the measurement noise covariance.

Comparison with other signal processing techniques
The proposed filtration technique is formulated in the time domain. Hence, a comparative study will be carried out with the use of the four most widely applied time-domain signal processing techniques: simple moving average, Gaussian, Savitzky-Golay and median filtering. The brief description of each method is given below. For each of the considered techniques the optimal parameters providing maximum value of SNR had been obtained by exhaustive search method. The reference signal and corrupted signal correspond to the data presented in Fig. 11 (a) for all applied filtration techniques.
The simple moving average method is the most popular technique due to its simplicity. In equation form, it can be written as [22]: is the input signal, M is the number of points in average. Results of filtration obtained with the optimal value of  14 M are presented in Fig. 14 (a). Processing of the input signal by Gaussian filter can be described using convolution [23]: where  M is the kernel shift which is typically equal to the half-width of the kernel M ( is the kernel. A normalized Gaussian kernel is defined as: Application of the Gaussian filtering with the optimal values of M =15 and  =12 which provide the highest value of SNR for this filtration technique is demonstrated in Fig. 14 (b). The basic idea of the Savitzky-Golay filtration technique is to fit input data set by polynomials of some degree [24][25]. Coefficients of polynomials should provide the minimal mean-square error n e : n M e p n z n (25) where M is the half-width of the approximation interval,      0 N k k k p n a n is the polynomial of the N th degree, k a is polynomial coefficients.
Tab. 3 summarizes values of the SNR for each considered signal processing technique calculated using the optimal parameters of the filter. It can be seen that Savitzky-Golay and median filters give the best results and the worst result is obtained by simple moving average. However, the developed Kalman-based filtration method gives higher values of the SNR even for the case when the applied model of temperature contrast evolution describes the reference data incorrectly (Fig. 12). Moreover, the analysis conducted in the previous section has shown that proper calibration of the R and Q parameters as well as a hypothesis on the approximate depth of the defect allows one to obtain a very high value of the SNR (near one thousand or more).

Signal processing technique SNR
Simple moving average 48

Gaussian filter 55
Savitzky-Golay filter 215 Median filter 191 Kalman-based filter 929 Table 3: Signal to noise ratio obtained for the various filtration techniques applied to the detection of the 8 mm defect located at the depth of 0.6 mm.

CONCLUSION
n this study, a numerical simulation of subsurface defect identification by pulsed thermography is presented. The object under investigation is a steel plate with artificial defects of various sizes located under the ceramic coating. Subsurface defects have been modelled as parts of the specimen filled with air. The simulation has been carried out on the base of the model which takes into account complex heat exchange by convection, conduction and radiation. The verification has shown that the proposed model can capture the main features of the temperature distribution obtained by a pulsed heating of the sample. The developed model has been applied to the investigation of the influence of the heating parameters on the maximum value of the temperature contrast and the peak contrast time. Due to the fact that the temperature contrast is often susceptible to surface noise of various nature the Kalman-based signal processing technique was developed. The results obtained by numerical simulation were used as reference noise-free signals. The proposed method of filtration is based on the analytical solution to the problem of an infinite plate heated by a Dirac pulse. This solution allows us to propose a universal filtration procedure for signals obtained by pulsed thermography which does not require specific information on the heating parameters. The technique gives satisfactory results even if this information is unknown. In addition, the proposed methodology has been compared to the main well-known and widely used techniques of signal reconstruction. The obtained results have shown that:  Specific features of the heat source (the shape, rise and fall times) substantially affect the temperature contrast value and simulation results. This means that heating of the object by heat sources with the same power can give different values of the temperature contrast.  There is a threshold size of the defect with the specific depth for which the temperature contrast is universal. For example, in the considered problem the size is equal to 6 mm for a depth of 0.6 mm.  Defects with smaller thickness attain the peak contrast time earlier. For instance, defects with the size of 2 mm and a thickness of 0.14 mm reach the peak at t=0.9 s, while defects with a size of 2 mm and a thickness of 3.5 mm have the peak at t=1 s.  Decrease in the heating time leads to the substantial decline in the maximum value of the temperature contrast which can make experimental defect detection difficult due to the presence of the noise.  The temperature contrast is non-linearly dependent on the coating thickness while the peak contrast time depends on it almost linearly. I  The proposed Kalman-based filtration technique provides better values of signal-to-noise ratio (near 1000 or more) in comparison to the simple mean average, Gaussian filter, Savitzky-Golay filter and median filter when proper calibration of the process noise covariance, measurement noise covariance and defect depth is carried out.  The efficiency of the developed technique can be raised by applying the more accurate analytical model which takes into account dependence of the temperature evolution on the specific defect size.