Undamped Free Vibration Analysis of Functionally Graded Beams: A Dynamic Finite Element Approach

Abstract

A Dynamic Finite Element (DFE) method for coupled axial–flexural undamped free vibration analysis of functionally graded beams is developed and subsequently used to investigate the system’s natural frequencies and mode shapes. The formulation is based on the Euler–Bernoulli beam theory and material grading is assumed to follow a power law variation through the thickness direction. Using the closed-form solutions to the uncoupled segments of the system’s governing differential equations as the basis functions of approximation space, the dynamic, frequency-dependent, trigonometric interpolation functions are developed. The interpolation functions are used with the weighted residual method to develop the DFE of the system. The resulting nonlinear eigenvalue problem is then solved to determine the coupled natural frequencies. Example elements using DFE, Finite Element Method (FEM) and the Dynamic Stiffness Method (DSM) are implemented in MATLAB for testing, verification, and validation. Good agreement was observed and the DFE formulation exhibited superior convergence performance compared to the FEM.

1. Introduction

Functionally graded materials (FGM) are a type of composite with continuously varying properties throughout its volume. The concept of FGM was originally developed in the 1980s in an attempt to create a reliable and durable heat shield for spacecraft re-entry [1]. FGM have seen applications in areas such as aerospace, automotive, electronic and biomedical engineering. Traditionally, the use of FGM has been limited due to the manufacturing challenges of creating smooth, continuous, varying properties. Recent advances in manufacturing techniques, such as electron beam freeform fabrication and additive manufacturing, are improving the manufacturability of FGM and could lead to wider use of and interest in the material [2,3,4]. Therefore, there is a growing interest and need for closed form and numerical analysis of FGM beams.

Most recent studies on FGM beams have generally been focused on damped beam models and damaged simulation cases such as [5,6,7,8,9]. Additional studies on undamped free vibrational analysis of FGM is still needed. Various methods and theories have been proposed to analyze the free vibration performance of FGM beams. Numerous studies have used direct analytical methods to derive and solve the governing differential equations of motion [10,11,12]. The finite element method (FEM) is the most popular computational method for solving structural mechanics and has also been explored to analyze the free-vibration characteristics of FGM beams based on Euler–Bernoulli and first-order shear deformation beam theories [13,14,15,16,17]. The Rayleigh–Ritz and Chebyshev collocation methods have also been explored and reported in [18,19,20]. The dynamic stiffness method (DSM) was also applied for free vibration analysis of FGM using both Euler–Bernoulli and Timoshenko beam theories by various authors [21,22,23,24,25].

Of the approaches mentioned above, the FEM and DSM are of particular importance for this study. The FEM is often used when analyzing structural mechanics and its greatest advantage is its robustness and generalizability. It can be used to model complex geometries and loadings. However, the accuracy of the method is dependent on the number of elements used and this becomes increasingly important when analyzing high frequency mode shapes, where increasing number of elements is required to produce accurate results. The DSM on the other hand is an analytical method similar to FEM, but the results of this method are exact within the limits of the theory. This method solves the system’s governing equations by combining the two coupled differential equations of motions into a single, higher-order ordinary differential equation. The general close-form solution is then used to obtain a frequency-dependent stiffness matrix [26,27]. A single DSM element can produce an infinite number of natural frequencies and mode shapes. Although the DSM method provides accurate results even at higher-frequency modes, it does not have the generality of the FEM and is often limited when dealing with complex geometry.

The Dynamic Finite Element (DFE) method is a hybrid method introduced that combines the strengths of DSM and FEM methods [28]. Like the FEM, the DFE formulation uses the general procedure of the weighted residual method. However, instead of using polynomial shape functions, the DFE shape functions are frequency-dependent, trigonometric expressions obtained from the basis functions of approximation space, satisfying the uncoupled portions of the governing differential equations. The solution of the DFE method is a frequency-dependent stiffness matrix like DSM, where the inertial and stiffness properties are combined into a single dynamic (frequency-dependent) stiffness matrix. When compared to the FEM, DFE has improved convergence and accuracy due to the use of Dynamic Trigonometric Shape Functions (DTSFs). The superior convergence performance is particularly apparent at higher frequencies.

The DFE method has been successfully extended to the vibration analysis of various intact and defective beam-like structures [29,30,31,32,33,34,35,36]. Formulations included various coupled/combined loading cases and layered sandwich beam models. In all cases, either homogenous material properties were assumed through the beam thickness or the homogenization method to evaluate apparent/equivalent properties. However, the DFE method has not been extended to free vibrational analysis of non-homogenous materials such as FGM beams.

In the present study, the Dynamic Finite Element (DFE) method is extended to the free vibration analysis of FGM Euler–Bernoulli beams, under various boundary conditions and material variation. In the following section, the material variation model and the mathematical procedure of the Euler–Bernoulli DFE formulation is presented. A one-element DFE model is then developed and validated against classical, FEM and DSM results for homogeneous and FGM beams. The efficiency and accuracy of a one-element DFE model is then demonstrated at various boundary conditions and with material variations. Materials and methods are presented in Section 2, the results are provided in Section 3, followed by final conclusions summarizing the most important achievements of the presented work.

2. Materials and Methods

The coordinate system and notation used is shown in Figure 1, where the beam has a length L, thickness h, and width b. Material properties are Young’s Modulus E and mass density $\mathsf{\rho }$.

The effective material property, P, for a FGM is evaluated using the following rule of mixtures, expressed as [9,21,22]:

$$P\left(z\right)=\left(P_t-P_b\right)V_t+P_b$$

(1)

where P_t and P_b are the material properties at the top and bottom surfaces of the beam, respectively. The effective material properties are assumed to vary according to a power law distribution. V_t is the volume fraction of the top constituent of the beam, defined as [9,21,22]:

$$V_t={\left(\frac{z}{h}+\frac{1}{2}\right)}^k$$

(2)

Variation of the volume fraction against thickness for different values, k, is represented in Figure 2, where k = 1 indicates a linear variation of material composition from the top to the bottom, k = 0 represents a beam composed entirely of the top material, and k = ∞ characterizes a beam composed entirely of the bottom material.

To derive the differential equations governing the free vibrations of Euler–Bernoulli FGM beams, Su & Banerjee established the expressions for the system’s strain and kinetic energy and applied Hamilton’s principle. The resulting coupled differential equations of motion are [21]:

$$-I_0\ddot{v}+{\mathrm{I}}_1\ddot{{\mathrm{w}}^{\prime }}+{\mathrm{A}}_0{\mathrm{v}}^{″}-{\mathrm{A}}_1{\mathrm{w}}^{‴}=0$$

