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

Authors: A. E. Alshorbagy; S. S. Alieldin; M. Shaat; F. F. Mahmoud
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

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.


---

## Body

## 1. 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 first-order 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.


## 2. 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.


### 2.1. 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
(1)u(x,y,z)=u0(x,y)+(z+h0)∅x(x,y),v(x,y,z)=v0(x,y)+(z+h0)∅y(x,y),ω(x,y,z)=ω0(x,y),
where (u0,v0,ω0,∅x,∅y) are unknown functions to be determined. h0 denotes the position of the neutral plane (see Figure 1) where for isotropic homogenous plates h0=0.Figure 1
Coordinate system used for a typical FG plate.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 (100-150) can be determined from (1) as follow:
(2)exx=u0,x+(z+h0)∅x,x,eyy=v0,y+(z+h0)∅y,y,γxy=u0,y+v0,x+(z+h0)(∅x,y+∅y,x),γxz=w0,x+∅x,γyz=w0,y+∅y,
where eij are the total strain components. The total strain components are the sum of the elastic strains εij, resulting from the applied mechanical loads and the thermal strains εTi produced due to temperature change. So, the total strains are given by
(3)eij=εij+εTi.
The total strain components could be divided into {eb} bending strain and {es} shear strain components as follows:
(4){exx(x,y,z)eyy(x,y,z)γxy(x,y,z)}=[100(z+h0)000100(z+h0)000100(z+h0)]×[u0,x(x,y)v0,y(x,y)u0,y(x,y)+v0,x(x,y)∅x,x(x,y)∅y,y(x,y)∅x,y(x,y)+∅y,x(x,y)]or{eb}=[Zb]{eb0},(5){γxz(x,y,z)γyz(x,y,z)}=[1001][ω0,x(x,y)+∅x(x,y)ω0,y(x,y)+∅y(x,y)]or{es}=[Zs]{es0},
where {eb0} are the nodal bending strains and {es0} 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(6)Π=(0.5∫A{eb0}T[DEb]{eb0}dA-∫A{eb0}T{DTb}dA)+(0.5∫A{es0}T[DEs]{es0}dA-∫A{es0}T{DTs}dA)-∑{P}{uo},
where {DTb} and {DTs} are given by
(7){DTb}=∫z[Zb]T[Db]{εTb}dz,{DTs}=∫z[Zs]T[Ds]{εTs}dz,where{εTb}={α(z)δT(z)α(z)δT(z)0},{εTs}={00},
where α(z) is the thermal coefficient of expansion and δT(z) 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(8)[DEb]=∫-h/2h/2[Zb]T[Db]T[Zb]dz,[DEs]=∫-h/2h/2[Zs]T[Ds]T[Zs]dz,
where [Db] and [Ds] are the bending and shear material matrices, respectively. These material matrices provide the stress-strain relations for FG plates as follows:
(9)[Db]=[Q-11Q-12Q-16Q-12Q-22Q-26Q-16Q-26Q-66],[Ds]=[Q-44Q-45Q-45Q-55].
For FG plates, the two equivalent material matrices can be written in the form
(10)[DEb]=[AijBijBijDij],[DEs]=[000Sij],
where Aij are the extensional stiffness components, Bij are the bending-extensional coupling stiffness components, and Dij 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
(11)Aij=∫-h/2h/2Q-ij(z)dzi,j=1,2,6,Bij=∫-h/2h/2(z+h0)Q-ij(z)dzi,j=1,2,6,Dij=∫-h/2h/2(z+h0)2Q-ij(z)dzi,j=1,2,6,Sij=Aij=∫-h/2h/2Q-ij(z)dzi,j=4,5.
Note that Q-ij(z) are the equivalent material property stiffnesses as a function of the material thickness direction z, 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(12)Q-11(z)=Q-22(z)=E(z)1-υ2,Q-12(z)=υQ-11(z),Q-66(z)=1-υ2Q-11(z),Q-44(z)=Q-55(z)=K1-υ2Q-11(z),Q-16(z)=Q-26(z)=Q-45(z)=0,
where K is the material shear correction factor, E(z) 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.


