Analytical Solution for Static and Dynamic Analysis of Graphene-Based Hybrid Flexoelectric Nanostructures

Abstract

Owing to their applications in devices such as in electromechanical sensors, actuators and nanogenerators, the consideration of size-dependent properties in the electromechanical response of composites is of great importance. In this study, a closed-form solution based on the linear piezoelectricity, Kirchhoff’s plate theory and Navier’s solution was developed, to envisage the electromechanical behaviors of hybrid graphene-reinforced piezoelectric composite (GRPC) plates, considering the flexoelectric effect. The governing equations and respective boundary conditions were obtained, using Hamilton’s variational principle for achieving static deflection and resonant frequency. Moreover, the different parameters considering aspect ratio, thickness of plate, different loadings (inline, point, uniformly distributed load (UDL), uniformly varying load (UVL)), the combination of different volume fraction of graphene and piezoelectric lead zirconate titanate are considered to attain the desired bending deflection and frequency response of GRPC. Different mode shapes and flexoelectric coefficients are also considered and the results reveal that the proper addition of graphene percentage and flexoelectric effect on the static and dynamic responses of GRPC plate is substantial. The obtained results expose that the flexoelectric effect on the piezoelastic response of the bending of nanocomposite plates are worth paying attention to, in order to develop a nanoelectromechanical system (NEMS). Our fundamental study sheds the possibility of evolving lightweight and high-performance NEMS applications over the existing piezoelectric materials.

1. Introduction

In recent years, the size-dependent properties, such as piezoelectric, flexoelectric and surface effects, find a unique niche in the various applications of engineering, such as energy harvesters, including nanoactuators, nanosensors, pulse power nanogenerators and transducers, and, most importantly, in the field of structural health monitoring, due to its extraordinary remarkable characteristics. The recent advancement in nanotechnologies has renewed the interest not only in size-dependent piezoelectricity but also in flexoelectricity, since the large strain gradients often present at the nanoscale that may lead to strong flexoelectric effects even in the non-piezoelectric materials. Piezoelectric effect was first noticed in the 1880 by P. Curie and J. Curie. Piezoelectric material exhibits the property to generate the electrical response to an applied homogeneous strain and vice versa. Moreover, the flexoelectric effect was observed by Mashkevich and Tolpygo [1] and defined as the spontaneous electric response to an applied strain gradient (inhomogeneous strain) and is allowed by symmetry in all materials [2,3,4,5]. Similar to piezoelectric effect, the flexoelectric effect is also characterized into two types: (i) direct flexoelectric effect, which is coupling between the electric response (in terms of polarization or electric field) and gradient of mechanical responses (stress gradient i.e., higher order stress or strain gradient); and (ii) converse flexoelectric effect, which refers to coupling between the mechanical responses (stress/strain) and gradient of electric response (electric field gradient/polarization gradient). Contrary to the piezoelectricity, flexoelectricity exits only in all dielectric materials. In addition to this, for the active control of lightweight smart structures subjected to static/dynamic loading conditions, these piezoelectric structures with flexoelectric effect plays significant role [6,7,8,9,10]. For instance, the most common and extensively used monolithic piezoelectric materials are lead zirconate titanate (PZT) and polyvinylidene fluoride (PVDF). According to recent advances in the composite material, the piezoelectric composites are emerged as new efficient materials that may be used as distributed sensors and actuators for the smart structures. In such piezoelectric composite material, the monolithic piezoelectric material (PZT/PVDF) is used as a fiber-reinforcing phase in the polymer matrix, such as epoxy.

If the 20th century was the era of “plastics”, the 21st century appears set to become the era of “graphene”. Such a ground-breaking discovery of graphene material was made mainly by Russian scientists, Novoselov and Geim [11] by using the “Scotch tape method” in 2004. Graphene is a one atom thick sheet of sp² hybridized carbon atoms with remarkably high elastic modulus, tensile strength and high specific surface area (~3 times that of CNTs), that has encouraged huge interests among research and industrial communities. It has been reported in many experimental studies that graphene-based composites show substantial enhancement in mechanical properties as that of CNT-based composites [12]. In this work, they reported that the use of small quantity of graphene in the conventional matrices is found to alter and improve the properties of resulting nanocomposite significantly. Very few research studies [13,14,15] have been dedicated to the introduction of graphene or its derivatives as the modifiers/interphase to the conventional bulk composites in order to improve their multifunctional properties. The thickness of the graphene is in the order of nanometer. Hence, the concepts of validation of continuum mechanics at such nanoscales are still doubtful. However, classical homogenization strategies for graphene-reinforced composites are commonly implemented in the scientific literature with relatively reasonable outcomes. Therefore, many researchers [16,17,18] considered graphene layers as continuum plate in different types of analysis. Hence, in this paper, we are concentrating on graphene as nanofillers in a piezoelectric composite which was not conveyed in scientific literatures earlier.

Here, a review is presented on different structural elements considering the piezoelectric and flexoelectric effect and is characterized with different modeling theories. The closed-form solution was derived by Ray et al. [6,19], for investigating the static behavior of simply supported piezoelectric plate subjected to the cylindrical bending, and they analyzed its electrostatic behavior in terms of deformation, stresses, electric potential and electric displacement. They supposed the laminate structure with a substrate of graphite/epoxy and coupled with PVDF. Ray and Pradhan [20] investigated the 1–3 piezoelectric vertically reinforced laminated composite to improve the performance of hybrid damping, using the active constrained layer-damping treatment. Afterwards, Majdoub et al. [21] conveyed the flexoelectric effect for energy harvester, such as piezoelectric nanobeam/ribbon due to substantial enrichment of the piezoelectric coefficient. Using the nonlocal differential constitutive relations of Eringen, Pradhan [22] reformulated higher-order shear deformation theory and developed the equations of motion of the nonlocal theories to investigate buckling characteristics of nanoplates, considering single-layer graphene sheet, using Navier’s approach. Song et al. [23] studied the buckling, free and forced vibration of multilayer-graphene-nanoplatelets-reinforced plates subjected to axial compression and transverse loadings. Recently, an analytical model considering the flexoelectric effect for unimorph piezoelectric energy harvesters at the nanoscale level was developed by Wang and Wang [24]. Their outcomes revealed that the flexoelectricity shows a significant effect for piezoelectric cantilever beams in the energy harvesting application. Feng et al. [25] examined the nonlinear bending behavior of a unique multilayer composite beams reinforced with graphene platelets that are non-uniformly dispersed in the thickness direction. They found that the bending performance of polymer matrix composite significantly improved by adding a small amount of graphene platelets. Srivastava and Kumar [26] investigated the postbuckling response of the graphene-sheet-reinforced plate comprising the effect of van der Waals bonding amongst the graphene sheet and matrix. By using the FEM-based homogenization technique, the effective elastic constants of the nanocomposite were evaluated. Recently, Lin et al. [27] demonstrated that the nonlinear aeroelastic response and vibration characteristics of the functionally graded (FG) multilayer composite plate reinforced with graphene nanoplatelets, as well as piezoelectric layers when subjected to electromechanical loadings. Naskar and co-authors [28,29,30,31,32] derived new micromechanical analysis for envisaging effective properties of fiber-reinforced composite and functionally graded material by using a novel stochastic representative volume element technique, considering spatial distribution with stochastic analysis and uncertainty quantifications. They have also carried out static and dynamic characteristic including sensitivity and frequency analysis of laminated composite and FG material by using probabilistic and non-probabilistic uncertainty quantification approach. Recently, Kundalwal and co-authors [10,33,34,35] studied the electromechanical response of graphene-based nanocomposite (GNC) and its different structural elements, such as beam, plate, wire and shell, by incorporating piezoelectric graphene nanofiber in polyimide matrix. They predicted the overall effective properties of GNC, using analytical and numerical models. They showed the substantial enrichment in structural response of GNC structures by accounting these size-dependent properties and also revealed that one cannot ignore these effects at nanoscale. They also presented a piezoelastic model which explains the effect of electric field on the effective bending rigidity and strain relaxation concept developed in the nanostructure. Kundalwal et al. [34] studied the stress-transfer characteristics and mechanical properties of multiscale composites, including nano- and micro-scale reinforcements via micromechanical pull out model and molecular dynamic (MD) simulations.

