Multiple Equilibria and Buckling of Functionally Graded Graphene Nanoplatelet-Reinforced Composite Arches with Pinned-Fixed End

: This paper presents an analytical study on the multiple equilibria and buckling of pinned-ﬁxed functionally graded graphene nanoplatelet-reinforced composite (FG-GPLRC) arches under central point load. It is assumed that graphene nanoplatelets (GPLs) in each GPLRC layer are uniformly distributed and randomly oriented with its concentration varying layer-wise along the thickness direction. The Halpin–Tsai micromechanics-based model is used to estimate the e ﬀ ective elastic modulus of GPLRC. The non-linear equilibrium path and buckling load of the pinned-ﬁxed FG-GPLRC arch are subsequently derived by employing the principle of virtual work. The e ﬀ ects of GPLs distribution, weight fraction, size and geometry on the buckling load are examined comprehensively. It is found that the buckling performances of FG-GPLRC arches can be signiﬁcantly improved by using GPLs as reinforcing nanoﬁllers. It is also found that the non-linear equilibrium path of the pinned-ﬁxed FG-GPLRC arch have multiple limit points and non-linear equilibrium branches when the arch is with a special geometric parameter. micromechanics-based equilibrium for equilibrium of


Introduction
The development of carbonaceous nanomaterials, such as graphene and carbon nanotubes, has been greatly advanced by rapidly developing material science and manufacturing techniques over the last decade [1][2][3][4]. Established research works have demonstrated that adding graphene and its derivatives into polymer matrix can form nanocomposites with excellent mechanical properties. Rafiee et al. [5,6] comprehensively examined the reinforcing effect of graphene nanoplatelets (GPLs) and carbon nanotubes on the mechanical properties of epoxy nanocomposites, and indicated the remarkable effect of graphene nanoplatelets. Liu et al. [7,8] studied the effects of GPLs on the microstructure and mechanical properties of Al 2 O 3 -based ceramic composites. The results shown that, by adding GPLs into Al 2 O 3 ceramic composites, the flexural strength and fracture toughness of the composite increased by 30.75% and 27.20%, respectively. Experimental results conducted by Tang et al. [9] showed that the epoxy nanocomposites with graphene reinforcements have better strength and fracture toughness when the graphene is highly dispersed in composites. Zhao et al. [10] reported that, at a graphene loading of 1.8 vol %, a 150% improvement of tensile strength and a nearly 10 times increase of Young's modulus of the graphene/poly (vinyl alcohol) (PVA) composites are obtained.
Except for the study on the material properties of the GPL-reinforced nanocomposite, the mechanical behaviors of structures made of these composites have been investigated widely. Yang and co-workers [11][12][13][14][15] carried out a series of studies on the vibration, dynamic stability, buckling and postbuckling of functionally graded graphene nanoplatelet-reinforced composite (FG-GPLRC) structures, which revealed that the performances of the FG-GPLRC structures can be remarkably improved by adding a low concentration of GPLs. Dong et al. [16][17][18] analytically investigated the non-linear dynamic behavior of graded graphene-reinforced composite cylindrical shells subjected to thermal loading, impact loading and harmonic loading. Shen et al. [19][20][21] studied the non-linear dynamic responses and instability of laminated plates consisting of graphene-reinforced composite (GRC) layers in thermal environments where the materials properties of GRC are estimated by an extended Halpin-Tsai model. Gao et al. [22] presented an investigation of the probabilistic stability characteristics of FG-GPLRC beams by taking into account the multidimensional probability distributions. Wang et al. [23][24][25] theoretically studied the buckling, postbuckling, nonlinear free vibration, nonlinear static and dynamic responses characteristics of GPL-reinforced Polyvinylidene fluoride (PVDF) composite beam with dielectric permittivity.
For the arch structure, the authors had performed in-depth studies on the nonlinear static and dynamic behaviors of FG-GPLRC arches with different boundary conditions, such as fixed ends, pinned ends, and elastic end restraints [26][27][28][29][30][31]. Liu et al. [32] presented an analytical solution for non-linear static responses and buckling of functionally graded porous (FGP) arches with GPLs reinforcements. Li et al. [33] studied the coupling effects of temperature and mechanical load on the stability of confined FGP arches. Up to now, no previous work has been carried out on the pinned-fixed FG-GPLRC arch for investigating its multiple equilibria and non-linear buckling behaviors.
Therefore, this paper aims to investigate the multiple equilibria and nonlinear buckling behaviors of a pinned-fixed FG-GPLRC arch. The effective elastic modulus of GPLRC is estimated by Halpin-Tsai micromechanics-based model. The principle of virtual work is used to establish the non-linear equilibrium equations of the pinned-fixed FG-GPLRC arch. The analytical solutions for the non-linear equilibrium path and buckling load are derived. The effects of GPLs distribution, weight fraction, size and geometry on the limit point buckling load are examined comprehensively.

