# An Investigation on the Nonlinear Free Vibration Analysis of Beams with Simply Supported Boundary Conditions Using Four Engineering Theories

Authors: R. A. Jafari-Talookolaei; M. H. Kargarnovin; M. T. Ahmadian; M. Abedi
Journal: Journal of Applied Mathematics
 (2011)

---

## Abstract

The objective of this study is to present a brief survey on the geometrically nonlinear free vibrations of
the Bernoulli-Euler, the Rayleigh, shear, and the Timoshenko beams with simple end conditions using the Homotopy
Analysis Method (HAM). Expressions for the natural frequencies, the transverse deflection, postbuckling
load-deflection relation to, and critical buckling load are presented. The results of nonlinear analysis are
validated with the published results, and excellent agreement is observed. The effects of some parameters,
such as slender ratio, the rotary inertia, and the shear deformation, are examined as other parameters are
fixed.


---

## Body

## 1. Introduction

Many structures, such as bridges, buildings, and spacecraft arms can be modeled as flexible beams. Most phenomena in these structures are essentially nonlinear and are described by nonlinear equations. Therefore, the study of nonlinear effects on vibration analysis of beam structures is of great importance for engineers working in the field of aerospace, mechanical, and civil engineering.Although the nonlinear free vibrations of classical Bernoulli-Euler’s beams have been investigated by many researchers [1–5], only a few publications have been devoted to include the effects of shear deformation and rotary inertia [6, 7].Azrar et al. [1] have developed a semianalytical approach to the nonlinear dynamic response problem based on Lagrange’s principle and the harmonic balance method. An analytical method for determining the vibration modes of geometrically nonlinear beams under various edge conditions has been presented by Qaisi [2]. Nonlinear modal analysis approach based on invariant manifold method is utilized to obtain the nonlinear normal modes of a clamped-clamped beam for large amplitude displacements by Xie et al. [3]. Singh et al. [4] have developed a spline-based differential quadrature method based on sextic cardinal spline functions. Guo and Zhong [5] have investigated nonlinear vibrations of thin beams based on sextic cardinal spline functions, a spline-based differential quadrature method.Rao et al. [6] have studied the large-amplitude free vibrations of slender beams using the continuum and finite element methods. The large-amplitude free-vibration behavior of the beams with central point-concentrated masses has investigated by Rao et al. [7].The main purpose of this survey is to present the analytical expressions for geometrically nonlinear vibration of Bernoulli-Euler’s, Rayleigh’s, shear, and Timoshenko’s beams with simple end conditions. First it is obtained nonlinear differential equations of motion on the basis of the four engineering beam theories. It is assumed that only the fundamental mode is excited and the governing nonlinear partial differential equation is reduced to a single nonlinear ordinary differential equation. These equations have been solved analytically in time domain using HAM. Finally, some numerical examples are shown for the nonslender and slender beams to signify the differences among the beam models.


## 2. Problem Formulation

### 2.1. Timoshenko’s Beam Model

Consider a straight Timoshenko’s beam of lengthL, a uniform cross-sectional area A, the mass per unit length of m, the moment of area of the cross-section I, Young's modulus E, and shear modulus G that is subjected to an axial force of magnitude N̂ as shown in Figure 1. We assume that the beam is made of a homogenous and isotropic material.Figure 1
A schematic of a simple ends’ beam.The kinetic energyT and the potential energy U of the vibrating beam can be written as(2.1)T=12∫0Lm(ŵ,t̂2+û,t̂2+r2ψ,t̂2)dx̂,U=12∫0L(Mψ,x̂+Vγ+N̂x̂ε)dx̂=12∫0L(EIψ,x̂2+kAGγ2+EAε2)dx̂,
where γ  represents the shear angle (γ=ŵ,x̂-ψ), û,ŵ, and ψ are axial, transverse displacements, and the cross-section rotation due to the bending moment, x̂ is the axial coordinate of the beam, r is the radius of gyration (I/A), and k is beam’s cross-sectional shape factor. Also M, V, and N̂x̂ are the bending moment, shear, and axial forces, respectively. Comma denotes differentiate with respect to x̂ or t̂.The strain-displacement relation of the Von-Karman type is used:(2.2)ε=û,x̂+ŵ,x̂22.
Applying Hamilton's principle, the governing equations of motion and boundary conditions are obtained as follows:(2.3)mŵ,t̂t̂-[kAG(ŵ,x̂-ψ)],x̂-[EA(û,x̂+ŵ,x̂22)ŵ,x̂],x̂=0,(2.4)mû,t̂t̂-[EA(û,x̂+ŵ,x̂22)ŵ,x̂],x̂=0,(2.5)mr2ψ,t̂t̂-kAG(ŵ,x̂-ψ)-(EIψ,x̂),x̂=0,(2.6)(Mδψ)0L=0,(N̂x̂δû)0L=0,(Vδŵ+N̂x̂ŵ,x̂δŵ)0L=0.
Upon neglecting the axial inertia, (2.3)-(2.4) can be combined to get the first type of differential equations which is coupled as follows: (2.7)mŵ,t̂t̂-[kAG(ŵ,x̂-ψ)],x̂+N̂x̂w,x̂x̂=0,(2.8)mr2ψ,t̂t̂-kAG(ŵ,x̂-ψ)-(EIψ,x̂),x̂=0,
where the axial force N̂x̂ is given by(2.9)N̂x̂=N̂-EA2L∫0Lŵ,x̂2dx̂
In the next step, (2.7)-(2.8) can be combined to arrive at the following decoupled equation of motion:(2.10)EIŵ,x̂x̂x̂x̂+mŵ,t̂t̂-(mr2+mEIkAG)ŵ,x̂x̂t̂t̂+mr2kAGŵ,t̂t̂t̂t̂+Nx̂(ŵ,x̂x̂+mr2kAGŵ,x̂x̂t̂t̂-EIkAGŵ,x̂x̂x̂x̂)+mr2kAGN̂x̂,t̂t̂ŵ,x̂x̂=0.
It may be noted that the equation of motion governing the rotation of the cross-section ψ  can be obtained using the same procedures to derive an equation similar to (2.10) with w being replaced byψ.


