Free and Forced Vibration Analyses of Functionally Graded Graphene-Nanoplatelet-Reinforced Beams Based on the Finite Element Method

The finite element method (FEM) is used to investigate the free and forced vibration characteristics of functionally graded graphene-nanoplatelet-reinforced composite (FG-GPLRC) beams. The weight fraction of graphene nanoplatelets (GPLs) is assumed to vary continuously along the beam thickness according to a linear, parabolic, or uniform pattern. For the FG-GPLRC beam, the modified Halpin–Tsai micromechanics model is used to calculate the effective Young’s modulus, and the rule of mixture is used to determine the effective Poisson’s ratio and mass density. Based on the principle of virtual work under the assumptions of the Euler–Bernoulli beam theory, finite element formulations are derived to analyze the free and forced vibration characteristics of FG-GPLRC beams. A two-node beam element with six degrees of freedom is adopted to discretize the beam, and the corresponding stiffness matrix and mass matrix containing information on the variation of material properties can be derived. On this basis, the natural frequencies and the response amplitudes under external forces are calculated by the FEM. The performance of the proposed FEM is assessed, with some numerical results obtained by layering method and available in published literature. The comparison results show that the proposed FEM is capable of analyzing an FG-GPLRC beam. A detailed parametric investigation is carried out to study the effects of GPL weight fraction, distribution pattern, and dimensions on the free and forced vibration responses of the beam. Numerical results show that the above-mentioned effects play an important role with respect to the vibration behaviors of the beam.


Introduction
Graphene is a single two-dimensional honeycomb lattice structure composed of carbon atoms and the basic constituent of other graphite materials with varied dimensions [1]. Lee et al. [2] first measured the elastic properties and internal fracture strength of monolayer graphene by the nanoindentation method under atomic force microscopy (AFM). They found that graphene has a Young's modulus of 1 ± 0.1 TPA and internal strength of 130 ± 10 GPa, which are superior to those of other materials. Graphene has been used by many researchers as a reinforcing material to improve the comprehensive properties of nanocomposites due to its thin thickness, strength, and thermal conductivity [3]. Rafiee et al. [4] carried out experimental studies on different types of epoxy nanocomposites and found that after 0.1% graphene nanoplatelets (GPLs) were added, the Young's modulus of GPLs/epoxy resin composites increased by 31%, and the tensile strength increased by 40% compared with that of the original epoxy resin matrix composites. In rotational speed on the natural frequencies. Based on the third-order shear deformation theory, Li et al. [27] investigated the primary and secondary resonances of FG-GPLRC beams. They examined the effects of GPL parameters on the nonlinear response of the primary, secondary, and combination resonances and found that the addition of a very low-weight fraction of GPLs significantly reduced the resonances of the beams. Based on Hamilton's principle and the nonlinear von Karman strain-displacement relationship, Wang et al. [28] studied the nonlinear transient dynamic response of FG-GPLRC doublecurved shallow shells under time-varying blast loads. According to their reported results, the nonlinear transient response of the shell is considerably affected by many factors, such as the temperature field, GPL parameters, shell parameters, etc. Many achievements have been made with respect to analysis of the mechanical performance of FG-GPLRC structures, but few studies have been reported on FG-GPLRC structures using the finite element method (FEM).
As a powerful numerical method, the FEM is flexible in terms of application to complex geometric configurations and various physical problems. The advantages of the FEM are well-suited for analysis of FG-GPLRC structures. In the present study, the free and forced vibration responses of FG-GPLRC beams were investigated using the FEM. We assumed that the weight fraction of GPLs varies continuously along the beam thickness according to a linear, parabolic, or uniform pattern. The effective Young's modulus was calculated according to the modified Halpin-Tsai micromechanics model, whereas the effective Poisson's ratio and mass density were determined by the rule of mixture. Finite element formulations for the analysis of the FG-GPLRC beams were derived based on the principle of virtual work under the assumptions of the Euler-Bernoulli beam theory. Thus, the corresponding stiffness matrix and mass matrix containing information on the variation of material properties can be derived. The proposed finite element model was verified by ABAQUS based on the layering method and previously published results, confirming agreement. Finally, the effects of the GPL weight fraction, distribution pattern, and dimensions on the free and forced vibration responses of the beam were investigated. Figure 1 represents an FG beam (length, l; width, b; thickness, h) in which GPLs are distributed in the matrix according to a certain law, with the weight fraction of GPLs (WGPL) changing continuously along the thickness of the beam. Four distribution patterns of GPLs were considered, as shown in Figure 2 [14,15]. (GPLRC) cylindrical shells based on first-order shear deformation theory. Using the appropriate technique of Navier solution, an analytical solution was provided to obtain the natural frequencies, as well as the number of modes. Niu et al. [26] investigated the free vibration characteristics of rotating pretwisted FG-GPLRC cylindrical shell plates under cantilever conditions using the Chebyshev-Ritz method and discussed the effects of factors such as GPL parameters, pretorsion angle, and rotational speed on the natural frequencies. Based on the third-order shear deformation theory, Li et al. [27] investigated the primary and secondary resonances of FG-GPLRC beams. They examined the effects of GPL parameters on the nonlinear response of the primary, secondary, and combination resonances and found that the addition of a very low-weight fraction of GPLs significantly reduced the resonances of the beams. Based on Hamilton's principle and the nonlinear von Karman strain-displacement relationship, Wang et al. [28] studied the nonlinear transient dynamic response of FG-GPLRC double-curved shallow shells under time-varying blast loads. According to their reported results, the nonlinear transient response of the shell is considerably affected by many factors, such as the temperature field, GPL parameters, shell parameters, etc. Many achievements have been made with respect to analysis of the mechanical performance of FG-GPLRC structures, but few studies have been reported on FG-GPLRC structures using the finite element method (FEM). As a powerful numerical method, the FEM is flexible in terms of application to complex geometric configurations and various physical problems. The advantages of the FEM are well-suited for analysis of FG-GPLRC structures. In the present study, the free and forced vibration responses of FG-GPLRC beams were investigated using the FEM. We assumed that the weight fraction of GPLs varies continuously along the beam thickness according to a linear, parabolic, or uniform pattern. The effective Young's modulus was calculated according to the modified Halpin-Tsai micromechanics model, whereas the effective Poisson's ratio and mass density were determined by the rule of mixture. Finite element formulations for the analysis of the FG-GPLRC beams were derived based on the principle of virtual work under the assumptions of the Euler-Bernoulli beam theory. Thus, the corresponding stiffness matrix and mass matrix containing information on the variation of material properties can be derived. The proposed finite element model was verified by ABAQUS based on the layering method and previously published results, confirming agreement. Finally, the effects of the GPL weight fraction, distribution pattern, and dimensions on the free and forced vibration responses of the beam were investigated. Figure 1 represents an FG beam (length, l; width, b; thickness, h) in which GPLs are distributed in the matrix according to a certain law, with the weight fraction of GPLs (WGPL) changing continuously along the thickness of the beam. Four distribution patterns of GPLs were considered, as shown in Figure 2 [14,15].       GPLs changes from the  maximum value on the top surface to the minimum value on the bottom surface,  along with the thickness, i.e.,

