Analytical Prediction for Nonlinear Buckling of Elastically Supported FG-GPLRC Arches under a Central Point Load

In this paper, we present an analytical prediction for nonlinear buckling of elastically supported functionally graded graphene platelet reinforced composite (FG-GPLRC) arches with asymmetrically distributed graphene platelets (GPLs). The effective material properties of the FG-GPLRC arch are formulated by the modified Halpin–Tsai micromechanical model. By using the principle of virtual work, analytical solutions are derived for the limit point buckling and bifurcation buckling of the FG-GPLRC arch subjected to a central point load (CPL). Subsequently, the buckling mode switching phenomenon of the FG-GPLRC arch is presented and discussed. We found that the buckling modes of the FG-GPLRC arch are governed by the GPL distribution pattern, rotational restraint stiffness, and arch geometry. In addition, the number of limit points in the nonlinear equilibrium path of the FG-GPLRC arch under a CPL can be determined according to the bounds of successive inflexion points. The effects of GPL distribution patterns, weight fractions, and geometric configurations on the nonlinear buckling behavior of elastically supported FG-GPLRC arches are also comprehensively discussed.


Introduction
Functionally graded material (FGM) structures, characterized by a continuous change in the material compositions along one or multiple directions, have attracted extensive attention from both research and industrial communities owing to their excellent stiffness and strength-to-weight properties as compared with homogeneous composite structures [1][2][3]. To better understand the performance of FGM structures in practical engineering applications, researchers have conducted a series of investigations on the structural behavior of FGM structures. For instance, Ke et al. [4] presented analytical solutions for the nonlinear vibration responses of FGM beams with different end supports and discussed the influence of bending-stretching coupling on the nonlinear vibration of FGM beams. Librescu et al. [5] analyzed the vibration and stability of FGM thin-walled beams under a high-temperature environment. Yan et al. reported analytical solutions for the dynamic instability [6] and dynamic responses [7] of FGM beams with open-edge cracks. They investigated the crack depth and location effects on the mechanical behavior of such beams. Nguyen et al. [8] studied the mechanical buckling of stiffened FGM plates by using the finite element method. They found that the addition of stiffeners to the FGM plate could significantly reduce the weight of the FGM plate. Chen et al. [9] investigated the imperfection sensitivity in the nonlinear vibration of initially stressed FGM plates. According to their results, the effects of the initial stress, geometric imperfection, and volume fraction index were quite significant on the nonlinear vibration behavior of the FGM plate. Hao et al. [10] employed an asymptotic perturbation method to analyze the nonlinear oscillations, bifurcations, and chaotic motions of FGM plates. More relevant studies to investigate the significant performance of FGM structures could be found from the open literature [11][12][13][14][15][16][17]. Graphene, as an emerging high-performance nanofillers, has attracted considerable attention in aerospace, mechanical, thermal, and electrical engineering fields. It has been demonstrated by researchers that the reinforced performance of graphene nanoplatelets (GPLs) is significantly superior to other reinforcement materials. Adding a low concentration of GPLs can improve the stiffness and strength of reinforced composites significantly [18][19][20][21][22][23]. By introducing GPLs to FGM materials, novel FG GPLs-reinforced composite (FG-GPLRC) structures have been developed recently and have since attracted extensive attention in both research and engineering communities [24]. Yang and his co-workers conducted pioneering studies on the mechanical behaviors of FG-GPLRC structures, such as beams [25][26][27][28], shells [29,30], and plates [31,32]. Focusing on the stability analysis of FG-GPLRC arches, the authors devoted extensive efforts to the investigation of such structures, including the characteristics of nonlinear static buckling, dynamic buckling, and free vibration for FG-GPLRC arches with different boundary conditions and external loads [33][34][35][36][37][38][39]. In addition, Liu et al. [40] analyzed the nonlinear behavior and stability of functionally graded porous (FGP) arches reinforced by GPLs and obtained the critical buckling load under a uniform load. Zhao et al. [41] discussed the linear buckling, fundamental frequency, and dynamic instability of porous arches using an analytical method. Li et al. [42] further investigated the mechanics of the confined porous arches by including the thermal effect. Hitherto, to the best of our knowledge, far too little attention has been paid to the nonlinear buckling behavior of elastically supported FG-GPLRC arches with asymmetric distributed GPLs under a central point load (CPL).
Therefore, to fill this research gap, the nonlinear buckling behavior of elastically supported FG-GPLRC arches under a CPL is investigated in this study. The effective material properties of the FG-GPLRC arch are formulated by the modified Halpin-Tsai micromechanical model, because the Young's modulus of the FG-GPLRC arch predicted by this model agreed well with the experimental results [19]. Subsequently, the principle of virtual work is employed to derive the nonlinear buckling load for the limit point buckling and bifurcation buckling modes, from which the critical geometric parameters to identify the buckling load and the number of limit points can also be determined. In addition, the influences of GPL distribution patterns, GPL weight fractions, geometric configurations, as well as the flexibility of the rotational constraints on the nonlinear buckling behavior of the arch are discussed in detail. The bifurcation stability criteria presented in this work can provide essential information for the structural/material design of FG-GPLRC arches that have great potential in various engineering applications, for example, arch-shaped micro-electromechanical devices as sensors and transducers [43], and dielectric elastomer actuators as lightweight speakers [44,45]. In addition, the presence of analytical solutions is useful to engineers and researchers for benchmarking the convergence and validity of numerical methods for arch buckling analysis.

