Semi-Analytical Solution for Thermo-Piezoelectric Bending of FG Porous Plates Reinforced with Graphene Platelets

Abstract

This research is devoted to investigating the thermo-piezoelectric bending of functionally graded (FG) porous piezoelectric plates reinforced with graphene platelets (GPLs). A refined four-variable shear deformation plate theory is utilized considering the transverse shear strain to describe the displacement components. The porous nanocomposite plate is composed of polymer piezoelectric material containing internal pores and reinforced with FG GPLs. In accordance with modified distribution laws, the porosity and GPLs volume fraction are presumed to vary continuously through the plate thickness. Four GPLs and porosity distribution types are presented. By applying the Halpin–Tsai model, the elastic properties of the nanocomposite plate are calculated. The governing equations are derived based on the present theory and the principle of virtual work. The deduced partial differential equations are converted to ordinary equations by employing Levy-type solution. These equations are numerically solved based on the differential quadrature method (DQM). In order to determine the minimum grid points sufficient to gain a converging solution, a convergence study is introduced. Moreover, the accuracy of the present formulations are examined by comparing the obtained results with those published in the literature. Additional parametric analyses are introduced to investigate the influences of the GPLs weight fraction, distribution types, side-to-thickness ratio, external electric voltage and temperature on the thermal bending of FG GPLs porous nanocomposite piezoelectric plates.

1. Introduction

Materials that can generate electric charges whenever they are deformed and subjected to a mechanical load and vice versa are called piezoelectric materials [1]. Composite structures with layers of such materials, known as active lightweight smart structures, have drawn considerable attention for variety different applications such as sensors, heat exchangers, automobiles, nuclear devices and transducers. Furthermore, piezoelectric sensors react to vibrations and produce an electric voltage; the resultant voltage can be prepared and intensified by a feedback gain and then applied to an actuator. Because of the converse piezoelectric impact, the actuator will create a control force. Therefore, the piezoelectric materials are extensively utilized for the manufacture of intelligent structures and systems due to their ability to suppress dynamic vibrations and to control shape [2]. Extraordinary composite structures with piezoelectric materials are characterized by electro-mechanical coupling properties, as well as the effective capability to convert energy types such as mechanical and electrical energy between each other [3]. Accordingly, several studies in the open literature have analyzed the behavior of such materials, which makes dealing with the FGMs more efficient, as discussed by Wu et al. [4]. Based on the DQM, Sobhy [3] discussed the axisymmetric bending response of sandwich annular and circular nanoplates with FG GPLs-reinforced face layers and an FG porous core, integrated with piezoelectric layers lying on elastic foundations. By using the finite element method and the Euler–Bernoulli theory, El Harti et al. [5] investigated the vibration control of FG porous beams with bonded piezoelectric materials under a thermal environment, in which the motion equations were obtained through Hamilton’s principle. Mallek et al. [6] presented the nonlinear dynamic behavior of piezolaminated FG carbon nanotube-reinforced composite shells based on the improved first-order shear deformation theory. Zenkour and Aljadani [7] examined the electro-mechanical buckling behavior of simply supported rectangular FG piezoelectric (FGP) plates under the impact of external electric voltage using a quasi-3D refined plate theory. Moreover, Moradi-Dastjerdi and Behdinan [8] analyzed the free vibration behavior of smart sandwich plates with piezoelectric face sheets and porous FG carbon nanotubes, employing Reddy’s third-order shear deformation theory.

Graphene guarantees efficient electro-mechanical and thermal characteristics, such as semi-perfect optical transparency and electrical conductivity [9]. The tensile strength of graphene is about 130.5 GPa, it has a Young’s modulus greater than 1 TPa and its electrical conductivity is 1000 times greater than copper in terms of electric current-carrying capacity. In the basic structure of graphene, carbon atoms are arranged in regular hexagonal pattern as in the case of graphite, but in a one-atom-thick sheet. It is very light, with a 1 m${}^2$ of a single-sheet weighing only 0.77 mg [10]. Moreover, the specific-surface-area of the graphene is around 2630 m${}^2$/g [11], whereas that of carbon nanotubes is in the range of 100–1000 m${}^2$/g. Graphene represents the best reinforcement for polymer-, metal- and ceramic-matrix composite structures to enhance their mechanical properties and piezoelectric characteristics, as well as their stiffness. As a consequence, many articles have investigated the properties and behavior of piezoelectric nanocomposite structures reinforced with graphene platelets (GPLs). Experimentally, Abolhasani et al. [12] prepared piezoelectric PVDF (polyvinylidene fluoride) composite nanofibers reinforced with graphene and investigated the polymorphism, crystallinity, electrical outputs and morphology of these composites. Mao et al. [13] discussed the small-scale effect on the frequencies of graphene nanoplatelet-reinforced FG piezoelectric composite microplates depending on the nonlocal constitutive relation in addition to von Karman geometric non-linearity. Furthermore, Sobhy and Al Mukahal [14] presented the free vibration of piezoelectromagnetic plates reinforced with FG graphene nanosheets (FG-GNSs) subjected to external electric and magnetic potentials. They found that the increase of the elastic foundation stiffness, graphene weight fraction and applied magnetic potential and electromagnetic properties of graphene can enhance plate stiffness. They also [15] studied wave propagation in a sandwich plate with GPLs-reinforced piezoelectromagnetic face layers and a honeycomb core based on higher-order shear deformation plate theory. Alazwari et al. [16] employed the DQM to investigate the critical buckling temperature of piezoelectric circular nanoplates reinforced with uniformly distributed GPLs resting on an elastic substrate and subjected to an external electric field.

Porous structures are attracting remarkable attention as advanced engineering materials in aerospace vehicles, the automotive industry and civil manufacturing due to their outstanding multi-functionality, such as low specific weight, reduced thermal and electrical conductivity, and efficient capacity for energy dissipation. Porous materials are now commercially available for manufacturing lightweight sandwich structures. Composite structures with porosity can resist bending and shear forces, alongside reducing damping vibration; as a result, they have been greatly studied by many authors. Therefore, it is important to consider the impact of porosity on the static and dynamic behavior of FGM plates [17,18]. Sahmani et al. [19] illustrated the nonlinear bending response of FG graphene-reinforced porous nanobeams. In addition, Amir et al. [20] have presented the free vibration of porous three-layered annular and circular microplates reinforced with FG carbon nanotubes. Furthermore, the free vibrational behavior of porous nanocomposite metal foam shells reinforced with GPLs has been demonstrated by Barati and Zenkour [21] considering different porosity distributions, in which their outcomes were performed with the help of first-order shear deformation theory together with Galerkin’s method. In their research, it was obvious that the distribution of porosity did, in fact, greatly affect vibrational frequencies. Furthermore, the hygrothermal buckling work of porous FGM microplates and microbeams was analyzed in an article of Sobhy [22], depending on a new four-variable shear deformation theory and the modified couple stress theory. Additional recent works in the literature have discussed the various behaviors of porous FGGPLs-reinforced composites, such as (Sahmani and Madyira [23], Zhao et al. [24], Ansari et al. [25], Zenkour and Aljadani [26], and Raza et al. [27]).

As shown in previous studies, several investigations about the behavior of porous polymer- or metal-GPLs-reinforced structures have been performed; nevertheless, there is little research on piezoelectric materials reinforced with GPLs. There are no research studies on the behavior of FG porous piezoelectric materials reinforced with GPLs. Therefore, a new investigation depending on the Levy procedure and the DQM is introduced to analyze the thermal bending of FG porous piezoelectric/GPLs plates. Moreover, increasing temperature has significant effects on the behavior of the structures. The increase of temperature negatively affects the mechanical properties of the panels. To reduce the temperature effects, various types of porosity are considered in the present composed plate. A refined four-variable shear deformation plate theory is presented to describe the displacement components. The current plate consists of a piezoelectric material containing internal pores and is supported by FG GPLs. Depending on the Halpin–Tsai model, the effective Young’s modulus of the nanocomposite plate is evaluated. According to the rule of the mixture, Poisson’s ratio, the thermal expansion coefficient and piezoelectric properties are computed for four FG GPLs and porosity-distribution types. In addition, according to modified distribution laws, the volume fraction of graphene and the porosity vary continuously throughout the thickness of the plate. The differential equations are deduced by using the principle of virtual work including thermal loads. Next, the accuracy of the present formulations and theory is validated by comparison with some examples from other research in the literature. Moreover, the influences of the GPLs volume fraction, distribution types, external electric voltage and other parameters on thermal bending are all investigated.

2. Mathematical Formulation

Consider a rectangular plate made of porous piezoelectric material reinforced with GPLs. The geometry of this plate is prescribed by the Cartesian coordinate system $x_1,x_2$ and $x_3$ placed in the mid-plane, with the origin positioned at the plate’s corner. As obvious in Figure 1, the related dimensions are the plate length A, width B and thickness H. Furthermore, the porosity and GPLs are continuously varied. Four different types of porous plates are considered in this paper, as shown in Figure 2, depending on how the interior pores are distributed throughout the thickness of the plate: (a) uniform porosity distribution—the pores in this type are evenly distributed (Type I); (b) porosity distribution Type II—the distribution of the pores is as small as possible at the top and bottom, and then gradually increases until it reaches its highest value at the middle of the plate, (c) porosity distribution Type III—in this type, we have the opposite of Type II; the distribution of the pores is maximum at the top and bottom of the plate, and is minimum in the middle; (d) porosity distribution Type IV—in the last case, the distribution of pores is maximum at the bottom and gradually decreases until it reaches its smallest value at the top. The internal porosities of Types II, III and IV vary in density along the plate thickness.

The effective mechanical and piezoelectric characteristics of porous nanocomposite plates as functions of $x_3$ are given as follows by the closed-cell cellular solids [28,29]:

$$\begin{array}{cc}\hfill E\left(x_3\right)& =\overline{E}[1-e_{0}K\left(x_3\right)],\hfill \\ \hfill \alpha \left(x_3\right)& =\overline{\alpha }[1-e_{1}K\left(x_3\right)],\hfill \\ \hfill {\zeta }_{ij}\left(x_3\right)& ={\overline{\zeta }}_{ij}[1-e_{0}K\left(x_3\right)],ij=1,2,3,4,5,\hfill \\ \hfill f_{ii}\left(x_3\right)& ={\overline{f}}_{ii}[1-e_{0}K\left(x_3\right)],i=1,2,3,\hfill \\ \hfill \nu \left(x_3\right)& =\overline{\nu }\left[1-1.21\check{P}+0.342{\check{P}}^2\right]+0.221\check{P},\hfill \end{array}$$

(1)

in which $\overline{E}=E_{max}$, $\overline{\alpha }={\alpha }_{max}$, ${\overline{\zeta }}_{ij}={\zeta }_{ij}^{max}$, ${\overline{f}}_{ii}=f_{ii}^{max}$ and $\overline{\nu }={\nu }_{max}$ denote Young’s modulus, the thermal expansion coefficient, the piezoelectric coefficients, the dielectric coefficients and Poisson’s ratio of the nanocomposite structures without porosity, respectively. Furthermore, according to Kitipornchai et al. [29], the mass density coefficient $e_1$ is presented in terms of the porosity coefficient $e_0$ as:

$$\begin{array}{c}\hfill e_1=\frac{\check{P}\left(x_3\right)}{K\left(x_3\right)},\check{P}\left(x_3\right)=1.121\{1-{[1-e_{0}K\left(x_3\right)]}^{\frac{1}{2.3}}\},\end{array}$$

(2)

and

$$K\left(x_3\right)=\left\{\begin{array}{cc}1\hfill & \mathrm{for}\mathrm{Type}\mathrm{I},\hfill \\ 1-{\left(2\frac{x_3}{H}\right)}^2\hfill & \mathrm{for}\mathrm{Type}\mathrm{II},\hfill \\ {\left(2\frac{x_3}{H}\right)}^2\hfill & \mathrm{for}\mathrm{Type}\mathrm{III},\hfill \\ \left(\frac{1}{2}-\frac{x_3}{H}\right)\hfill & \mathrm{for}\mathrm{Type}\mathrm{IV}.\hfill \end{array}\right.$$

(3)

The shear modulus of the plates with pores is

$$\begin{array}{c}\hfill G\left(x_3\right)=\frac{E\left(x_3\right)}{2[1+\nu \left(x_3\right)]}.\end{array}$$

(4)

Next, using the Halpin–Tsai micromechanics model and the mixture rule [30,31,32,33], it is possible to predict the effective material properties of the FG GPLs-reinforced plates:

$$\begin{array}{cc}\hfill \overline{E}& =\frac{E_m}{8}\left(3\frac{{\eta }_1+2{\xi }_L^G{v}_G}{{\eta }_1-v_G}+5\frac{{\eta }_2+2{\xi }_b^G{v}_G}{{\eta }_2-v_G}\right),\hfill \\ \hfill \overline{\rho }& =v_G{\rho }_G+(1-v_G){\rho }_m,\hfill \\ \hfill \overline{\nu }& =v_G{\nu }_G+(1-v_G){\nu }_m,\hfill \\ \hfill \overline{\alpha }& =v_G{\alpha }_G+(1-v_G){\alpha }_m,\hfill \\ \hfill {\overline{\zeta }}_{ij}& =v_G{\zeta }_{ij}^G+(1-v_G){\zeta }_{ij}^m,ij=1,2,3,4,5,\hfill \\ \hfill {\overline{f}}_{ii}& =v_G{f}_{ii}^G+(1-v_G)f_{ii}^m,i=1,2,3,\hfill \end{array}$$

(5)

where $E_G\left(E_m\right)$, ${\rho }_G\left({\rho }_m\right)$, ${\nu }_G\left({\nu }_m\right)$ and ${\alpha }_G\left({\alpha }_m\right)$ are Young’s modulus, mass density, Poisson’s ratio and thermal expansion coefficient of the GPLs (piezoelectric material), respectively. Moreover, the remaining properties, ${\zeta }_{ij}^G\left({\zeta }_{ij}^m\right)$ and $f_{ii}^G\left(f_{ii}^m\right)$ denote the piezoelectric coefficients and the dielectric of the GPLs (piezoelectric material), respectively. However, the values for the parameters ${\eta }_1,{\eta }_2,{\xi }_L^G$ and ${\xi }_b^G$ are provided as [29]:

$$\begin{array}{c}\hfill {\eta }_1=\frac{E_G+2{\xi }_L^G{E}_m}{E_G-E_m},{\xi }_L^G=\frac{L^G}{h^G},\\ \hfill {\eta }_2=\frac{E_G+2{\xi }_b^G{E}_m}{E_G-E_m},{\xi }_b^G=\frac{b^G}{h^G},\end{array}$$

(6)

in which $L^G,b^G$ and $h^G$ denote the length, width and thickness of the GPLs, respectively.

In addition, there are four different GPLs distribution patterns, created by altering the volume fraction of the GPLs through the thickness of the plates according to a cosine distribution rule. The GPLs are evenly dispersed throughout the thickness for Pattern A. In contrast, the GPLs for Patterns B, C and D are functionally graded through the thickness, as shown in Figure 3. As a result, the volume fraction $v_G\left(x_3\right)$ of GPLs can be expressed as follows [29]:

$$v_G\left(x_3\right)=\left\{\begin{array}{cc}{v}^1\hfill & \mathrm{for}\mathrm{Pattern}\mathrm{A},\hfill \\ v^{2}cos\left(\frac{\pi x_3}{H}\right)\hfill & \mathrm{for}\mathrm{Pattern}\mathrm{B},\hfill \\ v^3\left[1-cos\left(\frac{\pi x_3}{H}\right)\right]\hfill & \mathrm{for}\mathrm{Pattern}\mathrm{C},\hfill \\ v^4\left[1-cos(\frac{\pi x_3}{2H}+\frac{\pi }{4})\right]\hfill & \mathrm{for}\mathrm{Pattern}\mathrm{D},\hfill \end{array}\right.$$

(7)

where $v^i$ are defined as the top values of GPLs volume fraction, which can be calculated as [29]:

$${\int }_{-H/2}^{H/2}{v}_G\left(x_3\right)(1-\check{P})\mathrm{d}{x}_3=\frac{W_f}{W_f+{\rho }_G/{\rho }_m-W_f{\rho }_G/{\rho }_m}{\int }_{-H/2}^{H/2}(1-\check{P})\mathrm{d}{x}_3,$$

(8)

where $W_f$ is the GPLs weight fraction.

3. Constitutive Equations

The two-variable shear deformation plate theory was first established for isotropic panels by Shimpi [34]. This theory has been extended for FG structures by many authors [35,36], and contains four variables and various shape functions. This theory considers the transverse shear effects and parabolic distribution of the transverse shear stresses through the thickness. Since the transverse shear strains and stresses are correctly represented, it is unnecessary to use shear correction factors. In addition, this theory satisfies the traction-free boundary conditions at the panel faces. Accordingly, a modified refined four-variable shear deformation theory [14,34,37] for FG porous GPLs-reinforced plates is presented to describe the displacement components as follows:

$$\begin{array}{cc}\hfill U(x_1,x_2,x_3)& =u_0(x_1,x_2)-x_3\frac{\partial w_1}{\partial x_1}-f\left(x_3\right)\frac{\partial w_2}{\partial x_1},\hfill \\ \hfill V(x_1,x_2,x_3)& =v_0(x_1,x_2)-x_3\frac{\partial w_1}{\partial x_2}-f\left(x_3\right)\frac{\partial w_2}{\partial x_2},\hfill \\ \hfill W(x_1,x_2,x_3)& =w_1(x_1,x_2)+w_2(x_1,x_2),\hfill \end{array}$$

(9)

where $u_0(x_1,x_2)$ and $v_0(x_1,x_2)$ represent the components of the mid-plane displacements along the $x_1$- and $x_2$-axes, respectively. The deflection $W(x_1,x_2,x_3)$ is separated into two components, $w_1(x_1,x_2)$ and $w_2(x_1,x_2)$, which specify, respectively, the bending and the shear displacements, according to Shimpi’s theory [34] and $f\left(x_3\right)=x_3-F\left(x_3\right)$. In addition, the shape function that describes the configuration of the shear stress along the thickness of the plate can be demonstrated by [15,34]:

$$\begin{array}{c}\hfill F\left(x_3\right)=\frac{x_3}{1+{(x_3/H)}^2}-\frac{5}{8}\frac{x_3^3}{H^2}.\end{array}$$

(10)

The components of the strains can be calculated using the displacements field (9) as follows:

$$\begin{array}{cc}\hfill {\epsilon }_{11}& =\frac{\partial U}{\partial x_1}=\frac{\partial u_0}{\partial x_1}-x_3\frac{{\partial }^2{w}_1}{\partial x_1^2}-f\left(x_3\right)\frac{{\partial }^2{w}_2}{\partial x_1^2},\hfill \\ \hfill {\epsilon }_{22}& =\frac{\partial V}{\partial x_2}=\frac{\partial v_0}{\partial x_2}-x_3\frac{{\partial }^2{w}_1}{\partial x_2^2}-f\left(x_3\right)\frac{{\partial }^2{w}_2}{\partial x_2^2},\hfill \\ \hfill {\epsilon }_{33}& =\frac{\partial W}{\partial x_3}=0,\hfill \\ \hfill {\epsilon }_{12}& =\frac{\partial U}{\partial x_2}+\frac{\partial V}{\partial x_1}=\frac{\partial u_0}{\partial x_2}+\frac{\partial v_0}{\partial x_1}-2x_3\frac{{\partial }^2{w}_1}{\partial x_1\partial x_2}-2f\left(x_3\right)\frac{{\partial }^2{w}_2}{\partial x_1\partial x_2},\hfill \\ \hfill {\epsilon }_{13}& =\frac{\partial U}{\partial x_3}+\frac{\partial W}{\partial x_1}=F^{\prime }\left(x_3\right)\frac{\partial w_2}{\partial x_1},\hfill \\ \hfill {\epsilon }_{23}& =\frac{\partial V}{\partial x_3}+\frac{\partial W}{\partial x_2}=F^{\prime }\left(x_3\right)\frac{\partial w_2}{\partial x_2}.\hfill \end{array}$$

(11)

In terms of the piezoelasticity theory [3,38,39], the constitutive relationship for the components of the stresses $\sigma$ can be described as follows:

$$\left\{\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{23}\\ {\sigma }_{13}\\ {\sigma }_{12}\end{array}\right\}=\left[\begin{array}{ccccc}{R}_{11}& R_{12}& 0& 0& 0\\ R_{12}& R_{22}& 0& 0& 0\\ 0& 0& R_{44}& 0& 0\\ 0& 0& 0& R_{55}& 0\\ 0& 0& 0& 0& R_{66}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{11}-\alpha T\\ {\epsilon }_{22}-\alpha T\\ {\epsilon }_{23}\\ {\epsilon }_{13}\\ {\epsilon }_{12}\end{array}\right\}-\left[\begin{array}{ccc}0& 0& {\zeta }_{31}\\ 0& 0& {\zeta }_{32}\\ 0& {\zeta }_{24}& 0\\ {\zeta }_{15}& 0& 0\\ 0& 0& 0\end{array}\right]\left\{\begin{array}{c}{\widehat{E}}_1\\ {\widehat{E}}_2\\ {\widehat{E}}_3\end{array}\right\},$$

(12)

where $(R_{11},R_{12},R_{22},R_{44},R_{55},R_{66})$ are the composite plate’s elastic coefficients, which are written as:

$$R_{11}=R_{22}=\frac{E}{1-{\nu }^2},R_{12}=\frac{E\nu }{1-{\nu }^2},R_{44}=R_{55}=R_{66}=\frac{E}{2(1+\nu )},$$

(13)

and $T(x_1,x_2,x_3)$ denotes the applied temperature, which is supposed to be vary along the thickness of the plate according to the following law [32]:

$$\begin{array}{c}\hfill T(x_1,x_2,x_3)={\left(\frac{x_3}{H}+\frac{1}{2}\right)}^k\overline{T}(x_1,x_2),k⩾0,\end{array}$$

(14)

where k indicates the temperature exponent. Furthermore, the electric displacements $D_i$ can be written as [38]:

$$\left\{\begin{array}{c}{D}_1\\ D_2\\ D_3\end{array}\right\}=\left[\begin{array}{ccccc}0& 0& 0& {\zeta }_{15}& 0\\ 0& 0& {\zeta }_{24}& 0& 0\\ {\zeta }_{31}& {\zeta }_{32}& 0& 0& 0\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{11}-\alpha T\\ {\epsilon }_{22}-\alpha T\\ {\epsilon }_{23}\\ {\epsilon }_{13}\\ {\epsilon }_{12}\end{array}\right\}+\left[\begin{array}{ccc}{f}_{11}& 0& 0\\ 0& f_{22}& 0\\ 0& 0& f_{33}\end{array}\right]\left\{\begin{array}{c}{\widehat{E}}_1\\ {\widehat{E}}_2\\ {\widehat{E}}_3\end{array}\right\}.$$

(15)

According to [38], the electric field $\widehat{E}$ can be written as:

$$\begin{array}{c}\hfill \widehat{E}=-\nabla \mathsf{\Psi },\end{array}$$

(16)

where $\mathsf{\Psi }$ represents the electric potential of the FG-GPLs/piezoelectric plate, which can be described as [38]:

$$\begin{array}{c}\hfill \mathsf{\Psi }(x_1,x_2,x_3)=-\psi (x_1,x_2)cos\left(\frac{\pi x_3}{H}\right)+\frac{2x_3{\overline{V}}_0}{H},\end{array}$$

(17)

in which ${\overline{V}}_0$ is the external applied electric potential and $\psi (x_1,x_2)$ denotes the electric potential of the middle surface of the plate. Equation (17) is substituted into Equation (16) to provide the components of the electric field as follows:

$$\left\{\begin{array}{c}{\widehat{E}}_1\\ {\widehat{E}}_2\\ {\widehat{E}}_3\end{array}\right\}=\left\{\begin{array}{c}{\displaystyle \frac{\partial \psi }{\partial x_1}cos\left(\frac{\pi x_3}{H}\right)}\\ {\displaystyle \frac{\partial \psi }{\partial x_2}cos\left(\frac{\pi x_3}{H}\right)}\\ {\displaystyle -\frac{\pi }{H}\psi sin\left(\frac{\pi x_3}{H}\right)}\end{array}\right\}-\left\{\begin{array}{c}0\\ 0\\ {\displaystyle \frac{2{\overline{V}}_0}{H}}\end{array}\right\}.$$

(18)

4. Governing Equations

To conclude the governing equations, we will use the principal of virtual work, which can be expressed as:

$$\begin{array}{c}\hfill \delta {\mathcal{W}}_S-\delta {\mathcal{W}}_E-\delta {\mathcal{W}}_F=0,\end{array}$$

(19)

where the variation of the strain energy $\delta {\mathcal{W}}_S$ and the variation of the work done by the in-plane piezoelectric load $\delta {\mathcal{W}}_E$ and the external loads $\delta {\mathcal{W}}_F$ are given as follows:

$$\begin{array}{cc}\hfill \delta {\mathcal{W}}_S& ={\int }_V\left({\sigma }_{jk}\delta {\epsilon }_{jk}-D_j\delta {\widehat{E}}_j\right)\mathrm{d}V,\hfill \\ \hfill \delta {\mathcal{W}}_E& ={\int }_A{N}^E\left(\frac{{\partial }^{2}W}{\partial x_1^2}+\frac{{\partial }^{2}W}{\partial x_2^2}\right)\delta W\mathrm{d}A,\hfill \\ \hfill \delta {\mathcal{W}}_F& ={\int }_{A}q\delta W\mathrm{d}A,\hfill \end{array}$$

(20)

in which q denotes the vertical load applied to the plate and $N^E$ denotes the in-plane piezoelectric load, which is defined as [38]:

$$\begin{array}{c}\hfill N^E={\int }_{-H/2}^{H/2}{\zeta }_{31}\frac{2{\overline{V}}_0}{H}\mathrm{d}{x}_3.\end{array}$$

(21)

Inserting the equations of ${\epsilon }_{jk}$ and ${\widehat{E}}_j$ into the first equation of Equation (20) leads to

$$\begin{array}{cc}\hfill \delta {\mathcal{W}}_S=& {\int }_A[{\hat{C}}_{11}\delta \frac{\partial u_0}{\partial x_1}-{\widehat{L}}_{11}\delta \frac{{\partial }^2{w}_1}{\partial x_1^2}-{\widehat{O}}_{11}\delta \frac{{\partial }^2{w}_2}{\partial x_1^2}+{\widehat{C}}_{22}\delta \frac{\partial v_0}{\partial x_2}-{\widehat{L}}_{22}\delta \frac{{\partial }^2{w}_1}{\partial x_2^2}-{\widehat{O}}_{22}\delta \frac{{\partial }^2{w}_2}{\partial x_2^2}\hfill \\ \hfill & +{\widehat{C}}_{12}\delta \left(\frac{\partial u_0}{\partial x_2}+\frac{\partial v_0}{\partial x_1}\right)-2{\widehat{L}}_{12}\delta \frac{{\partial }^2{w}_1}{\partial x_1\partial x_2}-2{\widehat{O}}_{12}\delta \frac{{\partial }^2{w}_2}{\partial x_1\partial x_2}+Q_{13}\delta \frac{\partial w_2}{\partial x_1}+Q_{23}\delta \frac{\partial w_2}{\partial x_2}\hfill \\ \hfill & -S_1\delta \frac{\partial \psi }{\partial x_1}-S_2\delta \frac{\partial \psi }{\partial x_2}+S_3\delta \psi ]\mathrm{d}A,\hfill \end{array}$$

(22)

where

$$\begin{array}{cc}\hfill \{{\widehat{C}}_{11},{\widehat{C}}_{22},{\widehat{C}}_{12}\}& ={\int }_{-H/2}^{H/2}\{{\sigma }_{11},{\sigma }_{22},{\sigma }_{12}\}\mathrm{d}{x}_3,\hfill \\ \hfill \{{\widehat{L}}_{11},{\widehat{L}}_{22},{\widehat{L}}_{12}\}& ={\int }_{-H/2}^{H/2}{x}_3\{{\sigma }_{11},{\sigma }_{22},{\sigma }_{12}\}\mathrm{d}{x}_3,\hfill \\ \hfill \{{\widehat{O}}_{11},{\widehat{O}}_{22},{\widehat{O}}_{12}\}& ={\int }_{-H/2}^{H/2}f\left(x_3\right)\{{\sigma }_{11},{\sigma }_{22},{\sigma }_{12}\}\mathrm{d}{x}_3,\hfill \\ \hfill \left\{Q_{13},Q_{23}\right\}& ={\int }_{-H/2}^{H/2}{F}^{\prime }\left(x_3\right)\{{\sigma }_{13},{\sigma }_{23}\}\mathrm{d}{x}_3,\hfill \\ \hfill S_1& ={\int }_{-H/2}^{H/2}{D}_{1}cos\left(\frac{\pi x_3}{H}\right)\mathrm{d}{x}_3,\hfill \\ \hfill S_2& ={\int }_{-H/2}^{H/2}{D}_{2}cos\left(\frac{\pi x_3}{H}\right)\mathrm{d}{x}_3,\hfill \\ \hfill S_3& ={\int }_{-H/2}^{H/2}{D}_3\frac{\pi }{H}sin\left(\frac{\pi x_3}{H}\right)\mathrm{d}{x}_3.\hfill \end{array}$$

(23)

By substituting ${\sigma }_{jk}$ and $D_i$ into Equation (23), one obtains the following:

$$\begin{array}{cc}\hfill {\widehat{C}}_{11}& =N_{11}\frac{\partial u_0}{\partial x_1}-N_{12}\frac{{\partial }^2{w}_1}{\partial x_1^2}-N_{13}\frac{{\partial }^2{w}_2}{\partial x_1^2}+Z_{11}\frac{\partial v_0}{\partial x_2}-Z_{12}\frac{{\partial }^2{w}_1}{\partial x_2^2}-Z_{13}\frac{{\partial }^2{w}_2}{\partial x_2^2}\hfill \\ \hfill & -C_{11}^T+B_1\psi +C_{11}^E,\hfill \\ \hfill {\widehat{C}}_{22}& =Z_{11}\frac{\partial u_0}{\partial x_1}-Z_{12}\frac{{\partial }^2{w}_1}{\partial x_1^2}-Z_{13}\frac{{\partial }^2{w}_2}{\partial x_1^2}+N_{11}\frac{\partial v_0}{\partial x_2}-N_{12}\frac{{\partial }^2{w}_1}{\partial x_2^2}-N_{13}\frac{{\partial }^2{w}_2}{\partial x_2^2}\hfill \\ \hfill & -C_{11}^T+B_2\psi +C_{22}^E,\hfill \\ \hfill {\widehat{C}}_{12}& =A_{11}\left(\frac{\partial u_0}{\partial x_2}+\frac{\partial v_0}{\partial x_1}\right)-2A_{12}\frac{{\partial }^2{w}_1}{\partial x_1\partial x_2}-2A_{13}\frac{{\partial }^2{w}_2}{\partial x_1\partial x_2},\hfill \\ \hfill {\widehat{L}}_{11}& =N_{12}\frac{\partial u_0}{\partial x_1}-N_{22}\frac{{\partial }^2{w}_1}{\partial x_1^2}-N_{23}\frac{{\partial }^2{w}_2}{\partial x_1^2}+Z_{12}\frac{\partial v_0}{\partial x_2}-Z_{22}\frac{{\partial }^2{w}_1}{\partial x_2^2}-Z_{23}\frac{{\partial }^2{w}_2}{\partial x_2^2}\hfill \\ \hfill & -L_{11}^T+B_3\psi +L_{11}^E,\hfill \\ \hfill {\widehat{L}}_{22}& =Z_{12}\frac{\partial u_0}{\partial x_1}-Z_{22}\frac{{\partial }^2{w}_1}{\partial x_1^2}-Z_{23}\frac{{\partial }^2{w}_2}{\partial x_1^2}+N_{12}\frac{\partial v_0}{\partial x_2}-N_{22}\frac{{\partial }^2{w}_1}{\partial x_2^2}-N_{23}\frac{{\partial }^2{w}_2}{\partial x_2^2}\hfill \\ \hfill & -L_{11}^T+B_4\psi +L_{22}^E,\hfill \\ \hfill {\widehat{L}}_{12}& =A_{12}\left(\frac{\partial u_0}{\partial x_2}+\frac{\partial v_0}{\partial x_1}\right)-2A_{22}\frac{{\partial }^2{w}_1}{\partial x_1\partial x_2}-2A_{23}\frac{{\partial }^2{w}_2}{\partial x_1\partial x_2},\hfill \\ \hfill {\widehat{O}}_{11}& =N_{13}\frac{\partial u_0}{\partial x_1}-N_{23}\frac{{\partial }^2{w}_1}{\partial x_1^2}-N_{33}\frac{{\partial }^2{w}_2}{\partial x_1^2}+Z_{13}\frac{\partial v_0}{\partial x_2}-Z_{23}\frac{{\partial }^2{w}_1}{\partial x_2^2}-Z_{33}\frac{{\partial }^2{w}_2}{\partial x_2^2}\hfill \\ \hfill & -O_{11}^T+B_5\psi +O_{11}^E,\hfill \\ \hfill {\widehat{O}}_{22}& =Z_{13}\frac{\partial u_0}{\partial x_1}-Z_{23}\frac{{\partial }^2{w}_1}{\partial x_1^2}-Z_{33}\frac{{\partial }^2{w}_2}{\partial x_1^2}+N_{13}\frac{\partial v_0}{\partial x_2}-N_{23}\frac{{\partial }^2{w}_1}{\partial x_2^2}-N_{33}\frac{{\partial }^2{w}_2}{\partial x_2^2}\hfill \\ \hfill & -O_{11}^T+B_6\psi +O_{22}^E,\hfill \\ \hfill {\widehat{O}}_{12}& =A_{13}\left(\frac{\partial u_{\circ }}{\partial x_2}+\frac{\partial v_{\circ }}{\partial x_1}\right)-2A_{23}\frac{{\partial }^2{w}_1}{\partial x_1\partial x_2}-2A_{33}\frac{{\partial }^2{w}_2}{\partial x_1\partial x_2},\hfill \\ \hfill Q_{13}& =K_{13}\frac{\partial w_2}{\partial x_1}-Q_{13}^E\frac{\partial \psi }{\partial x_1},Q_{23}=K_{13}\frac{\partial w_2}{\partial x_2}-Q_{23}^E\frac{\partial \psi }{\partial x_2},\hfill \\ \hfill S_1& =Q_{13}^E\frac{\partial w_2}{\partial x_1}+S_{12}^E\frac{\partial \psi }{\partial x_1},S_2=Q_{23}^E\frac{\partial w_2}{\partial x_2}+S_{22}^E\frac{\partial \psi }{\partial x_2},\hfill \\ \hfill S_3& =B_1\frac{\partial u_0}{\partial x_1}-B_3\frac{{\partial }^2{w}_1}{\partial x_1^2}-B_5\frac{{\partial }^2{w}_2}{\partial x_1^2}+B_2\frac{\partial v_0}{\partial x_2}-B_4\frac{{\partial }^2{w}_1}{\partial x_2^2}-B_6\frac{{\partial }^2{w}_2}{\partial x_2^2}\hfill \\ \hfill & -S_3^T-B_7\psi -S_3^E,\hfill \end{array}$$

(24)

where

$$\begin{array}{cc}\hfill \left\{N_{11},N_{12},N_{13},N_{22},N_{23},N_{33}\right\}=& {\int }_{-H/2}^{H/2}{R}_{11}\left\{1,x_3,f,x_3^2,f{x}_3,f^2\right\}\mathrm{d}{x}_3,\hfill \\ \hfill \left\{Z_{11},Z_{12},Z_{13},Z_{22},Z_{23},Z_{33}\right\}=& {\int }_{-H/2}^{H/2}{R}_{12}\left\{1,x_3,f,x_3^2,f{x}_3,f^2\right\}\mathrm{d}{x}_3,\hfill \\ \hfill \left\{A_{11},A_{12},A_{13},A_{22},A_{23},A_{33}\right\}=& {\int }_{-H/2}^{H/2}{R}_{44}\left\{1,x_3,f,x_3^2,f{x}_3,f^2\right\}\mathrm{d}{x}_3,\hfill \\ \hfill \{C_{11}^T,L_{11}^T,O_{11}^T\}=& {\int }_{-H/2}^{H/2}\alpha (R_{11}+R_{12})T\left\{1,x_3,f\right\}\mathrm{d}{x}_3,\hfill \\ \hfill \left\{B_1,B_2,B_3,B_4,B_5,B_6,B_7\right\}=& {\int }_{-H/2}^{H/2}\frac{\pi }{H}sin\left(\frac{\pi x_3}{H}\right)\{{\zeta }_{31},{\zeta }_{32},{\zeta }_{31}{x}_3,{\zeta }_{32}{x}_3,\hfill \\ \hfill & {\zeta }_{31}f,{\zeta }_{32}f,f_{33}\frac{\pi }{H}sin\left(\frac{\pi x_3}{H}\right)\}\mathrm{d}{x}_3,\hfill \\ \hfill \left\{C_{11}^E,C_{22}^E,L_{11}^E,L_{22}^E,O_{11}^E,O_{22}^E\right\}=& {\int }_{-H/2}^{H/2}\frac{2{\overline{V}}_0}{H}\left\{{\zeta }_{31},{\zeta }_{32},{\zeta }_{31}{x}_3,{\zeta }_{32}{x}_3,{\zeta }_{31}f,{\zeta }_{32}f\right\}\mathrm{d}{x}_3,\hfill \\ \hfill \left\{K_{13},Q_{13}^E,Q_{23}^E\right\}=& {\int }_{-H/2}^{H/2}{F}^{\prime }\left\{F^{\prime }{R}_{44},{\zeta }_{15}cos\left(\frac{\pi x_3}{H}\right),{\zeta }_{24}cos\left(\frac{\pi x_3}{H}\right)\right\}\mathrm{d}{x}_3,\hfill \\ \hfill \left\{S_{12}^E,S_{22}^E\right\}=& {\int }_{-H/2}^{H/2}{cos}^2\left(\frac{\pi x_3}{H}\right)\left\{f_{11},f_{22}\right\},\hfill \\ \hfill S_3^T=& {\int }_{-H/2}^{H/2}\alpha \frac{\pi }{H}({\zeta }_{31}+{\zeta }_{32})Tsin\left(\frac{\pi x_3}{H}\right)\mathrm{d}{x}_3,\hfill \\ \hfill S_3^E=& {\int }_{-H/2}^{H/2}{f}_{33}\frac{\pi }{H}sin\left(\frac{\pi x_3}{H}\right)\frac{2{\overline{V}}_0}{H}\mathrm{d}{x}_3.\hfill \end{array}$$

(25)

By inserting Equations (20) and (22) into Equation (19) and applying integration by parts, one obtains the governing equations as follows:

$$\begin{array}{cc}\hfill & \frac{\partial {\widehat{C}}_{11}}{\partial x_1}+\frac{\partial {\widehat{C}}_{12}}{\partial x_2}=0,\hfill \\ \hfill & \frac{\partial {\widehat{C}}_{12}}{\partial x_1}+\frac{\partial {\widehat{C}}_{22}}{\partial x_2}=0,\hfill \\ \hfill & \frac{{\partial }^2{\widehat{L}}_{11}}{\partial x_1^2}+2\frac{{\partial }^2{\widehat{L}}_{12}}{\partial x_1\partial x_2}+\frac{{\partial }^2{\widehat{L}}_{22}}{\partial x_2^2}+q+N^E\left(\frac{{\partial }^2(w_1+w_2)}{\partial x_1^2}\right)+N^E\left(\frac{{\partial }^2(w_1+w_2)}{\partial x_2^2}\right)=0,\hfill \\ \hfill & \frac{{\partial }^2{\widehat{O}}_{11}}{\partial x_1^2}+2\frac{{\partial }^2{\widehat{O}}_{12}}{\partial x_1\partial x_2}+\frac{{\partial }^2{\widehat{O}}_{22}}{\partial x_2^2}+\frac{\partial Q_{13}}{\partial x_1}+\frac{\partial Q_{23}}{\partial x_2}+q+N^E\left(\frac{{\partial }^2(w_1+w_2)}{\partial x_1^2}\right)+N^E\left(\frac{{\partial }^2(w_1+w_2)}{\partial x_2^2}\right)=0,\hfill \\ \hfill & \frac{\partial S_1}{\partial x_1}+\frac{\partial S_2}{\partial x_2}+S_3=0.\hfill \end{array}$$

(26)

5. Solution Procedure

5.1. Levy-Type Solution

In this section, the Levy-type approach [32,40] is employed here with respect to the variety of boundary conditions to reduce governing Equation (26) to a system of ordinary differential equations that can be solved exactly or approximately. However, Levy-type solution can be developed for rectangular nanoplates with at least two opposite edges that are simply supported, while the remaining two edges can have any possible combination of boundary conditions: free, simply supported or clamped. In accordance with this approach, the edges parallel to the $x_2$-axis are assumed to be simply supported, given as:

$$v_0=w_1=w_2={\widehat{C}}_{11}={\widehat{L}}_{11}=0,\mathrm{at}{x}_1=0,A.$$

(27)

The boundary conditions of the edges parallel to the $x_1$-axis can have a combination of simply supported, clamped or free conditions, which can be expressed as [32]:

Simply supported (S):

$$u_0=w_1=w_2={\widehat{C}}_{22}={\widehat{L}}_{22}=0,\mathrm{at}{x}_2=0,B.$$

(28)

Clamped (C):

$$u_0=v_0=w_1=w_2=\psi =0,\mathrm{at}{x}_2=0,B.$$

(29)

Free (F):

$${\widehat{C}}_{22}={\widehat{L}}_{22}=2\frac{\partial {\widehat{L}}_{12}}{\partial x_1}+\frac{\partial {\widehat{L}}_{22}}{\partial x_2}=2\frac{\partial {\widehat{O}}_{12}}{\partial x_1}+\frac{{\partial }^2{\widehat{O}}_{22}}{\partial x_2}+Q_{23}=0,\mathrm{at}{x}_2=0,B.$$

(30)

The simply supported boundary conditions can be fulfilled by trigonometric functions in $x_1$, and the resulting ordinary differential equations in $x_2$ can be solved by applying the state-space method or the differential quadrature method. According to the Levy solution process, we consider the following representation of the displacements and loading as follows [32]:

$$\left\{\begin{array}{c}{v}_0(x_1,x_2)\\ w_1(x_1,x_2)\\ w_2(x_1,x_2)\\ \psi (x_1,x_2)\\ q(x_1,x_2)\\ \overline{T}(x_1,x_2)\end{array}\right\}=\sum _{n=1}^{\infty }\left\{\begin{array}{c}{v}_n\left(x_2\right)\\ w_{1n}\left(x_2\right)\\ w_{2n}\left(x_2\right)\\ {\psi }_n\left(x_2\right)\\ {\widehat{q}}_n\left(x_2\right)\\ {\widehat{T}}_n\left(x_2\right)\end{array}\right\}sin\left(\mu x_1\right),u_0(x_1,x_2)=\sum _{n=1}^{\infty }cos\left(\mu x_1\right)u_n\left(x_2\right)$$

(31)

where $\mu =n\pi /A$. For the sinusoidal load, the functions ${\widehat{q}}_n\left(x_2\right)$ and ${\overline{T}}_n\left(x_2\right)$ are expressed as:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\widehat{q}}_n\left(x_2\right)\\ {\widehat{T}}_n\left(x_2\right)\end{array}\right\}=\left\{\begin{array}{c}{q}_0\\ T_0\end{array}\right\}sin\left(\frac{\pi x_2}{B}\right).\end{array}$$

