Isogeometric Analysis of Graphene-Reinforced Functionally Gradient Piezoelectric Plates Resting on Winkler Elastic Foundations

In this paper, the refined plate theory (RPT), Hamilton’s principle, and isogeometric analysis (IGA) are applied to investigate the static bending, free vibration and buckling behaviors of functionally graded graphene-platelet-reinforced piezoelectric (FG-GRP) plates resting on a Winkler elastic foundation. The graphene platelets (GPLs) are distributed in polyvinylidene fluoride (PVDF) as a power function along the plate thickness direction to generate functionally gradient materials (FGMs). The modified Halpin–Tsai parallel model predicts the effective Young’s modulus of each graphene-reinforced piezoelectric composite plate layer, and the rule of the mixture can be used to calculate the effective Poisson’s ratio, mass density, and piezoelectric properties. Under different graphene distribution patterns and boundary conditions, the effects of a plate’s geometric dimensions, GPLs’ physical properties, GPLs’ geometric properties and the elastic coefficient of the Winkler elastic foundation on deflections, frequencies and bucking loads of the FG-GRP plates are investigated in depth. The convergence and computational efficiency of the present IGA are confirmed versus other studies. Furthermore, the results illustrate that a small amount of GPL reinforcements can improve the FG-GRP plates’ mechanical properties, i.e., GPLs can improve the system’s vibration and stability characteristics. The more GPL reinforcements spread into the surface layers, the more effective it is at enhancing the system’s stiffness.


Introduction
Functionally gradient materials (FGMs), first developed in 1987, are composite materials with two or more distinct phases. FGMs provide designers with a range of options to distribute strength and stiffness in an appropriate manner [1]. Since FGMs have great properties and match well with many engineering challenges, many researchers have focused on FG structures formed of FGMs. Additionally, their mechanical properties, such as bending, vibration, and buckling behaviors, have been explored and realized more extensively [2][3][4]. Novoselov et al. [5] succeeded in extracting graphene platelets (GPLs) in 2004. Since then, a number of theoretical and experimental studies have focused on understanding the behavior of GPLs. Because GPLs have high strength, stiffness as well as good thermal and electrical properties [6,7], many experimental and numerical results have shown that combining GPL reinforcements with FGMs can improve the mechanical, piezoelectric, and thermoelectric properties of the composite structures [8][9][10]. According to Lee et al. [11], when the graphene volume fraction reached 1.6%, the strength of the composites increased to roughly 80% compared with pure epoxy. As a result, functionally graded GPL-reinforced composites (FG-GPLRCs) have become one of the most popular composites today, with correction factor by introducing higher-order terms in the displacement field approximation. However, HSDT is computationally expensive. Relatively speaking, the refined plate theory (RPT) [40] offers apparent advantages. As a simplification of the third-order shear deformation theory (TSDT), RPT not only preserves the benefits of TSDT, but also eliminates an unknown variable relative to TSDT. On the other hand, RPT must match C 1 -continuity of the displacement fields, which presents challenges to the traditional finite element analysis (FEA) framework [41]. Furthermore, the generally employed numerical methods would result in disagreement between the geometric model and the analytical model when the geometries are complicated, unavoidably introducing a significant number of interactions between computer aided design (CAD) and computer aided engineering (CAE) data. To address such issues, Hughes et al. [42] originally presented the notion of isogeometric analysis (IGA) in 2005, which is a novelty computational approach that smoothly integrates CAD with CAE. The fundamental concept of IGA is to apply spline functions (non-uniform rational B-splines [43,44], T-splines and so on) to express the geometry and approximate the unknown fields. IGA not only avoids the approximation error caused by the interpolation function of polynomials, but it also solves the C 1 -continuity of displacement fields required by RPT. In recent years, researchers have frequently used IGA to analytically calculate the mechanical properties of functionally graded composites. On the basis of IGA and TSDT, Kiani et al. [45,46] investigated the large amplitude free vibration and post-buckling properties of GPL-reinforced laminates in a thermal environment. In the context of FSDT and HSDT, Li et al. [47] utilized IGA to investigate the static linear bending, natural frequency, and buckling behavior of GPL-reinforced FG porous plates. Liu et al. [48] combined IGA with simple first-order shear deformation theory (S-FSDT) to analyze the static bending, free vibration, and dynamic response of FG plates integrated with piezoelectric actuators and sensors. Phung-Van et al. [49] investigated the static bending, free vibration, and dynamic control of composite plates integrated with piezoelectric actuators and sensors based on the combination of HSDT and IGA.
To the best of the authors' knowledge, no literature is available for the analysis of functionally graded graphene-platelet-reinforced piezoelectric structures based on isogeometric analysis. As a result, the purpose of this work is to investigate the static bending, free vibration and buckling behaviors of functionally graded graphene-platelet-reinforced piezoelectric plates embedded on a Winkler elastic foundation using isogeometric analysis and refined plate theory. The object of this paper is to development isogeometric analysis of graphene-reinforced functionally gradient piezoelectric plates resting on Winkler elastic foundations. This paper is organized as follows. In Section 2, the graphene-reinforced functional gradient piezoelectric plate model is established by refined plate theory and Hamilton's principle. Isogeometric analysis formulations are established in Section 3. In Section 4, detailed parameter investigations are discussed to demonstrate the influences of a plate's geometric dimensions, GPLs' geometric properties, GPLs' volume fraction and the elastic coefficient of the Winkler elastic foundation on mechanical behaviors. Finally, some main conclusions are drawn in Section 5.