Even if the static and dynamic problem of conventional materials with their different structures were broadly studied by using theoretical, experimental or numerical approach, no previous study has been reported to investigate the electromechanical behavior of advanced layered hybrid piezoelectric composite structures. The effective properties of these advanced hybrid piezoelectric composites are predicted by using analytical and numerical models which are taken from Shingare and Naskar [36]. In this, they showed that electromechanical properties including elastic and piezoelectric properties of hybrid piezoelectric composite are enhanced due to addition of graphene nanofillers. Hence, the objective and the novelty of the current research work is to fill the research gap to study the electromechanical behavior of advanced hybrid graphene-based composite structures. In order to find the reliability of utilizing graphene as nanofiller in epoxy matrix with conventional PZT fibers for various energy harvesters and structural NEMS applications, the research of investigation on the analysis of graphene-based hybrid composite structures is required. Hence, the consideration of graphene nanofillers for modeling of hybrid composite structures is practically more relevant. To best of our knowledge, no work has appeared so far in that direction concerning the response of a hybrid piezoelectric composite plate when subjected to static and dynamic loading by taking flexoelectricity into account. Figure 1 shows the detailed framework and flowchart of our work, which was carried out to achieve the electromechanical behavior of graphene-reinforced piezoelectric composite (GRPC) plate. Hereafter, the paper is presented as per the following sections. In Section 2, an analytical model is derived for simply supported hybrid graphene-based piezoelectric composite plates with consideration of effects of flexoelectricity. In Section 3, the results for hybrid flexoelectric GRPC plates are presented, to investigate the static and modal/dynamic analysis, considering different parameters, such as aspect ratio, plate thickness, different loadings, the combination of different volume fraction of graphene and piezoelectric lead zirconate titanate (PZT) modes, as well as different flexoelectric coefficients. Section 4 offers concluding remarks and future direction.

2. Materials and Methods

In the present article, the mathematical formulation for simply supported (SS) piezoelectric GRPC plate subjected to the uniformly distributed load (UDL ${\mathrm{q}}_0$) is presented, to explore its static and modal analysis. Figure 2 shows the schematic of GRPC plate having thickness ($\mathrm{h}$), length ($\mathrm{a}$) and width (b). The principal material $(\mathrm{x}-\mathrm{y}-\mathrm{z})$ and problem $(1-2-3)$ coordinate system was used to define the composite plate with thickness along the z-axis and the mid-plane of the undeformed plate coincides with the $\mathrm{x}-\mathrm{y}$ plane.

2.1. Assumptions of Kirchhoff’s Plate Theory

• Straight lines normal to the mid-surface (transverse normal) before deformation remain straight after deformation.
• The transverse normal are inextensible.
• The thickness of the plate does not change during a deformation.
• The transverse normal rotate in such a way that they remain normal to the middle surface after deformation.

According to Kirchhoff’s plate theory, the displacement field are given by the following [37]:

$$\mathrm{u}(\mathrm{x},\mathrm{y},\mathrm{z},\mathrm{t})=-\mathrm{z}\frac{\partial \mathrm{w}(\mathrm{x},\mathrm{y},\mathrm{t})}{\partial \mathrm{x}},$$

(1)

$$\mathrm{v}(\mathrm{x},\mathrm{y},\mathrm{z},\mathrm{t})=-\mathrm{z}\frac{\partial \mathrm{w}(\mathrm{x},\mathrm{y},\mathrm{t})}{\partial \mathrm{y}},$$

(2)

$$\mathrm{w}(\mathrm{x},\mathrm{y},\mathrm{z},\mathrm{t})=\mathrm{w}(\mathrm{x},\mathrm{y},\mathrm{t}),$$

(3)

where u, v and w denote the in-plane displacement components corresponding to the x, y and z directions, respectively; $\mathrm{t}$ denotes the time; and $\mathrm{w}$ denotes the transverse displacement of a point on the mid-plane ($\mathrm{z}=0$). The displacement field (Equations (1)–(3)) indicates that straight lines perpendicular to the $\mathrm{x}-\mathrm{y}$ plane before deformation remain straight and normal to the mid-surface after deformation.

Subsequently, the nonzero strains are obtained as follows:

$${\mathsf{\epsilon }}_{\mathrm{x}\mathrm{x}}=-\mathrm{z}\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{x}}^2},{\mathsf{\epsilon }}_{\mathrm{y}\mathrm{y}}=-\mathrm{z}\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{y}}^2},{\mathsf{\epsilon }}_{\mathrm{x}\mathrm{y}}=-\mathrm{z}\frac{{\partial }^2\mathrm{w}}{\partial \mathrm{x}\partial \mathrm{y}},$$

(4)

Due to large ratio of in-plane dimension to thickness, it is considered that the electric field $({\mathrm{E}}_{\mathrm{z}})$ exists only in the $\mathrm{z}$-direction of the piezoelectric medium is constant along the plate thickness and ${\mathrm{E}}_{\mathrm{z}}$ components in the $\mathrm{x}-\mathrm{y}$ plane can be ignored when it is compared with those in the $\mathrm{z}$-direction [38].

The constitutive equations for the dielectric material considering the flexoelectric effect at nanoscale can be written as follows [38]:

$${\mathsf{\sigma }}_{\mathrm{i}\mathrm{j}}=\frac{\partial \mathrm{U}}{\partial {\mathsf{\epsilon }}_{\mathrm{i}\mathrm{j}}}={\mathrm{C}}_{\mathrm{i}\mathrm{j}\mathrm{k}\mathrm{l}}{\mathsf{\epsilon }}_{\mathrm{k}\mathrm{l}}-{\mathrm{e}}_{\mathrm{i}\mathrm{j}\mathrm{k}}{\mathrm{E}}_{\mathrm{k}},$$

(5)

$${\mathsf{\tau }}_{\mathrm{i}\mathrm{j}\mathrm{m}}=\frac{\partial \mathrm{U}}{\partial {\mathsf{\eta }}_{\mathrm{i}\mathrm{j}\mathrm{m}}}=-{\mathrm{f}}_{\mathrm{k}\mathrm{i}\mathrm{j}\mathrm{m}}{\mathrm{E}}_{\mathrm{k}},$$

(6)

$${\mathrm{D}}_{\mathrm{i}}=\frac{\partial \mathrm{U}}{\partial {\mathrm{E}}_{\mathrm{i}}}={\mathrm{e}}_{\mathrm{j}\mathrm{k}\mathrm{i}}{\mathsf{\epsilon }}_{\mathrm{j}\mathrm{k}}+{\mathsf{\chi }}_{\mathrm{i}\mathrm{j}}{\mathrm{E}}_{\mathrm{j}}+{\mathrm{f}}_{\mathrm{i}\mathrm{j}\mathrm{k}\mathrm{l}}{\mathsf{\eta }}_{\mathrm{j}\mathrm{k}\mathrm{l}},$$

(7)

