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

: A Dynamic Finite Element (DFE) method for coupled axial–ﬂexural 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, veriﬁcation, and validation. Good agreement was observed and the DFE formulation exhibited superior convergence performance compared to the FEM.


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 reentry [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 freevibration 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.

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 ρ. Figure 1. Coordinate system and notation for formulation [22].
The effective material property, P, for a FGM is evaluated using the following rule of mixtures, expressed as [9,21,22]: (1) where Pt and Pb 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. Vt is the volume fraction of the top constituent of the beam, defined as [9,21,22]: 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]: Figure 1. Coordinate system and notation for formulation [22].
The effective material property, P, for a FGM is evaluated using the following rule of mixtures, expressed as [9,21,22]: 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]: 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. The effective material property, P, for a FGM is evaluated using the following rule of mixtures, expressed as [9,21,22]: (1) where Pt and Pb 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. Vt is the volume fraction of the top constituent of the beam, defined as [9,21,22]: 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]: 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]: 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: Assuming simple harmonic motion, axial and lateral displacements, v(y, t) and w(y, t), respectively, both functions of space (y) and time (t), can be written as (i.e., separation of variables): v(y, t) = V(y)e iωt w(y, t) = W(y)e iωt (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: 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ω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: 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: A set of integration by parts is performed, leading to the following weak form of the above integral equations: The above expressions also satisfy the principle of virtual work, 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: 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: The above integral equations W k v and W k w , 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: Applying additional integration by parts to the above discretized (elemental) weakform equations, leads to the following (equivalent) form of the equations: 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: where C 1,2,3,4 and D 1,2 are constants and The interpolation functions, also referred to as shape functions of approximation space, are then obtained as follows. First, using the generalized parameters a , δa , b , and δb 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: The axial and flexural dynamic basis functions of approximation space, respectively, also reported in [28], are defined as: Replacing the generalized parameters a , δa , b , and δb with the nodal variables, and δW 1 δW 1 δW 2 δW 2 , respectively, expression (26) can be rewritten as: (29) where the matrices [P n ] v and [P n ] w are defined as: By combining (26), (30), and (31) the approximation in terms of nodal variables can be written as: where N(ξ) v and N(ξ) w are the frequency-dependent trigonometric shape functions for the axial and flexure, respectively. Expression (31) can be combined and rewritten as: 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]: where The following shape functions are used to state the approximate elemental bending displacement in terms of nodal displacements, along the domain [28,30]: Using expressions (20) and (21) and the shape functions (35) through (42), the element dynamic stiffness matrix, [K(ω)] k , is derived. This consists of uncoupled and coupled dynamic stiffness matrices, [K(ω)] k u and [K(ω)] k c , respectively. There are four coupled and four uncoupled matrix components. The element dynamic stiffness matrix, [K(ω)] k , is formulated by adding the eight coupled and uncoupled submatrices.
The system's global dynamic stiffness matrix, [K(ω)], is then obtained by assembling all the element matrices and applying the boundary conditions. The nonlinear eigenvalue problem, resulting from this method, is: where {U n } is the vector of global nodal displacements of the system. The natural frequencies of the system would be the values of ω, which yields a zero determinant, |K(ω)| = 0, for the global dynamic stiffness matrix [28].
It is worth noting that the basis functions are specifically chosen so that when the frequency, ω, 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 = ∞ (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: ρ t = 2700.0027 kg/m 3 , ρ b = 2700 kg/m 3 , E t = 70.00007 GPa, and E b = 70 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: 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 (Al 2 O 3 ), as also presented in Ref. [9]. The properties of the aluminum are: E b = 70 GPa, ρ b = 2700 kg/m 3 and the properties of the alumina are: E t = 380 GPa, ρ t = 3800 kg/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]: with I o and A o defined previously in (5). The natural frequencies, ω, 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 = 210GPa, ρ b = 7800 kg/m 3 and the properties of the alumina are: E t = 390 GPa, rho t = 3960 kg/m 3 , as also reported in [9]. A comparison is made between the convergence performance of DFE and FEM formulations.

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 (Tables 1-4).    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 Tables 1-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 4 non-dimensional frequencies for various FGM beams (k = 0.1, k = 1, and k = 5) and different boundary conditions are presented in Tables 6-9. 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.
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. 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 Tables 6-9. 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. 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

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 Figures 3 and 4 and Tables 1-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 Tables 6-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 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 frequencydependent 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.

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 Figures 3 and 4 and Tables 1-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 inhouse 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 Tables 6-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.

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.