*Corr. Author’s Address: Guilin University of Electronic Technology, Jinji Road 1#, Guilin, China, xlhuang@guet.edu.cn 571
Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 571-582 Received for review: 2022-07-08
© 2022 The Authors. CC BY 4.0 Int. Licensee: SV-JME Received revised form: 2022-09-02
DOI:10.5545/sv-jme.2022.274 Original Scientific Paper Accepted for publication: 2022-09-07
Nonlinear Vibration Analysis of Functionally Graded  
Porous Plates Reinforced by Graphene Platelets  
on Nonlinear Elastic Foundations
Huang, X.L. – Wang, C.Z. – Wang, J.H. – Wei, N.G.
Xiaolin Huang* – Chengzhe Wang – Jiaheng Wang – Nengguo Wei
Guilin University of Electronic Technology, School of Architecture and Transportation Engineering, China
This paper presents a nonlinear vibration analysis of functionally graded graphene platelet (GPL) reinforced plates on nonlinear elastic 
foundations. Uniformly or non-uniformly distributed internal pores were present in the plates. Based on the modified Halpin-Tsai 
micromechanics model and the extended rule of mixture, the material properties were evaluated. The governing equations, coupled with 
the effect of nonlinear foundations, were derived by using the higher-order plate theory and general von Kármán-type equations. A two-step 
perturbation technique was employed to obtain the nonlinear frequency and transient response. After the present method was verified, the 
effects of pores, GPLs, and elastic foundations were investigated in detail. A new finding is that the influence of the porosity coefficient on the 
natural frequency and dynamic response is relevant to foundation parameters. Moreover, the influence of the nonlinear foundation parameter 
can be negligible.
Keywords: functionally graded porous nanocomposites, graphene platelets, pores, nonlinear elastic foundation, nonlinear vibration, 
transient response
Highlights
• 	 The material properties model of functionally graded (FG) graphene platelets reinforced plates is modified. 
• 	 A two-step perturbation technique for the nonlinear vibration of porous plates on nonlinear elastic foundations is presented. 
• 	 Some interesting conclusions about the effects of pores and nonlinear elastic foundations are drawn.
0  INTRODUCTION
Due to the excellent load-carrying capacity 
with stronger bonding between the matrix and 
carbonaceous nanofillers [1] and [2], nano graphene 
platelets (GPL) are now increasingly used in many 
engineering fields, including aerospace, automobile, 
and civil engineering. It is necessary to study the 
dynamic characteristics of GPL-reinforced composite 
structures. 
In recent years, few results about the dynamic 
behaviour for functionally graded (FG) plates 
reinforced with GPL have been reported [3]. Song 
et al. [4] employed the first-order shear deformation 
plate theory to study the free and forced vibrations 
of FG multilayer GPL/polymer plates. Their results 
illustrated that a small amount of GPL can lead to 
higher natural frequency and lower dynamic response. 
By using the finite element method, Zhao et al. [5] 
studied investigated the bending and free vibration 
behaviours of FG trapezoidal plates. Their results also 
illustrated that the transient deflection was decreased 
by using a small number of graphene platelets. 
Gholami and Ansari [6] investigated the nonlinear 
harmonically excited vibration of FG graphene-
reinforced composite plates. They found that the 
nonlinear buckling was increased with the rise of 
GPLs weight fraction. Moreover, Wu et al. [7] studied 
the dynamic stability of FG nanocomposite plates 
subjected to thermal and mechanical loads.
Recently, a novel kind of graphene-reinforced 
nanocomposite was made by dispersing GPL into 
metal foams [8] to [10]. Due to the excellent mechanical 
and thermal properties, the novel composite materials 
have attracted some attention. Kittipornchai et al. 
[11] proposed a micromechanical model to estimate 
the typical properties of FG graphene reinforced 
nanocomposite, in which both GPL and internal pores 
are uniformly dispersed within each layer. Later, the 
model was employed to study the dynamic behaviour 
of graphene-reinforced nano-composite plates [12] 
and [13]. In these papers, the effect of internal pores 
has been discussed in detail, and the results showed 
that both pore volume fraction and the distribution of 
pores can affect the dynamic characteristics of porous 
structures. 
In practical engineering, plates are sometimes 
rested on elastic foundations. Thus, a suitable model 
is needed to evaluate foundation interaction. The 
simplest is the Winkler model, in which a series 
of springs are used to calculate the tensile and 
compressed forces of elastic foundations. This model 
was later developed into the Pasternak model, in which 
a shear spring is added to estimate the shear forces of 

Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 571-582
572 Huang, X.L. – Wang, C.Z. – Wang, J.H. – Wei, N.G.
the Winkler foundations. Thereafter, the Pasternak 
model has been widely used to study the static and 
dynamic behaviours of plates on elastic foundations 
[14] to [16]. To estimate the nonlinear interaction 
between foundations and plates, some researchers 
recently presented a nonlinear foundation model. For 
example, Nath et al. [17], Civalek [18] and Najafi et 
al. [19] studied the nonlinear dynamic behaviour of 
composite plates by using a nonlinear foundation 
model. Their results revealed that the nonlinear 
foundation parameter had a distinct influence on the 
dynamic characteristics. 
As mentioned above, the number of studies 
focused on the dynamic behaviour of FG porous plates 
reinforced with GPL is still rather scarce. According 
to the authors’ knowledge, no previous work has 
been done to study the nonlinear dynamic behaviour 
of FG graphene platelets reinforced porous plates 
on nonlinear elastic foundations. Hence, this paper 
attempts to study the nonlinear vibration of the plates. 
A modified material properties model is proposed, and 
the effects of internal pores and a nonlinear elastic 
foundation are discussed. 
1  A POROUS NANO-COMPOSITE PLATE
As depicted in Fig. 1, a functionally graded graphene-
reinforced porous plate (length a, width b, thickness 
h) on a nonlinear elastic foundation is taken into 
account. The origin of the coordinate system (X, Y, 
Z) is located at one corner of the middle plane of the 
plate. The Z-axis is perpendicular to the X-Y plane and 
points upwards. Three types of porosity distributions 
are considered (Fig. 2). “P-1” indicates that the 
largest-size pores are distributed in the middle. In 
contrast, “P-2” indicates that the largest-size pores are 
on the bottom and top surfaces. The symbol of even 
distribution is “P-3”.  
Unlike the material properties model presented 
by Kittipornchai et al. [11], the present model is based 
on the volume fraction of pores, which are assumed 
to be:
 VZ
