Dynamic Analysis of Functionally Graded Timoshenko Beams in Thermal Environment Using a Higher-Order Hierarchical Beam Element

A higher-order finite beam element for free and forced vibration analysis of functionally graded Timoshenko beams in thermal environment is formulated by using hierarchical functions to interpolate the kinematic variables. The shear strain is constrained to constant to improve the efficiency of the element. The effect of environmental temperature is taken into account in the element derivation by considering that the material properties are temperature-dependent and the temperature is nonlinear distribution in the beam thickness. The accuracy of the derived formulation is confirmed by comparing the results obtained in the present work with the published data. Numerical investigations show that the formulated element is efficient, and it is capable of giving accurate vibration characteristics by a small number of elements. A parametric study is carried out to highlight the effect of the material inhomogeneity, temperature rise, and loading parameter on the dynamic behaviour of the beams. The influence of the aspect ratio on the dynamic behaviour of the beam is also examined and highlighted.


Introduction
Functionally graded materials (FGMs), initiated by Japanese scientists in mid-1980s, have great potential to be used as structural materials in severe conditions.FGM structures in general and FGM beams in particular are increasingly used as structural components in aircraft and space vehicles.The static and vibration analyses of FGM beams by using different analytical and numerical methods are extensively reported in the literature [1][2][3][4][5][6][7][8][9][10][11][12][13]; contributions based on the finite element method are briefly discussed below.
Chakraborty et al. [14] derived a Timoshenko beam element for studying the static, free vibration and wave propagation problems in bimaterial beams fused with a FGM layer.The element, taking the effect of uniform temperature rise into account, is formulated by using the solution and static equilibrium equations of a Timoshenko beam segment to interpolate the kinematic variables.Bhangale and Ganesan [15] employed the finite element method to examine the influence of temperature on natural frequencies and loss factors of a functionally graded sandwich beam with constrained viscoelastic core.The third-order shear deformation beam theory was adopted by Kadoli et al. [16] to develop the stiffness matrix and load vector for stress analysis of FGM beams.Alshorbagy et al. [17] used the traditional Euler-Bernoulli beam element to compute the natural frequencies of beams with material properties to be graded in the thickness or longitudinal direction by a power-law distribution.The polynomials derived by Kosmatka [18] were employed by Shahba et al. [19] to derive the finite element formulation for studying the free vibration and buckling of tapered Timoshenko beams made of axially FGM.Kosmatka's polynomials were also employed by Nguyen et al. [20] in derivation of a finite element formulation for computing the dynamic response of FGM Timoshenko beams under a variable speed moving load.Eltaher et al. [21] considered the shift of the neutral axis in derivation of a finite beam element for the free vibration analysis of FGM macro-/nanobeams.The authors showed that the natural frequencies of the beams are slightly overestimated by ignoring the shift of the neutral axis position.In [22], Gan et al. presented a finite element procedure for the forced vibration analysis of axially FGM Timoshenko beams under multiple moving loads.In order to improve the accuracy and convergence of finite element solutions, Nguyen and Gan [23] and Nguyen and Tran [24] employed the exact solutions of static equilibrium equations of a beam element to formulate the corotational formulation for large displacement analysis of tapered FGM beams and FGM sandwich beams, respectively.Based on Euler-Bernoulli beam theory and differential quadrature method, Jin and Wang [25] derived a beam element for free vibration analysis of thin FGM beams.The element, respectively, used Lagrange and Hermite functions to interpolate the axial and transverse displacements, is accurate in evaluating the natural frequencies.Murin et al. [26] took into account the effects of shear force deformation and variation of longitudinal inertia and rotary inertia in derivation of a finite element formulation for modal analysis of FGM structures.The accuracy of the formulation derived in the work has been confirmed by comparing the obtained results with the ones using the commercial finite element package ANSYS.Based on the hierarchical beam theories, De Pietro et al. [27] formulated a number of finite beam elements for the thermal stress analysis of three-dimensional FGM beams.The elements derived via a unified formulation show high accuracy and efficiency.Adopting the higher-order shear deformation beam theory, Frikha et al. [28] formulated a two-node  0 beam element for analysis of FGM beams.Kahya and Turan [29] proposed a first-order finite element formulation for the free vibration and buckling analyses of FGM beams.The element with eleven degrees of freedom employs Lagrange polynomials to interpolate the displacements and rotation.Recently, Nguyen et al. [30] employed Kotsmatka's polynomials to formulate a beam element for studying the dynamic behaviour of bidirectional FGM Timoshenko beams excited by a moving load.
The beam elements in the above cited references employ Lagrange polynomials or polynomials derived from the static equilibrium equations to interpolate the displacements and rotation.Alternatively, a finite element can be derived by using hierarchical functions to interpolate the kinematic variables [31].In this regard, a Timoshenko beam element for dynamic analysis of FGM beams is in this paper by using the hierarchical functions to interpolate the displacements and rotation.The beam element based on the hierarchical functions, however, needs middle nodes, and this increases the number of degrees of freedom of the element.In order to improve the efficiency of the element, an additional procedure, in which the shear strain is constrained to constant, is applied to reduce the number of degrees of freedom.The effects of environmental temperature are considered in the element derivation by assuming that the material properties are temperature-dependent and the temperature is nonlinear distribution in the beam thickness.Numerical investigations are carried out to show the accuracy of the formulated element and to highlight the effects of the material inhomogeneity, temperature rise, and loading parameter on the vibration characteristics of the beams.