Effective Material Properties
In Figure 1, we consider an elastically restrained FG-GPLRC arch with multiple layers, N L , under a CPL Q. The rectangular cross section of the arch is b × h (width × thickness). It is assumed that the reinforcement of GPLs is uniformly distributed in the isotropic polymer matrix, so each individual GPLRC layer can be regarded as an isotropic and homogeneous material. The variation of GPLs is continuous along the arch thickness in accordance with the power law distribution. The GPL volume fraction of the k th GPLRC layer GPL k V is formulated as follows: where NL is the total number of the GPLRC layer, and n is the power law index that characterizes the distribution of GPLs. When n = 0, it corresponds to a uniform distribution of GPLs reinforcements in the thickness direction. * GPL V is the total GPL volume fraction which can be determined using: where WGPL is the GPLs weight fraction and GPL ρ and m ρ are the mass densities of GPLs and matrix, respectively. According to the modified Halpin-Tsai micromechanical model [43], the effective Young's modulus E k of the k th GPLRC layer is given by where EGPL and Em are the Young's moduli of GPLs and matrix, respectively. aGPL, bGPL, and tGPL are the length, width, and thickness of GPLs, respectively.
The Poisson's ratio k υ of each GPLRC layer is determined by the rule of mixture as follows: It is assumed that the reinforcement of GPLs is uniformly distributed in the isotropic polymer matrix, so each individual GPLRC layer can be regarded as an isotropic and homogeneous material. The variation of GPLs is continuous along the arch thickness in accordance with the power law distribution. The GPL volume fraction of the kth GPLRC layer V k GPL is formulated as follows: where N L is the total number of the GPLRC layer, and n is the power law index that characterizes the distribution of GPLs. When n = 0, it corresponds to a uniform distribution of GPLs reinforcements in the thickness direction. V * GPL is the total GPL volume fraction which can be determined using: where W GPL is the GPLs weight fraction and ρ GPL and ρ m are the mass densities of GPLs and matrix, respectively. According to the modified Halpin-Tsai micromechanical model [43], the effective Young's modulus E k of the kth GPLRC layer is given by with where E GPL and E m are the Young's moduli of GPLs and matrix, respectively. a GPL , b GPL , and t GPL are the length, width, and thickness of GPLs, respectively. The Poisson's ratio υ k of each GPLRC layer is determined by the rule of mixture as follows: where υ GPL and υ m are the Poisson's ratios of GPLs and matrix, respectively. Consider GPLs with the material properties E GPL = 1010 GPa, ρ GPL = 1062.5 kg/m 3 , υ GPL = 0.186, a GPL = 2.5 µm, b GPL = 1.5 µm, and t GPL = 1.5 nm, and epoxy (polymer matrix) with E m = 3 GPa, υ m = 0.34, and ρ m = 1200 kg/m 3 [25], the effective Young's moduli of the FG-GPLRC arch with 10 GPLRC layers under various GPLs distributions (different power law indices) are shown in Figure 2.
where GPL υ and m υ are the Poisson's ratios of GPLs and matrix, respectively.

