Nonlinear Buckling Analysis of Functionally Graded Graphene Reinforced Composite Shallow Arches with Elastic Rotational Constraints under Uniform Radial Load

The buckling behavior of functionally graded graphene platelet-reinforced composite (FG-GPLRC) shallow arches with elastic rotational constraints under uniform radial load is investigated in this paper. The nonlinear equilibrium equation of the FG-GPLRC shallow arch with elastic rotational constraints under uniform radial load is established using the Halpin-Tsai micromechanics model and the principle of virtual work, from which the critical buckling load of FG-GPLRC shallow arches with elastic rotational constraints can be obtained. This paper gives special attention to the effect of the GPL distribution pattern, weight fraction, geometric parameters, and the constraint stiffness on the buckling load. The numerical results show that all of the FG-GPLRC shallow arches with elastic rotational constraints have a higher buckling load-carrying capacity compared to the pure epoxy arch, and arches of the distribution pattern X have the highest buckling load among four distribution patterns. When the GPL weight fraction is constant, the thinner and larger GPL can provide the better reinforcing effect to the FG-GPLRC shallow arch. However, when the value of the aspect ratio is greater than 4, the flakiness ratio is greater than 103, and the effect of GPL’s dimensions on the buckling load of the FG-GPLRC shallow arch is less significant. In addition, the buckling model of FG-GPLRC shallow arch with elastic rotational constraints is changed as the GPL distribution patterns or the constraint stiffness changes. It is expected that the method and the results that are presented in this paper will be useful as a reference for the stability design of this type of arch in the future.


Introduction
Graphene-related research and application development are hot topics today, and graphene-based composites have become an important direction in application field of graphene [1][2][3][4]. Graphene and its derivatives have many advantageous material properties, such as high strength, lightweight, and high chemical stability. The tensile strength of graphene can reach 130 GPa and its theoretical elastic modulus is as high as 1 TPa [5]. As a result, graphene and its derivatives can enhance the special properties of a variety of composites. Rafiee et al. [6,7] added 0.1 wt % graphene nanoplatelets (GPLs) and carbon nanotubes (CNTs) into two composites with epoxy matrices. By comparing the mechanical properties of the two reinforced composites, they found that the GPL-reinforced composite had a higher tensile strength and fracture toughness than the CNT-reinforced composite. By thoroughly analyzing  1 where NL is the total number of GPLRC layers and * GPL V is the total volume fraction of the GPLRC, E E and L  and T  are the filler geometry factors defined as  The FG-GPLRC shallow arch with an even number of layers is studied in this paper. The volume fractions VGPL of the kth layer for the four distribution patterns that are shown in where NL is the total number of GPLRC layers and * GPL V is the total volume fraction of the GPLRC, which is defined as where WGPL is the GPL weight fraction and GPL  and m  are the mass densities of GPLs and polymer matrices, respectively.
The elastic modulus of composites with randomly oriented fillers can be calculated by the Halpin-Tsai model [32,33] (6) in which the parameters L  and T  are defined as and L  and T  are the filler geometry factors defined as The FG-GPLRC shallow arch with an even number of layers is studied in this paper. The volume fractions V GPL of the kth layer for the four distribution patterns that are shown in Figure 2 are given by where N L is the total number of GPLRC layers and V * GPL is the total volume fraction of the GPLRC, which is defined as where W GPL is the GPL weight fraction and ρ GPL and ρ m are the mass densities of GPLs and polymer matrices, respectively. The elastic modulus of composites with randomly oriented fillers can be calculated by the Halpin-Tsai model [32,33] (6) in which the parameters η L and η T are defined as and ξ L and ξ T are the filler geometry factors defined as where a GPL , b GPL , and t GPL are the length, width and thickness of GPLs, respectively. a GPL /b GPL , b GPL /t GPL are the width-to-length ratio and width-to-thickness ratio of GPLs, respectively.

