Free vibration analysis of functionally graded beams is carried out in various classical boundary conditions. Two separate finite element formulations, one based on Euler-Bernoulli beam theory and the other based on Timoshenko beam theory are developed. Principle of virtual work is used to obtain the finite element system of equations.Numerical results are provided to demonstrate the effect of transverse shear on the natural frequencies and mode shapes for different length-to-thickness ratios and volume fraction exponents of functionally graded material (FGM) beams for the boundary conditions considered. It is observed that transverse shear significantly affects the fundamental frequency and mode shape for lower length to thickness ratios of FGM beams. Further, the effect is observed to be more prominent at higher modes for all the volume fraction exponents of FGM beam.

Keywords:    Functionally graded material ,  finite element method,  free vibration,  Euler-Bernoulli beam theory,  Timoshenko beam theory 

Two or more constituent materials are combined in a functionally graded material (FGM). The FGM can be used in various engineering applications such as aerospace, nuclear, automotive and steel industry. Structures subjected to aero-thermal loads, especially base shroud, inter-stages and re-entry vehicle of missiles, supporting structures for antennas, structures subjected to severe thermal gradients in nuclear and steel industry can be the more specific examples of the use of FGM’s. The above mentioned structural systems contain many structural members. Some of these structural members can be idealised as beams, whereas the others can be idealised as plates and shells. It is always easy for an analyst to find out the dynamic behaviour in terms of natural frequencies and mode shapes. Once this information is available for all the structural members of a system, the dynamic behaviour of the entire assembly can be obtained by component mode synthesis. With this idea, an attempt is made in this study to evaluate the frequencies and mode shapes of one of the structural members, i.e. FGM beam.

Numerous literatures are available on vibration of uniform, homogeneous and isotropic beams^1-2. Over the years, many researchers have studied vibration of FGM structures. Alshorbagy³ studied free vibration characteristics of a functionally graded beam by finite element method where-in linear strain displacement relationship is considered. Birman⁴ studied the effect of damage on the free and forced vibrations of a functionally graded cantilever beam. Chakraborty⁵ developed a new beam finite element for analysis of functionally graded beams. They considered static, free vibration and wave propagation problems to highlight the behavioural difference of functionally graded material beams with pure metal or a pure ceramic beams. Li⁶ presented a unified approach for analyzing the static and dynamic behaviours of functionally graded Timoshenko and Euler-Bernoulli beams. Sina⁷ presented an analytical method for the free vibration analysis of functionally graded beams. They derived the governing equations of motion using Hamilton’s principle and investigated the effects of boundary conditions, volume fraction and shear deformation on the natural frequencies and mode shapes. Yang⁸studied the free vibration and buckling analyses of functionally graded beams with edge cracks. They conducted a detailed parametric study to show the influences of the location and total number of cracks, material properties, slenderness ratio and end supports on the flexural vibration and buckling characteristics of cracked FGM beams. Gholam⁹ studied free vibrations of functionally graded beams. Mahi¹⁰ proposed an analytical method for temperature dependent free vibration analysis of functionally graded beams with general boundary conditions. Khorramabadi¹¹ studied free vibration of functionally graded beams with piezo-electric layers subjected to axial load.

Although, some literature is available for the study of free vibration of FGM beams, the study pertaining to the effect of transverse shear on the frequencies and mode shapes, to classical boundary conditions, which are encountered during practical engineering use of structures, is not available. In this study, FGM beam with classical boundary conditions such as simply supported at both ends (S-S), clamped at both ends (C-C) and simply supported at one end and clamped at other end (S-C) is considered. The material properties are graded according to power law distribution in the thickness direction with ceramic on top face and metal on bottom face. Free vibration of uniform FGM beams is studied using classical Euler-Bernoulli beam theory and Timoshenko beam theory. The governing equations are obtained by using the principle of virtual work. A linear eigenvalue analysis is carried out for various material distributions and geometrical parameters. The main objective of studying this problem is that, generally an analyst will adopt slender beam analysis neglecting the effects of transverse shear on frequencies and mode shape. The authors have attempted to clearly bring out the significance of considering the transverse shear effects in dynamic analysis of beam type structures.