where ${\mathrm{C}}_{\mathrm{i}\mathrm{j}\mathrm{k}\mathrm{l}}$, ${\mathrm{e}}_{\mathrm{i}\mathrm{j}\mathrm{k}}$ and ${\mathsf{\chi }}_{\mathrm{k}\mathrm{l}}$ denote the elastic, piezoelectric and dielectric tensors, respectively; ${\mathrm{f}}_{\mathrm{i}\mathrm{j}\mathrm{k}\mathrm{l}}$ denotes the fourth-order tensor; ${\mathsf{\sigma }}_{\mathrm{i}\mathrm{j}}$ denotes the stress tensor; ${\mathrm{D}}_{\mathrm{i}}$, ${\mathrm{E}}_{\mathrm{i}}$, ${\mathsf{\tau }}_{\mathrm{i}\mathrm{j}\mathrm{m}}$ and ${\mathsf{\eta }}_{\mathrm{j}\mathrm{k}\mathrm{l}}$ denote the electric displacement, electric field vector, the higher-order stress tensor and strain gradient tensor, respectively.

The strain gradient tensor can be given as follows:

$${\mathsf{\eta }}_{\mathrm{j}\mathrm{k}\mathrm{l}}={\mathsf{\epsilon }}_{\mathrm{j}\mathrm{k},\mathrm{l}},$$

(8)

Making use of Equations (5)–(7), the constitutive relations are re-written as follows:

$${\mathsf{\sigma }}_{\mathrm{x}\mathrm{x}}={\mathrm{C}}_{11}{\mathsf{\epsilon }}_{\mathrm{x}\mathrm{x}}+{\mathrm{C}}_{12}{\mathsf{\epsilon }}_{\mathrm{y}\mathrm{y}}-{\mathrm{e}}_{31}{\mathrm{E}}_{\mathrm{z}},$$

(9)

$${\mathsf{\sigma }}_{\mathrm{y}\mathrm{y}}={\mathrm{C}}_{12}{\mathsf{\epsilon }}_{\mathrm{x}\mathrm{x}}+{\mathrm{C}}_{11}{\mathsf{\epsilon }}_{\mathrm{y}\mathrm{y}}-{\mathrm{e}}_{31}{\mathrm{E}}_{\mathrm{z}},$$

(10)

$${\mathsf{\tau }}_{\mathrm{x}\mathrm{y}}=2{\mathrm{C}}_{66}{\mathsf{\epsilon }}_{\mathrm{x}\mathrm{y}},$$

(11)

$${\mathsf{\tau }}_{\mathrm{x}\mathrm{x}\mathrm{z}}=-{\mathrm{f}}_{14}{\mathrm{E}}_{\mathrm{z}},$$

(12)

$${\mathsf{\tau }}_{\mathrm{y}\mathrm{y}\mathrm{z}}=-{\mathrm{f}}_{14}{\mathrm{E}}_{\mathrm{z}},$$

(13)

$${\mathrm{D}}_{\mathrm{z}}={\mathrm{e}}_{31}({\mathsf{\epsilon }}_{\mathrm{x}\mathrm{x}}+{\mathsf{\epsilon }}_{\mathrm{y}\mathrm{y}})+{\mathsf{\chi }}_{33}{\mathrm{E}}_{\mathrm{z}}+{\mathsf{\epsilon }}_{\mathrm{x}\mathrm{y}}({\mathsf{\eta }}_{\mathrm{x}\mathrm{x}\mathrm{z}}+{\mathsf{\eta }}_{\mathrm{y}\mathrm{y}\mathrm{z}}),$$

(14)

where ${\mathrm{f}}_{14}={\mathrm{f}}_{3113}={\mathrm{f}}_{3223}$ [39]. To simplify the mathematical formulation, other than ${\mathsf{\eta }}_{\mathrm{x}\mathrm{x}\mathrm{z}} \mathrm{a}\mathrm{n}\mathrm{d} {\mathsf{\eta }}_{\mathrm{y}\mathrm{y}\mathrm{z}}$ all strain gradients are considered to be zero, meanwhile the related strain gradients are much lesser that of those along the thickness direction of plate.

Because of the absence of external ${\mathrm{E}}_{\mathrm{z}}$, ${\mathrm{D}}_{\mathrm{z}}$ becomes equal to electric polarization. Therefore, it is undoubtedly seen that term ${\mathsf{\epsilon }}_{\mathrm{x}\mathrm{y}}({\mathsf{\eta }}_{\mathrm{x}\mathrm{x}\mathrm{z}}+{\mathsf{\eta }}_{\mathrm{y}\mathrm{y}\mathrm{z}})$ in Equation (14) signifies the polarization generated in the plate because of the inhomogeneous strain.

In the absence of free electric charge $({\mathsf{\rho }}_{\mathrm{s}})$, ${\mathrm{D}}_{\mathrm{z}}$ for the thin plate must satisfy Gauss’s law of electrostatics $(\nabla .\mathrm{D}={\mathsf{\rho }}_{\mathrm{s}})$ and can be expressed as follows:

$$\frac{\partial {\mathrm{D}}_{\mathrm{z}}}{\partial \mathrm{z}}=0,$$

(15)

For open-circuit case, on the surface of plate ${\mathrm{D}}_{\mathrm{z}}$ is zero. Hence, the internal ${\mathrm{E}}_{\mathrm{z}}$ can be derived from Equation (15) as follows:

$${\mathrm{E}}_{\mathrm{z}}=\frac{{\mathrm{e}}_{31}}{{\mathsf{\chi }}_{33}}\left(\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{x}}^2}+\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{y}}^2}\right)\mathrm{z}+\frac{{\mathrm{f}}_{14}}{{\mathsf{\chi }}_{33}}\left(\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{x}}^2}+\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{y}}^2}\right),$$

(16)

It is observed from Equation (16) that the term ${\mathrm{e}}_{31}/{\mathsf{\chi }}_{33}$ implies ${\mathrm{E}}_{\mathrm{z}}$ which is caused due to the applied of homogeneous strains while term ${\mathrm{f}}_{14}/{\mathsf{\chi }}_{33}$ implies ${\mathrm{E}}_{\mathrm{z}}$ which is caused due to the applied of inhomogeneous strains. In case of non-piezoelectric dielectrics, inhomogeneous strains-generated ${\mathrm{E}}_{\mathrm{z}}$ predicts structural response which is important for considerate basis of electromechanical coupling behaviors of dielectrics. Considering the flexoelectric effect, the piezoelectric effect accompanying internal ${\mathrm{E}}_{\mathrm{z}}$ in thickness direction of plate no longer remnants anti-symmetric corresponding to the middle plane. Above Equation (16) can be re-expressed as ${\mathrm{E}}_{\mathrm{z}}=({\mathrm{e}}_{31}\mathrm{z}+{\mathrm{f}}_{14})\mathrm{G}/{\mathsf{\chi }}_{33}$ with $\mathrm{G}=\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{x}}^2}+\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{y}}^2}$ by taking the summation of the curvatures at an any point in the plate. It can be noticed that the response of ${\mathrm{E}}_{\mathrm{z}}$ to G is deeply depends on the z-coordinate. Moreover, ${\mathrm{e}}_{31}\mathrm{z}+{\mathrm{f}}_{14}$ depends not only on the plate thickness but also on its flexoelectric coefficient. If the effect of piezoelectricity is not considered, then present solution attributes to the effect of flexoelectricity.

The governing equations for the plate problem can be derived, using dynamic version of Hamilton’s variational principle, as follows [40]:

$$\mathsf{\delta }{\displaystyle {\int }_0^{\mathrm{t}}\left(\mathrm{U}+\mathrm{W}-\mathrm{K}\right)}\mathrm{d}\mathrm{t}=0,$$

