Composite Structures 154 (2016) 219–230

Contents lists available at ScienceDirect

Composite Structures
journal homepage: www.elsevier.com/locate/compstruct

Linear and non-linear dynamic instability of functionally graded plate

subjected to non-uniform loading
Rajesh Kumar a, S.C. Dutta b, S.K. Panda b,⇑
a
Department of Civil Engineering, Indian Institute of Technology Kharagpur, Kharagpur 721302, India
b
Department of Civil Engineering, Indian School of Mines, Dhanbad 826004, India

a r t i c l e i n f o a b s t r a c t

Article history: The present paper deals with the linear and non-linear dynamic instability of functionally graded mate-
Received 8 February 2016 rials (FGMs) plate subjected to non-uniform loading. Due to this loading, the in-plane stress distributions
Revised 24 June 2016 within the FG plate in pre-buckling state are evaluated using plate membrane analysis. The FG plate is
Accepted 20 July 2016
modeled using higher order shear deformation theory (HSDT) considering von-Kármán geometric nonlin-
Available online 22 July 2016
earity. The Galerkin’s method is used to reduce the non-linear governing partial differential equations
into a set of ordinary differential equations describing the plate non-linear dynamic instability.
Keywords:
Bolotin’s method is used to obtain the boundaries of dynamic instability regions. These boundaries are
FG plate
HSDT
traced by the periodic solution of linear ordinary differential equation (Mathieu-Hill equation) with per-
Non-uniform loading iod T and 2T. Effects of power index, span-to-thickness ratio, aspect ratio, static load factor, boundary
Galerkin’s method conditions and various types of linearly varying loadings and parabolically distributed loadings on the
Non-linear dynamic instability linear and nonlinear dynamic instability are studied. The study of linear and non-linear response in stable
and unstable regions are carried out to identify the dynamic instability behaviour such as dependence of
forcing frequency, existence of beats, effect of nonlinearity on the response and influence of the initial
conditions.
Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Shen [9] presented free vibration and dynamic instability charac-
teristic of functionally graded cylindrical panels subjected to com-
The problems like debonding and sudden change in thermo- bined static and periodic axial forces and in thermal environment.
mechanical properties in conventionally used laminated compos- One dimensional differential quadrature approximation in con-
ites can be avoided by using functionally graded materials (FGMs) junction with Galerkin technique and Bolotin’s method is used to
where, the material properties are varied along the thickness direc- analyse the shell panel. Yang et al. [10] conducted dynamic stabil-
tion. FGMs are widely used in aerospace structures, defense related ity analysis of symmetrically laminated FGM rectangular plates
structures, nuclear reactors and other structures with high func- subjected to a uniaxial periodic in plane load and undergoing uni-
tional performance [1]. During their operational life, these struc- form temperature change. In this study authors concluded that the
tures are subjected to different types of non-uniform mechanical plate may even be unstable beyond a small range of dynamic loads
in-plane static and dynamic loadings. For proper design and perfor- when high level static compression was applied. Natarajan et al.
mances of structural systems, studies of dynamic instability due to [11] studied dynamic buckling of functionally graded spherical
non-uniform in-plane loadings are important. caps under suddenly applied loads by using a three-noded axisym-
Extensive research works [1–6] have been done on static and metric curved shell elements based on field consistency approach.
dynamic behaviour of FG plates and shell structures. Ng et al. Lanhe et al. [12] used moving least-square differential quadrature
[7,8] investigated on dynamic instability of FGM plates and cylin- method to study the dynamic instability behaviours of FG plates
drical shells subjected to harmonic in-plane loading. Authors subjected to in-plane periodic loads and thermal field. Pradyumna
graded the material properties in the direction of thickness accord- and Bandyopadhyay [13] reported dynamic instability behaviour
ing to a material volume fraction power law distribution. Yang and of FG shells subjected to in-plane periodic load and temperature
field using higher-order shear deformation theory in conjunction
⇑ Corresponding author. with finite element approach. In this paper authors assumed that
E-mail addresses: [email protected] (R. Kumar), scdind2000@gmail. the FG materials are temperature dependent and graded in thick-
com (S.C. Dutta), [email protected] (S.K. Panda).

http://dx.doi.org/10.1016/j.compstruct.2016.07.050
0263-8223/Ó 2016 Elsevier Ltd. All rights reserved.
220 R. Kumar et al. / Composite Structures 154 (2016) 219–230

ness direction according to power-law distribution in terms of vol- Ceramic Surface

ume fraction of the constituents. Chen et al. [14] studied the
Neutral Plane
dynamic stability characteristic of functionally graded plates under
C
arbitrary periodic loads. Authors assumed the properties of FGPs to
vary continuously across the plate thickness according to a simple
power law. Torki et al. [15] addressed flutter of FGM cylindrical
Metal Surface
shell under distributed axial follower forces. Authors used FSDT
to model the shell and the material properties are assumed to be Fig. 1a. The position of mid plane and neutral plane for functionally graded (FG)
graded in the thickness direction according to power law distribu- plate.
tion using the properties of two base material phases. Sofiyev [16]
investigated dynamic instability of FGM shells by using Bolotin’s
plate is calculated by taking the first moment of modulus of elas-
method in conjunction with Gaterkin’s method. Sofiyer and Kuruo-
ticity with respect to neutral plane being zero, following the liter-
glu [17,18] and [19] worked extensively on dynamic instability of
ature [26,27].
FGM shell structures and sandwich cylindrical shell with FG core
subjected to periodic axial loads. In the same study authors Z h=2

reduced the governing differential equation to Mathieu type equa- EðzÞðz
CÞdz ¼ 0 ð1Þ

h=2
tion using Galerkin’s method. Further, Mathieu equation is solved
by using Bolotin’s method. Asnafi and Abedi [20,21] studied ortho- where, z is the distance from neutral plane. The position of refer-
tropic sigmoid, power law and exponential functionally graded ence plane (neutral plane) from middle plane may be defined as,
plates under lateral stochastic loads. R h=2
It is observed from the above literature survey that, there is no
h=2
EðzÞz dz
C ¼ R h=2 ð2Þ
work available in the published literature on the dynamic instabil- EðzÞ dz

h=2
ity of FG plates subjected to dynamic non-uniform in-plane loads.
The linear and non-linear dynamic instability analysis of a plate Therefore, distance ðzÞ of the top surface and bottom surface from
   