2.1 Functionally Graded Material

The variation of any material property P such as Young’s modulus E, shear modulus G and density ρ across thickness is assumed to be governed by a power law distribution, as given in Eqn (1), with the co-ordinate z along thickness direction varying between -h/2 to h/2 where h is beam thickness.

$\begin{array}{l}P\left(z\right)=P_c{V}_c+P_m\left(1-V_c\right)\\ V_c={\left(0.5+\frac{z}{h}\right)}^n\end{array}$           (1)

The parameter, V is the volume fraction and subscripts, c and m refer to ceramic and metal respectively. The volume fraction exponent n can take any value between 0 to ∞, with n = 0 and n = ∞ corresponding to the two extremes of completely homogenous ceramic and metal beams respectively.

2.2 Finite Element Analysis

Finite element formulations based on Euler-Bernoulli beam theory and Timoshenko beam theory were developed. The details of governing differential equations and stress-strain relations for both the theories are given in Appendix-A and Appendix-B respectively. In Euler-Bernoulli beam theory, the effect of transverse shear is neglected. However, transverse shear is accounted in Timoshenko beam theory. This transverse shear makes the beam more flexible especially at lower length to thickness ratios.

An Euler-Bernoulli beam element with two nodes per element is considered. Three degrees of freedom are considered at each node. The degrees of freedom for the beam are written as,

$\begin{array}{l}{u}_0\left(x\right)={\displaystyle \sum _{j=1}^2{u}_j{\psi }_j}\\ w_0\left(x\right)={\displaystyle \sum _{j=1}^4{d}_j{\phi }_j}\end{array}\}$           (2)

where ψ_j are the Lagrange shape functions, φ_j are the Hermite shape functions, u₁, u₂ are the nodal axial displacements, d_j = [w₁ θ₁ w₂ θ₂] are the nodal transverse displacements and rotations and subscripts '1' and '2' refer to the nodes 1 and 2.

A Timoshenko beam element with two nodes per element and three degrees of freedom per node is considered. Using Lagrange linear interpolation functions for axial displacement u, lateral displacement w and rotation φ, the degrees of freedom vector for a Timoshenko beam element can be written as

$\begin{array}{l}{u}_0\left(x\right)={\displaystyle \sum _{j=1}^2{u}_j}{\psi }_j\\ w_0\left(x\right)={\displaystyle \sum _{j=1}^2{w}_j}{\psi }_j\\ \phi \left(x\right)={\displaystyle \sum _{j=1}^2{\phi }_j}{\psi }_j\end{array}\}$           (3)

By using the principle of virtual work, the weak forms of the governing equations of motion can be obtained. Substituting the expressions for degree of freedom vector in the weak forms and rearranging, the finite element system of equations can be expressed in a compact form as,

$\left[K\right]\left\{U\right\}={\omega }^2\left[M\right]\left\{U\right\}$           (4)

where {U} is degree of freedom vector, [M] is the mass matrix, [K] is the stiffness matrix and ω denotes circular frequency. This is a linear eigenvalue problem which is solved to obtain the frequencies and mode shapes. In the present study, the phenomenon of shear locking for Timoshenko beam element is alleviated by evaluating stiffness coefficients associated with transverse shear deformation using reduced integration and full integration is used for all other stiffness coefficients.

MATLAB codes are developed based on the above two formulations. A FGM beam with unit thickness, different lengths 'L', and material properties as given in Table 1 is considered.

A convergence study is carried out with different number of elements (NE). Results with 50 elements are found to give converged frequencies to the desired accuracy as given in Table 2. The subsequent study is carried out with 50 numbers of elements.

For validation, the non-dimensional fundamental frequency $\overline{\omega }=\omega L^2\sqrt{\frac{I_0}{E_m{h}^3}}\times {10}^3$ obtained from present Timoshenko beam formulation is compared with obtained from commercial finite element software ANSYS¹³ for completely homogenous ceramic beam. A linear beam element (BEAM3 of ANSYS package version 10.0) with three degrees of freedom (axial and transverse deformation, in-plane rotation) at each node is used. The structure is modeled as an assembly of 50 BEAM3 elements. A modal analysis is carried out with Block Lanczos mode extraction method. The results are given in Table 3, where-in an excellent match is observed for all the L/h ratios considered.

