1 INTRODUCTION
Elastic flexural members whose crosssectional profile changes partially or gradually along their length, known as nonprismatic beams, are widely spread in many engineering applications, due to their ability in improving both strength and stability, satisfaction aesthetic necessities and optimization weight of structures. By developing manufacturing process, nonprismatic beams are now adopted with different materials such as wood, steel, and composite materials. Functionally graded materials (FGMs), as a new class of advanced materials, are made up by changing the composition of two or more different materials gradually and smoothly in any desired direction. Thus, engineers can produce the structures with favorable stability and manage the distribution of material properties. In addition, functionally graded materials (FGMs) can overcome some disadvantages and weaknesses of laminated composites such as delamination and stress concentration, due to smooth variations in material properties. During the last two decades, the use of nonprismatic beams made of functionally graded (FG) materials has been increasing in complicated mechanical components such as turbine blades, rockets, aircraft wings, and space structures due to their conspicuous characteristics such as high strength, thermal resistance, and optimal distribution of weight.
The accurate evaluation of the critical buckling loads and natural frequencies of elastic members are considered as one of important issues in designing of different structures. So far, a large number of studies have been conducted regarding the stability and free vibration of beams with constant crosssection. Regarding axially graded EulerBernoulli beams with an arbitrarily variable crosssection, the stability and dynamic analyses become more complex due to the presence of variable coefficients in the governing fourth order differential equation. Thus, there are no closedform solutions for exact estimation of buckling loads and natural frequencies of nonprismatic beams made up of axially nonhomogeneous material. Therefore, different numerical techniques, especially based on the Finite Element Method (FEM), have been proposed for investigating the linear stability and free vibration behavior of nonuniform beams. For example, exact solution of special types of tapered columns presented by ^{Gere and Carter (1962}) for the first time. ^{Ermopulos (1977}) calculated the critical loads and the corresponding equivalent buckling lengths of framed nonuniform members based on slope deflection method. In addition, finite difference method was applied by ^{Iromenger (1980}) to determine critical buckling load of tapered and stepped columns. Further, a numerical technique based on Galerkin’s method was suggested by ^{Jategaonkar and Chehil (1989}) to analyze the free vibration behavior of beams with varying section properties. ^{Arbabi and Li (1991}) presented a semi analytical approach for measuring buckling load of columns with stepvarying profiles. ^{Sampaio and Hundhausen (1998}) presented a mathematical model and analytical solution for buckling of inclined beamcolumns. Only one type of boundary condition (pinnedends) was analyzed based on their method and its exact solution was shown. ^{Rahai and Kazemi (2008}) formulated a new approach for the linear stability problem of tapered column members. The exact buckling load was calculated by utilizing modified vibrational mode shape (MVM) and energy method. Coşkun and ^{Atay (2009}) used variational iteration method to evaluate the critical buckling load of elastic columns with variable crosssections. In another study, Atay (2009) applied a new approach, known as Homotopy perturbation method, to solve the buckling problem of nonprismatic columns. Further, ^{Okay et al. (2010}) derived buckling loads and mode shapes of a heavy column by applying the variational iteration method.
Regarding the abovementioned studies, all are merely related to the buckling and dynamic analyses of homogeneous EulerBernoulli beam. However, in the recent years, the investigation of linear stability and free vibration behavior of beams with gradation in material properties along the member axis and its application in the construction of members has attracted a lot of attention. Based on the stepreduction method, ^{Tong and Tabarrok (1995}) developed a new approach for the free and forced vibration of nonhomogeneous Timoshenko beams with nonuniform crosssections. In addition, an elasticity solution based on functionally graded EulerBernoulli beam under static transverse loads in which the elastic modulus of the beam varied exponentially in the thickness direction of the member has been proposed by ^{Sankar (2001}). Further, ^{Chakraborty et al. (2003}) introduced a new finite element solution based on the firstorder shear deformation theory to investigate the thermoelastic behavior of FG beam structures. ^{Wu et al. (2005}) adopted a semiinverse technique to obtain the solution to the dynamic equation of simply supported axially functionally graded beams. Furthermore, ^{Singh and Li (2009}) determined the natural frequencies of a nonuniform beams by using a newly developed numerical method. Based on a new beam theory different from the traditional firstorder shear deformation beam assumption, ^{Sina et al. (2009}) studied free vibration behavior of functionally graded beams by utilizing an analytical technique. ^{Simsek and Kocatürk (2009}) investigated free vibration characteristics and the dynamic behavior of a FG simply supported beam subjected to concentrated moving harmonic load. Additionally, Li et al. (2013) studied the free vibration behavior of axially functionally graded beams by assuming material properties of the beam including Young’s modulus and density varying exponentially. Moreover, the linear static and dynamic analyses of homogeneous and axially nonhomogeneous tapered EulerBernoulli beams were performed by ^{Shahba et al. (2013}a, b) by using a finite element model. They derived new shape functions in terms of basic displacement functions (BDFs), which are based on energy method. ^{Arefi (2013}) suggested an analytical method to survey the nonlinear thermoelastic analysis of a functionally graded piezoelectric (FGP) thickwalled cylinder under thermal, mechanical and electrical loads. ^{Pradhan and Chakraverty (2013}) used RayleighRitz method to analyze the free vibration of Euler and Timoshenko functionally graded beams. Arefi (2015) studied electromechanical stability of a functionally graded circular plate integrated with two functionally graded piezoelectric layers under radial compressive. ^{Khan et al (2016}) proposed a finite element solution for free vibration and static analyses of beam made up of functionally graded materials (FGMs), based on efficient zigzag theory (ZIGT). Based on the modified couple stress theory, mechanical behavior of nonuniform bidirectional functionally graded beam sensors was studied in detail by ^{Khaniki and Rajasekaran (2018}). Li et al. (2018) conducted instability analysis of a microscaled bidirectional functionally graded (FG) beam having rectangular crosssection by employing the Generalized Differential Quadrature Method (GDQM). Recently, ^{Soltani and Asgarian (2019}) performed the stability analysis of cantilever axially functionally graded nonprismatic Timoshenko beam through a new finite element model based on power series approximation.
The investigations of elastic buckling load and natural frequency of members resting on elastic foundations is also one of the complicated and significant problems in designing many components related to soilstructure interaction such as the foundation of buildings, the pipelines embedded in soil, highway pavements, and the like. In this regard, different types of elastic foundation models like Winkler, Pasternak and Vlasov were presented. The Winkler type of foundation is the most used mechanical model for solving the problems mentioned above. In this model, the soil is considered as the limiting case of an infinitely dense distribution of translational springs with linear behavior, which are independent of each other. However, modelling of the elastic foundation by Winkler’s theory was inadequate in several problems since this model overlooks the soil cohesion. In order to improve this weakness, various twoparameter elastic foundation models such as WinklerPasternak foundation were developed. In this model, an additional layer is considered in the widely used Winkler model in order to accomplish the effect of shear interactions between the springs.
During the last decades, the study of the stability and free vibration of beam resting on an elastic foundation has received a lot of attentions by different authors. Among the first investigations on this topic, the most important one is the study of ^{Eisenberger and Clastornik (1987}), in which the polynomial functions were adopted for buckling and free vibration analyses of elastic beams resting on variable Winkler elastic foundation. In order to study the free vibration behavior of beam resting on elastic foundation, a modified Vlasov model was applied by ^{Ayvaz and Ozgan (2002}). In addition, ^{Kim et al. (2005}) developed the dynamic stiffness matrix based on power series method to study the free vibration behavior of thinwalled beam with nonsymmetric crosssection resting on twoparameter elastic foundation by considering shear deformation. Static and dynamic stiffness matrices of members with variable crosssection resting on WinklerPasternak elastic foundations were derived by ^{Girgin and Girgin (2005}). This numerical approach was developed based on the wellknown Mohr method. Further, ^{Ruta (2006}) used Chebyshev polynomial method to solve the coupled system of equilibrium differential equations of Timoshenko beams with nonuniform crosssection resting on two parameter elastic foundation. Additionally, ^{Malekzadeh and Karami (2008}) proposed a mixed method based on finite element approach and differential quadrature method for stability and free vibration analyses of beam members on elastic foundation. A new finite element formulation for the nonlinear free and forced vibration analyses of nonprismatic Timoshenko beams resting on two parameter foundations was introduced by ^{Zhu and Leung (2009}). The proposed method is based on a hierarchical finite element method by accounting the effects of transverse shear deformation and rotatory inertia. He’s variational method is adopted by ^{Fallah and Aghdam (2011}) to obtain analytical expressions for large amplitude free vibration and postbuckling analysis of the functionally graded beams under an axial load and resting on nonlinear elastic foundation. ^{Mirzabeigy (2014}) presented a semianalytical technique based on the differential transformation method to obtain the dimensionless natural frequencies of nonuniform beams resting on an elastic foundation. Furthermore, ^{Tsiatas (2014}) suggested a new influential approach to exactly determine stiffness and mass matrices of nonuniform EulerBernoulli beam from inhomogeneous linearly elastic material resting on an elastic foundation. Adomian decomposition method (ADM) was applied by ^{Hassan and Nassar (2015}) in order to determine the critical buckling loads in the static analysis and the natural frequencies for free vibration behavior of Timoshenko beam resting on a two parameter elastic foundation. By contemplating the impact of elastic foundation and semirigid end conditions, buckling analysis of axially functionally graded EulerBernoulli beam having nonuniform crosssection was surveyed in detail by ^{Shvartsman and Majak (2016}). Based on nonlocal strain gradient elasticity theory, an investigation on the vibration characteristics of viscoelastic functionally graded EulerBernoulli beam embedded in a viscoelastic foundation by considering surface and thermal effects was accomplished by ^{Ebrahimi and Barati (2017}).
By considering the studies mentioned above, it seems that the researchers investigated the problem of buckling and natural frequency only for special types of members. However, in the present study, an improved and efficient hybrid numerical method to exactly evaluate the critical buckling loads and free transverse frequencies is proposed for any types of AFG members with linear, polynomial or exponential variation of mechanical properties and resting on uniform elastic foundation. The material properties of the nonprismatic beam are assumed to be graded smoothly along the beam axis by a powerlaw distribution of volume fractions of metal and ceramic, while the material features are constant in the direction of the thickness. In addition, the material behavior is linearelastic.
In fact, the current article aimed to calculate the critical buckling loads and natural frequencies concurrently for FGM members based on the power series expansions mixed with the RayleighRitz method. It is worth noting that EulerBernoulli beam theory is adopted in this study in order to analyze the beams with uniform or nonuniform crosssection, in which the effect of flexural deformation is consider while the influence of shear deformation and rotatory inertia is negligible. The followings are considered as the main steps for conducting the present study:
1 The power series expansions is implemented to facilitate the solution of the governing motion differential equation of nonprismatic member with variable coefficients. The variable terms in equation of motion including crosssectional area, moment of inertia, mechanical properties such as Young’s modulus of elasticity and density of material, as well as the vertical deflection function, is expanded into power series form. Then, the natural frequencies of the considered members is derived by imposing the boundary conditions and solving the eigenvalue problem.
2 Accordingly, the explicit expression of the deflected shape of the member is constructed by imposing both natural and geometric boundary conditions in polynomial solution related to the obtained fundamental circular frequency.
3 Finally, the critical buckling load of the considered beam with varying crosssection is evaluated based on the similarities between vibration and buckled deformation shapes of elastic members and adopting the RayleighRitz method with respect to the principle of stationary total potential energy along the beam axis.
Following the above mentioned steps, several numerical examples are performed to measure the accuracy and validity of the proposed procedure. In the case of homogenous beam with varying crosssection, the calculated buckling loads is validated by comparing with finite element simulations using ANSYS software and with the accessible results from the available numerical and analytical benchmarks. It should be pointed out that material properties of functionally graded tapered beam are supposed to vary through the longitudinal direction of the member according to simple powerlaw distribution (PFGM). Comments and conclusions are presented towards the end of the manuscript.
2 MOTION EQUATION FOR NONPRISMATIC BEAM RESTING ON AN ELASTIC FOUNDATION
An axially functionally graded beam with variable crosssection resting on WinklerPasternak foundation as depicted in Figure 1 is taken into account. It is worth noting that the EulerBernoulli beam hypothesis is considered, which is only applicable for long and slender beams. Based on this theory, the influences of shear deformation and rotatory inertia are negligible and only the effect of flexural deformation is taken into account. In the absence of damping, the differential equation of motion for nonuniform beam of length L is expressed as follows:
In the last formulation, E and ρ denote respectively Young's modulus and the material density which can be both arbitrary over the beam’s length (xaxis). w is the vertical displacement. k _{ w } and k _{ G } express Winkler elastic foundation constant and the second foundation parameter modulus in vertical direction, respectively. I and ω express the second moment of inertia and the natural frequency (circular), respectively.
For axially functionally graded nonprismatic beam lying on uniform twoparameter elastic foundation, the free vibration analysis becomes more complex due to the presence of variable coefficients in the governing fourthorder differential equation (Eq.(1)). Regarding this, the solution of the differential equation is not straightforward and only numerical procedures such as Galekin’s method, differential transformation method, finite difference method (FDM), RayleighRitz method and differential quadrature method (DQM), are possible. In this article, the power series method (PSM) is applied to solve Eq. (1). Based on this approach, all the variable terms in the differential equation including crosssectional area, moment of inertia, Young’s modulus of elasticity, and density of material should be presented in power series form as follows:
In order to facilitate the solution of the stability equation, a dimensionless variable (
Substituting equations (3) and the nondimensional variable
The vertical displacement
Then, one obtains:
Inserting Eqs. (5)(6) into equation (4), the following relationship can be derived:
in which:
After some simplifications, the following expression is obtained:
As mentioned before, the deformation shape of the beam with nonuniform crosssection was considered as follows:
Referring to the obtained recursive relationship (Eq. (9)) and using symbolic MATLAB software, the first few terms of the transverse displacement function (Eq. (10)) is derived in the following form:
According to the last expression and from mathematical point of view, it is concluded that all the b
_{
i
} coefficients can be obtained except for the first four (
where
In equation (13),
3 BOUNDARY CONDITIONS
It has to be notified that the four undetermined coefficients (
1) At
and
2) At
Then, one gets:
Referring to Figure 2 and using nondimensional coordinate, the natural end conditions at the right end (
3) At
and
4) At
Then, we have:
According to the boundary conditions presented above, a general system is derived for four linear equations with four undefined constants (
in which
where {Δ} represents the nodal displacement vector for the EulerBernoulli element. Note that all of the terms presented in the coefficient matrix ([V]) are real and [V] is an invertible matrix. Thus, the solution of previous simultaneous equation is derived as follows:
Regarding Eq. (11), it can be pointed out that
The four independent constants (
where
The mathematical expression describing the deformation shape of the axially functionally graded nonprismatic beam resting on two parameter elastic foundation under free vibration analysis can be obtained at any places in the member by replacing the obtained matrix {B} into Eq. (12):
In this study, the free vibration and liner buckling analyses are conducted for FGM nonprismatic beam having three different end supports including a beam with left end fixed and the other end free (fixedfree), a beam clamped at left end and hinged at right one (clampedhinged), and a beam hinged at both ends (hingedhinged). Accordingly, the vibration mode shape for each end conditions considered in the present study is determined in the followings.
In the case of cantilever members, vertical deflection and bending rotation at the left support (fixed one) are null. Further, the bending moment and the shear force at the free end equal to zero. Thus, the deformation shape of a clampedfree beam can be obtained in nondimensional coordinate as:
The explicit expression of vibrational shape for fixedhinged members is determined by preventing the transverse displacement and rotation at the fixed end (
For simply supported beams, the vertical displacement and bending moment at the both ends are equal to zero. Therefore, the vibrational mode shape becomes:
By observing the last three formulations and based on mathematical point of view, it is finally concluded that the approximate expression for the vertical displacement in the nondimensional coordinate can be derived in terms of the rotation of right end of the member (
According to the eigenvalue procedure, the approximate function of deformation shape of an elastic member not the exact magnitude of the amplitude of the deflection can be obtained. In the following, the RayleighRitz method is employed to obtain buckling load of nonprismatic columns on continuous linear elastic foundation. In order to estimate the linear buckling load using this technique which is based on variational calculus, an appropriate deformation shape of the element after flexural buckling satisfying both geometrical and natural boundary conditions of the system is essential. The expressions of deformation shapes are similar to both linear stability and free vibration analyses of an elastic member except for magnitude and amplitude. These constants are insignificant for evaluating the value of buckling load based on the RayleighRitz method. Therefore, the obtained vibrational mode (Eqs. (2325)) can be adopted as the deformed shape of column for the linear stability analysis.
4 THE RAYLEIGHRITZ METHOD
In this section, the value of critical buckling load of the beam with variable crosssection is determined by the RayleighRitz method based on the variation of total potential energy, which is:
In Eq. (27), δ represents a virtual variation. Ul and U0 denote the elastic strain energy and the strain energy due to the initial stresses on the considered element, respectively, and Uf is the energy related to continuous elastic foundation. The potential energy for nonuniform member in local coordinate can be written as:
It should be noted that the influences of axial stiffness (EA) and longitudinal displacement are not considered in the last expression, since these terms have no incidence on flexural buckling analysis of EulerBernoulli beam. In the current study, it is assumed that the displacement equation obtained in the previous section (Eqs. (2325)) is the approximate form of buckled shape of elastic beam. Thus, Eq. (27) can be replaced by:
Since
The smallest real root to the obtained algebraic equation (Eq. (30)) is considered as critical buckling load. It should be noted that all calculation procedures are performed with the aid of MATLAB software in order to compute the approximation displacement function based on power series expansions and the critical buckling load with respect to energy method.
5 NUMERICAL EXAMPLES
In order to clarify and demonstrate the performance and efficiency of proposed procedure in free vibration and stability analyses of AFG beams with variable crosssections whether resting on elastic foundation or not, several examples are conducted.
At the first, in the absence of elastic foundation and nonhomogeneous material, the validity of the formulation for analyzing the static buckling and free vibration is accomplished by comparing the results to those obtained by Ansys software and to other analytical and numerical solutions presented in the literature. This is due to the lack of similar study for analyzing axially functionally graded nonprismatic beams on continuous elastic foundation. The consistency in the results indicates the accuracy of the presented numerical method. Then, the static and dynamic analyses of nonprismatic EulerBernoulli beams in the presence of one or twoparameter elastic foundation and axially functionally graded materials are investigated.
To the best of our knowledge on Power Series Method (PSM) (^{Asgarian et al. (2013}), ^{Soltani et al. (2014}) and Soltani and Asgarian (2019)), the results of this semianalytical approach are extremely sensitive to the number of terms considered in the power series approximations. Therefore, it is important to estimate appropriate number of terms (n) in power series expansion in each case in order to obtain an explicit expression for deformation shape of the member and accordingly calculate the lowest vibration frequencies (circular) with the desired precision, as well as optimizing the symbolic mathematical procedures. One can check that the vibration modes and their relating frequencies ( ω ) are highly dependent on the number of terms in power series expansion. The values of free vibration frequencies of higher modes approach the exact solutions as the number of power series expansion increases. This can be observed in Soltani et al. (2014). The critical buckling loads are obtained by using the RayleighRitz method after determining the accurate expression for vibration mode shape of the member. It should be emphasized that the calculated natural frequencies and axial critical loads are considered as real values.
5.1 Example 1: Nonprismatic homogenous beam
The purpose of this example is to comprehensively study the exactness and performance of the proposed hybrid numerical method in calculating the critical buckling loads and natural frequencies of homogenous beams with variable crosssection. To this aim, as shown in Figure 3, three nonuniform members are investigated with different boundary conditions (fixedfree, hingedhinged and fixedhinged). In all of the presented cases, the crosssection is in the form of rectangle with width b
_{
0
} and depth h
_{
0
}
, which is assumed to be sufficiently small compared to the width. It is also assumed that the geometrical properties of the left end section of the member are constant. Each member is subjected to a concentrated compressive axial load. In the following, the geometrical parameters at the left support (x=0,
In the first case (Case A), the width of the beam’s section at the left support (b0) is linearly diminished to (
In the second case (Case B), the height of the beam’s section is (h0) at the left support and linearly decreased to (
In Case C, the height and breadth of the beam are concurrently allowed to vary linearly along the member’s length with the same tapering ratio of
In this example, the tapering parameter can change from
In which A0 and I0 are respectively crosssectional area and moment of inertia at the left support (x=0,
In order to simplify the illustration of the obtained results, the dimensionless natural frequency and buckling load parameters are used as follows:
where the mechanical properties at the left support (x=0) of the beam are indicated with the subscript 0. The first section of this example aims to define the required number of terms in power series expansions to obtain an acceptable accuracy on natural frequencies and critical buckling loads. Regarding this, the lowest value of normalized circular frequency (ω_{nor}) and nondimensional buckling load parameter (P _{ nor } ) of the tapered beam based on EulerBernoulli beam assumptions for two cases ( B and C ) and two different boundary conditions (pinnedended and fixedfree) are calculated with respect to the number of power series terms (n). Note that the convergence study is carried out for three taper ratios: β =0, 0.2, and 0.6. Table 1summarizes the dimensionless natural frequencies and critical loads for Case A. In addition, the effect of the number of power series terms (n) on convergence is displayed in Table 1. Similarly, the result of double tapered homogenous EulerBernoulli beam for Case C have been presented in Table 2. Besides, the elapsed time to perform numerical computations is displayed in this table. The obtained results by the proposed numerical technique have been compared with those of finite element method by using ANSYS software. The nonprismatic and prismatic members has been modeled using BEAM188 in ANSYS software. This member is a twonode element with six degrees of freedom three translational UX, UY, UZ and three rotary ROTX, ROTY, ROTZ of freedom at each node.
According to Tables 1 and 2, it can be easily observed that an increase in the number of polynomials has a great effect on the convergence rate of fundamental frequency and critical load parameters at each case. Table 2 shows that Central Processing Unit (CPU) needs an average of 8.85 seconds to accomplish power series method simulation for a pinnedpinned doubletapered beam (α=0.6) with the number terms in the series equals to 20. The normalized fundamental vibration frequency is 6.2346 (error Δ=0.64%). When the number of terms is increased to 40, the nondimensional vibration frequency is close to beam results of ANSYS software (ω_{nor}=6.2086, error Δ<0.5%). The needed CPU time is only 133.62s.
Case  End Conditions  β  Type of Analysis  Number of terms of power series (PSM)  ANSYS  
n=10  n=20  n=30  n=40  n=50  
B  HingedHinged  0  (ω_{nor})  9.6730  9.8696  9.8696  9.8696  9.8696  9.7980 
(P_{nor})  9.8820  9.8696  9.8696  9.8696  9.8696  9.8700  
0.2  (ω_{nor})  8.8944  8.8462  8.8462  8.8462  8.8462  8.8172  
(P_{nor})  7.0133  7.0456  7.0456  7.0456  7.0456  7.0556  
0.6  (ω_{nor})  6.8076  6.4754  6.4667  6.4666  6.4666  6.4542  
(P_{nor})  2.5324  2.6774  2.6789  2.6789  2.6789  2.6638  
ClampedFree  0  (ω_{nor})  3.5073  3.5160  3.5160  3.5160  3.5160  3.5040  
(P_{nor})  2.4669  2.4674  2.4674  2.4674  2.4674  2.4670  
0.2  (ω_{nor})  3.6273  3.6083  3.6083  3.6083  3.6083  3.6027  
(P_{nor})  2.0241  2.0233  2.0233  2.0233  2.0233  2.0203  
0.6  (ω_{nor})  4.7193  3.9660  3.9351  3.9343  3.9343  3.9557  
(P_{nor})  0.7016  1.1273  1.0999  1.0999  1.0999  1.0970 
Case  End Conditions  β  Type of Analysis  Number of terms of power series (PSM)  ANSYS  
n=10  n=20  n=30  n=40  n=50  
C  HingedHinged  0.2  Vibration  8,8933  8.8246  8.8246  8.8246  8.8246  8.7952 
(ω_{ nor } )  (3.14s)  (7.98s)  (17.92s)  (125.82s)  (3088.7s)  
Stability  6.2364  6.2668  6.2668  6.2668  6.2668  6.2846  
(P_{nor})  (4.88s)  (11.95s)  (25.81s)  (140.31s)  (3314.6s)  
0.6  Vibration  6.8936  6.2346  6.2092  6.2086  6.2085  6.1952  
(ω_{nor})  (2.89s)  (8.85s)  (22.64s)  (133.62s)  (3522.6s)  
Stability  1.4466  1.5933  1.5958  1.5958  1.5958  1.6006  
(P_{nor})  (4.51s)  (12.75s)  (31.28s)  (148.12s)  (3607.1s)  
ClampedFree  0.2  Vibration  3.8802  3.8551  3.8551  3.8551  3.8551  3.8166  
(ω_{nor})  (5.02s)  (11.56s)  (25.54s)  (142.79s)  (3504.8s)  
Stability  1.9164  1.8835  1.8835  1.8835  1.8835  1.8821  
(P_{nor})  (8.38s)  (17.07s)  (37.72s)  (160.11s)  (3690.4s)  
0.6  Vibration  6.1174  5.0720  5.0113  5.009  5.0090  5.0803  
(ω_{nor})  (5.04s)  (11.84s)  (25.80s)  (147.32s)  (3628.7s)  
Stability  0.8758  0.8163  0.7658  0.7624  0.7624  0.7605  
(P_{nor})  (7.11s)  (17.75s)  (34.99s)  (164.18s)  (3819.6s) 
For stability analysis, an accurate value of normalized buckling load parameter equal to 1.5933 (error Δ=0.46%) is obtained with a number terms in the series equals to 20. This method requires 12.75s CPU time. In order to improve the results, more CPU times are needed. With 40 terms the buckling parameter is 1.5958 (error Δ=0.3%) the CPU time is impressive (148.12 s). Note that when the number terms is increased over 40, the natural frequency and the buckling load become insensitive and the CPU time increases accordingly. The same statement is true for other cases illustrated in Table 2.
Finally, it can be concluded that for a high accurate solution, it is not required to take more than 30 terms of power series in the current technique. It is worth mentioning that 40 terms are necessary to calculate the value of vibration frequencies and critical buckling loads with an acceptable error rate for cases with taper ratios higher than 0.7. It is also observed that the proposed hybrid method is considerably capable of predicting the exact buckling loads and natural frequencies simultaneously.
In the following, the nondimensional free transverse frequencies (ω_{nor}) and dimensionless critical loads (P _{ nor } ), which play a significant role in designing and analyzing different structures, corresponding to case A to C, have been estimated by adopting 40 terms in the power series expansions. In the case of vibration analysis, the verifying computations are performed using Beam188 of ANSYS software. Note that, the obtained buckling loads are compared with the closedform solution introduced by ^{Wang et al. (2005}).
Table 3 consists of the results of free vibration analysis, while Table 4 summarizes the results of instability analysis. Furthermore, both tables include the relative errors (Δ), which are obtained by using:
As shown in Tables 3 and 4, the values of buckling loads and natural frequencies for nonuniform members derived from this method are respectively in excellent agreement with the results obtained by finite element method using ANSYS software and the exact solution presented in ^{Wang et al. (2005}).
Case  β  Free Vibration Analysis  
HingedHinged  ClampedFree  ClampedHinged  
ANSYS  Present Method 

ANSYS  Present Method 

ANSYS  Present Method 


A  0.1  9.792  9.868  0.776  3.618  3.631  0.359  15.540  15.527  0.084 
0.3  9.785  9.860  0.766  3.901  3.916  0.384  15.784  15.768  0.101  
0.5  9.752  9.825  0.748  4.297  4.315  0.418  16.065  16.044  0.131  
0.7  9.748  9.747  0.010  4.949  4.932  0.343  16.385  16.351  0.208  
0.9  9.545  9.563  0.188  6.091  6.084  0.115  16.720  16.724  0.024  
B  0.1  9.306  9.368  0.666  3.547  3.559  0.338  14.862  14.849  0.087 
0.3  8.253  8.302  0.594  3.655  3.667  0.328  13.658  13.640  0.132  
0.5  7.093  7.122  0.409  3.813  3.824  0.288  12.328  12.300  0.227  
0.7  5.736  5.745  0.157  4.134  4.082  1.258  10.788  10.737  0.473  
0.9  3.885  3.976  2.342  4.713  4.754  0.870  8.772  8.879  1.220  
C  0.1  9.301  9.362  0.656  3.661  3.674  0.355  14.969  14.955  0.094 
0.3  8.206  8.250  0.536  4.050  4.067  0.419  13.982  13.962  0.143  
0.5  6.926  6.957  0.448  4.610  4.625  0.325  12.883  12.850  0.256  
0.7  5.348  5.359  0.206  5.499  5.512  0.236  11.620  11.577  0.370  
0.9  3.045  3.063  0.591  7.267  7.300  0.454  10.097  10.159  0.617 
Case  β  Stability Analysis  
HingedHinged  ClampedFree  ClampedHinged  
Wang et al. (2005)  Present Method 

Wang et al. (2005)  Present Method 

Wang et al. (2005)  Present Method 


A  0.1  9.372  9.372  0.000  2.393  2.393  0.000  19.170  19.168  0.010 
0.3  8.343  8.343  0.000  2.235  2.235  0.000  17.030  17.035  0.029  
0.5  7.256  7.256  0.000  2.062  2.062  0.000  14.740  14.739  0.007  
0.7  6.069  6.069  0.000  1.865  1.879  0.751  12.180  12.177  0.025  
0.9  4.667  4.672  0.107  1.621  1.631  0.617  9.029  9.083  0.598  
B  0.1  8.436  8.434  0.024  2.246  2.246  0.000  17.250  17.252  0.012 
0.3  5.840  5.840  0.000  1.798  1.798  0.000  11.920  11.923  0.025  
0.5  3.628  3.628  0.000  1.336  1.336  0.000  7.362  7.362  0.000  
0.7  1.821  1.806  0.824  0.853  0.861  0.938  3.634  3.604  0.834  
0.9  0.471  0.474  0.637  0.321  0.325  1.246  0.875  0.881  0.686  
C  0.1  7.994  7.994  0.000  2.175  2.175  0.000  16.350  16.354  0.024 
0.3  4.836  4.836  0.000  1.595  1.595  0.000  9.893  9.893  0.000  
0.5  2.467  2.470  0.122  1.029  1.029  0.000  5.048  5.058  0.198  
0.7  0.888  0.894  0.676  0.498  0.502  0.803  1.817  1.832  0.825  
0.9  0.099  0.100  0.909  0.080  0.081  1.250  0.202  0.205  1.485 
Although in this study, the buckling loads are obtained on the basis of free vibration analysis, it is observed that even with 40 terms in the power series expansion, buckling loads determined by the present hybrid numerical method are competitive with those of other benchmarks (see Tables 34). Small difference in the results is observed for higher values of tapering ratio.
Further, the effect of tapering ratio on the fundamental frequency analysis of homogenous beam with varying crosssection is sensitive to the boundary conditions. Regarding the pinnedpinned beams, it is observed that natural frequencies decrease by increasing the taper parameter. However in the case of cantilever members, it is obviously seen that the fundamental frequency is increasing with an increase in nonuniformity parameter. In this regard, it can be stated that the free vibration behavior of nonuniform is not predictable.
Comparing the buckling results of prismatic beam tabulated in Table 1 with those of nonprismatic ones (Table 4), it can also be stated that the considered prismatic beam in this example has the highest critical buckling loads. In this regard, let us consider that the tapering parameter equals to 0.5, the differences between the values of buckling loads for cantilever nonprismatic beams of cases A, B and C are 16%, 46% and 58% of fixedfree prismatic beam, respectively. This can be explained by the fact that increasing in width and/or height taper ratios causes the reduction in crosssectional area and moment of inertia, and consequently the stiffness of the elastic member. Regarding Eq. (31), the power exponent in the expression of minor axis moment of inertia for the last formulation (Eq. (31a)) related to coincident variation of the height and width of the assumed crosssection equals to 4, while the powerindex for other expressions corresponding to the moment of inertia is 1 or 3. Therefore, the reduction of elastic stability capacity of nonprismatic member caused by increasing the width and height nonuniformity ratios (
5.2 Example 2 Nonprismatic member with AFG material
In this section, to study the effects of material nonhomogeneity on linear stability and free vibration, a nonprismatic member made up of axially functionally graded material is considered. Axially functionally graded EulerBernoulli beam is assumed with exponentially tapered crosssection which is shown in Figure 4. Therefore, the geometrical properties including the crosssection area and the moment of inertia vary along the beam axis with exponential functions. Thus, the variation of minor axis moment of inertia and crosssectional area in the local coordinate are described as follows:
The nonuniformity parameter (α) can change from zero (prismatic beam) to a range of [2 to 0.1] for nonuniform beams. Exponential variation of geometrical and/or material properties is regarded as one of the most special states of members that few numerical methods such as differential transformation method, generalized differential quadrature, method finite element method, and power series approach can solve its governing differential equation (^{Mohanty et al, 2012}; ^{Li et al., 2013}, 2018; ^{Ebrahimi and Mokhtari, 2015}; ^{Wang et al., 2016}; ^{Khaniki and Rajasekaran, 2018}; ^{Arefi et al., 2018}; ^{Soltani and Asgarian, 2019}). Regarding the power series method for solving the linear differential equation, all variable coefficients should be represented in a polynomial form. In the presence of exponential variation, the Maclaurin series expansion should be adopted. Accordingly, the explicit forms of Maclaurin series for crosssectional area and moment of inertia expressions (Eq. (35)) are defined as:
In this example, it is supposed that the beam is made of two different materials specifically zirconia ( ZrO2 ) and aluminum ( Al ) in the length direction with the following characteristics:
ZrO2: E_{0}=200GPa, ρ_{0} =5700 kg/m3;
Al: E_{1}=70GPa, ρ_{1} =2702 kg/m3;
The variation of Young’s modulus of elasticity and the material density along the beam axis are defined with the following powerlaw formulations:
In the last expression, m is the nonhomogeneity parameter of the material. Based on this formulation, it can be stated that the proportion of zirconia over the beam’s length increases by raising the powerlaw index ( m ). It should be noted that the material properties of the beam are assumed to be constant in the direction of the thickness. In addition, the Poisson’s ratio of the material remains constant in longitudinal direction. For convenience, m is assumed to be a positive integer. Thus, PSM could be easily applied to the problem.
Regarding the complicated process of calculating the buckling loads and free transverse frequencies of nonprismatic members with exponential variation of geometrical properties in the first stage of the current example, the first three nondimensional natural frequencies of simply supported homogenous beam for two different nonuniformity parameters (α), namely 0 and 1 are estimated according to the number of the terms of power series needed for convergence which are reported in Table 5. Then, the results are compared to the exact solution proposed by ^{Ece et al. (2007}).
α  Vibration Mode  Number of terms of power series (PSM)  Exact Solution  

n=10  n=20  n=30  n=40  n=50  Ece et al. (2007)  
0  1  9.673  9.8696  9.8696  9.8696  9.8696  9.8696 
2    39.41691  39.4784  39.4784  39.4784  39.4784  
3    149.0075  88.81072  88.8264  88,8264  88.8266  
1  1  9.65319  9.7729  9.7729  9.7729  9.7729  9.7729 
2    39.6464  39.5704  39.5704  39.5704  39.5704  
3    77.7059  88.9738  88.9705  88.9705  88.9705 
One can observe from Table 5 that there is an acceptable accuracy on free frequencies of the first two vibration modes when the taken number of terms of the power series exceeds over 20. It should be emphasized that 30 terms in the power series expansion is a good compromise for equivalent accuracy of the first three natural frequencies.
In the next stage of this numerical example, linear stability analysis for homogenous nonprismatic beam with different nonuniformity ratios (α) should be considered to estimate the required numbers of the terms in power series. Table 6 indicates the dimensionless parameter of critical load (
End Conditions  α  Number of terms of power series (PSM)  Wang et al. (2005) 


n=10  n=20  n=30  n=40  n=50  
HingedHinged  0.1  9.0847  9.3857  9.3857  9.3857  9.3857  9.380  0.061 
0.5  8.1909  7.6345  7.6345  7.6345  7.6345  7.634  0.006  
1.0  5.7210  5.8266  5.8265  5.8265  5.8265  5.827  0.009  
1.5  3.7071  4.3873  4.3885  4.3885  4.3885  4.389  0.011  
2.0  2.6917  3.2624  3.2636  3.2636  3.2636  3.264  0.012  
Clamped Free  0.1  2.3970  2.3945  2.3945  2.3945  2.3945  2.394  0.021 
0.5  2.0973  2.1121  2.1121  2.1121  2.1121  2.110  0.100  
1.0  1.7976  1.7821  1.7821  1.7821  1.7821  1.782  0.006  
1.5  1.7190  1.4804  1.4803  1.4803  1.4803  1.480  0.020  
2.0  1.8681  1.2092  1.2093  1.2093  1.2093  1.209  0.025  
Clamped Hinged  0.1  17.8937  19.2016  19.2018  19.2018  19.2018  19.200  0.009 
0.5  14.2300  15.6388  15.6399  15.6399  15.6399  15.640  0.001  
1.0  8.8990  11.9984  11.9884  11.9884  11.9884  11.990  0.013  
1.5  5.7520  9.0437  9.0980  9.0980  9.0980  9.098  0.000  
2.0  4.9715  6.7942  6.8397  6.8397  6.8397  6.839  0.010 
As can be seen in Table 6, the satisfactory results for engineering requirements can be reached by considering 30 terms of power series. Based on the definitions of the proposed numerical method presented earlier, it is observed that using 30 terms of power series is a good compromise to obtain a satisfactory accuracy on critical elastic buckling loads and natural frequencies. In the following sections, the number of the terms in power series (the maximum power of shape functions) is taken as N=30. Therefore, the relative error becomes sufficiently small.
In what follows, the influences of nonuniformity parameter (α) on dimensionless natural frequencies and critical loads of AFG nonprismatic beam by considering the nonhomogeneity indexes m=1, 2 and 3 for simply supported, fixedfree, and fixedpinned are respectively illustrated in Figures 5 and 6.
Regarding Eq. (37), it is obvious that by increasing the gradient index ( m ) from 1 to 3, the volume fraction of Zirconia is increased and as a result, a stiffer and also heavier beam is acquired. In addition, a decrease in nonuniformity parameter from 0.0 to 2 leads to a reduction in the moment of inertia and crosssectional area and accordingly a gradual decrease in both stiffness and mass of beam. Therefore, the beam easily becomes weaker and lighter. Since the free transverse vibration in the absence of damping relies on both the mentioned terms (mass and stiffness) predicting the variations of natural frequency with respect to nonhomogeneity index (m) and tapering parameter (α) concurrently is impossible. This fact can be easily observed in Figure 5.
Comparing the results of Figure 6, it can be concluded that when nonhomogeneity index is increased, the value of Young’s modulus increases, and as a result, a higher buckling load is achieved. For example in the case of simply supported AFG prismatic beam (α=0), the normalized critical buckling load increases from 6.379 to 7.785 and then to 8.497, when m increases from 1 to 3. It shows an increase by 22.04% and 33.2%, accordingly. As illustrated in Figure 6, regarding a constant powerlaw exponent the buckling resistance of uniform beam (α=0) and nonprismatic beam with α=2 are considered as the most and least amount, respectively, since a decrease in nonuniformity parameter leads to a reduction in crosssectional area and moment of inertia, and consequently in the stiffness of the elastic member.
5.3 Example 3 AFG Nonprismatic beam on continuous WinklerPasternak foundation
To study the effects of WinklerPasternak foundation on linear stability and free vibration, another example is considered including axially nonhomogeneous and homogeneous beams. In this regard, an elastic EulerBernoulli beam with rectangular crosssection whose width and thickness taper linearly along the member axis with similar taper ratio is taken into account. Therefore, the crosssectional area and moment of inertia are identical to the case C in the first example (Eq. (31c)). The convergence study of the results is carried out for two taper ratios (
Accordingly, estimating the required number of the terms in power series to obtain the natural frequency and the lowest buckling load with an excellent accuracy is important since there are no accessible results for linear stability analysis on AFGM tapered beams on a twoparameter foundation. Therefore, a converge study is conducted for linear stability and free vibration characteristics of the considered hingedhinged member with different foundation parameters:
Next, the variation of the frequencies of vibration of hingedhinged and fixedhinged EulerBernoulli beams resting on Winkler type of foundation versus the elastic foundation constant (
Finally, in order to investigate the influence of uniform WinklerPasternak foundation on the behavior of simply supported EulerBernoulli beam, the first transverse natural frequencies and the lowest normalized critical buckling loads for different combinations of elastic foundation parameters i.e.
Material  (β)  The Pasternak foundation constant (

Nondimensional bending frequency  Nondimensional critical buckling load  

The Winkler foundation constant (

The Winkler foundation constant (


0  20  40  60  80  0  20  40  60  80  
Homogeneous  0  0  9.8696  10.8340  11.7207  12.5449  13.3182  9.8694  11.8965  13.9236  15,9507  17,9778 
2  10.8217  11.7093  12.5343  13.3082  14.0395  11.8694  13.8965  15.9236  17,9507  19,9778  
4  11.6979  12.5236  13.2982  14.0300  14.7255  13.8694  15.8965  17.9236  19,9507  21,9778  
6  12.5130  13.2881  14.0205  14.7165  15.3810  15.8694  17.8965  19.9236  21,9507  23,9778  
8  13.2781  14.0109  14.7074  15.3723  16.0096  17.8694  19.8965  21.9236  23,9507  25,9778  
0.5  0  6.9566  9.2419  11.0495  12.5923  13.9584  2.4891  4.3334  6.0834  7,6699  9,0018  
2  9.2683  11.0618  12.5956  13.9558  15.1891  4.4891  6.3334  8.0834  9,6699  11,0018  
4  11.0660  12.5938  13.9499  15.1805  16.3142  6.4891  8.3334  10.0834  11,6699  13,0018  
6  12.5885  13.9418  15.1706  16.3031  17.3582  8.4891  10.3334  12.0834  13,6699  15,0018  
8  13.9322  15.1597  16.2915  17.3462  18.3373  10.4891  12.3334  14.0834  15,6699  17,0018  
Axially Functionally Graded m=2  0  0  9.5994  10.7618  11.8099  12.7718  13.6657  7.8308  9.8396  11.8455  13,8478  15,8453 
2  10.7502  11.7989  12.7614  13.6557  14.4947  9.8308  11.8396  13.8455  15,8478  17,8453  
4  11.7877  12.7507  13.6456  14.4849  15.2778  11.8308  13.8396  15.8455  17,8478  19,8453  
6  12.7399  13.6353  14.4750  15.2683  16.0220  13.8308  15.8396  17.8455  19,8478  21,8453  
8  13.6248  14.4650  15.2587  16.0128  16.7326  15.8308  17.8396  19.8455  21,8478  23,8453  
0.5  0  6.5127  9.3065  11.4053  13.1480  14.6634  1.7011  3.3688  4.7855  5,8456  6,5903  
2  9.3724  11.4325  13.1540  14.6577  16.0058  3.7011  5.3688  6.7855  7,8456  8,5903  
4  11.4441  13.1525  14.6487  15.9930  17.2215  5.7011  7.3688  8.7855  9,8456  10,5903  
6  13.1464  14.6377  15.9795  17.2072  18.3441  7.7011  9.3688  10.7855  11,8456  12,5903  
8  14.6252  15.9654  17.1928  18.3302  19.3938  9.7011  11.3688  12.7855  13,8456  14,5903 
6 CONCLUSIONS
In the present study, a hybrid numerical methodology to study the stability and free bending vibration in nonprismatic axially functionally graded beams based on EulerBernoulli theory is proposed. This new technique combines the power series method and RayleighRitz approach and applied to nonprismatic beams having generalized end conditions. In this regard, the power series approximation is implemented to solve the fourth order differential equation of motion since stiffness quantity is unstable in the presence of variable crosssection. Then, all variable parameters including displacement component, geometrical and material properties are developed based on power series of an identified degree. By solving the eigenvalue problem, one can acquire the natural frequencies.
Further, the explicit expression of vibrational shape function is evaluated based on the abovementioned method. The vibrational mode shapes of an elastic member are similar to the buckling ones. Therefore, the obtained deflected shapes of the considered nonprismatic beams under free transverse vibration analysis can be used as the deformation shape of the member for the linear buckling analysis. The critical buckling load of nonprismatic beam resting on an elastic foundation can be then estimated by adopting the principle of stationary total potential energy. The applicability and efficiency of mixed power series and RayleighRitz methods in free vibration and linear stability analyses of nonprismatic AFG beams on twoparameter elastic foundations with three different classical boundary conditions are demonstrated by providing several numerical examples in which the effects of uniform elastic foundation, mechanical variation and different boundary conditions are comprehensively investigated. The acquired results are contrasted with other analytical and numerical solutions presented in the literature and with finite element method by using ANSYS software. In all cases, the fast rate of convergence is demonstrated and the results are in consistent with those of other techniques even by considering 40 terms of power series. The impact of width and/or height tapering ratios, axial variation of material properties, boundary conditions and WinklerPasternak parameter on linear buckling resistance and the natural frequency of EulerBernoulli beam with varying crosssection are comprehensively surveyed. It should be pointed out that material properties of functionally graded beam with variable crosssection are supposed to vary through longitudinal direction of the constituents according to simple powerlaw distribution (PFGM). From the numerical examples, it can be concluded that the effects of width and/or height tapering ratios play important roles on the linear stability capacity and the fundamental frequency of AFG tapered beam lying on elastic foundation. It is also observed that, for ceramicmetal beams, buckling resistance increases with increase of functionally graded material content for powerlaw property distribution. In addition, the numerical outcomes reveal that the elastic foundation enhance the stability characteristics of axially nonhomogeneous and homogeneous beams with constant or variable crosssection.