(3)

$$-{\mathrm{I}}_0\ddot{\mathrm{w}}-{\mathrm{I}}_1\ddot{{\mathrm{v}}^{\prime }}+{\mathrm{I}}_2\ddot{{\mathrm{w}}^{″}}+{\mathrm{A}}_1{\mathrm{v}}^{‴}-{\mathrm{A}}_2{\mathrm{w}}^{″″}=0$$

(4)

where (′) stands for the spatial derivative with respect to the beam’s longitudinal axis, y, where 0 ≤ y ≤ L, and (˙) denotes the derivative with respect to t (time). The parameters I_i and A_i are defined as:

$$\begin{array}{c}{I}_i=\int z^i\mathsf{\rho }\left(z\right)dA\\ A_i=\int z^{i}E\left(z\right)dA\end{array}$$

(5)

Assuming simple harmonic motion, axial and lateral displacements, $v\left(y,t\right)$ and $w\left(y,t\right)$, respectively, both functions of space (y) and time (t), can be written as (i.e., separation of variables):

$$\begin{array}{c}v\left(y,t\right)=V\left(y\right)e^{i\mathsf{\omega }t}\\ w\left(y,t\right)=W\left(y\right)e^{i\mathsf{\omega }t}\end{array}$$

(6)

where ω denotes the frequency of vibration and V(y) and W(y) are the amplitudes of axial and flexural displacements, respectively. By substituting expressions (6) into Equations (3) and (4), the governing equations can be rewritten as:

$$I_{0}V{\mathsf{\omega }}^2-I_1{W}^{\prime }{\mathsf{\omega }}^2+A_0{V}^{″}-A_1{W}^{‴}=0$$

(7)

$$I_{0}W{\mathsf{\omega }}^2+I_1{V}^{\prime }{\mathsf{\omega }}^2-I_2{W}^{″}{\mathsf{\omega }}^2+A_1{V}^{‴}-A_2{W}^{″″}=0$$

(8)

where the space dependency (y) of amplitudes has been omitted for brevity, and both expressions have been divided by the non-zero term, $e^{i\mathsf{\omega }t}$.

The boundary conditions are produced as a by-product of the governing differential equation derivation. The loads or natural boundary conditions are the resultant axial force (F), shear force (S) and bending moment (M), defined as:

$$\begin{array}{c}F=-A_0{V}^{\prime }+A_1{W}^{″}\\ S=-I_{1}V{\mathsf{\omega }}^2-A_1{V}^{″}+I_2{W}^{\prime }{\mathsf{\omega }}^2+A_2{W}^{‴}\\ M=A_1{V}^{\prime }-A_2{W}^{″}\end{array}$$

(9)

The initial steps in a DFE formulation follow similar derivation as conventional FEM. First, by implementing the Galerkin-type method of weighted residuals, Equations (7) and (8) are formulated into residual equations weighted by the virtual axial and lateral displacements, δv and δw, respectively. The integral value of weighted residuals is then set to zero, resulting in the following integral form of the system’s governing equations:

$$W_v={{\displaystyle \int }}_0^L\mathsf{\delta }v\left(I_{0}V{\mathsf{\omega }}^2-I_1{W}^{\prime }{\mathsf{\omega }}^2+A_0{V}^{″}-A_1{W}^{‴}\right)\mathsf{\delta }y=0$$

(10)

$$W_w={{\displaystyle \int }}_0^L\mathsf{\delta }w\left(I_{0}W{\mathsf{\omega }}^2+I_1{V}^{\prime }{\mathsf{\omega }}^2-I_2{W}^{″}{\mathsf{\omega }}^2+A_1{V}^{‴}-A_2{W}^{″″}\right)\mathsf{\delta }y=0$$

(11)

A set of integration by parts is performed, leading to the following weak form of the above integral equations:

$$\begin{array}{l}{W}_v=\underset{\_}{\mathsf{\delta }v\left(-A_0{V}^{\prime }+A_1{W}^{″}\right)}{|}_0^L\\ +{{\displaystyle \int }}_0^L\left(\left(-{\mathsf{\omega }}^2\right)\mathsf{\delta }v\left(I_{0}V-I_1{W}^{\prime }\right)+\mathsf{\delta }{v}^{\prime }\left(A_0{V}^{\prime }+A_1{W}^{″}\right)\right)\mathsf{\delta }y=0\end{array}$$

(12)

$$\begin{array}{l}{W}_w={\underset{\_}{\left(\mathsf{\delta }w\left(-I_1{\mathsf{\omega }}^{2}V-A_1{V}^{″}+I_2{\mathsf{\omega }}^2{W}^{\prime }+A_2{W}^{‴}\right)+\mathsf{\delta }{w}^{\prime }\left(A_1{V}^{\prime }-A_2{W}^{″}\right)\right)|}}_0^L\\ +{{\displaystyle \int }}_0^L(\left(-{\mathsf{\omega }}^2\right)\mathsf{\delta }w\left(I_{0}W\right)+\left(-{\mathsf{\omega }}^2\right)\mathsf{\delta }{w}^{\prime }\left(I_2{W}^{\prime }-I_{1}V\right)\\ +\mathsf{\delta }{w}^{″}\left(A_2{W}^{″}-A_1{V}^{\prime }\right))\mathsf{\delta }y=0\end{array}$$

(13)

The above expressions also satisfy the principle of virtual work,

$$W=W_{int}-W_{ext}=0$$

(14)

where W represents the total virtual work, W_int stands for the internal virtual work, and W_ext is the external virtual work. The underlined terms in expressions (12) and (13) represent the virtual work generated by the load boundary conditions of the system, presented earlier as expressions (9). For the free vibration of a system, the total external work is null, W_ext = 0, and therefore:

$$\begin{array}{c}\mathsf{\delta }v\left(-A_0{V}^{\prime }+A_1{W}^{″}\right){|}_{x=0}=\mathsf{\delta }{W}_{ext}{|}_{x=0}\\ \mathsf{\delta }v\left(-A_0{V}^{\prime }+A_1{W}^{″}\right){|}_{x=L}={\mathsf{\delta }W}_{\mathit{ext}}{|}_{x=L}\\ \left(\mathsf{\delta }w\left(-I_1{\mathsf{\omega }}^{2}V–A_1{V}^{″}+I_2{\mathsf{\omega }}^2{W}^{\prime }+A_2{W}^{‴}\right)+\mathsf{\delta }{w}^{\prime }\left(A_1{V}^{\prime }–A_2{W}^{″}\right)\right){|}_{x=0}\\ ={\mathsf{\delta }W}_{\mathit{ext}}{|}_{x=0}\\ \left(\mathsf{\delta }w\left(-I_1{\mathsf{\omega }}^{2}V-A_1{V}^{″}+I_2{\mathsf{\omega }}^2{W}^{\prime }+A_2{W}^{‴}\right)+\mathsf{\delta }{w}^{\prime }\left(A_1{V}^{\prime }-A_2{W}^{″}\right)\right){|}_{x=L}\\ =\mathsf{\delta }{W}_{ext}{|}_{x=L}\end{array}$$