subjected to dynamic non-uniform (linearly varying and paraboli- the reference plane are zt ¼ 2h
C and zb ¼
2h
C , respectively.
cally varying) in-plane loads is studied considering higher order The effective material properties (modulus of elasticity and
shear deformation theory. Since the applied in-plane edge load is density) according to the Voigt rule of mixture (VRM) are given by
non-uniform, in the first step, the plane elasticity problem is solved
by minimizing the membrane strain energy to evaluate the stress EðzÞ ¼ Ec  V c ðzÞ þ Em  V m ðzÞ
ð3Þ
distribution (rx ; ry and sxy ) within the plate [22,23] and [24]. The qðzÞ ¼ qc  V c ðzÞ þ qm  V m ðzÞ
membrane strain energy of the plate is minimized using Ritz
where, subscripts c and m represents the ceramic and metal con-
method. Using the above stress distribution and using Hamilton’s
stituents, respectively. Here, the volume fraction of the ceramic
principle, the governing partial differential equations of plate
and metal phase is described by the following power-law
motion are derived. Adopting multiterm Galerkin’s approximation,
distribution,
the governing partial differential equations are converted into a set
 n
of ordinary differential equations (Mathieu type of equations) for zþC
studying and describing the dynamic instability of plate. Following Vc ¼ þ 0:5 ; Vm ¼ 1
Vc ð4Þ
h
Bolotin [25], the instability regions are determined from the
boundaries of instability, which represents the periodic solutions The volume fraction varies through the thickness based on the
of differential equation with period T and 2T. The study of linear power index n. In the present investigation, the Poisson’s ratio (m)
and non-linear time history response in stable and unstable is assumed to be constant across the thickness.
regions is carried out to identify the dynamic instability behaviour
such as existence of beats, effect of nonlinearity on the response, 2.1. Displacement field
dependence of forcing frequency and influence of the initial condi-
tions. The Newmark’s numerical integration method conjunction Based on higher order shear deformation theory [28], the dis-
with Newton-Raphson method is used to solve the set of non- placement fields may be written as,
linear ordinary differential equation to obtain the linear and non- 2
linear time history response of the FG plate. The industrial applica- u ¼ uo þ zu1 þ z3 ð4=3h Þ½
u1
wo;x 
tion of such study may be quite significant as mentioned in the v ¼ v o þ zu2 þ z3 ð4=3h2 Þ½
u2
wo;y  ð5Þ
abstract. w¼w o

The above displacement fields can be rearranged as proposed by

2. Formulation Soldatos [29]

u ¼ uo
zwo;x þ f ðzÞ/o1
The functionally graded (FG) plate of length a, width b and
thickness h is made by mixing two different materials ceramic v ¼ v o
zwo;y þ f ðzÞ/o2 ð6Þ
and metal in which top surface is ceramic and the bottom surface w ¼ wo
is metal. The coordinates x, y are along the in-plane directions and
where,
z is along the thickness direction (Fig. 1a). Due to unsymmetrical
material distribution (according power-law) along the thickness /01 ¼ u1 þ w0;x ; /02 ¼ u2 þ w0;y and f ðzÞ ¼ z½1
ð4=3Þðz=hÞ 
2
ð7Þ
direction, the middle plane does not coincide with the neutral
plane resulting coupling between stretching and bending deforma- Here u, v and w are displacement components along x, y, z direc-
tions. In the present case, the reference plane is taken as neutral tions respectively at a distance z away from reference plane and
plane for the analysis of FG plate using elementary theory to avoid u0, v0 and w0 are displacement component of a generic point on
the coupling between the stretching and bending deformation. The the reference surface. u1 and u2 are rotations of the cross sections
shift of neutral plane (C) from the middle plane within the FGMs initially perpendicular to the x and y axes respectively. h is the
 R. Kumar et al. / Composite Structures 154 (2016) 219–230 221

thickness of the plate and suffix x represents the differentiation In the above formulation, FG plate stiffnesses Aij ; Bij ; C ij ; Dij ; Eij ; F ij and
with respect to x. Hij are defined as,
Z h=2
C
2.2. Strain-displacement relations ðAij ; Bij ; Dij Þ ¼ Q ij ð1; z; z2 Þdz ði; j ¼ 1; 2; 6Þ

h=2
C
Z h=2
C
The non-linear strain-displacement relations at a distance ‘z’
away from the reference plane of a plate including von-Kármán ðC ij ; Eij ; F ij Þ ¼ Q ij ð1; z; f ðzÞÞf ðzÞdz ði; j ¼ 1; 2; 6Þ ð14Þ

h=2
C
type of nonlinearities are, Z h=2
C
0 0
ex ¼ eox
zwo;xx þ f ðzÞ/o1;x Hij ¼ Q ij f ðzÞf ðzÞdz ði; j ¼ 4; 5Þ

h=2
C
ey ¼ eoy
zwo;yy þ f ðzÞ/o2;y
where,
cxy ¼ coxy
2zwo;xy þ f ðzÞ/o1;y þ f ðzÞ/o2;x ð8a-eÞ 2 3
1 m 0
EðzÞ 6
cxz ¼ u;z þ w;x ¼ f 0 ðzÞ/o1 Q ij ¼ 1
m2 4
m 1 0 7
5 ði; j ¼ 1; 2; 6Þ
cyz ¼ v ;z þ w;y ¼ f 0 ðzÞ/o2 0 0 1
m
2 ð15a-bÞ
" #
1
m
and o o o
e e and c are reference surface strains and are defined as:
x, EðzÞ 0
y xy
Q ij ¼ 1
m2
2
1
m
ði; j ¼ 4; 5Þ
1  2 1  2 0 2
e ¼ þ wo;x ; eoy ¼ v o;y þ wo;y ;
o
x uo;x
2 2
coxy ¼ uo;y þ v o;x þ wo;x wo;y ð9a-cÞ 2.4. Governing equations of motion
0
where, f ðzÞ ¼ dz
d
f ðzÞ The equations of motion of the functionally graded plate are
obtained in terms of force, moment and shear resultants using
2.3. Constitutive relations Hamilton’s principle as
Z t1 
A plane stress constitutive law is used in the present study. The dð1Þ ðU
W
TÞ ¼ 0 ð16Þ
t0
stress and strain are related as,
U, W and T respectively are the strain energy, work done by the
frg ¼ ½Q feg ð10Þ external loads and the kinetic energy of the functionally graded plate
where, frgT ¼ frxx ryy sxz syz sxy g and fegT ¼ fexx eyy cxz cyz cxy g in the time interval t0 and t1 , respectively. dð1Þ stands for first variation.
The von-Kármán strain displacement relations are used to account for
Here, frgT and fegT are the cartesian components of stress ten-
the geometric nonlinearity. The partial differential equations govern-
sor and corresponding strain at a point respectively. The force and
ing the stability of functionally graded plate, can be written as,
moment resultants of the FG plate can be defined as
00 1 0 1 0 11 0 1 ^ xx;x þ N
N ^ xy;y ¼ q u;tt
Nxx M xx Maxx
rxx Z
1

