Finite Element Analysis of the Deformation of Functionally Graded Plates under Thermomechanical Loads

The first-order shear deformation plate model, accounting for the exact neutral plane position, is exploited to investigate the uncoupled thermomechanical behavior of functionally graded (FG) plates. Functionally graded materials are mainly constructed to operate in high temperature environments. Also, FG plates are used in many applications (such as mechanical, electrical, and magnetic), where an amount of heat may be generated into the FG plate whenever other forms of energy (electrical, magnetic, etc.) are converted into thermal energy. Several simulations are performed to study the behavior of FG plates, subjected to thermomechanical loadings, and focus the attention on the effect of the heat source intensity. Most of the previous studies have considered the midplane neutral one, while the actual position of neutral plane for functionally graded plates is shifted and should be firstly determined. A comparative study is performed to illustrate the effect of considering the neutral plane position.The volume fraction of the two constituent materials of the FG plate is varied smoothly and continuously, as a continuous power function of the material position, along the thickness of the plate.


Introduction
Functionally graded materials (FGMs) are microscopically inhomogeneous composite materials, in which the volume fraction of the two or more materials is varied smoothly and continuously as a continuous function of the material position along one or more dimension of the structure.These materials are mainly constructed to operate in high temperature environments.
In conventional laminated composite structures, homogeneous elastic laminae are bonded together to obtain enhanced mechanical and thermal properties.The main inconvenience of such an assembly is the creation of stress concentration sources along the interfaces and specifically when the structure is exposed to elevated temperatures.This can lead to many deficiencies such as delaminations, matrix cracks, and other damage mechanisms which may result from the abrupt change of the mechanical properties at the interface between the layers.One of the ways to overcome this problem is to use functionally graded materials (FGMs) with continuous material properties variations, which can lead to a continuity of the material properties.The concept of functionally graded material (FGM) was proposed in 1984 by the material scientists in Japan [1].Alieldin et al. [2] suggested three approaches to transform the laminated composite plate, with stepped material properties, to an equivalent functionally graded (FG) plate with a continuous property function across the plate thickness.Such transformations are used to determine the details of a functionally graded plate equivalent to the original laminated one.In addition it may provide an easy and efficient way to investigate the behavior of multilayer composite plates, with direct and less computational efforts.FGMs are usually made of a mixture of ceramic and metals.The ceramic constituent of the material provides a high temperature resistance due to its low thermal conductivity, while the ductile metal constituent, on the other hand, prevents the fracture caused by thermal stress due to high temperature gradient in a very short period of time.
The FGM is suitable for various applications, such as thermal coatings of barrier for ceramic engines, gas turbines, nuclear fusions, optical thin layers, and biomaterial electronics.Cheng and Batra [3] studied thermo-mechanical deformations of the FG plates.Tanigawa et al. [4] studied the linear thermal bending of FGM plate in steadystate condition.They also studied the stress in transient heat conduction with temperature-dependent material properties.Using the first-order shear deformation theory (FSDT), Praveen and Reddy [5] analyzed nonlinear static and dynamic responses of FG ceramic-metal plates in a steady temperature field.They used finite element method.Lanhe [6] used the FSDT and derived equilibrium and stability equations of a moderately thick rectangular plate made of FGM under thermal loads.He assumed that the material properties varied as a power law of thickness.Alibeigloo [7] derived an exact solution for thermoelastic response of functionally graded rectangular plates subjected to thermo-mechanical loads.A finite element analysis of thermoelastic field in a rotating FGM circular disk is studied by Afsar and Go [8].This study focuses on the finite element analysis of thermoelastic field in a thin circular functionally graded material disk subjected to a thermal load and an inertia force due to rotation of the disk.Tung and Duc [9] derived a simple analytical approach to investigate the nonlinear stability of functionally graded plates under mechanical and thermal loads.Equilibrium and compatibility equations for FG plates are derived by using the classical plate theory.The nonlinear behaviors of FGM plates under transverse distribution load are investigated by Singha et al. [10] using a high precision plate bending finite element.Material properties of the plate are assumed to be graded in the thickness direction according to the simple power law distribution.The formulation is developed based on the FSDT considering the exact natural surface position.A high-order control volume finite element method is propose by Chareonsuk and Vessakosol [11] to explore thermal stress analysis for functionally graded materials (FGMs) at steady state with unstructured mesh capability for arbitrary-shaped domain.This formulation, also known as cell-vertex finite volume formulation, is useful for material engineers and scientists in determining the thermal response and thermo deformation in FGM that subjected to thermal and mechanical loads.The heat conduction is considered for thermal analysis, whereas the plane elasticity is considered for stress analysis.Wang and Shen performed a nonlinear bending analysis for a FG plate [12] and also performed a nonlinear vibration, a nonlinear bending, and a postbuckling analyses for sandwich plate with FGM face sheets [13].The two plates are resting on an elastic foundation and subjected to thermal environments.
Functionally graded materials are used in many applications, owing to their stability in high thermal environments.To this aim, many approaches are developed to study the thermoelastic behavior of functionally graded materials.One of these approaches is the finite element analysis of such material type.
In this paper the boundary value problem of the uncoupled thermoelastic behavior of FG plate is formulated and solved.First, the temperature distribution is predicted to be used in the thermoelastic analysis of FG plate.Then, the firstorder shear deformation plate theory is proposed, accounting for the exact neutral plane position, for modeling the functionally graded plates.A parametric study is developed to investigate the effects of different distributions of the material properties on the response of the plates.Moreover, some numerical comparisons are performed to show the lack of accuracy while neglecting the effect of neutral plane position for FG plates.

