Temperature patterns obtained in thermoelastic stress test at different frequencies, a FEM approach

Abstract

Purpose

The goal of this work is to create a computational finite element model to perform thermoelastic stress analysis (TSA) with the usage of a non-ideal load frequency, containing the effects of the material thermal properties.

Design/methodology/approach

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.

Findings

This work had the main goal to create a model to perform TSA with the usage of non-ideal loading frequencies, considering the materials’ thermal properties. Loading frequencies out of the ideal range were applied and the model showed capable of good results. The created model reproduced acceptably the TSA, with the desired conditions.

Originality/value

This work creates a model to perform TSA with the usage of non-ideal loading frequencies, considering the materials’ thermal properties.

Citation

Duarte, G., Neves, A.M.A. and Ramos Silva, A. (2023), "Temperature patterns obtained in thermoelastic stress test at different frequencies, a FEM approach", International Journal of Structural Integrity, Vol. 14 No. 3, pp. 354-377. https://doi.org/10.1108/IJSI-10-2022-0126

Nomenclature

V̇

Flowrate (m³ s⁻¹)

ρ

Density (kg m⁻³)

τ

Period (s)

c_p

Specific heat capacity at constant pressure (J kg⁻¹K⁻¹)

E

Young’s modulus (GPa)

x′

Derivative of variable x with respect to time

(ξ, η)

Isoparametric coordinates

[B]

Strain matrix

[D]

Elasticity matrix (GPa)

[J]

Jacobian matrix

[K]

Stiffness matrix (N mm⁻¹)

[N]

Matrix of interpolation functions

[C]

Mass matrix (kg)

[Q]

External heat vector (W)

α

Thermal expansion coefficient (K⁻¹)

ν

Dynamic viscosity (m² s⁻¹)

υ

Poisson’s ratio

{t¯}

Surface traction forces (N mm⁻¹)

{F}

Load vector (N)

{u}

Displacement vector (mm)

{ϵ}

Strain vector

{σ}

Stress vector (MPa)

{b}

External forces applied over the element area (N mm⁻²)

a

Area (m²)

d

Crack height (m)

f

Frequency (Hz)

h

Convection coefficient (W m⁻²K⁻¹)

k

Thermal conductivity (W m⁻¹K⁻¹)

K_Q

Conversion coefficient

K_TSA

Thermoelastic coefficient

l

Crack length (m)

MNI

Maximum number of increments

n

Number of cycles

Nu

Nusselt number

ppc

Points per cycle

Pr

Prandtl number

Re

Reynolds number

T

Temperature (K)

t

Time (s)

V

Volume (m³)

v

Velocity (m s⁻¹)

w

Specimen thickness (m)

x_g

Variable x at Gauss point

1. Introduction

Structural integrity monitoring is a vital procedure to ensure the proper functioning and security of a mechanical system. Maintenance interventions are very common since these evaluate the state of the components and detect defects, preventing malfunctions and failures. There are several techniques to perform these tests, from measurement with strain gauges to linear ultrasonic evaluation to detect cracks. 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, and beyond structural integrity monitoring including failure detection, they are helpful in fatigue assessment or residual stresses 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 (La Rosa and Risitano, 2000; Luong, 1998; Wang et al., 2017). The usage of thermography has also been employed to characterize material properties that are not directly measurable (Laguela et al., 2011; Tian and Cole, 2012; Villière et al., 2013; Strzałkowski et al., 2015). 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 fast over the past several decades due to the progress in the digital camera and infrared camera performance and increased computational power (Ramos Silva et al., 2019, 2022). 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 (Sharpe, 2008; An et al., 2014; Silva et al., 2017; Farahani et al., 2020; Ancona et al., 2015).

One of the most interesting thermographic techniques is thermoelastic stress analysis (TSA). TSA is an experimental thermograhpy technique whose fundamentals rely on the thermoelasticity theory. In fact, it is one of the most promising techniques, as it provides full field measurement from the surface of the component, such as the temperature field induced by a time-varying load, thus calculating the respective map of surface stress. 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 anisotropy (An et al., 2014; Díaz et al., 2020, 2004b).

This work involves a computational numerical simulation based on the thermoelasticity theory. This theory combines the elasticity of the materials with the respective thermal effects and vice-versa, i.e. temperature changes in a body are correlated with changes in stress. To establish the thermoelasticity theory, thermal aspects such as temperature and heat flux have to be included in the classical elasticity theory. This implies the introduction of an energy equation applied to an elastic solid continuum, respecting the principle of conservation of energy. Considering reversible, adiabatic thermodynamic events, including anisotropy, mean stresses and temperature-variant elasticity, equation (1) was derived, where ρ represents the density of the material, c_ɛ is the specific heat for constant strain, T is the temperature and C_ijkl, α_kl and ɛ_kl are, respectively, the elasticity, thermal expansion coefficient and strain tensors (Sharpe, 2008; Sadd, 2014; Díaz et al., 2004a).

(1)ρ⋅cε⋅δTT=∂Cijkl∂T⋅(εkl−αkl⋅ΔT)−Cijkl⋅αkl+ΔT⋅∂αkl∂T⋅dεij

Perfect adiabatic conditions can only be obtained if there are no stress gradients on the specimen and the thermal conductivity of the material is zero (or very close) meaning there is no heat flux. However, the non-adiabatic effects can be minimized, and consequently the heat flux, by applying a load frequency that is high enough to minimize the thermal diffusion, normally above the range of 5–10 Hz. Combining this with only working within the elastic range of the material, the energy conversion between thermal and mechanical forms can be treated as fully reversible (Sharpe, 2008).

