**2. Problem formulation**

Derivation of governing equations related to stability analysis is given in detail in Timoshenko and Gere [5], Simitses and Hodges [6], and Wang et al. [8]. The reader can also refer to any textbook related to the subject. In this section, only the governing equation will be given for the related cases.

Consider the elastic columns given in Fig.1. The governing equation for the buckling of such columns is

$$\frac{d^2y}{d\mathfrak{x}^2}\bigg[EI(\mathfrak{x})\frac{d^2y}{d\mathfrak{x}^2}\bigg] + P\frac{d^2y}{d\mathfrak{x}^2} = 0\tag{1}$$

Elastic Stability Analysis of Euler Columns Using Analytical Approximate Techniques 117

**Figure 1.** Elastic columns with various end conditions

In the case of constant flexural rigidity (*i.e. EI* is constant), Eq.(1) becomes

$$\frac{d^4y}{dx^4} + \frac{P}{EI}\frac{d^2y}{dx^2} = 0\tag{2}$$

where *EI* is the flexural rigidity of the column, and *P* is the applied load. Both Eqs. (2) and (3) are solved due to end conditions of the column. Some of these conditions are shown in Fig.1. In this figure, letters are used for a simplification to describe the support conditions of the column. The first letter stands for the support at the bottom and the second letter for the top. Hence, CF is *Clamped-Fixed*, PP is *Pinned-Pinned*, C-P is *Clamped-Pinned* and C-S is *Clamped-Sliding Restraint*.

The governing equations (1) and (2) are both solved with respect to the problem's end conditions. The end conditions for the columns shown in Fig.1 are given below:

Pin support:

116 Advances in Computational Stability Analysis

foundations. [28-29]

discretized solution techniques.

related to the topic of the article.

**2. Problem formulation** 

be given for the related cases.

columns is

distribution method were analyzed by Kerekes [14]. The research of Siginer [15] was about the stability of a column whose flexural stiffness has a continuous linear variation along the column. Moreover, the analytical solutions of a multi-step bar with varying cross section were obtained by Li et al. [16-18]. The energy method was used by Sampaio et al. [19] to find the solution for the problem of buckling behavior of inclined beamcolumn. Some of the important researchers who studied the mechanical behavior of beamcolumns are Keller [20], Tadjbakhsh and Keller [21] and Taylor [22]. Later on, analytical approximate techniques were used for the stability analysis of elastic columns. Coşkun and Atay [23] and Atay and Coşkun [24] studied column buckling problems for the columns with variable flexural stiffness and for the columns with continuous elastic restraints by using the variational iteration method which produces analytical approximations. Coşkun [25, 26] used the homotopy perturbation method for buckling of Euler columns on elastic foundations and tilt-buckling of variable stiffness columns. Pnarbaş [27] also analyzed the stability of nonuniform rectangular beams using homotopy perturbation method. These techniques were also used successfully in the vibration analysis of Euler-Bernoulli beams and in the vibration of beams on elastic

Recently, by the emergence of new and innovative semi analytical approximation methods, research on this subject has gained momentum. Analytical approximate solution techniques are used widely to solve nonlinear ordinary or partial differential equations, integrodifferential equations, delay equations, etc. The main advantage of employing such techniques is that the problems are considered in a more realistic manner, and the solution obtained is a continuous function which is not the case for the solutions obtained by

The methods that will be used throughout this study are, Adomian Decomposition Method (ADM), Variational Iteration Method (VIM) and Homotopy Perturbation Method (HPM). Each technique will be explained first, and then all will be applied to a selected case study

Derivation of governing equations related to stability analysis is given in detail in Timoshenko and Gere [5], Simitses and Hodges [6], and Wang et al. [8]. The reader can also refer to any textbook related to the subject. In this section, only the governing equation will

Consider the elastic columns given in Fig.1. The governing equation for the buckling of such

(1)

2 22 2 22 ( ) <sup>0</sup> *dy dy dy EI x P dx dx dx* 

$$y = 0 \quad \text{and} \quad \frac{d^2y}{d\mathbf{x}^2} = 0\tag{3}$$

Clamped support:

$$
\delta y = 0 \quad \text{and} \quad \frac{dy}{d\mathbf{x}} = 0 \tag{4}
$$

Free end:

$$\text{and } \frac{d^3y}{d\mathbf{x}^3} + \frac{P}{EI}\frac{dy}{d\mathbf{x}} = 0 \tag{5}$$

Sliding restraint:

$$\frac{dy}{d\chi} = 0 \text{ and} \tag{6}$$

Elastic Stability Analysis of Euler Columns Using Analytical Approximate Techniques 119

*Lu Ru Nu g x* ( ) (10)

*Lu g x Ru Nu* ( ) (11)

1 11 *u L g x L Ru L Nu* ( ) (12)

1 1 *u f x L Ru L Nu* ( ) (13)

(14)

(9)

2 2 22 22 44

*S T ST T S S T S T*

2 ( ) sin sin 0

Although Eq.(8) is a linear equation with constant coefficients, obtaining a solution from Eq.(9) is not that easy. It is very interesting that, even with a software, one can not easily produce the buckling loads in a sequential order from Eq.(9). In view of this experience, an analytical solution for Eq.(7) is almost impossible to obtain except very limited *EI(x)*

Hence, analytical approximate techniques are efficient alternatives for solving these problems. By the use of these techniques, a solution which is continuous in the problem domain is possible for any variation in flexural rigidity. These techniques produce the buckling loads in a sequential order, and it is also very easy to obtain the buckling mode shapes from the solution provided by the method used. These are great advantages in the

**3. The methods used in the elastic stability analysis of Euler columns** 

where, *L* is the linear operator which is highest order derivative, *R* is the remainder of linear operator including derivatives of less order than *L*, *Nu* represents the nonlinear terms, and *g*