Mathematical Formulation of the Uncoupled Thermoelasticity System
In this section, a mathematical model is derived for the uncoupled thermoelasticity systems considering the exact neutral plane position.The temperature distribution and the effective material properties of the FGMs are determined firstly and used in the developed model as input data.The effect of neutral plane position is typically neglected in most previous studies, while the position of neutral plane for functionally graded plates must be predetermined.Modifications over the model formulated at [2] are presented to account for the thermal load effects on the finite element model and the neutral plane position effects on the material stiffnesses.

Governing Equations.
Based on the FSDT assumptions, the transverse normals would not remain perpendicular to the midsurface but remain straight after deformation.Thus, the transverse shear strains and consequently the shear stresses are constant throughout the laminate thickness.In practice a convenient shear correction factor, equals to 5/6, is assumed for the analysis of the plates [14,15].So, the displacement fields based on FSDT assumptions are where ( 0 , V 0 ,  0 , 0  , 0  ) are unknown functions to be determined.ℎ 0 denotes the position of the neutral plane (see Figure 1) where for isotropic homogenous plates ℎ 0 = 0.The total strain is the variation of the continuum deformation with respect to its volume, so the linear Green-Lagrange strains components for small deformations and moderate rotations (10 0 − 15 0 ) can be determined from (1) as follow: where   are the total strain components.The total strain components are the sum of the elastic strains   , resulting from the applied mechanical loads and the thermal strains   produced due to temperature change.So, the total strains are given by The total strain components could be divided into {  } bending strain and {  } shear strain components as follows: 0, (, ) + V 0, (, ) 0 , (, ) 0 , (, ) 0 , (, ) + 0 , (, ) ] where { 0  } are the nodal bending strains and { 0  } are the nodal shear strains.The nodal strain components will be determined in the coming sections.
The governing equations for the plate equilibrium are derived based on the principle of minimum total potential energy.So, the total potential energy takes the form where {  } and {  } are given by where where () is the thermal coefficient of expansion and () is the continuum temperature change through the plate thickness.
Based on the concept that the equivalent single-layer theories are built up, that a heterogeneous plate is treated as a statically equivalent, single layer having a complex constitutive behavior, reducing the 3D continuum problem to 2D problem, the equivalent layer of the FG plate can be obtained.By integrating for the plate material properties through the plate thickness the equivalent single-layer material matrix can be determined to be where [  ] and [  ] are the bending and shear material matrices, respectively.These material matrices provide the stress-strain relations for FG plates as follows: For FG plates, the two equivalent material matrices can be written in the form where   are the extensional stiffness components,   are the bending-extensional coupling stiffness components, and   are the bending stiffness components.
Prior to the determination of the material stiffnesses, the location of the neutral plane must be given.Clearly, due to varying young's modulus of the plate, the neutral plane is no longer at the midplane but shifted from the midplane unless for a plate with symmetrical young's modulus [16].
The material stiffnesses can be determined considering the exact neutral plane position to be in the form Note that   () are the equivalent material property stiffnesses as a function of the material thickness direction , where in FG plate, its material properties vary smoothly and continuously over the thickness of the structure (plate).
The equivalent material stiffnesses of isotropic FG plate are where  is the material shear correction factor, () is the effective young's modulus, and  is the effective Poisson's ratio of the material through the plate thickness.The convenient shear correction factor has been assumed to be given by 5/6 in our analysis of FG plates.So, by minimizing the total potential energy (6), the equilibrium equations for the FG plate can be obtained.