### 2.2. 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 thez coordinate, such that ∂T/∂x=∂T/∂y=0. Furthermore, a stress free reference temperature T0=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]
(13)∂∂z(k(z)∂T∂z)+q˙=0,
where k(z) is the thermal conductivity through the plate thickness, T(z) 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:
(14)T=TUatz=-h2,T=TLatz=h2,
by performing integration
(15)k(z)dTdz+q˙z=C1
by rearranging (15) and considering neutral plane position integration
(16)∫T(z)TUdT=∫zh/2C1-q˙(z+h0)k(z)dz
after solving integral
(17)TU-T(z)=C1∫zh/2dzk(z)-∫zh/2q˙(z+h0)dzk(z).
By imposing one of the previous boundary conditions, the constant in (15) may be evaluated as follows:
(18)C1=TU-TL∫-h/2h/2(dz/k(z))+q˙∫-h/2h/2(z+h0)(dz/k(z))∫-h/2h/2(dz/k(z)).
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 k(z) is given by
(19)T(z)=TU-C1∫zh/2dξk(ξ)+∫zh/2q˙zdξk(ξ)+h0q˙∫zh/2dξk(ξ).If the plate is in a steady state without any heat sources, (13) reduces to the one-dimensional steady state heat transfer equation [18]
(20)∂∂z(γ(z)∂T∂z)=0.
So, the temperature distribution through the plate thickness for any distribution of k(z), where in our upcoming thermal analyses a power law distribution (the linear rule of mixtures) of thermal conductivity k(z) is assumed, can be written as
(21)T(z)=TU-TU-TL∫-h/2h/2(dz/k(z))∫zh/2dξk(ξ),
where
(22)k(z)=kU+(kL-kU)(z+(h/2)h)n,
where kU and kL are the thermal conductivity of the upper and lower surfaces of the FG plate, respectively [19]. The grading material parameter n can be any nonnegative real number. In the case of a thermally homogenous plate; that is, when k does not depend on z, the temperature distribution is linear through the thickness. Excursions from the linear distribution are obtained by changing the grading parameter n.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 thermo-mechanical properties. Usually the linear rule of mixtures is being conventionally employed [20].


### 2.3. 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 newx-axis is placed at the neutral axis, which will be determined in the following (see Figure 1). Then we have [16]
(23)x=x1,z=z1-h0,
where h0 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(24)Nxx=B11φx,x.
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 h0, such that the axial force at the cross-section vanishes
(25)B11=∫-(h/2)+h0(h/2)+h0z1E(z1)1-υ2dz1=0.
By changing the interval of integral we have
(26)∫-h/2h/2(z+h0)E(z)1-υ2dz=0.
Then, assuming υ is constant across the thickness
(27)(∫-h/2h/2zE(z)dz)+(h0∫-h/2h/2E(z)dz)=0.
The position of neutral plane for FG plates can be determined from
(28)h0=-(∫-h/2h/2zE(z)dz∫-h/2h/2E(z)dz).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.


## 2.1. 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
(1)u(x,y,z)=u0(x,y)+(z+h0)∅x(x,y),v(x,y,z)=v0(x,y)+(z+h0)∅y(x,y),ω(x,y,z)=ω0(x,y),
where (u0,v0,ω0,∅x,∅y) are unknown functions to be determined. h0 denotes the position of the neutral plane (see Figure 1) where for isotropic homogenous plates h0=0.Figure 1
Coordinate system used for a typical FG plate.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 (100-150) can be determined from (1) as follow:
(2)exx=u0,x+(z+h0)∅x,x,eyy=v0,y+(z+h0)∅y,y,γxy=u0,y+v0,x+(z+h0)(∅x,y+∅y,x),γxz=w0,x+∅x,γyz=w0,y+∅y,
where eij are the total strain components. The total strain components are the sum of the elastic strains εij, resulting from the applied mechanical loads and the thermal strains εTi produced due to temperature change. So, the total strains are given by
(3)eij=εij+εTi.
The total strain components could be divided into {eb} bending strain and {es} shear strain components as follows:
(4){exx(x,y,z)eyy(x,y,z)γxy(x,y,z)}=[100(z+h0)000100(z+h0)000100(z+h0)]×[u0,x(x,y)v0,y(x,y)u0,y(x,y)+v0,x(x,y)∅x,x(x,y)∅y,y(x,y)∅x,y(x,y)+∅y,x(x,y)]or{eb}=[Zb]{eb0},(5){γxz(x,y,z)γyz(x,y,z)}=[1001][ω0,x(x,y)+∅x(x,y)ω0,y(x,y)+∅y(x,y)]or{es}=[Zs]{es0},
where {eb0} are the nodal bending strains and {es0} 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(6)Π=(0.5∫A{eb0}T[DEb]{eb0}dA-∫A{eb0}T{DTb}dA)+(0.5∫A{es0}T[DEs]{es0}dA-∫A{es0}T{DTs}dA)-∑{P}{uo},
where {DTb} and {DTs} are given by
(7){DTb}=∫z[Zb]T[Db]{εTb}dz,{DTs}=∫z[Zs]T[Ds]{εTs}dz,where{εTb}={α(z)δT(z)α(z)δT(z)0},{εTs}={00},
where α(z) is the thermal coefficient of expansion and δT(z) 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(8)[DEb]=∫-h/2h/2[Zb]T[Db]T[Zb]dz,[DEs]=∫-h/2h/2[Zs]T[Ds]T[Zs]dz,
where [Db] and [Ds] are the bending and shear material matrices, respectively. These material matrices provide the stress-strain relations for FG plates as follows:
(9)[Db]=[Q-11Q-12Q-16Q-12Q-22Q-26Q-16Q-26Q-66],[Ds]=[Q-44Q-45Q-45Q-55].
For FG plates, the two equivalent material matrices can be written in the form
(10)[DEb]=[AijBijBijDij],[DEs]=[000Sij],
where Aij are the extensional stiffness components, Bij are the bending-extensional coupling stiffness components, and Dij 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
(11)Aij=∫-h/2h/2Q-ij(z)dzi,j=1,2,6,Bij=∫-h/2h/2(z+h0)Q-ij(z)dzi,j=1,2,6,Dij=∫-h/2h/2(z+h0)2Q-ij(z)dzi,j=1,2,6,Sij=Aij=∫-h/2h/2Q-ij(z)dzi,j=4,5.
Note that Q-ij(z) are the equivalent material property stiffnesses as a function of the material thickness direction z, 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(12)Q-11(z)=Q-22(z)=E(z)1-υ2,Q-12(z)=υQ-11(z),Q-66(z)=1-υ2Q-11(z),Q-44(z)=Q-55(z)=K1-υ2Q-11(z),Q-16(z)=Q-26(z)=Q-45(z)=0,
where K is the material shear correction factor, E(z) 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.


