Document Type : Research Paper
Author
Department of Civil Engineering, University of Kashan, Kashan, Iran
Abstract
Keywords
Finite Element Modeling for Buckling Analysis of Tapered Axially Functionally Graded Timoshenko Beam on Elastic Foundation
M. Soltani
^{ }Department of Civil Engineering, University of Kashan, Kashan, Iran
KEYWORDS 

ABSTRACT 
Power Series Method Shape Functions Buckling Load Timoshenko Beam Functionally Graded Materials

In this study, an efficient finite element model with two degrees of freedom per node is developed for buckling analysis of axially functionally graded (AFG) tapered Timoshenko beams resting on Winkler elastic foundation. For this, the shape functions are exactly acquired through solving the system of equilibrium equations of the Timoshenko beam employing the power series expansions of displacement components. The element stiffness matrix is then formulated by applying the developed shape functions to the total potential energy along the element axis. It is demonstrated that the resulting shape functions, in comparison with Hermitian cubic interpolation functions, are proportional to the mechanical features of the beam element, including the geometrical properties, material characteristics, as well as the critical axial load. An exhaustive numerical example is implemented to clarify the efficiency and simplicity of the proposed mathematical methodology. Furthermore, the effects of end conditions, material gradient, Winkler parameter, tapering ratio, and aspect ratio on the critical buckling load of AFG tapered Timoshenko beam are studied in detail. The numerical outcomes reveal that the elastic foundation enhances the stability characteristics of axially nonhomogeneous and homogeneous beams with constant or variable crosssection. Moreover, the results show that the influence of nonuniformity in the crosssection and axially inhomogeneity in material characteristics play significant roles in the linear stability behavior of Timoshenko beams subjected to different boundary conditions. 
Knowing the stability characteristics and the specific buckling load of the axially loaded members are of significant importance in design considerations. Elastic tapered beams are a class of important structural components, which have wide applications in civil, mechanical, and aeronautical structures. This is because of their ability to increase both strength and stability, reduce the whole weight of the structure, and satisfy aesthetic necessities. The stability problem has been investigated by many researchers via the EulerBernoulli and the Timoshenko beam theories. Based on the EulerBernoulli theory (EBT), the influence of transverse shear deformation is neglected, and only the effect of flexural deformation is considered. The Euler–Bernoulli beam model is usually adopted for the bending of slender and long beams. To overcome the defects of EBT, especially when the beam is moderately deep with a small lengthtodepth ratio, for instance, towers, moveable arms, and antenna, researchers usually use Timoshenko beam assumptions, in which the effects of rotatory inertia, transverse shear, and bending deformations are taken into account. Furthermore, the use of functionally graded materials (FGMs) has been increasing in automotive, civil, electronic, optical, and mechanical industries due to their noticeable characteristics such as elimination or minimization of interfacial stress concentration, thermal resistance and optimal distribution of weight.
Therefore, until now, several types of research have been performed to assess the static, stability, and the vibration behaviors of FGMs and/or homogenous beams. In this context, Simsek [1] studied the free vibration of FG beams based on the firstorder and different higherorder shear deformation theories. The free vibration and linear buckling analyses of nonprismatic Timoshenko beams with axially varying materials subjected to various end conditions were perused by Shahba et al. [2] using the finite element procedure. Through the RayleighRitz method, the linear buckling resistance of tapered microcolumns having different taper ratios was studied by Akgoz and Civalek [3]. Arefi and Rahimi [4] employed the Adomian Decomposition Method (ADM) on the nonlinear analysis of the FG beam with variable thicknesses. The surface effects on the nonlinear free vibration of elastically restrained nonlocal beams with variable crosssections were examined by Malekzadeh and Shojaee [5]. An efficient finite element model for static and dynamic analysis of functionally graded piezoelectric beams was introduced by LezgyNazargah et al. [6]. Rajasekaran and Tochaei [7] applied the Differential Transformation Element Method (DTEM) and Differential Quadrature Element of Lowestorder (DQEL) to obtain the nondimensional natural frequencies of tapered Timoshenko beams made up of AFG materials. Ghasemi et al. [8] analyzed the stress and the strain in the pressure vessel made of functionally graded materials reinforced by laminated composite. Rahmani and Pedram [9] investigated the free vibration of FG nanobeams using the Timoshenko beam model and the nonlocal theory. Besides, Gan et al. [10] surveyed the dynamic behavior of nonuniform AFG Timoshenko beams subject to multiple moving point loads. Ebrahimi and Mokhtari [11] applied the Differential Transformation Method (DTM) to investigate vibration responses of rotating FG Timoshenko beams made of porous materials. Via the Generalized Differential Quadrature Method (GDQM), Ghasemi and Mohandes [1215] performed comprehensive investigations on free vibration analysis of laminated composite Timoshenko and EulerBernoulli beams affected by finite strain. The sinusoidal shear deformation theory was employed by Arefi and Zenkour [16] to expand a transient formulation for a threelayer curved nanobeam in thermo–magnetoelastic environments. Additionally, to analyze free vibration analysis of AFG Euler–Bernoulli beams with varying crosssections, Ghazaryan et al. [17] adopted the DTM. Li et al. [18] applied the GDQM to analyze the instability of a microscaled bidirectional functionally graded beam. Arefi and Takhayor Ardestani [19] assessed the electrothermoelastic analysis of functionally graded piezoelectric thickwalled spherical shell utilizing the division method. Based on the firstorder shear deformation theory (FSDT), Mirzaei et al. [2022] perused timedependent creep and thermoelastic analysis of a rotating cantilever tapered beam with axially varying materials subjected to mechanical, inertia and thermal loadings.
The problem of beams resting on an elastic foundation is also a crucial and technical topic in structural and geotechnical engineering. Practical examples of these are railroad structures, highway pavements, the foundations of buildings, and pipelines embedded in the soil. There are various types of foundation models for soilstructure interaction in static and dynamic analysis of structures on an elastic foundation. The Winkler model, which consists of infinitely closed spaced linear translational springs which are independent of each other, is extensively used in the solution of the problems mentioned above related to soilstructure interaction. Several types of research have been conducted on the mechanical behavior of beam elements lying on an elastic foundation.
Also, Zhu and Leung [23] suggested a new finite element formulation for the nonlinear free and forced vibration analysis of nonprismatic Timoshenko beams lying on twoparameter foundations. According to the nonlocal Timoshenko beam theory, stability analysis of nanotubes embedded in an elastic matrix was also performed by Wang et al. [24]. Mirzabeigy [25] investigated the free vibration behavior of nonuniform beams resting on an elastic foundation by presenting a semianalytical technique, based on the DTM. Tsiatas [26] developed a new influential approach to accurately determine stiffness and mass matrices of nonuniform EulerBernoulli beam from inhomogeneous linearly elastic material resting on an elastic foundation. Hassan and Nassar [27] assessed the linear buckling and the free vibration analysis of the Timoshenko beam resting on a twoparameter foundation by employing the ADM. Akgoz and Civalek [28] applied higherorder shear deformation microbeams and a modified strain gradient theory to analyze the static bending response of singlewalled carbon nanotubes embedded in an elastic medium. The stability analysis of the AFG nonprismatic beam on an elastic foundation was comprehensively examined by Shvartsman and Majak [29]. Based on the modified strain gradient theory and surface stress effects, Mohammadimehr et al. [30] exploited the sizedependent effect on the free vibration behavior of Timoshenko microbeams subjected to prestress loading embedded in an elastic medium. Mercan and Civalak [31] analyzed the stability of boron nitride nanotube on the elastic matrix by utilizing a discrete singular convolution technique. By considering the impact of the viscoelastic foundation, Calim [32] studied free and forced vibration of AFG Timoshenko beams. Soltani and Asgarian [33] combined the power series approximation and the RayleighRitz method to assess the free vibration and stability of AFG tapered beam resting on the WinklerPasternak foundation.
Therefore, the main objective of the present paper is to derive highly accurate static and buckling stiffness matrices of axially functionally graded Timoshenko tapered beams subjected to compressive axial concentrated load and supported by a uniform Winkler foundation. For this, a novel numerical methodology based on a general and straightforward procedure presented in [3437] is established. The superiority of the finite element method over the other semianalytical and mathematical techniques is its simplicity, excellent precision, and generality. This methodology is applicable to analyzing a vast range of problems subjected to various circumstances. Considering these facts, the majority of the structural engineering simulation software is commonly developed based on the finite element solution.
An overview of the important contents of this study is given below:
1 Coupled governing equations for the buckling of AFG tapered Timoshenko beams resting on a uniform Winkler foundation have been derived via the energy principle, and they are analytically solved using the power series method. According to the aforementioned method, the expressions of the vertical deflection and crosssectional rotation modes are also determined.
2 Next, the expressions of new shape functions extracted from the beam’s nodal displacements and the principle of virtual work along the beam axis are employed to determine the exact terms of 4*4 static and buckling stiffness matrices. Finally, one can acquire critical buckling loads by solving the eigenvalue problem.
3 To assure the precision and practical usefulness of this hybrid formulation, comparisons with existing results in the literature are provided for a particular case. Subsequently, the effects of the taper ratio, the material nonhomogeneity index, end restraints, and the Winkler parameter on critical buckling loads are investigated.
A significant point of departure of the present finite element solution from others is in the interpolating shape functions used for derivation the structural stiffness matrices. Unlike the Hermitian interpolation polynomials, the expressions of proposed shape functions are dependent on geometrical properties, Winker elastic foundation coefficient, material characteristics, and the compressive axial load. Besides, the accuracy of this method is improved by contemplating the influence of material gradient and varying crosssections in the calculation procedure of the terms of structural and buckling stiffness matrices. This numerical methodology is not restricted by any computational operations. It can be easily used for linear stability analysis of nonprismatic Timoshenko beam with axially varying materials subjected to different boundary conditions. It is also believed that the rate of convergence of the present formulation is faster than that of the conventional finite element technique.
2. Derivation of the Governing Equations
Consider a straight nonuniform beam element of length span L resting on Winkler’s elastic foundation (Fig. 1) and loaded by a constant axial compressive force P applied at both ends. We contemplate the right hand Cartesian coordinate system, with x the initial longitudinal axis measured from the left end of the beam, the yaxis in the lateral direction, and the zaxis along the thickness of the beam. The origin of these axes (O) is located at the centroid of the crosssection. The crosssection is in the form of a rectangle with breadth b and height h, which is assumed to be sufficiently small relative to the breadth. It should be noted that the Timoshenko beam assumptions are adopted here in order to take into consideration the influences of the shear deformation of the beam. Assuming that the deformation of the beam has taken place about the weak axis, in the xz plane, the total potential energy of the considered member can be expressed as [3840]:
(1) 
In Eq. (1), w is the transverse deflection of the centerline of the beam, q the slope due to bending, I the second moment of the area with respect to the yaxis, A the crosssectional area, k the shear correction factor. G and E are the shear and Young's modulus, which are variable along the beam’s length. k_{w} denotes Winkler’s foundation constant per unit length of the beam.
Fig. 1. AFG tapered beam on Winkler’s foundation and subjected to an axial load, Coordinate system and notation of displacement parameters
In the stationary state, the equilibrium equations for nonprismatic Timoshenko beam are derived from a variation of total potential energy which is
(2) 
d illustrates a virtual variation in the last formulation. In the case of the constant axial load (P) and according to Eq. (2) with respect to u_{0}, w and q, the equilibrium equations for a nonprismatic Timoshenko beam are derived as
(3a) 

