New hybrid approach for free vibration and stability analyses of axially functionally graded Euler-Bernoulli beams with variable cross-section resting on uniform Winkler-Pasternak foundation

Abstract

In this paper, a new hybrid approach is presented based on the combination of the power series expansions and the Rayleigh-Ritz method for stability and free vibration analyses of axially functionally graded non-uniform beams resting on constant Winkler-Pasternak elastic foundation. In the proposed novel technique, the power series approximation is first adopted to solve the motion equation. Regarding this numerical methodology, the transverse displacement and all mechanical properties are expanded in terms of power series of a known degree. By solving the eigenvalue problem, one can acquire the fundamental natural frequencies. According to aforementioned method, the expression of vibrational mode shape is also determined. Based on the similarities existing between the vibrational and buckling deformation shapes, Rayleigh-Ritz method is finally employed to construct eigenvalue problem for obtaining the critical loads. In order to illustrate the correctness and convergence of the method, several numerical examples of axially non-homogeneous and homogeneous beams are conducted. The obtained outcomes are compared to the results of Finite Element Analysis in terms of ANSYS software and those of other available numerical and analytical solutions. The accuracy of the method is then remarked.

Keywords:
Non-prismatic beams; Elastic foundation; Power series expansions; Rayleigh-Ritz method; Buckling load; natural frequency

1 INTRODUCTION

Elastic flexural members whose cross-sectional profile changes partially or gradually along their length, known as non-prismatic 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, non-prismatic 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 non-prismatic 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 cross-section. Regarding axially graded Euler-Bernoulli beams with an arbitrarily variable cross-section, 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 closed-form solutions for exact estimation of buckling loads and natural frequencies of non-prismatic beams made up of axially non-homogeneous 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 non-uniform 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 non-uniform 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 step-varying profiles. Sampaio and Hundhausen (1998) presented a mathematical model and analytical solution for buckling of inclined beam-columns. Only one type of boundary condition (pinned-ends) 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 cross-sections. In another study, Atay (2009) applied a new approach, known as Homotopy perturbation method, to solve the buckling problem of non-prismatic columns. Further, Okay et al. (2010) derived buckling loads and mode shapes of a heavy column by applying the variational iteration method.

Regarding the above-mentioned studies, all are merely related to the buckling and dynamic analyses of homogeneous Euler-Bernoulli 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 step-reduction method, Tong and Tabarrok (1995) developed a new approach for the free and forced vibration of non-homogeneous Timoshenko beams with non-uniform cross-sections. In addition, an elasticity solution based on functionally graded Euler-Bernoulli 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 first-order shear deformation theory to investigate the thermo-elastic behavior of FG beam structures. Wu et al. (2005) adopted a semi-inverse 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 non-uniform beams by using a newly developed numerical method. Based on a new beam theory different from the traditional first-order 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 non-homogeneous tapered Euler-Bernoulli beams were performed by Shahba et al. (2013a, 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 non-linear thermo-elastic analysis of a functionally graded piezoelectric (FGP) thick-walled cylinder under thermal, mechanical and electrical loads. Pradhan and Chakraverty (2013) used Rayleigh-Ritz 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 zig-zag theory (ZIGT). Based on the modified couple stress theory, mechanical behavior of non-uniform bi-directional functionally graded beam sensors was studied in detail by Khaniki and Rajasekaran (2018). Li et al. (2018) conducted instability analysis of a micro-scaled bi-directional functionally graded (FG) beam having rectangular cross-section by employing the Generalized Differential Quadrature Method (GDQM). Recently, Soltani and Asgarian (2019) performed the stability analysis of cantilever axially functionally graded non-prismatic 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 soil-structure 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 two-parameter elastic foundation models such as Winkler-Pasternak 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 thin-walled beam with non-symmetric cross-section resting on two-parameter elastic foundation by considering shear deformation. Static and dynamic stiffness matrices of members with variable cross-section resting on Winkler-Pasternak elastic foundations were derived by Girgin and Girgin (2005). This numerical approach was developed based on the well-known Mohr method. Further, Ruta (2006) used Chebyshev polynomial method to solve the coupled system of equilibrium differential equations of Timoshenko beams with non-uniform cross-section 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 non-linear free and forced vibration analyses of non-prismatic 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 post-buckling analysis of the functionally graded beams under an axial load and resting on nonlinear elastic foundation. Mirzabeigy (2014) presented a semi-analytical technique based on the differential transformation method to obtain the dimensionless natural frequencies of non-uniform beams resting on an elastic foundation. Furthermore, Tsiatas (2014) suggested a new influential approach to exactly determine stiffness and mass matrices of non-uniform Euler-Bernoulli 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 semi-rigid end conditions, buckling analysis of axially functionally graded Euler-Bernoulli beam having non-uniform cross-section 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 Euler-Bernoulli 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 non-prismatic beam are assumed to be graded smoothly along the beam axis by a power-law 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 linear-elastic.

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 Rayleigh-Ritz method. It is worth noting that Euler-Bernoulli beam theory is adopted in this study in order to analyze the beams with uniform or non-uniform cross-section, 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 non-prismatic member with variable coefficients. The variable terms in equation of motion including cross-sectional 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 cross-section is evaluated based on the similarities between vibration and buckled deformation shapes of elastic members and adopting the Rayleigh-Ritz 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 cross-section, 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 power-law distribution (P-FGM). Comments and conclusions are presented towards the end of the manuscript.

2 MOTION EQUATION FOR NON-PRISMATIC BEAM RESTING ON AN ELASTIC FOUNDATION

An axially functionally graded beam with variable cross-section resting on Winkler-Pasternak foundation as depicted in Figure 1 is taken into account. It is worth noting that the Euler-Bernoulli 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 non-uniform beam of length L is expressed as follows:

$\frac{d^2}{d{x}^2}\left[E(x)I(x)\frac{d^{2}w}{d{x}^2}\right]-k_G\frac{d^{2}w}{d{x}^2}+(k_w-\rho (x)A(x){\omega }^2)w(x)=0$ (1)

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 (x-axis). 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.

Figure 1:
Non-prismatic beam resting on a two-parameter elastic foundation

For axially functionally graded non-prismatic beam lying on uniform two-parameter elastic foundation, the free vibration analysis becomes more complex due to the presence of variable coefficients in the governing fourth-order 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), Rayleigh-Ritz 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 cross-sectional area, moment of inertia, Young’s modulus of elasticity, and density of material should be presented in power series form as follows:

$I(x)={\displaystyle \sum _{i=0}^{\infty }{I}_i{x}^i}\text{,}A(x)={\displaystyle \sum _{i=0}^{\infty }{A}_i{x}^i}\text{,}E(x)={\displaystyle \sum _{i=0}^{\infty }{E}_i{x}^i}\text{,}\rho (x)={\displaystyle \sum _{i=0}^{\infty }{\rho }_i{x}^i}$ (2)

In order to facilitate the solution of the stability equation, a dimensionless variable ($\xi =x/L$) is introduced. All abovementioned series can then be written in terms of $\xi$ as:

$I(\xi )={\displaystyle \sum _{i=0}^{\infty }{I}_i{L}^i{\xi }^i}\text{,}A(\xi )={\displaystyle \sum _{i=0}^{\infty }{A}_i{L}^i{\xi }^i}\text{,}E(\xi )={\displaystyle \sum _{i=0}^{\infty }{E}_i{L}^i{\xi }^i}\text{,}\rho (\xi )={\displaystyle \sum _{i=0}^{\infty }{\rho }_i{L}^i{\xi }^i}$ (3)