There are other important aspects to have in account when performing TSA, such as the equipment used, specimen preparation and calibration. Usually, it is used as an infrared (IR) detector camera that is capable to collect the IR radiation emitted by the specimen, and it’s necessary a signal-processing unit that extracts the required information from the noisy IR detector output signal. The component is often painted with a high emissivity layer to have a uniform surface emissivity and increase its value. TSA information is normally presented as a vector, where the modulus is a magnitude that is proportional to the temperature variation experienced by the specimen and the phase denotes the angular shift between the thermoelastic signal and the reference. Calibration is important because a small change in values such as emissivity or the thermoelastic coefficient can lead to big changes in the obtained stresses (Sharpe, 2008; Silva et al., 2017; Díaz et al., 2004a).

One of the main limitations of the TSA is the frequency of the applied load, more specifically when the values become out of a certain range, out of the ideal range of 10 Hz. In the presence of high-frequency vibrations, there are two parameters that are very important: frame rate and integration time. High frequencies usually are higher than the available thermal camera frame rate. This issue can be tested by using a sub-sampling strategy. For example, in (Molina-Viedma et al., 2021) criterion applied was to choose the frame rate for every resonance ensuring at least 10 points per cycle, meaning that the frame rate should be much higher than the frequency of interest. Thus, the frame rate was chosen so that the closest multiple of the frame rate to the excitation frequency must be as close as possible to the excitation frequency itself. As for the integration time, the problem is more severe. For an IR camera, during the integration time, the sensor is receiving infrared radiation. A pixel is intended to represent a specific point of the analyzed component, with the objective to obtain its temperature variation over time. So, a long integration time will come with the blurring of the information as a variation of the pixel intensity occurs due to the position change of the specimen within the field-of-view of the camera. On the other hand, during a short integration time, low levels of IR radiation are received by the camera sensor and the thermoelastic signal is polluted with noise. The integration time should be many times shorter than the vibration period, in order to collect information in as many points per cycle as possible, so choosing this parameter, based on the previous explanation, it’s a challenging task (Molina-Viedma et al., 2021).

With the rise of loading frequency values, the quality of the results falls gradually. This degradation is due to the fact that with higher frequencies a lower integration time is used, which decreases the amount of IR radiation received by the detector and increases the noise level. This decreases the signal-to-noise ratio (S/N ratio). On other hand, increasing the frequency increases the stress peak values, particularly at stress concentration areas. Not also that, but the stress profile is also affected by the rise of the frequency values (Silva et al., 2017; Molina-Viedma et al., 2021). As a consequence, in the computational simulation of this test, these effects must be taken into account and submitted for comparison with experimental work.

The majority of works around TSA involve its advantages and applications. There are a few that focus on developing the technique, by trying to acknowledge some of its limitations and better understand some effects behind this technique. Thus, one of the objectives of this work is to try and improve the knowledge of the loading frequency effect on the TSA, considering the few information about this topic and its importance described above.

The main goal of this work is to create a computational finite element model to perform TSA with the usage of a non-ideal load frequency, containing the effects of the material thermal properties. Throughout this document, the methodology of the model is presented first, followed by the procedure and results. The last part is reserved for results, discussion and conclusions.

2. Methodology

In this section the followed methodology is presented. To have a global perception of the structure, a sequence illustration is shown in Figure 1.

2.1 Mesh

The TSA simulation was carried out using the finite element method (FEM). This is a numerical technique that approximates continuous functions as discrete models for obtaining estimate solutions to boundary value problems. This discretization involves a subdivision of the domain into an assembly of sub-regions called elements, linked together in a finite number of nodal points or nodes, in what is called the mesh (Pires et al., 2021). This is step one of the sequence presented in Figure 1.

The smaller the element size is, more accurate and reliable are the approximate results obtained. However, a more refined mesh also implies a bigger number of nodes, which means a bigger number of unknowns and a greater simulation time. Ideally, a balance between accuracy and simulation time is achieved in applying the mesh to this simulation. Beyond this, the type of element is also an important parameter of the mesh. Accounting the computational expense that 6-node and 8/9-node elements present, the choice, in this simulation, relies in 3-node or 4-node elements. Considering that 4-node elements allow a structured mesh design and originate better results, these elements are the ones that are going to be used. Also, the dimensions of the meshes need to be considered. The specimens to analyze are designed to be thin enough to consider a plane stress state, so the meshes can be designed in 2D. Various meshes will be designed and to make sure that the right ones are used, a mesh convergence has to be made. The different meshes will be tested and its results will be compared. The meshes selected should be the ones that require less computational effort (the ones with less elements and nodes and with a better arrangement, preferably with a structured form) within the meshes that converge at the same results.

2.2 Elastic simulation

After the mesh design, its information such as number of nodes and number of elements is necessary to perform simulations with the help of FEM. All simulations were performed using MATLAB using the approach and approximations described by Erik G. Thompson in (Thompson, 2004).

The first simulation is the elastic and it is the second step in the sequence presented in Figure 1. The goal is to obtain a matrix of nodal stress sum, [σ_xx + σ_yy]. The intended results are the nodal displacements u_x,i and u_y,i. FEM leads to an algebraic set of equations in the form of equation (2), where [K] is the global stiffness matrix, {u} is the displacements vector of unknowns and {F} is the global external load vector. Equation (2) can also be used for each individual element and then assembled together (Pires et al., 2021).