(3b) 

(3c) 
The last two equilibrium equations (3b and 3c) are coupled differential equations due to the presence of vertical and rotation displacement components (w and θ) as well as shear rigidity (GA), while the axial stability equation (Eq. (3a)) is uncoupled from the others. It has no incidence on linear stability analysis of the Timoshenko beam.
In the following section, the application of the Power Series Method (PSM) in the linear stability analysis of nonhomogeneous Timoshenko beams with nonuniform crosssection is presented. According to this semianalytical method, all variable geometric and material properties of a beam and the displacement components are developed into the power series form.
3. Numerical Approach
In Eq. (3), due to the presence of a functionally graded Timoshenko beam with variable crosssection, the geometrical characteristics of the crosssection over the member length are variable (I(x), A(x)). Moreover, shear and young’s modulus of elasticity are not constant along the xaxis (G(x), E(x)). For these reasons, all the variable terms are presented in power series form, as follows:
(4) 
where , and are coefficients of power series at order i. In order to facilitate the solution of the equations (3b) and (3c), a nondimensional variable ( ) is introduced. Furthermore, the general solutions of two displacement parameters ( , ) should be presented by the following power series of the form:
(5) 
where a_{i} and b_{i} are unknown coefficients and the amplitude of the i component. Substituting Eqs. (4) and (5) and the nondimensional variable e into the system of stability equations, the following expressions are obtained:
(6a) 