eZ hP
eZ hP
P
P
 
  
  



 
 


0
0
1
12
13
cos/ ,
cos/ ,
,
*



 


. (1)
In the above equation, e
0
, e
0
*
 (0 ≤ e
0
(e
0
*
) < 1 and 
α denote the porosity coefficients for P – 1, P – 2, and 
P – 3 distributions, respectively . 
Yang’s elastic modulus E(Z) and shear elastic 
modulus G(Z) of the porous plate can be expressed by 
using the rule of mixture:
 
EZ EV
GZ GV
P
P
() () ,
() () ,


0
0
1
1 (2)
where E
0
 and G
0
 are the corresponding variations 
of the graphene-reinforced nanocomposites without 
internal pores.
Fig. 1.  A functionally graded GPL-reinforced porous plate on a 
nonlinear elastic foundation
 
a) 
b) 
c) 
Fig. 2.  Three types of porosity distributions;  
a) P – 1, b) P – 2, and c) P – 3
Because the mass density ρ(Z) cannot be 
calculated by using the rule of mixture [11], it is 
assumed to be

Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 571-582
573 Nonlinear Vibration Analysis of Functionally Graded Porous Plates Reinforced by Graphene Platelets on Nonlinear Elastic Foundations  
   

 ()
cos/ ,
(c os /,
*
Z
eZ hP
eZ hP 
  



 
 




 
01
01
11
11 2 
 








0
3
'
,
,
P
 (3)
where ρ
0
 is the mass density of the nanocomposites 
without pores. Coefficients e
1
, e
1
*
 and α′ are the mass 
density coefficients for the three types of porosity 
distributions. 
For open-cell metal foams, the relationship 
between Yang’s elastic modulus and mass density is  
Ez
E
Z () ()
00
2









 [12]. Hence, the coefficients e
1
, e
1
*
 
and α′ can be estimated as 
11 1
11
0














 







e
Z
h
e
Z
h
P
e
Z
h
1
1
*
cosc os ,
cos

  





 












 
  









11 2
3
0
e
Z
h
P
P
*
cos,
,


 

, (4)
It is assumed that the masses of all plates with 
different porosity distributions are equivalent. Hence,  
e
0
*
 and α can be calculated by 
11 1
0
2
2
0
2
2
















e
Z
h
dZ e
Z
h
dZ
h
h
h
h
*
/
/
/
/
(c os cos


 















h
h
h
h
e
Z
h
dZ
/
/
/
/
cos
,
2
2
0
2
2
1

 (5)
Based on the Halpin-Tsai micro-mechanical 
mode, the effective Young’s elastic modulus E
0
 can be 
expressed by
 E
V
V
E
V
V
LLGPL
LG PL
m
WWGPL
WG PL
0
3
8
1
1
5
8
1
1






















E E
m
, (6)
where E
m
 and V
GPL
 are Young’s elastic modulus of 
matrix and GPL volume fraction, respectively. The 
coefficients ξ
L
, ξ
W
, η
L
 and η
W
 are defined by
 




L
GPL
GPL
L
GPL
GPL
L
GPLm
GPLmL
W
GPL
a
h
b
h
EE
EE
EE





22
1
,,
/
/
,
/
m m
GPLmW
EE


1
/
,

 (7)
in which E
GPL
, a
GPL
, b
GPL
 and h
GPL
 are Young’s 
modulus, average length, width, and thickness of 
graphene platelets, respectively. 
Based on the rule of mixture, Poisson’s ratio ν, 
and mass density ρ
0
 can be calculated as
 
 
 
() () ,
() ,
ZV V
VV
GPLG PL mG PL
GPLG PL mG PL
 
 
1
1
0
 (8)
where ρ
GPL
 is the mass density of GPL, and ρ
m
 and ν
m 
are the mass density and Poisson’s ratio of the matrix. 
The shear modulus G
0
 can be obtained by
 G
E
0
0
0
21

  
. (9)
Fig. 3.  Three types of GPL dispersion patterns;  
a)	G	 ― 	 1,	b)	G	 ― 	 2,	and	c)	G	 ― 	 3
As shown in Fig. 3, three types of GPL dispersion 
patterns (G-1, G-2, G-3) are considered. V
GPL
 is 
expressed as 
 Vz
vZ h
vZ h
v
G
G
G
i
i
i
GPL
()
[c os(/ )]
cos( /) ] 
 




 
 

1
1
3
11
2
3


 
, (10)
where ν
i1
, ν
i2 
and ν
i3
 are the maximum value of 
V
GPL
, i (i = 1, 2, 3) indicate the three types of porosity 
distributions. ν
i1
, ν
i2 
and ν
i3
 can be reckoned as 
follows:
 V
W
WW
GPL
T GPLm
GPLm GPLG PL GPL



 
, (11)
V
Z
dZ
ZZ
H
dZ
v
GPL
T
i
h
h
i





 ()
()
cos
(
*
*
/
/
















v
1
2
2
1
1
Z ZZ
H
dZ
v
Z
dZ
h
h
i
h
h
)
cos
()
*
/
/
*
/
/



























2
2
3
2
2
h h
h
/
/
,
2
2

 (12)
in which V
GPL
T
 and W
GPL
 are the total volume and 
weight fractions of GPL, respectively. 
2  FORMULATIONS
2.1  Governing Equations
As noted by Civalek [18], the effect of nonlinear plate-
foundation interaction on the dynamic response of 

Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 571-582
574 Huang, X.L. – Wang, C.Z. – Wang, J.H. – Wei, N.G.
plates on elastic foundations must not be neglected. 
Therefore, the following three-parameter nonlinear 
foundation model is adopted:
RK WK
W
X
W
XY
W
Y
KW
f
















