**O R I G I NA L**

**Le Minh Thai** **·** **Doan Trac Luat** **·** **Van Binh Phung** **·**
**Phung Van Minh** **·** **Do Van Thom**

**Finite element modeling of mechanical behaviors** **of piezoelectric nanoplates with flexoelectric effects**

Received: 30 November 2020 / Accepted: 14 September 2021 / Published online: 17 November 2021

© The Author(s), under exclusive licence to Springer-Verlag GmbH Germany, part of Springer Nature 2021

**Abstract** This paper uses the finite element method to simulate the mechanical, electric, and polarization
behaviors of piezoelectric nanoplates resting on elastic foundations subjected to static loads, in which the
flexoelectric effect is taken into consideration. The finite element formulations are established by employing
a new type of shear deformation theory, which does not need any shear correction factors, but still accurately
describes the stress field of the plate. The numerical results show clearly that the flexoelectric effect has a strong
influence on the mechanical responses of the nanoplates. In particular, the normal stress distribution in the
thickness direction is no longer linear when the flexoelectric coefficient is large enough, and this phenomenon
differs completely from that of conventional plates. In addition, the distribution of the electric field and the
polarization also depend on boundary conditions, which were not investigated in the published works.

**Keywords** Nanoplates·Elastic foundation·Flexoelectric·Bending·Finite element method
**1 Introduction**

Nowadays, along with the development of science and technology, piezoelectric nanostructures have been used widely in engineering practices such as sensors, actuators, energy harvesters, and so on. For these structures, flexoelectricity is a common phenomenon, especially the electric polarization inside structures due to strain gradients. There have also been several studies showing that the mechanical behavior of the nanoplates takes into account the flexoelectric effect. Yan [1] used the classical plate theory to study static bending and free vibration of piezoelectric nanoplates, taking into account the flexoelectricity based on the analytical solution. Yang et al.

[2] introduced an analytical solution to show the static bending and free vibration response of nanoplates, in which the theory considered the effects of piezoelectric and flexoelectricity. The authors also relied on Kirchhoff’s plate theory (classical plate theory—CPT) to provide the solutions to be found in explicit form. Li et al. [3] studied static bending and free vibration of the circular microplate. The equations were established from the CPT, the solution was solved analytically, and the results pointed out the effect of flexoelectric effect on the mechanical response of this plate. Still, based on the classical plate theory, Wang et al. [4] employed the finite difference solution to show the static bending response of flexoelectric nanoplates. The authors also calculated the case where one edge of the plate was clamped. He and his co-workers [5] used the CPT combined with the exact solution to show clearly the effect of both the flexoelectricity and the piezoelectricity L. M. Thai

Faculty of Special Equipments, Le Quy Don Technical University, Hanoi City, Vietnam D. T. Luat·P. V. Minh·D. V. Thom (

### B

^{)}

Faculty of Mechanical Engineering, Le Quy Don Technical University, Hanoi City, Vietnam e-mail: thom.dovan@lqdtu.edu.vn

V. B. Phung

Faculty of Aerospace Engineering, Le Quy Don Technical University, Hanoi City, Vietnam

on the mechanical response of nanoplates. Ebrahimi1 and Barati [6] researched the bucking response of flexoelectric nanoplates based on the exact solution and the classical plate theory, where the plate rested on a two-parameter elastic foundation, and the impact of the stress surface was taken into consideration. Amin and his colleagues [7] employed the classical plate theory and the flexoelectric theory to establish the analytical solution to analyze the nonlinear free vibration of functionally graded (FG) flexoelectric nanoplates. Amir et al. [8] studied the free vibration of sandwich flexoelectric plates, in which the first-order shear deformation theory and Navier’s method were mixed to give a series of free vibration responses of the plate. Ghobadi et al.

[9] explored the nonlinear static bending of flexoelectric nanoplates. The analytic solution was incorporated into classical Kirchhoff’s plate theory to find the displacement and electric field of these plates subjected to the thermo-electromagnetic force. Based on published papers, Ghobadi and his team [10] investigated the static bending response of FG nanoplates, taking into account the effects of temperature, electric, and flexoelectric fields. Giannakopoulos et al. [11] introduced an antiplane dynamic flexoelectric issue, which was defined as a dielectric solid with electric polarization and flexoelectricity gradients owing to strain gradients. Qu and colleagues [12] investigated the torsion of a rectangular cross-sectional flexoelectric semiconductor rod. It was based on the macroscopic theory of flexoelectric semiconductors. Using double power series expansion of the coordinates inside the cross section, a one-dimensional model was established from the three-dimensional theory. Yang et al. [13] researched multilayer phononic crystal band structures with flexoelectricity. In this work, the classical plate theory of Kirchhoff was used to derive the finite element formulations.

From the summary of previously published works, it can be seen that these studies mainly relied on the classical plate theory to determine the mechanical response of nanoplates. It can be understood that the CPT is the simple theory to establish the equilibrium equation of structures. However, it cannot show all the behavior of nanoplates when taking into account the effects of the flexoelectric effect. Currently, there is no study showing the distribution of stress, electric field, and polarization of nanoplates based on shear deformation theories that correctly show the mechanical response of the plate, including the satisfaction of stress boundary conditions, has to be zero at the surfaces of the plate. As a result, this paper aims to use the finite element method to simulate the mechanical, electric, and polarization behaviors of piezoelectric nanoplates resting on an elastic foundation subjected to static load, taking into account the flexoelectric effect by employing the finite element method and the new type of shear deformation theory.