(32)

Incorporating Equation (31) into Equation (26) subject to Equation (24) provides the following ordinary differential equations

$$\begin{array}{cc}\hfill & e_{14}{u}_n+e_{12}{w}_{1n}+e_{13}{w}_{2n}+e_{17}{\psi }_n+e_{18}{\widehat{T}}_n\left(x_2\right)+e_{11}{v}_n^{\prime }+e_{21}{u}_n^{″}+(e_{15}+e_{19})w_{1n}^{″}+(e_{16}+e_{20})w_{2n}^{″}=0,\hfill \\ \hfill & e_{22}{v}_n+e_{66}{\widehat{T}}_n^{\prime }\left(x_2\right)+e_{23}{u}_n^{\prime }+e_{24}{w}_{1n}^{\prime }+e_{25}{w}_{2n}^{\prime }+e_{26}{\psi }_n^{\prime }+e_{27}{v}_n^{″}+e_{28}{w}_{1n}^{‴}+e_{29}{w}_{2n}^{‴}=0,\hfill \\ \hfill & e_{30}{u}_n+e_{31}{w}_{1n}+e_{32}{w}_{2n}+e_{33}{\psi }_n+e_{38}{\widehat{T}}_n\left(x_2\right)+e_{34}{v}_n^{\prime }+e_{35}{u}_n^{″}+e_{36}{w}_{1n}^{″}+e_{37}{w}_{2n}^{″}+e_{39}{\psi }_n^{″}+e_{40}{v}_n^{‴}\hfill \\ \hfill & +e_{41}{w}_{1n}^{⁗}+e_{42}{w}_{2n}^{⁗}+e_{422}{\widehat{T}}_n^{″}\left(x_2\right)+{\widehat{q}}_n\left(x_2\right)=0,\hfill \\ \hfill & e_{43}{u}_n+e_{44}{w}_{1n}+e_{45}{w}_{2n}+e_{46}{\psi }_n+e_{55}{\widehat{T}}_n\left(x_2\right)+e_{47}{v}_n^{\prime }+e_{48}{u}_n^{″}+e_{49}{w}_{1n}^{″}+e_{50}{w}_{2n}^{″}+e_{54}{\psi }_n^{″}+e_{65}{\widehat{T}}_n^{″}\left(x_2\right)\hfill \\ \hfill & +e_{53}{v}_n^{‴}+e_{51}{w}_{1n}^{⁗}+e_{52}{w}_{2n}^{⁗}+{\widehat{q}}_n\left(x_2\right)=0,\hfill \\ \hfill & e_{56}{u}_n+e_{57}{w}_{1n}+e_{58}{w}_{2n}+e_{59}{\psi }_n+e_{60}{v}_n^{\prime }+e_{61}{w}_{1n}^{″}+e_{62}{w}_{2n}^{″}+e_{63}{\psi }_n^{″}+e_{64}{\widehat{T}}_n\left(x_2\right)+e_{67}=0,\hfill \end{array}$$