Mathematical Formulation
A FGM beam with length , rectangular cross section ( × ℎ) as depicted in Figure 1, is considered.The Cartesian coordinate system in the figure is chosen as that the -axis is on the midplane, the -axis is perpendicular to the midplane, and the -axis is perpendicular to the (, ) plane.The beam is assumed to be formed from ceramic and metal phases whose volume fraction varies in the  direction according to where   and   are, respectively, the volume fraction of ceramic and metal and  is the nonnegative grading index.In (1) and hereafter, the subscripts "" and "" stand for ceramic and metal, respectively.The FGM beams, as mentioned above, are widely used in high-temperature environment.In such environment, the material properties may be changed with temperature, and they should be considered to be temperature-dependent.In this regard, a typical material property () is varied with environmental temperature () as [32] where  (K) is the environmental temperature;  −1 ,  0 ,  1 ,  2 ,  3 are the coefficients of the temperature , and they are unique to the constituent materials.
Based on Voigt model, the effective material property P(, ) is evaluated as P (, ) = P    + P    . (3) From ( 1) and (3), Young's modulus (), thermal expansion (), and mass density () are expressed in terms of the constituent properties as In the above equation, the mass density is assumed to be independent of the temperature since it is hardly changed by temperature [9].The temperature is considered to vary in the thickness direction only, and it is assumed that the temperature is imposed to a reference temperature on the bottom surface,  =   (K) at  = −ℎ/2, and a value on the top surface,  =   at  = ℎ/2.In this case, the temperature distribution can be obtained by solving the steady-state heat transfer equation: where  is the thermal conductivity, assumed to be independent of the temperature.The solution of ( 5) is of the following form [11,33]: As seen from the above equation, ( 6) results in a uniform temperature rise when   =   ; otherwise it leads to a nonlinear distribution of the temperature in the beam thickness.The uniform temperature rise is simple from computational point of view, and only the nonlinear temperature rise is considered herein.It is assumed that the temperature on the bottom surface is kept at room temperature,   =  0 = 300 K, and a temperature rise (Δ) is defined as the difference between the temperature on the top surface with that on the bottom surface, Δ =   −  0 .Based on Timoshenko beam theory, the axial displacement ( 1 ) and transverse displacement ( 3 ) at any point of the beam are given by  1 (, , ) =  (, ) −  (, ) ,  3 (, , ) =  (, ) , (7) where  is the distance from the midplane to the considering point; (, ) and (, ) are, respectively, the axial and transverse displacements of the corresponding point on the midplane, and (, ) is the cross-sectional rotation.
The normal strain (  ) and the shear strain (  ) resulting from ( 7) are of the following forms: where a subscript comma is used to indicate the derivative of the variable with respect to the spatial coordinate , (⋅) , = (⋅)/.
Based on Hooke's law, the constitutive equation is of the following form: where   and   are, respectively, the normal and shear stresses, and (, ) is the effective shear modulus.
The strain energy of the beam () resulting from ( 8) and ( 9) is as follows: where  is the cross-sectional area;  is the shear correction factor, taken by 5/6 for the beams with rectangular cross section considered herein;   ,   ,   , and   denote the strain energy resulting from the axial stretching, stretchingbending coupling, bending, and shear deformation, respectively;  11 ,  12 ,  22 , and  33 are, respectively, the extensional stiffness, bending-coupling stiffness, bending stiffness, and shear stiffness, and they are defined as follows: It should be noted that, due to the temperature,  is a nonlinear function of  and numerical integration should be employed in computing the above stiffness components.Simpson's 3/8 rule is adopted herein.Suppose the beam is initially stress free at temperature  0 .The beam is initially stressed by the temperature rise.The initial stress due to temperature rise is the following: where Young's modulus (, ) and the thermal expansion (, ) are calculated from (4).The strain energy resulting from the temperature rise is then given by the following [7,12]: where   is the axial force caused by the elevated temperature, defined as with Δ, as mentioned above, being the temperature rise.
The kinetic energy for the beam (T) resulting from ( 7) is of the following form: where an over dot is used to denote the differentiation with respect to the time variable ;

