Nonlinear stability analysis of FGM shallow arches under an arbitrary concentrated radial force

The nonlinear in-plane buckling behaviour of pinned-pinned shallow circular arches made of functionally graded material (FGM) is investigated using a one-dimensional Euler–Bernoulli model. The effect of the bending moment on the membrane strain is included. The material properties vary along the thickness of the arch. The external load is a radial concentrated force at an arbitrary position. The pre-buckling and buckled equilibrium equations are derived using the principle of virtual work. Analytical solutions are given for both bifurcation and limit point buckling. Comprehensive studies are performed to find the effect of various parameters on the buckling load and on the behaviour. The major findings: (1) arches have multiple stable and unstable equilibria; (2) when the load is at the crown, the lowest buckling load is related to bifurcation buckling for most geometries and material compositions but when the external force is replaced, only limit point buckling is possible; (3) the position of the load has huge influence on the buckling load; (4) such arches are sensitive to small loading imperfections when loaded in the vicinity of the crown. The paper intends to improve and extend the existing knowledge about the behaviour of pinned-pinned shallow FGM arches under arbitrary concentrated radial load.


Introduction
Curved members are frequently used in engineering structures, like in bridges or roofs, because of their favourable load carrying capabilities. In practise, arches may be supported by pins at the ends to the ground or to adjacent members. Such members are often made of functionally graded material (FGM) with material properties continuously varying over the cross section.
Comprehensive studies on the in-plane buckling of arches were published in the past (Timoshenko and Gere 1961;Szeidl 1975;Simitses 1976;Bradford et al. 2002;Simitses and Hodges 2006;Bazant and Cedolin 2010). It is therefore well-known that geometrically nonlinear model is needed when analysing the stability of shallow members because the deformations prior to buckling can be significant and the behaviour can become nonlinear. Since then, various circular arch problems with numerous methods have been solved. Articles Plaut (2015a, b) focus on the critical displacements instead of critical forces. The analytical solutions are given for limit point buckling with pinned boundary conditions. Paper Silveira et al. (2013) presents a new numerical strategy for the nonlinear equilibrium and stability analysis of slender curved elements under unilateral constraints. It clarifies the influence of the foundation position and its stiffness on the nonlinear behavior and stability. Study (Zhou et al. 2015) report about analytical solutions of the stability and remote, unconnected equilibria of arbitrary shallow arches with non-symmetric geometric imperfections. In Asgari et al. (2014), a procedure is presented to trace the primary equilibrium path of functionally graded, pinned shallow arches under inplane thermal load using the nonlinear Donnell shell theory. Investigations on the static buckling behaviour of steel arches with elastic horizontal and rotational end restraints with simulations and experiments can be found in Han et al. (2016). The authors of (Tsiatas and Babouskos 2017) develop an integral equation solution to the linear and nonlinear behaviour of thin shallow arches with varying cross section under a central concentrated load. The coupled differential equations are solved with the analog equation method. In Abdelgawad et al. (2013) the conditions of the static and dynamic snap-through of shallow arches resting on a two-parameter elastic foundation are studied. Article Rastgo et al. (2005) report on the thermal in-and out-of-plane buckling of FGM curved beams using the Galerkin method. Contributions (Liu and Jeffers 2017;Liu et al. 2017b;Liu and Jeffers 2018) present the isogeometric analysis of FG structural members, including snap-through instabilities in thermal environment.
Analytical solutions for pinned-pinned and fixedfixed shallow arches subjected to a vertical dead load at the crown point are provided in Bradford et al. (2002) and . Moreover, paper  generalizes the results achieved in Bradford et al. (2002) for the case when the ends of the arch are rotationally restrained (there are pins with linear torsional springs). Further studies are published in Pi and Bradford (2012a) under the assumption that the torsional springs at the ends are different. Similar investigations, though, with a refined model are presented in Kiss (2014), Kiss and Szeidl (2015a) and Kiss and Szeidl (2015b) valid for cross sectional inhomogeneity.
Articles (Ting Yan et al. 2017, 2018a report on the static stability of non-uniform shallow arches made of homogeneous material. Non-uniformity means three constant stiffness regions along the length of the arch. The Euler-Bernoulli hypothesis is used to find analytical solutions. Concentrated, distributed loads and various support conditions are investigated with geometrical imperfections also included. It is found that arches with stiffer center are favourable compared to arches with softer center. Studies Eslami 2014, 2015;Babaei et al. 2018) focus on the stability of FGM arches. Comprehensive analytical investigations are presented for both central concentrated and constant distributed loads using the classical shallow arch theory.
Despite concentrated crown load is really preferred in the literature, in practice, the load may act elsewhere. However, the influence of the load position has been under much less investigations. Early results for homogeneous members (Schreyer 1972;Plaut 1979) show the influence of the load position on the buckling load of shallow arches. The later article neglects the effect of axial compressive force and so overestimates the buckling load. Recently, some progress has been made in Pi et al. (2017) and Liu et al. (2017a). The authors have presented and evaluated an analytical Euler-Bernoulli model-the same as in Bradford et al. (2002)-to find the nonlinear equilibria and buckling load of fixed-fixed and pinnedfixed homogeneous shallow arches under arbitrary concentrated radial load. The authors have found that such arches have multiple stable equilibria and the load position has significant effects on the in-plane behaviour.
Based on the above literature research, in this paper it is intended to find how pinned-pinned FGM shallow arches behave when loaded by an arbitrary concentrated radial load. A refined single-layer Euler-Bernoulli theory is used and the equilibrium equations are gained from the principle of virtual work. Analytical solutions are provided for both limit point and bifurcation buckling to find the lowest buckling loads and nonlinear equilibria. Comprehensive studies are carried out to show the effect of certain geometric and material parameters on the buckling load and behaviour. The paper intends to improve and extend the existing results about pinned-pinned shallow arches under arbitrary concentrated radial load. The model introduced has further potential applications to dynamic problems, axially FG beams, non-shallow arches and to various structures and mechanisms with multiple stable configurations. Moreover, the analytical solutions can be used as benchmark for numerical solutions.