Thermal Analysis.
Knowledge of the temperature distribution within a body is basically important in many engineering problems.This information will be highly required in computing the capacity of the heat flow in or out the body.Further, if a body is not free to expand in all the directions, some stresses may be developed inside the body.The magnitude of these thermal stresses will influence dramatically on the design of devices such as boilers, steam turbines, and jet engines.The first step in calculating the thermal stresses is to determine the temperature distribution within the body.
FGMs are primarily used in situations where large temperature gradients are encountered.Also, FG plates are used in many applications (such as mechanical, electrical, and magnetic), where an amount of heat may be generated into the FG plate whenever other forms of energy (electrical, magnetic, etc.) are converted into thermal energy.Within our analysis, constant temperatures are imposed at the ceramic and metal surfaces, but a scalar temperature field is assumed to vary continuously along the  coordinate, such that / = / = 0. Furthermore, a stress free reference temperature  0 = 0 is considered.The one-dimensional differential equation governing the steady state heat conduction in a FGM with heat source strength, by assuming a FG plate (ceramic/metal) with a continuous material properties distribution along its thickness, can be written as [17]   ( ()   ) + q = 0, where () is the thermal conductivity through the plate thickness, () is the temperature distribution through the plate thickness, and q is the strength of a heat source inside the body (rate of heat generated per unit volume).
The temperature distribution along the thickness can be obtained by solving the one-dimensional steady state heat transfer equation with the presence of a heat source (13), subject to the following boundary conditions: by performing integration by rearranging (15) and considering neutral plane position integration after solving integral By imposing one of the previous boundary conditions, the constant in (15) may be evaluated as follows: So, in (17), the temperature distribution through the plate thickness, for a steady state FG plate with heat source of strength q , for any distribution of () is given by If the plate is in a steady state without any heat sources, (13) reduces to the one-dimensional steady state heat transfer equation [18] So, the temperature distribution through the plate thickness for any distribution of (), where in our upcoming thermal analyses a power law distribution (the linear rule of mixtures) of thermal conductivity () is assumed, can be written as where where   and   are the thermal conductivity of the upper and lower surfaces of the FG plate, respectively [19].The grading material parameter  can be any nonnegative real number.In the case of a thermally homogenous plate; that is, when  does not depend on , the temperature distribution is linear through the thickness.Excursions from the linear distribution are obtained by changing the grading parameter .
It is convenient here to mention that, in general, there are many approaches for homogenization of FGMs.The choice of the approach should be based on the gradient of gradation relative to the size of a typical representative volume element.One of such approaches is the approximation approach.The linear rule of mixtures and the modified rule of mixtures by Tamura are convenient methods for estimating the equivalent material properties of the FGMs based on the approximation approach.Previous studies predicted that the linear rule of mixtures cannot reflect the detailed constituent geometry and the microstructure and provides a highly questionable accuracy compared to the modified rule of mixtures.On the other hand, the modified rule of mixtures by Tamura provides a convenient accuracy for a wide range of volume fractions and loading conditions.But, the modified rule of mixtures is restricted to the Young's modulus, so any appropriate averaging method must be used to estimate the other thermomechanical properties.Usually the linear rule of mixtures is being conventionally employed [20].

