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

Authors: Dinh Kien Nguyen; Van Tuyen Bui
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

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.


---

## Body

## 1. 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–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 C0 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.


## 2. Mathematical Formulation

A FGM beam with lengthL, rectangular cross section (b×h) as depicted in Figure 1, is considered. The Cartesian coordinate system in the figure is chosen as that the x-axis is on the midplane, the z-axis is perpendicular to the midplane, and the y-axis is perpendicular to the (x,z) plane. The beam is assumed to be formed from ceramic and metal phases whose volume fraction varies in the z direction according to(1)Vc=zh+12n,Vm=1-Vc,where Vc and Vm are, respectively, the volume fraction of ceramic and metal and n is the nonnegative grading index. In (1) and hereafter, the subscripts “c” and “m” stand for ceramic and metal, respectively.Figure 1
A FGM beam in Cartesian coordinate system.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(P) is varied with environmental temperature (T) as [32](2)P=P0P-1T-1+1+P1T+P2T2+P3T3,where T (K) is the environmental temperature; P-1, P0, P1, P2, P3 are the coefficients of the temperature T, and they are unique to the constituent materials.Based on Voigt model, the effective material propertyP(z,T) is evaluated as(3)Pz,T=PcVc+PmVm.From (1) and (3), Young’s modulus (E), thermal expansion (α), and mass density (ρ) are expressed in terms of the constituent properties as(4)Ez,T=EcT-EmTzh+12n+EmT,αz,T=αcT-αmTzh+12n+αmT,ρz=ρc-ρmzh+12n+ρm.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,T=Tm(K) at z=-h/2, and a value on the top surface, T=Tcatz=h/2. In this case, the temperature distribution can be obtained by solving the steady-state heat transfer equation:(5)-ddzκzdTdz=0,where κ is the thermal conductivity, assumed to be independent of the temperature. The solution of (5) is of the following form [11, 33]:(6)Tz=Tm+Tc-Tm∫-h/2h/2dz/κz∫-h/2zdzκz.As seen from the above equation, (6) results in a uniform temperature rise when Tc=Tm; 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, Tm=T0=300 K, and a temperature rise (ΔT) is defined as the difference between the temperature on the top surface with that on the bottom surface, ΔT=Tc-T0.Based on Timoshenko beam theory, the axial displacement(u1) and transverse displacement (u3) at any point of the beam are given by(7)u1x,z,t=ux,t-zθx,t,u3x,z,t=wx,t,where z is the distance from the midplane to the considering point; u(x,t) and w(x,t) are, respectively, the axial and transverse displacements of the corresponding point on the midplane, and θ(x,t) is the cross-sectional rotation.The normal strain(εxx) and the shear strain (γxz) resulting from (7) are of the following forms:(8)εxx=u,x-zθ,x,γxz=w,x-θ,where a subscript comma is used to indicate the derivative of the variable with respect to the spatial coordinate x, (·),x=∂(·)/∂x.Based on Hooke’s law, the constitutive equation is of the following form:(9)σxxτxz=Ez,T00Gz,Tεxxγxz,where σxx and τxz are, respectively, the normal and shear stresses, and G(z,T) is the effective shear modulus.The strain energy of the beam(U) resulting from (8) and (9) is as follows:(10)U=12∫0L∫Aσxxεxx+τxzγxzdAdx=Ua+Uc+Ub+Us=12∫0LA11u,x2-2A12u,xθ,x+A22θ,x2+ψA33w,x-θ2dx,where A is the cross-sectional area; ψ is the shear correction factor, taken by 5/6 for the beams with rectangular cross section considered herein; Ua, Uc, Ub, and Us denote the strain energy resulting from the axial stretching, stretching-bending coupling, bending, and shear deformation, respectively; A11, A12, A22, and A33 are, respectively, the extensional stiffness, bending-coupling stiffness, bending stiffness, and shear stiffness, and they are defined as follows:(11)A11,A12,A22=∫AEz,T1,z,z2dA,A33=∫AGz,TdA.It should be noted that, due to the temperature, T is a nonlinear function of z 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 temperatureT0. The beam is initially stressed by the temperature rise. The initial stress due to temperature rise is the following:(12)σxxT=-Ez,Tαz,TΔTz,where Young’s modulus E(z,T) and the thermal expansion α(z,T) are calculated from (4).The strain energy resulting from the temperature rise is then given by the following [7, 12]:(13)UT=12∫0LNTw,x2dx,where NT is the axial force caused by the elevated temperature, defined as(14)NT=∫AσxxTdA=-∫AEz,Tαz,TΔTzdAwith ΔT, as mentioned above, being the temperature rise.The kinetic energy for the beam(T) resulting from (7) is of the following form:(15)T=12∫0LI11u˙2+I11w˙2-2I12u˙θ˙+I22θ˙2dx=Tu+Tw+Tc+Tθ,where an over dot is used to denote the differentiation with respect to the time variable t; Tu, Tw, Tc, and Tθ are, respectively, the kinetic energy stemming from the axial and transverse translations, axial translation-rotation coupling, and rotation of the cross section; I11, I12, and I22 are the mass moments which are defined as follows:(16)I11,I12,I22=∫Aρz1,z,z2dA.