Mathematical Modeling
As shown in Figure 1, an elastically supported FG-GPLRC arch with an angle 2Θ, a radius R, and an arc length S under the effect of a CPL is studied in this section. The stiffness of the elastic rotational constraints at both ends is k, the radial and axial displacement of the arch are v and w, respectively, and θ is the angular coordinate. By using the Donnell's shallow shell theory [39], the nonlinear strain-displacement relations for the FG-GPLRC arch are adopted as follows: According to the principle of virtual work, the governing equation is established as: is the Dirac delta function [27].

Mathematical Modeling
As shown in Figure 1, an elastically supported FG-GPLRC arch with an angle 2Θ, a radius R, and an arc length S under the effect of a CPL is studied in this section. The stiffness of the elastic rotational constraints at both ends is k, the radial and axial displacement of the arch are v and w, respectively, and θ is the angular coordinate. By using the Donnell's shallow shell theory [39], the nonlinear strain-displacement relations for the FG-GPLRC arch are adopted as follows: where ( ) = d( )/dθ. According to the principle of virtual work, the governing equation is established as: where δ D (θ) is the Dirac delta function [27]. Substituting Equation (6) into Equation (7), one obtains with the axial force As the present model is derived in accordance with the classical Euler-Bernoulli theory, so that A 11 , B 11 and D 11 are E k -based constants.
From Equation (8), the boundary conditions for the elastically supported FG-GPLRC arch are given by Substituting Equations (9) and (10) into Equation (8), obtains the equilibrium equations as follows: (NR) = 0 (13) Inserting Equations (9), (10) and (13) into Equation (14) yields the radial equilibrium equation as follows: v iv where κ is effective bending stiffness of the FG-GPLRC arch. Solving Equation (11) and considering the boundary conditions of the elastically support of the arch, the analytical solution of the dimensionless radial displacement can be determined as follows: where β = µΘ, B = B 11 /A 11 , and α is the flexibility of the rotational constraints defined as: The dimensionless CPL is defined as: and the parameters K 1 , K 2 , K 3 , and K 12 are given as: 2αβ cos β+sin β 2αβ cos β+sin β , K 12 = 2αβ 2αβ cos β+sin β (19) In Equation (16), H(θ) is a step function defined by [33] as: Substituting Equation (16) into Equation (9), obtains the nonlinear equilibrium relationship as: where A 1 , B 1 , and C 1 are determined as: 11 , and λ is the arch geometrical parameter as The parameters Ξ 1 , Ξ 2 , and Ξ 3 are given as

Limit Point Buckling
When the FG-GPLRC arch buckles in a limit point instability mode, the limit point loads may be either local maxima or local minima on the nonlinear equilibrium path, which are derived from Equation (21) as: Accordingly, the solutions of the limit point load and nonlinear equilibrium path of FG-GPLRC arches under a CPL can be determined by solving Equations (21) and (27).

