Stability of a half-sine shallow arch under sinusoidal and step loads in thermal environment

The complex structural behavior of shallow arches can be remarkably affected by many parameters. In this paper, the structural responses of a half-sine pin-ended shallow arch under sinusoidal and step loadings are accurately calculated. Additionally, the effects of environmental temperature changes are considered. Three types of sinusoidal loadings are separately investigated. Displacements, load-bearing capacity, the magnitude of the axial force and the locus of critical points (including limit and bifurcation points) are directly obtained without tracing the corresponding equilibrium path. Furthermore, the boundaries identifying the number of critical points are investigated. All mentioned structural responses are formulized based on the rise of the arch and the environmental temperature change, which are introduced in a dimensionless form. The proposed formulation is also developed for generalized sinusoidal loadings. Additionally, the structural behavior of the shallow arch under two types of step loadings is investigated. Finally, the accuracy of the suggested approach is examined by a non-linear finite element method.


INTRODUCTION
Shallow arches are widely used in structural, mechanical and aerospace engineering, and the investigation of structural stability has always been of the researchers' interest.The failure of such structures is in the form of material failures, structural instability or a combination of them.
The tendency of structure to return to the static state, after creating a perturbation in the degrees of freedom, is called stability (Thompson and Hunt, 1973;Khalil, 2002).In the analysis of structural stability, since the structure experiences the sudden deformations, the investigation of critical points (such as limit and bifurcation point) is crucial.Such deformations cause severe changes in strains and stresses.The geometry of the arch is an influential parameter on its load-bearing capacity (Cai et al., 2012;Bateni and Eslami, 2015;Bradford et al., 2015;Rezaiee-Pajand and Rajabzadeh-Safaei, 2016).In addition, various loadings (e.g., the sinusoidal (Plaut and Johnson, 1981), concentrated (Pi et al., 2008;Chandra et al., 2012;Tsiatas and Babouskos, 2017), distributed (Moghaddasie and Stanciulescu, 2013b) and end moment loads (Chen and Liao, 2005;Chen and Lin, 2005)), geometric imperfections (Virgin et al., 2014;Zhou et al., 2015a), and boundary conditions (Pi and Bradford, 2012;Pi and Bradford, 2013;Han et al., 2016) are other important factors in the structural design.
In most cases, shallow arches become elastically unstable when the lateral load reaches a critical value (Chen and Li, 2006).This means that a large deformation could be observed while the material remains elastic.Practical experiences also confirm this issue (Chen and Liao, 2005;Chen and Yang, 2007a;Chen and Ro, 2009).Consequently, the behavior of shallow arches can be explained by the non-linear theory of elastic stability.In some analyses, it is assumed that the displacements of the arch are limited to avoid a material failure (Pippard, 1990;Xu et al., 2002;Chen and Hung, 2012).In addition, the variation in the environment temperature can be influential on the stability of structures (Matsunaga, 1996;Hung and Chen, 2012;Stanciulescu et al., 2012;Kiani and Eslami, 2013).
Several approaches can be applied to investigate the structural behavior of shallow arches.Previously, both analytical and numerical methods are discussed in the literature (Plaut and Johnson, 1981;Reddy and Volpi, 1992;Pi et al., 2002;Xenidis et al., 2013).In some analytical techniques, the displacement field is replaced by a set of orthogonal functions to derive the non-linear equilibrium and buckling equations (Xu et al., 2002;Chen et al., 2009;Chen and Hung, 2012;Moghaddasie and Stanciulescu, 2013b;Zhou et al., 2015a).Using the principle of stationary potential energy is another robust analytical approach to investigate the equilibrium and stability of shallow arches (Moon et al., 2007;Pi et al., 2007;Pi et al., 2008;Pi et al., 2010).On the other hand, the non-linear finite element method has been widely applied by researchers to trace the equilibrium path (Chandra et al., 2012;Saffari et al., 2012;Stanciulescu et al., 2012;Zhou et al., 2015b).Identifying the corresponding critical point(s) and finding the relationship between imperfections and load-bearing capacity are the capability of this numerical technique (Eriksson et al., 1999;Moghaddasie and Stanciulescu, 2013b;Rezaiee-Pajand and Moghaddasie, 2014).
This paper provides an analytical method to find the exact response of the half-sine shallow arch under the sinusoidal and step loads.Furthermore, the effect of temperature change on the equilibrium paths is investigated.For this purpose, the displacements of the structure are rewritten in the form of the Fourier series.By the substitution of displacements into the governing equations of the arch, the initial and bifurcated equilibrium path are obtained.On the other hand, the critical (limit and bifurcation) points on the static paths are achieved when the stiffness matrix is singular.In this paper, the behavior of the shallow arch under five types of distributed loads are separately investigated by the suggested approach.
The advantages of the proposed method are: (1) obtaining the exact solution of displacement field, equilibrium paths and the locus of critical points, (2) performing one parametric analysis instead of multiple analyses with specified values, and (3) finding the critical points without tracing the equilibrium paths.On the other hand, some limitations of the supposed method can be listed as (1) the changes in environment temperature are gradual, (2) the theory of plane stress is applied, (3) the height of the arch is limited, and (4) the material remains elastic during the analysis.
In the following section, the governing equations of the half-sine shallow arch under an arbitrary load are provided and the relative equilibrium paths are obtained.Then, the way of finding the locus of critical (limit and bifurcation) points is proposed (Section 3).The behavior of the half-sine arch under a number of distributed loadings is investigated by using the suggested method in Section 4. Finally, concluding remarks are given.

