Temperature Patterns in TSA for Different Frequencies and Material Properties: A FEM Approach

Abstract

Thermography techniques are gaining popularity in structural integrity monitoring and analysis of mechanical systems’ behavior because they are contactless, non-intrusive, rapidly deployable, applicable to structures under harsh environments, and can be performed on-site. More so, the use of image optical techniques has grown quickly over the past several decades due to the progress in the digital camera, infrared camera, and computational power. This work focuses on thermoelastic stress analysis (TSA), and its main goal was to create a computational model based on the finite element method that simulates this technique, to evaluate and quantify how the changes in material properties, including orthotropic, affect the results of the stresses obtained with TSA. The numeric simulations were performed for two samples, compact and single lap joints. when comparing the numeric model developed with previous laboratory tests, the results showed a good representation of the stress test for both samples. The created model is applicable to various materials, including fiber-reinforced composites. This work also highlights the need to perform laboratory tests using anisotropic materials to better understand the TSA potential and improve the developed models.

1. Introduction

Structural integrity monitoring is a vital process of any mechanical system that analyzes its behavior in order to detect defects and prevent malfunctions and failures. With the major importance of structural integrity, there is interest in studying new and unconventional techniques with new characteristics and advantages. Thermography techniques are gaining popularity because they are contactless, real-time, non-intrusive, and can be performed on-site. The range of applications of these techniques is vast, going beyond structural integrity monitoring into failure detection, fatigue assessment, and residual stress measurement. There are some works that focus on fatigue assessment, determining fatigue limit, fatigue strength evaluation, damage location, and life prediction, showing the potential and advantages of thermography techniques [1,2,3]. The usage of thermography has also been employed to characterize material properties that are not directly measurable [4,5,6,7]. These measures can be used to obtain properties for numeric simulations or directly used in research or industrial applications. The use of image techniques has grown quickly over the past several decades due to the progress in the digital camera, infrared camera, and computational power [8,9]. Using image (or field) nature techniques enables the identification of strain concentrations and damage more accurately. This could be done for identifying the presence of defects or to quantify them. These analyses are more likely to be performed for high value and high-cost components, such as carbon-fiber-reinforced polymer (CFRP) [10,11]. One of the most interesting techniques in this field is thermoelastic stress analysis (TSA) [12,13,14,15,16].

TSA fundamentals are based on the thermoelastic effect, in which small temperature changes are induced in a solid component when a cyclic load is applied. The results of TSA depend on various parameters, including the load type (sinusoidal, square-wave, triangular, etc.), its amplitude/average, specimen mean temperature, and material type (isotropic, orthotropic or anisotropic). If adiabatic and reversible conditions are achieved, these temperature changes are proportional to the sum of the principal stresses. These temperature variations, which are in the order of tens of mK in metallic materials, can be measured at the surface of a component using infrared detectors. Therefore, by measuring the variations in temperature due to the applied cyclic load, the first stress invariant can be determined from the detector’s output. The stress fields obtained with computational modeling and the finite element method (FEM) are often used in comparisons, and the results obtained normally show a satisfying agreement [13,16,17,18,19,20].

Besides the typical applications of this technique, such as fatigue failure and crack analysis, a point of focus of TSA applications is their expansion to materials that are not isotropic, including composites, for example, which are of great interest to the engineering industry. Orthotropic materials have different property values in different directions, contrary to isotropic materials. This means that the considerations made in the classic thermoelastic theory are not applicable and the equations for these cases are different; the measured temperature variations are proportional to a linear combination of the principal stresses instead. The thermoelastic coefficient values change with direction also beacuse the coefficient of thermal expansion also varies. Besides this, the effects of the change on some thermal properties on the obtained stresses, such as the conductivity, are quite unclear. Regarding composite materials, TSA is more complicated, since the stresses depend on the constitutive material properties and the ply stacking sequence as well as the fiber orientation. Heat flux with the surroundings is usually negligible, but between plies of a laminate can be sufficient to negate the simpler thermoelasticity process of only measuring the surface stresses. Furthermore, the internal heat transfer effects can depend on the loading frequency. This influences the thermal patterns measured by TSA, in particular, in areas with great stress variation gradients. The TSA theory presented by William N.S. Jr. does not consider the influence of the stimulation frequency [12]. However, this has great influence in the TSA results [14] TSA allows one to analyze delamination failure in composites, the most common failure type in these materials. Estimation of the delamination crack length and crack growth can be carried out [21,22,23]. TSA can also evaluate fiber breakage and matrix cracking [24].

