Bending and free vibrational analysis of bi-directional functionally graded beams with circular cross-section

The bending and free vibrational behaviors of functionally graded (FG) cylindrical beams with radially and axially varying material inhomogeneities are investigated. Based on a high-order cylindrical beam model, where the shear deformation and rotary inertia are both considered, the two coupled governing differential motion equations for the deflection and rotation are established. The analytical bending solutions for various boundary conditions are derived. In the vibrational analysis of FG cylindrical beams, the two governing equations are firstly changed to a single equation by means of an auxiliary function, and then the vibration mode is expanded into shifted Chebyshev polynomials. Numerical examples are given to investigate the effects of the material gradient indices on the deflections, the stress distributions, and the eigenfrequencies of the cylindrical beams, respectively. By comparing the obtained numerical results with those obtained by the three-dimensional (3D) elasticity theory and the Timoshenko beam theory, the effectiveness of the present approach is verified.

If a Timoshenko or higher-order beam model is adopted to investigate the free vibrations of axially FG beams, two coupled governing differential equations with variable coefficients will be derived, which makes it difficult to obtain the exact solutions due to the arbitrary gradient changes. Many researchers have introduced different numerical methods to solve such problems. Based on the finite element method, Shahba et al. [11] discussed the free vibrations of arbitrarily tapered and axially FG Timoshenko beams under various boundary conditions. Huang et al. [12] introduced a unified method to analyze the free vibrations of Timoshenko beams with arbitrarily axial material parameters by transforming the governing equations into a system of linear algebraic equations under different boundary conditions. Based on the differential transformation element method, Rajasekaran and Tochaei [13] studied the free vibrations of axially FG Timoshenko beams. Tang et al. [14] studied the free vibrations of non-uniform FG Timoshenko beams, and derived the frequency equations in closed form, where the material properties were assumed to vary in a unified exponential law. Cao et al. [15] used the asymptotic development method to discuss the vibrational behaviors of uniform axially FG beams under different boundary conditions. Zhang et al. [16] presented an effective approximation to investigate the free vibrations of axially FG beams based on the Euler beam theory and the Timoshenko beam theory, respectively.
FG beams with bi-directional variations of material parameters have attracted increasing attention by many scholars recently [17][18] since bidirectional FG structures have better mechanical properties and can produce more effective resistance effects and realize the integrated design of materials and structures. With the symplectic method and the Hamiltonian state space approach, Zhao et al. [19] derived the elasticity bending solutions of the bi-directional FG beams under arbitrary lateral loadings. Based on the classical hairbrush hypothesis, Pydah and Sabale [20] obtained an analytical solution for the flexure of bi-directional FG circular beams subject to various tip loads. Armagan [21] presented a symmetric smoothed particle hydrodynamics method to derive the bending solutions of bi-directional FG beams under various boundary conditions with the Euler-Bernoulli, Timoshenko, and Reddy-Bickford beam theories, respectively. Li et al. [22] proposed a meshless total Lagrangian corrective smoothed particle method to study the static behavior of bi-directional FG beams with the exponential or power-law gradient assumption. Based on the Euler-Bernoulli and Timoshenko beam theories, Simsek [23] used the energy approach to analyze the free and forced vibrations of bi-directional FG beams subject to the action of moving loads. With the assumption of an exponential-law material distribution, Deng and Chen [24] introduced the variable substitution method and the dynamic stiffness matrix method to discuss the dynamic characteristics of bi-directional FG Timoshenko beams. Huynh et al. [25] used the isogeometric analysis method to study the free vibrations of bi-directional FG Timoshenko beams with power-and exponential-law models of material distributions, respectively. Nguyen et al. [26] introduced a finite element method to discuss the vibrations of bi-directional FG Timoshenko beams under the action of a moving concentrated loading. Based on the third-order shear deformation theory, Karamanli [27] investigated the free vibrations of bi-directional FG beams with exponentially graded material properties. Lal and Dangi [28] used the first-order shear deformation theory and Eringen's nonlocal elasticity theory to study the vibrational behavior of bi-directional FG non-uniform Timoshenko nanobeams, where the material properties varied in both the axial and the thickness directions according to the power-and exponential-law distributions, respectively.
It should be noted that most of the above mentioned papers are limited to FG beams with a rectangular cross-section. However, there are very few bending and eigenfrequency results reported for FG circular cylinders with radially or/and axially varying material inhomogeneities. Abadikhah and Folkow [29] adopted the three-dimensional (3D) elastodynamic theory to study the dynamic behaviors of simply supported FG cylinders with power-and exponential-law distributions of the material properties. Zhang et al. [30] introduced a higher-order shear deformation based spectral element model to calculate the stochastic natural frequency of radially FG beams with a circular cross-section. To the author's best knowledge, the bending and free vibrations of axially or two-dimensional (2D) FG cylindrical beams have not yet been reported. Solid circular beams are common structural elements in macro and micro scales. The analysis on the bending and free vibrations of FG cylindrical beams is useful for the precise design of circular cylindrical composite beams.
The purpose of this paper is to introduce a high-order circular beam model for studying the bending and free vibrations of 2D FG circular beams, where the shear-free surface condition is identically satisfied. The analytical bending solutions can be obtained for general cases with different boundary conditions. By use of a new auxiliary function, a single governing equation is furthermore derived in the free vibration analysis of 2D FG circular beams. The Chebyshev polynomials are used to get the characteristic equations. Illustrative examples are given to show the effects of the gradient parameters on the deflections, the stresses, and the natural frequencies, respectively.