(33)

where the derivation is denoted by ${\left(\right)}^{{}^{\prime }}=\mathrm{d}\left(\right)/\mathrm{d}{x}_2$ and

$$\begin{array}{cc}\hfill & e_{11}=\mu (A_{11}+Z_{11}),e_{12}=N_{12}{\mu }^3,e_{13}=N_{13}{\mu }^3,e_{14}=-N_{11}{\mu }^2,e_{15}=-2A_{12}\mu ,\hfill \\ \hfill & e_{16}=-2A_{13}\mu ,e_{17}=B_1\mu ,e_{18}=-C_{11}^T\mu ,e_{19}=-Z_{12}\mu ,e_{20}=-Z_{13}\mu ,e_{21}=A_{11},\hfill \\ \hfill & e_{22}=-A_{11}{\mu }^2,e_{23}=-\mu (A_{11}+Z_{11}),e_{24}=(2A_{12}+Z_{12}){\mu }^2,e_{25}=(2A_{13}+Z_{13}){\mu }^2,e_{26}=B_2,e_{27}=N_{11},\hfill \\ \hfill & e_{28}=-N_{12},e_{29}=-N_{13},e_{30}=N_{12}{\mu }^3,e_{31}=-(N_{22}{\mu }^2+N^E){\mu }^2,e_{32}=-(N_{23}{\mu }^2+N_E){\mu }^2,\hfill \\ \hfill & e_{33}=-B_3{\mu }^2,e_{34}=-(2A_{12}+Z_{12}){\mu }^2,e_{35}=-(2A_{12}+Z_{12})\mu ,e_{36}=4A_{22}{\mu }^2+2Z_{22}{\mu }^2+N^E,\hfill \\ \hfill & e_{37}=4A_{23}{\mu }^2+2Z_{23}{\mu }^2+N^E,e_{38}=L_{11}^T{\mu }^2,e_{39}=B_4,e_{40}=N_{12},e_{41}=-N_{22},e_{42}=-N_{23},e_{422}=-L_{11}^T,\hfill \\ \hfill & e_{43}=N_{13}{\mu }^3,e_{44}=-(N_{23}{\mu }^2+N^E){\mu }^2,e_{45}=-(N_{33}{\mu }^2+K_{13}+N^E){\mu }^2,e_{46}=-(-Q_{13}^E+B_5){\mu }^2,\hfill \\ \hfill & e_{47}=-(Z_{13}+2A_{13}){\mu }^2,e_{48}=-(Z_{13}+2A_{13})\mu ,e_{49}=4A_{23}{\mu }^2+2Z_{23}{\mu }^2+N^E,\hfill \\ \hfill & e_{50}=4A_{33}{\mu }^2+2Z_{33}{\mu }^2+K_{13}+N^E,e_{51}=-N_{23},e_{52}=-N_{33},e_{53}=N_{13},e_{54}=B_6-Q_{23}^E,\hfill \\ \hfill & e_{55}=O_{11}^T{\mu }^2,e_{56}=-B_1\mu ,e_{57}=B_3{\mu }^2,e_{58}=(-Q_{13}^E+B_5){\mu }^2,e_{59}=-S_{12}^E{\mu }^2-B_7,e_{60}=B_2,\hfill \\ \hfill & e_{61}=-B_4,e_{62}=-B_6+Q_{23}^E,e_{63}=S_{22}^E,e_{64}=-S_3^T,e_{65}=-O_{11}^T,e_{66}=-C_{11}^T,e_{67}=-A^E{S}_3^E,\hfill \\ \hfill & {\displaystyle A^E=\frac{{\int }_0^{A}sin\left(\mu x_1\right)\mathrm{d}{x}_1}{{\int }_0^A{sin}^2\left(\mu x_1\right)\mathrm{d}{x}_1}.}\hfill \end{array}$$

