# A Semianalytical Method for Nonlinear  Vibration of Euler-Bernoulli Beams with  General Boundary Conditions

Authors: Jian-She Peng; Yan Liu; Jie Yang
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

This paper presents a new semianalytical approach for geometrically nonlinear vibration analysis of Euler-Bernoulli beams with different boundary conditions. The method makes use of Linstedt-Poincaré perturbation technique to transform the nonlinear governing equations into a linear differential equation system, whose solutions are then sought through the use of differential quadrature approximation in space domain and an analytical series expansion in time domain. Validation of the present method is conducted in numerical examples through direct comparisons with existing solutions, showing that the proposed semianalytical method has excellent convergence and can give very accurate results at a long time interval.


---

## Body

## 1. Introduction

Geometrically  nonlinear vibration of beams with different boundary conditions has long been a subject receiving numerous research efforts, as evidenced by many analytical and numerical studies reported in the open literature. Woinowsky-Krieger [1] used the elliptic integral function to solve the nonlinear equation of motion for simply supported beams with immovable ends. Lewandowski [2] applied the Rayleigh-Ritz technique to the nonlinear vibration of beams. A comprehensive  review in this field was given by Sathyamoorthy [3]. Sze et al.  [4] applied incremental harmonic balance method for nonlinear vibration of axially moving beams. Gadagi and Benaroya  [5] studied the dynamic response of an axially loaded tendon of a tension leg platform. Ibrahim and Somnay [6] solved the nonlinear dynamic analysis of an elastic beam isolator sliding on frictional supports. Ozkaya and Tekin [7] analyzed the nonlinear vibrations of stepped beam system under different boundary conditions. Chen et al. [8] put forward the multidimensional Lindstedt-Poincare method for nonlinear vibration of axially moving beams. It is noted that many numerical investigations published so far employed step-by-step iterative time integration schemes to calculate the dynamic deflection response which requires a very small time step size in order to achieve the response with sufficient accuracy. This is inevitably computationally expensive. Another big concern associated with these iterative schemes is the error accumulation in each time step which may cause a huge loss of accuracy and sometimes numerical stability problem of the results in later time steps. To overcome this problem, several analytical and semianalytical approaches to the nonlinear vibration of beams which do not involve step-by-step time integration have been proposed; see, for example, those by Azrar et al. [9, 10], Leung and his coworkers [11–13], Lewandowski [14, 15], Nayfeh and his associates [16–18], and Ribeiro [19], to name just a few.The differential quadrature method (DQM) is an efficient numerical approach developed by Bellman et al.  [20] to solve linear and nonlinear differential equations. This method was later introduced to the structural analysis by Bert  and his coworkers  [21] and has been widely used in many structural engineering problems, including those in the nonlinear vibration of beams, plates, and shells; see, for example, the work by Zhong and Guo [22], Yang et al. [23, 24], Manaoach and Ribeiro [25], Hsu [26], and Tomasiello [27], among many others. Based on DQM approximation, this paper proposes a new semianalytic method for the geometrically nonlinear vibration of Euler-Bernoulli beams with different boundary conditions. The proposed method overcomes the disadvantage of error accumulation that occurs in conventional numerical methods involving iterative time integration and is therefore accurate and numerically stable even for a long time interval. The present paper is structured as follows: Section 2 briefly outlines the DQM rules. Section 3 gives a detailed description of mathematic formulations of the proposed semianalytical approach, starting with the transformation of the nonlinear governing equations into a group of linear differential equations by Linstedt-Poincare perturbation procedure which is followed by the semianalytical solution process for each perturbation equation by using DQM in space domain and an analytical series in the time domain. Section 4 presents some numerical results to validate the proposed method in both accuracy and convergence aspects through direct comparisons between the present results and existing solutions. Some concluding remarks are summarized in Section 4.


## 2. Mathematical Formulations

### 2.1. Perturbation Equations

Consider an isotropic, homogeneous slender beam of lengthL with uniform cross-section area subjected to a dynamic transverse load  q̅(x,t̅). Let u̅ and w̅ be the displacements parallel to the x- and y-axes, t̅   the time, and m denote  the mass density. Based on classical Euler-Bernoulli beam theory and von-Karman nonlinear displacement-strain relationship, the geometrically nonlinear governing equations for flexural vibration of the beam can be  derived as
(2.1)∂∂x[EA̅(∂u¯∂x+12(∂w¯∂x)2)]=0,EI∂4w̅∂x4-EA̅[∂u̅∂x+12(∂w̅∂x)2]∂2w̅∂x2+m∂2w̅∂t̅2=q̅(x,t̅),
where E is the elastic modulus of the beam, and A̅ and  I are the area and the second moment of the cross section. The present analysis considers general boundary conditions, that is, the beam end may be either simply supported or clamped, with the following boundary conditions: (2.2)u̅=0,w̅=0,M̅x=0(simplysupported),u̅=0,w̅=0,dw̅dx=0(clamped),
where M̅x is the bending moment of the beam. To facilitate the solution process of  (2.1)  by using the differential quadrature method in space domain, the following quantities are introduced: (2.3)u=uη¯,w=w¯ρ,ξ=xL,t=ωt¯,q=L4q̅ρEI,q=εF,p=ωωL,ωL=(1L)2EIm,ρ2=IA̅,
where ωL and ω refer to the linear and nonlinear frequencies of the beam, respectively, and ε is the small perturbation parameter. Substituting these quantities into (2.1)   leads to dimensionless nonlinear governing equations as below (2.4)∂2u∂ξ2+∂w∂ξ∂2w∂ξ2=0,∂4w∂ξ4-[∂u∂ξ+12(∂w∂ξ)2]∂2w∂ξ2+p2∂2w∂t2=εF.Deflectionw, axial displacement u, and frequency parameter p can be expanded into series forms in terms of a small perturbation  parameter ε(2.5)w=w1ε+w3ε3+w5ε5+⋯,(2.6)u=u2ε2+u4ε4+⋯,(2.7)p2=1+ε2σ2+ε4σ4+⋯.
Inserting   (2.5)–(2.7)  into   (2.4), and equating the terms of the same power of ε yield  a series of perturbation equations (2.8)ε1:∂4w1∂ξ4+∂2w1∂t2=F,(2.9)ε2:∂2u2∂ξ2+∂w1∂ξ∂2w1∂ξ2=0,(2.10)ε3:∂4w3∂ξ4-∂u2∂ξ∂2w1∂ξ2-12(∂w1∂ξ)2∂2w1∂ξ2+∂2w3∂t2+σ2∂2w1∂t2=0,(2.11)ε4:∂2u4∂ξ2+∂w1∂ξ∂2w3∂ξ2+∂w3∂ξ∂2w1∂ξ2=0.The initial conditions considered in the present study are(2.12)w1(ξ,0)=w10,∂w1(ξ,0)∂t=0,(2.13)w3(ξ,0)=0,∂w3(ξ,0)∂t=0,
where w10 refers to initial vibration amplitude. By making use of the dimensionless quantities and perturbation series defined above, the dimensionless form of the associated boundary conditions in each perturbation can also be obtained.