Material Properties
As illustrated in Figure 1, an FG-GRP plate of length a, width b, and thickness h on Winkler elastic foundation is established. The FG-GRP plate is supposed to be made of a multi-layer polymer composite reinforced by graphene platelets, where the GPL reinforcements are functionally graded along the thickness direction with a power law distribution. Furthermore, the thickness of each independent layer may be described as ∆h = h/N. In this research, three graphene distribution patterns are generated by smoothly varying the graphene volume fractions along the thickness direction, as shown in Figure 2: denoted as the U type, the X type, and the O type. The graphene volume fraction of the top and bottom surfaces is significantly larger than that of the mid-plane in the X type, indicating that the reinforcing effect of GPL in the top and bottom surfaces is more remarkable than that in the mid-plane. On the contrary, the graphene volume fraction of the top and bottom surfaces is significantly lower than that in the mid-plane in the O type, indicating that GPL reinforcements have opposite strengthening effects for both the O type and X type. The graphene volume fraction maintains constant along the thickness direction for the U type, resulting in a uniform distribution on the whole. Both the X and O types are functionally gradient distributions in which the graphene volume fraction increases or decreases linearly and symmetrically from the mid-plane to the top and bottom surfaces of the composite.  In this research, three graphene distribution patterns are generated by smoothly varying the graphene volume fractions along the thickness direction, as shown in Figure 2: denoted as the U type, the X type, and the O type. The graphene volume fraction of the top and bottom surfaces is significantly larger than that of the mid-plane in the X type, indicating that the reinforcing effect of GPL in the top and bottom surfaces is more remarkable than that in the mid-plane. On the contrary, the graphene volume fraction of the top and bottom surfaces is significantly lower than that in the mid-plane in the O type, indicating that GPL reinforcements have opposite strengthening effects for both the O type and X type. The graphene volume fraction maintains constant along the thickness direction for the U type, resulting in a uniform distribution on the whole. Both the X and O types are functionally gradient distributions in which the graphene volume fraction increases or decreases linearly and symmetrically from the mid-plane to the top and bottom surfaces of the composite. In this research, three graphene distribution patterns are generated by smoothly varying the graphene volume fractions along the thickness direction, as shown in Figure 2: denoted as the U type, the X type, and the O type. The graphene volume fraction of the top and bottom surfaces is significantly larger than that of the mid-plane in the X type, indicating that the reinforcing effect of GPL in the top and bottom surfaces is more remarkable than that in the mid-plane. On the contrary, the graphene volume fraction of the top and bottom surfaces is significantly lower than that in the mid-plane in the O type, indicating that GPL reinforcements have opposite strengthening effects for both the O type and X type. The graphene volume fraction maintains constant along the thickness direction for the U type, resulting in a uniform distribution on the whole. Both the X and O types are functionally gradient distributions in which the graphene volume fraction increases or decreases linearly and symmetrically from the mid-plane to the top and bottom surfaces of the composite. The total graphene volume fraction is expressed as V gpl , and that the minimum graphene volume fraction in each layer is expressed as V * for the X and O types, then The volume fraction in the k-th layer of the FG-GRP plate can be written as follows: The U type: The X type: The O type: As demonstrated by Layek et al. [28], GPL reinforcements are more easily dispersed in the matrix material, as the graphene volume fraction is less than 1%. As a result, the FG-GRP plate with a graphene volume fraction of less than 1% is the subject of this research. The effective Young's modulus of an FG-GRP plate is given for the k-th layer as where in which the subscripts G and M symbolize graphene platelet and matrix material, and a gpl , b gpl and h gpl represent the GPL's length, width and thickness, respectively. Furthermore, the effective material properties, such as Poisson's ratio, density, piezoelectric constant, and dielectric constant for the k-th layer may be predicted by