After validating the formulation for homogenous beam, free vibrations of FGM beam is considered. The trend of variation of natural frequencies with various boundary conditions, geometrical parameters and material distribution has been studied. Figure 1 shows the plot of non-dimensional fundamental frequency versus length-to-thickness ratio (L/h) for S-S FGM beam with different volume fraction exponents. The results are plotted for two beam formulations, namely Euler-Bernoulli and Timoshenko. As expected, for the same volume fraction exponent, the frequencies do not match for lower L/h ratios clearly showing the effect of transverse shear by increasing the flexibility thus reducing the frequency compared to that obtained from Euler theory. Further, the nature of vibration response is observed to be similar for all volume fraction exponents considered. The study is repeated with other two boundary conditions (C-C and S-C) as shown in Figs 2 and 3. For all the boundary conditions considered, there is a large difference in the non-dimensional frequency obtained from these two theories for L/h up to 10. For simply supported beam, the curve obtained with Timoshenko beam theory coincides with the curve obtained with Euler beam theory after L/h ≥ 25. For clamped beam, however, these two curves coincide only after L/h ≥ 50. For beam with one end simply supported and other end clamped, the response is observed to be in between that of S-S and C-C beam. Thus, for small values of L/h, transverse shear is more significant in reducing the fundamental frequency than for long and thin beams. As mentioned earlier, this can be attributed to the increase in flexibility of Timoshenko beam for small L/h.

The per cent difference in fundamental frequency obtained using the two formulations for various L/h and volume fraction exponents are tabulated in Tables 4 to 6 for various boundary conditions. As can be observed, for all the boundary conditions, the shear effects are less predominant after L/h ≥ 50. As L/h tends to 100, the difference between frequencies obtained by Euler and Timoshenko beam theory becomes negligible. For identical L/h, the maximum difference is observed for C-C beam as compared to S-S beam and S-C beam. For S-C beam, the difference is in between that of clamped beam and simply supported beam. Further, for a given L/h, the difference is observed to be almost similar for various volume fraction exponents.