## 3. Finite Element Formulation

The stiffness and mass matrices for an element with length ofl are derived in this section. In Timoshenko beam theory, the axial displacement u, transverse displacement w, 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 lower-order 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]:(17)N1=121-ξ,N2=121+ξ,N3=1-ξ2,N4=ξ1-ξ2,with(18)ξ=2xl-1being 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 w 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:(19)u=N1u1+N2u2,θ=N1θ1+N2θ2+N3θ3,w=N1w1+N2w2+N3w3+N4w4,where u1,u2,…,w4 are values of the variables at nodes. Figure 2 shows the shape functions Ni(i=1,…,4) in (17) and details of the nodal variables in (19).Figure 2
(a) Hierarchical shape functions; (b) detail on nodal displacements and rotations.
(a)
(b)A beam element for dynamic analysis can be formulated from nine nodal degrees of freedom in Figure2. However, a more efficient element with less degrees of freedom can be derived by constraining the shear strain γxz to constant [36]. Using (17) and (19) one can express the shear strain in (8) in the following form:(20)γxz=-6lw4+θ3ξ2-4lw3-12θ1+12θ2ξ+1lw2-w1+2w4-12θ1+θ2+2θ3.For γxz=const, we need(21)-6lw4+θ3=0,4lw3-12θ1+12θ2=0.The above equation gives(22)w3=l8θ1-θ2,w4=l6θ3.From (17) and (21), the kinematic variables in (19) can be written in the following forms:(23)u=121-ξu1+121+ξu2,θ=121-ξθ1+121+ξθ2+1-ξ2θ3,w=121-ξw1+121+ξw2+l81-ξ2θ1-θ2+l6ξ1-ξ2θ3.The shear strain in (20) is now of the following form:(24)γxz=1lw2-w1-12θ1+θ2-23θ3.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:(25)d=u1w1θ1θ3u2w2θ2T,where and hereafter the superscript ‘T’ 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:(26)U=12∑i=1nEdTka+kc+kb+ksd,where nE is the total number of elements; ka, kc, kb, and ks are, respectively, the element stiffness matrices stemming from the axial stretching, stretching-bending coupling, bending, and shear deformation, and they have the following forms:(27)ka=A11l1-1-11,kc=-A12l0100-1100-10000000-1001-10010,kb=A22l10-101630-101,ks=ψA33l1l22l3-1l2l2l24l23-l2l242l3l2349l2-2l3l23-1-l2-2l31-l2l2l24l23-l2l24.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(28)UT=12∑i=1nEdTkTd,where kT is the stiffness matrix resulting from the temperature rise with the following form:(29)kT=NTl100-100112l200-112l200445l200-100100-112l200112l2.The kinetic energy (15) can also be written as(30)T=12∑i=1nEd˙Tmu+mw+mc+mθd˙,with mu, mw, mc, 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:(31)mu=lI113112121,mw=lI1131l8-l3012-l8l8l2400l8-l240-l3002l2315l30012l8l301-l8-l8-l2400-l8l240,mc=-lI12301101210012010010012101120010,mθ=lI223111218511211.Zero entries are also needed to insert to the matrices in (29) and (31) to expand them to the (7×7) matrices.