Governing Equations for the Refined Plate Theory
The displacement fields (u, v, w) consisting of four unknown variables, according to the refined plate theory, may be described as where u 0 and v 0 represent the in-plane displacement components in the mid-plane, w b and w s denote the transverse bending and shear displacement components, and f (z) = −4z 3 /3h 2 .
The in-plane strain and transverse shear strains corresponding to the displacement fields are shown as According to the generalized Hooke's law, the k-th layer constitutive relationship is represented as where Q (k) For the FG-GRP plate, the constitutive equation is expressed as in which σ and ε represent the plate's stress and strain vectors, D is the electric displacement vector, c denotes the elastic stiffness matrix, e and g are the piezoelectric and dielectric constant matrices, respectively. It should be mentioned that the electric field vector E is given as According to Equations (9)-(13), the FG-GRP plate's mid-plane internal force N, bending moment M, high-order bending moment P, and transverse shear force R are written as follows: where and The internal force N E and bending moment M E , P E induced by an electric field are be defined as

Weak Form
The Galerkin weak form of the governing equations is calculated using the Hamilton principle [50].
in which u and . u represent the displacement and velocity vectors, respectively, φ denotes electric potential, f s and F p represent surface force and concentrated force, respectively, and q s denotes surface charge.

NURBS Basis Function
The main elements of the B-spline curve are a set of control points P i and a group of knot vectors Ξ formed of non-decreasing parameter values ξ i , i = 1, . . . , n + p. The definition formula of Ξ is as follows: where ξ i , n and p denote the ith knot, the number of basic functions, and the polynomial order, respectively. According to the Cox-de Boor recursive formula [42], the B-spline basis function is defined as It should be noted that the ratio 0/0 in the preceding Equation (20) is given as zero.
The B-spline basis function is mostly made up of specified spline basis function orders and related knot vectors, and its expression is expressed as follows: The B-spline function can be a particular case of the NURBS function when weights of control points are constant. The bi-variate NURBS basis functions are written as The NURBS spline surfaces are defined by supplying knot vectors and control points in two parameter directions, and the derivative of ξ and η in a bivariate NURBS basis function is expressed as Materials 2022, 15, 5727

Approximations of the Displacement Field
By adopting the NURBS basis functions, the approximate mechanical displacement field u(ξ, η) of the FG-GRP plates based on RPT can be described as in which m × n represents the total number of basic functions, and d I = u 0I v 0I w bI w sI T denotes the displacement vector of control point P I . Taking Equation (25) into Equation (9) yields where

Approximations of the Electric Potential
After approximating the potential field, the system can be discretized into several sublayer elements along the thickness direction [51]. The approximate potential field φ i (z) is written as follows: in which φ i represents electric potential of the upper and lower layers related to the ith sublayer, Here, assuming that each individual piezoelectric layer has the same potential at the same height, we may obtain the electric field of each sublayer unit by substituting Equation (28) into Equation (13). where

