Nonlinear free vibration analysis of functionally graded beams by using different shear deformation theories

Journal Pre-proof

Nonlinear free vibration analysis of functionally graded beams by using different shear deformation theories Ke Xie , Yuewu Wang , Xuanhua Fan , Tairan Fu PII: DOI: Reference:

S0307-904X(19)30554-2 https://doi.org/10.1016/j.apm.2019.09.024 APM 13028

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

19 March 2019 18 August 2019 12 September 2019

Please cite this article as: Ke Xie , Yuewu Wang , Xuanhua Fan , Tairan Fu , Nonlinear free vibration analysis of functionally graded beams by using different shear deformation theories, Applied Mathematical Modelling (2019), doi: https://doi.org/10.1016/j.apm.2019.09.024

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Inc.

Highlights 

A direct numerical integration method is proposed to obtain the nonlinear frequency of the FGM beam.



LCDA, the proposed deformation assumption, is more suitable for the nonlinear frequency analysis of beams than the GCDA.



The differences of the results, due to different beam theories, are more obvious for frequencies than frequency ratios.



The effects of power-law index, slenderness ratio and maximum deflection on the nonlinear frequency cannot be ignored.

Nonlinear free vibration analysis of functionally graded beams by using different shear deformation theories Ke Xie a,b,*, Yuewu Wang c, Xuanhua Fan a,b, Tairan Fu c a. Institute of Systems Engineering, China Academy of Engineering Physics, Mianyang, Sichuan, 621900, China b. Shock and Vibration of Engineering Material and Structures Key Laboratory of Sichuan Province, Mianyang, Sichuan, 621900, China c. Department of Energy and Power Engineering, Tsinghua University, Beijing, 100084, China *Corresponding author: [email protected] (K. Xie)

Abstract This study investigates the nonlinear free vibration of functionally graded material (FGM) beams by different shear deformation theories. The volume fractions of the material constituents and effective material properties are assumed to be changing in the thickness direction according to the power-law form. The von Kármán geometric nonlinearity has been considered in the formulation. The Ritz method and Lagrange equation are adopted to yield the discrete formulations. A direct numerical integration method for the motion equation in matrix form is developed to solve the nonlinear frequencies of FGM beams. Comparing with the global concordant deformation assumption (GCDA), a new deformation assumption named as local concordant deformation assumption (LCDA) is proposed in this study. The LCDA fits with the real deformation of the vibrating beam better, thus more accurate results of the nonlinear frequency can be expected. In numerical results, the comparison study of the GCDA and LCDA is carried out. In addition, the effects of power-law index, slenderness ratio and maximum deflection for different shear deformation theories and boundary conditions on the nonlinear frequency of the beam are discussed. Key words: Functionally graded material (FGM) beam, nonlinear vibration, different shear deformation beam theories, global concordant deformation assumption (GCDA), local concordant deformation assumption (LCDA)

1. Introduction The concept of functionally graded material (FGM) was originated in Japan in 1984 during the space-plane project. This material is a new type of material whose compositions are designed to vary smoothly and continuously within the solid. The FGMs have been widely used in engineering fields, such as aircraft, space vehicles, automotive industries, biomedical sectors, and nuclear reactors, owing to their

advantages over the conventional homogeneous and laminated composite materials [1,2]. The studies devoted to understanding the free vibration characteristics of FGM beams, plates and shells have drawn researchers’ attentions in the past decades. The linear free vibrations of FGM beams, plates and shells have been extensively investigated [3-14]. Compared with the linear free vibrations, the studies for nonlinear vibrations of FGM structures are relatively limited, thus have attracted increasing research efforts in recent years. Praveen and Reddy [15] analyzed the nonlinear dynamic response of FGM plates subjected to pressure loads and thickness varying temperature fields by using the first-order shear deformation plate theory and the finite element method. Woo et al. [16] studied the nonlinear free vibration behavior of thin rectangular functionally graded plates by an analytical method. Their results indicated that the coupling effects play a major role in dictating the fundamental frequency of FGM plates. Sundararajan et al. [17] studied the nonlinear free vibration of both rectangular and skew FGM plates by using finite element method. The variation of nonlinear frequency ratio with amplitude was highlighted considering various parameters such as gradient index, temperature, thickness and aspect ratios, and skew angle. Ke and his co-authors [18–21] devoted great efforts to investigate the nonlinear vibration behaviors of FGM beams by various beam theories and methods, and found that unlike isotropic homogeneous beams, the FGM beams show different vibration behaviors at positive and negative amplitudes due to the presence of bending–extension coupling in FGM beams. Talha and Singh [22] investigated the large amplitude free flexural vibrations of shear deformable FGM plates by using a higher order shear deformation theory with a special modification in the transverse displacement. Ding and Chen [23] studied the natural frequencies of nonlinear planar vibration of axially moving beams numerically via the fast Fourier transform. They reduced the governing equations of coupled planar of an axially moving beam to two nonlinear models of transverse vibration, and presented numerical schemes for the governing equations via the finite difference method with the simple support boundary condition. Shen [24] investigated the large amplitude vibration behavior of a shear deformable FGM cylindrical shell of finite length embedded in a large outer elastic medium and in thermal environments by using a higher order shear deformation shell theory and considering shell-foundation interaction. Yazdi [25] used the Homotopy perturbation method to analysis the geometrically nonlinear vibrations of thin rectangular laminated FGM plates, and investigated the effect of parameters such as initial deflection, aspect ratio and material properties on frequency ratio. By using differential transformation method, Wattanasakulpong and Ungbhakorn [26] investigated the linear and nonlinear vibration problems of elastically end restrained FGM beams. Ebrahimi and Zia [27] investigated the nonlinear vibration characteristics of FGM Timoshenko beams made of porous material, and pointed out that the vibration behavior of a FGM beams is significantly influenced by the effects of material distribution profile, porosity volume fraction, aspect ratio and mode number. Taeprasartsit [28] presented a finite element model of large amplitude free vibrations of thin FGM beams based on Euler-Bernoulli beam theory and von Kármán

geometric nonlinearity. The time response of each node was assumed to be a harmonic function, and the residual errors due to this assumption were minimized by employing the Galerkin method. Simsek [29] proposed a novel size-dependent beam model for nonlinear free vibration of a FGM nanobeam based on the nonlocal strain gradient theory and Euler-Bernoulli beam theory in conjunction with von Kármán geometric nonlinearity. The partial nonlinear differential equation describing the nonlinear vibration of FG nanobeam was reduced to an ordinary nonlinear differential equation with cubic nonlinearity via Galerkin’s approach, and a closed-form solution was obtained. Sheng and Wang [30] used multiple scales method to determine the primary resonance responses of FGM beams subjected to parametric excitation with consideration of the von Kármán geometric nonlinearity. They examined the effects of the damping coefficient, volume fraction exponent, amplitude and frequency of parametric excitations. Sınır et al [31] derived a mathematical model for the axially FGM Euler-Bernoulli beam with variable section including the von Kármán geometric nonlinearity, and used multiple time scales method and differential quadrature method simultaneously to solve the free and forced vibration problem of the axially FGM beam. Recently, Malekzadeh and his co-authors [32–35] published a series of literatures focusing on the nonlinear free or forced vibrations of FGM microbeams and plates. In the present work, the direct numerical integration method in conjunction with the local concordant deformation assumption (LCDA) is proposed to deal with the nonlinear vibration problem of the FGM beams based on different shear deformation theories. It is assumed that the volume fractions of the material constituents are changing in the thickness direction according to the power-law form, thus the effective material properties are changing in the same pattern. The von Kármán geometric nonlinearity is considered in the formulation. The Ritz method and Lagrange equation are adopted to derive the discrete formulations. A direct numerical integration method for the motion equation in matrix form is developed to yield the nonlinear frequency of the FGM beam. The first mode shape, which is used in solving the nonlinear frequency by the direct numerical integration method, is obtained via static analysis. Comparing with the traditional deformation pattern, named as global concordant deformation assumption (GCDA), a local concordant deformation assumption (LCDA), which fits with the real deformation of the vibrating beam better, is proposed in this study. In numerical study, the comparison of the GCDA and LCDA is performed, proving that the formulation derived by the LCDA is more advantageous than that by the GCDA, and can lead to results that are more accurate. In addition, the effects of power-law index, slenderness ratio and maximum deflection for different shear deformation theories and boundary conditions on the nonlinear frequencies of beams are studied.

2. Functionally graded material beam model 2.1 Theory and formulation A functionally graded material (FGM) beam with length L, width b, and height h is

considered. Fig. 1 shows the schematic of an FGM beam with simply supported (S–S) ends. The volume fractions of the constituents of FGM beam are assumed to vary from the top (z=h/2) to the bottom surface (z=-h/2) of the beam according to a power-law form which is expressed as the following formula

 z 1 V1     h 2