Functionally Graded Graphene-Nanoplatelet-Reinforced Beam
(b) Parabolic distribution pattern with rich GPLs in the surface (PDPRS): The weight fraction of GPLs is the highest at the top and bottom, lowest in the middle plane, and symmetric about the neutral axis, i.e., (c) Parabolic distribution pattern with rich GPLs in the middle plane (PDPRM): The weight fraction of GPLs is the lowest at the top and bottom, highest in the middle plane, and symmetric about the neutral axis, i.e., (d) Uniform distribution pattern (UDP): The weight fraction of GPLs is evenly distributed in the beam, forming an isotropic, homogeneous beam, i.e., where λ i (i = 1, 2, 3, 4) is a gradient indicator controlling the distribution of the weight fraction of GPLs along the thickness direction; and W 0 GPL is a specific weight fraction, which is taken as W 0 GPL = 1% in this paper. W t GPL represents the total weight fraction of GPLs in the whole beam and can be determined by integrating W GPL (z) along the thickness, i.e., W t GPL = h/2 −h/2 W GPL (z) dz. After integrating Equations (1)-(4), the relationship between W t GPL and λ i can be obtained, as listed in Table 1. The variation in GPL weight fraction along the thickness in the four distribution patterns is shown in Figure 3 when the overall weight fraction (W t GPL ) is 1%. Table 1. Relationship between gradient factor λ i and W t GPL .
The effective Young's modulus (E C ) of graphene-reinforced composites can be calculated according to the modified Halpin-Tsai micromechanics model [13], i.e., in which: where E M and E GPL are the Young's moduli of the polymer matrix and GPLs, respectively; V GPL is the volume fraction of GPLs; and l GPL , b GPL , and t GPL denote the length, width, and thickness of GPLs, respectively.  The effective Young's modulus ( C E ) of graphene-reinforced composites can be calculated according to the modified Halpin-Tsai micromechanics model [13], i.e., in which: where M E and GPL E are the Young's moduli of the polymer matrix and GPLs, respectively; GPL V is the volume fraction of GPLs; and lGPL, bGPL, and tGPL denote the length, width, and thickness of GPLs, respectively.
The effective Poisson's ratio ( C ν ) and mass density ( C ρ ) of graphene-reinforced composites are determined according to the rule of mixture [13], i.e., (1 ) The effective Poisson's ratio (ν C ) and mass density (ρ C ) of graphene-reinforced composites are determined according to the rule of mixture [13], i.e., in which: where ν GPL , ρ GPL , and W GPL are the Poisson's ratio, mass density, and weight fraction of GPLs, respectively; and ν M and ρ M are the Poisson's ratio and mass density of the matrix, respectively.

Governing Equations
Based on the Euler-Bernoulli beam theory, the axial and transverse displacements at any point in the beam can be expressed as [29]: where u 0 and w 0 are the axial displacement and transverse displacement of any point on the neutral axis of the beam, respectively, and t represents time. Then, the displacement vector, {γ s }, can be expressed as: Considering a small deformation, the displacement-strain relationship can be expressed as: where ε xx is the axial strain. Based on the linear elastic theory, the relationship between stress and strain follows Hooke's law, and the strain-stress constitutive equation can be expressed as: where E(z) is the Young's modulus, and σ xx is the axial stress.