Isogeometric FG-GRP Plate Formulation
The governing equation of motion of the FG-GRP plate can be found by inserting Equations (12), (26) and (30) into Equation (18): where where where f s denotes the transverse mechanical surface load, q E is charge surface density, and where U is the applied voltage to the upper and lower surfaces. The element stiffness matrix for buckling analysis can be represented as: where The edges of the FG-GRP plate are subjected to a uniaxial load and biaxial loads, as shown in Figure 3a,b. The connection between buckling forces and load parameter index is expressed as follows: where ξ 1 and ξ 2 are parameters that control load conditions.
The edges of the FG-GRP plate are subjected to a uniaxial load and biaxial loads, a shown in Figure 3a,b. The connection between buckling forces and load parameter index is expressed as follows: ...
Eliminating the potential and assembling elements, the governing equations of FG GRP are represented as: The governing equations for static analysis, free vibration and buckling behaviors ar obtained as follows To evaluate the impact of the Winkler elastic foundation, the additional matrix K e f of the elastic foundation is introduced: (51) where N = R P 1 I 4×4 R P 2 I 4×4 . . . R P N I 4×4 , k f is the elastic foundation coefficient. Eliminating the potential and assembling elements, the governing equations of FG-GRP are represented as: The governing equations for static analysis, free vibration and buckling behaviors are obtained as follows

Numerical Investigation
In this section, static analysis, free vibration and buckling behaviors of FG-GRP plates on Winkler elastic foundation are investigated by isogeometric analysis. Unless otherwise specified, an FG-GRP plate with geometric size a = b = 0.1 m, thickness h = 0.01 m is considered. Additionally, the graphene volume fraction is set as 0.5%. The GPLs' piezoelectric characteristics can be defined as e 31,G = αe 31 [30]. Furthermore, the relevant parameters of GPL and PVDF [28] are: The following dimensionless formula is used in this paper unless otherwise specified: Firstly, the convergence of the IGA technique is examined by adopting (p + 1)(q + 1) Gauss quadrature, where p and q indicate the orders of the basic functions along the directions of ξ and η, respectively. Furthermore, the effects of the geometric parameters of the plate, graphene volume fraction, the distribution pattern, and the Winkler elastic coefficient on the deflections, frequencies and bucking loads of FG-GRP plates resting in Winkler elastic foundation are discussed.

Verification
Firstly, the FG-GRP plates are modelled with 7 × 7, 11 × 11 and 19 × 19 control points as shown in Figure 4a-c. Then, the convergence and accuracy of the present solutions at different mesh levels are presented in Table 1 for an FG-GRP plate under CCCC and SSSS boundary conditions. Table 1 displays the dimensionless frequencies and critical buckling loads of the FG-GRP plates. It is observed from Table 1 that the present calculation results are highly consistent with those listed in studies [30,31], and the IGA method has reached the convergence state when the control points are greater than 11 × 11. As a result, the number of control points utilized in the following examples in this paper is 11 × 11. Lastly, the effects of Winkler elastic coefficient k l on the central deflection are shown in Table 2. It is seen that the normalized deflections in this paper match well with the results of Shojaee et al. [52]. Table 3 compares the dimensionless vibration frequency and percentage frequency change of the FG-GRP plate to the calculation results of Mao et al. [30] for different GPL distribution patterns. As shown in Table 3, the percentage change on frequency in terms of the X type is the highest one among three kinds of distribution patterns. Table 4 compares the dimensionless critical buckling loads P * cr = (P cr a 2 )/(E m h 3 ) of FG-GRP plates to the calculation findings of Mao et al. [31]. It can be seen that the calculation results are mostly consistent.       Figure 5a,b show the effect of the number of layers on the percentage deflection change of the FG-GRP plate under different graphene distribution patterns. For the U type, the number of layers has no effect on the percentage deflection change. Moreover, the influence of the layers on the percentage deflection change is extremely obvious for the X type and the O type. For the X type, with a rise in the number of layers N, more GPL reinforcements are spread on both sides of the plate to enhance the stiffness of the system. While for the other one, the O type plays opposite roles on the stiffness and dimensionless central deflections. Simultaneously, with increasing the number of layers N, the percentage deflection changes corresponding to the X and O types gradually converge when N ≥ 16. Therefore, in the following analysis, the number of layers N = 16 is employed for computation to improve calculation efficiency.  Figure 5a,b show the effect of the number of layers on the percentage deflection change of the FG-GRP plate under different graphene distribution patterns. For the U type, the number of layers has no effect on the percentage deflection change. Moreover the influence of the layers on the percentage deflection change is extremely obvious for the X type and the O type. For the X type, with a rise in the number of layers N, more GPL reinforcements are spread on both sides of the plate to enhance the stiffness of the system While for the other one, the O type plays opposite roles on the stiffness and dimensionless central deflections. Simultaneously, with increasing the number of layers N, the percentage deflection changes corresponding to the X and O types gradually converge when 16 N  . Therefore, in the following analysis, the number of layers 16 N  is employed for computation to improve calculation efficiency. The effective material properties for the composite are determined by the length and thickness of the graphene reinforcements but have nothing to do with the width in Halpin