BB C B C B a CC B C
h=2
C ^
Nxy;x þ N^ yy;y ¼ q v ;tt
@@ Nyy A; @ M yy A; @ Myy AA ¼ @ ryy Að1; z; f ðzÞÞdz 1

h=2
C ^ xx w; þ N
Mxx;xx þ 2Mxy;xy þ Myy;yy þ ðN ^ xy w; Þ;
Nxy Mxy M axy sxy x y x
ð17a-eÞ
^ ^
þðNxy w; þ Nyy w; Þ; ¼ q w;tt
ð11Þ x y y 1

Maxx;x þ M axy;y
Q axx ¼ q3 /1;tt

Z h=2
C
0 Maxy;x þ M ayy;y
Q ayy ¼ q3 /2;tt
ðQ axx ; Q ayy Þ ¼ ðsxz ; syz Þf ðzÞdz ð12Þ

h=2
C where
where N xx ; N yy ; Nxy and Mxx ; Myy ; Mxy are the force and moment Z h=2
C Z h=2
C
resultants; Maxx ; M ayy ; Maxy are additional moment resultants due to q1 ¼ qðzÞdz and q3 ¼ qðzÞz2 dz

h=2
C
h=2
C
additional changes of curvature /o1;x /o2;y /o2;x þ /o1;y ; Q axx ; Q ayy are ^ xx ¼ ½Nxx
nxx 
N ð17fÞ
transverse shear force resultants. The force, moment, additional 
and; ^
Nyy ¼ Nyy
nyy
moment and shear resultants are related to strains by

^ xy ¼ N xy
nxy
N

8 9 2 38 o 9 2 38 9 2 38 9
/o1;y >
<
w;xx >
o
>
< Nxx > = A11 A12 A16 > < exx >= B11 B12 B16 > = C 11 C 12 C 16 > < =
6 7 6 7 6 7
Nyy ¼ 4 A12 A22 A26 5 eoyy þ 4 B12 B22 B26 5
wo;yy þ 4 C 12 C 22 C 26 5 /o2;y
>
: >
; >
: co > ; >
:
2wo > ; >
: o >
;
Nxy A16 A26 A66 xy B16 B26 B66 xy C 16 C 26 C 66 /1;y þ /o2;x
8 9 2 38 o 9 2 38 9 2 38 9
/o1;x >
<
w;xx >
o
>
< Mxx > = B11 B12 B16 > < exx > = D11 D12 D16 > = E11 E12 E16 > < =
6 7 eo 6 7
wo 6 7 o
Myy ¼ 4 B12 B22 B26 5 yy þ 4 D12 D22 D26 5 ;yy þ 4 E12 E22 E26 5 /2;y
>
: >
; >
: co > ; >
:
2wo > ; >
: o >
;
Mxy B16 B26 B66 xy D16 D26 D66 xy E16 E26 E66 /1;y þ /o2;x ð13a-dÞ
8 a 9 2 38 o 9 2 38 9 2 38 9
/o1;x >
<
w;xx >
o
< Mxx >
> = C 11 C 12 C 16 > < exx >= E11 E12 E16 > = F 11 F 12 F 16 >< =
Myy ¼ 6
a 7 6 7
4 C 12 C 22 C 26 5 eyy þ 4 E12 E22 E26 5
w;yy
o o 6 7
þ 4 F 12 F 22 F 26 5 /o2;y
>
: a > ; >
: co > ; >
:
2wo > ; >
: o >
;
Mxy C 16 C 26 C 66 xy E16 E26 E66 xy F 16 F 26 F 66 /1;y þ /o2;x
( a ) 
o 
Q yy H44 H45 /2
¼
Q axx H45 H55 /o1
222 R. Kumar et al. / Composite Structures 154 (2016) 219–230

where, nxx ; nyy and nxy are the plate internal force resultants due to For x ¼ 0; a N  xy ¼ 0 N xx ¼ 4N
 0 y ð1
yÞ
applied parabolically distributed in-plane loading and Nxx ; Nyy ; Nxy
b b
ð20Þ
 xy ¼ 0 N
y ¼ 0; b N  yy ¼ 0
are the in-plane force resultants defined in Eq. 13(a). Thus,
N ^ yy ; N
^ xx ; N ^ xy are the resultant plate internal force resultants. For Following Timoshenko and Goodier [32,30], the stress function
biaxial loading along both principal direction, the biaxial loading is assumed in the form of a series as,
 

N
ratio is denoted as N b N b ¼ N yy
xx
. Here, prebuckling internal force U ¼ U0 þ a1 U1 þ a2 U2 þ a3 U3 þ a4 U4 þ . . . ð21Þ
resultants are the addition of force resultants developed due to both
the in plane edge loading i.e., N xx and N
 yy . Here, assumed stress functions U0 ; U1 ; U2 ; . . . to be such that the
boundary conditions are satisfied. The constant a1 ; a2 ; a3 ; . . . an are
2.5. In-plane stress distribution obtained with the help of Ritz method. Substituting the Eq. (21) in
(18) and carrying out integration, an expression in second degree
In the present investigation, parabolically varying and linearly in a1 ; a2 ; a3 ; . . . an are obtained, minimizing which
varying in-plane compressive dynamic loads as shown in Fig. 1b ai ði ¼ 1; 2; . . . ; nÞ are obtained. The stress function for the parabolic
are considered. The in-plane load distribution can be expressed in-plane edge loading is assumed as follows.
     
by f ðyÞ ¼ 4 by 1
by and f ðyÞ ¼ 1
k by for parabolic and linearly y2 y y2
0
U ¼ 2N
2 2