### 2.2. Shear Beam Model

This model neglects the effect of rotary inertia in Timoshenko’s beam model. Therefore, the kinetic energy can only be rewritten as(2.11)T=12∫0Lm(ŵ,t̂2+û,t̂2)dx̂.
Applying the same process mentioned in Timoshenko’s beam model, the decoupled equation of motion is given by(2.12)EIŵ,x̂x̂x̂x̂+mŵ,t̂t̂-mEIkAGŵ,x̂x̂t̂t̂+N̂x̂(ŵ,x̂x̂-EIkAGŵ,x̂x̂x̂x̂)=0,
and the boundary conditions are as listed in (2.6).


### 2.3. Rayleigh’s Beam Model

In the Rayleigh beam model, the effect of shear deformation is neglected in Timoshenko’s beam model. Unlike the Timoshenko and Shear beam models, there is only one dependent variable for the Rayleigh beam. Therefore, the potential and kinetic energies can be rewritten as(2.13)T=12∫0Lm(ŵ,t̂2+û,t̂2+r2ŵ,x̂t̂2)dx̂,U=12∫0L(Mŵ,x̂x̂+N̂x̂ε)dx̂=12∫0L(EIŵ,x̂x̂2+EAε2)dx̂.
Again following the same procedure outline before, the equation of motion in this case is obtained as(2.14)EIŵ,x̂x̂x̂x̂+mŵ,t̂t̂-mr2ŵ,x̂x̂t̂t̂+N̂x̂ŵx̂x̂=0.
The general form of any kind of boundary conditions can be listed as(2.15)(N̂x̂δû)0L=0,(Vδŵ+N̂x̂ŵ,x̂δŵ+Mδŵ,x̂)0L=0.


### 2.4. Bernoulli-Euler’s Beam Model

In this model both effects of shear deformation and rotary inertia are neglected in the Timoshenko beam model. Like the Rayleigh beam model, there is one dependent variable in this model. Based on above assumptions, the potential and kinetic energies in this case are(2.16)T=12∫0Lm(ŵ,t̂2+û,t̂2)dx̂,U=12∫0L(Mŵ,x̂x̂+N̂x̂ε)dx̂=12∫0L(EIŵ,x̂x̂2+EAε2)dx̂.
Similarly, the equation of motion is(2.17)EIŵ,x̂x̂x̂x̂+mŵ,t̂t̂+N̂x̂ŵ,x̂x̂=0,
and the boundary conditions are as listed in (2.15).


## 2.1. Timoshenko’s Beam Model

Consider a straight Timoshenko’s beam of lengthL, a uniform cross-sectional area A, the mass per unit length of m, the moment of area of the cross-section I, Young's modulus E, and shear modulus G that is subjected to an axial force of magnitude N̂ as shown in Figure 1. We assume that the beam is made of a homogenous and isotropic material.Figure 1
A schematic of a simple ends’ beam.The kinetic energyT and the potential energy U of the vibrating beam can be written as(2.1)T=12∫0Lm(ŵ,t̂2+û,t̂2+r2ψ,t̂2)dx̂,U=12∫0L(Mψ,x̂+Vγ+N̂x̂ε)dx̂=12∫0L(EIψ,x̂2+kAGγ2+EAε2)dx̂,
where γ  represents the shear angle (γ=ŵ,x̂-ψ), û,ŵ, and ψ are axial, transverse displacements, and the cross-section rotation due to the bending moment, x̂ is the axial coordinate of the beam, r is the radius of gyration (I/A), and k is beam’s cross-sectional shape factor. Also M, V, and N̂x̂ are the bending moment, shear, and axial forces, respectively. Comma denotes differentiate with respect to x̂ or t̂.The strain-displacement relation of the Von-Karman type is used:(2.2)ε=û,x̂+ŵ,x̂22.
Applying Hamilton's principle, the governing equations of motion and boundary conditions are obtained as follows:(2.3)mŵ,t̂t̂-[kAG(ŵ,x̂-ψ)],x̂-[EA(û,x̂+ŵ,x̂22)ŵ,x̂],x̂=0,(2.4)mû,t̂t̂-[EA(û,x̂+ŵ,x̂22)ŵ,x̂],x̂=0,(2.5)mr2ψ,t̂t̂-kAG(ŵ,x̂-ψ)-(EIψ,x̂),x̂=0,(2.6)(Mδψ)0L=0,(N̂x̂δû)0L=0,(Vδŵ+N̂x̂ŵ,x̂δŵ)0L=0.
Upon neglecting the axial inertia, (2.3)-(2.4) can be combined to get the first type of differential equations which is coupled as follows: (2.7)mŵ,t̂t̂-[kAG(ŵ,x̂-ψ)],x̂+N̂x̂w,x̂x̂=0,(2.8)mr2ψ,t̂t̂-kAG(ŵ,x̂-ψ)-(EIψ,x̂),x̂=0,
where the axial force N̂x̂ is given by(2.9)N̂x̂=N̂-EA2L∫0Lŵ,x̂2dx̂
In the next step, (2.7)-(2.8) can be combined to arrive at the following decoupled equation of motion:(2.10)EIŵ,x̂x̂x̂x̂+mŵ,t̂t̂-(mr2+mEIkAG)ŵ,x̂x̂t̂t̂+mr2kAGŵ,t̂t̂t̂t̂+Nx̂(ŵ,x̂x̂+mr2kAGŵ,x̂x̂t̂t̂-EIkAGŵ,x̂x̂x̂x̂)+mr2kAGN̂x̂,t̂t̂ŵ,x̂x̂=0.
It may be noted that the equation of motion governing the rotation of the cross-section ψ  can be obtained using the same procedures to derive an equation similar to (2.10) with w being replaced byψ.


