In-plane instability of shallow layered arches with interlayer slip

In this paper, a beam theory for predicting limit point buckling and bifurcation buckling of shallow arches composed of two layers flexibly bonded is presented. The flexibility of layer bond results in interlayer slip, which significantly affects the critical transverse loads. The presented theory is based on a layerwise assumption of the Euler–Bernoulli theory and a linear behavior of the interlayer. After establishing the equilibrium equations and boundary conditions, a numerical method for efficient solution of these equations is provided. In a first example, the presented theory is validated by comparative computations with a much more elaborate finite element analysis assuming a plane stress state. In several other examples, the effect of interlayer stiffness, load distribution and boundary conditions on the stable and unstable equilibrium paths of shallow arches with interlayer slip is investigated.


Introduction
When the critical transverse load is exceeded, an immovably supported shallow arch with fully constraint outof-plane displacement is subjected to in-plane buckling. This instability can take the form of limit point buckling (snap-through) or bifurcation buckling. While in a symmetric arch subjected to a symmetric transverse load the buckling shape of the snap-through mode is symmetric, in the bifurcation mode the shape is antisymmetric. The problem of in-plane buckling of shallow arches has been recognized for a long time and, accordingly, the results of numerous comprehensive studies have been published in the past (e.g., Kerr and El-Bayoumy [18], Lo and Conway [22], Virgin et al. [37], Bradford et al. [5], Tsiatas and Babouskos [36]).
The deformation of shallow arches can be substantial even before reaching the stability limit, and thus their response becomes nonlinear. The analysis of the stability and instability behavior of shallow arches should therefore be based on a geometrically nonlinear model [28]. If the shallow arch is elastically supported, the instability behavior may become much more complex (e.g., Pi et al. [27], Pi and Bradford [26], Han et al. [13]). For the estimation of the critical transverse load, it is generally sufficient to perform a static analysis. However, the deflection that occurs with further load increase is a dynamic process; thus, several papers consider the inertia terms to estimate the nonlinear dynamic response (e.g., Öz and Pakdemirli [25], Pi and Bradford [29], Keibolahi et al. [17], Zhong et al. [38]).
In most studies, homogeneous arches are considered, while limit point buckling and bifurcation buckling of layered and inhomogeneous curved structures have less frequently been analyzed so far (e.g., Heuer and Ziegler [16], Schultz et al. [34], Babaei and Eslami [4], Kiss [20], Chan et al. [9]). In particular, Kiss [19] presents a thorough analytical study of shallow arches made of functionally graded materials on the effects of various parameters on the buckling load and stable as well unstable equilibrium paths. However, to the best of the authors' knowledge, these studies refer solely to members whose layers are perfectly bonded.
In numerous applications, however, a perfect bond between the layers cannot be achieved because the fastener is flexible. In such a case, the longitudinal displacement of the fibers is not continuous over the height of the member, but exhibits a jump at the layer interface due to the flexible bonding, which is referred to as interlayer slip. The mechanical behavior of layered members with interlayer slip (e.g., Girhammar and Gopu [11], Heuer and Adam [15], Girhammar and Pan [12], Lorenzo et al. [23], Gahleitner and Schoeftner [10]) is much more complex than for homogeneous structures. Solutions for the geometric nonlinear behavior (e.g., Ranzi et al. [31], Adam et al. [2]) and buckling of columns with interlayer slip (e.g., KryŽanowski et al. [21], Schnabl and Planinc [32], Challamel and Girhammar [8], Schnabl et al. [33]) exist in the literature. To date, however, limit point buckling and bifurcation buckling of transversely layered shallow arches with interlayer slip subjected to transverse loads have not been investigated.
To fill this research gap, a nonlinear beam theory for shallow arches whose two layers are flexibly bonded is presented below to efficiently predict the stable and unstable deformation branches of these structures. After describing the strategy for numerically solving the differential equilibrium conditions, the effect of the interlayer slip on the snap-through and in-plane bifurcation on the response and the buckling loads is investigated in several application problems. To validate the presented theory, for one example problem the response is additionally determined with finite shell element analyses assuming plane stress and compared to the results of the presented theory.

