Closed-Form Exact Solution for Free Vibration Analysis of Symmetric Functionally Graded Beams

Abstract

The dynamic stiffness method is developed to analyze the natural vibration characteristics of functionally graded beams, where material properties change continuously across the beam thickness following a symmetric law distribution. The governing equations of motion and associated natural boundary conditions for free vibration analysis are derived using Hamilton’s principle and closed-form exact solutions are obtained for both Euler–Bernoulli and Timoshenko beam models. The dynamic stiffness matrix, which governs the relationship between force and displacements at the beam ends, is determined. Using the Wittrick–Williams algorithm, the dynamic stiffness matrix is employed to compute natural frequencies and mode shapes. The proposed procedure is validated by comparing the obtained frequencies with those given by approximated well-known formulas. Finally, a parametric investigation is conducted by varying the geometry of the structure and the characteristic mechanical parameters of the functionally graded material.

1. Introduction

Functionally graded materials (FGMs) are deliberately engineered to exhibit continuously varying properties in one or more directions, offering designers extensive versatility in distributing strength and stiffness as needed. Consequently, FGMs have found successful applications in various scientific and engineering fields, including the design of aircraft and spacecraft structures [1]. The study of the static and dynamic behavior of such materials, especially from a structural design perspective, is therefore essential. Beams, acting as load-bearing components, are principal candidates for fabrication from FGMs. Beams made of FGMs can be designed to have specific vibrational characteristics, thus enhancing the stability of structures. FGMs are also beneficial in terms of thermal resistance since the gradual transition in material properties can reduce thermal stresses and prevent thermal fatigue; experimental tests on FGM concrete beams have been performed in [2]. FGMs can also be used to model circular concrete columns confined by Fiber-Reinforced Polymers (FRP) in order to predict their collapse [3].

Consequently, there has been a significant focus on the free [4,5,6,7,8,9,10,11] and forced [9,10,12,13] vibration analysis of functionally graded beams (FGBs) in recent years. The solutions have been proposed for both Euler–Bernoulli [4,6,7,8,10,13] and Timoshenko [5,7,9,10] beam models. In these studies, material properties are usually assumed to vary continuously with a power law [4,5,6,7,12,13], an exponential law [9,10,11] or an arbitrary law [8] along one or more directions. A considerable portion of the study of free vibration analysis of FGBs has relied on finite element and other approximate methods [14,15,16]. Although these methods represent valuable advancements, their results are often contingent upon the number and quality of elements employed in the analysis, leading to potential unreliability, particularly at higher frequencies.

In addition to the predominantly employed numerical methods, further studies used analytical approaches to address the free vibration problem of FGMs. Some of these are based on the use of the dynamic stiffness method (DSM) [4,5,9]. Unlike many numerical methods, the DSM adopts accurate member theory, incorporating frequency-dependent shape functions derived from the solution of the governing differential equations of motion for the structural element undergoing free vibration. Consequently, the DSM yields exact results for all natural frequencies and mode shapes without resorting to any approximation. This independence from the number of elements used in the analysis makes the DSM notably attractive and accurate compared to finite element and other approximate methods.

In FGM, different laws of variation in the mechanical characteristics in the transversal and/or axial direction can be assumed. For example, Banerjee et al. [4,5] applied the DSM to functionally graded materials in which the material properties vary continuously through the beam thickness direction according to a power law distribution. To the author’s knowledge, the case of a symmetric variation in mechanical characteristics with respect to the centroid has never been studied. For this kind of variation, the results reported in [4,5] cannot be applied since some terms of the dynamic stiffness matrix would take an indeterminate form or tend towards infinity.

Symmetric FGMs could find an interesting technical application in modeling the retrofitting of reinforced concrete framed structures, where the reinforcement of the existing columns is realized by a symmetric section augmentation with high-performance concrete; this technique has several examples of application (see, for example, [17]).

With the aim of studying the free vibrations of beam structures in which the material properties are assumed to vary continuously along the beam thickness according to a symmetric distribution, in this paper, the dynamic stiffness matrix is determined, starting from the differential equations of motion derived from Hamilton’s principle.

Subsequently, the dynamic stiffness matrix is utilized in conjunction with the Wittrick–Williams algorithm [18] to compute the natural frequencies and mode shapes of some illustrative examples. In order to apply the Wittrick and Williams algorithm in conjunction with the dynamic stiffness matrix of symmetric functionally graded materials, the term J₀, representing the frequencies of vibration of clamped–clamped beams, has to be determined. For homogeneous materials, this term has a well-known expression provided in the scientific literature. In the case of FGM, this expression is not available, so an original derivation of J₀ is proposed here. The frequencies of vibration of selected structures composed of symmetric functionally graded materials have been validated through the comparison with the results obtained, evaluating the corresponding Rayleigh quotient. Finally, a parametric investigation is conducted, varying the geometry of the structure and the characteristic mechanical parameters of the functionally graded material. The numerical applications refer to the mechanical parameters of the outer material being either greater or smaller than those of the inner one, and the obtained results could be used, for example, in simulating the effects of retrofitting or material degradation, respectively.

2. The Considered Functionally Graded Beam

The elementary beam model considered in this study is a rectangular section of length L, height h and width b, with y as the beam axis and cross section in the x–z plane (Figure 1).

The beam is characterized by Young’s modulus E and density ρ varying through the height of the section in the z direction with a symmetric law with respect to the centroid. The variation could be parabolic, for example, as reported below:

$$E(\zeta )=E_2+(E_1-E_2){\zeta }^2, \rho (\zeta )={\rho }_2+({\rho }_1-{\rho }_2){\zeta }^2$$

(1)

where E₂ and ρ₂ are the properties of the beam at the center line of the section, E₁ and ρ₁ are the properties at the top and bottom surfaces of the beam, and ζ = z/(h/2) is the dimensionless abscissa along the vertical axis represented in Figure 1. In Figure 2, a clearer representation of Equation (1) is reported, where the variation in Young’s modulus E is shown as an example. Furthermore, the material properties are assumed to be constant along the horizontal x and y axes.

2.1. The Governing Differential Equations of Motion

In this subsection, the equations of motion in free vibrations of a beam made by FGM with symmetric variation in the mechanical properties and the associated closed-form solution are obtained. The formulations for the Euler–Bernoulli and Timoshenko models are briefly reported in the following subsubsections.

2.1.1. Euler–Bernoulli Beam Model

The displacement components u₁, v₁ and w₁, respectively, along the x, y and z axes, which characterize the Euler–Bernoulli beam model, can be assumed as

$$u_1=0, v_1(y,z,t)=v(y,t)-z\varphi (y,t), w_1(y,z,t)=w(y,t)$$

(2)

where $\varphi (y,t)={\displaystyle \frac{\partial w(y,t)}{\partial y}}$ is the flexural rotation in the y-z plane.

The potential and kinetic energies U_P and U_K of the FGB are, after some simplifications, given by [5]

$$\begin{array}{l}{U}_p={\displaystyle \frac{1}{2}}{\displaystyle \underset{0}{\overset{L}{\int }}(A_0{v^{\prime }}^2-2A_1{v}^{\prime }{w}^{″}+A_2{w^{″}}^2)dy}\\ U_K={\displaystyle \frac{1}{2}}{\displaystyle \underset{0}{\overset{L}{\int }}\left[B_0({\dot{v}}^2+{\dot{w}}^2)-2B_1\dot{v}{\dot{w}}^{\prime }+B_2{{\dot{w}}^{\prime }}^2\right]}dy\end{array}$$

(3)

where prime and over-dot denote differentiation with respect to space y and time t, respectively, and the parameters A and B, which consider the variation in material properties, are

$$A_i={\displaystyle \underset{A}{\int }{z}^{i}E(z)}dA, B_i={\displaystyle \underset{A}{\int }{z}^i\rho (z)}dA\hspace{1em}i=0,1,2$$