THE GOVERNING EQUATIONS OF THE SHALLOW ARCH
In this section, the governing equations of a half-sine shallow arch under an arbitrary loading * Q are ob- tained (Figure 1(a)).In this regard, a modified Bernoulli beam theory with large transversal displacement is applied.The material is assumed isotropic and variations in the temperature are gradual (no transversal temperature gradient is considered).E , A , I ,  and *  , respectively, represent the Young's modulus, area, moment of inertia, mass density and thermal expansion coefficient, which are assumed to be constant over the span L .The assumptions used for the analysis are as follows: (1) The axial force is constant over the span (Xu et al., 2002;Plaut, 2009;Chen and Hung, 2012); (2) The material is elastic (Chen and Li, 2006); (3) out-of-plane deflec-tions are neglected (Chen and Yang, 2007a); and (4) The range of displacements and curvatures of the arch is small in comparison with the length of the span Xu et al., 2002).Given the above assumptions, the equation of motion can be written as follows (Plaut and Johnson, 1981;Chen and Yang, 2007b;Chen et al., 2009): where, 0 y is the initial shape of the arch, the subscripts ",t " and ", x " show the partial differentiation with respect to the time and the longitudinal position of arch, respectively.The value of the axial force Here, 0 T  denotes environmental temperature changes.
Since the supports are fixed at the ends, the displacement of the arch in the direction x is neglected.In addi- tion, the temperature change can cause the axial force in the structure, which is denoted by the term Eqs. ( 1) and ( 2) could be rewritten in a dimensionless form: where, Here, r represents the radius of gyration of the cross section calculated by . The boundary conditions for Eq. ( 3) is as follows: By considering the boundary conditions, the initial and deformed shapes of the arch can be rewritten in the Fourier series form: In Eq. ( 7), h is the initial dimensionless rise of the arch.The external load Q based on Fourier series can be written as By substituting Eqs. ( 7)-( 9) into (3), a set of equations will be obtained: Latin American Journal of Solids and Structures, 2018, 15(8), e100 5/39 Here, and p is as follows: is equal to zero.Therefore, the Eq. ( 11) is written as: where, n R is the unbalanced force.The set of equations in (13) denotes the equilibrium state dependent on the external load Q .