Bifurcation Buckling
When the FG-GPLRC arch buckles in a bifurcation mode, the equilibrium equation for the arch can be derived by substituting the critical states The general solution of Equation (29) is solved as: where E 1 , E 2 , E 3 , and E 4 are unknown coefficients. Substituting the boundary conditions of the arch into Equation (30) yields the following characteristic equation for the coefficients E 1 , E 2 , E 3 , and E 4 as: From Equation (31), for the first case, it is true that the first term becomes zero, from which the critical axial force parameter β b for the bifurcation buckling is solved as: and the corresponding axial force is obtained as: Substituting Equation (32) into Equation (21), the bifurcation buckling equation is determined as: where from which the bifurcation load is solved as: The existence of real solutions in Equation (36) requires B 2 3 − 4A 3 C 3 ≥ 0, this condition gives a critical geometric parameter λ b1 that can trigger the bifurcation buckling of the arch.
Consider the second case of Equation (31) (i.e., the second term is zero), we may obtain the critical axial force parameter β sn as follows: β sn = η sn π (sn = s1, s2, s3, . . . , ) Following similar procedures, a critical geometric parameter λ sn that governs the number of inflection points in the nonlinear equilibrium path of the FG-GPLRC arch can also be determined. It is worth noting that the CPL corresponding to the first (as well as the lowest) axial force parameter β s1 is also the lowest buckling load of the arch whose geometric parameter is λ s1 . When the FG-GPLRC arch has λ ≤ λ s1 , the arch becomes a slightly curved beam and does not perform typical nonlinear buckling behavior.