## 2.2. 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 thez coordinate, such that ∂T/∂x=∂T/∂y=0. Furthermore, a stress free reference temperature T0=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]
(13)∂∂z(k(z)∂T∂z)+q˙=0,
where k(z) is the thermal conductivity through the plate thickness, T(z) 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:
(14)T=TUatz=-h2,T=TLatz=h2,
by performing integration
(15)k(z)dTdz+q˙z=C1
by rearranging (15) and considering neutral plane position integration
(16)∫T(z)TUdT=∫zh/2C1-q˙(z+h0)k(z)dz
after solving integral
(17)TU-T(z)=C1∫zh/2dzk(z)-∫zh/2q˙(z+h0)dzk(z).
By imposing one of the previous boundary conditions, the constant in (15) may be evaluated as follows:
(18)C1=TU-TL∫-h/2h/2(dz/k(z))+q˙∫-h/2h/2(z+h0)(dz/k(z))∫-h/2h/2(dz/k(z)).
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 k(z) is given by
(19)T(z)=TU-C1∫zh/2dξk(ξ)+∫zh/2q˙zdξk(ξ)+h0q˙∫zh/2dξk(ξ).If the plate is in a steady state without any heat sources, (13) reduces to the one-dimensional steady state heat transfer equation [18]
(20)∂∂z(γ(z)∂T∂z)=0.
So, the temperature distribution through the plate thickness for any distribution of k(z), where in our upcoming thermal analyses a power law distribution (the linear rule of mixtures) of thermal conductivity k(z) is assumed, can be written as
(21)T(z)=TU-TU-TL∫-h/2h/2(dz/k(z))∫zh/2dξk(ξ),
where
(22)k(z)=kU+(kL-kU)(z+(h/2)h)n,
where kU and kL are the thermal conductivity of the upper and lower surfaces of the FG plate, respectively [19]. The grading material parameter n can be any nonnegative real number. In the case of a thermally homogenous plate; that is, when k does not depend on z, the temperature distribution is linear through the thickness. Excursions from the linear distribution are obtained by changing the grading parameter n.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 thermo-mechanical properties. Usually the linear rule of mixtures is being conventionally employed [20].


## 2.3. 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 newx-axis is placed at the neutral axis, which will be determined in the following (see Figure 1). Then we have [16]
(23)x=x1,z=z1-h0,
where h0 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(24)Nxx=B11φx,x.
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 h0, such that the axial force at the cross-section vanishes
(25)B11=∫-(h/2)+h0(h/2)+h0z1E(z1)1-υ2dz1=0.
By changing the interval of integral we have
(26)∫-h/2h/2(z+h0)E(z)1-υ2dz=0.
Then, assuming υ is constant across the thickness
(27)(∫-h/2h/2zE(z)dz)+(h0∫-h/2h/2E(z)dz)=0.
The position of neutral plane for FG plates can be determined from
(28)h0=-(∫-h/2h/2zE(z)dz∫-h/2h/2E(z)dz).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.