Figure 4 shows the effect of volume fraction exponent n on non-dimensional fundamental frequency $\overline{\omega }$ for various L/h ratios for S-S case. For all the ratios of L/h considered, the fundamental frequency reduces with increase in volume fraction exponent. In the present study, since the Young’s modulus of metal (70 GPa) considered is much lower as compared to that of ceramic (151 GPa), higher values of volume fraction exponent will result in lower effective stiffness [Eqns (10 and (10)] of FGM at any cross section of beam, with n = ∞ corresponding to the completely homogenous metal beam. Thus, the fundamental frequency reduces with higher values of volume fraction exponent as the effective stiffness at any cross section reduces compared to completely homogenous ceramic beam (n = 0).

A study is carried out to observe the effect of transverse shear on the fundamental mode for various lengths to thickness ratios and boundary conditions. It is observed that the mode shape obtained using the two theories are identical for all L/h ratios of simply supported FGM beam. Typical mode shapes with volume fraction exponent n = 1.0 for L/h = 5.0 and L/h = 100.0 are shown in Fig. 5. These values of L/h are chosen as they typically represent the maximum and minimum effect of transverse shear on fundamental frequency as shown in Table 4.

Identical mode shapes are also observed for clamped FGM beam with high L/h. However, for smaller values of L/h, a difference is observed in mode shape for clamped FGM beam. Figure 6 shows comparison of fundamental mode obtained using two theories for n = 1.0 and L/h = 5.0, where-in the effect of transverse shear can be clearly observed as the Timoshenko beam shows more flexibility. From the mode shapes plotted in Figs. 5 and 6, it can be concluded that for lower lengths to thickness ratios (e.g. L/h = 5.0), the effect of transverse shear is observed to be predominant for clamped beam as compared to simply supported beam. However, for higher length to thickness ratios, the effect of transverse shear is insignificant on mode shape for all the boundary conditions considered in this analysis.

The study is further extended to higher frequencies obtained from both formulations with various L/h ratios, volume fraction exponents and boundary conditions. Figure 7 shows the variation of first five non-dimensional frequencies obtained from two theories for two typical L/h ratios and volume fraction exponents. It is observed that for short beams, the effect of transverse shear becomes more significant at higher frequencies even for simply supported FGM beam. At higher L/h ratios, however, transverse shear has insignificant effect on higher frequencies for all the boundary conditions considered in this study.

Figure 8 shows a comparison of second mode shape obtained from two theories for C-C FGM beam with volume fraction exponent n = 1 and L/h = 5. Lower value of L/h is chosen to clearly bring out the effect of transverse shear on the higher modes. It can be observed that, similar to the fundamental mode, the effect of transverse shear is to increase the flexibility of the beam as the effective length to thickness reduces at higher modes. This aspect is of more importance as the locations of nodes and anti-nodes predicted by the two theories can be different at higher modes.

The effect of transverse shear in reducing the frequencies, as compared to those obtained from Euler-Bernoulli beam theory, becomes more significant at higher modes. Further, a significant effect of transverse shear is observed on fundamental mode shape of short clamped FGM beam. The effect is observed to be more significant at higher mode shapes. Finally, it is concluded that for the real structures with lower length to thickness ratios, transverse shear effects should be included to predict accurate frequencies and mode shapes for homogenous as well as FGM beams.

1. Rao, S. S. Vibration of continuous system. John Wiley & Sons, Inc., 2007.

2. Thomson, W. T. Theory of vibration. Unwin Hyman Ltd., London, 1990.

3. Alshorbagy, A. E.; Eltaher, M. A.; & Mahmoud, F. F. Free vibration characteristics of a functionally graded beam by finite element method. Appl. Math. Modeling, 2011, 35 (1), 412-25.

4. Birman, V. & Byrd, L. W. Vibrations of damaged cantilever beams manufactured from functionally graded materials. AIAA Journal, 2007, 45(11), 2747-757.

5. Chakraborty, A.; Gopalkrishnan, S. & Reddy, J. N. A new beam finite element for the analysis of functionally graded materials. Int. J. Mechanical Sci., 2003, 45 (3), 519-39.

6. Li, X.-F. A unified approach for analyzing static and dynamic behaviours of functionally graded Timoshenko and Euler-Bernoulli beams. J. Sound Vibration, 2008, 318(4), 1210-229.

7. Sina, S.A.; Navazi, H.M. & Haddadpour, H. An analytical method for free vibration analysis of functionally graded beams. Materials Design, 2009, 30 (3), 741-47.

8. Yang, J. & Chen, Y. Free vibration and buckling analyses of functionally graded beams with edge cracks. Composite Structures, 2008, 83 (1), 48-60.

9.Gholam, R.K. Free vibration of functionally graded beams. World Academy Sci. Engg. Technol., 2011, 74, 366-69.

10. Mahi, A.; Adda Bedia, E.A.; Tounsi, A. & Mechab, I. An analytical method for temperature-dependent free vibration analysis of functionally graded beams with general boundary conditions. Composite Structures, 2010, 92 (8), 1877-887.

11. Khorramabadi, M.K. Free vibration of functionally graded beams with piezo-electric layers subjected to axial load. J. Solid Mechanics, 2009, 1 (1), 22-28.

12. Reddy, J.N. Mechanics of laminated composite plates and shells: theory and analysis. CRC Press, New York, 2003. pp. 131-166.

13. ANSYS Inc. ANSYS package version 10.0, Canons burgh, PA, USA, 2007. Ch. 14.3.

Euler-Bernoulli Beam

The differential equations governing the vibrations of a straight beam are¹²,

$\begin{array}{l}-\frac{\partial N_{xx}}{\partial x}+I_0\frac{{\partial }^2{u}_0}{\partial t^2}-I_1\frac{{\partial }^2}{\partial t^2}\left(\frac{\partial w_0}{\partial x}\right)=0\\ -\frac{{\partial }^2{M}_{xx}}{\partial x^2}-\frac{\partial }{\partial x}\left(N_{xx}\frac{\partial w_0}{\partial x}\right)-q\left(x\right)+I_0\frac{{\partial }^2{w}_0}{\partial t^2}\\ -I_2\frac{{\partial }^2}{\partial t^2}\left(\frac{{\partial }^2{w}_0}{\partial x^2}\right)+I_1\frac{{\partial }^2}{\partial t^2}\left(\frac{\partial u_0}{\partial x}\right)=0\end{array}$            (5)

where u₀ and w₀ are displacements in axial (x) and transverse (z) directions, 't' is time and q(x) is generalized transverse load. N_xx and M_xx are stress and moment resultants in axial direction. Mass moments of inertia terms I_i are given by,

$\left\{\begin{array}{c}{I}_0\\ I_1\\ I_2\end{array}\right\}={\displaystyle \underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}\left\{\begin{array}{c}1\\ z\\ z^2\end{array}\right\}\rho \left(z\right)dz}$            (6)