Dynamic Equations
The equations that govern the dynamic response of a structure can be obtained based on the principle of virtual work [30,31]. For the FG-GPLRC beam, the virtual work equation becomes: where W S and W C are the sums of the work done by the stress and damping force, respectively; W I is the sum of the work done by the inertia force; and W F is the sum of the work done by the external load. The work done by the stress on the virtual strain is: Substituting Equations (16) and (17) into Equation (19) yields: The work done by the damping forces on the virtual displacement is: where c d is the material damping parameter. Substituting Equation (15) into Equation (21) yields: The work done by the inertia force on the virtual displacement is: Substituting Equation (15) into Equation (23) yields: The work done by the external force on the virtual displacement is:

Beam Element for FG Beam Reinforced with GPLs
As shown in Figure 4, the FG-GPLRC beam is discretized by a two-node beam element with six degrees of freedom.  The components of element displacements at the mid plane of a beam can be expressed as [32]: where i U and i W denote the axial and transverse displacements of node i; i θ is the slope of node i;     The global stiffness matrix is expressed as: The components of element displacements at the mid plane of a beam can be expressed as [32]: where U i and W i denote the axial and transverse displacements of node i; θ i is the slope of node i; U j and W j denote the axial and transverse displacements of node j, respectively; θ j is the slope of node j; N 3i−2 and N 3j−2 are Lagrange's linear interpolation polynomials; and N 3i−1 , N 3j−1 , N 3i , and N 3j are first-order Hermite's interpolation polynomials, i.e., Because the virtual displacements are arbitrary, the combination of Equations (20), (22), (24), (25), and (26)-(33) yields: ..
where [M] is the global mass matrix, [K] is the global stiffness matrix, [C] is the damping matrix, F (t) is the nodal load vector, and r (t) is the displacement vector. The global stiffness matrix is expressed as: where numel denotes the total number of elements; and [K] e is the element stiffness matrix, which can be expressed as follows: (36) in which: Considering Equation (5), substituting Equations (37) and (38) into Equation (36) yields the integrand that contains functions of x. Thus, integration of the integrand in x over elemental length l e and the element stiffness matrix [K] e becomes: in which: where explicit forms of C jl and C jw (j = 1, 2, 3) are given in Appendix A. The global mass matrix is expressed as: where [M] e is the element mass matrix, which can be expressed as: in which: Similarly, considering Equation (11), substituting Equations (45) and (46) into Equation (44) yields the integrand that contains functions of x. Thus, integration of the integrand in x over elemental length l e and the element mass matrix [M] e becomes: in which: where explicit forms of C j (j = 1, 2, 3) are given in Appendix A.

Free Vibration Analysis
When the damping effect and external force are not considered, Equation (34) can be simplified as: [M] ..
This is the free vibration equation of the system, also known as the dynamic characteristic equation.
Suppose that: where {φ n } is the mode shape vector for the nth-order mode of vibration, ω n is the frequency of vibration of the φ n vector, t is the time variable, and t 0 is the time constant determined by the initial conditions. Substituting Equation (52) into Equation (51) yields: According to the theory of homogeneous equations, Equation (53) can be expressed as: Thus, ω n and {φ n } can be obtained by solving Equations (54) and (53), respectively.

Transient Analysis
Equation (34) is the dynamic equation of the system. The modal superposition method [33,34] is adopted to analyze the modal transient response of the system. By projecting the physical model into the eigenmodal coordinate system, a set of mutually independent single-degree-of-freedom equations of motion can be obtained, i.e., ..
The general form of the solution is: At this point, the dynamic response under modal coordinates can be obtained, and the physical coordinate response of the system can be obtained through modal superposition, i.e.,:

Numerical Results
The above framework offers the same procedures as those used in the standard FEM to analyze FG-GPLRC beams with elastic models and the procedures programmed with FORTRAN language through subroutine UEL of ABAQUS.
Unless otherwise stated, the following conditions are used to investigate the effects of factors, such as the total weight fraction, distribution patterns, and size of GPLs, on the free vibration and transient response of the FG-GPLRC beam. Epoxy resin is chosen as the composite matrix, and GPLs are used as material reinforcement for this study. The material-related parameters are shown in Table 2 [35]. The FG beam length (l = 0.2 m) and thickness (h) vary with the length-to-slenderness ratio. The support conditions are selected to be simply supported at both ends, as shown in Figure 5. The following quantities can be defined as follows: Table 2. Material-related parameters.
Unless otherwise stated, the following conditions are used to investigate the effects of factors, such as the total weight fraction, distribution patterns, and size of GPLs, on the free vibration and transient response of the FG-GPLRC beam. Epoxy resin is chosen as the composite matrix, and GPLs are used as material reinforcement for this study. The material-related parameters are shown in Table 2 [35]. The FG beam length (l = 0.2m) and thickness (h) vary with the length-to-slenderness ratio. The support conditions are selected to be simply supported at both ends, as shown in Figure 5. The following quantities can be defined as follows: Dimensionless fundamental frequency: where F0 is the maximum value of the dynamic load, and 3 12 bh I = .

Convergence and Validation
The two-node beam element with six degrees of freedom is used to discretize the FG-GPLRC beam with a length-to-thickness ratio of 25 and a total GPLs weight fraction of 1%. The number of elements varies from 5 to 60 in order to investigate the convergence. Table 3 illustrates the calculated first-order dimensionless free vibration frequencies for four distribution patterns of GPLs. As shown in Table 3, the frequencies decrease continuously and converge to a certain value with an increase in the number of elements from 5 to 20, and the frequency values remain unchanged after the number of elements exceeds 20. According to the results, the number of elements has an impact on the calculation accuracy, and the calculated values correspond with the increase in the number of elements. Considering the calculation efficiency and accuracy, the number of elements is taken as 20 for all cases in the following calculation. Dimensionless fundamental frequency: Dimensionless response vibration amplitude: where F 0 is the maximum value of the dynamic load, and I = bh 3 12 .