n

(1)

V1  V2  1

(2)

where V1 and V2 are the volume fractions of the constituent at the top and bottom surface, respectively; and n is the power-law index.

z E1 ,υ1 ,ρ1

I

O

h

E2 ,υ2 ,ρ2

x b

I

Section I-I

L Fig. 1. Schematic of a simply-supported FGM beam.

According to the rule of mixture, the effective properties of the FGM beam can be expressed as n

 z 1 P( z )  PV    P2 1 1  PV 2 2  ( P1  P2 )  h 2

(3)

where P1, P2 are the effective material properties (i.e. modulus of elastics E and mass density ρ) of the beam at the top and bottom surface, respectively. P(z) represents the effective material properties of the FGM beam. Poisson’s ratio υ(z) is assumed to be a constant. According to the shear deformation theory, the displacement field can be expressed as follows: w (4a) u  u0 ( x, t )  g ( z ) 0  f ( z ) ( x, t ) x

w  w0 ( x, t )

(4b)

where u and w denote the axial and transverse displacement of the beam at any point, the subscript “0” represents the mid-plane surface (z=0) of the beam, φ represents the shear strain at the mid-plane surface, and f(z) and g(z) represent the shape functions which characterize the transverse shear strain and stress distribution along the thickness of the beam. Various well-known shear deformation theories [36–40] and a new shear deformation theory proposed by the authors recently [41] can be obtained by choosing f(z) and g(z) as follows:

First-order shear deformation beam theory (FSDBT):  f ( z)  z   g ( z)   z

(5a)

Reissner [36]:

 5z  4 z 2   f ( z )  1  2  4  3h    g ( z)   z 

(5b)

Reddy [37]:

  4z2   f ( z )  z 1  2    3h   g ( z)   z 

(5c)

 h z   f ( z )  sin     h    g ( z)   z 

Touratier [38]:

Karama et al. [39]:

Soldatos [40]:

z   f ( z )  ze 2 h    g ( z )   z

(5d)

2

(5e)

 z 1  f ( z )  h sinh    z cosh   h  2   g ( z)   z 

(5f)

2   h2 16  z    f ( z )  z  2 2     Wang et al. [41]:   h +z 25  h     g ( z)   z

(5g)

It is worthnoting that except the shear deformation theories, recently, the 3D beam theory [44–48] has drawn more and more researches’ attentions for static, buckling and vibration analysis of the beams, since it is free of the assumption of shear deformation pattern and useful for arbitrary thickness-to-length ratio. However, the 3D beam theory is more complex in mathematic than the shear deformation theories. For the simplicity and uniformity of mathematical form, we focus on the shear deformation theories in the current paper. Using the von Kármán assumptions and the displacement field (4), the nonlinear strain–displacement relation of the FGM beam are

 xx 

u 1  w     x 2  x 

2

u  2 w0  1  w0   0 g f    2 x x x 2  x 

2

(6a)

u w  z x w   g   1 0  f  x

 xz 

(6b)

where εxx is the normal strain, and γxz is the shear strain. The stress–strain relation for the FGM beam is expressed as

 xx  Q11 ( z ) xx

(7a)

 xz  Q55 ( z ) xz

(7b)

where Q11 ( z )  E ( z ) , Q55 ( z )  E ( z) / 2(1   ) are the elastic constants that vary through the beam thickness, and E(z) is the Young’s modulus. The strain energy U and kinetic energy T of the FGM beam are given as

1 L h2 U    h  xx xx   xz xz bdzdx 2 0 2 T

(8a)

 u 2  w 2  1 L h2  ( z )     h  bdzdx 2 0  2  t   t  

(8b)

Substituting Eqs. (4), (6) and (7) into Eq. (8), the strain energy U and kinetic energy T can be expressed as 4 2   u0 2 1  u0  w0   w0  A  A  A  11   11  11     x  x    x  4  x   2  2 2 2 2  2 A u0  w0  A  w0  w0   A   w0    g g gg    x x 2 x 2  x  x 2   1 L  (9a) U   dx 2 2 0   w  w  u     2 0 0  2 Af 0   A55  A66     2 A56  x  x x x     2 2    2 w0    w0       Af   A ff      2 Agf x  x  x 2 x    x  2 2   u 2   2 w0   u0  2 w0  w0  0  I0   I gg      I0    2I g t xt xt   1 L   t   t   T   dx 2 2 2 0  u  w          2 I  0 0  2 I gf  I ff  f    t t xt t  t 

(9b)

The stiffness components and inertia related terms in Eq. (9) are defined as follows [ A11 Ag Agg A55 A66 A56 Af Agf Aff ] h    2h bQ11 1 g   2

g2

Q55 2  f  Q11

Q55 2  g   1 Q11

Q55 f   g   1 Q11

f

gf

(10a)  f 2 dz 

[I0

Ig

I gg

If

I gf

I ff ]

h

  2h b 1 g 

g2

f

f 2 dz

gf

(10b)

2

The total energy functional (Π) of the FGM beam for free vibration analysis is expressed as follows (11)  U T 2.2 Solution method 2.2.1 Discrete formulations Ritz method is adopted to obtain the approximate solutions of FGM beams, and the displacements and shear deformation at the mid-plane surface of the beam u0(x,t), w0(x,t) and φ(x,t) are expressed in the terms of space-dependent admissible functions and time-dependent generalized coordinates as u0  aT (t )α( x)

(12a)

w0  bT (t ) β ( x)

(12b)

  cT (t )θ ( x)

(12c)

α(x), β(x), and θ(x) are vectors composed of admissible functions of the beam, and they are expressed as

α( x)  1 ( x)  2 ( x)

 m ( x) 

(13a)

β ( x)   1 ( x) 2 ( x)

m ( x) 

(13b)

T

T

θ ( x)  1 ( x) 2 ( x)

 m ( x) 

T

(13c)

where m is the vector dimension. a(t), b(t), and c(t) are vectors composed of unknown generalized coordinates, and they are expressed as

a(t )   a1 (t ) a2 (t )

am (t )

(14a)

b(t )  b1 (t ) b2 (t )

bm (t )

(14b)

c(t )  c1 (t ) c2 (t )

cm (t )

T

T

T

(14c)

The polynomial type of admissible function is adopted as the Ritz approximation functions. The polynomial admissible functions for the FGM beam are written as

x  i ( x)    L

l

r

 x  x 1      L  L j

x 

i 1

k

x x

i ( x)     1      L  L  L

(15a) i 1

(15b)

p

q

x  x x i ( x)     1     L  L L

i 1

(15c)

where l, r, j, k, and p, q are the admissible function indices of axial, shear, and transverse deformation for various boundary conditions of the FGM beam at the left end (x=0) and right end (x=L), respectively. The values of the indices for different boundary conditions are listed in Table 1. Table 1. Admissible function indices for different boundary conditions.

Boundary condition C-C C-S S-S S-C C-F F-C

Left end (x=0) l 1 1 1 1 1 0

Right end (x=L)

j 2 2 1 1 2 0

p 1 1 0 0 1 0

r 1 1 1 1 0 1

k 2 1 1 2 0 2

q 1 0 0 1 0 1

C: Clamped; S: Simply-supported; F: Free. Substituting Eqs. (9) and (10) into the total energy functional (Π) in Eq. (11) for the nonlinear vibration analysis and then using the Lagrange equation as follows: d     0   dt  qi  qi

(16)

The system equation of motion can be obtained as

Mq  ( K L  K NL )q  0

(17)

where qi represents the unknown coefficients (ai, bi, and ci), the over-dot stands for the partial derivative with respect to time. M is the mass matrix of the FGM beam, KL is the linear part of the structural stiffness matrix, and KNL is the nonlinear part of the structural stiffness matrix. The mass matrix is

 M1  T M   M 2   T  M 3 

M2 M4

 M5 

T

M3   M5   M 6 

(18)

in which the submatices in the mass matrix are L

M1   I 0ααT dx 0

L

M2   I g α 0

L

β T dx x

M3   I f αθ T dx 0

(19a) (19b) (19c)

L

M 4   I 0 ββ T  I gg 0

L

M 5   I gf 0

β β T dx x x

(19d)

β T θ dx x

(19e)

L

M 6   I ff θθ T dx

(19f)

0

The linear part of the stiffness matrix is  KL K 2L K 3L   1  L L T L L  K   K2  K4 K5    L T L T L  K 3   K 5  K 6  where the submatrices of the linear part of stiffness matrix are L

K1L   A11 0

α αT dx x x

(21b)

α θ T dx x x

(21c)

2 β 2 βT β β T  A dx 66 x 2 x 2 x x

(21d)

0

L

K 3L   Af 0

L

0