(4)

By applying Hamilton’s principle to the displacement field, it is possible to obtain the differential equations of motion in free vibrations. In particular, it is assumed that v(y,t) and w(y,t) can be expressed in harmonic form:

$$v(y,t)=V(y)e^{i\omega t}, w(y,t)=W(y)e^{i\omega t}$$

(5)

where V(y) and W(y) are the mode shapes and ω is the natural frequency. Introducing the dimensionless abscissa ξ = y/L and the differential operator D = d/dξ, and considering that, according to the symmetric variation laws Equation (1), the terms with i = 1 in Equation (4) become zero, the equations of motion take the following form:

$$\begin{array}{l}(B_0{\omega }^2{L}^3+A_{0}L{D}^2)V(\xi )=0\\ (B_0{\omega }^2{L}^4-B_2{\omega }^2{L}^2{D}^2-A_2{D}^4)W(\xi )=0\end{array}$$

(6)

As it can be noticed, axial and bending contributions are decoupled and each displacement component can be easily determined from the individual equations.

The second equation in Equation (6), assuming $W(\xi )=e^{\lambda \xi }$, can be written in the form

$$(D^4+a{D}^2+b)e^{\lambda \xi }=0$$

(7)

where

$$a={\displaystyle \frac{B_2}{A_2}}{L}^2{\omega }^2, b=-{\displaystyle \frac{B_0}{A_2}}{L}^4{\omega }^2$$

(8)

Equation (7) can be simply solved, leading to the following expression for the transversal displacement component:

$$W(\xi )=Q_1{e}^{{\lambda }_1\xi }+Q_2{e}^{{\lambda }_2\xi }+Q_3{e}^{{\lambda }_3\xi }+Q_4{e}^{{\lambda }_4\xi }$$

(9)

with

$${\lambda }_{1,2,3,4}=\pm \sqrt{{\displaystyle \frac{-a\pm \sqrt{a^2-4b}}{2}}}$$

(10)

where $Q_j,\hspace{0.17em}\hspace{0.17em}j=1,\dots ,4$ are constants to be obtained from the boundary conditions.

With a similar procedure, it is possible to evaluate the axial displacements in the form

$$V(\xi )=P_1{e}^{{\eta }_1\xi }+P_2{e}^{{\eta }_2\xi }$$

(11)

where ${\eta }_{1,2}=\pm i\sqrt{c}$; $c={\displaystyle \frac{B_0}{A_0}}{L}^2{\omega }^2$ and $P_j,j=1,2$ are constants to be obtained from the boundary conditions.

2.1.2. Timoshenko Beam Model

For the Timoshenko beam, the displacement components can still be expressed in the form of Equation (2) but, while for the Euler–Bernoulli beam the rotation ϕ(y,t) is equal to the derivative of the transversal displacement, for the Timoshenko model, ϕ(y,t) is an independent variable related to the total rotation of the cross section ${\displaystyle \frac{\partial w\left(y,t\right)}{\partial y}}$ and the shear strain ψ(y,t) as

$$\psi \left(y,t\right)={\displaystyle \frac{\partial w\left(y,t\right)}{\partial y}}-\varphi \left(y,t\right)$$

(12)

The application of Hamilton’s principle to the Timoshenko beam therefore leads to a system of three differential equations of motion.

Assuming that ϕ(y,t) can also be expressed in harmonic form with amplitude Ф(y) and setting

$$A_3={\displaystyle \underset{A}{\int }G\left(z\right)}dA$$

(13)

where G(z) is the shear modulus of the beam varying through the height of the section according to a symmetric law formally identical to the ones reported in Equation (1), the differential equations of motion with respect to the dimensionless abscissa ξ, considering the property of symmetry in the variation in the mechanical parameters of the material, take the form

$$\begin{array}{l}(B_0{\omega }^2{L}^2+A_0{D}^2)V(\xi )=0\\ B_0{\omega }^2{L}^{2}W(\xi )+A_3{D}^{2}W(\xi )-A_{3}LD\Phi (\xi )=0\\ A_{3}LDW(\xi )+\left(B_2{\omega }^2-A_3\right)L^2\Phi (\xi )+A_2{D}^2\Phi (\xi )=0\end{array}$$

(14)

As it can be noticed, also for the Timoshenko model, the axial and bending problems are decoupled. Assuming $W(\xi )=e^{\lambda \xi }$ and $\Phi (\xi )=e^{\lambda \xi }$, obtaining Ф(ξ) from the second equation in Equation (14) and substituting in the third equation, the differential equation governing the transversal displacement can be written in the form

$$\left(D^4+a{D}^2+b\right)e^{\lambda \xi }=0$$

(15)

where

$$a={\displaystyle \frac{A_3{B}_2+A_2{B}_0}{A_2{A}_3}}{\omega }^2{L}^2, b={\displaystyle \frac{B_0\left(B_2{\omega }^2-A_3\right)}{A_2{A}_3}}{\omega }^2{L}^4$$

(16)

The same equation for Ф(ξ) would have been obtained if W(ξ) had been isolated in the third equation of Equation (14) and substituted in the second equation. Therefore, solutions for transversal displacement and flexural rotation are

$$\begin{array}{l}W(\xi )=Q_1{e}^{{\lambda }_1\xi }+Q_2{e}^{{\lambda }_2\xi }+Q_3{e}^{{\lambda }_3\xi }+Q_4{e}^{{\lambda }_4\xi }\\ \Phi (\xi )=R_1{e}^{{\lambda }_1\xi }+R_2{e}^{{\lambda }_2\xi }+R_3{e}^{{\lambda }_3\xi }+R_4{e}^{{\lambda }_4\xi }\end{array}$$

(17)

where λ_i has the same formal expression derived for the Euler–Bernoulli beam in Equation (10).

The constants Q_j and R_j, j = 1, …, 4 in Equation (17) can be related to each other by substituting the expressions for W(ξ) and Ф(ξ) in the second part of Equation (14):

$$Q_j={\beta }_j{R}_j$$

(18)

where

$${\beta }_j={\displaystyle \frac{A_{3}L{\lambda }_j}{B_0{\omega }^2{L}^2+A_3{\lambda }_j^2}}$$

(19)

The axial problem can be solved similarly to Euler–Bernoulli beam model.

2.2. The Dynamic Stiffness Matrix

In this subsection, the derivation of the dynamic stiffness matrix of a beam made of symmetric FGM is reported for the Euler–Bernoulli and Timoshenko beams models. To this aim, the nodal displacements W, V, and Ф and the nodal forces F, M, and S are evaluated according to the convention reported in Figure 3 and applying the boundary conditions at ξ = 0 and ξ = 1 reported in Figure 4.

2.2.1. Euler–Bernoulli Beam Model

From the expressions of the transversal and axial displacements in Equations (9) and (11), it is possible to obtain the following expressions of the flexural rotation Ф(ξ), axial force F(ξ), bending moment M(ξ) and shear force S(ξ):

$$\Phi (\xi )={\displaystyle \frac{W^{\prime }}{L}}={\displaystyle \frac{1}{L}}\left(Q_1{\lambda }_1{e}^{{\lambda }_1\xi }+Q_2{\lambda }_2{e}^{{\lambda }_2\xi }+Q_3{\lambda }_3{e}^{{\lambda }_3\xi }+Q_4{\lambda }_4{e}^{{\lambda }_4\xi }\right)$$

(20)

$$F(\xi )=-{\displaystyle \frac{A_0}{L}}{V}^{\prime }=-{\displaystyle \frac{A_0}{L}}\left(P_1{\eta }_1{e}^{{\eta }_1\xi }+P_2{\eta }_2{e}^{{\eta }_2\xi }\right)$$

(21)

