## Abstract

A buckled beam with shallow rise under lateral constraint is considered. The initial rise results from a prescribed end displacement. The beam is modeled as inextensible, and analytical solutions of the equilibria are obtained from a constrained energy minimization problem. For simplicity, the results are derived for the archetypal beam with pinned ends. It is found that there are an infinite number of zero lateral-load equilibria, each corresponding to an Euler buckling mode. A numerical model is used to verify the accuracy of the model and also to explore the effects of extensibility.

## 1 Introduction

Shallow arches are known to exhibit complex non-linear behavior when subjected to transverse loads. These relatively simple structures—which underlie more complicated civil, mechanical, aerospace, and biological systems—have a rich history of research (e.g., Refs. [16]), yet continue to yield new and exciting results for a variety of boundary and loading conditions (e.g., Refs. [715]). Recent interest has stemmed in part from “smart” devices, energy harvesting, and morphing systems that harness buckling and bistable post-buckling behavior [1618].

In a recent study [19], a relationship was discovered between the Euler buckling loads of an initially straight beam and the lateral force–displacement curve of a buckled beam under lateral loading. In particular, it was shown that the number of unstable equilibria could be determined purely based on the buckled beam’s geometry. The number of unstable equilibria was related to the number of “loops” in the lateral load–deflection relation, quantified by the number of equilibria where the transverse force was zero (called “zero-transverse-load crossings”). This number was determined from the “squash” load at which the extensible beam would flatten to conform to the prescribed end displacement. This paper seeks to prove that the number of loops, and correspondingly the number of zero-transverse-load crossings, is infinite for the inextensible beam. Furthermore, the previous study [19] relied on numerical (finite element method (FEM)) results to demonstrate the hypothesized relationship, whereas this paper derives the result analytically.

The following work considers an elastic, inextensible, initially straight beam that is buckled to give a pre-stressed, shallow arch [11], as opposed to a curved, stress-free arch [20]. It is hoped that this work can provide insight on the behavior post-buckled structures, approximate preliminary design equations, and motivate extensions for other boundary conditions. The analysis of the buckled beam is formulated in Sec. 2. Using a constrained energy minimization approach, the lateral load–deflection relation is derived, leading to a proof of the infinite number of zero lateral-load equilibria for the archetypal case of the pin-ended, inextensible, buckled beam. Some results are presented and discussed in Sec. 3, followed by a comparison to the extensible case in Sec. 4. Section 5 gives some concluding remarks.

## 2 Mathematical Formulation

Consider the initially straight, inextensible beam of length L, shown in Fig. 1(a). It is assumed homogeneous and uniform with flexural rigidity EI. When subjected to a prescribed end displacement of Δ0 (Fig. 1(b)), a shallow lateral deflection y0(x) results, with an associated rise of H at midspan. Now a lateral constraint is imposed at a point $x¯$ along the length of the buckled beam such that the lateral deflection $y(x¯)=A$ (Fig. 1(c)), where y(x) is the general deflected shape and A is the amplitude of the deflected shape at $x¯$. This constraint is taken here to be “knife-edge” (Fig. 1(c)), where the rotation at $x¯$ is not hindered.