2 þ ðx2
axÞ ðy2
ybÞ ða1 þ a2 x þ a3 y
varying loadings, respectively. Different types of in-plane load dis- 3 b 2b
tributions such as uniform, trapezoidal, triangular and partial ten-
þ a4 x2 þ a5 xy þ a6 y2 þ . . .Þ ð22Þ
sion are obtained by choosing various values of k = 0, 0.5, 1, and
1.5, respectively. It may be noted that both the static part (N s )
and the dynamic part (N t ) of the in-plane loading varies as the 2.6. Plate buckling analyses
above two non-uniform loading functions. In the case of linearly
varying in-plane load, the stress distribution within the plate coin- The critical buckling load of functionally graded rectangular
cides with the applied edge loading [30]. But, in the case of plates loaded by non-uniform in-plane loading (parabolically vary-
parabolically varying in-plane load the correct stress distribution ing and linearly varying) are obtained by using Galerkin’s method.
within the plate is the one which minimizes the membrane strain In the present investigation following four sets of boundary condi-
energy of the plate and satisfies the boundary condition. The mem- tions are used: SSSS, CSCS, SCSC and CCCC. Where S stands for sim-
brane strain energy of a plate of thickness h of FG plate is given by ply supported edge and C for clamped edge. The letters indicate the
8 9T 2 3
1 8 9 boundary conditions on the edge of the plate in the anti-clockwise
Z Z >< nxx >
= A11 A12 A16 >
< nxx >
= fashion starting from the left hand edge. The out-of-plane displace-
h 6 7 ð18Þ
V¼ nyy 4 A12 A22 A26 5 nyy dxdy ment field wðx; yÞ satisfying the boundary conditions of the plate is
2 A >
: >
; >
: >
;
nxy A16 A26 A66 nxy expressed as the product of beam function.

where, X
1 X
1

Z wðx; yÞ ¼ X m ðxÞY n ðyÞ ð23Þ

@2U @2U @2U h=2
C
m¼1 n¼1
nyy ¼ ; nxx ¼ 2 ; nxy ¼
; Aij ¼ Q ij dz ð19Þ
@x 2 @y @x@y
h=2
C
where X m ðxÞ and Y n ðyÞ are the eigenfunctions of the beam having
Here, U is the stress function. Ritz method [31] is adopted to the same boundary conditions as that of two opposite edges of
minimize the membrane strain energy of the plates in this study. the plate. In the present case following beam functions are assumed
The boundary conditions of the plate membrane problem are given for different edge conditions.
here for parabolically varying in-plane compressive edge loads as a. Simply supported at x = 0 and x = a
follows.

y y
N xx 4N 0 (1 ) b
b b
Parabolic Loading
x
a

0 0.5 1.0 1.5

Uniform Loading Trapezoidal Loading Triangular Loading Partial Tension

Fig. 1b. Parabolic and linearly varying (k ¼ 0; 0:5; 1:0 and 1.5) in-plane edge loading.
 R. Kumar et al. / Composite Structures 154 (2016) 219–230 223

mpx Following sets of displacement fields which satisfy the clamped

X ss
m ðxÞ ¼ sin ðm ¼ 1; 2; 3 . . .Þ ð24Þ
a boundary conditions are given below.
b. Clamped support along two opposite edges, i.e., at x = 0 and X
i X
j
   
x=a ~o ¼
u U mn cos mapx sin npb y
  m¼1 n¼1
x 1 sinðnm =2Þ
X cc
m ðxÞ ¼ cos nm
þ X
i X
a 2 sinhðnm =2Þ j
   
  v~ o ¼ V mn sin mapx cos npb y
x 1
 cosh nm
ðm ¼ 2; 4; 6 . . .Þ ð25aÞ m¼1 n¼1
a 2
X
i X
j

where, nm are the roots of the equation ~o ¼

w W mn X cc cc
m ðxÞX n ðyÞ ð28cÞ
m¼1 n¼1

X
i X
j
   
tanðnm =2Þ þ tanhðnm =2Þ ¼ 0 ð25bÞ ~o ¼
/ K mn sin mapx sin npb y
1
m¼1 n¼1
and
 
x 1 sinðnm =2Þ X
i X
j
   
X cc
m ðxÞ ¼ sin nm

~o ¼
/ Lmn sin mapx sin npb y
a 2 sinhðnm =2Þ 2
  m¼1 n¼1
x 1
 cosh nm
ðm ¼ 3; 5; 7 . . .Þ ð26aÞ Using Eqs. (27a-c) and (28a-c), and employing Galerkin’s method
a 2
the governing differential equation in displacement terms are con-
where, nn are obtained as roots of the equation, verted into a set of homogeneous algebraic equations for simply
tanðnm =2Þ
tanhðnm =2Þ ¼ 0 ð26bÞ supported and clamped boundary conditions, respectively. Solving
the associated eigenvalue problem, critical buckling loads are eval-
The function of Y n ðyÞ are similarly chosen based on the condition uated numerically.2.7. Dynamic instability analysis The applied
at y = 0 and y = b by replacing x by y and a by b and m by n in the non-uniform in-plane load (N x ¼ N s þ N t cos pt) has a static and a
above equations. m and n are, respectively, the number of nodal dynamic component. The static component (N s ) and the dynamic
lines along x and y directions. In the case of simple support only component (N t ) are assumed to vary accordingly either parabolic
normal in-plane displacements are allowed and in-plane tangential or linear distributions. The force, moment and shear resultants in
displacements and out of plane displacements are prevented. Eq. (17a-e) are expressed in terms of displacements and rotation
components. The von-Kármán strain displacement relations are
2.6.1. Boundary conditions at the simply supported edges are used to account for the geometric nonlinearity. The nonlinear dif-
Generalized force boundary conditions are, ferential equations for the case of functionally graded plate sub-
 xx ; Pxx ¼ M xx ¼ 0 at x ¼ 0; a
nxx
Nxx ¼
N jected to non-uniform dynamic loads in terms of displacements
ð27aÞ
 yy ; Pyy ¼ Myy ¼ 0 at y ¼ 0; b
nyy
Nyy ¼
N and rotations component are given in Appendix A. Adopting multi-
term Galerkin’s approximation, the governing partial differential
Generalized displacement boundary conditions are, equations are converted into a set of non-linear ordinary differ-
ential equations describing plate dynamic instability behaviour as,
v o
¼ wo ¼ /oy ¼ 0 at x ¼ 0; a
ð27bÞ
uo ¼ wo ¼ /ox ¼ 0 at y ¼ 0; b
Following sets of displacement fields satisfy the simply sup- ½Mf€dg þ ½½K L þ K NL 
ðNs þ Nt cos ptÞ½K G fdg ¼ f0g ð29aÞ
ported boundary conditions are,
where, ½M, ½K L , ½K NL  and ½K G  are respectively the mass, linear stiff-
X
i X
j
   
~o ¼
u U mn cos mapx sin npb y ness, nonlinear stiffness and geometric stiffness matrices and ‘p’ is
m¼1 n¼1 the excitation frequency. For the study of dynamic instability
X
i X regions using Bolotin [25] method, the nonlinear stiffness (K NL ) is
j
   