(34)

5.2. Differential Quadrature Solution

In this section, Equation (33) is solved using the differential quadrature approach. The DQM is a common numerical method that has been implemented to solve initial and boundary value problems by many authors (see, e.g., [33,41]). In particular, the DQM predicts precise results using significantly fewer grid points when compared to traditional-order finite difference and finite element analysis (FEM) [42,43]. In accordance with the DQM, the sth-order partial derivatives of a function $F\left(x_2\right)$ at a given discrete point i are approximated by putting a weighted linear sum of the function’s values over the entire domain; the nanocomposite plate is discretized by M mesh points. In particular, the sth-order derivative of a function $F\left(x_2\right)$ at a given discrete point i can be given for M discrete grid points as follows [42]:

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^{s}F\left(x_2\right)}{\mathrm{d}{x}_2^s}{|}_{x_2=x_{2i}}=\sum _{j=1}^M{P}_{ij}^{\left(s\right)}{F}_j,F_j=F\left(x_{2j}\right),i=1,2,\cdots M,\end{array}$$

(35)

in which $P_{ij}^{\left(s\right)}$ are the weighting coefficients associated with the sth derivative, which can be written for higher-order derivatives $(s>1)$ as [42]:

$$\begin{array}{c}\hfill P_{ij}^{\left(s\right)}=\sum _{k=1}^M{P}_{ik}^{\left(1\right)}{P}_{kj}^{(s-1)},i,j=1,2,\cdots M.\end{array}$$

(36)

In contrast, the weighting coefficients associated with the first-order derivative $(s=1)$ is written as [42]:

$$\begin{array}{cc}\hfill & P_{ij}^{\left(1\right)}=\frac{R\left(x_{2i}\right)}{(x_{2i}-x_{2j})R\left(x_{2j}\right)},i,j=1,2,\cdots M;i\ne j,\hfill \\ \hfill & P_{ii}^{\left(1\right)}=-\sum _{k=1}^M{P}_{ki}^{\left(1\right)},i=1,2,\cdots M;k\ne i,\hfill \\ \hfill & R\left(x_{2i}\right)=\prod _{k=1}^M(x_{2i}-x_{2k}),i\ne k.\hfill \end{array}$$

(37)

Furthermore, in accordance with the Gauss–Chebyshev–Lobatto method, the distributed grid points $x_{2i}$ are generated as follows [42]:

$$\begin{array}{c}\hfill x_{2i}=\frac{B}{2}\left[1-cos\left(\pi \frac{i-1}{M-1}\right)\right].\end{array}$$

(38)

By applying the DQ approximation formulas (35)–(38) into Equation (33), the governing equations can be discretized and converted to $(M-2)$ linear algebraic equations as follows:

$$\begin{array}{cc}\hfill & e_{14}{u}_{ni}+e_{12}{w}_{1ni}+e_{13}{w}_{2ni}+e_{17}{\psi }_{ni}+e_{18}{\widehat{T}}_{ni}+e_{11}\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{v}_{nj}+e_{21}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{u}_{nj}+(e_{15}+e_{19})\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{w}_{1nj}\hfill \\ \hfill & +(e_{16}+e_{20})\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{w}_{2nj}=0,\hfill \\ \hfill & e_{22}{v}_{ni}+e_{66}{\widehat{T}}_{ni}^{\prime }+e_{23}\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{u}_{nj}+e_{24}\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{w}_{1nj}+e_{25}\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{w}_{2nj}+e_{26}\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{\psi }_{nj}+e_{27}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{v}_{nj}\hfill \\ \hfill & +e_{28}\sum _{j=1}^M{P}_{ij}^{\left(3\right)}{w}_{1nj}+e_{29}\sum _{j=1}^M{P}_{ij}^{\left(3\right)}{w}_{2nj}=0,\hfill \\ \hfill & e_{30}{u}_{ni}+e_{31}{w}_{1ni}+e_{32}{w}_{2ni}+e_{33}{\psi }_{ni}+e_{38}{\widehat{T}}_{ni}+e_{34}\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{v}_{nj}+e_{35}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{u}_{nj}+e_{41}\sum _{j=1}^M{P}_{ij}^{\left(4\right)}{w}_{1nj}\hfill \\ \hfill & +e_{36}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{w}_{1nj}+e_{37}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{w}_{2nj}+e_{39}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{\psi }_{nj}+e_{40}\sum _{j=1}^M{P}_{ij}^{\left(3\right)}{v}_{nj}+e_{42}\sum _{j=1}^M{P}_{ij}^{\left(4\right)}{w}_{2nj}+e_{422}{\widehat{T}}_{ni}^{″}+{\widehat{q}}_{ni}=0,\hfill \\ \hfill & e_{43}{u}_{ni}+e_{44}{w}_{1ni}+e_{45}{w}_{2ni}+e_{46}{\psi }_{ni}+e_{55}{\widehat{T}}_{ni}+e_{47}\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{v}_{nj}+e_{48}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{u}_{nj}+e_{49}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{w}_{1nj}\hfill \\ \hfill & +e_{50}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{w}_{2nj}+e_{54}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{\psi }_{nj}+e_{65}{\widehat{T}}_{ni}^{″}+e_{53}\sum _{j=1}^M{P}_{ij}^{\left(3\right)}{v}_{nj}+e_{51}\sum _{j=1}^M{P}_{ij}^{\left(4\right)}{w}_{1nj}+e_{52}\sum _{j=1}^M{P}_{ij}^{\left(4\right)}{w}_{2nj}+{\widehat{q}}_{ni}=0,\hfill \\ \hfill & e_{56}{u}_{ni}+e_{57}{w}_{1ni}+e_{58}{w}_{2ni}+e_{59}{\psi }_{ni}+e_{60}\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{v}_{nj}+e_{61}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{w}_{1nj}+e_{62}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{w}_{2nj}\hfill \\ \hfill & +e_{63}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{\psi }_{nj}+e_{64}{\widehat{T}}_{ni}+e_{67}=0,i=2,3,\cdots ,(M-1).\hfill \end{array}$$