THE CRITICAL POINTS
In this section, the critical load of the shallow arch is addressed.If the external load is a function of an independent parameter  ( , the solution of Eq. ( 13) results in a relationship between the displacement u  and the load factor  .In the other words, the equilibrium states in the space of ( , ) u  represent a number of curves that are called equilibrium paths.An example of equilibrium paths is shown in Figure 2.Each point on the curves represents the position of an equilibrium state relative to the load factor.In some equilibrium states, sudden changes can be observed in the behavior of the structure.These such points, which are part of the equilibrium path, are called the critical points.The critical points are categorized into limit and bifurcation points.At limit points, the slope of the equilibrium path is zero (Point A in Figure 2), while bifurcation points are located at the intersection of equilibrium paths (Point B in Figure 2).One way to obtain critical points is equating the determinant of tangent stiffness matrix to zero.The (modal) tangent stiffness matrix is calculated by the derivation of the unbalanced force.This can be done by substituting the magnitude of p from Eq. ( 12) into Eq.( 13) and taking derivatives with respect to m  : Here, nm  is the Kronecker delta.By equating the determinant of nm K to zero, the critical condition illustrat- ing limit and bifurcation points is obtained: The magnitude of the determinant n m nm     in Eq. ( 15) can be calculated by the following procedure: By using algebraic operations, Eq. ( 17) becomes the determinant of an upper triangle matrix: Consequently, Eq. ( 15) is rewritten as

RESULTS FOR LOADING PATTERNS
In this section, the behavior of shallow arch under a number of distributed loads are separately investigated.The patterns of loadings, respectively, are half-sine

1(e)) and asymmetric step function
In this way, a new formulation is proposed to achieve the relationship between displacements and load parameter.The result is verified by a finite element method (FEM).
In addition, the equilibrium paths and the locus of critical points of the arch for each loading are drawn.

Half-sine loading
By considering the type of loading shown in Figure 1(b), the values of n q for 1,2,... n  are equal to 1 0, 2, 3,...
For 0 n q  , two types of structural responses are obtained.In the former type which is corresponding to the initial equilibrium path, the parameters n  are equal to zero for 2, 3,... n  , while there is a non-zero n  (for instance, j  ) in the latter type (bifurcation path).The parameters n  for the initial equilibrium path is obtained from Eqs. ( 13) and ( 21): By substituting n  into Eq.( 12), the value of p is calculated: From this equation,  is as follows: By substituting (22) into Eq.( 8), the displacement field of the shallow arch is obtained: This equation is compatible with the results previously presented in the literature for a pin-ended shallow arch under a half-sine distributed loading (Plaut and Johnson, 1981).Eqs. ( 24) and ( 25) reveals the equilibrium path in the space ( , ) u  for the different values of p .For example, the displacement in the middle of the span ) is equal to Figure 3 shows the equilibrium path for four different values of h and 0 T  .In this diagram, the solid curves are the obtained initial equilibrium paths from Eqs. ( 24) and ( 26).In order to verify the suggested method, a non-linear FEM procedure is applied.In Figure 3, the signs  represent the structural responses obtained by FEM.Here, 60 Timoshenko beam elements with large displacements are used for modeling the shallow arch (Reddy, 2004).Each element includes six DoFs in the space of ( , ) u  .The equilibrium paths are traced by an in- cremental-iterative procedure (Crisfield, 1991;Crisfield, 1997).In this way, a modified cylindrical arc-length method is applied to find the next static state from the previous one (Moghaddasie and Stanciulescu, 2013a;Rezaiee-Pajand and Moghaddasie, 2014).The procedure of non-linear FEM begins from the initial unloaded state ( 0   ) and traces the equilibrium path for both directions 0   and 0   .All calculations were performed with the software Wolfram Mathematica 10.  and Structures, 2018, 15(8), e100 8/39


The comparison between the formation of curves and the result of finite element method displays the performance of the proposed strategy.It is noteworthy that the non-linear FEM procedure obtains a number of discrete equilibrium points, while a continuous equilibrium curve is given by the suggested method.
As it is mentioned previously, there is a non-zero term in bifurcation paths ( 0 j   ).From Eq. ( 13), the value p is obtained ( 2 p j   ), and by considering Eq. ( 12), the magnitude of j  is calculated: Similar to the initial equilibrium path, the displacement field relative to the j th bifurcation path is obtained by substituting the coefficients n  into Eq.( 8): This equation shows a linear relationship between  and u  .By considering / 2    , the displacement of the midpoint in the j th bifurcation path is computed: The gray curves in Figure 4 represent the calculated bifurcation paths given by Eq. ( 29).
Latin American Journal of Solids and Structures, 2018, 15(8), e100 9/39  A particular case which can be of interest is the relationship between the external load parameter and the axial force along the equilibrium path.Figure 6