v~ o ¼ V mn sin mapx cos npb y neglected from Eq. (29a). After neglecting the non-linear stiffness,
m¼1 n¼1 the differential equation (Mathieu type of equation) is reduced to
X
i X
j the form as follows.
~o ¼
w W mn X ss ss
m ðxÞX n ðyÞ
ð27cÞ
½Mf€dg þ ½K L
ðNs þ Nt cos ptÞ½K G fdg ¼ f0g
m¼1 n¼1
ð29bÞ
X
i X
j
 
~o ¼
/ K mn cosðmapxÞ sin npb y
1
m¼1 n¼1 Here, Ns and Nt are varied as N s ¼ aN cr and N t ¼ bN cr ; where a
X
i X and b are static and dynamic load factors respectively and N cr is
j
mpx  
~o ¼
/ Lmn sin cos npb y the static buckling load. It may be noted that, the effect of non-
2 a
m¼1 n¼1 uniform in-plane load is reflected in the computation of ½K G 
matrix. The Eq. (29a) is a second order differential equation with
2.6.2. The boundary conditions for CCCC plate are periodic coefficients. The critical buckling load is evaluated from
Generalized force boundary conditions are, the solution of linear eigenvalue problem by neglecting the mass,
nonlinear stiffness and time dependant load terms. Similarly the
 xx ;
nxx
Nxx ¼
N at x ¼ 0; a solutions of the eigenvalue problem associated with the Eq. (29a)
 yy ; ð28aÞ
nyy
Nyy ¼
N at y ¼ 0; b neglecting terms containing N s and N t gives the natural frequen-
cies. The regions of instability are located by boundaries of instabil-
Generalized displacement boundary conditions are, ity and on the boundaries, the linear differential equation system
v o
¼ wo ¼ /ox ¼ /oy ¼ 0 at x ¼ 0; a (Eq. (29b)) has periodic solution with period T or 2T. Two solutions
ð28bÞ with same period confine the region of instability and two solu-
uo ¼ wo ¼ /ox ¼ /oy ¼ 0 at y ¼ 0; b
224 R. Kumar et al. / Composite Structures 154 (2016) 219–230

tions with different periods confine the region of stability. The Table 1
solution of the Eq. (29b) with period T and 2T, respectively, are The critical buckling loads a simply supported FG plate (a/b = 1) under uniform in-
plane loading for different span-to-thickness ratios (a/h).
assumed in the form of Fourier series as,
1 
X  Span-to-thickness a/h = 100 a/h = 10
kpt kpt ratio
dðtÞ ¼ b0 þ ak sin þ bk cos ð30Þ
k¼2;4;6
2 2 Power index Present Thai and Choi Present Thai and Choi
[33] [33]

X1   n=0 19.614 19.614 18.518 18.578

kpt kpt
dðtÞ ¼ ak sin þ bk cos ð31Þ n=1 9.777 9.777 9.314 9.339
k¼1;3;5
2 2 n=2 7.629 7.629 7.238 7.263

where, ak and bk are arbitrary constants. Substituting Eq. (30) or

(31) in (29) and equating the coefficients of identical sin kpt
2
and
cos kpt
2
leads to a system of homogeneous algebraic equations in ak
and bk . For a nontrivial solution the determinant of the coefficient
matrix of ak and bk must vanish. The size of the above determinant
is infinite as we have assumed the solution in the form of infinite
series. The determinants are shown to be belonging to a class of
converging determinant known as normal determinant. The first
order and second order approximation to boundaries of principal
regions of instability corresponding to period 2T is obtained by solv-
ing following two eigenvalue problems respectively.

jK   0:5bNcr K G
0:25Mp21 j ¼ 0
    
 K  0:5bNcr K G 0:5bNcr K G   0  ð32a-bÞ
 
p2  0:25M ¼0
 bN K 
K
2:25 Mp 2 2
0 0
cr G 1

Secondary instability region corresponding to period T are

determined from,
 
  2
K
Mp1  ¼ 0
     Fig. 2. Principle and secondary instability zones of simply supported FG plate
K
0:5bNcr K G  0 0 

p 2 ¼0 (a/b = 1, a/h = 100, n = 1, a ¼ 0) subjected to parabolic in-plane loading.

bN K K  1
0 M
cr G
    
K
0:5bNcr K G  M 0
 2 ¼0

0:5bN K 2 
p2  The principal and secondary instability zones of a simply sup-
cr G

K
4 Mp1 0 0 ported FG plate (a/b = 1, a/h = 100, n = 1, a = 0) under parabolic
    
K
0:5bNcr K G 0  0 0 0
    in-plane loading are presented in Fig. 2. Here, the principal zone
    

bNcr K G K
0:5bNcr K G 
p22  0 M 0  ¼ 0 of instability corresponds to period 2T, whereas, secondary zone
   
0
0:5bNcr K G K 
4Mp21  0 0 0 of instability corresponds to period T. It is observed from the figure
that the width of instability increases with the increase of dynamic
ð33a-dÞ
load factor. At a dynamic load factor b = 0.7, the width of instability
 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
where, K ¼ ½K L 
Ns ½K G  2
of principal zone (DP ¼ 3:14 Ec h =qc ) is approximately six times
The linear and non-linear time history responses are obtained
more than the width of secondary instability zone
by using Newmark’s numerical integration scheme [32] by solving qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
2
the set of non-linear and linear ordinary differential Eq. (29a) and (DP ¼ 0:53 Ec h =qc ). The principal zone has greater practical
Eq. (29b), respectively. The non-linear stiffness is considered only importance due to its maximum width of instability. Hence, prin-
for plotting the non-linear time history response and correspond- cipal instability zones with second-order approximation results
ing phase plot. are presented in subsequent figures.
At fixed value of dynamic load factor (b), width of instability is
3. Results and discussions defined by DX ¼ Xu
Xl where, Xu is upper critical values of exci-
tation frequency and Xl is lower critical values of excitation fre-
3.1. Validation of results quency. The unstable region is defined by Xl 6 X 6 Xu . The lower
stability region and upper stability region defined by X < Xl and
Following material properties are used for the functionally X > Xu , respectively. For dynamic load factor b = 0.4, the upper
graded material in the present example problem i.e., Em = 70 Gpa, critical values of excitation frequency (Xu ) is 9.98 and lower critical
Ec = 380 Gpa, qm = 2702 kg/m3, qc = 3800 kg/m3 and m = 0.3, where values of excitation frequency (Xl ) is 8.15 as shown in Fig. 2. Any
m and c represent the material properties corresponding to metal value of dimensionless excitation frequency value between 8.15
and ceramic, respectively. The dynamic instability regions are rep- and 9.98, the FGM plate is unstable for dynamic load factor of 0.40.
resented as a plot of dimensionless excitation frequency
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
2 3.2. Effect of power law and shear deformation
ðX ¼ p qc =Ec h Þ against dynamic load factor (b). To validate the
present formulation, the critical buckling loads of a simply sup-
The influence of power index of power law distribution on
ported FG plate (a/b = 1, a/h = 100) under uniform in-plane loading
the principal instability zones of a simply supported FG plate
for different span-to-thickness ratios are compared with Thai and
(a/b = 1, a/h = 100, a = 0) under parabolic in-plane loading is pre-
Choi [33]in Table 1. It is observed that the present results match
sented in Fig. 3. It is observed that the width of principal zone
well with the reference values.
increases with the decrease of power index (n). The influence of
 R. Kumar et al. / Composite Structures 154 (2016) 219–230 225