The remainder of this work is structured as follows: Sect.2develops finite element formulas to address the bending issue of nanoplates while accounting for the flexoelectric effect. Section3continues to provide comparisons to demonstrate the validity of the suggested theory and the calculation program coded in the MATLAB environment. Section4presents a set of findings from the static bending response of nanoplates, in which the impact of flexoelectric on the plate’s stress distribution, electric field, and polarization is clearly shown. Section5highlights the study’s major points.

**2 Finite element model of piezoelectric nanoplate in static bending problem**

Consider a mechanical structure of the piezoelectric nanoplate and the elastic foundation placed in the Cartesian
coordinate system with the length*a*, width*b*, and thickness*h*as shown in Fig.1.

In the analysis of plate structures, there are many different theories, such as classical deformation theory, first-order shear deformation theory, and higher-order shear deformation theory. The classical shear deforma- tion (CDT) theory is only suitable for thin structures because the effect of shear deformation is not taken into account, and the first-order shear deformation theory (FSDT) gives the acceptable results. However, it needs to add a shear correction factor into calculations. Higher-order shear deformation theories do not require a shear correction coefficient, but their finite element formulations are more complicated than those of the CDT and FSDT. As a result, Shimpi [14] developed a simple plate theory by isolating the transverse displacement component into bending and shear parts. The most amazing feature of Shimpi’s approach is that it has fewer unknown variables and governing equations than the CDT and FSDT, and does not need a shear correction factor. According to Shimpi’s work, some four unknown shear deformation theories were modified using some shape functions, such as polynomial functions [15], sinusoidal functions [16], and hyperbolic functions [17].

This work uses the simple form theory proposed by Shimpi with different shape functions such as sinusoidal functions as in works [18,19]. The main advantage of this theory is simplicity, no need for any shear correc- tions since it is satisfied that the shear stress boundary at the top and bottom surface of the plate is zero, and it accurately describes the mechanical response of the plate.

**Fig. 1** The model of a piezoelectric nanoplate resting on a two-parameter elastic foundation

According to the shear deformation theory based on hyperbolic sine functions [18,19], the displacements
*u*,*v*, and*w*in the*x-*,*y-*, and*z-*directions, respectively, depending on the coordinates of one point, have the
following forms:

⎧⎪

⎪⎪

⎪⎨

⎪⎪

⎪⎪

⎩

*u**x*(*x,y,z*) −*z∂w*^{b}(*x,y,*0)

*∂x* − *f*(*z*)*∂w*^{s}(*x,y,*0)

*∂x*
*v**y*(*x,y,z*) −*z∂w*^{b}(*x,y,*0)

*∂y* − *f*(*z*)*∂w*^{s}(*x,y,*0)

*∂y*
*w**z*(*x,y,z*)*w**b*(*x,y,*0) +*w**s*(*x,y,*0)

(1)

with *f*(*z*)*z*−*ζ*(*z*)* ,ζ*(

*z*)

*h.*sin

_{h}

^{z}−

*z.*cosh

^{1}

_{2}, where

*u*

*x*

*, v*

*y*, and

*w*

*z*are the displacements in the

*x-*,

*y-*, and

*z*-directions at one point within the plate;

*w*

*b*and

*w*

*s*are the bending displacement and shear displacement in the

*z*-direction. Therefore, to determine the longitudinal displacements and tangential displacements at any point with the coordinates (

*x*,

*y*,

*z*), two main components,

*w*

*b*and

*w*

*s*need to be defined.

Longitudinal and shear strain components are defined as follows:

⎧⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎨

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎪⎪

⎩
**ε**

⎧⎨

⎩
*ε**x x*

*ε**yy*

*γ**x y*

⎫⎬

⎭

⎧⎪

⎪⎪

⎪⎪

⎪⎪

⎨

⎪⎪

⎪⎪

⎪⎪

⎪⎩

−*z∂*^{2}*w*^{b}

*∂x*^{2} − *f*(*z*)*∂*^{2}*w*^{s}

*∂x*^{2} *ε**bx* +*zε**sx*

−*z∂*^{2}*w*^{b}

*∂y*^{2} − *f*(*z*)*∂*^{2}*w*^{s}

*∂y*^{2} *εby*+*zε**sy*

−*z∂*^{2}*w**b*

*∂x∂y* −2*f*(*z*)*∂*^{2}*w**s*

*∂x∂y* *γ**bx y*+*zγ**sx y*

⎫⎪

⎪⎪

⎪⎪

⎪⎪

⎬

⎪⎪

⎪⎪

⎪⎪

⎪⎭
*z*

⎧⎨

⎩
*ε**bx*

*ε**by*

*γ*_{bx y}

⎫⎬

⎭
*ε**b*

+*f*(*z*)

⎧⎨

⎩
*ε**sx*

*ε**sy*

*γ**sx y*

⎫⎬

⎭

*ε**s*

**γ***γ*^{x z}

*γ*^{yz}

*∂ζ (z)*

*∂z*

⎧⎪

⎪⎨

⎪⎪

⎩

*∂w**s*

*∂x*

*∂w**s*

*∂y*

⎫⎪

⎪⎬

⎪⎪

⎭ *∂ζ (z)*

*∂z* **γ**_{0}

(2)

In this work, only the strain gradients in the*x*- and*y*-axes are considered, and the strain gradient in the
*z*-axis is zero. It is obvious due to the displacement field chosen. This means that the strain gradients in the
*x*- and*y*-axes are much higher than that in the thickness direction. Then, the strain gradient is expressed as
follows:

**η**

⎧⎪

⎪⎨

⎪⎪

⎩

*η*^{x x z} −*∂*^{2}*w**b*

*∂x*^{2} −*∂f*(*z*)