(17)

where $\mathrm{U}$ is the electric Gibbs free energy density. For the simplification of formulation, we considered open-circuit condition $({\mathrm{D}}_{\mathrm{z}}=0)$ and hence, the relation of $\mathrm{U}$ can be reduced as follows:

$$\mathrm{U}=\frac{1}{2}{\mathsf{\sigma }}_{\mathrm{i}\mathrm{j}}{\mathsf{\epsilon }}_{\mathrm{i}\mathrm{j}}+\frac{1}{2}{\mathsf{\tau }}_{\mathrm{i}\mathrm{j}\mathrm{k}}{\mathsf{\eta }}_{\mathrm{i}\mathrm{j}\mathrm{k}},$$

(18)

By ignoring the vibration along $\mathrm{x}-\mathrm{y}$ plane, the kinetic energy $(\mathrm{K})$ is described as follows:

$$\mathrm{K}=\frac{1}{2}{\displaystyle \underset{\mathsf{\Omega }}{\int }\mathsf{\rho }{\left(\frac{\partial \mathrm{w}}{\partial \mathrm{t}}\right)}^2\mathrm{d}\mathsf{\Omega }},$$

(19)

where $\mathsf{\Omega }$ denotes the total volume occupied by the GRPC plate, and $\mathsf{\rho }$ is the mass density. Due to external load, the work done $(\mathrm{W})$ can be described as follows:

$$\mathrm{W}={\displaystyle \underset{0}{\overset{\mathrm{a}}{\int }}{\displaystyle \underset{0}{\overset{\mathrm{b}}{\int }}{\mathrm{q}}_0\mathrm{w}\mathrm{d}\mathrm{y}\mathrm{d}\mathrm{x}}},$$

(20)

It is significant to note that the variational principles are the fundamental basis of investigating and determining coupled multi-physics problems in multifaceted materials. Thus, the governing Eq. can be written as follows:

$$\frac{{\partial }^2{\mathrm{M}}_{\mathrm{x}\mathrm{x}}}{\partial {\mathrm{x}}^2}+\frac{{\partial }^2{\mathrm{M}}_{\mathrm{x}\mathrm{y}}}{\partial \mathrm{x}\partial \mathrm{y}}+\frac{{\partial }^2{\mathrm{M}}_{\mathrm{y}\mathrm{x}}}{\partial \mathrm{x}\partial \mathrm{y}}+\frac{{\partial }^2{\mathrm{M}}_{\mathrm{y}\mathrm{y}}}{\partial {\mathrm{y}}^2}+\frac{{\partial }^2{\mathrm{N}}_{\mathrm{x}\mathrm{x}\mathrm{z}}}{\partial {\mathrm{x}}^2}+\frac{{\partial }^2{\mathrm{N}}_{\mathrm{y}\mathrm{y}\mathrm{z}}}{\partial {\mathrm{y}}^2}-\mathsf{\rho }\mathrm{h}\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{t}}^2}+{\mathrm{q}}_0=0,$$

(21)

Boundary conditions prescribed on all four edges of SS rectangular GRPC plates can be given by the following [40]:

$\mathrm{a}\mathrm{t}\mathrm{x}=0\mathrm{and}\mathrm{x}=\mathrm{a}:$

$$\mathrm{w}=0,{\mathrm{M}}_{\mathrm{x}\mathrm{x}}=0$$

(22)

$\mathrm{a}\mathrm{t}\mathrm{y}=0\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{y}=\mathrm{b}:$

$$\mathrm{w}=0,{\mathrm{M}}_{\mathrm{y}\mathrm{y}}=0$$

(23)

where ${\mathrm{M}}_{\mathrm{x}\mathrm{x}},{\mathrm{M}}_{\mathrm{x}\mathrm{y}},{\mathrm{M}}_{\mathrm{y}\mathrm{x}}\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{M}}_{\mathrm{y}\mathrm{y}}$ indicate the bending moments, and ${\mathrm{N}}_{\mathrm{x}\mathrm{x}\mathrm{z}}\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{N}}_{\mathrm{y}\mathrm{y}\mathrm{z}}$ denote the axial higher-order stresses along the thickness, which can be obtained as follows:

$${\mathrm{M}}_{\mathrm{x}\mathrm{x}}={\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}{\mathsf{\sigma }}_{\mathrm{x}\mathrm{x}}}\mathrm{z}\mathrm{d}\mathrm{z},{\mathrm{M}}_{\mathrm{y}\mathrm{y}}={\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}{\mathsf{\sigma }}_{\mathrm{y}\mathrm{y}}}\mathrm{z}\mathrm{d}\mathrm{z},{\mathrm{M}}_{\mathrm{x}\mathrm{y}}={\mathrm{M}}_{\mathrm{y}\mathrm{x}}={\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}{\mathsf{\tau }}_{\mathrm{x}\mathrm{y}}}\mathrm{z}\mathrm{d}\mathrm{z},$$

(24)

$${\mathrm{N}}_{\mathrm{x}\mathrm{x}\mathrm{z}}={\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}{\mathsf{\tau }}_{\mathrm{x}\mathrm{x}\mathrm{z}}}\mathrm{d}\mathrm{z},{\mathrm{N}}_{\mathrm{y}\mathrm{y}\mathrm{z}}={\displaystyle \underset{-\mathrm{h}/2}{\overset{\mathrm{h}/2}{\int }}{\mathsf{\tau }}_{\mathrm{y}\mathrm{y}\mathrm{z}}}\mathrm{d}\mathrm{z},$$

(25)

By putting Equations (4) and (16) into Equations (9)–(13), one can get following explicit expressions for the stresses related to the transverse deflection.

$${\mathsf{\sigma }}_{\mathrm{x}\mathrm{x}}=-\left({\mathrm{C}}_{11}+\frac{{\mathrm{e}}_{31}^2}{{\mathsf{\chi }}_{33}}\right)\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{x}}^2}\mathrm{z}-\left({\mathrm{C}}_{12}+\frac{{\mathrm{e}}_{31}^2}{{\mathsf{\chi }}_{33}}\right)\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{y}}^2}\mathrm{z}-\frac{{\mathrm{e}}_{31}{\mathrm{f}}_{14}}{{\mathsf{\chi }}_{33}}\left(\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{x}}^2}+\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{y}}^2}\right),$$

(26)

$${\mathsf{\sigma }}_{\mathrm{y}\mathrm{y}}=-\left({\mathrm{C}}_{12}+\frac{{\mathrm{e}}_{31}^2}{{\mathsf{\chi }}_{33}}\right)\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{x}}^2}\mathrm{z}-\left({\mathrm{C}}_{11}+\frac{{\mathrm{e}}_{31}^2}{{\mathsf{\chi }}_{33}}\right)\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{y}}^2}\mathrm{z}-\frac{{\mathrm{e}}_{31}{\mathrm{f}}_{14}}{{\mathsf{\chi }}_{33}}\left(\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{x}}^2}+\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{y}}^2}\right),$$

(27)

$${\mathsf{\tau }}_{\mathrm{x}\mathrm{y}}=-2{\mathrm{C}}_{66}\frac{{\partial }^2\mathrm{w}}{\partial \mathrm{x}\partial \mathrm{y}}\mathrm{z},$$

(28)

$${\mathsf{\tau }}_{\mathrm{x}\mathrm{x}\mathrm{z}}={\mathsf{\tau }}_{\mathrm{y}\mathrm{y}\mathrm{z}}=-\frac{{\mathrm{e}}_{31}{\mathrm{f}}_{14}}{{\mathsf{\chi }}_{33}}\left(\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{x}}^2}+\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{y}}^2}\right)\mathrm{z}-\frac{{\mathrm{f}}_{14}^2}{{\mathsf{\chi }}_{33}}\left(\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{x}}^2}+\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{y}}^2}\right),$$

