Transient Analysis of a Functionally Graded Ceramic/Metal Layer considering Lord-Shulman Theory

The transient displacement, temperature, and stress fields in a functionally graded ceramic/metal layer under uniform thermal shock conditions at the upper surface are numerically studied based on the Lord-Shulmanmodel, employing a direct finite element method. The Newmark method is employed for the time integration of the problem. A Matlab finite element code is developed for the numerical analysis of the one-dimensional problem under consideration. The Voigt model (rule of mixture) is used for the estimation of the effective properties inside the functionally graded layer and the variation of the volume fraction of the materials follows the sigmoid function in terms of the introduced parameterp. Furthermore, a parametric studywith respect to the parameter p follows, where three different combinations of ceramic/metal materials are considered. It is concluded that the value p = 1, which corresponds to a linear variation of the properties, minimizes themaximum (tensile) stress applied at themiddle of the functionally graded layer.


Introduction
Ceramic materials are used as thermal barrier coatings for the thermal protection of metals in high-temperature environments. A common failure mechanism of those composite configurations is the spallation of the ceramic coating close to the interface with the metal substrate, mainly due to their thermal expansion mismatch [1]. A possible solution for this problem is the use of functionally graded materials (FGMs), which are advanced materials with gradually varying properties [2]. In FG thermal barrier coatings, an intermediate FGM layer connects the ceramic coating and the metal substrate. The thermomechanical properties of the FGM layer vary from the properties of the ceramic material to the properties of the metal material in a continuous way, thus eliminating the material discontinuities.
On the other hand, the classical theory of thermoelasticity predicts that the thermal disturbances propagate through a solid medium with infinite speed. During the last 50 years, generalized theories of thermoelasticity have been formulated that predict finite-speed thermal waves and overcome this physical paradox. The wave type heat propagation is frequently described as second sound. Although the effects of second sound are short-lived, they become important in thermal shock applications [3].
Ceramic/metal FG thermal barrier coatings are used as parts in machines that operate in high-temperature environments. Combustion chambers, exhaust pipes, power generators, aircraft engines, and space shuttles are typical examples of such machines. In these applications, the study of the thermomechanical response of a ceramic/metal FGM layer under thermal shock conditions using generalized thermoelasticity theories is of great importance.
The first generalized thermoelasticity theory was formulated by Lord and Shulman [4]. The Lord-Shulman theory modifies the classical Fourier's law of heat conduction by introducing a relaxation time, which can be interpreted as the time required to establish steady-state heat conduction in a volume element when a temperature gradient is suddenly imposed on that element [3,4]. Another important generalized theory is the Green-Lindsay theory [5], which introduces two relaxation times to modify the stress-strain relations and 2 Mathematical Problems in Engineering the entropy equation, respectively. Furthermore, Green and Naghdi proposed three models of thermoelasticity [6], which are labelled as models I, II, and III. The Green-Naghdi model II in particular assumes no dissipation of thermal energy. Bagri and Eslami [7] suggested a unified generalized thermoelasticity theory, which incorporates the theory of classical coupled thermoelasticity and the generalized thermoelasticity theories of Lord-Shulman, Green-Lindsay, and Green-Naghdi as special cases.
There are several studies available that analyze the stresses developed in FGMs by employing generalized theories of thermoelasticity, which regard different structural configurations such as disks [8,9], cylinders [10][11][12][13], and spheres [14,15]. Researches regarding generalized thermoelasticity analyses in FGM layered structures are rarer in literature. Bagri et al. [16] studied the problem of a FGM layer under thermal shock conditions in the context of Lord-Shulman theory, where the variation of the thermomechanical properties in the FGM layer followed a power law function. They used a transfinite element method to solve the coupled system of equations, where time is eliminated using the Laplace transform. The inverse Laplace transform was based on a numerical scheme. Youssef and El-Bary [17] studied the problem of a symmetric 3-layered strip under thermal shock conditions, in the context of Lord-Shulman theory. Each layer was made of homogeneous, isotropic, and thermoelastic material, where they also assumed that the thermal conductivity was a function of temperature. They also used the Laplace transform to solve the coupled system of equations. Hosseini Zad et al. [18] used a direct finite element method to study the problem of a 2-layered one-dimensional media under uniform thermomechanical loading at the free surfaces. Each layer consisted of homogeneous, isotropic and thermoelastic material. Nowruzpour Mehrian et al. [19] studied the dynamic thermoelastic response of a functionally graded plate subjected to thermal shock based on the Lord-Shulman theory, using the state space approach and the Laplace transform. Heydarpour and Aghdam [20] employed the Lord-Shulman theory to study the transient thermoelastic behavior of rotating functionally graded conical shells.
In most of the aforementioned studies the Laplace transform is employed for the solution of the coupled system of partial differential equations. However, especially in the case of FGMs, the inverse Laplace transform is based on numerical methods, which have numerical stability issues in problems with fast transient loading such as thermal shock applications [21]. Furthermore, in almost all of the above studies the variation of the properties in the FGM layer follows a power law function [8][9][10][11][12][13][14][15][16].
In the present paper the standard Galerkin finite element method [22] is used to study the transient fields of displacement, temperature, and stress in a ceramic/metal FGM layer under uniform thermal shock loading, based on the generalized theory of Lord-Shulman. The time marching scheme employed is based on the Newmark method [22]. The effective properties of the FGM layer are estimated based on the Voigt model (rule of mixture) [23] and the through thickness volume fraction of the materials is based on a sigmoid function in terms of parameter [24]. Moreover, a parametric study with respect to the parameter follows, where three different combinations of ceramic/metal materials are considered and the minimization of the maximum stress applied at the middle of the FGM is used as the optimization criterion. The properties of the materials are assumed to be temperature-independent. A Matlab finite element code [25] is developed for the numerical analysis of the problem under consideration.