Effective Elastic Modulus
Before the structural analysis, the effective elastic modulus of the composite needs to be determined. As follows, the Halpin-Tsai micromechanics-based model [34] is employed to calculate of the effective elastic modulus of GPLRC as: with, where E GPL and E m are the Young's modulus of GPLs and polymer matrix, respectively. V GPL is the volume fraction of GPLs, and a GPL , b GPL , t GPL , a GPL /b GPL , b GPL /t GPL are the length, width, thickness, aspect ratio, and width-to-thickness ratio of GPLs, respectively. As shown in Figure 1, three GPL distribution patterns are studied in this paper, namely U-GPLRC, X-GPLRC and O-GPLRC, where darker color represents a higher content of GPLs. In particular, the GPLs concentration in U-GPLRC is constant which represents an isotropic homogeneous arch. Crystals 2020, 10, x FOR PEER REVIEW 3 of 14 particular, the GPLs concentration in U-GPLRC is constant which represents an isotropic homogeneous arch. The GPL volume fraction of the kth layer for the above three GPL distribution patterns are formulated as: In which, NL is total number of layers, * GPL V and WGPL are the total volume fraction and weight fraction of GPLs, respectively.

Equilibrium Equations
The FG-GPLRC arch, shown in Figure 2, with pinned-fixed end, arc length S, central angle 2Θ , cross section b × h, subjected to central point load Q is investigated.  The GPL volume fraction of the kth layer for the above three GPL distribution patterns are formulated as: with, In which, N L is total number of layers, V * GPL and W GPL are the total volume fraction and weight fraction of GPLs, respectively.

Equilibrium Equations
The FG-GPLRC arch, shown in Figure 2, with pinned-fixed end, arc length S, central angle 2Θ, cross section b × h, subjected to central point load Q is investigated. The GPL volume fraction of the kth layer for the above three GPL distribution patterns are formulated as: In which, NL is total number of layers, * GPL V and WGPL are the total volume fraction and weight fraction of GPLs, respectively.

Equilibrium Equations
The FG-GPLRC arch, shown in Figure 2, with pinned-fixed end, arc length S, central angle 2Θ , cross section b × h, subjected to central point load Q is investigated.  According to Donnell's shallow shell theory, the normal strain ε in the arch is given as [35]: The non-linear equilibrium equation of FG-GPLRC arches is established by employing the principle of virtual work as: in which, δ D (θ) is the Dirac-delta function defined by: with the property: Substituting Equation (8) into Equation (9) leads to the virtual work statement as: with axial compressive force N and bending moment M as: with, where κ is the effective bending stiffness. Integrating Equation (13) by parts leads to the differential equations of equilibrium as: and the boundary conditions as: Substituting Equations (15) and (17) into Equation (18), the equilibrium equation in radial direction of the arch can be derived as: Crystals 2020, 10, 1003