(15)

Consequently, the net resultant virtual work caused by the external force and moment terms in (12) and (13) goes to zero. With the boundary conditions satisfied, the remaining parts of the resulting expressions can now be discretized, which means the beam length (L) is divided into several elements, leading to the following elemental form:

$$W_v^k={{\displaystyle \int }}_0^l\left(\left(-{\mathsf{\omega }}^2\right)\mathsf{\delta }v\left(I_{0}V-I_1{W}^{\prime }\right)+\mathsf{\delta }{v}^{\prime }\left(A_0{V}^{\prime }-A_1{W}^{″}\right)\right)\mathsf{\delta }y$$

(16)

$$\begin{array}{c}{W}_w^k={{\displaystyle \int }}_0^l(\left(-{\mathsf{\omega }}^2\right)\mathsf{\delta }w\left(I_{0}W\right)+\left(-{\mathsf{\omega }}^2\right)\left(\mathsf{\delta }{w}^{\prime }\right)\left(I_2{W}^{\prime }-I_{1}V\right)\\ +\mathsf{\delta }{w}^{″}\left(A_2{W}^{″}-A_1{V}^{\prime }\right))\mathsf{\delta }y\end{array}$$

(17)

The above integral equations $W_v^k$ and $W_w^k$, respectively, stand for the flexural and torsional contributions of each individual element ‘k’, of length ‘l’, in the total virtual work of the entire system, $W_v$ and $W_w.$At this point, the DFE formulation diverges from the conventional FEM derivation. Instead of using polynomial interpolation (shape) functions to express the field variables, V and W, in terms of nodal variables (i.e., conventional FEM), further manipulation is performed before DTSFs are applied, leading to the elemental dynamic stiffness matrix.

The closed-form solutions of the uncoupled governing differential equations are used as the basis functions of approximation space to derive the dynamic shape functions, which are then exploited to find the element (frequency-dependent) dynamic stiffness matrix.

The following term, ξ = y/l (0 ≤ ξ ≤ 1), is introduced into Equations (16) and (17), resulting in the following non-dimensional form:

$$W_v^k={{\displaystyle \int }}_0^1\mathsf{\delta }v\left(I_{0}lV{\mathsf{\omega }}^2-I_1{W}^{\prime }{\mathsf{\omega }}^2+\frac{A_0}{l}{V}^{″}-\frac{A_1}{l^2}{W}^{‴}\right)\mathsf{\delta }\mathsf{\xi }$$

(18)

$$W_w^k={{\displaystyle \int }}_0^l\mathsf{\delta }w\left(I_{0}lW{\mathsf{\omega }}^2+I_1{V}^{\prime }{\mathsf{\omega }}^2-\frac{I_2}{l}{W}^{″}{\mathsf{\omega }}^2+\frac{A_1}{l^2}{V}^{‴}-\frac{A_2}{l^3}{W}^{″″}\right)\mathsf{\delta }\mathsf{\xi }$$

(19)

Applying additional integration by parts to the above discretized (elemental) weak-form equations, leads to the following (equivalent) form of the equations:

$$\begin{array}{r}{W}_v^k=\left[\mathsf{\delta }{v}^{\prime }\left(\frac{A_0}{l}V\right)\right]{|}_0^1+\underset{*}{\underbrace{{{\displaystyle \int }}_0^{1}V\left(-{\mathsf{\omega }}^2{I}_{0}l\mathsf{\delta }v-\frac{A_0}{l}\mathsf{\delta }{v}^{″}\right)d\mathsf{\xi }}}\\ +\underset{\mathrm{Coupling}}{\underbrace{{\mathsf{\omega }}^2{{\displaystyle \int }}_0^1\mathsf{\delta }v\left(I_1{W}^{\prime }\right)d\mathsf{\xi }-{{\displaystyle \int }}_0^1\mathsf{\delta }{v}^{\prime }\left(\frac{A_1}{l^2}{W}^{″}\right)d\mathsf{\xi }}}\end{array}$$

(20)

$$\begin{array}{r}{W}_w^k=\left[\mathsf{\delta }{w}^{\prime }\left(-\frac{{\mathsf{\omega }}^2{I}_2}{l}W\right)+\mathsf{\delta }{w}^{″}\left(\frac{A_2}{l^3}{W}^{\prime }\right)-\mathsf{\delta }{w}^{‴}\left(\frac{A_2}{l^3}w\right)\right]{|}_0^1\\ +\underset{**}{\underbrace{{{\displaystyle \int }}_0^{1}W\left(\frac{A_2}{l^3}\mathsf{\delta }{w}^{″″}+\frac{{\mathsf{\omega }}^2{I}_2}{l}\mathsf{\delta }{w}^{″}-{\mathsf{\omega }}^2{I}_{0}l\mathsf{\delta }w\right)d\mathsf{\xi }}}\\ +\underset{\mathrm{Coupling}}{\underbrace{{\mathsf{\omega }}^2{{\displaystyle \int }}_0^1\mathsf{\delta }{w}^{\prime }\left(I_{1}V\right)d\mathsf{\xi }-{{\displaystyle \int }}_0^1\mathsf{\delta }{w}^{″}\left(\frac{A_1}{l^2}{V}^{\prime }\right)d\mathsf{\xi }}}\end{array}$$

(21)

each consisting of two integral (uncoupled and coupled) terms, and one boundary term. The coupling terms are separated and labeled here for clarity. The closed-form solutions of the uncoupled integral terms (*) and (**), respectively, can be written as:

$$\begin{array}{c}V\left(\mathsf{\xi }\right)=D_{1}cos\left(\mathsf{\gamma }\mathsf{\xi }\right)+D_{2}sin\left(\mathsf{\gamma }\mathsf{\xi }\right)\\ W\left(\mathsf{\xi }\right)=C_{1}cos\left(\mathsf{\alpha }\mathsf{\xi }\right)+C_{2}sin\left(\mathsf{\alpha }\mathsf{\xi }\right)+C_{3}cosh\left(\mathsf{\beta }\mathsf{\xi }\right)+C_{4}sinh\left(\mathsf{\beta }\mathsf{\xi }\right)\end{array}$$

