Transient thermal analysis of layered media based on thermal quadrupoles

The transient heat transfer within a layered structure can be modeled efficiently using the thermal quadrupole methods. Such numerical simulations would contribute to effective experimental design for thermography inspection applications of composite materials. Current research is part of such effort that focuses on the parametric study on both a two-layer model and a four-layer model. The variation of amplitude curves agrees with the lateral position of the interfacial discontinuity in general. The sudden increase in the phase value distinctively indicates the end point of the same discontinuity. Hence the lateral location of an embedded interfacial discontinuity can be readily identified from the Laplace surface temperature. The layer thickness does not have significant effect on the general trend observed in the parametric analysis presented in this work. The advantage of the numerical modeling of layered structures using thermal quadrupoles lies in that thermal quadrupoles can be easily extended to multiple layers containing internal anomaly between different layers.


INTRODUCTION
Fiber reinforced polymer composite materials (GFRP or CFRP) have widespread use in both fabrication of light-weight components and repair of damaged engineering structures. However new applications also present new challenges to structural integrity assessment and nondestructive evaluation (NDE). A visual examination in which no recordable indications of flaws are detected does not mean that no structurally significant flaws exist within the structure. More advanced nondestructive inspection (NDI) techniques such as ultrasonic scanning and X-ray computed tomography thus are sometimes applied to the damage assessment of CFRP components. The demand of continuous monitoring of progressive damage, however, presents some challenges to the applications of advanced NDI methods.
Composite structural elements under variable loading could produce heat due to thermal expansion effect and change of stress state. The former heating effect is likely caused by the irreversible thermoelastic plastic deformation. The latter is often resulted from the material discontinuity or embedded defects . Friction caused by two partially dis-bonded surfaces would also generate heat. The subsequent re-distribution of thermal stress field would act as internal heat sources that dissipate energy to the surface. Such inhomogeneous surface temperature may be detected by passive thermography. The recorded thermal images have gained some success on revealing embedded voids, de-bonding interface, or even progressive damages (Harizi et al., 2014;Zalameda and Jackson, 2020). Thermal analysis plays an important part in the effectiveness of thermography testing because of the nature of internal heat sources in a lot of cases. The detection of non-repetitive transient thermal signals can be quite difficult for passive thermography due to unwanted heat transfer and small temperature contrast and other issues.
The transient heat transfer within a layered structure can be modeled using a variety of numerical techniques, such as finite difference methods, finite element analysis and thermal quadrupole methods. The flexibility of thermal quadrupoles has been well documented, for example Winfree et al. (2018). In brief, such a technique allows the construction of an exact and direct model that is suitable for experimental design and inversion of thermal properties (Maillet et al., 2000). The solutions to the heat equations are expressed as linear matrices that allow numerical results be obtained more efficiently than finite element analysis. This is very appealing when experimental design requires simulations of the multi-layered structures such as CFRP components . Furthermore, the need for improving the quantitative analysis of internal defects based on the results of passive thermography also requires more efficient thermal modeling and simulations than those in the literature (Shiozawa et al., 2017;Palumbo et al., 2017).
The current research starts with the introduction to thermal quadrupoles. A two-layer model is analyzed that results in solution to the surface temperature in the Laplace domain. The effect of an interfacial defect is explored according to the number of layers and the locations of discontinuity at the interface between two adjacent layers. The simulated results for both the two-layer and four-layer models will be presented, followed by discussions and concluding remarks.

QUADRUPOLES FOR THERMAL ANALYSIS
The concept of quadrupoles was first introduced to thermal analysis by Maillet and co-workers in 1993. The temperature at any given time of a lumped body model can be expressed as the difference between input and output heat fluxes as, (1) The time derivative is removed after Laplace transform is applied to temperature T and heat flux Ф in Equation (1). The simplest case for such a model can then be written in a matrix formulation as, , and Ct is the heat capacity.
As an analogy to a discretized electrical network, Equation (2) is represented by the diagram shown as Fig. 1.
The capital letters inside the square indicate elements of the transfer matrix. The temperature variation along the z direction of the layer thickness is written as, λ is the heat conduction coefficient and S is the surface area perpendicular to the z direction. The heat flux can be determined as, (4) Let the heat transfer from z = 0 to z = e, the corresponding heat equation is written in the matrix form as,