Theory and formulation
An elastic cylindrical beam of the length L and the radius R is considered, subject to the action of the arbitrary transverse loading q(x) (see Fig. 1). The cylindrical beam is inhomogeneous, where the material properties are assumed to vary simultaneously along the length and radial directions. For such a structure, the Cartesian coordinates (x, y, z) and polar cylindrical coordinates (x, r, θ) are both used, where y = r cos θ, and z = r sin θ. With the high-order cylindrical beam model [31] , we can express the displacement field as where u x , u y , u z , and u r are the elastic displacement components in the x-, y-, z-, and rdirections, respectively. u 0 is the axial displacement of any point on the neutral axis. ϕ(x, t) and w(x, t) are, respectively, the cross-sectional rotation and transverse displacements of the beam.
Based on the small deformation assumption, the relations between the strain and displacement components can be derived as where ε and γ are the normal strain and the shear strain, respectively. According to the Saint-Venant principle, we can ignore the effects of σ yy , σ zz , and τ yz for a cylindrical beam [32] . As a result, the 3D constitutive equations change to where σ and τ are the normal stress and the shear stress, respectively. E(x, r) and G(x, r) are the Young's modulus and the shear modulus, respectively, which are 2D continuous functions dependent on both x and r. From Eq. (12), it can be easily found that the shear stress τ xr is zero at the circumference z 2 + y 2 = R 2 , which means that the shear-free surface condition is identically satisfied in this model. In this thesis, we mainly consider the bending and vibration in the vertical plane, while neglect the horizontal motion in the Oxy-plane. As a result, the equilibrium equations for the plane problem change to where ρ(x, r) is the radial and length-dependent mass density of the beam. Integrating both sides of Eq. (13) over the cross-sectional area A yields where is the axial normal force, and It should be noted that the following integral results have been used to derive Eq. (15): In the absence of the applied axial force N , Eqs. (15) and (16) change to Now, applying Eqs. (9)-(11) into Eq. (13), multiplying both sides of Eq. (13) by z, and integrating both sides over the cross-sectional area A yield where In the process of deriving the above formulae, the following integral results are used: Furthermore, integrating both sides of Eq. (14) over the cross-sectional area A yields Up to now, we finally get two coupled dynamic governing equations (22) and (28) for a 2D FG cylindrical beams in terms of the transverse deflection w and the rotation ϕ. The stress resultants of the shear force Q and the bending moment M are defined as In the Timoshenko beam theory, the formula of the shear force is Q = A κτ xz dA, where κ is the shear correction factor. Because the shear stress is assumed to be uniform in this theory, we naturally have to introduce a coefficient κ to relax the traction-free surface condition.