The solid black and gray curves are relative to the initial and bifurcated paths, respectively.Note that, negative values for the axial force p represent that the shallow arch is in compression.As it is seen, the arch is always in tension for the case (a).By increasing the parameters h and 0 T  , the magnitude of p becomes negative in some parts of initial equilibrium paths.All bifurcation paths happen when the arch is in compression.
In order to obtain limit points, Eq. ( 19) can be rewritten in a simpler form: By substituting the obtained values n  from Eq. ( 22) into Eq.( 30), the critical load is calculated: In this equation, cr  represents the critical load.Eqs. ( 23) and (31) display the locus of limit points in the space 0 , ) ( , Figure 7 shows the relationship between the magnitude of critical load and the values of h and  The projection of the surfaces displayed in Figure 7 on the plane of   0 , h T  draws a boundary which is identifying the number of limit points on the equilibrium path.Figure 8 illustrates the mentioned boundary and the locus of states (a)-(d) in Figure 3.As a result, the initial equilibrium paths corresponding to the upper side of the boundary (e.g.states (b)-(d)) include two limit points, while there is no limit point for the state (a).This issue would also be realized by the investigation of critical surfaces passing over the supposed h and 0 T  in Figure 7.   As previously mentioned, the bifurcation points have the following characteristics: Latin American Journal of Solids and Structures, 2018, 15(8), e100 13/39 Since bifurcation points are located on both initial and bifurcation paths, these points include all properties of both paths (especially, the condition 0 j   ).By considering Eq. ( 27), a relationship between h , 0 T  and cr  is obtained: The equation 0  In a similar way, the boundaries which are identifying the number of bifurcation points on the equilibrium path ( B B ), can be derived by projecting the surfaces on the plane of h and 0 T  .For this purpose, the con- straint 0 j cr c     should be satisfied.This constraint concludes cr h   .Consequently, the following formu- lation for the set of boundaries B B is obtained from (34): Figure 10 shows the boundaries B j B for different values of j . .By substituting the value of 2 Q in Eq. ( 10), the magnitude of n q is obtained as follows: Latin American Journal of Solids and Structures, 2018, 15(8), e100 15/39 2 1 0 , 3, 4,....
If the procedure, which is previously described in the Subsection 4.1, is applied, the values of n  and p will be calculated: Additionally, the dimensionless displacement field for the initial equilibrium path is obtained from Eq. ( 8): Consequently, the displacement of the midpoint is as follows: Since the procedure of FEM begins from the initial unloaded state (which is known a priori) and is capable to trace only the paths passing through this state, the secondary equilibrium paths cannot be achieved by FEM.
Eqs. ( 42) and ( 43), respectively, show the displacement field and displacement of the midpoint for the j th bifurcation state: In Figure 12, bifurcation equilibrium paths are shown by gray curves.Although the initial equilibrium path does not include any critical point, the second path has limit and bifurcation points.In order to obtain the locus of limit points, the values of n  given by Eq. ( 37) are substituted into Eq.( 30): Eqs. ( 38) and (44) reveal the location of limit points as a set of surfaces in the space of Latin American Journal of Solids and Structures, 2018, 15(8), e100 20/39 Figure 16 shows surfaces . The green, red and blue surfaces represent the magnitude of bifurcation points corresponding to the first, second and third bifurcated paths, respectively.
Furthermore, the calculated parameters p , u  and Mid u  are respectively given by Eqs. ( 49)-( 51):  In a similar way, the displacement field and displacement of midpoint for the j th bifurcation path can be cal- culated: By considering Eq. ( 12), j  is computed: In Figure 18, the bifurcation paths are denoted by gray curves.Figure 18 shows that the initial equilibrium path does not include any critical point, and all critical points (limit and bifurcation) are located on the second equilibrium path.By substituting the values of n  into Eq.( 30), the critical load is calculated: The locus of limit points in the interval   On the other hand, the bifurcation points should satisfy the following conditions: 2 0, 2, 4, 5,..., , 2, 4, 5,...
By substituting the value p from Eq. ( 56) into Eq.( 49), the location of bifurcation points in the space of 0 , ) ( , The locus of bifurcation points for the interval   ) on the structural behavior of the shallow arch is investigated.The parameter k describes the formation of the external loading.Based on this supposition, the magnitude of n q for 1,2, 3, n   can be calculated for different values of k : By considering Eqs. ( 13), ( 58) and ( 59), the coefficients n  are obtained in a generalized form: ( ) The States I and II are corresponding to 1 k  , while the others are relative to the condition 1 k  .Addi- tionally, the States I and III can describe the initial equilibrium path.For this purpose, the magnitude of axial force is calculated by substituting Eqs. ( 60) and ( 62) into (12): Subsequently, the dimensionless displacement field for the initial equilibrium path is achieved by considering Eq.( 8) for the States I and III: On the other hand, the States II and IV, are corresponding to the bifurcation equilibrium path.Similarly, the displacement field for the j th bifurcation path can be derived: The displacement of the midpoint is obtained when for the four mentioned states.The parameter j  for the States II and IV are calculated by substituting Eqs. ( 61) and ( 63) into ( 12) and considering In order to find limit points on the initial equilibrium path, Eq. ( 30) is applied: Furthermore, the condition (74) should be satisfied for the j th bifurcation point: By substituting Eq. ( 74) into Eqs.( 64) and ( 65), the locus of bifurcation points is computed: It is noteworthy that bifurcation points are located on both initial and bifurcated paths.
Latin American Journal of Solids and Structures, 2018, 15(8), e100 27/39 4.5 Symmetric step loading A type of symmetric step load is shown in Figure 1(e).This loading pattern can be defined in the following form:   0 0 0.25 0.7 and 5 0.25 0.75 By using the Fourier series, the values of n q for 1,2,... n  are obtained: Note that for even values of n , the magnitude of n q is equal to zero.If the procedure, which is previously de- scribed, is applied, the values of n  and p will be calculated: where, the function   13).By considering Eq. ( 12), 2 j  is achieved: Eqs. and (85), respectively, show the displacement field and displacement of the midpoint for the j th bifurcation state: In Figure 23, the bifurcation paths are given. into Eq.( 30), the critical load is calculated: The locus of limit points in the interval   The locus of bifurcation points for the interval   The values of n q for 1,2,... n  are achieved by using the Fourier series: It is noteworthy that the magnitude of n q is equal to zero when 4 n j  (for 1,2, j  ).Similar to Subsec- tion 4.5, the values of n  and p can be calculated: where, the function The initial equilibrium paths for different values of h and In the case of asymmetric loading, there can be a non-zero 4 j  for the j th bifurcation path.Consequently, the axial force p equals 2 16 j  according to Eq. ( 13).By considering Eq. ( 12), 4 j  will be obtained: Accordingly, the displacement field and displacement of the midpoint for the j th bifurcation state are as fol- lows: In the asymmetric step loading, there is only one bifurcation path which can be seen in the case of 9.5 h  and 0 8.0 T   (the gray solid curve in Figure 27).Similar to the previous subsection, the locus of limit points can be determined.In this way, the critical constraint (30) is rewritten in the following form by considering Eq. ( 90): The locus of limit points in the interval   As it is observed, based on the magnitude of h and 0 T  , the number of limit points could be zero, two, four, six or eight along the corresponding equilibrium path for this interval.