$$M(\xi )=-{\displaystyle \frac{A_2}{L^2}}{W}^{″}=-{\displaystyle \frac{A_2}{L^2}}\left(Q_1{\lambda }_1^2{e}^{{\lambda }_1\xi }+Q_2{\lambda }_2^2{e}^{{\lambda }_2\xi }+Q_3{\lambda }_3^2{e}^{{\lambda }_3\xi }+Q_4{\lambda }_4^2{e}^{{\lambda }_4\xi }\right)$$

(22)

$$\begin{array}{c}S(\xi )={\displaystyle \frac{A_2}{L^3}}\left(W^{‴}+{\displaystyle \frac{B_2{L}^2{\omega }^2}{A_2}}{W}^{\prime }\right)={\displaystyle \frac{A_2}{L^3}}[Q_1{\lambda }_1^3{e}^{{\lambda }_1\xi }+Q_2{\lambda }_2^3{e}^{{\lambda }_2\xi }+Q_3{\lambda }_3^3{e}^{{\lambda }_3\xi }+Q_4{\lambda }_4^3{e}^{{\lambda }_4\xi }+\\ +{\displaystyle \frac{B_2{L}^2{\omega }^2}{A_2}}(Q_1{\lambda }_1{e}^{{\lambda }_1\xi }+Q_2{\lambda }_2{e}^{{\lambda }_2\xi }+Q_3{\lambda }_3{e}^{{\lambda }_3\xi }+Q_4{\lambda }_4{e}^{{\lambda }_4\xi })]\end{array}$$

(23)

The nodal displacements and forces vectors can be defined as

$$\delta ={\left[\begin{array}{cccccc}{V}_1& W_1& {\Phi }_1& V_2& W_2& {\Phi }_2\end{array}\right]}^T, P={\left[\begin{array}{cccccc}{F}_1& S_1& M_1& F_2& S_2& M_2\end{array}\right]}^T$$

(24)

Using Equations (20)–(23), the vectors $\delta$ and $P$ can be expressed in matrix form as

$$\delta =B^{EB}R$$

(25)

$$P=A^{EB}R$$

(26)

where R = [P₁ P₂ Q₁ Q₂ Q₃ Q₄]^T. Obtaining R from Equation (25) and substituting into Equation (26), the dynamic stiffness matrix K^EB = A^EB(B^EB)⁻¹ that relates nodal displacements and forces can be obtained:

$$P=K^{EB}\delta$$

(27)

The explicit expressions of matrices A^EB and B^EB for the Euler–Bernoulli beam model are not reported here for the sake of shortness but can be found in Appendix A.

2.2.2. Timoshenko Beam Model

The shear force S(ξ) and the bending moment M(ξ) for the Timoshenko beam model assume the following expressions:

$$\begin{array}{c}S(\xi )={\displaystyle \frac{A_3}{L}}\left(-W^{\prime }+L\Phi \right)={\displaystyle \frac{A_3}{L}}[-(Q_1{\lambda }_1{e}^{{\lambda }_1\xi }+Q_2{\lambda }_2{e}^{{\lambda }_2\xi }+Q_3{\lambda }_3{e}^{{\lambda }_3\xi }+Q_4{\lambda }_4{e}^{{\lambda }_4\xi })+\\ +L(R_1{e}^{{\lambda }_1\xi }+R_2{e}^{{\lambda }_2\xi }+R_3{e}^{{\lambda }_3\xi }+R_4{e}^{{\lambda }_4\xi })]\end{array}$$

(28)

$$M(\xi )=-{\displaystyle \frac{A_2}{L}}{\Phi }^{\prime }=-{\displaystyle \frac{A_2}{L}}\left(R_1{\lambda }_1{e}^{{\lambda }_1\xi }+R_2{\lambda }_2{e}^{{\lambda }_2\xi }+R_3{\lambda }_3{e}^{{\lambda }_3\xi }+R_4{\lambda }_4{e}^{{\lambda }_4\xi }\right)$$

(29)

where the rotation Ф(ξ) is already reported in Equation (17). It has to be noted that the axial force does not depend on the beam model; therefore, the expression of F(ξ) in this case is the same as in Equation (21).

Considering the relationship between the constants Q_j = β_jR_j, it is possible to obtain nodal forces and displacements at ξ = 0 and ξ = 1 and collect them in the vectors reported in Equation (24). Analogously to the Euler–Bernoulli beam model, the dynamic stiffness matrix K^TIM = A^TIM(B^TIM)⁻¹ that relates nodal displacements and forces can be obtained:

$$P=K^{TIM}\delta$$

(30)

The explicit expressions of matrices A^TIM and B^TIM for the Timoshenko beam are reported in Appendix A.

3. Application of the Wittrick and Williams Algorithm

The exact frequencies of vibration of simple beams or framed structures may be obtained applying the Wittrick–Williams algorithm [18] in conjunction with the dynamic stiffness matrix of the considered structure. This algorithm allows the evaluation of the number J of vibration frequencies that are smaller than a trial value ω*, by means of an iterative procedure, to converge to any required accuracy. The number J is given by

$$J=J_K+J_0$$

(31)

where J_k is the number of negative eigenvalues of the dynamic stiffness matrix evaluated at the specified frequency value ω* and $J_0={\displaystyle {\sum }_{b=1}^{Nbeams}{J}_{0,b}}$ is the number of frequencies of vibration of the beams considered with both ends clamped which are lower than ω*.

The evaluation of J_k in the case of beams composed of symmetric functionally graded material can be obtained once the dynamic stiffness matrix of the structure is evaluated at the frequency value ω*. In order to compute J, the expression of J₀ for FGM has to be evaluated, and since this is not available in the scientific literature, it will be originally derived in the following.

The procedure is based on the consideration of a simply supported beam made of FGM, for which the natural frequencies of vibration, for the cases of both the Euler–Bernoulli and Timoshenko models, are given in [19]. In particular, for Timoshenko beams, the frequencies of vibration take the following expression:

$${\omega }_n^2={\displaystyle \frac{2{\widehat{E}}_2{\alpha }_n^4}{{\rho }_0}}{\displaystyle \frac{1}{{\zeta }_n+\sqrt{{\zeta }_n^2-\frac{4{\widehat{\rho }}_2{\widehat{E}}_2}{{\rho }_0{G}_0^s}{\alpha }_n^4}}}$$

(32)

where ${\widehat{E}}_2$ and ${\widehat{\rho }}_2$ are given by

$${\widehat{E}}_2=E_2-{\displaystyle \frac{E_1^2}{E_0}}, {\widehat{\rho }}_2={\rho }_2-{\displaystyle \frac{{\rho }_1^2}{{\rho }_0}}$$

(33)

and $E_i=A_i,\hspace{0.33em}{\rho }_i=B_i,\hspace{0.33em}i=0,1,2,\hspace{0.33em}{G}_0^s=A_3$ are terms that consider the variation of the material as in Equations (4)–(13) and α_n is given by

$${\alpha }_n={\displaystyle \frac{n\pi }{L}}$$

(34)

where L is the length of the beam and n is the frequency number. ζ_n has the following expression:

$${\zeta }_n=1+{\alpha }_n^2\left({\displaystyle \frac{{\widehat{\rho }}_2}{{\rho }_0}}+{\displaystyle \frac{{\widehat{E}}_2}{G_0^s}}\right)$$

(35)

In the case of Euler–Bernoulli, Equation (32) is simplified in

$${\omega }_n={\alpha }_n^2\sqrt{{\displaystyle \frac{{\widehat{E}}_2}{{\rho }_0+{\widehat{\rho }}_2{\alpha }_n^2}}}$$

(36)

Obtaining α_n from Equation (32) or Equation (36), depending on the beam model, it is possible to isolate n and, by setting ω_n = ω*, find the number of natural frequencies below the assigned frequency, thus obtaining J_b for a simply supported beam made of FGM.