K  L 5

(21a)

α  2 β T dx x x 2

L

K 2L   Ag

K 4L   Agg

(20)

 2 β θ T β T Agf  A56 θ dx 2 x x x

L

0

L

K 6L   A55θθ T  Aff 0

(21e)

θ θ T dx x x

(21f)

The nonlinear part of the stiffness matrix is

K NL

 0  T   K1NL     0

K1NL K 2NL (1)  K 2NL (2)  K 2NL (3)

K 

NL T 3

0   NL  K3   0 

(22)

where the submatrices of the nonlinear part of stiffness matrix are K1NL  

L

0

1 α T β β T A11 b dx 2 x x x

(23a)

K 2NL (1)

 1 β T  2 β β T 1  2 β T β β T b  Ag 2 b  A L 2 g x x 2 x 2 x x x   2 T 0  1 β T β  β 1 β T θ β T  A b  A c  g f x x x 2 2 x x x  2

(23c)

1 β β T T β β T A11 bb dx 2 x x x x

(23d)

1 β β T θ T Af b dx 2 x x x

(23e)

L

0

L

0

(23b)

1 β T α β T A11 a dx 2 x x x

K 2NL (2)   K 2NL (3)  

  dx   

K 3NL  

L

0

The mass matrix M and the linear part of the stiffness matrix KL are independent of the displacement, while the nonlinear part of the stiffness matrix KNL is dependent on it. Combining Eqs. (9a), (20) – (23), the strain energy U can be rewritten as 1 1 1 U  qT K L q  qT K aNL q  qT K bNL q (24) 2 3 4 where

K aNL

 0  T   K1NL     0

K

NL b

0  0 0

K1NL K 2NL (1)  K 2NL (2)

K 

NL T 3

0 K

NL (3) 2

0

0   NL  K3   0 

0 0  0 

Combining Eqs. (9b), (18), and (19), the kinetic energy T can be rewritten as 1 T  qT Mq 2

(25a)

(25b)

(26)

2.2.2 Linear vibration analysis For linear free vibration analysis, the time-dependent generalized vector q can be expressed as follows

q  qeiLt

(27)

where ωL is the natural frequency of the FGM beam, and q is the corresponding mode vector. Substituting Eq. (27) into Eq. (17), by neglecting the nonlinear part of stiffness matrix KNL(1) and KNL(2), the generalized eigenvalue problem for linear free vibration is expressed as follows ( K L  L2 M )q  0

(28)

By solving the eigenvalue equation (28), the linear natural frequency ωL can be obtained. 2.2.3 Nonlinear vibration analysis based on the global concordant deformation assumption For nonlinear free vibration analysis, due to the presence of bending–extension coupling, the FGM beams show different vibration behavior at positive and negative amplitudes [8], thus the maximum deformation generalized vector of the positive half cycle q is different from that of negative half cycle q . The maximum deformation generalized vectors q and q can be calculated in view of energy conservation, which will be discussed later. The time-dependent generalized vector q for the positive and negative half cycles is assumed as

q  q A(t )

for the positive half cycle

(29a)

q  q A(t )

for the negative half cycle

(29b)

where A(t) is the amplitude factor for the large-amplitude nonlinear free vibration analysis, of which the possible value range of A for the positive and negative half cycle is [0,1]. When A=1, the deformation of the FGM beam reaches to the maximum for the positive or negative half cycle. Eq. (29) implies that the axial, transverse, and shear deformations of beam are linearly related with each other. This assumption is named as the global concordant deformation assumption (GCDA) in this paper, comparing with so called the local concordant deformation assumption (LCDA) which will be discussed later. Substituting Eq. (29), (14), and (26) into Eq. (11), the total energy can be rewritten as    U 0 A2  U 1 A3  U 2 A4  T0 A2

for the positive half cycle

(30a)

   U 0 A2  U 1 A3  U 2 A4  T0 A2

for the negative half cycle

(30b)

where U 0 

1 T L q K q 2

1 U 1  qT K aNL q q  q  3 1 U 2  qT K bNL q q  q  4 1 T0  qT Mq 2 1 U 0  qT K L q 2 1 U 1  qT K aNL q q  q  3 1 U 2  qT K bNL q q  q  4

(31a) (31b) (31c) (31d) (31e) (31f) (31g)

1 (31h) T0  qT Mq 2 According to the energy conservation principle, the total energy of the positive half cycle is equal to that of the negative half cycle, i.e. Π+=Π-. When A=1, the deformation of the FGM beam reaches to the maximum for the positive and negative half cycle respectively, i.e. A  0 , consequently

U 0  U 1  U 2  U 0  U 1  U 2      

(32)

Combining Eqs. (30) and (32) leads to

U 0 (1  A2 )  U 1 (1  A3 )  U 2 (1  A4 ) T0

for the positive half cycle

(33a)

U 0 (1  A2 )  U 1 (1  A3 )  U 2 (1  A4 ) T0

for the negative half cycle

(33b)

for the positive half cycle

(34a)

for the negative half cycle

(34b)

A2 

A2 

Eq. (33) can be transformed into

dt 

dt 

T0 dA U 0 (1  A2 )  U 1 (1  A3 )  U 2 (1  A4 ) T0 dA U 0 (1  A2 )  U 1 (1  A3 )  U 2 (1  A4 )

Integrating the right side of Eqs. (34a) and (34b) from 0 to 1 give the half of the positive and negative half cycle period respectively

 NL 

1 T0 dA   NL 0 2 U 0 (1  A2 )  U 1 (1  A3 )  U 2 (1  A4 )

(35a)

 NL 

1 T0 dA   2NL 0 U 0 (1  A2 )  U 1 (1  A3 )  U 2 (1  A4 )

(35b)

where  NL and  NL are the positive and negative half cycle period; NL and NL are the nonlinear frequencies of the FGM beam for the positive and negative half cycle, which can be written as



NL  1

T0 dA

0

U 0 (1  A2 )  U 1 (1  A3 )  U 2 (1  A4 )

2



NL  1

T0 dA

0

U 0 (1  A2 )  U 1 (1  A3 )  U 2 (1  A4 )

2

The nonlinear frequency ωNL for the full cycle is obtained from

(36a)

(36b)

NL

2NLNL   NL 1 1   NL  NL NL 2



(37)



2.2.4 Nonlinear vibration analysis based on the local concordant deformation assumption In Eq. (29), based on the GCDA, the axial, transverse, and shear deformation are linearly related with each other. However, static results show that the axial deformation is quadratically related with the transverse or shear deformation of the beam. Thus, more precisely, time-dependent generalized vector q can be written as  A2 (t )a    q   A(t )b   A(t )c    

for the positive half cycle

(38a)

 A2 (t )a    q   A(t )b   A(t )c    

for the negative half cycle

(38b)

where a , b , c and a , b , c are the axial, transverse, and shear components of the maximum deformation generalized vectors q and q . This assumption is named as the local concordant deformation (LCDA) in this paper. Substituting Eq. (38), (14), and (26) into Eq. (11), the total energy can be rewritten as

   U 0 A2  U 1 A3  U 2 A4

for the positive half cycle

(39a)

for the negative half cycle

(39b)

T0 A2  T1 AA2  T2 A2 A2    U 0 A2  U 1 A3  U 2 A4 T0 A2  T1 AA2  T2 A2 A2 where

U 0

 0 1 T  q 0 2  0 

  0  1 T 3 L T U 1  q  K2  3 2 3 T   K 3L   2

0 K 4L

K 

L T 5

3 L K2 2 K 2NL (1)

K 

NL T 3

 0  K 5L  q  K 6L 

(40a)

3 L K3  2  NL  K3 q   0   q q

(40b)



U 2

 L  2 K1  1 T  4 NL T  q  K1  4 3  0   0 1 T T0  q 0 2  0

4 NL K1 3 4 NL (2) K2  K 2NL (3) 3 0

0 M4

 M5 

T

M2 0 0

0   M 5  q  M 6 

(40d)

M3   0  q  0 

(40e)

 M1 0 0  T2  2qT  0 0 0  q  0 0 0 

(40f)

 0 1 T  q 0 2  0 

 0  K 5L  q  K 6L 

(40g)

3 L K3  2  K 3NL  q   0   q q

(40h)

  0  T 1 3 U 1  qT   K 2L  3 2 3 T   K 3L   2

U 2

(40c)



 0  T T1  qT  M 2   T  M 3 

U 0

 0  0 q   0  q q

 L  2 K1  1 T  4 NL T  q  K1  4 3  0   0 1  T0  qT 0 2  0

0 K 4L

K 

L T 5

3 L K2 2 K 2NL (1)

K 

NL T 3



4 NL K1 3 4 NL (2) K2  K 2NL (3) 3 0

 0  0 q   0  q q