## 4. 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:(32)MD¨+KB+KTD=0,where D and D¨ are, respectively, the vectors of global nodal displacements and acceleration; M=∑i=1nE(mu+mw+mc+mθ), KB=∑i=1nE(ka+kc+kb+ks), and KT=∑i=1nEkT 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 Table1, the fundamental frequency parameter, ω¯=ω1L2/hI11/∫-h/2h/2E(z)dz (ω1 is the fundamental frequency), obtained by different number of the elements, is given for n=0.3, ΔT=0, and various values of the aspect ratio L/h. 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, L/h=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 n, 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 (Al2O3) and stainless steel (SUS304) with the temperature-dependent data listed in Table 3, and the dimensionless fundamental frequency is defined as follows:(33)ω∗=ω1L2hρmEm0,where ω1 is the fundamental frequency; Em0 and ρm are, respectively, Young’s modulus and mass density of SUS304 at room temperature.Table 1
Convergence and accuracy of the element in evaluating frequency parameter(ω-) of FGM beam (n=0.3, ΔT=0).
L
/
h
n
E
Ref. [4]
Ref. [37]
8
12
14
16
18
20
10
2.7021
2.7015
2.7014
2.7014
2.7013
2.7013
2.695
2.701
30
2.7386
2.7383
2.7382
2.7381
2.7381
2.7381
2.737
2.738
100
2.7429
2.7426
2.7425
2.7424
2.7424
2.7424
2.742
2.742Table 2
Comparison of dimensionless fundamental frequency(ω∗) of FGM beam in thermal environment.
Δ
T (K)
Source
n
=
0.1
n
=
0.2
n
=
0.5
n
=
1
20
Present
4.6950
4.4458
3.9781
3.6093
Ref. [12]
4.7018
4.4333
3.9353
3.5473
40
Present
4.5964
4.3471
3.8785
3.5079
Ref. [12]
4.6020
4.3278
3.8140
3.4112
80
Present
4.3931
4.1428
3.6704
3.2945
Ref. [12]
4.3956
4.1087
3.5590
3.1214Table 3
Temperature-dependent coefficients ofAl2O3 and SUS304 [12].
Material
Property
P
0
P
-
1
P
1
P
2
P
3
A
l
2
O
3
E (Pa)
349.55
e
+
9
0
-
3.853
e
-
4
4.027
e
-
7
-
1.673
e
-
10
ρ (kg/m3)
3800
0
0
0
0
α (K-1)
6.8269
e
-
6
0
1.838
e
-
4
0
0
κ (W/mK)
-
14.087
-
1123.6
-
6.227
e
-
3
0
0
SUS304
E (Pa)
201.04
e
+
9
0
3.079
e
-
4
-
6.534
e
-
7
0
ρ (kg/m3)
8166
0
0
0
0
α (K-1)
12.330
e
-
6
0
8.086
e
-
4
0
0
κ (W/mK)
15.379
0
-
1.264
e
-
3
2.092
e
-
6
-
7.223
e
-
10Figure3 illustrates the relation between the dimensionless frequency ω∗ with the grading index n of a FGM beam formed from Al2O3 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 L/h, as seen from Figure 3(a), the frequency parameter decreases by increasing the index n and the temperature rise. The decrease of the frequency by increasing the index n is due to the fact that the beam associated with a higher index n 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 n. 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 n contains more metal, and thus for a given temperature rise its stiffness decreases more significantly than the one with a lower index n. 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 n. 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 L/h=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 n=0.2, but this value is 5.20% for the beam with n=10. The corresponding values are 21.60% and 49.44% for the beam with l/h=25, which are much higher than that of the beam with L/h=10.Table 4
Dimensionless frequency of FGM beam with different aspect ratios and temperature rises.
L
/
h
Δ
T (K)
n
=
0.2
n
=
0.5
n
=
1
n
=
3
n
=
10
10
50
4.4184
3.9576
3.5941
3.1895
2.9432
100
4.3496
3.8889
3.5250
3.1185
2.8688
150
4.2803
3.8193
3.4543
3.0447
2.7902
15
50
4.3864
3.9204
3.5524
3.1417
2.8891
100
4.2408
3.7730
3.4021
2.9848
2.7230
150
4.0904
3.6193
3.2435
2.8164
2.5413
25
50
4.1663
3.6932
3.3171
2.8928
2.6232
100
3.7439
3.2564
2.8617
2.4014
2.0858
150
3.2663
2.7477
2.3123
1.7661
1.3263Figure 3
Relation between dimensionless frequencyω∗ and grading index n: (a) L/h=20 and ΔT is variable; (b) ΔT=80 K and L/h is variable.
(a)
(b)