Nonlinear Equilibrium Equation
The configuration and the coordinate system of the FG-GPLRC shallow arch are shown in Figure 1. The origin o of the axis system oxy is located at the centroidal axes of the arch. The axis oz coincides with a principle axis of the cross-section, its direction changes along the circumference of the arch and it is always towards the center of the arch, and the axis ox coincides with the centroidal axis. S, Θ, and R are the axial length, half-included angle, and radius of the arch, respectively. The stiffness of the elastic rotational constraints at both ends is k. The longitudinal strain of an arbitrary point Po at the cross-section can be expressed as where ε m is the membrane strain and ε b is the bending strain of the FG-GPLRC shallow arch. According to Donnell's shallow arch theory [18,22], the relationship between the strain and displacement are where ( ) ≡ d( )/dθ, v = v/R and w = w/R are the dimensionless axial and radial displacements, respectively. z is the coordinate of the point P at the cross-section of the axis oz.
According to the virtual work principle, the virtual work of the FG-GPLRC shallow arch with elastic rotational constraints under uniform radial load can be expressed as Substituting Equations (9)-(11) into Equation (12) leads to where N and M are the axial force and bending moment of the arch, which can be expressed as in which, A 11 , B 11 , and D 11 are the stiffness components of the FG-GPLRC shallow arch, which are defined as where z n and z n+1 are the coordinates of the bottom and upper surfaces of the kth GPLRC layer, respectively.
Integrating Equation (13) by parts leads At the same time, the following essential kinematic boundary conditions δ v = 0 and δ w = 0 at θ = ±Θ −M ± k v = 0 at θ = ±Θ (18) need to be satisfied. Hence, the last two terms of Equation (17) vanish. Because δ v = 0 and δ w = 0 are arbitrary by definition, to hold Equation (17), the differential equations of equilibrium have to be satisfied The essential kinematic boundary conditions are Substituting Equation (14) into Equation (15) leads to the bending moment of the FG-GPLRC elastic rotational constraint arch, as Substituting Equations (19) and (22) into Equation (20), the differential equation of equilibrium in the radial direction can be rewritten as v iv where, P is a dimensionless load, which is defined as µ is dimensionless axial force parameter, which is defined as A 11 (25) where κ is the equivalent stiffness of the FG-GPLRC shallow arch. The dimensionless radial displacement v can be obtained by solving Equation (23) under the boundary conditions given by Equation in which, β = µΘ and α is the ratio of the bending stiffness per unit length of the equivalent stiffness κ/S to the stiffness of elastic rotational end constraints k, which can be expressed as Substituting Equations (25) and (26) into (14), and integrating the function over the entire arch, leads to the nonlinear equilibrium equation between the axial force parameter µ and the dimensionless load P, expressed as where in which Γ = 1 + 2α 2α + tan β/β (32) λ is the geometric parameter of FG-GPLRC shallow arch defined as which increases as the arch length and the half-included angle increase but decreases as the thickness of cross section increases. Thus, Equation (28) establishes the relationship between the external load P and the axial compressive force parameter β. By solving Equation (28), the axial compressive force β under a different load P can be obtained. Then, by substituting β into Equation (26), the dimensionless radial displacement v can be obtained. Thus, the P-v curve can be drawn, and the nonlinear buckling behavior can be analyzed.

Limit Instability Buckling
The limit instability points are local maximum and minimum points on the nonlinear equilibrium path of an FG-GPLRC shallow arch, which can be determined by routing calculus. The localized load P can be expressed as an implicit function of the parameter β using Equation (28) as F(P, β) = 0. Hence, the limit loads can be obtained by setting which leads to the equation of equilibrium between the dimensionless load P and the axial force parameter β at the limit points, which is expressed as Solving Equations (28) and (35) simultaneously can yield the axial force parameter β and the limit point instability load P. Then, by substituting β and P into Equation (26), the dimensionless radial displacement v, corresponding to the limit instability loads P can be obtained.