(29)

The bending moments in terms of transverse deflection can be obtained by using Equations (24)–(27) in Equations (24) and (25), as follows:

$${\mathrm{M}}_{\mathrm{x}\mathrm{x}}=-\left({\mathrm{C}}_{11}+\frac{{\mathrm{e}}_{31}^2}{{\mathsf{\chi }}_{33}}\right)\frac{{\mathrm{h}}^3}{12}\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{x}}^2}\mathrm{z}-\left({\mathrm{C}}_{12}+\frac{{\mathrm{e}}_{31}^2}{{\mathsf{\chi }}_{33}}\right)\frac{{\mathrm{h}}^3}{12}\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{y}}^2},$$

(30)

$${\mathrm{M}}_{\mathrm{y}\mathrm{y}}=-\left({\mathrm{C}}_{12}+\frac{{\mathrm{e}}_{31}^2}{{\mathsf{\chi }}_{33}}\right)\frac{{\mathrm{h}}^3}{12}\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{x}}^2}-\left({\mathrm{C}}_{11}+\frac{{\mathrm{e}}_{31}^2}{{\mathsf{\chi }}_{33}}\right)\frac{{\mathrm{h}}^3}{12}\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{y}}^2},$$

(31)

$${\mathrm{M}}_{\mathrm{x}\mathrm{y}}={\mathrm{M}}_{\mathrm{x}\mathrm{y}}=-2{\mathrm{C}}_{66}\frac{{\mathrm{h}}^3}{12}\frac{{\partial }^2\mathrm{w}}{\partial \mathrm{x}\partial \mathrm{y}}\mathrm{z},$$

(32)

$${\mathrm{N}}_{\mathrm{x}\mathrm{x}\mathrm{z}}={\mathrm{N}}_{\mathrm{y}\mathrm{y}\mathrm{z}}=-\frac{{\mathrm{f}}_{14}^2}{{\mathsf{\chi }}_{33}}\left(\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{x}}^2}+\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{y}}^2}\right)\mathrm{h},$$

(33)

By observing Equations (26)–(33), it is obvious that the flexoelectric effect has significant influence on the distribution of stress and higher-order stresses. Therefore, the higher-order stress disappears as the effect of flexoelectricity is not considered, while it does not have any influence on conservative bending moments. Furthermore, the higher-order stresses are raised by the introduction of flexoelectric effect.

Making use of Equations (30)–(33) into Equation (21), one can find the following governing Equation (34) in terms of transverse deflection.

$${\mathrm{D}}_{11}\left(\frac{{\partial }^4\mathrm{w}}{\partial {\mathrm{x}}^4}+\frac{{\partial }^4\mathrm{w}}{\partial {\mathrm{y}}^4}\right)+2\left({\mathrm{D}}_{12}+2{\mathrm{D}}_{66}\right)\frac{{\partial }^4\mathrm{w}}{\partial {\mathrm{x}}^2\partial {\mathrm{y}}^2}+\mathsf{\rho }\mathrm{h}\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{t}}^2}={\mathrm{q}}_0,$$

(34)

where