## 3. The Finite Element Model

The displacements and normal rotations at any point into a finite elemente may be expressed, in terms of the n nodes of the element, as follow
(29){u0(x,y)v0(x,y)ω0(x,y)∅x(x,y)∅y(x,y)}=∑i=1n[ψie00000ψie00000ψie00000ψie00000ψie]{ujvjωjSj1Sj2},
where ψie is the Lagrange interpolation function at node i. The Lagrange interpolation functions for nine-node rectangular element (see Figure 2) are given by in terms of the natural coordinates [15]
(30){ψ1ψ2ψ3ψ4ψ5ψ6ψ7ψ8ψ9}=14[(1-ξ)(1-η)(-ξ-η-1)+(1-ξ2)(1-η2)(1+ξ)(1-η)(ξ-η-1)+(1-ξ2)(1-η2)(1+ξ)(1+η)(ξ+η-1)+(1-ξ2)(1-η2)(1-ξ)(1+η)(-ξ+η-1)+(1-ξ2)(1-η2)2(1-ξ2)(1-η)-(1-ξ2)(1-η2)2(1+ξ)(1-η2)-(1-ξ2)(1-η2)2(1-ξ2)(1+η)-(1-ξ2)(1-η2)2(1-ξ)(1-η2)-(1-ξ2)(1-η2)4(1-ξ2)(1-η2)].
The nodal bending strain can be written as follows:
(31){eb0}={u0,x(x,y)v0,y(x,y)u0,y(x,y)+v0,x(x,y)∅x,x(x,y)∅y,y(x,y)∅x,y(x,y)+∅y,x(x,y)}=∑j=19[ψj,xe00000ψj,ye000ψj,yeψj,xe000000ψj,xe00000ψj,ye000ψj,yeψj,xe]{ujvjωjSj1Sj2}or{eb0}=∑j=1n=9[Bbj]{uj0},(32){es0}={ω0,x(x,y)+∅x(x,y)ω0,y(x,y)+∅y(x,y)}=∑j=19[00ψj,xeψje000ψj,ye0ψje]{ujvjωjSj1Sj2}or{es0}=∑j=1n=9[Bsj]{uj0},
where [Bbj] is the curvature-displacement matrix, [Bsj] is the shear strain-displacement matrix, and {uj0}={ujvjωjSj1Sj2}T are the nodal degrees of freedom.Figure 2
Nine-node quadratic Lagrange rectangular element.So, the total potential energy can be obtained, and (6) can be written as follow:
(33)Π=(0.5∫A{u0}T[Bb]T[DEb][Bb]{u0}dA-∫A{u0}T[Bb]T{DTb}dA)+(0.5∫A{u0}T[Bs]T[DEs][Bs]{u0}dA-∫A{u0}T[Bs]T{DTs}dA)-∑{uo}T{p}.
The minimum potential energy principle states that
(34)δΠ=(∫A[Bb]T[DEb][Bb]dA){u0}+(∫A[Bs]T[DEs][Bs]dA){u0}-∫A{DTb}dA-∫A{DTs}dA-∑i=1n[ψie]T{p}=0.
In another form
(35)[[Kb]+[Ks]]{u0}={FT}+{P},
where [Kb], [Ks] are the element bending and shear stiffness matrices, respectively, defined as
(36)[Kb]=∫A[Bb]T[DEb][Bb]dA,[Ks]=∫A[Bs]T[DEs][Bs]dA
and {FT}, {P} are the element thermal and mechanical load vectors, respectively, defined as
(37){P}=∫A[ψie]T{p}dA,{FTb}=∫A[Bb]T{DTb}dA,{FTs}=∫A[Bs]T{DTs}dA.
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:
(38){σxxσyyσxy}=[Q-11(z)Q-12(z)0Q-12(z)Q-22(z)000Q-66(z)]×{{exxeyyγxy}-{α(z)δT(z)α(z)δT(z)0}}or{σb}=[Db]{εb},
and the shear stress components are
(39){σyzσxz}=[Q-44(z)00Q-55(z)]{{γyzγxz}-{00}}or{σs}=[Ds]{εs}.


## 4. 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 ofn=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.


