Homotopy perturbation study of nonlinear vibration of Von Karman rectangular plates

Computers and Structures 106–107 (2012) 46–55

Contents lists available at SciVerse ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Homotopy perturbation study of nonlinear vibration of Von Karman rectangular plates M.M. Rashidi a,b, A. Shooshtari b, O. Anwar Bég c,⇑ a

Mechanical Engineering, Engineering Department, Université de Sherbrooke, 2500 boul. de l’Université Sherbrooke, (Québec), Canada J1K 2R1 Mechanical Engineering, Engineering Department, Engineering Faculty, Bu-Ali Sina university, Hamedan, Iran c Biomechanics and Astronautics, Aerospace Engineering, Department of Engineering and Mathematics, Sheaf Building, Sheaf Street, Sheffield Hallam University, Sheffield S11WB, United Kingdom b

a r t i c l e

i n f o

Article history: Received 13 November 2011 Accepted 9 April 2012 Available online 1 May 2012 Keywords: Nonlinear vibration Plates Homotopy perturbation method Time function

a b s t r a c t Based on the Von Karman theory, the equations of motion for a rectangular isotropic plate, considering the effect of shear deformation and rotary inertia, have been derived. For the nonlinear vibration of the plate, a nonlinear coupled equation is obtained with an Airy stress function. Using the Galerkin method, this equation is separated into position and time functions. The Homotopy Perturbation Method (HPM) is employed to solve the nonlinear time function. It is shown that the obtained results demonstrate excellent agreement with numerical solutions obtained using the fourth-order Runge–Kutta method. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction Plates are one of the most popular continuous structural systems used in many engineering applications including aircraft structures, nuclear vessels, bridges, power hydraulics etc. Such systems are often exposed to dynamic forces. In this regard the nonlinear vibration of plates is an important area of applied mechanics and has stimulated considerable interest among researchers. There are several methods such as numerical, analytical and approximation methods which can be employed to study the nonlinear vibration of plates. Sathyamoorphy [1] has provided a lucid review on this subject. Bikri et al. [2] investigated the free vibration of thin laminated isotropic plates with a numerical method. Ribeiro and Petyt [3] used the hierarchical finite-element method to investigate the free vibration of an isotropic plate and obtained internal resonances. Analytical methods are however very attractive since they yield closed-form solutions for the dynamical behavior of plates and these solutions can show directly the effects of various system parameters and may be implemented in design applications. One of the most powerful analytical methods in nonlinear problems is the perturbation method [4]. Shooshtari and Khadem [5] used the multiple scale method for study nonlinear vibrations of a rectangular plate. However in traditional perturbation methods of this type, such as the Lindested–Poincare method, multiple time scales method and the generalized averaging methods, the assumption of weak nonlinearity and a small parameter is essential to obtain a solution. ⇑ Corresponding author. E-mail addresses: [email protected], [email protected] (O. Anwar Bég). 0045-7949/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compstruc.2012.04.004

In recent years, a powerful and robust analytical technique, namely the Homotopy Perturbation Method (HPM) [6–8] has been developed and implemented for solving strongly nonlinear problems. Many problems in, for example, fluid dynamics and heat transfer have been solved with this method [9,10]. This method and its variants, has also been employed successfully in some nonlinear vibration problems. Hoseini et al. [11] used the Homotopy Analysis Method (HAM) to analyze nonlinear free vibration of a conservative oscillator with inertia and cubic nonlinear stiffness, presenting an analytical solution valid for a wide range of system parameters. Zand and Ahmadian [12], used this method to derive an accurate expression for dynamic Pull-in instability in micromechanical systems. Pirbodaghi et al. used the HAM approach for studying the non-linear vibration of Euler–Bernoulli beams subjected to axial loads [13], they obtained some analytical relations for nonlinear frequency and buckling load which demonstrated good agreement with the numerical results. In the present article, HPM is used for studying the nonlinear vibration of a rectangular plate. It is shown that the obtained analytical solutions yield very good agreement with the numerical results. Also the effects of various plate parameters on the amplitude of oscillation are investigated in detail. 2. Definition of the problem Considering an isotropic, elastic, rectangular plate with dimensions a and b and thickness h as shown in Fig. 1. In the analysis, the effect of shear deformations and rotary inertia are considered. The plate experiences large deflection and the displacements of the

M.M. Rashidi et al. / Computers and Structures 106–107 (2012) 46–55

47

Fig. 1. Rectangular isotropic plate and its displacement after deflection.

plate are unknown. The displacement relations based on the first order deformation and Von Karman theory, are as follows

uðx; y; z; tÞ ¼ u0 ðx; y; tÞ þ zaðx; y; tÞ;

ð1Þ

v ðx; y; z; tÞ ¼ v 0 ðx; y; tÞ þ zbðx; y; tÞ;

ð2Þ

wðx; y; z; tÞ ¼ w0 ðx; y; tÞ;

ð3Þ

which u, v, w are the displacement in directions of x, y, z respectively, u0, v0, w0 are the displacements of the mid-plane. Also a and b are the angles between the normal to mid-plane before and after deformation. The components of stresses are as follows:

Nx Mx Nx 12Mx þ ¼ þ 3 ; h I h h N M N 12M ry ¼ x þ y ¼ y þ 3 y ; h I h h Nxy M xy Nxy 12Mxy rxy ¼ ; þ ¼ þ 3 h I h! h 2 3 Qx Z sxz ¼ 14 2 ; 2 h h ! 3 Qy Z2 syz ¼ 14 2 ; 2 h h

rx ¼

ð4Þ ð5Þ ð6Þ ð7Þ ð8Þ

where Nx, Ny, Nxy are the membrane in-plane forces, Mx, My, Mxy are the bending moments and Qx and Qy are the internal shear forces. These parameters are shown in Fig. 2.

Considering the first order shear deformation theory and variation calculus, the force–displacement relations for an elastic isotropic plate will become as

 2 @u0 1 @w0 Nx Ny  þ þm ¼ 0; @x 2 @x Eh Eh  2 @ v 0 1 @w0 Ny Nx  þ þm ¼ 0; @y 2 @y Eh Eh @u0 @ v 0 @w0 @w0 ð1 þ mÞ Nxy ¼ 0; þ þ 2 Eh @y @x @x @y @ a 12Mx 12mMy þ ¼ 0;  3 3 @x Eh Eh @b 12My 12mM x  þ ¼ 0; 3 3 @y Eh Eh @ a @b 24ð1 þ mÞ M xy ¼ 0;  þ 3 @y @x Eh @w 12ð1 þ mÞ aþ  Q x ¼ 0; @x 5Eh @w 12ð1 þ mÞ bþ  Q y ¼ 0: @y 5Eh

ð9Þ ð10Þ ð11Þ ð12Þ ð13Þ ð14Þ ð15Þ ð16Þ

Computing the kinetic and potential energy and by the Hamilton principle for conservative systems, the following equations are obtained: 3

@M x @M xy qh att ; þ  Qx ¼ @x @y 12

ð17Þ

3

@M y @M xy qh þ  Qy ¼ b ; @y @x 12 tt @Nx @Nxy þ ¼ qhutt ; @x @y @Ny @Nxy þ ¼ qhv tt ; @y @x

ð18Þ ð19Þ ð20Þ

@Nx @ 2 w @Ny @ 2 w @Nxy @Nx y wx wx þ Nx 2 þ wy þ Ny 2 þ wy þ @x @y @y @x @y @x þ 2Nxy Fig. 2. Internal forces and bending moment in the plate.

@ 2 w @Q x @Q y þ þ ¼ qhwtt : @x@y @x @y

ð21Þ

48

M.M. Rashidi et al. / Computers and Structures 106–107 (2012) 46–55

The general boundary conditions which are obtained from Hamilton principle can be described as: at x = 0 and x = a:

u ¼ 0 or Nx ¼ 0;

ð22Þ

v ¼ 0 or Nxy ¼ 0;

ð23Þ

w ¼ 0 or Q x þ N x

@w @w þ Nxy ¼ 0; @x @y

ð24Þ

a ¼ 0 or Mx ¼ 0;

ð25Þ

b ¼ 0 or Mxy ¼ 0;

ð26Þ

and at y = 0 and y = b:

u ¼ 0 or Nxy ¼ 0;

ð27Þ

v ¼ 0 or Ny ¼ 0;

ð28Þ

w ¼ 0 or Q y þ N xy

@w @w þ Ny ¼ 0; @x @y

ð30Þ

b ¼ 0 or My ¼ 0;

ð31Þ

Eqs. (17)–(21) are derived for the nonlinear vibration of a rectangular isotropic elastic plate under large amplitude incorporating the shear deformation and rotary inertia phenomena. Assuming principally transverse motion occurs, which is based on the results of previous investigations that the in-plane natural frequencies are very far away from transverse natural frequencies, an Airy stress function u may be introduced as follows: 2

@ u Nx ¼ 2 ; @y

ð32Þ

@2u Ny ¼ 2 ; @x @2u Nxy ¼  : @x@y

ð33Þ ð34Þ

 @w qh3  att ¼ 0; @x 12   @ 2 b Dð1 mÞ @ 2 b Dð1þ mÞ @ 2 a Eh @w qh3 D 2þ  þ 5 bþ b ¼ 0; @y 2 @x2 2 12ð1þ mÞ @y @x@y 12 tt   @2 u @2 w @2 u @2 w @2u 5Eh @ a @b þ 2 2 þ þ þ r2 w  qhwtt ¼ 0: 2 2 2 @x @y @y @x @x@y 12ð1þ v Þ dx dy @ 2 a Dð1 mÞ @ 2 a Dð1þ mÞ @ 2 b 5Eh  þ þ 2 2 @x@y 12ð1þ mÞ @x2 @y2



aþ

@2w r u ¼ Eh4 @x@y

!2

3 @ 2 w @ 2 w5 :  2 @x @y2

where u is the Airy stress function and is defined in Eqs. (32)–(34) which must satisfy the Eq. (38). Comparing the Eqs. (39) and (40), it is evident that the presence of shear deformation and rotary inertia will complicate the equation of motion and will generate more nonlinear terms. Also the general boundary condition of the plate (Eqs. (22)–(33)) can be presented based on the definition of Airy stress function as: u0 ¼

Z

þa=2

"

a=2

v0 ¼

Z

þb=2