12 3
4
4
4
22
4
4
3
2, (13)
where 
1
K
 ,
2
K
, and 
3
K
 are Winkler, Pasternak, and 
nonlinear foundation parameters. 
According to Reddy’s higher-order thick plate 
theory [20], the displacements of the thick composite 
plate are assumed to be 
uU XYtZ
Z
h
W
X
uV XY
1
2
2
4
3
 
























(,,) ,
(,,

11
ttZ
Z
h
W
Y
uW XYt
yy
),
(,,) ,


























4
3
2
3
 (14)
in which U , V , W , ψ
1
 and ψ
2
 are the displacements 
and rotations of a point (X, Y) on the mid-plane.
The von Karman strains associated with the 
displacement field in Eq. (14) can be stated as



1
2
1
3
2
1
2
2
2
1
2
4
3





























U
X
W
X
Z
X
Z
h X
W
Y
,
V V
Y
W
Y
Z
Y
Z
h Y
W
Y 

























1
2
4
3
0
2
3
2
2
2
34


22

,
,  












22
W
Y
Z
h
W
Y
4
2
2
,



5
2
2 1
6
1
4



























1
W
X
Z
h
W
X
U
Y
V
X
W
X
W
Y
Z
Y
,
 

























2
X
Z
h YX
W
XY
4
3
2
3
2
12
2
. (15)
Based on Hook’s law, the relationship between 
stresses and strains can be expressed as 
 






1
2
6
11 12
12 22
66
1
2
0
0
00























QQ
QQ
Q



6 6
4
5
44
55
4
5
0
0































,
,




Q
Q


 (16)
in which
    
QQ
E
Q
E
QQ QQQ
E
 
  
11 22 2 12 2
16 26 44 55 66
11
0
2




 



,,
,
(1 1)
. (17)
The in-plane forces N
i
, bending moment M
i
, 
higher-order bending moment P
i
, shear forces Q
i
, 
and higher-order shear forces R
i
 are 
  
NMPZ ZZi
QR Z
ii ii
h
h
,, ,, ,( ,,) ,
,,
/
/











11 26
1
3
2
2
22 4
2
d
  d
 d
Z
QR ZZ
h
h
h
h
,
,, .
/
/
/
/







2
2
11 5
2
2
2
1  (18)
By using the Hamilton principle, the equations of 
motion for the plate can be derived as 
   














N
X
N
Y
IU II
h
W
X
N
X
N
Y
IV I
1 6
12 142
6 2
12
4
3




 ,
 
24 2
4
3



I
h
W
Y
,
   


































Q
X
Q
YX
N
W
X
N
W
Y
Y
N
W
X
N
W
Y
12
16
62
4 4
2
12
h
R
X
R
Y











   
















 




4
3
2
4
3
2
2
1
2
2
6
2
2
2
17 2
h
P
X
P
YY
P
Y
R
qI WI
h
f

 














2
2
2
2
2
 
W
X
W
Y
,
   





 












M
X
M
Y
Q
h
R
h
P
X
P
Y
IU II
x
1 6
1 2 1 2
1 6
23 5
44
3



4 4
3
2
h
W
X



,
   





 











 
M
X
M
Y
Q
h
R
h
P
X
P
Y
IV II
y
6 1
2 2 2 2
6 2
23 5
44
3



4 4
3
2
h
W
Y



. (19)
In Eq. (19), the constants I
j
 and I
j
 were given by 
Reddy [20]. The superposed dots indicate the 
differentiation with respect to time. 
The strain compatibility equation is 
 






















2
1
0
2
2
2
0
2
2
6
0 2
2
2
2
2
2
 
YXXY
W
XY
W
X
W
Y
. (20)

Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 571-582
575 Nonlinear Vibration Analysis of Functionally Graded Porous Plates Reinforced by Graphene Platelets on Nonlinear Elastic Foundations  
The in-plane forces N
i
 can be expressed by stress 
function FXYt (,,) :
 N
F
Y
N
F
X
N
F
XY
1
2
2 2
2
2 6
2









,, . (21)
By substituting Eqs. (18) and (21) into Eqs. 
(19) and (20), the governing equations of nonlinear 
vibration for the plate can be derived as follows:
lW lllF R
lWFl WI
X
f 11 12 11 32 14
17 8
1
() () () ()
(,)( )

 






 



 
I
Y
q
lF lllW lWW
l
8
2
21 22 12 32 24
31
1
2



,
() () () () (,),
(W Wl ll FI
W
X
I
lW l
)( )( )( ),
() ()





32 13 32 34 91 01
41 42 1
 



 


 ll FI
W
Y
I
43 24 4910 2
() () , 


 (22)
where the constants I
j
 (  j = 8, 9, 10), linear operators l
ij
 
and nonlinear operator l were given by Shen [21] and 
Huang and Zheng [22].
The four edges of the plate are assumed to be 
simply supported. The boundary conditions are 
expressed as
 
XW MPN
YW MPN
  
  
00
00
11 6
22 26
,,
,.
a:
b:
1

 (23)
2.2  Solution Procedure
To solve the nonlinear equations, Eqs. (22) and 
(23), we first introduce the following dimensionless 
parameters:
x
X
a
y
Y
b
z
Z
h
a
b
W
W
DDAA
F
F
DD
 


 ,, ,,
[]
,
[
****/
*
11 22 11 22
14
11 22
* */ ****/
]
,(,)
(, )
[]
,
(, )
12 12
12
11 22 11 22
14
4
4






a
DDAA
KK
a
12
,, (, )
,,
*
** *
*
a
DKK
K
Ka DDAA
D
q
q
2
2
11
11 22 11 22
11











12
3
3
4*
4
a a
DD DAA
t
a
E
m
m
4
4
11 11 22 11 22
14




* ****/
[]
,, 
x
X
a
y
Y
b
z
Z
h
a
b
W
W
DDAA
F
F
DD
 


 ,, ,,