## 2.2. Shear Beam Model

This model neglects the effect of rotary inertia in Timoshenko’s beam model. Therefore, the kinetic energy can only be rewritten as(2.11)T=12∫0Lm(ŵ,t̂2+û,t̂2)dx̂.
Applying the same process mentioned in Timoshenko’s beam model, the decoupled equation of motion is given by(2.12)EIŵ,x̂x̂x̂x̂+mŵ,t̂t̂-mEIkAGŵ,x̂x̂t̂t̂+N̂x̂(ŵ,x̂x̂-EIkAGŵ,x̂x̂x̂x̂)=0,
and the boundary conditions are as listed in (2.6).


## 2.3. Rayleigh’s Beam Model

In the Rayleigh beam model, the effect of shear deformation is neglected in Timoshenko’s beam model. Unlike the Timoshenko and Shear beam models, there is only one dependent variable for the Rayleigh beam. Therefore, the potential and kinetic energies can be rewritten as(2.13)T=12∫0Lm(ŵ,t̂2+û,t̂2+r2ŵ,x̂t̂2)dx̂,U=12∫0L(Mŵ,x̂x̂+N̂x̂ε)dx̂=12∫0L(EIŵ,x̂x̂2+EAε2)dx̂.
Again following the same procedure outline before, the equation of motion in this case is obtained as(2.14)EIŵ,x̂x̂x̂x̂+mŵ,t̂t̂-mr2ŵ,x̂x̂t̂t̂+N̂x̂ŵx̂x̂=0.
The general form of any kind of boundary conditions can be listed as(2.15)(N̂x̂δû)0L=0,(Vδŵ+N̂x̂ŵ,x̂δŵ+Mδŵ,x̂)0L=0.


## 2.4. Bernoulli-Euler’s Beam Model

In this model both effects of shear deformation and rotary inertia are neglected in the Timoshenko beam model. Like the Rayleigh beam model, there is one dependent variable in this model. Based on above assumptions, the potential and kinetic energies in this case are(2.16)T=12∫0Lm(ŵ,t̂2+û,t̂2)dx̂,U=12∫0L(Mŵ,x̂x̂+N̂x̂ε)dx̂=12∫0L(EIŵ,x̂x̂2+EAε2)dx̂.
Similarly, the equation of motion is(2.17)EIŵ,x̂x̂x̂x̂+mŵ,t̂t̂+N̂x̂ŵ,x̂x̂=0,
and the boundary conditions are as listed in (2.15).


## 3. Method of Solution

### 3.1. Homotopy Analysis Method (HAM)

For convenience of the reader, we will first present a brief description of the HAM [8]. Consider the following nonlinear homogeneous differential equations: (3.1)N[W(t)]=0,
where N are nonlinear operators, t denotes the independent variable, and W(t) are unknown functions.Liao constructed the so-called zero-order deformation equation as follows [8]:(3.2)(1-q)L[ϕ(t;q)-W0(t)]=qℏh(t)N[ϕ(t;q)],
where q∈[0,1] is an embedding parameter,(ℏ,h(t)) are nonzero auxiliary parameter and function; respectively, L is an auxiliary linear operator, W0(t) are initial guesses of W(t), and ϕ(t;q) are unknown functions. As q increases from 0 to 1, the ϕ(t;q) varies from the initial approximation to the exact solution. In other words,ϕ(t;0)=W0(t) andϕ(t;1)=W(t).Differentiating once more from (3.2) with respect to the embedding parameter q and then setting q=0, we obtain the first-order deformation equation as(3.3)L[W1(t)]=ℏh(t)N[ϕ(t;q)]q=0,
which gives the first-order approximation for the W(t) [8]. The higher-order approximation of the solution can be obtained by calculating the mth-order (m>1) deformation equation which can be calculated by differentiating (3.3) m times with respect to the q, then setting q=0, and finally dividing them by m! [8].To solve (3.2), we employ Taylor’s expansion series for ϕ(t;q) as(3.4)ϕ(t;q)=ϕ(t;0)+∑m=1∞1m!∂mϕ(t;q)∂qm|q=0qm=W0(t)+∑m=1∞Wm(t)qm,
where Wm(t) is called the mth-order derivative of unknown function ϕ(t;q).


### 3.2. Application of the HAM in Beam Vibrations

#### 3.2.1. Timoshenko’s Beam Model