Finite Element Formulation
The stiffness and mass matrices for an element with length of  are derived in this section.In Timoshenko beam theory, the axial displacement , transverse displacement , and rotation  are independent variables, and thus the interpolation functions for these variables can be chosen separately.Linear functions can be adopted for all these variables, but the element based on the linear functions encounters the shear-locking problem, and some technique such as the reduced integration must be applied to overcome this problem [34].Alternatively, the shear locking can be avoided by using appropriate interpolation functions for the independent variables.Standard polynomial-based shape functions can be adopted to approximate the displacement field of Timoshenko beam.The finite element formulation derived from the standard shape functions, however, has a serious drawback.Since the coefficients of the polynomials are determined from the element boundary conditions, related to nodal values of the variables, totally new shape functions have to be redetermined whenever the element refinement is made [35].The finite element method using the hierarchical functions, in which the higher-order shape functions contain the lowerorder ones, to interpolate the displacement field is able to overcome this problem.For one-dimensional element, the linear, quadratic, and cubic forms of the hierarchical shape functions are as follows [31]: with being the natural coordinate.For a Timoshenko beam element, a quadratic variation of the rotation should be chosen to represent linearly varying bending moments along the element.In addition, with the shear strain given by (8), the shape function for  should be one order higher than that of .In this regard, the displacements and rotation can be interpolated by using the hierarchical shape functions as follows: where  1 ,  2 , . . .,  4 are values of the variables at nodes. Figure 2 shows the shape functions   ( = 1, . . ., 4) in (17) and details of the nodal variables in (19).
A beam element for dynamic analysis can be formulated from nine nodal degrees of freedom in Figure 2.However, a more efficient element with less degrees of freedom can be derived by constraining the shear strain   to constant [36].Using (17) and ( 19) one can express the shear strain in (8) in the following form: For   = const, we need The above equation gives From ( 17) and ( 21), the kinematic variables in (19) can be written in the following forms: The shear strain in ( 20) is now of the following form: The beam element in the present work is formulated from the displacement field in (23) and the shear strain (24).The element vector of nodal displacements (d) has seven components as follows: where and hereafter the superscript '' is used to indicate the transpose of a vector or a matrix.Using ( 23) and (24), one can rewrite the strain energy (10) in the following form: where  is the total number of elements; k  , k  , k  , and k  are, respectively, the element stiffness matrices stemming from the axial stretching, stretching-bending coupling, bending, and shear deformation, and they have the following forms: It should be noted that zero coefficients corresponding to the nodal variables in (25) must be added to the above matrices to expand them to (7 × 7) matrices.
Similarly, the strain energy due to temperature rise (13) can be written as where k  is the stiffness matrix resulting from the temperature rise with the following form: The kinetic energy (15) can also be written as with m  , m  , m  , and m  are, respectively, the element mass matrices stemming from the axial and transverse translations, axial translation-rotation coupling, and cross-sectional rotation, and they have the following forms: Zero entries are also needed to insert to the matrices in (29) and (31) to expand them to the (7 × 7) matrices.

