Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

54, 4, pp. 1219-1230, Warsaw 2016
DOI: 10.15632/jtam-pl.54.4.1219

HOMOTOPY ANALYSIS OF A FORCED NONLINEAR BEAM MODEL

WITH QUADRATIC AND CUBIC NONLINEARITIES

Shahram Shahlaei-Far, Airton Nabarrete, José Manoel Balthazar

Aeronautics Institute of Technology, ITA, São José dos Campos, Brazil

e-mail: shazzz 85@yahoo.de

This study investigates forced nonlinear vibrations of a simply supported Euler-Bernoulli
beam on a nonlinear elastic foundation with quadratic and cubic nonlinearities. Applying
the homotopy analysis method (HAM) to the spatially discretized governing equation, we
derivenovel analytical solutions anddiscuss their convergence topresentnonlinear frequency
responses with varying contributions of the nonlinearity coefficients. A comparisonwith nu-
merical solutions is conducted and nonlinear time responses and phase planes are compared
to the results from linear beam theory. The study demonstrates that apart from nonlinear
problems of free vibrations, HAM is equally capable of solving strongly nonlinear problems
of forced vibrations.

Keywords: forced nonlinear vibration, HAM, quadratic and cubic nonlinearities, Galerkin
method

1. Introduction

Nonlinear vibrations of a simply supported Euler-Bernoulli beam on a nonlinear elastic founda-
tion with distributed quadratic and cubic nonlinearities subject to harmonic excitation are, in
nondimensional form, governed by (Abe, 2006)

∂2w

∂t2
+
∂4w

∂x4
+2µ
∂w

∂t
+α2w

2+α3w
3 =F(x)cos(Ωt) (1.1)

subject to the boundary conditions

w=
∂2w

∂x2
=0 at x=0,1 (1.2)

wherew is the displacement of the beam,µ is the viscous damping coefficient,α2 andα3 are the
quadratic and cubic nonlinearity coefficients, respectively, and F(x) andΩ are the distribution
and frequency of the harmonic load, respectively.
Abe (2006) providednumerical solutions for the cases of primary resonance and subharmonic

resonance of the order one-half by means of the finite difference and shooting method, the first
applied directly toEqs. (1.1) and (1.2) and the latter to their spatially discretized formusing the
Galerkinmethod.Heconcludedthat theGalerkindiscretizationyieldsmoreaccurate results than
the direct approach. Considering only small amplitude vibrations, he also obtained perturbative
approximate solutions using the method of multiple scales (MMS) (Nayfeh and Mook, 1979)
applied directly to Eq. (1.1)modeled as aweakly nonlinear system involving small perturbation
parameters. However, Abe et al. (1998,a,b,c, 2000) showed that for the direct approach it is
unlikely to obtain highly accurate solutions as it is not possible to define a detuning parameter
in the quadratic form between the natural and the excitation frequency.
The present work provides accurate analytical solutions to Eqs. (1.1) and (1.2) using the

homotopy analysis method (HAM) introduced by Liao (Liao, 1992, 1995, 2003, 2004, 2009,



1220 S. Shahlaei-Far et al.

2012; Liao and Cheung, 1998; Liao and Tan, 2007) to investigate the case when the excitation
frequency is close to the natural frequency of the fundamental mode.Making no use of small or
large parameters, thismethod allows us to consider the general nonlinear distributed-parameter
system for small aswell as large amplitude vibrations.Thus it overcomes the requirement for the
MMS tomodel the problemas aweakly nonlinear system involving only small finite amplitudes.
To ensureaccurate results,we reduce thegoverning equation to anordinarynonlineardifferential
equation with the Galerkin method before applying HAM. For the convergence analysis of the
obtained analytical solutions, we plot the so-called h-curves for higher-order approximations
and achieve the optimal value of h byminimizing the square residual of the governing equation
for a chosen order of the approximation. We verify our higher-order solutions by comparison
with numerical solutions given by the fourth-order Runge-Kutta method. Demonstrating how
powerful a first-order approximation of HAMalready is, we derive explicit closed-form solutions
of the mean of motion, amplitude and phase of the generalized coordinate to present frequency
response curves for various values of the quadratic and cubic nonlinearity coefficients addressing
the cases of softening- and hardening-type behavior. Moreover, we compare the nonlinear time
responses and their respective phase planes for different values of the quadratic nonlinearities
to the results from linear beam theory.
HAMhas seen extensive successful applications to highly nonlinear problems in science and