[]
,
[
****/
*
11 22 11 22
14
11 22
* */ ****/
]
,(,)
(, )
[]
,
(, )
12 12
12
11 22 11 22
14
4
4






a
DDAA
KK
a
12
,, (, )
,,
*
** *
*
a
DKK
K
Ka DDAA
D
q
q
2
2
11
11 22 11 22
11











12
3
3
4*
4
a a
DD DAA
t
a
E
m
m
4
4
11 11 22 11 22
14




* ****/
[]
,,  (24)
where the stiffness constants A
ij
, B
ij
, and D
ij
 are 
defined in the standard way [20].
The dimensionless form of the equation (22) can 
be rewritten as
lW ll lF KW KW
KW lWF
11 12 11 32 11 4
2
3
1
2
() () () ()
(,)
 
 
 

12
3
l lW
xy
lF ll
q 17 2
1
2
2
21 22 21 22 3
() ,
() () (

 













 
2 22 42 42
2
31 32 13 32 13 4
1
2
)( )( ,) ,
() () () (
 

 
 
lW lWW
lW ll lF ) ), 


  
45 1


W
x
 
lW ll lF KW KW
KW lWF
11 12 11 32 11 4
2
3
1
2
() () () ()
(,)
 
 
 

12
3
l lW
xy
lF ll
q 17 2
1
2
2
21 22 21 22 3
() ,
() () (

 













 
2 22 42 42
2
31 32 13 32 13 4
1
2
)( )( ,) ,
() () () (
 

 
 
lW lWW
lW ll lF ) ), 


  
45 1


W
x
 (25)
where the constants γ
i
(i = 1, 2, ..., 5), the dimensionless 
operators l
ij 
and l can be seen in the previous work 
[22].
The dimensionless form of the boundary 
conditions in Eq. (23) are also rewritten as:
 
x:
:
1
  
  
00
00
11 6
22 26
,,
,.


WM PN
yWMPN  (26)
A two-step perturbation technique is employed to 
solve the nonlinear governing equations, Eq. (25). As 
the essence of this procedure, the asymptotic solution 
is supposed to be
 
Wxyw xy
Fxyf xy
i
i
i
i
i
i
(,,,)( ,,),
(,,,)( ,,),
 
 








1
1
   
  
x
i
i
i
y
i
i
xy xy
xy xy
(,,,)( ,,),
(,,,)( ,,)






 1
1
2
i i
q
i
i
i
xy xy





1
,
(,,,)( ,,).   

0
 (27)
In Eq. (27), the time parameter τ

 ( τ

= ετ) is used 
to improve the perturbation procedure. Substituting 
Eq. (27) into Eqs. (25) and (26), then solving the 
perturbation equations step by step, the displacements 
W, ψ
x
, ψ
y
 and stress function F can be obtained. The 
dimensionless transverse load λ
q
 can be derived as
  

q
xy gw gw mx ny
wg
(,,) [( )( )]sins in
(( )) (


41 14 31
1
2
441

c cosc os )
[( )] (sin sin) () .
22
442
42 1
33 3
mx gn y
gw mx ny o

   (28)
In Eq. (28), it should be noted that τ

 is replaced 
by τ. Multiplying Eq. (28) by (sin mx  sin ny) and 
integrating over the plate area, the following nonlinear 
ordinary differential equation can be obtained:
g
dw
d
gwgw gw q
1
2
1
2 21 31
2
1
3
()
() () () () ,


    
4

(29)
in which
   


 

q
q
xy mx nydxdy () (,,) sins in . 
 
4
2
0 0
 (30)

Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 571-582
576 Huang, X.L. – Wang, C.Z. – Wang, J.H. – Wei, N.G.
Table 1.  Comparison of dimensionless fundamental frequencies Ω

 for GPL reinforced plates resting on Winkler-Pasternak elastic foundations
(k
1
, k
2
) Method P – 1, G – 1 P – 2, G – 1 P – 1, G – 3 P – 2, G – 3
(100, 100)
Ref. [13] 0.0637 0.0591 0.0745 0.0696
present 0.0652 0.0603 0.0761 0.0710
Error [%] 2.3 2.0 2.1 2.0
The nonlinear ordinary equation, Eq. (29) can be 
solved by using the Runge-Kutta iteration Scheme 
[23]. For the free vibration problem (λ
q
(τ) = 0), the 
approximate nonlinear frequency can be derived as
 
NL
g
g
gg g
gg
A 
 





41
43
42 41 44
2
41 43
2
12
91 0
12
/
, (31)
where AW h =
max
/ is the dimensionless vibration 
amplitude. If A=0, the dimensionless natural frequency 
is 
L
gg 
41 43
/ .
3  RESULTS AND DISCUSSION
In the section, the several dimensionless parameters 
are used as follows:
 
k
Ka
D
k
Ka
D
k
Kah
D
D
Eh a
h
w
m
w
mm
m
m
m
L m
1
4
2
2
3
3
42
3
2
2
12 1




,, ,
()
,

 

E E
h
E
m
L
m
m
,. 



3.1  Comparison Studies
To verify the accuracy and effectiveness of the present 
method, two numerical examples are presented in this 
subsection.
Example 1. The dimensionless fundamental 
frequencies of GPL-reinforced porous plates resting 
on elastic foundations are calculated and listed in 
Table 1. The material properties and dimensions 
of GPL are E
GPL
 = 1.01 TPa, ρ
GPL
 = 1062.5 kg/m
3
,  
ν
GPL
 = 0.186, a
GPL
 = 2.5 μm, b
GPL
 = 1.5 μm, h
GPL
 = 
1.5 nm. The material properties of the Matrix are E
m
 = 
200 GPa, ρ
m
 = 8908 kg/m
3
, ν
m
 = 0.31.
The geometrical parameters of the plate h = 
0.05 m, a = b = 1.0 m. The GPL weight fraction and 
porosity coefficient are W
GPL
 = 5 %, e
0
 = 0.4. It can 