Convergence and Validation
The two-node beam element with six degrees of freedom is used to discretize the FG-GPLRC beam with a length-to-thickness ratio of 25 and a total GPLs weight fraction of 1%. The number of elements varies from 5 to 60 in order to investigate the convergence. Table 3 illustrates the calculated first-order dimensionless free vibration frequencies for four distribution patterns of GPLs. As shown in Table 3, the frequencies decrease continuously and converge to a certain value with an increase in the number of elements from 5 to 20, and the frequency values remain unchanged after the number of elements exceeds 20. According to the results, the number of elements has an impact on the calculation accuracy, and the calculated values correspond with the increase in the number of elements. Considering the calculation efficiency and accuracy, the number of elements is taken as 20 for all cases in the following calculation. Owing to the lack of suitable data available in the published literature for validation analysis, the solutions generated by finite element software ABAQUS for FG-GPLRC beams are chosen as the reference solutions to verify the proposed method. To solve the FG material inhomogeneity problem, a layering method [36][37][38] is used to divide the FG material into a certain number of homogeneous layers with the same material parameters in the same layer but with different material parameters in each layer. Figure 6 represents a schematic diagram of FEM mesh based on the layering method for FG-GPLRC beams. In the finite element model, the number of homogeneous layers is taken as 20, and the layers are discretized by the four-node S4R shell element with a size ratio of a: b: t = 25:10:1.  Owing to the lack of suitable data available in the published literature for validatio analysis, the solutions generated by finite element software ABAQUS for FG-GPLR beams are chosen as the reference solutions to verify the proposed method. To solve th FG material inhomogeneity problem, a layering method [36][37][38] is used to divide the F material into a certain number of homogeneous layers with the same material paramete in the same layer but with different material parameters in each layer. Figure 6 represen a schematic diagram of FEM mesh based on the layering method for FG-GPLRC beam In the finite element model, the number of homogeneous layers is taken as 20, and th layers are discretized by the four-node S4R shell element with a size ratio of a: b: t = 25:10: Both the proposed method based on a two-node beam element with six degrees freedom and the ABAQUS solution based on a layering method are used to solve th problem. In addition, further verification results according to a numerical model are give in reference [35]. Under different GPL weight fractions, slenderness ratios, and GPL di tribution patterns for a simply supported FG-GPLRC beam, the dimensionless fundamen tal frequencies obtained by the proposed method, ABAQUS, and Wang [35] are presente in Table 4. A comparison of the results shows that the solutions of the proposed metho are in agreement with those of ABAQUS and Wang [35]. The maximum deviation of th presented solutions with those of Wang [35]    Both the proposed method based on a two-node beam element with six degrees of freedom and the ABAQUS solution based on a layering method are used to solve this problem. In addition, further verification results according to a numerical model are given in reference [35]. Under different GPL weight fractions, slenderness ratios, and GPL distribution patterns for a simply supported FG-GPLRC beam, the dimensionless fundamental frequencies obtained by the proposed method, ABAQUS, and Wang [35] are presented in Table 4. A comparison of the results shows that the solutions of the proposed method are in agreement with those of ABAQUS and Wang [35]. The maximum deviation of the presented solutions with those of Wang [35] is 0.576% (l/h = 25, PDPRS, W t GPL = 1%) and the maximum deviation of the presented solutions with those of ABAQUS is 0.920% (l/h = 20, PDPRS, W t GPL = 1%), which indicates that the proposed method highly accurate and is trustworthy. In addition, fewer elements are used to model the beam with the proposed method than the ABAQUS solution based on a laying method, indicating the efficiently of the proposed method. The investigation of transient response is based on the modal superposition method, which first required determination of the free vibration frequency. Because no similar results have been published in the literature with respect to the transient response of the problem, in this study, we directly compare the amplitude of dimensionless transient response in the middle span of FG-GPLRC beams with the ABAQUS results. When W t GPL = 0.5%, the dynamic load applied in the span of the beam adopts the following form: As shown in Figure 7, the results of the proposed method agree with the ABAQUS data under four distribution patterns of GPLs. Therefore, the reliability of the proposed method is further verified.