For convenience, the following nondimensional variables are used:(3.5)x=x̂L,w=ŵr,t=t̂EI(mL4).
In order to obtain the nonlinear natural frequency, we start to apply the Galerkin method on which the deflection for the beam with simply supported boundary conditions is expressed as(3.6)w(x,t)=W(t)sin⁡⁡(πx)
As a result (2.10) can be written as follows:(3.7)W,tttt+[η4λ+η2π2(1+1λ-N)+π44]W,tt+η4π2[π2λ-N(π2+η2λ)]W+η2π44[η2λ+π2]W3+π42(WW,t2+W2W,tt)=0,
in which(3.8)η=Lr,λ=EkG,N=N̂EA.
The beam centroid is subjected to the following initial conditions:(3.9)W(0)=a,W,t(0)=0,W,tt(0)=0,W,ttt(0)=0.
Note that the free vibration of a system without damping is expressed by a kind of periodic function which has to satisfy the initial conditions. For our case this function simply is(3.10)cos⁡⁡(nωNLt),n=1,2,3,….
It is well known that for a Timoshenko’s beam with simple ends there are two natural frequencies. In other words, there are two sinusoidal modes of different frequencies corresponding to the same spatial mode. This is the most important difference between Timoshenko’s model and the other beam models with simple end conditions. This phenomenon is also reported by many researchers [9–11]. Thus, the initial guess for W(t) is as follows:(3.11)W0(t)=A1cos⁡⁡(ωNL1t)+A2cos⁡⁡(ωNL2t).
To satisfy the initial conditions given in (3.9), we apply those conditions into (3.11) which yields the following expressions for the A1 and A2 as(3.12)A1=a1-(ωNL1/ωNL2)2,A2=a1-(ωNL2/ωNL1)2.
Now we use (3.11) in conjunction with HAM to get the natural frequencies and the beam deflection. To construct the homotopy function, the auxiliary linear operator is selected as(3.13)L[ϕ(t;q)]=∂4ϕ(t;q)∂t4+(ωNL12+ωNL22)∂2ϕ(t;q)∂t2+(ωNL1ωNL2)2.
For convenience and simplify, from (3.7), the nonlinear operator is defined as(3.14)N[ϕ(t;q)]=∂4ϕ(t;q)∂t4+α1∂2ϕ(t;q)∂t2+α2ϕ(t;q)+α3ϕ(t;q)3+α4(ϕ(t;q)2∂2ϕ(t;q)∂t2+ϕ(t;q)(∂ϕ(t;q)∂t)2),
where(3.15)α1=η4λ+η2π2(1+1λ-N)+π44,α2=η4π2[π2λ-N(π2+η2λ)],α3=η2π44[η2λ+π2],α4=π42.
In our case, to obtain the first-order approximation, (3.3) for the function of ϕ(t;q) can be expressed as(3.16)L[W1(t)]=ℏh(t)N[ϕ(t;q)]|q=0,W1(0)=0,W1,t(0)=0,W1,tt(0)=0,W1,ttt(0)=0.
Assuming ℏ=-1andh(t)=1 and substituting (3.11) in (3.16), (3.16) becomes as(3.17)W1,tttt+(ωNL12+ωNL22)W1,tt+(ωNL1ωNL2)2=B1cos⁡⁡(3ωNL1t)+B2cos⁡⁡(3ωNL2t)+B3cos⁡⁡((ωNL1-2ωNL2)t)+B4cos⁡⁡((ωNL1+2ωNL2)t)+B5cos⁡⁡((2ωNL1-ωNL2)t)+B6cos⁡⁡((2ωNL1+ωNL2)t)+B7cos⁡⁡(ωNL1t)+B8cos⁡⁡(ωNL2t),
where B’s are given in the appendix, equation (A.1).These solutions should obey the general form of the base function. This means that the coefficients of the secular terms (cos(ωNL1t),cos(ωNL2t)) must be zero, that is, B7=0andB8=0. This provides two algebric equations which yield the nonlinear frequencies: (3.18)(ωNL14-α1ωNL12+α2)(ωNL24-ωNL14)2+34α3a2(ωNL24+2ωNL14)-12α4a2(ωNL16+ωNL14ωNL22+ωNL12ωNL24)=0,(ωNL24-α1ωNL22+α2)(ωNL24-ωNL14)2+34α3a2(ωNL14+2ωNL24)-12α4a2(ωNL26+ωNL14ωNL22+ωNL12ωNL24)=0,
and solving (3.17), the W1(t) is obtained as follows:(3.19)W1(t)=Acos⁡(ωNL1t)+Bcos⁡⁡(ωNL2t)+D1cos⁡⁡(3ωNL1t)+D2cos⁡⁡(3ωNL2t)+D3cos⁡⁡((ωNL1-2ωNL2)t)+D4cos⁡⁡((ωNL1+2ωNL2)t)+D5cos⁡⁡((2ωNL1-ωNL2)t)+D6cos⁡⁡((2ωNL1+ωNL2)t),
where A and B are obtained by applying initial conditions. Coefficients A, B and D’s have been presented in the appendix, equation (A.2).Thus, the first-order approximation of theW(t) becomes as(3.20)W(t)=W0(t)+W1(t).


#### 3.2.2. Shear, Rayleigh’s, and Bernoulli-Euler’s Beam Models

The second-order differential equations of motion for either of shear (S), Rayleigh’s (R), or Bernoulli- Euler’s (B) beams using Galerkin’s method are the same and is(3.21)W,tt+αi1W+αi3W3=0,i=S,R,B.
Under the following initial conditions(3.22)w(0)=a,w,t(0)=0.
Using the same method described for Timoshenko’s beam model, we obtain the nonlinear natural frequency and the deflection of the beam as follows.Nonlinear natural frequency:(3.23)ωNL=123αi3a3+4αi1.Deflection of the beam,(3.24)W(t)=(αi3a332ωNL2+a)cos⁡⁡(ωNLt)-αi3a332ωNL2cos⁡⁡(3ωNLt),
where αij,(i=S,R,B,j=1,3) are given in the appendix equation (A.3).


## 3.1. Homotopy Analysis Method (HAM)