be observed that the present results are close to those 
given by Gao et al. [13]. The maximum error is about 
2.3 %. This is because Gao et al. [13] employed the 
classic plate theory and differential quadrature method 
to calculate the fundamental frequencies, which is 
different from the present method.
Example 2. The dynamic response of a FG 
GPL reinforced plate is investigated in this example. 
The dispersion pattern of GPL is G – 2. The plate is 
subjected to the explore load:
 qXYt
Pt tt t
tt t
mp p
p
(, ,)
(/ ),
,
, 
 






10
00 and 
 (32)
where the peak pulse p
m
 is 500 kPa and the loading 
time t
p
 is 0.01 s. The material properties of the matrix 
are E
m
 = 3.0 GPa, ρ
m
 = 1290 kg/m³ and ν
m
 = 0.34. 
The corresponding parameters of GPLs are E
GPL
= 
1.01 TPa, ρ
GPL
 = 1062.5 kg/m³ and ν
GPL
 = 0.186. The 
geometric parameters of GPLs and the plate are a = 
b = 0.45, h = 0.045 m, a
GPL
 = 2.5 μm, b
GPL
 = 1.5 μm 
and h
GPL
 = 1.5 nm. The curves of central transient 
deflection versus time are depicted in Fig. 4. It is 
found that the present results are agreement with those 
given by Song et al. [4].
0.0 0.4 0.8 1.2 1.6 2.0 2.4 2.8 3.2
-0.02
0.00
0.02
0.04
0.06
Present
Song et al.[4]
t [ms]
W/h
Fig. 4.  Comparison of the central transient deflection  
for an FG GPL-reinforced plate
3.2 Parametric Studies
In what follows, the effects of material properties 
and foundation parameters are investigated in detail. 
The material and geometric parameters are the same 
as those in Example 1. Unless specially stated, the 
weight fraction W
GPL
 and pore coefficient e
0
 are 5 % 
and 0.2. The dimensionless foundation parameters (k
1
, 
k
2
, k
3
) are (50, 50, 50). 

Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 571-582
577 Nonlinear Vibration Analysis of Functionally Graded Porous Plates Reinforced by Graphene Platelets on Nonlinear Elastic Foundations  
Tables 2 to 4 list the fundamental natural 
frequencies of the plate with different GPL dispersion 
pattern, porosity distribution, weight fractions W
GPL
, 
porosity coefficients e
0
 and foundation parameters 
(k
1
, k
2
). It can be seen that the natural frequency is 
increased with the rising weight fraction W
GPL
 and 
foundation parameters (k
1
, k
2
). The frequency for 
G – 1 is higher than those for G – 2 and G – 3. This 
illustrates that G – 1 can strengthen the plate more 
effectively than the other two patterns. Also, it can 
be seen that the frequency for P – 1 is lower than 
those for P – 2 and P – 3. The fact demonstrated that 
P – 1 can weaken the plate more seriously. In past 
studies [11] to [13], a conclusion was drawn that the 
natural frequency monotonously decreases with 
the increasing porosity coefficient. However, in the 
present case, the conclusion was correct only on P – 2. 
On P – 1 and P – 3, the effect of the porosity coefficient 
Table 2.  Dimensionless fundamental frequencies Ω for porosity distribution P – 1
GPL e
0
(k
1
, k
2
) =(0, 0) (k
1
, k
2
) =(50, 0) (k
1
, k
2
) =(50, 50)
[wt%] [wt%] [wt%]
2.0 5.0 8.0 2.0 5.0 8.0 2.0 5.0 8.0
G – 1
0.0 9.026 12.234 14.731 9.313 12.488 14.977 13.813 16.739 19.202
0.2 8.970 12.149 14.623 9.280 12.422 14.889 14.055 16.952 19.400
0.4 8.928 12.076 14.530 9.265 12.376 14.820 14.378 17.247 19.683
G – 2
0.0 7.359 9.072 10.500 7.709 9.411 10.842 12.783 14.580 16.172
0.2 7.311 8.997 10.404 7.687 9.363 10.773 13.056 14.848 16.440
0.4 7.283 8.946 10.337 7.672 9.345 10.739 13.415 15.212 16.812
G – 3
0.0 8.014 10.364 12.260 8.336 10.663 12.554 13.172 15.420 17.371
0.2 7.982 10.323 12.210 8.328 10.643 12.526 13.444 15.690 17.643
0.4 7.973 10.310 12.198 8.348 10.559 12.440 13.803 16.056 18.021
Table 3.  Dimensionless fundamental frequencies Ω for porosity distribution P – 2
GPL e
0
(k
1
, k
2
) =(0, 0) (k
1
, k
2
) =(50, 0) (k
1
, k
2
) =(50, 50)
[wt%] [wt%] [wt%]
2.0 5.0 8.0 2.0 5.0 8.0 2.0 5.0 8.0
G – 1
0.0 9.026 12.234 14.731 9.313 12.488 14.977 13.813 16.739 19.202
0.2 8.283 11.328 13.526 8.617 11.534 13.822 13.629 16.313 18.595
0.4 7.397 10.035 12.088 7.800 10.394 12.435 13.484 15.889 17.959
G – 2
0.0 7.359 9.072 10.500 7.709 9.411 10.842 12.783 14.580 16.172
0.2 6.775 8.382 9.719 7.180 8.774 10.112 12.765 14.486 16.017
0.4 6.110 7.606 8.845 6.593 8.072 9.312 12.720 14.468 15.942
G – 3
0.0 8.014 10.364 12.260 8.336 10.663 12.554 13.172 15.420 17.371
0.2 7.337 9.488 11.223 7.712 9.836 11.566 13.072 15.155 16.977
0.4 6.553 8.473 10.024 7.004 8.894 10.438 13.037 14.945 16.630
Table 4.  Dimensionless fundamental frequencies Ω for porosity distribution P – 3
GPL e
0
(k
1
, k
2
) =(0, 0) (k
1
, k
2
) =(50, 0) (k
1
, k
2
) =(50, 50)
[wt%] [wt%] [wt%]
2.0 5.0 8.0 2.0 5.0 8.0 2.0 5.0 8.0
G – 1
0.0 9.026 12.234 14.731 9.313 12.488 14.977 13.813 16.739 19.202
0.2 9.026 12.234 14.731 9.313 12.488 14.977 13.813 16.739 19.202
0.4 9.026 12.234 14.731 9.313 12.488 14.977 13.813 16.739 19.202
G – 2
0.0 9.026 12.234 14.731 9.313 12.488 14.977 13.813 16.739 19.202
0.2 9.026 12.234 14.731 9.313 12.488 14.977 13.813 16.739 19.202
0.4 9.026 12.234 14.731 9.313 12.488 14.977 13.813 16.739 19.202
G – 3
0.0 9.026 12.234 14.731 9.313 12.488 14.977 13.813 16.739 19.202
0.2 8.721 11.821 14.233 9.039 12.103 14.506 13.898 16.719 19.317
0.4 8.372 11.348 13.664 8.731 11.666 13.972 14.041 16.747 19.354

Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 571-582
578 Huang, X.L. – Wang, C.Z. – Wang, J.H. – Wei, N.G.
a) 
0.0 0.2 0.4 0.6 0.8 1.0
1.00
1.05
1.10
1.15
1.20
1.25
G-1
G-2
G-3
ω
NL
/ω
L
W
max
  /h