### 4.1. 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 ceramic-metal. The lower plate surface is assumed to be aluminum, while the top surface is assumed to be zirconia. Material properties parametern=[0,2] is considered. Physical material properties are given in Table 1. An all edges simply supported square plate, of side a=0.2 m and thickness h=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).Table 1
Material properties.
Property
Aluminum
Zirconia
Young’s modulus
E
L = 70 Gpa
E
U = 151 Gpa
Poisson’s ratio
ν
L = 0.3
ν
U = 0.3
Thermal conductivity
k
L
=
204 W/mK
k
U
=
2.09 W/mK
Thermal expansion
α
L
=
23
×
10
-
6
/C
α
U
=
10
×
10
-
6
/CFigure 3
A simply supported FG plate subjected to a uniformly distributed mechanical load and thermal loading.To investigate the thermo-mechanical behavior of the plate, the temperature distribution through the FG plate thickness should be firstly determined. Figure4 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 n 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.Figure 4
Temperature distribution through the FG plate thickness for various values of grading parametern (without heat source).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 parametern. The following nondimensional parameters have been manipulated throughout the numerical simulations:central deflectionω-=ω/h;nondimensional load intensityP-=a4p/ELh4;thicknessh-=h/a,wherep denotes the intensity of the applied mechanical load and EL is the Young’s modulus of the aluminum bottom face.Figure5, 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 (n), 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 (n). Table 2 shows the numerical results of a set of numerical simulations with different values of (n). 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].Table 2
Nondimensional deflection of midpoint versus nondimensional load intensity (P-).
Nondimensional central deflection (ω-)
P
-
Metal
n
=
0.5
n
=
1
n
=
2
Ceramic
Midplane
Neutral plane
Midplane
Neutral plane
Midplane
Neutral plane
1
0.0452
0.0317
0.0328
0.0287
0.03
0.0267
0.0277
0.021
2
0.0905 (0.089)*
0.0633 (0.061)
0.0656
0.0573 (0.056)
0.06
0.0534 (0.052)
0.0554
0.042 (0.0413)
3
0.1357
0.095
0.0985
0.086
0.09
0.0801
0.0831
0.0629
4
0.181 (0.1804)
0.1266 (0.124)
0.1313
0.1146 (0.113)
0.119
0.1067 (0.103)
0.1108
0.0839 (0.083)
5
0.2262
0.1583
0.1641
0.1433
0.149
0.1334
0.1385
0.1049
6
0.2714 (0.269)
0.1899 (0.185)
0.1969
0.172 (0.1695)
0.179
0.1601 (0.154)
0.1662
0.1258 (0.126)
7
0.3167
0.2216
0.2298
0.2006
0.209
0.1868
0.1939
0.1468
8
0.3619 (0.361)
0.2532 (0.248)
0.2626
0.2293 (0.228)
0.239
0.2135 (0.208)
0.2216
0.1678 (0.166)
9
0.4072
0.2849
0.2954
0.2579
0.269
0.2402
0.2493
0.1888
10
0.4524 (0.452)
0.3166 (0.308)
0.3282
0.2866 (0.285)
0.299
0.2669 (0.261)
0.277
0.2097 (0.206)
11
0.4979
0.3482
0.361
0.3153
0.329
0.2935
0.3047
0.2307
12
0.5428 (0.54)
0.3799 (0.365)
0.3939
0.3439 (0.343)
0.359
0.3202 (0.313)
0.3324
0.2517 (0.252)

*The consistent results of Croce and Venini [18].Nondimensional center deflection of P-FGM plate versus nondimensional load intensity (mechanical load): (a) neglecting neutral plane position, (b) considering neutral plane position.
(a)
(b)Figure6 shows the central deflection ω-=ω/h due to a sequence of mechanical loads for different values of n. The nondimensional load P- 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].Table 3
Nondimensional deflection of midpoint versus nondimensional load intensity (P-) (thermo-mechanical load without heat source).
(ω-) Nondimensional central deflection for
P
-
Metal
n
=
1
n
=
2
Ceramic
Midplane
Neutral plane
Midplane
Neutral plane
1
−0.205 (−0.2)*
−0.0665 (−0.0654)
−0.0338
−0.0547 (−0.0595)
−0.0338
−0.0879 (−0.083)
2
−0.1598 (−0.1547)
−0.0379 (−0.0345)
−0.0038
−0.028 (−0.0298)
−0.0061
−0.0669 (−0.0654)
3
−0.1146 (−0.1083)
−0.0092 (−0.0095)
0.0262
−0.0013 (−0.0015)
0.0216
−0.0459 (−0.0452)
4
−0.069 (−0.0654)
0.0195 (0.0189)
0.0562
0.0254 (0.0226)
0.0493
−0.025 (−0.0238)
5
−0.024 (−0.0214)
0.0481 (0.0476)
0.0862
0.0521 (0.0571)
0.077
−0.004 (−0.0038)
6
0.021
0.0768
0.1161
0.0788
0.1047
0.017
7
0.066
0.1054
0.1461
0.1054
0.1324
0.038
8
0.1115
0.1341
0.1761
0.1321
0.1601
0.0589
9
0.1568
0.1628
0.2061
0.1588
0.1878
0.0799
10
0.202
0.1914
0.2361
0.1855
0.2155
0.1009
11
0.2473
0.2201
0.2661
0.2122
0.2432
0.1219
12
0.2925
0.2487
0.296
0.2389
0.2709
0.1428