*∂z*

*∂*^{2}*w**s*

*∂x*^{2}
*η*^{yyz} −*∂*^{2}*w**b*

*∂y*^{2} −*∂f*(*z*)

*∂z*

*∂*^{2}*w**s*

*∂y*^{2}

⎫⎪

⎪⎬

⎪⎪

⎭

−^{∂}_{∂}^{2}_{x}^{w}2^{b}

−^{∂}_{∂}^{2}_{y}^{w}2^{b}

+*∂f*(*z*)

*∂z*

−^{∂}_{∂}^{2}_{x}^{w}2^{s}

−^{∂}_{∂}^{2}_{y}^{w}2^{s}

**η***b*+*∂f*(*z*)

*∂z* **η***s* (3)

When the flexoelectric effect is taken into consideration, the stress components and electric displacement vector for a nanoscale dielectric material are expressed as follows:

*T**i j* *c**i j kl**ε**kl*−*e**ki j**E**k*

^{i j m} −*f**ki j m**E**k*

*P**i* *c**i j k**ε**j k*+*κ**i j**E**k*+ *f**i j kl**η**j kl* (4)

in which*c**i j kl*,*e**ki j*, *f**ki j m*, and*κ*^{i j}are the components of elastic, piezoelectric, flexoelectric, and permittivity
constant tensors; they are the material parameters. *T**i j* is the stress tensor, which is similar to that of the
traditional elastic foundation. *P**i* is the electric displacement vector, and^{i j m}is the moment stress tensor or
the higher-order stress tensor.

From the expressions of the strain components, we have specific expressions of stress and electric dis- placement vectors as follows:

**T**

⎡

⎣*c*_{11}*c*_{12} 0
*c*_{12}*c*_{11} 0
0 0 *c*_{66}

⎤

⎦

⎧⎨

⎩
*ε**x*

*ε**y*

*γ**x y*

⎫⎬

⎭−*e*_{31}

⎧⎨

⎩ 1 1 0

⎫⎬

⎭*E**z***C***b** ε*− ˜

*(5)*

**E**

**S***s**x z*

*s**yz*

*c*_{66} 0
0 *c*_{66}

**γ****C***s** γ* (6)

**Ψ***x x z*

*yyz*

−*f*_{14}
1

1

*E**z* (7)

*P*_{z}^{0}*e*_{31}

*ε*^{x x}+*ε*^{yy}

+*κ*33*E**z*+ *f*_{14}

*η*^{x x z}+*η*^{yyz}

(8)
where*f*_{14}*f*_{3113}and*f*_{14}*f*_{3223}[20].

Herein, some assumptions are made as follows: there is no external electric field acting on the plate, the
electric displacement equals the electric polarization. Hence, the last element in Eq. (8) *f*_{14}

*η**x x z*+*η**yyz*

is the
polarization, which is caused by the strain gradients in this plate; and this element depends on the coordinate
*z*of the plate due to the derivative of the function*f (z)*in the expression of the strain gradient.

The electrical field is calculated from the partial derivative of electrical potential as follows:

*E**z* −*∂ϕ*

*∂z* (9)

According to Gaussian’s law in electrostatics, the electric displacement is given by the formula:

*∂P*_{z}^{0}

*∂z* 0 (10)

From Eqs. (8), (9), and (10), we have:

*ϕ* *e*_{31}
2*a*_{33}

*∂*^{2}*w*

*∂x*^{2} +*∂*^{2}*w*

*∂y*^{2}

*z*^{2}+*c*_{1}+*c*_{2}*z* (11)

where*c*_{1}and*c*_{2}are the coefficients to be found. With the open-circuit condition, we have the following
conditions:

*P*_{z}^{0}

±*h*
2

0 (12)

Then, the following coefficient can be found:

*c*_{2} −*f*_{14}
*a*_{33}

*∂*^{2}*w*

*∂x*^{2} +*∂*^{2}*w*

*∂y*^{2}

(13) Substituting Eqs. (11) and (13) into Eq. (9), we obtain the expression of the internal electric field as follows:

*E**z* −*e*_{31}
*κ*33

*ε**x x*+*ε**yy*

+ *f*_{14}
*κ*33

*∂*^{2}*w**b*

*∂x*^{2} +*∂*^{2}*w**b*

*∂y*^{2}

+ *f*_{14}
*κ*33

*∂*^{2}*w**s*

*∂x*^{2} +*∂*^{2}*w**s*

*∂y*^{2}

*∂f*(*z*)

*∂z*
*e*_{31}

*κ*33

*∂*^{2}*w*^{b}

*∂x*^{2} +*∂*^{2}*w*^{b}

*∂y*^{2}

*z*+*e*_{31}
*κ*33

*∂*^{2}*w*^{s}

*∂x*^{2} +*∂*^{2}*w*^{s}

*∂y*^{2}

*f*(*z*)
+ *f*_{14}

*κ*33

*∂*^{2}*w**b*

*∂x*^{2} +*∂*^{2}*w**b*

*∂y*^{2}

+ *f*_{14}
*κ*33

*∂*^{2}*w**s*

*∂x*^{2} +*∂*^{2}*w**s*

*∂y*^{2}

*∂f*(*z*)

*∂z* (14)

Note that the expression of the electric field*E**z*depends on the coordinate*z*, function*f*(*z*), and the derivative
of the function*f (z)*. Therefore,*E**z*may convert roughly to linear or nonlinear form depending on the value of
*f*_{14}; it means the flexoelectric effect will significantly affect the distribution of*E**z*over the thickness direction,
which will be specified in the example below. This is a very interesting problem. However, there are no
published works dealing with this issue.