Fundamental equations
The considered shallow arch shown in Fig. 1 is composed of two elastic homogeneous layers with constant cross-section along the span, which are elastically bonded in circumferential direction. In radial direction the layer connection is rigid. The parameters of the upper layer are identified by the subscript "1", and those of the lower layer by the subscript "2". The cross-sectional area is symmetrical in the vertical direction. The member axis corresponds to the line connecting the elastic centers of gravity of the arch with rigidly bonded layers. The member is immovably supported at both ends.
The shallow arch is referenced to the curvilinear coordinate x following the member axis, z is the in-plane coordinate perpendicular to x, and y is the out-of-plane coordinate perpendicular to x and z, see Fig. 1. The origin of the reference coordinate system is located in the member axis at the left support. The member axis is the line connecting the elastic centers of the composite cross section of the beam with rigidly bonded layers, whose distance d 1 (d 2 ) from the local centroid of the upper (lower) layer is determined as (see Fig. 1) In these relations, d denotes the distance between both layerwise centroids, E A 1 = E 1 A 1 is the product of Young's modulus E 1 and cross-sectional area A 1 of the upper layer, correspondingly for the lower layer In addition, the local coordinates (ζ 1 , η 1 ) and (ζ 2 , η 2 ) parallel to the z-and y-coordinates, respectively, are defined for the two layers, whose origins are in the respective local centroid of the considered layer. The variable R(x) represents the radius of principal curvature, which is in the present case of a shallow arch very large compared to the member length, and k(x) = R −1 denotes the principal curvature. The member is subjected to the transverse load per unit length p(x), see Fig. 1. Let u (∞) (x) and w(x) denote the components of the displacement at the member axis in x-and z-direction, and u(x) the interlayer slip, i.e., the relative displacement in circumferential direction between the layers at the interface. In a shallow arch with non-compressible cross section, it is reasonable to assume that the radial displacement w(x) is the same for each fiber and small compared to the principal radius of curvature R(x). Moreover, the two layers are slender and stiff, and thus, their shear deformation is negligible. Therefore, the validity of the kinematic assumptions of the Euler-Bernoulli theory can be assumed for each layer separately. Consequently, the displacement field of the shallow arch is fully described by the three kinematic variables w, u (∞) and u. The Euler-Bernoulli theory implies that the cross-sectional rotation χ is the same for both  layers. It is composed of a part due to the tangential displacement and a part due to the radial displacement difference, However, in several studies (e.g., Pi et al. [28], Bradford et al. [6]) it has been shown that for immovably supported nonlinear shallow curved arches, the effect of axial deformations u (∞) k on cross-sectional rotation and curvature is negligible [26]. Assuming that the member axis is in the lower layer, the longitudinal displacements u 2 (x) of the local layer axes can be expressed as a function of the governing kinematic variables w, u (∞) and u as follows, see Fig. 2, Then, the longitudinal displacements u 1 (x, ζ 1 ) (u 2 (x, ζ 2 )) in the radial distance ζ 1 (ζ 2 ) from the central axis of the upper (lower) layer follow from Since the instability of shallow arches is associated with moderately large displacements and membrane stresses, the longitudinal strains in the layer-wise central axes need to be formulated nonlinearly. For an arch structure, the full expression for the membrane strains with all nonlinear terms reads [7] However, in the case of an immovably supported shallow arch, only the nonlinearity 1 2 w 2 ,x is of significant order, while the other nonlinear terms are negligible, and therefore, the membrane strains in the layer axes can be approximated by the following expression (e.g., Heuer [14], Adam [1]): and further after inserting Eq. (3), The expression for the bending strain at distance ζ i from the ith local layer axis can also be simplified (e.g., Pi and Bradford [26]), because on the one hand the contribution of (u (∞) k) ,x to the curvature is negligible and on the other hand The total strain at the distance ζ i from the ith local layer axis is then Next, the stress resultants are expressed as a function of the governing kinematic variables and their derivatives. Multiplying the strains Eq. (9) by Young's modulus E 1 (E 2 ) gives the normal stresses, which are integrated over the cross-sectional area A 1 (A 2 ) to obtain the normal force in the upper (lower) layer, The sum of the layerwise axial forces yields the overall axial force as When calculating the layer-wise moments, the normal stresses are multiplied by ζ i before layerwise integration with the result (e.g., Ziegler [39]) where E J i = E i J i is the bending stiffness and J i the area moment of inertia about the η i -axis of the ith layer. The overall bending moment M is related to the layerwise stress resultants M i and N i and further with Eqs. (12) and (10) to the kinematic variables as with E J ∞ denoting the bending stiffness of the member with rigidly bonded layers, where E J 0 is the bending stiffness of the member with unbonded layers. The interlaminar shear traction, i.e., the longitudinal force transferred between the two layers in the interface, is proportional to the interlayer slip u, K s denotes the slip modulus, see, e.g., Girhammar and Pan [12],