b) 
0.0 0.2 0.4 0.6 0.8 1.0
1.00
1.05
1.10
1.15
1.20
1.25
G-1
G-2
G-3
ω
NL
/ω
L
W
max
  /h
Fig. 5. Effect of GPL dispersion pattern on the frequency ratio;  
a) P – 1,and b) P – 2
a) 
0.0 0.2 0.4 0.6 0.8 1.0
1.00
1.05
1.10
1.15
1.20
P-1
P-2
P-3
ω
NL
/ω
L
W
max
  /h
b) 
0.0 0.2 0.4 0.6 0.8 1.0
1.00
1.05
1.10
1.15
1.20
1.25
W
max
  /h
P-1
P-2
P-3
ω
NL
/ω
L
Fig. 6.  Effect of porosity distribution on the frequency ratio,  
a) G – 1, and b) G – 2
a) 
0.0 0.2 0.4 0.6 0.8 1.0
1.00
1.05
1.10
1.15
1.20
W
GPL
= 0
W
GPL
= 2%
W
GPL
= 5%
ω
NL
/ω
L
W
max
  /h
 b)
0.0 0.2 0.4 0.6 0.8 1.0
1.00
1.05
1.10
1.15
1.20
1.25
W
GPL
= 0
W
GPL
= 2%
W
GPL
= 5%
ω
NL
/ω
L
W
max
  /h
Fig. 7.  Effect of GPL weight fractions on the frequency ratio;  
a) G – 1/P-3, and b) G – 2/P – 3
a) 
0.0 0.2 0.4 0.6 0.8 1.0
1.00
1.05
1.10
1.15
1.20
1.25
e
0
= 0
e
0
= 0.2
e
0
= 0.4
ω
NL
/ω
L
W
max
  /h
b) 
0.0 0.2 0.4 0.6 0.8 1.0
1.00
1.05
1.10
1.15
1.20
1.25
e
0
= 0
e
0
= 0.2
e
0
= 0.4
ω
NL
/ω
L
W
max
  /h
Fig. 8.  Effect of porosity coefficient on the frequency ratio;  
a) G – 3/P – 1, and b) G – 3/P – 2

Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 571-582
579 Nonlinear Vibration Analysis of Functionally Graded Porous Plates Reinforced by Graphene Platelets on Nonlinear Elastic Foundations  
is relevant to the values of foundation parameters (k
1
, 
k
2
). If foundation parameters are (50, 0) or (0, 0), the 
natural frequency is decreased. In contrast, the natural 
frequency for (50, 50) is increased.
Figs. 5 to 8 reveal the influences of GPL 
dispersion pattern, porosity distribution, GPL weight 
fraction W
GPL
, and pore coefficient e
0
 on the nonlinear 
to linear frequency ratio ω
NL
 / ω
L
. As can be observed, 
the frequency ratio for G – 2 / P – 2 is higher than those 
for other GPL and pore distributions. The frequency 
ratio is decreased as pore coefficient e
0
 rises but rose 
with the increase of GPL fraction W
GPL
. 
0.0 0.2 0.4 0.6 0.8 1.0
1.0
1.1
1.2
1.3
1.4
1.5
(k
1
,k
2
,k
3
) = (0,0,0)
(k
1
,k
2
,k
3
) = (50,0,0)
(k
1
,k
2
,k
3
) = (50,50,0)
(k
1
,k
2
,k
3
) = (50,50,50)
ω
NL
/ω
L
W
max
  /h
 
Fig. 9.  Effect of elastic foundation on the frequency ratio  
for G – 3/P – 3
0.0 0.2 0.4 0.6 0.8 1.0
1.00
1.05
1.10
1.15
1.20
1.25
k
3
= 0
k
3
= 50
k
3
= 100
ω
NL
/ω
L
W
max
  /h
Fig. 10.  Effect foundation parameter k3 on the frequency ratio  
for G – 3/P – 3
The effect of nonlinear elastic foundation 
parameters on frequency ratio is shown in Figs. 9 and 
10. The two figures demonstrated that the frequency 
ratio reduces with the increasing   parameters k
1
 and 
k
2
. However, the ratio rises as the nonlinear parameter 
k
3
 increases.
The curves of central transient deflections and 
bending moments for different distributions of GPLs 
and pores are depicted in Figs. 11 and 12. It is found 
a)
0 5 10 15
-3
0
3
6
G-1
G-2
G-3
W [mm]
t [ms]
b)
0 5 10 15
-40
-20
0
20
40
60
80
G-1
G-2
G-3
M
x
[kN.m/m]
t [ms]
Fig. 11.  Effect of GPL dispersion patterns on the central transient 
responses for P – 3; a) dynamic deflection, and  
b) dynamic bending moment
 