The relation between the thermoelastic theory and the material type is complex, and many studies focus on it. The main goal of this works was to take in account various effects that influence the obtained results, such as the mean stress at a crack tip, residual stresses, and non-isotropy. To take in account these effects, the authors propose new TSA models and equations that allow materials to undergo any loading conditions and orthotropic/anisotropic conditions too [25,26,27,28,29].

The goal for this work was to create a finite element-based computational model that simulates a material’s properties’ influences on TSA (including the frequency influence) that is applicable to orthotropic materials, including composites. The created model was tested in order to evaluate and quantify how the changes in material properties and load frequency affect the results of the obtained stresses with TSA. To have a better comprehension of this work’s goals and planifications, a workflow was designed, as shown in Figure 1. Throughout this document, the methodology of the model is presented first, followed by the procedure and results. The last part is reserved to results discussion and conclusions.

2. Methodology

This section describes the methodology used to create a computational model that simulates the TSA technique. A model was created using MATLAB software using the finite element method (FEM). The model and its approximations had a major influence on the results that are described by Erik G. Thompson in [30]. TSA relates the principal stresses’ sum and the temperature changes in a body due to the application of a cyclic load. Thus, the first step towards this goal is to perform a static elastic simulation to obtain a stress-sum distribution resulting from the application of a static load with a value equal to the amplitude value of the desired cyclic load. Next, the obtained results from this simulation are transformed into an external heat source to transpose the results to the next stage, a thermal simulation. Finally, with this information, and after data processing with the help of a fast Fourier transform (FFT), the signal amplitude is extracted, and the TSA equations can be applied. Then, a full-field stress-sum distribution map is obtained.

2.1. Elastic Simulation

Consider plane elasticity and its governing equation, Equation (1), where K is the global stiffness matrix, {u} is the unknown displacement’s vector, and {F} is the global external load vector. Then, Equation (1) can be applied to each element and then assembled to represent the entire domain [31].

$$\begin{array}{c}\hfill \left[K\right]·\left\{u\right\}=\left\{F\right\}\end{array}$$

(1)

The balance between size and accuracy was found considering only 2D analysis, more specifically in a plane-stress state. In the simulations, 4-node elements were used to balance the mesh detailing and computational cost [31]. The stiffness matrix for each element was calculated using Equation (2), where B, D, and J represent the strain, elasticity, and Jacobian matrices, and the following considerations were made:

• Usage of isoparametric elements, because it is suitable for numerical quadrature and systematic definition of the interpolation functions. This means a change from x and y coordinates to $\xi$ and $\eta$, with the relationship between both being a matrix of interpolation functions, N. The interpolation function $N_i$ is relative to node i, and its value is one at that node and null in the finite elements that do not share this node;
• Application of numerical quadrature, where the obtained integral for the calculation of K according to the FEM can be replaced by a summation of the integrated function applied at several Gauss points, $n_{GP}$ (with (${\xi }_g$, ${\eta }_g$) coordinates), multiplied by the respective weight, $w_{GP}$. The number of Gauss points chosen was four, meaning full integration (instead of 1 GP, reduced integration) because shear locking is not a problem in this study, and to prevent the hourglass effect [31].