With respect to the determination of J_k,b, the DSM for the simply supported beam can be obtained, according to the boundary conditions, by imposing M₁ = M₂ = v₁ = v₂ = 0; therefore, the matrix system is as follows:

$$\left[\begin{array}{c}0\\ 0\end{array}\right]=\left[\begin{array}{cc}{K}_{2,2}& K_{2,4}\\ K_{4,2}& K_{4,4}\end{array}\right]\left[\begin{array}{c}{\phi }_1\\ {\phi }_2\end{array}\right]$$

(37)

The term J_k,b is the number of negative eigenvalues of the DSM in Equation (37) evaluated at the specified frequency value ω*.

Finally, the term, J_0,b (which, as already said, refers to the number of frequencies of vibration of the beams considered with both ends clamped which are lower than ω*) can be evaluated by applying Equation (31) to a simply supported beam composed of FGM and calculating J_0,b = J_b − J_k,b.

4. Results and Discussion

In this section, with the aim of investigating the influence of material variability on the eigenproperties of structures, the proposed procedure has been applied both to single beams with different boundary conditions and to framed structures. Firstly, in order to validate the obtained results, the frequencies of vibration of some single beams composed of symmetric functionally graded materials have been compared with those obtained by evaluating the corresponding Rayleigh quotient, showing a very good correspondence. Furthermore, some results reported in the literature with reference to homogeneous material have been re-obtained by means of the dynamic stiffness method. Finally, some parametric analyses by varying the E₁/E₂ ratio on beams and on a simple portal frame have been conducted.

4.1. Validation of the Proposed Procedure

4.1.1. Rayleigh’s Quotient

In this subsection, the fundamental natural frequency of vibration obtained with the proposed procedure is compared to the one obtained through Rayleigh’s quotient, defined as

$${\omega }^2={\displaystyle \frac{{\displaystyle {\int }_{L}k{\left({\varphi }^{″}\right)}^2\hspace{0.33em}dy}}{{\displaystyle {\int }_{L}m{\varphi }^2\hspace{0.33em}dy}}}$$

(38)

where $k$ and m are the distributed stiffness and mass of the beam, respectively, ω² represents the frequency of vibration, and ϕ is the mode shape of the beam.

In Equation (38), the numerator is a quantity proportional to the potential energy, while the denominator is a measure of the kinetic energy. Thus, this quotient may also be derived by equating the maximum value of kinetic energy to the maximum value of potential energy. If potential and kinetic energies are calculated by imposing the exact mode shape of the beam under particular boundary conditions, Rayleigh’s quotient leads to the exact value of the fundamental natural frequency of the beam.

The kinetic and potential energies for an Euler–Bernoulli beam composed of symmetric FGM have been reported in Equation (3). Assuming a harmonic variability with time for the transversal displacement,

$$w\left(y,t\right)=\psi \left(y\right)\cdot cos\left({\omega }_{n}t\right)$$

(39)

The maximum values of U_k and U_p assume the form

$$\left.\begin{array}{l}{U}_{K,\mathit{max}}={\displaystyle \frac{1}{2}}{\displaystyle \underset{0}{\overset{L}{\int }}\left(B_0{\psi }^2\left(y\right){{\omega }_n}^2+B_2{{\psi }^{\prime }}^2\left(y\right){{\omega }_n}^2\right)dy}\\ U_{p,\mathit{max}}={\displaystyle \frac{1}{2}}{\displaystyle \underset{0}{\overset{L}{\int }}{A}_2{{\psi }^{″}}^2\left(y\right)dy}\end{array}\right\}\Rightarrow U_{K,\mathit{max}}=U_{p,\mathit{max}}$$

(40)

Rayleigh’s quotient can be therefore written as

$${{\omega }_n}^2={\displaystyle \frac{{\displaystyle \underset{0}{\overset{L}{\int }}{A}_2{{\psi }^{″}}^2\left(y\right)dy}}{{\displaystyle \underset{0}{\overset{L}{\int }}\left(B_0{\psi }^2\left(y\right)+B_2{{\psi }^{\prime }}^2\left(y\right)\right)dy}}}$$

(41)

With reference, for example, to an Euler–Bernoulli simply supported beam, the Rayleigh quotient can be exactly evaluated, assuming the shape of the vibration mode as

$$\psi \left(y\right)=Asin\left({\displaystyle \frac{\pi y}{L}}\right)$$

(42)

and therefore

$${{\omega }_n}^2={\displaystyle \frac{A_2{\displaystyle \underset{0}{\overset{L}{\int }}\left[\frac{{\pi }^4}{L^4}{A}^{2}si{n}^2\left(\frac{\pi y}{L}\right)\right]dy}}{B_0{\displaystyle \underset{0}{\overset{L}{\int }}\left[A^{2}si{n}^2\left(\frac{\pi y}{L}\right)\right]dy}+B_2{\displaystyle \underset{0}{\overset{L}{\int }}\left[\frac{{\pi }^2}{L^2}{A}^{2}si{n}^2\left(\frac{\pi y}{L}\right)\right]dy}}}={\displaystyle \frac{A_2\frac{{\pi }^4}{L^4}}{B_0+B_2\frac{{\pi }^2}{L^2}}}$$

(43)

The comparison of the frequencies of vibration obtained by means of the dynamic stiffness matrix and the Rayleigh quotient of simply supported beams is reported in

Table 1. Comparison of frequencies of vibration in rad/s for FGM simply supported beam.

L/h	$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}=\frac{9}{5}$		$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}=\frac{7}{5}$		$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}=\frac{3}{5}$		$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}=\frac{1}{5}$
	Proposed Procedure	Rayleigh Quotient	Proposed Procedure	Rayleigh Quotient	Proposed Procedure	Rayleigh Quotient	Proposed Procedure	Rayleigh Quotient
10	267.3849	267.3849	244.7467	244.7467	191.6077	191.6077	158.4922	158.4922
15	119.1080	119.1080	109.0237	109.0237	85.3526	85.3526	70.6011	70.6011
20	67.0517	67.0517	61.3747	61.3747	48.0491	48.0491	39.7448	39.7448

. A rectangular cross section with dimensions 30 × 50 cm² has been considered and the mechanical properties of the inner material have been assumed to be E₂ = 30 GPa and ρ = 2000 kg/m³. The ratio between the length of the beam L and the height of the cross section h has been assumed to be variable, with values 10, 15, and 20. The ratio between the Young moduli of the outer and inner materials has also been assumed variable, with values 9/5, 7/5, 3/5, and 1/5.

The reported comparison clearly shows that the proposed procedure allows for evaluating the exact frequencies of vibration of a simply supported beam composed of symmetric FGM. Other boundary conditions, not reported for the sake of brevity, have also been analyzed, and the results confirm the accuracy of the approach.

4.1.2. Homogeneous Beams with Different Boundary Conditions

A second validation of the reliability of the proposed procedure can be obtained by evaluating the non-dimensional frequency parameter ${\lambda }_i={\omega }_i\sqrt{{\displaystyle \frac{\rho A{L}^4}{EI}}}$ for single beams with different boundary conditions and different length/cross-section height ratios L/h when the material properties are constant in the cross section. The results obtained by means of the proposed procedure considering the Euler–Bernoulli beam model have been compared to those reported by Lee and Lee in [20] for homogeneous beams.

Table 2. Comparison of frequency parameter for the Euler–Bernoulli beam model.

L/h	C-F			C-C
	λ1	λ2	λ3	λ1	λ2	λ3
10	3.509	21.743	59.801	22.259	60.522	116.21	Proposed procedure
	3.509	21.743	59.802	22.259	60.522	116.21	Lee and Lee
20	3.514	21.960	61.206	22.345	61.379	119.68	Proposed procedure
	3.514	21.961	61.207	22.345	61.379	119.68	Lee and Lee