With the open-circuit condition, the electric Gibbs free energy has the following form:

*U* 1
2

*V*

**ε**^{T}* T* +

**γ**^{T}

*+*

**S**

**η**^{T}

**Ψ**d*V* 1

2

*V*

⎛

⎝**ε**^{T}* Cε* +

**γ**^{T}

*−*

**S**

**ε**^{T}

*e*

_{31}

⎧⎨

⎩ 1 1 0

⎫⎬

⎭*E**z*−**η**^{T} *f*_{14}
1

1

*E**z*

⎞

⎠d*V*

−1 2

*V*

⎛

⎝**ε**^{T}*e*_{31}*.e*_{31}
*k*_{33}

⎧⎨

⎩ 1 1 0

⎫⎬

⎭
*∂*^{2}*w**s*

*∂x*^{2} +*∂*^{2}*w**s*

*∂y*^{2}

*f*(*z*)

⎞

⎠d*V*

−1 2

*V*

⎛

⎝**ε**^{T}*e*_{31}

⎧⎨

⎩ 1 1 0

⎫⎬

⎭
*f*14

*k*_{33}
*∂*^{2}*w*^{b}

*∂x*^{2} +*∂*^{2}*w*^{b}

*∂y*^{2}

⎞⎠d*V* −1
2

*V*

⎛

⎝**ε**^{T}*e*_{31}

⎧⎨

⎩ 1 1 0

⎫⎬

⎭
*f*14

*k*_{33}
*∂*^{2}*w*^{s}

*∂x*^{2} +*∂*^{2}*w*^{s}

*∂y*^{2}

*∂f*(*z*)

*∂z*

⎞

⎠d*V*

−1 2

*V*

**η**^{T} *f*14

1 1

*e*_{31}
*k*_{33}

*∂*^{2}*w**b*

*∂x*^{2} +*∂*^{2}*w**b*

*∂y*^{2}

*.z*

d*V* −1
2

*V*

**η**^{T} *f*14

1 1

*e*_{31}
*k*_{33}

*∂*^{2}*w**s*

*∂x*^{2} +*∂*^{2}*w**s*

*∂y*^{2}

*f*(*z*)

d*V*

−1 2

*V*

**η**^{T} *f*_{14}

1 1

*f*_{14}
*k*_{33}

*∂*^{2}*w**b*

*∂x*^{2} +*∂*^{2}*w**b*

*∂y*^{2}

d*V* −1
2

*V*

**η**^{T} *f*_{14}

1 1

*f*_{14}
*k*_{33}

*∂*^{2}*w**s*

*∂x*^{2} +*∂*^{2}*w**s*

*∂y*^{2}

*∂f*(*z*)

*∂z*

d*V*
(15)
Since the plate is rested on an elastic foundation, the potential energy expression of the plate, taking into
account the effect of the elastic foundation, has the following form:

*U*^{found} 1
2

*S*

"

*k*_{w}*w*^{2}+*k**s*

"

*∂w*

*∂x*
2

+
*∂w*

*∂y*
2##

d*S* (16)

where*k**w*and*k**s*are the two coefficients of the elastic foundation.

In this paper, a four-node element is used. Each node has six degrees of freedom:

**q**_{e}

$4
*i*1

*w**bi**, w**si**,*
*∂w**b*

*∂x*

*i*

*,*
*∂w**s*

*∂x*

*i*

*,*
*∂w**b*

*∂y*

*i*

*,*
*∂w**s*

*∂y*

*i*

*T*

(17) where

⎧⎪

⎪⎪

⎪⎪

⎨

⎪⎪

⎪⎪

⎪⎩
*w*^{b}

$4
*i*1

*H**i**w*^{bi}+*H*_{i+1}
*∂w**b*

*∂x*

*i*

+*H*_{i+2}
*∂w**b*

*∂y*

*i*

**H***b***q**_{e}
*w**s*

$4
*i*1

*H**i**w**si* +*H**i*+1

*∂w**s*

*∂x*

*i*

+*H*_{i+2}
*∂w**s*

*∂y*

*i*

**H***s***q**_{e}

(18)

in which,*H**j*is the Hermit interpolation function.

Then, the displacement vector at any point in the element is interpolated through the element’s nodal displacement vector as follows:

**u**

*w**b**, w**s**,*
*∂w**b*

*∂x*

*,*
*∂w**s*

*∂x*

*,*
*∂w**b*

*∂y*

*,*
*∂w**s*

*∂y*
*T*

**H**.**q**_{e} (19)
According to this definition, the strain vectors are computed through the nodal displacement vector as
follows:

**ε***b*

⎡

⎢⎢

⎣

−^{∂}^{2}**H***b*

*∂**x*^{2}

−^{∂}^{2}**H***b*

*∂**y*^{2}

−2^{∂}^{2}**H***b*

*∂**x y*

⎤

⎥⎥

⎦**q**_{e} **B**_{1}**q**_{e}; **ε***s*

⎡

⎢⎢

⎣

−^{∂}^{2}**H***s*

*∂**x*^{2}

−^{∂}^{2}**H***s*

*∂**y*^{2}

−2^{∂}^{2}**H***s*

*∂**x y*

⎤

⎥⎥

⎦**q**_{e} **B**_{2}**q**_{e} (20)

* γ*0
'

_{∂}

**H***s*

*∂**x*

*∂***H***s*

*∂**y*