(6b) 
in which:
(7) 
After some simplifications, the following expressions are obtained:
(8a) 

(8b) 
or:
(9a) 

(9b) 
To satisfy these equations for every value of ewe must have:
(10a) 

(10b) 
Based on Eqs. (10a) and (10b), the following recurrence formulas are obtained:
(11a) 

(11b) 
According to the above recurrence formulations and from a mathematical point of view, it is culminated that all the a_{k} and b_{k}coefficients can be obtained except for the first two ( and ), which can be derived by imposing the natural boundary conditions of Timoshenko beam. Note that the terms of a_{k+2} and b_{k+2} converge to zero as .
Applying Eqs. (11a) and (11b) together with Eq. (5), the fundamental solutions of the coupled system of differential equations (3b)–(3c) can be thus determined explicitly in terms of the four constants ( and ). Then, the general solution of equilibrium equations can be expressed in the following matrix equation of the form:
(12) 
in which
(13a) 

(13b) 

(13c) 
In equation (12), [B] is a matrix including the fundamental solutions of equilibrium equations for linear stability ( and ), and {A} represents the column vector of four unknown parameters. All terms of and are derived with the aid of the symbolic software MATLAB [44] and the expressions of displacement functions ( and ) are shown in Appendix A.
4. Boundary Conditions and Shape Functions
It has to be noted that the four undefined coefficients ( ) are functions of the displacements of the degree of freedom (DOF), then all the rest of coefficients are also functions of the displacements of DOF. The expression of the angle of rotation and vertical displacement ( and ) can thus be derived as a function of the displacement of DOF. These mentioned unknown parameters can be obtained by imposing the right and left end boundary conditions of the element (two at each end).
In the current finite element model, there are two nodes with two degrees of freedom per node (DOF) for each element. The two nodes by which the finite element can be assembled into the structure are located at its ends. The nodal displacements of the beam element in the local coordinate at and are illustrated in Fig. 2.
Fig. 2. The nodal displacement vectors of a tapered Timoshenko beam element of length L_{e}
According to Fig. 2, the following boundary conditions in the local coordinate must be satisfied:
(14) 
where
(15a) 

(15b) 
The considered degrees of freedom at the left and right nods of each element are: (the transverse displacement in the zdirection), (the angle of rotation). The DOF vector of the element is given in Eq. (15). From (14), one gets
(16) 
Based on Eq. (15a) and knowing that and (see Appendix A), the inverse of matrix [V] is acquired as
(17) 

Subsequently, the system (12) changes into
(18) 
Multiplying the three matrices and expanding Eq. (18), we obtain
(19) 
The terms of the element stiffness matrix can be found from the derivation of the interpolation functions. The shape functions define the deformation shape of the element from applying unit translation or rotation at each of the four degrees of freedom while constraining the other three nodal displacements. With these four interpolation functions, the exact deformed shape of the beam element can be expressed in terms of its nodal displacements. Therefore, the corresponding shape functions must satisfy the following boundary conditions:
(20) 
Referring to these four boundary conditions corresponding to each shape function and using Eq. (19), four sets of interpolation functions can be derived. These shape functions could take any arbitrary shapes which satisfy the boundary conditions of element and internal continuity requirements. In other words, the general displacement of the considered beam is related to four evaluated shape functions for each displacement parameter. Using equations (18) to (20), the general displacement expressions in terms of the shape functions can thus be expressed as
(21) 
where
(22) 
The equivalent shape functions ( ) are then given as
(23a) 

(23b) 

(23c) 

(23d) 
After noticing the resulting interpolating functions (Eq. (23)) and the symbolic expressions presented in Appendix A which are related to the fundamental solutions of the equilibrium equations (Eq. (3)), it can be stated that the new shape functions, on the contrary to Hermitian cubic interpolation (cubic spline) functions, are proportional to the mechanical features of beam element including the geometrical properties, material characteristics, Winkler coefficient, as well as the applied compressive axial load (P). In the remaining part of the current paper, the finite element formulation using the acquired interpolation shape functions in the powers series form (Eq. (23)) is developed.
5. Finite Element Formulation
In section 2, the governing equilibrium equations have already been achieved from the principle of stationary total potential energy. While, in the present part, the finite element formulation is outlined based on the principle of minimum potential energy. The terms of the elemental stiffness matrix of AFG nonprismatic Timoshenko beam in the nondimensional coordinate are carried out by substituting the acquired interpolation shape functions (Eq. (23)) into Eq. (2)
(24a) 