engineering (Wen andCao, 2007; Hoseini et al., 2008; Pirbodaghi et al., 2009; Abbasbandy and
Shivanian, 2011; Mastroberardino, 2011; Mehrizi et al., 2012; Mohammadpour et al., 2012; Mu-
stafa et al., 2012; Sedighi et al., 2012; Ray andSahoo, 2015). Previous applications ofHAMhave
beenmainly devoted to nonlinear differential equations of free vibratory continuous systems.As
seen inEq. (1.1), the important case of forced nonlinear vibrations of a damped system involving
cubic and quadratic nonlinearities, where the mean of motion cannot be disregarded, poses a
greater challenge with respect to finding analytical solutions and understanding nonlinear beha-
vior.Our results provide an example, applicable to other beammodels andboundary conditions,
of a forced nonlinear vibratory system in which HAM yields accurate convergent solutions for
all values of the relevant parameters.

2. Governing differential equation

The model is discretized by the Galerkin method with w(x,t) = W(t)φ(x), where
φ(x) =

√
2sin(πx) is the normalized eigenfunction of the fundamental mode and W(t), the

corresponding time-dependent amplitude, is the generalized coordinate. Thus, governing equ-
ation (1.1) is reduced to the nonlinear ordinary differential equation

Ẅ +2µẆ +ω2W +
8
√
2

3π
α2W

2+
3

2
α3W

3 = f cos(Ωt−ϕ) (2.1)

and without loss of generality subject to the initial conditions

W(0)= δ+A Ẇ(0)= 0 (2.2)

whereω=π2 is the normalized natural frequency of the fundamentalmode, f =
∫1
0 F(x)φ(x)dx

is the first modal force, A is an unknown amplitude and δ = (1/T)
∫T
0 W(t)dt is the mean

of motion being generally nonzero for oscillations with the quadratic nonlinearity. The dot
represents differentiation with respect to time t. Note that, for convenience and without loss of
generality, we introduce the phase angleϕ in the expression of the harmonic load as a quantity
to be determined.
Defining the variables

τ =Ωt W(t)= δ+AV (τ) (2.3)



Homotopy analysis of a forced nonlinear beam model... 1221

and inserting them into Eqs. (2.1) and (2.2), we obtain

Ω2A
∂2V (τ)

∂τ2
+Ω2µA

∂V (τ)

∂τ
+ω2(δ+AV (τ))+

8
√
2

3π
α2(δ+AV (τ))

2

+
3

2
α3(δ+AV (τ))

3 = f1cosτ+f2 sinτ

(2.4)

subject to the initial conditions

V (0)= 1
∂V (0)

∂τ
=0 (2.5)

with the constants f1 and f2 satisfying

f21 +f
2
2 = f

2 ϕ=arctan
f2
f1

(2.6)

3. Homotopy analysis method

The homotopy analysis method is a nonperturbative analytical technique for solving nonlinear
differential equations. Bymeans of an embedding parameter ranging from zero to one, it trans-
forms a nonlinear differential equation into an infinite number of linear differential equations
and derives a family of solution series.
The periodic solution to Eq. (2.4) can be expressed by a set of base functions

{sin(mτ),cos(mτ) |m=1,2,3, . . .} (3.1)

such that

V (τ)=
infty
∑

k=1

(αk sin(kτ)+βk cos(kτ)) (3.2)