!  2 # 1 @2 u @2 u 1 @w0 @2u  dx ¼ 0; or N x ¼ 2 ¼ 0;  v 2 2 Eh @y 2 @x @x @y

"

b=2

ð41Þ

!  2 # 1 @2u @2 u 1 @w0 @2 u dx ¼ 0; or N xy ¼ v 2 þ 2  ¼ 0; ð42Þ Eh 2 @y @y @x @x@y

w0 ¼ 0; or KGh

a ¼ 0; or Mx ¼



 @w0 @ 2 u @w0 @ 2 u @w0 þa þ 2  ¼ 0; @x @y @x @x@y @y

ð43Þ

@a @b þ v ¼ 0; @y @x

ð44Þ

@ a @b þ ¼ 0: @y @x

ð45Þ

b ¼ 0; or M xy ¼

At y = b/2 and y = b/2: u0 ¼

Z

þa=2

"

v0 ¼

Z

þb=2

!  2 # 1 @2 u @2 u 1 @w0 @2 u  dx ¼ 0; or N xy ¼  v ¼ 0; Eh @y2 2 @x @x2 @x@y

"

b=2

!  2 # 1 @2u @2 u 1 @w0 @2 u dx ¼ 0; or N y ¼ 2 ¼ 0; v 2 þ 2  Eh 2 @y @y @x @x

w0 ¼ 0; or KGh



 @w0 @ 2 u @w0 @ 2 u @w0 þb   2 ¼ 0; @y @x@y @x @x @y

@ a @b þ ¼ 0; @y @x

ð46Þ

ð47Þ

ð48Þ

a ¼ 0; or Mxy ¼

ð36Þ

b ¼ 0; or M y ¼ v

ð37Þ

The procedure to obtain the Eqs. (41)–(50) has been presented in appendix.

ð38Þ

With eliminating a and b from equations, the following equation is obtained:

" 3 D2 ð1  v Þ 4 Dqh ð3  v Þ 2 @ 2 5DEhð3  v Þ 2 r  r 2 r 2 24 24ð1 þ v Þ @t # 2 4 25E2 h 5qEh @2 þ þ 144ð1 þ v Þ2 72ð1 þ v Þ @t 2 " # @2u @2w @2u @2w @2u @2w 5Eh 2 € þ þ 2 r w  qhw  @x2 @y2 @y2 @x2 @x@y @x@y 12ð1 þ v Þ " # 2 4 5Eh 5DEh ð1  v Þ 4 25E2 h 5qEh @2 2 r  r  þ 12ð1 þ v Þ 24 ð1 þ v Þ 144ð1 þ v Þ @t 2 144ð1 þ v Þ2  w ¼ 0:

ð40Þ

ð35Þ

The proper compatibility equation must be considered for the middle surface strains, which may be stated as follows [27]:

2

# @2 @2u @2 @2u @2 @2u @2 w ¼ 0; Dr þ qh 2  2  2 @x @x2 @y2 @y2 @x@y @x@y @t 4

a=2

Substituting the above functions in the equations of motion, the Eqs. (19) and (20) have been satisfied automatically and the other equations take the form:

4

"

ð29Þ

a ¼ 0 or Mxy ¼ 0;

D

Eq. (39) is the equation of motion in the transverse direction for a Eh3 rectangular isotropic elastic plate, where D ¼ 12ð1 m2 Þ : Eq. (40) is derived for the classical plate theory (CPT) without considering the effect of shear deformation and rotary inertia as:

ð39Þ

@ a @b þ ¼ 0: @x @y

ð49Þ ð50Þ

3. Analytical solution method In order to calculate the transverse mode shapes of the system, Eq. (39) should be solved. The boundary conditions of problem have been selected as movable simply supported boundary conditions in all four edge of the plate. So the selected boundary condition equations for this boundary condition are: At x = a/2 and x = a/2: u0 ¼

Z

þa=2

a=2

v0 ¼

Z

þb=2

b=2

"

!  2 # 1 @2 u @2 u 1 @w0 @2u  dx ¼ 0; or N x ¼ 2 ¼ 0;  v Eh @y2 2 @x @x2 @y

"

!  2 # 1 @2u @2 u 1 @w0 @2 u dx ¼ 0; or N xy ¼ v 2 þ 2  ¼ 0; ð52Þ Eh 2 @y @y @x @x@y

w0 ¼ 0; M xy ¼

ð51Þ

@ a @b þ ¼ 0; @y @x

b ¼ 0;

and at y = b/2 and y = b/2:

ð53Þ ð54Þ ð55Þ

49

M.M. Rashidi et al. / Computers and Structures 106–107 (2012) 46–55

"

!

2 #

 1 @2 u @2 u 1 @w0 @2 u dx ¼ 0; or N xy ¼ v 2  ¼ 0; 2 Eh @y 2 @x @x @x@y a=2 " ! #  2 Z þb=2 1 @2 u @2 u 1 @w0 @2 u dx ¼ 0; or N y ¼ 2 ¼ 0; v 2 þ 2  v0 ¼ Eh 2 @y @y @x @x b=2 u0 ¼

Z

þa=2

w0 ¼ 0; a ¼ 0; @ a @b M y ¼ v þ ¼ 0: @x @y

0.0004

ð56Þ ð57Þ

0.0003