(2)[K]⋅{u}={F}

As is usually done in FEM, since is more convenient because is suitable for numerical quadrature and for systematic definition on the interpolation functions, all elements will be considered isoparametric, with the x and y coordinates being replaced by ξ and η that vary between −1 and 1. The relation between the x and y coordinates and ξ and η ones is given by equation (3), where [N] is the matrix of interpolation functions. The interpolation functions are now dependent on the isoparametric coordinates and are the same for all elements (Pires et al., 2021).

(3){u}e(ξ,η)=[N]e⋅{u}e(x,y)

The isoparametric element implies a change in the integration domain by changing from the element physical domain, Ω_e to the natural domain, A, as shown in equation (4), where (•) is the integrand function and [J] is the Jacobian matrix, which is given by equation (5) (Pires et al., 2021).

(4)∫Ωe(•)dxdy=∫A(•)detJdξdη

(5)J(ξ,η)=∂N1∂ξ∂N2∂ξ⋯∂Nn∂ξ∂N1∂η∂N2∂η⋯∂Nn∂η⋅x1ey1ex2ey2e⋯⋯xneyne

The stiffness matrix and the load vector are given by the integrals represented in the left-hand side of equations (6) and (7), respectively, for the case of plane elasticity. In order to facilitate the calculus, numerical quadrature are applied, so that the integral can be replaced by a summation of the integrand function applied at a number of Gauss integration points, n_GP, multiplied by the respective weight, w_GP, as shown in the right-hand side of equations (6) and (7) (Pires et al., 2021).

(6)[K]e=∫ABTDBdetJdξdη=∑g=1nGPBT(ξg,ηg)DB(ξg,ηg)detJ(ξg,ηg)wg

(7){F}e=∫ANTbdetJdξdη+∫AσNTt¯detJσdξσ=∑g=1nGPNT(ξg,ηg)bdetJ(ξg,ηg)wg+∑g=1nGPNT(ξg,ηg)t¯detJσ(ξg,ηg)wg

To calculate the elements stiffness matrix, applying equation (6), it’s necessary to define the strain matrix, B, the elasticity matrix, D and the number of Gauss integration points. The strain matrix is defined as shown in equation (8). The elasticity matrix is defined as represented in equation (9), considering a plane stress state, as already explained, where E is the material Young’s modulus and υ is its Poisson’s coefficient. Regarding the number of integration points, the choice relies in full integration, which means 4 Gauss points. The choice of full integration instead of reduced integration is due to the fact that shear locking is not a problem in this study case and to prevent the hourglass effect (Pires et al., 2021).

(8)[B]=∂N1∂x∂N2∂x⋯∂Nn∂x∂N1∂y∂N2∂y⋯∂Nn∂y=J−1(ξ,η)⋅∂N1∂ξ∂N2∂ξ⋯∂Nn∂ξ∂N1∂η∂N2∂η⋯∂Nn∂η

(9)[D]=E(1−υ2)⋅1υ0υ10001−υ2

In order to calculate the external load vector by applying equation (7), both b and t¯ have to be known, which represent, respectively, the external forces applied over the element area and surface traction applied to the element sides. Thus, boundary conditions have to be known. Besides providing the necessary data to calculate the external load vector, the boundary conditions are also applied in the global stiffness matrix, obtained after assembling all elements stiffness matrices and load vectors. To apply the boundary conditions, the stiffness value is multiplied by a “big number” (in the order of 1e10 or bigger), at the correspondent nodal direction that should be restricted. This overly excessive stiffness prevents displacements in that direction (Pires et al., 2021). With all matrices and vectors calculated, the nodal displacements can be calculated by implementing equation (2). Following the calculus of the displacements, the next step toward obtaining the nodal stress sum matrix is to calculate the nodal strains. The strains values at each node can be reached by using equation (10), with the strain matrix here being defined as equation (11).

(10){ε}e=[B]e⋅{u}e

(11)[B]=∂N1∂x0∂N2∂x0⋯∂Nn∂x00∂N1∂y0∂N2∂y⋯0∂Nn∂y∂N1∂y∂N1∂x∂N2∂y∂N2∂x⋯∂Nn∂y∂Nn∂x

Having the strain vectors calculated for all nodes, the nodal stresses can be easily obtained by applying equation (12). By summing the principal stress in the x direction and in the y direction, the nodal matrix of stress sum is obtained as desired, and the elastic simulation is concluded. This nodal matrix of stress sum represents the stress sum that the component presents when submitted to a static load with the value of the amplitude of the cyclic load.

(12){σ}e=[D]⋅{ε}e→[σxx+σyy]

2.3 Nodal heat source

In order to transpose the obtained results in the elastic simulation to the next thermal simulation, in between step is needed, step 3 of Figure 1. The stress sum matrix will be transformed in to a heat source, applied at each node, proportional to the nodal stresses obtained previously, with the help of an unknown, K_Q, that correlates the elastic simulation and the thermal one. This relation is presented in equation (13). The value of K_Q can be estimated by doing the reverse process of the one that is explained in this section, for a case that the results are known. In order to insert the cyclic aspect of the load, a sine function is applied to the value of Qv* obtained, function of the load frequency, as shown in equation (14).

(13)[Qv*]=KQ⋅[σxx+σyy]

(14)[Qv]=[Qv*]Δtsin(2π⋅f⋅t)