### 4.2. 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:(34)MD¨+KB+KTD=Rex,where Rex 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 (Al2O3) with the material data in Table 3 to a harmonic point load F=F0sin(Ωt) 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 L/h=20 are shown for Ω=10 rad/s and various values of the grading index n and the temperature rise ΔT. The deflection in the figure is normalized by the static mid-span deflection of the pure steel beam, that is,(35)w0=F0L348Em0Iwith I being the moment of inertia and Em0 being Young’s modulus of SUS304 at room temperature.Figure 4
Dynamic response of FGM beam withL/h=20 to a harmonic load: (a) ΔT=0 and n is variable; (b) n=0.3 and ΔT is variable (Ω=10 rad/s).
(a)
(b)The small red circles in Figure4(a) are the analytical solution of [38] for a homogeneous beam (pure alumina beam), which is of the following form:(36)w=2F0L3ρA∑i=1∞siniπc/Lsiniπx/Li4a4π2-Ω2L4sin⁡Ωt-2ΩF0L5ρAπ2a∑i=1∞siniπc/Lsiniπx/Li2i4a4π2-Ω2L4sin⁡i2π2atL2,where c=L/2, and a=EI/ρA. 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 Figure4. 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 n. The influence of the temperature rise, as seen from Figure 4(b), is similar to that of the grading index n, 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 n 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 n, temperature rise ΔT, and aspect ratio L/h. The maximum mid-span deflection in the table steadily increases by the increase of the index n, 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.Table 5
Maximum normalized mid-span deflection of FGM beam under harmonic load with different aspect ratios and temperature rises (Ω=10 Hz).
L
/
h
Δ
T (K)
n
=
0.2
n
=
0.5
n
=
1
n
=
3
n
=
10
10
50
0.7866
0.8515
0.9095
0.9869
1.0628
100
0.8135
0.8831
0.9476
1.0327
1.1179
150
0.8401
0.9163
0.9893
1.0858
1.1791
15
50
0.8384
0.9315
1.0058
1.1130
1.2084
100
0.9118
1.0007
1.1002
1.2397
1.3815
150
0.9795
1.0831
1.2244
1.3992
1.6050
20
50
0.9759
1.0843
1.2074
1.3765
1.5265
100
1.0653
1.2672
1.3892
1.7410
2.0633
150
1.3084
1.4657
1.8049
2.3815
3.0587In the next numerical investigation, the forced vibration of a simply supported FGM beam withL/h=20 excited by a moving load is studied. The beam is also assumed to be formed from SUS304 and Al2O3, 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 Table6, 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 n, 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].Table 6
Comparison of maximum normalized mid-span deflection for FGM beam under a moving load (ΔT=0).
Source
A
l
2
O
3
n
=
0.2
n
=
0.5
n
=
1
n
=
2
SUS304
(
252
m
/
s
)
∗
(222 m/s)
(198 m/s)
(179 m/s)
(164 m/s)
(132 m/s)
Present
0.9382
1.0306
1.1509
1.2569
1.3450
1.7422
Ref. [30]
0.9380
1.0402
1.1505
1.2566
1.3446
1.7420
Ref. [39]
0.9328
1.0344
1.1444
1.2503
1.3376
1.7324
Ref. [40]
0.9317
1.0333
1.1429
1.2486
1.3359
1.7301
Note. ∗Moving load speed.Figure5 shows the effect of temperature rise on the time histories for the normalized mid-span deflection of the FGM beam for n=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 n and the temperature rise. The beam tends to execute less vibration cycles by increasing the index n 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 n and the temperature rise ΔT. The peak of the curves in Figure 6 increases and it is attained at a lower moving load speed when increasing the index n and the temperature rise. The numerical results depicted in Figures 5 and 6 have been computed for the beam with an aspect ratio L/h=20, and also 500 time steps have been used for the Newmark method.Figure 5
Effect of temperature rise on time histories for maximum normalized mid-span deflection of FGM beam withn=0.5.
(a)
v
=
30 m/s
(b)
v
=
60 m/sFigure 6
Relation between maximum normalized deflection with moving load speed: (a)ΔT=80 K and n is variable; (b) n=1 and ΔT is variable.
(a)
(b)


## 4.1. Free Vibration