a)  
0 5 10 15
-2.5
0.0
2.5
5.0
P-1
P-2
P-3
t [ms]
W [mm]
b)
0 5 10 15
-40
0
40
80
t [ms]
M
x
[kN.m/m]
P-1
P-2
P-3
Fig. 12.  Effect of porosity distributions on the transient responses 
for G – 3; a) dynamic deflection, and b) dynamic bending moment

Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 571-582
580 Huang, X.L. – Wang, C.Z. – Wang, J.H. – Wei, N.G.
that the maximum dynamic deflection for G – 2 / P – 2 
is the largest among all patterns and distributions. 
In contrast, the amplitude of the dynamic bending 
moment for G – 2 / P – 2 is the smallest.
The effects of GPL weight fraction W
GPL
 and 
pore coefficient e
0
 on the transient responses is 
illustrated in Figs. 13 and 14. It is discerned that 
the rise of W
GPL
 reduces the amplitude of transient 
deflection, but increases the amplitude of bending 
moment. The maximum deflection increases by about 
8 % as the porosity coefficient e
0
 rises from 0.0 to 0.4. 
Therefore, a conclusion may be made that the effect of 
the porosity coefficient on the dynamic response can 
be negligible.
a)
0.000 0.005 0.010 0.015
-3
0
3
6
W
GPL
= 0
W
GPL
= 2%
W
GPL
= 5%
W [mm]
t [ms]
b)
0 5 10 15
-40
0
40
80
W
GPL
= 0
W
GPL
= 2%
W
GPL
= 5%
M
x
[kN.m/m]
t [ms]
Fig. 13.  Effect of GPL weight fraction on the transient responses 
for P – 1 / G – 1; a) dynamic deflection, and b) dynamic bending 
moment
The effect of foundation parameters (k
1
, k
2
, k
3
) 
on dynamic responses is presented in Fig. 15. As 
expected, Winkler and Pasternak elastic foundation 
parameters reduce the dynamic responses. The 
dynamic response for Pasternak elastic foundations 
(50, 50, 0) is very close to those for the nonlinear 
elastic foundations (50, 50, 50). Hence, a conclusion 
may be drawn that the effect of nonlinear foundation 
parameter k
3
 on dynamic responses may be neglected, 
which is different from the statement mentioned by 
a) 
0 5 10 15
-2
0
2
4
e
0
= 0
e
0
= 0.2
e
0
= 0.4
W [mm]
t [ms]
b)
0 5 10 15
-40
0
40
80
e
0
= 0
e
0
= 0.2
e
0
= 0.4
M
x
[kN.m/m]
t [ms]
Fig. 14.  Effect of porosity coefficient on the transient responses for 
P – 1 / G – 1; a) dynamic deflection, and b) dynamic bending moment
a)
0 5 10 15
-4
0
4
8
(k
1
, k
2
, k
3
) = (0,0,0)
(k
1
, k
2
, k
3
) = (50,0,0)
(k
1
, k
2
, k
3
) = (50,50,0)
(k
1
, k
2
, k
3
) = (50,50,50)
W [mm]
t [ms]
b)
0 5 10 15
-100
-50
0
50
100
150
200
(k
1
, k
2
, k
3
) = (0,0,0)
(k
1
, k
2
, k
3
) = (50,0,0)
(k
1
, k
2
, k
3
) = (50,50,0)
(k
1
, k
2
, k
3
) = (50,50,50)
M
x
[kN.m/m]
t [ms]
Fig. 15.  Effect of elastic foundation on the transient responses for 
P – 1 / G – 1; a) dynamic deflection, and b) dynamic bending moment

Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 571-582
581 Nonlinear Vibration Analysis of Functionally Graded Porous Plates Reinforced by Graphene Platelets on Nonlinear Elastic Foundations  
Civalek [18] and Najafi et al. [19]. They deemed that 
the nonlinear foundation parameter k
3 
had a significant 
effect on the dynamic responses of laminated and 
FGM plates.
4  CONCLUSIONS
The present work presents a reliable and effective 
method to investigate the nonlinear free and forced 
vibrations of functionally graded GPL reinforced 
porous plates on nonlinear elastic foundations. The 
effects of internal pores, GPLs, and nonlinear elastic 
foundations are discussed in detail. Some interesting 
conclusions can be drawn from the numerical results:
1. Both GPL dispersion patterns and porosity 
distribution can affect the nonlinear vibrations 
and responses for porous plates. Furthermore, 
the effect of a GPL dispersion pattern is more 
significant than that of a porosity distribution. 
2. The GPL weight fraction and foundation 
parameters k
1
 and k
2
 increase the natural 