*The consistent results of Croce and Venini [18].Nondimensional center deflection of P-FGM in steady state without heat source present versus nondimensional load intensity (P-) (thermo-mechanical load): (a) neglecting neutral plane position, (b) considering neutral plane position.
(a)
(b)The variation of the axial stressσxx at the center of the FG plate and along the thickness direction, for different values of the material parameter n, is shown in Figure 7. The stress of power law FGM (P-FGM) plate can be represented as a cubic function of z for material parameter n=2. The maximum tensile stress along the thickness of the FG plate is located at the bottom edge (z=h/2) and increases as the ratio EU/EL increases. However, the maximum compressive stress is occurred at the top surface (z=-h/2) and is high for small EU/EL. For the ratio EU/EL=1 in which the FGM plate becomes a homogenous isotropic plate, the stress distribution is a linear function of z, and the maximum stress value is occurred at the top and bottom surfaces of the plate. The neutral plane—where zero axial stress is located—is shifted from the midplane to the top surface, (Figure 7(b)), by variation of grading parameter n, 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.Nondimensional axial stress distribution through the thickness of FG plate subjected to mechanical loads: (a) neglecting neutral plane position, (b) considering neutral plane position (*nondimensional axial stressσxx-=σxx(h2/a2p), *nondimensional coordinate across the thickness h-=z/h).
(a)
(b)Figure8 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].Nondimensional axial stress distribution, through the plate thickness, for steady state thermo-mechanical loads: (a) neglecting actual neutral plane position. (b) Considering the location of the neutral plane position (*nondimensional axial stressσxx-=σxx(h2/a2p), *nondimensional coordinate across the thickness h-=z/h).
(a)
(b)


### 4.2. 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×106 watt/m3. 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.Figure9 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 n, the FG plate becomes more sensitive for the amount of heat delivered from the heat source.Figure 9
Temperature distributions through the plate thickness for various values of grading parametern (with heat source present).Figure10 and Table 4 show the central deflection ω-=ω/h due to a sequence of mechanical loads for different values of n.Table 4
Nondimensional deflection of midpoint versus nondimensional load intensity (P-) (thermo-mechanical load with heat source present).
P
-
(ω-) Non dimensional central deflection for
Metal
n
=
1
n
=
2
Ceramic
1
−0.2051
−0.0387
−0.0471
−0.0879
2
−0.1599
−0.0087
−0.0194
−0.0669
3
−0.1146
0.0213
0.0083
−0.0459
4
−0.0694
0.0513
0.036
−0.025
5
−0.0241
0.0813
0.0637
−0.004
6
0.0211
0.1112
0.0914
0.017
7
0.0663
0.1412
0.1192
0.038
8
0.1116
0.1712
0.1469
0.0589
9
0.1568
0.2012
0.1746
0.0799
10
0.2021
0.2312
0.2023
0.1009
11
0.2473
0.2612
0.2300
0.1219
12
0.2925
0.2911
0.2577
0.1428Figure 10
Nondimensional center deflection of P-FG plate in steady state with heat source present versus nondimensional load intensity (P-) (thermo-mechanical load).Figure11, shows the maximum axial stress σxx distribution through the plate thickness for different material parameter n. The plate is expected to be a steady state heat conductive plate subjected to an amount of generated heat of strength (q˙=100×106 W/m3) 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×106 W/m3) and only thermal load. Figure 13 shows the nondimensional deflection of the FG steady state plate versus the heat source strength.Figure 11
Nondimensional axial stress distributions through the thickness of a steady state plate, with heat source present, for thermo-mechanical loads.Figure 12
Nondimensional axial stress distributions through the thickness of a steady state plate, with heat source present, for thermal loads.Figure 13
Non dimensional deflection of the central point of a steady state plate versus the heat source strength (P-=2).Figures6 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.Figures11 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.