(22)

where C_1,2,3,4 and D_1,2 are constants and

$$\mathsf{\gamma }=\sqrt{\frac{-{\mathsf{\omega }}^2{I}_0{l}^2}{A_0}}$$

(23)

$$\mathsf{\alpha },\mathsf{\beta }=\frac{1}{\sqrt{2\mathrm{A}}}\sqrt{\left\{-\mathrm{B}\pm \sqrt{{\mathrm{B}}^2–4\mathrm{AC}}\right\}}$$

(24)

with

$$\begin{array}{c}A=A_2\\ B={\mathsf{\omega }}^2{I}_2\\ C=-{\mathsf{\omega }}^2{I}_0\end{array}$$

(25)

The interpolation functions, also referred to as shape functions of approximation space, are then obtained as follows. First, using the generalized parameters $\langle a\rangle$, $\langle \delta a\rangle$, $\langle b\rangle$, and $\langle \delta b\rangle$ for the field variables (i.e., solution functions), V and W, and virtual displacements (i.e., test functions), δV and δW, can be written as:

$$\begin{array}{c}V=\langle P{\left(\mathsf{\xi }\right)\rangle }_v\left\{a\right\}\\ \mathsf{\delta }V=\langle P{\left(\mathsf{\xi }\right)\rangle }_v\left\{\mathsf{\delta }a\right\}\\ \mathrm{W}={\langle \mathrm{P}\left(\mathsf{\xi }\right)\rangle }_{\mathrm{w}}\left\{\mathrm{b}\right\}\\ \mathsf{\delta }\mathrm{W}=\langle \mathrm{P}{\left(\mathsf{\xi }\right)\rangle }_{\mathrm{w}}\left\{\mathsf{\delta }\mathrm{b}\right\}\end{array}$$

(26)

The axial and flexural dynamic basis functions of approximation space, respectively, also reported in [28], are defined as:

$$V_n\left(\mathsf{\xi }\right)=\langle cos\left(\mathsf{\gamma }\mathsf{\xi }\right) \frac{sin\left(\mathsf{\gamma }\mathsf{\xi }\right)}{\mathsf{\gamma }}\rangle$$

(27)

$$W_n\left(\mathsf{\xi }\right)=\langle cos\left(\mathsf{\alpha }\mathsf{\xi }\right) \frac{sin\left(\mathsf{\alpha }\mathsf{\xi }\right)}{\mathsf{\alpha }} \frac{cosh\left(\mathsf{\beta }\mathsf{\xi }\right)-cos\left(\mathsf{\alpha }\mathsf{\xi }\right)}{{\mathsf{\alpha }}^2+{\mathsf{\beta }}^2} \frac{sinh\left(\mathsf{\beta }\mathsf{\xi }\right)-sin\left(\mathsf{\alpha }\mathsf{\xi }\right)}{{\mathsf{\alpha }}^3+{\mathsf{\beta }}^3}\rangle$$

(28)

Replacing the generalized parameters $\langle a\rangle$, $\langle \delta a\rangle$, $\langle b\rangle$, and $\langle \delta b\rangle$ with the nodal variables, $\begin{array}{cc}\langle V_1& V_2\rangle \end{array}$, $\begin{array}{cc}\langle V_1^{\prime }& V_2^{\prime }\rangle \end{array}$, $\begin{array}{cccc}\langle W_1& W_1^{\prime }& W_2& W_2^{\prime }\rangle \end{array}$, and $\begin{array}{cccc}\langle \delta W_1& \delta W_1^{\prime }& \delta W_2& \delta W_2^{\prime }\rangle \end{array}$, respectively, expression (26) can be rewritten as:

$$\begin{array}{c}\left\{V_N\right\}={\left[P_n\right]}_v\left\{a\right\}\\ \left\{\mathsf{\delta }{V}_n\right\}={\left[P_n\right]}_v\left\{\mathsf{\delta }a\right\}\\ \left\{W_n\right\}={\left[P_n\right]}_w\left\{b\right\}\\ \left\{\mathsf{\delta }{W}_n\right\}={\left[P_n\right]}_w\left\{\mathsf{\delta }b\right\}\end{array}$$

(29)

where the matrices [P_n]_v and [P_n]_w are defined as:

$${\left\{P_n\right\}}_v={\left\{\begin{array}{c}{P}_1\left(\xi =0\right)\\ P_2\left(\xi =1\right)\end{array}\right\}}_v=\left[\begin{array}{cc}1& 0\\ \mathit{cos}\left(\gamma \right)& \frac{\mathit{sin}\left(\gamma \right)}{\gamma }\end{array}\right]$$

(30)

$$\begin{array}{l}{\left\{P_n\right\}}_v={\left\{\begin{array}{c}{P}_1\left(\xi =0\right)\\ {P_1}^{\prime }\left(\xi =0\right)\\ P_2\left(\xi =1\right)\\ {P_2}^{\prime }\left(\xi =1\right)\end{array}\right\}}_w\\ =\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& \frac{\beta -\alpha }{{\alpha }^3+{\beta }^3}\\ \mathit{cos}\left(\gamma \right)& \frac{\mathit{sin}\left(\gamma \right)}{\gamma }& \frac{\mathit{cosh}\left(\beta \right)-\mathit{cos}\left(\alpha \right)}{a_i^2+{\beta }^2}& \frac{\mathit{sinh}\left(\beta \right)-\mathit{sin}\left(\alpha \right)}{a_i^3+{\beta }^3}\\ -\alpha \mathit{sin}\left(\alpha \right)& \mathit{cos}\left(\alpha \right)& \frac{\beta \mathit{sinh}\left(\beta \right)+\alpha \mathit{sin}\left(\alpha \right)}{{\alpha }^2+{\beta }^2}& \frac{\beta \mathit{cosh}\left(\beta \right)-\alpha \mathit{cos}\left(\alpha \right)}{{\alpha }^3+{\beta }^3}\end{array}\right]\end{array}$$

(31)

By combining (26), (30), and (31) the approximation in terms of nodal variables can be written as:

$$\begin{array}{c}V\left(\xi \right)=\langle P{\left(\xi \right)\rangle }_v{\left[P_n\right]}_v^{-1}\left\{V_n\right\}=\langle N{\left\{\xi \right\}\rangle }_v\left\{V_n\right\}\\ W\left(\xi \right)=\langle P{\left(\xi \right)\rangle }_w{\left[P_n\right]}_w^{-1}\left\{W_n\right\}=\langle N{\left\{\xi \right\}\rangle }_w\left\{W_n\right\}\end{array}$$