Equilibrium equations
The differential equations of nonlinear equilibrium of the shallow arch can be found, for instance, with the minimum total potential energy principle δ = 0 (e.g., Ziegler [39]). The potential energy of the system reads as After inserting the strains Eqs. (9), (7) and (8) and some algebra, the potential energy as a function of the governing kinematic variables is obtained as follows Applying the principles of calculus of variations and partial integration to the potential and using the internal forces, we obtain By definition, the coefficients δu (∞) , δ u, and δw are arbitrary, and therefore Eq. (18) is met only if the following differential equations of equilibrium: and the boundary conditions are satisfied. The subscript b denotes a boundary at x = 0 and x = l, respectively. From Eq. (19) it follows that the overall axial force is constant along the arch. Therefore, in Eq.
Substituting the internal forces into Eqs. (19)-(21) leads to the differential equilibrium conditions in the governing kinematic variables, Equation (20) (respectively Eq. (24)) implies that the free-body diagram of an infinitesimal element of the upper layer is in equilibrium. Likewise, a free-body diagram of an infinitesimal element of the lower layer must be in equilibrium, i.e., N 2,x − t s = 0, or expressed in kinematic variables, In this relation neither the longitudinal displacement u (∞) nor the nonlinear terms appear. Therefore, in the following, this equilibrium equation is used instead of Eq. (24) to solve the present boundary value problem more efficiently.

Boundary conditions
As can be seen from Eq. (22), for the solution of the equilibrium Eqs. (27), (23), (25), a total of four boundary conditions per boundary must be defined. Here, shallow arches with three different boundary conditions are examined, i.e., soft-hinged support, hard-hinged support, and clamped support. Common to all of these supports is that the member axis is immovable at these points in all directions, see also Eq. (22), [2],

Soft-hinged support
In the case of a soft-hinged support (sh), the beam axis is hinged, i.e., (w ,x ) b = 0, resulting in zero total bending moment, M b = 0, see Eq. (22). Expressed in kinematic variables according to Eq. (13), it follows that [3] − Since the interlayer slip is not constrained at the support, i.e., u = 0, the axial force in the upper layer is therefore zero at this end, (N 1 ) b = 0 (see Eq. (22)) [3], or with the first of Eq. (10)

Hard-hinged support
The boundary condition Eq. (30) also applies at a hard hinged support (hh). In this case, however, a rigid plate at the member end prevents the interlayer slip (see, e.g., Adam et al. [2]),

Clamped support
At a clamped end (cl), the rotation of the cross-section is constrained (see, e.g., Adam et al. [2]), and furthermore the interlayer slip is constrained, see Eq. (32).

Analysis
The solution of the present boundary value problem is based on the Ritz approach for the deflection [30] w where the J − i a + 1 shape functions φ i (x) are polynomials. The initial value i a of the summation depends on the geometric boundary conditions in w to be satisfied at x = 0, and the exponent i b on the geometric boundary conditions to be satisfied at x = l. For the boundary conditions considered here, i a and i b are [30] soft-hinged (sh) and hard-hinged support (hh): The Ritz approach is substituted into the two equilibrium conditions Eqs. (23) and (27). These two ordinary differential equations are solved together with the corresponding boundary conditions (Eqs. (29); (31) for sh respectively Eq. (31) for hh and cl) for u and u (∞) analytically, which are then also a function of the unknown coefficients q i , denoted as u * and u (∞) * .
The coupled J − i a + 1 nonlinear equations in the coefficients q i , i = i a , ..., J , are found by applying the Galerkin method [39] to the third equilibrium condition Eq. (25). Therefore, Eq. (25) is consecutively multiplied by the J − i a + 1 shape functions φ i and integrated over l, yielding the following set of equations: N * is the approximation of the overall axial force obtained by substituting w * , u * and u (∞) * into Eq. (11) and evaluated at any x (for instance at x = l/2). Since the chosen Ritz approach violates the static boundary condition M b = 0 in the case of a hinged support, the work of this boundary moment was added in the above expression. Integration of this equation partially twice with respect to x, the order of the derivatives to x is reduced by two. The simultaneously appearing negative work of the boundary moment cancels out with the positive work of the boundary moment, and thus, Eq. (37) becomes The solution of the nonlinear coupled set of J − i a + 1 equations for q i , i = i a , ..., J , resulting from evaluation of this integral is found by numerical standard solvers. The approximation of the deformation of the shallow arch follows by substituting the found coefficients q i into w * (x), u * (x) and u (∞) * (x). It should be noted that with increasing number of shape functions not only the deformation but also the static boundary conditions are better approximated.