(39)

Furthermore, the boundary conditions can be expressed as follows in a discretization form:

Simply supported:

$$\begin{array}{cc}\hfill & u_{ni}=w_{1ni}=w_{2ni}=0,\hfill \\ \hfill & Z_{12}{\mu }^2{w}_{1ni}+Z_{13}{\mu }^2{w}_{2ni}-Z_{11}\mu u_{ni}+B_2{\psi }_{ni}-C_{11}^T{\widehat{T}}_{ni}+N_{11}\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{v}_{nj}\hfill \\ \hfill & -N_{12}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{w}_{1nj}-N_{13}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{w}_{2nj}+C_{22}^E{A}^E=0,\hfill \\ \hfill & Z_{22}{\mu }^2{w}_{1ni}+Z_{23}{\mu }^2{w}_{2ni}-Z_{12}\mu u_{ni}+B_4{\psi }_{ni}-L_{11}^T{\widehat{T}}_{ni}+N_{12}\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{v}_{nj}\hfill \\ \hfill & -N_{22}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{w}_{1nj}-N_{23}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{w}_{2nj}+L_{22}^E{A}^E=0,i=1,M.\hfill \end{array}$$

(40)

Clamped:

$$\begin{array}{c}\hfill u_{ni}=v_{ni}=w_{1in}=w_{1ni}={\psi }_{ni}=0,i=1,M.\end{array}$$

(41)

Free:

$$\begin{array}{cc}\hfill & Z_{12}{\mu }^2{w}_{1ni}+Z_{13}{\mu }^2{w}_{2ni}-Z_{11}\mu u_{ni}+B_2{\psi }_{ni}-C_{11}^T{\widehat{T}}_{ni}+N_{11}\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{v}_{nj}\hfill \\ \hfill & -N_{12}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{w}_{1nj}-N_{13}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{w}_{2nj}+C_{22}^E{A}^E=0,\hfill \\ \hfill & Z_{22}{\mu }^2{w}_{1ni}+Z_{23}{\mu }^2{w}_{2ni}-Z_{12}\mu u_{ni}+B_4{\psi }_{ni}-L_{11}^T{\widehat{T}}_{ni}+N_{12}\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{v}_{nj}\hfill \\ \hfill & -N_{22}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{w}_{1nj}-N_{23}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{w}_{2nj}+L_{22}^E{A}^E=0,\hfill \\ \hfill & Z_{23}{\mu }^2{w}_{1ni}+Z_{33}{\mu }^2{w}_{2ni}-Z_{13}\mu u_{ni}+B_6{\psi }_{ni}-O_{11}^T{\widehat{T}}_{ni}+N_{13}\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{v}_{nj}\hfill \\ \hfill & -N_{23}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{w}_{1nj}-N_{33}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{w}_{2nj}+O_{22}^E{A}^E=0,\hfill \\ \hfill & 2A_{12}\left(-{\mu }^2{v}_{ni}-\mu \sum _{j=1}^M{P}_{ij}^{\left(1\right)}{u}_{nj}\right)+4A_{22}{\mu }^2\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{w}_{1nj}+4A_{23}{\mu }^2\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{w}_{2nj}\hfill \\ \hfill & -Z_{12}\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{u}_{nj}+Z_{22}{\mu }^2\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{w}_{1nj}+Z_{23}{\mu }^2\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{w}_{2nj}+N_{12}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{v}_{nj}\hfill \\ \hfill & -N_{22}\sum _{j=1}^M{P}_{ij}^{\left(3\right)}{w}_{1nj}-N_{23}\sum _{j=1}^M{P}_{ij}^{\left(3\right)}{w}_{2nj}-L_{11}^T{\widehat{T}}_{ni}^{\prime }+B_4\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{\psi }_{nj}=0,\hfill \\ \hfill & 2A_{13}\left(-{\mu }^2{v}_{ni}-\mu \sum _{j=1}^M{P}_{ij}^{\left(1\right)}{u}_{nj}\right)+4A_{23}{\mu }^2\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{w}_{1nj}+4A_{33}{\mu }^2\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{w}_{2nj}\hfill \\ \hfill & -Z_{13}\mu \sum _{j=1}^M{P}_{ij}^{\left(1\right)}{u}_{nj}+Z_{23}{\mu }^2\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{w}_{1nj}+Z_{33}{\mu }^2\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{w}_{2nj}+N_{13}\sum _{j=1}^M{P}_{ij}^{\left(2\right)}{v}_{nj}\hfill \\ \hfill & -N_{23}\sum _{j=1}^M{P}_{ij}^{\left(3\right)}{w}_{1nj}-N_{33}\sum _{j=1}^M{P}_{ij}^{\left(3\right)}{w}_{2nj}-O_{11}^T{\widehat{T}}_{ni}^{\prime }+B_6\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{\psi }_{nj}\hfill \\ \hfill & +K_{13}\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{w}_{2nj}-Q_{23}^E\sum _{j=1}^M{P}_{ij}^{\left(1\right)}{\psi }_{nj}=0,i=1,M.\hfill \end{array}$$

(42)

6. Numerical Results

Various numerical examples are introduced in the current section to illustrate the influences of different parameters such as porosity factor, porosity distribution types, GPLs weight fraction, GPLs distribution patterns, side-to-thickness ratio, temperature parameter, temperature exponent, external electric voltage and boundary conditions on the deflection and stresses of the FG porous piezoelectric plate reinforced with GPLs. For this purpose, the following quantities are defined:

$$\begin{array}{cc}\hfill & W^{\ast }=\frac{{10}^{-3}W{E}_m}{q_0},{\sigma }_{11}^{\ast }=\frac{{10}^{2}H}{A^2{E}_m{q}_0}{\sigma }_{11}\left(\frac{A}{2},\frac{B}{2},Z\right),\hfill \\ & {\sigma }_{13}^{\ast }=-\frac{{10}^{4}H}{A^2{E}_m{q}_0}{\sigma }_{13}\left(0,\frac{B}{2},Z\right),{\sigma }_{12}^{\ast }=-\frac{{10}^{2}H}{A^2{E}_m{q}_0}{\sigma }_{12}\left(0,0,Z\right),\hfill \\ \hfill & V_0=\frac{{\overline{V}}_0}{H^2{E}_G},Z=\frac{x_3}{H}.\hfill \end{array}$$

The material properties of the piezoelectric matrix and GPLs are given in

Table 1. Material properties of the piezoelectric matrix and GPLs [44].

Materials	Piezoelectric	GPLs
E(GPa)	1.4	1010
$\nu$	0.29	0.186
$\rho$(g/cm${}^3$)	1.92	1.06
$\alpha$(${10}^{-6}$K${}^{-1}$)	60	5
${\zeta }_{31}={\zeta }_{32}$ (${10}^{-3}{\mathrm{C}/\mathrm{m}}^2$)	50.535	50.535 $S_0$
${\zeta }_{24}={\zeta }_{15}$ (${10}^{-3}{\mathrm{C}/\mathrm{m}}^2$)	−15.93	−15.93 $S_0$
$f_{11}=f_{22}$ (${10}^{-9}\mathrm{C}/\mathrm{Vm}$)	0.5385	0.5385 $S_0$
$f_{33}$ (${10}^{-9}\mathrm{C}/\mathrm{Vm}$)	0.59571	0.59571 $S_0$

, noting that the parameter $S_0$ is named as the piezoelectric multiple [44].

The following fixed data (unless otherwise declared) are used in the numerical examples: $A/H=10,B/A=1,S_0=100,e_0=0.2,n=1,W_f=0.1,k=1,V_0=0.1,$$q_0=10$ Pa, $H=0.002$ m, $L^G=2.5\mathsf{\mu }$m, $b^G=1.5\mathsf{\mu }$m, $h^G=1.5$nm.

Firstly,

Table 2. Convergence of the results of FG porous nanocomposite piezoelectric square plates.

M	${\mathit{W}}^{\ast }$	${\mathit{\sigma }}_{11}^{\ast }$	${\mathit{\sigma }}_{13}^{\ast }$	${\mathit{\sigma }}_{12}^{\ast }$
7	7.047371446	3.769956328	12.185145430	5.930577935
9	7.047291563	3.770929050	12.159459440	5.931145760
11	7.047291007	3.771028796	12.160105040	5.931197603
13	7.047291003	3.771028832	12.160092250	5.931196924
15	7.047291003	3.771028839	12.160092460	5.931196935
17	7.047291003	3.771028839	12.160092460	5.931196931

displays a convergence study for the present results of the FG porous piezoelectric plate reinforced with GPLs. In this table, we determine the minimum number of mesh points required for converged solutions in the DQM. It is notable that 15 mesh points are sufficient to achieve a converged solution.

Secondly, in order to check the accuracy of the results obtained by the current formulations, the present deflection ($\overline{W}=10W{E}_c{H}^3/\left(q_0{A}^4\right)$) of an FG square plate is compared with those provided by Thai and Kim [45], as shown in

Table 3. Comparison of the central deflection $\overline{W}$ of an FGM square plate.

$\mathit{A}/\mathit{H}$	Source	$\mathit{\kappa }=0$	$\mathit{\kappa }=1$	$\mathit{\kappa }=2$	$\mathit{\kappa }=5$	$\mathit{\kappa }=10$
5	Ref. [45]	0.3433	0.6688	0.8671	1.0885	1.2276
	Present	0.3511	0.6844	0.8872	1.1122	1.2536
10	Ref. [45]	0.2961	0.5890	0.7573	0.9114	1.0087
	Present	0.3038	0.6046	0.7773	0.9350	1.0347
20	Ref. [45]	0.2842	0.5689	0.7298	0.8669	0.9538
	Present	0.2920	0.5846	0.7498	0.8906	0.9798
100	Ref. [45]	0.2804	0.5625	0.7209	0.8527	0.9362
	Present	0.2882	0.5782	0.7410	0.8764	0.9622

. The FG plate is composed of ceramic (c) and metal (m) with the following properties: $E_c=380$ GPa, $E_m=70$ GPa, ${\nu }_c={\nu }_m=0.3$. The effective Young’s modulus $E_{eff}$ is calculated as:

$$E_{eff}=E_m+(E_c-E_m){\left(\frac{x_3}{H}+\frac{1}{2}\right)}^{\kappa }.$$