Fig. 6. Principal instability zones of FG plates (a/b = 1, a/h = 100, n = 1, a ¼ 0)

Fig. 3. Principal instability zones of simply supported FG plates (a/b = 1, a/h = 100, subjected to parabolic in-plane loading for different boundary conditions.
a ¼ 0) subjected to parabolic in-plane loading for different power law (n).

Fig. 4. Principal instability zones of the simply supported FG plates (a/b = 1, n = 1, Fig. 7. Principal instability zones of simply supported FG plates (a/h = 100, n = 1,
a ¼ 0) subjected to parabolic in-plane loading for different span to thickness (a/h) a ¼ 0) subjected to parabolic in-plane loading for different aspect (a/b) ratio.
ratio.

Fig. 8. Principal and secondary instability zones of simply supported FG plates (a/
Fig. 5. Principal instability zones of simply supported FG plates (a/b = 1, a/h = 100, b = 1, a/h = 100, n = 1) subjected to parabolic in-plane loading for different static
n = 1, a ¼ 0:25) for different types of linearly varying in-plane loadings. load factor (a).
226 R. Kumar et al. / Composite Structures 154 (2016) 219–230

shear deformation on principal instability zone of a simply sup-

ported FG plate (a/b = 1, n = 1, a = 0) subjected to parabolic in-
plane loading is studied in Fig. 4. It is observed from the figure that
the influence of shear deformation on the frequency becomes pro-
nounced when the ratio of side to thickness decreases, resulting
decrease in non-dimensional natural frequency. It is also observed
that as the side to thickness ratio decreases, the width of the insta-
bility zone decreases.

3.3. Effect of different type of in-plane loading

The principal instability zone of a simply supported FG plate

(a/b = 1, a/h = 100, n = 1, a = 0.25) under parabolic and different
types of linearly varying in-plane periodic loadings are presented
in Fig. 5. The static and dynamic parts of different types of in-
plane loads are denoted as N s ¼ aN cr f ðyÞ and N t ¼ bN cr f ðyÞ,
where aN cr and bN cr are the intensity of static and dynamic com-
ponent of compressive load at the edge x = 0, a, and f ðyÞ is in-
plane load distribution for different types of in-plane periodic
Fig. 9. Principal instability zones of simply supported FG plates (a/b = 1, a/h = 100,
n = 1, a ¼ 0) subjected to biaxial parabolic in-plane loading.
loadings, as mentioned earlier. The critical buckling load (N cr )
and dynamic load factor (b) of FG plate under uniform loading

Fig. 10. (a) Linear response and corresponding phase plot, and (b) Non-linear response and corresponding phase plot of simply supported FG plate (a/b = 1, a/h = 100, n = 1,
a ¼ 0 and b ¼ 0:4) under parabolic in-plane loading with excitation frequency parameter (X ¼ 9:10 i.e. Xl < X < Xu) from unstable regions.
 R. Kumar et al. / Composite Structures 154 (2016) 219–230 227

are considered. It is concluded that the width of principal insta- observed from figure that the principal instability zones increases
bility region is the highest for uniform loading and the minimum as the slenderness of the FG plate increases.
for partial tension.
3.5. Effect of static load factor and bi-axial loading
3.4. Effect of edge restraint and aspect ratio
Influence of static load factor on the principal instability zone of
The effect of boundary conditions on the principal instability FG plate (a/b = 1, a/h = 100, n = 1) under parabolic in-plane loading
zones of FG plate (a/b = 1, a/h = 100, n = 1, a = 0) subjected to para- is shown in the Fig. 8. It is clear from the figure that the width of
bolic in-plane loading is presented in Fig. 6. In the present analysis, the primary instability zone increases with increase of static load
four different boundary conditions considered as follows: all edges factor (a). Fig. 9, shows the effect of biaxial parabolic in-plane load-
simply supported (SSSS), loaded edges simply supported and other ing on the principal zone of a simply supported FG plate (a/b = 1,
two edges clamped (SCSC), loaded edges clamped and other two a/h = 100, n = 1, a = 0). The biaxial edge loading is denoted by the
 yy
edges simply supported (CSCS) and all edges clamped (CCCC). It load ratio N b ð¼
N
 xx Þ
N
which is the ratio of compressive (positive) or
shows that the frequency of FG plate decreases under the parabolic
tensile (negative) edge load in the y-direction ðN yy Þ. to compres-
in-plane load in the order: CCCC > CSCS > SCSC > SSSS at the same
sive edge load in x-direction (Nxx). It is clear from the figure that
dynamic load factor ðbÞ ¼ N t =N cr of SSSS FG plate. The width of the width of the principal zone increases with the increase of bi-
instability zone of CCCC plate is the lowest among all boundary axial load factor (N b ).
conditions due to its higher value of natural frequency. The effect
of aspect ratio on the principal instability zone of FG plate 3.6. Response and corresponding phase plot
(a/h = 100, n = 1, a = 0) subjected to parabolic in-plane loading is
presented in Fig. 7. Here, length ‘a’ is kept constant and width ‘b’ Fig. 10 shows (a) the linear response and corresponding phase
is varied to obtain different aspect ratio of the FG plate. It is plot and (b) the non-linear response and corresponding phase plot

Fig. 11. (a) Linear response and corresponding phase plot, and (b) Non-linear response and corresponding phase plot of simply supported four layered FG plate (a/b = 1,
a/h = 100, n = 1, a ¼ 0 and b ¼ 0:4) under parabolic in-plane loading with dimensionless excitation frequency parameter (X ¼ 8:00 i.e. X < Xl) from lower stability region.
228 R. Kumar et al. / Composite Structures 154 (2016) 219–230

Fig. 12. (a) Linear response and corresponding phase plot, and (b) Non-linear response and corresponding phase plot of simply supported FG plate (a/b = 1, a/h = 100, n = 1,
a ¼ 0 and b ¼ 0:4) under parabolic in-plane loading with dimensionless excitation frequency parameter (X ¼ 10:15 i.e. X > Xu) from upper stability region.

