Bending and Elastic Vibration of a Novel Functionally Graded Polymer Nanocomposite Beam Reinforced by Graphene Nanoplatelets

A novel functionally graded (FG) polymer-based nanocomposite reinforced by graphene nanoplatelets is proposed based on a new distribution law, which is constructed by the error function and contains a gradient index. The variation of the gradient index can result in a continuous variation of the weight fraction of graphene nanoplatelets (GPLs), which forms a sandwich structure with graded mechanical properties. The modified Halpin–Tsai micromechanics model is used to evaluate the effective Young’s modulus of the novel functionally graded graphene nanoplatelets reinforced composites (FG-GPLRCs). The bending and elastic vibration behaviors of the novel nanocomposite beams are investigated. An improved third order shear deformation theory (TSDT), which is proven to have a higher accuracy, is implemented to derive the governing equations related to the bending and vibrations. The Chebyshev–Ritz method is applied to describe various boundary conditions of the beams. The bending displacement, stress state, and vibration frequency of the proposed FG polymer-based nanocomposite beams under uniformly distributed loads are provided in detail. The numerical results show that the proposed distributions of GPL nanofillers can lead to a more effective pattern of improving the mechanical properties of GPL-reinforced composites than the common ones.

The dispersion of nanofillers in the polymer matrix based on a reasonable graded distribution can effectively make use of reinforcements. To address the effects of nanofiller distributions on the mechanical behaviors of FG polymer-based nanocomposites, different types of distributions, such as the uniform distribution (UD), FG-V shape, FG-O shape and FG-X shape, were introduced and employed in many reports, e.g., references [21][22][23][24][25][26]. In addition, the distribution laws in forms of general polynomials have also been implemented [27][28][29]. All existing functions to describe the nanofiller distribution law are not adjustable since no adjustable parameter is included. Until now, the distribution law having adjustable parameters, which can lead to continuously graded mechanical properties for the nanofiller-reinforced polymer nanocomposites, has not been reported. It is believed that an adjustable distribution has great potential for introducing a novel type of nanocomposite, and also can be used to optimize the mechanical performances of the nanofiller-reinforced structures. The first contribution of the current work is to propose a new distribution law with an adjustable parameter as the gradient index.
In recent years, the graphene nanoplatelet (GPL) has been considered as an ideal reinforcement because of its excellent material properties, including high stiffness and strength. A polymer with a low concentration of GPLs exhibits distinct improvement in mechanical properties. Yang and his coauthors devoted great efforts to the bending, buckling, and vibration behaviors of FG-GPLRC structures [30][31][32]. Shen et al. [33,34] discussed various results of the nonlinear bending, vibration, and buckling behaviors of FG-GPLRC plates, panels, and shells. Mao and Zhang [35] examined the linear and nonlinear vibration responses of graphene-reinforced piezoelectric composite plates under external voltage excitation. Wang et al. [36] studied the vibration and bending behaviors of functionally graded nanocomposite doubly curved shallow shells reinforced by GPLs. Gholami et al. [37] investigated the nonlinear, harmonically excited vibration of a third-order shear deformable FG-GPLRC rectangular plate. Although the studies concerning on the mechanical behaviors of FG-GPLRC structures are massive, few studies focus on improving the distribution of the GPLs. Furthermore, the bending and vibration performances are the most basic and critical properties for structural components. Lezgy-Nazargah and Salahshuran [38] proposed a novel mixed-field theory with relatively low number of unknown variables for investigating bending and vibration analysis of multi-layered composite plates. Wattanasakulpong and Ungbhakorn [39] studied the bending, buckling and vibration behaviors of carbon-nanotube-reinforced composite (CNTRC) beams. Ghannadpour et al. [40] investigated the bending, buckling and vibration analyses of nonlocal Euler beams. Zhao et al. [41] analyzed the bending and vibrations of functionally graded trapezoidal nanocomposite plates reinforced with GPLs. The second contribution of our work is to introduce a novel FG-GPLRC beam using the proposed distribution law and investigate its bending and vibration behaviors.
In present work, a novel functionally graded graphene nanoplatelet reinforced composite (FG-GPLRC) is proposed by developing a new distribution law. The variation of the gradient index results in a continuous variation of the weight fraction of graphene nanoplatelets (GPLs) and forms graded mechanical properties within composite according to the thickness direction. The modified Halpin-Tsai micromechanics model is employed to evaluate the effective Young's modulus. The bending and elastic vibration behaviors of the novel FG-GPLRC beams are investigated. An improved third-order shear deformation theory (TSDT), which is proved to have a higher accuracy, is implemented to derive the governing equations. Chebyshev-Ritz method is applied to describe various boundary conditions of the beam. Parametric studies are performed and present the bending displacement, stress state, and vibration frequency of the new FG-GPLRC beam under the uniformly distributed loads. The effects of gradient index of the proposed GPL distributions on the bending and vibration performance are addressed.