Mechanical model
A pinned-pinned circular shallow arch is considered. The (E-weighted) centroidal axis is shown in Fig. 1, where the initial radius of curvature is q, the included angle is h ¼ 2#, the length of the arch is S while u and s are the angle and arc coordinates. The radial concentrated load P is applied at u ¼ a. The arch has uniform, symmetric cross section made of FGM with the material composition being independent of n. The material properties, like the Young modulus E, are assumed to follow a power-law distribution along the thickness as (Babaei et al. 2018) with subscripts c and m noting the ceramic and metallic constituents while b is the height of the cross section and k is the volume fraction index. Axis g is a major principal axis. Using the Euler-Bernoulli hypothesis, the longitudinal normal strain at an arbitrary point is (Kiss and Szeidl 2015b) Here e L is the strain on the centroidal axis according to the classical linear theory (Timoshenko and Gere 1961), j is the curvature and w is the mid-plane rotation about g. Therefore, the membrane strain at an arbitrary point on the centroidal axis (f ¼ 0) is given by (Kiss and Szeidl 2015a;Pi and Bradford 2012b) where u and w are the circumferential and lateral displacements. The third term on the right side is the source of nonlinearity-moderately large rotations are assumed. The E-weighted first moment of the cross section to g on the centroidal axis vanishes: Q e ¼ 0see Eq. (7). The material behaviour is assumed to be linearly elastic and isotropic, thus according to the Hooke law, the axial force N and the bending moment M can be expressed as (Kiss and Szeidl 2015b) Thus, the effect of the bending moment on the membrane strain is accounted-see ''Inner forces in terms of displacements'' section in ''Appendix'' for the details. In the recent articles Liu et al. 2017a;Ting Yan et al. 2018d) the following approximations can be found: In the above equations are the E-weighted area, E-weighted first moment-and moment of inertia to major principal axis g: 3 Pre and post-buckling equilibrium Since the pre-buckling deformations are significant for shallow arches (Dimopoulos and Gantes 2008a, b), these effects must be included. The principle of virtual work is used to find the nonlinear pre-buckling equilibrium configuration of pinned-pinned shallow arches under an arbitrary concentrated radial load. This principle requires that Z V r n de n dV À Pdwj u¼a ¼ 0 to be satisfied for every kinematically admissible virtual field dðÞ. From this principle-with the details gathered in ''Pre-buckling equilibrium equations'' section in ''Appendix''-the following equilibrium equations are found In terms of strains and displacements-see ''Prebuckling equilibrium equations'' section in ''Appendix''-these become All of Eq. (12) depend on the material composition and geometry, and k is the (E-weighted) modified slenderness of the arch (Liu et al. 2017a). So, the effect of the geometry and material distribution is incorporated into the model through these parameters, with r being the (E-weighted) radius of gyration of the cross section about its major principal axis g. Equation (11) 1 expresses that the membrane strain is constant while Eq. (11) 2 is comparable to the likes of Ting Yan et al. 2018d). The major reasons of the differences: the simplification q þ f ' q is not implied throughout this model, the effect of the circumferential displacement field u is kept in w when deriving Eq. (11) 2 and the inequality N ) M=q is not considered.
Differential equations (11) are paired with boundary and (dis-)continuity conditions. The displacements and bending moment are zero at the supports, while all fields are continuous at the tip of the external force except the shear force distribution which has a discontinuity with magnitude P. Recalling Eqs. (5), all these conditions can be given in terms of the dimensionless lateral displacement as where P is the dimensionless load defined by P ¼ Pq 2 # 2I e . To find the post-buckling equilibrium state, let us decompose the typical quantities. Let, e.g., the total lateral displacement be w Ã ¼ w þ w b , where w is the pre-buckling-and w b is the buckling displacement. With this rule, similarly as in ''Inner forces in terms of displacements'' section in ''Appendix'', the following expressions are found for the buckling strain, axial force and bending moment: Applying the principle of virtual work for the buckled equilibrium configuration Z has to be fulfilled. Following the integration by parts theorem, as detailed in ''Buckled equilibrium configuration'' section in ''Appendix'', it is found that the inner forces should satisfy equilibrium equations As shown in ''Buckled equilibrium configuration'' section in ''Appendix'', these relations can be given in terms of the kinematical quantities as These results are again comparable to the achievements of Bateni and Eslami 2014). The related boundary and continuity conditions are with all the fields being continuous at u ¼ a as there is no load increment.