CONCLUSIONS
The stability behavior of shallow arches is always being of the researchers' interest.In this paper, an analytical method to find the exact solution of a half-sinusoidal elastic shallow arch in the thermal environment under sinusoidal and step loads is proposed.For this purpose, the structural displacement is rewritten in a form of Fourier series, and subsequently, both initial and bifurcated equilibrium paths are obtained by substituting the transformed displacements into the governing equations of the arch.In addition, the critical points (such as limit and bifurcation points) are calculated by equating the determinant of stiffness matrix to zero.Furthermore, a new generalized formulation for various types of sinusoidal loadings is proposed.
In this research, the stability behavior of a half-sine shallow arch under three types of sinusoidal and two types of step function loads is separately investigated.Simultaneously, a non-linear finite element method is applied to show the accuracy and robustness of the suggested approach.In some cases, FEM becomes divergent during the path following procedure, while the proposed method is able to obtain the equilibrium path(s) com-   .6)

Figure 2 :
Figure 2: Primary (black) and bifurcation (gray) equilibrium paths Stability of a half-sine shallow arch under sinusoidal and step loads in thermal environment Latin American Journal of Solids

Figure 3 :
Figure 3: Initial equilibrium paths for four different values of h and

Figure 4 :
Figure 4: Initial and bifurcation paths for different values of h and