The Jacobian matrix is obtained by multiplying a matrix of the $N_i$ interpolation functions’ partial derivatives of each node x and y coordinates. The strain matrix is obtained by multiplying the inverse of the Jacobian matrix by the same partial derivative matrix used before.

$$\begin{array}{c}\hfill \left[K\right]=\sum _{g=1}^{n_{GP}}\left[B^T({\xi }_g,{\eta }_g)DB({\xi }_g,{\eta }_g)detJ({\xi }_g,{\eta }_g)w_g\right]\end{array}$$

(2)

Since a plane stress state is being considered, and orthotropic materials, the elasticity matrix is given by Equation (3), where $E_i$ represents the Young’s modulus and $v_{ij}$ the Poisson’s coefficient in i and j directions. The shear modulus, $G_{xy}$, is given by Equation (4) [32,33].

$$\begin{array}{c}\hfill \left[D\right]=\left[\begin{array}{ccc}\frac{E_x}{1-{\upsilon }_{xy}·{\upsilon }_{yx}}& \frac{{\upsilon }_{xy}·E_y}{1-{\upsilon }_{xy}·{\upsilon }_{yx}}& 0\\ \frac{{\upsilon }_{xy}·E_y}{1-{\upsilon }_{xy}·{\upsilon }_{yx}}& \frac{E_y}{1-{\upsilon }_{xy}·{\upsilon }_{yx}}& 0\\ 0& 0& G_{xy}\end{array}\right]\end{array}$$

(3)

$$\begin{array}{c}\hfill G_{xy}=\frac{E_x·E_y}{E_x+E_y+2·E_y·{\upsilon }_{xy}}\end{array}$$

(4)

The external loading vector was calculated by applying Equation (5). As one can see, both external forces applied over the element area and surface traction, respectively, b and $\overline{t}$, take part in this relation. Included in the boundary conditions are also restricted movements, which can be applied in the global stiffness matrix by multiplying the value by ${10}^{12}$, in the corresponding restricted nodal direction, originating an overly excessive stiffness and preventing movement in that direction. With both K and {F} calculated, the nodal displacements are obtained with the use of Equation (1). With these values known, both deformations and stress vectors can be easily calculated by applying Equation (6).

$$\begin{array}{c}\hfill \left\{F\right\}=\sum _{g=1}^{n_{GP}}\left[N^T({\xi }_g,{\eta }_g)bdetJ({\xi }_g,{\eta }_g)w_g\right]+\sum _{g=1}^{n_{GP}}\left[N^T({\xi }_g,{\eta }_g)\overline{t}det{J}_{\sigma }({\xi }_g,{\eta }_g)w_g\right]\end{array}$$

(5)

$$\begin{array}{c}\hfill \left\{\sigma \right\}=\left[\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\tau }_{xy}\end{array}\right]=\left[D\right]·\left\{\epsilon \right\}=\left[B\right]·\left\{u\right\}\end{array}$$

(6)

2.2. Thermal Simulation

With the FEM approximations, the problem comes in the form of Equation (7), where C is the mass matrix, {$T\left(t\right)$} is the temperature with respect to time, {$T^{\prime }\left(t\right)$} is the temperature’s time derivative, and {Q} is the external heat vector. This analysis is divided into two steps: (1) solving Equation (7) to extract the temperature derivative with respect to time, for $t=0$ s; and (2) using that result to calculate the temperature at the next time increment [34,35].

$$\begin{array}{c}\hfill \left[K\right]·\left\{T\left(t\right)\right\}+\left[C\right]\left\{T^{\prime }\left(t\right)\right\}=\left\{Q\right\}\end{array}$$

(7)

In order to calculate {$T^{\prime }\left(t\right)$}, we need $t=0$ s and {$T\left(0\right)$}, and is assumed that the analyzed components are initially at an ambient temperature of $T_{\infty }=297.5$ K. The matrices and vectors are given by Equations (8) and (9). In the stiffness matrix equation, {k} is the conductivity vector containing the conductivity values in both x and y directions, since the materials considered are orthotropic. Regarding the mass matrix equation, $\rho$ represents the density of the material, and $c_p$ the specific heat capacity at constant pressure.

