Accepted Manuscript

Application of sinusoidal shear deformation theory and physical neutral surface to

analysis of functionally graded piezoelectric plate

Mohammad Arefi, Elyas Mohammad-Rezaei Bidgoli, Rossana Dimitri, Michele

Bacciocchi, Francesco Tornabene

PII: S1359-8368(18)31269-1
DOI: 10.1016/[Link].2018.05.050
Reference: JCOMB 5720

To appear in: Composites Part B

Received Date: 24 April 2018

Revised Date: 24 May 2018
Accepted Date: 29 May 2018

Please cite this article as: Arefi M, Mohammad-Rezaei Bidgoli E, Dimitri R, Bacciocchi M, Tornabene F,
Application of sinusoidal shear deformation theory and physical neutral surface to analysis of functionally
graded piezoelectric plate, Composites Part B (2018), doi: 10.1016/[Link].2018.05.050.

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to
our customers we are providing this early version of the manuscript. The manuscript will undergo
copyediting, typesetting, and review of the resulting proof before it is published in its final form. 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.
 ACCEPTED MANUSCRIPT

APPLICATION OF SINUSOIDAL SHEAR DEFORMATION THEORY

AND PHYSICAL NEUTRAL SURFACE TO ANALYSIS OF
FUNCTIONALLY GRADED PIEZOELECTRIC PLATE

PT
Mohammad Arefi1, Elyas Mohammad-Rezaei Bidgoli1, Rossana Dimitri2, Michele Bacciocchi3,

Francesco Tornabene3*

RI
ABSTRACT. The concept of neutral surface for a Functionally Graded Piezoelectric (FGP)

SC
plate is developed in this paper. The electro-elastic analysis of a FGP plate resting on Winkler-
Pasternak foundation is performed in the theoretical framework provided by a two-variable
sinusoidal shear deformation theory, including the aforementioned concept of neutral surface.

U
First, the location of neutral surface is defined with respect to the position of the middle surface
AN
and then the governing bending equations are derived using the principle of virtual work. An
analytical method is presented to investigate the influence on the displacement and stress
components of the main parameters of the model, which are volume fraction exponent of the
M

constituents, the electric potential and foundation parameters. The numerical results are validated
through the comparison with available references. The numerical results prove that the value of
D

transverse bending deflection is greater than the corresponding shear deflection. In addition, it
TE

can be observed that the volume fraction exponent has peculiar influence on the distribution of
displacements and stresses.
EP

KEYWORDS: Neutral Surface; Functionally Graded Piezoelectric materials; Two-variable

C

sinusoidal shear deformation plate theory; Electric potential; Pasternak elastic foundation.
AC

1
Department of Solid Mechanics, Faculty of Mechanical Engineering, University of Kashan, Kashan 87317-51167, I.R. Iran
2
Department of Innovation Engineering, Università del Salento, 73100 Lecce, Italy
3
DICAM - Department, School of Engineering and Architecture, University of Bologna, Italy. [Link]@[Link];
web page: [Link]

1
 ACCEPTED MANUSCRIPT

1. INTRODUCTION

It is well-known that Functionally Graded Materials (FGMs) are characterized by a smooth

variation of the mechanical properties within a certain medium [1-6]. They were introduced in
the late 1980s by a Japanese group during a space plane project. In general, these advanced

PT
components are made of a combination of different materials, typically a mixture of ceramic and
metal constituents [7-12]. Their increasing use is due to their ability to bear high temperature
gradients without losing the structural integrity [13-15]. In addition, the continuous variation of

RI
the mechanical properties reduces the stress concentrations at the interface of different layers in
laminated composite structures. As proven by many papers in the literature, FGMs still represent

SC
an interesting topic which draws the attention of many researchers. It should be mentioned that
the same idea of smooth variation is employed nowadays to characterize the mechanical

U
properties of nanocomposites [16-19].
AN
Since FGMs are characterized by variable mechanical properties from a surface to another one,
the location of the neutral surface can differ from the middle surface of the structure, if the
variation is applied along the thickness direction. This difference could lead to an inaccurate
M

analysis of FGM structures, in which the thickness coordinate is measured from their middle
surface. A brief literature review related to the concept of neutral surface and FGMs is presented
D

below for the sake of completeness. Reddy et al. [20] presented the bending analysis of
TE

Functionally Graded (FG) solid and annular circular plates based on First-order Shear
Deformation Theory (FSDT). They studied the influence of various percentages of ceramic-
metal volume fractions on the stress and displacement distributions. Reddy [21] employed the
EP

Third-order Shear Deformation Theory (TSDT) and nonlinear Von-Karman relations for the
theoretical and analytical solution of a FG plate subjected to thermo-mechanical loads. Reddy et
C

al. [22] presented a three-dimensional thermo-elastic analysis of a FG rectangular plate based on

AC

an asymptotic method. The effective material properties were estimated using the Mori-Tanaka
scheme. Reddy’s higher order shear deformation plate theory was employed by Yang and Shen
for the nonlinear bending analysis of a FG plate subjected to thermal and mechanical loads [23].
In particular, temperature dependency was considered for evaluating the material properties and
various numerical approaches, such as perturbation technique, one-dimensional differential
quadrature approximation and Galerkin method, were employed to obtain the solution of the
problem. The FSDT was used by Lanhe [24] to study the buckling behavior of a FG plate due to

2
 ACCEPTED MANUSCRIPT

thermal loads. The influence of various material and geometric parameters was investigated. Kim
[25] used the TSDT for the vibration analysis of a FG plate using the Rayleigh-Ritz method.
Qian et al. [26] presented the static, free and forced vibration analyses of FG plates based on
Higher-order Shear Deformation Theories (HSDTs) using the meshless local Petrov-Galerkin
method. The free vibration, buckling, and static analyses of a FG plate were studied also by

PT
Abrate [27]. A generalized shear deformation theory was employed by Zenkour [28] for the
static analysis of a simply-supported FG plate subjected to a transverse uniform load. The FG

RI
plate was graded along the thickness direction by means of a simple power-law distribution.
Zenkour [29] investigated also the influence of relevant parameters, such as aspect ratio, side-to-

SC
thickness ratio and the exponentially graded parameter, on the strain and stress distributions of
the FG plate. Chi and Chung [30] used the Classical Plate Theory (CPT) for the elastic analysis

U
of a FG plate subjected to transverse loads using a Fourier-based approach. Prakash et al. [31]
employed the Finite Element Method (FEM) to study the nonlinear behavior of FG skew plates.
AN
In particular, the influence of neutral surface, boundary conditions and geometric parameters
were analyzed. Zhang and Zhou [32] used the concept of physical neutral surface for the analysis
M

of a FG thin plate. They assumed variable mechanical properties along the thickness direction, as
far as the Young’s modulus and the mass density are concerned. Singha et al. [33] studied the
D

nonlinear behavior of FG plates through the FEM considering the concept of physical neutral
surface in the theoretical framework provided by the FSDT. Alinia and Ghannadpour [34]
TE

investigated the nonlinear behavior of FG plates subjected to external loads including Von-
Karman nonlinearities. Their results were obtained by minimizing the total potential energy with
EP

respect to the unknown displacements. Kulikov and Plotnikova [35] presented a three-
dimensional thermo-electro-elastic analysis of FG piezoelectric laminated plates in thermal
C

environment. Jandaghian and Rahmani [36] investigated the thermo-electro-mechanical vibration

problem of FGP rectangular nano-plates subjected to applied electric voltage, biaxial forces and
AC

uniform temperature increase, by means of the Kirchhoff plate theory. Matsunaga [37] employed
a HSDT for the two-dimensional analysis of a FG plate subjected to thermal and mechanical
loads. The free vibration analysis of an annular FG plate integrated with two piezoelectric layers
placed at the top and the bottom of the structure was performed by Ebrahimi and Rastgoo [38] by
using the Kirchhoff plate theory. The free and forced vibrations of a FGP rectangular plate with
exponent-law dependent material properties were investigated by Zhong and Yu [39] using a

3
 ACCEPTED MANUSCRIPT

Fourier-based approach. In the papers [40-42], the results of the free vibration and bending
analyses of FG rectangular and circular plates resting on elastic foundation based on new refined
plate theory were presented in terms of mechanical indices and geometric parameters. Arefi and
Zenkour [43-47] presented some interesting results related to the electro-elastic analysis of FG
structures. Further details can be found also in the works [48-56].

PT
To the best of authors’ knowledge, there is no systematic published work concerning the effect
of neutral surface on the electro-elastic analysis of FGP plate. For the sake of completeness, it

RI
should be recalled that a FGP plate is made of two piezoelectric materials, whose electrical and
mechanical properties are graded along the thickness direction of the structure. A power-law

SC
distribution is used to this aim. The reader can find different through-the-thickness distributions
for FGMs in the papers [57-61]. In the current work, the FGP plate is subjected to transverse