Moreover, the results in this table are given for different values of the side-to-thickness ratio $A/H$ and power law index $\kappa$. From Table 3, the present deflection $\overline{W}$ is in good agreement with that obtained by Thai and Kim [45].

Table 4. Dimensionless central deflection $W^{\ast }$ of FG porous nanocomposite piezoelectric square plates under various boundary conditions and for different values of side-to-thickness ratio $A/H$ and GPLs weight fraction $W_f$.

$\mathit{A}/\mathit{H}$	${\mathit{W}}_{\mathit{f}}$	SS	SC	CC	FS	FC
5	0.1	1.9502024	1.9502248	1.9502472	2.2611658	2.2612573
	0.2	1.6322019	1.6322216	1.6322414	1.8975372	1.8976191
	0.3	1.3590695	1.3590865	1.3591035	1.5833144	1.5833859
	0.4	1.1226046	1.1226191	1.1226335	1.3101761	1.3102377
10	0.1	7.0473344	7.0473355	7.0473365	7.7182566	7.7182636
	0.2	5.9247160	5.9247169	5.9247178	6.5231102	6.5231166
	0.3	4.9475774	4.9475782	4.9475790	5.4679584	5.4679640
	0.4	4.0963349	4.0963356	4.0963363	4.5414866	4.5414915
15	0.1	13.6580812	13.6580816	13.6580819	13.8244905	13.8244918
	0.2	11.5497255	11.5497258	11.5497261	11.7910596	11.7910611
	0.3	9.6816468	9.6816470	9.6816473	9.9431761	9.9431777
	0.4	8.0407740	8.0407742	8.0407744	8.2987022	8.2987037
20	0.1	20.3341739	20.3341741	20.3341743	18.8196737	18.8196646
	0.2	17.2976822	17.2976824	17.2976825	16.2053831	16.2053769
	0.3	14.5566201	14.5566203	14.5566204	13.7525487	13.7525442
	0.4	12.1282823	12.1282824	12.1282826	11.5374732	11.5374700
30	0.1	31.2422929	31.2422930	31.2422931	23.7936452	23.7934750
	0.2	26.8380929	26.8380930	26.8380931	20.8718941	20.8717673
	0.3	22.7327667	22.7327668	22.7327669	17.9352086	17.9351111
	0.4	19.0428742	19.0428742	19.0428743	15.2002292	15.2001539

,

Table 5. Dimensionless normal stress ${\sigma }_{11}^{\ast }$ of FG porous nanocomposite piezoelectric square plates under various boundary conditions and for different values of side-to-thickness ratio $A/H$ and GPLs weight fraction $W_f$ ($x_3/H=0.5$).

$\mathit{A}/\mathit{H}$	${\mathit{W}}_{\mathit{f}}$	SS	SC	CC	FS	FC
5	0.1	17.9083147	14.2299878	8.6026809	22.0421057	17.6169790
	0.2	39.7862064	32.7299747	21.9766406	46.2524481	37.9476666
	0.3	66.0203904	55.6527912	39.9056717	73.4591357	61.5069862
	0.4	95.0825801	81.4973138	60.9227184	102.4762900	87.0968504
10	0.1	3.7706838	2.8510300	1.4441312	4.4323964	3.3258352
	0.2	8.8800785	7.1159003	4.4274463	9.9192326	7.8425622
	0.3	15.2775451	12.6854970	8.7485686	16.4571412	13.4685088
	0.4	22.5133983	19.1169208	13.9731110	23.6512968	19.8057824
15	0.1	1.2687137	0.8599784	0.2346896	1.3957948	0.9039876
	0.2	3.3256340	2.5415543	1.3466853	3.5225628	2.5995945
	0.3	6.0703433	4.9183212	3.1685746	6.2782875	4.9500017
	0.4	9.2648491	7.7553028	5.4691644	9.4364205	7.7272979
20	0.1	0.4823713	0.2524577	−0.0992672	0.4856278	0.2089857
	0.2	1.5136873	1.0726425	0.4005286	1.5134553	0.9942848
	0.3	2.9977062	2.3496937	1.3654612	2.9799047	2.2327430
	0.4	4.7794906	3.9303708	2.6444179	4.7316853	3.7703029
30	0.1	0.0464381	−0.0557457	−0.2120679	0.0200638	−0.1028896
	0.2	0.4094121	0.2133922	−0.0853251	0.3621829	0.1314385
	0.3	1.0215842	0.7335787	0.2961420	0.9549210	0.6228470
	0.4	1.7994355	1.4220489	0.8505143	1.7139098	1.2866266

,

Table 6. Dimensionless transverse shear stress ${\sigma }_{13}^{\ast }$ of FG porous nanocomposite piezoelectric square plates under various boundary conditions and for different values of side-to-thickness ratio $A/H$ and GPLs weight fraction $W_f$ ($x_3/H=0$).

$\mathit{A}/\mathit{H}$	${\mathit{W}}_{\mathit{f}}$	SS	SC	CC	FS	FC
5	0.1	27.5725468	27.5728630	27.5731795	35.5991272	35.6004401
	0.2	41.8047883	41.8052929	41.8057982	54.0959779	54.0981093
	0.3	48.3696904	48.3702947	48.3708997	62.7064278	62.7090189
	0.4	49.7835309	49.7841707	49.7848114	64.6427662	64.6455467
10	0.1	12.4811878	12.4811896	12.4811915	15.6183845	15.6183973
	0.2	19.0079651	19.0079681	19.0079710	23.8609151	23.8609361
	0.3	22.0559114	22.0559150	22.0559185	27.7532155	27.7532413
	0.4	22.7531791	22.7531828	22.7531866	28.6890691	28.6890970
15	0.1	7.1699769	7.1699770	7.1699772	8.6438176	8.6438183
	0.2	10.9833193	10.9833196	10.9833199	13.2912633	13.2912648
	0.3	12.7930095	12.7930099	12.7930102	15.5239358	15.5239380
	0.4	13.2383284	13.2383288	13.2383292	16.1019428	16.1019454
20	0.1	4.5039926	4.5039926	4.5039926	5.2462239	5.2462218
	0.2	6.9405200	6.9405200	6.9405201	8.1157913	8.1157887
	0.3	8.1156866	8.1156866	8.1156867	9.5164542	9.5164516
	0.4	8.4250986	8.4250987	8.4250988	9.9026979	9.9026956
30	0.1	2.0506008	2.0506008	2.0506008	2.2739293	2.2739175
	0.2	3.1909701	3.1909701	3.1909701	3.5497144	3.5496984
	0.3	3.7556360	3.7556360	3.7556360	4.1872969	4.1872798
	0.4	3.9198921	3.9198922	3.9198922	4.3788605	4.3788441

and

Table 7. Dimensionless in-plane shear stress ${\sigma }_{12}^{\ast }$ of FG porous nanocomposite piezoelectric square plates under various boundary conditions and for different values of side-to-thickness ratio $A/H$ and GPLs weight fraction $W_f$ ($x_3/H=0.5$).

$\mathit{A}/\mathit{H}$	${\mathit{W}}_{\mathit{f}}$	SS	SC	CC	FS	FC
5	0.1	25.3738484	19.5453075	31.9011320	11.1960158	11.1960758
	0.2	38.5561770	27.5053073	51.5049609	18.1304595	18.1305722
	0.3	42.7533121	26.6810333	62.3163630	22.0082296	22.0083827
	0.4	40.3460882	19.4731472	66.5907895	23.5899193	23.5900995
10	0.1	5.9309953	4.4738289	7.5626158	3.0122594	3.0122599
	0.2	8.9985189	6.2357483	12.2353733	4.8493441	4.8493449
	0.3	9.9328869	5.9147504	14.8232207	5.8603925	5.8603936
	0.4	9.2961913	4.0778823	15.8568939	6.2585055	6.2585067
15	0.1	2.3982223	1.7505927	3.1233861	1.4290626	1.4290626
	0.2	3.6262604	2.3983621	5.0648607	2.2915821	2.2915822
	0.3	3.9716420	2.1858032	6.1451217	2.7601589	2.7601590
	0.4	3.6656777	1.3464290	6.5815433	2.9392490	2.9392491
20	0.1	1.2139283	0.8496366	1.6218329	0.8331441	0.8331440
	0.2	1.8253245	1.1346317	2.6345371	1.3343452	1.3343451
	0.3	1.9774569	0.9729225	3.2000390	1.6048102	1.6048101
	0.4	1.7903302	0.4857528	3.4305045	1.7065881	1.7065880
30	0.1	0.4414378	0.2795304	0.6227287	0.3670740	0.3670736
	0.2	0.6530614	0.3460869	1.0127115	0.5888630	0.5888625
	0.3	0.6876045	0.2411448	1.2309744	0.7083674	0.7083669
	0.4	0.5914927	0.0116805	1.3204590	0.7532405	0.7532400

depict the effects of the side-to-thickness ratio $A/H$ and GPLs weight fraction $W_f$ on the dimensionless central deflection $W^{\ast }$, the normal stress ${\sigma }_{11}^{\ast }$, the transverse shear stress ${\sigma }_{13}^{\ast }$ and the in-plane shear stress ${\sigma }_{12}^{\ast }$ of FG porous nanocomposite piezoelectric square plates under various boundary conditions. It is observed that regardless of the type of boundary conditions, central deflection $W^{\ast }$ gradually increases with increasing side-to-thickness ratio $A/H$. On the other hand, it decreases as the GPLs weight fraction $W_f$ increases. Furthermore, the sensitivity of the deflection to variations of the boundary conditions is very noticeable. It is notable that the plates with clamped edges have minimum values of normal and transverses stresses, while plates with free and clamped edges suffer more deflection and stresses. It can be also observed that the stresses ${\sigma }_{11}^{\ast }$, ${\sigma }_{13}^{\ast }$ and ${\sigma }_{12}^{\ast }$ decrease as the side-to-thickness ratio $A/H$ increases. They increase as the GPLs weight fraction $W_f$ increases, except the in-plane shear stress ${\sigma }_{12}^{\ast }$, which no longer increases as $W_f$ increases.

The effect of the porosity factor $e_0$ on the central deflection $W^{\ast }$ of FG porous nanocomposite piezoelectric square plates is plotted in Figure 4. Four different porosity distribution types—I, II, III and IV—are illustrated in the corresponding figures (a), (b), (c) and (d), respectively. The effect of the porosity factor $e_0$ on the central deflection $W^{\ast }$ is more pronounced for small values of $q_0$, especially for Type I and IV. In contrast, the effect of $e_0$ on $W^{\ast }$ is minimum for large values of $q_0$, especially for porosity Type II, since the thermal conductivity of the porosity is much lower than that of the composite materials. Therefore, increasing the porosity reduces the thermal effects on the plate. Subsequently, a noticeable reduction in the deflection occurs as the porosity factor increases, especially for Types I and III. However, for Type IV, the deflection behaves in an opposite manner to the variation of the porosity factor. Note that, at the top surface of the plate (Type IV), there is no porosity and the temperature is maximum (see, Equation (14)); then, the porosity linearly increases through the thickness; while the temperature linearly decreases ($k=1$) in the same direction. Therefore, the role of pores in reducing the effects of the temperature on the plate stiffness are very weak or may be missed because the temperature naturally decreases in the direction of increasing pores. Subsequently, the deflection increases as the porosity factor increases. This also explains the weak role of pores to reduce the deflection for Type II.