Bending of 2D FG cylindrical beams
For the bending of 2D FG circular beams, by removing all time-dependent terms, the governing equations (22) and (28) can be simplified as Integrating both sides of Eq. (32) with respect to x from 0 to x yields Substituting Eq. (32) into Eq. (31) and then integrating both sides of Eq. (31) with respect to x from 0 to x yield Derivating both sides of Eq. (33) yields Solving the above linear equations (34) and (35), we can easily get the solution of d 2 w dx 2 as d 2 w where Then, integrating two times from x = 0 to x on both sides of Eq. (36) yields where A i (i = 1, 2, 3, 4) are unknown integration constants in the general solution, which can be determined by means of the boundary conditions. Inserting Eq. (37) into Eq. (33) yields In addition, applying the above expressions of w and ϕ into Eqs. (29) and (30) yields Substituting Eqs. (37) and (38) into the normal stress σ xx and the shear stresses τ xz and τ xz yields It is obvious that the above physical quantities will be uniquely determined when the four constants A j are obtained. For instance, we consider a simply 2D FG cylindrical beam, where the corresponding boundary conditions can be stated as In view of w = 0 and M = 0 at x = 0 and L, from Eqs. (37) and (40), one can derive four linear equations for A i as Therefore, A 2 = A 4 = 0, and A 1 and A 3 can be obtained by solving the linear equations of Eqs. (47) and (48) as For other boundary conditions, we can use the similar method to determine the four constants which are omitted here.

Free vibrations of 2D FG cylindrical beams
In this section, we will study the free vibrations of 2D FG cylindrical beams. To this end, applying q = 0 into Eqs. (22) and (28) yields the governing equations as Since the governing equations (51) and (52) are two coupled higher-order differential equations, it is almost impossible to find the closed-form solutions of Eqs. (51) and (52) for different inhomogeneities. Therefore, it is much desired to find a numerical method for dealing with the free vibrations of such beams effectively. If we can transform two coupled governing differential equations (51) and (52) into a single governing equation, it is going to be easier for our subsequent analysis. For this purpose, the deflection w and the rotation ϕ are assumed to be where F (x, t) is a new auxiliary function. Apply Eqs. (53) and (54) 53) and (54) into Eq. (51), we can immediately get a single fourth-order partial differential governing equation as where The advantage of this treatment is that the physical quantities of interest can be expressed by the auxiliary function F (x, t). For example, the deflection w and the rotation ϕ can be determined by means of Eqs. (53) and (54), respectively. Keeping Eqs. (29) and (30) in mind, the bending moment M and the shear force Q can be derived in terms of F as In order to analyze the free vibrations of the 2D FG circular cylindrical beams, the auxiliary function F can be expressed as where ω is the circular frequency. Plugging Eq. (58) into Eq. (55) yields The key of the problem lies in exactly calculating the natural frequencies w by solving the fourth-order differential equation (59) associated with the corresponding boundary conditions.

Solution of the resulting eigenproblem
We know that Chebyshev polynomials play an important role in numerical calculation applications. In the following, a simple approach called the Chebyshev polynomial expansion method will be proposed to determine the eigenvalues, which are related to the natural frequencies. The solution f (x) of Eq. (59) can be expressed approximately in terms of the shifted Chebyshev polynomials with respect to the length coordinate as where a i and N are unknown coefficients and the orthogonal polynomial order, respectively. T i (x) are the first kind of Chebyshev polynomials over the interval [0, L], which are defined by the following recurrence relations: It should be noted that the Chebyshev polynomial expansion (60) must satisfy the governing equation (59) and the boundary conditions at the ends of the beam simultaneously. Now, putting the Chebyshev expansion (60) into the governing equation (59) for each case, we have where T (k) Multiplying both sides of Eq. (63) by factors T j (x) (j = 0, 1, 2, · · · , N −4) and then integrating both sides from 0 to 1 with respect to x yield N − 3 linear equations of unknown coefficients a i as where Next, take clamped-clamped (CC) 2D FG circular cylindrical beams as an example. Substitute the expansion (60) to the clamped conditions w = ϕ = 0 at x = 0. Then, the other two linear equations of a 0 , a 1 , · · · , a N will be obtained as On the other hand, applying the clamped conditions at From Eqs. (64)-(68), we can find that a series of of linear algebraic equations for unknown coefficients a i have been derived, which forms a system as To make sure that the resulting linear system has a non-zero solution, the determinant of the coefficient matrix of the system should be zero, i.e., It is clear that the obtained equation (70) is actually a polynomial in the natural frequency ω. With the help of scientific software, we can easily obtain the multi-roots of its positive solutions, which correspond to different orders of natural frequencies. For other relevant familiar end conditions, such as simply-simply (SS), clamped-simply (CS), and clamped-free (CF), we can also apply the Chebyshev expansion (60) to the corresponding boundary conditions for getting the last four linear equations, and making the first N − 3 linear equations as the same form as Eq. (64), where the details are omitted here. This treatment could give rise to a great simplification on studying the free vibrations of 2D FG circular cylindrical beams with different end supports.