U
static loads and applied electric potential. A two-variable sinusoidal shear deformation theory is
used to develop the governing equations of the problem, which are solved by means of an
AN
analytical approach. The Navier’s method is employed to this aim.
M

2. FORMULATION
The electro-elastic analysis of a FGP plate resting on Pasternak foundation is presented in this
D

paper. The concept of neutral axis and a sinusoidal shear deformation theory are employed for
this purpose. The FGP plate is subjected to transverse loads. In addition, an electric potential is
TE

applied at the top of the structure. Let us to consider a plate with length a , width b and
thickness h . A Cartesian coordinate system x, y , z is introduced in order to define the various
EP

structural parameters. The schematic representation of this model is depicted in Figure 1. As

illustrated in Figure 1, the model is of two piezoelectric materials at top and bottom of the
C

structure.
The coordinate zms is considered from mid-surface of plate (the subscript “ ms ” stands for
AC

middle surface). On the other hand, the coordinate zns is considered from the neutral surface

(denoted by “ ns ”). The volume fractions Vc1 ( z ) ,Vc2 ( z ) of FGP materials (respectively Ceramic 1

and Ceramic 2 specified in Figure 1) are defined as follows, according to a power-law

distribution [48]

4
 ACCEPTED MANUSCRIPT

1  z +c 1
r r
z
Vc2 ( z ) =  ms +  =  ns +  (1)
 h 2  h 2

Vc1 ( z ) = 1 − Vc2 ( z ) (2)

where c define the position of the neutral surface with respect to the middle surface, whereas r

PT
stands for the exponent of the volume fractions. The mechanical characterization of the
composite is provided in terms of Young’s modulus E ( z ) , which is computed as follows

E ( z ) = Ec1Vc1 ( z ) + Ec2 Vc2 ( z )

RI
(3)

in which Ec1 , Ec2 are the Young’s moduli of the two constituents. It can be noted that the top of the plate

SC
is entirely made of Ceramic 2, whereas at the bottom of the plate the material has the mechanical feature
of the Ceramic 1. As specified in the paper by Singha et al. [33], the following relation can be
used to compute the position of the neutral surface c

∫
h

U
E ( zms )( zms − c ) dzms = 0
AN
2
−h (4)
2

from which one gets

M

∫ 2
−h E ( zms ) zms dzms
c= 2
h
(5)
∫ E ( zms ) dzms
D

2
−h
2
TE

It is obvious that for homogeneous materials, the neutral surface coincides with the middle
surface of the plate. On the other hand, for non-homogeneous materials c must be calculated
since it assumes a different value. In particular, its value increases if the difference between Ec1
EP

and Ec2 is bigger.

Once the position of the neutral surface is defined, it is possible to introduce the displacement
C

field, as well as the constitutive relations. A two-variable sinusoidal shear deformation plate
AC

theory is used in this paper. According to this approach [44], the three displacement components
u ( x, y, zns ) , v ( x, y, zns ) , w( x, y, zns ) are defined as follows

∂wb ( x, y ) ∂ws ( x, y )
u ( x, y, zns ) = u0 ( x, y ) − zns − f ( zns ) (6)
∂x ∂x
∂wb ( x, y ) ∂ws ( x, y )
v ( x, y, zns ) = v0 ( x, y ) − zns − f ( zns ) (7)
∂y ∂y

5
 ACCEPTED MANUSCRIPT

w ( x, y, zns ) = wb ( x, y ) + ws ( x, y ) (8)

where wb , ws are the plate deflections due to bending and shear, respectively. On the other hand,
u0 , v0 are the translational displacements along x, y , respectively. Since a sinusoidal shear

deformation theory is used, the function f ( z ns ) assumes the following aspect

PT
π  π ( zns + c ) 
f ( zns ) = zns + c − sin   (9)
h  h 

RI
Once the displacement field is defined, the kinematic equations can be carried out. In particular,
the three-dimensional strain components are defined as follows, assuming ε zz = 0

SC
ε xx = ε x0 + zns k xb + f ( zns ) k xs (10)

ε yy = ε y0 + zns k yb + f ( zns ) k ys (11)

U
γ xy = γ xy0 + zns k xyb + f ( zns ) k xys (12)
AN
γ yz = g ( zns ) γ yzs (13)

γ xz = g ( zns ) γ xzs (14)

M

in which
df ( zns )
g ( zns ) = 1 − (15)
D

dzns

The following definitions are required as well

TE

∂u 0 ( x , y ) ∂v0 ( x, y )
ε x0 = , ε y0 = (16)
∂x ∂y
EP

∂ 2 wb ( x, y ) ∂ 2 ws ( x, y ) ∂ 2 wb ( x, y )
k xb = − , k xs = − , k xyb = −2 (17)
∂x 2 ∂x 2 ∂x∂y

∂ 2 wb ( x, y ) ∂ 2 ws ( x, y ) ∂ 2 ws ( x, y )
C

k yb = − , k ys = − , k xys = −2 (18)
∂y 2 ∂y 2 ∂x∂y
AC

∂u0 ( x, y ) ∂v0 ( x, y )
γ xy0 = + (19)
∂y ∂x

∂ws ( x, y ) ∂ws ( x, y )
γ yzs = , γ xzs = (20)
∂y ∂x

As far as the constitutive laws are concerned, the following electro-elastic relations allow to
compute the stress components

6
 ACCEPTED MANUSCRIPT

σ xx   c c12 0 0 0  ε xx   0 0 e31 
   11  ε   0
σ yy   c12 c22 0 0 0   yy   0 e32   Ex 
   
τ yz  =  0 0 c44 0 0  γ yz  −  0 e24 0   Ey  (21)
  0  
0  γ xz   e15

0   Ez 
τ xz  
0 0 c55 0
 
τ xy   0 0 0 0 c66  γ xy   0 0 0 
   

PT
where cij are the stiffness coefficients and eij stand for the piezoelectric constants. The electric

field components are denoted by Ei , for i = 1, 2, 3 . Both the stiffness coefficients and

RI
piezoelectric constants depend on the thickness coordinate as shown below
cij = cij ,c1Vc1 ( z ) + cij ,c 2Vc 2 ( z ) (22)

SC
eij = eij ,c1Vc1 ( z ) + eij ,c 2Vc 2 ( z ) (23)

in which cij ,c1 , cij ,c 2 and eij ,c1 , eij ,c 2 are the stiffness coefficients and piezoelectric constants

U
associated with the two constituents of the plate. At this point, the electric potential Φ ( x, y, zns )
AN
must be introduced. It can be assumed as specified in the paper by Arefi and Zenkour [47]
2 znsφ0 πz 
Φ ( x, y, zns ) = − cos  ns  φ ( x, y )
M

(24)
h  h 
where φ0 is the applied voltage, whereas φ ( x, y ) stands for the two-dimensional electric potential
D

distribution. Consequently, the electric field components can be evaluated as follows

∂Φ ( x, y, zns )  π z  ∂φ ( x, y )
TE

Ex = − = cos  ns  (25)
∂x  h  ∂x
∂Φ ( x, y, zns )  π z  ∂φ ( x, y )
EP

Ey = − = cos  ns  (26)
∂y  h  ∂y
∂Φ ( x, y, zns ) π πz  2φ0
Ez = − = − sin  ns  φ ( x, y ) − (27)
C

∂zns h  h  h
AC

The electric displacements for piezoelectric materials are given by

ε xx 
 
 Dx   0 0 0 e15 0  ε yy  κ11 0 0   Ex 
       
 Dy  =  0 0 e24 0 0 γ yz  + 0 κ 22 0   E y 
 (28)
  
  e 0  γ xz   0 0 κ 33   Ez 
 Dz   31 e32 0 0
 
γ xy 
 

7
 ACCEPTED MANUSCRIPT

in which κ ii represent the dielectric coefficients that must be evaluated according to the volume
fraction law.
In this stage, the equations that govern the electro-elastic behavior of the structure can be
accomplished. For this purpose, the principle of virtual works is applied as follows
δΠ = δ U − δ V = 0

PT
(29)
where U is strain energy and V is energy contribution due to external work. The following
quantity defines the variation of the strain energy δ U

RI
h
−c

∫ ∫ (σ δε xx + σ yyδε yy + τ xyδγ xy + τ yzδγ yz + τ xzδγ xz − Dxδ Ex − Dyδ E y − Dzδ Ez ) dAdzns

2
δU = (30)

SC
xx
−h
−c A
2

in which A represents the area of the middle surface. On the other hand, the variation of the

U
energy due to external work δ V includes the contribution of both elastic foundation Ff and
AN
external forces q

δ V = − ∫ (q − Ff )δ wdA = 0 (31)
A
M

The Winkler-Pasternak foundation effect is taken into account as follows

 ∂ 2 w ∂ 2 wb ∂ 2 ws ∂ 2 ws 
Ff = k1 ( wb + ws ) − k2  2b + + 2 + 2 
D

(32)
 ∂x ∂y 2 ∂x ∂y 
TE