For convenience of the reader, we will first present a brief description of the HAM [8]. Consider the following nonlinear homogeneous differential equations: (3.1)N[W(t)]=0,
where N are nonlinear operators, t denotes the independent variable, and W(t) are unknown functions.Liao constructed the so-called zero-order deformation equation as follows [8]:(3.2)(1-q)L[ϕ(t;q)-W0(t)]=qℏh(t)N[ϕ(t;q)],
where q∈[0,1] is an embedding parameter,(ℏ,h(t)) are nonzero auxiliary parameter and function; respectively, L is an auxiliary linear operator, W0(t) are initial guesses of W(t), and ϕ(t;q) are unknown functions. As q increases from 0 to 1, the ϕ(t;q) varies from the initial approximation to the exact solution. In other words,ϕ(t;0)=W0(t) andϕ(t;1)=W(t).Differentiating once more from (3.2) with respect to the embedding parameter q and then setting q=0, we obtain the first-order deformation equation as(3.3)L[W1(t)]=ℏh(t)N[ϕ(t;q)]q=0,
which gives the first-order approximation for the W(t) [8]. The higher-order approximation of the solution can be obtained by calculating the mth-order (m>1) deformation equation which can be calculated by differentiating (3.3) m times with respect to the q, then setting q=0, and finally dividing them by m! [8].To solve (3.2), we employ Taylor’s expansion series for ϕ(t;q) as(3.4)ϕ(t;q)=ϕ(t;0)+∑m=1∞1m!∂mϕ(t;q)∂qm|q=0qm=W0(t)+∑m=1∞Wm(t)qm,
where Wm(t) is called the mth-order derivative of unknown function ϕ(t;q).


## 3.2. Application of the HAM in Beam Vibrations

### 3.2.1. Timoshenko’s Beam Model

For convenience, the following nondimensional variables are used:(3.5)x=x̂L,w=ŵr,t=t̂EI(mL4).
In order to obtain the nonlinear natural frequency, we start to apply the Galerkin method on which the deflection for the beam with simply supported boundary conditions is expressed as(3.6)w(x,t)=W(t)sin⁡⁡(πx)
As a result (2.10) can be written as follows:(3.7)W,tttt+[η4λ+η2π2(1+1λ-N)+π44]W,tt+η4π2[π2λ-N(π2+η2λ)]W+η2π44[η2λ+π2]W3+π42(WW,t2+W2W,tt)=0,
in which(3.8)η=Lr,λ=EkG,N=N̂EA.
The beam centroid is subjected to the following initial conditions:(3.9)W(0)=a,W,t(0)=0,W,tt(0)=0,W,ttt(0)=0.
Note that the free vibration of a system without damping is expressed by a kind of periodic function which has to satisfy the initial conditions. For our case this function simply is(3.10)cos⁡⁡(nωNLt),n=1,2,3,….
It is well known that for a Timoshenko’s beam with simple ends there are two natural frequencies. In other words, there are two sinusoidal modes of different frequencies corresponding to the same spatial mode. This is the most important difference between Timoshenko’s model and the other beam models with simple end conditions. This phenomenon is also reported by many researchers [9–11]. Thus, the initial guess for W(t) is as follows:(3.11)W0(t)=A1cos⁡⁡(ωNL1t)+A2cos⁡⁡(ωNL2t).
To satisfy the initial conditions given in (3.9), we apply those conditions into (3.11) which yields the following expressions for the A1 and A2 as(3.12)A1=a1-(ωNL1/ωNL2)2,A2=a1-(ωNL2/ωNL1)2.
Now we use (3.11) in conjunction with HAM to get the natural frequencies and the beam deflection. To construct the homotopy function, the auxiliary linear operator is selected as(3.13)L[ϕ(t;q)]=∂4ϕ(t;q)∂t4+(ωNL12+ωNL22)∂2ϕ(t;q)∂t2+(ωNL1ωNL2)2.
For convenience and simplify, from (3.7), the nonlinear operator is defined as(3.14)N[ϕ(t;q)]=∂4ϕ(t;q)∂t4+α1∂2ϕ(t;q)∂t2+α2ϕ(t;q)+α3ϕ(t;q)3+α4(ϕ(t;q)2∂2ϕ(t;q)∂t2+ϕ(t;q)(∂ϕ(t;q)∂t)2),
where(3.15)α1=η4λ+η2π2(1+1λ-N)+π44,α2=η4π2[π2λ-N(π2+η2λ)],α3=η2π44[η2λ+π2],α4=π42.
In our case, to obtain the first-order approximation, (3.3) for the function of ϕ(t;q) can be expressed as(3.16)L[W1(t)]=ℏh(t)N[ϕ(t;q)]|q=0,W1(0)=0,W1,t(0)=0,W1,tt(0)=0,W1,ttt(0)=0.
Assuming ℏ=-1andh(t)=1 and substituting (3.11) in (3.16), (3.16) becomes as(3.17)W1,tttt+(ωNL12+ωNL22)W1,tt+(ωNL1ωNL2)2=B1cos⁡⁡(3ωNL1t)+B2cos⁡⁡(3ωNL2t)+B3cos⁡⁡((ωNL1-2ωNL2)t)+B4cos⁡⁡((ωNL1+2ωNL2)t)+B5cos⁡⁡((2ωNL1-ωNL2)t)+B6cos⁡⁡((2ωNL1+ωNL2)t)+B7cos⁡⁡(ωNL1t)+B8cos⁡⁡(ωNL2t),
where B’s are given in the appendix, equation (A.1).These solutions should obey the general form of the base function. This means that the coefficients of the secular terms (cos(ωNL1t),cos(ωNL2t)) must be zero, that is, B7=0andB8=0. This provides two algebric equations which yield the nonlinear frequencies: (3.18)(ωNL14-α1ωNL12+α2)(ωNL24-ωNL14)2+34α3a2(ωNL24+2ωNL14)-12α4a2(ωNL16+ωNL14ωNL22+ωNL12ωNL24)=0,(ωNL24-α1ωNL22+α2)(ωNL24-ωNL14)2+34α3a2(ωNL14+2ωNL24)-12α4a2(ωNL26+ωNL14ωNL22+ωNL12ωNL24)=0,
and solving (3.17), the W1(t) is obtained as follows:(3.19)W1(t)=Acos⁡(ωNL1t)+Bcos⁡⁡(ωNL2t)+D1cos⁡⁡(3ωNL1t)+D2cos⁡⁡(3ωNL2t)+D3cos⁡⁡((ωNL1-2ωNL2)t)+D4cos⁡⁡((ωNL1+2ωNL2)t)+D5cos⁡⁡((2ωNL1-ωNL2)t)+D6cos⁡⁡((2ωNL1+ωNL2)t),
where A and B are obtained by applying initial conditions. Coefficients A, B and D’s have been presented in the appendix, equation (A.2).Thus, the first-order approximation of theW(t) becomes as(3.20)W(t)=W0(t)+W1(t).