(

**q**_{e} * B*3

**q**_{e};

**η***b*

'−^{∂}^{2}**H***b*

*∂**x*^{2}

−^{∂}^{2}**H***b*

*∂**y*^{2}

(

**q**_{e} * B*4

**q**_{e};

**η***s*

'−^{∂}^{2}**H***s*

*∂**x*^{2}

−^{∂}^{2}**H***s*

*∂**y*^{2}

(

**q**_{e} **B**_{5}**q**_{e} (21)

So Eq. (15) is rewritten as follows:

*U**e* 1
2**q**_{e}^{T}

*V*_{e}

)

*z B*

^{T}

_{1}+

*f*(

*z*)

**B**^{T}

_{2}

**C***b*

*z B*

_{1}+

*f*(

*z*)

**B**_{2}

+ **B**^{T}_{3}**C***s***B**_{3}

*
d*V q*

_{e}

−1
2**q**^{T}_{e}

*V*_{e}

*z B*

_{1}

^{T}+

*f*(

*z*)

**B**^{T}

_{2}

*e*_{31}**I**_{3}*e*_{31}
*k*_{33}

*∂*^{2}**H***b*

*∂x*^{2} +*∂*^{2}**H***b*

*∂y*^{2}

*.z*

d*V q*

_{e}

−1
2**q**^{T}_{e}

*V*_{e}

*z B*

_{1}

^{T}+

*f*(

*z*)

**B**^{T}

_{2}

*e*

_{31}

*.e*

_{31}

*k*

_{33}

**I**_{3}

*∂*^{2}**H***s*

*∂x*^{2} +*∂*^{2}**H***s*

*∂y*^{2}

*f*(*z*)

d*V q*

_{e}

−1
2**q**^{T}_{e}

*V*_{e}

*z B*

_{1}

^{T}+

*f*(

*z*)

**B**^{T}

_{2}

*e*_{31}**I**_{3}*f*14

*k*_{33}
*∂*^{2}**H***b*

*∂x*^{2} +*∂*^{2}**H***b*

*∂y*^{2}

d*V q*

_{e}

−1
2**q**^{T}_{e}

*V*_{e}

*z B*

^{T}

_{1}+

*f*(

*z*)

**B**_{2}

^{T}

*e*_{31}**I**_{3} *f*14

*k*_{33}
*∂*^{2}**H***s*

*∂x*^{2} +*∂*^{2}**H***s*

*∂y*^{2}

*∂f*(*z*)

*∂z*

d*V q*

_{e}

−1
2**q**^{T}_{e}

*V*_{e}

**B**^{T}_{4} +*∂f*(*z*)

*∂z* **B**^{T}_{5}

*f*_{14}**I**_{2}*e*31

*k*_{33}
*∂*^{2}**H***b*

*∂x*^{2} +*∂*^{2}**H***b*

*∂y*^{2}

*.z*

d*V q*

_{e}

−1
2**q**^{T}_{e}

*V*_{e}

**B**^{T}_{4} +*∂f*(*z*)

*∂z* **B**^{T}_{5}

*f*_{14}**I**_{2}*e*31

*k*_{33}
*∂*^{2}**H***s*

*∂x*^{2} +*∂*^{2}**H***s*

*∂y*^{2}

*f*(*z*)

d*V q*

_{e}

−1
2**q**^{T}_{e}

*V*_{e}

**B**^{T}_{4} +*∂f*(*z*)

*∂z* **B**^{T}_{5}

*f*_{14}**I**_{2} *f*14

*k*_{33}
*∂*^{2}**H***b*

*∂x*^{2} +*∂*^{2}**H***b*

*∂y*^{2} d*V q*

_{e}

−1
2**q**^{T}_{e}

*V*_{e}

**B**^{T}_{4} +*∂f*(*z*)

*∂z* **B**^{T}_{5}

*f*_{14}**I**_{2} *f*_{14}
*k*_{33}

*∂*^{2}**H***s*

*∂x*^{2} +*∂*^{2}**H***s*

*∂y*^{2}

*∂f*(*z*)

*∂z*

d*V q*

_{e}(22)

with**I**_{3}{1,1,0}^{T}và**I**_{2}{1,1}^{T}.

Equation (22) can be written in short form as follows:

*U**e* 1

2**q**^{T}_{e} **K***e***q**_{e} (23)

**Fig. 2** Boundary conditions of piezoelectric nanoplates

Equation (16) relates to the element nodal displacement vector as follows:

*U*_{e}^{found} 1
2**q**_{e}^{T}

*S*_{e}

⎛

⎜⎜

⎜⎝*k*_{w}*( H*

*b*+

**H***s*

*)*

^{T}

*(*

**H***b*+

**H***s*

*)*+

*k*

*s*

⎛

⎜⎜

⎜⎝
*∂ H*

*b*

*∂x* +*∂ H*

*s*

*∂x*

*T**∂ H*

*b*

*∂x* +*∂ H*

*s*

*∂x*

+
*∂ H*

*b*

*∂y* +*∂ H*

*s*

*∂y*

*T**∂ H*

*b*

*∂y* +*∂ H*

*s*

*∂y*

⎞

⎟⎟

⎟⎠

⎞

⎟⎟

⎟⎠d*S q*

_{e}

1

2**q**_{e}^{T}**K***e*^{f}**q**_{e} (24)

The external force exerted on the plate has the following formula:

*W**e*

*S**e*

*(w**b*+*w**s**)*^{T}**F***S*d*S q*

^{T}

_{e}

*S**e*

*( H*

*b*+

**H***s*

*)*

^{T}

**F***S*d

*S*

*F*_{e}

**q**^{T}_{e}**F***e* (25)

So the static equilibrium equation of the piezoelectric nanoplate the flexoelectricity is obtained from the
minimal equation*U**e*,*U*_{e}^{found}, and*W**e*is as follows:

**K***e*+**K***e*^{f}

**q**_{e} **F***e* (25)

By solving Eq. (26), the displacement and stress fields will be obtained, then substituting them into Eqs. (8)
and (14),*E**z*and*P**z*will be obtained due to strain gradients of this plate.

In this paper, some boundary conditions are defined as follows [18]:

• As one edge is simply supported:

*w*^{b}0*, w*^{s}0*,* *∂w**b*

*∂y* 0*,* *∂w**s*

*∂y* 0 at *x* 0*,a*
*w**b*0*, w**s*0*,* *∂w*^{b}

*∂x* 0*,* *∂w*^{s}

*∂x* 0 at *y*0*,b*

(27a)

• As one edge is clamped:

*w*^{b}0*, w*^{s}0*,* *∂w**b*

*∂x* 0*,* *∂w**b*

*∂y* 0*,* *∂w**s*

*∂x* 0*,* *∂w**s*

*∂y* 0 (27b)

and*S*represents for a simply supported boundary condition,*C*represents a clamped supported boundary
condition. Some boundary conditions used in this paper are shown in Fig.2.

**Table 1** Non-dimensional deflection*w*of the plate resting on the two-parameter elastic foundation
*K*_{w}^{∗} *K*_{s}^{∗} *a*/*h*10

[21] [22] This work

8×8 elements

10×10 elements 12×12 elements 14×14 elements

1 5 3.3455 3.3455 3.3797 3.3643 3.3560 3.3512

10 2.7505 2.7504 2.7743 2.7635 2.7578 2.7545

15 2.3331 2.3331 2.3508 2.3428 2.3386 2.3361

20 2.0244 2.0244 2.0382 2.0320 2.0287 2.0268

81 5 2.8422 2.8421 2.8667 2.8557 2.8499 2.8464

10 2.3983 2.3983 2.4163 2.4083 2.4040 2.4015

15 2.0730 2.0730 2.0868 2.0806 2.0774 2.0755

20 1.8245 1.8244 1.8355 1.8306 1.8280 1.8265

625 5 1.3785 1.3785 1.3835 1.3816 1.3804 1.3797

10 1.2615 1.2615 1.2658 1.2642 1.2632 1.2626

15 1.1627 1.1627 1.1665 1.1650 1.1642 1.1637

20 1.0782 1.0782 1.0815 1.0802 1.0795 1.0791

**3 Verification study**

In this section, some comparison problems will be carried out to verify the accuracy of the proposed theory and mathematical model. The numerical results of the deflection and stress in this work are compared with those of exact published data.

* Example 1* Consider a plate resting on an elastic foundation with dimensions

*ab*0.2 m,

*ha*/10 and

*a*/200. The material properties are

*E*320.24 GPa, Poisson’s ratio of 0.26. The plate is

*SSSS*and under a uniformly distributed load

*q*

_{0}. The non-dimensional elastic foundation parameters are calculated as:

⎧⎪

⎪⎨

⎪⎪

⎩

*K*_{w}^{∗} *K*_{w}*a*^{4}
*D*
*K*_{s}^{∗} *K**s**a*^{2}

*D*

(28)

with

*D* *E h*^{3}
12

1−*ν*^{2} (29)

The non-dimensional deflection at the center point of the plate is defined by the following formula:

*w* 10^{3}*D*
*q*_{0}*a*^{4} *w*

*a*
2*,b*

2

(30) The comparative results obtained by this work, the differential quadrature method [21], and the analytical method [22] are presented in Table1with the increase in the mesh size, and one can see that with the mesh size of 12×12, the accuracy is acceptable. Therefore, in all the following investigations, this mesh size will be used.

* Example 2* Consider a square plate with

*a*/

*b*1,

*a*/

*h*10,

*E*380 GPa, and

*ν*0.3. The plate is fully simply supported, and the load is bi-sinusoidal type as follows:

*qq*_{0}sin*πx*
*a*

sin*πy*

*b*

(31) Non-dimensional stresses are calculated as:

⎧⎪

⎪⎨

⎪⎪

⎩

*σ**x*^{∗∗}*, σ**y*^{∗∗}

*σ**x*

*q*0

*a*
2*,b*

2*,z*

*,* *σ**y*

*q*0

*a*
2*,b*

2*,z*
*τ**x y*^{∗∗} *τ*^{x y}*(*0*,*0*,*z*)*

*q*_{0}

(32)

**a** **b**

**c** **d**

**Fig. 3** The stress distributions in the thickness direction of the square plate under bi-sinusoidal load

For different volume fraction index*n*, the stress distributions in the thickness direction are presented in
Fig.3. Herein, the numerical results of this work are compared with those of Hiroyuki’s exact solution [23]. It
can be seen that the stresses are distributed nonlinearly in the thickness direction, and they reach their maximum
values at the top and bottom surfaces.

* Example 3* Finally, non-dimensional maximum defection of the

*SSSS*nanoplate with

*h*20 nm,

*ab*50

*h*, and material properties

*c*

_{11}102 GPa;

*c*

_{12}31 GPa;

*c*

_{33}35.50 GPa;

*e*

_{31}−17.05 C/m

^{2};

*k*

_{33}1.76.10

^{–8}C/(Vm);

*f*

_{14}10

^{–7}C/m is considered. The plate is under a uniformly distributed load of

*q*

_{0}0.05 MPa. The numerical results of this work compared with those of Yang et al.’s analytical solution [2] are presented in Fig.4. It can see that the data has good agreement.

**4 Numerical results**

In this section, a piezoelectric nanoplate with*h*20 nm,*ab*50*h*, and the material properties*c*_{11}
102 GPa;*c*_{12}31 GPa;*c*_{33}35.50 GPa;*e*_{31} −17.05 C/m^{2};*k*_{33}1.76.10^{–8}C/(Vm) is considered; the
uniformly distributed load is*q*_{0}0.05 MPa. The plate is rested on a two-parameter elastic foundation with
*k**w*and*k**s*. The boundary conditions are as follows:

• Fully simply supported plate: SSSS

• Fully clamped plate: CCCC

• Two opposite edges are clamped, the other edges are simply supported: CSCS

• Two adjacent edges are clamped, the other edges are simply supported: CCSS

**Fig. 4** Non-dimensional deflection at*y**b*/2 taking into account the flexoelectric effect

Non-dimensional parameters are calculated as follows:

*w*∗ 10^{3}*c*11*h*^{3}

12*q*_{0}*a*^{4} *w*; *f*_{14}^{∗} *f*14

*f*_{14}^{0}; *K*_{w}^{∗} *k*_{w}*a*^{4}
*D**f*

; *K*_{s}^{∗} *k**s**a*^{2}
*D**f*

; *D**f* *c*11*h*^{3}
) 12

*σ**x*^{∗}*, τ**x y*^{∗}*, τ**x z*^{∗}

* 1
*q*0

)*σ**x**, τ**x y**, τ**x z*