Pre-buckling state
The general solution satisfying equation (11) 2 is Plugging it to conditions (13)-(14), a system of linear algebraic equations can be set for the constants X i . The closed-form solutions are given in ''Coefficients in the pre-buckling displacement field and the pre-buckling average strain'' section in ''Appendix''.
Since the membrane strain is constant on the centroidal axis, it is therefore equal to the average membrane strain: Recalling Eq. (4), the nonlinear equilibrium condition for pinned-pinned shallow arches is found to be Here, the effect of the circumferential displacement on the rotation field was dropped. This assumption is generally accepted for shallow arches (Ting Yan et al. 2018a). The constants K 0 , K 1 and K 2 are given in ''Coefficients in the pre-buckling displacement field and the pre-buckling average strain'' section in ''Appendix''.

Buckled equilibrium
It is well known (Bradford et al. 2002) that for certain shallow arches it is possible that the pre-buckling equilibrium path changes to an infinitesimally close orthogonal path with e mb ¼ 0. Based on Eq. (23), this phenomenon is characterized by the homogeneous equation whose general solution is Equation (31) and boundary conditions (24) determine a homogeneous linear system of algebraic equations. The lowest physically possible nontrivial solution is sin v# ¼ 0 ! v# ¼ p (Bradford et al. 2002). Therefore, the (critical) buckling strain for bifurcation buckling is with the associated antisymmetric buckled arch shape being W b ðuÞ ¼ F 3 sinðpu=#Þ. Plugging Eq. (32) to (29) yields the buckling load for bifurcation buckling.
From (17), the average buckling strain for shallow arches is 1 2# which has to vanish for bifurcation buckling. The first term on the right is zero because of the supports. Since the pre-buckling and buckling lateral displacements are orthogonal for bifurcation buckling, the third term is as well zero, meaning that the second term should also vanish. This buckling displacement, however, vanishes only if W b ðuÞ ¼ ÀW b ðuÞ, i.e., when the buckling displacement is antisymmetric in u. It is only possible for radial loads exerted at the crown. Hence, if the position of the load is arbitrary ðu ¼ a 6 ¼ 0Þ, bifurcation buckling can not occur for pinned-pinned shallow arches. However, when there is strain increment after buckling, equation (23) governs the problem. The general solution to this inhomogeneous differential equation is With the boundary-(24) and continuity conditions (25-(26), determination of the constants Y i is possible with the results provided in ''Coefficients in the buckling displacement and average buckling strain increment'' section in ''Appendix''. Furthermore, using kinematical equations (17)-(18) with (22) leads to which is the nonlinear relationship between the strain parameter v and the dimensionless load P. Thus, a quadratic equation is achieved with constants L 0 ; L 1 and L 2 defined in ''Coefficients in the post-buckling average strain'' section in ''Appendix''. Solutions for limit point buckling can be evaluated in terms of the geometry and material distribution by solving nonlinear Eqs. (29) and (36) simultaneously.