Bifurcation Buckling
The FG-GPLRC shallow arches with elastic rotational constraints may also buckle in a bifurcation mode. The critical condition for buckling may be stated that the second variation of the total potential of the arch and load system is equal to zero for any admissible infinitesimal deformation variations. Thus, the differential equation for bifurcation buckling is [22]. v iv The general solution of Equation (37) is where E 1 , E 2 , E 3 , and E 4 are the undetermined coefficients. The boundary conditions for deformations are The existence of a nontrivial solution to Equation (40) for E 1 -E 4 requires that their determinant of the coefficient matrix vanishes, which leads to the characteristic equation: When the first term of Equation (41) is equal to zero, there is Solving Equation (42) leads to Substituting Equation (43) into Equation (25), the axial compressive force N b of the FG-GPLRC shallow arch with elastic rotational constraints corresponding to an antisymmetric buckling is obtained as Substituting Equation (43) into Equation (28), the critical load P b of the FG-GPLRC shallow arch with elastic rotational constraints corresponding to an antisymmetric buckling is obtained, as The existence of real solutions for P b require that expressions in the square root of Equation (45) is non-negative, i.e., From Equations (29)- (31), it can be seen that, when β = η b π and the GPL distribution model is given, the Equation (45) is a function of λ and α. By solving Equation (46), the minimum geometric parameter λ b of the FG-GPLRC shallow arch, at which antisymmetric buckling is possible can be solved.
When the second factor of Equation (41) is equal to zero, there is From Equation (47), the solution for β can be obtained as which is a parameter related to the value of α. Substituting Equation (48) into Equation (25), the axial compressive force N s of the FG-GPLRC shallow arch with elastic rotational constraints is obtained as Substituting the solution given by Equation (47) into Equation (28) leads to By solving Equations (24), (49) and (50), the lowest buckling load qR can be obtained as which indicates that the normal axial compressive force qR during buckling is equal to the axial compressive force N s .

Comparison and Discussion
In order to meet the requirement of Donnell's shallow shell theory, in the following numerical FG-GPLRC shallow arch models, the width of section is chosen as b = 0.3 m, the total thickness is chosen as h = 0.25 m, and the axial length of the arch is chosen as S = 20 m. The geometric parameters of GPL are a GPL = 2.5 µm, b GPL = 1.5 µm and t GPL = 1.5 nm. Young's modulus of polymer matrix and GPL are 3 GPa and 1010 GPa, respectively. The density of the polymer matrix and GPL are 1200 kg/m 3 and 1062.5 kg/m 3 , respectively. The material properties of GPLs that are used in this paper are referred from Ref. [6].

Effect of the Number of Nanocomposite Layers on the Numerical Results
When a multilayer graphene-reinforced nanocomposite is used to simulate functionally graded properties, the number of layers of the nanocomposite (N L ) will significantly affect the simulation results. Hence, the elastically constrained FG-GPLRC shallow arch was simulated using GPLRCs with varying numbers of layers to examine the effect of the number of nanocomposite layers on the critical buckling load of the arch. The calculated results are summarized in Table 1. Note that the limit point buckling loads of FG-GPLRC shallow arches are obtained by solving Equations (28) and (35) simultaneously, while the bifurcation loads are obtained by solving Equation (45).
As demonstrated in Table 1, there is a relatively insignificant difference between the buckling loads of the elastically constrained functionally graded X-GPLRC arch, with N L = 20 and N L = 1000. This suggests that the layers number of the nanocomposite N L = 20 is sufficient to facilitate an accurate simulation of functionally graded property, whose thickness is continuously changing. Therefore, N L = 20 was used in the subsequent numerical simulation and analysis.
Additionally, the proposed analytical solutions of the stable equilibrium path and the buckling load of the elastically constrained FG-GPLRC shallow arch are validated by finite element (FE) solutions, which were obtained using the commercial FE software ANSYS. The FG-GPLRC shallow arch model was constructed using SHELL181 layered elements and the elastic constraints at the ends were established using COMBIN14 elements. The circular arch was meshed into 100 SHELL181 elements and six COMBIN14 elements. A comparison of the analytical and FE solutions is shown in Figure 3. solutions, which were obtained using the commercial FE software ANSYS. The FG-GPLRC shallow arch model was constructed using SHELL181 layered elements and the elastic constraints at the ends were established using COMBIN14 elements. The circular arch was meshed into 100 SHELL181 elements and six COMBIN14 elements. A comparison of the analytical and FE solutions is shown in Figure 3. In Figure 3, vc and f are the radial displacement and rise of the FG-GPLRC shallow arch. The analytical solutions are highly consistent with the FE simulation results, which indicates that the analytical method proposed in this study can be used to obtain a sufficiently accurate stable equilibrium path and critical buckling load of the elastically constrained FG-GPLRC shallow arch.