Position of Neutral Surface.
For the analyses of the flexural behavior of a functionally graded plate, subjected only to transversal applied load, we have to determine the location of the neutral plane before solving the equilibrium equation of the plate.Clearly, due to the varying of Young's modulus of the FG plate through the thickness, the neutral plane is no longer located at the midplane but shifted from it.To determine the position of the neutral plane, we construct a new coordinate system such that the new -axis is placed at the neutral axis, which will be determined in the following (see Figure 1).Then we have [16] where ℎ 0 is the distance of the neutral plane from the midplane of the plate.
In this case, similar to the usual treatment in the FSDT, the axial force for an infinitely wide FG plate subjected to transverse mechanical load is given by By putting the axial force equal to zero to determine the position of the neutral surface, where the position of the neutral plane can be determined by choosing ℎ 0 , such that the axial force at the cross-section vanishes By changing the interval of integral we have Then, assuming  is constant across the thickness The position of neutral plane for FG plates can be determined from ) .
(28) Equation ( 28) provides an applicable way to manage and control the position of neutral plane for FG plates.For design considerations, sometimes we have to adapt the neutral plane position with the required design constraints.Changing the grading continuous function of the FG plate controls the position of the neutral plane.

The Finite Element Model
The displacements and normal rotations at any point into a finite element  may be expressed, in terms of the  nodes of the element, as follow where    is the Lagrange interpolation function at node .The Lagrange interpolation functions for nine-node rectangular element (see Figure 2) are given by in terms of the natural coordinates [15] The nodal bending strain can be written as follows: V 0, (, )  0, (, ) + V 0, (, ) 0 , (, ) 0 , (, ) 0 , (, ) + 0 , (, ) where [  ] is the curvature-displacement matrix, [  ] is the shear strain-displacement matrix, and are the nodal degrees of freedom.
So, the total potential energy can be obtained, and ( 6) can be written as follow: The minimum potential energy principle states that In another form where [  ], [  ] are the element bending and shear stiffness matrices, respectively, defined as and {  }, {} are the element thermal and mechanical load vectors, respectively, defined as So by substitution of ( 31) and (32), into (4) and ( 5) the bending and shear strain vectors can be obtained.The normal stress components can be determined as follow: and the shear stress components are

Numerical Result
In this section we present several numerical simulations in order to assess the behavior of functionally graded plates subjected to thermo-mechanical loads.A simple supported plate is considered for the investigation.The plate is made up of a ceramic material at the top and a metallic at the bottom.The simple power law with different values of  = 0 : 2 is used for the through-the-thickness variation.The plate is subjected to mechanical loadings in addition to a temperature gradient through its thickness.