ð58Þ ð59Þ

0.0002

f

ð60Þ 0.0001

For solution of Eq. (39), there are many different methods which can be employed including the numerical methods (e.g. finite-element method) and analytical methods (e.g. perturbation methods). The method adopted in the present study is the HPM. To solve the above nonlinear equations, first by the Galerkin method, one discretizes the equation by writing the solution as: M X N X wðx; y; tÞ ¼ wij ðx; yÞf ðtÞ:

0

-0.0001

ð61Þ

i¼1 j¼1

-0.0002

Next, by selecting wij(x, y) which satisfies the boundary conditions, we use the HPM to find the time function f(t). For this purpose without any lack of generality, the case of a square plate with all sides hinged and movable boundary condition will be considered in the first mode. Therefore we have:

w11 ðx; yÞ ¼ h sin

Numerical HPM

px a

sin

py ; a

0

5

10

w11 ðx; y; tÞ ¼ hf ðtÞ sin

a

sin

a

1E-07

p4 f 2 2

2a2 b



cos

    2p x 2p y þ cos : a b

5E-08

ð64Þ

Considering movable boundary condition (Eqs. (51)–(60)) and solving the above equation, the function u will be obtained as:

 3 Eh 2p2 2px 2p y 2 2 : þ cos u ¼ f ðtÞ þ y Þ þ cos ðx a a 32 a2 ð1  v 2 Þ 2

Error

3

Eh

x x ¼ ; a

y y ¼ ; a

ð65Þ -1E-07

ð67Þ

where

p8 r6 ð3  mÞ 5p6 ð3  mÞ2 r4 25p4 r 2 ð3  mÞ þ þ ; 2304 ð1  m2 Þ2 4608ð1  m2 Þ2 4608ð1 þ m2 Þð1  mÞ 5p6 r4 25p4 r 2 a2 ¼ þ ; 2 1728ð1 þ mÞ ð1  m2 Þ 1728ð1 þ mÞ2 ð1  m2 Þ p4 r4 5p2 r 2 ð23  11mÞ 25 ; a3 ¼ þ þ 2 2Þ 2 3456ð1 þ m Þð1  m 576ð1 þ mÞ2 1728ð1 þ mÞ ð1  m Þ p6 r6 ð3  mÞ2 5p4 r 4 ð3  mÞ a4 ¼ þ ; 2 1536ð1  mÞð1  m Þ 768ð1  m2 Þ

a1 ¼

a5 ¼

p6 r6 ð3  mÞ2 5p4 r 4 ð3  mÞ þ ; 768ð1  mÞð1  m2 Þ 384ð1  m2 Þ

-1.5E-07

0

5

ð68Þ ð69Þ ð70Þ ð71Þ

10

15

20

25

30

t

ð66Þ

and substituting Eqs. (63) and (65) in Eq. (39) and multiplying the



left-hand side by sin pax sin pay and integrating over the plate area leads to a non-dimensional equation of the form:

a1 f 3 þ a2 f þ a3 f 00 þ a4 f 2 f 00 þ a5 f ðf 0 Þ2 ¼ 0;

0

-5E-08

Defining the non-dimensional parameters (i.e. time and distances) as:

sffiffiffiffiffiffiffiffi E s¼t ; qa2

Numerical-HPM

ð63Þ

Substituting the above equation in the compatibility Eq. (38) that equation became as:

r4 u ¼ 

30

2E-07

1.5E-07

:

25

Fig. 3. Comparison of amplitude oscillation for the HPM and numerical solution (f(0) = 0.0002, f0 (0) = 0, m = 0.3, r = 0.1).

ð62Þ

py

20

t

so

px

15

Fig. 4. Comparison of amplitude oscillation error between the HPM and numerical solution (f(0) = 0.0002, f0 (0) = 0, m = 0.3, r = 0.1).

and r = h/a. Comparing Eq. (67) with the equation which is derived by Nayfeh [14,15], we observe that the first term of these equation represents the nonlinearity in the stiffness and the last two terms are the nonlinear terms in inertia. So, this nonlinear system possesses both nonlinearity in stiffness and inertia. Now Eq. (67) is rewritten as:

f 00 þ b1 f þ b2 f 3 þ b3 f 2 f 00 þ b4 f ðf 0 Þ2 ¼ 0;

ð74Þ

where

b1 ¼ a2 =a3 ¼ Oð1Þ;

ð75Þ

b2 ¼ a1 =a3 ¼ Oðr2 Þ;

ð76Þ

b3 ¼ a4 =a3 ¼ Oðr4 Þ;

ð77Þ

4

b4 ¼ a5 =a3 ¼ Oðr Þ;

ð78Þ

with definition of the parameter r = h/a as a small parameter

ð73Þ

pffiffiffi

e ¼ r ¼ h=a:

ð79Þ

50

M.M. Rashidi et al. / Computers and Structures 106–107 (2012) 46–55

0.004

0.04 Numerical HPM

0.003

Numerical HPM

0.03

f

0.02

f

0.002

0.001

0.01

0

0

-0.001

-0.01

-0.002

0

5

10

15

20

25

-0.02

30

0

5

10

t

25

30

Fig. 7. Comparison of amplitude oscillation for the HPM and numerical solution (f(0) = 0.02, f0 (0) = 0, m = 0.3, r = 0.1).