Limit Point Buckling Analysis
The limit point buckling loads of the elastically constrained arch when reinforced with GPLs in various distribution patterns are compared in Table 2. GPLs of various distribution patterns (U-, X-, O-, and A-types) significantly increase the critical buckling load of the elastically constrained pure epoxy resin arch, as demonstrated in Table 2. In addition, GPLs of the X-type distribution pattern increase the critical buckling load of the elastically constrained pure epoxy resin arch to the highest extent, followed by GPLs of the U-, A-, and O-type distribution patterns, in that order. This is because the elastically constrained composite arch reinforced with GPLs of the X-type distribution pattern contains a relatively large amount of reinforcing GPL filler in the top and bottom layers and is capable of resisting a relatively high bending stress; consequently, the critical buckling load of this type of reinforced composite arch is relatively high. GPLs with the O-type distribution pattern, the opposite of the X-type distribution pattern, increase the critical buckling load of the elastically constrained composite arch relatively insignificant. This is because most of the GPL fillers in the O-type distribution are concentrated near the neutral axis of the section, leading to a minor increase in cross sectional bending stiffness. In Figure 3, v c and f are the radial displacement and rise of the FG-GPLRC shallow arch. The analytical solutions are highly consistent with the FE simulation results, which indicates that the analytical method proposed in this study can be used to obtain a sufficiently accurate stable equilibrium path and critical buckling load of the elastically constrained FG-GPLRC shallow arch.