30	3.515	22.001	61.478	22.361	61.542	120.35	Proposed procedure
	3.515	22.002	61.478	22.361	61.542	120.35	Lee and Lee
L/h	S-S			C-S
	λ1	λ2	λ3	λ1	λ2	λ3
10	9.829	38.845	85.711	15.345	49.095	100.39	Proposed procedure
	9.829	38.845	85.711	15.345	49.095	100.39	Lee and Lee
20	9.860	39.317	88.016	15.400	49.743	103.24	Proposed procedure
	9.860	39.317	88.016	15.400	49.743	103.24	Lee and Lee
30	9.865	39.407	88.463	15.410	49.866	103.80	Proposed procedure
	9.865	39.407	88.464	15.410	49.866	103.80	Lee and Lee

shows the comparison and highlights the accuracy of the obtained results. Comparisons were made for clamped–free (C-F), clamped–clamped (C-C), simply supported (S-S) and clamped–supported (C-S) end conditions. The cross section of the beams has been assumed to have dimensions 30 × 50 cm² and the mechanical properties are E = 70 GPa and ρ = 2702 kg/m³.

In Figure 5, the first three modes of vibration of the differently constrained beams have been reported.

With the aim of validating the proposed procedure, the case shear deformability is taken into account; the fundamental frequency parameter of simply supported homogeneous Timoshenko beams ${\lambda }_1={\displaystyle \frac{{\omega }_1{L}^2}{h}}\sqrt{{\displaystyle \frac{\rho }{E}}}$ has been evaluated for different length/cross-section height ratios L/h and compared in

Table 3. Comparison of fundamental frequency parameter for Timoshenko beam model for simply supported beam.

L/h	Proposed Approach	Su and Banerjee	Sina et al. [21]
10	2.8023	2.8023	2.797
30	2.8438	2.8439	2.843
100	2.8486	2.8496	2.848

to corresponding results reported in [5,21]. The mechanical properties of the beam are E = 70 GPa, ρ = 2702 kg/m³ and Poisson’s ratio υ = 0.3; a rectangular cross section has been considered. For the case of C-C and C-F beams, the results are compared in

Table 4. Comparison of fundamental frequency parameter for Timoshenko beam model for C-C and C-F beams.

L/h	C-C	C-F
5	5.1946	0.9843	Proposed procedure
	5.2138	0.9845	Şimşek
20	6.3495	1.0130	Proposed procedure
	6.3512	1.0130	Şimşek

with those reported by Şimşek in [22]. For the latter comparison, the considered mechanical properties are E = 70 GPa, ρ = 2707 kg/m³ and υ = 0.3. The reported values show a very good agreement with the results available in the literature.

4.2. Eigenproperties of Single Beams Composed of Symmetric Functionally Graded Materials

In this subsection, beams with different boundary conditions are analyzed, assuming a symmetric variation in the mechanical properties in the cross section. In particular, the cross section is rectangular with dimensions 30 × 50 cm², and the values for the mechanical properties of the inner material have been assumed equal to E₂ = 30 GPa and ρ = 2000 kg/m³. Both increasing and decreasing values of the ratio E₁/E₂ with respect to unity have been considered, assuming greater or lower stiffness of the outer material with respect to the inner one. In general, both the Young modulus and density can vary. The influence of the variation in these two parameters on the natural frequencies of vibration can be represented in terms of surfaces, as shown, for example, in Figure 6, for a C-C beam with slenderness ratio L/h = 10. The correspondent surfaces for other constraint conditions and slenderness ratios have not been reported for the sake of brevity. In order to observe in detail the numerical variations in natural frequencies, sections of these surfaces can be reported in terms of tables, assuming constant ρ or E₁/E₂.

Table 5. Results for frequencies of vibration in rad/s for symmetric FGM for Euler–Bernoulli beam model.

L/h	C-F			C-C			$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}$
	ω1	ω2	ω3	ω1	ω2	ω3
3	1040.5150	5837.4920	14,278.2317	6406.4823	15,594.6431	26,398.0259	9/5
	997.4404	5595.8349	13,687.1499	6141.2705	14,949.0652	25,305.2159	8/5
	952.4196	5343.2597	13,069.3627	5864.0763	14,274.3198	24,163.0324	7/5
	905.1624	5078.1375	12,420.8861	5573.1122	13,566.0556	22,964.1092	6/5
	855.2980	4798.3889	11,736.6342	5266.0960	12,818.7176	21,699.0436	1
	802.3407	4501.2878	11,009.9388	4940.0359	12,025.0230	20,355.5072	4/5
	745.6315	4183.1384	10,231.7605	4590.8761	11,175.0990	18,916.7877	3/5
	684.2384	3838.7111	9389.3073	4212.8768	10,254.9741	17,359.2349	2/5
	616.7642	3460.1674	8463.4073	3797.4358	9243.7087	15,647.4029	1/5
10	95.4612	591.4606	1626.7698	605.5223	1646.3690	3161.1801	9/5
	91.5094	566.9757	1559.4258	580.4553	1578.2136	3030.3154	8/5
	87.3790	541.3845	1489.0391	554.2557	1506.9789	2893.5382	7/5
	83.0434	514.5221	1415.1558	526.7546	1432.2055	2749.9664	6/5
	78.4686	486.1777	1337.1966	497.7363	1353.3070	2598.4740	1
	73.6101	456.0751	1254.4016	466.9180	1269.5145	2437.5847	4/5
	68.4074	423.8399	1165.7410	433.9164	1179.7857	2265.2972	3/5
	62.7749	388.9421	1069.7573	398.1890	1082.6456	2078.7792	2/5
	56.5845	350.5877	964.2662	358.9227	975.8836	1873.7863	1/5
15	42.4728	264.8157	735.4772	269.8832	739.3562	1435.6814	9/5
	40.7145	253.8530	705.0303	258.7107	708.7488	1376.2479	8/5
	38.8768	242.3950	673.2079	247.0335	676.7585	1314.1292	7/5
	36.9478	230.3679	639.8046	234.7761	643.1790	1248.9246	6/5
	34.9124	217.6772	604.5585	221.8426	607.7471	1180.1228	1
	32.7507	204.1993	567.1262	208.1068	570.1173	1107.0533	4/5
	30.4359	189.7666	527.0419	193.3979	529.8216	1028.8072	3/5
	27.9299	174.1417	483.6468	177.4741	486.1976	944.0983	2/5
	25.1757	156.9692	435.9534	159.9730	438.2526	850.9987	1/5
20	23.8999	149.3472	416.2481	151.9602	417.4227	813.8886	9/5
	22.9105	143.1646	399.0165	145.6694	400.1424	780.1957	8/5
	21.8764	136.7027	381.0064	139.0944	382.0815	744.9806	7/5
	20.7910	129.9197	362.1016	132.1928	363.1233	708.0161	6/5
	19.6456	122.7626	342.1538	124.9105	343.1193	669.0123	1
	18.4292	115.1615	320.9687	117.1764	321.8744	627.5892	4/5
	17.1266	107.0220	298.2828	108.8944	299.1245	583.2314	3/5
	15.7165	98.2101	273.7230	99.9284	274.4954	535.2098	2/5
	14.1667	88.5254	246.7306	90.0742	247.4268	482.4316	1/5
L/h	S-S			C-S			${\displaystyle \frac{{\mathit{E}}_{\mathbf{1}}}{{\mathit{E}}_{\mathbf{2}}}}$
	ω1	ω2	ω3	ω1	ω2	ω3
