# Stability Analysis of Nonuniform Rectangular Beams Using Homotopy Perturbation Method

Authors: Seval Pinarbasi
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

The design of slender beams, that is, beams with large laterally unsupported lengths, is commonly controlled by stability limit states. Beam buckling, also called “lateral torsional buckling,” is different from column buckling in that a beam not only displaces laterally but also twists about its axis during buckling. The coupling between twist and lateral displacement makes stability analysis of beams more complex than that of columns. For this reason, most of the analytical studies in the literature on beam stability are concentrated on simple cases: uniform beams with ideal boundary conditions and simple loadings. This paper shows that complex beam stability problems, such as lateral torsional buckling of rectangular beams with variable cross-sections, can successfully be solved using homotopy perturbation method (HPM).


---

## Body

## 1. Introduction

A beam is a structural element which spans large distances between supports and which primarily carries transverse loads with negligible axial loads. If a beam has sufficient lateral bracing, it can easily be designed by selecting the most economical “compact” cross-section satisfying the strength and serviceability limit states. However, just like slender columns which buckle under compressive loads much smaller than their “stable” load carrying capacities, a “laterally unbraced” slender beam can also buckle under transverse loads. For this reason, the design of slender beams has to consider stability limit states as well.Beam buckling, which is also called “lateral torsional buckling,” differs from column buckling in that a beam not only displaces laterally but also rotates about its axis during buckling. The coupling between twist and outward lateral displacement makes stability analysis of beams more complex than that of columns. For this reason, most of the analytical studies in the literature are concentrated on simple cases: uniform beams with ideal boundary conditions and simple loadings. For exact solutions to simple beam buckling problems, one can refer to one of the well-known structural stability books, such as [1–4].However, in an attempt to construct ever-stronger and ever-lighter structures, many engineers currently design light slender members with variable cross-sections. Unfortunately, design engineers are lack of sufficient guidance on design of nonuniform structural elements since most of the provisions in design specifications are developed for uniform elements. Consequently, there is a need for a practical tool to analyze complex beam stability problems.In recent years, many analytical approaches such as homotopy perturbation method (HPM), Adomian decomposition method (ADM), and variational iteration method (VIM), are proposed for the solution of nonlinear equations, and many researchers have shown that complex engineering problems, which do not have exact closed-form solutions, can easily be solved using these techniques. A review of some recently developed nonlinear analytical techniques is given in [5]. A kind of nonlinear analytical technique which was proposed by He [6] in 1999, homotopy perturbation method (HPM) has many successful applications to various kinds of nonlinear problems. For a review of the state-of-the-art of HPM, the work by He [7] can be referred to. Very recently, HPM is also applied to stability problems of columns. Coşkun [8] and Coşkun [9] and Atay [10] analyzed the elastic stability of Euler columns with variable cross-sections under different loading and boundary conditions using HPM and verified that HPM is a very efficient and powerful technique in buckling analysis of columns with variable cross-sections.In this paper, this powerful analytical technique is applied to two fundamental beam stability problems: lateral torsional buckling of (i) simply supported rectangular beams under pure bending and (ii) cantilever rectangular beams subjected to a concentrated load at their free ends. In the analyses, two different types of stiffness variations, linear and exponential variations, are considered. Exact solutions to these problems, some of which are considerably complex, are available in literature only for uniform beams and some particular cases of linearly tapered beams. For this reason, before studying beams with variable cross-sections, uniform beams with constant cross-sections are analyzed and HPM solutions are compared with the exact solutions. After verifying the effectiveness of HPM in solving lateral buckling problems, HPM is applied to more complex beam buckling problems.


## 2. Lateral Torsional Buckling of Rectangular Beams

### 2.1.  Basic Theory

Consider a narrow rectangular beam subjected to an arbitrary loading iny-z plane causing its bending about its strong axis (x). Locate x, y, z coordinate system to define the undeformed configuration of the beam as shown in Figure 1(a). Similarly, locate ξ, η, ζ coordinate system at the centroid of the cross section at an arbitrary section of the beam along its length to define the deformed configuration of the beam as shown in Figure 1(b).Undeformed and deformed shapes of a narrow rectangular beam loaded to bend about its major axis.
(a)
undeformed shape(b)
deformed (bended) and buckled shapeThe deformation of the beam can be defined by lateral (u) and vertical (v) displacements of the centroid of the beam and angle of twist (ϕ) of the cross section (Figure 1(b)). Assume u and v are positive in the positive directions of x and y, respectively. Then, obeying the right-hand rule, ϕ is positive about positive z axis. Hence, while the twist shown in Figure 1(b) is positive, the displacements are both negative.For small deformations, the cosines of the angles between axes are as listed in Table1. Also, the curvatures in xz and yz planes can be taken as d2u/dz2 and d2v/dz2, respectively. Since one can realistically take the “warping rigidity” of a narrow rectangular beam as zero, the equilibrium equations for the buckled (deformed) beam can be written [1] as follows:
(2.1)EIξd2vdz2=Mξ,EIηd2udz2=Mη,GItdϕdz=Mζ,
representing, respectively, the major-axis bending, minor-axis bending, and twisting of the beam. In (2.1), EIξ and EIη denote, respectively, the strong-axis and weak-axis flexural stiffnesses of the beam and GIt denotes the torsional stiffness of the beam. Positive directions of internal moments are defined in Figure 2.Table 1
Cosines of angles between axes [1].
xyzξ1ϕ-du/dzη−ϕ1-dv/dzζdu/dzdv/dz1Figure 2
Positive directions for internal moments.


### 2.2.  Lateral Buckling of Simply Supported Beams under Pure Bending