Static Analysis
Tsai's parallel model [28]. Assuming  The effective material properties for the composite are determined by the length and thickness of the graphene reinforcements but have nothing to do with the width in Halpin Tsai's parallel model [28]. Assuming b gpl = 1.5 µm, Figure 6a-d illustrate the percentage change on deflection under a different graphene length-to-thickness ratio a gpl /h gpl . It is clearly discerned from this figure that the percentage deflection change increases with the increase in a gpl /h gpl , indicating that smaller and thinner graphene reinforcements can more effectively enhance the stiffness of the FG-GRP plate; this is due to the fact that the greater the contact area between the graphene platelet and the matrix material, the better the bonding ability. Furthermore, compared to Figure 6a-d, it can be shown that the percentage deflection change for a given pattern under different boundary conditions is almost the same. As a result, the following analysis only takes into account the percentage deflection changes under the SSSS boundary condition.
The impact of the Winkler elastic coefficients k l on the normalized central deflections w r of the FG-GRP plate under varied graphene volume fractions V gpl is depicted in Figure 7a. The dimensionless central deflections of the FG-GRP plate decrease as the Winkler elastic coefficients increase, as seen in this figure. Simultaneously, it is discerned that the normalized central deflections decrease with increasing GPL volume fractions. This is due to the fact that graphene reinforcements can effectively enhance the system stiffness, and the higher the graphene volume fraction (less than 1%), the stronger the impact. Figure 7b displays  Simultaneously, it discerned that the normalized central deflections decrease with increasing GPL volum fractions. This is due to the fact that graphene reinforcements can effectively enhance t system stiffness, and the higher the graphene volume fraction (less than 1%), the strong the impact. Figure

Vibration Analysis
In what follows, the influence of detailed parameters on free vibration analysis for a FG-GRP plate with Figure 8 illustrates the percentage chang on frequency under a different graphene length-to-thickness ratio / gpl gpl a h , graphen

Vibration Analysis
In what follows, the influence of detailed parameters on free vibration analysis for an FG-GRP plate with b gpl = 1.5 µm is analyzed. Figure 8 illustrates the percentage change on frequency under a different graphene length-to-thickness ratio a gpl /h gpl , graphene volume fraction and graphene distribution patterns (the U type, the X type and the O type). It can be clearly seen from this figure that the percentage frequency change increases with the increase in a gpl /h gpl , indicating that an increasing graphene length-to-thickness ratio can enhance the stiffness of the FG-GRP plate. Furthermore, this figure demonstrates that, regardless of graphene distribution patterns, the percentage change on frequency increases with an increasing graphene volume fraction.