Evaluation of Effective Mechanical Properties
It is well known that the Eshelby-Mori-Tanaka theory and the Halpin-Tsai micromechanics model are the most common methods to determine the effective elastic moduli of nanocomposite materials. Much literature suggest that the Eshelby-Mori-Tanaka theory is able to yield a more accurate prediction of the nanocomposite mechanical properties [42][43][44][45][46], and this conclusion was supported by experiments for CNT-reinforced nanocomposite.
The Halpin-Tsai model is based on the generalized self-consistent method. Although this model was initially developed for traditional fiber composites, successful efforts were undertaken to apply this model to nanofillers-reinforced nanocomposites by introducing orientation and shape factors that account for the geometry of the filler phase. The Halpin-Tsai model was used to calculate the effective Young's modulus of a polymer nanocomposite reinforced by graphene nanoplatelets (GPLs), and the comparisons between the theoretical predictions and experimental results were also performed [47,48]. Due to the simple form in mathematics, the Halpin-Tsai micromechanics model has been widely employed to estimate the effective Young's modulus of functionally graded graphene nanoplatelets reinforced composites [27][28][29][30][31][32][33][34][35][36][37]. The main objective of the current work is to propose an adjustable distribution law to find a more effective way to use the GPL reinforcements. It is believed that the effective material properties resulted from the Halpin-Tsai model can provide the results with some reference values. Thus, the Halpin-Tsai micromechanics model is selected.
A beam composed of FG-GPLRCs is considered herein. The length, width, and thickness of the beam are L, b, and h, respectively. The origin of the coordinate system is fixed at the center at the left end of the beam. The matrix of the nanocomposite is a polymer, and the GPL nanofillers are uniformly or non-uniformly dispersed across the thickness direction of the beam. A modified Halpin-Tsai model is used to calculate the effective Young's modulus of the GPL/polymer composite. We assume that GPLs are effective rectangular solid fillers dispersed in a polymer matrix, and the effective Young's modulus E C of the GPL/polymer composite is approximated by Voigt-Reuss model as follows [18]: where longitudinal modulus E L and transverse modulus E T are determined by the Halpin-Tsai model as follows: Substituting Equations (2) and (3) into Equation (1), we obtain the following expression: where where E M and E GPL are Young's moduli of the polymer matrix and GPLs, respectively; V GPL is the volume fraction of the graphene nanoplatelets; ξ L and ξ W characterize the geometry and size of the GPL nanofillers, respectively, and are defined as follows: where l GPL , w GPL , and h GPL are the average length, width, and thickness of GPLs, respectively. Mass density ρ C and Poisson's ratio v C of the GPL/polymer nanocomposite are calculated by the rule of mixtures as follows: where V M is the volume fraction of the polymer matrix; subscripts "GPL", "M", and "C" denote the GPLs, polymer matrix, and GPL/polymer nanocomposite, respectively. The volume fraction of GPLs is: where g GPL is the total weight fraction of the GPLs in the nanocomposite.

A New GPL Distribution
Previous works proved that different reinforcement distributions significantly affect the mechanical properties of the composites. FG-X, FG-O and FG-V are the most common types of distributions for the GPLs or CNTs, where the total weight fractions or volume fractions of the GPLs or CNTs are constant regardless how the reinforcements dispersed in the polymer matrix.
The motivation of the present work is to find a GPL distribution with an adjustable parameter that controls the weight fraction of GPLs in a preferred direction and maintain a constant total weight fraction. This will be helpful to address the impacts of the GPL distributions on the structural responses and provide a new distribution for designing polymer nanocomposites and optimizing the mechanical properties of FG-GPLRC composites in different applications.
The new distribution function is constructed by using an error function: where Erf(r) = 2 √ π r 0 e −t 2 dt is the error function; r is the gradient index and takes only positive values (r > 0); g 0 GPL is the total weight fraction of the GPLs. Although the value of r cannot be zero, when the selected r approaches 0, the value of Equation (8) approaches 1 as its limit, which represents the uniform GPL distribution pattern. From Equation (8), the following features of the proposed GPL distribution are obtained: 1.
The total weight fraction of GPLs g 0 GPL remains constant with the r variations.