Consider a simply supported rectangular beam with variable flexural and torsional stiffnessesEIξ(z), EIη(z), and GIt(z) along its length L (Figure 3). Under pure bending, the beam is subjected to equal end moments Mo about x-axis. The bending and twisting moments at any cross section can be found by determining the components of Mo about ξ, η, ζ axes. Considering the sign convention defined in Figure 2 and using Table 1, these components can be written as
(2.2)Mξ=Mo,Mη=ϕMo,Mζ=(-dudz)Mo.
Substituting (2.2) into (2.1) yields
(2.3)[EIξ(z)]d2vdz2=Mo,[EIη(z)]d2udz2=ϕMo,[GIt(z)]dϕdz=(-dudz)Mo.
It is apparent from (2.3) that v is independent from u and ϕ. Thus, in this problem, it is sufficient to consider only the coupled equations between u and ϕ. Differentiating the last equation in (2.3) with respect to z and using the resulting equation to eliminate u in the second equation in (2.3) give the following second-order differential equation for the angle of twist (ϕ) of the beam:
(2.4)d2ϕdz2+d[GIt(z)]dz1[GIt(z)]dϕdz+Mo2[GIt(z)][EIη(z)]ϕ=0.
The boundary conditions for (2.4) can be written from the end conditions of the beams. Since the ends of the beam are restrained against rotation about z axis, ϕ=0 at both z=0 and z=L.Figure 3
Simply supported rectangular beam under pure bending.


#### 2.2.1.  Beams with Constant Stiffnesses

If the minor-axis flexural and torsional stiffnesses of the beam are constant, that is,EIη(z)=EIη and GIt(z)=GIt, then (2.4) reduces to the following simpler equation:
(2.5)d2ϕdz2+Mo2GItEIηϕ=0.
For easier computations, the nondimensional form of (2.5) can be written as follows:
(2.6)(ϕ̅)′′+α(ϕ̅)=0
with
(2.7)α=Mo2L2GItEIη,
where z̅=z/L, ϕ̅=ϕ, prime denotes differentiation with respect to z̅ and α is the “nondimensional critical moment.” The boundary conditions for this buckling problem can also be written in nondimensional form as
(2.8)ϕ̅(0)=0ϕ̅(1)=0.
It is to be noted that (2.8) is also applicable to beams with variable stiffnesses.


#### 2.2.2.  Beams with Linearly Varying Stiffnesses

If both the minor axis flexural and torsional stiffnesses of the beam changes in linear form, that is, if(2.9)GIt(z)=GIt(1+bzL),EIη(z)=EIη(1+bzL),
where b is a constant determining the “sharpness” of stiffness changes along the length of the beam, then the buckling equation (2.4) becomes
(2.10)d2ϕdz2+bL+bzdϕdz+Mo2GItEIηL2(1+bz/L)2ϕ=0,
the nondimensional form of which can be written as
(2.11)(ϕ̅)′′+b1+bz̅(ϕ̅)′+α1(1+bz̅)2(ϕ̅)=0.


#### 2.2.3.  Beams with Exponentially Varying Stiffnesses

If the beam stiffnesses change in the following exponential form:(2.12)GIt(z)=GIte-a(z/L),EIη(z)=EIηe-a(z/L).
where a is a positive constant, then the nondimensional form of the buckling equation can be written as
(2.13)(ϕ̅)′′-a(ϕ̅)′+αe2az̅(ϕ̅)=0,


### 2.3.  Lateral Buckling of Cantilever Beams with Vertical End Load

Consider a narrow rectangular cantilever beam of lengthL (Figure 4(a)) with flexural and torsional stiffnesses EIξ(z), EIη(z), and GIt(z). When subjected to a vertical load P passing through its centroid at its free end (Figure 4(a)), the beam deforms as shown in Figure 4(b). u1 is the lateral displacement of the loaded end of the beam. The components of the moments of the load at an arbitrary section m-n about x, y, z axes are
(2.14)Mx=-P(L-z),My=0,Mz=P(-u1+u).
Using the sign convention defined in Figure 2, the bending and twisting moments at this arbitrary section can be written as
(2.15)Mξ=-P(L-z),Mη=-ϕP(L-z),Mζ=(dudz)P(L-z)-P(u1-u).
Then, the equilibrium equations for the buckled beam become
(2.16)[EIξ(z)]d2vdz2=-P(L-z),[EIη(z)]d2udz2=-ϕP(L-z),[GIt(z)]dϕdz=(dudz)P(L-z)-P(u1-u).
Similar to the pure bending case, v is independent from u and ϕ. Differentiating the last equation in (2.16) with respect to z and using the resulting equation to eliminate u in the second equation in (2.16), the following second order differential equation is obtained for ϕ:
(2.17)d2ϕdz2+d[GIt(z)]dz1[GIt(z)]dϕdz+P2[GIt(z)][EIη(z)](L-z)2ϕ=0.
Since the fixed end of the beam is restrained against rotation and since the twisting moment at the free end is known to be zero, the boundary conditions for this problem are ϕ=0 at z=0 and dϕ/dz=0 at z=L.Lateral buckling of a narrow rectangular cantilever beam carrying concentrated load at its free end.
(a)
undeformed shape(b)
buckled shape


#### 2.3.1.  Beams with Constant Stiffnesses

IfEIη(z)=EIη and GIt(z)=GIt, then (2.17) takes the following simpler form:
(2.18)d2ϕdz2+P2GItEIη(L-z)2ϕ=0.
Equation. (2.18) can be rewritten in nondimensional form as
(2.19)(ϕ̅)′′+β(1-z̅)2(ϕ̅)=0,
where the “nondimensional critical load” β is defined as
(2.20)β=P2L4GItEIη.
The boundary conditions for this buckling problem can be written in nondimensional form as(2.21)ϕ̅(0)=0,dϕ̅dz̅(1)=0,
which are also applicable to the beams with variable stiffnesses.


#### 2.3.2.  Beams with Linearly Varying Stiffnesses

If both stiffnesses of the beam change in linear form, that is, if(2.22)GIt(z)=GIt(1-bzL),EIη(z)=EIη(1-bzL),
where b is a positive constant that can take values between zero and one, then (2.17) becomes(2.23)d2ϕdz2-bL-bzdϕdz+P2L2GItEIη(1-z/L)2(1-bz/L)2ϕ=0,
which, when written in nondimensional form, takes the following simpler form:
(2.24)(ϕ̅)′′-b1-bz̅(ϕ̅)′+β(1-z̅)2(1-bz̅)2(ϕ̅)=0


#### 2.3.3.  Beams with Exponentially Varying Stiffnesses