*

(33)
with *f*_{14}^{0} 10^{–7}C/m.

Since this study calculates solely in the case of disregarding external voltage, only an evenly distributed
mechanical load*q*_{0}is applied to the plate. Therefore, in order to clearly see the response of the electric field
and the polarization (due to strain gradient*P**z* *f*_{14}

*η**x x z*+*η**yyz*

) corresponding to this mechanical load, this work gives the normalized parameters of the electric field and polarization as follows:

*E*_{z}^{∗} *E**z*

*q*_{0}
*V m*

*N*

; *P*_{z}^{∗} *P**z*

*q*_{0}
*C*

*N*

; (34)

4.1 Influence of the flexoelectric effect

Consider an*SSSS* plate, to investigate the influence of the flexoelectric effect with the values of*f*_{14}are 0
and 10^{–7}C/m (*f*_{14}^{∗} 1,*K*_{w}^{∗} 100,*K*_{s}^{∗}10). This means the plate with and without the flexoelectric effect
is considered. The non-dimensional deflection*w**at*yb*/2 is presented in Fig.5. The distributions of the
electric fields*E**z*,*P**z*, and the stress components in the thickness direction at the center point of the plate are
shown in Figs.6,7,8,9, and10. Some discussions are as follows:

• Fig. 5points out that when taking into account the flexoelectric effect*,*the non-dimensional maximum
deflection*w*of the plate decreases. However, the normal stress*σ**x* and shear stress*τ**x z*increase (Figs.8and
10). This represents a complete difference from conventional plates and is also something to keep in mind
when testing nanoplates in the presence of the flexoelectric effect.