## 4.1. 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 ceramic-metal. The lower plate surface is assumed to be aluminum, while the top surface is assumed to be zirconia. Material properties parametern=[0,2] is considered. Physical material properties are given in Table 1. An all edges simply supported square plate, of side a=0.2 m and thickness h=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).Table 1
Material properties.
Property
Aluminum
Zirconia
Young’s modulus
E
L = 70 Gpa
E
U = 151 Gpa
Poisson’s ratio
ν
L = 0.3
ν
U = 0.3
Thermal conductivity
k
L
=
204 W/mK
k
U
=
2.09 W/mK
Thermal expansion
α
L
=
23
×
10
-
6
/C
α
U
=
10
×
10
-
6
/CFigure 3
A simply supported FG plate subjected to a uniformly distributed mechanical load and thermal loading.To investigate the thermo-mechanical behavior of the plate, the temperature distribution through the FG plate thickness should be firstly determined. Figure4 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 n 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.Figure 4
Temperature distribution through the FG plate thickness for various values of grading parametern (without heat source).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 parametern. The following nondimensional parameters have been manipulated throughout the numerical simulations:central deflectionω-=ω/h;nondimensional load intensityP-=a4p/ELh4;thicknessh-=h/a,wherep denotes the intensity of the applied mechanical load and EL is the Young’s modulus of the aluminum bottom face.Figure5, 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 (n), 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 (n). Table 2 shows the numerical results of a set of numerical simulations with different values of (n). 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].Table 2
Nondimensional deflection of midpoint versus nondimensional load intensity (P-).
Nondimensional central deflection (ω-)
P
-
Metal
n
=
0.5
n
=
1
n
=
2
Ceramic
Midplane
Neutral plane
Midplane
Neutral plane
Midplane
Neutral plane
1
0.0452
0.0317
0.0328
0.0287
0.03
0.0267
0.0277
0.021
2
0.0905 (0.089)*
0.0633 (0.061)
0.0656
0.0573 (0.056)
0.06
0.0534 (0.052)
0.0554
0.042 (0.0413)
3
0.1357
0.095
0.0985
0.086
0.09
0.0801
0.0831
0.0629
4
0.181 (0.1804)
0.1266 (0.124)
0.1313
0.1146 (0.113)
0.119
0.1067 (0.103)
0.1108
0.0839 (0.083)
5
0.2262
0.1583
0.1641
0.1433
0.149
0.1334
0.1385
0.1049
6
0.2714 (0.269)
0.1899 (0.185)
0.1969
0.172 (0.1695)
0.179
0.1601 (0.154)
0.1662
0.1258 (0.126)
7
0.3167
0.2216
0.2298
0.2006
0.209
0.1868
0.1939
0.1468
8
0.3619 (0.361)
0.2532 (0.248)
0.2626
0.2293 (0.228)
0.239
0.2135 (0.208)
0.2216
0.1678 (0.166)
9
0.4072
0.2849
0.2954
0.2579
0.269
0.2402
0.2493
0.1888
10
0.4524 (0.452)
0.3166 (0.308)
0.3282
0.2866 (0.285)
0.299
0.2669 (0.261)
0.277
0.2097 (0.206)
11
0.4979
0.3482
0.361
0.3153
0.329
0.2935
0.3047
0.2307
12
0.5428 (0.54)
0.3799 (0.365)
0.3939
0.3439 (0.343)
0.359
0.3202 (0.313)
0.3324
0.2517 (0.252)

*The consistent results of Croce and Venini [18].Nondimensional center deflection of P-FGM plate versus nondimensional load intensity (mechanical load): (a) neglecting neutral plane position, (b) considering neutral plane position.
(a)
(b)Figure6 shows the central deflection ω-=ω/h due to a sequence of mechanical loads for different values of n. The nondimensional load P- 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].Table 3
Nondimensional deflection of midpoint versus nondimensional load intensity (P-) (thermo-mechanical load without heat source).
(ω-) Nondimensional central deflection for
P
-
Metal
n
=
1
n
=
2
Ceramic
Midplane
Neutral plane
Midplane
Neutral plane
1
−0.205 (−0.2)*
−0.0665 (−0.0654)
−0.0338
−0.0547 (−0.0595)
−0.0338
−0.0879 (−0.083)
2
−0.1598 (−0.1547)
−0.0379 (−0.0345)
−0.0038
−0.028 (−0.0298)
−0.0061
−0.0669 (−0.0654)
3
−0.1146 (−0.1083)
−0.0092 (−0.0095)
0.0262
−0.0013 (−0.0015)
0.0216
−0.0459 (−0.0452)
4
−0.069 (−0.0654)
0.0195 (0.0189)
0.0562
0.0254 (0.0226)
0.0493
−0.025 (−0.0238)
5
−0.024 (−0.0214)
0.0481 (0.0476)
0.0862
0.0521 (0.0571)
0.077
−0.004 (−0.0038)
6
0.021
0.0768
0.1161
0.0788
0.1047
0.017
7
0.066
0.1054
0.1461
0.1054
0.1324
0.038
8
0.1115
0.1341
0.1761
0.1321
0.1601
0.0589
9
0.1568
0.1628
0.2061
0.1588
0.1878
0.0799
10
0.202
0.1914
0.2361
0.1855
0.2155
0.1009
11
0.2473
0.2201
0.2661
0.2122
0.2432
0.1219
12
0.2925
0.2487
0.296
0.2389
0.2709
0.1428