(24b) 

(24c) 

(24d) 
where , and are respectively, the terms of the flexural stiffness matrix, shear stiffness matrix, foundation stiffness matrix, and the geometric stiffness matrix, which account for the secondary effect of the axial force (P) on the elemental stiffness matrix. L_{e} is also the length of each segment.
Afterwards, by assembling each element stiffness matrices based on its nodal displacements, the stiffness matrix of the whole structure can be achieved. In most finite element method textbooks [41, 42], one can find the description of the process of assemblage in detail. It should be pointed out that the linear analysis is under consideration in the present study. In the linear buckling analysis, the geometric stiffness matrix is proportional to the initial stress forces. In order to perform stability analysis and evaluate the critical buckling loads, the following eigenvalue problem should be solved:
(25) 
in which and are respectively the eigenvalues and their related eigenvectors, which associate with the total Degree of Freedom of each element. It is well known that for a system with n Degrees of Freedom, there exist n buckling modes, but in practice, only the lowest ones are of interest.
In order to clarify the numerical procedure, a general algorithm is illustrated in Fig. 3 for determination of fundamental solutions, interpolating shape functions, and structural matrices for AFG nonprismatic Timoshenko beams on Winkler elastic foundation. The presented algorithm is used for computer applications of the method.
6. Numerical Examples
This section aims to measure the accuracy and check the validity of the present numerical procedure in the buckling analysis of AFG beams with variable crosssections resting on an elastic foundation. In order to achieve this goal, two illustrative examples are carried out. The influences of end conditions, material gradient index, slenderness ratio,elastic foundation modulus, and tapering parameter on the buckling loads of beam are examined.
Fig. 3. The algorithm used for stability analysis of AFG nonprismatic Timoshenko beams resting on elastic foundation
Through this section, linear buckling analysis is performed for a double tapered beam with rectangular crosssection whose height (h_{0}) and breadth (b_{0}) both are concurrently allowed to vary linearly along the member’s length with the same tapering ratio. Therefore, the crosssectional area A(x) and the area moment of inertia I(x) vary along the beam as
(26) 
where b is breadth and height taper ratios. Note that the tapering parameter can change from b=0 (prismatic beam) to b=0.10.9 (nonuniform ones). In Eq. (26), A_{0} and I_{0} are respectively crosssectional area and moment of inertia at the left support (x=0). They are defined as: and .
In the following benchmark examples, it is also 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; Al: E_{1}=70GPa. The variation of Young’s modulus of elasticity along the beam axis is defined with the following powerlaw formulation:
(27) 
In the last expression, m signifies the material nonhomogeneity parameter. By notifying this formulation, it can be stated that by descending the gradient index (m), the proportion of zirconia over the beam’s length increases. It should be noted that the material properties of the beam are assumed to be constant in the direction of the thickness. Poisson’s ratio of the material also remains constant in a longitudinal direction. Further, Poisson’s ratio and the shear correction factor are assumed 0.3 and 5/6, respectively. In the numerical computation, the nondimensional forms of buckling load and elastic foundation parameter are introduced as

(28a, b) 
6.1 Example 1: Tapered beam from functionally graded materials
To study the effect of the Winkler foundation on the buckling capacity, this comprehensive example is considered including axially nonhomogeneous and homogeneous beams with nonuniform crosssection. The following parameters are used: length L =1m, the ratio between moment of inertia, and crosssectional area I_{0}/A_{0} = 0.01.
At first, the validation of the present formulation for buckling analysis of AFG tapered Timoshenko beam without elastic foundation is checked by comparing the obtained results with those available in the literature when possible.
Note that it is an important task to evaluate the convergence characteristics of any numerical methods, especially the finite element solution to guarantee the exactness and performance of the adopted technique to different engineering problems. Accordingly, estimating the required number of segments in the assemblage of beam element stiffness matrices to determine acceptable buckling results is essential before the presentation of numerical results of buckling loading of AFG Timoshenko columns with varying crosssection.
According to the authors’ knowledge of power series approximation [3337], for most nonprismatic members under various loading conditions, it is not practicle to use more than 10 terms in power series expansions to derive the exact shape functions. Therefore, the taken number in terms of the power series (the maximum power of shape functions) is equal to 10 in the following cases.
In this regard, the lowest value of the nondimensional buckling load of cantilever tapered beam with two different taper ratios (b = 0.2 and 0.5) is acquired with respect to the number of meshes (n) considered in the procedure for the assemblage of the stiffness matrices of the whole structure. The convergence study is carried out for homogenous beams, as well as axially functionally graded ones by contemplating the gradient parameter (m) equals to two. The obtained results by the proposed numerical technique compared with those acquired via new finite element modeling introduced by Soltani and Asgarian [36]. The percentage of relative errors (D (%)) between the results of the present study and the mentioned numerical method is also calculated by the following expression:
(29) 
In the following, the graphic illustrations of the variation of relative errors with the number of considered segments (n) in the FE discretization phase are provided in Fig. 4.
After noticing the results represented in Fig. 4, the following outcomes could be expressed:
Fig. 4: Timoshenko beams with variable crosssection: variation of the relative errors (D) versus the number of elements (n) along the beam’s length
1 An outstanding compatibility between the elastic buckling loads acquired by the current study and those computed from another available benchmark solution is noteworthy
2 Even by applying 4 segments in the beam’s length according to the suggested finite element solution, the elastic buckling loads can be precisely calculated bellow the acceptable relative error (1%).
3 It is not required to use more than 6 segments in the finite element approach, in order to obtain a satisfactory accuracy on critical elastic buckling loads.
Following the procedure mentioned above, the first nondimensional buckling load parameters for tapered Timoshenko beam from homogenous materials and axially functionally ones with different gradient indexes (m=1, 2, and 3) are derived using the proposed finite element solution by dividing the beam into 6 equal segments. The dimensionless buckling loads of the member for various tapering ratios and different boundary conditions are arranged in Table 1 and compared with Soltani et al. [43]. As presented in Table 1, an excellent agreement is observed between the critical buckling loads for different values of nonuniformity ratios acquired by the present study and those computed from the other benchmark solutions.
Comparing the results of the prismatic beams presented in Table 1 with those related to nonprismatic ones, it can be culminated that the considered prismatic beam in this example had the highest critical buckling loads. This can be explained by the fact that increasing taper ratios causes the reduction in crosssectional area and moment of inertia and consequently stiffness of the elastic member. Since the linear stability behavior is directly proportional to the member’s stiffness, an increase in tapering ratio causes a decrement in the dimensionless axial buckling load. Furthermore, it can be argued from Table 1 that the increase in the gradientindices leads to an increase in the critical buckling loads. For example, in the case of the simply supported FG prismatic beam ( ), the critical buckling load increases from 4.6618 to 5.5390 and then to 6.0492, when m increases from 1 to 3. It shows a rise of 18.82% and 29.76%, accordingly. Regarding Eq. (27), it is clear that when the nonhomogeneity index is increased from 1 to 3, the volume fraction of Zirconia and, consequently, the value of Young’s modulus are increased. As a result, a higher buckling load is achieved.
Afterward, the lowest buckling loads variation versus the taper ratio (b) and the gradient index (m) for simply supported and fixedfree beams are presented in Fig. 5. Each of the depictions of Fig. 5a and b present six different plots relating to m = 0.5, 1, 1.5, 2, 2.5, and 3.
Table 1. Critical axial load parameter ( ) for tapered Timoshenko beams with different material nonhomogeneity indexes (m)
End Conditions 
(b) 
Homogenous 
m=1 
m=2 
m=3 