of 13
Solving Equation (20) by incorporating the boundary conditions given by Equation (19), the explicit expression of radial displacement for FG-GPLRC arch with pinned-fixed boundary condition is determined as: with, and with, and H(θ) is a step function as: By substituting Equation (21) into Equation (14) and integrating along the entire FG-GPLRC arch leads to the equilibrium equation as: with, and where λ = RΘ 2 /h is the dimensionless geometrical parameter of the FG-GPLRC arch. Accordingly, by solving Equations (21) and (25), the solutions of the nonlinear equilibrium path of the FG-GPLRC arch with pinned-fixed boundary condition can be determined.

Limit Point Buckling
When the applied point load reaches a maximum load, namely limit point load, the limit point buckling may occur to the FG-GPLRC arch, and the limit point buckling load is the local maximum point on the non-linear equilibrium path defined in Equation (25). Therefore, by setting Equation (25) as Crystals 2020, 10, 1003 6 of 13 an implicit function F(Q, β) = 0 and performing dQ/dβ, the extrema equation for limit point buckling load is derived as: In response to this, the limit point buckling load can be obtained by solving Equations (25) and (30) simultaneously.

Multiple Equilibria and Inflection Point
According to the previous research [26], when the arch buckle is in a bifurcation mode, the arch is in an equilibrium state v, w, N, M, Q associated with the pre-buckling condition and also in an infinitesimally closed bifurcation equilibrium state as: Substituting Equation (31) into Equation (18) and neglecting the second-and higher-order terms, the differential equation of the bifurcation equilibrium state can be formulated as: If a pinned-fixed arch can buckle in a bifurcation buckling mode, the bifurcation equilibrium state N b associated with Equation (14) should become zero, and Equation (32) is simplified as: However, as pointed out by Pi and Bradford [36], the arch with pinned-fixed end cannot buckle in a bifurcation mode due to its asymmetric deformation, and Equation (33) is false. Even so, Equation (33) may hold at an inflection point for a specific arch with a certain value of the arch geometry parameter λ sb1 (or λ sb2 ), and the value of λ sb1 (or λ sb2 ) can be derived from solving Equation (33) as follows.
Solving Equation (33), its general solution is obtained as: where the coefficients D 1 , D 2 , D 3 , and D 4 can be determined from specified boundary conditions of the pinned-fixed FG-GPLRC arch.
Substituting the boundary conditions of the pinned-fixed FG-GPLRC arch into Equation (34), leads to the following four linear homogeneous algebraic equations: The existence of a non-trivial solutions to Equation (35) for coefficients D 1 -D 4 requires the vanishing of their determinant of the coefficient matrix, which leads to the characteristic equation as: The first two solutions of Equation (36) are: From the definitions of µ and β, the axial compressive forces corresponding to β 1 and β 2 are: The dimensionless load P in inflexion point can be solved from Equation (25) in β = β 1 or β 2 as: It should be mentioned that the existence of real solutions for Equation (39) requires that: By substituting β 1 or β 2 into Equation (39), respectively, a special geometric parameter λ sb1 (or λ sb2 ) making Equation (33) true can be solved, which is also defined as the minimum geometric parameter for distinguishing the number of limit points.