respectively of a simply supported FG plate (a/b = 1, a/h = 100, 4. Conclusion

a = 0.0 and b = 0.4) subjected to parabolic in-plane loading. The lin-
ear response of plate is evaluated in unstable region with dimen- In this paper, linear and non-linear dynamic instability of FG
sionless excitation frequency X = 9.10 and initial condition w(0)/ plates subjected to parabolic in-plane loading is studied. The study
h = 0.3. It is observed that the linear response grows exponentially shows that the width of principal instability region is more com-
and the plate becomes unstable, and the non-linear response in the pared to all other instability regions and hence has greater practi-
unstable regions is bounded and also beats are observed. The cal importance. For higher dynamic load factor, second order
amplitude of response is slightly higher than the initial amplitude approximation predicts accurate results. The width of instability
assumed as w(0)/h = 0.3. The linear and non-linear response, and obtained due to uniform load is the maximum compared to all
phase plots corresponding to the dimensionless excitation fre- other types of different in-plane loadings. The width of instability
quency, X = 8.00 ðX < XL Þ and X = 10.15 ðX > Xu Þ are plotted in decreases with the decrease of span to thickness ratio and
Figs. 11 and 12, respectively. It is observed from Fig. 11 that the increases with the increase of aspect ratios. The width of instability
maximum displacement obtained in the linear response is slightly region is the maximum for simply supported boundary conditions
more than non-linear response and the frequency of the beats is in all edges compared to all other combinations of boundary con-
different for both the responses. It is noted from Fig. 12 that in ditions. Also the width of instability region increases with the
the linear response the maximum amplitude is higher than the increase of static and dynamic loads. For values of dimensionless
non-linear response and the frequency of the beats are different excitation frequency, which lies inside the instability region, the
for both the responses. It is also observed from above figures that linear response grows exponentially but the non-linear response
the amplitude of both linear response and non-linear response is bounded. The maximum displacement obtained by linear and
for X = 10.15 are more than that of amplitude of responses for non-linear dynamic response for dimensionless excitation fre-
X = 8.00. quency of X = 10.15 (above the upper instability envelop) is
 R. Kumar et al. / Composite Structures 154 (2016) 219–230 229

approximately three times the initial amplitude of w(0)/h = 0.3. [10] Yang J, Liew KM, Kitipornchai S. Dynamic stability of laminated FGM plates
based on higher-order shear deformation theory. Comput Mech 2004;33
Here, the frequency of beats is different for both linear and nonlin-
(4):305–15.
ear responses. However, for dimensionless excitation frequency of [11] Natarajan S, Thiruvengadam P, Manickam G. Dynamic buckling of functionally
X = 8.00 (below the lower instability envelop), the maximum dis- graded spherical caps. AIAA J 2006;44(5):1097–102.
placement obtained by the linear and nonlinear dynamic responses [12] Lanhe W, Hongjun W, Daobin W. Dynamic stability analysis of FGM plates by
the moving least squares differential quadrature method. Compos Struct
are almost the same as initial amplitude of w(0)/h = 0.3. 2007;77(3):383–94.
[13] Pradyumna S, Bandyopadhyay JN. Dynamic instability of functionally graded
shells using higher-order theory. J Eng Mech 2009;136(5):551–61.
Appendix A. [14] Chen CS, Chen CW, Chen WR. Dynamic stability characteristics of functionally
graded plates under arbitrary periodic loads. Int J Struct Stab Dyn 2013;13(06).
Nonlinear governing partial differential equations of FGM plate http://dx.doi.org/10.1142/S0219455413500260.
[15] Torki ME, Kazemi MT, Reddy JN, Haddadpoud H, Mahmoudkhani S. Dynamic
in displacement variables are
stability of functionally graded cantilever cylindrical shells under distributed
axial follower forces. J Sound Vib 2014;333(3):801–17.

A11 uo;xx þ A66 uo;yy þ ðA12 þ A66 Þv o;xy
fB11 wo;xxx þ ðB12 þ 2B66 Þwo;xyy g þ ðC 11 /o1;xx þ C 66 /o1;yy Þ
ðA-1Þ
þðC 12 þ C 66 Þ/o2;xy þ ðA11 wo;xx þ A66 wo;yy Þwo;x þ ðA12 þ A66 Þwo;y wo;xy
ðnxx;x þ nxy;y Þ ¼ q1 uott
ðA12 þ A66 Þuo;xy þ A66 v o;xx þ A22 v o;yy
fB22 wo;yyy þ ðB12 þ 2B66 Þwo;xxy g þ ðC 12 þ C 66 Þ/o1;xy
ðA-2Þ
þðC 22 /o2;yy þ C 66 /o2;xx Þ þ ðA66 wo;xx þ A22 wo;yy Þwo;y þ ðA12 þ A66 Þðwo;x wo;xy Þ
ðnyy;y þ nxy;x Þ ¼ q1 v ott

B11 fuo;xxx þ wo;x wo;xxy þ ðwo;xx Þ2 g þ B12 fv o;xxy þ uo;xyy þ wo;x wo;xyy þ wo;y wo;xxy þ 2ðwo;xy Þ2 g
þB22 fv o;yyy þ wo;y wo;yyy þ ðwo;yy Þ2 g þ 2B66 fuo;xyy þ v o;xxy þ wo;x wo;xyy þ wo;yy wo;xx þ wo;y wo;xxy þ ðwo;xy Þ2 g

D11 wo;xxxx
2D12 wo;xxyy
D22 wo;yyyy
4D66 wo;xxyy þ E11 /o1;xxx þ E12 ð/o1;xyy þ /o2;xxy Þ þ E22 /o2;yyy
þ2E66 ð/o1;xyy þ /o2;xxy Þ þ wo;xx ðA11 p1 þ A12 q1 þ B11 p2 þ B12 q2 þ C 11 p3 þ C 12 q3 Þ ðA-3Þ
þwo;x ðA11 p1;x þ A12 q1;x þ B11 p2;x þ B12 q2;x þ C 11 p3;x þ C 12 q3;x Þ
nxx wo;xx
nxx;x wo;x
þ2wo;xy ðA66 r 1 þ B66 r 3 þ C 66 r 3 Þ þ wo;x ðA66 r1;y þ B66 r 3;y þ C 66 r 3;y Þ þ wo;y ðA66 r 1;x þ B66 r3;x þ C 66 r 3;x Þ