Vibration Analysis
In what follows, the influence of detailed parameters on free vibration analysis for an FG-GRP plate with 1.5μm gpl b  is analyzed. Figure 8 illustrates the percentage change on frequency under a different graphene length-to-thickness ratio / gpl gpl a h , graphene volume fraction and graphene distribution patterns (the U type, the X type and the O type). It can be clearly seen from this figure that the percentage frequency change increases with the increase in / gpl gpl a h , indicating that an increasing graphene length-tothickness ratio can enhance the stiffness of the FG-GRP plate. Furthermore, this figure demonstrates that, regardless of graphene distribution patterns, the percentage change on frequency increases with an increasing graphene volume fraction. As shown in Figure 9, the effect of plate geometric dimension on the free vibration characteristics is discussed, including aspect ratio / a b and length-to-thickness ratio / a h . The FG-GRP plate with =0.1m a is applied in the following parametric studies. Because of the uniformity, the geometry dimension of the plate has no effect on the vibration characteristics corresponding to the U type. Since the surface area of the square plate in this paper is larger than that of the rectangular plate, the GPL reinforcements of the X type are more distributed at the top and bottom surfaces of the square plate when the same GPL reinforcements are used. As a result, the square plate presents more remarkable reinforcement effects than the rectangle plate in terms of the X-type plate. As shown in Figure 9, the effect of plate geometric dimension on the free vibration characteristics is discussed, including aspect ratio a/b and length-to-thickness ratio a/h. The FG-GRP plate with a = 0.1 m is applied in the following parametric studies. Because of the uniformity, the geometry dimension of the plate has no effect on the vibration characteristics corresponding to the U type. Since the surface area of the square plate in this paper is larger than that of the rectangular plate, the GPL reinforcements of the X type are more distributed at the top and bottom surfaces of the square plate when the same GPL reinforcements are used. As a result, the square plate presents more remarkable reinforcement effects than the rectangle plate in terms of the X-type plate. The effect of the Winkler elastic coefficient on the dimensionless vibration frequencies of the FG-GRP plate is shown in Figure 10a,b. Figure 10a presents the influences of the Winkler elastic coefficient l k on the dimensionless vibration frequency  Figure 9. Variation of the percentage frequency change along the side-to-width ratio and side-tothickness ratio for an SSSS FG-GRP plate.
The effect of the Winkler elastic coefficient on the dimensionless vibration frequencies of the FG-GRP plate is shown in Figure 10a,b. Figure 10a presents the influences of the Winkler elastic coefficient k l on the dimensionless vibration frequency ω l under different graphene volume fractions V gpl . It is observed that increasing Winkler elastic coefficients leads to increased dimensionless vibration frequencies. Figure 10b examines the dimensionless vibration frequency of the FG-GRP plate versus elastic coefficients under different graphene distribution patterns. It is clearly demonstrated that the increasing frequencies increase with increasing elastic coefficients under different graphene distribution patterns, which shows that the increase in Winkler elastic coefficient can effectively enhance the stiffness of the FG-GRP plate. Figure 9. Variation of the percentage frequency change along the side-to-width ratio and side-t thickness ratio for an SSSS FG-GRP plate.
The effect of the Winkler elastic coefficient on the dimensionless vibratio frequencies of the FG-GRP plate is shown in Figure 10a

Buckling Analysis
To explore the impact of detailed parameters on the buckling behaviors, the uniaxi or biaxial in-plane forces must be applied to the piezoelectric plates.

Buckling Analysis
To explore the impact of detailed parameters on the buckling behaviors, the uniaxial or biaxial in-plane forces must be applied to the piezoelectric plates. Figure 11a-d show the effect of graphene volume fraction V gpl on the normalized critical buckling load P * cr of the FG-GRP plates for different boundary conditions and graphene distribution patterns. The normalized critical buckling load P * cr subjected to uniaxial load is approximately twice as large as that subjected to biaxial load, owing to the fact that the biaxial load in x and y directions is loaded from four directions, limiting the load it can carry. Next, the normalized critical buckling load P * cr is the largest one under CCCC, followed by CCSS, CSSS, and SSSS boundary conditions. For the uniaxial and biaxial loads, the normalized critical buckling load P * cr increases with the increasing graphene volume fraction V gpl . When the graphene volume fraction is the same, the X-type plate has the highest critical buckling load, while the O-type plate has the lowest one. Hence, the following study exclusively analyzes the buckling behavior of the X-type plate under the SSSS boundary condition for the sake of brevity.
The impact of the length-to-width ratio a/b and the graphene length-to-thickness ratio a gpl /h gpl on the normalized critical buckling load of the X-type plate under uniaxial load is examined in Figure 12. When a gpl /h gpl increases from 0 to 2, the normalized critical buckling load P * cr increases significantly, and when a gpl /h gpl continues to rise, P * cr rises slowly, indicating that the longer and thinner GPL reinforcements can reinforce the buckling resistance. Furthermore, under certain a gpl /h gpl , the normalized critical buckling load P * cr increases with length-to-width ratio a/b increasing. Under two loading conditions, Figure 13 shows the effect of the Winkler elastic foundation's elastic coefficient k l on the normalized critical buckling load P * cr of the X type. The normalized critical buckling load increases as the elastic coefficient grows, showing that increasing the elastic coefficient enhances system stiffness.
the load it can carry. Next, the normalized critical buckling load P cr is the largest one under CCCC, followed by CCSS, CSSS, and SSSS boundary conditions. For the uniaxia and biaxial loads, the normalized critical buckling load P cr  increases with the increasing graphene volume fraction V gpl . When the graphene volume fraction is the same, the X type plate has the highest critical buckling load, while the O-type plate has the lowest one Hence, the following study exclusively analyzes the buckling behavior of the X-type plat under the SSSS boundary condition for the sake of brevity. The impact of the length-to-width ratio / a b and the graphene length-to-thicknes ratio / gpl gpl a h on the normalized critical buckling load of the X-type plate under uniaxia load is examined in Figure 12 Figure 13 shows the effect of the Winkler elastic foundation's elastic coefficien l k on the normalized critical buckling load P cr  of the X type. The normalized critica buckling load increases as the elastic coefficient grows, showing that increasing the elasti coefficient enhances system stiffness.