Fig. 1
Fig. 1
Close modal
For a general (unconstrained) deflected shape y(x) and end displacement Δ, the arc length S of the beam midline is given by
$S[y(x)]=∫0L−Δ1+[y′(x)]2dx≈∫0L−Δ{1+12[y′(x)]2}dx=L−Δ+12∫0L−Δ[y′(x)]2dx$
(1)
If the midline axial deformation of the beam is defined as δLS (hence shortening is positive), then the end displacement Δ, which is assumed small, is given by
$Δ[y(x)]=δ+12∫0L−Δ[y′(x)]2dx≈δ+12∫0L[y′(x)]2dx$
(2)
For an inextensional beam, there is no midline shortening (i.e., δ = 0), and the end displacement simplifies to [6]
$Δ[y(x)]=12∫0L[y′(x)]2dx$
(3)
The total elastic strain energy U in the beam for the deflected shape y(x) can be expressed as [6]
$U[y(x)]=12EI∫0L[y″(x)]2{1−[y′(x)]2}−1/2dx≈12EI∫0L[y″(x)]2dx$
(4)
Since the beam buckles to a shallow rise H (≪L) under small end displacement Δ (≪L), only the leading terms of the expansions of the integrands in Eqs. (1) and (4) need to be retained as a suitable approximation.
Stable equilibria exist when the minima of the strain energy occur, and this minimization is constrained by both the end displacement and the lateral constraint, which can be expressed as follows:
${miny(x)U[y(x)]=12EI∫0L[y″(x)]2dxsubjectto12∫0L[y′(x)]2dx=Δ0y(x¯)=A$
(5)
where Δ0 = Δ[y0(x)] = const., as the beam is inextensible and the support will be held stationary after the initial end displacement Δ0 is applied. The second constraint equation (i.e., $y(x¯)=A$) is for the knife-edge lateral constraint assumed in this study, but other constraint conditions could alternatively be considered by modifying this constraint equation (e.g., $y(x¯)⩽A$ for unilateral constraint [20,21]) or appending additional equations (e.g., $y′(x¯)=0$ for an infinitesimally long clamped constraint).
The general expression for the deflection curve can be expressed by the series
$y(x)=∑n=1∞Qnψn(x)$
(6)
where ψn(x) = a shape function consistent with the boundary conditions and Qn = the associated coefficient. Substituting this form of y(x) into Eq. (5), the minimization problem becomes
${minQ1,Q2,…U[Q1,Q2,…]=12EI∫0L(∑n=1∞Qnψn″(x))2dxsubjectto12∫0L(∑n=1∞Qnψn′(x))2dx=Δ0∑n=1∞Qnψn(x¯)=A$
(7)
Introducing the notation
$kij=EI∫0Lψi″(x)ψj″(x)dx$
(8)
$aij=∫0Lψi′(x)ψj′(x)dx$
(9)
Equation (7) can be rewritten as
${minQ1,Q2,…U[Q1,Q2,…]=12∑i=1∞∑j=1∞kijQiQjsubjectto12∑i=1∞∑j=1∞aijQiQj=Δ0∑j=1∞Qjψj(x¯)=A$
(10)
A solution to the minimization problem (10) is sought using a multiplier method, wherein multipliers μ and λ are used to enforce the constraints. The end-displacement equation is multiplied by μ and the lateral-constraint equation by λ. These equations are then added to the strain energy equation to provide an augmented strain energy function, UaUa(Q1, Q2, …, λ, μ), which will be maximized over the multipliers and minimized over the coefficients, as described by
$minQ1,Q2,…maxλ,μUa=12∑i=1∞∑j=1∞kijQiQj−μ(Δ0−12∑i=1∞∑j=1∞aijQiQj)+λ(A−∑j=1∞Qjψj(x¯))$
(11)
The necessary conditions for optimality are given by
$∂Ua∂Qn=0:∑j=1∞knjQj+(∑j=1∞anjQj)μ−ψn(x¯)λ=0,n=1,2,…$
(12a)
$∂Ua∂μ=0:Δ0−12∑i=1∞∑j=1∞aijQiQj=0$
(12b)
$∂Ua∂λ=0:A−∑j=1∞Qjψj(x¯)=0$
(12c)
The roots of these ∞ + 2 equations in terms of the coefficients (Q1, Q2, …) and multipliers (μ, λ) are sought. For the following discussion, it is convenient to write Eq. (12a) in a matrix form as follows:
$[k11+a11μk12+a12μ⋯k21+a21μk22+a22μ⋯⋮⋮⋱]{Q1Q2⋮}={ψ1(x¯)ψ2(x¯)⋮}λ$
(13)
A test of stability is performed by analyzing the eigenvalues of the bordered Hessian matrix:
$H(Ua)=[00−∑j=1∞a1jQj−∑j=1∞a2jQj⋯00−ψ1(x¯)−ψ2(x¯)⋯−∑j=1∞a1jQj−ψ1(x¯)k11+a11μk12+a12μ⋯−∑j=1∞a2jQj−ψ2(x¯)k21+a21μk22+a22μ⋯⋮⋮⋮⋮⋱]$
(14)
A stable system will have two positive eigenvalues corresponding to the dual variables (i.e., multipliers μ and λ), with the rest of the eigenvalues negative.

Finally, note that the physical interpretations of μ and λ are the axial load required to induce Δ = Δ0 and the lateral load required to induce $y(x¯)=A$, respectively. Because of the way the constraints were augmented to U (see Eq. (11)), the axial load μ is positive in tension and the lateral load λ is positive upward (i.e., opposite as drawn in Fig. 1(c)).

### 2.1 The Pin-Ended Beam.

Hereinafter the case of a pin-ended beam is considered. For this case, the following expansion is assumed:
$y(x)=∑n=1∞QnsinnπxL$
(15)
for which kij = aij = 0 for ij, knn = EI()4/2L3, and ann = ()2/2L. Hence, Eq. (13) simplifies to
$[k11+a11μ0⋯0k22+a22μ⋱⋮⋱⋱]{Q1Q2⋮}={ψ1(x¯)ψ2(x¯)⋮}λ$
(16)
and the bordered Hessian (14) simplifies as follows:
$H(Ua)=[00−a11Q1−a22Q2⋯00−ψ1(x¯)−ψ2(x¯)⋯−a11Q1−ψ1(x¯)k11+a11μ0⋯−a22Q2−ψ2(x¯)0k22+a22μ⋱⋮⋮⋮⋱⋱]$
(17)
Moreover, the end-displacement constraint (12b) simplifies as follows:
$Δ0=12∑n=1∞annQn2$
(18)
To establish the initial end displacement, the unconstrained case will be considered first.

#### 2.1.1 No Lateral Constraint.

In the absence of a lateral constraint, the multiplier λ ≡ 0 and the right side of Eq. (16) goes to zero. The solution to the homogeneous system of equations is found from the eigenvalue problem, i.e., the determinant of the matrix is zero, giving the characteristic equation in terms of μ:
$∏n=1∞(knn+annμ)=0$
(19)
The roots of the characteristic equation are given by
$μn=−knnann=−(nπ)4EI/2L3(nπ)2/2L=−(nπ)2EIL2,n=1,2,…$
(20)
which are recognized as the Euler buckling loads for a pin-ended beam. The eigenvector corresponding to the nth buckling load is Qi = 0 for in, with Qn arbitrary. The first buckling mode will be assumed by the beam, as it corresponds to the lowest load; therefore, this will be used as the initial condition to define Δ0 as follows:
$Δ0≡12a11H2=π2H24L$
(21)
where it is recalled that H is the initial rise. Substituting this expressions for Δ0 into Eq. (18), the following equation is recovered:
$12∑n=1∞(nπ)22LQn2=π2H24L⇒∑n=1∞n2Qn2=H2$
(22)
This equation represents a hyper-ellipsoid in the space of Qn (n = 1, 2, …), which satisfies the end-displacement condition.

#### 2.1.2 Solution With Lateral Constraint at $x¯$⁠.

The buckling loads given by Eq. (20) are for the system without lateral constraint. Now, a solution is sought for the case in which the beam is laterally constrained at $x¯$. Solving Eq. (16) for the coefficients in terms of μ and λ gives
$Qn=ψn(x¯)knn+annμλ,n=1,2,…$
(23)
Substituting this expression into Eqs. (22) and (12c), noting that the multiplier λ can be factored out of the summations, gives
$(∑n=1∞n2[ψn(x¯)]2(knn+annμ)2)λ2=H2$
(24)
and
$(∑n=1∞[ψn(x¯)]2knn+annμ)λ=A$
(25)
Solving Eq. (25) for λ, and substituting this expression into Eq. (24) gives
$(∑n=1∞n2[ψn(x¯)]2(knn+annμ)2)(∑n=1∞[ψn(x¯)]2knn+annμ)−2A2=H2$
(26)
Rearranging terms in this equation gives
$(A/H)2=(∑n=1∞n2[ψn(x¯)]2(knn+annμ)2)−1(∑n=1∞[ψn(x¯)]2knn+annμ)2$
(27)
Equation (27) provides a straightforward relationship to determine the deflection A at $x¯$ for a given μ. Then, the corresponding lateral constraint force λ can be found from Eq. (25), which can then be used in Eq. (23) along with μ to find the coefficients Qn (n = 1, 2, …). It should be noted that, for a given μ, there are two solutions that differ only in the signs of A, λ, and Qn. The stability of these solutions (μ, ± A, ± λ, ± Q1, ± Q2, …) is determined by looking at the number of positive eigenvalues of the bordered Hessian (17). As noted previously, a stable system will have two positive eigenvalues corresponding to the dual variables (i.e., multipliers μ and λ), with the rest of the eigenvalues negative.

A zero lateral-load crossing (ZLLC) corresponds to an equilibrium where the Lagrange multiplier λ = 0 [19], which is consistent with the case without lateral constraint (Sec. 2.1.1). The number of ZLLCs can be observed by inspection of Eq. (16). By setting λ = 0, Eq. (16) simplifies to a homogeneous system of equations, from which the characteristic equation (Eq. (19)) is recovered whose roots are the Euler buckling loads given by Eq. (20). Hence, this proves that there are an infinite number of ZLLCs for the inextensible pin-ended buckled beam with shallow rise, which correspond to the Euler buckling loads. In the case of a different choice of shape functions (potentially non-orthogonal) or boundary conditions other than pinned ends, the ZLLCs are found from Eq. (13) by setting λ = 0 and solving the characteristic equation
$|k11+a11μk12+a12μ⋯k21+a21μk22+a22μ⋯⋮⋮⋱|=0$
(28)
There are an infinite number of unique roots (i.e., ZLLCs) to this equation, which correspond to the buckling loads. As will be shown in the following section, there are in general four stable ZLLCs—first and second buckling modes and their mirror images—and the rest unstable, but under certain conditions (i.e., six or more terms retained in approximation of y(x) and a constraint close the beam’s end) the number of stable ZLLCs decreases to two.

## 3 Results

In this section, some results for the pin-ended beam are shown to demonstrate the theory developed in Sec 2. The results are presented primarily in the form of normalized lateral load–deflection (λA) relations. Additionally, the influence of increasing the number of modes retained in the expansion [Eq. (15)] on the lateral load λ and axial load μ is presented.

### 3.1 Effect of Lateral Constraint Location $x¯$

Figure 2 presents the normalized lateral load–deflection (λA) relation for varying lateral constraint location $x¯=$ (a) 0.499L, (b) 0.49L, (c) 0.45L, and (d) 0.33L, assuming a three-term approximation of y(x). The lateral load λ is normalized by the critical load [6]
$PC=3π4EIH/2L3$
(29)
and the deflection A at $x=x¯$ is normalized by the initial rise H. Detached equilibria are observed (distinguished by line color), corresponding to the sign of A, λ, and Qn, as discussed before.
Fig. 2
Fig. 2
Close modal

Representative equilibrium configurations are shown for seven equilibria (i–vii) along one of the detached equilibrium paths (gray) for $x¯=0.45L$ (Fig. 2(c)). The initial equilibrium configuration (ii) corresponds to the first Euler buckling mode, for which λ = 0 and A/H is slightly less than 1 because A is measured at $x¯=0.45L$, whereas H is at x = L/2. Configuration (i) corresponds to a slight increase in A accompanied by a sharp increase in load, which continues to infinity. Alternatively, if A is decreased starting from configuration (ii), the lateral load λ decreases (positive stiffness) until a minimum is reached at configuration (iii)—a limit point, corresponding to snap-through under load-control—at which point the load begins to increase (negative stiffness) until reaching zero at configuration (iv). Configuration (iv) corresponds to the second Euler buckling mode, with A slightly less than zero (for similar reason as before); configuration (iv) also corresponds to the point of snap-through under unilateral displacement control [9,11,21]. Continuing to decrease A, the load increases up to a critical point—corresponding to configuration (v)—at which stability is lost. This critical point coincides with a vertical tangent. Following the equilibrium path further along the unstable branch, with increasing A, the load increases and then decreases, returning to zero load at configuration (vi). Configuration (vi) corresponds to the third Euler buckling mode. Finally, if A is further increased, passing through configuration (vii), the load blows up to negative infinity.

It is worth noting that, experimentally, the loss of stability at configuration (v) would propagate a jump—sometimes called “snap-back” [9,11]—to a neighboring stable equilibrium, i.e., the black line in Fig. 2(c). Similar equilibrium configurations exist on this path (black line), but are not shown because they are simply mirror images of those shown for the other path (gray line). As discussed in Sec. 3.2, retaining more terms in the expansion of y(x) results in more “loops” in the equilibrium paths and more ZLLCs, but the two detached equilibrium paths will never connect for the case of the inextensible beam.

With only three terms, there are six equilibria with λ = 0 (i.e., ZLLCs): $A/H≃$ ±1, 0 (double root), and ±1/3. The ZLLCs at A/H ≃ ±1 and 0 (e.g., configurations (ii) and (iv), respectively, in Fig. 2(c)) correspond to stable solutions, whereas the ZLLCs at A/H ≃ ±1/3 (e.g., configuration (vi) in Fig. 2(c)) are unstable. As the lateral constraint moves away from midspan ($x¯=L/2$), the stable ZLLCs near A/H ≃ 0 (i.e., close to the horizontal, the chord connecting the arch ends) progressively move away from A/H = 0, corresponding to a more gradual transition from positive to negative stiffnesses. This positive-to-negative transition is not pronounced for large eccentricities, such as $x¯=0.33L$ (Fig. 2(d)). Furthermore, the critical points at which stability is lost (e.g., configuration (v) in Fig. 2(c)) progressively move away from the nearby stable branch, resulting in larger jumps when stability is lost. There is a corresponding increase in the critical lateral load at which this jump occurs, discussed in Sec. 3.1.2.

It should be noted that not all parts of the equilibrium paths, such as those in Fig. 2, will be experimentally observable, at least not readily. Of course, at each value of λ, the equilibria occur where the total potential has zero-valued gradient, but these equilibria may be valleys (stable), hilltops (unstable), or saddles (unstable) in the potential energy hypersurface. The unstable equilibria will both be difficult to “find” experimentally [22], and it would also be difficult to hold the structures in such configurations under the unavoidable perturbations and imperfections that occur in experiments. Furthermore, many regions of the potential energy landscape are “steep” (i.e., this is a stiff system), meaning that divergence away from unstable equilibria after a small perturbation may be quite violent. Another subtle issue is that displacement control provides a constraint that restricts certain routes through the potential energy landscape; this may stabilize some, but not all unstable equilibria. The simplest example of this is that displacement control at or near the midpoint stabilizes the primary bifurcation path, but not paths with higher spatial wave numbers. Additional constraints could potentially be applied to stabilize more of the equilibrium path [23]; however, without feedback control these constraints may inadvertently change the force pattern to be something other than what is intended.

#### 3.1.2 Critical Equilibria.

From the preceding discussion, three critical equilibria can be identified from the lateral load–deflection curves: (1) a limit point resulting in snap-through under load-control at a critical lateral load, denoted $λLP*$; (2) a snap-through instability under unilateral displacement-control at a critical displacement, denoted $AUL*$; and (3) a snap-back instability at a critical lateral load, denoted $λSB*$. These three critical conditions are depicted in Fig. 3 and are recognized as corresponding to configurations of the type (iii), (iv), and (v), respectively, illustrated in Fig. 2(c). All three of these critical equilibria correspond to the onset of snapping behavior depending on the experimental loading protocol used; namely, types (1) and (2) for load-control and unilateral displacement-control, respectively, whereas type (3) occurs for both load- and (bilateral) displacement-control tests.

Fig. 3
Fig. 3
Close modal

Figure 4 portrays the effect of the lateral constraint location $x¯$ on these three critical values when three terms are used in the approximation of y(x). The values reported are for equilibrium paths corresponding to the gray lines in Fig. 2. For a midpoint lateral constraint ($x¯=L/2$), the critical values are $λLP*$ = PC, $AUL*$ = 0, and $λSB*$ = −PC. As the lateral constraint moves away from midspan, the critical lateral loads $λLP*$ and $λSB*$ both increase, whereas the critical displacement $AUL*$ decreases as previously discussed. The critical limit-point lateral load $λLP*$ increases to about −PC/2 near $x¯=L/3$ and 2L/3, before decreasing and going to negative infinite as $x¯→0$ and L; therefore, it is physically the easiest (i.e., requiring the smallest absolute load) to initiate snap-through under load-control for a constraint applied close to the third points of the buckled beam. The critical displacement $AUL*$ decreases to about −H/2 at $x¯=L/4$ and 3L/4, before returning to 0 at $x¯=0$ and L; therefore, under unilateral (and bilateral) displacement-control, the constraint location point is able to pass below the horizontal ($AUL*$ < 0) before the arch snaps downward for all constraint locations except $x¯=0,$L/2, and L. The critical snap-back lateral load $λSB*$ peaks at $x¯=L/3$ (i.e., the node line of the third buckling mode), reaching nearly 2.5 times PC; this large load corresponds to a bifurcation, which can be seen in Fig. 2(d) for $x¯=0.33L$. The lowest $λSB*$ is half of PC for a constraint location close to $x¯=L/9$ and 8L/9, but this minimum critical lateral load is dependent on the number of terms used to approximate y(x) (see Sec. 3.2).

Fig. 4
Fig. 4
Close modal

#### 3.1.3 Axial Load μ Relation.

To better understand these results, Fig. 5 portrays the axial load μ versus deflection A, lateral load λ, and coefficients Qi for the same cases considered in Fig. 2. The axial load μ is normalized by the first Euler buckling load for a pin-ended beam (i.e., Pe = π2EI/L2); note that, as defined, μ is negative for compression, so the sign is reversed in Fig. 5. For the nearly midspan load [Fig. 5(a)], the deflection A, lateral load λ, and coefficients Qi (i = 1, 2, 3) gradually vary up to an axial load of −4π2EI/L2 (i.e., the second Euler buckling load), at which point the deflection A and lateral load λ abruptly change signs. This transition is realized through the growth of the second coefficient (Q2), while the first and third coefficients (Q1 and Q3, respectively) go to zero. This is more clearly shown by the coefficient for $x¯=0.45L$ in Fig. 5(c), where the coefficients are distinguished by the line weight—Q1 (thickest) to Q3 (thinnest)—and are normalized by $∑j=13|Qj|$. A similar effect is observed at an axial load of −9π2EI/L2 (i.e., the third Euler buckling load), where only the third coefficient (Q3) is present. Finally, it should be noted that the ZLLCs directly correspond to the Euler buckling loads [i.e., −μ/(π2EI/L2) = 1, 4, and 9].

Fig. 5
Fig. 5
Close modal

### 3.2 Effect of Number of Terms Used to Represent y(x).

Figure 2(b) showed the normalized lateral load–deflection relation for a three-term approximation of the deflected shape y(x) for a constraint location $x¯=0.49L$. Figure 6 shows the same relation, but for varying numbers of terms used in the approximation: (a) 5, (b) 8, (c) 11, and (d) 22 terms. As more terms are retained in the expansion of y(x), additional “loops” in the load–deflection curves are observed. A similar effect has been observed when the rise of an extensible arch is increased [19], resulting in progressively more ZLLCs. For the system considered in Ref. [19], however, the number of loops was determined to be limited based on the axial load required to “flatten” the beam. In the formulation presented here, progressively more “loops” and ZLLCs are recovered as more terms are retained in the expansion, with the relationship holding up to an infinite number of ZLLCs for the infinite series (Eq. (15)). The ZLLCs for an approximately midspan constraint ($x¯≈0.5L$) correspond to A/H ≃ ±1, 0, ±1/3, 0, ±1/5, …, where each 0 is a double root.

Fig. 6
Fig. 6
Close modal

It is worth noting that the additional loops all correspond to unstable branches of the load–deflection relations, yet the stable branches remain largely unchanged with more terms being retained. As a result, the limit-point lateral load $λLP*$ (not shown) was unaffected (within numerical precision) by the number of terms retained in the approximation of y(x), whereas the snap-back lateral load $λSB*$ (Fig. 7(a)) and critical displacement $AUL*$ (Fig. 7(b)) showed varying degrees of sensitivity to the number of terms retained. In particular, for constraint locations near midspan ($0.45L⩽x¯⩽0.55L$), the snap-back lateral load $λSB*$ (Fig. 7(a)) is relatively insensitive to the number of terms. This critical load is more sensitive to the number of terms when the constraint is located near the beam’s third points ($x¯≃L/3$ and 2L/3) or ends ($0 and $5L/6⩽x¯) where the critical load peaks. In the case of the former (third points), retaining more terms reduces the peak critical load, with less effect when more than four terms are retained. In the case of the latter (close to beam ends), the effect of the number of terms is more dramatic, resulting in a wide range of critical loads, including both positive and negative loads. For three, four, or five terms, $λSB*$ has a distinct minimum that approaches $x¯=0$ and L as more terms are retained, but these minimum critical loads never go below zero in these cases. Conversely, for six or more terms, $λSB*$ crosses zero for load locations close to 0.05L and 0.95L and asymptotes to negative infinity as $x¯→0$ and L. If $λSB*$ is less than zero, then the stable branch of the lateral load–deflection path (e.g., solid lines in Figs. 2 and 6) would not reach zero load before going unstable, and hence instability under unilateral displacement-control testing would occur at a lateral load of $λSB*$ instead of at the critical displacement $AUL*$ which corresponds to a ZLLC (i.e., λ = 0). Therefore, two of the ZLLCs that are typically stable would be unstable in this case. This behavior was experimentally observed in Ref. [21]; the portion of the stable branch between configurations (iv) and (v) in Fig. 2(c), called the “secondary branch” in Ref. [21], could not be found experimentally for constraint locations near the supports ($x¯<0.152L$ and $x¯>0.848L$). Furthermore, the critical displacement $AUL*$ as defined in Sec. 3.1.2 is not well defined in this case (i.e., it is unstable), indicated in Fig. 7(b) by the termination of the lines. Otherwise, the critical displacement $AUL*$ is independent of the number of terms retained in the approximation of y(x). Note that, in Fig. 7(b), there are in fact seven lines, but they fall on top of each other, obstructing the cases with fewer terms.

Fig. 7
Fig. 7
Close modal

Figure 8(a) shows the full range of the lateral load–deflection relation for the case of an 11-term approximation of y(x), corresponding to Fig. 6(c). It can be observed that the lateral load λ increases with each successive loop that is added. This is illustrated in Fig. 8(b); as the axial load μ increases, the amplitude of the lateral load λ increases with each successive ZLLC. Note that in Fig. 8(b) the square root of the normalized axial load is plotted to stretch the μ-axis. In doing so, the ZLLCs occur at integer values (n = 1, 2, …), corresponding to the Euler buckling loads n2π2EI/L2. At each buckling load, the corresponding coefficient Qn is present with all other coefficient going to zero.

Fig. 8
Fig. 8
Close modal

## 4 Effect of Extensibility

The formulation in Sec. 2 relies on the assumption of inextensibility, which is applied through the constraint in Eq. (5). This section explores the validity of this assumption through a comparison with an in-house co-rotational finite element model (cFEM). The details of the co-rotational model formulation can be found in Refs. [24,25]. The cFEM additionally provides a validation of the small deformation assumptions used in the formulation as it accurately captures not only the effect of extensibility, but also arbitrarily large global deformations, and is limited only by the fact that it assumes small strains.

For shallow extensible beams with axial rigidity AE, the axial thrust can be closely approximated by
$P=AELδ$
(30)
where δ is the midline axial deformation (shortening) of the beam as introduced in Sec. 2. This result for the axial load is a constitutive expression, which cannot be used for the inextensible case. Extensibility limits the maximum static compressive force that can be sustained for a given Δ0 value. As noted in Ref. [19], this maximum axial force corresponds with the case where the beam is flattened between the two supports. From Fig. 1, it can be seen that this squash force is (identically) given by
$Psq=AELΔ0$
(31)
However, for the extensible case, the expression for Δ0 in Eq. (21) must be modified as follows. First, substituting the post-buckled non-transversely loaded deformed shape y0(x) into Eq. (2) makes the result specific to the midline shortening value corresponding to the initial post-buckled configuration, which will be denoted by δ0. As hinted by Fig. 1(b), this buckled shape can be approximated by a half sine wave, i.e.,
$y0(x)=Hsin(πxL−Δ0)$
(32)
In fact, as this is the first buckling mode, it is the identical deformed shape in the small deformation limit. The result of the substitution of y0(x) into Eq. (2) is
$Δ0=δ0+π2H24(L−Δ0)≈δ0+π2H24L$
(33)
which collapses to Eq. (21) for the inextensible case (i.e., δ = δ0 = 0). Finally, making the approximation that the initial axial load in the post-buckled (but not transversely loaded) configuration in Fig. 1(b) will be very close to the first Euler buckling load giving
$AELδ0≈π2EIL2$
(34)
This approximation relies on the fact that the classical supercritical pitchfork post-buckling path of axially loaded struts exhibits only a small increase in axial load for small transverse deformations.
The equations presented above can be used to approximate the thrust force for any deformed shape, and agree with the result for thrust force in  Appendix and Appendix B of Ref. [11] (though different labels and normalization are used). Furthermore, combining Eqs. (31), (32), and (34) also provides a simple but accurate estimate of the squash load as a function of the buckled rise and beam properties:
$Psq=(1+H24r2)π2EIL2≡n~2π2EIL2$
(35)
where $r=I/A$ is the radius of gyration and $n~$ is defined as the term in parentheses. Hence, the squash load is dependent on H/r, often called the dimensionless rise [26] (though a factor of 1/2 is also often applied). Additionally, the value $n~$, or more appropriately, the integer floor($n~$), gives the highest buckling mode number that can be achieved when the beam is flattened between supports. It should be noted that this squash limit is discussed in Ref. [19], but it is not explicitly derived as a function of the dimensionless rise. With increasing interest in exploiting bistable beams for smart structures and energy harvesting applications, Eq. (35), and ultimately extensions thereof to clamped boundary conditions, may prove useful in preliminary design steps. For example, the rise needed to transition from symmetric limit point to asymmetric bifurcation snap-through can be readily estimated without the need for the FEM simulation.

The effect of the compression limit is demonstrated in Fig. 9, which shows a comparison of the results for extensible (cFEM) and inextensible (analytical) cases. The different results shown for extensible beams correspond to scenarios with identical normalized forces and deformations, but differing axial rigidities (see the simulation parameters in the figure caption). This is readily done in an FEM simulation, as the cross-sectional area can be increased without changing the bending stiffness. Figure 9(a) shows that the normalized lateral load–deflection curves approach the inextensible case as the axial stiffness is increased. The axial load plots shown in Fig. 9(b) for successively (axially) stiffer beams are perhaps more informative. They show how the squash limit (denoted by horizontal dash-dotted lines), from Eq. (35), serves as a “short circuit” that allows the equilibrium path to begin winding down in axial force toward the snapped-through configuration. The superimposed inextensible results, on the other hand, are unable to pass through this flattened configuration and hence two disconnected paths both extend as vertical asymptotes. While adding additional modes to the inextensible solution increases the number of loops that appear, the two paths will never connect.

Fig. 9
Fig. 9
Close modal

The peak load (which is a ZLLC) corresponds with the completely flattened beam. This peak load is not (necessarily) identically equal to a buckling load, unless $n~$ in Eq. (35) is coincidentally an integer value. The approximate results for the four cases in Fig. 9(b), for example, are $n~=6.08,7.81,11.00$, and 15.52, respectively. It should be noted that the curves in Fig. 9(b) come from the cFEM, which is exact, even for large deformation. Hence, the approximate expression for the squash load from Eq. (35) is surprisingly accurate at predicting both the squash load and the number of loops that will occur in the load–deflection curve. Additionally, $n~$ passing through increasing integer values corresponds with bifurcations resulting in additional loops in the force deformation curves. The transition through $n~=2$, in particular, provides a useful estimate of the system parameters at which a bifurcation buckling may occur as opposed to limit point buckling, though the bifurcation path initially appears (generates) on the down-sloping part of the equilibrium paths.

## 5 Concluding Remarks

An elastic, inextensible, initially straight beam that is buckled to give a pre-stressed, shallow arch was analyzed. The buckled beam, subject to both a prescribed end displacement and a lateral constraint, was modeled analytically using a constrained energy approach, and equilibria were characterized. These analytical results prove that, for inextensible beams, there are an infinite number of equilibrium configurations that require zero lateral load. Curves of the lateral constraint force versus the displacement at the constraint location were presented, which exhibited an increasing number of loops and zero lateral-load crossings with the number of terms retained in the expression for the deflected shape. Furthermore, it was shown that the zero lateral-load crossings corresponded to the Euler buckling modes for the flat beam. These results were verified numerically by adjusting the axial rigidity of a buckled beam modeled within a co-rotational finite element formulation. It was shown that the primary (stable) equilibrium path was insensitive to the number of terms retained in the analytical model, as well as the numerical model.  Appendix gives an extension of the energy minimization formulation presented in Sec. 2 to include extensibility. Other possible extensions of this work are the generalization to arbitrary boundary conditions and curved beams (i.e., stress-free arches).

## Acknowledgment

This material is based upon work supported by the National Science Foundation under Grant Nos. NSF-CMMI-1663376 and NSF-CMMI-1943917 and the Air Force Office of Scientific Research Grant No. FA9550-19-1-0031. This support is greatly appreciated.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The authors attest that all data for this study are included in the paper.

### Appendix: Relationship Between the Extensible and Inextensible Cases

For the case of an extensible beam, the strain energy due to axial deformation δ is given by
$Uaxial[y(x)]=12AELδ2$
(A1)
where AE/L is the axial stiffness of the beam. By Eq. (2) and noting that Δ[y(x)] = Δ0 = const., as the pin supports are not permitted to move, the axial strain energy can be rewritten as follows:
$Uaxial[y(x)]=12AEL(Δ[y(x)]−12∫0L[y′(x)]2dx)2=12AEL(Δ0−12∫0L[y′(x)]2dx)2$
(A2)
The strain energy due to bending is given by Eq. (4), and the total strain energy of the beam is then given by
$U[y(x)]=12EI∫0L[y″(x)]2dx+12AEL(Δ0−12∫0L[y′(x)]2dx)2$
(A3)
Proceeding in the usual way with the calculus of variations:
$0=∂U[y(x)+εw(x)]∂ε|ε=0=∫0L[EIv⁗(x)+AEL{Δ0−12∫0L[y′(x)]2dx}y″(x)]w(x)dx+b.c.′s.$
(A4)
where ɛw(x) is the variation of y(x). Noting that w(x) is arbitrary, the following equation is recovered:
$EIv⁗(x)+P*y″(x)=0$
(A5)
where
$P*=AEL{Δ0−12∫0L[y′(x)]2dx}$
(A6)
By Eqs. (31) and (35), this expression for P* can be written as follows:
$P*=(1+H24r2)π2EIL2−AE2L∫0L[y′(x)]2dx$
(A7)
which is equivalent to Eq. (B.1) in Ref. [11]. Furthermore, P* is recognized as the axial thrust given in Eq. (30).

## References

1.
Fung
,
Y. C.
, and
Kaplan
,
A.
, “
Buckling of Low Arches or Curved Beams of Small Curvature
,”
Technical Note number 2840
,
,
1952
.
2.
Hoff
,
N. J.
, and
Bruce
,
V. G.
,
1953
, “
Dynamic Analysis of the Buckling of Laterally Loaded Flat Arches
,”
J. Math. Phys.
,
32
(
1–4
), pp.
276
288
. 10.1002/sapm1953321276
3.
Lock
,
M. H.
,
1966
, “
Snapping of a Shallow Sinusoidal Arch Under a Step Pressure Load
,”
AIAA J.
,
4
(
7
), pp.
1249
1256
. 10.2514/3.3656
4.
Walker
,
A. C.
,
1969
, “
A Non-linear Finite Element Analysis of Shallow Circular Arches
,”
Int. J. Solids. Struct.
,
5
(
2
), pp.
97
107
. 10.1016/0020-7683(69)90023-7
5.
Plaut
,
R. H.
,
1979
, “
Influence of Load Position on the Stability of Shallow Arches
,”
J. Appl. Math. Phys. (ZAMP)
,
30
, pp.
548
552
. 10.1007/BF01588902
6.
Thompson
,
J. M. T.
, and
Hunt
,
G. W.
,
1983
, “
On the Buckling and Imperfection-Sensitivity of Arches With and Without Prestress
,”
Int. J. Solids. Struct.
,
19
(
5
), pp.
445
459
. 10.1016/0020-7683(83)90055-0
7.
Chen
,
J.-S.
,
Ro
,
W.-C.
, and
Lin
,
J.-S.
,
2009
, “
Exact Static and Dynamic Critical Loads of a Sinusoidal Arch Under a Point Force at the Midpoint
,”
Int. J. Non-Linear Mech.
,
44
(
1
), pp.
66
70
. 10.1016/j.ijnonlinmec.2008.08.006
8.
Virgin
,
L. N.
,
Wiebe
,
R.
,
Spottswood
,
S. M.
, and
Eason
,
T. G.
,
2014
, “
Sensitivity in the Structural Behavior of Shallow Arches
,”
Int. J. Non-Linear Mech.
,
58
, pp.
212
221
. 10.1016/j.ijnonlinmec.2013.10.003
9.
Plaut
,
R. H.
,
2015
, “
Snap-through of Arches and Buckled Beams Under Unilateral Displacement Control
,”
Int. J. Solids. Struct.
,
63
, pp.
109
113
. 10.1016/j.ijsolstr.2015.02.044
10.
Zhou
,
Y.
,
Chang
,
W.
, and
Stanciulescu
,
I.
,
2015
, “
Non-Linear Stability and Remote Unconnected Equilibria of Shallow Arches with Asymmetric Geometric Imperfections
,”
Int. J. Non-Linear Mech.
,
77
, pp.
1
11
. 10.1016/j.ijnonlinmec.2015.06.015
11.
Plaut
,
R. H.
, and
Virgin
,
L. N.
,
2017
, “
Snap-Through Under Unilateral Displacement Control With Constant Velocity
,”
Int. J. Non-Linear Mech.
,
94
, pp.
292
299
. 10.1016/j.ijnonlinmec.2017.01.015
12.
Sano
,
T. G.
, and
,
H.
,
2018
, “
Snap-Buckling in Asymmetrically Constrained Elastic Strips
,”
Phys. Rev. E
,
97
(
1
), p.
013002
. 10.1103/PhysRevE.97.013002
13.
Gao
,
R.
,
Li
,
M.
,
Wang
,
Q.
,
Zhao
,
J.
, and
Liu
,
S.
,
2018
, “
A Novel Design Method of Bistable Structures With Required Snap-Through Properties
,”
Sens. Actuators., A.
,
272
, pp.
295
300
. 10.1016/j.sna.2017.12.019
14.
Yan
,
S.-T.
,
Shen
,
X.
,
Chen
,
Z.
, and
Jin
,
Z.
,
2018
, “
Collapse Behavior of Non-uniform Shallow Arch Under a Concentrated Load for Fixed and Pinned Boundary Conditions
,”
Int. J. Mech. Sci.
,
137
, pp.
46
67
. 10.1016/j.ijmecsci.2018.01.005
15.
Zhao
,
J.
,
Zhang
,
J.
,
Wang
,
K. W.
,
Cheng
,
K.
,
Wang
,
H.
,
Huang
,
Y.
, and
Liu
,
P.
,
2020
, “
On the Nonlinear Snap-Through of Arch-Shaped Clamped-Clamped Bistable Beams
,”
ASME J. Appl. Mech.
,
87
(
2
), p.
024502
. 10.1115/1.4045593
16.
Harne
,
R. L.
, and
Wang
,
K. W.
,
2013
, “
A Review of the Recent Research on Vibration Energy Harvesting Via Bistable Systems
,”
Smart Mater. Struct.
,
22
(
2
), p.
023001
. 10.1088/0964-1726/22/2/023001
17.
Hu
,
N.
, and
Burgueño
,
R.
,
2015
, “
Buckling-Induced Smart Applications: Recent Advances and Trends
,”
Smart Mater. Struct.
,
24
(
6
), p.
063001
. 10.1088/0964-1726/24/6/063001
18.
Zhang
,
Z.
,
Li
,
Y.
,
Yu
,
X.
,
Li
,
X.
,
Wu
,
H.
,
Wu
,
H.
,
Jiang
,
S.
, and
Chai
,
G.
,
2019
, “
Bistable Morphing Composite Structures: A Review
,”
Thin-Walled Struct.
,
142
, pp.
74
97
. 10.1016/j.tws.2019.04.040
19.
Nistor
,
M.
,
Wiebe
,
R.
, and
Stanciulescu
,
I.
,
2017
, “
Relationship Between Euler Buckling and Unstable Equilibria of Buckled Beams
,”
Int. J. Non-Linear Mech.
,
95
, pp.
151
161
. 10.1016/j.ijnonlinmec.2017.06.016
20.
Plaut
,
R. H.
,
2015
, “
Snap-Through of Shallow Extensible Arches Under Unilateral Displacement Control
,”
ASME J. Appl. Mech.
,
82
(
9
), p.
094503
. 10.1115/1.4030741
21.
Harvey, Jr.
,
P. S.
, and
Virgin
,
L. N.
,
2015
, “
Coexisting Equilibria and Stability of a Shallow Arch: Unilateral Displacement-Control Experiments and Theory
,”
Int. J. Solids. Struct.
,
54
, pp.
1
11
. 10.1016/j.ijsolstr.2014.11.016
22.
van Iderstein
,
T.
, and
Wiebe
,
R.
,
2019
, “Experimental Path Following of Unstable Static Equilibria for Snap-Through Buckling,”
Nonlinear Dynamics
, Vol.
1
,
Kerschen
,
G.
, ed.,
Springer International Publishing
,
Cham
, pp.
17
22
.
23.
Neville
,
R. M.
,
Groh
,
R. M. J.
,
Pirrera
,
A.
, and
Schenk
,
M.
,
2020
, “
Beyond the Fold: Experimentally Traversing Limit Points in Nonlinear Structures
,”
Proce. R. Soc. A: Math., Phys. Eng. Sci.
,
476
(
2233
), p.
20190576
. 10.1098/rspa.2019.0576
24.
Masashi
,
I.
,
1994
, “
Effects of Coordinate System on the Accuracy of Corotational Formulation for Bernoulli-Euler’s Beam
,”
Int. J. Solids. Struct.
,
31
(
20
), pp.
2793
2806
. 10.1016/0020-7683(94)90069-8
25.
Battini
,
J.-M.
,
2002
, “
Co-Rotational Beam Elements in Instability Problems
,” Ph.D. Thesis,
KTH Royal Institute of Technology
,
Stockholm
.
26.
Zhou
,
Y.
,
Yi
,
Z.
, and
Stanciulescu
,
I.
,
2019
, “
Nonlinear Buckling and Postbuckling of Shallow Arches With Vertical Elastic Supports
,”
ASME J. Appl. Mech.
,
86
(
6
), p.
061001
. 10.1115/1.4042572