It should be noted that the K_Q parameter is a computational unknown only, since this model begins with an elastic simulation and a conversion of results to a thermal domain is needed. At an experimental level, the temperature measurements are made directly from a components surface and this variable is not taken into account. Also, this parameter was determined for the flat bar (without any stress concentration) and, as already stated, is used to convert the simulated mechanical stress to the heat used in the thermal simulation. Q_V is a heat per area as such, changing the element size does not change the heat “density” and consequently does not depend on the mesh division.

2.4 Thermal simulation

The next steps on this methodology are toward implementing the varying aspect to the simulation. This is done in a cyclic thermal simulation. The main goal of this thermal simulation is to obtain the nodal matrix of temperature variation with time [T(t)], and it is the fourth step in the methodology sequences presented in Figure 1. Differently from the previous elastic simulation, this one is time-dependent. Thus, the governing equation is different too. The governing equation of heat transfer is given by equation (15), where k is the material’s thermal conductivity, ρ its density and c_p its specific heat capacity at constant pressure. With the FEM approximations, the governing equation comes to equation (16), where {T’(t)} represents the nodal values of the derivative of the temperature with respect to time (Thompson, 2004; Elsevier B.V., 2022a).

(15)∂∂xk∂T∂x+∂∂yk∂T∂y+Qext=ρcp∂T∂t

(16)[K]⋅{T(t)}+[C]{T′(t)}={Q}

Following a similar approach to the elastic simulation, the global system can be solved by resolving equation (16) for all elements and then assembling them together. As the temperature changes with time, there are two steps to obtain the desired results. Firstly, equation (16) is rearranged into equation (17) in order to obtain {T’}, where [C] represents the mass matrix. As one can see, the temperature values at an initial state have to be known, and it will be considered that the ambient temperature, T_∞, is 297.5 K.

(17){T′(t)}={Q}−[K]⋅{T(t)}[C]

As it was done for solving the previous problem, considering isoparametric elements and with the help of numerical quadrature, the stiffness matrix [K] can be calculated with equation (18) and the mass matrix with equation (19). When it comes to the external heat load vector, {Q}, similar to the elastic approach, the boundary conditions have to be known, as can be observed in equation (20). Accounting in {Q} are two parcels: the heat flux to the exterior due to convection, {Q_c}, and a nodal heat source, {Q_v}, presented previously.

(18)[K]e=∫ABTkBdetJdξdη=∑g=1nGPBT(ξg,ηg)kB(ξg,ηg)detJ(ξg,ηg)wg

(19)[C]e=∫ANTρcpdetJdξdη=∑g=1nGPNT(ξg,ηg)ρcpdetJ(ξg,ηg)wg

(20){Q}e=∫ANTQdetJdξdη=∑g=1nGPNT(ξg,ηg)QdetJ(ξg,ηg)wg

Regarding the convection, two types of convection exist and will be considered: natural and forced convection. Both types of convection follow the same equation to perform the heat flux calculus, equation (21), only differing in the convection coefficient, h. When it comes to the natural convection coefficient, h_n, its value will be considered to be 5 W/m².K, as is a normal value for this kind of situation (Elsevier B.V., 2022b). Forced convection could be not considered, simplifying the problem. However, considering there is some forced air movement, especially due to crack opening and closure at high frequency rates, and in order to evaluate this phenomenon effect on the heat transfer and consequently on the obtained stresses, it will be considered.

(21){Qc}=h⋅a⋅(T(t)−T∞)

To obtain the values for the forced convection coefficient, h_f, a thermodynamic calculus is performed, with the used relations presented in equation (22). According to thermodynamics, the forced convection coefficient can be calculated applying the first formula presented, where Nu is the local Nusselt number, k_air is the air’s thermal conductivity and x is the node coordinate. The local Nusselt’s number is given by the second relation presented, where Re is the local Reynolds number, and Pr is the air’s Prandtl number. The Reynolds number is calculated with the help of the last formula presented in equation (22), where the velocity is accounted, as well as the node coordinate and the air’s dynamic viscosity, ν. With the help of these three equations, the forced convection coefficient can be calculated, with the only unknown being the air’s velocity, which is dependent of the analyzed part and loading frequency. The value of this coefficient is also dependent on the temperature. However, as explains the thermoelasticity theory, the temperature changes are so small that the dependency can be neglected. So, to extract the necessary air properties, it was considered a constant 293 K (≈20°C). The values of these properties are presented in Table 1 (Frank et al., 2008).

(22)hf=Nu⋅kairx,Nu=0.0332⋅Pr1/3⋅Re0.5,Re=v⋅xν

With both heat fluxes contributions calculated, the vector of external heat sources can be obtained and so can be the nodal derivative of temperature with respect to time, considering equation (17). Thus, it is necessary to calculate the nodal temperatures at the next time increment. Here, the forward difference method is applied, which is a method used in finite differences approximations. Using Eulerian integration in equation (16), equation (23) is obtained (Thompson, 2004). As it happens in the finite differences’ formulation, this method becomes unstable as Δt increases in value (Thompson, 2004). So, Δt is a major parameter in this simulation. It is vital that this value is not high enough, and considering this, the calculus of this value is presented in equation (24). The maximum number of increments, MNI, is established as 500, since the number of cycles chosen to simulate was 25, 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.

(23){T(t+Δt)}={T(t)}+{T′(t)}⋅Δt

(24)Δt=τ⋅nMNI