If beam stiffnesses change as in (2.12), the nondimensional form of (2.17) can be written as(2.25)(ϕ̅)′′-a(ϕ̅)′+βe2az̅(1-z̅2)(ϕ̅)=0.


## 2.1.  Basic Theory

Consider a narrow rectangular beam subjected to an arbitrary loading iny-z plane causing its bending about its strong axis (x). Locate x, y, z coordinate system to define the undeformed configuration of the beam as shown in Figure 1(a). Similarly, locate ξ, η, ζ coordinate system at the centroid of the cross section at an arbitrary section of the beam along its length to define the deformed configuration of the beam as shown in Figure 1(b).Undeformed and deformed shapes of a narrow rectangular beam loaded to bend about its major axis.
(a)
undeformed shape(b)
deformed (bended) and buckled shapeThe deformation of the beam can be defined by lateral (u) and vertical (v) displacements of the centroid of the beam and angle of twist (ϕ) of the cross section (Figure 1(b)). Assume u and v are positive in the positive directions of x and y, respectively. Then, obeying the right-hand rule, ϕ is positive about positive z axis. Hence, while the twist shown in Figure 1(b) is positive, the displacements are both negative.For small deformations, the cosines of the angles between axes are as listed in Table1. Also, the curvatures in xz and yz planes can be taken as d2u/dz2 and d2v/dz2, respectively. Since one can realistically take the “warping rigidity” of a narrow rectangular beam as zero, the equilibrium equations for the buckled (deformed) beam can be written [1] as follows:
(2.1)EIξd2vdz2=Mξ,EIηd2udz2=Mη,GItdϕdz=Mζ,
representing, respectively, the major-axis bending, minor-axis bending, and twisting of the beam. In (2.1), EIξ and EIη denote, respectively, the strong-axis and weak-axis flexural stiffnesses of the beam and GIt denotes the torsional stiffness of the beam. Positive directions of internal moments are defined in Figure 2.Table 1
Cosines of angles between axes [1].
xyzξ1ϕ-du/dzη−ϕ1-dv/dzζdu/dzdv/dz1Figure 2
Positive directions for internal moments.


## 2.2.  Lateral Buckling of Simply Supported Beams under Pure Bending

Consider a simply supported rectangular beam with variable flexural and torsional stiffnessesEIξ(z), EIη(z), and GIt(z) along its length L (Figure 3). Under pure bending, the beam is subjected to equal end moments Mo about x-axis. The bending and twisting moments at any cross section can be found by determining the components of Mo about ξ, η, ζ axes. Considering the sign convention defined in Figure 2 and using Table 1, these components can be written as
(2.2)Mξ=Mo,Mη=ϕMo,Mζ=(-dudz)Mo.
Substituting (2.2) into (2.1) yields
(2.3)[EIξ(z)]d2vdz2=Mo,[EIη(z)]d2udz2=ϕMo,[GIt(z)]dϕdz=(-dudz)Mo.
It is apparent from (2.3) that v is independent from u and ϕ. Thus, in this problem, it is sufficient to consider only the coupled equations between u and ϕ. Differentiating the last equation in (2.3) with respect to z and using the resulting equation to eliminate u in the second equation in (2.3) give the following second-order differential equation for the angle of twist (ϕ) of the beam:
(2.4)d2ϕdz2+d[GIt(z)]dz1[GIt(z)]dϕdz+Mo2[GIt(z)][EIη(z)]ϕ=0.
The boundary conditions for (2.4) can be written from the end conditions of the beams. Since the ends of the beam are restrained against rotation about z axis, ϕ=0 at both z=0 and z=L.Figure 3
Simply supported rectangular beam under pure bending.


### 2.2.1.  Beams with Constant Stiffnesses

If the minor-axis flexural and torsional stiffnesses of the beam are constant, that is,EIη(z)=EIη and GIt(z)=GIt, then (2.4) reduces to the following simpler equation:
(2.5)d2ϕdz2+Mo2GItEIηϕ=0.
For easier computations, the nondimensional form of (2.5) can be written as follows:
(2.6)(ϕ̅)′′+α(ϕ̅)=0
with
(2.7)α=Mo2L2GItEIη,
where z̅=z/L, ϕ̅=ϕ, prime denotes differentiation with respect to z̅ and α is the “nondimensional critical moment.” The boundary conditions for this buckling problem can also be written in nondimensional form as
(2.8)ϕ̅(0)=0ϕ̅(1)=0.
It is to be noted that (2.8) is also applicable to beams with variable stiffnesses.


### 2.2.2.  Beams with Linearly Varying Stiffnesses

If both the minor axis flexural and torsional stiffnesses of the beam changes in linear form, that is, if(2.9)GIt(z)=GIt(1+bzL),EIη(z)=EIη(1+bzL),
where b is a constant determining the “sharpness” of stiffness changes along the length of the beam, then the buckling equation (2.4) becomes
(2.10)d2ϕdz2+bL+bzdϕdz+Mo2GItEIηL2(1+bz/L)2ϕ=0,
the nondimensional form of which can be written as
(2.11)(ϕ̅)′′+b1+bz̅(ϕ̅)′+α1(1+bz̅)2(ϕ̅)=0.


### 2.2.3.  Beams with Exponentially Varying Stiffnesses

If the beam stiffnesses change in the following exponential form:(2.12)GIt(z)=GIte-a(z/L),EIη(z)=EIηe-a(z/L).
where a is a positive constant, then the nondimensional form of the buckling equation can be written as
(2.13)(ϕ̅)′′-a(ϕ̅)′+αe2az̅(ϕ̅)=0,


## 2.2.1.  Beams with Constant Stiffnesses

If the minor-axis flexural and torsional stiffnesses of the beam are constant, that is,EIη(z)=EIη and GIt(z)=GIt, then (2.4) reduces to the following simpler equation:
(2.5)d2ϕdz2+Mo2GItEIηϕ=0.
For easier computations, the nondimensional form of (2.5) can be written as follows:
(2.6)(ϕ̅)′′+α(ϕ̅)=0
with
(2.7)α=Mo2L2GItEIη,
where z̅=z/L, ϕ̅=ϕ, prime denotes differentiation with respect to z̅ and α is the “nondimensional critical moment.” The boundary conditions for this buckling problem can also be written in nondimensional form as
(2.8)ϕ̅(0)=0ϕ̅(1)=0.
It is to be noted that (2.8) is also applicable to beams with variable stiffnesses.