Present method 
Soltani et al. [43] 
Present method 
Soltani et al. [43] 
Present method 
Soltani et al. [43] 
Present method 
Soltani et al. [43] 

HingedHinged 
0.0 
7.546 
4.546 
4.662 
4.667 
5.539 
5.556 
6.049 
6.075 
0.1 
6.274 
6.248 
3.756 
3.749 
4.484 
4.467 
4.897 
4.902 

0.2 
5.055 
5.036 
2.930 
2.926 
3.501 
3.487 
3.833 
3.839 

0.3 
3.936 
3.924 
2.205 
2.203 
2.634 
2.623 
2.890 
2.896 

0.4 
2.934 
2.928 
1.584 
1.584 
1.879 
1.882 
2.076 
2.082 

0.5 
2.062 
2.062 
1.070 
1.071 
1.264 
1.267 
1.398 
1.405 

0.6 
1.333 
1.337 
0.661 
0.663 
0.776 
0.780 
0.858 
0.865 

0.7 
0.753 
0.765 
0.356 
0.359 
0.414 
0.418 
0.457 
0.463 

0.8 
0.341 
0.350 
0.150 
0.154 
0.171 
0.176 
0.188 
0.194 

0.9 
0.087 
0.085 
0.035 
0.032 
0.039 
0.036 
0.042 
0.039 

FixedFree 
0.0 
2.291 
2.291 
1.710 
1.712 
1.982 
1.986 
2.097 
2.104 
0.1 
2.016 
2.016 
1.469 
1.471 
1.712 
1.716 
1.821 
1.827 

0.2 
1.743 
1.743 
1.233 
1.235 
1.446 
1.450 
1.547 
1.553 

0.3 
1.472 
1.472 
1.006 
1.008 
1.186 
1.190 
1.277 
1.283 

0.4 
1.204 
1.205 
0.789 
0.791 
0.935 
0.939 
1.013 
1.020 

0.5 
0.942 
0.944 
0.586 
0.588 
0.696 
0.701 
0.760 
0.767 

0.6 
0.688 
0.691 
0.401 
0.404 
0.476 
0.481 
0.524 
0.530 

0.7 
0.449 
0.453 
0.241 
0.244 
0.285 
0.289 
0.314 
0.320 

0.8 
0.236 
0.241 
0.114 
0.117 
0.133 
0.136 
0.146 
0.151 