(32)

where $\langle N{\left(\xi \right)\rangle }_v$ and $\langle N{\left(\xi \right)\rangle }_w$ are the frequency–dependent trigonometric shape functions for the axial and flexure, respectively. Expression (31) can be combined and rewritten as:

$$\left\{\begin{array}{c}V\left(\xi \right)\\ W\left(\xi \right)\end{array}\right\}=\left[N\right]\left\{u_n\right\}$$

(33)

where

$$\begin{array}{c}\left[N\right]=\left[\begin{array}{cccccc}{N}_{1v}\left(\omega \right)& 0& 0& N_{2v}\left(\omega \right)& 0& 0\\ 0& N_{1w}\left(\omega \right)& N_{2w}\left(\omega \right)& 0& N_{3w}\left(\omega \right)& N_{4w}\left(\omega \right)\end{array}\right]\\ \left\{u_n\right\}={\begin{array}{cccccc}{V}_1& W_1& W_1^{\prime }& V_2& W_2& W_2^{\prime }\end{array}}^T\end{array}$$

(34)

The following shape functions, also presented in [30], are used to express the approximate elemental axial displacement in terms of nodal displacements, and along the element domain k, of length l [28,32]:

$$N_{v_1}=\mathit{cos}\left(\gamma \xi \right)-\frac{\mathit{cos}\left(\gamma \right)\mathit{sin}\left(\gamma \xi \right)}{D_t}$$

(35)

$$N_{v_2}=\frac{\mathit{sin}\left(\gamma \xi \right)}{D_t}$$

(36)

where

$$D_t=\mathit{sin}\left(\gamma \right)$$

(37)

The following shape functions are used to state the approximate elemental bending displacement in terms of nodal displacements, along the domain [28,30]:

$$\begin{array}{c}{N}_{w_1}=\frac{\alpha \beta }{D_f}\{-\mathit{cos}\left(\alpha \xi \right)+\mathit{cos}\left(\alpha \left(1-\xi \right)\right)\mathit{cosh}\left(\beta \right)+\mathit{cos}\left(\alpha \right)\mathit{cosh}\left(\beta \left(1-\xi \right)\right)\\ -\mathit{cosh}\left(\beta \xi \right)-\frac{\beta }{\alpha }\mathit{sin}\left(\alpha \left(1-\xi \right)\right)\mathit{sinh}\left(\beta \right)+\frac{\alpha }{\beta }\mathit{sin}\left(\alpha \right)\mathit{sinh}\left(\beta \left(1-\xi \right)\right)\}\end{array}$$

(38)

$$\begin{array}{c}{N}_{w_2}=\frac{1}{D_f}\{\beta \left[\mathit{cosh}\left(\beta \left(1-\xi \right)\mathit{sin}\left(\alpha \right)\right)-\mathit{cosh}\left(\beta \right)\mathit{sin}\left(\alpha \left(1-\xi \right)\right)-\mathit{sin}\left(\alpha \xi \right)\right]\\ +\alpha \left[\mathit{cos}\left(\alpha \left(1-\xi \right)\right)\mathit{sinh}\left(\beta \right)-\mathit{cos}\left(\alpha \right)\mathit{sinh}\left(\beta \left(1-\xi \right)\right)-\mathit{sinh}\left(\beta \xi \right)\right]\}\end{array}$$

(39)

$$\begin{array}{c}{N}_{w_3}=\frac{\alpha \beta }{D_f}\{-\mathit{cos}\left(\alpha \left(\xi -l\right)\right)+\mathit{cos}\left(\alpha \xi \right)\mathit{cosh}\left(\beta \right)-\mathit{cosh}\left(\beta \left(1-\xi \right)\right)\\ +\mathit{cos}\left(\alpha \right)\mathit{cosh}\left(\beta \xi \right)-\frac{\beta }{\alpha }\mathit{sin}\left(\alpha \xi \right)\mathit{sinh}\left(\beta \right)+\frac{\alpha }{\beta }\mathit{sin}\left(\alpha \right)\mathit{sinh}\left(\beta \xi \right)\}\end{array}$$

(40)

$$\begin{array}{c}{N}_{w_4}=\frac{1}{D_f}\{\beta \left[-\mathit{cosh}\left(\beta \xi \right)\mathit{sin}\left(\alpha \right)+\mathit{sin}\left(\alpha \left(1-\xi \right)\right)+\mathit{cosh}\left(\beta \right)\mathit{sin}\left(\alpha \xi \right)\right]\\ +\alpha \left[-\mathit{cos}\left(\alpha \xi \right)\mathit{sinh}\left(\beta \right)+\mathit{sinh}\left(\beta \left(1-\xi \right)\right)+\mathit{cos}\left(\alpha \right)\mathit{sinh}\left(\beta \xi \right)\right]\}\end{array}$$

(41)

where

$$D_f=\alpha \beta \left\{-2\left(1-\mathit{cos}\left(\alpha \right)\mathit{cosh}\left(\beta \right)\right)+\left(\frac{{\alpha }^2-{\beta }^2}{\alpha \beta }\right)\mathit{sin}\left(\alpha \right)\mathit{sinh}\left(\beta \right)\right\}$$

(42)

Using expressions (20) and (21) and the shape functions (35) through (42), the element dynamic stiffness matrix, ${\left[K\left(\omega \right)\right]}^k$, is derived. This consists of uncoupled and coupled dynamic stiffness matrices, ${\left[K\left(\omega \right)\right]}_u^k$ and ${\left[K\left(\omega \right)\right]}_c^k$, respectively. There are four coupled and four uncoupled matrix components. The element dynamic stiffness matrix, ${\left[K\left(\omega \right)\right]}^k$, is formulated by adding the eight coupled and uncoupled submatrices.

The system’s global dynamic stiffness matrix, $\left[K\left(\omega \right)\right]$, is then obtained by assembling all the element matrices and applying the boundary conditions. The nonlinear eigenvalue problem, resulting from this method, is:

$$\left[K\left(\omega \right)\right]\left\{U_n\right\}=\left\{0\right\}$$

(43)

where $\left\{U_n\right\}$ is the vector of global nodal displacements of the system. The natural frequencies of the system would be the values of $\mathsf{\omega }$, which yields a zero determinant, $\left|K\left(\omega \right)\right|=0$, for the global dynamic stiffness matrix [28].

