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

1.     Introduction

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 Euler-Bernoulli and the Timoshenko beam theories. Based on the Euler-Bernoulli 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 length-to-depth 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 first-order and different higher-order shear deformation theories. The free vibration and linear buckling analyses of non-prismatic 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 Rayleigh-Ritz method, the linear buckling resistance of tapered micro-columns 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 non-local beams with variable cross-sections 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 Lezgy-Nazargah et al. [6]. Rajasekaran and Tochaei [7] applied the Differential Transformation Element Method (DTEM) and Differential Quadrature Element of Lowest-order (DQEL) to obtain the non-dimensional 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 non-uniform 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 [12-15] performed comprehensive investigations on free vibration analysis of laminated composite Timoshenko and Euler-Bernoulli beams affected by finite strain. The sinusoidal shear deformation theory was employed by Arefi and Zenkour [16] to expand a transient formulation for a three-layer curved nanobeam in thermo–magneto-elastic environments. Additionally, to analyze free vibration analysis of AFG Euler–Bernoulli beams with varying cross-sections, Ghazaryan et al. [17] adopted the DTM. Li et al. [18] applied the GDQM to analyze the instability of a micro-scaled bi-directional functionally graded beam. Arefi and Takhayor Ardestani [19] assessed the electro-thermo-elastic analysis of functionally graded piezoelectric thick-walled spherical shell utilizing the division method. Based on the first-order shear deformation theory (FSDT), Mirzaei et al. [20-22] perused time-dependent creep and thermo-elastic 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 soil-structure 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 soil-structure 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 non-linear free and forced vibration analysis of non-prismatic Timoshenko beams lying on two-parameter 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 non-uniform beams resting on an elastic foundation by presenting a semi-analytical technique, based on the DTM. Tsiatas [26] developed a new influential approach to accurately determine stiffness and mass matrices of non-uniform Euler-Bernoulli 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 two-parameter foundation by employing the ADM. Akgoz and Civalek [28] applied higher-order shear deformation microbeams and a modified strain gradient theory to analyze the static bending response of single-walled carbon nanotubes embedded in an elastic medium. The stability analysis of the AFG non-prismatic 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 size-dependent effect on the free vibration behavior of Timoshenko microbeams subjected to pre-stress 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 Rayleigh-Ritz method to assess the free vibration and stability of AFG tapered beam resting on the Winkler-Pasternak 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 [34-37] is established. The superiority of the finite element method over the other semi-analytical 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 cross-sectional 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 non-homogeneity 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 cross-sections 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 non-prismatic 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 non-uniform 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 y-axis in the lateral direction, and the z-axis along the thickness of the beam. The origin of these axes (O) is located at the centroid of the cross-section. The cross-section 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 x-z plane, the total potential energy of the considered member can be expressed as [38-40]:

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 y-axis, A the cross-sectional 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 non-prismatic Timoshenko beam are derived from a variation of total potential energy which is

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₀, w and q, the equilibrium equations for a non-prismatic Timoshenko beam are derived as

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 non-homogeneous Timoshenko beams with non-uniform cross-section is presented. According to this semi-analytical 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 cross-section, the geometrical characteristics of the cross-section over the member length are variable (I(x), A(x)). Moreover, shear and young’s modulus of elasticity are not constant along the x-axis (G(x), E(x)). For these reasons, all the variable terms are presented in power series form, as follows:

where   ,  and  are coefficients of power series at order i. In order to facilitate the solution of the equations (3b) and (3c), a non-dimensional variable ( ) is introduced. Furthermore, the general solutions of two displacement parameters ( , ) should be presented by the following power series of the form:

where a_i and b_i are unknown coefficients and the amplitude of the i component. Substituting Eqs. (4) and (5) and the non-dimensional variable e into the system of stability equations, the following expressions are obtained:

in which:

After some simplifications, the following expressions are obtained:

or:

To satisfy these equations for every value of ewe must have:

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

According to the above recurrence formulations and from a mathematical point of view, it is culminated that all the a_k and b_kcoefficients 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:

 in which

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 co-ordinate must be satisfied:

where

The considered degrees of freedom at the left and right nods of each element are:  (the transverse displacement in the z-direction),  (the angle of rotation). The DOF vector of the element is given in Eq. (15). From (14), one gets