3	2855.5125	10,211.2928	19,887.7552	4429.7108	12,775.5312	23,040.8390	9/5
	2737.3016	9788.5717	19,064.4535	4246.3323	12,246.6573	22,087.0079	8/5
	2613.7500	9346.7519	18,203.9548	4054.6686	11,693.8886	21,090.0823	7/5
	2484.0608	8882.9841	17,300.7097	3853.4837	11,113.6604	20,043.6330	6/5
	2347.2168	8393.6310	16,347.6340	3641.1999	10,501.4220	18,939.4530	1
	2201.8845	7873.9238	15,335.4401	3415.7482	9851.2070	17,766.7818	4/5
	2046.2562	7317.3979	14,251.5369	3174.3245	9154.9274	16,511.0323	3/5
	1877.7734	6714.9048	13,078.1072	2912.9599	8401.1376	15,151.5624	2/5
	1692.6021	6052.7334	11,788.4466	2625.7066	7572.6831	13,657.4338	1/5
10	267.3849	1056.6879	2331.5864	417.4227	1335.5276	2730.8517	9/5
	256.3159	1012.9438	2235.0648	400.1424	1280.2402	2617.8015	8/5
	244.7467	967.2232	2134.1823	382.0815	1222.4549	2499.6437	7/5
	232.6029	919.2315	2028.2882	363.1233	1161.7991	2375.6162	6/5
	219.7891	868.5921	1916.5522	343.1193	1097.7969	2244.7463	1
	206.1804	814.8116	1797.8854	321.8744	1029.8248	2105.7587	4/5
	191.6077	757.2210	1670.8115	299.1245	957.0372	1956.9245	3/5
	175.8313	694.8737	1533.2418	274.4954	878.2375	1795.7970	2/5
	158.4922	626.3507	1382.0455	247.4268	791.6326	1618.7096	1/5
15	119.1080	473.8501	1056.6879	186.0130	599.3424	1238.9827	9/5
	114.1772	454.2340	1012.9438	178.3125	574.5312	1187.6920	8/5
	109.0237	433.7315	967.2232	170.2641	548.5990	1134.0840	7/5
	103.6141	412.2106	919.2315	161.8160	521.3786	1077.8129	6/5
	97.9061	389.5024	868.5921	152.9017	492.6565	1018.4375	1
	91.8441	365.3857	814.8116	143.4345	462.1527	955.3791	4/5
	85.3526	339.5603	757.2210	133.2966	429.4879	887.8532	3/5
	78.3249	311.6019	694.8737	122.3214	394.1252	814.7500	2/5
	70.6011	280.8742	626.3507	110.2590	355.2596	734.4057	1/5
20	67.0517	267.3849	598.5720	104.7296	338.2895	702.1165	9/5
	64.2759	256.3159	573.7927	100.3941	324.2852	673.0507	8/5
	61.3747	244.7467	547.8938	95.8627	309.6482	642.6717	7/5
	58.3294	232.6029	520.7084	91.1061	294.2840	610.7836	6/5
	55.1161	219.7891	492.0232	86.0872	278.0723	577.1362	1
	51.7035	206.1804	461.5587	80.7570	260.8549	541.4018	4/5
	48.0491	191.6077	428.9359	75.0491	242.4178	503.1357	3/5
	44.0929	175.8313	393.6186	68.8698	222.4578	461.7090	2/5
	39.7448	158.4922	354.8030	62.0784	200.5208	416.1789	1/5

and

Table 6. Results for frequencies of vibration in rad/s for symmetric FGM for Timoshenko beam model.

L/h	C-F			C-C			$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}$
	ω1	ω2	ω3	ω1	ω2	ω3
3	976.7450	4397.4773	9584.6797	4186.2677	8648.0538	13,854.5002	9/5
	938.0058	4243.4098	9264.3486	4049.1202	8379.5514	13,429.7922	8/5
	897.4863	4082.3999	8930.1297	3906.0296	8100.3366	12,988.3476	7/5
	854.9102	3913.2982	8579.7513	3756.0063	7808.7243	12,527.5718	6/5
	809.9217	3734.5893	8210.2086	3597.7339	7502.4769	12,044.0356	1
	762.0500	3544.2094	7817.3864	3429.3991	7178.5043	11,533.0319	4/5
	710.6499	3339.2333	7395.3979	3248.3926	6832.3253	10,987.7858	3/5
	654.7982	3115.3012	6935.3464	3050.7483	6457.0195	10,397.9306	2/5
	593.0950	2865.4806	6422.7824	2829.9967	6040.9824	9746.2551	1/5
10	94.8493	567.0480	1486.8251	567.9982	1449.9738	2611.4853	9/5
	90.9403	544.2368	1428.8213	545.4725	1394.6966	2515.5404	8/5
	86.8542	520.3822	1368.1448	521.9122	1336.8656	2415.1783	7/5
	82.5646	495.3244	1304.3756	497.1566	1276.0725	2309.6782	6/5
	78.0376	468.8581	1236.9710	470.9992	1211.7878	2198.0977	1
	73.2287	440.7122	1165.2098	443.1662	1143.3040	2079.1675	4/5
	68.0776	410.5171	1088.0998	413.2827	1069.6417	1951.1129	3/5
	62.4986	377.7453	1004.2138	380.8122	989.3830	1811.3341	2/5
	56.3647	341.6006	911.3762	344.9408	900.3504	1655.7839	1/5
15	42.3507	259.6980	703.7623	261.9847	694.8156	1299.8421	9/5
	40.6010	249.0905	675.4856	251.3573	667.2282	1249.4363	8/5
	38.7721	238.0008	645.9152	240.2453	638.3740	1196.7053	7/5
	36.8523	226.3559	614.8525	228.5751	608.0567	1141.2822	6/5
	34.8264	214.0625	582.0428	216.2521	576.0235	1082.6919	1
	32.6747	200.9979	547.1496	203.1519	541.9398	1020.3027	4/5
	30.3702	186.9951	509.7133	189.1049	505.3473	953.2439	3/5
	27.8749	171.8174	469.0787	173.8700	465.5906	880.2612	2/5
	25.1317	155.1082	424.2557	157.0839	421.6756	799.4388	1/5
20	23.8611	147.6924	405.6701	149.4019	402.5527	766.7621	9/5
	22.8745	141.6252	389.1699	143.2889	386.2950	736.2677	8/5
	21.8432	135.2829	371.9181	136.8983	369.2947	704.3738	7/5
	20.7606	128.6240	353.8008	130.1879	351.4385	670.8636	6/5
	19.6183	121.5958	334.6718	123.1044	332.5809	635.4585	1
	18.4051	114.1287	314.3388	115.5770	312.5301	597.7895	4/5
	17.1058	106.1285	292.5401	107.5101	291.0244	557.3518	3/5
	15.6990	97.4613	268.9035	98.7676	267.6912	513.4231	2/5
	14.1527	87.9264	242.8687	89.1451	241.9686	464.9091	1/5
L/h	S-S			C-S			${\displaystyle \frac{{\mathit{E}}_{\mathbf{1}}}{{\mathit{E}}_{\mathbf{2}}}}$
	ω1	ω2	ω3	ω1	ω2	ω3
3	2540.7814	7688.0999	13,383.1858	3355.9825	8221.7511	13,606.4807	9/5
	2443.3654	7422.2724	12,953.3475	3238.2762	7953.5982	13,180.4554	8/5
	2341.4467	7144.3359	12,504.8210	3115.2314	7673.9842	12,736.8627	7/5
	2234.3082	6852.2464	12,034.3978	2985.9610	7381.0235	12,272.8404	6/5
	2121.0202	6543.3005	11,537.7691	2849.2965	7072.2107	11,784.5740	1
	2000.3425	6213.8038	11,008.9315	2703.6474	6744.0942	11,266.7847	4/5
	1870.5614	5858.5037	10,439.1357	2546.7628	6391.7012	10,711.8189	3/5
	1729.2003	5469.5539	9814.9013	2375.2948	6007.4509	10,107.9024	2/5
	1572.4632	5034.4637	9113.9457	2183.9351	5578.9064	9435.4708	1/5