Applying the inverse operator *L-1* to both sides of Eq.(11) and employing given conditions;

After integrating source term and combining it with the terms arising from given conditions

The nonlinear operator *Nu F u* ( ) is represented by an infinite series of specially generated (Adomian) polynomials for the specific nonlinearity. Assuming *Nu* is analytical, we write

> 0 ( ) *<sup>k</sup> k Fu A*

 

2 2 ( ) 2 cos cos ( ) ( )

choices.

we obtain

solution of such problems.

**3.1. Adomian Decomposition Method (ADM)** 

is the source term. Eq.(10) can be rearranged as

of the problem, a function *f(x)* is defined in the equation as

In the ADM a differential equation of the following form is considered

*ST S T T S*

The governing equation given in Eq.(2) is a fourth order differential equation with constant coefficients which makes it possible to obtain analytical solutions easily. However, Eq.(1) includes variable coefficients due to variable flexural rigidity. For this type of differential equations, analytical solutions are limited for the special cases of *EI(x)* only. It is not possible to obtain a solution for any form of the function *EI(x)*.

In some problems, obtaining analytical solutions is very difficult even for a constant coefficient governing equation. Consider the buckling of a column on an elastic foundation shown in Fig.2.

**Figure 2.** Column with continuous elastic restraints.

The governing equation for the column in Fig.2 is

$$
\left[\frac{d^2y}{d\mathbf{x}^2}\right]EI(\mathbf{x})\frac{d^2y}{d\mathbf{x}^2}\Big] + P\frac{d^2y}{d\mathbf{x}^2} + ky = 0\tag{7}
$$

which, for the constant *EI* becomes

$$\frac{d^4y}{dx^4} + \frac{P}{EI}\frac{d^2y}{dx^2} + \frac{k}{EI}y = 0\tag{8}$$

In Eqs.(7) and (8), *k* is the stiffness parameter for the elastic restraint. The solution of Eq.(8) for the *CF* column is given in [8] as

Elastic Stability Analysis of Euler Columns Using Analytical Approximate Techniques 119

$$\begin{aligned} \left[\alpha(\mathbb{S}^2 + T^2) - 2\mathbb{S}^2 T^2\right] \cos T \cos S - \alpha(\mathbb{S}^2 + T^2) + (\mathbb{S}^4 + T^4) + \\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad \begin{aligned} \text{ST}\left[2a - (\mathbb{S}^2 + T^2)\right] \sin T \sin S &= 0 \end{aligned} \tag{9}$$

Although Eq.(8) is a linear equation with constant coefficients, obtaining a solution from Eq.(9) is not that easy. It is very interesting that, even with a software, one can not easily produce the buckling loads in a sequential order from Eq.(9). In view of this experience, an analytical solution for Eq.(7) is almost impossible to obtain except very limited *EI(x)* choices.

Hence, analytical approximate techniques are efficient alternatives for solving these problems. By the use of these techniques, a solution which is continuous in the problem domain is possible for any variation in flexural rigidity. These techniques produce the buckling loads in a sequential order, and it is also very easy to obtain the buckling mode shapes from the solution provided by the method used. These are great advantages in the solution of such problems.

### **3. The methods used in the elastic stability analysis of Euler columns**

#### **3.1. Adomian Decomposition Method (ADM)**

118 Advances in Computational Stability Analysis

to obtain a solution for any form of the function *EI(x)*.

**Figure 2.** Column with continuous elastic restraints.

which, for the constant *EI* becomes

for the *CF* column is given in [8] as

The governing equation for the column in Fig.2 is

2 22

4 2

2 22 ( ) <sup>0</sup> *dy dy dy EI x P ky dx dx dx*

4 2 <sup>0</sup> *dy dy P k <sup>y</sup> dx dx EI EI*

In Eqs.(7) and (8), *k* is the stiffness parameter for the elastic restraint. The solution of Eq.(8)

<sup>0</sup> *dy*

The governing equation given in Eq.(2) is a fourth order differential equation with constant coefficients which makes it possible to obtain analytical solutions easily. However, Eq.(1) includes variable coefficients due to variable flexural rigidity. For this type of differential equations, analytical solutions are limited for the special cases of *EI(x)* only. It is not possible

In some problems, obtaining analytical solutions is very difficult even for a constant coefficient governing equation. Consider the buckling of a column on an elastic foundation

*dx* and (6)

(7)

(8)

Sliding restraint:

shown in Fig.2.

In the ADM a differential equation of the following form is considered

$$L\mu + R\mu + N\mu = \mathfrak{g}(\mathfrak{x})\tag{10}$$

where, *L* is the linear operator which is highest order derivative, *R* is the remainder of linear operator including derivatives of less order than *L*, *Nu* represents the nonlinear terms, and *g* is the source term. Eq.(10) can be rearranged as

$$L\boldsymbol{\mu} = \boldsymbol{g}(\mathbf{x}) - \mathbf{R}\boldsymbol{\mu} - \mathbf{N}\boldsymbol{\mu} \tag{11}$$

Applying the inverse operator *L-1* to both sides of Eq.(11) and employing given conditions; we obtain

$$u = L^{-1}\left(g(\mathbf{x})\right) - L^{-1}\left(Ru\right) - L^{-1}\left(Nu\right) \tag{12}$$

After integrating source term and combining it with the terms arising from given conditions of the problem, a function *f(x)* is defined in the equation as

$$\mu = f(\mathbf{x}) - L^{-1} \left( R\mu \right) - L^{-1} \left( N\mu \right) \tag{13}$$

The nonlinear operator *Nu F u* ( ) is represented by an infinite series of specially generated (Adomian) polynomials for the specific nonlinearity. Assuming *Nu* is analytical, we write