A Functionally Graded Plate Subjected to a Steady State
Thermomechanical Load.The analysis of FG plates is performed for a combination of materials of type ceramicmetal.The lower plate surface is assumed to be aluminum, while the top surface is assumed to be zirconia.Material properties parameter  = [0,2] is considered.Physical material properties are given in Table 1.An all edges simply  supported square plate, of side  = 0.2 m and thickness ℎ = 0.01 m, is taken for study.The plate is subjected to a uniformly distributed mechanical transverse load on the top surface, in addition to a temperature gradient along the thickness, where the temperature of the ceramic rich top surface is 300 ∘ C and the temperature of the metal rich bottom surface is 20 ∘ C (see Figure 3).To investigate the thermo-mechanical behavior of the plate, the temperature distribution through the FG plate thickness should be firstly determined.Figure 4 shows the steady state temperature distribution given by ( 21), for the aforementioned FG plate.The results are highly consistent with that given at [18].It is obviously noticed that the resulted temperature distribution within FG plates of ceramic and metallic constituents is usually smaller than those of homogeneous plates, either purely ceramic or purely metallic plates.The material parameter  has dominant effect on the thermal behavior of the FG plates, where the temperature distribution depends basically on the equivalent thermal conductivity of the constituents and the temperatures of both the upper and lower surfaces.
To investigate the elastostatic behavior of FG plates, with aluminum and zirconia material constituents, several numerical simulations are performed, for different values of the grading parameter .The following nondimensional parameters have been manipulated throughout the numerical simulations: where  denotes the intensity of the applied mechanical load and   is the Young's modulus of the aluminum bottom face.
Figure 5, shows the nondimensional central deflection of a uniform temperature FG plate, subjected to a mechanical load and free of any thermal excitation.From Figure 5, it can be noticed that the central deflection of the FG plate (of the previous example) decreases noticeably by raising the value of the material index (), because of the increasing of the material rigidity.The rigidity of the material increases by increasing the volume fraction of ceramic, as a consequence of raising the value of ().Table 2 shows the numerical results of a set of numerical simulations with different values of ().The table shows a difference in central deflection while considering and neglecting neutral plane position.The plate provides a higher central deflection when considering neutral plane position rather than the midplane of the plate.The results of Figure 5(a) are consistent with those presented in [18].
Figure 6 shows the central deflection  = /ℎ due to a sequence of mechanical loads for different values of .The nondimensional load  takes values in the interval [0, 12].In addition to the uniformly distributed mechanical transverse load on the top surface, the plate is subjected to a thermal field where the ceramic rich top surface is held at 300 ∘ C and the metal rich bottom surface is held at 20 ∘ C (see Figure 1).One may see that all the plates with intermediate material properties experience intermediate values of deflection.This was expected since the metallic plate is the one with the lowest stiffness and the ceramic plate is the one with the highest stiffness.Also, the figure shows that the contribution of the thermal load is higher than that of the mechanical load, which reduces with increasing of the mechanical load intensity.Table 3 shows the numerical results of the simulations.The results of Figure 6(a) are in excellent agreement with those presented in [18,21].
The variation of the axial stress   at the center of the FG plate and along the thickness direction, for different values of the material parameter , is shown in Figure 7.The stress of power law FGM (P-FGM) plate can be represented as a cubic function of  for material parameter  = 2.The maximum tensile stress along the thickness of the FG plate is located at the bottom edge ( = ℎ/2) and increases as the ratio   /  increases.However, the maximum compressive stress is occurred at the top surface ( = −ℎ/2) and is high for small   /  .For the ratio   /  = 1 in which  the FGM plate becomes a homogenous isotropic plate, the stress distribution is a linear function of , and the maximum stress value is occurred at the top and bottom surfaces of the plate.The neutral plane-where zero axial stress is locatedis shifted from the midplane to the top surface, (Figure 7(b)), by variation of grading parameter , while for isotropic homogenous plate, the zero stress is located at the midplane.The results of Figure 7(a) are consistent with those presented by Croce and Venini [18], where they assumed that neutral plane coincides with the midplane.
Figure 8 illustrates the through-the-thickness maximum axial stress along the center line of the plate on the centroidal axis (midplane).We observe in Figure 8 that the stresses are compressive both at the top and at the bottom surfaces.The profiles of the compressive stress for the graded plates are close to each other, and their magnitude is less than the one of the homogeneous plates.All the results of Figure 8(a) are in agreement with those presented in [18,21].