where αk and βk are coefficients to be determined.We choose the initial guess as

V0 =cosτ (3.3)

which satisfies the initial conditions of Eq. (2.5). To ensure the rule of solution expression given
by Eq. (3.2), we choose the linear operator to be

L[φ(τ,q)] =Ω2
(∂2φ(τ,q)

∂τ2
+φ(τ,q)

)

(3.4)

with the property

L[C1 sinτ+C2cosτ] = 0 (3.5)

where C1 and C2 are constants of integration. According to Eq. (2.4), we define the nonlinear
operator

N [φ(τ,q),Λ(q),∆(q)] =Ω2Λ(q)∂
2φ(τ,q)

∂τ2
+2µΩΛ(q)

∂φ(τ,q)

∂τ

+ω2(∆(q)+Λ(q)φ(τ,q))+
8
√
2

3π
α2(∆(q)+Λ(q)φ(τ,q))

2

+
3

2
α3(∆(q)+Λ(q)φ(τ,q))

3 −f1cosτ−f2 sinτ

(3.6)



1222 S. Shahlaei-Far et al.

where q∈ [0,1] is the embedding parameter, φ(τ,q) is a function of τ and q,Λ(q) and∆(q) are
functions of q. The zeroth-order deformation equation is given by

(1−q)L[φ(τ,q)−V0(τ)] = qhH(τ)N [φ(τ,q),Λ(q),∆(q)] (3.7)

where h 6= 0 is the convergence-control parameter and H(τ) a nonzero auxiliary function. For
q=0 and q=1we have

φ(τ,0)=V0(τ) φ(τ,1)=V (τ) Λ(0)=A0

Λ(1)=A ∆(0)= δ0 ∆(1)= δ
(3.8)

Thus, the function φ(τ,q) varies from the initial guess V0(τ) to the desired solution as q varies
from 0 to 1. The Taylor expansions of φ(τ,q), Λ(q) and∆(q) with respect to q are

φ(τ,q) =V0(τ)+
infty
∑

m=1

Vm(τ)q
m Λ(q)=A0+

infty
∑

m=1

Amq
m

∆(q)= δ0+
infty
∑

m=1

δmq
m

(3.9)

where

Vm(τ)=
1

m!

∂mφ(τ,q)

∂qm

∣

∣

∣

q=0
Am =

1

m!

∂mΛ(q)

∂qm

∣

∣

∣

q=0
δm =

1

m!

∂m∆(q)

∂qm

∣

∣

∣

q=0

(3.10)

Choosing properly the auxiliary function H(τ) and the convergence-control parameter h, the
series in Eqs. (3.9) converge when q=1, such that

V (τ)=V0(τ)+
infty
∑

m=1

Vm(τ) A=A0+
infty
∑

m=1

Am δ= δ0+
∞
∑

m=1

δm (3.11)

Differentiating zeroth-order equation (3.7) m times with respect to q, dividing it by m! and
setting q=0, them-th order deformation equation is obtained as

L[Vm(τ)−χmVm−1(τ)] =hH(τ)Rm(Vm−1,Am−1,δm−1) (3.12)

subject to the initial conditions

Vm(0)=
∂Vm(0)

∂τ
=0 (3.13)

where

χm =