0.9 
0.072 
0.075 
0.031 
0.031 
0.034 
0.036 
0.037 
0.035 
As shown in Fig. 5, for any value of the powerlaw exponent, the stability of prismatic beam (b=0) and tapered beam with b=0.9 is most and least, respectively. It is observable that increasing the gradient index leads to the enlargement of the dimensionless buckling load for all values of the tapering ratio. The reason is the higher portion of the ZrO2 phase as the value of the gradient index rises. Also, it is found that for , the nondimensional critical loads increase sharply whereas, for m > 1.5, the buckling resistance increases slightly and approaches maximum magnitude.
In the following, the linear stability problem for simply supported and clampedfree beams made up of axially functionally graded materials in the presence of the Winkler foundation is investigated. For this purpose, the nondimensional buckling load is carried out for four different taper ratios (b), namely 0, 0.2, 0.5, and 0.8; where the first one is a prismatic beam, while the others are tapered ones. Note that the FEM results are obtained by discretizing the beam into 6 elements.
The variation of the lowest buckling load parameters for hingedhinged and fixedfree Timoshenko beams resting on the Winkler type of foundation versus the elastic foundation constant ( ) are respectively presented in Figs. 6 and 7.
As can be seen, the variation of the Winkler elastic foundation parameter has a significant influence on the linear stability behavior of both beams under different circumstances. It is found from these figures that the critical buckling load parameters corresponding to the first mode are increased as the stiffness of the elastic foundation increases. In other words, the numerical outcomes show that the elastic foundation has a stabilizing effect on the stability characteristics of axially nonhomogeneous and homogeneous Timoshenko beams with constant or variable crosssection. By pondering Fig. 6, one can remark that for homogeneous prismatic simply supported beam, the buckling load increases linearly due to the effect of uniform Winkler foundation. At the same time, this variation for doubletapered members is nonlinear.
It can be stated that the variation of stability behavior for the cantilever is similar to those for the hingedhinged beam, but the latter beam appears to be more unstable than the former. Moreover, for any value of the Winkler’s parameters, the corresponding buckling load for the beam having constant material properties and a uniform crosssection is the highest and that for tapered beams having a material property according to powerlaw with index m = 1 is the lowest.
Fig. 5. Effects of the gradient index (m) on the normalized buckling load of tapered Timoshenko beam with different tapering ratios: (a) simply supported, (b) fixedfree
Fig. 6. Influence of Winkler type elastic foundation modulus on the critical axial load parameters of simply supported Timoshenko beams: (a) , (b) , (c) , (d)
Fig. 7. Influence of Winkler type elastic foundation modulus on the critical axial load parameters of fixedfree Timoshenko beams: (a) , (b) , (c) , (d) .
6.2 Example 2: Tapered Timoshenko beam versus tapered EulerBernoulli beam
To depict the difference between the critical buckling loads of the EulerBernoulli beam theory and those of the Timoshenko beam model, this illustrative example is also considered. In this regard, for different values of slenderness ratio (L/b_{0}), the lowest critical loads are calculated for hingedhinged beam, and the results are listed in Table 2. The critical axial loads are also carried out for two cases: axially nonhomogeneous and homogeneous beams. In the case of axially FG members, the distribution of modulus of elasticity is contemplated to vary in the longitudinal direction with a powerlaw formulation as expressed in Eq. (27). In this case, the material nonhomogeneity parameter (m) is assumed to be equal to 1. Moreover, the buckling load is acquired for four different tapering parameters: b=0, 0.2, 0.5, and 0.8.
It should be noted here that the static stability responses of EBT are independent of the aspect ratio. Also, it is worth mentioning that the values of the taper parameter and the type of boundary conditions are selected in such a way that make comparison possible with available reference in the case of the tapered EulerBernoulli beam [33].
From Table 2, one can see that the elastic buckling loads calculated by employing the EulerBernoulli theory are overestimated for all values of thickness ratio (L/b_{0}). The difference between the EBT and shear deformation theory (Timoshenko beam model) is significant for deep members (L/b_{0}<10). Note that for slender and long beams (L/b_{0}>50), this difference is small and negligible (see Table 2). According to EBT, the influences of shear deformation and rotatory inertia are ignorable, and only the effect of flexural deformation is taken into account. In general, the effect of transverse shear deformation is to increase the deflection, which leads to a noticeable decrease in the value of the stiffness and rigidity of the member, and consequently, a weaker member being obtained. Therefore, a significant decrease in the elastic buckling loads of the member is observed. The shear deformation loses its effect on the transverse deflections, as the value of the thickness ratio increases. It is also clear that the exclusion of the shear deformation decreases the deflections and increases the buckling loads. Moreover, it is reasonable that the elastic critical load increases noticeably, when considering the Winkler foundation effects. This fact can be easily observed in Table 2.
7. Conclusions
In this paper, a hybrid numerical methodology to acquire the axial critical buckling loads of nonprismatic functionally graded beams based on the Timoshenko theory was proposed. This new technique combined the power series method and finite element solution and applied them to tapered beams lying on a uniform Winkler foundation. In this regard, the power series expansions method was used to solve the system of secondorder differential equations with variable coefficients and determine the exact four sets shape functions for the nonuniform element with four degrees of freedom. In turn, based on the principle of internal virtual work along the element axis, the exact element stiffness matrices for the buckling analysis of the nonuniform Timoshenko beam resting on an elastic foundation were established.
Table 2. The dimensionless buckling load of simply supported Timoshenko beam and EulerBernoulli beam
Material 
b 
L/b_{0} 
Timoshenko Beam Theory 
EulerBernoulli Beam Theory Soltani and Asgarian [33] 

Homogenous

0.0 
5 
8.9687 
13.0088 
17.0567 
9.8694 
13.9236 
17.9778 
10 
9.6289 
13.6775 
17.7285 

50 
9.8607 
13.9127 
17.9652 

100 
9.8681 
13.9203 
17.9728 

0.2 
5 
5.8342 
9.8243 
13.7806 
6.3328 
10.3542 
14.3598 

10 
6.1894 
10.1952 
14.1805 

50 
6.3120 
10.3228 
14.3151 

100 
6.3159 
10.3269 
14.3194 

0.5 
5 
2.3178 
5.7425 
8.0535 
2.4891 
6.0834 
9.0018 

10 
2.4288 
5.9595 
8.6897 

50 
2.4664 
6.0269 
8.8817 

100 
2.4676 
6.0291 
8.8876 

0.8 
5 
0.3761 
1.4630 
1.8132 
0.4189 
1.5901 
2.0081 

10 
0.3911 
1.5449 
1.9263 

50 
0.3963 
1.5729 
1.9647 

100 
0.3973 
1.5742 
1.9659 

m= 1 
0.0 
5 
5.7001 
9.6451 
13.4935 
6.3961 
10.3999 
14.3801 
10 
6.1959 
10.1828 
14.1351 

50 
6.3715 
10.3696 
14.3422 

100 
6.3772 
10.3755 
14.3487 

0.2 
5 
3.5125 
7.2245 
10.3464 
3.9030 
7.7599 
11.4402 

