Free Vibration Analysis of Functionally Graded Beams

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.

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.

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.

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)=\sum _{j=1}^{2}{u}_{j}{\psi }_{j}\\ {w}_{0}\left(x\right)=\sum _{j=1}^{4}{d}_{j}{\phi }_{j}\end{array}\right\}$           (2)

where ψj are the Lagrange shape functions, φj are the Hermite shape functions, u1, u2 are the nodal axial displacements, dj = [w1 θ1 w2 θ2] 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)=\sum _{j=1}^{2}{u}_{j}{\psi }_{j}\\ {w}_{0}\left(x\right)=\sum _{j=1}^{2}{w}_{j}{\psi }_{j}\\ \phi \left(x\right)=\sum _{j=1}^{2}{\phi }_{j}{\psi }_{j}\end{array}\right\}$           (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.

Table 1. Material properties of ceramic and metal.

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.

Table 2. Convergence study (S-S, n = 0, L/h = 100, Timoshenko theory)

For validation, the non-dimensional fundamental frequency $\overline{\omega }=\omega {L}^{2}\sqrt{\frac{{I}_{0}}{{E}_{m}{h}^{3}}}×{10}^{3}$ obtained from present Timoshenko beam formulation is compared with obtained from commercial finite element software ANSYS13 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.

Table 3. Comparison of $\overline{\omega }$ (Timoshenko beam formulation)

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.

Figure 1.Effect of L/h on non-dimensional fundamental frequency (S-S).

Figure 2.Effect of L/h on non-dimensional fundamental frequency (C-C).

Figure 3.Effect of L/h on non-dimensional fundamental frequency (S-C).

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.

Table 4.Effect of geometrical parameters and material distribution on (S-S)

Table 5.Effect of geometrical parameters and material distribution on (C-C)

Table 6.Effect of geometrical parameters and material distribution on (S-C).

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).

Figure 4.Effect of volume fraction exponent n on non-dimensional fundamental frequency (S-S).

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.

Figure 5. Normalised first mode for S-S FGM beam (n = 1.0).

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.

Figure 6. Normalised first mode for C-C FGM beam (n = 1, L/h = 5).

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 7. Variation of non-dimensional frequencies for different theories.

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.

Figure 8. Normalized second mode for C-C FGM beam (n = 1, L/h = 5).

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.

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 are12,

$\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 u0 and w0 are displacements in axial (x) and transverse (z) directions, 't' is time and q(x) is generalized transverse load. Nxx and Mxx are stress and moment resultants in axial direction. Mass moments of inertia terms Ii are given by,

$\left\{\begin{array}{c}{I}_{0}\\ {I}_{1}\\ {I}_{2}\end{array}\right\}=\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 'Axx', 'Bxx' and 'Dxx 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\}=\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 by12,

$\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\left(x\right)+{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, Qx 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 Ks is shear correction factor (5/6). The shear stiffness is given by,

${A}_{55}=\underset{-h/2}{\overset{h/2}{\int }}G\left(z\right)dz$           (17)