It is worth noting that the basis functions are specifically chosen so that when the frequency, $\mathsf{\omega }$, approaches zero, the roots, α, β, and γ of the characteristic equations (22) also approach zero. Subsequently, the dynamic trigonometric shape functions and basis functions become identical to those generated from the Hermite approximations, which are commonly used in conventional beam FEM. The form of the DFE basis functions were formulated to exhibit this property to ensure a complete solution at all frequencies. In other words, if this was not considered, a static deformation solution would not be possible.

Two illustrative examples are presented in this study to validate the presented DFE formulation through free vibration analysis. First, the DFE Euler–Bernoulli formulation is validated against the well-defined solution for solid beams. In this example, the beam is made of aluminum and results are compared against classical methods, FEM and the DSM method for Euler–Bernoulli beams. A pure aluminum beam can be considered a special case of an FGM beam. Using the current model, a solid aluminum beam would require an unattainable value of $k=\infty$ (refer to Figure 2). To avoid the use of infinity as a value two methods have been proposed, either setting the top and bottom properties within ±0.0001% or using a value of $k={10}^6$ [21,22]. It was observed that either method would produce the same results, and for the purposes of this investigation the first method was used. The material properties for the solid aluminum beam are: ${\rho }_t=2700.0027 \mathrm{kg}/{\mathrm{m}}^3$, ${\rho }_b=2700 \mathrm{kg}/{\mathrm{m}}^3$, $E_t=70.00007 \mathrm{GPa}$, and $E_b=70 \mathrm{GPa}$. An in-house DSM code was created in MATLAB to validate against reference [21], where various material grading, boundary conditions (clamped–free, simply supported, clamped–clamped, clamped–pinned) and slender ratios were investigated.

The following non-dimensional parameter was used to normalize and compare results reported in different sources:

$${\lambda }_i={\omega }_i{L}^2\sqrt{\frac{{\rho }_{b}A}{E_{b}I}}$$

(44)

To further validate the proposed DFE formulation, with coupling terms in effect, a second study is done with an FGM beam ($k=0.3)$composed of aluminum and alumina ($A{l}_2{O}_3)$, as also presented in Ref. [9]. The properties of the aluminum are: $E_b=70 \mathrm{GPa}$, ${\rho }_b=2700 \mathrm{kg}/{\mathrm{m}}^3$ and the properties of the alumina are: $E_t=380 \mathrm{GPa}$, ${\rho }_t=3800 \mathrm{kg}/{\mathrm{m}}^3$ [31]. Simsek [9] used a governing equation developed based on the third-order shear deformation theory, and the solutions for both Timoshenko and Euler–Bernoulli beams were presented. The frequency results were non-dimensionalized using the following parameter [9]:

$${\lambda }_i=\frac{{\omega }_i{L}^2}{h}\sqrt{\frac{I_o}{A_o}}$$

(45)

with $I_o$ and $A_o$defined previously in (5). The natural frequencies, $\omega$, were calculated by manipulating Equation (45), and then using Equation (44) converted into a consistent non-dimensional form.

Finally, free vibration analysis was performed to find the first four natural frequencies and mode shapes of a series of FGM beams. The beams are composed of steel and alumina for various $k$values and at a slenderness ratio of $L/h=100$. The properties of the steel are: $E_b=210\mathrm{GPa}$, ${\rho }_b=7800 \mathrm{kg}/{\mathrm{m}}^3$and the properties of the alumina are: $E_t=390 \mathrm{GPa}$, $rh{o}_t=3960 \mathrm{kg}/{\mathrm{m}}^3$, as also reported in [9]. A comparison is made between the convergence performance of DFE and FEM formulations.

3. Results

The convergence test results and performance for the in-house FEM code is presented in Figure 3. For the 3rd mode of a pure aluminum beam, the conventional FEM converges at around 10–12 elements. In contrast, DFE results show exact solutions at a single element (

Table 1. The first three non-dimensional frequencies of a pure aluminum beam for a range of slenderness ratios and simply supported boundary conditions.

Frequency No. i	L/h	$\mathbf{Non-Dimensional}\mathbf{Fundamental}\mathbf{Natural}\mathbf{Frequency}\left({\mathit{\lambda }}_{\mathit{i}}\right)$
		Classical Results [27]	DSM Ref [13]	DSM Code	DFE (1 ELE)	DFE (3 ELE)	FEM (10 ELE)
1	10	9.8696	9.8293	9.8293	9.8293	9.8293	9.8293
	30	9.8696	9.8651	9.8651	9.8651	9.8651	9.8652
	100	9.8696	9.8692	9.8692	9.8692	9.8692	9.8693
2	10	39.478	38.845	38.845	38.845	38.845	38.849
	30	39.478	39.406	39.406	39.406	39.406	39.411
	100	39.478	39.472	39.472	39.472	39.472	39.476
3	10	88.826	85.711	85.711	85.711	85.711	85.757
	30	88.826	88.463	88.463	88.463	88.462	88.511
	100	88.826	88.794	88.794	88.794	88.794	88.841

,

Table 2. The first three non-dimensional frequencies of a pure aluminum beam for a range of slenderness ratios and clamped–clamped boundary conditions.

Frequency No. i	L/h	$\mathbf{Non-Dimensional}\mathbf{Fundamental}\mathbf{Natural}\mathbf{Frequency}\left({\mathit{\lambda }}_{\mathit{i}}\right)$
		Classical Results [27]	DSM Ref [12]	DSM Code	DFE (1 ELE)	DFE (3 ELE)	FEM (10 ELE)
1	10	22.373	22.259	22.259	22.259	22.259	22.260
	30	22.373	22.361	22.361	22.361	22.361	22.361
	100	22.373	22.372	22.372	22.373	22.372	22.373
2	10	61.673	60.522	60.522	60.522	60.522	60.538
	30	61.673	61.542	61.542	61.542	61.542	61.558
	100	61.673	61.673	61.661	61.661	61.660	61.661
3	10	120.90	116.21	116.21	116.21	116.21	116.32
	30	120.90	120.35	120.35	120.35	120.35	120.47
	100	120.90	120.85	120.85	120.86	120.85	120.97

,

Table 3. The first three non-dimensional frequencies of a pure aluminum beam for a range of slenderness ratios and clamped–free boundary conditions.

Frequency No. i	L/h	$\mathbf{Non-Dimensional}\mathbf{Fundamental}\mathbf{Natural}\mathbf{Frequency}\left({\mathit{\lambda }}_{\mathit{i}}\right)$
		Classical Results [27]	DSM Ref [12]	DSM Code	DFE (1 ELE)	DFE (3 ELE)	FEM (10 ELE)
1	10	3.5160	3.5092	3.5092	3.5092	3.5092	3.5092
	30	3.5160	3.5153	3.5153	3.5153	3.5153	3.5153
	100	3.5160	3.5160	3.5159	3.5159	3.5159	3.5159
2	10	22.035	21.743	21.743	21.743	21.743	21.743
	30	22.035	22.002	22.002	22.001	22.001	22.002
	100	22.035	22.032	22.032	22.032	22.032	22.032
3	10	61.677	59.801	59.801	59.801	59.801	59.816
	30	61.677	61.478	61.478	61.478	61.478	61.493
	100	61.677	61.677	61.677	61.677	61.677	61.693

and

Table 4. The first three non-dimensional frequencies of a pure aluminum beam for a range of slenderness ratios and clamped–pinned boundary conditions.

Frequency No. i	L/h	$\mathbf{Non-Dimensional}\mathbf{Fundamental}\mathbf{Natural}\mathbf{Frequency}\left({\mathit{\lambda }}_{\mathit{i}}\right)$
		Classical Results [27]	DSM Ref [12]	DSM Code	DFE (1 ELE)	DFE (3 ELE)	FEM (10 ELE)
1	10	15.418	15.345	15.345	15.345	15.345	15.345
	30	15.418	15.410	15.410	15.410	15.410	15.410
	100	15.418	15.418	15.417	15.417	15.417	15.418
2	10	49.965	49.095	49.095	49.095	49.095	49.103
	30	49.965	49.866	49.866	49.866	49.866	48.875
	100	49.965	49.956	49.956	49.956	49.956	49.965
3	10	104.25	100.39	100.39	100.39	100.39	100.46
	30	104.25	103.80	103.80	103.80	103.80	103.87
	100	104.25	104.21	104.21	104.21	104.21	104.28

).