10 
3.7852 
7.5997 
11.1608 

50 
3.8804 
7.7214 
11.3714 

100 
3.8834 
7.7252 
11.3777 

0.5 
5 
1.2538 
3.7693 
4.7225 
1.3790 
4.3439 
5.8474 

10 
1.3369 
4.1660 
5.4850 

50 
1.3654 
4.2892 
5.7569 

100 
1.3663 
4.2930 
5.7656 

0.8 
5 
0.1704 
0.7620 
0.9490 
0.1941 
0.8995 
1.1354 

10 
0.1794 
0.8272 
1.0411 

50 
0.1825 
0.8500 
1.0734 

100 
0.1827 
0.8508 
1.0746 
The applicability and efficiency of the proposed finite element formulation in linear stability analysis of tapered AFG beams on the uniform elastic foundation were demonstrated by providing two illustrative examples in which the effects of the Winkler parameter, mechanical variation, and different boundary conditions were investigated. The acquired outcomes were contrasted with other analytical and numerical solutions presented in the literature. In most cases, it was concluded that by considering up to 10 terms of power series and only 6 elements in mesh process, the buckling loads related to the lowest buckling loads of axially functionally graded Timoshenko beams resting on the elastic foundation and having variable crosssections could be determined with high competency. Moreover, the numerical outcomes illustrated that Winkler’s elastic foundation had a stabilizing effect on the stability characteristics of axially nonhomogeneous and homogeneous beams with different tapering ratios. It was also observed that the results obtained based on the EulerBernoulli model were greater than those predicted by the shear deformation beam model.
Appendix A
By using symbolic MATLAB software [44], the displacements functions and are derived. The first few terms are expressed in the following forms:
(A.1) 


(A.2) 

(A.3) 



(A.4) 

(A.5) 

(A.6) 




(A.7) 


(A.8) 