Shallow arch 1 subjected to distributed load
A first example is a shallow arch (hereafter referred to as shallow arch 1), fixed on both sides and soft-hinged supported (boundary conditions sh − sh), and curved in the shape of a sine half-wave against the positive z-direction,ŵ i.e., k(x) ≡ −ŵ ,x x = −ŵ 0 π 2 /l 2 sin π x l . The rise of this structure is 3 cm, i.e.,ŵ 0 = 0.03 m. The arch is composed of two elastically bonded layers and has a length l = 1.0 m. The cross-sectional dimensions of the upper layer with Young's modulus of E 1 = 7 · 10 10 N/m 2 are h 1 = 4 mm and b 1 = 0.1 m. The lower layer with Young's modulus E 2 = 1 · 10 10 N/m 2 has the same width as the upper layer, b 2 = b 1 , and thickness of αl = 14.97 and thus to a moderate interaction of the two layers [12]. A uniformly distributed transverse load is applied to the shallow arch, i.e., p(x) = p 0 , which leads to an unstable response above a certain critical load.
For the structural response analysis described in the previous section, the deflection w is approximated by 11 shape functions φ i , i = 1, ..., 11, according to the Ritz approach Eq. (34), which has been shown to be sufficient in a convergence study. Consequently, 11 nonlinear coupled equations in q i , i = 1, ..., 11 result from the Galerkin method Eq. (38). Not only the critical loads are of interest, but the entire response path as a function of load p 0 in both the stable and unstable branches. Therefore, the load p 0 is increased stepwise, and at each load step the real solutions of this system of equations are computed. These computations were performed with the software package Mathematica Mat [24] using the function "FindRoot" to find the real roots of this system of equations.
To validate the presented theory, the response of this shallow arch is additionally analyzed with the finite element (FE) method in the software suite Abaqus [35] assuming a plane stress state. In the corresponding numerical model, the two layers are discretized with the 8-noded quadrilateral plane stress elements of type CPS8R. A Poisson's ratio of 0.3 is assumed. The intermediate layer is modeled with 4-noded linear cohesive elements of type COH2D4. In the FE model, the interlayer has a certain thickness, which in this case is chosen as 0.002 mm. The stiffness normal to the interlayer, which is infinite in the beam theory, is chosen to be very large, i.e., 10,000 K s . In total, the structure is discretized into 18,500 finite elements, corresponding to 109,132 degrees of freedom, which is a multiple of the 11 degrees of freedom of the beam model. The soft-hinged supports are modeled as kinematic coupling of the outer surface of the lower layer at an additional node. To determine the complete branch of the response, the Riks arc-length algorithm implemented in Abaqus is used to solve the resulting system of equations. Figure 3 shows with thick lines the response of this shallow arch in terms of the kinematic variables found with the described beam theory. While the left column shows the non-dimensional external force p = p 0 l 3 /E J ∞ over the considered non-dimensional response variable at the indicated location of the beam axis, the right column displays the distribution of this response variable over the beam axis x/l at four discrete load levels, denoted A, B, C and D in the left column. The first row contains the normalized deflection w/l and the second row the normalized interlayer slip u/l. As can be seen from Fig. 3a, in the primary stable deformation branch the slope of the deflection at midspan w(0.5l)/l decreases with increasing load until the critical load p cr = 2.46 is reached at limit point A, where the deflection snaps through to the remote stable equilibrium path at point B. At load reversal, at the critical load p cr = 1.60 (limit point C), the arch snaps back to the primary stable deformation branch. As observed, the limit points are local extrema on the stable deformation branches. The stable deformation branches are illustrated by solid black lines, the unstable deformation branch between the limit points A and C is depicted by a dashed black line. The solution of the FE plane stress analysis represented by a thin red line with markers is virtually identical to the beam solution, confirming the assumptions of the proposed beam theory. The deflection shape shown in Fig. 3b is similar for the four load levels A, B, C, D given in Fig. 3a. The deflection shape depicted with green color (load level D) is located in the unstable deformation branch of the shallow arch.
The load p plotted against the interlayer slip u(l)/l at the right support has the form of a loop, see Fig. 3c. The largest value of the interlayer slip is found in the unstable response branch. When the structure snaps through from A to B, the interlayer slip at the supports becomes smaller in magnitude, but in the interior of the structure it is larger than in the other three discrete load levels B, C, D, as Fig. 3d shows. The interlayer slip from the beam theory and the FE analysis agree well.
To illustrate the importance of considering the interlayer flexibility when estimating the critical loads, Fig. 3a also shows the response of the shallow arch with rigidly bonded layers as thin black lines. As observed, on the primary stable deformation branch, the critical load p cr = 3.47, which is 41% larger than for the arch with elastically bonded layers. The critical load p on the remote stable deformation branch decreases from 1.60 to 1.12.
The internal forces shown in Fig. 4 also demonstrate the influence of the interlayer slip on the nonlinear stable and unstable response of the shallow arch. The total normal force N divided by E A e is much larger for the structure with rigid bond than that of the arch with flexibly bonded layers, see Fig. 4a. In Fig. 4b, the layer-wise normal forces N 1 /E A e and N 2 /E A e are plotted as a function of x/l for load levels A, B, C and D in addition to the total normal force N /E A e . It can be seen that at the supports the boundary condition (N 1 ) b = 0 (Eq. 31) is satisfied, while the total axial force is fully transferred to the second layer, i.e., (N 2 ) b = N . Figure 4c shows p versus the normalized total moment Ml/E J ∞ at midspan. Next to it, Fig. 4d shows the normalized