### 3.2.2. Shear, Rayleigh’s, and Bernoulli-Euler’s Beam Models

The second-order differential equations of motion for either of shear (S), Rayleigh’s (R), or Bernoulli- Euler’s (B) beams using Galerkin’s method are the same and is(3.21)W,tt+αi1W+αi3W3=0,i=S,R,B.
Under the following initial conditions(3.22)w(0)=a,w,t(0)=0.
Using the same method described for Timoshenko’s beam model, we obtain the nonlinear natural frequency and the deflection of the beam as follows.Nonlinear natural frequency:(3.23)ωNL=123αi3a3+4αi1.Deflection of the beam,(3.24)W(t)=(αi3a332ωNL2+a)cos⁡⁡(ωNLt)-αi3a332ωNL2cos⁡⁡(3ωNLt),
where αij,(i=S,R,B,j=1,3) are given in the appendix equation (A.3).


## 3.2.1. Timoshenko’s Beam Model

For convenience, the following nondimensional variables are used:(3.5)x=x̂L,w=ŵr,t=t̂EI(mL4).
In order to obtain the nonlinear natural frequency, we start to apply the Galerkin method on which the deflection for the beam with simply supported boundary conditions is expressed as(3.6)w(x,t)=W(t)sin⁡⁡(πx)
As a result (2.10) can be written as follows:(3.7)W,tttt+[η4λ+η2π2(1+1λ-N)+π44]W,tt+η4π2[π2λ-N(π2+η2λ)]W+η2π44[η2λ+π2]W3+π42(WW,t2+W2W,tt)=0,
in which(3.8)η=Lr,λ=EkG,N=N̂EA.
The beam centroid is subjected to the following initial conditions:(3.9)W(0)=a,W,t(0)=0,W,tt(0)=0,W,ttt(0)=0.
Note that the free vibration of a system without damping is expressed by a kind of periodic function which has to satisfy the initial conditions. For our case this function simply is(3.10)cos⁡⁡(nωNLt),n=1,2,3,….
It is well known that for a Timoshenko’s beam with simple ends there are two natural frequencies. In other words, there are two sinusoidal modes of different frequencies corresponding to the same spatial mode. This is the most important difference between Timoshenko’s model and the other beam models with simple end conditions. This phenomenon is also reported by many researchers [9–11]. Thus, the initial guess for W(t) is as follows:(3.11)W0(t)=A1cos⁡⁡(ωNL1t)+A2cos⁡⁡(ωNL2t).
To satisfy the initial conditions given in (3.9), we apply those conditions into (3.11) which yields the following expressions for the A1 and A2 as(3.12)A1=a1-(ωNL1/ωNL2)2,A2=a1-(ωNL2/ωNL1)2.
Now we use (3.11) in conjunction with HAM to get the natural frequencies and the beam deflection. To construct the homotopy function, the auxiliary linear operator is selected as(3.13)L[ϕ(t;q)]=∂4ϕ(t;q)∂t4+(ωNL12+ωNL22)∂2ϕ(t;q)∂t2+(ωNL1ωNL2)2.
For convenience and simplify, from (3.7), the nonlinear operator is defined as(3.14)N[ϕ(t;q)]=∂4ϕ(t;q)∂t4+α1∂2ϕ(t;q)∂t2+α2ϕ(t;q)+α3ϕ(t;q)3+α4(ϕ(t;q)2∂2ϕ(t;q)∂t2+ϕ(t;q)(∂ϕ(t;q)∂t)2),
where(3.15)α1=η4λ+η2π2(1+1λ-N)+π44,α2=η4π2[π2λ-N(π2+η2λ)],α3=η2π44[η2λ+π2],α4=π42.
In our case, to obtain the first-order approximation, (3.3) for the function of ϕ(t;q) can be expressed as(3.16)L[W1(t)]=ℏh(t)N[ϕ(t;q)]|q=0,W1(0)=0,W1,t(0)=0,W1,tt(0)=0,W1,ttt(0)=0.
Assuming ℏ=-1andh(t)=1 and substituting (3.11) in (3.16), (3.16) becomes as(3.17)W1,tttt+(ωNL12+ωNL22)W1,tt+(ωNL1ωNL2)2=B1cos⁡⁡(3ωNL1t)+B2cos⁡⁡(3ωNL2t)+B3cos⁡⁡((ωNL1-2ωNL2)t)+B4cos⁡⁡((ωNL1+2ωNL2)t)+B5cos⁡⁡((2ωNL1-ωNL2)t)+B6cos⁡⁡((2ωNL1+ωNL2)t)+B7cos⁡⁡(ωNL1t)+B8cos⁡⁡(ωNL2t),
where B’s are given in the appendix, equation (A.1).These solutions should obey the general form of the base function. This means that the coefficients of the secular terms (cos(ωNL1t),cos(ωNL2t)) must be zero, that is, B7=0andB8=0. This provides two algebric equations which yield the nonlinear frequencies: (3.18)(ωNL14-α1ωNL12+α2)(ωNL24-ωNL14)2+34α3a2(ωNL24+2ωNL14)-12α4a2(ωNL16+ωNL14ωNL22+ωNL12ωNL24)=0,(ωNL24-α1ωNL22+α2)(ωNL24-ωNL14)2+34α3a2(ωNL14+2ωNL24)-12α4a2(ωNL26+ωNL14ωNL22+ωNL12ωNL24)=0,
and solving (3.17), the W1(t) is obtained as follows:(3.19)W1(t)=Acos⁡(ωNL1t)+Bcos⁡⁡(ωNL2t)+D1cos⁡⁡(3ωNL1t)+D2cos⁡⁡(3ωNL2t)+D3cos⁡⁡((ωNL1-2ωNL2)t)+D4cos⁡⁡((ωNL1+2ωNL2)t)+D5cos⁡⁡((2ωNL1-ωNL2)t)+D6cos⁡⁡((2ωNL1+ωNL2)t),
where A and B are obtained by applying initial conditions. Coefficients A, B and D’s have been presented in the appendix, equation (A.2).Thus, the first-order approximation of theW(t) becomes as(3.20)W(t)=W0(t)+W1(t).