10	263.8169	1005.6455	2111.9527	402.5527	1225.1508	2364.1911	9/5
	252.9963	965.3857	2030.0611	386.2950	1177.2484	2274.9703	8/5
	241.6846	923.2815	1944.3898	369.2947	1127.1392	2181.6289	7/5
	229.8079	879.0493	1854.3418	351.4385	1074.4768	2083.5041	6/5
	217.2717	832.3241	1759.1449	332.5809	1018.8142	1979.7343	1
	203.9517	782.6240	1657.7716	312.5301	959.5574	1869.1654	4/5
	189.6792	729.2903	1548.8062	291.0244	895.8897	1750.1959	3/5
	174.2147	671.3826	1430.2094	267.6912	826.6376	1620.4995	2/5
	157.1986	607.4789	1298.8667	241.9686	750.0132	1476.4898	1/5
15	118.3873	462.9097	1005.6455	182.9546	575.0322	1151.1363	9/5
	113.5069	444.0511	965.3857	175.4667	551.8873	1105.7626	8/5
	108.4056	424.3341	923.2815	167.6388	527.6841	1058.3019	7/5
	103.0502	403.6286	879.0493	159.4193	502.2605	1008.4288	6/5
	97.3985	381.7681	832.3241	150.7427	475.4082	955.7238	1
	91.3949	358.5334	782.6240	141.5227	446.8528	899.6307	4/5
	84.9641	333.6263	729.2903	131.6420	416.2193	839.3863	3/5
	77.9995	306.6230	671.3826	120.9340	382.9734	773.8964	2/5
	70.3410	276.8856	607.4789	109.1485	346.3082	701.4983	1/5
20	66.8218	263.8169	581.3592	103.7477	330.2736	672.0839	9/5
	64.0622	252.9963	557.7678	99.4807	316.8241	645.0739	8/5
	61.1777	241.6846	533.1009	95.0203	302.7624	616.8282	7/5
	58.1497	229.8079	507.1950	90.3375	287.9955	587.1573	6/5
	54.9543	217.2717	479.8404	85.3951	272.4047	555.8177	1
	51.5604	203.9517	450.7610	80.1444	255.8336	522.4882	4/5
	47.9254	189.6792	419.5808	74.5192	238.0691	486.7300	3/5
	43.9893	174.2147	385.7650	68.4258	218.8089	447.9175	2/5
	39.6620	157.1986	348.5075	61.7233	197.5975	405.1056	1/5

report the frequencies of vibration for different values of the ratio E₁/E₂ both for Euler–Bernoulli and Timoshenko beam models. Only the frequencies for different values of E₁/E₂ assuming constant density have been reported because these are the results with a more significant frequency variation.

Slender beams with L/h = 10, 15, and 20 have been considered, but, in order to highlight the influence of shear deformation, the case of a squat beam with L/h = 3 has also been reported.

The results reported in Table 5 and Table 6 and Figure 6 show that, as expected, for decreasing values of the ratio E₁/E₂, the FGM becomes more deformable, and therefore, the natural frequencies decrease.

From the values reported in Table 6, the influence of the shear deformability consists in a reduction in each frequency of vibration, and this influence is more accentuated for lower values of the slenderness ratios L/h.

It is interesting to point out that for all the support conditions and all the considered values of the slenderness ratio L/h, the ratio ${\displaystyle \frac{{\omega }_i\left(E_1/E_2=1\right)-{\omega }_i\left(E_1/E_2=1/5\right)}{{\omega }_i\left(E_1/E_2=1\right)}}$, where ω_i is the i-th frequency of vibration referring to the value of E₁/E₂ considered, remains constant for the Euler–Bernoulli beam model and equal to 27.89%. The corresponding ratio for the Timoshenko beam model varies between 27.86% and 24.67%. In particular, the value of the ratio decreases when ω increases (i.e., for example, when decreasing the value of L/h or when considering a boundary condition with fewer degrees of freedom). Therefore, when shear deformability is considered, the reduction in the frequencies of vibration due to the decreasing E₁/E₂ is higher for more deformable beams, whereas the frequency reduction is not affected by the boundary conditions or the slenderness of the beam when shear deformability is neglected.

4.3. Eigenproperties of a Framed Structure Composed of Symmetric Functionally Graded Materials

In this subsection, frequencies of vibration and mode shapes of a framed structure with elements composed of symmetric FGM are obtained. Columns of the frame have height H and the beam has length L, as shown in Figure 7.

Each element of the structure has a rectangular cross section with base b = 0.3 m and height h = 0.6 m. The mechanical properties of the inner material are E₂ = 30 GPa and ρ₂ = 2000 kg/m³. The results are obtained for different values of H/L = 1, 1/2, and 1/3, with H = 3 m assumed constant. For the Euler–Bernoulli beam model, frequencies for different values of the E₁/E₂ ratio by decreasing the Young modulus of the outer material only and assuming ρ₁/ρ₂ = 1 have been reported in

Table 7. Frequencies for framed structure with varying Young modulus for Euler–Bernoulli beam model.

H/L	$\mathbf{Frequencies}\left[\frac{\mathbf{r}\mathbf{a}\mathbf{d}}{\mathbf{s}}\right]$	$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}=\frac{9}{5}$	$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}=\frac{8}{5}$	$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}=\frac{7}{5}$	$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}=\frac{6}{5}$	$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}=1$	$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}=\frac{4}{5}$	$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}=\frac{3}{5}$	$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}=\frac{2}{5}$	$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}=\frac{1}{5}$
1	ω1	289.9012	277.9001	265.3567	252.1902	238.2974	223.5427	207.7428	190.6379	171.8387
	ω2	1126.9970	1080.3422	1031.5796	980.3946	926.3858	869.0269	807.6044	741.1086	668.0263
	ω3	1822.2656	1746.8285	1667.9832	1585.2211	1497.8931	1405.1483	1305.8329	1198.3145	1080.1461
1/2	ω1	207.5357	198.9443	189.9647	180.5390	170.5933	160.0307	148.7198	136.4746	123.0166
	ω2	371.7867	356.3957	340.3093	323.4238	305.6068	286.6846	266.4218	244.4854	220.3762
	ω3	1027.8928	985.3407	940.8661	894.1821	844.9227	792.6077	736.5865	675.9381	609.2824
1/3	ω1	163.2659	156.5071	149.4430	142.0279	134.2037	125.8943	116.9961	107.3630	96.7757
	ω2	178.6683	171.2719	163.5413	155.4267	146.8644	137.7710	128.0334	117.4915	105.9054
	ω3	510.2879	489.1633	467.0842	443.9084	419.4540	393.4828	365.6716	335.5632	302.4726

, whereas frequencies for different values of ρ₁/ρ₂ by decreasing the density of the outer material only and assuming E₁/E₂ = 1 have been reported in

Table 8. Frequencies for framed structure with varying density for Euler–Bernoulli beam model.

H/L	$\mathbf{Frequencies}\left[\frac{\mathbf{r}\mathbf{a}\mathbf{d}}{\mathbf{s}}\right]$	$\frac{{\mathit{\rho }}_1}{{\mathit{\rho }}_2}=\frac{9}{5}$	$\frac{{\mathit{\rho }}_1}{{\mathit{\rho }}_2}=\frac{8}{5}$	$\frac{{\mathit{\rho }}_1}{{\mathit{\rho }}_2}=\frac{7}{5}$	$\frac{{\mathit{\rho }}_1}{{\mathit{\rho }}_2}=\frac{6}{5}$	$\frac{{\mathit{\rho }}_1}{{\mathit{\rho }}_2}=1$	$\frac{{\mathit{\rho }}_1}{{\mathit{\rho }}_2}=\frac{4}{5}$	$\frac{{\mathit{\rho }}_1}{{\mathit{\rho }}_2}=\frac{3}{5}$	$\frac{{\mathit{\rho }}_1}{{\mathit{\rho }}_2}=\frac{2}{5}$	$\frac{{\mathit{\rho }}_1}{{\mathit{\rho }}_2}=\frac{1}{5}$
1	ω1	211.6497	217.4671	223.7923	230.7035	238.2974	246.6941	256.0459	266.5488	278.4603
	ω2	820.7426	843.7390	868.7835	896.1992	926.3858	959.8437	997.2099	1039.3099	1087.2358
	ω3	1325.2258	1362.7514	1403.6565	1448.4809	1497.8931	1552.7333	1614.0739	1683.3096	1762.2936