Shallow arch 2 subjected to distributed load
After the comparative FE analysis of the first example has confirmed the presented beam theory, the investigations are continued on a second shallow arch where the thickness of the two layers is smaller. The thickness of the upper layer is now h 1 = 2 mm and the thickness of the lower layer h 2 = 5 mm. The total height of the cross-section is therefore 7 mm. The rise of the shallow arch isŵ 0 = 25 mm. The slip modulus is reduced by one power of ten and is K s = 10 8 N/m 2 , and thus the composite action parameter according to Eq. (40) becomes αl = 10.41. All other parameters as well as the boundary conditions (sh − sh) of this structure referred to as shallow arch 2 are the same as in the previous example. Figure 5 shows the response of this structure in analogy to Fig. 3. The response of the symmetric shallow arch subjected to the symmetric transverse load p 0 based on the presented beam theory is shown in the figures of the left column with a thin red line with rectangular markers. Such an analysis can predict limit point buckling, but not bifurcation buckling. Therefore, a second analysis is performed in which the transverse load has a slight perturbation, i.e., the external load is slightly asymmetric in the form with H denoting the unit step-function. The total resultant of the transverse load thus remains unchanged, but its application point shifts slightly to the right of midspan. The result of this analysis is shown in the figures of the left column with thick black lines. Comparison of the outcomes of the two computations shows that this very slender shallow arch does indeed become unstable due to bifurcation buckling. As observed, a bifurcation point (point A) is located on the primary stable equilibrium path. At this point, the equilibrium path changes to an orthogonal path on which the total axial force remains constant until bifurcation point C is reached on the remote stable equilibrium path. The actual critical load at bifurcation point A is p cr = 3.61, compared to . It is noticeable that the response on the stable deformation branches up to the bifurcation point is identical from both computations. The response in the unstable branches is however much more complex for bifurcation buckling than for limit point buckling. In the configuration leading to bifurcation, the unstable branch of the limit point analysis are preserved, but two additional load paths are obtained. This is especially evident in the case of interlayer slip (Fig. 5c). It should be emphasized here that in reality there are always imperfections in the load as well as in the structure, which means that the critical load of a perfect configuration with subsequent limit point bucking is never reached. The response shown in red is therefore purely hypothetical. In reality, however, a bifurcation problem always exists for the considered configuration of shallow arch 2. Figure 5b shows the distribution of w/l over x/l at the four discrete load levels A, B, C, D (indicated in Fig. 5a) in the case of bifurcation buckling. It can be seen that the deflection becomes asymmetric at the bifurcation point A. In contrast, at the opposite point B on the second stable response branch, far away from an unstable deformation branch, the deflection is virtually symmetrically distributed. At bifurcation point C, the deflection becomes asymmetric again. In the unstable branch at point D, the deflection is also asymmetric. For the interlayer slip (Fig. 5d), a strong deviation from the antimetric response is observed near the bifurcation points. Figure 6a proves that during bifurcation buckling in the unstable response branches between bifurcation points A and C, the overall normal force N is virtually constant. It can also be seen that, in contrast to shallow arch 1, there is a change from compressive to tensile normal force during the transition from A to B. The asymmetric distribution of the layerwise normal forces N 1 and N 2 over x/l at the onset of bifurcation buckling is apparent from Fig. 6b. In contrast to shallow arch 1, the plot of the external force p against the moment Ml/E J ∞ exhibits a loop, as can be seen in Fig. 6c. In addition, Fig. 6d shows that a sign change in the normalized moment Ml/E J ∞ along the beam axis occurs at the onset of bifurcation buckling.