2.
The GPLs symmetrically disperse in the matrix about the mid-plane of the beam.

3.
When r increases, many GPLs are increasingly dispersed to the upper and lower surfaces of the beam.

4.
This GPL graded distribution is continuous along the thickness and makes the mechanical properties of the FG-GPLRCs continuously vary. To clearly illustrate the mentioned characteristics for the proposed distribution law, Figure 1 presents the variation of GPL weight fraction with respect to r. The total weight fraction of the GPLs in Figure 1 is assumed to be 0.5%. Nanomaterials 2019, 9, x FOR PEER REVIEW 5 of 21 4. This GPL graded distribution is continuous along the thickness and makes the mechanical properties of the FG-GPLRCs continuously vary.
To clearly illustrate the mentioned characteristics for the proposed distribution law, Figure 1 presents the variation of GPL weight fraction with respect to r. The total weight fraction of the GPLs in Figure 1 is assumed to be 0.5%.  In fact, a functionally graded GPL-reinforced composite (FG-GPLRC) structure is ideal in combining the advantages of both FGMs and GPLs. However, the fabrication of such functionally graded structures with a continuous and smooth variation of GPLs across the thickness is extremely difficult due to the constraint of manufacture technology. For overcoming this problem, the functionally graded GPL reinforced multilayer structures are introduced, because a functionally graded GPL reinforced multilayer nanocomposite structure in which each individual layer is made from a mixture of uniformly distributed GPL reinforcements and polymer matrix with GPL concentration incrementally varying layer by layer is much easier to fabricate. The published results have provided authentic evidence that such a multilayer structure is an excellent approximation to the ideal functionally graded structure with a continuous and smooth variation of GPLs across the thickness direction when the total number of layers is sufficiently large [18,31,32,[49][50][51].

Theory and Formulations
In the present work, an accurate and efficient modeling is established based on an improved third-order shear deformation theory, and the governing equations associated with the static and   4. This GPL graded distribution is continuous along the thickness and makes the mechanical properties of the FG-GPLRCs continuously vary.
To clearly illustrate the mentioned characteristics for the proposed distribution law, Figure 1 presents the variation of GPL weight fraction with respect to r. The total weight fraction of the GPLs in Figure 1 is assumed to be 0.5%.  In fact, a functionally graded GPL-reinforced composite (FG-GPLRC) structure is ideal in combining the advantages of both FGMs and GPLs. However, the fabrication of such functionally graded structures with a continuous and smooth variation of GPLs across the thickness is extremely difficult due to the constraint of manufacture technology. For overcoming this problem, the functionally graded GPL reinforced multilayer structures are introduced, because a functionally graded GPL reinforced multilayer nanocomposite structure in which each individual layer is made from a mixture of uniformly distributed GPL reinforcements and polymer matrix with GPL concentration incrementally varying layer by layer is much easier to fabricate. The published results have provided authentic evidence that such a multilayer structure is an excellent approximation to the ideal functionally graded structure with a continuous and smooth variation of GPLs across the thickness direction when the total number of layers is sufficiently large [18,31,32,[49][50][51].

Theory and Formulations
In the present work, an accurate and efficient modeling is established based on an improved third-order shear deformation theory, and the governing equations associated with the static and In fact, a functionally graded GPL-reinforced composite (FG-GPLRC) structure is ideal in combining the advantages of both FGMs and GPLs. However, the fabrication of such functionally graded structures with a continuous and smooth variation of GPLs across the thickness is extremely difficult due to the constraint of manufacture technology. For overcoming this problem, the functionally graded GPL reinforced multilayer structures are introduced, because a functionally graded GPL reinforced multilayer nanocomposite structure in which each individual layer is made from a mixture of uniformly distributed GPL reinforcements and polymer matrix with GPL concentration incrementally varying layer by layer is much easier to fabricate. The published results have provided authentic evidence that such a multilayer structure is an excellent approximation to the ideal functionally graded structure with a continuous and smooth variation of GPLs across the thickness direction when the total number of layers is sufficiently large [18,31,32,[49][50][51].

Theory and Formulations
In the present work, an accurate and efficient modeling is established based on an improved third-order shear deformation theory, and the governing equations associated with the static and Nanomaterials 2019, 9, 1690 6 of 21 vibration problems of the beams with various boundary conditions are derived using the minimum energy methodology and Chebyshev-Ritz method.

An Improved Third Order Shear Deformation Theory
The improved third-order shear deformation theory (TSDT) was originally proposed by Shi [52] based on rigorous kinematic of displacements. The published literature has proven that the results obtained by the improved TSDT are more reliable and accurate than many other theories because the kinematics of displacements is derived from elasticity theory instead of the hypothesis of displacement like other existing approaches.
The displacement fields of the improved third-order shear deformation theory (TSDT) are expressed as follows [52,53]: where u 0 and w 0 define the generalized displacements at the mid-plane of the beam in the x and z directions, and ϕ x is the rotation of the beam. From the above displacement fields, the small normal strain ε xx and transverse shear strain γ xz can be written as Then, the stresses are obtained under the assumption of Hook's law: where Q 11 (z) and Q 55 (z) are the elastic constants that continuously vary along the beam thickness and expressed as The strain energy U of the beam is Substituting Equations (9)-(12) into Equation (13), we obtain the strain energy expression as a function of the material stiffness and strain components as The kinetic energy of the composite beam is where ρ(z) is the mass density of the beam, which varies in the thickness direction. Substituting Equation (9) into Equation (16), we can rewrite the kinetic energy as where The work done by a uniform distributed load q is: The total energy function (Π) of the FG-GPLRC beam for the bending problem can be written as: The total energy functional (Π) of FG-GPLRC beams for the free-vibration analysis is expressed as:

Chebyshev-Ritz Method
Ritz method is known as an effective tool to analyze the structural behavior. Since the adoption of trial functions only depends on the essential type of boundary condition [54], various functions may be selected as the admissible functions. In the present work, each displacement function can be written in the form of triplicate series of Chebyshev polynomials multiplied by a boundary function, which ensures that the displacement component satisfies the essential geometric boundary condition of the beam, i.e., Nanomaterials 2019, 9, 1690 8 of 21 where N is the truncated number of Chebyshev polynomials. B Ξ (x) (Ξ = u, w and φ) are the boundary functions; U n , W n , V n denote the unknown coefficients corresponding to time and are expressed as where ω n is the vibration frequency. P i (x) in Equation (21) is the i-th Chebyshev polynomial of the first kind, which is commonly known as "the most optimal expansion" [55,56], and defined in the interval of [−1, 1] as: The recursive relationship is There are two distinct advantages of selecting Chebyshev polynomial series as the admissible functions for each displacement component [57,58]: (1) P i (x) is a set of complete and orthogonal series in the interval of [−1, 1] and has more rapid convergence and better numerical stability in computation than other polynomials; (2) P i (x) can be expressed in a simple and unified form of cosine function as shown in Equation (22), which reduces the coding efforts.
The boundary functions B Ξ (x) (Ξ = u, w and φ) that correspond to u, w and φ are provided in the following uniform formula: where L Θ and R Θ are indices from the following essential geometric boundary conditions: (1) Hinged-Hinged (H-H) Table 1 lists the values of the indices for various boundary conditions.

Boundary Conditions
Substituting Equation (21) into the total energy functional (Π) for the bending in Equation (19) and free vibration in Equation (20) and taking the derivative with respect to the unknown coefficients in the procedure of finding minimization requires The aforementioned procedure produces a system of simultaneous equations with an equal number of unknown coefficients (U n , W n , V n ).
The equation system of solving the static bending problem of the FG-GPLRC beam under a distributed load q can be given as: where [K] is the structural stiffness matrix; F is a column vector associated with the external load from Equation (17); and ∆ is the solution column vector for the static problem.
The generalized eigenvalue problem for free vibration is expressed as follows: where [K] is the structural stiffness matrix, [M] is the mass matrix, and ω is the natural frequency. Vector ∆ is the eigenvector from the displacement functions, which represents the modal shapes of the structures. The dimensions of the aforementioned matrices are 3 N × 3 N. The dimension of ∆ is 3 N × 1.

Convergence and Validation Studies
The validation and accuracy of the current model developed based on the improved third-order shear deformation theory and Chebyshev-Ritz method are verified by comparing the transverse deflections, stresses and vibration frequencies of the FG and FG sandwich beams reported in the existing literature.
A hinged-hinged (H-H) sandwich beam with power-law type FG face sheets and a homogeneous hardcore, which is subjected to a uniform distributed load, is selected for the first example. The transverse deflection (w) at mid-span, normal stresses (σ xx ) at the top surface of the mid-span, transverse shear stress (σ xz ) at mid-plane at the left end, and natural frequencies of the sandwich beam are calculated and compared with the results from [59] in the following dimensionless form: Table 2 tabulates the dimensionless results, which show excellent consistency between the present results and the published results based on a high-order shear deformation theory. Meanwhile, 6 is a reasonable number of Chebyshev polynomial terms, which is sufficiently large to obtain accurate results. In addition, Tables 3 and 4 list the dimensionless mid-span displacements and fundamental frequencies of the sandwich beams with the layer thickness ratio of 1:2:1. The H-H and C-C boundary conditions are used, and various power-law indices (p) are considered. The results are compared with the available ones obtained by the first-order shear deformation theory (FSDT), third-order shear deformation theory (TSDT) and quasi-3D beam theory from the literature [60,61]. The present results well match the published results. It is important to note that the numerical results from the present model are closer to the results from the quasi-3D theory than the other beam theories, which indicates that the present model is accurate and efficient.

Results and Discussion
To facilitate the presentation, the following dimensionless parameters are introduced as follows: Dimensionless normal stress : σ xx = σ xx L 2 , z h qL Dimensionless shear stress : σ xz = σ xz (0, z) h qL Dimensionless frequencies : where ρ M and E M are the density and elastic modulus of the polymer matrix, respectively.

Bending Deflection
In this section, the effects of the gradient indices and GPL weight fractions on the bending behaviors of the new FG-GPLRC beams are investigated. Tables 5-7 present the numerical results of the maximum dimensionless displacement w max of the FG-GPLRC beams with various boundary conditions.  Figure 3 shows the relationship between maximum dimensionless displacement w max and GPL weight fractions. There is a nonlinear variation between the values of w max with respect to the GPL weight fractions. Regardless of the gradient index, w max decreases when the total weight fraction of the GPLs g 0 GPL increases. The curves in Figure 3 clearly indicate that the decrease in w max increases when the GPL content is relatively low, while the decrease tends to be slight when the GPL content increases to a relatively high level, which suggests that the bending response of the beam is more sensitive to the variations of the GPL contents when g 0 GPL is relatively low.
the GPL weight fractions. Regardless of the gradient index, max w decreases when the total weight fraction of the GPLs 0 GPL g increases. The curves in Figure 3 clearly indicate that the decrease in max w increases when the GPL content is relatively low, while the decrease tends to be slight when the GPL content increases to a relatively high level, which suggests that the bending response of the beam is more sensitive to the variations of the GPL contents when 0 GPL g is relatively low.  Figure 4 indicate that parameter r has remarkable effects on the bending performances of the FG-GPLRC beam because the variation of r changes the variation of GPLs in the thickness direction of the beam when the weight fraction is fixed and affects the beam stiffness. For the thick beam (L/h = 5), as shown in Figure 4a, the regulation between r and max w becomes complex and interesting. When the total weight fraction 0 GPL g is very small, e.g., 0.25%, the increase in r will result in a decrease of max w implies that dispersing more reinforcements near the surfaces can improve the bending resistance of the beam. However, when the total weight fraction 0 GPL g is relatively high, e.g., 1.25%, when gradient index r increases, max w first reduces Another concerning point is the effect of the gradient index on the bending displacements of the beam. Figure 4 presents the relationship between the maximum dimensionless displacements w max and the values of gradient index r. L/h = 5 and L/h = 20 are considered. The curves in Figure 4 indicate that parameter r has remarkable effects on the bending performances of the FG-GPLRC beam because the variation of r changes the variation of GPLs in the thickness direction of the beam when the weight fraction is fixed and affects the beam stiffness.  Figure 4b, the maximum dimensionless displacements max w decrease when gradient index r increases regardless of the geometric parameters and boundary conditions. As mentioned, when the total GPL weight fractions remain constant, the increase in r implies that more GPLs are dispersed in the top and bottom surfaces of the beam. Therefore, Figure 4b indicates that dispersing more GPL nanofillers near the top and bottom surfaces of the beam is most effective in enhancing the beam stiffness. The numerical results in Tables 5-7 also support this conclusion.  In addition, for the thin beam, Figure 4b demonstrates that there is a nonlinear descent in the maximum dimensionless displacements max w versus the increase in gradient index r. For comparison, the maximum dimensionless displacement max w of the polymer composite beam of dispersing the GPLs in a FG-X type and UD type, are also displayed in Figure 4b. Figure 4b shows an intersection point of the two curves, which implies that the proposed FG-GPLRC beam and conventional FG-X composite beams have identical bending deflections. The corresponding value of r for this point will be of significance because when r is beyond this point, the new FG-GPLRC beams can provide better deformed-resistant capability than the conventional FG-X ones. In addition, we find that the value in Figure 4b is equal to 1.6 (r = 1.6), and this r at the intersection point is constant regardless of the boundary conditions, geometric parameters and total GPL weight fractions. This result can be mathematically explained that when r = 1.6, the values of the function from Equation (8)  ) at a given z-position, i.e., the GPL distribution along the thickness direction is similar to the GPL distribution from the FG-X For the thick beam (L/h = 5), as shown in Figure 4a, the regulation between r and w max becomes complex and interesting. When the total weight fraction g 0 GPL is very small, e.g., 0.25%, the increase in r will result in a decrease of w max implies that dispersing more reinforcements near the surfaces can improve the bending resistance of the beam. However, when the total weight fraction g 0 GPL is relatively high, e.g., 1.25%, when gradient index r increases, w max first reduces and subsequently increases, which suggests that there is a gradient index value that can minimize w max . This r value can produce the best GPL distribution that provides the greatest improvement to the FG-GPLRC beam.
For the thin beam (L/h = 20), as shown in Figure 4b, the maximum dimensionless displacements w max decrease when gradient index r increases regardless of the geometric parameters and boundary conditions. As mentioned, when the total GPL weight fractions remain constant, the increase in r implies that more GPLs are dispersed in the top and bottom surfaces of the beam. Therefore, Figure 4b indicates that dispersing more GPL nanofillers near the top and bottom surfaces of the beam is most effective in enhancing the beam stiffness. The numerical results in Tables 5-7 also support this conclusion.
In addition, for the thin beam, Figure 4b demonstrates that there is a nonlinear descent in the maximum dimensionless displacements w max versus the increase in gradient index r. For comparison, the maximum dimensionless displacement w max of the polymer composite beam of dispersing the GPLs in a FG-X type and UD type, are also displayed in Figure 4b. Figure 4b shows an intersection point of the two curves, which implies that the proposed FG-GPLRC beam and conventional FG-X composite beams have identical bending deflections. The corresponding value of r for this point will be of significance because when r is beyond this point, the new FG-GPLRC beams can provide better deformed-resistant capability than the conventional FG-X ones. In addition, we find that the value in Figure 4b is equal to 1.6 (r = 1.6), and this r at the intersection point is constant regardless of the boundary conditions, geometric parameters and total GPL weight fractions. This result can be mathematically explained that when r = 1.6, the values of the function from Equation (8) are closer to those from the FG-X function (g(z) = 4 z h ) at a given z-position, i.e., the GPL distribution along the thickness direction is similar to the GPL distribution from the FG-X one when r = 1.6.

Stress State
A hinged-hinged FG-GPLRC beam with a length-to-height ratio of 5 is used to evaluate the normal and shear stresses in the thickness direction, which reveals the effects of the gradient index and GPL weight fraction on the stress distributions of the beam.
The effect of the gradient index on the normal stress distributions is displayed in Figure 5. Figure 5 shows that the normal stress linearly varies with the thickness coordinates when r approaches 0 because in this situation, the GPLs are uniformly dispersed in the polymer matrix, and the mechanical properties of the FG-GPLRCs have no change in the thickness direction. When r increases, nonlinear change occurs and become more remarkable because a greater r indicates that more GPLs are concentrated on the top and bottom surfaces, which makes a significantly uneven GPL distribution and produces the graded mechanical properties of the FG-GPLRCs. In addition, the increasingly concentrated dispersion of GPLs to the surfaces increases the normal stresses near the top and bottom of the beam when the gradient index increases but reduces those in the middle part of the beam. As mentioned, when r increases, the new FG-GPLRC beam can become a sandwich structure, and the stress distributions in Figure 5 prove the presence of the smooth and continuous variation of stresses throughout the entire cross-section. increases, nonlinear change occurs and become more remarkable because a greater r indicates that more GPLs are concentrated on the top and bottom surfaces, which makes a significantly uneven GPL distribution and produces the graded mechanical properties of the FG-GPLRCs. In addition, the increasingly concentrated dispersion of GPLs to the surfaces increases the normal stresses near the top and bottom of the beam when the gradient index increases but reduces those in the middle part of the beam. As mentioned, when r increases, the new FG-GPLRC beam can become a sandwich structure, and the stress distributions in Figure 5 prove the presence of the smooth and continuous variation of stresses throughout the entire cross-section. Meanwhile, the effects of the total weight fractions of the GPLs on the normal stress distributions are presented in Figure 5. The cases of r = 1 and r = 4 are considered in Figure 6a Meanwhile, the effects of the total weight fractions of the GPLs on the normal stress distributions are presented in Figure 5. The cases of r = 1 and r = 4 are considered in Figure 6a,b, respectively. The variations of normal stresses are linear when g 0 GPL = 0 for a pure polymer beam, but they nonlinearly change when g 0 GPL > 0 because the uneven distribution of GPLs produces the gradient variation of material properties. Further addition of GPLs increases the normal stresses near the top and bottom surfaces of the beam and reduces those in the internal part far away from the two surfaces. Meanwhile, the effects of the total weight fractions of the GPLs on the normal stress distributions are presented in Figure 5. The cases of r = 1 and r = 4 are considered in Figure 6a    maximum shear stress appears near the mid-surface of the beam. The curves in Figure 7 show that the change in gradient index significantly affects the shear stress distribution. When the gradient index increases, the maximum shear stress moves to the beam surfaces.      Figure 8 shows the effects of the total weight fraction on the shear stress distributions. As observed in Figure 8, when the gradient index is fixed, the increase in GPL weight fraction increases the shear stresses near the top and bottom surfaces and reduces those in the center portion of the beam.

Elastic Vibration Behavior
Tables 8-10 show the dimensionless fundamental frequencies of FG-GPLRC beams with different gradient indices and weight fractions. The total weight fraction of the GPL is 0-1.5% with an increment of 0.25%. Three immovable boundary conditions are considered.
The numerical results in Tables 8-10 indicate that the change in gradient index can significantly affect the fundamental frequency and consequently affect the elastic vibration characteristics of the composite beam. The effect becomes increasingly obvious when the total weight fraction increases. For a fixed concentration of GPLs in the FG-GPLRC beam, when the gradient index increases, the dimensionless frequency increases. As mentioned, the increase in gradient index implies that more GPLs are dispersed near the top and bottom surfaces of the beam. Therefore, the above phenomenon proves that dispersing more GPL nanofillers near the top and bottom surfaces of the beam is most effective in enhancing the beam stiffness.

Elastic Vibration Behavior
Tables 8-10 show the dimensionless fundamental frequencies of FG-GPLRC beams with different gradient indices and weight fractions. The total weight fraction of the GPL is 0-1.5% with an increment of 0.25%. Three immovable boundary conditions are considered. The numerical results in Tables 8-10 indicate that the change in gradient index can significantly affect the fundamental frequency and consequently affect the elastic vibration characteristics of the composite beam. The effect becomes increasingly obvious when the total weight fraction increases. For a fixed concentration of GPLs in the FG-GPLRC beam, when the gradient index increases, the dimensionless frequency increases. As mentioned, the increase in gradient index implies that more GPLs are dispersed near the top and bottom surfaces of the beam. Therefore, the above phenomenon proves that dispersing more GPL nanofillers near the top and bottom surfaces of the beam is most effective in enhancing the beam stiffness.
The curves in Figure 9 reflect the relationship between the fundamental frequencies and the weight fractions. The results in Figure 9 indicate that the fundamental frequencies increase when the GPL concentration increases regardless of the gradient index. The aforementioned phenomenon clearly indicates that low-content GPL additions significantly improve the stiffness of the beam.  Figure 10 shows the change in dimensionless fundamental frequencies with respect to the gradient index. The total weight fraction of GPLs is 1.25%. Figure 10a shows the results of the thick beam (L/h = 5), and Figure 10b is related to the thin one (L/h = 20). We can reach a contrary conclusion with that of the bending behaviors: for the thick beam, when the gradient index increases, the fundamental frequencies first increase and subsequently decrease. There is a suitable r value that can result in a peak fundamental frequency, which indicates that the FG-GPLRC beam with this r value is the stiffest. However, this value is not constant and related to the geometric, weight fraction and boundary conditions.
For the thin beam (L/h = 20), the fundamental frequency increases with the increase in gradient index. A remarkable variation is observed when r is relatively small, and the variation becomes slight when r increases. This result demonstrates that the frequency is more sensitive to the gradient index when the gradient index is relatively small.  which indicates that the FG-GPLRC beam with this r value is the stiffest. However, this value is not constant and related to the geometric, weight fraction and boundary conditions. For the thin beam (L/h = 20), the fundamental frequency increases with the increase in gradient index. A remarkable variation is observed when r is relatively small, and the variation becomes slight when r increases. This result demonstrates that the frequency is more sensitive to the gradient index when the gradient index is relatively small. As mentioned before, the multilayer structures are introduced to solve the manufacturing problems caused by the distribution laws in practice. It is expected that as a total layer number increase, the functionally graded graphene reinforced multilayer beam will have the some static and vibration response to the monolithic FG-GPLRC beams. Therefore, we calculated the static displacements and fundamental frequencies of the multilayered nanocomposite beams with various layer numbers and made a comparison with those from the monolithic ones. It is worth to note that  . We can reach a contrary conclusion with that of the bending behaviors: for the thick beam, when the gradient index increases, the fundamental frequencies first increase and subsequently decrease. There is a suitable r value that can result in a peak fundamental frequency, which indicates that the FG-GPLRC beam with this r value is the stiffest. However, this value is not constant and related to the geometric, weight fraction and boundary conditions.
For the thin beam (L/h = 20), the fundamental frequency increases with the increase in gradient index. A remarkable variation is observed when r is relatively small, and the variation becomes slight when r increases. This result demonstrates that the frequency is more sensitive to the gradient index when the gradient index is relatively small.
As mentioned before, the multilayer structures are introduced to solve the manufacturing problems caused by the distribution laws in practice. It is expected that as a total layer number increase, the functionally graded graphene reinforced multilayer beam will have the some static and vibration response to the monolithic FG-GPLRC beams. Therefore, we calculated the static displacements and fundamental frequencies of the multilayered nanocomposite beams with various layer numbers and made a comparison with those from the monolithic ones. It is worth to note that the weight fractions of each layer GPL-reinforced nanocomposites are determined by substituting the corresponding z-coordinate of mid-plane of each layer into Equation (8). In this example, the length-to-height ratio is set to be 5, and the total weight fraction of GPLs is 1.0%. The clamped-clamped boundary condition is taken into account. Table 11 lists the maximum displacements and fundamental frequencies of GPL-reinforced multilayer beams for different gradient index r. One can observe that when the total layer number is beyond 12, the results of multilayer beams are very closed to the results of monolithic ones, implying that a multilayer GPLRC beam with 12 or more layers is an excellent approximation for an ideal functionally graded beam structure with a continuous and smooth variation in both material composition and properties. In addition, this result also demonstrates that the new FG-GPLRC structures proposed in present work can be realized easily using the multilayer alternatives in the manufacturing processing.

Conclusions
A new distribution law based on the error function has been proposed and is further employed to introduce a novel FG-GPLRC. The GPL distributions in the thickness direction can be adjusted by a gradient index to maintain a constant total GPL weight fraction. The modified Halpin-Tsai micromechanics model is used to evaluate the effective Young's modulus of the FG-GPLRCs.
A computational modeling based on an improved third-order shear deformation theory and Chebyshev-Ritz methodology is developed. The comparisons prove that the developed modeling is more accurate and efficient. The bending deflections, stresses, and natural frequencies of a novel FG-GPLRC beam are investigated using the proposed modeling. The parameter studies are performed, and some conclusions are summarized based the numerical results: (1) With respect to the GPL distribution patterns, low-content GPL additions can significantly improve the stiffness and bending resistance of the beam. (2) For the thick FG-GPLRC beam, the new distribution law can produce the stiffest FG-GPLRC beams with the lowest bending displacements and highest fundamental frequencies. In other words, by adjusting the gradient index, the most optimized distribution law for the new FG-GPLRCs can be found, which results in an FG-GPLRC beam with the greatest bending resistance and vibration stiffness. (3) For the thin FG-GPLRC beam, the increase in gradient index reduces the bending displacements and increases the fundamental frequencies. In addition, if the gradient index is beyond 1.6, the new FG-GPLRCs exhibit better capabilities in the static and vibration analysis than the common FG-X ones. (4) A multilayer GPLRC beam with 12 or more layers is an ideal alternative structure for fabricating the new FG-GPLRC structures.