• Figs.6,7,8,9and10show that when not taking the flexoelectric effect into the calculation (and*P**z*0)*,*the
electric field*E**z*and the normal stress*σ**x*symmetrically across the middle plane of the plate. However, when
*f*_{14}is not equal to zero,*E**z*, and the stress*σ**x*are no longer symmetrical across the middle plane of the plate,
and maximum values of*P**z*and*σ**x*increase. In particular, this is completely different from traditional plates.

That is the maximum displacement of*w*decreases, the stress*σ**x* increases. On the other hand, maximum
values of*σ**x* and*τ**x z*increase when the value of*f*_{14}is not equal to zero. However, the value of*τ**x y*decreases
when taking into account the flexoelectric effect.

**Fig. 5** The deflection of the plate at*y**b*/2 in the cases of with and without the flexoelectric effect

**Fig. 6** The distribution of the electric field*E**z*along the thickness direction in the cases of with and without the flexoelectric effect

**Fig. 7** The distribution of*P**z*across the thickness in the cases of with and without the flexoelectric effect

**Fig. 8** The distribution of the stress*σ**x*^{∗}across the thickness in the cases with and without the flexoelectric effect

**Fig. 9** The distribution of the stress*τ**x y*^{∗} across the thickness in the cases with and without the flexoelectric effect

**Fig. 10** The distribution of the stress*τ**x z*^{∗} across the thickness in the cases with and without the flexoelectric effect

**Fig. 11** The deflection of the plate at*y**b*/2 with different values of*f*14

4.2**Influence of the value of f**_{14}

The appearance of the influence of the flexoelectricity effect is shown in the coefficient*f*_{14}, and this will make
the mechanical response of the nanoplate show a difference compared to conventional structures. Since*E**z*

depends on the coefficient*f*_{14}as in Eq. (14) and makes the polarization-dependent on*f*_{14}as in expression (8),
the stress depends on*f*_{14}as in expression (4), so the energy of the plate changes. Thus, the response of the
nanoplate also changes. The coefficient*f*_{14}will vary depending on the material, and this work investigates
how*f*_{14}varies in a certain range to see clearly how the flexoelectricity effect affects the deflection response,
stress, electric field, and polarization of nanoplates. The computed data can be used as a reference in design,
manufacturing, and practical applications.

To investigate the influence of the value of*f*_{14}on the behavior of the plate (*SSSS*,*K*_{w}^{∗} 100,*K*_{s}^{∗}10),
the value of *f*_{14}^{∗} increases gradually from 1 to 10. The deflection of the plate at*yb*/2 is shown in Fig.11;