Limit Point Buckling Analysis
The limit point buckling loads of the elastically constrained arch when reinforced with GPLs in various distribution patterns are compared in Table 2. GPLs of various distribution patterns (U-, X-, O-, and A-types) significantly increase the critical buckling load of the elastically constrained pure epoxy resin arch, as demonstrated in Table 2. In addition, GPLs of the X-type distribution pattern increase the critical buckling load of the elastically constrained pure epoxy resin arch to the highest extent, followed by GPLs of the U-, A-, and O-type distribution patterns, in that order. This is because the elastically constrained composite arch reinforced with GPLs of the X-type distribution pattern contains a relatively large amount of reinforcing GPL filler in the top and bottom layers and is capable of resisting a relatively high bending stress; consequently, the critical buckling load of this type of reinforced composite arch is relatively high. GPLs with the O-type distribution pattern, the opposite of the X-type distribution pattern, increase the critical buckling load of the elastically constrained composite arch relatively insignificant. This is because most of the GPL fillers in the O-type distribution are concentrated near the neutral axis of the section, leading to a minor increase in cross sectional bending stiffness.  Figure 4 shows the nonlinear equilibrium paths and buckling loads for the elastically constrained FG-GPLRC shallow arch with various mass fractions of GPLs in the X-type distribution pattern (W GPL = 0.0% signifies the pure epoxy resin arch). As shown in Figure 4, as the mass fraction of GPLs increases, the critical load for the elastically constrained FG-GPLRC shallow arch to undergo limit point buckling increases, and the stability of the elastically constrained FG-GPLRC shallow arch increases. In addition, when compared to the elastically constrained pure epoxy resin arch, the enhancing effects of GPLs increase significantly as the GPL content increases.   Figure 4 shows the nonlinear equilibrium paths and buckling loads for the elastically constrained FG-GPLRC shallow arch with various mass fractions of GPLs in the X-type distribution pattern (WGPL = 0.0% signifies the pure epoxy resin arch). As shown in Figure 4, as the mass fraction of GPLs increases, the critical load for the elastically constrained FG-GPLRC shallow arch to undergo limit point buckling increases, and the stability of the elastically constrained FG-GPLRC shallow arch increases. In addition, when compared to the elastically constrained pure epoxy resin arch, the enhancing effects of GPLs increase significantly as the GPL content increases.   Figure 5 shows the effects of the geometric dimensions of GPLs on the limit point buckling load of the elastically constrained X-type FG-GPLRC shallow arch when the GPL content is fixed. As demonstrated in Figure 5, when the height/width (aGPL/bGPL) ratio of GPLs is fixed, the buckling load for the elastically constrained X-type FG-GPLRC shallow arch increases as the width/thickness (bGPL/tGPL) ratio increases. When the bGPL/tGPL ratio is greater than 10 3 , the effect of the bGPL/tGPL ratio on the buckling load is no longer significant. When the bGPL/tGPL ratio is fixed, the buckling load increases as the aGPL/bGPL ratio increases. However, when the aGPL/bGPL ratio is greater than 4, the effect of the aGPL/bGPL ratio on the buckling load gradually decreases as the aGPL/bGPL ratio increases. Specifically, the greater the aGPL/bGPL and bGPL/tGPL ratios of GPLs, the larger the surface area of GPLs, and the smaller the thickness of GPLs. Therefore, within a certain range, the larger and thinner GPLs have more significant enhancing effect on the stable bearing capacity of the elastically constrained composite shallow arch.  Figure 5 shows the effects of the geometric dimensions of GPLs on the limit point buckling load of the elastically constrained X-type FG-GPLRC shallow arch when the GPL content is fixed. As demonstrated in Figure 5, when the height/width (a GPL /b GPL ) ratio of GPLs is fixed, the buckling load for the elastically constrained X-type FG-GPLRC shallow arch increases as the width/thickness (b GPL /t GPL ) ratio increases. When the b GPL /t GPL ratio is greater than 10 3 , the effect of the b GPL /t GPL ratio on the buckling load is no longer significant. When the b GPL /t GPL ratio is fixed, the buckling load increases as the a GPL /b GPL ratio increases. However, when the a GPL /b GPL ratio is greater than 4, the effect of the a GPL /b GPL ratio on the buckling load gradually decreases as the a GPL /b GPL ratio increases. Specifically, the greater the a GPL /b GPL and b GPL /t GPL ratios of GPLs, the larger the surface area of GPLs, and the smaller the thickness of GPLs. Therefore, within a certain range, the larger and thinner GPLs have more significant enhancing effect on the stable bearing capacity of the elastically constrained composite shallow arch. Materials 2018, 11, x FOR PEER REVIEW 11 of 15 Figure 5. Effect of GPLs geometric dimensions on the limit point buckling load. Figure 6a shows the limit point buckling loads for the elastically constrained X-type FG-GPLRC shallow arch under various constraint stiffness conditions. The geometric parameter λ of the elastically constrained FG-GPLRC shallow arch varies between 4 and 60. Specifically, the constraint stiffness parameters a = 0 and   signify a completely clamped FG-GPLRC shallow arch and a completely hinged FG-GPLRC shallow arch, respectively. As shown in Figure 6a, as α gradually increases from 0 to ∞, the FG-GPLRC shallow arch gradually transitions from being a completely clamped arch to being a completely hinged arch, and the buckling load of the arch gradually decreases. When α = 1, the critical buckling load of the elastically constrained FG-GPLRC shallow arch is very close to that of a completely hinged FG-GPLRC shallow arch, thus suggesting that a certain low stiffness can facilitate a relatively satisfactory simulation of a completely hinged FG-GPLRC shallow arch.
In addition, the limit point buckling load of the elastically constrained X-type FG-GPLRC shallow arch increases as the geometric parameter λ increases. When λ is relatively low, the elastically constrained FG-GPLRC shallow arch is equivalent to a curved beam that undergoes nonlinear bending under loading; it has no extreme points for limit point buckling, and thus will not undergo buckling failure, as shown in Figure 6b.