1/2	ω1	151.5350	155.6962	160.2202	165.1630	170.5933	176.5971	183.2829	190.7906	199.3034
	ω2	271.3307	278.8103	286.9446	295.8351	305.6068	316.4154	328.4585	341.9905	357.3457
	ω3	748.6100	769.5766	792.4097	817.4038	844.9227	875.4222	909.4826	947.8553	991.5346
1/3	ω1	119.2165	122.4890	126.0467	129.9336	134.2037	138.9246	144.1815	150.0842	156.7767
	ω2	130.4482	134.0323	137.9290	142.1866	146.8644	152.0366	157.7968	164.2656	171.6012
	ω3	372.2256	382.5257	393.7311	405.9826	419.4540	434.3621	450.9818	469.6679	490.8874

.

Figure 8 summarizes the results reported in Table 7 and Table 8 for H/L = 1 in terms of 3D surfaces in which the first three frequencies of vibration are shown with the variation in the ratios E₁/E₂ and ρ₁/ρ₂. As can be observed, the relationships between the frequencies of vibration and the mechanical property ratios are almost linear. It can also be observed that, as expected, the frequencies increase by increasing E₁ (stiffer material), whereas they decrease by increasing ρ₁ (higher values of mass).

The first three mode shapes of the framed structure for H/L = 1 and E₁/E₂ = 4/5 are shown in Figure 9 as an example.

Analogous results for the Timoshenko beam model are reported in

Table 9. Frequencies for framed structure with varying Young modulus for Timoshenko beam model.

H/L	$\mathbf{Frequencies}\left[\frac{\mathbf{r}\mathbf{a}\mathbf{d}}{\mathbf{s}}\right]$	$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}=\frac{9}{5}$	$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}=\frac{8}{5}$	$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}=\frac{7}{5}$	$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}=\frac{6}{5}$	$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}=1$	$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}=\frac{4}{5}$	$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}=\frac{3}{5}$	$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}=\frac{2}{5}$	$\frac{{\mathit{E}}_1}{{\mathit{E}}_2}=\frac{1}{5}$
1	ω1	273.2436	262.3707	250.9990	239.0513	226.4284	212.9990	198.5833	182.9241	165.6319
	ω2	1026.1581	986.2055	944.4121	900.4884	854.0600	804.6293	751.5111	693.7194	629.7500
	ω3	1556.4546	1498.2359	1437.3257	1373.2883	1305.5563	1233.3686	1155.6661	1070.9072	976.7086
1/2	ω1	199.0289	191.0202	182.6453	173.8480	164.5564	154.6753	144.0749	132.5704	119.8823
	ω2	352.8565	338.7469	323.9914	308.4903	292.1159	274.6992	256.0090	235.7153	213.3186
	ω3	951.2203	913.8004	874.6600	833.5304	790.0654	743.8055	694.1199	640.1035	580.3797
1/3	ω1	157.9586	151.5656	144.8807	137.8597	130.4453	122.5625	114.1089	104.9385	94.8316
	ω2	172.7205	165.7303	158.4213	150.7449	142.6388	134.0210	124.7794	114.7549	103.7071
	ω3	484.9226	465.5170	445.2231	423.9042	401.3845	377.4319	351.7290	323.8225	293.0266

and

Table 10. Frequencies for framed structure with varying density for Timoshenko beam model.

H/L	$\mathbf{Frequencies}\left[\frac{\mathbf{r}\mathbf{a}\mathbf{d}}{\mathbf{s}}\right]$	$\frac{{\mathit{\rho }}_1}{{\mathit{\rho }}_2}=\frac{9}{5}$	$\frac{{\mathit{\rho }}_1}{{\mathit{\rho }}_2}=\frac{8}{5}$	$\frac{{\mathit{\rho }}_1}{{\mathit{\rho }}_2}=\frac{7}{5}$	$\frac{{\mathit{\rho }}_1}{{\mathit{\rho }}_2}=\frac{6}{5}$	$\frac{{\mathit{\rho }}_1}{{\mathit{\rho }}_2}=1$	$\frac{{\mathit{\rho }}_1}{{\mathit{\rho }}_2}=\frac{4}{5}$	$\frac{{\mathit{\rho }}_1}{{\mathit{\rho }}_2}=\frac{3}{5}$	$\frac{{\mathit{\rho }}_1}{{\mathit{\rho }}_2}=\frac{2}{5}$	$\frac{{\mathit{\rho }}_1}{{\mathit{\rho }}_2}=\frac{1}{5}$
1	ω1	201.1193	206.6449	212.6525	219.2165	226.4284	234.4025	243.2830	253.2559	264.5652
	ω2	757.1785	778.2845	801.2598	826.3976	854.0600	884.7000	918.8931	957.3837	1001.1544
	ω3	1156.9441	1189.3036	1224.5392	1263.1036	1305.5563	1352.5977	1405.1186	1464.2721	1531.5826
1/2	ω1	146.1774	150.1904	154.5532	159.3198	164.5564	170.3458	176.7926	184.0316	192.2395
	ω2	259.3760	266.5211	274.2912	282.7831	292.1159	302.4383	313.9384	326.8588	341.5181
	ω3	700.4027	719.9347	741.1974	764.4625	790.0654	818.4258	850.0771	885.7092	926.2332
1/3	ω1	115.8806	119.0609	122.5184	126.2956	130.4453	135.0329	140.1412	145.8769	152.3799
	ω2	126.6978	130.1782	133.9621	138.0964	142.6388	147.6610	153.2541	159.5350	166.6574
	ω3	356.2458	366.0919	376.8023	388.5114	401.3845	415.6284	431.5049	449.3519	469.6138

.

As in the case of the beam with different boundary conditions of the previous paragraph, changing the shear deformability reduces the stiffness and, consequently, the frequencies of vibration.

5. Conclusions

Inhomogeneous beams in which the mechanical properties vary along the beam thickness according to a symmetric distribution have been considered. In particular, the Young modulus and the mass density of the inner material have been assumed to be constant, while, with a symmetric parabolic variation, an outer material with variable characteristics has been introduced. The eigenproperties of single beams with different support conditions and slenderness ratios have been evaluated by means of the Wittrick and Williams algorithm applied in conjunction with the dynamic stiffness matrix for both Euler–Bernoulli and Timoshenko beam models. As expected, the reduction in mechanical characteristics of the outer material involves a reduction in vibration frequencies as the structures become more deformable. Considerations of the percentages of reduction for the considered beams have been reported, allowing us to differentiate the behavior of Euler–Bernoulli beams from Timoshenko beams and also to identify which support condition and slenderness ratio mainly affect the frequency reduction. The proposed approach can also be applied to framed structures, and a simple example of an elementary frame has been reported in this paper. The developed study aims to make a contribution to the design of structures using FGMs, since these materials allow for planning advanced structures tailored to specific dynamic conditions.

It must be pointed out that, although the proposed approach provides an exact mathematical solution within the adopted classical beam theories, the fundamental hypotheses of these theories could not be suitable for considering the presence of complex materials which could require the adoption of more sophisticated approaches based on the use of solid FEM elements.