Substituting equations (3) and the non-dimensional variable $\xi$ into equation (1), the following expression is obtained:

$\frac{d^2}{d{\xi }^2}\left[({\displaystyle \sum _{i=0}^{\infty }{E}_i{L}^i{\xi }^i})({\displaystyle \sum _{j=0}^{\infty }{I}_j{L}^j{\xi }^j})\frac{d^{2}w}{d{\xi }^2}\right]-L^2{k}_G\frac{d^{2}w}{d{\xi }^2}+L^4\left[\left(k_w-{\omega }^2({\displaystyle \sum _{i=0}^{\infty }{\rho }_i{L}^i{\xi }^i})({\displaystyle \sum _{j=0}^{\infty }{A}_j{L}^j{\xi }^j})\right)w(\xi )\right]=0$ (4)

The vertical displacement $w(\xi )$ is also presented by a power series of the form:

$w(\xi )={\displaystyle \sum _{k=0}^{\infty }{b}_k{\xi }^k}$ (5)

Then, one obtains:

$\begin{array}{l}\frac{dw(\xi )}{d\xi }={\displaystyle \sum _{k=1}^{\infty }k{b}_k{\xi }^{k-1}}={\displaystyle \sum _{k=0}^{\infty }(k+1)b_{k+1}{\xi }^k}\\ \frac{d^{2}w(\xi )}{d{\xi }^2}={\displaystyle \sum _{k=2}^{\infty }k(k-1)b_k{\xi }^{k-2}}={\displaystyle \sum _{k=0}^{\infty }(k+1)(k+2)b_{k+2}{\xi }^k}\end{array}$ (6a,b)

Inserting Eqs. (5)-(6) into equation (4), the following relationship can be derived:

$\begin{array}{l}\frac{d^2}{d{\xi }^2}\left[\left\{{\displaystyle \sum _{i=0}^{\infty }{E}_i^{*}{\xi }^i}\right\}\left\{{\displaystyle \sum _{j=0}^{\infty }{I}_j^{*}{\xi }^j}\right\}\left\{{\displaystyle \sum _{k=0}^{\infty }{b}_{k+2}(k+2)(k+1){\xi }^k}\right\}\right]\\ -L^2{k}_G\left\{{\displaystyle \sum _{k=0}^{\infty }{b}_{k+2}(k+2)(k+1){\xi }^k}\right\}\\ +\left(L^4\right)\left[\left\{k_w-{\omega }^2({\displaystyle \sum _{i=0}^{\infty }{\rho }_i^{*}{\xi }^i})({\displaystyle \sum _{j=0}^{\infty }{A}_j^{*}{\xi }^j})\right\}\left\{{\displaystyle \sum _{k=0}^{\infty }{b}_k{\xi }^k}\right\}\right]=0\end{array}$ (7)

in which:

$I_i^{*}=I_i{L}^i\text{,}{A}_i^{*}=A_i{L}^i\text{,}{E}_i^{*}=E_i{L}^i\text{,}{\rho }_i^{*}={\rho }_i{L}^i$ (8)

After some simplifications, the following expression is obtained:

$\begin{array}{l}{\displaystyle \sum _{k=0}^{\infty }[{\displaystyle \sum _{j=0}^{k+2}{\displaystyle \sum _{i=0}^j{E}_i^{*}{I}_{j-i}^{*}(k-j+4)(k-j+3)(k+2)(k+1)b_{k-j+4}}}-L^2{k}_G(k+2)(k+1)b_{k+2}}\\ +L^4{k}_w{b}_k-{\omega }^2{L}^4{\displaystyle \sum _{j=0}^k{\displaystyle \sum _{i=0}^j{\rho }_i^{*}{A}_{j-i}^{*}{b}_{k-j}}}]{\xi }^k=0\end{array}$ (9)

As mentioned before, the deformation shape of the beam with non-uniform cross-section was considered as follows:

$w\left(\xi \right)={\displaystyle \sum _{i=0}^{\infty }{b}_i{\xi }^i}=b_0+b_1\xi +b_2{\xi }^2+b_3{\xi }^3+b_4{\xi }^4+\mathrm{...}$ (10)

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:

$\begin{array}{l}w(\xi )=b_0\{1-\frac{{\xi }^4}{24E_0^{\ast }{I}_0^{\ast }}\left(L^4(k_w-{\omega }^2{\rho }_0^{\ast }{A}_0^{\ast })\right)-\frac{{\xi }^5}{120E_0^{\ast }{I}_0^{\ast }}\left(L^4{\omega }^2({\rho }_1^{\ast }{A}_0^{\ast }+{\rho }_0^{\ast }{A}_1^{\ast })+3L^4(k_0^{\ast }-{\omega }^2{\rho }_0^{\ast }{A}_0^{\ast })(\frac{I_1^{\ast }}{I_0^{\ast }}+\frac{E_1^{\ast }}{E_0^{\ast }})\right)\\ -\frac{{\xi }^6}{360E_0^{\ast }{I}_0^{\ast }}(L^4{\omega }^2({\rho }_0^{\ast }{A}_2^{\ast }+{\rho }_1^{\ast }{A}_1^{\ast }+{\rho }_2^{\ast }{A}_0^{\ast })+L^4(k_w-{\omega }^2{\rho }_0^{\ast }{A}_0^{\ast })\left(6\frac{I_2^{\ast }}{I_1^{\ast }}+6\frac{E_2^{\ast }}{E_1^{\ast }}+\frac{12E_1^{\ast }{I}_1^{\ast }+L^2{k}_G}{2E_0^{\ast }{I}_0^{\ast }}\right)\\ -2L^4(\frac{I_1^{\ast }}{I_0^{\ast }}+\frac{E_1^{\ast }}{E_0^{\ast }})\left({\omega }^2({\rho }_0^{\ast }{A}_1^{\ast }+{\rho }_1^{\ast }{A}_0^{\ast })+3(\frac{I_1^{\ast }}{I_0^{\ast }}+\frac{E_1^{\ast }}{E_0^{\ast }})(k_w-{\omega }^2{\rho }_0^{\ast }{A}_0^{\ast })\right))+\mathrm{...}\}\\ +b_1\left\{\xi -\frac{{\xi }^5}{120E_0^{\ast }{I}_0^{\ast }}\left(L^4(k_w-{\omega }^2{\rho }_0^{\ast }{A}_0^{\ast })\right)-\frac{{\xi }^6}{360E_0^{\ast }{I}_0^{\ast }}\left({\omega }^2{L}^4({\rho }_0^{\ast }{A}_1^{\ast }+{\rho }_1^{\ast }{A}_0^{\ast })+2L^4(\frac{I_1^{*}}{I_0^{*}}+\frac{E_1^{*}}{E_0^{*}})(k_w-{\omega }^2{\rho }_0^{\ast }{A}_0^{\ast })\right)+\mathrm{...}\right\}\\ +b_2\{{\xi }^2-\frac{{\xi }^4}{24E_0^{\ast }{I}_0^{\ast }}\left(4(E_0^{\ast }{I}_2^{\ast }+E_2^{\ast }{I}_0^{\ast }+E_1^{\ast }{I}_1^{\ast })-2L^2{k}_G\right)\\ -\frac{{\xi }^5}{120E_0^{\ast }{I}_0^{\ast }}\left(12(E_0^{\ast }{I}_3^{\ast }+E_1^{\ast }{I}_2^{\ast }+E_2^{\ast }{I}_1^{\ast }+E_3^{\ast }{I}_0^{\ast })-3(\frac{I_1^{\ast }}{I_0^{\ast }}+\frac{E_1^{\ast }}{E_0^{\ast }})(4(E_0^{\ast }{I}_2^{\ast }+E_2^{\ast }{I}_0^{\ast }+E_1^{\ast }{I}_1^{\ast })-2L^2{k}_G\right)\\ -\frac{{\xi }^6}{360E_0^{\ast }{I}_0^{\ast }}(24(E_0^{\ast }{I}_4^{\ast }+E_1^{\ast }{I}_3^{\ast }+E_3^{\ast }{I}_1^{\ast }+E_2^{\ast }{I}_2^{\ast }+E_4^{\ast }{I}_0^{\ast })-\left(6(\frac{I_1^{\ast }}{I_0^{\ast }}+\frac{E_1^{\ast }}{E_0^{\ast }})-\frac{L^2{k}_G-12E_1^{\ast }{I}_1^{\ast }}{2E_0^{\ast }{I}_0^{\ast }}\right)(4(E_0^{\ast }{I}_2^{\ast }+E_2^{\ast }{I}_0^{\ast }+E_1^{\ast }{I}_1^{\ast })-2L^2{k}_G)\\ \text{-2(}\frac{I_1^{\ast }}{I_0^{\ast }}+\frac{E_1^{\ast }}{E_0^{\ast }})(12(E_0^{\ast }{I}_3^{\ast }+E_1^{\ast }{I}_2^{\ast }+E_2^{\ast }{I}_1^{\ast }+E_3^{\ast }{I}_0^{\ast })-3\text{(}\frac{I_1^{\ast }}{I_0^{\ast }}+\frac{E_1^{\ast }}{E_0^{\ast }})(4(E_0^{\ast }{I}_2^{\ast }+E_2^{\ast }{I}_0^{\ast }+E_1^{\ast }{I}_1^{\ast })-2L^2{k}_G))+\mathrm{...}\}\\ +b_3\{{\xi }^3-\frac{{\xi }^4}{24E_0^{\ast }{I}_0^{\ast }}\left(12(E_0^{\ast }{I}_1^{\ast }+E_1^{\ast }{I}_0^{\ast })\right)-\frac{{\xi }^5}{120E_0^{\ast }{I}_0^{\ast }}\left(36(E_0^{\ast }{I}_2^{\ast }+E_1^{\ast }{I}_1^{\ast }+E_2^{\ast }{I}_0^{\ast })-36(\frac{I_1^{\ast }}{I_0^{\ast }}+\frac{E_1^{\ast }}{E_0^{\ast }})(E_0^{\ast }{I}_1^{\ast }+E_1^{\ast }{I}_0^{\ast })-6L^2{k}_G\right)\\ -\frac{{\xi }^6}{360E_0^{\ast }{I}_0^{\ast }}(72(E_0^{\ast }{I}_3^{\ast }+E_1^{\ast }{I}_2^{\ast }+E_2^{\ast }{I}_1^{\ast }+E_3^{\ast }{I}_0^{\ast })-72(\frac{I_2^{\ast }}{I_0^{\ast }}+\frac{E_2^{\ast }}{E_0^{\ast }}-\frac{L^2{k}_G-12E_1^{\ast }{I}_1^{\ast }}{12E_0^{\ast }{I}_0^{\ast }})(E_0^{\ast }{I}_1^{\ast }+E_1^{\ast }{I}_0^{\ast })\\ -2\text{(}\frac{I_1^{\ast }}{I_0^{\ast }}+\frac{E_1^{\ast }}{E_0^{\ast }})(36(E_0^{\ast }{I}_2^{\ast }+E_1^{\ast }{I}_1^{\ast }+E_2^{\ast }{I}_0^{\ast })-36(\frac{I_1^{\ast }}{I_0^{\ast }}+\frac{E_1^{\ast }}{E_0^{\ast }})(E_0^{\ast }{I}_1^{\ast }+E_1^{\ast }{I}_0^{\ast })-6L^2{k}_G))+\mathrm{...}\}\end{array}$ (11)

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 ($b_0,b_1,b_2,b_3$). Therefore, the fundamental solution of the motion equation for functionally graded beam with variable cross-section can be expressed in the matrix form:

$w(\xi )=\langle B\rangle \left\{A\right\}$ (12)

where

$\langle B\rangle =\langle \begin{array}{cccc}{w}_0\left(\xi \right)& w_1\left(\xi \right)& w_2\left(\xi \right)& w_3\left(\xi \right)\end{array}\rangle \left\{A\right\}=\left\{\begin{array}{c}{b}_0\\ b_1\\ b_2\\ b_3\end{array}\right\}$ (13 a, b)

In equation (13), $\langle B\rangle$ indicates the row vector including the fundamental solutions of the motion equation for free vibration ($w_i(\xi ),i=\mathrm{0,1,2,3}$), and {A} represents the column vector of four unknown parameters. Note that all terms of $w_i,i=\mathrm{0,1,2,3}$ are determined with the help of the symbolic software MATLAB. Finally, the transverse vibration characteristics are computed by implementing four boundary conditions corresponding to a single span Euler-Bernoulli beam (two at each end of the beam) and solving the eigenvalue problem. For most cases, both geometric boundary conditions (deflection and slope) and natural boundary conditions (shear force and bending moment) need to be specified. In the next stage, the approximate deflected shape of the member is derived by imposing the natural end conditions and real non-zero eigenvalues with respect to each vibrational mode. It is worth mentioning that the outcomes of the power series methodology are very sensitive to the number of terms considered in the power series approximations. An iterative procedure is thus an essential for estimating the required terms of power series expansions to obtain the natural frequencies related to the first vibration modes with the desired precision.

3 BOUNDARY CONDITIONS

It has to be notified that the four undetermined coefficients ($b_0,b_1,b_2,b_3$) are considered as the functions of the displacements of Degree of Freedom (DOF). These unknown parameters can be estimated by imposing the right- and left-end boundary conditions of beam element (two at each end of the beam). The four natural boundary conditions including displacements and rotations of the beam in the local coordinate are illustrated in Figure 2. As shown, the boundary conditions at the left end of the beam ($\xi =0$) including the vertical displacement (${\delta }_1$) and the bending rotation (${\theta }_1$) can be respectively written as:

1) At $x=0$ ($\xi =0$), then

$w(0)={\delta }_1\to b_0{w}_0(0)+b_1{w}_1(0)+b_2{w}_2(0)+b_3{w}_3(0)={\delta }_1$ (14a)

and

2) At $x=0$($\xi =0$), then $\frac{dw}{dx}={\theta }_1$

Then, one gets:

${\frac{dw(x)}{dx}|}_{x=0}={\theta }_1={\frac{1}{L}\frac{dw(\xi )}{d\xi }|}_{\xi =0}\to b_0\frac{{w^{\prime }}_0(0)}{L}+b_1\frac{{w^{\prime }}_1(0)}{L}+b_2\frac{{w^{\prime }}_2(0)}{L}+b_3\frac{{w^{\prime }}_3(0)}{L}={\theta }_1$ (14b)

Figure 2:
Boundary conditions of a beam element in global and local coordinate systems.

Referring to Figure 2 and using non-dimensional coordinate, the natural end conditions at the right end ($\xi =1$) of the elastic member are:

3) At $x=L$ ($\xi =1$), then