Bifurcation Buckling Analysis
According to Equation (46), when the elastically constrained FG-GPLRC shallow arch undergoes antisymmetric bifurcation buckling, its geometric parameter must meet the following condition:   Figure 6a shows the limit point buckling loads for the elastically constrained X-type FG-GPLRC shallow arch under various constraint stiffness conditions. The geometric parameter λ of the elastically constrained FG-GPLRC shallow arch varies between 4 and 60. Specifically, the constraint stiffness parameters α = 0 and α = ∞ signify a completely clamped FG-GPLRC shallow arch and a completely hinged FG-GPLRC shallow arch, respectively. As shown in Figure 6a, as α gradually increases from 0 to ∞, the FG-GPLRC shallow arch gradually transitions from being a completely clamped arch to being a completely hinged arch, and the buckling load of the arch gradually decreases. When α = 1, the critical buckling load of the elastically constrained FG-GPLRC shallow arch is very close to that of a completely hinged FG-GPLRC shallow arch, thus suggesting that a certain low stiffness can facilitate a relatively satisfactory simulation of a completely hinged FG-GPLRC shallow arch. Figure 5. Effect of GPLs geometric dimensions on the limit point buckling load. Figure 6a shows the limit point buckling loads for the elastically constrained X-type FG-GPLRC shallow arch under various constraint stiffness conditions. The geometric parameter λ of the elastically constrained FG-GPLRC shallow arch varies between 4 and 60. Specifically, the constraint stiffness parameters a = 0 and   signify a completely clamped FG-GPLRC shallow arch and a completely hinged FG-GPLRC shallow arch, respectively. As shown in Figure 6a, as α gradually increases from 0 to ∞, the FG-GPLRC shallow arch gradually transitions from being a completely clamped arch to being a completely hinged arch, and the buckling load of the arch gradually decreases. When α = 1, the critical buckling load of the elastically constrained FG-GPLRC shallow arch is very close to that of a completely hinged FG-GPLRC shallow arch, thus suggesting that a certain low stiffness can facilitate a relatively satisfactory simulation of a completely hinged FG-GPLRC shallow arch.
In addition, the limit point buckling load of the elastically constrained X-type FG-GPLRC shallow arch increases as the geometric parameter λ increases. When λ is relatively low, the elastically constrained FG-GPLRC shallow arch is equivalent to a curved beam that undergoes nonlinear bending under loading; it has no extreme points for limit point buckling, and thus will not undergo buckling failure, as shown in Figure 6b.

Bifurcation Buckling Analysis
According to Equation (46), when the elastically constrained FG-GPLRC shallow arch undergoes antisymmetric bifurcation buckling, its geometric parameter must meet the following condition:  In addition, the limit point buckling load of the elastically constrained X-type FG-GPLRC shallow arch increases as the geometric parameter λ increases. When λ is relatively low, the elastically constrained FG-GPLRC shallow arch is equivalent to a curved beam that undergoes nonlinear bending under loading; it has no extreme points for limit point buckling, and thus will not undergo buckling failure, as shown in Figure 6b.