4E-07

2E-07 Numerical-HPM

1.5E-07

1E-07

2E-07

5E-08

1E-07

0

0

-5E-08

-1E-07

-1E-07

-2E-07

0

5

10

15

20

Numerical-HPM

3E-07

Error

Error

20

t

Fig. 5. Comparison of amplitude oscillation for the HPM and numerical solution (f(0) = 0.002, f0 (0) = 0, m = 0.3, r = 0.1).

-1.5E-07

15

25

-3E-07

30

0

5

t

10

15

20

25

30

t

Fig. 6. Comparison of amplitude oscillation error between the HPM and numerical solution (f(0) = 0.002, f0 (0) = 0, m = 0.3, r = 0.1).

Fig. 8. Comparison of amplitude oscillation error between the HPM and numerical solution (f(0) = 0.02, f0 (0) = 0, m = 0.3, r = 0.1).

Eq. (74) will be written as

f 00 þ x2 f þ Gðf ; f 0 ; f 00 Þ ¼ 0;

ð80Þ

and 





Gðf ; f 0 ; f 00 Þ ¼ b2 f 3 þ b3 f 2 f 00 þ b4 f ðf 0 Þ2 :

ð81Þ

where A is a general differential operator, B is a boundary operator, f(r) is a known analytic function and C is the boundary of the domain X. The operator A can, generally speaking, be divided into two parts L and N, where L is linear, and N is non-linear. Therefore Eq. (82) can be written as follows:

Now, Eq. (80) can be solved by the HPM.

LðuÞ þ NðuÞ  f ðrÞ ¼ 0:

4. Explicit analytical solution with the HPM

By the homotopy technique, one can construct a homotopy tðr; pÞ : X  ½0; 1 ! R which satisfies

Now, for convenience, consider the following general non-linear differential equation

Hðt; pÞ ¼ ð1  pÞ½LðtÞ  Lðu0 Þ þ p½AðtÞ  f ðrÞ ¼ 0;

AðuÞ  f ðrÞ ¼ 0;

r 2 X;

ð82Þ

with boundary condition:

  @u ¼ 0; B u; @n

r 2 C;

ð83Þ

ð84Þ

ð85Þ

or

Hðt; pÞ ¼ LðtÞ  Lðu0 Þ þ pLðu0 Þ þ p½NðtÞ  f ðrÞ ¼ 0;

ð86Þ

where p 2 [0, 1] is an embedding parameter and u0 is the initial approximation of Eq. (84) which satisfies the boundary conditions. Clearly, we have

51

M.M. Rashidi et al. / Computers and Structures 106–107 (2012) 46–55

4E-05 Numerical HPM

3E-05

f

2E-05

1E-05

0

-1E-05

-2E-05

0

5

10

15

20

25

30

t

Fig. 11. The second mode shape of oscillation for different time (f(0) = 0.0002, f0 (0) = 0, m = 0.3, r = 0.1, a = 10, h = 1, t = 2).

Fig. 9. Comparison of amplitude oscillation for the HPM and numerical solution (f(0) = 0.00002, f0 (0) = 0, m = 0.3, r = 0.04).

1E-07 Numerical-HPM

5E-08

Error

0

-5E-08

-1E-07

Fig. 12. The second mode shape of oscillation for different time (f(0) = 0.0002, f0 (0) = 0, m = 0.3, r = 0.1, a = 10, h = 1, t = 10).

0

5

10

15

20

25

30

t Fig. 10. Comparison of amplitude oscillation error between the HPM and numerical solution (f(0) = 0.00002, f0 (0) = 0, m = 0.3, r = 0.04).

Hðt; 0Þ ¼ LðtÞ  Lðu0 Þ ¼ 0;

ð87Þ

Hðt; 1Þ ¼ AðtÞ  f ðrÞ ¼ 0:

ð88Þ

The changing process of p from zero to unity is just that of t(r, p) changing from u0(r) to u(r); this is referred to as deformation, and also, L(t)  L(u0) and A(t)  f(r) are called homotopic in topology. Furthermore if the embedding parameter p(0 6 p 6 1) is considered as a ‘‘small parameter’’, applying the classical perturbation technique, we can naturally assume that the solution of Eqs. (87) and (88) can be given as a power series in p, i.e.,

t ¼ t0 þ pt1 þ p2 t2 þ    ;

ð89Þ

and setting p = 1 results in the approximate solution of Eq. (62) as

u ¼ limt ¼ t0 þ t1 þ t2 þ    : p!1

ð90Þ

The convergence of series (90) has been proved by He [16]. It is pertinent to mention that the major advantage of He’s homotopy

perturbation method is that the perturbation equation can be freely constructed in many ways (therefore is problem dependent) by homotopy in topology and the initial approximation can also be freely selected. Moreover, the construction of the homotopy for the perturbation problem plays a very important role in obtaining the desired accuracy [17]. According to Eq. (85), a homotopy can be constructed as follows

ð1  pÞ½f 00 þ b1 f  þ p½f 00 þ b1 f þ b2 f 3 þ b3 f 2 f 00 þ b4 f ðf 0 Þ2  ¼ 0:

ð91Þ

Assume the solution of Eq. (91) in the form