## 2.2.2.  Beams with Linearly Varying Stiffnesses

If both the minor axis flexural and torsional stiffnesses of the beam changes in linear form, that is, if(2.9)GIt(z)=GIt(1+bzL),EIη(z)=EIη(1+bzL),
where b is a constant determining the “sharpness” of stiffness changes along the length of the beam, then the buckling equation (2.4) becomes
(2.10)d2ϕdz2+bL+bzdϕdz+Mo2GItEIηL2(1+bz/L)2ϕ=0,
the nondimensional form of which can be written as
(2.11)(ϕ̅)′′+b1+bz̅(ϕ̅)′+α1(1+bz̅)2(ϕ̅)=0.


## 2.2.3.  Beams with Exponentially Varying Stiffnesses

If the beam stiffnesses change in the following exponential form:(2.12)GIt(z)=GIte-a(z/L),EIη(z)=EIηe-a(z/L).
where a is a positive constant, then the nondimensional form of the buckling equation can be written as
(2.13)(ϕ̅)′′-a(ϕ̅)′+αe2az̅(ϕ̅)=0,


## 2.3.  Lateral Buckling of Cantilever Beams with Vertical End Load

Consider a narrow rectangular cantilever beam of lengthL (Figure 4(a)) with flexural and torsional stiffnesses EIξ(z), EIη(z), and GIt(z). When subjected to a vertical load P passing through its centroid at its free end (Figure 4(a)), the beam deforms as shown in Figure 4(b). u1 is the lateral displacement of the loaded end of the beam. The components of the moments of the load at an arbitrary section m-n about x, y, z axes are
(2.14)Mx=-P(L-z),My=0,Mz=P(-u1+u).
Using the sign convention defined in Figure 2, the bending and twisting moments at this arbitrary section can be written as
(2.15)Mξ=-P(L-z),Mη=-ϕP(L-z),Mζ=(dudz)P(L-z)-P(u1-u).
Then, the equilibrium equations for the buckled beam become
(2.16)[EIξ(z)]d2vdz2=-P(L-z),[EIη(z)]d2udz2=-ϕP(L-z),[GIt(z)]dϕdz=(dudz)P(L-z)-P(u1-u).
Similar to the pure bending case, v is independent from u and ϕ. Differentiating the last equation in (2.16) with respect to z and using the resulting equation to eliminate u in the second equation in (2.16), the following second order differential equation is obtained for ϕ:
(2.17)d2ϕdz2+d[GIt(z)]dz1[GIt(z)]dϕdz+P2[GIt(z)][EIη(z)](L-z)2ϕ=0.
Since the fixed end of the beam is restrained against rotation and since the twisting moment at the free end is known to be zero, the boundary conditions for this problem are ϕ=0 at z=0 and dϕ/dz=0 at z=L.Lateral buckling of a narrow rectangular cantilever beam carrying concentrated load at its free end.
(a)
undeformed shape(b)
buckled shape


### 2.3.1.  Beams with Constant Stiffnesses

IfEIη(z)=EIη and GIt(z)=GIt, then (2.17) takes the following simpler form:
(2.18)d2ϕdz2+P2GItEIη(L-z)2ϕ=0.
Equation. (2.18) can be rewritten in nondimensional form as
(2.19)(ϕ̅)′′+β(1-z̅)2(ϕ̅)=0,
where the “nondimensional critical load” β is defined as
(2.20)β=P2L4GItEIη.
The boundary conditions for this buckling problem can be written in nondimensional form as(2.21)ϕ̅(0)=0,dϕ̅dz̅(1)=0,
which are also applicable to the beams with variable stiffnesses.


### 2.3.2.  Beams with Linearly Varying Stiffnesses

If both stiffnesses of the beam change in linear form, that is, if(2.22)GIt(z)=GIt(1-bzL),EIη(z)=EIη(1-bzL),
where b is a positive constant that can take values between zero and one, then (2.17) becomes(2.23)d2ϕdz2-bL-bzdϕdz+P2L2GItEIη(1-z/L)2(1-bz/L)2ϕ=0,
which, when written in nondimensional form, takes the following simpler form:
(2.24)(ϕ̅)′′-b1-bz̅(ϕ̅)′+β(1-z̅)2(1-bz̅)2(ϕ̅)=0


### 2.3.3.  Beams with Exponentially Varying Stiffnesses

If beam stiffnesses change as in (2.12), the nondimensional form of (2.17) can be written as(2.25)(ϕ̅)′′-a(ϕ̅)′+βe2az̅(1-z̅2)(ϕ̅)=0.


## 2.3.1.  Beams with Constant Stiffnesses

IfEIη(z)=EIη and GIt(z)=GIt, then (2.17) takes the following simpler form:
(2.18)d2ϕdz2+P2GItEIη(L-z)2ϕ=0.
Equation. (2.18) can be rewritten in nondimensional form as
(2.19)(ϕ̅)′′+β(1-z̅)2(ϕ̅)=0,
where the “nondimensional critical load” β is defined as
(2.20)β=P2L4GItEIη.
The boundary conditions for this buckling problem can be written in nondimensional form as(2.21)ϕ̅(0)=0,dϕ̅dz̅(1)=0,
which are also applicable to the beams with variable stiffnesses.


## 2.3.2.  Beams with Linearly Varying Stiffnesses

If both stiffnesses of the beam change in linear form, that is, if(2.22)GIt(z)=GIt(1-bzL),EIη(z)=EIη(1-bzL),
where b is a positive constant that can take values between zero and one, then (2.17) becomes(2.23)d2ϕdz2-bL-bzdϕdz+P2L2GItEIη(1-z/L)2(1-bz/L)2ϕ=0,
which, when written in nondimensional form, takes the following simpler form:
(2.24)(ϕ̅)′′-b1-bz̅(ϕ̅)′+β(1-z̅)2(1-bz̅)2(ϕ̅)=0


## 2.3.3.  Beams with Exponentially Varying Stiffnesses