2.5 Signal amplitude and phase

Next in the TSA simulation is the signal processing, in order to extract the temperature amplitudes and phase. To help with the processing, a Fast Fourier Transform (FFT) function, one of the most useful tools in signal processing, is applied to the script. The temperature signal amplitude can be calculated with the help of the FFT function and taking in consideration that the vector length is equal to the MNI and that the number of saved points per cycle, ppc, is 20. To calculate the temperature signal phase, the angle function will be added to the script in order to determine the phase angle of the complex vector.

2.6 TSA stress

With the temperature signal amplitude calculated, the TSA governing equation (1) can now be applied. However, this equation is complex and with some considerations it can be simplified. Considering an adiabatic and reversible system and an isotropic material under plane stress, within the elastic range, the previous relation can be simplified to equation (25) (Sharpe, 2008). K_TSA is normally given and its value depends on which material is being analyzed. However, it can also be calculated using the last term of equation (25).

(25)[σxx+σyy]TSA={ΔT}KTSA⋅{T},KTSA=αρ⋅cp

Since this equation is only applicable within the elastic range of the material, a problem presents if this range is surpassed. As the plasticity constitutive law in the FEM theory is very complex and of great difficulty to apply in the created model, another solution was found, consisting in applying a different energy, frequency dependent, only at regions that present plastic behavior. The heat applied at those zones nodes is still given by equation (14), but the value of [Qv*] is now given by equation (26), where f represents the frequency.

(26)[Qv*]=562449⋅f0.1106

The existence of a concentrated stress results in the existence of an area under the plastic conditions/domain. Depending on the scale (zoom) used in the thermal images, this plastic domain could be visible. Several hypothesis were tested that could lead to mechanical stress variations due to the load frequency. From the results obtained, and presented later on the text, none of the initial hypothesis validated the stress variation observed experimentally. Thus, it was hypothesized that the thermal elastic model (equation (1)) is incomplete and it was performed an iterative curve fit between the experimental and the numeric simulations result. From this, equation (26) was obtained.

3. Procedure

3.1 Components

In Silva et al. (2017), two different parts were used to perform TSA experimentally. In order to have a good comparison, the above explained method will be applied to these two parts: (1) a simple bar and (2) compact tension specimen (CT specimen). The analysis of the bar provides a validation of the method, as the pretended results are known. The analysis of the CT specimen pretends to evaluate the TSA around the crack tip and the parts behavior when working with non-ideal load frequencies. The dimensions of both components are illustrated in Figure 2, with the CT specimen following the ASTM E 647-08 standard (Dassault Systemes, 2011). Both parts are considered to be made of aluminum AA2219-T851, just like in the tests performed in Silva et al. (2017). Aluminum properties necessary to the simulations are presented in Table 2.

3.2 Meshes

Taking in consideration the parameters chosen in the methodology, six different meshes were designed, three for each part, with the help of the Gmsh software. The meshes illustrated in Figure 3 (called mesh 1, 2 and 3, accordingly) are relative to the bar, whereby the first two consider a symmetry along the x-axis – half bar considered – and the third one is relative to only a sixth of the bar. The meshes become more refined than the other. Meshes 4, 5 and 6, presented in Figure 4, regard the CT specimen. These three also consider symmetry along the x-axis and become more refined than the other. One of the main advantages of mesh 6, besides having smaller elements than the other ones, is the structured form that presents, contrasting the random type of mesh 4 and 5. Between the three, the main concern was the area around the crack tip, which is the point of interest. It’s noticeable that the element size around the tip is much smaller than for the rest of the specimen, as that area is the most critical. Information about each mesh is presented in Table 3.

3.3 Elastic simulation

In the elastic simulations, besides the material’s properties, presented in Table 2, boundary conditions for each case have to be known. Boundary conditions for both cases are illustrated in Figure 5. Those will be applied by differing different types of nodes:

1. Node type 0 – Free movement and no load applied

2. Node type 1 – Free movement and load applied

3. Node type 2 – Constrained movement in y direction and no load applied

4. Node type 3 – Constrained movement in both directions and no load applied

In both cases a fixed node (node type 3) is necessary to prevent the parts free movement. In the case of the bar, the distributed load is applied equally along the top surface of the component. Regarding the CT specimen, this distributed load assumes a parabola-type of behavior, as can be observed in Figure 5.

3.4 Thermal simulation

Similar to the elastic simulation, the material properties and boundary conditions have to be known. The different conditions were applied by considering different types of nodes, namely:

1. Node type 0 – No convection

2. Node type 1 – Natural convection

3. Node type 2 – Forced Convection

Node type 0 was applied at every interior node, and node type 1 at every boundary, except along the crack in the CT specimen, where node type 2 was considered. As it can be concluded, majorly natural convection is considered, although the parts, when realizing the TSA, are in movement. This is to accentuate the effect of the forced convection near the crack area; thus, this is the only region where forced convection is considered.

It is also necessary to calculate the forced convection coefficient values. As this parameter is only applied at the crack, observe Figure 6 to a further detail on the crack area, as well as to better understand the used variables in the below mentioned equations. In this figure only two nodes, besides the crack tip, are illustrated, in order to explain that both crack length and height change depending on which node the equations are applied. To obtain the forced convection coefficient values, the velocities at each node that belongs to the crack have to be calculated, as function of the load frequency. In order to accomplish this, the crack volume is calculated by using equation (27), where w is the specimen thickness, l is the crack length and d is the height. The crack height can be extracted from the nodal displacement value at that specific node, obtained from the previous elastic simulation. Dividing the crack volume by the period of the frequency applied, the air’s flowrate is obtained, as in equation (28). With the air’s flowrate calculated, if divided by the cross-section area, as explained in equation (29), the nodal velocity, function of frequency, is obtained. Again, to calculate the area, the nodal displacements along the crack are already known due to the elastic simulation.