A Functionally Graded Plate Subjected to a Steady State
Thermomechanical Load and Including a Heat Source.The plate considered at the previous section is used to simulate a FG plate in a steady state with heat source strength (rate of heat generated per unit volume) q = 100 × 10 6 watt/m 3 .In addition to a uniformly distributed mechanical transverse load on the top surface, the plate subjected to a thermal field where the ceramic rich top surface is held at 300 ∘ C and the metal rich bottom surface is held at 20 ∘ C.
Figure 9 shows the temperature distribution for the considered FG plate.It is observed that temperature within the FGM plates of ceramic and metallic constituents is always smaller than that corresponding to a purely ceramic or metallic plate.Figure 4 shows the temperature distribution for a steady state behavior of the FG plate, without heat source, where the temperature distribution for the isotropic material is a linear function of the plate thickness coordinate.However, in Figure 9, the temperature distribution for metal through the plate thickness is nearly a linear function, although it is exactly not a linear function.The reason for the extreme difference of the temperature distributions for isotropic materials (zirconia and aluminum) is due to the extreme difference in the thermal conductivity for the two materials, where the temperature distribution depends on the thermal conductivity of the material.While increasing the value of material parameter , the FG plate becomes more sensitive for the amount of heat delivered from the heat source.Figure 10 and Table 4 show the central deflection  = /ℎ due to a sequence of mechanical loads for different values of .
Figure 11, shows the maximum axial stress   distribution through the plate thickness for different material  parameter .The plate is expected to be a steady state heat conductive plate subjected to an amount of generated heat of strength ( q = 100×10 6 W/m 3 ) and thermo-mechanical loads.Figure 12 shows the axial stress distribution through the plate thickness for a steady state plate in which an amount of heat is generated from heat source ( q = 100 × 10 6 W/m 3 ) and only thermal load.Figure 13 shows the nondimensional deflection of the FG steady state plate versus the heat source strength.Figures 6 and 10 show that the isotropic homogenous materials experience nearly the same deflection, independent on the amount of strength of the heat source.However, for FGMs, it is observed that they are more sensitive for the heat strength change.Also, Figure 13 proves that isotropic material provides a very lack sense for heat strength change.But FGMs provide a changeable deflection by changing the amount of heat source strength.
Figures 11 and 12 show that the stress distribution through the plate thickness depends mainly on the material properties more than any other parameter.The reason for the extreme difference of the stress profiles for isotropic materials (zirconia and aluminum) is due to the extreme difference in the thermal conductivity for the two materials, where the temperature distribution and consequently the thermal strains depend on the thermal conductivity of the material.Figures 11 and 12 show that the ceramic material is subjected to the highest thermal stresses because of its low thermal conductivity.Metal material is subjected to thermal stresses that are less than the ceramic material due to its high thermal conductivity, while the FGMs are subjected to thermal stresses that are less than the isotropic materials due to their moderate thermal conductivity.Also, Figures 11 and 12 show the high ability of FGMs to withstand thermal stresses, which reflects its ability to operate at elevated temperatures.

Conclusions
In this study, a finite element model based on the firstorder shear deformation plate (FSDT) theory is developed for the investigation of thermo-mechanical behavior of functionally graded plates.Different numerical simulations have been developed to investigate the thermoelastic behavior of a simply supported FG plate, with different material distributions along the thickness.The numerical results lead to the following conclusions: (1) there is a difference in plate deflection while considering the effect of shifting the neutral plane position; (2) the neutral plane of the FG plate is shifted towards the surface with the higher young's modulus.Also, the position of the neutral plane depends mainly on the ratio of the young's modulus of the two plate constituents; (3) FG plates provide a high ability to withstand thermal stresses, which reflects its ability to operate at elevated temperatures; (4) the FGMs are more sensitive to the variation of the intensity of the heat flow, in or out of the structure, than that may be happened in the case of the isotropic material structures.The FGMs provide a highly stable response for the thermal loading comparing to that of the isotropic materials; (5) due to the continuity of the material properties distribution along the thickness of the plates, the strains and stresses are varied smoothly without any sort of singularities and on contrary to what may be happened in the conventional laminated plates.

Figure 3 : 1 Figure 4 :
Figure 3: A simply supported FG plate subjected to a uniformly distributed mechanical load and thermal loading.

Figure 6 :
Figure 6: Nondimensional center deflection of P-FGM in steady state without heat source present versus nondimensional load intensity () (thermo-mechanical load): (a) neglecting neutral plane position, (b) considering neutral plane position.

Figure 10 :
Figure 10: Nondimensional center deflection of P-FG plate in steady state with heat source present versus nondimensional load intensity () (thermo-mechanical load).

Figure 11 :
Figure 11: Nondimensional axial stress distributions through the thickness of a steady state plate, with heat source present, for thermo-mechanical loads.

Figure 12 :Figure 13 :
Figure 12: Nondimensional axial stress distributions through the thickness of a steady state plate, with heat source present, for thermal loads.