(40i)



0 M4

 M5 

T

0   M 5  q  M 6 

(40j)

 0  T T1  qT  M 2   T  M 3 

M2 0 0

M3   0  q  0 

 M1 0 0  T2  2q  0 0 0  q  0 0 0  T 

(40k)

(40l)

Apply the energy conservation principle as mentioned in the section 2.2.3, we have Π+=Π-. When the FGM beam deforms to the maximum in the positive or negative cycle, one can write A=1 and A  0 , thus U 0  U 1  U 2  U 0  U 1  U 2      

(41)

Combining Eqs. (39) and (41) leads to

U 0 (1  A2 )  U 1 (1  A3 )  U 2 (1  A4 ) A  T0  AT1  A2T2

for the positive half cycle

(42a)

U 0 (1  A2 )  U 1 (1  A3 )  U 2 (1  A4 ) T0  AT1  A2T2

for the negative half cycle

(42b)

for the positive half cycle

(43a)

for the negative half cycle

(43b)

2

A2 

Eq. (42) can be transformed into

dt 

dt 

T0  AT1  A2T2 dA U 0 (1  A2 )  U 1 (1  A3 )  U 2 (1  A4 ) T0  AT1  A2T2 dA U 0 (1  A2 )  U 1 (1  A3 )  U 2 (1  A4 )

Integrating the right side of Eqs. (43a) and (43b) from 0 to 1 give the half of the positive and negative half cycle period respectively 1 T0  AT1  A T2 dA   NL 0 2 U 0 (1  A2 )  U 1 (1  A3 )  U 2 (1  A4 )

(44a)

1 T0  AT1  A T2 dA    NL 0 2 U 0 (1  A2 )  U 1 (1  A3 )  U 2 (1  A4 )

(44b)

 NL  

2

2

NL 

The nonlinear frequencies of the FGM beam for the positive and negative half cycle NL , NL can be written as



NL  1

T0  AT1  A2T2 dA

0

U 0 (1  A2 )  U 1 (1  A3 )  U 2 (1  A4 )

2

(45a)



NL  1

T0  AT1  A2T2 dA

0

U 0 (1  A2 )  U 1 (1  A3 )  U 2 (1  A4 )

2

(45b)

The nonlinear frequency ωNL for the full cycle can be obtained by Eq. (37). 2.2.4 Solution procedure for nonlinear vibration analysis Until now, the only problem left is to obtain the maximum deformation generalized vectors q and q . In this study, the static analysis is adopted to obtain the maximum deformation generalized vectors. When the deformation of the FGM beam reaches to the maximum for the positive or negative half cycle, the kinetic energy all transfers to the strain energy. It is assumed that the deformation of the FGM beam, at the instant A=1 for the positive and negative half cycle, is identical with the deformation when a uniform pressure of certain value p+ or p- is applied along the beam, respectively. For a given maximum deflection of the beam in the positive cycle wmax, the value of p+ can be determined by calculating the nonlinear static equation with Newton-Raphson (N-R) iterative method and adjusting the uniform load pressure until the maximum deflection of the beam approximates to wmax within a selected tolerance. Specifically, for a given initial positive load pressure p=p+, calculate the following equation by N-R iterative method ( K L  K NL )q  F (46) where the load vector F is written as

 0m1   L  F    pβdx   0   0  m  1  

(47)

When the generalized vector q is converged within a selected tolerance, the maximum deflection of the beam w+ is obtained. Scale up the load pressure p+ by multiplying wmax/ w+, and repeat calculation for the nonlinear static equation (46) until |1–wmax/ w+|<10-6, one can write q  q . The strain energy of the FGM beam U+0+ U+1+ U+2 (by the GCDA, or U 0  U 1  U 2 by the LCDA) for the positive deflection wmax is obtained by calculating Eqs. (31a) – (31c) (or Eqs. (40a) – (40c)) with q  q . Also the energy components U+0, U+1, U+2, T+0 (by the GCDA, or U 0 , U 1 , U 2 , T0 , T1 , T2 by the LCDA) for the positive maximum deformation can be obtained. Then

give an initial negative load pressure p-, calculate Eq. (46) with N-R iterative method until convergence, and calculate the strain energy of the FGM beam U-0+ U-1+ U-2 (by the GCDA, or U 0  U 1  U 2 by the LCDA). Scale up the load pressure p- by

multiplying

U

0

U 0  U 1  U 2  / U 0  U 1  U 2 

the

GCDA,

or

 U 1  U 2  / U 0  U 1  U 2  by the LCDA), and repeat the calculation for

Eq. (36) until 1 

1