Figure 5 :
Figure 5: Equilibrium states corresponding to the points a-g shown in Figure 4(b) draws this relationship by considering Eq. (23).In this figure, four cases corresponding to the values of h and 0 T  given in Figure 3(a)-(d) are shown.

Figure 6 :
Figure 6: Relationship between load parameter and axial force for four different values of h and    .As it can be seen, depending on the values of

Figure 7 :
Figure 7: The locus of limit points in the space of   0 , , cr

Figure 8 :
Figure 8: The boundary identifying the number of limit points in the space of   0 , h T  describe the value of critical loads cr  corresponding to bifurcation points on the equilibrium path.According to the values of h and 0 T  , the number of bifurcation points can be zero, two, four, six and eight.In Figure9, the green, red, blue and yellow surfaces represent the magnitude of critical points corresponding to the first, second, third and fourth bifurcation paths, respectively.

Figure 9 :
Figure 9: The locus of bifurcation points in the space of   0 , , cr

Figure 10 :
Figure 10: The boundaries identifying the number of bifurcation points in the space of   0 , h T 

Figure 11
Figure 11 displays the equilibrium paths for four different values of h and

39 Figure 11 :
Figure 11: Initial and secondary equilibrium paths for four different values of h and

Figure 12 :
Figure 12: Initial, secondary and bifurcation paths for different values of h and

Figure 13 :
Figure 13: Relationship between load parameter and axial force for four different values of h and

Figure 14 :
Figure 14: Equilibrium states corresponding to the points a-g shown in Figure 12(a) figure, each surface specifies a couple of critical points corresponding to the supposed h and

Figure 15 :
Figure 15: The locus of limit points in the space of   0 , , cr

Figure 16 :
Figure 16: The locus of bifurcation points in the space of   0 , , cr Figure17shows the equilibrium paths for four different values of h and

Figure 17 :
Figure 17: Initial and secondary equilibrium paths for four different values of h and

Figure 18 :
Figure 18: Initial, secondary and bifurcation paths for different values of h and

Figure 20 :
Figure 20: The locus of limit points in the space of   0 , , cr

Figure 21 .
Figure 21.The green, red and blue surfaces represent the magnitude of bifurcation points corresponding to the first, second and third bifurcated paths, respectively.

Figure 21 :
Figure 21: The locus of bifurcation points in the space of   0 , , cr

Figure 22 39 Figure 22 :
Figure 22 displays the equilibrium paths for four different values of h and

Figure 23 :
Figure 23: Initial and bifurcation paths for different values of h and Figure 24.In this figure, each surface demonstrates a couple of limit points corresponding to the specific h and 0 T  .

Figure 24 :
Figure 24: The locus of limit points in the space of   0 , , cr

Figure 25 .
Figure 25.The green and red surfaces represent the magnitude of bifurcation points corresponding to the first and second bifurcated paths, respectively.

Figure 25 :
Figure 25: The locus of bifurcation points in the space of   0 , , cr

39 Figure 26 :
Figure 26: Initial equilibrium paths for four different values of h and

Figure
Figure 27: Initial and bifurcation paths for 9.5 h  and

Figure 28 :Figure 29 :
Figure 28: The locus of limit points in the space of   0 , , cr  are defined and their magnitudes are given for some values of i :