The first three natural frequencies for an aluminum beam using classical Euler–Bernoulli beam solution [37], DSM [21], in-house DSM and FEM codes, and DFE for various boundary conditions are presented in Table 1, Table 2, Table 3 and Table 4. It was observed that, as the slenderness ratio, $L/h$, is increased, the results of the element-based methods would converge towards the classical solution [37]. There is a slight deviation between some results of the element-based methods with the classical solution, particularly at higher frequencies.

The first (fundamental) natural frequency values obtained using DFE, DSM, FEM and reference results for Euler–Bernoulli and Timoshenko beam theories are presented in

Table 5. The first (fundamental) non-dimensional frequency of a FGM beam with $k=0.3$ under simply supported boundary conditions.

L/h	$$\mathbf{Non-Dimensional}\mathbf{Fundamental}\mathbf{Natural}\mathbf{Frequency}\left({\mathit{\lambda }}_{\mathit{i}}\right)$$
	Ref [25]		DSM Ref [13]	DSM (1 ELE)	DFE (1 ELE)	DFE (3 ELE)	FEM (3 ELE)
	Euler-Bernoulli	Timoshenko
10	17.329	17.138	17.614	17.328	17.378	17.350	17.365
30	17.392	17.373	17.676	17.395	17.447	17.420	17.433
100	17.405	17.398	17.684	17.402	17.455	17.430	17.441

. Excellent agreement was found for reference Euler–Bernoulli results and the three-element DFE solution, and deviation was found to be very small, i.e., 0.121%, 0.160%, and 0.145% for $L/h$ratios of 10, 30, and 100, respectively. For the same slenderness ratios, DSM method performed better, with a difference of 0.004%, 0.0164% and 0.0163%. Deviation between the DFE results and the Timoshenko beam values are 1.24%, 0.269%, and 0.181%. Again, DSM outperforms other methods, with a difference of 1.107%, 0.126%, and 0.0202%. The FEM showed the greatest discrepancy between results and reference values. When using the same numbers of elements, errors of 0.210%, 0.235% and 0.208% for Euler–Bernoulli and 1.323%, 0.345%, and 0.208% for Timoshenko beam theories were found.

In Figure 4, the convergence performance of the DFE and FEM for k = 1 and cantilevered (clamped–free) boundary condition is presented. The DFE formulation error decreases more at fewer elements when compared to FEM.

The first 4 non-dimensional frequencies for various FGM beams (k = 0.1, k = 1, and k = 5) and different boundary conditions are presented in

Table 6. The first 4 non-dimensional frequencies for a FGM beam for (k = 0.1, 1, and 5), under clamped–free boundary conditions.

Frequency No. i	$\mathbf{Non-Dimensional}\mathbf{Fundamental}\mathbf{Natural}\mathbf{Frequency}\left({\mathit{\lambda }}_{\mathit{i}}\right)$
	k = 0.1			k = 1			k = 5
	DSM Code (1 ELE)	DFE (10 ELE)	FEM (10 ELE)	DSM Code (1 ELE)	DFE (10 ELE)	FEM (10 ELE)	DSM Code (1 ELE)	DFE (10 ELE)	FEM (10 ELE)
1	6.2673	6.2673	6.2673	4.7669	4.7672	4.7672	4.0493	4.0495	4.0495
2	39.272	39.272	39.274	29.870	29.883	29.883	25.374	25.381	25.382
3	109.92	109.95	109.97	83.599	83.703	83.722	71.014	71.085	71.102
4	215.38	215.40	215.60	163.81	164.11	164.26	139.16	139.34	139.47

,

Table 7. The first 4 non-dimensional frequencies for a FGM beam for (k = 0.1, 1, 5), under simply supported boundary conditions.

Frequency No. i	$\mathbf{Non-Dimensional}\mathbf{Fundamental}\mathbf{Natural}\mathbf{Frequency}\left({\mathit{\lambda }}_{\mathit{i}}\right)$
	k = 0.1			k = 1			k = 5
	DSM Code (1 ELE)	DFE (10 ELE)	FEM (10 ELE)	DSM Code (1 ELE)	DFE (10 ELE)	FEM (10 ELE)	DSM Code (1 ELE)	DFE (10 ELE)	FEM (10 ELE)
1	17.592	17.592	17.592	13.381	13.382	13.382	11.366	11.367	11.367
2	70.360	70.361	70.367	53.515	53.542	53.547	45.459	45.476	45.481
3	158.28	158.28	158.37	120.38	120.51	120.58	102.26	102.34	102.40
4	281.30	281.32	281.78	213.94	214.36	214.69	181.74	182.00	182.29

,

Table 8. The first 4 non-dimensional frequencies for a FGM beam for (k = 0.1, 1, and 5), under clamped–clamped boundary conditions.