MODELLING OF HEAT TRANSFER IN LAYERED MEDIA
To construct a thermal quadrupole model for a layered media, the electrical circuit analogy is again applied to describing the relation between temperature and flux that go through two contacting layers and the interface between them, as shown in Fig. 2.
Given that the surface temperature is the key feature for most experimental study using infrared thermography, we will only focus on the thermal modeling in the z direction of a layered media. The temperature θ (0, s) and heat flux ϕ (0, s) at the surface are determined by inverting the transfer matrix of thermal quadrupoles as shown in Equation (6), Similarly, the equations for a two-layer system are . 2. Equivalent diagram for a two-layer system (after Maillet et al., 2000) Fig. 3. Schematics of a two-layer system containing an interfacial discontinuity, shown in red The surface temperature θ (0, s) is obtained as, The advantage of thermal quadrupoles is now clear that a more complicated system can be readily analyzed by extending Equation (7) to model multiple layers. The added complexity of solutions in the Laplace domain is certainly worth the effort. Furthermore, the applications to numerical simulation are relatively straightforward, as can be seen in the following section.
The symbols used and their units are listed in the Appendix as Table A1.

MODELLING RESULTS OF TWO-LAYER AND FOUR-LAYER SYSTEMS
The aforementioned analysis forms the basis of numerical simulation that has been successfully implemented using a computer code developed by the authors. The Laplace surface temperature (0, ) can be inverted to the time domain using the Stehfest numerical algorithm (Maillet et al., 2000). We now present the results for various configurations in a two-layer system (Fig. 3). The exact dimensions of each configuration are listed in Table 1. L is 2 mm in the lateral direction of the model.
The data shown in Fig. 4 to Fig. 6 are computed based on a thermal conductivity of 6 W/mK and the thermal diffusivity in z direction = 6.493X10 -7 m 2 /s and that in x direction = 3.250X10 -6 m 2 /s, respectively. The time history of the simulated thermal response to the periodical load at the interface is shown in Fig. 4(a) to Fig. 7(a) for  Fig. 4(b). The amplitude of the Laplace surface temperature at various positions along the x direction on the surface corresponding to the peak values are indicated by the blue dots in Fig. 4(a). The time difference of the blue dots with respect to the maximum peak is the phase change with respect to the maximum peak. The normalized phase change is plotted as the red curves in Fig. 4(b). The area in brown indicates the position of interface discontinuity along the x direction.  The variation of amplitude curves agrees with the lateral location of the interfacial discontinuity in general. On the other hand, the sudden increase in the phase value distinctively indicates the end point of the same discontinuity. Hence the lateral location of the discontinuity can be identified in such a way. The drawback of the phase variation approach is that the farther away from the discontinuity in the lateral direction the smaller the Laplace amplitude. This makes it more difficult to determine the accurate phase value. Further, the results obtained from a thinner upper layer appear to be less stable as shown in Fig.  6.
The results of six configurations are presented in Fig. 7 and Fig. 8 for a four-layer model. A thickness of 0.05 mm is selected for all four layers. The exact dimensions for the configurations analyzed are listed in Table 2. In general, the results are in agreement with those in the same configuration of a two-layer model. Difference in the phase variation does exist between F-1 and D-1. Such deviation can be expected since the farther away from the discontinuity in the lateral direction the smaller the Laplace amplitude. Same rationale applies to the difference between F-6 and D-9. Such deviations are of no concern to the objective of finding the lateral location of the discontinuity.

CONCLUDING REMARKS
The thermal analysis based on thermal quadrupoles has been successfully applied to both a two-layer model and a four-layer model. The variation of Laplace temperature amplitude curves agrees with the lateral position of the interfacial discontinuity in general. The sudden increase in the phase value distinctively indicates the end point of the same discontinuity. Hence the lateral location of the discontinuity can be readily identified from the surface. The layer thickness does not have significant effect on the general trend observed in the parametric analysis presented in this work. The advantage of the numerical modeling of layered structures using thermal quadrupoles lies in that thermal quadrupoles can be easily extended to multiple layers containing internal anomaly between different layers.
The proposed analysis would be valuable in applying passive infrared thermography to CFRP components under loading conditions. We will report in a separate study about the outcome of the experimental results and detailed analysis.