$$\begin{array}{c}\hfill \left[K\right]=\sum _{g=1}^{n_{GP}}\left[B^T({\xi }_g,{\eta }_g)\left\{k\right\}B({\xi }_g,{\eta }_g)detJ({\xi }_g,{\eta }_g)w_g\right]\end{array}$$

(8)

$$\begin{array}{c}\hfill \left[C\right]=\sum _{g=1}^{n_{GP}}\left[N^T({\xi }_g,{\eta }_g)\rho c_{p}detJ({\xi }_g,{\eta }_g)w_g\right]\end{array}$$

(9)

For the external heat, Equation (10) is used, with Q being the effect of natural convection. An external heat source, proportional to the stress from the elastic simulation, is applied. In this relation (Equation (11)), h is the natural convection coefficient, with a chosen value of $h=5$ Hz (a normal value for this situation [36]); $K_Q$ is a coefficient to relate both simulations, with the value of $K_Q=1758$ (value obtained in previous work done by the same authors [20]); and f is the frequency of the applied load. The external heat source is applied at all nodes, and the convection is only applied to the nodes on the boundaries.

$$\begin{array}{c}\hfill \left\{Q\right\}=\sum _{g=1}^{n_{GP}}\left[N^T({\xi }_g,{\eta }_g)QdetJ({\xi }_g,{\eta }_g)w_g\right]\end{array}$$

(10)

$$\begin{array}{c}\hfill \left\{Q\right\}=h·area·(T\left(t\right)-T_{\infty })+\frac{K_Q·[{\sigma }_{xx}+{\sigma }_{yy}]}{\mathrm{\Delta }t}sin(2\pi ·f·t)\end{array}$$

(11)

The time increment is obtained using the forward difference method. This method becomes unstable if the value of $\mathrm{\Delta }t$ is too high; see Equation (12). Here, $\tau$ represents the signal period and n is the total number of cycles. The maximum number of increments, $MNI$, was established as 500, for 25 cycles, avoiding a computationally expensive simulation, and the number of iterations per cycle chosen was 20, which is a number big enough to accurately describe a cycle. Finally, {$T(t+\mathrm{\Delta }t)$} is given by Equation (13) [30]. In the end, a matrix of the nodal temperatures at each time increment is obtained, $T\left(t\right)$.

$$\begin{array}{c}\hfill \mathrm{\Delta }t=\frac{\tau ·n}{MNI}\end{array}$$

(12)

$$\begin{array}{c}\hfill \left\{T(t+\mathrm{\Delta }t)\right\}=\left\{T\left(t\right)\right\}+\left\{T^{\prime }\left(t\right)\right\}·\mathrm{\Delta }t\end{array}$$

(13)

2.3. TSA Stress

To extract the temperature signal amplitude, $\mathrm{\Delta }T$, a fast Fourier transform (FFT) function is applied. The TSA stresses are obtained using Equation (14), with ${\alpha }_i$ being the thermal expansion coefficient for the i direction. This equation is a simplification of the general TSA equation, applied to an orthotropic solid [21,22,23]. Throught this work, and in the majority of the cases, a uniform $\alpha$ was considered. In these cases, Equation (14) can be simplified, with $K_{TSA}$ being the thermoelastic coefficient.

$$\begin{array}{c}\hfill \rho c_p\frac{\mathrm{\Delta }T}{T}=({\alpha }_x{\sigma }_{xx}+{\alpha }_y{\sigma }_{yy})=\alpha [{\sigma }_{xx}+{\sigma }_{yy}]\equiv [{\sigma }_{xx}+{\sigma }_{yy}]=\frac{\left\{\mathrm{\Delta }T\right\}}{K_{TSA}·\left\{T\right\}}\end{array}$$

(14)