Modified shallow arch 2 subjected to distributed load
In a third example (referred to as modified shallow arch 2), the slip modulus is five times larger than in shallow arch 2 (i.e., K s = 5 · 10 8 N/m 2 and further αl = 23.28), otherwise all parameters and the load distribution remain unchanged. Figure 7 shows that also in this structure instability is determined by bifurcation buckling This result again proves that capturing the flexibility of the interlayer is essential for predicting the nonlinear stable and unstable response. Moreover, the response behavior is much more complex than for shallow arch 2, i.e., the number of unstable equilibrium paths increases. The external force p versus normal force N /E A e diagram has two loops in this case, as can be seen in Fig. 7c. Instead of two, this structure has four bifurcation points. However, only the two with the smaller critical loads are of practical importance because they indicate the transition from stable to unstable deformation branches. As can be seen, in the associated unstable response branches the normal force is again virtually constant (at N /E A e = −1.02 · 10 −4 and N /E A e = −2.99 · 10 −4 , respectively).

Shallow arch 2 subjected to single force
In the fourth example, shallow arch 2 is considered again, but this time it is loaded by a single force. In a first computation, this single load is applied exactly in the center, where δ is the Dirac delta function. In a second computation, this force is shifted by 1% to the left, to capture the influence of a small load imperfection. In Fig. 8, the normalized single force p = P 0 l 2 /E J ∞ is plotted over the response of this structure: in Fig. 8a the normalized deflection w/l at x/l = 0.5, in Fig. 8b the normalized interlayer slip u/l at the right support, in Fig. 8c the normalized total normal force N /E A e , and in Fig. 8d the normalized total bending moment Ml/E J ∞ at x/l = 0.5. As observed, also in this example instability is governed by bifurcation buckling. Comparison with Fig. 5 shows that a single force is more critical than a distributed load, since the critical loads are smaller with p cr = 2.29 and p cr = −0.45, respectively (compared to p cr = 3.61 and p cr = −0.83, respectively, for a distributed load). The response exhibits more unstable deformation branches than for a distributed load. Four bifurcation points are observed as in modified shallow arch 2, as can be seen particularly from Fig. 8c.