{

0 m¬ 1
1 m> 1

Vm−1 = [V0(τ),V1(τ), . . . ,Vm−1(τ)]

Am−1 = [A0,A1, . . . ,Am−1] δm−1 = [δ0,δ1, . . . ,δm−1]

(3.14)



Homotopy analysis of a forced nonlinear beam model... 1223

and

Rm(Vm−1,Am−1,δm−1)=
1

(m−1)!
dm−1N [φ(τ,q),Λ(q),∆(q)]

dqm−1

∣

∣

∣

q=0

=Ω2Am−1
∂2Vm−1
∂τ2

+2µΩAm−1
∂Vm−1
∂τ

+ω2(δm−1+Am−1Vm−1)

+
8
√
2

3π
α2

m−1
∑

n=0

(δn+AnVn)(δm−1−n+Am−1−nVm−1−n)

+
3

2
α3

m−1
∑

n=0

(

n
∑

j=0

(δj +AjVj)(δn−j +An−jVn−j)

)

· (δm−1−n+Am−1−nVm−1−n)− (1−χm)(f1cosτ+f2 sinτ)

=Cm,0+

l(m)
∑

k=1

(

cm,k(Vm−1,Am−1,δm−1)cos(kτ)+dm,k(Vm−1,Am−1,δm−1)sin(kτ)
)

(3.15)

For the nonzero auxiliary function to obey the rule of solution expression and the rule of coeffi-
cient ergodicity, we choose it to be

H(τ)= cos(2τ) (3.16)

where κ is an integer. It can be determined uniquely asH(τ)= 1.
According to the property of the linear operator, if the terms sinτ and cosτ exist in Rm,

the secular terms τ cosτ and τ sinτ will appear in the final solution, therefore cm,1 and dm,1
have to be equal to zero.Moreover, ifCm,0 6=0, a constant termwill appear in the final solution
violating the rule of solution expression, thus it must be set to zero.
The general solution to Eq. (3.12) for m­ 1 is derived to be

Vm(τ)=χmVm−1(τ)+
h

Ω2

l(m)
∑

k=2

cm,k cos(kτ)+dm,k sin(kτ)

1−k2 +C1 sinτ+C2cosτ (3.17)

whereC1 andC2 need to be determined by the initial conditions in Eq. (3.13).
For the first-order approximation (m=1) we obtain from Eq. (3.15)

R1(V0,A0,δ0)=−Ω2A0cosτ−2µΩA0 sinτ +ω2(δ0+A0cosτ)

+
8
√
2

3π
α2(δ0+A0cosτ)

2+
3

2
α3(δ0+A0cosτ)

3−f1cosτ

−f2 sinτ =
8
√
2

6π
α2A

2
0+

(

ω2+
9

4
α3
)

δ0+
8
√
2

3π
α2δ
2
0 +
3

2
α3δ
3
0

+
(

−Ω2A0+ω2A0+
16
√
2

3π
α2δ0A0+

9

2
α3δ
2
0A0+

9

8
α3A

3
0−f1

)

cosτ

+(−2µΩA0−f2)sinτ+
(8
√
2

6π
α2A

2
0+
9

8
α3δ0A

2
0

)

cos(2τ)+
3

8
α3A

3
0cos(3τ)

=C1,0+

l(1)=3
∑

k=1

(

c1,k(V0,A0,δ0)cos(kτ)+d1,k(V0,A0,δ0)sin(kτ)
)

(3.18)

with d1,2 = d1,3 = 0. Enforcing the constant term and the coefficients of sinτ and cosτ to be
equal to zero, and using the conditions inEq. (2.6), we obtain steady state solutions of themean
of motion δ0, the amplitudeA0 and the phaseϕ, respectively, as

δ0 =−
16
√
2α2

27α3π
− K1
27α3

3
√
2π

3

√

K2+
√

K22 +4K
3
1

+
1

54α3
3
√
2π

3

√

K2+
√

K22 +4K
3
1 (3.19)



1224 S. Shahlaei-Far et al.

with

K1 =−2048α22 +162α3π(9A20α3π+4π5) K2 =−131072
√
2α32+62208

√
2α2α3π

6

and

(

(ω2−Ω2)A0+
16
√
2

3π
α2δ0A0+

9

2
α3δ
2
0A0+

9

8
α3A

3
0

)2
+(−2µΩA0)2 = f2

ϕ=arctan
2µΩ

Ω2−ω2− 16
√
2

3π
α2δ0− 92α3δ20 −

9
8
α3B

B=
ω2δ0+

8
√
2
3π
α2δ
2
0 +

3
2
α3δ
3
0

9
4
α3δ0+

8
√
2
6π
α2

(3.20)

Finally, solving the first-order deformation equation of Eq. (3.12), the general solution is

V1(τ)=
h

Ω2

l(1)=3
∑

k=2

c1,k
1−k2 cos(kτ)+C1 sinτ+C2cosτ

=
h

Ω2

(

(4
√
2

9π
α2+

3

4
α3δ0

)

A20(cosτ− cos(2τ))+
3

64
α3A

3
0(cosτ − cos(3τ))

)

(3.21)

where the constants C1 andC2 are obtained from the initial conditions in Eq. (3.13).

Thus, with Eqs. (3.11), the first-order approximation ofW(t) becomes

W(t)= δ+AV (τ)≈ δ0+A0(V0(τ)+V1(τ))

= δ0+
h

Ω2

(

(Ω2

h
A0+

4
√
2

9π
α2A

3
0+
3

4
α3δ0A

3
0+
3

64
α3A

4
0

)

cos(Ωt)

+
(4
√
2

9π
α2A

3
0+
3

4
α3δ0A

3
0

)

cos(2Ωt)+
3

64
α3A

4
0cos(3Ωt)

)

(3.22)

4. Convergence of HAM solution

Applying HAM to a nonlinear problem results in a family of solution series which depend
on the convergence-control parameter h. In order to ensure the convergence and rate of the
approximation for the HAM solution, valid convergence regions for the auxiliary parameter h
need to be obtained. By means of plotting h-curves this can be achieved and thus a convergent
solution series is guaranteed. Since a valid region comprises a range of possible values of h, the
optimal choice is obtained by minimizing the square residual of the governing equation for a
given order of the approximation. For this, we consider the Nth-order approximations of Eqs.
(3.11) given by

VN(τ)=V0(τ)+
N
∑

m=1

Vm(τ) A=A0+
N
∑

m=1

Am δ= δ0+
N
∑

m=1

δm (4.1)

Inserting Eqs. (4.1)-(4.3) into Eq. (3.6) with q = 1, we can define the square residual error for
theN-th order approximation as

eN(h)=

1
∫

0

(N [VN(τ),AN ,δN])2 dτ (4.2)



Homotopy analysis of a forced nonlinear beam model... 1225

The solution series is convergent when eN(h)→ 0 asN →∞. The optimal value ofh for a given
orderN of the approximation is obtained by the solution of the algebraic equation

deN
dh
=0 (4.3)

It is to be noted that the calculation for each order is done separately and not iteratively.

Fig. 1. h-curves of the amplitudeA for α2 =5: (a) α3 =0.5, (b) α3 =1, (c) α3 =1.5

Figure 1 shows the effect of the auxiliary parameter on the solution convergence for higher-
-order approximations. The valid region is characterized by the flat portion which is common
to all h-curves displayed. The corresponding optimal values of h for varying values of the cubic
nonlinearity coefficient α3 are presented in Table 1 using the software Mathematica.

5. Discussion of results

Asimply supportedEuler-Bernoulli beam resting on a nonlinear elastic foundationwith quadra-
tic and cubic nonlinearities is considered. First, higher-order approximations from the general
solution in Eq. (3.17), obtained with the software Mathematica, are compared to numerical re-
sults. Secondly, considering thefirst-order approximation ofHAM, the frequency response curves
of the amplitude obtained in Eq. (3.20)1 are presented for different values of the quadratic and
cubic nonlinearity coefficients. Moreover, the nonlinear time response curves and phase planes
are compared to the results from linear beam theory showing the effect of these nonlinearities
upon the distributed-parameter system.
In Fig. 2, the accuracy of nonlinear time responses obtained by a sixth-order HAM ap-

proximation for α2 = 5 is validated by comparison with the numerical results achieved by the
fourth-orderRunge-Kuttamethod for varying values ofα3 and the corresponding optimal values
of the auxiliary parameter h. There is accurate agreement between the analytical and numerical
results.



1226 S. Shahlaei-Far et al.

Table 1. Optimal values of h and minimum values of eN for µ = 0.05, F(x) =
√
2sin(πx),

α2 =5

α3 N Optimal value of h Minimum value of eN

0.5 2 −0.3567 8.374 ·10−4
4 −0.3324 1.705 ·10−7
6 −0.3286 4.119 ·10−9
8 −0.3045 7.805 ·10−11

1 2 −0.2943 5.096 ·10−2
4 −0.2755 1.643 ·10−5
6 −0.2619 4.295 ·10−6
8 −0.2578 9.349 ·10−8

1.5 2 −0.2581 6.912 ·10−1
4 −0.2319 3.782 ·10−2
6 −0.2247 2.186 ·10−4
8 −0.2169 6.237 ·10−6

Fig. 2. Comparison of sixth-order HAM solutions (solid line) with numerical results by the fourth-order
Runge-Kutta method (symbols)

Fixing α2 =5 andα3 =0.5 and using the optimal values of h fromTable 1, approximations
of order 2, 4, 6 and 8 are compared to numerical results (fourth-order Runge-Kutta) in Table 2.
The results show that low-order approximations by HAM agree well with numerical solutions,
although by increasing the order of HAM iterations, the accuracy increases.

For the rest of this Section, we consider the first-order approximation byHAM. InFig. 3, we
investigate the influence of the nonlinearity coefficients on the frequency response curves, that
is, we show the variation of the response curves with α2 for different values of α3. In Fig. 3a,
as the excitation frequency Ω nears the fundamental frequency ω, the response of the system
exhibits hardening-spring nonlinear characteristics due to the cubic nonlinearity. Increasing the
value ofα3, this hardening-type behavior is further increased. In the presence of bothα2 andα3,
in Fig. 3b, the response is almost linear (for α3 = 0.5) suggesting that the magnitude of the
nonlinearities cancel each other out. Further increase of α3 results, as in the previous case, in
hardening-type behavior. Figure 3c predicts softening- as well as hardening-type responses for
increasing values of α3 demonstrating the softening effect of the quadratic nonlinearity. This
transition is also evident in Fig. 3d, whereby the softening-spring effect is more pronounced.
For higher values of the quadratic nonlinearity, the softening effect on the system becomesmore
distinctive.



Homotopy analysis of a forced nonlinear beam model... 1227

Table 2. Comparison of higher-order HAM solutions with numerical results for µ = 0.05,
F(x)=

√
2sin(πx), α2 =5, α3 =0.5

Time HAM solution Numerical
t 2-nd order 4-th order 6-th order 8-th order results

0 2.70874 2.71538 2.71796 2.71923 2.71925

0.1 1.34795 1.35269 1.35344 1.35416 1.35418

0.2 −1.50326 −1.51083 −1.51278 −1.51490 −1.51492
0.3 −3.21418 −3.21964 −3.22461 −3.22791 −3.22793
0.4 −2.15817 −2.16227 −2.16435 −2.16520 −2.16521
0.5 0.65104 0.66129 0.66483 0.66659 0.66662

0.6 2.61598 2.62004 2.62135 2.62337 2.62339

0.7 1.92956 1.93118 1.93207 1.93479 1.93480

Fig. 3. Amplitude-frequency curves for µ=0.05,F(x)=
√
2sin(πx), α3 =0.5 (solid line),

α3 =1 (dashed line), α3 =1.5 (dotted line): (a)α2 =5, (b) α2 =7.4, (c) α2 =10.5, (d) α2 =12

Investigating the impact of the quadratic nonlinearity on the distributed-parameter system
we compare, in Fig. 4, the nonlinear time responses obtained by HAM with those from linear
beam theory maintaining a fixed value for α3 and varying values of the quadratic nonlinearity
coefficientα2. Theoptimal value ofh is obtainedby solvingEq. (4.5) forN =1. It is evident that
for lower values of α2 there is more agreement between the linear and nonlinear time response,
whereas by increasing the value of α2, the difference becomes considerable.

In Fig. 5, the phase planes for both nonlinear and linear responses are presented for different
values ofα2 while fixingα3 =0.5.With an increase inα2, the phase planes significantly diverge
from their linear counterparts emphasizing the effect of the quadratic nonlinearity term on



1228 S. Shahlaei-Far et al.

Fig. 4. Linear time response (dashed line) versus nonlinear time response (solid line) for µ=0.05,
F(x)=

√
2sin(πx), α3 =0.5, h=−0.3719: (a) α2 =5 (Ω/ω=0.9806), (b) α2 =7.4 (Ω/ω=0.9910),

(c) α2 =10.5 (Ω/ω=1.0075), (d) α2 =12 (Ω/ω=1.0160)

Fig. 5. Phase plane of linear response (thin line) versus nonlinear response (thick line) for µ=0.05,
F(x)=

√
2sin(πx), α3 =0.5, h=−0.3719: (a) α2 =5 (Ω/ω=0.9806), (b) α2 =7.4 (Ω/ω=0.9910),

(c) α2 =10.5 (Ω/ω=1.0075), (d) α2 =12 (Ω/ω=1.0160)



Homotopy analysis of a forced nonlinear beam model... 1229

the system when the excitation frequency is close to the natural frequency of the fundamental
mode.

It should be noted that the analysis presented in this study can be expanded to predict
beam responses for different boundary conditions. To this end, the mode shape φ(x) satisfying
the boundary conditions at both ends should be inserted intow(x,t) =W(t)φ(x) for the spatial
discretization by the Galerkin method. The solutions can then be derived analogously with the
methodology described in Section 3.

6. Conclusion

The present study provides analytical solutions to forced nonlinear vibrations of a simply sup-
ported Euler-Bernoulli beam resting on a nonlinear elastic foundation with quadratic and cubic
nonlinearities using the homotopy analysis method. The convergence of the solution has been
investigated by optimizing the value of the auxiliary convergence-control parameterh. Using the
optimal values of h, higher-order solutions by HAM have been compared to numerical results
demonstrating the effectiveness of the method for low-order approximations and varying values
of the cubic nonlinearity. The derived closed-form solution of the amplitude yields frequency
response curves for various values of the quadratic and cubic nonlinearity coefficients presenting
their softening-/hardening-type effect on the distributed-parameter system. Phase planes and
nonlinear time response curves illustrate the considerable difference with respect to the results
from linear beam theory for various values of the quadratic nonlinearity coefficient. The findings
reveal that HAM is a general solution method that can successfully address highly nonlinear
problems of forced vibrations.

Acknowledgement

The authors acknowledge the financial support received during this work from CAPES and the

MCT/CNPq/FAPEMIG-National Institute of Science and Technology on Smart Structures in Engine-

ering, grant number 574001/2008-5.

References

1. Abbasbandy S., Shivanian E., 2011, Predictor homotopy analysis method and its application
to some nonlinear problems,Communications in Nonlinear Science andNumerical Simulation, 16,
2456-2468

2. Abe A., 2006, On non-linear vibration analyses of continuous systems with quadratic and cubic
non-linearities, International Journal of Non-Linear Mechanics, 41, 873-879

3. Abe A., Kobayashi Y., Yamada G., 1998a, Analysis of subharmonic resonance of moderately
thick antisymmetric angle-ply laminated plates by using method of multiple scales, Journal of
Sound and Vibration, 217, 467-484

4. Abe A., Kobayashi Y., Yamada G., 1998b, Internal resonance of rectangular laminated plates
with degenerate modes, JSME International Journal, Series C, 41, 718-726

5. AbeA.,KobayashiY.,YamadaG., 1998c,Two-mode response of simply supported, rectangular
laminated plates, International Journal of Non-Linear Mechanics, 33, 675-690

6. Abe A., Kobayashi Y., Yamada G., 2000, Non-linear vibration characteristics of clamped la-
minated shallow shells, Journal of Sound and Vibration, 234, 405-426

7. Hoseini S.H.,Pirpodaghi T., AsghariM., FarrahiG.H.,AhmadianM.T., 2008,Nonlinear
free vibration of conservative oscillators with inertia and static type cubic nonlinearities using
homotopy analysis method, Journal of Sound and Vibration, 316, 263-273



1230 S. Shahlaei-Far et al.

8. Liao S.J., 1992, Proposed homotopy analysis techniques for the solution of nonlinear problems,
Ph.D. Dissertation, Shanghai Jiao Tong University, Shanghai

9. LiaoS.J., 1995,Anapproximate solution techniquewhichdoesnotdependupon small parameters:
a special example, International Journal of Non-Linear Mechanics, 30, 371-380

10. Liao S.J., 2003, Beyond Perturbation: Introduction to Homotopy Analysis Method, Chapman &
Hall/CRCPress, Boca Raton

11. Liao S.J., 2004, On the homotopy analysis method for nonlinear problems,Applied Mathematics
and Computation, 147, 499-513

12. Liao S.J., 2009, Notes on the homotopy analysis method: some definitions and theorems, Com-
munications in Nonlinear Science and Simulation, 14, 983-997

13. Liao S.J., 2012,Homotopy Analysis Method in Nonlinear Differential Equations, Springer Educa-
tion Press, Heidelberg

14. LiaoS.J.,CheungA.T., 1998,Applicationofhomotopyanalysismethod innonlinearoscillations,
Journal of Applied Mechanics, 65, 914-922

15. Liao S.J., Tan Y., 2007, A general approach to obtain series solutions of nonlinear differential
equations, Studies in Applied Mathematics, 119, 297-355

16. Mastroberardino A., 2011, Homotopy analysis method applied to electrohydrodynamic flow,
Communications in Nonlinear Science and Numerical Simulation, 16, 2730-2736

17. Mehrizi A.A., Vazifeshenas Y., Domairry G., 2012, New analysis of natural convection bo-
undary layer flow on a horizontal plate with variable wall temperature, Journal of Theoretical and
Applied Mechanics, 50, 4, 1001-1010

18. MohammadpourA., Rokni E., FooladiM., KimiaeifarA., 2012,Bernoulli-Euler beams un-
der different boundary conditions with non-linearWinkler type foundation, Journal of Theoretical
and Applied Mechanics, 50, 2, 339-355

19. Mustafa M., HayatT., Hendi A.A., 2012, Influence ofmelting heat transfer in the stagnation-
point flow of a Jeffrey fluid in the presence of viscous dissipation, Journal of Applied Mechanics,
79, 2, 4501-4505

20. Nayfeh A.H., Mook D.T., 1979,Nonlinear Oscillations, Wiley, NewYork

21. Pirbodaghi T., Ahmadian M.T., Fesanghary M., 2009, On the homotopy analysis method
for non-linear vibration of beams,Mechanics Research Communications, 36, 143-148

22. Ray S.S., Sahoo S., 2015, Traveling wave solutions to Riesz time-fractional Camassa-Holm equ-
ation inmodeling for shallow-waterwaves,Journal of Computational andNonlinear Dynamics,10,
6, 1026-1030

23. Sedighi H.M., Shirazi K.H., Zare J., 2012, An analytic solution of transversal oscillation of
quintic non-linear beam with homotopy analysis method, International Journal of Non-Linear
Mechanics, 47, 777-784

24. Wen J., Cao Z., 2007, Sub-harmonic resonances of nonlinear oscillationswith parametric excita-
tion bymeans of the homotopy analysis method,Physics Letters A, 371, 427-431

Manuscript received August 15, 2015; accepted for print February 18, 2016