where k1 , k2 are the stiffness and shear parameters of the foundation. Having in mind the
definitions of the strain components and the electric field, the principle of virtual works provides
EP

the following result after some manipulations

∫ ( N δε + N y δε y0 + N xy δγ xy0 + M xbδ k xb + M ybδ k yb + M xyb δ k xyb + M xs δ k xs +

0
x x
A
C

Dx Dy 
+ M ys δ k ys + M xys δ k xys + S yzs δγ yzs + S xzs δγ xzs + δf+ δ f + Dz δ f  dA + (33)

AC

x y 

− ∫ ( q − Ff ) (δ w b + δ ws ) = 0
A

where the following stress resultants are introduced

8
 ACCEPTED MANUSCRIPT

 Nx Ny N xy  h2 − c σ xx 
 b   
M x M yb M xyb  = ∫ σ yy  1 zns f ( zns )  dzns (34)
 M xs M s
M xys  − h2 − c  τ xy 
 y

h
−c
 S xzs  2 τ xz 
 s  = ∫   g ( zns ) dzns (35)

PT
 S yz  − h − c τ yz 
2

 Dx  h2 − c  Dx 

RI
      π zns   π zns  π  π zns  
 Dy  = ∫  Dy  cos  h  cos  h  h sin  h   dzns (36)
 Dz  − h2 − c  Dz        

SC
These quantities can be written as a function of the displacement components and the electric
potential as follows, by using the definitions of the stress components and the electric

U
displacements introduced before
∂u0 ∂ 2 wb s ∂ ws ∂v0 ∂ 2 wb s ∂ ws
AN
2 2
N x = A11 − B11 2 − B11 2 + A12 − B12 − B12 + P11φ + P31 (37)
∂x ∂x ∂x ∂y ∂y 2 ∂y 2

∂u0 ∂ 2 wb s ∂ ws
2
∂v0 ∂ 2 wb s ∂ ws
2
N y = A12 − B12 − + − − + P12φ + P32
M

B A B B (38)
∂x ∂x 2 ∂x 2 ∂y ∂y 2 ∂y 2
12 11 11 11

 ∂u ∂v  ∂ 2 wb ∂ 2 ws
N xy = A66  0 + 0  − 2 B66 − 2 B66s
D

(39)
 ∂y ∂x  ∂x∂y ∂x∂y
TE

∂u0 ∂2w ∂ 2 ws ∂v0 ∂ 2 wb s ∂ ws

2
M xb = B11 − D11 2b − D11s + B − D − D + P13φ + P33 (40)
∂x ∂x ∂x 2 ∂y ∂y 2 ∂y 2
12 12 12

∂u0 ∂ 2 wb s ∂ ws ∂v0 ∂ 2 wb s ∂ ws
EP

2 2
M yb = B12 − D12 − D + B − D − D + P14φ + P34 (41)
∂x ∂x 2 ∂x 2 ∂y ∂y 2 ∂y 2
12 11 11 11

 ∂u0 ∂v0  ∂ 2 wb s ∂ ws
2
C

M = B66 
b
+  − 2 D66 − 2 D66 (42)
 ∂y ∂x  ∂x∂y ∂x∂y
xy
AC

∂u0 ∂ 2 wb s ∂ ws
2
s ∂v0 s ∂ wb
2
s ∂ ws
2
M xs = B11s − D11s − H + B − D − H + P15φ + P35 (43)
∂x ∂x 2 ∂x 2 ∂y ∂y 2 ∂y 2
11 12 12 12

∂u0 s ∂ wb
2
s ∂ ws
2
s ∂v0 s ∂ wb
2
s ∂ ws
2
M =B
s
− D12
s
− H12 2 + B11 − D11 2 − H11 2 + P16φ + P36 (44)
∂x ∂x 2 ∂x ∂y ∂y ∂y
y 12

 ∂u ∂v  ∂ 2 wb ∂ 2 ws
M xys = B66s  0 + 0  − 2 D66s − 2 H 66s (45)
 ∂y ∂x  ∂x∂y ∂x∂y

9
 ACCEPTED MANUSCRIPT

∂ws ∂φ
S xzs = A55s − P17 (46)
∂x ∂x
∂ws ∂φ
S yzs = A44s − P18 (47)
∂y ∂y
∂ws ∂φ
Dx = P17 + P20 (48)

PT
∂x ∂x
∂ws ∂φ
Dy = P18 + P22 (49)
∂y ∂y

RI
∂u0 ∂2w ∂2w ∂v ∂ 2 wb ∂ 2 ws
Dz = P11 − P13 2b − P15 2s + P12 0 − P14 − P − P30φ − P29 (50)
∂x ∂x ∂x ∂y ∂y 2 ∂y 2
16

SC
where the electro-elastic properties of the composites are defined in the Appendix. Finally, a
system of five governing equations is accomplished by means of the principle of virtual works

U
∂N x ∂N xy
+ =0 (51)
AN
∂x ∂y

∂N xy ∂N y
+ =0 (52)
∂x ∂y
M

∂ 2 M xb ∂ 2 M xyb ∂ 2 M yb
+ 2 + + q − Ff = 0 (53)
∂x 2 ∂x∂y ∂y 2
D

∂ 2 M xs ∂ 2 M xys ∂ 2 M ys ∂S xzs ∂S yzs

TE

+2 + + + + q − Ff = 0 (54)
∂x 2 ∂x∂y ∂y 2 ∂x ∂y

∂Dx ∂Dy
+ + Dz = 0 (55)
EP

∂x ∂y
The fundamental equations can be obtained by inserting definitions (37)-(50) into the
C

equilibrium equations
∂ 2 u0 ∂ 2u 0 ∂ 2v0 ∂ 3wb ∂ 3 wb
AC

A11 + A + ( A + A ) − B − ( B + 2 B ) +
∂x 2 ∂y 2 ∂x∂y ∂x3 ∂x∂y 2
66 12 66 11 12 66

(56)
∂ 3 ws s ∂ ws
3
∂φ ∂P31
− ( B + 2B )
s s
− B + P11 + =0
∂x∂y ∂x ∂x ∂x
12 66 2 11 3

10
 ACCEPTED MANUSCRIPT

∂ 2v0 ∂ 2v0 ∂ 2 u0 ∂ 3 wb ∂ 3 wb
A11 + A + ( A + A ) − B − ( B + 2 B ) +
∂y 2 ∂x 2 ∂x∂y ∂y 3 ∂x 2∂y
66 12 66 11 12 66

(57)
∂ 3w ∂ 3 ws ∂φ ∂P32
− ( B12s + 2 B66s ) 2 s − B11s + P12 + =0
∂x ∂y ∂y 3
∂y ∂y

∂ 3u0 ∂ 3u0 ∂ 3v0 ∂ 3v0 ∂ 4 wb

+ ( + ) + ( + ) + − +

PT
B11 B 2 B B 2 B B D
∂x3 ∂x∂y 2 ∂x 2∂y ∂y 3 ∂x 4
12 66 12 66 11 11

∂ 4 wb ∂ 4 wb s ∂ ws
4
∂ 4 ws
− 2 ( D12 + 2 D66 )
∂x 2∂y 2
+ − D11
∂y 4
− D11
∂x 4
− 2 ( D s
12 + 2 D s
66 ) ∂x 2∂y 2
+

RI
(58)
∂ 4 ws ∂ 2φ ∂ 2φ ∂ 2 P33 ∂ 2 P34
−D s
+ P + P + + + q − k1 ( wb + ws ) +
∂y 4 ∂x 2 ∂y 2 ∂x 2 ∂y 2
11 13 14

SC
 ∂ 2 w ∂ 2 wb ∂ 2 ws ∂ 2 ws 
+ k2  2b + + + =0
 ∂x ∂y 2 ∂x 2 ∂y 2 

U
∂ 3u0 ∂ 3u0 ∂ 3v0 s ∂ v0
3
s ∂ wb
4
B11s
∂x3
+ ( 12 66 ) ∂x∂y 2 ( 12 66 ) ∂x 2∂y 11 ∂y 3 11 ∂x 4 +
B s
+ 2 B s
+ B s
+ 2 B s
+ B − D
AN
∂ 4 wb s ∂ wb
4
s ∂ ws
4
∂ 4 ws
− 2 ( D12s + 2 D66s ) − D − H − 2 ( 12
H s
+ 2 H 66 )
s
+
∂x 2∂y 2 ∂y 4 ∂x 4 ∂x 2∂y 2
11 11

(59)
∂ 4 ws s ∂ ws s ∂ ws ∂ 2φ ∂ 2φ
M

2 2
−H s
+ A + A + ( P − P ) + ( P − P ) +
∂y 4 ∂x 2 ∂y 2 ∂x 2 ∂y 2
11 55 44 15 17 16 18

∂ 2 P35 ∂ 2 P36  ∂ 2 wb ∂ 2 wb ∂ 2 ws ∂ 2 ws 
( )
D

+ + + q − k w + w + k 2 + + + =0
∂x 2 ∂y 2  ∂x ∂y 2 ∂x 2 ∂y 2 
1 b s 2
TE

∂u0 ∂v0 ∂ 2 wb ∂ 2 wb ∂ 2 ws ∂ 2 ws
− P11 − P12 + P13 2 + P14 + ( P15 − P17 ) 2 + ( P16 − P18 ) 2 +
∂x ∂y ∂x ∂y 2 ∂x ∂y
(60)
∂ 2φ ∂ 2φ
− P20 2 − P22 2 + P30φ + P29 = 0
EP

∂x ∂y
C

3. SOLUTION PROCEDURE
AC

In this section, an analytical method for the solution of the electro-elastic problem is presented.
In particular, a Navier-type solution is carried out as shown in the paper by Zenkour et al. [49].
For this purpose, the FGP plate is assumed as simply-supported along its external edges and it is
subjected to homogeneous electrical boundary conditions. The following definitions are required

11
 ACCEPTED MANUSCRIPT

  mπ x   nπ y  
 u ⋅ cos   ⋅ sin  
  a   b 
  mπ x   nπ y  
 uo   v ⋅ sin   ⋅ cos   
    a   b  
 vo  ∞ ∞  
 wb  = ∑∑  wb ⋅ sin  mπ x  ⋅ sin  nπ y   (61)

PT
  n =1 m =1   a   b 
 ws  
φ   mπ x   nπ y  
   w s ⋅ sin   ⋅ sin  
  a   b 

RI
 mπ x   nπ y  
 φ ⋅ sin   ⋅ sin   
  a   b  

SC
where u , v , wb , ws ,φ are the maximum values of the displacements and electric potential,
respectively. Analogously, the transverse load can be taken in the following form given by the

U
Fourier series
∞ ∞
 mπ x   nπ y 
AN
q ( x, y ) = ∑∑ qmn ⋅ sin   ⋅ sin   (62)
n =1 m =1  a   b 
For the sake of conciseness, the five variables u , v , wb , ws ,φ can be included in the following
M

vector
T
x = u v wb ws φ  (63)
D

which collects the unknowns of the problem in hand. Consequently, the system of governing
TE

equations assumes the aspect shown below

Kx = F (64)
EP

where K is the stiffness matrix whose size is 5 ⋅ 5 and F is the load vector, which is given by
T
 ∂P ∂P ∂ 2 P33 ∂ 2 P34 ∂ 2 P35 ∂ 2 P36 16 
F =  − 31 − 32 − q0 − − −q0 − − P
2 29 
(65)
 ∂x ∂y ∂x 2 ∂y 2 ∂x 2 ∂y 2 π
C


As far as the stiffness operator is concerned, only fifteen elements are required for its complete
AC

characterization due to its symmetric form. In particular, its elements K ij are shown below

A11m 2π 2 A66 n 2π 2
K11 = − − (66)
a2 b2

K12 = −
( A12 + A66 ) mπ 2 n (67)
ab

12
 ACCEPTED MANUSCRIPT

B11m3π 3 ( B12 + 2 B66 ) mπ n

3 2

K13 = + (68)
a3 ab 2

K14 =
(B s
12 + 2 B66s ) mπ 3n 2
+
B11s m3π 3
(69)
ab 2 a3
P11mπ

PT
K15 = (70)
a
A11n 2π 2 A66 m 2π 2
K 22 = − − (71)

RI
b2 a2
B11n3π 3 ( B12 + 2 B66 )m 2π 3n
K 23 = + (72)

SC
b3 a 2b

K 24 =
(B s
12 + 2 B66s ) m 2π 3n
+
B11s n3π 3
(73)
a 2b b3

U
K 25 =
P12 nπ
(74)
AN
b
D11m 4π 4 2 ( D12 + 2 D66 ) m π n D11n 4π 4  m2π 2 n 2π 2 
2 4 2

K 33 = − − − − k1 + k 2− − 2  (75)
a4 a 2b 2 b4  a
2
b 
M

D11s m4π 4 2 ( D12 + 2 D66 ) m π n D11s n 4π 4  m 2π 2 n 2π 2 

s s 2 4 2

K 34 = − − − − k + k 2 − − 2  (76)
D

1
a4 a 2b 2 b4  a
2
b 

P13m 2π 2 P14 n 2π 2
TE

K 35 = − − (77)
a2 b2

 m 4π 4 n 4π 4  2 ( H12 + 2 H 66 ) m π n
s s 2 4 2

K 44 = − H  4 + 4  −
s
+
EP

11
 a b  a 2b 2
(78)
As m 2π 2 A44s n 2π 2  m 2π 2 n 2π 2 
− 55 2 − − k1 + k 2 − − 2 
a b2  a
2
b 
C

P15m 2π 2 P16 n 2π 2 P17 m 2π 2 P18n 2π 2

AC

K 45 = − − + + (79)
a2 b2 a2 b2
P20m 2π 2 P22 n 2π 2
K 55 = + + P30 (80)
a2 b2

13
 ACCEPTED MANUSCRIPT

4. NUMERICAL RESULTS AND DISCUSSION

The numerical results are discussed in the present section. In particular, the structural response of
a FGP plate of thickness h = 0.01m resting on Pasternak foundation ( k1 = 109 N m3 ,

k2 = 107 N m ) due to transverse loads q = 1000 N m 2 and applied electric potential is

PT
investigated. The size of the plate is given by a = b = 0.4 m . For the sake of conciseness, the
mechanical and electrical properties of the constituents are listed in Table 1. These values can be

RI
found also in the papers by Giannakopoulos and Suresh [50] and Tabib-Azar [51]. It should be
specified that the top and the bottom surfaces of the plate have the mechanical properties of the

SC
Lead zirconate titanate constituents, named PZT-4 and PZT-5A, respectively. Before presenting
the numerical results, a comparison with previous works is performed for validating the present
formulation and methodology.

U
AN
4.1 Comparison and validation
The results presented in the papers by Zenkour [28] and Singha et al. [33] are taken as references
and used for the validation process. To this aim, a functionally graded square plate made of
M

Alumina and Aluminum with a h = 10 is considered. A sinusoidal transverse load of maximum

magnitude q0 is applied to the structure. In order to perform the comparison, the following
D

dimensionless quantities must be introduced

TE

h a b h h a b h h  h
σx = σ x  , , , σ y = σ y  , ,  , τ xy = τ xy  0,0, −  (81)
aq0  2 2 2  aq0  2 2 2  aq0  3
EP

h a h h  b  10h3 Ec  a b 
τ yz = τ yz  , 0,  , τ xz = τ xz  0, , 0  , w= w , ,c  (82)
aq0 2 6 aq0  2  a 4 q0 2 2 
where Ec is the Young’s modulus of the ceramic constituent. The results of the comparison tests
C

are shown in Tables 2-3. A good agreement can be observed between the proposed approach and
AC

the reference solutions [28, 33]. In particular, it can be noted that the dimensionless central
displacement and in-plane stresses increase significantly for higher values of the volume fraction
exponent r . This variation is due to the decrease of the stiffness of the FGP plate. In fact, the
flexural rigidity of the FGP plate decreases noticeably and the central displacement and in-plane
stresses increase as a consequence.

14
 ACCEPTED MANUSCRIPT

4.2 Electro-mechanical bending analysis

The results for a FGP plate subjected to applied electric potential and transverse loads are
presented in this section. The structure is made of a combination of PZT-4 and PZT-5A. Table 4
shows the maximum dimensionless and shear stresses σ x q , σ y q , τ xz q , τ yz q and τ xy q

evaluated at the top of the plate, as a function of different values of r and applied potential φ0 .

PT
The maximum stresses in hand occur at
a b a b a   b

RI
σ x = σ x  ,  , σ y = σ y  ,  , τ xy = τ xy ( 0,0 ) , τ yz = τ yz  , 0  , τ xz = τ xz  0,  (83)
 
2 2 
2 2  2   
2
The numerical results presented in Table 4 shows that by increasing the volume fraction

SC
exponents r , the stresses σ x q , σ y q , τ xy q decrease for r < 2 , whereas the same quantities

increase for r ≥ 2 . This behavior is due to the non-homogeneities of the material properties. In

U
addition, it can be noted that the maximum in-plane normal stresses σ x q , σ y q decrease by
AN
increasing the value of applied electric potential φ0 . On the other hand, the in-plane shear stress

τ xy q assumes higher values. The variation of the electric potential has no considerable effects
M

for the extreme cases ( r = 0, r = ∞ ). Tables 5 and 6 show the same dimensionless stresses in

terms of various volume fraction exponents r and applied electric potentials φ0 at neutral and
D

lower surfaces of FGP plate, respectively. In particular, the results in Table 5 prove that τ xy q is
TE

equal to zero at the neutral surface, whereas τ xz q and τ yz q reach their maximum values. In

addition, it can be noted that σ x q and σ y q increase for r < 2 and decrease for r ≥ 2 . On the
EP

other hand, the values of τ xz q and τ yz q are lower if the exponent r is bigger. The results in

Table 6, instead, show that the normal stresses σ x q , σ y q and the shear stress τ xy q at the
C

lower surface of the plate decrease for higher values of r .

AC

Table 7 shows the dimensionless bending deflection 105 wb a , shear deflection 108 ws a and

maximum electric potential φ of the FGP plate, as a function of various volume fraction
exponents r and applied electric potential φ0 . Results in Table 7 prove that the dimensionless
deflections slightly increase for higher values of r , due to the stiffness reduction of the structure.
In addition, the maximum electric potential φ decreases if the value of r is increased. As far as

15
 ACCEPTED MANUSCRIPT

the influence of the applied electric potential φ0 is concerned, it can be observed that the

dimensionless deflections increase for bigger values of φ0 . On the other hand, the maximum

electric potential φ decreases significantly by increasing the value of φ0 . It can be deducted that

the influence of linear term 2 znsφ0 h in equation (24) is greater than the effect given by the

PT
cosine term. Consequently, its effect can be neglected.
The variation c h is shown in Figure 2 in terms of various volume fraction exponents r . The

RI
graph shows that c h initially increases for bigger value of r and reaches its maximum value for
r ≅ 1.6 . On the other hand, the location of the neutral surface is closer to the middle surface for

SC
r > 1.6 , whereas it coincides with the middle surface for r = ∞ .
Figures 3 and 4 present the variations of the dimensionless bending deflection 105 wb a and

U
shear deflection 108 ws a as a function of the volume fraction exponent r , for different values
AN
of the applied electric potential φ0 . The graphs prove that the deflections at issue increases

significantly for higher values of φ0 . In addition, the influence of the volume fraction exponent
M

r on the bending and shear deflections strongly depends on the applied electric potentials φ0 .
According to the results shown in Figure 3, it can be deduced that the maximum bending
D

deflection occurs for r ≅ 1.6 when the applied electric potential is taken as φ0 > 4 V . If the
TE

applied electric potential is assumed as φ0 ≤ 4 V , the maximum bending deflection increases

uniformly. On the other hand, the graphs in Figure 4 show that the maximum shear deflection is
reached for 2 ≤ r ≤ 5 . Figure 5 shows the variation of the maximum electric potential φ in terms
EP

of the volume fraction exponent r , for different values of the applied electric potential φ0 .

Analogously, its maximum value is reached for r ≅ 1.6 .

C

Figures 6-8 present respectively the variations of σ x q = σ y q , τ xz q = τ yz q , and τ xy q in

AC

terms of zns h for different values of applied electric potential φ0 , assuming r = 1 . It can be
noted that in-plane shear stress is not affected by the applied electric potential. Furthermore, the
minimum value of the in-plane shear stress occurs at the neutral surface, whereas its maximum
value is obtained at the inner and outer surfaces. On the other hand, the maximum value of the
out-of-plane shear stress is gotten at the neutral surface, whereas its minimum value occurs at the
inner and outer surfaces.

16
 ACCEPTED MANUSCRIPT

Figures 9-11 present respectively the variations of σ x q = σ y q , τ xz q = τ yz q , and τ xy q in

terms of zns h for different values of applied electric potential φ0 , assuming r = 5 . The same
considerations are still valid in this circumstance.
Finally, the last graphs depicted in Figures 12-17 show the influence of the Pasternak foundation

PT
parameters k1 , k2 on the structural response. Figures 12 and 13 present the two-dimensional

variation of the dimensionless bending and shear deflections ( 105 wb a , 108 ws a ) in terms of

RI
k1 , k2 , respectively. The values of the exponent is assumed as r = 1.6 , whereas the applied

electric potentials is given by φ0 = 1V . It can be noted that the displacements at issue decrease

SC
significantly for higher values of k1 , k2 . Figure 14 shows the two-dimensional variation of the

maximum electric potential φ as a function of k1 , k2 . It can be easily deducted that its value is

U
lower by increasing the foundation parameters. Finally, Figures 15-17 present the variations of
AN
the dimensionless stresses σ x q = σ y q , τ xz q = τ yz q , and τ xy q in terms of the Pasternak

foundation parameters k1 , k2 . It can be observed that all the stress components decrease
M

significantly for higher values of the foundation parameters.

D

5. CONCLUSIONS
TE

In this paper, the concept of neutral surface is employed for the electro-elastic analysis of a FGP
plate resting on Pasternak foundation. A two-variable sinusoidal shear deformation plate theory
is used for the description of the displacement field. The governing equations are obtained by
EP

means of the principle of virtual works. Several investigations are carried out and the solutions
are obtained analytically. The following remarks can be introduced:
C

- The increase of the volume fraction exponent r leads to important changes in terms of
stress components. These variations are not uniform and depend strongly on the value of
AC

r . In particular, the in-plane normal and shear stresses decrease for r < 2 , whereas they
increase for r ≥ 2 . The out-of-plane shear stresses uniformly decrease for higher value of
r . In addition, both maximum bending and shear deflections and maximum electric
potential increase noticeably for greater values of the volume fraction exponent r .
- The increase of the applied electric potential φ0 leads to a decrease of the maximum
values of the in-plane normal stresses and in-plane shear stress. The effect of the applied

17
 ACCEPTED MANUSCRIPT

electric potential can be neglected for r = 0 and r = ∞ . In addition, the bending and
shear deflections increase for higher value of φ0 . Finally, the maximum electric potential

φ decreases significantly by increasing the value of φ0 .

PT
APPENDIX

RI
 A11 B11 D11 B11s D11s H11s  h2 −c  c11 
 
 A12 B12 D12 B12s D12s H12s  = ∫  c12  1 zns zns2 f ( zns ) zns f ( zns ) f 2 ( zns )  dzns
 A66 B66 D66 B66s D66s H 66s  − h2 −c c66 


SC
h
−c
2
s
A44 = A55s = ∫ g 2 ( zns ) dzns

U
h
− −c
2

 P11 P12  2h − c  1 
AN
P π πz  
 13 P14  = ∫ sin  ns   zns  [ e31 e32 ] dzns
h  h 
 P15 P16  − h2 −c  f ( zns ) 
M

h
−c
2
πz 
[ P17 P18 ] = ∫ g ( z ) cos  h  [e
ns
ns
15 e24 ] dzns
D

h
− −c
2
h
−c
TE

πz 
2
[ P20 P22 ] = ∫ cos 2  ns  [κ11 κ 22 ] dzns
h
− −c
 h 
2
EP

h
−c
2π πz 
2
P29 = ∫ 2
κ 33 sin  ns  φ0 dzns
h
− −c
h  h 
2
C

h
−c 2
π  π z 
2
P30 = ∫ κ 33  sin  ns   φ0 dzns
AC

h
− −c
h  h 
2

 P31 P32  h2 −c  1 
P  2  
 33 P34  = ∫ φ0  zns  [ e31 e32 ] dzns
h
 P35 P36  − 2h −c  f ( zns ) 

18
 ACCEPTED MANUSCRIPT

REFERENCES

[1] Tornabene F., Viola E., Free Vibration Analysis of Functionally Graded Panels and Shells
of Revolution. Meccanica, 2009, 44, 255-281.
[2] Viola E., Tornabene F., Free Vibrations of Three Parameter Functionally Graded Parabolic
Panels of Revolution. Mech. Res. Commun., 2009, 36, 587-594.
[3] Tornabene F., Viola E., Free Vibrations of Four-Parameter Functionally Graded Parabolic

PT
Panels and Shells of Revolution. Eur. J. Mech. A/Solids, 2009, 28, 991-1013.
[4] Tornabene F., Free Vibration Analysis of Functionally Graded Conical, Cylindrical Shell
and Annular Plate Structures with a Four-parameter Power-Law Distribution. Comput.

RI
Method. Appl. M., 2009, 198, 2911-2935.
[5] Tornabene F., Viola E., Inman D.J., 2-D Differential Quadrature Solution for Vibration
Analysis of Functionally Graded Conical, Cylindrical Shell and Annular Plate Structures.

SC
J. Sound Vib., 2009. 328, 259-290.
[6] Tornabene F., Liverani A., Caligiana G., FGM and Laminated Doubly-Curved Shells and
Panels of Revolution with a Free-Form Meridian: a 2-D GDQ Solution for Free Vibrations.
Int. J. Mech. Sci., 2011, 53(6), 446-470.

U
[7] Viola E., Rossetti L., Fantuzzi N., Numerical Investigation of Functionally Graded
Cylindrical Shells and Panels Using the Generalized Unconstrained Third Order Theory
AN
Coupled with the Stress Recovery. Compos. Struct., 2012, 94, 3736-3758.
[8] Tornabene F., Viola E., Static Analysis of Functionally Graded Doubly-Curved Shells and
Panels of Revolution. Meccanica, 2013, 48, 901-930.
M

[9] Tornabene F., Reddy J.N., FGM and Laminated Doubly-Curved and degenerate Shells
resting on Nonlinear Elastic Foundation: a GDQ Solution for Static Analysis with a
Posteriori Stress and strain Recovery. J. Indian I. Sci, 2013, 93, 635-688.
[10] Viola E., Rossetti L., Fantuzzi N., Tornabene F., Static Analysis of Functionally Graded
D

Conical Shells and Panels Using the Generalized Unconstrained Third Order Theory
Coupled with the Stress Recovery. Compos. Struct., 2014, 112, 44-65.
TE

[11] Brischetto S., Tornabene F., Fantuzzi N., Viola E., 3D Exact and 2D Generalized
Differential Quadrature Models for Free Vibration Analysis of Functionally Graded Plates
and Cylinders. Meccanica, 2016, 51, 2059-2098.
EP

[12] Fantuzzi N., Brischetto S., Tornabene F., Viola E., 2D and 3D Shell Models for the Free
Vibration Investigation of Functionally Graded Cylindrical and Spherical Panels. Compos.
Struct., 2016, 154, 573-590.
[13] Kandasamy R., Dimitri R., Tornabene F., Numerical Study on the Free Vibration and
C

Thermal Buckling Behaviour of Moderately Thick Functionally Graded Structures in

Thermal Environments. Compos. Struct., 2016, 157, 207-221.
AC

[14] Viola E., Rossetti L., Fantuzzi N., Tornabene F., Generalized Stress-Strain Recovery
Formulation Applied to Functionally Graded Spherical Shells and Panels Under Static
Loading. Compos. Struct., 2016, 156, 145-164.
[15] Tornabene F., Brischetto S., Fantuzzi N., Bacciocchi M., Boundary Conditions in 2D
Numerical and 3D Exact Models for Cylindrical Bending Analysis of Functionally Graded
Structures. Shock Vib., 2016, Vol. 2016, Article ID 2373862, 1-17.
[16] Tornabene F., Fantuzzi N., Bacciocchi M., Viola E., Effect of Agglomeration on the
Natural Frequencies of Functionally Graded Carbon Nanotube-Reinforced Laminated
Composite Doubly-Curved Shells. Compos. Part B-Eng., 2016, 89, 187-218.

19
 ACCEPTED MANUSCRIPT

[17] Tornabene F., Fantuzzi N., Bacciocchi M., Linear Static Response of Nanocomposite
Plates and Shells Reinforced by Agglomerated Carbon Nanotubes, Compos. Part B-Eng.,
2017, 115(1), 449-476.
[18] Fantuzzi N., Tornabene F., Bacciocchi M., Dimitri R., Free Vibration Analysis of
Arbitrarily Shaped Functionally Graded Carbon Nanotube-Reinforced Plates. Compos.
Part B-Eng., 2017, 115, 384-408.
[19] Banić D., Bacciocchi M., Tornabene F., Ferreira A.J.M., Influence of Winkler-Pasternak

PT
Foundation on the Vibrational Behavior of Plates and Shells Reinforced by Agglomerated
Carbon Nanotubes. Appl. Sci., 2017, 7, 1228, 1-55.
[20] Reddy J.N., Wang C.M., Kitipornchai S., Axisymmetric bending of functionally graded

RI
circular and annular plates. Euro. J. Mech. A/Solids, 1999, 18, 185-199.
[21] Reddy J.N., Analysis of functionally graded plates. Int. J. Numer. Meth. Engng., 2000, 47,
663–684.

SC
[22] Reddy J.N., Cheng Z.Q., Three-dimensional thermomechanical deformations of
functionally graded rectangular plates. Eur. J. Mech. A/Solids, 2001, 20, 841–855.
[23] Yang J., Shen H.S., Nonlinear bending analysis of shear deformable functionally graded
plates subjected to thermo-mechanical loads under various boundary conditions. Compos.

U
Part B Eng., 2003, 34, 103–15.
[24] Lanhe W., Thermal buckling of a simply supported moderately thick rectangular FGM
AN
plate. Compos. Struct., 2004, 64, 211–218.
[25] Kim Y.W., Temperature dependent vibration analysis of functionally graded rectangular
plates. J. Sound. Vib., 2005, 284, 531–549.
[26] Qian L.F., Batra R.C., Chen L.M., Static and dynamic deformations of thick functionally
M

graded elastic plates by using higher-order shear and normal deformable plate theory and
meshless local Petrov–Galerkin method. Compos. Part B Eng., 2004, 35, 685–97.
[27] Abrate S., Free vibration, buckling, and static deflections of functionally graded plates.
D

Compos. Sci. Technol., 2006, 66, 2383–94.

[28] Zenkour A.M., Generalized shear deformation theory for bending analysis of functionally
TE

graded plates. Appl. Math. Model., 2006, 30, 67–84.

[29] Zenkour A.M.: Benchmark trigonometric and 3-D elasticity solutions for an exponentially
graded thick rectangular plate. Arch. Appl. Mech., 2007, 77, 197–214.
EP

[30] Chi S.H., Chung Y.L., Mechanical behavior of functionally graded material plates under
transverse load-Part I: Analysis. Int. J. Solids Struct., 2006, 43, 3657–74.
[31] Prakash T., Singha M.K., Ganapathi M., Influence of neutral surface position on the
nonlinear stability behavior of functionally graded plates. Comput. Mech., 2009, 43, 341–
C

50.
[32] Zhang D.G., Zhou Y.H., A theoretical analysis of FGM thin plates based on physical
AC

neutral surface. Comput. Mater. Sci., 2008, 44, 716–20.

[33] Singha M.K., Prakash T., Ganapathi M., Finite element analysis of functionally graded
plates under transverse load. Finite Elem. Anal. Des., 2011, 47, 453–60.
[34] Alinia M.M., Ghannadpour S.A.M., Nonlinear analysis of pressure loaded FGM plates.
Compos. Struct., 2009, 88, 354–359.
[35] Kulikov G.M., Plotnikova S.V., An analytical approach to three-dimensional coupled
thermoelectroelastic analysis of functionally graded piezoelectric plates. J. Intel. Mater.
Sys. Struct., 2017, 28, 435-450.

20
 ACCEPTED MANUSCRIPT

[36] Jandaghian A.A., Rahmani O., Size-dependent free vibration analysis of functionally
graded piezoelectric plate subjected to thermo-electro-mechanical loading. J. Intel. Mater.
Sys. Struct., 2017, In Press.
[37] Matsunaga H., Stress analysis of functionally graded plates subjected to thermal and
mechanical loadings. Compos. Struct., 2009, 87, 344–57.
[38] Ebrahimi F., Rastgoo A., Free vibration analysis of smart annular FGM plates integrated
with piezoelectric layers. Smart Mater. Struct., 2008, 17, 15044.

PT
[39] Zhong Z., Yu T., Vibration of a simply supported functionally graded piezoelectric
rectangular plate. Smart Mater. Struct., 2006, 15, 1404–1412.
[40] Xu T.F., Xing Y.F., Closed-form solutions for free vibration of rectangular FGM thin

RI
plates resting on elastic foundation. Acta Mech. Sin., 2016, 32, 1088–1103.
[41] Mechab I., Atmane H.A., Tounsi A., Belhadj H.A., Bedia E.A.A., A two variable refined
plate theory for the bending analysis of functionally graded plates. Acta Mech. Sin., 2010,

SC
26, 941–9.
[42] Wang Y., Xu R.Q., Ding H.J., Free axisymmetric vibration of FGM circular plates. Appl.
Math. Mech. -Engl. Ed., 2009, 30, 1077–1082.
[43] Arefi M., Zenkour A.M., A simplified shear and normal deformations nonlocal theory for

U
bending of functionally graded piezomagnetic sandwich nanobeams in magnetthermo-
electric environment. J. Sand. Struct. Mater., 2016, 18, 624–651.
AN
[44] Arefi M., Zenkour A.M., Nonlocal electro-thermo mechanical analysis of a sandwich
nanoplate containing a Kelvin–Voigt viscoelastic nanoplate and two piezoelectric layers.
Acta Mech., 2017, 228, 475–493.
[45] Arefi M., Zenkour A.M., Vibration and bending analysis of a sandwich microbeam with
M

two integrated piezo-magnetic face sheets. Compos. Struct., 2017, 159, 479–490.
[46] Arefi M., Zenkour A.M., Employing sinusoidal shear deformation plate theory for transient
analysis of three layers sandwich nanoplate integrated with piezo-magnetic face-sheets.
D

Smart Mater. Struct., 2016, 25, 115040.

[47] Arefi M., Zenkour A.M., Thermo-electro-mechanical bending behavior of sandwich
TE

nanoplate integrated with piezoelectric face-sheets based on trigonometric plate theory.

Compos. Struct., 2017, 162, 108-122.
[48] Ebrahimi F., Salari E., Thermal buckling and free vibration analysis of size dependent
EP

Timoshenko FG nanobeams in thermal environments. Compos. Struct., 2015, 128, 363–

380.
[49] Zenkour A.M., Allam M.N.M., Radwan A.F., El-Mekawy H.F., Thermo-Mechanical
Bending Response of Exponentially Graded Thick Plates Resting on Elastic Foundations.
C

Int. J. Appl. Mech., 2015, 7, 1–24.

[50] Giannakopoulos A.E., Suresh S., Theory of indentation of piezoelectric materials. Acta
AC

Mater., 1999, 47, 2153–2164.

[51] Tabib-Azar M., Microactuators: Electrical, Magnetic, Thermal, Optical, Mechanical,
Chemical and Smart Structures. Springer Science & Business Media, 1998.
[52] Liew K.M., He X.Q., Tan M.J., Lim H.K., Dynamic analysis of laminated composite plates
with piezoelectric sensor/actuator patches using the FSDT mesh-free method. Int. J. Mech.
Sci., 2004, 46, 411-431-
[53] Gürses M, Civalek O, Korkmaz A, Ersoy H., Free vibration analysis of symmetric
laminated skew plates by discrete singular convolution technique based on first-order shear
deformation theory, Int. J. Numer. Meth. Eng., 2009, 79, 290-313.

21
 ACCEPTED MANUSCRIPT

[54] Baltacioglu AK, Akgoz B, Civalek O., Nonlinear static response of laminated composite
plates by discrete singular convolution method. Compos. Struct., 2010, 93, 153-161.
[55] Fereidoon A, Asghardokht SM, Mohyeddin A., Bending analysis of thin functionally
graded plates using generalized differential quadrature method, Arch. Appl. Mech., 2011,
81, 1523-1539.
[56] Demir C, Mercan K, Civalek, O., Determination of critical buckling loads of isotropic,
FGM and laminated truncated conical panel, Comps. Part B-Eng., 2016, 94, 1-10.

PT
[57] Tornabene F., Ceruti A., Mixed Static and Dynamic Optimization of Four-Parameter
Functionally Graded Completely Doubly-Curved and Degenerate Shells and Panels Using
GDQ Method. Math. Probl. Eng., 2013, Vol. 2013, Article ID 867079, 1-33.

RI
[58] Tornabene F., Fantuzzi N., Bacciocchi M., Viola E., Reddy J.N., A Numerical
Investigation on the Natural Frequencies of FGM Sandwich Shells with Variable
Thickness by the Local Generalized Differential Quadrature Method. Appl. Sci., 2017, 7,

SC
131, 1-39.
[59] Tornabene F., Fantuzzi N., Bacciocchi M., Free Vibrations of Free-Form Doubly-Curved
Shells Made of Functionally Graded Materials Using Higher-Order Equivalent Single
Layer Theories. Compos. Part B-Eng., 2014, 67, 490-509.

U
[60] Tornabene F., Fantuzzi N., Viola E., Batra R.C., Stress and Strain Recovery for
Functionally Graded Free-Form and Doubly-Curved Sandwich Shells Using Higher-Order
AN
Equivalent Single Layer Theory. Compos. Struct., 2015, 119, 67-89.
[61] Fantuzzi N., Tornabene F., Viola E., Four-Parameter Functionally Graded Cracked Plates
of Arbitrary Shape: a GDQFEM Solution for Free Vibrations. Mech. Adv. Mater. Struct.,
2016, 23, 89-107.
M
D
TE
C EP
AC

22
 ACCEPTED MANUSCRIPT

PT
RI
SC
Figure 1 - Schematic representation of a FGP plate resting on Pasternak foundation.

U
AN
M
D
TE
C EP
AC

Figure 2 – Variation of c h as a function of r .

23
 ACCEPTED MANUSCRIPT

PT
RI
U SC
AN
Figure 3 – Variation of the dimensionless bending deflection 105 wb a as a function of r , for different values of
M

applied electric potential φ0 .

D
TE
C EP
AC

24
 ACCEPTED MANUSCRIPT

PT
RI
U SC
AN
Figure 4 – Variation of the dimensionless shear deflection 108 ws a as a function of r , for different values of
applied electric potential φ0 .
M
D
TE
C EP
AC

25
 ACCEPTED MANUSCRIPT

PT
RI
U SC
AN
Figure 5 – Variation of the maximum electric potential φ as a function of r , for different values of applied electric
potential φ0 .
M
D
TE
C EP
AC

26
 ACCEPTED MANUSCRIPT

PT
RI
U SC
AN
Figure 6 – Variation of the dimensionless stresses σ x q = σ y q as a function of zns h , for different values of
applied electric potential φ0 , assuming r = 1 .
M
D
TE
C EP
AC

27
 ACCEPTED MANUSCRIPT

PT
RI
U SC
AN
Figure 7 – Variation of the dimensionless stresses τ xz q = τ yz q as a function of zns h , for different values of
applied electric potential φ0 , assuming r = 1 .
M
D
TE
C EP
AC

28
 ACCEPTED MANUSCRIPT

PT
RI
U SC
AN
Figure 8 – Variation of the dimensionless stress τ xy q as a function of zns h , for different values of applied
electric potential φ0 , assuming r = 1 .
M
D
TE
C EP
AC

29
 ACCEPTED MANUSCRIPT

PT
RI
U SC
AN
Figure 9 – Variation of the dimensionless stresses σ x q = σ y q as a function of zns h , for different values of
applied electric potential φ0 , assuming r = 5 .
M
D
TE
C EP
AC

30
 ACCEPTED MANUSCRIPT

PT
RI
U SC
AN
Figure 10 – Variation of the dimensionless stresses τ xz q = τ yz q as a function of zns h , for different values of
applied electric potential φ0 , assuming r = 5 .
M
D
TE
C EP
AC

31
 ACCEPTED MANUSCRIPT

PT
RI
U SC
AN
Figure 11 – Variation of the dimensionless stress τ xy q as a function of zns h , for different values of applied
electric potential φ0 , assuming r = 5 .
M
D
TE
C EP
AC

32
 ACCEPTED MANUSCRIPT

PT
RI
U SC
AN
Figure 12 – Two-dimensional variation of the dimensionless bending deflection 105 wb a in terms of the
M

foundation parameters k1 , k2 .
D
TE
C EP
AC

33
 ACCEPTED MANUSCRIPT

PT
RI
U SC
AN
Figure 13 – Two-dimensional variation of the dimensionless shear deflection 108 ws a in terms of the foundation
parameters k1 , k2 .
M
D
TE
C EP
AC

34
 ACCEPTED MANUSCRIPT

PT
RI
U SC
AN
Figure 14 – Two-dimensional variation of the electric potential φ in terms of the foundation parameters k1 , k2 .
M
D
TE
C EP
AC

35
 ACCEPTED MANUSCRIPT

PT
RI
U SC
AN
Figure 15 – Two-dimensional variation of the dimensionless stresses σ x q = σ y q in terms of the foundation
parameters k1 , k2 .
M
D
TE
C EP
AC

36
 ACCEPTED MANUSCRIPT

PT
RI
U SC
AN
Figure 16 – Two-dimensional variation of the dimensionless stresses τ xz q = τ yz q in terms of the foundation
parameters k1 , k2 .
M
D
TE
C EP
AC

37
 ACCEPTED MANUSCRIPT

PT
RI
U SC
AN
Figure 17 – Two-dimensional variation of the dimensionless stress τ xy q in terms of the foundation parameters
k1 , k2 .
M
D
TE
C EP
AC

38
 ACCEPTED MANUSCRIPT

Table 1 – Mechanical and electrical properties of the constituents [50].

Properties Lead zirconate titanate Lead zirconate titanate

PZT-4 PZT-5A
Young’s modulus (GPa) 75 61
Poisson’s ratio 0.3 0.3
Elastic stiffness coefficients (GPa)

PT
c11 = c22 139 121

c12 77.8 75.4

c44 = c55 25.6 21.1

RI
c66 30.6 22.8

Piezoelectric coefficients ( C m 2 )

SC
e31 = e32 -5.2 -5.4

e15 = e24 12.7 12.3

Dielectric constants ( 10-9 F m )

U
κ11 = κ 22 6.461 8.1066

κ 33 5.620 7.3455
AN
M
D
TE
C EP
AC

39
 ACCEPTED MANUSCRIPT

Table 2 - Comparison of central displacement w and in-plane stresses σ x , σ y for various volume fraction
exponents r .

Central displacement w In-plane stress σ x In-plane stress σ y

r Zenkour Singha et Zenkour Singha et Zenkour
Present Present Present
[28] al. [33] [28] al. [33] [28]

PT
Ceramic 0.2960 0.2960 0.2961 1.9955 1.9955 1.9679 1.3121 1.3121
1 0.5889 0.5889 0.5891 3.0870 3.0870 3.0389 1.4894 1.4894
2 0.7573 0.7573 0.7582 3.6094 3.6094 3.5456 1.3954 1.3954

RI
3 0.8377 0.8377 - 3.8742 3.8742 - 1.2748 1.2748
4 0.8819 0.8819 0.8831 4.0693 4.0693 3.9813 1.1783 1.1783

SC
5 0.9118 0.9118 - 4.2488 4.2488 - 1.1029 1.1029
6 0.9356 0.9356 - 4.4244 4.4244 - 1.0417 1.0417

U
7 0.9562 0.9562 - 4.5971 4.5971 - 0.9903 0.9903
8 0.9750 0.9750 - 4.7661 4.7661 - 0.9466 0.9466
AN
9 0.9925 0.9925 - 4.9303 4.9303 - 0.9092 0.9092
10 1.0089 1.0089 - 5.0890 5.0890 - 0.8775 0.8775
M

Metal 1.6070 1.6070 1.6072 1.9955 1.9955 1.9679 1.3121 1.3121

D
TE
C EP
AC

40
 ACCEPTED MANUSCRIPT

Table 3 - Comparison of shear stresses τ yz ,τ xz ,τ xy for various volume fraction exponents r .

Shear stress τ yz Shear stress τ xz Shear stress τ xy

r
Present Zenkour [28] Present Zenkour [28] Present Zenkour [28]
Ceramic 0.2132 0.2132 0.2462 0.2462 0.7065 0.7065
1 0.2622 0.2622 0.2462 0.2462 0.6110 0.6110

PT
2 0.2763 0.2763 0.2265 0.2265 0.5441 0.5441
3 0.2715 0.2715 0.2107 0.2107 0.5525 0.5525

RI
4 0.2580 0.2580 0.2029 0.2029 0.5667 0.5667
5 0.2429 0.2429 0.2017 0.2017 0.5755 0.5755
6 0.2296 0.2296 0.2041 0.2041 0.5803 0.5803

SC
7 0.2194 0.2194 0.2081 0.2081 0.5834 0.5834
8 0.2121 0.2121 0.2124 0.2124 0.5856 0.5856
9 0.2072 0.2072 0.2164 0.2164 0.5875 0.5875

U
10 0.2041 0.2041 0.2198 0.2198 0.5894 0.5894
AN
Metal 0.2132 0.2132 0.2462 0.2462 0.7065 0.7065
M
D
TE
C EP
AC

41
 ACCEPTED MANUSCRIPT

Table 4 – Dimensionless normal and shear stresses σ x q , σ y q , τ xz q , τ yz q , τ xy q at the top of the FGP plate
( zns = h 2 − c ) for various volume fraction exponents r and applied electric potential.

r φ0 σx q , σy q τ xz q , τ yz q τ xy q
0 28.5488 0 -7.7756
r =0
1 27.5088 0 -7.7756
(PZT-4)
2 26.4688 0 -7.7756

PT
0 27.9617 0 -7.6429
1 1 26.6353 0 -7.6451
2 25.3088 0 -7.6473

RI
0 27.9369 0 -7.6404
1.6 1 26.6015 0 -7.6425
2 25.2661 0 -7.6445
0 27.9485 0 -7.6454

SC
2 1 26.6202 0 -7.6475
2 25.2919 0 -7.6495
0 28.0051 0 -7.6641
3 1 26.7066 0 -7.6662

U
2 25.4081 0 -7.6683
0 28.0661 0 -7.6831
AN
4 1 26.7981 0 -7.6852
2 25.5302 0 -7.6872
0 28.1201 0 -7.6997
5 1 26.8785 0 -7.7017
2 25.6369 0 -7.7037
M

r=∞ 0 28.6150 0 -7.8504

(PZT- 1 27.5750 0 -7.8504
5A) 2 26.5350 0 -7.8504
D
TE
C EP
AC

42
 ACCEPTED MANUSCRIPT

Table 5 – Dimensionless normal and shear stresses σ x q , σ y q , τ xz q , τ yz q , τ xy q at the neutral surface of the
FGP plate ( zns = 0 m ) for various volume fraction exponents r and applied electric potential.

r φ0 σx q , σy q τ xz q , τ yz q τ xy q
0 0 0.7400 0
r=0
1 -1.0400 0.7400 0
(PZT-4)
2 -2.0800 0.7400 0

PT
0 0.3251 0.6947 0
1 1 -0.7441 0.6876 0
2 -1.8133 0.6806 0

RI
0 0.3409 0.6852 0
1.6 1 -0.7358 0.6778 0
2 -1.8124 0.6705 0
0 0.3356 0.6807 0

SC
2 1 -0.7439 0.6735 0
2 -1.8235 0.6663 0
0 0.3068 0.6732 0
3 1 -0.7762 0.6666 0

U
2 -1.8591 0.6601 0
0 0.2755 0.6690 0
AN
4 1 -0.8083 0.6632 0
2 -1.8922 0.6574 0
0 0.2480 0.6666 0
5 1 -0.8359 0.6615 0
2 -1.9198 0.6564 0
M

r=∞ 0 0 0.6491 0
(PZT- 1 -1.0800 0.6491 0
5A) 2 -2.1600 0.6491 0
D
TE
C EP
AC

43
 ACCEPTED MANUSCRIPT

Table 6 – Dimensionless normal and shear stresses σ x q , σ y q , τ xz q , τ yz q , τ xy q at the lower surface of the
FGP plate ( zns = − h 2 − c ) for various volume fraction exponents r and applied electric potential.

r φ0 σx q , σy q τ xz q , τ yz q τ xy q
0 -28.5488 0 7.7756
r=0
1 -29.5888 0 7.7756
(PZT-4)
2 -30.6288 0 7.7756

PT
0 -26.5469 0 5.9498
1 1 -27.3482 0 5.9560
2 -28.1495 0 5.9622

RI
0 -26.5690 0 5.9588
1.6 1 -27.3631 0 5.9654
2 -28.1571 0 5.9719
0 -26.5583 0 5.9578

SC
2 1 -27.3590 0 5.9642
2 -28.1598 0 5.9706
0 -26.5081 0 5.9485
3 1 -27.3362 0 5.9541

U
2 -28.1643 0 5.9597
0 -26.4555 0 5.9378
AN
4 1 -27.3118 0 5.9426
2 -28.1681 0 5.9473
0 -26.4097 0 5.9283
5 1 -27.2906 0 5.9324
2 -28.1715 0 5.9365
M

0 -26.0292 0 5.8493
r=∞
1 -27.1092 0 5.8493
(PZT-5A)
2 -28.1892 0 5.8493
D
TE
C EP
AC

44
 ACCEPTED MANUSCRIPT

Table 7 – Dimensionless bending deflection 105 wb a , shear deflection 108 ws a and maximum electric potential
φ for various volume fraction exponents r and applied electric potential φ0 .

r φ0 105 wb a 108 ws a φ
0 0.10268 0.8439 -0.6145
0.5 0.10268 0.8439 -0.6145
r =0

PT
1 0.10268 0.8439 -0.6145
(PZT-4)
2 0.10268 0.8439 -0.6145
3 0.10268 0.8439 -0.6145
0 0.10306 0.8550 -0.5619

RI
0.5 0.10308 0.8607 -0.4896
0.5 1 0.10311 0.8664 -0.4172
2 0.10316 0.8777 -0.2726
3 0.10321 0.8891 -0.1280

SC
0 0.10318 0.8726 -0.5437
0.5 0.10321 0.8799 -0.4536
1 1 0.10324 0.8871 -0.3634
2 0.10331 0.9016 -0.1830

U
3 0.10337 0.9160 -0.0026
0 0.10323 0.8902 -0.5340
AN
0.5 0.10327 0.8978 -0.4412
1.6 1 0.10330 0.9053 -0.3484
2 0.10337 0.9205 -0.1628
3 0.10344 0.9356 0.0227
M

0 0.10326 0.8990 -0.5301

0.5 0.10329 0.9065 -0.4396
2 1 0.10332 0.9139 -0.3490
2 0.10339 0.9288 -0.1678
D

3 0.10345 0.9438 0.0134

0 0.10329 0.9132 -0.5241
0.5 0.10332 0.9200 -0.4426
TE

3 1 0.10335 0.9269 -0.3612

2 0.10341 0.9405 -0.1982
3 0.10347 0.9541 -0.0353
0 0.10332 0.9205 -0.5201
EP

0.5 0.10335 0.9266 -0.4480

4 1 0.10337 0.9326 -0.3760
2 0.10342 0.9448 -0.2318
3 0.10347 0.9569 -0.0877
C

0 0.10334 0.9240 -0.5171

0.5 0.10337 0.9294 -0.4531
AC

5 1 0.10339 0.9349 -0.3892

2 0.10344 0.9457 -0.2613
3 0.10348 0.9565 -0.1334
0 0.10365 0.9074 -0.4929
0.5 0.10365 0.9074 -0.4929
r=∞
1 0.10365 0.9074 -0.4929
(PZT-5A)
2 0.10365 0.9074 -0.4929
3 0.10365 0.9074 -0.4929

45