Results and discussion
To verify the new model, first, the commercial finite element software Abaqus 6.13 (Abaqus 2017) was used. A square cross section was considered with r being 0.108 m and the constant Young modulus E of 200 GPa. In the one-dimensional, geometrically nonlinear beam model eighty B21 elements and the Static, Riks step were used. The boundary and load conditions were assigned in polar coordinate system. The equilibrium paths are shown by diamond symbols on Figs. 5, 6 and a very good correlation is found between the current model and the Abaqus computations. Some further comparisons were also carried out with focus on the buckling loads. The typical values are gathered in Table 1. Some additional comparisons were performed with the analytical results of (Bateni and Eslami 2014) valid for FG arches under crown load. The arch with rectangular cross section is made of Silicon nitride Si 3 N 4 and stainless steel SUS304. The Young modulus of the ceramic and metallic parts are E c ¼ 348:43 GPa and E m ¼ 201:04 GPa. Four power law indices were selected. A very good agreement was found as shown in Table 2. Using the same geometry and materials, further 3D Abaqus computations were also performed using C3D8R elements. The body was composed of 25 perfectly tied layers along the thickness with the effective material parameters being constant in each layer. Results are shown in Table 3.
Evaluation of the new model reveal that when the load is central, there are three typical kinds of behaviour. The geometrical limits for these are almost constant in the slenderness as shown in Fig. 2a. When kðl ¼ 10 3 . . .10 6 Þ % k\3:91, there is no buckling. After that, there is possibly limit point buckling with symmetric buckled shape and then, for most arches (k [ % 9:81), bifurcation buckling occurs first with antisymmetric buckled shape. When the external load is placed elsewhere, there is either no buckling or limit point buckling. The slenderness limits against the load position between these two ranges are shown in Fig. 2b for three magnitudes of l.
To demonstrate the significant effect of the modified slenderness and the load position on the nonlinear in-plane behaviour of shallow arches, some the typical equilibrium paths for S=r ¼ 200 are shown in Figs. 3, 4, 5 and 6 with W a ¼ Wðu ¼ aÞ=ð1 À cos #Þ as the dimensionless displacement at the tip of the load and e str: ¼ 2pr=S ð Þ 2 being the second mode buckling strain of a pin-ended column under compression with the same cross section and material distribution as the selected arch (Bradford et al. 2002). Structural members with k ¼ 3:1 show no buckling behaviour. There is nonlinear bending only, and the load position already has visible influence on the equilibrium path (Fig. 3). When the slenderness is set to 6.1, one upper and one lower limit point appear (Fig. 4) with two stable and one unstable branch. The buckling load increases with the selected load positions and the critical (buckling) strain is always less than for a straight member.
Increasing the slenderness to 10.1 (Fig. 5) reveals that the number of limit points increase with k.
Although there is still one upper and one lower limit point when the load is at the crown, there are two upper and two lower limits for any other load position. When a ¼ 0, there are two further bifurcation points, one on the primary stable and one on the remote stable branch. The critical strain for bifurcation buckling is lower and thus, it occurs first. However, for other load position parameters, there is only limit point buckling. The number of equilibrium branches also increases when a 6 ¼ 0: two stable ones (primary and remote) and three unstable ones can be found. One typical load is related to two stable and three unstable equilibrium points. It is also noticed that meanwhile the lowest buckling loads for a ¼ 0 and 0.6 are almost identical, the in-plane behaviour is completely different. At e m =e str: ¼ 1 there is an intersection point of the typical curves in Fig. 5b and the lowest critical strain is always less than e str . Further upper and lower limit points and equilibrium branches appear when k ¼ 14:3 with bifurcation points located on the stable branches for crown loads as shown in Fig. 6. The total number of limit points for a [ 0 is now six, with additional equilibrium branches. From the previous figures it is understood that the modified slenderness and the load position both have significant   Some related arch shapes in the pre-buckling equilibrium are shown in Fig. 7 with normalized displacements W against the arc coordinate s. The The effect of the load position parameter on the buckling load is evaluated in Figs. 8 and 9. Because the end supports are identical, only positive parameters are selected. When l ¼ 10 3 and # ¼ 0:4 the buckling load gradually increases with a. Otherwise, for any other selected parameter set, replacing the load from the crown causes the lowest buckling load decrease for a while. This decrease is more remarkable when l is smaller and # is greater. Overall, it can be seen, the load position can almost halve the buckling load compared to Pða ¼ 0Þ ¼ Pð0Þ. After reaching this lower limit, there is an increase and the load bearing capabilities of arches become even better when the load is placed in the direct vicinity of the supports. The effect of the load position and included angle have the greatest influence when l ¼ 10 3 and least when l ¼ 10 5 : Next, the lowest critical loads are shown in Figs. 10 and 11 in terms of k and # for arches with S=r ¼ 100 and 200. Initially, the critical load increases steeply and after reaching the peak, there is a slow decrease. It is also visible that a small imperfection in the load position-e.g., a=# ¼ 0:1 instead of 0-has quite a huge impact for most parameter sets. Increasing the ratio S / r shifts the typical curves to the left (buckling occurs for lower angles). It can be concluded that, independently of the load position parameter, arches with the same slenderness (or included angle) have better load carrying capabilities when S / r is greater.
In Figs. 12, 13 and 14 it is investigated how the dimensionless buckling load varies with the semi-vertex angle and the modified slenderness. Three magnitudes of parameter l are selected: 10 3 ; 10 4 ; 10 5 and four load positions are chosen for the illustrations: a=# ¼ 0; 0:1; 0:3; 0:6 f g : There are many common findings in these diagrams. For lowest k; # parameters, there is no buckling. When the load is at the crown, limit point buckling occurs first for small angles, then bifurcation type is dominant. When the load is replaced from the crown, only limit point buckling can occur. Due to a small perturbation in the load position parameter around the crown, (e.g., a=# ¼ 0:1-it is only 5% of the included angle), for a short while, the typical buckling load is actually unchanged but with # or k increasing, the effect of this small imperfection becomes really considerable. Depending on the parameters, the dimensionless buckling load can decrease by up to 20% compared to crown load. When the load position parameter is changed from zero, the critical load decreases for a while but after that there is an increase. This finding is in accord with Figs. 8 and 9. There is always an intersection between the curves for a ¼ 0 and a 6 ¼ 0 in the figures. Thus, for certain arches, the buckling loads might be identical but with completely different buckled shapes and nonlinear behaviour.
When comparing Figs. 12, 13 and 14 it is visible that increasing l makes the curves shift to the left: buckling behaviour appears for smaller angles and the related buckling load increases for any given # and a. When the semi-vertex angle and the load position are both prescribed, we find greater buckling loads when l is greater. When a ¼ 0 the findings correspond with the results in Kiss and Szeidl (2015b) valid for pinnedpinned shallow arches under a central concentrated load.
The next example attempts to demonstrate the possible effect of the material distribution on the buckling load. The arch with rectangular cross section is made of Silicon nitride Si 3 N 4 and stainless steel SUS304. The material distribution follows the Voigt rule of Eq. (1). If k ¼ 0, the arch is made of ceramic, while if k ¼ 10 5 (tends to infinity) the material is pure steel with one thin ceramic layer at the top. The power law index k is varied continuously between 0 and 5. In Fig. 15a it is shown how the ratios A e ðk 6 ¼ 0Þ=A e ðk ¼ 10 5 Þ; I e ðk 6 ¼ 0Þ=I e ðk ¼ 10 5 Þ-defined by Eq. (7)-depend on the power law index. The maximum value is 1.77 times of the homogeneous steel for both, and remain considerably above 1 between k ¼ 0. . .5. In Fig. 15b it is shown how the radius of the centroidal axis changes with the material distribution. The plotted ratio depends on q=b, but overall, this effect is not relevant. Next, it is investigated how the material distribution affects the buckling load. For k ¼ 10 5 , S / r was set to 200 and r ¼ 0:108 m, as formerly. Various central angles, load positions and radius/cross sectional height ratios were selected, and it is found that the results plotted in Fig. 16, i.e., the critical load ratio Pðk 6 ¼ 0Þ=Pðk ¼ 10 5 Þ, are independent of these with a good accuracy. It is also clear that the power law index has a huge impact on the lowest critical load.  The nonlinear in-plane buckling behaviour of FGM shallow circular arches under an arbitrary concentrated radial load was investigated using a refined beam model. The effect of the bending moment on the membrane strain was accounted. The pre-buckling and buckled equilibrium equations were obtained from the principle of virtual work. Analytical solutions to these equations were presented for both bifurcation and limit point buckling. Parametric studies revealed that arches have multiple stable and unstable equilibrium   branches. When the load is at the crown, the lowest buckling load is generally related to bifurcation buckling but for most arches and any other load positions, there is only limit point buckling. Changing the position of the load can almost halve the critical load compared to a force exerted at the crown. Moving the load close to the support causes the buckling load become even multiple times of crown load. Most arches are sensitive to small loading imperfections when the load is in the vicinity of the crown. This phenomenon requires further thorough investigations. The effect of the load position can become quite considerable as the load carrying capabilities decrease rapidly. The effects of various geometrical and material parameters on the behaviour and buckling load were shown through comprehensive studies.
Acknowledgements Open access funding provided by University of Miskolc (ME). This research was supported by the National Research, Development and Innovation Office-NKFIH, K115701.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix
Inner forces in terms of displacements Comparing Eqs. (2)-(4) and (7) yields the axial force and bending moment as in Eq. (5): Therefore It can easily be shown that similar relations can be set up for the increments N b and M b as given by Eq. (19).

Pre-buckling equilibrium equations
With dV ¼ 1 þ f=q ð Þ ds dA being the infinitesimal volume element, the first term in the principle of virtual work (8) can be expanded by recalling Eqs.
Coefficients in the pre-buckling displacement field and the pre-buckling average strain The coefficients in the pre-buckling displacement field (27) are given by with Closed-form solutions to the constants K 0 , K 1 and K 2 in (29) can be found after evaluating the integrals If a ¼ 0, the above solutions are identical to those in Kiss and Szeidl (2015b) set up for a concentrated load at the crown.

Buckled equilibrium configuration
Expanding Eq. (20), the principle of virtual work for the buckled configuration is Z V r n de nb dV þ Z V r n b de nb dV À Pdw b j s¼qa ¼ 0 : The pre-buckling equilibrium state is already known, thus the variation of such quantities is identically zero. Based on Eq. (16), the variation of the total axial strain is ð56Þ So the principle of virtual work is Integrating the first three terms by parts while using kinematical Eqs. (16)-(19) result the following relations Z V r n de 2 and finally, the third term in (57) is the same as (40) but for the increments. Altogether, because of the arbitrariness of the virtual quantities, the buckled equilibrium is governed by Eq. (21) After substituting (5) and (19) where the second term can be dropped, so e 0 mb ¼ 0. Furthermore, Eq. (21) simplifies to lead to that is identical to (23).