References
[1] Şimşek, M., 2010. Fundamental frequency analysis of functionally graded beams by using different higherorder beam theories. Nuclear Engineering and Design, 240(4), pp.697705.
[2] Shahba, A., Attarnejad, R., Marvi, M.T., and Hajilar, S., 2011. Free vibration and stability analysis of axially functionally graded tapered Timoshenko beams with classical and nonclassical boundary conditions. Composites Part B: Engineering, 42(4), pp.801808.
[3] Akgoz, B. and Civalek, O., 2013. Buckling analysis of linearly tapered microcolumns based on strain gradient elasticity. Structural Engineering and Mechanics, 48(2), pp.195205.
[4] Arefi, M. and Rahimi, G.H., 2013. Nonlinear analysis of a functionally graded beam with variable thickness. Scientific Research and Essays, 8(6), pp.256264.
[5] Malekzadeh, P. and Shojaee, M., 2013. Surface and nonlocal effects on the nonlinear free vibration of nonuniform nanobeams. Composites Part B: Engineering, 52, pp.8492.
[6] LezgyNazargah, M., Vidal, P. and Polit, O., 2013. An efficient finite element model for static and dynamic analyses of functionally graded piezoelectric beams. Composite Structures, 104, pp.7184.
[7] Rajasekaran, S. and Tochaei, E.N., 2014. Free vibration analysis of axially functionally graded tapered Timoshenko beams using differential transformation element method and differential quadrature element method of lowestorder. Meccanica, 49(4), pp.9951009.
[8] Ghasemi, A.R., Kazemian, A. and Moradi, M., 2014. Analytical and numerical investigation of FGM pressure vessel reinforced by laminated composite materials.
[9] Rahmani, O. and Pedram, O., 2014. Analysis and modeling the size effect on vibration of functionally graded nanobeams based on nonlocal Timoshenko beam theory. International Journal of Engineering Science, 77, pp.5570.
[10] Gan, B.S., Trinh, T.H., Le, T.H. and Nguyen, D.K., 2015. Dynamic response of nonuniform Timoshenko beams made of axially FGM subjected to multiple moving point loads. Structural Engineering and Mechanics, 53(5), pp.981995.
[11] Ebrahimi, F. and Mokhtari, M., 2015. Transverse vibration analysis of rotating porous beam with functionally graded microstructure using the differential transform method. Journal of the Brazilian Society of Mechanical Sciences and Engineering, 37(4), pp.14351444.
[12] Ghasemi, A.R. and Mohandes, M., 2016. The effect of finite strain on the nonlinear free vibration of a unidirectional composite Timoshenko beam using GDQM. Advances in aircraft and spacecraft science, 3(4), p.379.
[13] Mohandes, M. and Ghasemi, A.R., 2016. Finite strain analysis of nonlinear vibrations of symmetric laminated composite Timoshenko beams using generalized differential quadrature method. Journal of Vibration and Control, 22(4), pp.940954.
[14] Ghasemi, A.R. and Mohandes, M., 2017. Nonlinear free vibration of laminated composite EulerBernoulli beams based on finite strain using generalized differential quadrature method. Mechanics of Advanced Materials and Structures, 24(11), pp.917923.
[15] Mohandes, M. and Ghasemi, A.R., 2017. Modified couple stress theory and finite strain assumption for nonlinear free vibration and bending of micro/nanolaminated composite Euler–Bernoulli beam under thermal loading. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 231(21), pp.40444056.
[16] Arefi, M. and Zenkour, A.M., 2017. Transient sinusoidal shear deformation formulation of a sizedependent threelayer piezomagnetic curved nanobeam. Acta Mechanica, 228(10), pp.36573674.
[17] Ghazaryan, D., Burlayenko, V.N., Avetisyan, A. and Bhaskar, A., 2018. Free vibration analysis of functionally graded beams with nonuniform crosssection using the differential transform method. Journal of Engineering Mathematics, 110(1), pp.97121.
[18] Li, X., Li, L.I. and Hu, Y., 2018. Instability of functionally graded microbeams via microstructuredependent beam theory. Applied Mathematics and Mechanics, 39(7), pp.923952.
[19] Arefi, M. and Takhayor Ardestani, M., 2019. The Effect of Temperature Dependency on the ThermoElectroElastic Analysis of Functionally Graded Piezoelectric Spherical Shell. Mechanics of Advanced Composite Structures, 6(1), pp.18.
[20] Mirzaei, M.M.H., Loghman, A. and Arefi, M., 2019. Timedependent creep analysis of a functionally graded beam with trapezoidal cross section using firstorder shear deformation theory. Steel and Composite Structures, 30(6), pp.567576.
[21] Mirzaei, M.M.H., Loghman, A. and Arefi, M., 2019. Effect of Temperature Dependency on Thermoelastic Behavior of Rotating Variable Thickness FGM Cantilever Beam. Journal of Solid Mechanics, 11(3), pp.657669.
[22] Arefi, M., Loghman, A. and Mohammad Hosseini Mirzaei, M., 2019. Thermoelastic analysis of a functionally graded simple blade using firstorder shear deformation theory. Mechanics of Advanced Composite Structures.
[23] Zhu, B. and Leung, A.Y.T., 2009. Linear and nonlinear vibration of nonuniform beams on twoparameter foundations using pelements. Computers and Geotechnics, 36(5), pp.743750.
[24] Wang, B.L., Hoffman, M. and Yu, A.B., 2012. Buckling analysis of embedded nanotubes using gradient continuum theory. Mechanics of Materials, 45, pp.5260.
[25] Mirzabeigy, A., 2014. Semianalytical approach for free vibration analysis of variable crosssection beams resting on elastic foundation and under axial force. International Journal of Engineering, Transactions C: Aspects, 27(3), pp.455463.
[26] Tsiatas, G.C., 2014. A new efficient method to evaluate exact stiffness and mass matrices of nonuniform beams resting on an elastic foundation. Archive of Applied Mechanics, 84(5), pp.615623.
[27] Hassan, M.T. and Nassar, M., 2015. Analysis of stressed Timoshenko beams on two parameter foundations. KSCE Journal of Civil Engineering, 19(1), pp.173179.
[28] Akgöz, B. and Civalek, Ö., 2016. Bending analysis of embedded carbon nanotubes resting on an elastic foundation using strain gradient theory. Acta Astronautica, 119, pp.112.
[29] Shvartsman B. and Majak J. Numerical method for stability analysis of functionally graded beams on elastic foundation. Applied Mathematical Modelling, 2015; 44: 3713–3719.
[30] Mohammadimehr, M., Mohammadi Hooyeh, H., Afshari, H. and Salarkia, M.R., 2016. Sizedependent Effects on the Vibration Behavior of a Timoshenko Microbeam subjected to Prestress Loading based on DQM. Mechanics of Advanced Composite Structures, 3(2), pp.99112.
[31] Mercan, K. and Civalek, Ö., 2016. DSC method for buckling analysis of boron nitride nanotube (BNNT) surrounded by an elastic matrix. Composite Structures, 143, pp.300309.
[32] Calim, F.F., 2016. Free and forced vibration analysis of axially functionally graded Timoshenko beams on twoparameter viscoelastic foundation. Composites Part B: Engineering, 103, pp.98112.
[33] Soltani, M. and Asgarian, B., 2019. New hybrid approach for free vibration and stability analyses of axially functionally graded EulerBernoulli beams with variable crosssection resting on uniform WinklerPasternak foundation. Latin American Journal of Solids and Structures, 16(3).
[34] Soltani, M., Asgarian, B. and Mohri, F., 2014. Finite element method for stability and free vibration analyses of nonprismatic thinwalled beams. ThinWalled Structures, 82, pp.245261.
[35] Soltani, M. and Mohri, F., 2016. Stability and vibration analyses of tapered columns resting on one or twoparameter elastic foundations. Journal of Numerical Methods in Civil Engineering, 1(2), pp.5766.
[36] Soltani, M. and Asgarian, B., 2019. Finite Element Formulation for Linear Stability Analysis of Axially Functionally Graded Nonprismatic Timoshenko Beam. International Journal of Structural Stability and dynamics, 19(02), p.1950002.
[37] Soltani M, Asgarian B, Mohri F. Improved finite element model for lateral stability analysis of axially functionally graded nonprismatic Ibeams. International Journal of Structural Stability and dynamics, 2019; 19(9): 30 pages.
[38] Arefi, M. and Zenkour, A.M., 2016. A simplified shear and normal deformations nonlocal theory for bending of functionally graded piezomagnetic sandwich nanobeams in magnetothermoelectric environment. Journal of Sandwich Structures & Materials, 18(5), pp.624651.
[39] Arefi, M. and Zenkour, A.M., 2017. Vibration and bending analysis of a sandwich microbeam with two integrated piezomagnetic facesheets. Composite Structures, 159, pp.479490.
[40] Arefi, M. and Zenkour, A.M., 2017. Sizedependent vibration and bending analyses of the piezomagnetic threelayer nanobeams. Applied Physics A, 123(3), p.202.
[41] Zienkiewicz, O.C. and Taylor, R.L., 2005. The finite element method for solid and structural mechanics. Elsevier.
[42] Logan, D.L., 2011. A first course in the finite element method. Cengage Learning.
[43] Soltani M, Asgarian B, Jafarzadeh F. Finite difference method for buckling analysis of tapered Timoshenko beam made of functionally graded material. AUT Journal of Civil Engineering, 2019; DOI: 10.22060/AJCE.2019.15195.5525.
[44] MATLAB Version 7.6. MathWorks Inc, USA, 2008.