(27)V=l⋅d2⋅w

(28)V̇=Vτ

(29)v=V̇a=V̇w⋅d

3.5 Coefficients

In order to transpose the elastic simulation results to the thermal one, it is necessary to know the value of a coefficient, K_Q, as represented in equation (13). To know this value, a reverse process was done, since the expected results for the bar are known. An uniform stress distribution of 22.5 MPa is expected, in the case of the experiment performed in (Silva et al., 2017). However, to correlate the stress and the temperature changes, another coefficient is needed, K_TSA, as shown in equation (25). This value depends on the parts’ material, and knowing both components are made of aluminum, as already said, this value is known (Silva et al., 2017). Both coefficient values are presented in Table 4.

4. Results

4.1 Results for the bar sample

The bar part was submitted to the presented TSA simulation. The three meshes showed pretty similar results regarding stress values under adiabatic conditions, so mesh 3, as it was more computationally expensive, was discarded. When it came to the convection effect, mesh 1 showed coarse results, so mesh 2 was the one used. The results obtained using mesh 2 are presented in Figure 7a and b. As the main goal of this work relies on the TSA dependency on frequency, and as it was done in (Silva et al., 2017), the simulation was performed considering different frequencies varying from 2.5 Hz up to 70 Hz. The stresses values obtained were of 22.4921 MPa for all the different loading frequencies simulated.

Besides the simulations, a test was also performed to see what was the effect of the convection on the results, and also to validate the chosen natural convection coefficient. The test results are illustrated in Figure 8a–c. The bar is zoomed in so the difference between the different pictures can be more visible. The case for f = 5 Hz and h = 0 W/m².K is not presented since the values of the stress and of the temperature amplitude do not change with the loading frequency in a adiabatic system, for this simple bar case, meaning that the distribution map is pretty much equal to (Figure 8a).

4.2 CT specimen

After analyzing the bar, the CT specimen was submitted to the simulations. As this case is more complicated than the bar case, mesh comparison is vital. So, the stress profile along the x-axis, resultant from the elastic simulation, was plotted for the three different meshes, presented in the graphic of Figure 9. In the top corner of the figure is presented a plot with a zoom around the crack area, with the crosses representing each node from the different meshes. The crack tip is the node with the bigger stress value for all the cases, as expected, located at x = 18 mm.

Besides the mesh analysis, the forced convection coefficient values need to be calculated. Using the above presented equations and the displacements obtained in the elastic simulation (using mesh 5), these values were calculated and are plotted in the graph of Figure 10a. The determined coefficients were fine-tuned to better represent/mimic the convection effect observed experimentally, in particular at the node of the crack tip. Thus, the values implemented were the ones shown in Figure 10b.