Bifurcation Buckling Analysis
According to Equation (46), when the elastically constrained FG-GPLRC shallow arch undergoes antisymmetric bifurcation buckling, its geometric parameter must meet the following condition: The values of the critical geometric parameter λ b required for the arch to undergo bifurcation buckling of the FG-GPLRC arches under a α = 1.5 elastic constraint condition are summarized Table 3. Table 3. Critical geometrical parameter λ b (S/h = 80, a GPL /b GPL = 4, b GPL /t GPL = 10 3 , W GPL = 0.5%). As demonstrated in Table 3, when α is fixed, the value of λ b for the O-type FG-GPLRC shallow arch to undergo bifurcation buckling is the smallest, suggesting that the elastically constrained FG-GPLRC shallow arch with GPLs in the O-type distribution pattern is the most prone to bifurcation buckling, followed by the FG-GPLRC shallow arches with GPLs of the U-, A-, and X-type distribution patterns, in that order.

Patterns
shallow arch to undergo bifurcation buckling is the smallest, suggesting that the elastically constrained FG-GPLRC shallow arch with GPLs in the O-type distribution pattern is the most prone to bifurcation buckling, followed by the FG-GPLRC shallow arches with GPLs of the U-, A-, and Xtype distribution patterns, in that order.  Figure 7 shows the buckling behavior of the FG-GPLRC shallow arch with GPLs in various distribution patterns under various constraint stiffness conditions. As shown in Figure 7a, when 1.5   , the bifurcation buckling point of the X-type FG-GPLRC shallow arch occurs after its limit point buckling point, indicating that limit point buckling will occur before bifurcation buckling for this type of arch. Under the same condition, the bifurcation buckling point of the O-type FG-GPLRC shallow arch occurs before its limit point buckling point, indicating that bifurcation buckling will occur before limit point buckling for this type arch (Figure 7b). This suggests that the distribution pattern of GPLs can alter the buckling mode of elastically constrained FG-GPLRC shallow arches. In addition, the stable equilibrium path for the O-type FG-GPLRC shallow arch when   Figure 7 shows the buckling behavior of the FG-GPLRC shallow arch with GPLs in various distribution patterns under various constraint stiffness conditions. As shown in Figure 7a, when α = 1.5, the bifurcation buckling point of the X-type FG-GPLRC shallow arch occurs after its limit point buckling point, indicating that limit point buckling will occur before bifurcation buckling for this type of arch. Under the same condition, the bifurcation buckling point of the O-type FG-GPLRC shallow arch occurs before its limit point buckling point, indicating that bifurcation buckling will occur before limit point buckling for this type arch (Figure 7b). This suggests that the distribution pattern of GPLs can alter the buckling mode of elastically constrained FG-GPLRC shallow arches. In addition, the stable equilibrium path for the O-type FG-GPLRC shallow arch when α = 0.1 was analysed, as shown in Figure 7c. A comparison of Figure 7b,c shows that, under the same conditions, altering the constraint stiffness can result in a change in the buckling mode of the FG-GPLRC shallow arch. Figure 8 shows the critical bifurcation buckling loads for the FG-GPLRC shallow arch under various constraint stiffness conditions (i.e., various values of α). As demonstrated in Figure 8, as α gradually increases, the constraint at the ends of the FG-GPLRC shallow arch gradually transitions from being a clamp constraint to being a hinge constraint, and the critical bifurcation buckling load gradually decreases. In addition, for the FG-GPLRC shallow arch with GPLs of various distribution patterns, the critical bifurcation buckling load of the FG-GPLRC shallow arch with GPLs of the X-type distribution pattern is the largest, followed by those of the FG-GPLRC shallow arch with GPLs of the U-, A-, and O-type distribution patterns, in that order. This finding is consistent with the conclusion that is derived from Table 2.
Materials 2018, 11, x FOR PEER REVIEW 13 of 15 Figure 8 shows the critical bifurcation buckling loads for the FG-GPLRC shallow arch under various constraint stiffness conditions (i.e., various values of α). As demonstrated in Figure 8, as α gradually increases, the constraint at the ends of the FG-GPLRC shallow arch gradually transitions from being a clamp constraint to being a hinge constraint, and the critical bifurcation buckling load gradually decreases. In addition, for the FG-GPLRC shallow arch with GPLs of various distribution patterns, the critical bifurcation buckling load of the FG-GPLRC shallow arch with GPLs of the Xtype distribution pattern is the largest, followed by those of the FG-GPLRC shallow arch with GPLs of the U-, A-, and O-type distribution patterns, in that order. This finding is consistent with the conclusion that is derived from Table 2.