(by

U

0

U 0  U 1  U 2  / U 0  U 1  U 2 

 106 (by the GCDA, or

 U 1  U 2  / U 0  U 1  U 2   106 by the LCDA), one can write

q  q . Thus the energy components U-0, U-1, U-2, T-0 (by the GCDA, or U 0 , U 1 , U 2 , T0 , T1 , T2 by the LCDA) for the negative maximum deformation are

obtained. The nonlinear frequencies of the FGM beam for the positive and negative half cycle NL , NL can be calculated by Eq. (36) or Eq. (45), depending on which deformation assumption is selected; and the nonlinear frequency ωNL of the beam can be obtained by Eq. (37). The flowchart in Fig. 2 provides a good understanding of the whole procedure of the nonlinear frequency analysis.

Given an initial positive load pressure p+.

Starting Given the maximum deflection in the positive cycle wmax.

End

Using N-R iterative method, solve the equation ( K L  K NL )q  F

Calculate  NL by

NL  Update p+ by multiplying wmax/ w+

2NLNL NL  NL

Calculate the value of w+.

Calculate NL by 

NL  1

T0 dA

0

U0 (1  A )  U1 (1  A3 )  U2 (1  A4 )

1

T0  AT1  A2T2 dA

2

No

2

or |1–wmax/ w+|<10-6

NL  2



U0 (1  A )  U1 (1  A3 )  U2 (1  A4 )

0

2

Yes Calculate NL by Calculate the values of U+0, U+1, U+2 and T+0 (or U 0, U 1, U 2, T0 , T1 and T2 ) with q  q



NL  1

T0 dA

0

U0 (1  A2 )  U1 (1  A3 )  U2 (1  A4 )

1

T0  AT1  A2T2 dA

2

or NL  2

Given an initial negative load pressure p-.

Using N-R iterative method, solve the equation ( K L  K NL )q  F



U0 (1  A2 )  U1 (1  A3 )  U2 (1  A4 )

0

Calculate the values of U-0, U-1, U-2 and T-0 (or U 0, U 1, U 2, T0 , T1 and T2 ) with q  q

Update p- by multiplying U 0  U 1  U 2 or U 0  U 1  U 2 U 0  U 1  U 2 U 0 U 6 1  U 2 1  U  / U  10

No

1

Calculate the value of U-0+U-1+U-2 or U 0  U 1  U 2

or 1 

U 0  U 1  U 2  106 U 0  U 1  U 2

Yes

U 0  U 1  U 2  106 U 0  U 1  U 2

Fig. 2. Solution procedure of the nonlinear frequency analysis.

3. Results and discussion 3.1 Convergence and validation study The convergence study of the present method, for the boundary conditions S-S, S-C, C-C and C-F, is carried out at first. In the convergence study, the Reddy’s shear deformation theory and GCDA are adopted. For a given maximum deflection wmax=r ( r  I / A , I=bh3/12, A=bh), the dimensionless nonlinear-to-linear frequency ratios ωNL/ωL for the isotropic beam with different slenderness ratio L/h and different

polynomial terms number m, are presented in table 2. In table 2, one can see that eight is a suitable number of the polynomial term, which is sufficiently large to obtain accurate results. Thus eight will be taken as the polynomial term number in the following study. Table 2. Convergence of polynomial term number m for isotropic beams with wmax=r.

Boundary condition

m

L/h

S-S

4 6 8 10

5 1.0902 1.0908 1.0908 1.0908

20 1.0891 1.0896 1.0896 1.0896

100 1.0890 1.0895 1.0895 1.0895

S-C

4 6 8 10

1.0495 1.0481 1.0481 1.0481

1.0482 1.0466 1.0466 1.0466

1.0481 1.0465 1.0465 1.0465

C-C

4 6 8 10

1.0296 1.0281 1.0281 1.0281

1.0276 1.0259 1.0259 1.0259

1.0274 1.0257 1.0257 1.0257

C-F

4 1.0020 1.0011 1.0010 6 1.0022 1.0010 1.0009 8 1.0023 1.0010 1.0009 10 1.0023 1.0010 1.0009 Tables 3–5 and Fig. 3 give the nonlinear-to-linear frequency ratios ωNL/ωL for the isotropic beam with different maximum deflection values wmax. The FSDBT is adopted, and the boundary conditions of S-S, S-C, C-C are selected in this example. The results by GCDA and LCDA are presented for comparison. The shear correction factor is considered as ks=5/6 for FSDBT. The slenderness ratio of the beam is L/h=100. The published results by Gunda et al. [42] and Rao et al. [43] are also listed in table 3 and plotted in Fig. 3, for the purpose of validation of the proposed method. Table 3–5 and Fig. 3 show the results in present work demonstrate good agreement with the published results. Moreover, the results by LCDA are slightly larger than the results by GCDA, and are closer to the published results when wmax/r is larger than 2.0. It implies that the results by the GCDA will slightly underestimate the nonlinear frequencies of the beams. Table 3. Nonlinear-to-linear frequency ratios ωNL/ωL for the isotropic beam with different maximum deflection values wmax. L/h=100. S-S boundary condition.

wmax/r 0.0 0.2

Present by GCDA 1.0000 1.0044

Present by LCDA 1.0000 1.0045

Gunda et al. [42] 1.0000 1.0038

Rao et al. [43] 1.0000 1.0038

0.4 0.6 0.8 1.0 2.0 3.0 4.0

1.0155 1.0337 1.0585 1.0896 1.3188 1.6304 1.9882

1.0158 1.0342 1.0595 1.0910 1.3238 1.6401 2.0032

1.0149 1.0333 1.0585 1.0900 1.3237 1.6409 2.0023

1.0149 1.0333 1.0585 1.0900 1.3237 1.6409 2.0023

Table 4. Nonlinear-to-linear frequency ratios ωNL/ωL for the isotropic beam with different maximum deflection values wmax. L/h=100. S-C boundary condition.

wmax/r 0.0 0.2 0.4 0.6 0.8 1.0 2.0 3.0 4.0

Present by GCDA 1.0000 1.0040 1.0094 1.0184 1.0308 1.0465 1.1679 1.3430 1.5525

Present by LCDA 1.0000 1.0040 1.0096 1.0187 1.0314 1.0475 1.1711 1.3492 1.5618

Gunda et al. [42] 1.0000 1.0020 1.0079 1.0176 1.0311 1.0482 1.1810 1.3741 1.6060

Rao et al. [43] 1.0000 1.0018 1.0073 1.0164 1.0289 1.0448 1.1689 1.3506 1.5701

Table 5. Nonlinear-to-linear frequency ratios ωNL/ωL for the isotropic beam with different maximum deflection values wmax. L/h=100. C-C boundary condition.

wmax/r 0.0 0.2 0.4 0.6 0.8 1.0 2.0 3.0 4.0

Present by GCDA 1.0000 1.0043 1.0071 1.0115 1.0178 1.0257 1.0894 1.1869 1.3100

Present by LCDA 1.0000 1.0044 1.0071 1.0117 1.0180 1.0261 1.0906 1.1894 1.3139

Gunda et al. [42] 1.0000 1.0009 1.0036 1.0081 1.0144 1.0224 1.0867 1.1863 1.3130

Rao et al. [43] 1.0000 1.0009 1.0038 1.0085 1.0151 1.0235 1.0908 1.1946 1.3264

Fig. 3. Nonlinear-to-linear frequency ratios ωNL/ωL for the isotropic beam with different maximum deflection values wmax. L/h=100.

Table 6 give the displacements related ratio u w for the isotropic beam at the point x=3/4L subjected to different pressures at the top surface along the whole beam. u and w are the axial and transverse dimensionless displacements normalized by r, which are written as u  u / r and w  w / r . The pressure applied on top surface of the beam is measured by the times of p0, where p0  EI I / A / bL4 . The Reissner’s beam theory is adopted. Table 6 shows that displacements related ratio u w stay almost constant with the variation of p for any given boundary condition, which implies the axial deformation is quadratically related with the transverse deformation. These results demonstrate that the formulation derived by the LCDA is more advantageous than that by the GCDA, since the LCDA is a more proper assumption than the GCDA for the deformation of vibrating beams. The results by the LCDA are more accurate than that by the GCDA in the nonlinear vibration analysis of beam. Thus, the LCDA is adopted in the following study. Table 6. Dimensionless displacements related ratio

u w , for the isotropic beam at the

point x=3/4L, with different pressures. L/h=100. i  1 .

p/p0 1.0 2.0 3.0 4.0 5.0 6.0

S-S 0.047990 0.047990 0.047990 0.047991 0.047991 0.047991

S-C 0.048499 0.048499 0.048499 0.048499 0.048499 0.048499

C-C 0.014634 0.014634 0.014634 0.014634 0.014634 0.014634

C-F 0.041265i 0.041265i 0.041265i 0.041265i 0.041265i 0.041265i

3.2 Nonlinear frequency analysis for FGM beams An FGM beam made of Tungsten (W) and Copper (Cu) is considered. W dominates

the topmost layer and Cu dominates the bottommost layer. Young’s modulus of W and Cu are 411GPa and 120GPa, respectively. The density of W and Cu is 19250 kg/m3 and 8960 kg/m3, respectively. The slenderness ratio is L/h=8. The maximum deflection of the positive cycle is wmax=r. Tables 7–10 present the nonlinear-to-linear frequency ratios ωNL/ωL for the FGM beam with different power-law indexes n by different beam theories, for the boundary conditions of S-S, S-C, C-C and C-F. The results by Reddy’s beam theory for different boundary conditions are plotted in Fig. 4. From tables 7–10, one can see that the results by different high-order beam theories except the FSDBT are almost identical with each other. However, the results by the FSDBT is slightly smaller than the results by other high-order beam theories, although the difference is very small. It implies that the differences of the results, due to different beam theories, are more obvious for linear [10] or nonlinear frequencies than the frequency ratios ωNL/ωL. This conclusion is valid for the following study for the nonlinear frequency analysis. One can expect that the linear frequency of the FGM beam decreases as the power-law index n increases for any boundary condition, since the natural frequency of the beam made of Cu is lower than that of W, and the increase in power-law index n results in the more material proportion of Cu. However, the frequency ratio ωNL/ωL does not follow the same variation trend. Viewing Fig. 4, it is evident that for any given maximum deflection wmax, the frequency ratio ωNL/ωL drops in accordance with the order of boundary condition S-S, S-C, C-C and C-F. Also the sensibility of the frequency ratio ωNL/ωL to the power-law index n drops in accordance with the order of boundary condition S-S, S-C, C-C and C-F. This phenomenon can be explained that the bending–stretching coupling effect, which plays an important role in the nonlinear vibration of FGM beams [18], recedes in accordance with the order of S-S, S-C, C-C and C-F. Moreover, the frequency ratio ωNL/ωL increases at first, then drops, and increases at last as the increase of the power-law index n for the boundary condition of S-S. The frequency ratio ωNL/ωL increases at first, and then drops slowly as the increase of the power-law index n for the boundary condition of S-C and C-C. The second order variation trends of the frequency ratio curves are concordant with each other for the boundary condition S-S, S-C and C-C. However, the frequency ratio ωNL/ωL drops at first, and then increase slowly as the increase of the power-law index n for the boundary condition of C-F, and the second order variation trend of the frequency ratio curve for C-F is different from those for S-S, S-C and C-C. These results indicate that the sensitivity of the nonlinear frequency to the power-law index n is different from that of the linear frequency to the power-law index n for any given boundary condition. In other words, the sensitivity of the nonlinear frequency to the power-law index n depends on the maximum deflection of the positive cycle wmax. Table 7. Nonlinear-to-linear frequency ratios ωNL/ωL for the FGM beam with different power-law indexes n by different beam theories, for the boundary condition of S-S. L/h=8.

wmax=r. n

Beam theory

0.0 0.1 0.2 0.5 1.0 2.0 5.0 10.0

FSDBT

Reissner

Reddy

Touratier

Karama et al.

Soldatos

Wang et al.

1.094547 1.098056 1.099483 1.096237 1.083644 1.068518 1.067388 1.075772

1.094541 1.098024 1.099433 1.096168 1.083638 1.068744 1.067902 1.076161

1.094541 1.098024 1.099433 1.096168 1.083638 1.068744 1.067902 1.076161

1.094535 1.098018 1.099426 1.096160 1.083630 1.068748 1.067918 1.076155

1.094516 1.097999 1.099407 1.096139 1.083609 1.068737 1.067917 1.076131

1.094542 1.098024 1.099433 1.096168 1.083638 1.068743 1.067900 1.076161

1.094527 1.098010 1.099418 1.096151 1.083622 1.068740 1.067909 1.076146

Table 8. Nonlinear-to-linear frequency ratios ωNL/ωL for the FGM beam with different power-law indexes n by different beam theories, for the boundary condition of S-C. L/h=8.

wmax=r. n

0.0 0.1 0.2 0.5 1.0 2.0 5.0 10.0

Beam theory FSDBT

Reissner

Reddy

Touratier

Karama et al.

Soldatos

Wang et al.

1.051291 1.052537 1.053343 1.053936 1.052085 1.047673 1.043781 1.044215

1.051279 1.052498 1.053287 1.053866 1.052073 1.047877 1.044307 1.044636

1.051279 1.052498 1.053287 1.053866 1.052073 1.047877 1.044307 1.044636

1.051271 1.052490 1.053278 1.053856 1.052064 1.047878 1.044323 1.044630

1.051250 1.052469 1.053258 1.053834 1.052041 1.047863 1.044320 1.044607

1.051279 1.052499 1.053287 1.053866 1.052073 1.047876 1.044305 1.044635

1.051260 1.052481 1.053273 1.053853 1.052064 1.047872 1.044315 1.044629

Table 9. Nonlinear-to-linear frequency ratios ωNL/ωL for the FGM beam with different power-law indexes n by different beam theories, for the boundary condition of C-C. L/h=8.

wmax=r. n

0.0 0.1 0.2 0.5 1.0 2.0 5.0 10.0

Beam theory FSDBT

Reissner

Reddy

Touratier

Karama et al.

Soldatos

Wang et al.

1.030555 1.030954 1.031197 1.031285 1.030396 1.028334 1.026421 1.026682

1.030527 1.030900 1.031127 1.031201 1.030367 1.028526 1.026925 1.027064

1.030527 1.030900 1.031127 1.031201 1.030367 1.028526 1.026925 1.027064

1.030516 1.030889 1.031116 1.031188 1.030355 1.028523 1.026934 1.027052

1.030491 1.030865 1.031092 1.031163 1.030328 1.028503 1.026924 1.027022

1.030528 1.030901 1.031128 1.031202 1.030368 1.028526 1.026923 1.027064

1.030509 1.030886 1.031108 1.031185 1.030347 1.028514 1.026925 1.027041

Table 10. Nonlinear-to-linear frequency ratios ωNL/ωL for the FGM beam with different power-law indexes n by different beam theories, for the boundary condition of C-F. L/h=8.

wmax=r.

n

0.0 0.1 0.2 0.5 1.0 2.0 5.0 10.0

Beam theory FSDBT

Reissner

Reddy

Touratier

Karama et al.

Soldatos

Wang et al.

1.004340 1.004245 1.004031 1.003072 1.001556 1.000382 1.001474 1.002943

1.004335 1.004236 1.004020 1.003059 1.001551 1.000406 1.001546 1.003000

1.004335 1.004236 1.004020 1.003059 1.001551 1.000406 1.001546 1.003000

1.004333 1.004235 1.004018 1.003057 1.001549 1.000405 1.001547 1.002998

1.004329 1.004231 1.004014 1.003054 1.001545 1.000402 1.001545 1.002993

1.004335 1.004237 1.004020 1.003060 1.001551 1.000406 1.001546 1.003000

1.004332 1.004234 1.004017 1.003057 1.001548 1.000404 1.001546 1.002996

Fig. 4. Nonlinear-to-linear frequency ratios ωNL/ωL for the FGM beam with different power-law indexes n for different boundary conditions. L/h=8. wmax=r.

Tables 11–14 present the nonlinear-to-linear frequency ratios ωNL/ωL for the FGM beam with different slenderness ratios L/h by different beam theories, for the boundary conditions of S-S, S-C, C-C and C-F. The maximum deflection of the positive cycle is wmax=r; and the power-law index is n=1.0. In linear vibration analysis, one can expect the natural frequency of the beam increases as the increase of the slenderness ratio L/h [10], because the increase of the slenderness ratio L/h results in the decrease of the shear deformation effect of the beam, which positively correlates to the flexibility of the beam. However, the relationship of the nonlinear-to-linear frequency ratio ωNL/ωL and the slenderness ratio L/h of the beam does not follow the same pattern. The results by Reddy’s beam theory for different boundary conditions are plotted in Fig. 5. It is evident that the nonlinear-to-linear frequency ratio ωNL/ωL drops as the increase of slenderness ratio L/h for all boundary conditions, and converges at about L/h=25. These results indicate that the nonlinear frequency of the beam is less sensitive to the slenderness ratio than the linear frequency of the beam. The linear vibration analysis will overestimate the sensibility of the frequency to the

slenderness ratio for large deflection vibration. It is worthnoting that the nonlinear-to-linear frequency ratio ωNL/ωL is insensitive to the slenderness ratio L/h for the boundary condition of C-F, and the nonlinear-to-linear frequency ratio curve of C-F is almost flat. Table 11. Nonlinear-to-linear frequency ratios ωNL/ωL for the FGM beam with different slenderness ratios L/h by different beam theories, for the boundary condition of S-S. wmax=r. n=1.0. L/h

5 10 15 25 50 100

Beam theory FSDBT

Reissner

Reddy

Touratier

Karama et al.

Soldatos

Wang et al.

1.089970 1.082164 1.080693 1.079936 1.079616 1.079536

1.089933 1.082161 1.080692 1.079936 1.079616 1.079536

1.089933 1.082161 1.080692 1.079936 1.079616 1.079536

1.089910 1.082157 1.080690 1.079935 1.079616 1.079536

1.089853 1.082144 1.080685 1.079933 1.079615 1.079536

1.089934 1.082161 1.080692 1.079936 1.079616 1.079536

1.089895 1.082153 1.080688 1.079934 1.079616 1.079536

Table 12. Nonlinear-to-linear frequency ratios ωNL/ωL for the FGM beam with different slenderness ratios L/h by different beam theories, for the boundary condition of S-C. wmax=r. n=1.0. L/h

5 10 15 25 50 100

Beam theory FSDBT

Reissner

Reddy

Touratier

Karama et al.

Soldatos

Wang et al.

1.057967 1.050758 1.049463 1.048808 1.048532 1.048464

1.057894 1.050753 1.049462 1.048807 1.048532 1.048464

1.057894 1.050753 1.049462 1.048807 1.048532 1.048464

1.057864 1.050748 1.049460 1.048807 1.048532 1.048464

1.057801 1.050734 1.049454 1.048805 1.048532 1.048464

1.057895 1.050753 1.049462 1.048807 1.048532 1.048464

1.057866 1.050745 1.049458 1.048806 1.048532 1.048464

Table 13. Nonlinear-to-linear frequency ratios ωNL/ωL for the FGM beam with different slenderness ratios L/h by different beam theories, for the boundary condition of C-C. wmax=r. n=1.0. L/h

5 10 15 25 50 100

Beam theory FSDBT

Reissner

Reddy

Touratier

Karama et al.

Soldatos

Wang et al.

1.036720 1.028842 1.027258 1.026424 1.026067 1.025978

1.036577 1.028829 1.027255 1.026424 1.026067 1.025978

1.036577 1.028829 1.027255 1.026424 1.026067 1.025978

1.036535 1.028822 1.027252 1.026423 1.026067 1.025978

1.036460 1.028806 1.027246 1.026421 1.026067 1.025977

1.036580 1.028830 1.027255 1.026424 1.026067 1.025978

1.036511 1.028816 1.027251 1.026422 1.026067 1.025978

Table 14. Nonlinear-to-linear frequency ratios ωNL/ωL for the FGM beam with different

slenderness ratios L/h by different beam theories, for the boundary condition of C-F. wmax=r. n=1.0. L/h

5 10 15 25 50 100

Beam theory FSDBT

Reissner

Reddy

Touratier

Karama et al.

Soldatos

Wang et al.

1.001802 1.001500 1.001444 1.001415 1.001403 1.001400

1.001770 1.001497 1.001443 1.001415 1.001403 1.001400

1.001770 1.001497 1.001443 1.001415 1.001403 1.001400

1.001762 1.001496 1.001443 1.001415 1.001403 1.001400

1.001749 1.001494 1.001442 1.001415 1.001403 1.001400

1.001771 1.001497 1.001443 1.001415 1.001403 1.001400

1.001759 1.001496 1.001443 1.001415 1.001403 1.001400

Fig. 5. Nonlinear-to-linear frequency ratios ωNL/ωL for the FGM beam with different slenderness ratios L/h for different boundary conditions. wmax=r. n=1.0.

Tables 15–18 present the nonlinear-to-linear frequency ratios ωNL/ωL for the FGM beam with different maximum deflection values of the positive cycle wmax by different beam theories, for the boundary conditions of S-S, S-C, C-C and C-F. The slenderness ratio is L/h=8; and the power-law index is n=1.0. The results by Reddy’s beam theory for different boundary conditions are plotted in Fig. 6. It is evident that the nonlinear-to-linear frequency ratio ωNL/ωL increases as the increase of the maximum deflection wmax for the boundary conditions of S-S, S-C and C-C. However, the nonlinear-to-linear frequency ratio ωNL/ωL increases at the very first, and then drops slightly as the increase of the maximum deflection wmax for the boundary condition of C-F. Moreover, the nonlinear-to-linear frequency ratio ωNL/ωL drops below 1.0 for C-F boundary condition when wmax/r is more than about 3.0, which implies that the nonlinear frequency is less than the linear frequency in that situation. This phenomenon is because of the lack of pull tension, which contributes to the acceleration of vibration, for the C-F boundary condition. Table 15. Nonlinear-to-linear frequency ratios ωNL/ωL for the FGM beam with different

maximum deflection values of the positive cycle wmax by different beam theories, for the boundary condition of S-S. L/h=8. n=1.0. wmax/r

0.0 0.2 0.4 0.6 0.8 1.0 2.0 3.0 4.0

Beam theory FSDBT

Reissner

Reddy

Touratier

Karama et al.

Soldatos

Wang et al.

1.000000 1.003063 1.011241 1.026542 1.050426 1.083644 1.359137 1.721656 2.116012

1.000000 1.003061 1.011238 1.026539 1.050422 1.083638 1.359117 1.721614 2.115943

1.000000 1.003061 1.011238 1.026539 1.050422 1.083638 1.359117 1.721614 2.115943

1.000000 1.003060 1.011237 1.026536 1.050417 1.083630 1.359093 1.721573 2.115884

1.000000 1.003058 1.011232 1.026528 1.050403 1.083609 1.359030 1.721464 2.115729

1.000000 1.003061 1.011238 1.026539 1.050422 1.083638 1.359117 1.721615 2.115944

1.000000 1.003060 1.011235 1.026533 1.050416 1.083628 1.359087 1.721559 2.115867

Table 16. Nonlinear-to-linear frequency ratios ωNL/ωL for the FGM beam with different maximum deflection values of the positive cycle wmax by different beam theories, for the boundary condition of S-C. L/h=8. n=1.0. wmax/r

0.0 0.2 0.4 0.6 0.8 1.0 2.0 3.0 4.0

Beam theory FSDBT

Reissner

Reddy

Touratier

Karama et al.

Soldatos

Wang et al.

1.000000 1.002508 1.008589 1.018867 1.033384 1.052085 1.199892 1.408807 1.649826

1.000000 1.002493 1.008574 1.018853 1.033371 1.052073 1.199894 1.408840 1.649899

1.000000 1.002493 1.008574 1.018853 1.033371 1.052073 1.199894 1.408840 1.649899

1.000000 1.002488 1.008569 1.018847 1.033363 1.052064 1.199877 1.408815 1.649868

1.000000 1.002479 1.008558 1.018833 1.033346 1.052041 1.199820 1.408722 1.649737

1.000000 1.002493 1.008574 1.018853 1.033371 1.052073 1.199894 1.408839 1.649897

1.000000 1.002485 1.008565 1.018844 1.033361 1.052060 1.199869 1.408806 1.649857

Table 17. Nonlinear-to-linear frequency ratios ωNL/ωL for the FGM beam with different maximum deflection values of the positive cycle wmax by different beam theories, for the boundary condition of C-C. L/h=8. n=1.0. wmax/r

0.0 0.2 0.4 0.6 0.8 1.0 2.0

Beam theory FSDBT

Reissner

Reddy

Touratier

Karama et al.

Soldatos

Wang et al.

1.000000 1.003832 1.007201 1.012784 1.020536 1.030396 1.108430

1.000000 1.003785 1.007156 1.012743 1.020501 1.030367 1.108456

1.000000 1.003785 1.007156 1.012743 1.020501 1.030367 1.108456

1.000000 1.003773 1.007144 1.012731 1.020488 1.030355 1.108445

1.000000 1.003754 1.007124 1.012709 1.020464 1.030328 1.108398

1.000000 1.003786 1.007157 1.012744 1.020501 1.030368 1.108455

1.000000 1.003771 1.007141 1.012726 1.020478 1.030345 1.108433

3.0 4.0

1.225629 1.370784

1.225736 1.370980

1.225736 1.370980

1.225727 1.370975

1.225655 1.370870

1.225733 1.370976

1.225697 1.370946

Table 18. Nonlinear-to-linear frequency ratios ωNL/ωL for the FGM beam with different maximum deflection values of the positive cycle wmax by different beam theories, for the boundary condition of C-F. L/h=8. n=1.0. wmax/r

0.0 0.2 0.4 0.6 0.8 1.0 2.0 3.0 4.0

Beam theory FSDBT

Reissner

Reddy

Touratier

Karama et al.

Soldatos

Wang et al.

1.000000 1.001903 1.001860 1.001787 1.001686 1.001787 1.001686 0.998684 0.996195

1.000000 1.001898 1.001854 1.001782 1.001681 1.001782 1.001681 0.998679 0.996189

1.000000 1.001898 1.001854 1.001782 1.001681 1.001782 1.001681 0.998679 0.996189

1.000000 1.001896 1.001852 1.001780 1.001679 1.001780 1.001679 0.998677 0.996186

1.000000 1.001892 1.001849 1.001776 1.001675 1.001776 1.001675 0.998672 0.996182

1.000000 1.001898 1.001854 1.001782 1.001681 1.001782 1.001681 0.998678 0.996189

1.000000 1.001895 1.001852 1.001779 1.001678 1.001779 1.001678 0.998675 0.996186

Fig. 6. Nonlinear-to-linear frequency ratios ωNL/ωL for the FGM beam with different maximum deflection values of the positive cycle wmax for different boundary conditions. L/h=8. n=1.0.

4. Conclusions A direct numerical integration method for dealing with the motion equation in matrix form was proposed to solve the nonlinear frequency of the FGM beam. The first mode shape that obtained from the static analysis was used in solving the nonlinear frequency by the proposed direct numerical integration method. Comparing with the traditional deformation assumption GCDA, a new deformation assumption named as LCDA was proposed in this study. The LCDA fitted with the real deformation of the vibrating beam better, and was able to yield more accurate results

of the nonlinear frequencies of the FGM beams. The comparison study between the GCDA and LCDA was carried out. Further, the effects of power-law index, slenderness ratio and maximum deflection for different beam theories and boundary conditions on the nonlinear frequency of the beam were investigated. According to the numerical results, the following conclusions were reached: 1) The displacements related ratio u w stayed almost constant with the variation of p for any given boundary condition, suggesting that the axial deformation was quadratically related with the transverse deformation. These results demonstrated that the LCDA was a more proper assumption than the GCDA for the free vibration of the beam, and the formulation derived by the LCDA was more advantageous than that by the GCDA. 2) The differences of the results, due to different beam theories, were more obvious for linear or nonlinear frequencies than the frequency ratios ωNL/ωL. 3) For any given maximum deflection wmax, the frequency ratio ωNL/ωL and the sensibility of the frequency ratio dropped in accordance with the order of boundary condition S-S, S-C, C-C and C-F. 4) The sensitivity of the nonlinear frequency to the power-law index n depended on the maximum deflection of the positive cycle wmax. 5) The nonlinear-to-linear frequency ratio ωNL/ωL dropped as the increase of slenderness ratio L/h for all boundary conditions, and converged at about L/h=25. The nonlinear-to-linear frequency ratio ωNL/ωL was insensitive to the slenderness ratio L/h for the boundary condition of C-F. 6) The nonlinear-to-linear frequency ratio ωNL/ωL increased as the increase of the maximum deflection wmax for the boundary conditions of S-S, S-C and C-C. However, the nonlinear-to-linear frequency ratio ωNL/ωL did not follow the same pattern for the boundary condition of C-F as of S-S, S-C and C-C. 7) The nonlinear-to-linear frequency ratio ωNL/ωL dropped below 1.0 for C-F boundary condition when wmax/r was more than about 3.0. Thus the nonlinear frequency was less than the linear frequency in that situation.