Free Vibration
The relationship between the first three-order dimensionless frequencies and the total weight fraction of GPLs is shown in Figure 8 for simply supported FG-GPLRC beams with four different GPL distribution patterns and a slenderness ratio of l/h = 20. The dimensionless frequency increases significantly after adding a small amount of graphene. In order to explore the effect of the weight fraction and distribution pattern of GPLs on the frequency, Table 5 lists the specific values of the first three dimensionless frequencies of FG-GPLRC beams with 0%, 0.5%, and 1% weight fractions of GPLs, respectively.
For the LDP of GPLs, the first three-order dimensionless frequencies increase by 54.006%, 53.815%, and 53.129% with the addition of a 0.5% weight fraction of GPLs and by 89.248%, 88.881%, and 87.740% with the addition of a 1.0% weight fraction of GPLs. Like the LDP of GPLs, PDPRS, PDPRM, and UDP are subject to a similar effect of the weight fraction of GPLs on the frequency. The addition of a small amount of graphene can greatly increase the frequency because graphene with a high Young's modulus increases the stiffness of the beam. Under the same GPL distribution pattern, when the weight fraction of GPLs is the same, the growth rate of different order frequencies is almost the same, which indicates that there is no obvious difference in terms of the effect of the GPL weight fraction on different order frequencies.
As shown in Figure 7, the results of the proposed method agree with the ABAQUS data under four distribution patterns of GPLs. Therefore, the reliability of the proposed method is further verified.   tal weight fraction of GPLs is shown in Figure 8 for simply supported FG-GPLRC beams with four different GPL distribution patterns and a slenderness ratio of l/h = 20. The dimensionless frequency increases significantly after adding a small amount of graphene. In order to explore the effect of the weight fraction and distribution pattern of GPLs on the frequency, Table 5 lists the specific values of the first three dimensionless frequencies of FG-GPLRC beams with 0%, 0.5%, and 1% weight fractions of GPLs, respectively. For the LDP of GPLs, the first three-order dimensionless frequencies increase by 54.006%, 53.815%, and 53.129% with the addition of a 0.5% weight fraction of GPLs and by 89.248%, 88.881%, and 87.740% with the addition of a 1.0% weight fraction of GPLs. Like the LDP of GPLs, PDPRS, PDPRM, and UDP are subject to a similar effect of the weight fraction of GPLs on the frequency. The addition of a small amount of graphene Under the same GPL weight increment, such as a W t GPL increase from 0.0% to 0.5%, the growth rate of different order frequencies is about 103% for the PDPRS of GPLs, 65% for the UDP of GPLs, 53% for the LDP of GPLs, and 43% for the PDPRM of GPLs. The distribution pattern of GPLs considerably affects the frequency of FG-GPLRC beams. When the addition of GPL weight fraction is the same, the PDPRS of GPLs has the greatest influence on the frequency, whereas the PDPRM of GPLs has the least influence on the frequency. This means that distributing more GPLs on the upper and lower surfaces of the beam is the most effective method to increase the flexural stiffness and the frequency of the beam.
In order to study the effect of the geometric shape and size of GPLs on the dynamic performance of beams, three sizes of GPLs are considered, as shown in Figure 9. The free vibration characteristics of simply supported FG-GPLRC with a slenderness ratio of l/h = 20 are studied by adding GPLs of varying sizes, with the GPL weight fraction kept constant (W t GPL = 1%). In this example, the length of GPLs (l GPL ) remains unchanged, and the aspect ratio l GPL /b GPL = 1 corresponds to square GPLs, whereas the aspect ratios l GPL /b GPL = 2 and l GPL /b GPL = 3 correspond to two rectangular GPLs of different sizes. A higher length-thickness ratio (l GPL /t GPL ) indicates that the GPLs are thinner and consist of fewer monolayers of graphene.  Figure 10 shows the variation of the first-order dimensionless fundament quency of simply supported FG-GPLRC beams with varying sizes of GPLs. Figu shows that the fundamental frequency of the beam increases substantially whe length-thickness ratio (lGPL/tGPL) of the GPLs increases to 1000 and then increase slightly as the length-thickness ratio increases further. Whereas the total weight fr of GPLs is kept constant, thinner GPLs result in a larger contact area between the and the matrix, so the GPLs are fully in contact with the matrix material, highlightin high elastic modulus properties and considerably increasing the frequency of the accordingly. However, as GPLs become thinner, they are prone to curling beh changing from a planar configuration to a curved configuration, which cannot fu ploit their own properties and those of the composite material. Figure 10 also show under the four different distribution patterns of GPLs, the frequency of a square enhanced FG beam is higher than that of a rectangular-GPL-enhanced FG beam be square GPLs have a larger specific surface area, resulting in superior performan hancement.  Figure 9. Schematic diagram of graphene nanosheets. Figure 10 shows the variation of the first-order dimensionless fundamental frequency of simply supported FG-GPLRC beams with varying sizes of GPLs. Figure 10 shows that the fundamental frequency of the beam increases substantially when the length-thickness ratio (l GPL /t GPL ) of the GPLs increases to 1000 and then increases only slightly as the length-thickness ratio increases further. Whereas the total weight fraction of GPLs is kept constant, thinner GPLs result in a larger contact area between the GPLs and the matrix, so the GPLs are fully in contact with the matrix material, highlighting their high elastic modulus properties and considerably increasing the frequency of the beam accordingly. However, as GPLs become thinner, they are prone to curling behavior, changing from a planar configuration to a curved configuration, which cannot fully exploit their own properties and those of the composite material. Figure 10 also shows that under the four different distribution patterns of GPLs, the frequency of a square-GPL-enhanced FG beam is higher than that of a rectangular-GPL-enhanced FG beam because square GPLs have a larger specific surface area, resulting in superior performance enhancement. Figure 11 represents a simply supported FG-GPLRC beam with a rectangular cross section applied by a dynamic load over a period of time at mid span of the beam. Four dynamic load types are selected to investigate the effects of the graphene distribution pattern, mass, and size on the transient response, as shown in Figure 12, with a maximum value of dynamic load of F 0 = 1 KN, T = 0.2 s. Fu and Chung [39] tested the loss tangent of epoxy resin under low-frequency external forces, with an obtained value of 0.030 ± 0.007. The same result was obtained by Fereidoon through experimental testing as 0.036 [40]. Thus, in the following numerical study, the low modal damping ratio (ζ) is assumed to be 0.03. The amplitude of the response to dynamic loads in the mid span of an FG-GPLRC beam is studied as an eigenvalue. Figure 13 compares the effects of four different GPL distribution patterns on the amplitude of the spanwise dynamic response of a simply supported FG-GPLRC beam under four dynamic loads. Due to the presence of damping in the system, the response amplitude becomes progressively smaller with time. Figure 13 shows that under the PDPRS of GPLs, the beam has the lowest dynamic response amplitude and the shortest response vibration period, whereas under the PDPRM of GPLs, the beam has the highest dynamic response amplitude and the longest response vibration period because PDPRS can distribute GPLs more on the upper and lower surfaces of the beam to effectively improve the flexural stiffness, which results in a higher vibration frequency and a shorter vibration period. changing from a planar configuration to a curved configuration, which cannot fully exploit their own properties and those of the composite material. Figure 10 also shows that under the four different distribution patterns of GPLs, the frequency of a square-GPLenhanced FG beam is higher than that of a rectangular-GPL-enhanced FG beam because square GPLs have a larger specific surface area, resulting in superior performance enhancement.      Figure 11 represents a simply supported FG-GPLRC beam with a rectangular cross section applied by a dynamic load over a period of time at mid span of the beam. Four dynamic load types are selected to investigate the effects of the graphene distribution pattern, mass, and size on the transient response, as shown in Figure 12, with a maximum value of dynamic load of F0 = 1KN, T = 0.2s. Fu and Chung [39] tested the loss tangent of epoxy resin under low-frequency external forces, with an obtained value of 0.030 ± 0.007. The same result was obtained by Fereidoon through experimental testing as 0.036 [40]. Thus, in the following numerical study, the low modal damping ratio (ζ) is assumed to be 0.03. The amplitude of the response to dynamic loads in the mid span of an FG-GPLRC beam is studied as an eigenvalue.    Figure 11 represents a simply supported FG-GPLRC beam with a rectangular cross section applied by a dynamic load over a period of time at mid span of the beam. Four dynamic load types are selected to investigate the effects of the graphene distribution pattern, mass, and size on the transient response, as shown in Figure 12, with a maximum value of dynamic load of F0 = 1KN, T = 0.2s. Fu and Chung [39] tested the loss tangent of epoxy resin under low-frequency external forces, with an obtained value of 0.030 ± 0.007. The same result was obtained by Fereidoon through experimental testing as 0.036 [40]. Thus, in the following numerical study, the low modal damping ratio (ζ) is assumed to be 0.03. The amplitude of the response to dynamic loads in the mid span of an FG-GPLRC beam is studied as an eigenvalue.   Figure 13 compares the effects of four different GPL distribution patterns o plitude of the spanwise dynamic response of a simply supported FG-GPLRC be four dynamic loads. Due to the presence of damping in the system, the respon tude becomes progressively smaller with time. Figure 13 shows that under the GPLs, the beam has the lowest dynamic response amplitude and the shortest vibration period, whereas under the PDPRM of GPLs, the beam has the highest response amplitude and the longest response vibration period because PDPR tribute GPLs more on the upper and lower surfaces of the beam to effectively im flexural stiffness, which results in a higher vibration frequency and a shorter period. In order to investigate the effects of total weight fraction of the GPLs on the dynamic response, a simply supported FG-GPLRC beam with the PDPRS of GPLs is considered under four kinds of dynamic load. The total weight fraction of the GPLs ranges from 0.0% to 2.0%, with an increment of 0.5%. Figure 14 represents the relationships between the dynamic responses and total weight fraction of the GPLs. Figure 14 shows that when W t GPL = 0.0%, i.e., a pure epoxy beam, the amplitudes are largest, and the vibration periods are longest. With an increase in the total weight fraction of the GPL, the amplitudes of the dynamic responses decrease, and the vibration period gradually shortens. As described in Section 5.2, adding a small amount of GPLs can increase the flexural stiffness and the frequency of the beam, meaning that the amplitude of the dynamic response of the beam under dynamic load can be suppressed by adding a small amount of GPLs.