The CT specimen was then submitted to the TSA simulation, with different loading frequencies just like the bar was. As it was quite unclear at the moment what node should be studied, there were two points of focus during the simulations: one at the crack tip (x = 18 mm) and one around the crack (x = 17.5 mm, at an elastic-domain area where results were nearly the same for the different meshes. The uncertainty of the analysis point is due to the fact that at x = 17.5 mm the heat transfer effect is practically null and so the stress results behavior with the increase of the loading frequency is not well represented, considering the experimental results of (Silva et al., 2017). However, at the crack tip, the node is now on the plastic domain and its analysis is much more complex, although the heat transfer effect is present.

Due to its element size around the crack tip, mesh 6 was the one used. Temperature variation for both points of analysis and for different frequencies are presented in Figure 11a–d. The stress distribution map was plotted and it’s illustrated in Figure 12, where the visible discontinuity around the holes is probably due to the mesh distribution around it. The stress values obtained, for both focus points and the different loading frequencies, are presented in Table 5.

As this first simulation was performed, it was noticed that the region around the crack tip surpassed the elastic range of the material. With this being the case, equation (26) was applied at this area. After this change, the simulations were ran again. The obtained stress results were plotted to display the behavior dependent of the loading frequency, illustrated in the graph of Figure (13). A curve fitting was also plotted, with the respective equation and R² coefficient value.

The respective stress values, with also the ones extracted in (Silva et al., 2017) for comparison purposes (first row of the table, the fitting curve obtained in the respective mentioned work), are presented in Table 6. Also, the stress and phase profiles for the different loading frequencies were plotted, along the symmetry axis, presented in Figure 14a and b, with zooms in areas of interest. In the case of the stress profile, the zoom is around the crack tip area. Regarding the phase profile, the zoom is exactly at the crack tip, showing a portion of the crack. Temperature variation and stress distribution maps are not shown since there is no visual difference from the ones presented before.

The final step was to apply the created model to the usage of non-ideal loading frequencies, as it was the main goal of this work. A similar procedure to the one just done was followed. The loading frequencies applied varied from 50 to 2000 Hz, and the plot of the obtained stress and respective loading frequencies is illustrated in the graph of Figure 15. The correspondent values are present in Table 7. As was done in the case before, stress and phase profiles along the symmetry axis were plotted and are shown in Figure 16a and b.

5. Discussion

5.1 Bar

Analyzing the obtained results presented in Figure 7, it can be observed that they corresponded to what was expected. Mesh 2 showed good results and there was no need to apply mesh 3 that is more computationally expensive. Temperature variation, observed in Figure 7a, has the desired behavior and the chosen 20 points per cycle, clearly seen represented by the points along the curve, are enough to properly describe the temperatures behavior. Stress plot, observed in Figure 7b, shows a uniform distribution, which was also expected. Finally, analyzing the stress values, a stress around 22.5 MPa was obtained for each loading frequency tested. Besides the stress value being the one expected, the constant, independent of the loading frequency, behavior also was the one expected, according to (Silva et al., 2017).

Regarding the test realized to acknowledge the natural convection effect on the TSA simulation, it can be seen that in the surfaces of the bar the temperature amplitude signal is different from the inside of the bar, although the changes are of small magnitude. More so, the chosen natural convection coefficient of 5 W/m².K is proven to be an adequate value. Comparing the adiabatic case presented in Figure 8a with the ones where convection is present, it can be clearly seen that with the increasing of the loading frequency, Figure 8b and c, this effect is becoming less accentuated and the system behavior approximates an adiabatic system. This behavior corresponds to what is explained in the TSA theory, which states, as said above, that above a certain loading frequency, the system can be considered adiabatic.

5.2 CT specimen

Different from the bar case, the CT specimen meshes were compared, due to the more complex design. As it can be seen in the graph of Figure 9, the three meshes present very similar results only away from the crack tip. This is because the area around the crack tip (where the colors are less accentuated in Figure 9) represents a FEM mesh singularity, meaning that refining the mesh, i.e. decreasing the element size, will lead to stress results that tend to infinity. Because of this, this area needs particular attention to be studied with the FEM, More so, this leaves an uncertainty on the whereabouts of the study node, since the crack tip belongs to the singularity. The choice relies on mesh 6 to perform the simulations, since it is the one that presents more accuracy when analyzing the crack tip results.

When it comes to first TSA simulation, for the case of 5 Hz of loading frequency, the stress distribution map, Figure 12, it’s what was expected. It can be clearly seen the stress intensity factor around the crack tip, with a heart-shaped distribution around it, as expected. Also, the compression zones, like at the specimen beginning (x = 0 mm), are also visible.

Regarding the plots of the temperature variation with time, the two focus points are studied. Observing Figure 11a–d, the amplitude differences are correct for both crack tip analysis nodes, being that at the crack tip nodes the signal amplitudes are bigger than at the rest of the specimen, proportional to the stress, as it should be. However, analyzing the temperature variation plotted in (Figure 11a), results are not so good, since the temperature behavior seems to follow a constant cycle pattern. Ideally, some amplitude attenuation with time, at the crack tip, should be visible, due to the conduction that occurs from the high stress zones (and higher temperatures) to the low stress ones. This effect should be reduced with the increment of the loading frequency, since the heat has less and less time to transfer. Analyzing the plots (Figure 11b) to (Figure 11d), this effect is present. Thus, the crack tip should be studied at x = 18 mm. The downside of this is, as already explained, this node represents a FEM singularity, meaning that extra care should be taken when analyzing this node.

Another aspect that must be pointed out is the fact that although the conduction effect is noticeable, the signal amplitude values practically don’t change, as it can be concluded by looking at the obtained stress values present in Table 5, where the stress values practically don’t change with the increment of loading frequencies. The temperature attenuation occurs, but its values go down below the ambient temperature, meaning that at some points in time the specimen is compressing at the crack tip, which it shouldn’t. The tests performed in (Silva et al., 2017) prevent that by applying a load that guarantees that the specimen is submitted to traction at all times. Some simulations were ran in order to change the range of the sine used in the thermal simulation (change it so it doesn’t vary from −1 to 1, bur from 0 to 2 or even never be zero) but this problem was not corrected.

In order to try and fix the above-mentioned problem, the methodology was changed at the thermal aspect of the simulation. The heat applied at the nodes that are outside the elastic range of the material is now given by equation (26), as it was explained better in the methodology. The equation used is an attempt to try and solve the independent behavior of frequency that the stress results present, but also to correct the behavior at the plastic zone by trying to eliminate the tendency to infinity that it has. Obviously, other solutions exist but this was the one used. Analyzing the graph of Figure 13, it is clear that the frequency now has an effect on the obtained stresses. Also, the curve fitting is pretty good as it presents a R² value of 1, and the obtained values with this fitting curve and the ones extracted from the simulation are almost the same, as it can be seen it Table 6. The behavior is what was expected and similar to the one obtained in (Silva et al., 2017). The curve fitting obtained in that work is compared to the one present in this simulation, in Table 6. Although the stress values evolutions with the increment of loading frequency are very similar, the magnitudes of the values are not. This is due to difference in pixel and element size. When performing the TSA, a camera with a certain resolution is used. The one used in (Silva et al., 2017) had a 320 × 100 resolution, with all pixels being the same size, different from the mesh used that becomes more refined close to the crack tip. The camera captures radiation coming from each pixel, integrated to the pixel area, when the simulation can be more precise and measure the stress at the precise crack tip node. In order to become more visual, observe Figure 17 to compare element and pixel size, where one square represents one pixel.

The stress and phase profiles, along the symmetry axis, for different loading frequencies, were also plotted. The stress profile, Figure 14a, somewhat corresponds to what was expected. It can be seen that the only difference between theirselves occurs in the area where plasticity happens, different from the ones obtained in (Silva et al., 2017) where a noticeable change in profiles happens also at the crack. This could be due to the lack of the forced convection effect. More so, in the mentioned work, the stress peak moves toward the crack propagation direction, while in the obtained profiles in this simulation the stress peak remains constant at the crack tip node, which is more correct. The change in stress peak location should be because of crack propagation of some sort. Analyzing the phase profiles globally, Figure 14b, a change of 180° at the point where the stresses change signal was expected. However, the values of around 25° and 155° seem a bit excessive since these values should be as close as possible to 0° and 180°, if the system were to be considered adiabatic. Looking at the zoom figure, the forced convection can be detected. Although the h_f coefficient increases with the frequency, it is at the low frequencies that these effect is more noticeable, probably due to the lower stress values at those frequencies.

Finally, the model was submitted to the use of non-ideal loading frequencies, as it was the main goal of this work. The obtained stresses are plotted in the graph of Figure 15 and a new curve fitting was plotted, with also a R² value of 1. Both the equations that relate stresses and loading frequencies show a great correlation, even between themselves, being almost the same. Some kind of stabilization was expected as the frequency rises to high values, since a value in which the heat has no time to transfer should exist. However, that behavior was not observable. The rest of the presented data was as expected, proving that this model can be used with high frequencies.

6. Conclusion

This work had the main goal to create a computational finite element model to perform TSA simulations with the usage of various loading frequencies, considering the materials’ thermal properties. The methodology based on FEM was presented first, with the model being validated by running the simulations for the case of a simple bar. The model was used to simulate a CT specimen, where the results helped to update the model in order to simulate the TSA properly. Finally, loading frequencies out of the ideal range were applied and the model could perform those simulations.

The results for the flat bar showed great similarity to the laboratory tests, including the temperature profiles at the borders. When simulating the CT specimen for different frequencies, the stress amplitude presented various amplitude values, mainly near the crack tip. This was also observed in the laboratory tests published previously. The stresses obtained for the CT specimen revealed a heart-shaped pattern, as predicted in the literature. When comparing the load frequency vs stress curve obtained in this work with the laboratory tests previously published, they are overlapping.

As a final remark, the created model reproduced acceptably the TSA, with the desired conditions. However, there’s still room for improvement and to perfect this model.

Figures

Figure 1

Sequence of procedures followed in this methodology

Figure 2

Dimensions of the bar (left) and the CT specimen (right)

Figure 3

Designed meshes for the bar: mesh 1, mesh 2 and mesh 3, respectively

Figure 4

Designed meshes for the CT specimen: mesh 4, mesh 5 and mesh 6, respectively

Figure 5

Boundary conditions applied to the bar (left) and to the CT specimen (right)

Figure 6

Zoom around the crack area

Figure 7

Bar TSA simulation for f = 20 Hz

Figure 8

Bar signal amplitude distribution maps for different cases

Figure 9

Elastic stress profile comparison for mesh 4, mesh 5 and mesh 6

Figure 10

Forced convection coefficient at the crack

Figure 11

Temperature variation with time at the specimen for different cases

Figure 12

CT specimen stress distribution map for f = 5 Hz

Figure 13

Stress vs loading frequency for ideal frequencies

Figure 14

Profiles along the x-axis (symmetry axis) for different ideal loading frequencies

Figure 15

Stress vs loading frequency for non-ideal loading frequencies

Figure 16

Profiles along the x-axis (symmetry axis) for different non-ideal loading frequencies

Figure 17

Element and pixel size

Table 1

Air properties at 20°C

Thermal	Prandtl	Dynamic
Conductivity (W/m.K)	Number (−)	Viscosity (m2/s)
0.0257	0.713	15.11e−6

Source(s): Frank et al. (2008)

Table 2

Aluminum AA2219-T851 properties

Yield	Young’s	Poisson’s	Thermal	Density	Specific heat
Strength (MPa)	Modulus (GPa)	Coefficient (−)	Conductivity (W/m.K)	(kg/m3)	Capacity (J/kg.K)
352	73.1	0.33	120	2,840	864

Source(s): A. S. M. Aerospace Specification Metals (2022)

Table 3

Information about the meshes

Meshes	Mesh 1	Mesh 2	Mesh 3	Mesh 4	Mesh 5	Mesh 6
Number of nodes	185	1,001	2,821	617	693	5,479
Number of elements	144	900	2,700	556	630	5,244

Table 4

Coefficient values

KQ	KTSA
1758	7.7E–6

Table 5

Stress at two different nodes for different loading frequencies

Frequency (Hz)	2.5	10	20	50	70
Stress (MPa) for x = 17.5 mm	167	167	167	167	167
Stress (MPa) for x = 18 mm	703	703	703	703	703

Table 6

Stresses at the crack tip for different cases and ideal frequencies

Frequency (Hz)	2.5	10	20	50	70
Stress (MPa) for y = 68.911 ⋅ x0.3029	91	138	171	225	250
Stress (MPa) obtained in simulation	354	413	445	493	512
Stress (MPa) for y = 319.88 ⋅ x0.1106	354	413	446	493	512

Table 7

Stresses at the crack tip for different cases and non-ideal frequencies

Frequency (Hz)	50	500	1,000	1,500	2000
Stress (MPa) obtained in simulation	493	636	686	718	742
Stress (MPa) for y = 319.83 ⋅ x0.1106	493	636	687	718	741