$w(L)=w(1)={\delta }_2\to b_0{w}_0(1)+b_1{w}_1(1)+b_2{w}_2(1)+b_3{w}_3(1)={\delta }_2$ (15a)

and

4) At $x=L$ ($\xi =1$), then $\frac{dw}{dx}={\theta }_2$

Then, we have:

${\frac{dw(x)}{dx}|}_{x=L}={\theta }_2={\frac{1}{L}\frac{dw(\xi )}{d\xi }|}_{\xi =1}\to b_0\frac{{w^{\prime }}_0(1)}{L}+b_1\frac{{w^{\prime }}_1(1)}{L}+b_2\frac{{w^{\prime }}_2(1)}{L}+b_3\frac{{w^{\prime }}_3(1)}{L}={\theta }_2$ (15b)

According to the boundary conditions presented above, a general system is derived for four linear equations with four undefined constants ($b_0,b_1,b_2,b_3$). The last four algebraic equations are equivalent to a matrix equation of the form:

$\left\{\Delta \right\}=\left[V\right]\left\{B\right\}$ (16)

in which

$\left\{\Delta \right\}=\left\{\begin{array}{c}{\delta }_1\\ {\theta }_1\\ {\delta }_2\\ {\theta }_2\end{array}\right\}[V]=\left[\begin{array}{cccc}{w}_0(0)& w_1(0)& w_2(0)& w_3(0)\\ \frac{{w^{\prime }}_0(0)}{L}& \frac{{w^{\prime }}_1(0)}{L}& \frac{{w^{\prime }}_2(0)}{L}& \frac{{w^{\prime }}_3(0)}{L}\\ w_0(1)& w_1(1)& w_2(1)& w_3(1)\\ \frac{{w^{\prime }}_0(1)}{L}& \frac{{w^{\prime }}_1(1)}{L}& \frac{{w^{\prime }}_2(1)}{L}& \frac{{w^{\prime }}_3(1)}{L}\end{array}\right]$ (17a, b)

where {Δ} represents the nodal displacement vector for the Euler-Bernoulli 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:

$\left\{B\right\}={\left[V\right]}^{-1}\left\{\Delta \right\}$ (18)

Regarding Eq. (11), it can be pointed out that $w_0(0)={w^{\prime }}_1(0)=1$ and $w_1(0)=w_2(0)=w_3(0)={w^{\prime }}_0(0)={w^{\prime }}_2(0)={w^{\prime }}_3(0)=0$; therefore, the inverse of matrix [V] is acquired as:

${\left[V\right]}^{-1}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& L& 0& 0\\ \frac{(w_0(1){w^{\prime }}_3(1)-w_3(1){w^{\prime }}_0(1))}{\left(w_3(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_3(1)\right)}& \frac{L(w_1(1){w^{\prime }}_3(1)-w_3(1){w^{\prime }}_1(1))}{\left(w_3(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_3(1)\right)}& -\frac{{w^{\prime }}_3(1)}{\left(w_3(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_3(1)\right)}& \frac{L\left(w_3(1)\right)}{\left(w_3(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_3(1)\right)}\\ -\frac{(w_0(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_0(1))}{\left(w_3(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_3(1)\right)}& -\frac{L(w_1(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_1(1))}{\left(w_3(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_3(1)\right)}& \frac{{w^{\prime }}_2(1)}{\left(w_3(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_3(1)\right)}& -\frac{L\left(w_2(1)\right)}{\left(w_3(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_3(1)\right)}\end{array}\right]$ (19)

The four independent constants ($b_0,b_1,b_2,b_3$) are then obtained by substituting Eq. (19) into Eq. (18) as:

$\left\{B\right\}=\left\{\begin{array}{c}{\delta }_1\\ L{\theta }_1\\ \frac{1}{C}\left((w_0(1){w^{\prime }}_3(1)-w_3(1){w^{\prime }}_0(1)){\delta }_1+L(w_1(1){w^{\prime }}_3(1)-w_3(1){w^{\prime }}_1(1)){\theta }_1-{w^{\prime }}_3(1){\delta }_2+L{w}_3(1){\theta }_2\right)\\ -\frac{1}{C}\left((w_0(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_0(1)){\delta }_1+L(w_1(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_1(1)){\theta }_1-{w^{\prime }}_2(1){\delta }_2+L{w}_2(1){\theta }_2\right)\end{array}\right\}$ (20)

where

$C=w_3(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_3(1)$ (21)

The mathematical expression describing the deformation shape of the axially functionally graded non-prismatic 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):

$\begin{array}{l}w(\epsilon )=\left\{w_0(\xi )+\left(\frac{(w_0(1){w^{\prime }}_3(1)-w_3(1){w^{\prime }}_0(1))}{\left(w_3(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_3(1)\right)}\right)w_2(\xi )+\left(-\frac{(w_0(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_0(1)}{\left(w_3(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_3(1)\right)}\right)w_3(\xi )\right\}{\delta }_1\\ +\left\{L{w}_1(\xi )+\left(\frac{L(w_1(1){w^{\prime }}_3(1)-w_3(1){w^{\prime }}_1(1))}{\left(w_3(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_3(1)\right)}\right)w_2(\xi )+\left(-\frac{L(w_1(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_1(1))}{\left(w_3(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_3(1)\right)}\right)w_3(\xi )\right\}{\theta }_1\\ +\left\{\left(-\frac{{w^{\prime }}_3(1)}{\left(w_3(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_3(1)\right)}\right)w_2(\xi )+\left(\frac{{w^{\prime }}_2{}_2(1)}{\left(w_3(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_3(1)\right)}\right)w_3(\xi )\right\}{\delta }_2\\ +\left\{\left(\frac{L\left(w_3(1)\right)}{\left(w_3(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_3(1)\right)}\right)w_2(\xi )+\left(-\frac{L\left(w_2(1)\right)}{\left(w_3(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_3(1)\right)}\right)w_3(\xi )\right\}{\theta }_2\end{array}$ (22)

In this study, the free vibration and liner buckling analyses are conducted for FGM non-prismatic beam having three different end supports including a beam with left end fixed and the other end free (fixed-free), a beam clamped at left end and hinged at right one (clamped-hinged), and a beam hinged at both ends (hinged-hinged). 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 clamped-free beam can be obtained in non-dimensional coordinate as:

$\begin{array}{l}w(\xi )=\{\left(\frac{L\left(w_2(1){w^{\prime }}_3(1){w^{‴}}_2(1)-w_2(1){w^{\prime }}_2(1){w^{‴}}_3(1)\right)}{\left({w^{‴}}_2(1){w^{\prime }}_3(1){w^{\prime }}_2(1)-{w^{\prime }}_2(1){w^{\prime }}_2(1){w^{‴}}_3(1)\right)}\right)\left(-{w^{\prime }}_3(1){w^{\prime }}_2(1)w_2(\xi )+{w^{\prime }}_2(1){w^{\prime }}_2(1)w_3(\xi )\right)\\ +L\left(w_2(1){w^{\prime }}_3(1)w_2(\xi )-w_2(1){w^{\prime }}_2(1)w_3(\xi )\right)\}\frac{{\theta }_2}{{w^{\prime }}_2(1)}\end{array}$ (23)

The explicit expression of vibrational shape for fixed-hinged members is determined by preventing the transverse displacement and rotation at the fixed end ($x=0;\xi =0$), as well as the vertical displacement and bending moment at the pinned end ($x=L;\xi =1$) as follows:

$w(\xi )=\left(\frac{L}{w_3(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_3(1)}\right)\left\{\frac{w_2(1){w^{″}}_3(1)}{{w^{″}}_2(1)}{w}_2(\xi )-w_2(1)w_3(\xi )\right\}{\theta }_2$ (24)

For simply supported beams, the vertical displacement and bending moment at the both ends are equal to zero. Therefore, the vibrational mode shape becomes:

$\begin{array}{l}w(\xi )=\frac{L}{\left(w_3(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_3(1)\right)\left({w^{\prime }}_1(1){w^{″}}_3(1)-{w^{″}}_1(1){w^{\prime }}_3(1)\right)}\{({w^{″}}_1(1){w^{\prime }}_3(1)-{w^{\prime }}_1(1){w^{″}}_3(1))w_2(1)w_3(\xi )\\ +\left((w_3(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_3(1)){w^{″}}_3(1)w_1(\xi )+(w_3(1){w^{″}}_1(1){w^{\prime }}_2(1)-w_2(1){w^{\prime }}_1(1){w^{″}}_3(1))w_3(\xi )\right)\}{\theta }_2\end{array}$ (25)

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 non-dimensional coordinate can be derived in terms of the rotation of right end of the member (${\theta }_2$), which can be written as bellow:

$w(\xi )=\varphi (\xi ){\theta }_2$ (26)

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 Rayleigh-Ritz method is employed to obtain buckling load of non-prismatic 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 Rayleigh-Ritz method. Therefore, the obtained vibrational mode (Eqs. (23-25)) can be adopted as the deformed shape of column for the linear stability analysis.

4 THE RAYLEIGH-RITZ METHOD

In this section, the value of critical buckling load of the beam with variable cross-section is determined by the Rayleigh-Ritz method based on the variation of total potential energy, which is:

$\delta \Pi =\delta (U_l+U_0+U_f)=0$ (27)

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 non-uniform member in local coordinate can be written as:

$\begin{array}{l}\prod =\frac{1}{2L^3}{\displaystyle {\int }_0^{1}E(\xi )I(\xi )}{\left(\frac{d^{2}w}{d{\xi }^2}\right)}^{2}d\xi +\frac{P}{2L}{\displaystyle {\int }_0^1{\left(\frac{dw}{d\xi }\right)}^{2}d\xi }+k_w\frac{L}{2}{\displaystyle {\int }_0^1{w}^2(\xi )d\xi }\\ +\frac{k_G}{2L}{\displaystyle {\int }_0^1{\left(\frac{dw}{d\xi }\right)}^{2}d\xi }\end{array}$ (28)

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 Euler-Bernoulli beam. In the current study, it is assumed that the displacement equation obtained in the previous section (Eqs. (23-25)) is the approximate form of buckled shape of elastic beam. Thus, Eq. (27) can be replaced by:

$\frac{d(U_l+U_0+U_f)}{d{\theta }_2}\delta {\theta }_2=0$ (29)

Since $\delta {\theta }_2$ is an arbitrary virtual displacement in terms of bending rotation, the above formulation can also be expressed as:

$\frac{d(U_l+U_0+U_f)}{d{\theta }_2}=0$ (30)

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 cross-sections whether resting on elastic foundation or not, several examples are conducted.

At the first, in the absence of elastic foundation and non-homogeneous 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 non-prismatic 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 non-prismatic Euler-Bernoulli beams in the presence of one or two-parameter 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 semi-analytical 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 Rayleigh-Ritz 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: Non-prismatic 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 cross-section. To this aim, as shown in Figure 3, three non-uniform members are investigated with different boundary conditions (fixed-free, hinged-hinged and fixed-hinged). In all of the presented cases, the cross-section is in the form of rectangle with width b ₀ and depth h ₀ , 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,$\xi =0$) and the right one (x=L, $\xi =1$) of the beam are respectively indicated with the subscripts 0 and 1, relating to the dimensionless coordinate system ($\xi =x/L$).

In the first case (Case A), the width of the beam’s section at the left support (b0) is linearly diminished to ($b_1=(1-\beta )b_0$) at the other end. The height of the beam remains constant along the member’s length. Therefore, the obtained section at the right end becomes a rectangle with a depth of (h0) and width of (b1).

In the second case (Case B), the height of the beam’s section is (h0) at the left support and linearly decreased to ($h_1=(1-\beta )h_0$) at the other end. However, the width of the beam (b0) remains constant.

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 $b_1/b_0=h_1/h_0=1-\beta$.

Figure 3:
Members with various end conditions; left and right sections corresponding to cases A to C.

In this example, the tapering parameter can change from $\beta =0$ (prismatic beam) to $\beta =0.1-0.9$ (non-uniform ones). For mentioned three cases, the minor axis moment of inertia and the cross-sectional area can be represented in the following forms:

$CaseA:I\left(\xi \right)=I_0\left(1-\beta \xi \right);A\left(\xi \right)=A_0\left(1-\beta \xi \right)$ (31a)

$CaseB:I\left(\xi \right)=I_0{\left(1-\beta \xi \right)}^3;A\left(\xi \right)=A_0\left(1-\beta \xi \right)$ (31b)

$CaseC:I\left(\xi \right)=I_0{\left(1-\beta \xi \right)}^4;A\left(\xi \right)=A_0{\left(1-\beta \xi \right)}^2$ (31c)

In which A0 and I0 are respectively cross-sectional area and moment of inertia at the left support (x=0, $\xi =0$). They are defined as:

$I_0=\frac{b_0{h}_0^3}{12}and{A}_0=b_0{h}_0$ (32)

In order to simplify the illustration of the obtained results, the dimensionless natural frequency and buckling load parameters are used as follows:

$P_{nor}=P_{cr}\times \frac{L^2}{E_0{I}_0}$ (33a)

${\omega }_{nor}=\omega \times \sqrt{\frac{{\rho }_0{A}_0{L}^4}{E_0{I}_0}}$ (33b)

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 non-dimensional buckling load parameter (P _nor ) of the tapered beam based on Euler-Bernoulli beam assumptions for two cases (B and C ) and two different boundary conditions (pinned-ended and fixed-free) 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 Euler-Bernoulli 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 non-prismatic and prismatic members has been modeled using BEAM188 in ANSYS software. This member is a two-node 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 pinned-pinned double-tapered 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 non-dimensional vibration frequency is close to beam results of ANSYS software (ω_nor=6.2086, error Δ<0.5%). The needed CPU time is only 133.62s.

Table 1:
Convergence of the proposed numerical technique in determination of the first non-dimensional vibration frequency (ω_nor) and normalized buckling load (P_nor) of tapered beam (Case B) with two different end conditions.

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	Hinged-Hinged	0	(ωnor)	9.6730	9.8696	9.8696	9.8696	9.8696	9.7980
			(Pnor)	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
			(Pnor)	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
			(Pnor)	2.5324	2.6774	2.6789	2.6789	2.6789	2.6638
	Clamped-Free	0	(ωnor)	3.5073	3.5160	3.5160	3.5160	3.5160	3.5040
			(Pnor)	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
			(Pnor)	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
			(Pnor)	0.7016	1.1273	1.0999	1.0999	1.0999	1.0970

Table 2:
Effect of power series terms on the non-dimensional natural frequency (ω_nor) and buckling load (P_nor) of double tapered beam with two different end conditions and CPU time comparisons

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	Hinged-Hinged	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
			(Pnor)	(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
			(Pnor)	(4.51s)	(12.75s)	(31.28s)	(148.12s)	(3607.1s)
	Clamped-Free	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
			(Pnor)	(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
			(Pnor)	(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 non-dimensional 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 closed-form 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:

$\Delta =\frac{\left|P_{cr}^{PSM}-P_{cr}^{\mathrm{Re}f}\right|}{P_{cr}^{\mathrm{Re}f}}\times 100$ (34)

As shown in Tables 3 and 4, the values of buckling loads and natural frequencies for non-uniform 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).

Table 3:
Normalized natural frequencies for columns of different cases (A, B and C) with various boundary conditions

Case	β	Free Vibration Analysis
		Hinged-Hinged			Clamped-Free			Clamped-Hinged
		ANSYS	Present Method	$\Delta (\%)$	ANSYS	Present Method	$\Delta (\%)$	ANSYS	Present Method	$\Delta (\%)$
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

Table 4:
Normalized critical buckling loads for columns of different cases (A, B and C) with various boundary conditions

Case	β	Stability Analysis
		Hinged-Hinged			Clamped-Free			Clamped-Hinged
		Wang et al. (2005 )	Present Method	$\Delta (\%)$	Wang et al. (2005 )	Present Method	$\Delta (\%)$	Wang et al. (2005 )	Present Method	$\Delta (\%)$
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 3-4). 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 cross-section is sensitive to the boundary conditions. Regarding the pinned-pinned 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 non-uniformity parameter. In this regard, it can be stated that the free vibration behavior of non-uniform is not predictable.

Comparing the buckling results of prismatic beam tabulated in Table 1 with those of non-prismatic 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 non-prismatic beams of cases A, B and C are 16%, 46% and 58% of fixed-free prismatic beam, respectively. This can be explained by the fact that increasing in width and/or height taper ratios causes the reduction in cross-sectional 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 cross-section equals to 4, while the power-index for other expressions corresponding to the moment of inertia is 1 or 3. Therefore, the reduction of elastic stability capacity of non-prismatic member caused by increasing the width and height non-uniformity ratios ($\beta$) simultaneously is more noticeable than that of others. This statement can be observed in Table 4.

5.2 Example 2- Non-prismatic member with AFG material

In this section, to study the effects of material non-homogeneity on linear stability and free vibration, a non-prismatic member made up of axially functionally graded material is considered. Axially functionally graded Euler-Bernoulli beam is assumed with exponentially tapered cross-section which is shown in Figure 4. Therefore, the geometrical properties including the cross-section area and the moment of inertia vary along the beam axis with exponential functions. Thus, the variation of minor axis moment of inertia and cross-sectional area in the local coordinate are described as follows:

$I\text{(ξ)}=I_0{e}^{\alpha \xi }$ (35a)

$A\text{(ξ)}=A_0{e}^{\alpha \xi }$ (35b)

Figure 4:
(a) Schematic representation of AFG beam with exponentially tapered cross-section, (b) Exponentially varying width along member, (c) Constant thickness through beam’s length.

The non-uniformity parameter (α) can change from zero (prismatic beam) to a range of [-2 to -0.1] for non-uniform 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 cross-sectional area and moment of inertia expressions (Eq. (35)) are defined as:

$I(\xi )=I_0{\displaystyle \sum _{i=0}^{\infty }\frac{{\alpha }^i}{i!}{\xi }^i}$ (36a)

$A(\xi )=A_0{\displaystyle \sum _{i=0}^{\infty }\frac{{\alpha }^i}{i!}{\xi }^i}$ (36b)

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₀=200GPa, ρ₀ =5700 kg/m3;

Al: E₁=70GPa, ρ₁ =2702 kg/m3;

The variation of Young’s modulus of elasticity and the material density along the beam axis are defined with the following power-law formulations:

$E(\xi )=E_0+(E_1-E_0){\xi }^m$ (37a)

$\rho (\xi )={\rho }_0+({\rho }_1-{\rho }_0){\xi }^m$ (37b)

In the last expression, m is the non-homogeneity 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 power-law 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 non-prismatic members with exponential variation of geometrical properties in the first stage of the current example, the first three non-dimensional natural frequencies of simply supported homogenous beam for two different non-uniformity 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).

Table 5:
Natural frequencies comparison of power series method and finite element results for simply supported beam

α	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 non-prismatic beam with different non-uniformity ratios (α) should be considered to estimate the required numbers of the terms in power series. Table 6 indicates the dimensionless parameter of critical load ($P_{nor}$) evaluated by the proposed numerical method and compared with those of Wang et al. (2005). Further, Table 6 demonstrates the effect of the degree of the power series on convergence. Besides, the relative errors between the closed-form solution suggested by Wang et al. (2005) and the proposed methodology are arranged in the mentioned table.

Table 6:
The effect of number of power series terms (P _nor ) on precision of normalized critical axial load parameter for homogenous exponentially tapered beams with different boundary conditions

End Conditions	α	Number of terms of power series (PSM)				Wang et al. (2005 )		$\Delta (\%)$
		n=10	n=20	n=30	n=40	n=50
Hinged-Hinged	-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 non-uniformity parameter (α) on dimensionless natural frequencies and critical loads of AFG non-prismatic beam by considering the non-homogeneity indexes m=1, 2 and 3 for simply supported, fixed-free, and fixed-pinned 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 non-uniformity parameter from 0.0 to -2 leads to a reduction in the moment of inertia and cross-sectional 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 non-homogeneity index (m) and tapering parameter (α) concurrently is impossible. This fact can be easily observed in Figure 5.

Figure 5:
Variations of the non-dimensional free transverse frequency of AFG exponentially tapered beam with respect to non-uniformity parameter: (a) Simply supported, (b) Cantilever, (c) Clamped-Pinned

Comparing the results of Figure 6, it can be concluded that when non-homogeneity 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 power-law exponent the buckling resistance of uniform beam (α=0) and non-prismatic beam with α=-2 are considered as the most and least amount, respectively, since a decrease in non-uniformity parameter leads to a reduction in cross-sectional area and moment of inertia, and consequently in the stiffness of the elastic member.

Figure 6:
Variations of the non-dimensional buckling load parameter of AFG exponentially tapered beam with respect to non-uniformity parameter: (a) Simply supported, (b) Cantilever, (c) Clamped-Pinned

5.3 Example 3- AFG Non-prismatic beam on continuous Winkler-Pasternak foundation

To study the effects of Winkler-Pasternak foundation on linear stability and free vibration, another example is considered including axially non-homogeneous and homogeneous beams. In this regard, an elastic Euler-Bernoulli beam with rectangular cross-section whose width and thickness taper linearly along the member axis with similar taper ratio is taken into account. Therefore, the cross-sectional 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 ($\beta$) namely 0 and 0.5; where the first one is a prismatic beam, while the latter is a non-prismatic one. Additionally, the natural frequencies and critical buckling loads are conducted for two cases: axially non-homogeneous and homogeneous beams. In the case of axially FG member, the distribution of modulus of elasticity and density of material are assumed to vary in the longitudinal direction with a power law formulation as expressed in Eq. (37). In this case, the material non-homogeneity parameter (m ) is equal to 2. In the current example, the linear stability and free vibration problems for simply supported and clamped-pinned beams are investigated. The following non-dimensional parameters related to Winkler and Pasternak constants are respectively introduced as:

${\overline{k}}_w=k_w\frac{L^4}{E_0{I}_0}$ (38a)

${\overline{k}}_G=k_G\frac{L^2}{E_0{I}_0}$ (38b)

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 two-parameter foundation. Therefore, a converge study is conducted for linear stability and free vibration characteristics of the considered hinged-hinged member with different foundation parameters: ${\overline{k}}_w$=0 and 30, as well as ${\overline{k}}_G$=0 and 3. The convergent performance of the results relating to vibration and stability analyses which are obtained by the present method versus the number of terms (n) considered in the power series is respectively plotted in Figures 7 and 8. Moreover, the needed Central Processing Unit (CPU) time to accomplish the computational procedure is shown for axially functionally graded double tapered beam lying on Winkler-Pasternak foundation. As presented, the value of normalized free vibration frequencies and buckling loads gradually approach the acceptable solutions when the taken number of terms of the power series is exceeded over 20. Note that 30 terms is a good compromise for equivalent accuracy of different cases.

Next, the variation of the frequencies of vibration of hinged-hinged and fixed-hinged Euler-Bernoulli beams resting on Winkler type of foundation versus the elastic foundation constant (${\overline{k}}_w$) are respectively presented in Figure 9 a-b. Further, Figure 10 shows the variation of the lowest buckling load parameters with respect to the elastic foundation constant for simply supported and clamped-pinned beams. As can be seen, the variation of Winkler elastic foundation parameter has a significant influence on the bending vibration and linear stability behavior of both beams under different circumstance. As can be seen in Figure 9, the natural frequency parameters are increased due to the influence of Winkler type of foundation. Furthermore, as displayed in Figure 9, the interaction curves are interchanged for different values of elastic foundation parameters (${\overline{k}}_w$). Like the case of the natural frequencies, the critical buckling load parameters corresponding to the first mode are also increased as the stiffness of the elastic foundation increases. In the other words, based on the numerical outcomes the elastic foundation plays a stabilizing effect on the stability characteristics of axially non-homogeneous and homogeneous beams with constant or variable cross-section. One can also remark that for axially non-homogeneous and homogeneous prismatic beams, the buckling load increases linearly due to the effect of uniform Winkler foundation although this variation is non-linear for double-tapered members. It is found from Figure 10 that the variation of stability behavior for the fixed-hinged beam is similar to that of the hinged-hinged beam, but the former beam is obviously more stable than the latter one. Further, the corresponding buckling load for the beam having constant material properties and uniform cross-section is the highest for any value of Winkler’s parameters while it is the lowest for non-prismatic beam having property according to power-law with index m = 2. Regarding AFG non-prismatic beam, the first buckling load increases slightly as the non-dimensional spring constant (${\overline{k}}_w$) increases over 50

Figure 7:
Convergence of non-dimensional natural frequency of simply supported Euler-Bernoulli beam resting on one or two parameter elastic foundation: (a) homogeneous beam with constant cross-section, (b) AFG beam with constant cross-section, (c) homogeneous tapered beam, (d) AFG tapered beam.

Figure 8:
Convergence of normalized buckling load parameter of simply supported Euler-Bernoulli beam resting on one or two parameter elastic foundation: (a) homogeneous beam with constant cross-section, (b) AFG beam with constant cross-section, (c) homogeneous tapered beam, (d) AFG tapered beam.

Figure 9:
Influence of Winkler type elastic foundation modulus on the natural frequency parameters of (a) pinned-ended beams, (b) clamped-pinned beams

Figure 10:
Influence of Winkler type elastic foundation modulus on the critical buckling load parameters of (a) pinned-ended columns, (b) clamped-pinned columns

Finally, in order to investigate the influence of uniform Winkler-Pasternak foundation on the behavior of simply supported Euler-Bernoulli beam, the first transverse natural frequencies and the lowest normalized critical buckling loads for different combinations of elastic foundation parameters i.e. ${\overline{k}}_w$and ${\overline{k}}_G$are calculated through the proposed novel technique and the results are shown in Table 7.

Table 7:
Natural frequency and Critical axial load parameters (${\omega }_{nor},P_{nor}$) for simply supported beams resting on two-parameter elastic foundation.

Material	(β)	The Pasternak foundation constant (${\overline{k}}_G$)	Non-dimensional bending frequency					Non-dimensional critical buckling load
			The Winkler foundation constant (${\overline{k}}_w$)					The Winkler foundation constant (${\overline{k}}_w$)
			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 non-prismatic axially functionally graded beams based on Euler-Bernoulli theory is proposed. This new technique combines the power series method and Rayleigh-Ritz approach and applied to non-prismatic 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 cross-section. 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 above-mentioned method. The vibrational mode shapes of an elastic member are similar to the buckling ones. Therefore, the obtained deflected shapes of the considered non-prismatic 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 non-prismatic 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 Rayleigh-Ritz methods in free vibration and linear stability analyses of non-prismatic AFG beams on two-parameter 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 Winkler-Pasternak parameter on linear buckling resistance and the natural frequency of Euler-Bernoulli beam with varying cross-section are comprehensively surveyed. It should be pointed out that material properties of functionally graded beam with variable cross-section are supposed to vary through longitudinal direction of the constituents according to simple power-law distribution (P-FGM). 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 ceramic-metal beams, buckling resistance increases with increase of functionally graded material content for power-law property distribution. In addition, the numerical outcomes reveal that the elastic foundation enhance the stability characteristics of axially non-homogeneous and homogeneous beams with constant or variable cross-section.

References

• ANSYS, Version 5.4, Swanson Analysis System, Inc, 2007.
• Arbabi F, Li F. (1991) “Buckling of variable cross-section columns: integral equation approach”. Journal of Structural Engineering, 117 (8): 2426-2441.
• Arefi M. (2013). “Nonlinear thermoelastic analysis of thick-walled functionally graded piezoelectric cylinder”. Acta Mechanica, 224(11): 2771-2783.
• Arefi M. (2015). “Nonlinear electromechanical stability of functionally graded circular plate integrated with functionally graded piezoelectric layers under radial compressive”, Latin American Journal of Solids and Structures, 12(9): 1653-1665.
• Arefi M, Zamani M.H., Kiani M. (2018) “Size-dependent free vibration analysis of three layered exponentially graded nanoplate with piezomagnetic face-sheets resting on Pasternak’s foundation”. Journal of Intelligent Material Systems and Structures 29(5): 774-786.
• Asgarian B., Soltani M., Mohri F. (2013) “Lateral-torsional buckling of tapered thin-walled beams with arbitrary cross-sections”. Thin-Walled Structures, 62: 96-108.
• Atay M.T. (2009) “Determination of critical buckling loads for variable stiffness Euler columns using homotopy perturbation method”. International Journal of Nonlinear Sciences and Numerical Simulation, 10(2): 199-206.
• Ayvaz Y. and Ozgan K. (2002). “Application of modified Vlasov model to free vibration analysis of beam resting on elastic foundation”, Journal of Sound and Vibration, 255: 111-127.
• Chakraborty A, Gopalakrishnan S, Reddy JN. (2003). “A new beam finite element for the analysis of functionally graded materials”. International Journal Mechanical Science, 45: 519-39.
• Coşkun S.B, Atay M.T. (2009). “Determination of critical buckling load for elastic columns of constant and variable cross-sections using variational iteration method”. Computers and Mathematic with Applications, 58(11-12): 2260-2266.
• Ebrahimi F., Mokhtari M. (2015) “Vibration analysis of spinning exponentially functionally graded Timoshenko beams based on differential transform method, Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering; 229(14): 2559-2571.
• Ebrahimi F. and Barati M.R. (2017) “Vibration analysis of viscoelastic inhomogeneous nanobeams resting on a viscoelastic foundation based on nonlocal strain gradient theory incorporating surface and thermal effects”. Acta Mechanica, 228(3): 1197-1210.
• Ece M.C, Aydogdu M, Taskin V. (2007). “Vibration of a variable cross-section beam”. Mechanics Research Communications, 34: 78-84.
• Eisenberger M and Clastornik J. (1987) “Vibrations and buckling of a beam on a variable Winkler elastic foundation”, Journal of Engineering Mechanics, Journal of Sound and Vibration; 115(2): 233-241.
• Ermopulos J. (1977) “Equivalent buckling length of non-uniform members”. Journal of Constructional Steel Research 1977; 42(4):141-158.
• Fallah A and Aghdam M.M. (2011). “Nonlinear free vibration and post-buckling analysis of functionally graded beams on nonlinear elastic foundation”, European Journal of Mechanics A/Solids, 30: 571-583.
• Gere J.M, Carter W.O. (1962) “Critical buckling loads for tapered columns”. Journal of Structural Engineering ASC, 88(1): 1-11.
• Girgin Z.C and Girgin K. (2005). “A numerical method for static or dynamic stiffness matrix of non-uniform members resting on variable elastic foundations”, Engineering Structures, 27: 1373-1384.
• Hassan M. T. and Nassar M. (2015) “Analysis of stressed Timoshenko beams on two parameter foundations”, KSCE Journal of Civil Engineering, 19(1): 173-179.
• Iromenger M.J. (1980). “Finite difference buckling analysis of non-uniform columns”. Computers & Structures, 12(5): 741-748.
• Jategaonkar R and Chehil DS. (1989). “Natural frequencies of a beam with varying section properties”, Journal of Sound and Vibration, 133: 303-322.
• Khan A.A, Alam M.N, Rahman N, Wajid M. (2016). “Finite element modelling for static and free vibration response of functionally graded beam”, Latin American Journal of Solids and Structures, 13: 690-714.
• Khaniki H.B. and Rajasekaran S. (2018). “Mechanical analysis of non-uniform bi-directional functionally graded intelligent micro-beams using modified couple stress theory”. Materials Research Express, 5(5), 055703.
• Kim N-II, Lee J-H and Kim M-Y. (2005). “Exact dynamic stiffness matrix of non-symmetric thin-walled beams on elastic foundation using power series method”, Advances in Engineering Software, 36: 518-532.
• Li XF, Kang YA, Wu JX. (2013). “Exact frequency equations of free vibration of exponentially graded beams”. Applied Acoustic, 74: 413-420.
• Li X., Li L., Hu Y. (2018). “Instability of functionally graded micro-beams via micro-structure-dependent beam theory”. Applied Mathematics and Mechanics (English Edition), 39: 923-952.
• Malekzadeh P and Karami G. (2008). “A mixed differential quadrature and finite element free vibration and buckling analysis of thick beams on two-parameter elastic foundations”, Applied Mathematical Modelling, 32: 1381-1394.
• MATLAB 7.6, The MathWorks, Inc., Natick, Massachusetts, United States.
• Mirzabeigy A. (2014) “Semi-analytical approach for free vibration analysis of variable cross-section beams resting on an elastic foundation and under axial force”, International Journal of Engineering, 27(3): 385-394.
• Mohanty S.C., Dash R.R., Rout T. (2012) “Static and Dynamic Stability Analysis of a Functionally Graded Timoshenko Beam”, International Journal of Structural Stability and Dynamics. 12(4) 1250025.
• Okay F, Atay M.T, Coçkun S.B. (2010). “Determination of buckling loads and mode shapes of a heavy vertical column under its own weight using the variational iteration method”. International Journal of Nonlinear Sciences and Numerical Simulation, 11(10): 851-857.
• Pradhan K.K, Chakraverty S. (2013). “Free vibration of Euler and Timoshenko functionally graded beams by Rayleigh-Ritz method”. Composite Part B, 51: 175-184.
• Rahai A.R and Kazemi S. (2008). “Buckling analysis of non-prismatic column based on modified vibration method”. Communications in Nonlinear Science and Numerical Simulation, 13: 1721-1735.
• Ruta P. (2006). “The application of Chebyshev Polynomials to the solution of the non-prismatic Timoshenko beam vibration problem”, Journal of Sound and Vibration, 296(1-2): 243-263.
• Sampaio J.H.B and Hundhausen J.R. (1998). “A mathematical model and analytical solution for buckling of inclined beam columns”. Applied Mathematical Modeling, 22(6): 405-421.
• Sankar BV. (2001). “An elasticity solution for functionally graded beams”. Composite Science and Technology, 61: 689-696.
• Shahba A., Attarnejad R., Hajilar S. (2013). “A mechanical-based solution for axially functionally graded tapered Euler-Bernoulli beams”. Mechanics of Advanced Materials and Structures, 20: 696-707.
• Shahba A., Attarnejad R., Zarrinzadeh H. (2013). “Free Vibration Analysis of Centrifugally Stiffened Tapered Functionally Graded Beams”. Mechanics of Advanced Materials and Structures, 20(5): 331-338.
• Shvartsman B. and Majak J. (2016) “Numerical method for stability analysis of functionally graded beams on elastic foundation”, Applied Mathematical Modelling, 44: 3713-3719.
• Simsek M, Kocatürk T. (2009). “Free and forced vibration of a functionally graded beam subjected to a concentrated moving harmonic load”. Composite Structures, 90: 465-73.
• Sina S.A, Navazi H.M, Haddadpour H. (2009). “An analytical method for free vibration analysis of functionally graded beams”. Material Design, 30:741-7.
• Singh K.V. and Li G. (2009). “Buckling of functionally graded and elastically restrained non-uniform columns”. Composites Part B, 40(5): 393-403.
• Soltani M., Asgarian B., Mohri F. (2014) “Elastic instability and free vibration analyses of tapered thin-walled beams by the power series method”, Journal of Constructional Steel Research, 96:106-126.
• Soltani M. and Asgarian B. (2019). “Finite element formulation for linear stability analysis of axially functionally graded non-prismatic Timoshenko beam”. International Journal of Structural Stability and dynamics, 19(2): 30 pages.
• Tong X. and Tabarrok B. (1995). “Vibration analysis of Timoshenko beams with non-homogeneity and varying cross-section”. Journal of Sound and Vibration, 186:821-35.
• Tsiatas G.C. (2014) “A new efficient method to evaluate exact stiffness and mass matrices of non-uniform beams resting on an elastic foundation”, Archive of Applied Mechanic. 84: 615-623.
• Wang C.M, Wang, C.Y and Reddy J.N. (2005). “Exact Solutions for Buckling of Structural Members”. CRC Press LLC, Florida.
• Wang ZH, Wang XH, Xu GD, Cheng S and Zeng T. (2016) “Free vibration of two-directional functionally graded beams”, Composite Structures; 135:191-198.
• Wu L, Wang QS, Elishakoff I. (2005). “Semi-inverse method for axially functionally graded beams with an anti-symmetric vibration mode”. Journal of Sound and Vibration, 284(3): 1190-202.
• Zhu B and Leung A.Y.T. (2009) “Linear and nonlinear vibration of non-uniform beams on two-parameter foundations using p-elements”, Computers and Geotechnics, 36: 743-750.

• Available online: February 21, 2019

Publication Dates

• Publication in this collection
  08 Apr 2019
• Date of issue
  2019

History

• Received
  11 Nov 2017
• Reviewed
  03 Jan 2019
• Accepted
  18 Feb 2019