Transient Response
Consider dimensionless response vibration amplitude (A MAX ) at mid span of simply supported FG-GPLRC beams with different weight fractions and distribution patterns of GPLs. The total weight fraction of the GPLs ranges from 0.0% to 2.5%, with an increment of 0.25%. Figure 15 represents the relationships between A MAX and W t GPL under different distribution patterns of GPLs and different dynamic loads. The curves in Figure 15 clearly show that adding a small amount of GPLs can significantly reduce the maximum response amplitudes and enhance the bending performance of the beams. In the initial stage of the GPL weight fraction increase (W t GPL from 0.0% to 1.0%), the amplitude (A MAX ) decreases obviously. As the GPL weight fraction continues to increase, the amplitude (A MAX ) continues to decrease, although not obviously. For the same GPL weight fraction, A MAX of the beam with the PDPRS of GPLs decreases more than that of the other three distribution patterns (i.e., LDP, PDPRM, and UDP) under four kinds of dynamic load. It can be concluded that the PDPRS distributing more GPLs on the upper and lower surfaces of the beam is the most effective way to increase the stiffness of the beam, reduce the maximum response amplitude, and enhance the bending performance.
Considering the PDPRS of GPLs achieves excellent performance with respect to transient response, the PDPRS is selected to study the effect of GPL size on transient response. An FG-GPLRC beam with a GPL weight fraction of W t GPL = 1% is subjected to four different types of dynamic forces. Figure 16 shows the relationship between the dynamic response and the aspect ratio of GPLs. As shown in Figure 16, the response amplitude of the beam increases with increased aspect ratio, but the change in the response period is not obvious. The values of GPL length and width are taken from Table 2, i.e., l GPL × b GPL = 2.5 µm × 1.5 µm, and the thickness of GPLs is changed to further investigate the effect of thickness on the dynamic response of the beam. Figure 17 shows the relationship between the dynamic response and the length-to-thickness ratio of GPLs. As the GPLs become thinner, the dynamic response period of the beam becomes shorter, and the maximum amplitude becomes lower. The greatest decrease is observed under rectangular loading, when l GPL /b GPL changes from 500 to 2000, and the maximum amplitude decreases from 0.3384 to 0.2676-a decrease of 26%. This means that adding thinner GPLs can reduce the vibration amplitude of the beam structure under external action. T Time(s) Figure 12. Four types of dynamic loads. Figure 13 compares the effects of four different GPL distribution patterns on the am plitude of the spanwise dynamic response of a simply supported FG-GPLRC beam unde four dynamic loads. Due to the presence of damping in the system, the response ampl tude becomes progressively smaller with time. Figure 13 shows that under the PDPRS o GPLs, the beam has the lowest dynamic response amplitude and the shortest respons vibration period, whereas under the PDPRM of GPLs, the beam has the highest dynam response amplitude and the longest response vibration period because PDPRS can dis tribute GPLs more on the upper and lower surfaces of the beam to effectively improve th flexural stiffness, which results in a higher vibration frequency and a shorter vibratio period. In order to investigate the effects of total weight fraction of the GPLs on the dynam response, a simply supported FG-GPLRC beam with the PDPRS of GPLs is considere under four kinds of dynamic load. The total weight fraction of the GPLs ranges from 0.0% to 2.0%, with an increment of 0.5%. Figure 14 represents the relationships between th dynamic responses and total weight fraction of the GPLs. Figure 14 shows that whe GPL = 0.0%, i.e., a pure epoxy beam, the amplitudes are largest, and the vibration peri ods are longest. With an increase in the total weight fraction of the GPL, the amplitude of the dynamic responses decrease, and the vibration period gradually shortens. As de scribed in Section 5.2, adding a small amount of GPLs can increase the flexural stiffnes and the frequency of the beam, meaning that the amplitude of the dynamic response o the beam under dynamic load can be suppressed by adding a small amount of GPLs.      As a further illustration of the effects of GPL geometry on transient response, the changes in the aspect ratio (l GPL /b GPL ) and the length-thickness ratio (l GPL /t GPL ) are simultaneously considered. Figure 18 shows the variation curves of the maximum dimensionless response vibration amplitude (A MAX ) at mid span of the beam with different sizes and thicknesses of GPLs. The curves in Figure 18 show that in the process of increasing the length-thickness ratio to 500, A MAX declines continuously, with a wide range of decline. However, as the length-thickness ratio exceeds 500 and increases further, the decline in A MAX tends occurs gradually. When the l GPL /t GPL ratio is increased from 500 to 2000 for all operating conditions, the maximum dimensionless response amplitude decreases in the range of about 75% to 85%. When the total weight fraction of GPLs remains unchanged, relatively thin GPLs have a larger contact area between the GPLs and the matrix, which can transfer the load and stress better. When the GPLs become thinner and reach a certain thickness, they so thin that curling behavior occurs, and their advantages cannot be fully exploited to transfer the stress. Therefore, relatively thin GPLs can fully reduce the vibration amplitude of the beam, whereas excessively thin GPLs modulate the reduction effect to a certain extent.
clearly show that adding a small amount of GPLs can significantly reduce the maximum response amplitudes and enhance the bending performance of the beams. In the initia stage of the GPL weight fraction increase ( t GPL W from 0.0% to 1.0%), the amplitude ( M AX A ) decreases obviously. As the GPL weight fraction continues to increase, the amplitude M AX A ) continues to decrease, although not obviously. For the same GPL weight fraction