Governing Equations
Consider a ceramic/metal FGM layer of total thickness , as shown in Figure 1, which is subjected to uniform thermomechanical loading at the upper and the lower surfaces. The other two dimensions of the layer are assumed to be infinite compared to its thickness, and the properties of the FGM layer depend only on the depth , where 0 ≤ ≤ . At = 0 the FGM is fully ceramic, while at = the FGM is fully metal. The material is considered to be isotropic and linear thermoelastic. Initially, the layer is stress-free, undeformed, and at uniform temperature 0 .
Due to the geometry and the loading conditions the problem becomes one-dimensional. The governing equations of the Lord-Shulman theory are the equations of motion and conservation of energy, which in the absence of body forces and heat sources take the following form [3]: where = ( , ) is the displacement, = ( , ) is the temperature change with respect to the reference temperature 0 , = ( ) and = ( ) are the Lamé constants, = ( ) is the thermal constant, = ( ) is the density, = ( ) is the specific heat capacity, 0 = 0 ( ) is the relaxation time of the Lord-Shulman theory, denotes time, and a dot above a variable indicates differentiation with respect to time.
The normal stress and the heat flow across the thickness of the layer are expressed as Mathematical Problems in Engineering 3 For normalization purposes the following dimensionless variables are introduced: where is a characteristic length, V is a characteristic speed, is the normalized speed of the thermal wave, = ( ) is a dimensionless coupling parameter, and the index indicates reference to the metal material. Hence (1) take the following normalized form: where ∼ is the normalized speed of the elastic dilatational wave in the metal material, ∼ /( ∼ 0 ) 1/2 is the normalized speed of the thermal wave in the metal material, ∼ 0 is the normalized relaxation time of the metal material, and a dot above a variable indicates differentiation with respect to normalized time ∼ . Furthermore, the normalized forms for the stress ∼ and the heat flow ∼ are given by

Finite Element Method
In order to solve numerically (4), the Galerkin finite element method is applied [22]. Linear Lagrange shape functions are used for both the displacement and the temperature variables. By integrating (4) in the domain of the th finite element, the variational formulation of the problem is obtained: where 0 ≤ ∼ ≤ 1 is the natural normalized coordinate, ℎ ∼ is the normalized length of the th finite element, is the vector of the linear Lagrange shape functions, and u ∼ , ∼ , ∼ , and q ∼ are the normalized vectors of displacement, temperature change, normal stress, and heat flow, respectively, at the nodes of the th finite element. The Gauss method [22] is employed for the calculation of the integrals that appear in (6).

Mathematical Problems in Engineering
Equations (6) can also be written in matrix form as where Summing up (7) for all the finite elements and adding up the boundary terms lead to the construction of the global matrix equation, which can be written as where M, C, and K are the global matrices, F is the global force vector which contains only boundary terms, and Λ is the unknown vector that incorporates the displacement and the temperature change at each node. A finite element code is developed in Matlab for the numerical solution of the problem, based on (9), where the Newmark method [22] is employed for the time integration. If the parameters of the Newmark method are appropriately chosen, then the method is unconditionally stable.

Applications
In the following applications the characteristic length and the characteristic speed are selected so that the parameters ∼ (normalized speed of the elastic dilatational wave in the metal material) and ∼ are set to unity. Hence it is obtained:

Verification of the Finite Element Code.
Due to the absence of analytical results for the Lord-Shulman theory, the finite element code shall be verified based on known numerical data in the literature. As a verification example, the results for a homogeneous, isotropic and thermoelastic layer which is subjected to thermal shock loading at the upper surface are compared to the numerical results derived by Bagri et al. [26]. The normalized form of the initial-boundary value problem (IBVP) considered is described by the following equations: Mathematical Problems in Engineering where ∼ = / is the normalized thickness of the layer and the material index has been dropped since the layer is homogeneous.
The normalized parameters considered in [26] are ∼ = 1, ∼ = 1, = 0.05, and ∼ 0 = 1. In order to get the results for the homogeneous layer, the properties of the metal and the ceramic materials are simply set equal to each other. Figures  2-4 show the variations of temperature change, displacement, and stress, respectively, at several times, including some numerical data based on results obtained in [26]. As it is shown in these figures, the results of the present work are in close agreement with the numerical results obtained in [26]. Moreover, the results in this paper are numerically more stable compared to the results in [26]. In this analysis, 1000 linear finite elements and 1250 time steps were used.

A Zirconia/Titanium Alloy FGM Layer under Thermal
Shock Conditions. In the present example, a zirconia/titanium alloy (ZrO 2 /Ti-6Al-4V) FGM layer is considered. The FGM layer is subjected to thermal shock loading at the upper surface, while the temperature of the lower surface is kept steady. The normalized form of the IBVP under consideration is described by the following equations: The normalized thickness of the FGM layer is ∼ = 1 and the initial uniform temperature is 0 = 300 K. The thermomechanical properties of the materials at temperature = 300 K are shown in Table 1 [27,28]. Since in Table 1 the modulus elasticity , the Poisson coefficient V, and the thermal expansion coefficient are given, it is mentioned that the Lamé constants and and the thermal constant are given by [3] = ] (1 + ]) (1 − 2]) , The volume fraction of the metal material in the FGM layer is assumed to follow the sigmoid function in terms of parameter [24]. Based on the Voigt model (rule of mixture) [23], a property = ( ∼ ) in the FGM layer is estimated by the following equation: where and are the relevant properties of the ceramic and the metal material, respectively, and > 0 is the parameter of the sigmoid function. For example, the variation of the thermal conductivity in the FGM layer for several values of is shown in Figure 5. The sigmoid function has a basic characteristic: regardless of the value of the parameter , the total amounts of the ceramic and the metal materials in the FGM layer are equal. As → 0 + the FGM layer tends to become a homogeneous layer with properties equal to the mean value of the ceramic and the metal properties, while as → ∞ the FGM layer tends to become a composite layer consisting of a fully ceramic layer and a fully metal layer. In other words, as increases the ceramic material tends to concentrate at the area 0 ≤ ∼ ≤ 0.5.
According to (3) and the values in Table 1, the normalized parameters for zirconia are ∼ = 1.11, ∼ = 0.49, and = 0.005 and for titanium alloy is ∼ = 1.00, ∼ = 1.00, and = 0.005. Furthermore, the normalized values for the relaxation times are assumed to be ∼ 0 = 1.5625 for zirconia and ∼ 0 = 0.64 for titanium alloy [9].  The elastic wave speed is ∼ = 1.11 in the ceramic material and ∼ = 1.00 in the metal material. Due to the similar elastic wave speeds for the two materials, the speed of the elastic wave does not change significantly as it propagates in the FGM layer ( Figure 6). On the other hand, the thermal wave speed is ∼ /( ∼ 0 ) 1/2 = 0.40 in the ceramic material and ∼ /( 0 ) 1/2 = 1.25 in the metal material. As it is shown in Figure 7, the speed of the thermal wave increases for larger values of ∼ due to the increase of the volume fraction of the metal material. Time ∼ = 1.56 shows the reflection of the thermal wave at the boundary at ∼ = 1. The propagation of the stress wave is shown in Figure 8. Initially, the stress wave is compressive and it consists of two wavefronts. The slower wavefront is related to the thermal wave and the faster wavefront is related to the elastic wave. Moreover, it is observed that the magnitude of the stress wave increases for larger values of ∼ due to the increase of the volume fraction of the metal constituent.

Dependence on the Mesh Density.
In order to ensure that the numerical results are correct, in this section the dependence of the analysis to the mesh density is investigated. It is anticipated that the mesh should be more dense for higher values of , since higher values of lead to more rapid change of the material properties inside the FGM layer (see also Figure 5). Figures 9-11 show the fields of displacement, temperature change, and stress at time ∼ = 1.00, respectively, for = 1 with variable density of the mesh. In these figures symbolizes the number of the finite elements and symbolizes the number of time steps used. As it is expected, the results become more accurate as the number of finite elements and time steps used increases. Nevertheless, the dependency of the results on the density of the mesh is rather small and the extracted fields can be considered as accurate even when the mesh used is sparse. From Figures 12-14, where the value = 10 is considered, it is implied that this holds true even for higher values of .
Consequently, when the total normalized time is ∼ = 1.00, it is concluded that using 2000 linear finite elements and 2000 time steps for the numerical analysis would lead to accurate results, even for large values of .