Parameters Studies
In order to verify the proposed results, the commercial finite element (FE) package ANSYS 17.0 (ANSYS, Inc., Cannonsburg, PA, USA) were used to perform nonlinear analysis for the nonlinear buckling of a pinned-fixed FG-GPLRC arch. Table 1 lists the material properties of the FG-GPLRC arch for the following parameters studies. In ANSYS simulation, the shell element SHELL181 with six degrees of freedom (DOF) at each node is used to develop the FE model of the FG-GPLRC arch, and the arch is meshed into 100 elements for non-linear buckling analysis. The layered composite is defined by using the section commands to specify the thickness, material and orientation of each layer. In addition, the equilibrium path and structural displacement of pinned-fixed FG-GPLRC arches under the external load are traced by using the arc-length method in the FE analysis.
The finite element result and present result for the non-linear equilibrium path of the FG-GPLRC arch are shown together in Figure 3 for comparison, where N E0 is the N E2 given by Equation (37) for the pure epoxy arch, and v c /f is the dimensionless displacement of the arch crown with f being the rise of the arch. It is found that the present result is almost identical to the finite element result. Figure 4 shows the multiple equilibria and buckling of the pinned-fixed FG-X-GPLRC arch in the form of load-displacement curve and load-axial force curve. It can be seen that, when the pinned-fixed FG-X-GPLRC arch has λ ≤ λ sb1 = 2.1672 (Figure 4a,b), there is no limit point on its non-linear equilibrium path, which means that the pinned-fixed FG-GPLRC arch cannot buckle under the central point load. When the arch has λ sb1 < λ ≤ λ sb2 = 4.9504 (Figure 4c,d), there are two limit points (b and c) and three branches (ab, bc and cd) on the non-linear equilibrium paths. When the arch has λ = 8 > λ sb2 (Figure 4e,f), there are four limit points (b, c, d and e) and five branches (ab, bc, cd, de and ef ) on the non-linear equilibrium paths. Hence, it can be concluded that the special arch geometric parameter λ sb1 and λ sb2 can be used as a switch between arches with two or four limit points and those without limit point. It is worth pointing out that no matter how the λ increases, if λ > λ sb2 , there are four and only four limit points on the non-linear equilibrium path. It also can be observed from Figure 4c-f that the curve on both sides of the upper limit point is concave downwards with the slopes changing from positive to negative, whereas the curve on both sides of the lower limit point is concave upwards with the slopes changing from negative to positive. The deformed shape of the arch under limit point load has a greater change in curvature.

GPL length
2.5 µm GPL width 1.5 µm GPL thickness 1.5 nm Young's modulus (epoxy) 3 GPa Young's modulus (GPL) 1010 GPa In ANSYS simulation, the shell element SHELL181 with six degrees of freedom (DOF) at each node is used to develop the FE model of the FG-GPLRC arch, and the arch is meshed into 100 elements for non-linear buckling analysis. The layered composite is defined by using the section commands to specify the thickness, material and orientation of each layer. In addition, the equilibrium path and structural displacement of pinned-fixed FG-GPLRC arches under the external load are traced by using the arc-length method in the FE analysis.
The finite element result and present result for the non-linear equilibrium path of the FG-GPLRC arch are shown together in Figure 3 for comparison, where NE0 is the NE2 given by Equation (37) for the pure epoxy arch, and vc/f is the dimensionless displacement of the arch crown with f being the rise of the arch. It is found that the present result is almost identical to the finite element result.      Taking the first upper limit point as the limit point buckling load, the limit point buckling loads of FG-X-GPLRC arches with different geometry are listed in Table 2. It is found that the FG-X-GPLRC arch has the highest limit point buckling load, and the cracked FG-X-GPLRC beam has the smallest one. This is because the FG-X-GPLRC arch with more GPLs filler dispersed in the top and bottom layers possesses a higher bending stiffness and consequently has a higher buckling load. In addition,   Taking the first upper limit point as the limit point buckling load, the limit point buckling loads of FG-X-GPLRC arches with different geometry are listed in Table 2. It is found that the FG-X-GPLRC arch has the highest limit point buckling load, and the cracked FG-X-GPLRC beam has the smallest one. This is because the FG-X-GPLRC arch with more GPLs filler dispersed in the top and bottom layers possesses a higher bending stiffness and consequently has a higher buckling load. In addition, Figure 6 plots the buckling load of FG-GPLRC arches with varying weight fraction WGPL, and shows that the buckling load of FG-GPLRC arches are much higher than that of pure epoxy arch (WGPL = 0.0%) in all cases, which indicates the remarkable reinforcing effect of GPLs. Taking the first upper limit point as the limit point buckling load, the limit point buckling loads of FG-X-GPLRC arches with different geometry are listed in Table 2. It is found that the FG-X-GPLRC arch has the highest limit point buckling load, and the cracked FG-X-GPLRC beam has the smallest one. This is because the FG-X-GPLRC arch with more GPLs filler dispersed in the top and bottom layers possesses a higher bending stiffness and consequently has a higher buckling load. In addition, Figure 6 plots the buckling load of FG-GPLRC arches with varying weight fraction W GPL , and shows that the buckling load of FG-GPLRC arches are much higher than that of pure epoxy arch (W GPL = 0.0%) in all cases, which indicates the remarkable reinforcing effect of GPLs.     Figure 7 presents the influences of GPL geometry on the limit point buckling load of the FG-X-GPLRC arch with the GPL width-to-thickness ratio b GPL /t GPL increasing from 10 to 5000. According to the Equations (1)-(3), GPL with higher value of width-to-thickness ratio and aspect ratio have a more significant enhancing effect on the effective elastic modulus of GPLRC. Therefore, the buckling load increases when GPL width-to-thickness ratio or aspect ratio increases, but the effect of GPL geometry tends to be insignificant when b GPL /t GPL > 1000 and b GPL /t GPL > 4, respectively. Figure 8 investigates the effect of the arch's central angle on limit point buckling loads of the FG-X-GPLRC arch where the arch's central angle varies from 10 • to 120 • . It can be seen that the pinned-fixed FG-X-GPLRC arch with a bigger central angle has higher limit point buckling load.
GPLRC arch with the GPL width-to-thickness ratio bGPL/tGPL increasing from 10 to 5000. According to the Equations (1)-(3), GPL with higher value of width-to-thickness ratio and aspect ratio have a more significant enhancing effect on the effective elastic modulus of GPLRC. Therefore, the buckling load increases when GPL width-to-thickness ratio or aspect ratio increases, but the effect of GPL geometry tends to be insignificant when bGPL/tGPL > 1000 and bGPL/tGPL > 4, respectively.