The linear strain-displacement relation is given by,

$\begin{array}{l}{\epsilon }_{xx}={\epsilon }_{xx}^0+z{k}_{xx}\\ {\epsilon }_{xx}^0=\frac{\partial u_0}{\partial x}\\ k_{xx}=-\frac{{\partial }^2{w}_0}{\partial x^2}\end{array}$           (7)

The stress in the axial direction can be expressed as,

${\sigma }_{xx}=E\left(z\right){\epsilon }_{xx}$           (8)

The stress and moment resultants in the axial direction can be expressed as,

$\begin{array}{l}{N}_{xx}=A_{xx}{\epsilon }_{xx}^0+B_{xx}{\kappa }_{xx}\\ M_{xx}=B_{xx}{\epsilon }_{xx}^0+D_{xx}{\kappa }_{xx}\end{array}$           (9)

where 'A_xx', 'B_xx' and 'D_xx are extensional stiffness, bending-extension coupling stiffness and bending stiffness respectively given by,

$\left\{\begin{array}{c}{A}_{xx}\\ B_{xx}\\ D_{xx}\end{array}\right\}={\displaystyle \underset{-\frac{h}{2}}{\overset{\frac{h}{2}}{\int }}E\left(z\right)}\left\{\begin{array}{c}1\\ z\\ z^2\end{array}\right\}dz$          (10)

Timoshenko Beam

The differential equations governing the vibrations of shear flexible beam based on first order shear deformation theory are given by¹²,

$\begin{array}{l}-\frac{\partial N_{xx}}{\partial x}+I_0\frac{{\partial }^2{u}_0}{\partial t^2}+I_1\frac{{\partial }^2{\phi }_x}{\partial t^2}=0\\ -\frac{\partial Q_x}{\partial x}-\frac{\partial }{\partial x}\left(N_{xx}\frac{\partial w_0}{\partial x}\right)-q(x)+I_0\frac{{\partial }^2{w}_0}{\partial t^2}=0\\ -\frac{\partial M_{xx}}{\partial x}+Q_x+I_2\frac{{\partial }^2{\phi }_x}{\partial t^2}+I_1\frac{{\partial }^2{u}_0}{\partial t^2}=0\end{array}$            (11)

where φ_x is rotation, Q_x is transverse shear and the other terms are as explained in Appendix-A. The linear strain-displacement relations can be written as,

${\epsilon }_{xx}=\frac{\partial u_0}{\partial x}-z\frac{\partial {\phi }_x}{\partial x}$            (12)

and the transverse shear strain is given by,

${\gamma }_{xz}={\phi }_x+\frac{\partial w_0}{\partial x}$            (13)

The axial stress 'σ _xx' and axial strain 'ε_xx' are related by Eqn,

${\sigma }_{xx}=E\left(z\right){\epsilon }_{xx}$            (14)

The shear stress 'τ_xz 'and shear strain 'γ_xz ' are related by Eqn

${\tau }_{xz}=G\left(z\right){\gamma }_{xz}$            (15)

The transverse shear is,

$Q_x=K_s{A}_{55}\left(\frac{d{w}_0}{dx}+{\phi }_x\right)$           (16)

where K_s is shear correction factor (5/6). The shear stiffness is given by,

$A_{55}={\displaystyle \underset{-h/2}{\overset{h/2}{\int }}G(z)dz}$           (17)