Numerical Investigation
Numerical investigation is carried out in this section to show the accuracy of the proposed beam element and the effects of the material inhomogeneity, temperature rise, and loading parameter on the vibration characteristics of FGM beams.
Only beams with simply supported ends are considered as follows.
4.1.Free Vibration.The equation for free vibration analysis of a FGM beam can be written as follows: where D and D are, respectively, the vectors of global nodal displacements and acceleration; =1 k  are the global mass matrix, global stiffness matrix due to the beam deformation, and the global stiffness matrix due to the temperature rise, respectively.The summation symbol is understood herewith as the assembly of the element matrices over the total number of elements.By assuming a harmonic form for the displacement vector D, (32) leads to an eigenvalue problem, and the standard solution method gives the natural frequencies and vibration modes [34].
The convergence and accuracy of the derived formulation in evaluating natural frequencies are firstly examined.In Table 1, the fundamental frequency parameter,  =  1  2 /ℎ √  11 / ∫ ℎ/2 −ℎ/2 () ( 1 is the fundamental frequency), obtained by different number of the elements, is given for  = 0.3, Δ = 0, and various values of the aspect ratio /ℎ.For the comparison purpose, the results based on an analytical method by Sina et al. [4] and a semianalytical method by Simsek [37] are also given in the table.The table shows an excellent agreement between the frequency parameter of the present paper with that of [4,37].The convergence of the present solution, as seen from Table 1, is fast, and it is capable of giving accurate frequency parameter with just sixteen elements.The rate of convergence of the present element is comparable with that of the beam element using Kosmatka shape functions as reported in [30], where the convergence of the element in evaluating the fundamental frequency of transverse FGM beams is achieved by using fourteen elements.In addition, the beam element derived in the present work is free of the shear locking, and it is able to evaluate accurately the frequency of the beam having a high aspect ratio, /ℎ = 100.In Table 2, the nondimensionless fundamental frequencies of a FGM beam are listed for various values of the temperature rise and the grading index , where the result obtained by an analytical method for a FGM Euler-Bernoulli beam theory in [12] is also given.A good agreement between the finite element solution of the present work with the cited reference is noted.The result in Table 2 was obtained for the beam formed from alumina (Al 2 O 3 ) and stainless steel (SUS304) with the temperature-dependent data listed in Table 3, and the dimensionless fundamental frequency is defined as follows: where  1 is the fundamental frequency;  0  and   are, respectively, Young's modulus and mass density of SUS304 at room temperature.
Figure 3 illustrates the relation between the dimensionless frequency  * with the grading index  of a FGM beam formed from Al 2 O 3 and SUS304 for various values of the temperature rise and the aspect ratio.The effects of temperature rise, grading index, and aspect ratio on the frequency are clearly seen from the figure.For a given value of the aspect ratio /ℎ, as seen from Figure 3(a), the frequency parameter decreases by increasing the index  and the temperature rise.The decrease of the frequency by increasing the index  is due to the fact that the beam associated with a higher index  contains more metal, and thus its stiffness, defined by (11), is lower.In a higher temperature environment, Young's modulus of the constituent materials decreases, and this also leads to a decrease in the beam stiffness.Figure 3(a) also shows that the effect of temperature rise on the frequency is more significant for the beam with a higher index .This can be explained by the fact that the metal is more sensitive with temperature than the ceramic.As mentioned above, the beam with a higher index  contains more metal, and thus for a given temperature rise its stiffness decreases more significantly than the one with a lower index .The effect of the aspect ratio on the dimensionless frequency, as seen from Figure 3(b), is similar to that of the temperature rise, and the effect is also more significant for the beam associated with a higher index .The effect of the temperature rise and the aspect ratio on the fundamental frequency can also be seen clearly from Table 4 where the dimensionless frequency is given for different values of the temperature rise and the aspect ratio.A careful examination of the table shows that, with /ℎ = 10, a decrease of 3.13% in the frequency when increasing the temperature rise from 50 K to 150 K is attained  from the beam with  = 0.2, but this value is 5.20% for the beam with  = 10.The corresponding values are 21.60% and 49.44% for the beam with /ℎ = 25, which are much higher than that of the beam with /ℎ = 10.