## 3.2.2. Shear, Rayleigh’s, and Bernoulli-Euler’s Beam Models

The second-order differential equations of motion for either of shear (S), Rayleigh’s (R), or Bernoulli- Euler’s (B) beams using Galerkin’s method are the same and is(3.21)W,tt+αi1W+αi3W3=0,i=S,R,B.
Under the following initial conditions(3.22)w(0)=a,w,t(0)=0.
Using the same method described for Timoshenko’s beam model, we obtain the nonlinear natural frequency and the deflection of the beam as follows.Nonlinear natural frequency:(3.23)ωNL=123αi3a3+4αi1.Deflection of the beam,(3.24)W(t)=(αi3a332ωNL2+a)cos⁡⁡(ωNLt)-αi3a332ωNL2cos⁡⁡(3ωNLt),
where αij,(i=S,R,B,j=1,3) are given in the appendix equation (A.3).


## 4. Postbuckling Load-Deflection Relation and Critical Buckling Load

In order to obtain the postbuckling load-deflection relation one can set all time derivative terms in (3.7) and (3.21) equal to zero which yields the following.Rayleigh’s and Bernoulli-Euler’s beam models,(4.1)N=π24η2[4+W2].Shear and Timoshenko’s beam models,(4.2)N=π24η2(η2+λπ2)[4η2+(η2+λπ2)W2].
Neglecting the contribution of W  in the above equations, the critical buckling load can be determined as follows:Rayleigh and Bernoulli-Euler Beam Models:(4.3)Ncr=(πη)2.Shear and Timoshenko’s beam models:(4.4)Ncr=π2η2+λπ2.
As it is seen the postbuckling and critical buckling load predicted by Timoshenko model is the same as the one from shear model. The same trend exists between the Rayleigh and the Bernoulli-Euler beam models.


## 5. Results and Discussions