f ¼ f0 þ pf1 þ p2 f2 þ    ;

ð92Þ

substituting Eq. (92) into Eq. (91) and equating coefficients of like powers of p yields the following equations

p0 : f000 þ b1 f0 ¼ 0;

ð93Þ

p1 : f100 þ b1 f1 þ b2 f03 þ b3 f02 f000 þ b4 f0 f002 ¼ 0; 2

p :

f200

þ b1 f2 þ

3b2 f02 f1

þ

b3 ðf02 f100

þ

2f 0 f1 f000 Þ

ð94Þ þ

b4 f1 f002

¼ 0; . . . :

ð95Þ

For example, the solution of Eqs. (93)–(95) with the initial conditions f(0) = 0.0002, f0 (0) = 0, m = 0.3, r = 0.1 can be written as follows:

52

M.M. Rashidi et al. / Computers and Structures 106–107 (2012) 46–55

Fig. 13. The second mode shape of oscillation for different time (f(0) = 0.0002, f0 (0) = 0, m = 0.3, r = 0.1, a = 10, h = 1, t = 15).

Fig. 16. Other mode shape of oscillation for f(0) = 0.0002, f0 (0) = 0, m = 0.3, r = 0.1, a = 10, h = 1, t = 1.

Fig. 14. The second mode shape of oscillation for different time (f(0) = 0.0002, f0 (0) = 0, m = 0.3, r = 0.1, a = 10, h = 1, t = 25).

Fig. 17. Other mode shape of oscillation for f(0) = 0.0002, f0 (0) = 0, m = 0.3, r = 0.1, a = 10, h = 1, t = 5.

Fig. 15. Other mode shape of oscillation for f(0) = 0.0002, f0 (0) = 0, m = 0.3, r = 0.1, a = 10, h = 1, t = 0.1.

Fig. 18. Other mode shape of oscillation for f(0) = 0.0002, f0 (0) = 0, m = 0.3, r = 0.1, a = 10, h = 1, t = 13.6.

M.M. Rashidi et al. / Computers and Structures 106–107 (2012) 46–55

53

Fig. 19. Other mode shape of oscillation for f(0) = 0.0002, f0 (0) = 0, m = 0.3, r = 0.04, a = 10, h = 1, t = 25.

Fig. 22. Other mode shape of oscillation for f(0) = 0.0002, f0 (0) = 0, m = 0.3, r = 0.1, a = 15, h = 1, t = 25.

Fig. 20. Other mode shape of oscillation for f(0) = 0.0002, f0 (0) = 0, m = 0.3, r = 0.1, a = 5, h = 1, t = 25.

Fig. 23. Other mode shape of oscillation for f(0) = 0.0002, f0 (0) = 0, m = 0.3, r = 0.1, a = 10, h = 2, t = 25.

Fig. 21. Other mode shape of oscillation for f(0) = 0.0002, f0 (0) = 0, m = 0.3, r = 0.1, a = 7, h = 1, t = 25.

Fig. 24. Other mode shape of oscillation for f(0) = 0.0002, f0 (0) = 0, m = 0.3, r = 0.1, a = 10, h = 3, t = 25.

54

M.M. Rashidi et al. / Computers and Structures 106–107 (2012) 46–55

Table 1 Comparison of amplitude oscillation for the HPM and the numerical solution and error between the HPM and numerical solution (for f(0) = 0.002, f0 (0) = 0, m = 0.3, r = 0.1). t

Numerical

HPM

Error (Numerical-HPM)

5 10 15 20 25 30

0.0019346 0.00174274 0.00143697 0.00103726 0.000569759 0.0000650041

0.00193463 0.00174279 0.00143703 0.00103732 0.000569807 0.000065044

2.58193  108 4.84862  108 5.84358  108 5.85421  108 4.80009  108 3.9926  108

f0 ðtÞ ¼ 0:0002 cosð0:577041 tÞ; 12

f1 ðtÞ ¼ 1:38526  10

 1:38526  10

cosð1:73112 tÞ;

2

8 r 9 > > < x =

E

6

1

v

6v 1 ry ¼ 24 > ; 1v 0 0 rxy

38 > > ex 7< 7 0 5 ey > > 1v : c

9 > > =

0

xy

2

ðA1Þ

> > ;

and based on the Von Karman nonlinear strain relation which is stated as:

ð96Þ ð97Þ

f2 ðtÞ ¼ 5:70548  1020 cosð0:577041 tÞ þ 2:87842

9 8 08 91 3 > @u0 þ 1 @w0 2 @a > > > @x 2 @x > > > > > @x > > > > B> =C = <  2 7B< C E 6 @b 6 v 1 0 7B @ v 0 þ 1 @w0 C: ry ¼ þ z @y 5B> @y 2 @y C 24 > > > > > : > > ; 1v > > A @ > > > > 1v @b > 0 0 2 rxy : @u0 @ v 0 @w0 @w0 > ; : @@ya þ @x ; þ þ @y @x @x @y 2

8 r 9 > < x > =

1

v

0

ðA2Þ

 1020 cosð1:73112 tÞ þ 2:82706  1020  cosð2:88521 tÞ:

Hooke’s law describes the relations between stress and strain for an isotropic material as [27]:

> :

cosð0:577041 tÞ

12

Appendix A

ð98Þ