the distributions along the thickness direction at the center point of the plate of*E**z*,*P**z*, and stress components
are presented in Figs.12,13,14, and15. Some discussions are as follows:

• When increasing the value of*f*_{14}, maximum deflection of the plate*w*decreases (Fig.11).

• By looking at Figs.11,12,13, and14, when the value of*f*_{14}is small,*E**z*,*P**z**,*and*σ**x*distribute linearly in the
thickness of the plate. However, when the value of*f*_{14}is high, these components change nonlinearly along
the plate thickness. This is also a complete difference from plates made of a homogeneous material whose
mechanical properties do not vary in the thickness direction. It can be explained that in the expressions of
*E**z*,*P**z**,*and*σ**x* including the multiplication of*f*_{14}and the function*g*(*z*) vary nonlinearly in the thickness
direction. Therefore, when the value of*f*_{14}is small, it has a slight influence on the function*g*(*z*). However,
when the value of*f*_{14}is high, its effect on the function*g*(*z*) becomes stronger. Hence, the values of*E**z*,*P**z**,*
and*σ**x* change nonlinearly along the thickness direction. In addition, one can see that the increase in*f*_{14}
will make the law of transformation of the maximum values of*E**z*,*P**z*, and*σ**x*completely different from the
transformation of the maximum deflection*w*_{max}. And the surface position of the plate is always the position
where*E**z*,*P**z**,*and*σ**x* reach their maximum value.

• The stress component*τ**x y* decreases gradually with increasing the value of *f*_{14}. It is proportional to the
decrease in deflection*w*(see Fig.15), as is the case of conventional plates.

4.3 Influence of boundary condition

In this subsection, to see more clearly the influence of the boundary condition on the distribution of the electric
field*E**z*and the polarization*P**z*along the*x*and*y*edges of the plate, four boundary conditions, SSSS, CCCC,
CSCS, and CCSS (*f*_{14}^{∗} 3, *K*_{w}^{∗} 100,*K*_{s}^{∗}10) are considered. The distribution surfaces of*E**z*

*z* −^{h}_{2}
and*P**z*

*z* −^{h}_{2}

along the*x*and*y*edges are presented in Fig.16. The value lines at*yb*/2 are shown in
Figs.17and18. The minimum and maximum values of*w*,*E**z*, and*P**z*are listed in Table2. One can see that:

• Boundary conditions affect strongly not only the shapes of the distribution surfaces of*E**z*and*P**z*along the
*x-*and*y*-axes, but also to the maximum values of these components.

**Fig. 12** The distribution of the electric field*E*_{z}along the thickness direction with different values of*f*_{14}

**Fig. 13** The distribution of*P**z*along the thickness direction with different values of*f*14

**Fig. 14** The distribution of*σ**x*^{∗}along the thickness direction with different values of*f*14

**Fig. 15** The distribution of*τ**x y*^{∗} along the thickness direction with different values of*f*14

**Table 2** The minimum and maximum values of*w, E**z**,*and*P**z*with different boundary conditions,*z* −^{h}_{2}

Value Boundary condition

*SSSS* *CCCC* *SCSC* *CCSS*

*w*max 1.292 0.5403 0.7455 0.8213

*E*_{z(max)}[Vm/N] 7.845 5.595 7.576 7.446

*E**z*(min)[Vm/N] 0.146 4.344 5.271 5.873

*P**z*(max)[C/N] 8.146×10^{–8} 4.171×10^{–8} 6.184×10^{–8} 6.098×10^{–8}

*P**z*(min)[C/N] 2.260×10^{–10} 4.107×10^{–8} 5.170×10^{–8} 5.857×10^{–8}

• For the SSSS plate, the maximum values of*E**z*and*P**z*appear in the center of the plate, and they are much
greater than the minimum values at the edge of the plate.

• For the CCCC plate, the maximum values of*E**z*and*P**z*obtain the maximum values in the center of the plate,
and the minimum values at the edge of the plate; their magnitudes are not much different than in the case of
the*SSSS*plate.

• For the CSCS and CCSS plates, the maximum values of*E**z*and*P**z*do not appear in the center of the plate
and tend to move to the simply supported edges; their magnitudes are not much different than in the case of
the SSSS plate.

• Figures17 and 18 also show clearly that when more than one edge of the plate is fixed, there will be
positions, where the values of*E**z*and*P**z*equal to zero. In the other worlds, there is a position, in which there
is no electric field and polarization. In contrast, for the*SSSS*plate, this phenomenon will not occur. This is
suggested for experimental measurements to determine the electric field or polarization. It is necessary to
determine the measurement location to obtain the desired results.

4.4 Effect of elastic foundation

Consider a piezoelectric nanoplate under the*SSSS*boundary condition (*f*_{14}^{∗} 3). Changing the parameters of
the elastic foundation between 0 and 100, the maximum displacement, stress, electric field, and polarization
in the center of the plate are listed in Table3. It can be seen that when increasing the parameters of the elastic
foundation, the deflection, stress, electric field, and polarization of the plate reduce significantly.

4.5 Effect of the plate thickness

Finally, to see more clearly the maximum deflection of the plate with and without the flexoelectric effect in
the case of changing the plate thickness*h*, this subsection uses the maximum deflection ratio as follows:

*R*_{w} *w*^{∗}max*(*without flexoelectric effect*)*
*w*_{max}^{∗} *(*with flexoelectric effect*)*

*SSSS* *CCCC*

*CSCS* *CCSS*

**Fig. 16** The distributions of*P*_{z}^{∗}and*E*^{∗}_{z}along the edges of the plate with different boundary conditions,*z* −^{h}_{2}