Frequency No. i	$\mathbf{Non-Dimensional}\mathbf{Fundamental}\mathbf{Natural}\mathbf{Frequency}\left({\mathit{\lambda }}_{\mathit{i}}\right)$
	k = 0.1			k = 1			k = 5
	DSM Code (1 ELE)	DFE (10 ELE)	FEM (10 ELE)	DSM Code (1 ELE)	DFE (10 ELE)	FEM (10 ELE)	DSM Code (1 ELE)	DFE (10 ELE)	FEM (10 ELE)
1	39.872	39.872	39.881	30.326	30.333	30.352	25.761	25.766	25.778
2	109.89	109.88	109.95	83.576	83.627	83.734	70.997	71.027	71.103
3	215.35	215.36	215.66	163.79	164.00	164.38	139.13	139.27	139.56
4	355.83	355.85	356.99	270.63	271.21	272.40	229.90	230.27	231.19

and

Table 9. The first 4 non-dimensional frequencies for a FGM beam for (k = 0.1, 1, and 5), under clamped–pinned boundary conditions.

Frequency No. i	$\mathbf{Non-Dimensional}\mathbf{Fundamental}\mathbf{Natural}\mathbf{Frequency}\left({\mathit{\lambda }}_{\mathit{i}}\right)$
	k = 0.1			k = 1			k = 5
	DSM Code (1 ELE)	DFE (10 ELE)	FEM (10 ELE)	DSM Code (1 ELE)	DFE (10 ELE)	FEM (10 ELE)	DSM Code (1 ELE)	DFE (10 ELE)	FEM (10 ELE)
1	27.485	27.485	27.485	20.947	20.954	20.955	17.783	17.788	17.788
2	89.050	89.053	89.069	67.770	67.832	67.842	57.560	57.597	57.607
3	185.75	185.77	185.90	141.32	141.56	141.66	120.04	120.19	120.27
4	317.55	317.59	318.25	241.55	242.21	242.70	205.20	205.60	206.02

. For all boundary conditions, the results from DFE and FEM showed greater deviation from the baseline DSM results at higher frequency modes. In regard to material grading, the highest deviation between DFE and FEM from DSM results occurred in the linear variation (k = 1) case. For all cases, DFE performed as well or better than FEM, particularly at higher modes. The mode shapes for k = 1 and each boundary condition are presented in Figure 5. The mode shapes for varied material grading (k = 0.1 and k = 5) are presented in Figure 6.

4. Discussion

Excellent agreement is achieved for DFE in all conditions. The classical results used generic closed-form solutions to the 4th-order Euler–Bernoulli governing equation solution that is often found in standard texts (see, e.g., [37]). Due to the highly convergent nature of DSM and DFE formulations, their results converge at higher frequencies, while FEM results deviate. For FEM, additional elements would be required to match the accuracy of DFE and DSM.

Comparison between the data in Figure 3 and Figure 4 and Table 1, Table 2, Table 3 and Table 4 further illustrates the superior convergence performance of DFE over the conventional FEM. The 3rd mode of an aluminum beam converges at around 10–12 elements using FEM. In contrast, DFE results show that a single element is capable of producing the exact solution. This is due to the fact that the DFE method uses the solution of the uncoupled governing equations to create the DTSFs. Since this validation case uses no material grading, the coupling terms approach zero and the DFE solution is exact within the limits of the theory, i.e., attaining the same results as the DSM. When material grading is introduced, a single DFE element does not produce the exact results like the DSM method, i.e., due to coupling effects. Additional DFE elements are required to converge to the solution. However, as seen in Figure 4, convergence occurs with fewer DFE elements than FEM.

Note that the in-house DSM code and the reference DSM results deviate by about 1.6%. The authors were unable to resolve the difference between the two results; however, there is an excellent agreement between the results of the in-house DSM code and reference [25]. Furthermore, the results of the in-house DSM, DFE and FEM converge to the same solution at higher element numbers. Therefore, for validation and comparison purposes, the in-house DSM results were used. Since the formulations in this paper are based on Euler–Bernoulli assumptions, a greater discrepancy is expected when comparing to reference Timoshenko beam results.

The results of Table 6, Table 7, Table 8 and Table 9 show that, when using the same number of elements, the DFE model outperforms the FEM. This becomes more apparent at higher frequencies. For example, in the case of a cantilever beam with $k=1$, the 4th frequency shows a discrepancy of 0.274% and 0.183%, for FEM and DFE, respectively, when compared with DSM. Now referring to the first frequency results, the error is found to be 0.00629% for both FEM and DFE when compared with DSM. To accurately model the 4th mode using FEM, one would need at least 4 elements to express the mode shape. Since DFE uses frequency-dependent interpolation functions, a single element can produce an infinite number of frequencies. Figure 4 shows that even with a single element, the DFE produces less errors than 5 FEM elements. At higher frequencies, fewer DFE elements can be used to develop results that are better or comparable to FEM. For the 4th mode, DFE and FEM begin to converge at 10 elements.

The mode shapes for a FGM beam of $k=1$, and clamped–free, simply supported, clamped–clamped, and clamped–pinned boundary conditions, are shown in Figure 5. Bending was found to be the dominant mode and therefore the axial mode shapes are not shown. As can be observed in Figure 6, when varying $k$, the dominant mode shape remained unchanged. Therefore, for the current theory, the through thickness variance of material properties does not affect the mode shape. This is expected due to the fact that the Euler–Bernoulli bending theory assumes that the cross section remains a plane, and the beam deforms only on the neutral axis.

5. Conclusions

A Dynamic Finite Element (DFE) formulation for the free vibration analysis of Functionally Graded Beams was presented. Analysis on the natural frequency and mode shapes of functionally graded Euler–Bernoulli beams, using DFE and other methods was performed. Results from the presented DFE element were compared against FEM, exact results of DSM, and other data found in the open literature. Boundary conditions, slenderness ($L/h$) ratio, and $k$ values were varied to investigate their effects on the beam vibration response. The convergence of the DFE formulation outperformed conventional FEM, particularly at higher frequencies. Since the formulation is based on Euler–Bernoulli assumption, the mode shapes were not affected by the through the thickness variation in properties. In the absence of material grading, the coupling effects, and the corresponding terms in the element’s DFE matrix, reduces to zero and a one element DFE model produces exact results within the limits of the theory.

This study was limited to Euler–Bernoulli theory and in future studies the DFE model could be extended to higher-order beam theories such as Timoshenko beams. Also, this study was limited to solving a simple beam element, whereas a future study could demonstrate the generality of DFE in solving structures such as frames.