Forced Vibration.
The equations of motion for the forced vibration analysis in the context of finite element analysis can be written in the following form: where R ex is the consistent vector of external loads, which can be easily derived by using the shape functions.The equations of motion (34) can be solved by the direct integration Newmark method.The average acceleration method which ensures the unconditional stability [34] is adopted herein.The harmonic response of a simply supported FGM beam composed of stainless steel (SUS304) and alumina (Al 2 O 3 ) with the material data in Table 3 to a harmonic point load  =  0 sin(Ω) acting at the mid-span of the beam is firstly studied.In Figure 4, the time histories for the normalized mid-span deflections of the beam with /ℎ = 20 are shown for Ω = 10 rad/s and various values of the grading index  and the temperature rise Δ.The deflection in the figure is normalized by the static mid-span deflection of the pure steel beam, that is, with  being the moment of inertia and  0  being Young's modulus of SUS304 at room temperature.
The small red circles in Figure 4(a) are the analytical solution of [38] for a homogeneous beam (pure alumina beam), which is of the following form: where  = /2, and  = √/.A good agreement between the finite element solution of the present work with the analytical solution of [38] is noted.The effects of the material inhomogeneity and the temperature rise are clearly seen from Figure 4.At the given value of the excitation frequency, the maximum mid-span deflection in Figure 4(a) is increased with an increase of the grading index .The influence of the temperature rise, as seen from Figure 4(b), is similar to that of the grading index , and the maximum mid-span deflection gradually increases by increasing the temperature rise.The increase of the maximum mid-span deflection by increasing the index  and the temperature rise resulted from the decrease of the beam stiffness, as already explained above.The effects of material inhomogeneity and temperature rise on the harmonic response of the beam can also be seen from Table 5, where the maximum normalized mid-span deflections are  given for Ω = 10 rad/s and various values of the index , temperature rise Δ, and aspect ratio /ℎ.The maximum mid-span deflection in the table steadily increases by the increase of the index , regardless of the temperature rise and the aspect ratio.The temperature rise also increases the deflection, and, as in the case of the fundamental frequency, the effect of temperature rise is more significant for the beam with a higher aspect ratio.It should be noted that 500 time steps have been used for the Newmark method, and the numerical results in Figure 4 and Table 5 are converged by using sixteen elements.
In the next numerical investigation, the forced vibration of a simply supported FGM beam with /ℎ = 20 excited by a moving load is studied.The beam is also assumed to be formed from SUS304 and Al 2 O 3 , and the moving load speed V is considered to be constant.The problem has been previously investigated by several authors, including Simsek and Kocatürk [39] and Khalili et al. [40], for a FGM beam in room temperature.
In Table 6, the maximum normalized mid-span deflections of the FGM beam in room temperature subjected to the moving load are compared to the result obtained in several papers.Regardless of the grading index , the table shows a good agreement between the result of the present work with that obtained by using a semianalytical method in [39] and the different transformation method in [40].In addition, the deflections evaluated by the present beam element are in excellent agreement with the finite element solution in [30], where the beam element based on Kosmatka's shape functions was used.Note that the result in Table 6 has been obtained by using the geometric and material data stated in [39].Figure 5 shows the effect of temperature rise on the time histories for the normalized mid-span deflection of the FGM beam for  = 0.5 and two values of the moving speed, V = 30 m/s and V = 60 m/s.The mid-span deflection of the beam significantly increases by the increase of the index  and the temperature rise.The beam tends to execute less vibration cycles by increasing the index  and the temperature rise, regardless of the moving load speed.The effect of the moving load speed, temperature rise, and material inhomogeneity on the dynamic response of the beam can also be seen from Figure 6, where the relation between the maximum normalized mid-span deflection with the moving load speed of the beam is depicted for different values of the grading index  and the temperature rise Δ.The peak of the curves in Figure 6 increases and it is attained at a lower moving load speed when increasing the index  and the temperature rise.The numerical results depicted in Figures 5 and 6 have been computed for the beam with an aspect ratio /ℎ = 20, and also 500 time steps have been used for the Newmark method.

Conclusions
A higher-order beam element for dynamic analysis of FGM Timoshenko beams in thermal environment has been formulated by using the hierarchical functions to interpolate the displacement and rotation variables.The shear strain was constrained to constant for reducing the number of degrees of freedom of the element.The material properties are considered to be temperature-dependent and the temperature is assumed to be nonlinear distribution in the beam thickness.Explicit expressions for the element stiffness and mass matrices have been given in detail.Numerical examples were presented to show the accuracy and efficiency of the derived element.The numerical results showed that the convergence of the formulated element is fast, and it is capable of giving accurate vibration characteristics by using a small number of elements.The parametric study reveals that the material inhomogeneity, temperature rise, and loading parameters have significant influence on the dynamic behaviour of the FGM beams.The effect of the temperature rise on the vibration characteristics of the FGM beams is more significant for the beam with a higher aspect ratio.Though the numerical investigations have been carried out for simply supported beams under a harmonic load and a moving point load, the beam element derived in the present work can be applied for the beam with other boundary conditions and other dynamic loads as well.

Figure 1 :
Figure 1: A FGM beam in Cartesian coordinate system.

Figure 5 :
Figure 5: Effect of temperature rise on time histories for maximum normalized mid-span deflection of FGM beam with  = 0.5.

Figure 6 :
Figure 6: Relation between maximum normalized deflection with moving load speed: (a) Δ = 80 K and  is variable; (b)  = 1 and Δ is variable.

Table 2 :
Comparison of dimensionless fundamental frequency ( * ) of FGM beam in thermal environment.

Table 4 :
Dimensionless frequency of FGM beam with different aspect ratios and temperature rises.

Table 5 :
Maximum normalized mid-span deflection of FGM beam under harmonic load with different aspect ratios and temperature rises (Ω = 10 Hz).

Table 6 :
Comparison of maximum normalized mid-span deflection for FGM beam under a moving load (Δ = 0).