Conclusions
The multiple equilibria and buckling of the pinned-fixed FG-GPLRC arch under central point load has been comprehensively studied in this paper. The analysis demonstrates that there exists a special geometric parameter for distinguishing the number of multiple limit points and non-linear equilibrium branches of the FG-GPLRC arch. The analytical solutions for the special geometric parameter, limit point buckling load, structural displacement and axial force are derived. According to numerical studies, it is found that the pinned-fixed FG-GPLRC arch under central point load can buckle in unsymmetrical limit point instability mode, and there are four and only four limit points on the non-linear equilibrium path for the arch possessing a geometric parameter higher than a special value. It is also found that distributing more GPLs at the surface of the arch can considerably increase the buckling load of the FG-GPLRC arch.

Conclusions
The multiple equilibria and buckling of the pinned-fixed FG-GPLRC arch under central point load has been comprehensively studied in this paper. The analysis demonstrates that there exists a special geometric parameter for distinguishing the number of multiple limit points and non-linear equilibrium branches of the FG-GPLRC arch. The analytical solutions for the special geometric parameter, limit point buckling load, structural displacement and axial force are derived. According to numerical studies, it is found that the pinned-fixed FG-GPLRC arch under central point load can buckle in unsymmetrical limit point instability mode, and there are four and only four limit points on the non-linear equilibrium path for the arch possessing a geometric parameter higher than a special value. It is also found that distributing more GPLs at the surface of the arch can considerably increase the buckling load of the FG-GPLRC arch.
Author Contributions: All authors contributed to the main idea of this paper. Z.Y., H.L. and A.L. carried out the theoretical derivation. Z.Y., J.X., J.L. and J.F. analysed the numerical results. Z.Y. and H.L. wrote the article. All authors have read and agreed to the published version of the manuscript.