The amplitude and phase information is not related to the first stress invariant, and so it is very difficult to extract information regarding the stress values in both principal directions to use in Equation (14). To solve this problem, the chosen approach was to consider a stress sum vector with (${\sigma }_{xx}$, ${\sigma }_{yy}$) coordinates, and we considered both stresses as versors, $\overrightarrow{{\sigma }_{xx}}$ and $\overrightarrow{{\sigma }_{yy}}$. With this approach, Equation (15) was the one used to extract the first stress invariant with this TSA simulation.

$$\begin{array}{c}\hfill \rho c_p\frac{\mathrm{\Delta }T}{T}=[{\alpha }_x\overrightarrow{{\sigma }_{xx}}({\sigma }_{xx}+{\sigma }_{yy})+{\alpha }_y\overrightarrow{{\sigma }_{yy}}({\sigma }_{xx}+{\sigma }_{yy})]\end{array}$$

(15)

3. Procedure

The methodology presented in the previous section was applied to two different components: a compact tension (CT) specimen and single lap joints (SLJ). The tests performed on the CT specimen had the intention of explaining the effects of the changes in material properties on the obtained results, and guaranteeing that orthotropy was well established in the TSA model created. The simulations performed on the SLJ had the intention of demonstrating a TSA simulation on a composite material, more specifically, with the purpose of trying to see if TSA is applicable for detecting internal delamination failures.

3.1. CT Specimen

The CT specimen’s dimensions following the ASTM E 647-08 standard and an illustrative drawing, with the respective dimension values, are presented in Figure 2. With respect to this part, a mesh was designed, which is shown in Figure 3a. Although this mesh does not possess a structured form, the obtained results should be accurate enough for this purpose, even more so due to the small dimensions of the elements near the crack tip area, which is the area of most interest. The mesh refinement was performed using the same method as in other studies [9,37]. As the orthotropic characteristic is one of the points of focus, the analysis is centered a lot around two specific nodes near the crack tip, node 11 and node 130, presented in Figure 3b. Both nodes are at the same distance from the crack tip, $0.1$ mm, so there are no major differences when comparing the results.

To perform the elastic simulations, boundary conditions have to be known. The boundary conditions considered are illustrated in Figure 4. Those were applied by differing various types of nodes:

• Node type 0—free movement and no load applied.
• Node type 1—free movement and load applied.
• Node type 2—constrained movement in the y direction and no load applied.
• Node type 3—constrained movement in both directions and no load applied.

A fixed node (node type 3) is necessary to prevent the parts free movement, that being the bottom left corner one. It is worth noting that the distributed load applied assumes a parabola-type behavior, as can be observed in Figure 4. Regarding the thermal simulation, the boundary conditions applied relate to the convection, and that was only applied at the boundary nodes, as already explained.

3.2. Single Lap Joint

Regarding the SLJ, the dimensions, boundary conditions, and material properties were based on the work done in [38]. All the dimensions used in this SLJ specimens are presented in the illustration of Figure 5; the x and y directions being as represented in Figure 6, and the z direction represents the direction of the thickness. According to the figure shown, $0\le x\le 240$ mm, $0\le y\le 6.6$ mm, and $0\le z\le 25$ mm. Differently from [38], we added a defect to this part, more specifically, an air bubble in the resin layer that connected both carbon fiber reinforced plastic (CFRP) plates, in the middle of the specimen, in order to study the potential of the TSA technique for detecting failures in composite materials. When designing the mesh, the air bubble was taken in account, and so, three identical meshes were designed, changing only in the length of the defect, called mesh 1, mesh 2, and mesh 3, respectively. The meshes were designed considering a middle plane at $z=12.5$ mm. Mesh 1 has no defect present, for comparison purposes. The bubble was considered an ellipse, and since the layer of resin was only $0.2$ mm thick, the radius in the y direction was kept constant ($r_y=0.1$ mm). The value of the radius in the x direction: $r_x=2.5$ mm for mesh 2 and $r_x=7.5$ mm for mesh 3. Since the three meshes were identical, only mesh 1 is fully presented at the bottom of Figure 6. At the top of the figure, the zooms around the defect for mesh 2 and mesh 3 are presented. The meshes, away from the air bubble, had a structured form, and near it the, element was small enough and should produce good results.