Acknowledgements This work was supported by the National Natural Science Foundation of China (Grant no.11872059), the Science Challenge Project (Grant no.TZ2018002), and the Defense Industrial Technology Development Program (Grant no. C1520110002).

References [1] S. Suresh, A. Mortensen, Fundamentals of Functionally Graded Materials, IOM Communications Ltd., London, 1998. [2] S. Suresh, A. Mortensen, Modeling and design of multi-layered and graded materials, Prog. Mater. Sci. 42 (1997) 243–251. [3] Z.Q. Cheng, R.C. Batra, Exact correspondence between eigenvalues of membranes and functionally graded simply supported polygonal plates, J. Sound Vib. 229 (4) (2000) 879–895. [4] A. Chakraborty, S. Gopalakrishnan, J.N. Reddy, A new beam finite element for the

analysis of functionally graded materials, Int. J. Mech. Sci. 45 (2003) 519–539. [5] C.F. Lu, W.Q. Chen, Free vibration of orthotropic functionally graded beams with various end conditions, Struct. Eng. Mech. 13 (2005) 1430–7. [6] L. Wu, Q.S. Wang, I. Elishakoff, Semi-inverse method for axially functionally graded beams with an anti-symmetric vibration mode, J. Sound Vib. 284 (2005) 1190–202. [7] M. Aydogdu, V. Taksin, Free vibration analysis of functionally graded beams with simply supported edges, Mater. Design. 28 (2007) 1651–1656. [8] H.J. Xiang, J. Yang, Free and forced vibration of a laminated FGM Timoshenko beam of variable thickness under heat conduction, Composites: Part B. 39 (2008) 292–303. [9] S. Kapuria, M. Bhattacharyya, A.N. Kumar, Bending and free vibration response of layered functionally graded beams: a theoretical model and its experimental validation, Compos. Struct. 82 (3) (2008), 390–402. [10] M. Simsek, Fundamental frequency analysis of functionally graded beams by using different higher-order beam theories, Nucl. Eng. Des. 240 (2010) 697–705. [11] A.E. Alshorbagy, M.A. Eltaher, F.F. Mahmoud, Free vibration characteristics of a functionally graded beam by finite element method, Appl. Math. Model. 35 (2011) 412– 425. [12] S.H. Arshad, M.N. Naeem, N. Sultana, A.G. Shah, Z. Iqbal, Vibration analysis of bi-layered FGM cylindrical shells, Arch. Appl. Mech. 81 (2011) 319–343. [13] D.S. Mashat, E. Carrera, A.M. Zenkour, S.A. Al Khateeb, M. Filippi, Free vibration of FGM layered beams by various theories and finite elements, Composites: Part B. 59 (2014) 269–278. [14] C.M. Wang, L.L. Ke, A.N. Roy Chowdhury, J. Yang, S. Kitipornchai, D. Fernando, Critical examination of midplane and neutral plane formulations for vibration analysis of FGM beams, Eng. Struct. 130 (2017) 275–281. [15] G.N. Praveen, J.N. Reddy, Nonlinear transient thermoelastic analysis of functionally graded ceramic–metal plates, Int. J. Solids Struct. 35 (1998) 4457–4476. [16] J. Woo, S.A. Meguid, L.S. Ong, Nonlinear free vibration behavior of functionally graded plates, J. Sound Vib. 289 (2006) 595–611. [17] N. Sundaraja, T. Prakash, M. Ganapathi, Nonlinear free flexural vibration of functionally graded rectangular and skew plates under thermal environments, Finite Elem. Anal. Des. 42 (2005) 152–168. [18] S. Kitipornchai, L.L. Ke, J. Yang, Y. Xiang, Nonlinear vibration of edge cracked functionally graded Timoshenko beams, J. Sound Vib. 324 (2009) 962–982. [19] L.L. Ke, J. Yang, S. Kitipornchai, An analytical study on the nonlinear vibration of functionally graded beams, Meccanica, 45 (2010) 743–752. [20] L.L. Ke, J. Yang, S. Kitipornchai, Nonlinear free vibration of functionally graded carbon nanotube-reinforced composite beams, Compos. Struct. 92 (2010) 676–683. [21] L.L. Ke, Y.S. Wang, J. Yang, S. Kitipornchai, Nonlinear free vibration of size-dependent functionally graded microbeams, Int. J. Engineering Sci. 50 (2012) 256–267. [22] M. Talha, B.N. Singh, Large amplitude free flexural vibration analysis of shear deformable FGM plates using nonlinear finite element method, Finite Elem. Anal. Des. 47 (2011) 394–401. [23] H. Ding, L.Q. Chen, Natural frequencies of nonlinear vibration of axially moving beams,