If beam stiffnesses change as in (2.12), the nondimensional form of (2.17) can be written as(2.25)(ϕ̅)′′-a(ϕ̅)′+βe2az̅(1-z̅2)(ϕ̅)=0.


## 3. Formulations of the Studied Buckling Problems Using HPM

### 3.1. Brief Review of HPM

Consider a general nonlinear differential equation,(3.1)A(u)-f(r)=0,r∈Ω,
with boundary conditions
(3.2)B(u,∂u∂n)=0,r∈Γ,
where A is a general differential operator, B is a boundary operator, f(r) is a known analytic function, and Γ is the boundary of the domain Ω [6]. Dividing the operator A into linear (L) and nonlinear (N) parts, the differential equation can be written as follows:
(3.3)L(u)+N(u)-f(r)=0.
The basic idea of homotopy perturbation technique (HPM) is to construct a homotopy v(r,p)=Ω×[0,1]→ℝ which satisfies
(3.4)H(v,p)=(1-p)[L(v)-L(u0)]+p[L(v)+N(v)-f(r)]=0,
where p∈[0,1] is an embedding parameter, u0 is an initial approximation satisfying the boundary conditions. Equation (3.4) can be rearranged in the following form:
(3.5)H(v,p)=L(v)-L(u0)+pL(u0)+p[N(v)-f(r)]=0.
From (3.5), it is obvious that
(3.6)H(v,0)=L(v)-L(u0)=0,H(v,1)=N(v)+L(v)-f(r)=0.
In other words, as p changes from zero to unity, v(r,p) changes from u0 to u(r). Using the embedding parameter as a small parameter, HPM defines the solution of (3.5) as
(3.7)v=v0+pv1+p2v2+p3v3+⋯.
Thus, the approximate solution of (3.1) or (3.3) can be obtained from(3.8)u=lim⁡p→1v=v0+v1+v2+v3+⋯.


### 3.2. HPM Formulations of the Studied Buckling Equations

The nondimensional forms of the buckling equations derived for the studied stability problems are presented in (2.6), (2.11), (2.13), (2.19), (2.24), and (2.25). One can see that all of these equations can be written in the following form:(3.9)(ϕ̅)′′+λ1(ϕ̅)′+λ2(ϕ̅)=0,
where λ1 and λ2 are coefficient functions which depend on stiffness variations, end conditions and loading of the beam. For example, for the buckling problem of a cantilever beam with linearly varying stiffnesses along its length and carrying concentrated load at its free end, these functions are(3.10)λ1=-b1-bz̅,λ2=β(1-z̅)2(1-bz̅)2.
As it can also be inferred from (3.10) that, for particular values of a or b, λ1 is a function of z̅ only, while λ2 is function of both z̅ and the nondimensional critical moment α or load β.The linear and nonlinear parts of (3.9) can be taken as
(3.11)L(ϕ̅)=(ϕ̅)′′,N(ϕ̅)=λ1(z̅)(ϕ̅)′+λ2(z̅)(ϕ̅)=0,
with
(3.12)f(r)=0.
Substituting (3.7) into (3.5), in view of (3.9), (3.11) and (3.12), and equating the terms with similar powers of the embedding parameter p, the following iteration equations are obtained:
(3.13)p0:(v0)′′=(ϕ̅0)′′,p1:(v1)′′=-(ϕ̅0)′′-λ1(v0)′-λ2(v0),p2:(v2)′′=-λ1(v1)′-λ2(v1),p3:(v3)′′=-λ1(v2)′-λ2(v2),⋮pn:(vn)′′=-λ1(vn-1)′-λ2(vn-1).
For all cases considered in the study, the solution of the linear part of (3.9), that is, L(ϕ̅)=(ϕ̅)′′=0, can be taken as an initial guess ϕ̅0. Thus,
(3.14)ϕ̅0(z̅)=Az̅+B,
where A and B are unknown coefficients to be determined from the boundary conditions of the problems. Substituting (3.14), into the equations given in (3.13), vi (i: 0–n) can be obtained with n successive iterations. Finally, the approximate solution can be obtained from
(3.15)ϕ̅=lim⁡p→1v≅∑i=0nvi.
For each particular case of the studied problems, substituting the approximate solution to the related boundary conditions, two homogeneous equations are obtained in terms of the unknown coefficients A and B. These equations can be put into the following matrix form:
(3.16)[M(αorβ)]{AB}={00}.
Thus, each problem reduces to an eigenvalue problem. For a nontrivial solution, the determinant of the coefficient matrix has to be zero, that is, |M(αorβ)|=0. The smallest possible real root of the characteristic equation gives the nondimensional buckling moment or load (αorβ) in the first buckling mode.


## 3.1. Brief Review of HPM

Consider a general nonlinear differential equation,(3.1)A(u)-f(r)=0,r∈Ω,
with boundary conditions
(3.2)B(u,∂u∂n)=0,r∈Γ,
where A is a general differential operator, B is a boundary operator, f(r) is a known analytic function, and Γ is the boundary of the domain Ω [6]. Dividing the operator A into linear (L) and nonlinear (N) parts, the differential equation can be written as follows:
(3.3)L(u)+N(u)-f(r)=0.
The basic idea of homotopy perturbation technique (HPM) is to construct a homotopy v(r,p)=Ω×[0,1]→ℝ which satisfies
(3.4)H(v,p)=(1-p)[L(v)-L(u0)]+p[L(v)+N(v)-f(r)]=0,
where p∈[0,1] is an embedding parameter, u0 is an initial approximation satisfying the boundary conditions. Equation (3.4) can be rearranged in the following form:
(3.5)H(v,p)=L(v)-L(u0)+pL(u0)+p[N(v)-f(r)]=0.
From (3.5), it is obvious that
(3.6)H(v,0)=L(v)-L(u0)=0,H(v,1)=N(v)+L(v)-f(r)=0.
In other words, as p changes from zero to unity, v(r,p) changes from u0 to u(r). Using the embedding parameter as a small parameter, HPM defines the solution of (3.5) as
(3.7)v=v0+pv1+p2v2+p3v3+⋯.
Thus, the approximate solution of (3.1) or (3.3) can be obtained from(3.8)u=lim⁡p→1v=v0+v1+v2+v3+⋯.