Model validation and convergence studies
In order to verify the introduced cylindrical beam model, we first present some results applied to static and vibrational problems on homogeneous cylindrical beams. In general, however, it is rather difficult to get the exact bending elasticity solutions of the cylindrical beam for complex loads under different boundary conditions [32] , which can be completely solved only for some certain special cases. For example, consider a homogeneous cantilevered circular beam subjected to a concentrated force P at the free end, where the boundary conditions are x = L : M = 0, Q = P.
In this case, one can use the Saint-Venant semi-inverse method to get the exact solutions [33][34] : τ xy = − (2ν + 1)P yz 4(ν + 1)I , It is easily found that the expression of the normal stress σ xx in Eq. (76) is identical to the exact solution (73). In addition, the distributions are similar between two results of the shear stresses. If we do not take the effect of Poisson's ratio into account, the present solutions (77) and (78) are equivalent to the exact results. Moreover, it should be noted that the Timoshenko beam theory cannot give the distribution of the shear stress at the cross-section, but only presents the shear force resultant at the cross-section. For examining the convergence and effectiveness of the Chebyshev polynomial expansion method, next, we will study the free vibration of a homogeneous cylindrical beam with a uniform circular cross-section. The convergence and verification studies are performed by using different numbers of terms in the Chebyshev polynomial expansion. By using the proposed method, we calculate the dimensionless natural frequencies Ω n = ω n R ρ/G for CF cylindrical beams with various aspect ratios and truncation orders N , where the results are noted in Table 1. For comparison, the 3D solutions calculated by the Ritz method [35] and the Chebyshev-Ritz method [36] are also listed in this table. It can be found that the presented approach has a fast convergence rate. With an increase in the aspect ratio L/R, there is good agreement between our numerical results and the 3D solutions.
We further use the introduced model to analyze the free vibrations of linearly tapered cylindrical beams, where the Young's modulus and the mass density of the beams keep constants, and the radius R of the circular cross-section is assumed to be where R 0 and R L are the corresponding radii of the cross-section at the ends x = 0 and x = L, respectively. Based on the 3D theory of elasticity, Kang and Leissa used the Ritz method to investigate the free vibration frequencies and mode shapes of such tapered beams with a circular cross-section for isotropic materials [37] . For different ratio values of L/(R 0 + R L ), the first three dimensionless natural frequencies Ω n = ω n L ρ/G (n = 1, 2, 3) are evaluated for free-free (FF) boundary conditions, which are displayed in Table 2. We set Poisson's ratio ν = 0.3 and the aspect ratio R L /R 0 = 3 in this example. Compared with the results derived previously by the Euler-Bernoulli beam theory [37] , it is obvious that our results agree well with the 3D solutions. This indicates that the suggested high-order model and numerical method can effectively handle the dynamic analysis of cylindrical beams.

Bending of bi-directional FG cylindrical beams
In this section, we will discuss the bending behavior of bi-directional FG cylindrical beams with a uniform circular cross-section. A typical model called power-law distribution is analyzed, where the material properties are assumed to be where p Set p x = 0, which means that the Young's modulus, the shear modulus, and the mass density of the beam depend only on r. In this case, the cylindrical beam becomes the radial-dependent FG beam. Keep p r = 0, which means that the beam becomes the axially FG cylindrical beam. If p x = 0 and p r = 0, the beam becomes a homogeneous cylindrical beam with a uniform cross-section. Subject to a uniformly transverse loading, i.e., q(x) = q 0 , the dimensionless transverse deflections W (x) = 100EmR 3 q0L 3 ω(x) of bi-directional FG cylindrical beams are calculated for various values of the material gradient indices (p x , p r ). The maximum results of the transverse deflection W (L/2) under the SS boundary condition are listed in Table 3. For further comparison, the results of the dimensionless deflection obtained by the Timoshenko beam theory are also listed in Table 3 with the shear correction factor κ = 6/7. We can find that the results calculated by these two beam models appear to have good consistency. It can also be observed that an increase in the power-law gradient parameter p x or p r implies a stiffness-hardening effect of the FG beam, which can control the dimensionless maximum transverse deflection to increase monotonically. The results of the dimensionless maximum transverse deflection with CC supported ends are tabulated in Table 4.  The variations of the dimensionless deflection 100EmR 3 q0L 3 ω(x) with the dimensionless coordinate x/L are illustrated in Figs. 2 and 3 for CS and CF bi-directional FG cylindrical beams, respectively, in which p r = 1 and L/h = 10. It is observed from Figs. 2 and 3 that the maximum values of the transverse deflections for CS and CF beams are not located in the middle positions of the beams because of the asymmetrical boundary conditions.  The results of the normal stress distributions under CF boundary conditions have the same variation tendencies, as shown in Figs. 6 and 7, where p r = 1 and L/h = 10. We derive the explicit expression of the shear stress under the SS and CF boundary conditions based on the high-order cylindrical beam model as