frequency but decrease the nonlinear to linear 
frequency ratio and transient deflection. 
3. The increase of porosity coefficient does not 
always lead to the rise of natural frequency and 
transient responses.
4. Nonlinear foundation parameters have 
insignificant effects on the nonlinear to linear 
frequency ratio and transient response.
5  ACKNOWLEDGEMENTS
The authors are thankful for the financial support of 
the National Natural Science Foundation of China 
[12162010] and the Natural Science Foundation of 
Guangxi [2021GXNSFAA 220087].
6  REFERENCES
[1] Yusufoglu, E., Avey, M. (2021). Nonlinear dynamic behavior 
of hyperbolic paraboloidal shells reinforced by carbon 
nanotubes with various distributions. Journal of Applied and 
Computational Mechanics, vol. 7, p. 913-921, DOI:10.22055/
JACM.2021.36043.2783.
[2] Avey, M., Yusufoglu, E. (2020). On the solution of large-
amplitude vibration of carbon nanotube-based double-
curved shallow shells. Mathematical Methods in the Applied 
Sciences, vol. 8, p. 1-13, DOI:10.1002/mma.6820.
[3] Zhao, S.Y., Zhao, Z., Yang, Z.C., Ke, L.L., Kitipornchai, S., Yang, 
Y. (2020). Functionally graded graphene reinforced composite 
structures: A review. Engineering Structures, vol. 210, art. ID 
110339, DOI:10.1016/j.engstruct.2020.110339.
[4] Song, M., Kittipornchai S., Yang, J., (2017). Free and forced 
vibrations of functionally graded polymer composite plates 
reinforced with graphene nanoplatelets. Composite Structures, 
vol. 159, p. 579-588, DOI:10.1016/j.compstruct.2016.09.070.
[5] Zhao, Z., Feng, C., Wang, Y., Yang, J. (2017). Bending and 
vibration of functionally graded trapezoidal nanocomposite 
plates reinforced with graphene nanoplatelets (GPLs). 
Composite Structures, vol. 180, p. 799-808, DOI:10.1016/j.
compstruct.2017.08.044.
[6] Gholami R., Ansari R. (2018). Nonlinear harmonically 
excited vibration of third-order shear deformable functionally 
graded graphene platelet- reinforced composite rectangular 
plates. Engineering and Structures, vol. 156, p. 197-209, 
DOI:10.1016/j.engstruct.2017.11.019.
[7] Wu, H., Yang, J., Kittipornchai, S. (2018). Parametric instability 
of thermo-mechanically loaded functionally graded graphene 
reinforced nanaocomposite plates. International Journal of 
Mechanical Science, vol. 135, p. 431-440, DOI:10.1016/j.
ijmecsci.2017.11.039.
[8] Bartolucci, S.F., Paras, J., Rafiee, M.A., Rafiee, J., Lee, S., 
Kapoor, D. (2011). Graphene-aluminum Nanocomposites. 
Material Science and Engineering: A, vol. 528, no. 27, p. 
7933-7937, DOI:10.1016/j.msea.2011.07.043.
[9] Rashad, M., Pan, F., Tang, A., Asif, M. (2014). Effect of 
graphene nanoplatelets addition on mechanical properties 
of pure aluminum using a semi-powder method. Progress of 
Natural Science: Materials International, vol. 24, p. 101-108, 
DOI:10.1016/j.pnsc.2014.03.012.
[10] Wu, H.L, Yang, J., Kitipornchai, S. (2020). Mechanical analysis 
of functionally graded porous structures: A review. International 
Journal of Structural Stability and Dynamics, vol. 20, no. 13, 
art. ID 2041015, DOI:10.1142/S0219455420410151.
[11] Kittipornchai, S., Chen, D., Yang, J. (2017). Free vibration 
and elastic buckling of functionally graded porous beams 
reinforced with graphene platelets. Material Design, vol. 116, 
p. 656-665, DOI:10.1016/j.matdes.2016.12.061.
[12] Yang, J., Chen, D., Kittipornchai, S. (2018). Buckling and free 
vibration analysis of functionally graded graphene reinforced 
porous nanocomposite plates on Chebyshev-Ritz method. 
Composite Structures, vol. 193, p. 281-294, DOI:10.1016/j.
compstruct.2018.03.090.
[13] Gao, K., Gao, W., Chen, D., Yang, J. (2018). Nonlinear 
vibration of functionally graded graphene platelets reinforced 
porous nanocomposite plates resting on elastic foundation. 
Composite Structures, vol. 204, p. 831-846, DOI:10.1016/j.
compstruct.2018.08.013.
[14] Saidi, A.R., Bahaadini, R., Majidi-Mozafari, K. (2019). 
On vibration and stability of porous plates reinforced by 
graphene platelets under aerodynamic loading. Composite 
Part B: Engineering, vol. 163, p. 778-799, DOI:10.1016/j.
compositesb.2019.01.074.
[15] Duc, N.D., Lee, J., Nguyen-Thoi, T., Thang, P.T. (2017). Static 
response and vibration of functionally graded carbon 
nanotube-reinforced composite rectangular plates resting 
on Winkler-Pasternak elastic foundations. Aerospace 
Science and Technology, vol. 68, p. 391-402, DOI:10.1016/j.
ast.2017.05.032.
[16] Haciyev, V.C., Sofiyev, A.H., Kuruoglu, N. (2018). Free bending 
vibration analysis of thin bidirectionally exponentially graded 
orthortropic rectangular plates resting on two-parameter 

Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 571-582
582 Huang, X.L. – Wang, C.Z. – Wang, J.H. – Wei, N.G.
elastic foundations. Composite Structures, vol. 184, p. 372-
377, DOI:10.1016/j.compstruct.2017.10.014.
[17] Nath, Y., Prithviraju, M., Mufti, A.A. (2006). Nonlinear 
statics and dynamics of antisymmetric composite 
laminated square plates supported on nonlinear elastic 
subgrade. Communication in Nonlinear Science and 
Numerical Simulation, vol. 11, p. 340-354, DOI:10.1016/j.
cnsns.2004.11.003.
[18] Civalek, Ö. (2013). Nonlinear dynamic response of laminated 
plates resting on nonlinear elastic foundations by the 
discrete singular convolution-differential quadrature coupled 
approaches. Composite Part B: Engineering, vol. 50, p. 171-
179, DOI:10.1016/j.compositesb.2013.01.027.
[19] Najafi, F., Shojaeefard, M.H., Googarchin, H.S. (2016). Low-
velocity impact response of functionally graded plates with 
nonlinear elastic foundation in thermal fields. Composite 
Part B: Engineering, vol. 10, p. 123-140, DOI:10.1016/j.
compositesb.2016.09.070.
[20] Reddy, J.N. (1984). A refined nonlinear theory of plates with 
transverse shear deformation. International Journal of Solids 
and Structures, vol. 20, p. 881-896, DOI:10.1016/0020-
7683(84)90056-8.
[21] Shen, H. (1997). Kármán-type equations for a higher-order 
shear deformation plate theory and its use in the thermal 
postbuckling analysis. Applied Mathematics and Mechanics, 
vol. 18, p. 1137-1152, DOI:10.1007/BF00713716.
[22] Huang, X.-L., Zheng, J.-J. (2003). Nonlinear vibration and 
dynamic response of simply supported shear deformable 
laminated plates on elastic foundations. Engineering 
Structures, vol. 25, p. 1107-1119, DOI:10.1016/S0141-
0296(03)00064-6.
[23] Pearson, C.E. (1986). Numerical Methods in Engineering & 
Science. Van Nostrand Reinhold Company, New York.