Based on Eq. (15a) and knowing that  and (see Appendix A), the inverse of matrix [V] is acquired as

Subsequently, the system (12) changes into

Multiplying the three matrices and expanding Eq. (18), we obtain

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:

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

where

The equivalent shape functions ( ) are then given as

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 non-prismatic Timoshenko beam in the non-dimensional coordinate are carried out by substituting the acquired interpolation shape functions (Eq. (23)) into Eq. (2)

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:

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 non-prismatic 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 cross-sections 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 non-prismatic Timoshenko beams resting on elastic foundation

Through this section, linear buckling analysis is performed for a double tapered beam with rectangular cross-section whose height (h₀) and breadth (b₀) both are concurrently allowed to vary linearly along the member’s length with the same tapering ratio. Therefore, the cross-sectional area A(x) and the area moment of inertia I(x) vary along the beam as

where b is breadth and height taper ratios. Note that the tapering parameter can change from b=0 (prismatic beam) to b=0.1-0.9 (non-uniform ones). In Eq. (26), A₀ and I₀ are respectively cross-sectional 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₀=200GPa; Al: E₁=70GPa. The variation of Young’s modulus of elasticity along the beam axis is defined with the following power-law formulation:

In the last expression, m signifies the material non-homogeneity 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 non-dimensional forms of buckling load and elastic foundation parameter are introduced as

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 non-homogeneous and homogeneous beams with non-uniform cross-section. The following parameters are used: length L =1m, the ratio between moment of inertia, and cross-sectional area I₀/A₀ = 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 cross-section.

According to the authors’ knowledge of power series approximation [33-37], for most non-prismatic 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 non-dimensional 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:

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 cross-section: 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 non-dimensional 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 non-uniformity 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 non-prismatic 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 cross-sectional 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 non-homogeneity 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 fixed-free 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 non-homogeneity indexes (m)

As shown in Fig. 5, for any value of the power-law 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 non-dimensional 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 clamped-free beams made up of axially functionally graded materials in the presence of the Winkler foundation is investigated. For this purpose, the non-dimensional 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 hinged-hinged and fixed-free 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 non-homogeneous and homogeneous Timoshenko beams with constant or variable cross-section. 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 double-tapered members is non-linear.

It can be stated that the variation of stability behavior for the cantilever is similar to those for the hinged-hinged 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 cross-section is the highest and that for tapered beams having a material property according to power-law 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) fixed-free

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 fixed-free Timoshenko beams: (a) , (b) , (c) , (d) .

6.2 Example 2: Tapered Timoshenko beam versus tapered Euler-Bernoulli beam

To depict the difference between the critical buckling loads of the Euler-Bernoulli 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₀), the lowest critical loads are calculated for hinged-hinged beam, and the results are listed in Table 2. The critical axial loads are also carried out for two cases: axially non-homogeneous 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 power-law formulation as expressed in Eq. (27). In this case, the material non-homogeneity 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 Euler-Bernoulli beam [33].