The equation for free vibration analysis of a FGM beam can be written as follows:(32)MD¨+KB+KTD=0,where D and D¨ are, respectively, the vectors of global nodal displacements and acceleration; M=∑i=1nE(mu+mw+mc+mθ), KB=∑i=1nE(ka+kc+kb+ks), and KT=∑i=1nEkT 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 Table1, the fundamental frequency parameter, ω¯=ω1L2/hI11/∫-h/2h/2E(z)dz (ω1 is the fundamental frequency), obtained by different number of the elements, is given for n=0.3, ΔT=0, and various values of the aspect ratio L/h. 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, L/h=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 n, 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 (Al2O3) and stainless steel (SUS304) with the temperature-dependent data listed in Table 3, and the dimensionless fundamental frequency is defined as follows:(33)ω∗=ω1L2hρmEm0,where ω1 is the fundamental frequency; Em0 and ρm are, respectively, Young’s modulus and mass density of SUS304 at room temperature.Table 1
Convergence and accuracy of the element in evaluating frequency parameter(ω-) of FGM beam (n=0.3, ΔT=0).
L
/
h
n
E
Ref. [4]
Ref. [37]
8
12
14
16
18
20
10
2.7021
2.7015
2.7014
2.7014
2.7013
2.7013
2.695
2.701
30
2.7386
2.7383
2.7382
2.7381
2.7381
2.7381
2.737
2.738
100
2.7429
2.7426
2.7425
2.7424
2.7424
2.7424
2.742
2.742Table 2
Comparison of dimensionless fundamental frequency(ω∗) of FGM beam in thermal environment.
Δ
T (K)
Source
n
=
0.1
n
=
0.2
n
=
0.5
n
=
1
20
Present
4.6950
4.4458
3.9781
3.6093
Ref. [12]
4.7018
4.4333
3.9353
3.5473
40
Present
4.5964
4.3471
3.8785
3.5079
Ref. [12]
4.6020
4.3278
3.8140
3.4112
80
Present
4.3931
4.1428
3.6704
3.2945
Ref. [12]
4.3956
4.1087
3.5590
3.1214Table 3
Temperature-dependent coefficients ofAl2O3 and SUS304 [12].
Material
Property
P
0
P
-
1
P
1
P
2
P
3
A
l
2
O
3
E (Pa)
349.55
e
+
9
0
-
3.853
e
-
4
4.027
e
-
7
-
1.673
e
-
10
ρ (kg/m3)
3800
0
0
0
0
α (K-1)
6.8269
e
-
6
0
1.838
e
-
4
0
0
κ (W/mK)
-
14.087
-
1123.6
-
6.227
e
-
3
0
0
SUS304
E (Pa)
201.04
e
+
9
0
3.079
e
-
4
-
6.534
e
-
7
0
ρ (kg/m3)
8166
0
0
0
0
α (K-1)
12.330
e
-
6
0
8.086
e
-
4
0
0
κ (W/mK)
15.379
0
-
1.264
e
-
3
2.092
e
-
6
-
7.223
e
-
10Figure3 illustrates the relation between the dimensionless frequency ω∗ with the grading index n of a FGM beam formed from Al2O3 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 L/h, as seen from Figure 3(a), the frequency parameter decreases by increasing the index n and the temperature rise. The decrease of the frequency by increasing the index n is due to the fact that the beam associated with a higher index n 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 n. 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 n contains more metal, and thus for a given temperature rise its stiffness decreases more significantly than the one with a lower index n. 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 n. 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 L/h=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 n=0.2, but this value is 5.20% for the beam with n=10. The corresponding values are 21.60% and 49.44% for the beam with l/h=25, which are much higher than that of the beam with L/h=10.Table 4
Dimensionless frequency of FGM beam with different aspect ratios and temperature rises.
L
/
h
Δ
T (K)
n
=
0.2
n
=
0.5
n
=
1
n
=
3
n
=
10
10
50
4.4184
3.9576
3.5941
3.1895
2.9432
100
4.3496
3.8889
3.5250
3.1185
2.8688
150
4.2803
3.8193
3.4543
3.0447
2.7902
15
50
4.3864
3.9204
3.5524
3.1417
2.8891
100
4.2408
3.7730
3.4021
2.9848
2.7230
150
4.0904
3.6193
3.2435
2.8164
2.5413
25
50
4.1663
3.6932
3.3171
2.8928
2.6232
100
3.7439
3.2564
2.8617
2.4014
2.0858
150
3.2663
2.7477
2.3123
1.7661
1.3263Figure 3
Relation between dimensionless frequencyω∗ and grading index n: (a) L/h=20 and ΔT is variable; (b) ΔT=80 K and L/h is variable.
(a)
(b)