## 3.2. HPM Formulations of the Studied Buckling Equations

The nondimensional forms of the buckling equations derived for the studied stability problems are presented in (2.6), (2.11), (2.13), (2.19), (2.24), and (2.25). One can see that all of these equations can be written in the following form:(3.9)(ϕ̅)′′+λ1(ϕ̅)′+λ2(ϕ̅)=0,
where λ1 and λ2 are coefficient functions which depend on stiffness variations, end conditions and loading of the beam. For example, for the buckling problem of a cantilever beam with linearly varying stiffnesses along its length and carrying concentrated load at its free end, these functions are(3.10)λ1=-b1-bz̅,λ2=β(1-z̅)2(1-bz̅)2.
As it can also be inferred from (3.10) that, for particular values of a or b, λ1 is a function of z̅ only, while λ2 is function of both z̅ and the nondimensional critical moment α or load β.The linear and nonlinear parts of (3.9) can be taken as
(3.11)L(ϕ̅)=(ϕ̅)′′,N(ϕ̅)=λ1(z̅)(ϕ̅)′+λ2(z̅)(ϕ̅)=0,
with
(3.12)f(r)=0.
Substituting (3.7) into (3.5), in view of (3.9), (3.11) and (3.12), and equating the terms with similar powers of the embedding parameter p, the following iteration equations are obtained:
(3.13)p0:(v0)′′=(ϕ̅0)′′,p1:(v1)′′=-(ϕ̅0)′′-λ1(v0)′-λ2(v0),p2:(v2)′′=-λ1(v1)′-λ2(v1),p3:(v3)′′=-λ1(v2)′-λ2(v2),⋮pn:(vn)′′=-λ1(vn-1)′-λ2(vn-1).
For all cases considered in the study, the solution of the linear part of (3.9), that is, L(ϕ̅)=(ϕ̅)′′=0, can be taken as an initial guess ϕ̅0. Thus,
(3.14)ϕ̅0(z̅)=Az̅+B,
where A and B are unknown coefficients to be determined from the boundary conditions of the problems. Substituting (3.14), into the equations given in (3.13), vi (i: 0–n) can be obtained with n successive iterations. Finally, the approximate solution can be obtained from
(3.15)ϕ̅=lim⁡p→1v≅∑i=0nvi.
For each particular case of the studied problems, substituting the approximate solution to the related boundary conditions, two homogeneous equations are obtained in terms of the unknown coefficients A and B. These equations can be put into the following matrix form:
(3.16)[M(αorβ)]{AB}={00}.
Thus, each problem reduces to an eigenvalue problem. For a nontrivial solution, the determinant of the coefficient matrix has to be zero, that is, |M(αorβ)|=0. The smallest possible real root of the characteristic equation gives the nondimensional buckling moment or load (αorβ) in the first buckling mode.


## 4. HPM Solutions to the Studied Stability Problems

### 4.1. Critical Moments for Pure Bending Cases

Exact solution to (2.5) is given in [1] as
(4.1)Mcr=πGItEIηL,
where Mcr is the critical moment of the beam in the first buckling mode. In view of (2.7), this result corresponds to a nondimensional critical moment of
(4.2)α=π2≅9.8696.
In order to show how HPM is applied to the studied buckling problem and how the approximate solutions converge to the exact solution as the number of iterations increases, (2.6) is solved using different number of iterations defined in (3.13) with the initial guess given in (3.14). As an example, the terms obtained for the first five iterations are given below:
(4.3)v0=Az̅+B,v1=-16A(z̅)3α-16B(z̅)2α,v2=1120A(z̅)5α2+124B(z̅)4α2,v3=-15040A(z̅)7α3-1720B(z̅)6α3,v4=1362880A(z̅)9α4+140320B(z̅)8α4,v5=-139916800A(z̅)11α5-13628800B(z̅)10α5.
In fact, even five iterations are sufficient to obtain almost exact result when beam stiffnesses are constant along beam length as shown in Figure 5, where the convergence of HPM solutions to the exact one with increasing number of iterations is shown. Error is only 0.03% when n=5.Figure 5
Pure bending case: Beams with constant stiffnesses: convergence of HPM solution to the exact solution as the number of iterations increase.Exact solution to (2.10) is also available in the literature ([3]):
(4.4)Mcr=πbln⁡(1+b)GItEIηL,
which in view of (2.7) equals
(4.5)α=π2b2[ln⁡(1+b)]2
The normalized buckling moments for two particular values of b, 0.3 and 0.7, are computed using HPM for different numbers of iterations, and the convergences of the approximate results to the exact ones are shown in Figure 6. To simplify the integration processes, variable coefficients in the iteration integrals, that is, λ1 and λ2, are expanded in series using nine terms. As it is seen from Figure 6(a), for b=0.3, HPM solutions converge to the exact result as the number of iterations increases, and to obtain the exact result, it is sufficient to perform only eight iterations. On the other hand, when b=0.7, there remains some small error, not more than 1%, even when nine iterations are performed. This is due to the fact that as b increases, that is, as the nonlinearity in λ1 and λ2 increases, it becomes necessary to expand these coefficients in series using more terms in iteration integrals. As given in Table 2, as the number of terms in series is increased, HPM results converge to the exact result (α=17.1757).Table 2
Pure bending base: Beams with linearly varying stiffnesses (b=0.7). HPM solutions for nondimensional buckling moment (α) for different numbers of expansion of variable coefficients in the iteration integrals in series.
n9 terms11 terms15 terms17 terms21 terms25 terms29 terms17.558737.503227.470037.465697.462927.462367.46225212.843212.688012.590712.577312.568512.566612.5663314.663314.571914.517814.510814.506414.505514.5053418.297718.289818.296818.299118.301218.301718.3019517.151017.091317.060117.056617.054617.054317.0542617.277517.222717.195217.192317.190617.190417.1903717.262917.207417.179317.176317.174717.174417.1743817.264217.208817.180817.177917.176217.175917.1759917.264117.208717.180717.177717.176117.175817.1757Pure bending case: beams with linearly varying stiffnesses: convergence of HPM solutions to the exact solutions as the number of iterations increases.
(a)
b=0.3(b)
b=0.7To investigate the effects of exponential stiffness variations on buckling moment of a simply supported rectangular beam under pure bending, (2.13) is solved using HPM for various values of a and the smallest α values in the first buckling modes are obtained. For all values of a, the variable coefficient in the iteration integrals is expanded in series using seventeen terms and nine iterations are conducted. To the best knowledge of the author, there is no exact solution in the literature for this case of the problem. Critical moments of nonuniform beams normalized to that of the uniform beam (a=0) are plotted in Figure 7, which shows severe decrease in buckling moment as a increases.Figure 7
Pure bending case: beams with exponentially varying stiffnesses: variation of nondimensional critical moment ratio (α/αa=0) with “a” values.