From Table 2, one can see that the elastic buckling loads calculated by employing the Euler-Bernoulli theory are overestimated for all values of thickness ratio (L/b₀). The difference between the EBT and shear deformation theory (Timoshenko beam model) is significant for deep members (L/b₀<10). Note that for slender and long beams (L/b₀>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 non-prismatic 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 second-order differential equations with variable coefficients and determine the exact four sets shape functions for the non-uniform 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 non-uniform Timoshenko beam resting on an elastic foundation were established.

Table 2. The dimensionless buckling load of simply supported Timoshenko beam and Euler-Bernoulli beam

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 cross-sections 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 non-homogeneous and homogeneous beams with different tapering ratios. It was also observed that the results obtained based on the Euler-Bernoulli 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:

References

[1] Şimşek, M., 2010. Fundamental frequency analysis of functionally graded beams by using different higher-order beam theories. Nuclear Engineering and Design, 240(4), pp.697-705.

[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 non-classical boundary conditions. Composites Part B: Engineering, 42(4), pp.801-808.

[3] Akgoz, B. and Civalek, O., 2013. Buckling analysis of linearly tapered micro-columns based on strain gradient elasticity. Structural Engineering and Mechanics, 48(2), pp.195-205.

[4] Arefi, M. and Rahimi, G.H., 2013. Non-linear analysis of a functionally graded beam with variable thickness. Scientific Research and Essays, 8(6), pp.256-264.

[5] Malekzadeh, P. and Shojaee, M., 2013. Surface and nonlocal effects on the nonlinear free vibration of non-uniform nanobeams. Composites Part B: Engineering, 52, pp.84-92.

[6] Lezgy-Nazargah, 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.71-84.

[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 lowest-order. Meccanica, 49(4), pp.995-1009.

[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.55-70.

[10] Gan, B.S., Trinh, T.H., Le, T.H. and Nguyen, D.K., 2015. Dynamic response of non-uniform Timoshenko beams made of axially FGM subjected to multiple moving point loads. Structural Engineering and Mechanics, 53(5), pp.981-995.

[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.1435-1444.

[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.940-954.

[14] Ghasemi, A.R. and Mohandes, M., 2017. Nonlinear free vibration of laminated composite Euler-Bernoulli beams based on finite strain using generalized differential quadrature method. Mechanics of Advanced Materials and Structures, 24(11), pp.917-923.

[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.4044-4056.

[16] Arefi, M. and Zenkour, A.M., 2017. Transient sinusoidal shear deformation formulation of a size-dependent three-layer piezo-magnetic curved nanobeam. Acta Mechanica, 228(10), pp.3657-3674.

[17] Ghazaryan, D., Burlayenko, V.N., Avetisyan, A. and Bhaskar, A., 2018. Free vibration analysis of functionally graded beams with non-uniform cross-section using the differential transform method. Journal of Engineering Mathematics, 110(1), pp.97-121.

[18] Li, X., Li, L.I. and Hu, Y., 2018. Instability of functionally graded micro-beams via micro-structure-dependent beam theory. Applied Mathematics and Mechanics, 39(7), pp.923-952.

[19] Arefi, M. and Takhayor Ardestani, M., 2019. The Effect of Temperature Dependency on the Thermo-Electro-Elastic Analysis of Functionally Graded Piezoelectric Spherical Shell. Mechanics of Advanced Composite Structures‎, 6(1), pp.1-8.

[20] Mirzaei, M.M.H., Loghman, A. and Arefi, M., 2019. Time-dependent creep analysis of a functionally graded beam with trapezoidal cross section using first-order shear deformation theory. Steel and Composite Structures, 30(6), pp.567-576.

[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.657-669.

[22] Arefi, M., Loghman, A. and Mohammad Hosseini Mirzaei, M., 2019. Thermoelastic analysis of a functionally graded simple blade using first-order shear deformation theory. Mechanics of Advanced Composite Structures‎.

[23] 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(5), pp.743-750.

[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.52-60.

[25] Mirzabeigy, A., 2014. Semi-analytical approach for free vibration analysis of variable cross-section beams resting on elastic foundation and under axial force. International Journal of Engineering, Transactions C: Aspects, 27(3), pp.455-463.

[26] 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 Mechanics, 84(5), pp.615-623.

[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.173-179.

[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.1-12.

[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. Size-dependent Effects on the Vibration Behavior of a Ti-moshenko Microbeam subjected to Pre-stress Loading based on DQM. Mechanics of Advanced Composite Structures‎, 3(2), pp.99-112.

[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.300-309.

[32] Calim, F.F., 2016. Free and forced vibration analysis of axially functionally graded Timoshenko beams on two-parameter viscoelastic foundation. Composites Part B: Engineering, 103, pp.98-112.

[33] Soltani, M. and Asgarian, B., 2019. 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. 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 non-prismatic thin-walled beams. Thin-Walled Structures, 82, pp.245-261.

[35] Soltani, M. and Mohri, F., 2016. Stability and vibration analyses of tapered columns resting on one or two-parameter elastic foundations. Journal of Numerical Methods in Civil Engineering, 1(2), pp.57-66.

[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 non-prismatic I-beams. 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 magneto-thermo-electric environment. Journal of Sandwich Structures & Materials, 18(5), pp.624-651.

[39] Arefi, M. and Zenkour, A.M., 2017. Vibration and bending analysis of a sandwich microbeam with two integrated piezo-magnetic face-sheets. Composite Structures, 159, pp.479-490.

[40] Arefi, M. and Zenkour, A.M., 2017. Size-dependent vibration and bending analyses of the piezomagnetic three-layer 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.