### 2.2. Differential Quadrature Method

According to the DQM rule, an unknown functiong(x) and its nth order derivative with respect to x can be approximated as the weighted linear sums of its function values at a number of sampling points in the x-axis as (2.14)g(x)=∑i=1Ng(xi)Li(x),(2.15)dng(x)dxn|x=xi=∑j=1NCij(n)g(xj)(n=1,…,N-1),
where N is the total number of sampling points, Li(x) are Lagrange interpolation polynomials, and the weighting coefficients Cij(n) can be calculated from recursive formulae [20, 21].


### 2.3. Semianalytical Solution of the 1st-Order Perturbation Equation

This study is focused on the nonlinear vibration due to an initial deflection. In such a case, the right-hand-side term in   (2.8)  F=0. Hence, the first-order perturbation   becomes(2.16)∂4w1∂ξ4+∂2w1∂t2=0.
The solution of   (2.16) can be readily obtained by separation of variables(2.17)w1(ξ,t)=∑m=0NsT1m(t)ϕ1m(ξ),
where Ns is the number of terms in the truncated time series,                (2.18)T1m(t)=cos(pmt)
is the solution in the time domain for the first-order perturbation equation, and  pm and δ1m(ξ) are the frequency parameter and the associated mode shape function.The transient deflection at an arbitrary sampling point “i’’ can be expressed as a cosine series with a supplementary term as (2.19)w1(ξi,t)=δ0i+∑m=1NsT1m(t)δ1mi(i=1,…,N),
where δ1mi=ϕ1m(ξi) needs to be determined. The frequency parameter pm of beams with different boundary conditions can be written in a general form as (2.20)pm=[(m+j)π]2
in which j=0 for a simply supported beam,  j=1/2 for a clamped beam,  j=-1/2 for a beam clamped at one end and free at the other end, and  j=1/4  for a beam clamped at one end but simply supported at the other end.Sincew1(ξi,t) needs to satisfy the initial conditions in (2.12), one has (2.21)δ0i=w10i-∑m=1Nδ1mi.
The time domain function in   (2.19) can be readily obtained as(2.22)w1(ξi,t)=w10i+∑m=1Ns(cospm1t-1)δ1mi(i=1,…,N).
Differentiation of   (2.22) with respect to time twice gives(2.23)ẅ1(ξi,t)=-∑m=1Nspm2(cospm1t)δ1mi(i=1,…,N).The deflection response of the 1st-order perturbation equation can be determined by(2.24)w1(ξ,t)=∑i=1Nw1(ξi,t)Li(ξ).Substitution of   (2.23) and application of DQM approximation (2.15) to   (2.16) result  in a semianalytical algebraic equation at an arbitrary sampling point(2.25)∑j=1NCij(4)∑m=1Ns(cospmtk-1)δ1mj-p12∑m=1Nspm2cos(pmtk)δ1mi=-∑j=1NCij(4)w10j(i=1,…,N;k=1,…,Ns).
This linear equation system can be written in a matrix form as(2.26)[D1]{δ1}={f1},
where [D1] and {δ1} are the coefficient matrix and unknown vector to be solved, and {f1} is the generalized load vector due to initial deflection. Once {δ1} is determined from  (2.26), the deflections at all of the sampling points can be calculated from   (2.22) and the 1st-order deflection solution of the beam is then obtained from   (2.24).


### 2.4. Semianalytical Solution of the 2nd-Order Perturbation Equation

To solve the 2nd-order perturbation equation (2.9), the unknown axial displacement u2 is expressed as(2.27)u2(ξ,t)=∑m=1NsT2m(t)ϕ2m(ξ).
Substituting the 1st-order deflection solution and   (2.27) into   (2.9) yields (2.28)∑m=1NsT2m(t)d2ϕ2mdξ2=-∑m=1NS∑n=1NST1m(t)T1n(t)dϕ1mdξd2ϕ1ndξ2.
This gives (2.29)T2m=-∑n=1NsT1mT1n(dϕ1m/dξ)(d2ϕ1n/dξ2)d2ϕ2m/dξ2.
Therefore,   (2.27) becomes(2.30)u2(ξ,t)=∑m=1Ns∑n=1NsT1mT1nψ2mn(ξ)
with(2.31)ψ2mn(ξ)=-(dϕ1m/dξ)(d2ϕ1n/dξ2)d2ϕ2m/dξ2ϕ2m.
From   (2.30), the axial displacement u2 at an arbitrary sampling point “i’’ is    (2.32)u2(ξi,t)=∑m=1Ns∑n=1NsT1m(t)T1n(t)δ2mni.
Putting   (2.30) into   (2.9) and then applying DQM approximation, one has (2.33)∑j=1NCij(2)∑m=1Ns∑n=1NsT1mT1nδ2mnj=-∑j=1NCij(1)∑m=1NsT1mδ1mj∑k=1NCik(2)∑n=1NsT1nδ1nk(i=1,…,N)
in a matrix form (2.34)[D2]·{δ2}={f2}.
Note that the right-hand-side term {f2} is the pseudodynamic force vector totally dependent on the 1st-order deflection solution already obtained. After the unknown vector {δ2} composed of δ2mnj is determined from   (2.34), u2 can be calculated from the relationship (2.35)ψ2mn(ξ)=∑i=1Nδ2mniLi(ξ)
and   (2.30).