### 4.2. Critical Loads for Cantilever Cases

Exact solution to (2.18) is given in Timoshenko and Gere [1] as
(4.6)J-1/4(PL2/2GItEIη)=0,
where J-1/4 represent the Bessel function of the first kind of order −1/4. The smallest root of this equation yields the first mode critical load Pcr(4.7)Pcr=4.0126GItEIηL2.
In view of (2.20), this result corresponds to a nondimensional critical load of
(4.8)β=(4.0126)2=16.1010.
Even though it may be rather difficult to obtain the roots of a Bessel function, one can solve (4.6) to obtain the critical loads in higher modes. The results extremely depend on the initial guess and while deciding which root corresponds to which mode, one should be very careful. After some trial and errors, the larger two roots of (4.6) are obtained, which correspond to nondimensional buckling loads of 104.9830 and 272.775 in the second and third modes, respectively.Equation (2.19) is solved using HPM to evaluate the effectiveness of HPM in determining buckling loads in higher modes. Seventeen iterations are performed to get the higher mode values. Unlike the exact solutions, the roots of HPM results are much easier to determine since the characteristic equation obtained using HPM is a polynomial. This is one of the advantages of using HPM in this problem, even in the case of constant stiffnesses. Figure 8 shows how HPM solutions converge to the exact solutions as the numbers of iterations increases. While it is sufficient to execute six iterations to achieve the exact result for the first mode, iteration numbers have to be increased for higher mode values; ten iterations for the second mode and fifteen iterations for the third mode.Figure 8
Cantilever case: beams with constant stiffnesses: convergences of HPM solutions to the exact solutions as the number of iterations increases for nondimensional buckling loads in the first three buckling modes.To the best knowledge of author, there are no exact solutions in literature for lateral buckling of nonuniform cantilever beams supporting a concentrated load at its free end when the beam stiffnesses vary along its length linearly or exponentially due to the complex buckling equations ((2.24) and (2.25)) to be solved.Approximate solutions to (2.24) and (2.25) are obtained using HPM with nine iterations and presented in Figures 9 and 10. Variable coefficients in the iteration integrals are expanded in series using seventeen terms. Since the buckling load values of nonuniform beams plotted in Figures 9 and 10 are normalized with respect to that of uniform beams (a=0 or b=0), one can easily see how the buckling capacities of nonuniform beams decrease as the beam stiffnesses decrease. When b=1, the capacity drops 35% of its uniform capacity, whereas when a=2, the capacity drop is almost quarter.Figure 9
Cantilever case: beams with linearly varying stiffnesses: variation of nondimensional critical moment ratio (β/βb=0) with “b” values.Figure 10
Cantilever case: beams with exponentially varying stiffnesses: variation of nondimensional critical moment ratio (β/βa=0) with “a” values.


## 4.1. Critical Moments for Pure Bending Cases

Exact solution to (2.5) is given in [1] as
(4.1)Mcr=πGItEIηL,
where Mcr is the critical moment of the beam in the first buckling mode. In view of (2.7), this result corresponds to a nondimensional critical moment of
(4.2)α=π2≅9.8696.
In order to show how HPM is applied to the studied buckling problem and how the approximate solutions converge to the exact solution as the number of iterations increases, (2.6) is solved using different number of iterations defined in (3.13) with the initial guess given in (3.14). As an example, the terms obtained for the first five iterations are given below:
(4.3)v0=Az̅+B,v1=-16A(z̅)3α-16B(z̅)2α,v2=1120A(z̅)5α2+124B(z̅)4α2,v3=-15040A(z̅)7α3-1720B(z̅)6α3,v4=1362880A(z̅)9α4+140320B(z̅)8α4,v5=-139916800A(z̅)11α5-13628800B(z̅)10α5.
In fact, even five iterations are sufficient to obtain almost exact result when beam stiffnesses are constant along beam length as shown in Figure 5, where the convergence of HPM solutions to the exact one with increasing number of iterations is shown. Error is only 0.03% when n=5.Figure 5
Pure bending case: Beams with constant stiffnesses: convergence of HPM solution to the exact solution as the number of iterations increase.Exact solution to (2.10) is also available in the literature ([3]):
(4.4)Mcr=πbln⁡(1+b)GItEIηL,
which in view of (2.7) equals
(4.5)α=π2b2[ln⁡(1+b)]2
The normalized buckling moments for two particular values of b, 0.3 and 0.7, are computed using HPM for different numbers of iterations, and the convergences of the approximate results to the exact ones are shown in Figure 6. To simplify the integration processes, variable coefficients in the iteration integrals, that is, λ1 and λ2, are expanded in series using nine terms. As it is seen from Figure 6(a), for b=0.3, HPM solutions converge to the exact result as the number of iterations increases, and to obtain the exact result, it is sufficient to perform only eight iterations. On the other hand, when b=0.7, there remains some small error, not more than 1%, even when nine iterations are performed. This is due to the fact that as b increases, that is, as the nonlinearity in λ1 and λ2 increases, it becomes necessary to expand these coefficients in series using more terms in iteration integrals. As given in Table 2, as the number of terms in series is increased, HPM results converge to the exact result (α=17.1757).Table 2
Pure bending base: Beams with linearly varying stiffnesses (b=0.7). HPM solutions for nondimensional buckling moment (α) for different numbers of expansion of variable coefficients in the iteration integrals in series.
n9 terms11 terms15 terms17 terms21 terms25 terms29 terms17.558737.503227.470037.465697.462927.462367.46225212.843212.688012.590712.577312.568512.566612.5663314.663314.571914.517814.510814.506414.505514.5053418.297718.289818.296818.299118.301218.301718.3019517.151017.091317.060117.056617.054617.054317.0542617.277517.222717.195217.192317.190617.190417.1903717.262917.207417.179317.176317.174717.174417.1743817.264217.208817.180817.177917.176217.175917.1759917.264117.208717.180717.177717.176117.175817.1757Pure bending case: beams with linearly varying stiffnesses: convergence of HPM solutions to the exact solutions as the number of iterations increases.
(a)
b=0.3(b)
b=0.7To investigate the effects of exponential stiffness variations on buckling moment of a simply supported rectangular beam under pure bending, (2.13) is solved using HPM for various values of a and the smallest α values in the first buckling modes are obtained. For all values of a, the variable coefficient in the iteration integrals is expanded in series using seventeen terms and nine iterations are conducted. To the best knowledge of the author, there is no exact solution in the literature for this case of the problem. Critical moments of nonuniform beams normalized to that of the uniform beam (a=0) are plotted in Figure 7, which shows severe decrease in buckling moment as a increases.Figure 7
Pure bending case: beams with exponentially varying stiffnesses: variation of nondimensional critical moment ratio (α/αa=0) with “a” values.