In all computational cases, the shear coefficient is taken ask=5/6 and the axial force N̂=0 unless mentioned otherwise.In order to demonstrate the effectiveness of the HAM, the natural frequencies of a simply supported beam are compared with other existing results in literatures. To do this the nondimensional nonlinear frequencies are calculated, and the results are listed for differentη and a in the Table 1.Table 1
Comparison of the nondimensional nonlinear frequency(ωNL).
aηReference [6]EBBM2RBM3SBM4TBM5(ωNL1)(ωNL2)0.4108.5458(4)110.0179.55618.82008.545267.412259.7198(3)9.93849.79069.7200369.04509.9391(3)9.99689.95859.93931442.81009.9968(3)10.01210.0029.99705737.30.8109.0018(6)10.4459.96489.30379.010367.4382510.1494(5)10.36410.22810.155369.055010.3650(4)10.42410.38910.3691442.810010.4217(4)10.44010.43110.4265737.31.0109.3200(7)10.75510.2619.65069.343667.4582510.4563(5)10.67110.54510.469369.055010.6702(5)10.73410.70110.6801442.810010.7265(5)10.75010.74210.7365737.32.01011.4868(16)13.05612.45612.16311.74767.6232512.6458(8)12.95412.88412.790369.085012.8689(7)13.03113.01212.9871442.910012.9279(7)13.05013.04513.0395737.33.01014.1558(47)16.18015.43615.46814.87667.9072515.4532(14)16.05416.04115.922369.125015.7178(11)16.14816.14416.1131442.910015.7889(11)16.17216.17116.1635737.31Numbers in parentheses indicate the number of iterations required to achieve an accuracy of 10-6 in the evaluation of the nonlinear- to-linear frequency ratio.2EBBM: Bernoulli-Euler beam model.3RBM: Rayleigh beam model.4SBM: shear beam model.5TBM: Timoshenko beam model.It can be observed that there is an excellent agreement between the results obtained from HAM and those reported by Rao et al. [6]. It should be noted that only the first natural frequency of a Timoshenko’s beam is reported in [6] and in a rather numerical way using FE method. It is clear that Bernoulli-Euler’s beam predicts the nonlinear natural frequency higher than the one obtained by other theories. Also the nonlinear natural frequency obtained by Timoshenko’s beam theory is lower than the one obtained by other theories.The convergence and accuracy of the presented analytical solution is validated against variation of nondimensional amplitude ratio of the beam center versust using the Runge-Kutta method. Figure 2 illustrates comparison between these results. They are found to be in a very good agreement.Figure 2
Variation of the nondimensional amplitude of the beam center versust.Now, after being satisfied in obtaining results given by this method, in the next step we try to investigate the effect of different parameters on the natural frequencies of four above-mentioned beam models.The variation of the nonlinear frequency versus nondimensional amplitude for all different models and for two different slenderness ratios (η) using HAM is shown in Figure 3. For thin beams and based on the three employed beam models, namely, Bernoulli-Euler’s, Rayleigh’s and shear beams, almost no difference is seen for different nondimensional amplitude. However, due to the nature of Timoshenko’s simply supported beam, there are two ωNL variations. In the region 0<a<5  as it is seen from Figure 3(a), the ωNL1 represents a linear fashion variation whereas the ωNL2 variation remains almost unchanged. In the region a>20 as it is seen from the same figure, variation of ωNL1 follows a rather constant trend and ωNL2 follows a linear fashion variation. For the intermediate region, that is, 5<a<20 frequencies of both ωNL1 and ωNL2 keep changing in a rather nonlinear fashion. As can be seen from Figure 3(a), for thick beams this interval occurs in small amplitude. In other words, for thin and thick beams though in each amplitude both nonlinear frequencies (ωNL1,ωNL2) are involved in beam's response; however, one of them is mostly affected by variation of the nondimensional amplitude.Nonlinear frequency versus nondimensional amplitude (- - -) Bernoulli-Euler’s,(⋯) Rayleigh’s, (– - –) shear, and (—) Timoshenko’s Beam Models.
(a)(b)Figure4 shows nondimensional natural frequency ((ωNL1) for Timoshenko’s beam model) versus small nondimensional amplitude for different slenderness ratio and different beam models. It can be seen that the natural nonlinear frequency obtained from four different beam models yields the same value provided that the beam's slenderness ratio at least 100. This finding is also reported by [11] in which only the linear analysis was conducted.Nonlinear frequency versus nondimensional amplitude for different slenderness ratios (- - -) Bernoulli-Euler’s,(⋯) Rayleigh’s, (-·-) shear, and (—) Timoshenko’s beam models.
(a)(b)(c)(d)The variations of the nonlinear frequency ratios versus nondimensional amplitude of Timoshenko beam for different slenderness ratios are depicted in Figure5. As it is seen from these figures the point of shift from flat to nonflat and vise versa of ωNL  increases by increasing the beam slenderness ratio.Nonlinear frequency versus nondimensional amplitude for different beam slenderness ratios  (- - -)η=10, (⋯)η=20, (-·-)  η=35, and (—) η=50.
(a)(b)Variation of the nondimensional postbuckling load(N) versus nondimensional amplitude is shown in Figure 6. As it is seen the postbuckling load predicted by all four different models are the same for rather a slender beam (η=50). On the other hand for the thick beams (η=10), the differences between Timoshenko and/or Shear beam models with Rayleigh’s and/or Bernoulli-Euler’s beam models are more pronounced. Furthermore, for both cases the postbuckling load increases as amplitude (a) increases.Variation of the nondimensional buckling load versus nondimensional amplitude for different beam models and different slenderness ratios (- - -) Bernoulli-Euler’s and Rayleigh’s, (—) shear and Timoshenko’s beam models.
(a)(b)Variation of the critical buckling load(Ncr) versus slenderness ratio for the problem under investigation is shown in Figure 7. It can be observed that for η>12, that is, thin beams, all four models give the same answer for the critical buckling load.Figure 7
Variation of the critical buckling load(Ncr) versus slenderness ratio for different beam models (- - -) Bernoulli-Euler’s and Rayleigh’s, (—) shear and Timoshenko’s beam models.


## 6. Conclusion

The natural frequencies, deflection, postbuckling, and critical buckling loads of four approximate models of Bernoulli-Euler’s, Rayleigh’s, shear, and Timoshenko’s beams having nonlinear behavior and simply supported end conditions have been discussed in the present survey in detail. A first-order approximation of HAM is used primarily to get the analytic expressions for the nonlinear frequencies and consequently the related deflection of neutral axis of such beam. Based on derived results, the followings are concluded.(1)
For a wide range of cases having large amplitude vibration, HAM seems to be very efficient in handling such problems.(2)
For thin and thick beams with simple end conditions, there are two nonlinear natural frequencies. Depending on the value of amplitude, only one of these two nonlinear natural frequencies is mostly affected by variation of the amplitude. There is a primary interval of amplitude in which the lower frequency varies in a linear fashion, and the higher frequency remains almost constant. Contrary to this trend, at the ending interval the frequency variations change their character; that is, the lower frequency is almost constant and the higher frequency varies in a linear fashion. Moreover, in the intermediate interval both frequencies increase in a nonlinear fashion. It is also seen that for the thick beams this intermediate range shifts to the left; that is, the primary interval shrinks further. In addition, it is shown that the point of shift from flat to nonflat part ofωNL curves and vise versa increases by increasing the beam slenderness ratio.(3)
The Bernoulli-Euler beam model predicts the nonlinear natural frequency higher than the other theories. Whereas, the nonlinear natural frequency obtained by Timoshenko’s beam theory is lower than the other theories.(4)
It can be seen that the natural nonlinear frequency obtained from four different beam models yields the same value for the beams with slenderness ratio of at least 100.(5)
The shear effect is always more dominant than the rotary effect for a given geometry and material in predicting postbuckling and critical buckling loads. Therefore, it is highly recommended to use either of the shear or the Timoshenko models for a thick beam. Nonetheless, in the postbuckling and critical buckling loads analysis and in a more general view the shear beam model conveys less complexity than the Timoshenko beam model.(6)
The postbuckling and critical buckling loads predicted by all four different models are the same for rather a slender beam. On the other hand, for the thick beams the differences between Timoshenko’s and/or shear beam models with Rayleigh’s and/or Euler-Bernoulli’s beam models are more pronounced. Furthermore, for both cases the postbuckling load(N) increases as amplitude increases.



---

*Source: 842805-2011-12-22.xml*