M AX
A of the beam with the PDPRS of GPLs decreases more than that of the other thre distribution patterns (i.e., LDP, PDPRM, and UDP) under four kinds of dynamic load. I can be concluded that the PDPRS distributing more GPLs on the upper and lower surface of the beam is the most effective way to increase the stiffness of the beam, reduce th maximum response amplitude, and enhance the bending performance. Considering the PDPRS of GPLs achieves excellent performance with respect to tran sient response, the PDPRS is selected to study the effect of GPL size on transient respons An FG-GPLRC beam with a GPL weight fraction of t GPL W = 1% is subjected to four diffe ent types of dynamic forces. Figure 16 shows the relationship between the dynamic re sponse and the aspect ratio of GPLs. As shown in Figure 16, the response amplitude of th beam increases with increased aspect ratio, but the change in the response period is no obvious. The values of GPL length and width are taken from Table 2, i.e., lGPL × bGPL = 2 μm × 1.5 μm, and the thickness of GPLs is changed to further investigate the effect o thickness on the dynamic response of the beam. Figure 17 shows the relationship betwee the dynamic response and the length-to-thickness ratio of GPLs. As the GPLs becom the dynamic response and the length-to-thickness ratio of GPLs. As the GPLs become thinner, the dynamic response period of the beam becomes shorter, and the maximum amplitude becomes lower. The greatest decrease is observed under rectangular loading when lGPL/bGPL changes from 500 to 2000, and the maximum amplitude decreases from 0.3384 to 0.2676-a decrease of 26%. This means that adding thinner GPLs can reduce the vibration amplitude of the beam structure under external action.        Figure 18 also shows that the vibration amplitude of square-GPL-enhanced FG beams is lower than that of rectangular-GPL-enhanced beams for all four types of dynamic forces because square GPLs have a larger specific surface area, so they can transfer stress more effectively, resulting in an improved enhancement effect. This result provides a basis for future preparation of graphene-reinforced materials.