When it comes to the boundary conditions applied, the method applied was the same as to the CT specimen, i.e, considering the same different node types. A schematic illustration of the applied boundary conditions is presented in Figure 7. The restriction in the y direction that can be observed at the end of the plates corresponds to an alignment tab.

To perform the simulations, the material properties were needed, and so simplifications were made. It was to consider the part as one and made from the same material: an CFRP/epoxy composite (only the values for the Young’s modulus and Poisson’s coefficient are relative to a ply of CFRP) with the properties extracted from [5,38,39], which are presented in

Table 1. CFRP’s mechanical properties.

Yield Strength (MPa)	${\mathit{E}}_{\mathit{x}}$ (GPa)	${\mathit{E}}_{\mathit{y}}$ (GPa)	Poisson’s Coefficient (-)	Density (kg/m3)
1230	109	8.819	0.342	1410

and

Table 2. CFRP’s thermal properties.

${\mathit{k}}_{\mathit{x}}$ (W/m·K)	${\mathit{k}}_{\mathit{y}}$ (W/m·K)	Specific Heat Capacity (J/kg·K)	${\mathit{\alpha }}_{\mathit{x}}$ $({10}^{-6} {\mathbf{K}}^{-1})$	${\mathit{\alpha }}_{\mathit{y}}$ $({10}^{-6} {\mathbf{K}}^{-1})$
6.3	0.6	1130	21.3	67.6

. It was considered that the Poisson’s ratio was the same in all directions and the yield strength used was an average value, which was considered constant as well in all directions.

4. Results

4.1. CT Specimen

The first test that was performed consisted in considering an isotropic material (in this case, aluminum 2219-T851). This means that the properties were the same in both x and y directions, and as presented in

Table 3. Aluminum 2219-T851’s properties [40].

Yield Strength (MPa)	Young’s Modulus (GPa)	Poisson’s Coefficient (-)	Thermal Conductivity (W/m·K)	Density (kg/m ${}^3$)	Specific Heat Capacity (J/Kg·K)
352	73.1	0.33	120	2840	864

. Additionally, this means that Equation (14) was used with $K_{TSA}=7.7\times {10}^{-6}$ for the case of aluminum.

The goal of this first analysis was to understand the effect of conductivity, and so the thermal conductivity value was changed in each simulation performed to one the values presented in

Table 4. Thermal conductivity values for the performed test.

Parameter	Value 1	Value 2	Value 3	Value 4	Value 5
Thermal Conductivity (W/m·K)	0.12	12	120	1200	12,000

. The magnitude of differences between the chosen values was large in order to really see the effect of this parameter. After these first simulations were done, the temperature variation with time was plotted for both above-referenced nodes. As the lower conductivity values produced nearly no effect on the temperature variation, remaining almost constant, only the plots of the highest conductivity values are shown in Figure 8a,b.

Besides these graphics, the stress profile along the symmetry axis is also plotted, see Figure 9, with a zoom around the crack tip area.

The second experiment that was conducted had the objectives of: (a) seeing if the model was property simulating the different values in the two directions; (b) analyzing the effect of having different conductivity’s in x and y directions. To realize these simulations, all properties used were the same as before; only the conductivity in the x direction, $k_x$, was changed, to the values presented in Table 4. As before, the temperature variation with time was plotted for different cases, as shown in Figure 10a,b. As explained above, only the plots containing relevant information are shown. The obtained stress as a function of the loading frequency was also plotted—presented in Figure 11. This graphic was plotted for both node 11 and 130, since the values should be different for the two cases.