### 2.5. Semianalytical Solution of the 3rd-Order Perturbation Equation

The 3rd-order perturbation equation (2.10) can be rewritten as (2.36)∂4w3∂ξ4+∂2w3∂t2=f(ξ,t),
where the right-hand-side term(2.37)f(ξ,t)={[∂u2∂ξ+12(∂w1∂ξ)2]}∂2w1∂ξ2-σ2∂2w1∂t2.Substituting   (2.17) and (2.30) into   (2.37) gives(2.38)f(ξ,t)=∑l=1Ns∑m=1NS∑n=1NsT1lT1mT1nφlmn+∑n=1Nspn2σ2T1nϕ1n,
where(2.39)φlmn(ξ)={dψ2lmdξ+12[dϕ1ldξdϕ1mdξ]}d2ϕ1ndξ2.
The solution of   (2.36) can be readily obtained by the method of separation of variables (2.40)w3(ξ,t)=∑r=1NsT3r(t)ϕ3r(ξ).
It is evident from   (2.16) and (2.36) that ϕ3m(ξ) is the same as ϕ1m(ξ). Therefore,(2.41)w3(ξ,t)=∑r=1NsT3r(t)ϕ1r(ξ),
and   (2.36) becomes(2.42)∑r=1NsT3r(t)∂4ϕ1r∂ξ4+∑r=1Ns∂2T3r∂t2ϕ1r=f(ξ,r).
From   (2.16)–(2.18), one has (2.43)∂4ϕ1r∂ξ4=Pr2ϕ1r.
Expanding  f(ξ,t) into Fourier series in function ϕ1m,   (2.42) becomes (2.44)∑r=1Ns[T3r′′+T3rpr2]ϕ1r(ξ)=∑r=1Nsfr(ξ,t)ϕ1r(ξ)
which yields (2.45)T3r′′+T3rpr2=fr(ξ,t),
where right-hand-side term (2.46)fr(t)=1N3r2∫01f(ξ,t)ϕ1r(ξ)dξ=∑l=1Ns∑m=1Ns∑n=1NsalmnrT1lT1mT1n+∑n=1Nsσ2bnrpn2T1n,N3r2=∫01ϕ1r2(ξ)dξ,almnr=1N3r2∫01φlmn(ξ)ϕ1r(ξ)dξ,bnr=1N3r2∫01ϕ1n(ξ)ϕ1r(ξ)dξ.
Since ϕ1r are orthogonal functions, we have(2.47)bnr={0,n≠r,1N3r2∫01ϕ1r·ϕ1r(ξ)dξ=1,n=r.
To remove the secular terms in   (2.43), f3m(t) is expanded as(2.48)fr(t)={∑l=1Ns∑m=1Ns∑n=1Nscos(pl+pm+pn)talmnr+∑l=1Ns∑m=1Ns∑n=1Nsl≠n,n≠mcos(pl+pm-pn)talmnr+∑l=1Ns∑m=1Ns∑n=1Nsl≠m,m≠ncos(pl-pm+pn)talmnr+∑l=1Ns∑m=1Ns∑n=1Nsl≠m,l≠ncos(pl-pm-pn)talmnr}+{2∑l=1Ns∑m=1Nsl=ncospmtalmlr+2∑l=1Ns∑m=1Nsm=ncospltalmmr+2∑l=1Ns∑n=1Nsl=mcospntallnr+3∑l=1Nsl=m=ncospltalllr}+pr2σ2rcosprt.
The following condition must be satisfied in order to remove the last two terms in   (2.48),  (2.49){2∑l=1Ns∑m=1Nsl=ncospmtalmlr+2∑l=1Ns∑m=1Nsm=ncospltalmmr+2∑l=1Ns∑n=1Nsl=mcospntallnr+3∑l=1Nsl=m=ncospltalllr}+pr2σ2rcosprt=0.
The sum of the coefficients for cos(prt) must be zero, which gives the expressions of σ2 for different values of m as (2.50)σ2r=-1pr2{2∑m=1Ns[amrmr+armmr+ammrr]+3arrrr}.
Thus the term f3m(t) without secular terms is(2.51)fr(t)={∑l=1Ns∑m=1Ns∑n=1Nscos(pl+pm+pn)talmnr+∑l=1Ns∑m=1Ns∑n=1Nsl≠n,n≠mcos(pl+pm-pn)talmnr+∑l=1Ns∑m=1Ns∑n=1Nsl≠m,m≠ncos(pl-pm+pn)talmnr+∑l=1Ns∑m=1Ns∑n=1Nsl≠m,l≠ncos(pl-pm-pn)t·almnr}.
Equation (2.45) can then be rewritten as(2.52)T3r′′+T3rpr2=fr(t).
From   (2.13), the initial condition for   (2.52) can be derived by using separation of variable method  (2.53)T3r(0)=dT3r(0)dt=0.
The solution of   (2.52) under the given initial condition (2.53) can be obtained through Laplace transform as(2.54)T3r(t)=∑l=1M∑m=1M∑n=1Mcos(prt)-cos((pl+pm+pn)t)(pl+pm+pn)2-pr2γlmnr+∑l=1M∑m=1M∑n=1Ml≠n,n≠mcos(prt)-cos((pl+pm-pn)t)(pl+pm-pn)2-pr2γlmnr+∑l=1M∑m=1M∑n=1Ml≠m,m≠ncos(prt)-cos((pl-pm+pn)t)(pl-pm+pn)2-pr2γlmnr+∑l=1M∑m=1M∑n=1Ml≠m,l≠ncos(prt)-cos((pl-pm-pn)t)(pl-pm-pn)2-pr2γlmnr(whereγlmnr=14almnr).
Finally, we have(2.55)w3(ξ,t)=∑r=1NsT3r(t)ϕ1r(ξ).
The deflections response of nonlinear free vibration of the beam is then calculated from (2.56)w=w1ε+w3ε3+⋯
while the nonlinear frequency response, with higher-order terms being neglected, is given by   (2.50) and (2.7) as(2.57)pr2=1+ε2σ2r,
and the nonlinear frequency response can be obtained from(2.58)(ωrωA)2=1-ε21pr2{2∑m=1Ns[amrmr+armmr+ammrr]+3arrrr}.


## 2.1. Perturbation Equations

Consider an isotropic, homogeneous slender beam of lengthL with uniform cross-section area subjected to a dynamic transverse load  q̅(x,t̅). Let u̅ and w̅ be the displacements parallel to the x- and y-axes, t̅   the time, and m denote  the mass density. Based on classical Euler-Bernoulli beam theory and von-Karman nonlinear displacement-strain relationship, the geometrically nonlinear governing equations for flexural vibration of the beam can be  derived as
(2.1)∂∂x[EA̅(∂u¯∂x+12(∂w¯∂x)2)]=0,EI∂4w̅∂x4-EA̅[∂u̅∂x+12(∂w̅∂x)2]∂2w̅∂x2+m∂2w̅∂t̅2=q̅(x,t̅),
where E is the elastic modulus of the beam, and A̅ and  I are the area and the second moment of the cross section. The present analysis considers general boundary conditions, that is, the beam end may be either simply supported or clamped, with the following boundary conditions: (2.2)u̅=0,w̅=0,M̅x=0(simplysupported),u̅=0,w̅=0,dw̅dx=0(clamped),
where M̅x is the bending moment of the beam. To facilitate the solution process of  (2.1)  by using the differential quadrature method in space domain, the following quantities are introduced: (2.3)u=uη¯,w=w¯ρ,ξ=xL,t=ωt¯,q=L4q̅ρEI,q=εF,p=ωωL,ωL=(1L)2EIm,ρ2=IA̅,
where ωL and ω refer to the linear and nonlinear frequencies of the beam, respectively, and ε is the small perturbation parameter. Substituting these quantities into (2.1)   leads to dimensionless nonlinear governing equations as below (2.4)∂2u∂ξ2+∂w∂ξ∂2w∂ξ2=0,∂4w∂ξ4-[∂u∂ξ+12(∂w∂ξ)2]∂2w∂ξ2+p2∂2w∂t2=εF.Deflectionw, axial displacement u, and frequency parameter p can be expanded into series forms in terms of a small perturbation  parameter ε(2.5)w=w1ε+w3ε3+w5ε5+⋯,(2.6)u=u2ε2+u4ε4+⋯,(2.7)p2=1+ε2σ2+ε4σ4+⋯.
Inserting   (2.5)–(2.7)  into   (2.4), and equating the terms of the same power of ε yield  a series of perturbation equations (2.8)ε1:∂4w1∂ξ4+∂2w1∂t2=F,(2.9)ε2:∂2u2∂ξ2+∂w1∂ξ∂2w1∂ξ2=0,(2.10)ε3:∂4w3∂ξ4-∂u2∂ξ∂2w1∂ξ2-12(∂w1∂ξ)2∂2w1∂ξ2+∂2w3∂t2+σ2∂2w1∂t2=0,(2.11)ε4:∂2u4∂ξ2+∂w1∂ξ∂2w3∂ξ2+∂w3∂ξ∂2w1∂ξ2=0.The initial conditions considered in the present study are(2.12)w1(ξ,0)=w10,∂w1(ξ,0)∂t=0,(2.13)w3(ξ,0)=0,∂w3(ξ,0)∂t=0,
where w10 refers to initial vibration amplitude. By making use of the dimensionless quantities and perturbation series defined above, the dimensionless form of the associated boundary conditions in each perturbation can also be obtained.


## 2.2. Differential Quadrature Method

According to the DQM rule, an unknown functiong(x) and its nth order derivative with respect to x can be approximated as the weighted linear sums of its function values at a number of sampling points in the x-axis as (2.14)g(x)=∑i=1Ng(xi)Li(x),(2.15)dng(x)dxn|x=xi=∑j=1NCij(n)g(xj)(n=1,…,N-1),
where N is the total number of sampling points, Li(x) are Lagrange interpolation polynomials, and the weighting coefficients Cij(n) can be calculated from recursive formulae [20, 21].


## 2.3. Semianalytical Solution of the 1st-Order Perturbation Equation

This study is focused on the nonlinear vibration due to an initial deflection. In such a case, the right-hand-side term in   (2.8)  F=0. Hence, the first-order perturbation   becomes(2.16)∂4w1∂ξ4+∂2w1∂t2=0.
The solution of   (2.16) can be readily obtained by separation of variables(2.17)w1(ξ,t)=∑m=0NsT1m(t)ϕ1m(ξ),
where Ns is the number of terms in the truncated time series,                (2.18)T1m(t)=cos(pmt)
is the solution in the time domain for the first-order perturbation equation, and  pm and δ1m(ξ) are the frequency parameter and the associated mode shape function.The transient deflection at an arbitrary sampling point “i’’ can be expressed as a cosine series with a supplementary term as (2.19)w1(ξi,t)=δ0i+∑m=1NsT1m(t)δ1mi(i=1,…,N),
where δ1mi=ϕ1m(ξi) needs to be determined. The frequency parameter pm of beams with different boundary conditions can be written in a general form as (2.20)pm=[(m+j)π]2
in which j=0 for a simply supported beam,  j=1/2 for a clamped beam,  j=-1/2 for a beam clamped at one end and free at the other end, and  j=1/4  for a beam clamped at one end but simply supported at the other end.Sincew1(ξi,t) needs to satisfy the initial conditions in (2.12), one has (2.21)δ0i=w10i-∑m=1Nδ1mi.
The time domain function in   (2.19) can be readily obtained as(2.22)w1(ξi,t)=w10i+∑m=1Ns(cospm1t-1)δ1mi(i=1,…,N).
Differentiation of   (2.22) with respect to time twice gives(2.23)ẅ1(ξi,t)=-∑m=1Nspm2(cospm1t)δ1mi(i=1,…,N).The deflection response of the 1st-order perturbation equation can be determined by(2.24)w1(ξ,t)=∑i=1Nw1(ξi,t)Li(ξ).Substitution of   (2.23) and application of DQM approximation (2.15) to   (2.16) result  in a semianalytical algebraic equation at an arbitrary sampling point(2.25)∑j=1NCij(4)∑m=1Ns(cospmtk-1)δ1mj-p12∑m=1Nspm2cos(pmtk)δ1mi=-∑j=1NCij(4)w10j(i=1,…,N;k=1,…,Ns).
This linear equation system can be written in a matrix form as(2.26)[D1]{δ1}={f1},
where [D1] and {δ1} are the coefficient matrix and unknown vector to be solved, and {f1} is the generalized load vector due to initial deflection. Once {δ1} is determined from  (2.26), the deflections at all of the sampling points can be calculated from   (2.22) and the 1st-order deflection solution of the beam is then obtained from   (2.24).


## 2.4. Semianalytical Solution of the 2nd-Order Perturbation Equation

To solve the 2nd-order perturbation equation (2.9), the unknown axial displacement u2 is expressed as(2.27)u2(ξ,t)=∑m=1NsT2m(t)ϕ2m(ξ).
Substituting the 1st-order deflection solution and   (2.27) into   (2.9) yields (2.28)∑m=1NsT2m(t)d2ϕ2mdξ2=-∑m=1NS∑n=1NST1m(t)T1n(t)dϕ1mdξd2ϕ1ndξ2.
This gives (2.29)T2m=-∑n=1NsT1mT1n(dϕ1m/dξ)(d2ϕ1n/dξ2)d2ϕ2m/dξ2.
Therefore,   (2.27) becomes(2.30)u2(ξ,t)=∑m=1Ns∑n=1NsT1mT1nψ2mn(ξ)
with(2.31)ψ2mn(ξ)=-(dϕ1m/dξ)(d2ϕ1n/dξ2)d2ϕ2m/dξ2ϕ2m.
From   (2.30), the axial displacement u2 at an arbitrary sampling point “i’’ is    (2.32)u2(ξi,t)=∑m=1Ns∑n=1NsT1m(t)T1n(t)δ2mni.
Putting   (2.30) into   (2.9) and then applying DQM approximation, one has (2.33)∑j=1NCij(2)∑m=1Ns∑n=1NsT1mT1nδ2mnj=-∑j=1NCij(1)∑m=1NsT1mδ1mj∑k=1NCik(2)∑n=1NsT1nδ1nk(i=1,…,N)
in a matrix form (2.34)[D2]·{δ2}={f2}.
Note that the right-hand-side term {f2} is the pseudodynamic force vector totally dependent on the 1st-order deflection solution already obtained. After the unknown vector {δ2} composed of δ2mnj is determined from   (2.34), u2 can be calculated from the relationship (2.35)ψ2mn(ξ)=∑i=1Nδ2mniLi(ξ)
and   (2.30).


## 2.5. Semianalytical Solution of the 3rd-Order Perturbation Equation

The 3rd-order perturbation equation (2.10) can be rewritten as (2.36)∂4w3∂ξ4+∂2w3∂t2=f(ξ,t),
where the right-hand-side term(2.37)f(ξ,t)={[∂u2∂ξ+12(∂w1∂ξ)2]}∂2w1∂ξ2-σ2∂2w1∂t2.Substituting   (2.17) and (2.30) into   (2.37) gives(2.38)f(ξ,t)=∑l=1Ns∑m=1NS∑n=1NsT1lT1mT1nφlmn+∑n=1Nspn2σ2T1nϕ1n,
where(2.39)φlmn(ξ)={dψ2lmdξ+12[dϕ1ldξdϕ1mdξ]}d2ϕ1ndξ2.
The solution of   (2.36) can be readily obtained by the method of separation of variables (2.40)w3(ξ,t)=∑r=1NsT3r(t)ϕ3r(ξ).
It is evident from   (2.16) and (2.36) that ϕ3m(ξ) is the same as ϕ1m(ξ). Therefore,(2.41)w3(ξ,t)=∑r=1NsT3r(t)ϕ1r(ξ),
and   (2.36) becomes(2.42)∑r=1NsT3r(t)∂4ϕ1r∂ξ4+∑r=1Ns∂2T3r∂t2ϕ1r=f(ξ,r).
From   (2.16)–(2.18), one has (2.43)∂4ϕ1r∂ξ4=Pr2ϕ1r.
Expanding  f(ξ,t) into Fourier series in function ϕ1m,   (2.42) becomes (2.44)∑r=1Ns[T3r′′+T3rpr2]ϕ1r(ξ)=∑r=1Nsfr(ξ,t)ϕ1r(ξ)
which yields (2.45)T3r′′+T3rpr2=fr(ξ,t),
where right-hand-side term (2.46)fr(t)=1N3r2∫01f(ξ,t)ϕ1r(ξ)dξ=∑l=1Ns∑m=1Ns∑n=1NsalmnrT1lT1mT1n+∑n=1Nsσ2bnrpn2T1n,N3r2=∫01ϕ1r2(ξ)dξ,almnr=1N3r2∫01φlmn(ξ)ϕ1r(ξ)dξ,bnr=1N3r2∫01ϕ1n(ξ)ϕ1r(ξ)dξ.
Since ϕ1r are orthogonal functions, we have(2.47)bnr={0,n≠r,1N3r2∫01ϕ1r·ϕ1r(ξ)dξ=1,n=r.
To remove the secular terms in   (2.43), f3m(t) is expanded as(2.48)fr(t)={∑l=1Ns∑m=1Ns∑n=1Nscos(pl+pm+pn)talmnr+∑l=1Ns∑m=1Ns∑n=1Nsl≠n,n≠mcos(pl+pm-pn)talmnr+∑l=1Ns∑m=1Ns∑n=1Nsl≠m,m≠ncos(pl-pm+pn)talmnr+∑l=1Ns∑m=1Ns∑n=1Nsl≠m,l≠ncos(pl-pm-pn)talmnr}+{2∑l=1Ns∑m=1Nsl=ncospmtalmlr+2∑l=1Ns∑m=1Nsm=ncospltalmmr+2∑l=1Ns∑n=1Nsl=mcospntallnr+3∑l=1Nsl=m=ncospltalllr}+pr2σ2rcosprt.
The following condition must be satisfied in order to remove the last two terms in   (2.48),  (2.49){2∑l=1Ns∑m=1Nsl=ncospmtalmlr+2∑l=1Ns∑m=1Nsm=ncospltalmmr+2∑l=1Ns∑n=1Nsl=mcospntallnr+3∑l=1Nsl=m=ncospltalllr}+pr2σ2rcosprt=0.
The sum of the coefficients for cos(prt) must be zero, which gives the expressions of σ2 for different values of m as (2.50)σ2r=-1pr2{2∑m=1Ns[amrmr+armmr+ammrr]+3arrrr}.
Thus the term f3m(t) without secular terms is(2.51)fr(t)={∑l=1Ns∑m=1Ns∑n=1Nscos(pl+pm+pn)talmnr+∑l=1Ns∑m=1Ns∑n=1Nsl≠n,n≠mcos(pl+pm-pn)talmnr+∑l=1Ns∑m=1Ns∑n=1Nsl≠m,m≠ncos(pl-pm+pn)talmnr+∑l=1Ns∑m=1Ns∑n=1Nsl≠m,l≠ncos(pl-pm-pn)t·almnr}.
Equation (2.45) can then be rewritten as(2.52)T3r′′+T3rpr2=fr(t).
From   (2.13), the initial condition for   (2.52) can be derived by using separation of variable method  (2.53)T3r(0)=dT3r(0)dt=0.
The solution of   (2.52) under the given initial condition (2.53) can be obtained through Laplace transform as(2.54)T3r(t)=∑l=1M∑m=1M∑n=1Mcos(prt)-cos((pl+pm+pn)t)(pl+pm+pn)2-pr2γlmnr+∑l=1M∑m=1M∑n=1Ml≠n,n≠mcos(prt)-cos((pl+pm-pn)t)(pl+pm-pn)2-pr2γlmnr+∑l=1M∑m=1M∑n=1Ml≠m,m≠ncos(prt)-cos((pl-pm+pn)t)(pl-pm+pn)2-pr2γlmnr+∑l=1M∑m=1M∑n=1Ml≠m,l≠ncos(prt)-cos((pl-pm-pn)t)(pl-pm-pn)2-pr2γlmnr(whereγlmnr=14almnr).
Finally, we have(2.55)w3(ξ,t)=∑r=1NsT3r(t)ϕ1r(ξ).
The deflections response of nonlinear free vibration of the beam is then calculated from (2.56)w=w1ε+w3ε3+⋯
while the nonlinear frequency response, with higher-order terms being neglected, is given by   (2.50) and (2.7) as(2.57)pr2=1+ε2σ2r,
and the nonlinear frequency response can be obtained from(2.58)(ωrωA)2=1-ε21pr2{2∑m=1Ns[amrmr+armmr+ammrr]+3arrrr}.


## 3. Numerical Results

To illustrate its efficiency and accuracy, the proposed semianalytical approach is used to study the nonlinear frequency and dynamic response of Euler-Bernoulli beams of rectangular cross section(b×h=0.01m×0.02m) under various boundary conditions with a given initial vibration amplitude and zero-valued initial velocity. The dimensionless vibration amplitude is defined as A=W/ρ (or A=6W/3h for a rectangular beam) where W is the vibration amplitude at the midpoint of the beam. The parameters used herein are L=0.5m, E=150 GPa, and  m=8×103kg/m3. The results are presented in both tabular and graphical forms in terms of the normalized frequency ratio ωn/ωL and dimensionless dynamic deflections at the midpoint of the beam.Among the sampling point distribution schemes available, a nonuniform distribution is chosen in the present analysis due to its excellent convergence [21, 26], that is,(3.1)xk=12(1-coskπN-1)(k=0,1,…,N-1).In what follows, the total number of sampling points in the space domain isN = 11, unless stated otherwise.In Figures1–3, the present solutions and exact solutions are represented by solid lines and circular dots, respectively.Figure 1
Nonlinear fundamental frequency versus dimensionless central amplitude relationship of beams with different boundary conditions.Figure 2
Nondimensional central deflection response of a simply supported beam under different initial vibration amplitudes.Figure 3
Nondimensional central deflection response of a simply supported beam.


### 3.1. Nonlinear Vibration Frequency

Tables1, 2, and 3 list, respectively, the normalized fundamental frequency ratios ωn/ωL for a simply supported beam, a clamped beam, and a beam simply supported at one end but clamped at the other end at varying dimensionless vibration amplitudes A=0.1~2.0.  All beams exhibit typical “hard-spring” behavior, that is, the frequency ratio increases with an increase in vibration amplitude. Comparisons between the results with varying total number of sampling points N=9,11,13 show that the present method converges very well to produce very accurate results when N≥11. Our results are compared with the existing results obtained by using elliptic function [1] and the finite element method [28]. The relative errors, defined as(3.2)(PresentSolution-ExistingSolution)ExistingSolution×100%,
are also provided. Excellent agreement can be observed in these tables. The relationship between nonlinear fundamental frequency and dimensionless vibration amplitude (A=0.0~6.0) of these three beams are shown in Figure 1 where curves a, b, and c are for a simply supported beam, a beam simply supported at one end but clamped at the other end, and a clamped beam. The scatter points represent the elliptical function results using the formulae provided by Woinowsky-Krieger [1]. As can be seen, the present results and the elliptical function results are almost identical. The clamped beam has the highest fundamental frequency with the lowest deflection while the simply supported beam has the smallest fundamental frequency and the largest deflection. This is because the end support of the clamped beam is much stronger than its simply supported counterpart.Table 1
Normalized frequency ratioωn/ωL for a simply supported beam.
AN=9N=11N=13Woinowsky-Krieger [1]Mei [28]0.11.00091.00091.00091.00081.00090.21.00371.00371.00371.00381.00370.41.01491.01491.01491.01501.01480.61.03311.03311.03341.03801.03390.81.05801.05811.05881.05801.05781.01.08911.08931.09081.08901.08891.51.18931.19101.19321.19001.19022.01.31401.32011.33001.31601.3022Table 2
Normalized frequency ratioωn/ωL for a clamped beam.
AN=9N=11N=13Mei [28]0.11.00031.00031.00031.00030.21.00121.00121.00121.00120.41.00501.00471.00481.00480.61.01121.01061.01071.01070.81.01981.01871.01901.01901.01.03071.02901.02951.02961.51.06761.06401.06501.06532.01.11691.11001.11251.1135Table 3
Normalized frequency ratioωn/ωL for a beam simply supported at one end and clamped at the other end.
AN=9N=11N=13Mei [28]0.11.00061.00061.00061.00060.21.00261.00261.00261.00260.41.01031.01031.01031.01060.61.02301.02311.02311.02380.81.04041.04061.04071.04181.01.06241.06271.06281.06471.51.13461.13531.13551.14042.01.22721.22831.22861.2385


### 3.2. Nonlinear Dynamic Response

Convergence study on the nonlinear dynamic response is conducted in both Tables4 and 5 where central deflections of a simply supported beam at different time (in sec), calculated by using different total number of sampling points N in the space domain and different number of truncated series terms Ns in the time domain, are given and compared with the results calculated by using the formulae given by Woinowsky-Krieger [1]. The dimensionless vibration amplitude is A=1.0. Error 1, Error 2, and Error 3 in Table 4 indicate the relative error percentages between the exact solutions and the present results with N=9 and Ns=1, N=9 and Ns=3, N=9 and Ns=5 while Error 4, Error 5, and Error 6 in Table 5 are the relative error percentages between the exact solutions and the present results with N=9 and Ns=3, N=11 and Ns=3, N=13 and Ns=3. Our semianalytical solutions converge monotonically as Ns increases but nonmonotonically with an increasing N. However, the present method is capable of giving very accurate results with relative error less than 1.5% when Ns≥3 and N≥11. We have calculated the nonlinear dynamic response of beams with other end supports as well and found that the results exhibit the same convergence characteristics as observed in Tables 4 and 5. These results are, therefore, not detailed here.Table 4
Nonlinear central deflection response of a simply supported beam: solutions with varying number of terms in time domain seriesNs.
Time (s)Ns=1Ns=3Ns=5Woinowsky-Krieger [1]Error 1Error 2Error 30.00500.86220.85660.85530.85331.05050.38950.24050.01000.48110.47200.47000.46653.13601.18280.75400.0150-0.0420-0.0398-0.0395-0.04043.7991-1.6178-2.20360.0300-0.9965-0.9939-0.9942-0.99630.0134-0.2437-0.21460.0400-0.4060-0.3952-0.3938-0.39283.35880.62300.27170.04500.12590.12400.12310.12104.06922.47801.72170.05500.93350.93000.92970.92860.52940.15430.11770.06500.76720.75960.75740.75451.68040.66770.38040.0750-0.2086-0.2022-0.2010-0.20093.82590.66910.06130.0850-0.9596-0.9547-0.9549-0.95660.3065-0.2001-0.17670.0900-0.9698-0.9655-0.9657-0.96760.2261-0.2173-0.19600.0950-0.7113-0.6999-0.6983-0.69741.98720.34960.12450.1000-0.2477-0.2404-0.2390-0.23873.76150.70060.1317Table 5
Nonlinear central deflection response of a simply supported beam: solutions with varying total number of sampling pointsN.
Time (s)N=9N=11N=13Woinowsky-Krieger [1]Error 4Error 5Error 60.00500.85660.85620.85620.85330.38950.34580.34840.01000.47200.47110.47120.46651.18280.99221.00260.0150-0.0398-0.0410-0.0410-0.0404-1.61781.46241.42280.0300-0.9939-0.9967-0.9967-0.9963-0.24370.04110.03780.0400-0.3952-0.3969-0.3970-0.39280.62301.06061.06570.04500.12400.12230.12230.12102.47801.06511.08390.05000.61060.60970.60970.60480.95130.79840.80540.06500.75960.75850.75860.75450.66770.53000.53350.0750-0.2022-0.2033-0.2033-0.20090.66911.19241.19030.0850-0.9547-0.9580-0.9579-0.9566-0.20010.13910.13770.0900-0.9655-0.9687-0.9687-0.9676-0.21730.11400.11220.0950-0.6999-0.7017-0.7017-0.69740.34960.60710.60930.1000-0.2404-0.2415-0.2415-0.23870.70061.17161.1717Figure2 depicts the nonlinear dynamic central deflection response of a simply supported beam with dimensionless initial deflection  w10=0.5sin(πξ), 0.8sin(πξ),1.0sin(πξ),  and 1.4sin(πξ). The time interval in this example is relatively short, that is, t=0~0.1 sec. The result within a much longer time interval t=0~1 sec is displayed in Figure 3  where  w10=1.0sin(πξ). N=11 and Ns=3 are used in these two examples. Again, the present results agree very well with the elliptical function solutions by Woinowsky-Krieger [1]. It is also seen that a bigger initial deflection magnitude results in larger vibration amplitude.It is worthy of noting that unlike the conventional numerical integration schemes commonly used to calculate dynamic response whose accuracy is quickly degraded at later time steps due to the accumulation of numerical errors, the proposed technique does not involve iterative time integration and accumulative errors and is therefore able to give results with excellent accuracy even for a very long time interval, as demonstrated in Figure3.


## 3.1. Nonlinear Vibration Frequency

Tables1, 2, and 3 list, respectively, the normalized fundamental frequency ratios ωn/ωL for a simply supported beam, a clamped beam, and a beam simply supported at one end but clamped at the other end at varying dimensionless vibration amplitudes A=0.1~2.0.  All beams exhibit typical “hard-spring” behavior, that is, the frequency ratio increases with an increase in vibration amplitude. Comparisons between the results with varying total number of sampling points N=9,11,13 show that the present method converges very well to produce very accurate results when N≥11. Our results are compared with the existing results obtained by using elliptic function [1] and the finite element method [28]. The relative errors, defined as(3.2)(PresentSolution-ExistingSolution)ExistingSolution×100%,
are also provided. Excellent agreement can be observed in these tables. The relationship between nonlinear fundamental frequency and dimensionless vibration amplitude (A=0.0~6.0) of these three beams are shown in Figure 1 where curves a, b, and c are for a simply supported beam, a beam simply supported at one end but clamped at the other end, and a clamped beam. The scatter points represent the elliptical function results using the formulae provided by Woinowsky-Krieger [1]. As can be seen, the present results and the elliptical function results are almost identical. The clamped beam has the highest fundamental frequency with the lowest deflection while the simply supported beam has the smallest fundamental frequency and the largest deflection. This is because the end support of the clamped beam is much stronger than its simply supported counterpart.Table 1
Normalized frequency ratioωn/ωL for a simply supported beam.
AN=9N=11N=13Woinowsky-Krieger [1]Mei [28]0.11.00091.00091.00091.00081.00090.21.00371.00371.00371.00381.00370.41.01491.01491.01491.01501.01480.61.03311.03311.03341.03801.03390.81.05801.05811.05881.05801.05781.01.08911.08931.09081.08901.08891.51.18931.19101.19321.19001.19022.01.31401.32011.33001.31601.3022Table 2
Normalized frequency ratioωn/ωL for a clamped beam.
AN=9N=11N=13Mei [28]0.11.00031.00031.00031.00030.21.00121.00121.00121.00120.41.00501.00471.00481.00480.61.01121.01061.01071.01070.81.01981.01871.01901.01901.01.03071.02901.02951.02961.51.06761.06401.06501.06532.01.11691.11001.11251.1135Table 3
Normalized frequency ratioωn/ωL for a beam simply supported at one end and clamped at the other end.
AN=9N=11N=13Mei [28]0.11.00061.00061.00061.00060.21.00261.00261.00261.00260.41.01031.01031.01031.01060.61.02301.02311.02311.02380.81.04041.04061.04071.04181.01.06241.06271.06281.06471.51.13461.13531.13551.14042.01.22721.22831.22861.2385


## 3.2. Nonlinear Dynamic Response

Convergence study on the nonlinear dynamic response is conducted in both Tables4 and 5 where central deflections of a simply supported beam at different time (in sec), calculated by using different total number of sampling points N in the space domain and different number of truncated series terms Ns in the time domain, are given and compared with the results calculated by using the formulae given by Woinowsky-Krieger [1]. The dimensionless vibration amplitude is A=1.0. Error 1, Error 2, and Error 3 in Table 4 indicate the relative error percentages between the exact solutions and the present results with N=9 and Ns=1, N=9 and Ns=3, N=9 and Ns=5 while Error 4, Error 5, and Error 6 in Table 5 are the relative error percentages between the exact solutions and the present results with N=9 and Ns=3, N=11 and Ns=3, N=13 and Ns=3. Our semianalytical solutions converge monotonically as Ns increases but nonmonotonically with an increasing N. However, the present method is capable of giving very accurate results with relative error less than 1.5% when Ns≥3 and N≥11. We have calculated the nonlinear dynamic response of beams with other end supports as well and found that the results exhibit the same convergence characteristics as observed in Tables 4 and 5. These results are, therefore, not detailed here.Table 4
Nonlinear central deflection response of a simply supported beam: solutions with varying number of terms in time domain seriesNs.
Time (s)Ns=1Ns=3Ns=5Woinowsky-Krieger [1]Error 1Error 2Error 30.00500.86220.85660.85530.85331.05050.38950.24050.01000.48110.47200.47000.46653.13601.18280.75400.0150-0.0420-0.0398-0.0395-0.04043.7991-1.6178-2.20360.0300-0.9965-0.9939-0.9942-0.99630.0134-0.2437-0.21460.0400-0.4060-0.3952-0.3938-0.39283.35880.62300.27170.04500.12590.12400.12310.12104.06922.47801.72170.05500.93350.93000.92970.92860.52940.15430.11770.06500.76720.75960.75740.75451.68040.66770.38040.0750-0.2086-0.2022-0.2010-0.20093.82590.66910.06130.0850-0.9596-0.9547-0.9549-0.95660.3065-0.2001-0.17670.0900-0.9698-0.9655-0.9657-0.96760.2261-0.2173-0.19600.0950-0.7113-0.6999-0.6983-0.69741.98720.34960.12450.1000-0.2477-0.2404-0.2390-0.23873.76150.70060.1317Table 5
Nonlinear central deflection response of a simply supported beam: solutions with varying total number of sampling pointsN.
Time (s)N=9N=11N=13Woinowsky-Krieger [1]Error 4Error 5Error 60.00500.85660.85620.85620.85330.38950.34580.34840.01000.47200.47110.47120.46651.18280.99221.00260.0150-0.0398-0.0410-0.0410-0.0404-1.61781.46241.42280.0300-0.9939-0.9967-0.9967-0.9963-0.24370.04110.03780.0400-0.3952-0.3969-0.3970-0.39280.62301.06061.06570.04500.12400.12230.12230.12102.47801.06511.08390.05000.61060.60970.60970.60480.95130.79840.80540.06500.75960.75850.75860.75450.66770.53000.53350.0750-0.2022-0.2033-0.2033-0.20090.66911.19241.19030.0850-0.9547-0.9580-0.9579-0.9566-0.20010.13910.13770.0900-0.9655-0.9687-0.9687-0.9676-0.21730.11400.11220.0950-0.6999-0.7017-0.7017-0.69740.34960.60710.60930.1000-0.2404-0.2415-0.2415-0.23870.70061.17161.1717Figure2 depicts the nonlinear dynamic central deflection response of a simply supported beam with dimensionless initial deflection  w10=0.5sin(πξ), 0.8sin(πξ),1.0sin(πξ),  and 1.4sin(πξ). The time interval in this example is relatively short, that is, t=0~0.1 sec. The result within a much longer time interval t=0~1 sec is displayed in Figure 3  where  w10=1.0sin(πξ). N=11 and Ns=3 are used in these two examples. Again, the present results agree very well with the elliptical function solutions by Woinowsky-Krieger [1]. It is also seen that a bigger initial deflection magnitude results in larger vibration amplitude.It is worthy of noting that unlike the conventional numerical integration schemes commonly used to calculate dynamic response whose accuracy is quickly degraded at later time steps due to the accumulation of numerical errors, the proposed technique does not involve iterative time integration and accumulative errors and is therefore able to give results with excellent accuracy even for a very long time interval, as demonstrated in Figure3.


## 4. Conclusions

A new semianalytical method for geometrically nonlinear vibration analysis of Euler-Bernoulli beams with different boundary conditions is presented in this paper. The proposed method is based on the perturbation technique, differential quadrature approximation, and an analytical series expansion in the time domain. Numerical results show that this method has excellent accuracy and convergence characteristics. Compared with other numerical approaches that use step-by-step time integration, the present method has a unique advantage of being capable of producing results with very good and stable accuracy at a long time interval because the error accumulation is avoided due to the use of the analytical time series.



---

*Source: 591786-2010-07-21.xml*