2nxy wo;xy
nxy;y wo;x
nxy;x wo;y þ wo;yy ðA12 p1 þ A22 q1 þ B12 p2 þ B22 q2 þ C 12 p3 þ C 22 q3 Þ
þwo;y ðA12 p1;y þ A22 q1;y þ B12 p2;y þ B22 q2;y þ C 12 p3;y þ C 22 q3;y Þ
nyy wo;yy
nyy;y wo;y ¼ q1 wott
ðC 12 þ C 66 Þuo;xy þ ðC 22 v o;yy þ C 66 v o;xx Þ þ þðC 12 þ C 66 Þwo;x wo;xy þ ðC 22 wo;yy þ C 66 wo;xx Þwo;y
ðA-4Þ

E22 wo;yyy
ðE12 þ 2E66 Þwo;xxy þ F 66 /o2;xx þ F 22 /o2;yy þ ðF 12 þ F 66 Þ/o1;xy
H44 /o2 ¼ q3 /o1;tt
ðC 11 uo;xx þ C 66 uo;yy Þ þ ðC 12 þ C 66 Þv o;xy þ ðC 11 wo;xx þ C 66 wo;yy Þwo;x þ ðC 12 þ C 66 Þwo;y wo;xy
ðA-5Þ

E11 wo;xxx
ðE12 þ 2E66 Þwo;xyy þ F 11 /o1;xx þ F 66 /o1;yy þ ðF 12 þ F 66 Þ/o2;xy
H55 /o1 ¼ q3 /o2;tt

p1 ¼ uo;x þ 0:5ðwo;x Þ2 ; q1 ¼ v o;x þ 0:5ðwo;y Þ2 ; r 1;x ¼ uo þ v o;x þ wo;x wo;y ; [16] Sofiyev AH. Influences of shear stresses on the dynamic instability of
exponentially graded sandwich cylindrical shells. Compos B Eng
p2 ¼
wo;xx ; q2 ¼
wo;yy ; r 2 ¼ 2wo;xy ; 2015;77:349–62.
[17] Sofiyev AH, Kuruoglu N. Parametric instability of shear deformable sandwich
p3 ¼ /o1;x ; q3 ¼ /o2;y ; r 3 ¼ /o1;y þ /o2;x ; cylindrical shells containing an FGM core under static and time dependent
ðA-6Þ periodic axial loads. Int J Mech Sci 2015;101:114–23.
[18] Sofiyev AH, Kuruoglu N. Dynamic instability of three-layered cylindrical shells
containing an FGM interlayer. Thin Walled Struct 2015;93:10–21.
[19] Sofiyev AH, Kuruoglu N. Domains of dynamic instability of FGM conical shells
References under time dependent periodic loads. Compos Struct 2016;136:139–48.
[20] Asnafi A, Abedi M. A comparison between the dynamic stability of three types
[1] Chi Shyang-Ho, Chung Yen-Ling. Mechanical behavior of functionally graded of nonlinear orthotropic functionally graded plates under random lateral
material plates under transverse load—part I: analysis. Int J Solids Struct loads. J Vib Control 2015. 1077546315617857.
2006;43(13):3657–74. [21] Asnafi A, Abedi M. A complete analogical study on the dynamic stability
[2] Feldman E, Aboudi J. Buckling analysis of functionally graded plates subjected analysis of isotropic functionally graded plates subjected to lateral stochastic
to uniaxial loading. Compos Struct 1997;38(1):29–36. loads. Acta Mech 2015;226(7):2347–63.
[3] Praveen GN, Reddy JN. Nonlinear transient thermoelastic analysis of [22] Dey T, Ramachandra LS. Buckling and postbuckling response of sandwich
functionally graded ceramic-metal plates. Int J Solids Struct 1998;35 panels under non-uniform mechanical edge loadings. Compos B Eng
(33):4457–76. 2014;60:537–45.
[4] Reddy JN. Analysis of functionally graded plates. Int J Numer Methods Eng [23] Panda SK, Ramachandra LS. Postbuckling analysis of cross-ply laminated
2000;47(1–3):663–84. cylindrical shell panels under parabolic mechanical edge loading. Thin Walled
[5] Loy CT, Lam KY, Reddy JN. Vibration of functionally graded CYL shells. Int J Struct 2010;48(8):660–7.
Mech Sci 1999;41:309–24. [24] Kumar R, Kumar A, Panda SK. Parametric resonance of composite skew plate
[6] Pradhan SC, Loy CT, Lam KY, Reddy JN. Vibration characteristics of functionally under non-uniform in-plane loading. Struct Eng Mech 2015;55(2):435–59.
graded cylindrical shells under various boundary conditions. Appl Acoust [25] Bolotin VV. The dynamic stability of elastic systems. San Francisco: Holden-
2000;61(1):111–29. day; 1964.
[7] Ng TY, Lam KY, Liew KM. Effects of FGM materials on the parametric resonance [26] Morimoto T, Tanigawa Y, Kawamura R. Thermal buckling of functionally
of plate structures. Comput Methods Appl Mech Eng 2000;190(8):953–62. graded rectangular plates subjected to partial heating. Int J Mech Sci 2006;48
[8] Ng TY, Lam KY, Liew KM, Reddy JN. Dynamic stability analysis of functionally (9):926–37.
graded cylindrical shells under periodic axial loading. Int J Solids Struct [27] Zhang DG, Zhou YH. A theoretical analysis of FGM thin plates based on
2001;38(8):1295–309. physical neutral surface. Comput Mater Sci 2008;44(2):716–20.
[9] Yang J, Shen HS. Free vibration and parametric resonance of shear deformable [28] Reddy JN, Liu CF. A higher-order shear deformation theory of laminated elastic
functionally graded cylindrical panels. J Sound Vib 2003;261(5):871–93. shells. Int J Eng Sci 1985;23(3):319–30.
230 R. Kumar et al. / Composite Structures 154 (2016) 219–230

[29] Soldatos KP. A refined laminated plate and shell theory with applications. J [32] Subbaraj K, Dokainish MA. A survey of direct time-integration methods in
Sound Vib 1991;144(1):109–29. computational structural dynamics—II. Implicit methods. Comput Struct
[30] Zhong H, Gu C. Buckling of simply supported rectangular Reissner-Mindlin 1989;32(6):1387–401.
plates subjected to linearly varying in-plane loading. J Eng Mech 2006;132 [33] Thai HT, Choi DH. An efficient and simple refined theory for buckling analysis
(5):578–81. of functionally graded plates. Appl Math Model 2012;36(3):1008–22.
[31] Timoshenko SP, Goodier JN. Theory of elastic stability. New York: McGraw-
Hill; 1970.