A third test was performed considering only orthotropical elastic properties. The previous experiments focused only on the conductivity. Thus, the main goal of this test was to access the effect that variation in the elastic properties has on the simulations. To achieve this goal, the simulations were run considering the same properties presented in Table 3; however, the Young’s modulus in the x direction, $E_x$, was changed in each simulation. The values chosen are presented in

Table 5. Young’s modulus values for the tests.

Parameter	Value 1	Value 2	Value 3	Value 4	Value 5
Young’s Modulus (GPa)	10	70	130	200	300

. Just as before, the differences in the selected values were large in order to accentuate the changes from each simulation.

The temperature’s variation with time was once again plotted, for the cases of interest, and is shown in Figure 12a–c. Similarly to what was shown before, stress profiles along the symmetry axis were also plotted; see Figure 13. Finally, a stress distribution map was traced for each case, to make more noticeable the effect of the change in elastic properties. These maps are illustrated in Figure 14a–c.

The fourth and final experiment realized in the simulations using the CT specimen combined changes in conductivity and elastic properties, to see how these two aspects would affect the TSA simulation together. All properties remained constant with the same values as before, while varying only $k_x$ and $E_x$ using the values presented in Table 4 and Table 5, respectively.

The temperature’s variation with time was plotted for every case tested. However, due to there being too much information and much irrelevant data, since the effects of both could already be seen from the previous tests, only two plots are presented for the most extreme cases simulated; see Figure 15a,b. Other data, such as stress plots and distribution maps, are not shown either, due to the same reason.

4.2. Single Lap Joint

The SLJ specimen was submitted to TSA simulation with a loading frequency of 20 Hz. The three meshes were tested. The overall stress distribution maps are very similar; only in the area around the defect changes. Thus, the global TSA stress distribution map is only presented for mesh 1, but the stress distribution around the defect is present for both mesh 2 and mesh 3; see Figure 16.

As the TSA technique only allows surface measurements, the stress analysis on one of the two plies surfaces was of most interest. Therefore, the stress profile along the bottom surface of the bottom ply was plotted for the three different meshes, as presented in Figure 17. The vertical lines correspond to the beginning and end of the defect, for the respective mesh, with the center of the defect being at $x=120$ mm.

5. Discussion

When considering an orthotropic material, and comparing the temperature variations obtained for two different frequencies, 5 and 50 Hz, in Figure 8a and Figure 8b, respectively, it is noticeable that loading frequency has an effect on the conductivity, as shown by the more clear temperature signal attenuation at low frequencies, just as was expected, since with higher frequencies, heat has less time to propagate. Besides this, a further look at these two plots can reveal that nodes 11 and 130 have the same attenuation, although the temperature values are a little lower in the vertical direction. When it comes to the conductivity value itself, a little attenuation is noticeable at the $120$ W/m·K, growing bigger when increasing this value. When analyzing the stress profiles plotted in Figure 9, it can be seen that the change in the conductivity values barely affects the stress values obtained, because although the attenuation of the temperature signal is altered, its amplitude remains nearly the same. Thus, this is not responsible for the thermal patterns observed by different frequencies.

In the second test, the conductivity value remained constant in the y direction, and it was changed at each simulation for the x direction. As for the case of $k_x=\mathrm{12,000}$ W/m·K, the effects are more noticeable; that was the case used for analysis—see Figure 10a,b. For this case, it was expected that at node 11, the conductivity effect would be more accentuated than at node 130. However, it seems that was not the case. Actually, when comparing to the isotropic case shown before, the results are very similar. Regarding the stress plot for different loading frequencies, Figure 11, it is observable that the different values and cases have no clear effect on the stress behavior with the loading frequency. This raises the question, what is the influence of the frequency on the TSA results? This issue was reported in previous works by the authors [14] and requires deeper research. This was the case not only for mode I but also for mode II and mode III. Looking at the TSA theory, the load frequency is not considered an influencing parameter. The results presented here, and the experimental ones obtained previously, indicate that this consideration is incomplete. As such, it is imperative that future research tackles this issue in full detail, especially for when the TSA is used with high stress values and high stress gradients.