For different porosity distribution Types I, II, III and IV, the impact of the porosity factor $e_0$ on the normal stress ${\sigma }_{11}^{\ast }$, transverse shear stress ${\sigma }_{13}^{\ast }$ and in-plane shear stress ${\sigma }_{12}^{\ast }$ through the thickness of FG porous nanocomposite piezoelectric square plates is illustrated in Figure 5, Figure 6 and Figure 7, respectively. It is obvious that, generally, for all porosity distribution types, the impact of the porosity factor $e_0$ on the normal stress ${\sigma }_{11}^{\ast }$, transverse shear stress ${\sigma }_{13}^{\ast }$ and in-plane shear stress ${\sigma }_{12}^{\ast }$ is significantly dependent on the porosity type. Moreover, for Type I, the porosity is uniformly distributed through the thickness; therefore, the plate becomes an isotropic structure. Accordingly, the normal stress and in-plane shear stress are linearly varied through the thickness of the plate, whereas the transverse shear stress is parabolically changed. It is found that the maximum normal stress ${\sigma }_{11}^{\ast }$ decreases as the porosity factor $e_0$ increases. Note that for porosity distribution Types I and II, the maximum stress ${\sigma }_{11}^{\ast }$ occurs at the bottom surface of the plate, while it occurs near the bottom surface for Types III and IV, as shown in Figure 5.

In Figure 6a,c, the maximum transverse shear stress ${\sigma }_{13}^{\ast }$ occurs at the mid-plane of the plate for all values of porosity factor $e_0$. However, for Type IV, it occurs near the mid-plane of the plate, as shown in Figure 6d, because this type has an asymmetric porosity distribution. Since the volume fraction of porosity in Type II is maximum at the mid-plane ($x_3=0$), the shear stress ${\sigma }_{13}^{\ast }$ at $x_3=0$ is very sensitive to the variation of the porosity factor $e_0$. Therefore, the stress ${\sigma }_{13}^{\ast }$ has non-parabolic shapes for the largest values of $e_0$. Note that for porosity distribution Types I, II and III, the maximum stress ${\sigma }_{13}^{\ast }$ decreases with increasing $e_0$, while for Type IV, it increases with the increase of $e_0$.

It can be observed that for porosity distribution Types I, II and IV, the maximum in-plane shear stress ${\sigma }_{12}^{\ast }$ occurs at the top surface of the plate, as shown in Figure 7, while for Type III, the maximum ${\sigma }_{12}^{\ast }$ occurs near the top surface of the plate. Furthermore, for Types I, II and III, the maximum stress ${\sigma }_{12}^{\ast }$ is decreased as the porosity factor $e_0$ increases. However, this trend is reversed for Type IV; increasing $e_0$ leads to the increase of the maximum stress ${\sigma }_{12}^{\ast }$.

To explain the impact of the GPLs weight fraction $W_f$ on the obtained results, the central deflection $W^{\ast }$ of FG porous nanocomposite piezoelectric square plates against transverse load $q_0$ is plotted in Figure 8 for various values of $W_f$ and for different graphene distribution patterns. In agreement with the review of the literature, the proportion of graphene in the plates greatly improves the mechanical properties of the plates and enhances their stiffness. Therefore, we notice a successive decreasing of the central deflection $W^{\ast }$ as the weight fraction $W_f$ increases.

The Influence of the GPLs weight fraction $W_f$ on the normal stress ${\sigma }_{11}^{\ast }$, transverse shear stress ${\sigma }_{13}^{\ast }$ and in-plane shear stress ${\sigma }_{12}^{\ast }$ through the thickness of FG porous nanocomposite piezoelectric square plates is plotted in Figure 9, Figure 10 and Figure 11, respectively, for different graphene distribution patterns. We can clearly observe from Figure 9 that the maximum stress ${\sigma }_{11}^{\ast }$ increases as the GPLs weight fraction $W_f$ increases. For more clarity, ${\sigma }_{11}^{\ast }$ is linearly varied through the thickness of the plate for GPLs Pattern A (see Figure 9a), indicating the even distribution of GPLs through the plate thickness. In particular, for GPLs pattern B (see Figure 9b), the normal stress ${\sigma }_{11}^{\ast }$ is nearly zero at the top and bottom surfaces of the plate, and the maximum ${\sigma }_{11}^{\ast }$ occurs nearly at the mid-plane of the plate, because the GPLs weight fraction is maximum at the middle plane of the plate and equals to zero at the upper and lower faces. It is also notable that the normal stress ${\sigma }_{11}^{\ast }$ in the upper part of the plate of Pattern B has extremums at the large values of $W_f$ (as shown in Figure 9b) due to the sensitivity of graphene to temperature, which has a maximum value at the top surface (see Equation (14)). Furthermore, since the temperature has small values in the lower part of the plate, the stress ${\sigma }_{11}^{\ast }$ has no extremums in this part. Moreover, due to the asymmetric distribution of GPLs through the thickness of the plate of Pattern D (see, Figure 9d), the normal stress ${\sigma }_{11}^{\ast }$ is zero at the bottom surface of the plate and is maximum at the top surface of the plate.

It can be seen from Figure 10 that the stress ${\sigma }_{13}^{\ast }$ of all GPL distribution patterns has the same behavior with the variation of $W_f$. It increases to reach its maximum and then decreases as $W_f$ increases. These maximums depend on the GPL distribution pattern. This means that, adding more graphene to the structures makes the stress ${\sigma }_{13}^{\ast }$ behave in an opposite sense, especially with considering the thermal load, since graphene possesses high thermal conductivity.

In addition, owing to the same reasons discussed in Figure 9b, the maximum stress ${\sigma }_{12}^{\ast }$ in the upper part of the plate has the same distribution with the variation of $W_f$ for all patterns (see Figure 11); it no longer increases as the weight fraction $W_f$ increases, while, the maximum ${\sigma }_{12}^{\ast }$ in the lower part of the plate increases monotonically as $W_f$ increases. Note that the positive and negative signs indicate the tensile and compressive stresses, respectively.

It is noteworthy that the stresses ${\sigma }_{11}^{\ast }$, ${\sigma }_{13}^{\ast }$ and ${\sigma }_{12}^{\ast }$ significantly depend on the volume fraction of GPLs $v_G$, they equal zero when $v_G$ is equal to zero. Note that the volume fractions of GPLs patterns B, C and D are equal to zero at $x_3=\pm H/2$, $x_3=0$ and $x_3=-H/2$, respectively.

Figure 12 depicts the variation of the temperature parameter $T_0$ on the central deflection $W^{\ast }$, normal stress ${\sigma }_{11}^{\ast }$, transverse shear stress ${\sigma }_{13}^{\ast }$ and in-plane shear stress ${\sigma }_{12}^{\ast }$ of FG porous nanocomposite piezoelectric square plates. As is well known, an increase in temperature weakens the structures, so the deflection $W^{\ast }$ and the shear stresses ${\sigma }_{13}^{\ast }$ and ${\sigma }_{12}^{\ast }$ increase monotonically as the temperature parameter $T_0$ increases. Moreover, the normal and shear stresses may be equal to zero at the top and bottom surfaces of the plate, indicating the significant dependence of the stresses on the GPLs volume fraction. Since the normal stress ${\sigma }_{11}^{\ast }$ is directly dependent on the temperature (see Equation (12)), it is more sensitive to variation of the temperature and behaves differently from the shear stresses. Note that the temperature is maximum at the top surface of the plate and then decreases gradually to its minimum at the lower surface. Therefore, for small values of temperature ($T_0$ = 50,100 or in the lower part of the plate), the maximum tensile stress ${\sigma }_{11}^{\ast }$ decreases as the temperature parameter increases, while for large values of temperature (that occur in the upper part of the plate), the maximum compressive stress ${\sigma }_{11}^{\ast }$ increases with the increase in temperature.

In Figure 13, the impact of the temperature exponent k on the central deflection $W^{\ast }$, normal stress ${\sigma }_{11}^{\ast }$, transverse shear stress ${\sigma }_{13}^{\ast }$ and in-plane shear stress ${\sigma }_{12}^{\ast }$ of FG porous nanocomposite piezoelectric square plates is investigated. It can be observed that the temperature exponent k has a hardening effect. Accordingly, the deflection $W^{\ast }$ is decreased as k and $q_0$ increase. Furthermore, it is seen that for linear variation of the temperature ($k=1$), the normal stress ${\sigma }_{11}^{\ast }$ is also linearly changed through the thickness of the plate. The in-plane shear stress ${\sigma }_{12}^{\ast }$ is linearly varied with respect to Z for all values of k. The shear stresses ${\sigma }_{13}^{\ast }$ and ${\sigma }_{12}^{\ast }$ decrease with increasing k, while this trend may be reversed for the normal stress, especially for $-0.35\le Z\le 0.35$.

The influence of the applied electric voltage $V_0$ on the central deflection $W^{\ast }$, the normal stress ${\sigma }_{11}^{\ast }$, the transverse shear stress ${\sigma }_{13}^{\ast }$ and the in-plane shear stress ${\sigma }_{12}^{\ast }$ of FG porous nanocomposite piezoelectric square plates is illustrated in Figure 14. The deflection $W^{\ast }$ decreases gradually as the electric voltage $V_0$ increases. This is reversed for the stresses; they increase as $V_0$ increases.

7. Conclusions

Levy-type solution and the DQM are utilized here to analyze the thermo-piezoelectric bending of FG porous piezoelectric plates reinforced with GPLs. The displacement field is formulated based on a refined four-variable shear deformation plate theory. To reduce the temperature effects, various types of porosity are considered in the present composite plate. This plate is composed of a piezoelectric material containing internal pores and is supported by FG GPLs. The differential equations are deduced by using the principle of virtual work including thermal, mechanical and electrical loads. These equations are converted to ordinary systems by using the Levy procedure and then discretized by the DQ method. The present formulations and theory are examined by comparing them to some examples from other research in the literature. Additional numerical results are presented to explain the effects of various parameters on the bending of the reinforced composite plates. Moreover, it is concluded that a considerable effect on the deflection and stresses is seen with variation of the edge boundary conditions. The deflection and stresses monotonically increase as the temperature parameter increases because elevated temperature leads to a reduction in the plate stiffness. The deflection decreases as the porosity factor, temperature exponent, transverse load, electric voltage and GPLs weight fraction increase. In addition, the maximum normal stresses decrease as the porosity factor increases. Increasing the electric voltage and the GPLs weight fraction leads to an increase in the shear stresses. The variation of the stresses through the thickness of the plate depends considerably on the porosity and GPLs distribution type. Furthermore, it can be concluded that with presence of a thermal load, adding a certain amount of graphene to structures enhances their toughness; consequently, a decrease in the deflection occurs. By adding more graphene, the structures become more flexible; accordingly, the deflection increases. To reduce the deformation of engineering structures caused by temperature, it is recommended to use structures reinforced with GPLs containing porosity, especially porosity Types I and III. However, it is not recommended to use porosity Type II or IV to reduce the thermal effect. The present results can be used as a benchmark for future investigations.