Conclusions
To investigate the buckling behavior of an elastically constrained FG-GPLRC shallow arch under uniform radial load, this study establishes the nonlinear stable equilibrium equation for this type of arch while using the Halpin-Tsai micromechanical model of graphene and the principle of virtual work. The equation is then solved to determine the critical load of the elastically constrained FG-GPLRC shallow arch to undergo symmetric limit point buckling and antisymmetric bifurcation buckling under uniform radial load. Next, the proposed method that is employed in this study is validated through a finite element analysis. A parametric analysis was also carried out using the proposed theoretical method, and the results show that GPLs in various distribution patterns (U-, X-, O-, and A-types) all increase to some extent the critical buckling load of the elastically constrained pure epoxy resin arch. The critical buckling load of the elastically FG-GPLRC shallow arch reinforced with GPLs of the X-type distribution pattern, is the highest, so this type arch is the most stable. When the total graphene content is fixed, the enhancing effects increase as the aGPL/bGPL and bGPL/tGPL ratios of graphene increase. When the aGPL/bGPL ratio is greater than 4 and the bGPL/tGPL ratio is greater than 10 3 , the enhancing effects of the geometric shape of graphene on the buckling load of the arch are no longer significant. As α gradually increases, the constraint at the ends of the FG-GPLRC shallow arch gradually transitions from a clamp constraint to a hinge constraint, and the critical buckling load decreases. Altering the distribution pattern of graphene or the stiffness of the constraint at the ends of the FG-GPLRC shallow arch can result in a change in the buckling mode of the arch. Functionally graded graphene-reinforced composites have advantageous material properties, and one can therefore expect that it will be more and more commonly applied in engineering as the development of material manufacturing technology. The method and results that are presented in this paper would

Conclusions
To investigate the buckling behavior of an elastically constrained FG-GPLRC shallow arch under uniform radial load, this study establishes the nonlinear stable equilibrium equation for this type of arch while using the Halpin-Tsai micromechanical model of graphene and the principle of virtual work. The equation is then solved to determine the critical load of the elastically constrained FG-GPLRC shallow arch to undergo symmetric limit point buckling and antisymmetric bifurcation buckling under uniform radial load. Next, the proposed method that is employed in this study is validated through a finite element analysis. A parametric analysis was also carried out using the proposed theoretical method, and the results show that GPLs in various distribution patterns (U-, X-, O-, and A-types) all increase to some extent the critical buckling load of the elastically constrained pure epoxy resin arch. The critical buckling load of the elastically FG-GPLRC shallow arch reinforced with GPLs of the X-type distribution pattern, is the highest, so this type arch is the most stable. When the total graphene content is fixed, the enhancing effects increase as the a GPL /b GPL and b GPL /t GPL ratios of graphene increase. When the a GPL /b GPL ratio is greater than 4 and the b GPL /t GPL ratio is greater than 10 3 , the enhancing effects of the geometric shape of graphene on the buckling load of the arch are no longer significant. As α gradually increases, the constraint at the ends of the FG-GPLRC shallow arch gradually transitions from a clamp constraint to a hinge constraint, and the critical buckling load decreases. Altering the distribution pattern of graphene or the stiffness of the constraint at the ends of the FG-GPLRC shallow arch can result in a change in the buckling mode of the arch. Functionally graded graphene-reinforced composites have advantageous material properties, and one can therefore expect that it will be more and more commonly applied in engineering as the development of material manufacturing technology. The method and results that are presented in this paper would be useful as a reference for the stability design of the FG-GPLRC shallow arches in the future.