Numerical Studies and Discussion
In this section, the nonlinear buckling behavior of the FG-GPLRC arch with a cross section of b × h = 0.03 m × 0.025 m and 10 GPLRC layers is investigated. The material properties of GPLs and matrix are adopted as the same in Section 2. To verify the present solutions, a finite element (FE) analysis conducted by ANSYS 17.0 is used to obtain the nonlinear equilibrium path of the FG-GPLRC arch. In numerical modeling, the SHELL181 element is used to establish a multi-layer FG-GPLRC arch model. The COMBIN14 element is adopted to model elastic end restraints whose rotational restraint stiffness k is defined by Equation (17). The section commands are used to define the layered structures which can provide the input options for specifying the thickness, material, and orientation of each layer. The FG-GPLRC arch model is meshed by 100 elements along the arch length direction, and 10 layers are considered to be in the thickness direction, as depicted in Figure 3.
As a large deflection and rotation may occur during the whole deformation of the arch, the NLGEOM command is used for nonlinear analysis when considering the effect of geometric nonlinearities. In addition, the arc-length method activated by the ARCLEN command is, then, used to trace the nonlinear equilibrium path of the FG-GPLRC arch with a load step of 200, as specified by the NSUBST command. By setting the minimum and maximum multipliers for the arc-length radius, the subsequent displacement and load proportional factors and the increment size can be adjusted and computed automatically using the arc-length radius.
For comparison, the FE and present results for the limit point buckling path and bifurcation buckling path are plotted in Figure 4, in which N E0 = the critical bifurcation axial force N E2 of the pure epoxy arch, and v c /f = the dimensionless displacement of the arch crown with f being the rise of the arch. It should be mention that, for triggering the bifurcation buckling behavior of the FG-GPLRC arch, an antisymmetric geometric imperfection of 0.1% arch length is introduced to the arch in the bifurcation buckling analysis. It is observed that the present analytical solutions agree very well with the FE results.   Next, we discuss the buckling load of the elastically supported arch, i.e., the limit point buckling load (the first upper limit point load) for brevity, unless stated otherwise. Figure 5 shows the influence of the power law index on the buckling load and the effective bending stiffness of the arch for 2 α = . In Figure 5b, "D110" is the bending stiffness of the pure epoxy arch. It is found that the buckling load and the effective bending stiffness significantly decrease as the power law index increases, but the effect of the power law index tends to be less pronounced when the power law index is higher than three. This is be-    Next, we discuss the buckling load of the elastically supported arch, i.e., the limit point buckling load (the first upper limit point load) for brevity, unless stated otherwise. Figure 5 shows the influence of the power law index on the buckling load and the effective bending stiffness of the arch for 2 α = . In Figure 5b, "D110" is the bending stiffness of the pure epoxy arch. It is found that the buckling load and the effective bending stiffness significantly decrease as the power law index increases, but the effect of the power law index tends to be less pronounced when the power law index is higher than three. This is be-  Next, we discuss the buckling load of the elastically supported arch, i.e., the limit point buckling load (the first upper limit point load) for brevity, unless stated otherwise. Figure 5 shows the influence of the power law index on the buckling load and the effective bending stiffness of the arch for α = 2. In Figure 5b, "D 110 " is the bending stiffness of the pure epoxy arch. It is found that the buckling load and the effective bending stiffness significantly decrease as the power law index increases, but the effect of the power law index tends to be less pronounced when the power law index is higher than three. This is because a higher content of GPLs would be concentrated on the bottom and less on other layers when the FG-GPLRC arch has a higher power law index, thereby leading to a lower bending stiffness of the arch.  Table 1 lists the buckling load of the arch having different GPL weight fractions. It is found that, by comparing with a pure epoxy arch (i.e., WGPL = 0.0%), GPLs have a remarkably reinforced effect on improving the buckling bearing capacity of the FG-GPLRC arch and the buckling load of the arch increases with an increase in WGPL. In addition, the effects of the geometry and size of GPLs on the buckling load are depicted in Figure  6. For bGPL/tGPL = 10 3 , the reinforced effect of GPLs on improving the buckling bearing capacity of the FG-GPLRC arch becomes stable under different ratios of aGPL/bGPL. Figure 7 presents the flexibility of the rotational constraints on the buckling load of the FG-GPLRC arch. It is noted that, in all cases, the buckling load decreases as the flexibility of the rotational constraints increases, and the fixed FG-GPLRC arch has the highest one, as expected.   Table 1 lists the buckling load of the arch having different GPL weight fractions. It is found that, by comparing with a pure epoxy arch (i.e., WGPL = 0.0%), GPLs have a remarkably reinforced effect on improving the buckling bearing capacity of the FG-GPLRC arch and the buckling load of the arch increases with an increase in WGPL. In addition, the effects of the geometry and size of GPLs on the buckling load are depicted in Figure 6. For b GPL /t GPL = 10 3 , the reinforced effect of GPLs on improving the buckling bearing capacity of the FG-GPLRC arch becomes stable under different ratios of a GPL /b GPL . Figure 7 presents the flexibility of the rotational constraints on the buckling load of the FG-GPLRC arch. It is noted that, in all cases, the buckling load decreases as the flexibility of the rotational constraints increases, and the fixed FG-GPLRC arch has the highest one, as expected.    ), the bifurcation point is located behind the limit point, so the buckling mode of the arch is the limit point instability mode. However, when the FG-GPLRC arch has a higher geometric parameter (e.g., 5.5 λ = ), its bifurcation point is located before the limit point, as shown in Figure 8b, leading to a bifurcation buckling mode of the arch instead of a limit point buckling mode. Accordingly, a critical geometric parameter 2 λ b is defined to identify the buckling mode of the arch, whose value is determined by solving Equations (27) and (34) Figure 8a that, when the arch has a smaller geometric parameter (e.g., λ = 2.6), the bifurcation point is located behind the limit point, so the buckling mode of the arch is the limit point instability mode. However, when the FG-GPLRC arch has a higher geometric parameter (e.g., λ = 5.5), its bifurcation point is located before the limit point, as shown in Figure 8b, leading to a bifurcation buckling mode of the arch instead of a limit point buckling mode. Accordingly, a critical geometric parameter λ b2 is defined to identify the buckling mode of the arch, whose value is determined by solving Equations (27) and (34) at β = β b simultaneously. When the FG-GPLRC arch has λ > λ b2 , the buckling mode of the arch is bifurcation mode. higher geometric parameter (e.g., 5.5 λ = ), its bifurcation point is located before the limit point, as shown in Figure 8b, leading to a bifurcation buckling mode of the arch instead of a limit point buckling mode. Accordingly, a critical geometric parameter 2 λ b is defined to identify the buckling mode of the arch, whose value is determined by solving Equations The effects of the rotational constraints and the power law index on the critical geometric parameters to distinguish the buckling mode are presented in Figure 9. As shown in this figure, in Region 1 (λ > λ b2 ), the buckling mode of the arch is bifurcation mode. In Region 2 (λ b1 ≤ λ ≤ λ b2 ), the buckling mode of the arch is either limit point instability mode or bifurcation mode, depending on which mode would occur first. In Region 3 (λ s1 ≤ λ < λ b1 ), the buckling mode of the arch is limit point instability mode. In Region 4 (λ < λ s1 ), no buckling occurs in the arch. It should be mentioned that there is no solution to λ b2 for α = 0 which refer to the fixed support, it implies that the FG-GPLRC arch with fixed ends cannot buckle in a bifurcation mode.   The nonlinear equilibrium paths of the FG-GPLRC arch, having the geometric parameters λ = λ s2 = 4.0861 and λ = λ s3 = 5.7552, are given in Figure 10. The number of limit points of the arch, in which the flexibility of the rotational constraints is set as two, is indicated in the figure. Note that we can determine the number of limit points within the bounds of successive inflexion points. When λ < λ s1 , there is no limit point in the nonlinear equilibrium path. When λ s1 < λ < λ s2 , there are two limit points in the nonlinear equilibrium path. When λ s2 < λ < λ s3 , there are four limit points in the nonlinear equilibrium path, and so on. Hence, the first 18 limit points and the corresponding critical geometric parameters λ sn of the arch with the power law index n = 1.5 and n = 3 are plotted in Figure 11 for further clarification. linear equilibrium path. When 1 2 λ λ λ < < s s , there are two limit points in the nonlinear equilibrium path. When 2 3 λ λ λ < < s s , there are four limit points in the nonlinear equilibrium path, and so on. Hence, the first 18 limit points and the corresponding critical geometric parameters λ sn of the arch with the power law index n = 1.5 and n = 3 are plotted in