$$F(\mu) = \sum\_{k=0}^{\infty} A\_k \tag{14}$$

The polynomials *Ak*'s are generated for all kinds of nonlinearity, so that they depend only on *uo* to *uk* components and can be produced by the following algorithm.

$$A\_0 = F(\mu\_0) \tag{15}$$

Elastic Stability Analysis of Euler Columns Using Analytical Approximate Techniques 121

*Lu Nu f r r* ( ) ( ) ( ) , (24)

*vrp* ( , ) [0,1] *R* (25)

<sup>0</sup> *Hv Lv Lu* ( ,0) ( ) ( ) 0 (28)

*Hv Lv Nv f r* ( ,1) ( ) ( ) ( ) 0 (29)

01 2 3 *v v pv p v p v* ... (30)

(31)

<sup>0</sup> *Hvp p Lv Lu pLv Nv f r* ( , ) (1 )[ ( ) ( )] [ ( ) ( ) ( )] 0 (26)

0 0 *H v p L v L u pL u p N v f r* ( , ) ( ) ( ) ( ) [ ( ) ( )] 0 (27)

HPM provides an analytical approximate solution for problems at hand as the other previously explained techniques. Brief theoretical steps for the equation of following type

with boundary conditions *Bu u n* (, ) 0 . In Eq.(24) *L* is a linear operator, *N* is a nonlinear operator, *B* is a boundary operator, and *f(r)* is a known analytic function. HPM defines

where *r* and *p*[0,1] is an imbedding parameter, *u0* is an initial approximation which

As *p* is changing from zero to unity, so is that of *vrp* (, ) from 0 *u* to *u r*( ) . In topology, this deformation 0 *Lv Lu* () ( ) and *Lv Nv f r* () () () are called homotopic. The basic assumption is that the solutions of Eq.(34) and Eq.(35) can be expressed as a power series in *p* such that:

> <sup>0123</sup> <sup>1</sup> lim ... *<sup>p</sup> u vv v v v*

The convergence of the series in Eq.(31) has been proved in [33]. The method is described in

The governing equation for this case was previously given in Eq.(1). ADM, VIM and HPM will be applied to this equation in order to compute the buckling loads for the clamped-

2 3

satisfies the boundary conditions. Obviously, from Eq.(26) and Eq.(27) , we have :

The approximate solution of *Lu Nu f r r* ( ) ( ) ( ) , can be obtained as :

**3.3. Homotopy Perturbation Method (HPM)** 

which satisfies the following inequalities:

detail in references [33-36].

**4.1. Buckling of a clamped-pinned column** 

**4. Case studies** 

can be given as

homotopy as

or

$$A\_1 = \mu\_1 F'(\mu\_0) \tag{16}$$

$$A\_2 = \mu\_2 F'(\mu\_0) + \frac{1}{2!} \mu\_1^2 F''(\mu\_0) \tag{17}$$