Regarding the changes in Young’s modulus, its effects can be clearly seen by looking at the obtained stress distribution maps. Compared to the isotropic case, Figure 14b, where the stress distribution has a heart shape, the decrease in Young’s modulus in the x direction tightens in this form, Figure 14a, and the increase in this value stretches across this stress distribution, Figure 14c. These changes are also visible when analyzing the stress profiles; see Figure 13. However, is also noticeable that the stresses near the crack tip tend practically to the same values. This was expected, since the mesh used for the simulations was the same, and so it is normal that the plots tend to have the same value. When looking at the temperature variations, Figure 12a–c, it is clear that these changes only affect the amplitude of the signal, as expected. Still, the differences between the temperature plot for $E_x=70$ GPa and for $E_x=300$ GPa are almost nonexistent for node 11, although the difference in values is big. This also proves the tendency for the same stress values around the crack tip.

The fourth and last test, varying both conductivity and Young’s modulus in the x direction, was more proof of what was concluded from the previous tests. The temperature variations observed in Figure 15a,b show the expected behavior considering the above-presented discussion.

Finally, the created TSA model was tested with a composite material and applied to the SLJ specimen. The main goal of this simulation was to evaluate if the TSA technique could be used to detect the presence of defects, such as delamination, in these types of materials. The stress distribution maps obtained, Figure 16, were as expected, with the stress concentrations being bigger as the defect grew bigger. However, there are some areas, specifically around the defect tips, where the stress distribution has sudden stress changes. This is due to the discretization process inherent in the FEM. Different meshes will produce different results, but with similar patterns. For this reason, the TSA measures and high-stress analysis should be performed in areas without stress concentration. An alternative and future work will be to perform the simulations using meshless techniques. Since the technique allows surface measurement, the stress profile on the surface of the bottom ply was plotted in Figure 17. By analyzing these plots overall, it can be seen that the defect presence is noted in the surface stress profiles obtained with the TSA simulation. The presence of the defect deviates from the stress values when compared to the simulation without the defect. It can be concluded that the bigger the defect, the more detectable the presence of the defect. The exact locations of the beginning and the end of the defect are not precise, but a good estimation can be made. The detection of it is noticeable, meaning that TSA is capable of detecting the presence of a defect in composite materials.

The model conceived in this work was based on experimental tests performed with a FLIR SC7500. This camera has an NEDT (noise equivalent differential temperature) of less than 25 mK. As such, and at the current time, the authors cannot guarantee that errors in the simulations are smaller than 25 mK. Assuming the behaviors such as other applications of lock-in processing and similar test conditions, the errors between the simulations and the inferred laboratory conditions are expected to be smaller than 1.0 mK for both CT and SLJ. This is based on the authors’ experience.

6. Conclusions

This work had the main goal of creating a computational model based on the finite element method that simulates TSA applied to orthotropic materials, in order to study the effects of the material’s properties on the results and to evaluate the application of this technique to different materials, including composites. The methodology based on FEM was presented firstly. Then, the model was used on a CT specimen; the results helped to better understand the effects of the material’s properties and changes in its values. Finally, an SLJ specimen was simulated by the model. It was composed of a composite orthotropic material, so the potential of the application of the present model was tested. The simulations done showed that defect detection is possible with this technique.

As a final remark, the created model acceptably reproduced the TSA technique. The simulations helped to better understand material’s properties’ effects on the obtained results and to prove this technique’s potential for applications for orthotropic materials, including composites. However, there is still room for improvement and to perfect this model, especially considering the difficulty of the TSA analysis of these materials. Future works should consider anisotropic materials and the inherent concentration of stresses in certain areas, along with these areas’ effects on the material models and the TSA readings. Both numeric simulations and experimental testing should be involved.