Conclusions
In this paper, we studied the nonlinear buckling behavior of elastically supported FG-GPLRC arches with asymmetric distributed GPLs under a CPL. Analytical solutions for the limit point buckling, bifurcation buckling, and mode switching were derived. A finite element analysis was also employed to verify the present results, showing good accuracy of the present method for predicting the nonlinear buckling of the FG-GPLRC arch. According to the numerical studies, the influences of the distribution of GPLs, weight fractions, geometric configurations, as well as boundary restraints are discussed. It is found that the nonlinear buckling load of the FG-GPLRC arch decreases as the power law index or the flexibility of the rotational constraints increases. Switching of the buckling mode of the FG-GPLRC arch is quite sensitive to the power law index, the flexibility of rotational constraints, and the arch geometry. It is also found that, according to the present analytical method, the number of limit points in the nonlinear equilibrium path of the FG-GPLRC arch under a CPL can be easily determined.

Conclusions
In this paper, we studied the nonlinear buckling behavior of elastically supported FG-GPLRC arches with asymmetric distributed GPLs under a CPL. Analytical solutions for the limit point buckling, bifurcation buckling, and mode switching were derived. A finite element analysis was also employed to verify the present results, showing good accuracy of the present method for predicting the nonlinear buckling of the FG-GPLRC arch. According to the numerical studies, the influences of the distribution of GPLs, weight fractions, geometric configurations, as well as boundary restraints are discussed. It is found that the nonlinear buckling load of the FG-GPLRC arch decreases as the power law index or the flexibility of the rotational constraints increases. Switching of the buckling mode of the FG-GPLRC arch is quite sensitive to the power law index, the flexibility of rotational constraints, and the arch geometry. It is also found that, according to the present analytical method, the number of limit points in the nonlinear equilibrium path of the FG-GPLRC arch under a CPL can be easily determined.