$$\{\begin{array}{c}{\mathrm{D}}_{11}=\left({\mathrm{C}}_{11}+\frac{{\mathrm{e}}_{31}^2}{{\mathsf{\chi }}_{33}}\right)\frac{{\mathrm{h}}^3}{12}+\frac{{\mathrm{f}}_{14}^2}{{\mathsf{\chi }}_{33}}\mathrm{h}\\ {\mathrm{D}}_{12}=\left({\mathrm{C}}_{12}+\frac{{\mathrm{e}}_{31}^2}{{\mathsf{\chi }}_{33}}\right)\frac{{\mathrm{h}}^3}{12}+\frac{{\mathrm{f}}_{14}^2}{{\mathsf{\chi }}_{33}}\mathrm{h}\\ {\mathrm{D}}_{66}={\mathrm{C}}_{66}\frac{{\mathrm{h}}^3}{12}\end{array},$$

(35)

From Equation (35), it is obvious that the bending rigidity of plate is substantially influenced by effect of flexoelectricity and thickness of plate.

2.2. Closed-Form Solution for Static Analysis of GRPC Plates

To obtain the closed-form solution for static deflection of hybrid GRPC plates, the governing Equation (34) can be re-expressed as follows [40]:

$${\mathrm{D}}_{11}\left(\frac{{\partial }^4\mathrm{w}}{\partial {\mathrm{x}}^4}+\frac{{\partial }^4\mathrm{w}}{\partial {\mathrm{y}}^4}\right)+2\left({\mathrm{D}}_{12}+2{\mathrm{D}}_{66}\right)\frac{{\partial }^4\mathrm{w}}{\partial {\mathrm{x}}^2\partial {\mathrm{y}}^2}={\mathrm{q}}_0,$$

(36)

It can be noticed that the above equation follows the conventional classical Kirchhoff plate theory, considering piezoelectric effect only if we set the flexoelectric effect to zero.

To solve Equation (36) of the SS GRPC plate, the following Fourier series is used to compute $\mathrm{w}(\mathrm{x},\mathrm{y})$ according to conventional classical Kirchhoff plate theory.

$$\mathrm{w}(\mathrm{x},\mathrm{y})={\displaystyle \sum _{\mathrm{m}=1}^{\infty }{\displaystyle \sum _{\mathrm{n}=1}^{\infty }{\mathrm{A}}_{\mathrm{m}\mathrm{n}}}}\sin \mathsf{\alpha }\mathrm{x}+\sin \mathsf{\beta }\mathrm{y},$$

(37)

where $\mathsf{\alpha }=\frac{\mathrm{m}\mathsf{\pi }}{\mathrm{a}},\mathsf{\beta }=\frac{\mathrm{n}\mathsf{\pi }}{\mathrm{b}}\mathrm{a}\mathrm{n}\mathrm{d}{\mathrm{A}}_{\mathrm{m}\mathrm{n}}$ denote the coefficients to be evaluated for each $\mathrm{m}$ and $\mathrm{n}.$ It was already verified that Equation (37) satisfies the boundary conditions in Equations (22) and (23). The uniformly distributed load ${\mathrm{q}}_0(\mathrm{x},\mathrm{y})$ can be represented as a Fourier series:

$$\mathrm{q}(\mathrm{x},\mathrm{y})={\displaystyle \sum _{\mathrm{m}=1}^{\infty }{\displaystyle \sum _{\mathrm{n}=1}^{\infty }{\mathrm{Q}}_{\mathrm{m}\mathrm{n}}}}\sin \mathsf{\alpha }\mathrm{x}+\sin \mathsf{\beta }\mathrm{y},$$

(38)

$$\mathrm{w}(\mathrm{x},\mathrm{y})=\frac{{\mathrm{Q}}_{\mathrm{m}\mathrm{n}}}{{\mathrm{d}}_{\mathrm{m}\mathrm{n}}}\mathrm{a}\mathrm{n}\mathrm{d},$$

(39)

where ${\mathrm{Q}}_{\mathrm{m}\mathrm{n}}$ denotes the load coefficient for respective loading and can be expressed as follows:

For (i) uniform distributed load:

$${\mathrm{Q}}_{\mathrm{m}\mathrm{n}}=\frac{16{\mathrm{q}}_0}{\mathrm{m}\mathrm{n}{\mathsf{\pi }}^2}\mathrm{m},\mathrm{n}=1,3,5,\dots ,$$

(40)

(ii) Linearly varying load:

$${\mathrm{Q}}_{\mathrm{m}\mathrm{n}}=\frac{8{\mathrm{q}}_0\cos \mathrm{m}\mathsf{\pi }}{\mathrm{m}\mathrm{n}{\mathsf{\pi }}^2}\mathrm{m},\mathrm{n}=1,3,5,\dots ,$$

(41)

(iii) Point load:

$${\mathrm{Q}}_{\mathrm{m}\mathrm{n}}=\frac{4\mathrm{P}}{\mathrm{a}\mathrm{b}}\sin \frac{\mathrm{m}\mathsf{\pi }{\mathrm{x}}_0}{\mathrm{a}}\sin \frac{\mathrm{m}\mathsf{\pi }{\mathrm{y}}_0}{\mathrm{b}}\mathrm{m},\mathrm{n}=1,2,3,\dots ,$$

(42)

(iv) Line load:

$${\mathrm{Q}}_{\mathrm{m}\mathrm{n}}=\frac{8{\mathrm{q}}_0}{\mathsf{\pi }\mathrm{a}\mathrm{n}}\sin \frac{\mathrm{m}\mathsf{\pi }{\mathrm{x}}_0}{\mathrm{a}}\mathrm{m},=1,3,5,\dots ;\mathrm{n}=1,2,3,\dots ,$$

(43)

$${\mathrm{d}}_{\mathrm{m}\mathrm{n}}={\left(\frac{\mathsf{\pi }}{\mathrm{b}}\right)}^4\left\{{\mathrm{D}}_{11}{\left(\frac{\mathrm{m}}{\mathrm{a}}\right)}^4+{\mathrm{D}}_{12}{\left(\frac{\mathrm{n}}{\mathrm{b}}\right)}^4+2({\mathrm{D}}_{12}+{\mathrm{D}}_{66}){\left(\frac{\mathrm{m}\mathrm{n}}{\mathrm{a}\mathrm{b}}\right)}^2\right\},$$

(44)

Making use of Equations (37) and (38) into Equation (36), one can derive the transverse deflection of the SS plate.

$$\begin{array}{l}\mathrm{w}(\mathrm{x},\mathrm{y})=\frac{16{\mathrm{q}}_0}{{\mathsf{\pi }}^6}\\ ={\displaystyle \sum _{\mathrm{m}=1,3,5\dots }^{\infty }{\displaystyle \sum _{\mathrm{n}=1,3,5\dots }^{\infty }\frac{\sin \mathsf{\alpha }\mathrm{x}\sin \mathsf{\beta }\mathrm{y}}{\mathrm{m}\mathrm{n}\left\{{\mathrm{D}}_{11}{\left(\frac{\mathrm{m}}{\mathrm{a}}\right)}^4+{\mathrm{D}}_{12}{\left(\frac{\mathrm{n}}{\mathrm{b}}\right)}^4+2({\mathrm{D}}_{12}+{\mathrm{D}}_{66}){\left(\frac{\mathrm{m}\mathrm{n}}{\mathrm{a}\mathrm{b}}\right)}^2\right\}}}}\end{array}$$

(45)

2.3. Closed-Form Solution for Modal Analysis Considering Free Vibration of GRPC Plates

For the free vibration of GRPC plate, the governing equation can be obtained by using Equation (34):

$${\mathrm{D}}_{11}\left(\frac{{\partial }^4\mathrm{w}}{\partial {\mathrm{x}}^4}+\frac{{\partial }^4\mathrm{w}}{\partial {\mathrm{y}}^4}\right)+2\left({\mathrm{D}}_{12}+2{\mathrm{D}}_{66}\right)\frac{{\partial }^4\mathrm{w}}{\partial {\mathrm{x}}^2\partial {\mathrm{y}}^2}+\mathsf{\rho }\mathrm{h}\frac{{\partial }^2\mathrm{w}}{\partial {\mathrm{t}}^2}=0,$$

(46)

Similar to classical Kirchhoff plate theory, the harmonic solution for $\mathrm{w}(\mathrm{x},\mathrm{y},\mathrm{t})$ can be obtained as follows:

$$\mathrm{w}(\mathrm{x},\mathrm{y},\mathrm{t})={\displaystyle \sum _{\mathrm{m}=1}^{\infty }{\displaystyle \sum _{\mathrm{n}=1}^{\infty }{\mathrm{B}}_{\mathrm{m}\mathrm{n}}}}\sin \frac{\mathrm{m}\mathsf{\pi }\mathrm{x}}{\mathrm{a}}\sin \frac{\mathrm{m}\mathsf{\pi }\mathrm{x}}{\mathrm{b}}{\mathrm{e}}^{\mathrm{i}{\mathrm{w}}_{\mathrm{m}\mathrm{n}}\mathrm{t}},$$

(47)

where ${\mathrm{B}}_{\mathrm{m}\mathrm{n}}$ denotes a constant signifying the mode shape; $\mathrm{m}$ and $\mathrm{n}$ denote the half wave numbers; ${\mathrm{w}}_{\mathrm{m}\mathrm{n}}$ denotes the resonant frequency; and $\mathrm{i}=\sqrt{-1}$.

By substituting Equation (47) into Equation (46), we get the resonant frequency of plate.

$$\left\{{\mathrm{D}}_{11}\left[{\left(\frac{\mathrm{m}\mathsf{\pi }}{\mathrm{a}}\right)}^4+{\left(\frac{\mathrm{n}\mathsf{\pi }}{\mathrm{b}}\right)}^4\right]+2({\mathrm{D}}_{12}+2{\mathrm{D}}_{66}){\left(\frac{\mathrm{m}\mathsf{\pi }}{\mathrm{a}}\right)}^2{\left(\frac{\mathrm{n}\mathsf{\pi }}{\mathrm{b}}\right)}^2\right\}-\mathsf{\rho }\mathrm{h}{\mathrm{w}}_{\mathrm{m}\mathrm{n}}^2=0,$$

(48)

Hereafter, for different order numbers $\mathrm{m}$ and $\mathrm{n}$, the resonant frequency can be obtained as follows:

$${\mathrm{w}}_{\mathrm{m}\mathrm{n}}={\left(\frac{\mathsf{\pi }}{\mathsf{\rho }\mathrm{h}}\right)}^2\sqrt{{\mathrm{D}}_{11}\left[{\left(\frac{\mathrm{m}}{\mathrm{a}}\right)}^4+{\left(\frac{\mathrm{n}}{\mathrm{b}}\right)}^4\right]+2({\mathrm{D}}_{12}+2{\mathrm{D}}_{66}){\left(\frac{\mathrm{m}\mathsf{\pi }}{\mathrm{a}}\right)}^2{\left(\frac{\mathrm{n}\mathsf{\pi }}{\mathrm{b}}\right)}^2},$$

(49)

From above expressions, it is clearly seen that the resonant frequency can be evaluated using the traditional piezoelectric theory if the effect of flexoelectricity is eliminated. To demonstrate the effect of flexoelectricity on the free vibration characteristics, we can consider the different value of modes with respect to the aspect ratio and thickness of plate, as presented in the subsequent section, Section 3.

3. Results

In this article, hybrid graphene reinforced piezoelectric composite (GRPC) plate with combination of PZT and graphene $({\mathrm{v}}_{\mathrm{g}}=0.2{\mathrm{v}}_{\mathrm{p}})$ into epoxy matrix was selected to investigate how the flexoelectricity influences the static and modal analysis, in which ${\mathrm{v}}_{\mathrm{g}}$ and ${\mathrm{v}}_{\mathrm{p}}$ denote the volume fraction of graphene nanofillers and PZT fiber. The elastic, piezoelectric and dielectric properties of GRPC are taken from Shingare and Naskar [36] and are summarized in

Table 1. Effective properties of hybrid GRPC $({\mathrm{v}}_{\mathrm{g}}=0.2{\mathrm{v}}_{\mathrm{p}})$.

Material	${\mathbf{C}}_{11}(\mathbf{G}\mathbf{P}\mathbf{a})$	${\mathbf{C}}_{12}(\mathbf{G}\mathbf{P}\mathbf{a})$	${\mathbf{C}}_{66}(\mathbf{G}\mathbf{P}\mathbf{a})$	${\mathbf{e}}_{31}(\mathbf{C}/{\mathbf{m}}^2)$	${\mathbf{\chi }}_{33}\times {10}^{-9}(\mathbf{F}/\mathbf{m})$
GRPC	112.43	3.34	2.03	−6.9337	3.264

. The mass density $(\mathsf{\rho })$ of the hybrid GRPC is considered as 3798.4 $\mathrm{k}\mathrm{g}/{\mathrm{m}}^3$. All the results are estimated for load $\mathrm{q}=0.1\mathrm{M}\mathrm{P}\mathrm{a}.$ From the above mathematical formulation, it is observed that the transverse deflection of hybrid GRPC plate is directly proportional to the applied load and independent of the mass density, while the dynamic/modal response is independent of the applied load but proportional to the mass density of the material.

3.1. Static Deflection of Hybrid Flexoelectric GRPC Plate

In case of static analysis of hybrid GRPC plate, we considered the bending rigidity, center deflection and maximum deflection of hybrid plate with respect to different parameters, such as different graphene contents, flexoelectric coefficient, thickness and aspect ratio of plate. Figure 2 shows the variation of center deflection against the aspect ratio of hybrid GRPC plate, considering the different flexoelectric coefficient for $\mathrm{m}=1$ and $\mathrm{n}=1$. From Figure 3, it can be observed that the center deflection of plate decreases as the value of flexoelectric coefficient increases. Due to lack of atomistic calculations or experiments, the flexoelectric coefficients of hybrid GRPC are not readily available. However, in experimental measurement, it is observed that, for certain crystals and ceramics, the flexoelectric coefficient is estimated to have the order of $\approx {10}^{-7} \mathrm{C}/\mathrm{m}$ [41,42]. Therefore, the flexoelectric coefficient for hybrid GRPC is considered as ${10}^{-7} \mathrm{C}/\mathrm{m}$ for further calculation.

.

Figure 4 and

Table 2. Three-dimensional representation of deflection of hybrid flexoelectric plate, with and without considering flexoelectricity for different values of m and n.

(m, n)	Without Flexoelectric Effect	With Flexoelectric Effect
$\mathrm{m}=1$; $\mathrm{n}=1$
$\mathrm{m}=3$; $\mathrm{n}=3$

illustrate the effect of different values of $\mathrm{m}$ and $\mathrm{n}$ on the deflection of hybrid GRPC plate, with and without considering flexoelectricity. From this, it can be observed that the deflection of the hybrid GRPC plate decreases as the value of $\mathrm{m}$ and $\mathrm{n}$ increases. Additionally, Table 2 shows the 3D representation of the deflection of the hybrid GRPC plate for a clear visualization of the deflection of the plate.

Figure 5 shows the variation of center deflection against the aspect ratio of hybrid GRPC plate, considering the different combination of graphene and PZT for $\mathrm{m}=3$ and $\mathrm{n}=3$. From this, in Figure 5, it can be observed that the center deflection of the plate decreases as the value of graphene content in the hybrid composite plate increases. According to the experimental studies, it is clear that the appropriate graphene volume fraction in the composite should be less than 10%. Therefore, in this study, the volume fraction of graphene is considered as $0.2$ times the volume fraction of PZT (i.e., 8% in composite). In view of experimental possibilities, we considered $({\mathrm{v}}_{\mathrm{g}}=0.2{\mathrm{v}}_{\mathrm{p}})$ for further calculation of static and modal analysis, which consists of a graphene content of less than 10% volume fraction in the hybrid GRPC plate. Figure 6 summarizes the values of maximum deflections of hybrid GRPC plates for different in-plane dimensions $(\mathrm{a}=\mathrm{b}=45\mathrm{h};55\mathrm{h})$. It can be clearly seen that the influence of flexoelectricity on the maximum deflection of a plate diminishes as its thickness increases, indicating that the flexoelectric effect is size dependent.

However, it can be concluded from the previous figures that the deflection of the plate with the flexoelectric effect is lower that of the conventional plate theory, without considering the flexoelectric effect for both cases of in-plane dimension. The variation between both cases is significant, which indicates the importance of flexoelectricity in the static bending of the hybrid GRPC plate.

So far, to explore the effect of flexoelectricity and value of different values of m and n on the static deflection of the hybrid GRPC plate, four discrete types of loading conditions, namely point load (P), inline, uniformly distributed load (UDL) and uniformly varying load (UVL, i.e., hydrostatic load), are used. These cases represent the practical situation of different types of loadings on the thin plate. We considered the equivalent magnitude of loading in all situations. Our selection of in-plane dimensions of the plate is based on the fact that the theory of Kirchhoff plate provides better results when the aspect ratio of a plate is in the range of 5–80 [43]. As expected, the maximum deflection of the plate occurs at its center, irrespective of the type of loading in both the cases (with and without flexoelectricity). It is evident from this table that the consideration of the flexoelectric effect results in the lowering of the deflection of the hybrid GRPC plates, as compared to that of conventional plates. The significant reductions in static deflection of hybrid GRPC plates with respect to the value of m and n are observed from

Table 3. Deflection hybrid flexoelectric plate, with and without considering flexoelectricity for different loadings $(\mathrm{m}=1,3;\mathrm{n}=1,3).$

Load	$\mathrm{m}=1;\mathrm{n}=1$	$\mathrm{m}=3;\mathrm{n}=3$
Point
In-line
UDL
UVL

and

Table 4. Deflection hybrid flexoelectric plate, with and without considering flexoelectricity for different loadings $(\mathrm{m}=5,7;\mathrm{n}=5,7).$

Load	$\mathrm{m}=5;\mathrm{n}=5$	$\mathrm{m}=7;\mathrm{n}=7$
Point
In-line
UDL
UVL

. From Table 3 and Table 4, it can be clearly seen that there is a significant influence of flexoelectricity on the center deflection of the plate. In this, the legend or graph with the blue color indicates the results for the hybrid flexoelectric plate, while the graph with the red color indicates the results for the conventional plate. Here, we considered different values of m and n, such as $\mathrm{m}=1,3,5,7\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{n}=1,3,5,7$, and predicted the results for different loading conditions. It is evident from these tables that the consideration of the flexoelectric effect and different values of m and n results in the lowering the deflection of hybrid GRPC plates, as compared to that of conventional plates (without considering flexoelectricity). For example, it may be observed from this table that the reductions in static deflections of hybrid GRPC plates are found to be ∼6.0%, corresponding to 4 nm thickness of nanoplate in all loading cases. In the case of the deflection of plate, point load condition gives maximum deflection of the SS hybrid GRPC plate when compared with other loadings (point load > In-line load > UDL > UVL), accounting for both cases (with and without flexoelectricity).

Figure 7 illustrates the variation of normalized bending rigidity $({\mathrm{D}}_{11}/{\mathrm{D}}_{11}^0)$ of the hybrid flexoelectric GRPC plate with respect to its thickness $(\mathrm{h})$, where ${\mathrm{D}}_{11}$ and ${\mathrm{D}}_{11}^0$ are the bending rigidities of hybrid flexoelectric plate, with and without effect of flexoelectricity, respectively. Note that the value of ${\mathrm{D}}_{11}$ depends only on the thickness of plate and is independent of its in-plane dimensions, as seen from Equation (32). It is clearly observed from Figure 6 that the effect of flexoelectricity on the normalized bending rigidity is size-dependent, and it is more eminent for the plate with a smaller thickness. For example, in regard to the bending rigidity, the effect of flexoelectricity is ~3.8 times the bending rigidity, without considering the effect of flexoelectricity when the plate thickness is considered as 10 nm. If we consider both flexoelectric and surface effects, then the bending rigidity becomes ~4.5 times the bending rigidity, without considering the effect of flexoelectricity when the plate thickness is considered as 10 nm. Effective bending rigidity of the plate, considering both flexoelectric and surface effects, is determined by referring to Reference [8]. However, for the sake of brevity, results with surface effects are not considered for the remaining results. This difference is noticeable, and it cannot be ignored for predicting the electromechanical coupling responses of nanostructures accurately. As the thickness of the plate increases, the effect of flexoelectricity diminishes promptly, and this effect vanishes as the normalized bending rigidity approaches the unity. It can be noted from the above discussion that the strain gradient from the constitutive equations related to flexoelectricity has an influence on the static deflection, with respect to the thickness of hybrid flexoelectric plates. This may be due to the fact that, since the flexoelectric effect is size dependent, i.e., the smaller the plate’s thickness, better the result.

3.2. Modal Analysis of Hybrid Flexoelectric GRPC Plate

In this, the investigations were carried out to study the effect of flexoelectricity on the dynamic response of hybrid GRPC plates with respect to different parameters, such as different graphene contents, thickness and aspect ratio of the plate. Similar to Figure 5, Figure 8 shows the variation of the resonant frequency against the thickness of hybrid GRPC plate, considering the different combination of graphene and PZT for mode $(1,1).$ From this Figure 8, it can be observed that the resonant frequency of plate increases as the value of graphene content in the hybrid composite plate increases. From this figure, it is concluded that effect of flexoelectricity decreases as the thickness of the hybrid plate increases, because the flexoelectricity is a size-dependent phenomenon.

Figure 9 illustrates the results of the effects of flexoelectricity on the dynamic response of piezoelectric hybrid GRPC plates for the normalized fundamental frequency of mode $(1,1)$ for the conventional plate, and the hybrid flexoelectric plate is plotted against the plate thickness. From Figure 9, it can be observed that the hybrid flexoelectric plate envisages a larger resonant frequency over the conventional plate without flexoelectricity. Hence, the numerical results for plate resonant frequency can be enhanced by the flexoelectricity effect. Figure 10 demonstrates the results of the effects of flexoelectricity on the dynamic response of hybrid flexoelectric plates for resonant frequency of mode $(1,1)$, and $(2,2)$ is plotted against the plate aspect ratio. In Figure 10, the legend or graph with the blue color indicates the results for hybrid flexoelectric plate, while the graph with the red color indicates the results for conventional plate. We kept the in-plane dimensions of the hybrid plate constant $(\mathrm{a}=\mathrm{b}=1200\mathrm{n}\mathrm{m})$ and varied their modes. From Figure 10, it can be observed that the resonant frequency increases as the mode number increases for the same in-plane dimension of the plate $(\mathrm{a}=\mathrm{b}=1200\mathrm{n}\mathrm{m}).$ The effect of flexoelectricity is noteworthy in the case of the thin plate. It can be also concluded that the influence of flexoelectricity on the resonant frequency of a plate diminishes as its thickness increases, which shows that the flexoelectric effect is size-dependent. It is obvious that the resonant frequency is higher for the hybrid flexoelectric plate over that of the conventional plate when the plate thickness is less than $30\mathrm{n}\mathrm{m}$. The flexoelectricity does not much influence the resonant frequencies of plates having larger thickness $(>40\mathrm{n}\mathrm{m})$, and this is because of the effect of size-dependent flexoelectricity diminishes as the thickness of plate increases, despite the fact that the resonant frequency at mode $(1,1)$ is less, as compared to mode $(2,2).$ The same is true for another set of modes, $(2,2)$ and $(3,3)$, and so on. This is attributed to the fact that mode $(1,1)$ is a fundamental mode of vibration, which is the lowest natural frequency of the particular system. Normally, only the first few modes are vital from the practical application point of view; although, the results for the higher modes are also presented here. Modes are numbered according to the number of half and full waves (crest and trough) in the vibration. Differences between the numerical outcomes of all modes may be due to the bending, contracting and expanding of plates under the different modes while the specific amount of energy is stored in each mode. It means that, if mode number increases, then the crest and trough in the vibration get increased, which results in more energy being stored in the structures. Therefore, this resonant frequency is directly proportional to the mode numbers. Moreover, from Equation (44) it is clearly seen that the resonant frequency due to piezoelectric/flexoelectric effect strongly depends on the mode numbers $\mathrm{m}$ and $\mathrm{n}.$ Figure 9 and Figure 10 also reveal that the resonant frequency is largely depends on the in-plane dimensions of plate; resonant frequency of the plate diminishes as its in-plane dimensions increase. Here, the results of resonant frequencies of hybrid flexoelectric plates are presented to investigate the effect of flexoelectricity and different modes by varying in-plane dimension with respect to plate thickness. Figure 11 demonstrates the effect of flexoelectricity on the resonant frequency of hybrid flexoelectric GRPC plate for different modes, such as $(1,1),(2,2),(3,3),(4,4)\mathrm{a}\mathrm{n}\mathrm{d}(5,5)$, with respect to plate thickness. The results shown in Figure 9 and Figure 10 are significant, thus indicating that the flexoelectricity plays an important role in the dynamics of thin plates and needs to be accounted properly. It can be observed that the resonant frequency decreases as the aspect ratio of plate increases. It is observed from looking at Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11 that, as the thickness of the plate increases, the effect of flexoelectricity starts diminishing.

4. Conclusions

In this novel article, static and modal analyses of hybrid graphene-reinforced piezoelectric composite (GRPC) plates with flexoelectric effect were investigated for the first time in the literature. Finally, we derived the closed-form solutions for hybrid flexoelectric GRPC plate based on Kirchhoff’s plate theory, Navier’s solution and extended linear theory of piezoelectricity. In this, the governing equations and respective boundary conditions were obtained, using Hamilton’s variational principle for achieving static deflection and natural frequency. Moreover, the different parameters, such as aspect ratio, thickness of plate, different loadings (uniformly distributed, uniformly varying, inline and point load), the combination of different volume fraction of graphene and piezoelectric lead zirconate titanate (PZT), were considered. In addition to this, different mode shapes and flexoelectric coefficients were also considered. It is found that the bending rigidity of hybrid flexoelectric plates having a thickness less than 50 nm increases significantly, and such an effect cannot be neglected for predicting the static response of thin structures. Similarly, the dynamics response of hybrid GRPC plates is enhanced due to the flexoelectric effect as the plate thickness reduces. In addition to this, our results clearly indicate that the consideration of adding the proper percentage of graphene and flexoelectricity plays an important role in the reducing static bending deflection (overall highest displacement) and increasing modal analysis in terms of resonant frequency of thin plates and cannot be neglected while modeling 2D nanostructures.

The current study might provoke the investigators for building this novel hybrid GRPC and aid the purpose of validating the experimental approximations. Meanwhile one may use this closed-form solution for prediction of behavior of hybrid composite structures. Some of the improvements outlined in the current paper will guide the path for studying the structural health monitoring and energy harvester in MEMS/NEMS application in the future.