Effect of the Parameter of the Sigmoid Function.
The effect of the variation of the thermomechanical properties in the FGM layer is examined in the following analysis. The same problem is considered, as described by (12). Figures  15-17 show the time evolution of displacement, temperature, and stress, respectively, at ∼ = 0.5 for = 1/10, = 1, and = 10. In this analysis and for each value of , 2000 linear finite elements and 12000 time steps were used, since the total normalized time of the analysis is ∼ = 6.00.
From Figures 15 and 16 it is concluded that for smaller values of the propagation speed of the thermoelastic disturbances increases. Furthermore, in Figure 16 the dissipation of the thermal wave is clearly shown, as the temperature quickly stabilizes at its final value. The final temperature at ∼ = 0.5 decreases for larger values of due to the increase of the volume fraction of the ceramic material in the area 0 ≤ ∼ ≤ 0.5, which leads to the decrease of the thermal conductivity and consequently to the increase of the temperature gradient in the same area. It is also observed that for larger values of the maximum temperature and the maximum displacement at ∼ = 0.5 decrease. Finally, from Figure 17 it is deduced that the value = 1 leads to significant decrease of the magnitude of the stress wave. This is an important result, since the minimization of stresses can be used as an optimization criterion for the selection of the most appropriate value for .
For the derivation of more general results, two more combinations of ceramic/metal materials are examined: silicon nitride/stainless steel (Si 3 NO 4 /SUS304) and alumina/nickel (Al 2 O 3 /Ni). The properties of those materials at temperature = 300 K are shown in Tables 2 [27,28] and 3 [29]. Table 4 shows the normalized speeds of the thermal wave and the elastic wave in the materials of the three ceramic/metal combinations considered. Generally, the elastic wave propagates faster in the ceramic material, while the thermal wave propagates faster in the metal material.
The maximum values of the normalized temperature change and normalized stress for the three different ceramic/ metal material combinations considered and for several values  of are shown in Table 5. Based on the results presented in Table 5, there are two basic conclusions. The first conclusion is that the maximum temperature at ∼ = 0.5 decreases for larger values of . This is related to the increase of the volume fraction of the ceramic material in the area 0 ≤ ∼ ≤ 0.5, as it has already been discussed. The combination of Si 3 NO 4 /SUS304 is an exception to this result, due to   the fact that the difference of the thermal wave speeds in the ceramic and the metal materials of this combination is small ( Table 4). The second conclusion regards the maximum tensile stress at ∼ = 0.5, which is minimized for = 1 for all three combinations of materials. This result suggests that, in the context of the Lord-Shulman theory, the stresses in a ceramic/metal FGM layer subjected to the boundary conditions described in (12) are minimized for = 1. Nevertheless, this conclusion is valid only for the assumed normalized relaxation times of the materials and cannot be generalized.

Conclusions
In the present paper, the transient fields of displacement, temperature, and stress in a ceramic/metal FGM layer under thermal shock loading are studied in the context of the generalized Lord-Shulman theory. The numerical analysis is based on a direct finite element method. The Newmark method is employed for the time integration. The verification example demonstrates that the present numerical method produces reliable results for thermal shock applications, while other methods used in the literature (such as the Laplace transform) have numerical stability issues. Furthermore, the properties of the FGM layer are assumed to follow the sigmoid function in terms of parameter , while in the literature a power law distribution is usually used.  Initially, the method is applied to a homogeneous layer and the results are found to be in agreement with those reported in the literature. Then, a zirconia/titanium alloy FGM layer with linear variation of the properties ( = 1) is considered and the results are presented. Finally, a parametric analysis with regard to the parameter follows, where the minimization of the transient stress applied at the middle of the FGM layer is used as an optimization criterion. Three different combinations of ceramic/metal materials are considered for this analysis, namely, zirconia/titanium alloy, silicon nitride/stainless steel, and alumina/nickel. It is concluded that for larger values of the temperature at the middle of the FGM layer generally decreases due to the increase of the volume fraction of the ceramic material in the area 0 ≤ ∼ ≤ 0.5. On the other hand, the stress applied at the middle of the FGM layer is minimized for = 1. It should be mentioned that these results are valid for the specific boundary conditions considered and for the assumed values for the relaxation times of the materials.

Conflicts of Interest
The authors declare that they have no conflicts of interest.