Conclusions
Isogeometric analysis integrated with the refined plate theory is utilized in this pape to analyze the static bending, free vibration and buckling behaviors of functionally grade graphene-platelet-reinforced piezoelectric plates on Winkler elastic foundations. Whe Hamilton's principle is combined with isogeometric analysis, the isogeometric finite ele ment governing equations of the functionally graded graphene-platelet-reinforced piezo electric plate are obtained on the basis of the refined plate theory. Comparisons with othe studies validate the accuracy and convergence of the present isogeometric analysis. Fur thermore, the influences of a plate's geometrical parameters, the elastic coefficient, an GPLs' physical and geometric properties on the percentage change of central deflection frequencies and critical buckling loads of a functionally graded graphene-platelet-rein forced piezoelectric plate are thoroughly investigated under different boundary cond tions and graphene distribution patterns.
The results of this paper indicate that adding a small amount of GPL not only de creases the central deflections, but also increases the vibration frequencies and critica buckling loads of functionally graded graphene-platelet-reinforced piezoelectric plate Additionally, the X type represents the optimum distribution form, i.e., the higher th graphene volume fractions spread into the surface layer, the more effective it is to enhanc the stiffness of the system. Additionally, increasing the graphene length-to-thickness rati can significantly enhance system stiffness, that is, the longer and thinner the GPLs are, th more obvious the reinforcement is. In terms of a plate's geometrical parameters, the effec of a functionally graded graphene-platelet-reinforced piezoelectric plate's aspect ratio o

Conclusions
Isogeometric analysis integrated with the refined plate theory is utilized in this paper to analyze the static bending, free vibration and buckling behaviors of functionally graded graphene-platelet-reinforced piezoelectric plates on Winkler elastic foundations. When Hamilton's principle is combined with isogeometric analysis, the isogeometric finite element governing equations of the functionally graded graphene-platelet-reinforced piezoelectric plate are obtained on the basis of the refined plate theory. Comparisons with other studies validate the accuracy and convergence of the present isogeometric analysis. Furthermore, the influences of a plate's geometrical parameters, the elastic coefficient, and GPLs' physical and geometric properties on the percentage change of central deflections, frequencies and critical buckling loads of a functionally graded graphene-platelet-reinforced piezoelectric plate are thoroughly investigated under different boundary conditions and graphene distribution patterns.
The results of this paper indicate that adding a small amount of GPL not only decreases the central deflections, but also increases the vibration frequencies and critical buckling loads of functionally graded graphene-platelet-reinforced piezoelectric plates. Additionally, the X type represents the optimum distribution form, i.e., the higher the graphene volume fractions spread into the surface layer, the more effective it is to enhance the stiffness of the system. Additionally, increasing the graphene length-to-thickness ratio can significantly enhance system stiffness, that is, the longer and thinner the GPLs are, the more obvious the reinforcement is. In terms of a plate's geometrical parameters, the effect of a functionally graded graphene-platelet-reinforced piezoelectric plate's aspect ratio on central deflections, vibration frequencies and critical buckling loads should be handled differently depending on the graphene distribution pattern. For the X type, increasing the aspect ratio or plate surface area can improve system stiffness. In terms of the impacts of boundary conditions on the system, the stiffness of the functionally graded graphene-platelet-reinforced piezoelectric plate under the CCCC boundary condition is the highest one, followed by CCSS, CSSS, and SSSS. Furthermore, the dimensionless critical buckling load subjected to the uniaxial force is approximately twice as large as that subjected to the biaxial force. Effectively increasing the elastic coefficient of the Winkler elastic foundation can obtain a remarkable reinforcement effect on the plate stiffness.