Time(s) Time(s)
(c) (d)      to a certain extent. Figure 18 also shows that the vibration amplitude of square-GPL-enhanced FG beams is lower than that of rectangular-GPL-enhanced beams for all four types of dynamic forces because square GPLs have a larger specific surface area, so they can transfer stress more effectively, resulting in an improved enhancement effect. This result provides a basis for future preparation of graphene-reinforced materials.

Conclusions
Free vibration characteristics and the dynamic behaviors of FG-GPLRC beams were investigated based on the finite element method. The GPL weight fraction varies continuously along the beam thickness according to a linear, parabolic, or uniform pattern. For GPL-reinforced composites, the effective Young's modulus (E C ) is calculated according to the modified Halpin-Tsai micromechanics model, and the effective Poisson's ratio (ν C ) and mass density (ρ C ) are determined by the rule of mixture. Finite element formulations for the vibration analyses of FG-GPLRC beams are derived according to the principle of virtual work under the assumptions of the Euler-Bernoulli beam theory. A two-node beam element with six degrees of freedom is used to discretize the beam. Thus, the corresponding stiffness matrix and mass matrix, which contain information on the variation of material properties along the beam thickness, can be derived, and the equations of motion of the system are obtained for FG-GPLRC beams. On this basis, the natural frequencies of FG-GPLRC beams and the response amplitudes over time under the action of external forces are obtained using the FEM. The effects of four GPL distribution patterns, weight fractions, and dimensions of GPLs on the free vibration characteristics and transient response of the beams are investigated. Based on the numerical results, the following conclusions can be drawn: The performance of the proposed FEM was assessed relative to some numerical results obtained by layering method and available in the published literature. Comparison results indicate that the proposed FEM is highly accuracy and efficiency and is therefore capable of analyzing FG-GPLRC beams.
Numerical results show that adding a small amount of GPLs can considerably increase the natural frequency of the beam and significantly reduce the dynamic amplitude of the beam under forced vibration. The maximum increase in the intrinsic frequency of beams is 170% when the weight fraction of GPLs is increased from 0% to 1%. Figure 14 shows that the dynamic amplitude of the beam under forced vibration can be reduced significantly as the weight fraction of GPLs is increased from 0% to 2%. Because graphene has very high Young's modulus, adding a small amount of GPLs can considerably increase the stiffness of the FG-GPLRC beam and naturally improve its bending performance.
The distribution of GPLs in the beam has a considerable influence on the dynamic performance of the beam. Numerical results reveal that an FG-GPLRC beam with the PDPRS of GPLs has a higher natural frequency and lower dynamic amplitude than the other three distribution patterns, i.e., LDP, PDPRM, and UDP. Under the same GPL weight increment, such as a W t GPL increase from 0.0% to 0.5%, the growth rate of the frequency is about 103% for the PDPRS of GPLs, 65% for the UDP of GPLs, 53% for the LDP of GPLs, and 43% for the PDPRM of GPLs. Figure 13 shows that under the PDPRS of GPLs, the beam has the lowest dynamic response amplitude, which indicates that distributing more GPLs on the upper and lower surfaces of the beam is the most effective method to increase the stiffness of the beam and improve its bending performance.
The geometrical size and shape of GPLs can affect the dynamic performance of the beam. Figure 10 indicates that under the condition a constant weight fraction, an FG-GPLRC beam with square GPLs and thinner GPLs has a higher natural frequency. Similarly, Figure 18 reveals that an FG-GPLRC beam with square GPLs and thinner GPLs less dynamic amplitude than a beam with rectangular GPLs and thicker GPLs. Because square GPLs and thinner GPLs have a larger specific surface areas, they can come into full contact with the matrix material and provide better load and stress transfer. Thus, for an FG-GPLRC beam, square GPLs and thinner GPLs have better potential to improve the flexural performance of the beam than rectangular GPLs and thicker GPLs.
The findings presented herein contribute to improved understanding of the influence of GPL weight fraction, distribution pattern, and dimensions on the mechanical properties of FG-GPLRC beams, which is beneficial with respect to the applications of FG-GPLRC beams in many fields, such as aerospace engineering, chemical engineering, civil engineering, etc. For example, the mechanical properties of beams play an important role in structural design and structural safety assessment. Thus, understanding how to add graphene to reinforced concrete beams for improved mechanical properties is helpful with respect to the application of research results the construction field.

Conflicts of Interest:
The authors declare no conflict of interest.