## 4.2. 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:(34)MD¨+KB+KTD=Rex,where Rex 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 (Al2O3) with the material data in Table 3 to a harmonic point load F=F0sin(Ωt) 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 L/h=20 are shown for Ω=10 rad/s and various values of the grading index n and the temperature rise ΔT. The deflection in the figure is normalized by the static mid-span deflection of the pure steel beam, that is,(35)w0=F0L348Em0Iwith I being the moment of inertia and Em0 being Young’s modulus of SUS304 at room temperature.Figure 4
Dynamic response of FGM beam withL/h=20 to a harmonic load: (a) ΔT=0 and n is variable; (b) n=0.3 and ΔT is variable (Ω=10 rad/s).
(a)
(b)The small red circles in Figure4(a) are the analytical solution of [38] for a homogeneous beam (pure alumina beam), which is of the following form:(36)w=2F0L3ρA∑i=1∞siniπc/Lsiniπx/Li4a4π2-Ω2L4sin⁡Ωt-2ΩF0L5ρAπ2a∑i=1∞siniπc/Lsiniπx/Li2i4a4π2-Ω2L4sin⁡i2π2atL2,where c=L/2, and a=EI/ρA. 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 Figure4. 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 n. The influence of the temperature rise, as seen from Figure 4(b), is similar to that of the grading index n, 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 n 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 n, temperature rise ΔT, and aspect ratio L/h. The maximum mid-span deflection in the table steadily increases by the increase of the index n, 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.Table 5
Maximum normalized mid-span deflection of FGM beam under harmonic load with different aspect ratios and temperature rises (Ω=10 Hz).
L
/
h
Δ
T (K)
n
=
0.2
n
=
0.5
n
=
1
n
=
3
n
=
10
10
50
0.7866
0.8515
0.9095
0.9869
1.0628
100
0.8135
0.8831
0.9476
1.0327
1.1179
150
0.8401
0.9163
0.9893
1.0858
1.1791
15
50
0.8384
0.9315
1.0058
1.1130
1.2084
100
0.9118
1.0007
1.1002
1.2397
1.3815
150
0.9795
1.0831
1.2244
1.3992
1.6050
20
50
0.9759
1.0843
1.2074
1.3765
1.5265
100
1.0653
1.2672
1.3892
1.7410
2.0633
150
1.3084
1.4657
1.8049
2.3815
3.0587In the next numerical investigation, the forced vibration of a simply supported FGM beam withL/h=20 excited by a moving load is studied. The beam is also assumed to be formed from SUS304 and Al2O3, 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 Table6, 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 n, 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].Table 6
Comparison of maximum normalized mid-span deflection for FGM beam under a moving load (ΔT=0).
Source
A
l
2
O
3
n
=
0.2
n
=
0.5
n
=
1
n
=
2
SUS304
(
252
m
/
s
)
∗
(222 m/s)
(198 m/s)
(179 m/s)
(164 m/s)
(132 m/s)
Present
0.9382
1.0306
1.1509
1.2569
1.3450
1.7422
Ref. [30]
0.9380
1.0402
1.1505
1.2566
1.3446
1.7420
Ref. [39]
0.9328
1.0344
1.1444
1.2503
1.3376
1.7324
Ref. [40]
0.9317
1.0333
1.1429
1.2486
1.3359
1.7301
Note. ∗Moving load speed.Figure5 shows the effect of temperature rise on the time histories for the normalized mid-span deflection of the FGM beam for n=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 n and the temperature rise. The beam tends to execute less vibration cycles by increasing the index n 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 n and the temperature rise ΔT. The peak of the curves in Figure 6 increases and it is attained at a lower moving load speed when increasing the index n and the temperature rise. The numerical results depicted in Figures 5 and 6 have been computed for the beam with an aspect ratio L/h=20, and also 500 time steps have been used for the Newmark method.Figure 5
Effect of temperature rise on time histories for maximum normalized mid-span deflection of FGM beam withn=0.5.
(a)
v
=
30 m/s
(b)
v
=
60 m/sFigure 6
Relation between maximum normalized deflection with moving load speed: (a)ΔT=80 K and n is variable; (b) n=1 and ΔT is variable.
(a)
(b)


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



---

*Source: 7025750-2017-09-20.xml*