$$A\_3 = \mu\_3 F(\mu\_0) + \mu\_1 \mu\_2 F''(\mu\_0) + \frac{1}{3!} \mu\_1^3 F'''(\mu\_0) \tag{18}$$

The reader can refer to [30, 31] for the algorithms used in formulating Adomian polynomials. The solution *u(x)* is defined by the following series

$$\mu = \sum\_{k=0}^{\infty} u\_k \tag{19}$$

where, the components of the series are determined recursively as follows:

$$
\mu\_0 = f(\mathbf{x}) \tag{20}
$$

$$
\mu\_{k+1} = -L^{-1} \left( R \mu\_k \right) - L^{-1} \left( A\_k \right), \quad k \ge 0 \tag{21}
$$

#### **3.2. Variational Iteration Method (VIM)**

According to VIM, the following differential equation may be considered:

$$
\Delta \mathbf{u} + \mathbf{N} \mathbf{u} = \mathbf{g}(\mathbf{x}) \tag{22}
$$

where *L* is a linear operator, *N* is a nonlinear operator, and *g(x)* is an inhomogeneous source term. Based on VIM, a correct functional can be constructed as follows:

$$\mu\_{n+1} = \mu\_n + \int\_0^\chi \mathcal{J}(\xi) \left\{ L\mu\_n(\xi) + N\tilde{\mu}\_n(\xi) - g(\xi) \right\} \, d\xi \tag{23}$$

where is a general Lagrangian multiplier which can be identified optimally via the variational theory. The subscript *n* denotes the *nth*-order approximation, *u* is considered as a restricted variation *i.e.* 0 *u* . By solving the differential equation for obtained from Eq.(23) in view of 0 *u* with respect to its boundary conditions, Lagrangian multiplier can be obtained. For further details of the method the reader can refer to [32].

#### **3.3. Homotopy Perturbation Method (HPM)**

HPM provides an analytical approximate solution for problems at hand as the other previously explained techniques. Brief theoretical steps for the equation of following type can be given as

$$L(\mu) + N(\mu) = f(r) \quad , \ r \in \Omega \tag{24}$$

with boundary conditions *Bu u n* (, ) 0 . In Eq.(24) *L* is a linear operator, *N* is a nonlinear operator, *B* is a boundary operator, and *f(r)* is a known analytic function. HPM defines homotopy as

$$w(r, p) = \Omega \times [0, 1] \to \mathbb{R} \tag{25}$$

which satisfies the following inequalities:

$$H(\upsilon, p) = (1 - p)[L(\upsilon) - L(\mu\_0)] + p[L(\upsilon) + N(\upsilon) - f(r)] = 0 \tag{26}$$

or

120 Advances in Computational Stability Analysis

The polynomials *Ak*'s are generated for all kinds of nonlinearity, so that they depend only on

220 1 0 <sup>1</sup> () () 2!

3 3 0 12 0 1 0 <sup>1</sup> () () () 3!

 The reader can refer to [30, 31] for the algorithms used in formulating Adomian

> 0 *k k u u*

 1 1 <sup>1</sup> , 0 *k kk u L Ru L A k*

where *L* is a linear operator, *N* is a nonlinear operator, and *g(x)* is an inhomogeneous source

<sup>1</sup>

 

variational theory. The subscript *n* denotes the *nth*-order approximation, *u* is considered as a

*nn n n u u Lu Nu g d*

( ) ( ) ( ) ( )

is a general Lagrangian multiplier which can be identified optimally via the

*u* . By solving the differential equation for

 

*u* with respect to its boundary conditions, Lagrangian multiplier

 (23)

2

0 0 *A F u*( ) (15)

11 0 *A uF u*( ) (16)

(19)

<sup>0</sup> *u fx* ( ) (20)

*Lu Nu g x* ( ) (22)

obtained from

(21)

*A uF u uF u* (17)

3

*A uF u uu F u u F u* (18)

*uo* to *uk* components and can be produced by the following algorithm.

polynomials. The solution *u(x)* is defined by the following series

**3.2. Variational Iteration Method (VIM)** 

where

restricted variation *i.e.* 0

Eq.(23) in view of 0

where, the components of the series are determined recursively as follows:

According to VIM, the following differential equation may be considered:

term. Based on VIM, a correct functional can be constructed as follows:

0

*x*

 

can be obtained. For further details of the method the reader can refer to [32].

$$H(\upsilon, p) = L(\upsilon) - L(\mu\_0) + pL(\mu\_0) + p[N(\upsilon) - f(r)] = 0 \tag{27}$$

where *r* and *p*[0,1] is an imbedding parameter, *u0* is an initial approximation which satisfies the boundary conditions. Obviously, from Eq.(26) and Eq.(27) , we have :

$$LH(\upsilon,0) = L(\upsilon) - L(\mu\_0) = 0\tag{28}$$

$$H(\upsilon, 1) = L(\upsilon) + N(\upsilon) - f(r) = 0 \tag{29}$$

As *p* is changing from zero to unity, so is that of *vrp* (, ) from 0 *u* to *u r*( ) . In topology, this deformation 0 *Lv Lu* () ( ) and *Lv Nv f r* () () () are called homotopic. The basic assumption is that the solutions of Eq.(34) and Eq.(35) can be expressed as a power series in *p* such that:

$$
\Delta v = v\_0 + pv\_1 + p^2v\_2 + p^3v\_3 + \dots \tag{30}
$$

The approximate solution of *Lu Nu f r r* ( ) ( ) ( ) , can be obtained as :

$$\mu = \lim\_{p \to 1} v = v\_0 + v\_1 + v\_2 + v\_3 + \dots \tag{31}$$

The convergence of the series in Eq.(31) has been proved in [33]. The method is described in detail in references [33-36].

#### **4. Case studies**

#### **4.1. Buckling of a clamped-pinned column**

The governing equation for this case was previously given in Eq.(1). ADM, VIM and HPM will be applied to this equation in order to compute the buckling loads for the clamped-

pinned column with constant flexural stiffness, *i.e.* constant *EI*, and variable flexural stiffness, *i.e.* variable *EI* with its corresponding mode shapes. To achieve this aim, a circular rod is defined with an exponentially varying radius. The case is given in Fig.3, and rod and its associated boundary conditions are also provided in Eqs.(3-4). As a case study, first the formulations for constant stiffness column by using ADM, VIM and HPM are given, and then applied to the governing equation of the problem. Afterwards, a variable flexural rigidity will be defined for the same column, and the same techniques will be used for the analysis.

**Figure 3.** *CP* column with constant and variable flexural rigidity

#### **4.2. Formulation of the algorithms for uniform column**

#### *4.2.1. ADM*

The linear operator and its inverse operator for Eq.(2) is

$$L(\cdot) = \frac{d^4}{d\alpha^4}(\cdot) \tag{32}$$

Elastic Stability Analysis of Euler Columns Using Analytical Approximate Techniques 123

1 ( ) 2! 3!

(35)

(36)

<sup>0123</sup> *YY Y Y Y* ... (37)

(38)

*iv LY Y* (40)

01 2 3 *Y Y pY p Y p Y* ... (42)

0 0 <sup>0</sup> *iv iv Y y* (43)

(41)

would be obtained for

(34)

*LY Y* 

*x x Y A Bx C D L Y*

1 <sup>1</sup> ( ), 0 *Y LY n n n* 

Based on the formulation given previously, Lagrange multiplier,

Finally, the solution is defined by

the governing equation, *i.e.* Eq.(2), as

equation into the formulation given in Eq.(31) as

*4.2.2. VIM* 

where 

*4.2.3. HPM* 

2 3

<sup>3</sup>

( ) 3! *x*

An iterative algorithm can be constructed inserting Lagrange multiplier and governing

<sup>1</sup>

the algorithm is chosen as the solution of 0 *LY* which is a cubic polynomial with four

*NY Y* 

2 3

*iv YY Y Y d nn n n* 

( ) ( ) ( )

 (39)

 

is normalized buckling load for the column considered. Initial approximation for

 

0

unknowns to be determined by the end conditions of the column.

Based on the formulation, Eq.(2) can be divided into two parts as

The solution can be expressed as a power series in *p* such that

Inserting Eq.(50) into Eq.(35) provides a solution algorithm as

*x*

$$L^{-1}(\cdot) = \bigcap\_{0 \le 0} \bigcap\_{0 \le 0} \{ \cdot \} \text{ dɔ x } d\mathbf{x} \text{ dɔ} \tag{33}$$

To keep the formulation a general one for all configurations to be considered, the boundary conditions are chosen as *Y A* (0) , *Y B* (0) , *Y C* (0) and *Y D* (0) . Suitable values should be replaced in the formulation with these constants. In this case, 0 *A* and 0 *B* should be inserted for the *CP* column. Hence, the equation to be solved and the recursive algorithm can be given as

Elastic Stability Analysis of Euler Columns Using Analytical Approximate Techniques 123

$$LY = \zeta Y'' \tag{34}$$

$$Y = A + B\mathbf{x} + C\frac{\mathbf{x}^2}{2!} + D\frac{\mathbf{x}^3}{3!} + L^{-1}(\zeta Y'') \tag{35}$$

$$Y\_{n+1} = L^{-1} \langle \zeta \Upsilon\_n'' \rangle, \quad n \ge 0 \tag{36}$$

Finally, the solution is defined by

$$Y = Y\_0 + Y\_1 + Y\_2 + Y\_3 + \dots \tag{37}$$

#### *4.2.2. VIM*

122 Advances in Computational Stability Analysis

*P*

pinned column with constant flexural stiffness, *i.e.* constant *EI*, and variable flexural stiffness, *i.e.* variable *EI* with its corresponding mode shapes. To achieve this aim, a circular rod is defined with an exponentially varying radius. The case is given in Fig.3, and rod and its associated boundary conditions are also provided in Eqs.(3-4). As a case study, first the formulations for constant stiffness column by using ADM, VIM and HPM are given, and then applied to the governing equation of the problem. Afterwards, a variable flexural rigidity will

be defined for the same column, and the same techniques will be used for the analysis.

**Figure 3.** *CP* column with constant and variable flexural rigidity

The linear operator and its inverse operator for Eq.(2) is

*4.2.1. ADM* 

algorithm can be given as

**4.2. Formulation of the algorithms for uniform column** 

*x*

1

0000 ( ) ( ) *xxxx*

To keep the formulation a general one for all configurations to be considered, the boundary conditions are chosen as *Y A* (0) , *Y B* (0) , *Y C* (0) and *Y D* (0) . Suitable values should be replaced in the formulation with these constants. In this case, 0 *A* and 0 *B* should be inserted for the *CP* column. Hence, the equation to be solved and the recursive

4 <sup>4</sup> () () *<sup>d</sup> <sup>L</sup>*

*dx* (32)

*x* 

*<sup>L</sup> dx dx dx dx* (33)

*P* 

Based on the formulation given previously, Lagrange multiplier, would be obtained for the governing equation, *i.e.* Eq.(2), as

$$\mathcal{A}(\xi) = \frac{\left(\xi - \alpha\right)^{\mathfrak{Z}}}{\mathfrak{Z}!} \tag{38}$$

An iterative algorithm can be constructed inserting Lagrange multiplier and governing equation into the formulation given in Eq.(31) as

$$Y\_{n+1} = Y\_n + \int\_0^\chi \lambda(\xi) \left\{ Y\_n^{iv}(\xi) + \zeta \tilde{Y}\_n''(\xi) \right\} d\xi \tag{39}$$

where is normalized buckling load for the column considered. Initial approximation for the algorithm is chosen as the solution of 0 *LY* which is a cubic polynomial with four unknowns to be determined by the end conditions of the column.

#### *4.2.3. HPM*

Based on the formulation, Eq.(2) can be divided into two parts as

$$LY = Y^{iv} \tag{40}$$

$$NYY = \zeta Y'' \tag{41}$$

The solution can be expressed as a power series in *p* such that

$$Y = Y\_0 + pY\_1 + p^2Y\_2 + p^3Y\_3 + \dots \tag{42}$$

Inserting Eq.(50) into Eq.(35) provides a solution algorithm as

$$Y\_0^{i\upsilon} - y\_0^{i\upsilon} = \mathbf{0} \tag{43}$$

$$Y\_1^{iv} + y\_0^{iv} + \zeta^4 Y\_0'' = 0 \tag{44}$$

Elastic Stability Analysis of Euler Columns Using Analytical Approximate Techniques 125

Twenty iterations are conducted for each method, and the computed values are compared with the corresponding exact values for the first four modes of buckling in the following table.

From the table it can be seen that the computed values are highly accurate which show that the techniques used in the analysis are very effective. Only a few iterations are enough to obtain the critical buckling load which is Mode 1. Additional modes require additional iterations. The table also shows that additional two or three iterations will produce an excellent agreement for Mode 4. Even with twenty iterations, the error is less than 0.014%

The buckling mode shapes of uniform column for the first four modes are depicted in Fig.4. To prevent a possible confusion to the reader, the exact mode shapes and the computed ones are not shown separately in the figure since the obtained mode shapes coincide with the

*Mode Exact ADM VIM HPM*  20.19072856 20.19072856 20.19072856 20.19072856 59.67951594 59.67951594 59.67951594 59.67951594 118.89986916 118.89986857 118.89986857 118.89986868 197.85781119 197.88525697 197.88522951 197.88520511

**Table 1.** Comparison of normalized buckling loads ( <sup>2</sup> *PL EI* / ) for the *CP* column

for all the methods used in the analyses.

**Figure 4.** Buckling modes of *CP* column.

exact ones.

$$\mathcal{Y}\_n^{iv} + \mathcal{J}\mathcal{Y}\_{n-1}^{\prime} = 0, \quad n \ge 2 \tag{45}$$

Hence, an approximate solution would be obtained as

$$Y = Y\_0 + Y\_1 + Y\_2 + Y\_3 + \dots \tag{46}$$

Initial guess is very important for the convergence of solution in HPM. A cubic polynomial with four unknown coefficients can be chosen as an initial guess which was shown previously to be an effective one in problems related to Euler beams and columns [23-29].

#### **4.3. Computation of buckling loads**

By the use of described techniques, an iterative procedure is constructed and a polynomial including the unknown coefficients resulting from the initial guess is produced as the solution to the governing equation. Besides two unknowns from the initial guess, an additional unknown also exists in the solution. Applying far end boundary conditions to the solution produces a linear algebraic system of equations which can be defined in a matrix form as

$$\begin{bmatrix} M(\zeta) \end{bmatrix} \begin{bmatrix} \alpha \end{bmatrix} = \begin{Bmatrix} 0 \end{Bmatrix} \tag{47}$$

where , *<sup>T</sup> A B* . For a nontrivial solution, determinant of coefficient matrix must be zero. Determinant of matrix *M*( ) yields a characteristic equation in terms of . Positive real roots of this equation are the normalized buckling loads for the Clamped-Pinned column.

#### **4.4. Determination of buckling mode shapes**

Buckling mode shapes for the column can also be obtained from the polynomial approximations by the methods considered in this study. Introducing, the buckling loads into the solution, normalized polynomial eigen functions for the mode shapes are obtained from

$$\overline{Y}\_{j} = \frac{Y\_{N}\left(\mathbf{x}, \boldsymbol{\zeta}\_{j}\right)}{\left[\int\_{0}^{1} \left|Y\_{N}\left(\mathbf{x}, \boldsymbol{\zeta}\_{j}\right)\right|^{2} d\mathbf{x}\right]^{1/2}} \quad j = 1, 2, 3, \dots \tag{48}$$

The same approach can also be employed to predict mode shapes for the cases including variable flexural stiffness.

#### **4.5. Analysis of a uniform column**

After applying the procedures explained in the text, the following results are obtained for the buckling loads. Comparison with the exact solutions is also provided in order that one can observe an excellent agreement between the exact results and the computed results.


Twenty iterations are conducted for each method, and the computed values are compared with the corresponding exact values for the first four modes of buckling in the following table.

**Table 1.** Comparison of normalized buckling loads ( <sup>2</sup> *PL EI* / ) for the *CP* column

124 Advances in Computational Stability Analysis

Hence, an approximate solution would be obtained as

**4.4. Determination of buckling mode shapes** 

**4.3. Computation of buckling loads** 

unknown

where , *<sup>T</sup>* 

Determinant of matrix *M*( )

variable flexural stiffness.

**4.5. Analysis of a uniform column** 

4 10 0 <sup>0</sup> *iv iv Yy Y* 

<sup>1</sup> 0, 2 *iv YY n n n*

Initial guess is very important for the convergence of solution in HPM. A cubic polynomial with four unknown coefficients can be chosen as an initial guess which was shown previously to be an effective one in problems related to Euler beams and columns [23-29].

By the use of described techniques, an iterative procedure is constructed and a polynomial including the unknown coefficients resulting from the initial guess is produced as the solution to the governing equation. Besides two unknowns from the initial guess, an additional

produces a linear algebraic system of equations which can be defined in a matrix form as

roots of this equation are the normalized buckling loads for the Clamped-Pinned column.

 

The same approach can also be employed to predict mode shapes for the cases including

After applying the procedures explained in the text, the following results are obtained for the buckling loads. Comparison with the exact solutions is also provided in order that one can observe an excellent agreement between the exact results and the computed results.

*N j*

 1/2 <sup>1</sup> <sup>2</sup>

,

*Y x Y j Y x dx*

*N j*

0

*j*

Buckling mode shapes for the column can also be obtained from the polynomial approximations by the methods considered in this study. Introducing, the buckling loads into the solution, normalized polynomial eigen functions for the mode shapes are obtained from

 *M*() 0 

also exists in the solution. Applying far end boundary conditions to the solution

*A B* . For a nontrivial solution, determinant of coefficient matrix must be zero.

yields a characteristic equation in terms of

, , 1,2,3,...

(44)

(47)

(48)

. Positive real

(45)

<sup>0123</sup> *YY Y Y Y* ... (46)

From the table it can be seen that the computed values are highly accurate which show that the techniques used in the analysis are very effective. Only a few iterations are enough to obtain the critical buckling load which is Mode 1. Additional modes require additional iterations. The table also shows that additional two or three iterations will produce an excellent agreement for Mode 4. Even with twenty iterations, the error is less than 0.014% for all the methods used in the analyses.

The buckling mode shapes of uniform column for the first four modes are depicted in Fig.4. To prevent a possible confusion to the reader, the exact mode shapes and the computed ones are not shown separately in the figure since the obtained mode shapes coincide with the exact ones.

**Figure 4.** Buckling modes of *CP* column.

#### **4.6. Buckling of a rod with variable cross-section**

A circular rod having a radius changing exponentially is considered in this case. Such a rod is shown below in Fig.5. The function representing the radius would be as

$$R(\mathbf{x}) = R\_0 \mathbf{e}^{-a\mathbf{x}} \tag{49}$$

Elastic Stability Analysis of Euler Columns Using Analytical Approximate Techniques 127

 (56)

*x* (55)

is provided by ADM, buckling

<sup>4</sup> () e *ax*

ADM gives the following formulation with the previously defined fourth order linear

Lagrange multiplier is the same as used in the uniform column case due to the fourth order

*iv Y Y Y aY a x Y d nn n n <sup>n</sup>*

Application of HPM produces the following set of recursive equations as the solution

 <sup>2</sup> 10 0 <sup>0</sup> 8 16 ( ) 0 *iv iv Y y aY a x Y*

 <sup>2</sup> 1 1 8 16 ( ) 0, 2 *Y aY a x Y n n n <sup>n</sup>* 

The proposed formulations are applied for two different variations, *i.e.* 0.1 *aL* and 0.2 *aL* . Twenty iterations are conducted for each method, and the computed normalized buckling

<sup>0</sup> *PL EI* / values are given for the first four modes of buckling in Tables 2 and 3.

*Mode ADM VIM HPM*  16.47361380 16.47361380 16.47361380 48.69674135 48.69674135 48.69674135 97.02096924 97.02096916 97.02096921 161.45155447 161.45151518 161.45150000

 

<sup>0</sup> *PL EI* / ) for *aL* 0.1

derivative in Eq.(38). Hence, an algorithm by using VIM can be constructed as

<sup>0</sup> *PL EI* / . Once

 2 3 1 2 8 16 ( ) 2! 3! *x x Y A B L aY a x Y*

<sup>2</sup>

 

(57)

 

0 0 <sup>0</sup> *iv iv Y y* (58)

(59)

(60)

( ) 8 16 ( )

is normalized buckling load <sup>2</sup>

1

**Table 2.** Comparison of normalized buckling loads ( <sup>2</sup>

0

*x*

 

mode shapes for the rod can also be easily produced from the solution.

where

and, where

operator.

*4.6.1.2. VIM* 

*4.6.1.3. HPM* 

algorithm.

load <sup>2</sup>

*4.6.2. Results of the analyses* 

where *Ro* is the radius at the bottom end, *L* is the length of the rod and 1 *aL* .

**Figure 5.** Circular rod with a variable cross-section

Employing Eq.(49), cross-sectional area and moment of inertia for a section at an arbitrary point *x* becomes:

$$A(\mathbf{x}) = A\_0 \mathbf{e}^{-2a\mathbf{x}} \tag{50}$$

$$I(\mathbf{x}) = I\_0 \mathbf{e}^{-4a\mathbf{x}} \tag{51}$$

where

$$A\_0 = \pi R\_0^2\tag{52}$$

$$I\_0 = \frac{\pi R\_0^4}{4} \tag{53}$$

Governing equation for the rod was previously given in Eq.(1) as

$$\left. \frac{d^2 y}{d\alpha^2} \right| EI(\alpha) \frac{d^2 y}{d\alpha^2} \left| + P \frac{d^2 y}{d\alpha^2} = 0 \right|$$

#### *4.6.1. Formulation of the algorithms*

#### *4.6.1.1. ADM*

Application of ADM leads to the following

$$(Y^{\dot{\upsilon}\upsilon} - 8aY'' + \left(16a^2 + \zeta'\right)\wp(\chi)Y'' = 0\tag{54}$$

where

126 Advances in Computational Stability Analysis

**4.6. Buckling of a rod with variable cross-section** 

**Figure 5.** Circular rod with a variable cross-section

*4.6.1. Formulation of the algorithms* 

Application of ADM leads to the following

point *x* becomes:

*R<sup>o</sup>*

where

*4.6.1.1. ADM* 

A circular rod having a radius changing exponentially is considered in this case. Such a rod

Employing Eq.(49), cross-sectional area and moment of inertia for a section at an arbitrary

*L*

<sup>2</sup> *A R* 0 0 

2 22 2 22 ( ) <sup>0</sup> *dy dy dy EI x P dx dx dx* 

 <sup>2</sup> 8 16 ( ) 0 *iv Y aY a x Y* 

Governing equation for the rod was previously given in Eq.(1) as

2

4

4 0 <sup>0</sup> 4 *<sup>R</sup> <sup>I</sup>* 

<sup>0</sup> () e *ax Rx R* (49)

*x* 

<sup>0</sup> () e *ax Ax A* (50)

<sup>0</sup> () e *ax Ix I* (51)

(52)

(54)

(53)

is shown below in Fig.5. The function representing the radius would be as

where *Ro* is the radius at the bottom end, *L* is the length of the rod and 1 *aL* .

$$\boldsymbol{\wp}(\boldsymbol{\chi}) = \mathbf{e}^{4a\boldsymbol{\chi}} \tag{55}$$

and, where is normalized buckling load <sup>2</sup> <sup>0</sup> *PL EI* / . Once is provided by ADM, buckling mode shapes for the rod can also be easily produced from the solution.

ADM gives the following formulation with the previously defined fourth order linear operator.

$$Y = A \frac{\mathbf{x}^2}{2!} + B \frac{\mathbf{x}^3}{3!} + L^{-1} \left( 8aY''' - \left( 16a^2 + \zeta' \right) \wp(\mathbf{x}) Y'' \right) \tag{56}$$

*4.6.1.2. VIM* 

Lagrange multiplier is the same as used in the uniform column case due to the fourth order derivative in Eq.(38). Hence, an algorithm by using VIM can be constructed as

$$Y\_{n+1} = Y\_n + \int\_0^\varepsilon \mathbb{A}(\xi) \left\{ Y\_n^{iv} - 8a\tilde{Y}\_n{}^\mathbf{w} + \left( 16a^2 + \zeta \right) \wp(\mathbf{x}) \tilde{Y}\_n{}^\mathbf{w} \right\} d\xi \tag{57}$$

#### *4.6.1.3. HPM*

Application of HPM produces the following set of recursive equations as the solution algorithm.

$$Y\_0^{iv} - y\_0^{iv} = 0 \tag{58}$$

$$Y\_1^{iv} + y\_0^{iv} - 8aY\_0''' + \left(16a^2 + \zeta\right)\wp(\mathbf{x})Y\_0'' = 0\tag{59}$$

$$(Y\_n - 8aY\_{n-1} \stackrel{\text{'''}}{+} + \left(16a^2 + \zeta\right)\wp(\mathbf{x})Y\_{n-1} \stackrel{\text{'''}}{=} 0, \quad n \ge 2\tag{60}$$

#### *4.6.2. Results of the analyses*

The proposed formulations are applied for two different variations, *i.e.* 0.1 *aL* and 0.2 *aL* . Twenty iterations are conducted for each method, and the computed normalized buckling load <sup>2</sup> <sup>0</sup> *PL EI* / values are given for the first four modes of buckling in Tables 2 and 3.


**Table 2.** Comparison of normalized buckling loads ( <sup>2</sup> <sup>0</sup> *PL EI* / ) for *aL* 0.1


Elastic Stability Analysis of Euler Columns Using Analytical Approximate Techniques 129

**Figure 7.** Comparison of buckling modes for *CP* rod (Mode 2)

**Figure 8.** Comparison of buckling modes for *CP* rod (Mode 3)

**Table 3.** Comparison of normalized buckling loads ( <sup>2</sup> <sup>0</sup> *PL EI* / ) for *aL* 0.2

The buckling mode shapes of the rod for the first four modes are depicted in between Figs.5- 9. To demonstrate the effect of variable cross-section in the results, a comparison is made with normalized mode shapes for a uniform rod which are given in Fig.4. Constant flexural rigidity is defined as 0.1 *aL* in these figures.

**Figure 6.** Comparison of buckling modes for *CP* rod (Mode 1)

**Figure 7.** Comparison of buckling modes for *CP* rod (Mode 2)

**Table 3.** Comparison of normalized buckling loads ( <sup>2</sup>

rigidity is defined as 0.1 *aL* in these figures.

**Figure 6.** Comparison of buckling modes for *CP* rod (Mode 1)

*Mode ADM VIM HPM*  13.35006457 13.35006457 13.35006457 39.47004813 39.47004813 39.47004813 78.64155457 78.64155458 78.64155466 130.86858532 130.86856343 130.86853842

The buckling mode shapes of the rod for the first four modes are depicted in between Figs.5- 9. To demonstrate the effect of variable cross-section in the results, a comparison is made with normalized mode shapes for a uniform rod which are given in Fig.4. Constant flexural

<sup>0</sup> *PL EI* / ) for *aL* 0.2

**Figure 8.** Comparison of buckling modes for *CP* rod (Mode 3)

Elastic Stability Analysis of Euler Columns Using Analytical Approximate Techniques 131

[1] Vaziri H. H., and Xie. J., Buckling of columns under variably distributed axial loads,

[3] Karman, T.R, Biot, M.A., Mathematical Methods in Engineering, McGraw Hill,

[5] Timoshenko, S.P., Gere, J.M., Theory of Elastic Stability, McGraw Hill, NewYork, 1961. [6] G.J. Simitses , D. H. Hodges, Fundamentals of structural Stability, Elsevier-

[7] Iyengar NGR, Structural Stability of Columns and Plates. New York, John Wiley and

[8] Wang, C.M., Wang, C.Y., Reddy, J.N., Exact Solutions for Buckling of Structural

[9] Ermopoulos, J. C., Buckling of tapered bars under stepped axial loading, J. structural

[10] Li Q. S., exact solutions for buckling of non-uniform columns under axial concentrated

[11] Lee S.Y. and Kuo Y. H., Elastic stability of non-uniform columns, J. Sound Vibration,

[12] Gere, J.M., Carter, W.O., Critical buckling loads for tapered columns, J. Struct. Eng.-

[13] Arbabi. F., Li, F., Buckling of variable cross-section columns – integral equation

[14] Kerekes F., Hulsbos CL., Elastic stability of the top chord of a three-span continuous

[15] Siginer, A., Buckling of columns of variable flexural stiffness, J. Eng. Mech.-ASCE,

[16] Li, Q.S., Cao, H., Li, G., Stability analysis of bars with multi-segments of varying cross

[17] Li, Q.S., Cao, H., Li, G., Stability analysis of bars with varying cross section, Int. J. Solids

[18] Li, Q.S., Cao, H., Li, G.,Static and dynamic analysis of straight bars with variable cross-

[19] Sampaio Jr., J.H.B, Hundhausen, J.R.,A mathematical model and analytical solution for

[20] Keller, J.B.,The shape of the strongest column, Archive for Rational Mechanics and

[21] Tadjbakhsh, I., Keller., J.B., Strongest columns and isoperimetric inequalities for

[22] Taylor, J.E.,The strongest column – an energy approach, J. Appl. Mech. ASME, 34

buckling of inclined beam columns, Appl. Math. Model., 22 (1998), 405-421.

and distributed loading, Eur. J. Mech. A / Solids, 20(2001), 485-500.

approach, J. Struct. Eng.-ASCE, 117(8) (1991), 2426-2441.

pony truss bridge , Iowa Eng. : Expt. Sta. Bull., (1954), 177.

section, Comput. Struct., 53 (1994), 1085-1089.

section, Comput. Struct., 59 (1996), 1185-1191.

eigenvalues, J. Appl. Mech. ASME, 29 (1962), 159-164.

[2] Dinnik, A.N., Design of columns of varying cross-section, Trans. ASME 51, 1929.

[4] Morley A., Critical loads for long tapering struts, Engineering, 104-295, 1917.

**6. References** 

NewYork, 1940.

Sons, 1988

148(1) (1991), 11-24

ASCE, 88(1962), 1-11.

118(1992), 543-640.

Struuct., 32 (1995), 3217-3228.

Analysis, 5 (1960), 275-285.

(1967), 486-487.

Comput. Struct., 45(3) (1992), 505-509.

Butterworth-Heinemann Publishing, 2005.

Members, CRC Press LLC, Florida, 2005.

Eng. ASCE, 112(6) (1986), 1346-1354.

**Figure 9.** Comparison of buckling modes for *CP* rod (Mode 4)