Natural frequencies of 2D FG cylindrical beams
Till now, the free vibrations of FG circular cylinders have seldom been discussed. Abadikhah and Folkow [29] adopted the Fourier and power series expansions to calculate the eigenfrequencies of simply supported cylinders by using the 3D elastodynamic theory for radially varying material inhomogeneities. Zhang et al. [30] proposed a higher-order shear deformation based spectral element model to calculate the natural frequency of FG cylindrical beams, where the material properties were assumed to vary along the r-axis. However, there has been no published work on discussing the dynamic behaviors of the axially or bi-directional FG cylindrical beam.
In order to investigate the effects of the power-law indeices p x and p r in Eqs. (80)-(82) on the natural frequency of 2D FG cylindrical beams, the results of the first three dimensionless natural frequencies Ω n = ω n R ρ m /G m under the SS and CF boundary conditions are presented with N = 12 in Tables 5 and 6, respectively. Some natural frequency results cited from Ref. [30] are also listed in Tables 5 and 6 with p x = 0. Very good agreement between the present results and the existing numerical results [30] can be observed for the SS and CF 2D FG cylindrical beams. For a given value of the index p r , one can find from Table 5 that the first three natural frequencies of the SS beam decrease with the increase in the gradient index p x . When p r increases, the first two dimensionless natural frequencies of the SS beam for a given value of the index p x increase, reach the maximum values, and then decrease.
As the last example, we consider an exponential-law gradient model, where the material properties can be represented in an exponential form as  where E m , G m , and ρ m denote the Young's modulus, the shear modulus, and the mass density value at the reference point (0, 0, 0), respectively. α x and α r are the gradation indices.
To examine the effects of the material gradient indices α x and α r on the vibrational behaviors of 2D FG cylindrical beams, the first three dimensionless natural frequencies Ω n = ω n R ρ m /G m (n = 1, 2, 3) are evaluated for SS and CC boundary conditions, where the results are tabulated in Tables 7 and 8, respectively. We set ν = 0.3 and L/h = 10 in this example. It is obvious from Table 7 that if the exponential-law parameter α x increases, the  natural frequencies decrease, but the increase in α r can result in the increase in the natural frequencies for SS beams. For CC beams, the natural frequencies gradually increase when the index α x or α r increases.

Conclusions
Based on the high-order circular beam theory, where the shear deformation and rotary inertia are both considered without introducing the shear correction factor, the bending and free vibrations of 2D FG cylindrical beams are discussed. For any radial/axial nonhomogeneity of material properties, the coupled governing differential equations of the deflection and rotation are derived. The analytic bending solutions are derived in a closed form for different boundary conditions. By expanding the auxiliary function into shifted Chebyshev polynomials, the characteristic polynomial equations in natural frequencies are obtained. Higher accuracy can be achieved through increasing the order of Chebyshev polynomials. By comparing with the exact 3D solutions, we present some static and vibrational results of homogeneous cylindrical beams to verify the introduced cylindrical beam model. The effects of gradient indices on the transverse deflections, the stresses, and the natural frequencies of 2D FG cylindrical beams are studied. Compared with the results calculated by the Timoshenko beam theory and other exiting numerical results, the accuracy and effectiveness of the introduced approach can be confirmed.