Also Eq. (67) is solved numerically. Figs. 3–10 show the amplitude of oscillation with respect to time. Very good agreement has been found between the analytical and numerical results; however the amount of error is shown to increase when the initial amplitude is increased or the ratio of thickness to dimension is decreased. Also it is apparent from inspection of the Figs. 3–10, that when the plate thickness decreases, the period of oscillation is enhanced and the frequency is decreased. Figs. 11–14 show the second mode shape of oscillation for different values of time. In Figs. 15–24 other mode shapes of oscillation for different physical parameters are presented. In Table 1 comparison of amplitude oscillation for the HPM and the numerical solution is presented. These figures demonstrate that HPM clearly captures the dynamic behavior of the plate.

And definition of in plane and shear forces and bending moment by the following relations:

8 9 8 9 > @2 u > 2 > > 1 N @y2 > > > > x < = < = E 6 2 v Ny ¼ @@xu2 ¼ 4 2 > > > : > ; > > 1v > 0 Nxy :  @2 u > ; @x@y

v 1 0

8 3> @u0 þ 1 @w0 2 > > < @x 2  @x 2 v 0 1 @w0 0 7 5 @@y þ 2 @y > 1v > > : @u0 v 0 @w0 2 þ @@x þ @x @y 0

@w0 @y

9 > > > = ; > > > ;

ðA3Þ (

Qx Qy

)

8 9 < @w0 þ a = @x ¼ KGh @w ; : 0þb;

ðA4Þ

@y

2 8 9 1 M > 3 < x > = 6 Eh 6 v My ¼ 2 4 > > : ; 12ð1  v Þ 0 Mxy

v 1 0

38 @ a > > > @x 7< @b 0 7 5> @y > 1v > : @a 0

2

9 > > > = : > > @b > ;

ðA5Þ

þ @x @y

5. Conclusions In this paper, the equations of motion for an isotropic elastic plate are derived incorporating the effects of shear deformation and rotary inertia. With an Airy stress function, these equations are converted to a single coupled equation and a compatibility equation. By the Galerkin method, a nonlinear differential equation with respect to time is obtained for the dynamic behavior of the plate. This equation has nonlinearities in stiffness and inertia. HPM is used to solve this equation and the responses of the system are obtained. It is shown that the analytical results exhibit very good correlation with the numerical results. The advantage of the present solution is that the effect of nonlinearities can be determined accurately. HPM therefore holds significant promise in nonlinear structural dynamics and it is hoped that the present article will encourage further investigations by researchers in this regard. HPM has recently been applied to many nonlinear multi-physical mechanics problems (with magnetic fields) as described by Anwar Bég et al. [18]. It has also proved extremely accurate in peristaltic hydrodynamics [19] and is currently being employed to study carbon nanotubes [20]. As a special case of the homotopy analysis method (HAM) [21,22], HPM is still very versatile and will provide a very useful tool for engineering mechanics research in the 21st century as a parallel to finite element [23], network simulation [24] or finite difference methods [25,26]. Acknowledgments We express our gratitude to the anonymous referees for their constructive reviews of the manuscript and for helpful comments.

One can write the boundary condition for movable boundary condition for which in-plane forces, shear forces and bending moments are zero, by Eqs. (A3)–(A5), one can write boundary condition of the problems as: At x = a/2 and x = a/2: !  2 # 1 @2 u @2 u 1 @w0 @2u  dx ¼ 0; or N x ¼ 2 ¼ 0;  v ðA5Þ Eh @y2 2 @x @x2 @y a=2 " ! #  2 Z þb=2 1 @2u @2 u 1 @w0 @2 u dx ¼ 0; or N xy ¼ v 2 þ 2  v0 ¼ ¼ 0; ðA6Þ Eh 2 @y @y @x @x@y b=2   @w0 @ 2 u @w0 @ 2 u @w0 w0 ¼ 0; or KGh þa þ 2   ¼ 0; ðA7Þ @x @y @x @x@y @y u0 ¼

Z

þa=2

"

a ¼ 0; or Mx ¼

@a @b þ v ¼ 0; @y @x

ðA8Þ

@ a @b þ ¼ 0; @y @x

ðA9Þ

b ¼ 0; or M xy ¼

at y = b/2 and y = b/2: !  2 # 1 @2 u @2 u 1 @w0 @2 u  dx ¼ 0; or N xy ¼  v ¼ 0; 2 2 Eh 2 @y @x @x @x@y a=2 " ! #  2 Z þb=2 1 @2u @2 u 1 @w0 @2 u dx ¼ 0; or N y ¼ 2 ¼ 0; v 2 þ 2  v0 ¼ Eh 2 @y @x @y @x b=2   2 2 @w0 @ u @w0 @ u @w0 þb  þ 2 ¼ 0; w0 ¼ 0; or KGh @y @x@y @x @x @y @ a @b a ¼ 0; or Mxy ¼ þ ¼ 0; @y @x @ a @b b ¼ 0; or M y ¼ v þ ¼ 0; @x @y u0 ¼

Z

þa=2

"

ðA10Þ ðA11Þ ðA12Þ ðA13Þ ðA14Þ

M.M. Rashidi et al. / Computers and Structures 106–107 (2012) 46–55