Effect of the boundary conditions
In the next example, the influence of the boundary conditions on the critical loads and the nonlinear response is studied. For this investigation, the support conditions in shallow arch 2 are changed at the left end (x/l = 0). In the first case, the left support is hard-hinged (sh − hh), in the second case the left support is clamped (sh − cl). The right end (x/l = 1) remains a soft-hinged support in both cases. The load is uniformly distributed p(x) = p 0 over x. Figure 9a shows the normalized external load p = p 0 l 3 /E J ∞ versus normalized deflection w/l at x/l = 0.5 for the two arches with modified support (green line: case 1 hh − sh; blue line: case 2 cl − sh) and, in addition, the response of the soft-hinged structure on both sides (sh − sh) subjected to an imperfectly distributed load already shown in Fig. 5a. As can be seen from this figure, at the end of the primary stable deformation path the critical load is the same for the two arches with modified boundary conditions at the left support, with p cr = 3.44, and slightly smaller than for bifurcation buckling of the symmetrically supported structure (where p cr = 3.61). While the unstable deformation branch of the structure with boundary conditions hh − sh subsequently follows the course of the load-deflection path of case sh − sh (bifurcation buckling), the shallow arch clamped on the left end shows a completely different deformation pattern. This is also reflected in the critical load on the remote stable deformation branch, which is p cr = −0.17 for the structure with boundary conditions hh − sh and p cr = 2.07 for the boundary conditions cl − sh.
In Fig. 9b, p is plotted versus the normalized interlayer slip u/l at the right end. In the first response path (loading), the interlayer slip is virtually the same for all three structures with different boundary conditions. The deviations become significant only in the unstable deformation branch. It can be seen that the deformation pattern has only one loop when the left boundary conditions are modified. This loop is much smaller in area for boundary conditions hh − sh than for the other two cases. The largest interlayer slip at the member end x/l = 1 in a stable branch occurs for the boundary conditions cl − sh immediately before snap-back to the primary stable deformation branch.

Variation of the rise
Finally, the influence of the riseŵ 0 of the shallow arch in the form of a sinusoidal half-wave on the critical load p cr in the first load path is investigated. To this end, in Fig. 10 the normalized critical load p cr = p 0cr l 3 /E J ∞ on the primary stable deformation branch is plotted against the normalized riseŵ 0 /l in the range from 0 to 0.05 for variations of the shallow arch 2 under uniform load p 0 with interlayer stiffness K s = 10 8 N/m 2 (αl = 10.41), K s = 5 · 10 8 N/m 2 (αl = 23.28), and with rigidly bonded layers (αl = ∞), respectively. The computations were performed with both symmetric loading and imperfectly distributed loading according to Eq. (41). This figure shows that, as expected, the critical load increases with increasing normalized riseŵ 0 /l as well as with increasing interlayer stiffness. The minimum value ofŵ 0 /l above which a stability problem occurs can also be read. For small values ofŵ 0 /l after this limit, the structure becomes unstable due to limit point buckling. For largerŵ 0 /l values, instability is governed by bifurcation buckling, which can be seen from the fact that the corresponding critical load is smaller than that for hypothetical limit point buckling. The limit values ofŵ 0 /l at the transition from limit point buckling to bifurcation buckling again depends on the interlayer stiffness.

Summary and conclusions
Since in-plane limit point buckling (snap-through) and in-plane bifurcation buckling of shallow arches composed of flexibly bonded layers have not been studied in the literature to date, this paper introduced a beam theory for the stability and instability analysis of shallow arches with interlayer slip. In a first example, a comparative computation on a much more elaborate finite shell element model based on a plane stress state showed that the presented theory is both efficient and suitable to predict both the stable and unstable nonlinear response of slender shallow arches with interlayer slip. In the example problems, deflection, interlayer slip, total and layer-wise normal forces, and total bending moment and layer-wise moments are shown graphically as governing response variables.
In order to detect a possible bifurcation buckling of a symmetric shallow arch in a numerical computation, it is essential to impose a small imperfection (see, e.g. Kiss [19]). In the present examples, this was realized by applying a load with extremely small eccentricity. For the more slender arches considered, it was shown that otherwise limit point buckling instead of the actual unstable behavior (bifurcation buckling) is predicted and consequently the critical load is overestimated. Compared to limit point buckling, bifurcation buckling is accompanied by additional unstable branches in the load-deformation diagrams.
With increasing slenderness of the shallow arch as well as increasing interlayer stiffness, the number of unstable branches in the load-deformation diagrams increases. Likewise, the influence of the load distribution on the response was shown. In several examples, multiple unstable deformation branches were predicted, depending on various parameters such as the stiffness of the interlayer and the loading conditions. Selected examples have demonstrated that neglecting the flexibility of the interlayer leads to a considerable overestimation of the critical load.
The presented theory can be used to perform parametric studies due to its simplicity and numerical efficiency. Thus, a tool is available to comprehensively investigate and better understand the stability problem of transversely loaded shallow arches with interlayer slip. In future studies, inertia terms should be added to the theory in order to investigate the dynamic nature of the snap-through phenomenon.
Funding Open access funding provided by University of Innsbruck and Medical University of Innsbruck.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Conflict of interests
The authors have no relevant financial or non-financial interests to disclose.