*The consistent results of Croce and Venini [18].Nondimensional center deflection of P-FGM in steady state without heat source present versus nondimensional load intensity (P-) (thermo-mechanical load): (a) neglecting neutral plane position, (b) considering neutral plane position.
(a)
(b)The variation of the axial stressσxx at the center of the FG plate and along the thickness direction, for different values of the material parameter n, is shown in Figure 7. The stress of power law FGM (P-FGM) plate can be represented as a cubic function of z for material parameter n=2. The maximum tensile stress along the thickness of the FG plate is located at the bottom edge (z=h/2) and increases as the ratio EU/EL increases. However, the maximum compressive stress is occurred at the top surface (z=-h/2) and is high for small EU/EL. For the ratio EU/EL=1 in which the FGM plate becomes a homogenous isotropic plate, the stress distribution is a linear function of z, and the maximum stress value is occurred at the top and bottom surfaces of the plate. The neutral plane—where zero axial stress is located—is shifted from the midplane to the top surface, (Figure 7(b)), by variation of grading parameter n, 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.Nondimensional axial stress distribution through the thickness of FG plate subjected to mechanical loads: (a) neglecting neutral plane position, (b) considering neutral plane position (*nondimensional axial stressσxx-=σxx(h2/a2p), *nondimensional coordinate across the thickness h-=z/h).
(a)
(b)Figure8 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].Nondimensional axial stress distribution, through the plate thickness, for steady state thermo-mechanical loads: (a) neglecting actual neutral plane position. (b) Considering the location of the neutral plane position (*nondimensional axial stressσxx-=σxx(h2/a2p), *nondimensional coordinate across the thickness h-=z/h).
(a)
(b)


## 4.2. 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×106 watt/m3. 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.Figure9 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 n, the FG plate becomes more sensitive for the amount of heat delivered from the heat source.Figure 9
Temperature distributions through the plate thickness for various values of grading parametern (with heat source present).Figure10 and Table 4 show the central deflection ω-=ω/h due to a sequence of mechanical loads for different values of n.Table 4
Nondimensional deflection of midpoint versus nondimensional load intensity (P-) (thermo-mechanical load with heat source present).
P
-
(ω-) Non dimensional central deflection for
Metal
n
=
1
n
=
2
Ceramic
1
−0.2051
−0.0387
−0.0471
−0.0879
2
−0.1599
−0.0087
−0.0194
−0.0669
3
−0.1146
0.0213
0.0083
−0.0459
4
−0.0694
0.0513
0.036
−0.025
5
−0.0241
0.0813
0.0637
−0.004
6
0.0211
0.1112
0.0914
0.017
7
0.0663
0.1412
0.1192
0.038
8
0.1116
0.1712
0.1469
0.0589
9
0.1568
0.2012
0.1746
0.0799
10
0.2021
0.2312
0.2023
0.1009
11
0.2473
0.2612
0.2300
0.1219
12
0.2925
0.2911
0.2577
0.1428Figure 10
Nondimensional center deflection of P-FG plate in steady state with heat source present versus nondimensional load intensity (P-) (thermo-mechanical load).Figure11, shows the maximum axial stress σxx distribution through the plate thickness for different material parameter n. The plate is expected to be a steady state heat conductive plate subjected to an amount of generated heat of strength (q˙=100×106 W/m3) 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×106 W/m3) and only thermal load. Figure 13 shows the nondimensional deflection of the FG steady state plate versus the heat source strength.Figure 11
Nondimensional axial stress distributions through the thickness of a steady state plate, with heat source present, for thermo-mechanical loads.Figure 12
Nondimensional axial stress distributions through the thickness of a steady state plate, with heat source present, for thermal loads.Figure 13
Non dimensional deflection of the central point of a steady state plate versus the heat source strength (P-=2).Figures6 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.Figures11 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.


## 5. Conclusions

In this study, a finite element model based on the first-order 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.



---

*Source: 569781-2013-05-13.xml*