References [1] Sathyamoorphy M. Nonlinear vibrations of plates: an update of recent research developments. Appl Mech Rev 1996;49:55–62. [2] Bikri KEI, Benamar R, Bennouna M. Geometrically non-linear free vibrations of clamped simply supported rectangular plates. Part I: The effects of large vibration amplitudes on the fundamental mode shapes. Comput Struct 2003;81:2029–43. [3] Ribeiro P, Petyt M. Non-linear free vibration of isotropic plates with internal resonance. Int J Nonlinear Mech 2000;36:263–78. [4] Nayfeh AH, Mook DT. Nonlinear oscillations. USA: John Wiley & Sons, Inc.; 1979. pp. 500-519. [5] Shooshtari A, Khadem SE. A multiple scale method solution for the nonlinear vibration of rectangular plates. Scientia Iranica 2007;14:64–71. [6] He JH. A coupling method for homotopy technique and perturbation technique for nonlinear problem. Int J Non-Linear Mech 2000;35:37–43. [7] He JH. The homotopy perturbation method for non-linear oscillators with discontinuities. Appl Math Comput 2004;151:287–92. [8] He JH. Application of homotopy perturbation method to nonlinear wave equations. Chaos Solitons Fract 2005;26:695–700. [9] Rashidi MM, Shahmohamadi H, Dinarvand S. Analytic approximate solutions for unsteady two-dimensional and axisymmetric squeezing flows between parallel plates. Math Probl Eng 2008:1–13. [10] Rashidi MM, Ganji DD, Dinarvand S. Explicit analytical solutions of the generalized Burger and Burger–Fisher equations by homotopy perturbation method. Numer Methods Partial Diff Equat 2009;25:409–17. [11] Hoseini SH, Pirbodaghi T, Asghari M, Farrahi GH, Ahmadian MT. Nonlinear free vibration of conservative oscillators with inertia and static type cubic nonlinearities using homotopy analysis method. J Sound Vib 2008;316:263–73. [12] Zand MM, Ahmadian MT. Application of homotopy analysis method in studying dynamic pull-in instability of microsystems. Mech Res Commun 2009;36:851–8. [13] Pirbodaghi T, Ahmadian MT, Fesanghary M. On the homotopy analysis method for non-linear vibration of beams. Mech Res Commun 2009;36:143–8. [14] Nayfeh AH, Chin C. Nonlinear normal modes of a cantilever beams. J Vib Acoust 1995;117:474–81.

55

[15] Nayfeh AH, Chin C, Nayfeh SA. On nonlinear normal modes of systems with internal resonance. J Vib Acoust 1996;118:340–5. [16] He JH. Homotopy perturbation technique. Comput Methods Appl Mech Eng 1999;178:257–62. [17] Özisß T, Yıldırım A. A comparative study of He’s homotopy perturbation method for determining frequency–amplitude relation of a non-linear oscillator with discontinuities. Commun Nonlinear Sci Numer Simulation 2007;8(2):243–8. [18] Anwar Bég O, Ghosh SK, Bég TA. Applied magneto-fluid dynamics: modelling and computation. Germany: Lambert Academic Publishing; 2011. pp. 443, ISBN 978-3-8465-0865-7. [19] Tripathi D, Anwar Bég O, Curiel-Sosa J.L. Homotopy semi-numerical simulation of peristaltic flow of generalized oldroyd-b fluids with slip effect. Comput Methods Biomech Biomed Eng , in press. [20] Anwar Bég O, Islam M.N. HPM simulation of double-walled carbon nanotubes in bending. Technical Report, Aero-786, Dept. Engineering and Mathematics, Sheffield Hallam University, July (2011). [21] Mehmood Ahmer, Ali Asif, Takhar HS, Anwar Bég O, Islam MN, Wilson LS. Unsteady Von Kármán swirling flow: analytical study using the homotopy Method. Int J Appl Math Mech 2010;6(2):67–84. [22] Rashidi MM, Anwar Bég O, Rastegari MT. A study of non-Newtonian flow and heat transfer over a non-isothermal wedge using the homotopy analysis method. Chem Eng Commun 2012;199(2):231–56. [23] Rawat S, Bhargava R, Bhargava Renu, Anwar Bég O. Transient magnetomicropolar free convection heat and mass transfer through a non-Darcy porous medium channel with variable thermal conductivity and heat source effects. Proc IMechE Part C- J Mech Eng Sci 2009;223:2341–55. [24] Zueco J, Anwar Bég O. Network numerical analysis of hydromagnetic squeeze film flow dynamics between two parallel rotating disks with induced magnetic field effects. Tribology Int 2010;43:532–43. [25] Anwar Bég O, Prasad VR, Vasu B, Bhaskar Reddy N, Li Q, Bhargava R. Free convection heat and mass transfer from an isothermal sphere to a micropolar regime with Soret/Dufour effects. Int J Heat Mass Transfer 2011;54:9–18. [26] Rashidi MM, Moradi Bastani M, Anwar Bég O. Numerical simulation of axisymmetric supersonic viscous flow over a blunt cone with a diagonal fourth-order finite difference method. Proc Institution Mech Eng Part G: J Aerospace Eng 2012;226(3):310–26. [27] Chia Ch. Nonlinear analysis of plates. New York: McGraw- Hill; 1980.