Nonlinear Dyn. 63 (2011) 125–134. [24] H.S. Shen, Nonlinear vibration of shear deformable FGM cylindrical shells surrounded by an elastic medium, Compos. Struct. 94 (2012) 1144–1154. [25] A.A. Yazdi, Homotopy perturbation method for nonlinear vibration analysis of functionally graded plate, ASME J. Vibr. Acoust. 135 (2013) 021012. [26] N. Wattanasakulpong, V. Ungbhakorn, Linear and nonlinear vibration analysis of elastically restrained ends FGM beams with porosities, J. Aerosp. Sci. Technol. 32 (2014) 111–120. [27] F. Ebrahimi, M. Zia, Large amplitude nonlinear vibration analysis of functionally graded Timoshenko beams with porosities, Acta Astronautica, 116 (2015) 117–125. [28] S. Taeprasartsit, Nonlinear free vibration of thin functionally graded beams using the finite element method, Journal of Vibration and Control, 21 (2015) 29–46. [29] M. Simsek, Nonlinear free vibration of a functionally graded nanobeam using nonlocal strain gradient theory and a novel Hamiltonian approach, Int. J. Engineering Sci. 105 (2016) 12–27. [30] G.G. Sheng, X. Wang, Nonlinear vibration of FG beams subjected to parametric and external excitations, Eur. J. Mech. Solid. 71 (2018) 224–223. [31] S. Sınır, M. Çevik, B.G. Sınır, Nonlinear free and forced vibration analyses of axially functionally graded Euler-Bernoulli beams with non-uniform cross-section, Composites: Part B. 148 (2018) 123–131. [32] P. Malekzadeh, S.M. Monajjemzadeh, Nonlinear response of functionally graded plates under moving load, Thin-Walled Struct. 96 (2015) 120–129. [33] P. Malekzadeh, A. Alibeygi Beni, Nonlinear free vibration of in-plane functionally graded rectangular plates, Mech. Adv. Mater. Struct. 22 (2015) 633-640. [34] A.G. Shenas, S. Ziaee, P. Malekzadeh, Nonlinear vibration analysis of pre-twisted functionally graded microbeams in thermal environment, Thin-Walled Struct. 118 (2017) 87-104. [35] A.G. Shenas, S. Ziaee, P. Malekzadeh, Post-buckling and vibration of post-buckled rotating pre-twisted FG microbeams in thermal environment, Thin-Walled Struct. 138 (2019) 335-360. [36] E. Reissner, On transverse bending of plates, including the effect of transverse shear deformation, Int. J. Solids Struct. 11 (1974) 569–573. [37] J.N. Reddy, A simple higher-order theory for laminated composite plates, J. Appl. Mech. 51 (1984) 745–752. [38] M. Touratier, An efficient standard plate theory, Int. J. Eng. Sci. 29 (1991) 901–16. [39] M. Karama, K.S. Afaq, S. Mistou, Mechanical behavior of laminated composite beam by the new multilayered laminated composite structures model with transverse shear stress continuity, Int. J. Solid Struct. 40 (2003) 1525–1546. [40] K.P. Soldatos, A transverse shear deformation theory for homogeneous monoclinic plates, Acta Mech. 94 (1992) 195–220. [41] Y.W. Wang, K. Xie, T.R. Fu, C.L. Shi, Vibration response of a functionally graded graphene nanoplatelet reinforced composite beam under two successive moving masses, Compos. Struct. 209 (2019) 928–939. [42] J.B. Gunda, R.K. Gupta, G.R. Janardhan, G.V. Rao, Large amplitude free vibration