## 4.2. Critical Loads for Cantilever Cases

Exact solution to (2.18) is given in Timoshenko and Gere [1] as
(4.6)J-1/4(PL2/2GItEIη)=0,
where J-1/4 represent the Bessel function of the first kind of order −1/4. The smallest root of this equation yields the first mode critical load Pcr(4.7)Pcr=4.0126GItEIηL2.
In view of (2.20), this result corresponds to a nondimensional critical load of
(4.8)β=(4.0126)2=16.1010.
Even though it may be rather difficult to obtain the roots of a Bessel function, one can solve (4.6) to obtain the critical loads in higher modes. The results extremely depend on the initial guess and while deciding which root corresponds to which mode, one should be very careful. After some trial and errors, the larger two roots of (4.6) are obtained, which correspond to nondimensional buckling loads of 104.9830 and 272.775 in the second and third modes, respectively.Equation (2.19) is solved using HPM to evaluate the effectiveness of HPM in determining buckling loads in higher modes. Seventeen iterations are performed to get the higher mode values. Unlike the exact solutions, the roots of HPM results are much easier to determine since the characteristic equation obtained using HPM is a polynomial. This is one of the advantages of using HPM in this problem, even in the case of constant stiffnesses. Figure 8 shows how HPM solutions converge to the exact solutions as the numbers of iterations increases. While it is sufficient to execute six iterations to achieve the exact result for the first mode, iteration numbers have to be increased for higher mode values; ten iterations for the second mode and fifteen iterations for the third mode.Figure 8
Cantilever case: beams with constant stiffnesses: convergences of HPM solutions to the exact solutions as the number of iterations increases for nondimensional buckling loads in the first three buckling modes.To the best knowledge of author, there are no exact solutions in literature for lateral buckling of nonuniform cantilever beams supporting a concentrated load at its free end when the beam stiffnesses vary along its length linearly or exponentially due to the complex buckling equations ((2.24) and (2.25)) to be solved.Approximate solutions to (2.24) and (2.25) are obtained using HPM with nine iterations and presented in Figures 9 and 10. Variable coefficients in the iteration integrals are expanded in series using seventeen terms. Since the buckling load values of nonuniform beams plotted in Figures 9 and 10 are normalized with respect to that of uniform beams (a=0 or b=0), one can easily see how the buckling capacities of nonuniform beams decrease as the beam stiffnesses decrease. When b=1, the capacity drops 35% of its uniform capacity, whereas when a=2, the capacity drop is almost quarter.Figure 9
Cantilever case: beams with linearly varying stiffnesses: variation of nondimensional critical moment ratio (β/βb=0) with “b” values.Figure 10
Cantilever case: beams with exponentially varying stiffnesses: variation of nondimensional critical moment ratio (β/βa=0) with “a” values.


## 5. Conclusions

The design of slender beams with large laterally unsupported lengths is usually governed by their lateral torsional buckling capacities. In this limit state, structural deformation of the beam suddenly changes from in-plane deformation (strong-axis bending) to a combination of out-of plane deformation (weak-axis bending) and twisting. If the slenderness of the beam is considerably large, the lateral buckling capacity of the beam can be much smaller than its strong-axis bending capacity.In an attempt to construct ever-stronger and ever-lighter structures, many engineers currently design light slender members with variable cross-sections, which are especially prone to this type of buckling. Unfortunately, design engineers are lack of sufficient guidance on design of such nonuniform structural elements since most of the provisions in design specifications are developed for uniform elements. For nonuniform members, buckling equations usually become so complex that it becomes impractical and sometimes even impossible to obtain exact closed-form solutions to these equations. However, approximate solutions can easily be obtained to these complex problems using recently developed nonlinear analytical techniques, such as homotopy perturbation method (HPM).In this paper, two fundamental beam buckling problems, lateral torsional buckling of (i) simply supported rectangular beams subjected to pure bending and (ii) rectangular cantilever beams carrying concentrated load at their free ends, is studied using HPM. Exact solutions to these problems are available in literature only for uniform beams and some particular cases of linearly tapered beams. In order to verify the effectiveness of HPM on solving beam stability problems and to show the application of the method, first the lateral buckling of uniform beams are studied. The excellent match of the HPM results with the exact results verifies the efficiency of the technique in the analysis of lateral torsional buckling problems. Then, beams with variable minor-axis flexural and torsional stiffnesses along their lengths are studied. Both linear and exponential variations are considered in nonuniform beams. The stability analyses of nonuniform beams lead to differential equations with variable coefficients, for which it can be rather difficult to derive exact solutions. However, as shown in the paper, it is relatively easy to write HPM algorithms to these complex differential equations, which give buckling moment/load of the beam after a few iterations.



---

*Source: 197483-2012-02-27.xml*