analysis of Timoshenko beams using a relatively simple finite element formulation, Int. J. Mech. Sci. 52 (2010) 1597–1604. [43] G.V. Rao, K.M. Saheb, Concept of coupled displacement field for large amplitude free vibrations of shear flexible beams, ASME J. Vibr. Acoust. 128 (2006) 251–255. [44] A. Bakhshandeh, B. Navayi Neya, P. Nateghi Babagi, Benchmark solution for free vibration analysis of transversely isotropic thick rectangular plates, Acta Mech. 228 (2017) 3977–3995. [45] A. Moslemi, B. Navayi Neya, J. Vaseghi Amiri, Benchmark solution for buckling of thick rectangular transversely isotropic plates under biaxial load, Int. J. Mech. Sci. 131– 132 (2017) 356–367. [46] S. Brischetto, A general exact elastic shell solution for bending analysis of functionally graded structures, Compos. Struct. 175 (2017) 70–85. [47] Y. Gholami, R. Ansari, R. Gholami, H. Rouhi, Nonlinear bending analysis of nanoplates made of FGMs based on the most general strain gradient model and 3D elasticity theory, Eur. Phys. J. Plus, Eur. Phys. J. Plus (2019) 134: 167. [48] Z. Vafakhah, B. Navayi Neya, An exact three dimensional solution for bending of thick rectangular FGM plate, Compos. Part B. 156 (2019) 72–87.