Dynamic Behavior of Bistable Shallow Arches: From Intrawell to Chaotic Motion

Historically, arches were considered static structural elements, designed to resist deformation and minimize deflection [1]. Recently, an interest in multistable mechanical metamaterials inverted this perception and found ways to take advantage of elastic instabilities [2–4]. In particular, the “snap-through” phenomenon has spurred promising new applications for the bistable arch as a building block in multistable mechanical metamaterials [5–8]. Snap-through is the sudden transition between stable configurations resulting in large displacement and potential energy release [9]. Given that the snap-through instability is scale independent, bistable arches are functional for numerous applications across length scales ranging from the micro to macro scale. On the micro-scale, micro-electromechanical systems utilize bistable arches in electrical systems as microactuators [10], relays, switches [11], and memory cells [12]. Macro-scale bistable arches are used for energy production [13], energy dissipation [14], and passive sensing devices [15]. The switch-like response of bistable elements has promoted the study of structures, such as morphing aircraft structures, that can reconfigure themselves to optimize performances while adjusting physical properties [16,17]. Logic gates and mechanical computers have shown promising avenues to leverage bistability to encode bits of information and perform logic operations [18,19].

Previous studies on bistable arches have imposed constant velocity displacement to study the dynamic response and approximate the growth rate of the phenomenon [20]. The snapping dynamics of an elastic arch have been shown to manifest a critical slowing down on the edge of the snap-through transition [21]. Rotational boundary actuation and asymptotic analysis have been used to characterize the dynamics of arches near a shape transition [22]. Differences in the snapping behavior between elastically and plastically deformed arches have also been investigated [23].

In terms of actuation mechanisms, magnetic forces [24–26], photosensitive materials [27,28], heat activation [29,30], and electric loads [31] have all been leveraged to trigger the snap-through instability in bistable arches. The bistable arch has also been used to control fluid flow while being actuated from the fluid itself [32–34]. Each activation method is tailored to a specific application but in all cases, energy is added to the system to overcome the energy barrier between the energy wells and trigger a snap-through transition. Typically, quasi-static loading conditions are assumed to determine the critical load or displacement needed to overcome the energy barrier to initiate the shape transition [20]. However, under oscillatory loads, the necessary conditions for snap-through will rely both on the amplitude and the frequency of the external force. With a single-degree-of-freedom (SDOF) representation, the bistable arch can be modeled as a Duffing oscillator [35], a prototypical example of a bistable oscillator which has been extensively studied to predict the onset of chaos [36], the nonlinear hardening and softening behavior [37], and the criteria for interwell motion (i.e., motion between both wells) [38].

Here, we explore the dynamic behavior under oscillatory loads of two varieties of bistable shallow arches: the elastically deformed arch (i.e., buckled beam) and the plastically deformed arch with initial curvature. First, we analytically derive an expression describing the relationship between the first natural frequency and the rise of the elastically deformed arch. We investigate both the linear and nonlinear frequency response for both varieties of bistable arches. Then, we classify the response regime using the potential energy function and interpret the motion with respect to the two energy wells. We define a critical energy criterion needed for vacillation between the two energy wells. Finally, we explore the parameter space for a single-mode approximation of an arch as well as a three-mode approximation.

These findings have the potential to set the foundation for the design of next generation multistable metamaterials with energy harvesting and reconfiguration capabilities.

We consider two classes of bistable shallow arches: elastically deformed buckled beams and plastically deformed beams. The motion of both shallow arches is described by Eq. (1), but differ in their expression for end shortening and initial curvature:

• Elastically Deformed Arch. For an elastically deformed shallow arch (a buckled beam), the initial unstressed position of the midsurface is flat with zero initial curvature. An axial force is applied to buckle the beam into a deformed shape with rise b (Fig. 1(a)). Removing the external axial force would permit the arch to straighten back out to a flat beam. The end shortening ΔL can be found imposing a constant arc length between the undeformed and deformed beam configuration. Thus, the elastically deformed arch is described by

  $w0(x)=0$

  (3a)

  $ΔL≈1+(b2)2$

  (3b)
• Plastically Deformed Arch. A plastically deformed arch is manufactured by permanently deforming the arch with an initial curvature. We impose an initial curvature of a half sine wave with rise e. When the plastically deformed arch is removed from the supports, the initial curvature remains (Fig. 1(e)). There is no end shortening in the manufacturing process. The arch is manufactured to have the same length in its final configuration. The plastically deformed arch is expressed using

  $w0(x)=esinx$

  (4a)

  $ΔL=0$

  (4b)

To classify the dynamic behavior of the arch, we turn our attention to the potential energy of the structure. The potential energy, which is a scalar function of the mode shapes, governs the behavior of the unforced, undamped system. The stability of fixed points, energy wells, escape velocities, and domains of attraction are all dependent on the potential energy. Compared to attempts to classify the motion of the arch by tracking points along the length, the potential energy conveniently encompasses the entire shape of the arch. Furthermore, the energy function can easily be extended to higher dimensional models.

For N = 3, the inclusion of a₃ preserves the two stable equilibrium points with symmetric (elastic) or asymmetric (plastic) energy wells and two more unstable equilibrium points are generated. However, the associated energies of these new unstable equilibrium points are higher than the energy of the unstable equilibrium points generated by the second mode. Thus, the energy wells of the three-dimensional model are still defined by the energy barrier from the unstable equilibrium points found for N = 2 via Eqs. (24) and (26) for the elastic and plastic arches, respectively. The wells for the N = 3 system are visualized by plotting the boundaries of the wells in the a₁, a₂, a₃ coordinate system for the elastic (Fig. 1(d)) and plastic arches (Fig. 1(f)). It is seen that in the elastic case, the energy wells remain symmetric about the plane a₁ = 0, while for the plastic case, the energy wells are asymmetric and separated by the plane described by Eq. (20).

For any combination of a₁, a₂, …, a_N, we can determine which region of the energy landscape the structure is in. First, we calculate the potential energy for a given a using Eq. (11). If V(a) > V_cr, the structure is not in either well. Otherwise, the current configuration’s energy well can be located by comparing the component of a₁ to the unstable equilibrium point from N = 1 described by Eq. (14) for an elastic arch or Eq. (20) for a plastic arch.

To introduce oscillatory loads, we begin with a discussion of the natural frequencies of the shallow arch. Analytical formulations for the natural frequencies of arches have been derived for various boundary conditions [41,42]. As shown in Ref. [43], the first natural frequency of a hinged–hinged arch about an elastically buckled configuration is a function of the initial rise of the arch while higher frequencies are independent of the initial rise. Following this derivation, we develop closed-form expressions for the natural frequencies of the elastic deformed shallow arch. From these expressions, the initial rise can be chosen to investigate sub or super harmonics between the first and the n-th frequencies. In addition, we identify integer internal harmonics between pairs of n-th and j-th natural frequencies for hinged–hinged elastically deformed arch. Finally, we demonstrate the nonlinear effects of the bistable arch using frequency response curves.

To verify the analytical expression for the natural frequencies of the shallow arches, we implemented frequency sweep analyses using the software auto [44] based on a pseudo-arclength continuation technique. The frequency responses of the following three cases are investigated: (i) elastic arch about a stable equilibrium point, (ii) plastic arch about its lower energy stable equilibrium point, and (iii) plastic arch about its higher energy stable equilibrium point.

Using a seven degrees-of-freedom (DOFs) representation of the arch (N = 7), the systems of differential equations described by Eq. (6) are cast into first-order form and used for a continuation analysis. An initial rise b = 40 and e = 40 are chosen for the elastic and plastic arches, respectively. We excite the arches at η = 1/8 with a force amplitude Q = 0.1 and we assume light modal damping β = 0.1. Figures 2(a), 2(b), and 2(c) report the maximum amplitude modal response of the first five modes against the normalized excitation frequency Ω/ω₁ for the three cases. The analytical results for natural frequencies described by Eqs. (29) and (34) are plotted with vertical dashed lines. It is seen that peaks in the modal response from the elastically deformed arch precisely align with analytical values. The results for both cases of the plastically deformed arch include a slight shift in natural frequencies with respect to the elastic arch. When excited in the lower energy well, the plastically deformed arch presents a slight shift to the right of all the natural frequencies. Indeed, this suggests that the system is stiffer in this configuration, slightly increasing the natural frequencies. Similarly, a shift to the left is observed when exciting the plastic arch about its higher energy well, suggesting a less stiff system.

For sufficiently large forcing amplitudes Q, nonlinear phenomena can be observed in the frequency response of the shallow arches. For a small initial rise b = 5 (elastic) and e = 5 (plastic), we determine the dynamic response of each case for increasing forcing amplitudes (Q = 0.25, Q = 0.5, and Q = 0.75). Using auto’s pseudo-arc length continuation method, we observe a nonlinear softening behavior in the elastic arch case (Fig. 3(a)) as well as in the plastically deformed arch when excited both about its lower (Fig. 3(b)) and higher (Fig. 3(c)) energy wells. Despite the same initial rise, different amplitude responses are found between cases. The plastic arch excited about its lower energy well has the smallest amplitude response, consistent with the linear frequency response (Fig. 2(b)), suggesting the system requires more energy to oscillate around this well. Concurrently, when excited about the higher energy well, the plastic arch has the largest nonlinear response, demonstrating that the system requires less energy for oscillating around this well as also shown for linear vibrations (Fig. 2(c)). For each force amplitude Q, the nonlinear response of the elastic arch has vibration amplitudes in between the two previous plastic arches cases. In all the cases, an increase in the forcing amplitude results in an increased nonlinear effect, such as the nonlinear softening behavior and a larger range of unstable solutions. Arches with larger initial rises also exhibit nonlinear phenomena. An elastic arch with rise b = 25 subjected to an external force (Q = 0.50 and Q = 1.25) located at the midpoint of the arch η = 0.5 presents similar nonlinearities (Fig. 3(d)). Since the arch is excited about its midpoint, only the odd number modes are activated. The frequency response plot demonstrates a 1 : 2 internal resonance (ω₁ ≅ 2ω₃) which is found for Ω ≅ 0.5 ω₁.

We first introduce an energy-based criterion to characterize the motion of the arch over ranges of forcing amplitude Q and frequency Ω. In the first approximation, we explore the design space leveraging the Duffing oscillator as the archetypal SDOF model of the bistable arch. Then, we investigate the parameter space for a three-DOF model of the bistable elastic and plastic arches.

We classify the response of the structure into one of five regimes (Fig. 4): intrawell motion, switching motion, reverting motion, periodic vacillation, and aperiodic vacillation [46]. Intrawell motion is constrained to a single well over the entire time period N_C (intrawell in Fig. 4). If the motion ever escapes the initial well, we examine the behavior in N_S. If during N_S the structure is constrained to a single well, then the motion is either “switching” or “reverting.” Switching motion is classified as motion constrained to the well opposite the initial well. In this case, the structure “switched” from the initial energy well to the other well (switching in Fig. 4). Reverting motion is constrained to the original well over N_S. Reverting is different from intrawell motion since for reverting over the period N_T, the structure must have visited the other well at least once. In this case, although the opposite well was visited, the structure reverted to its original position (reverting in Fig. 4). The motion is described as vacillating when the structure is not constrained to a single well over the steady-state response N_S. Vacillating motion can be one of two types: periodic vacillation (periodic in Fig. 4) or aperiodic (chaotic) vacillation (aperiodic in Fig. 4).

The algorithm implemented to classify the behavior is represented as a flowchart in Fig. 4 together with schematics of the time response and corresponding positions in the energy landscapes. First, the system is integrated over the entire range N_C using the fourth-order Runge–Kutta time integration scheme. By examining the time response, we check if the motion is constrained to a single well. If true, the algorithm exits and the system is classified as intrawell. Otherwise, we check which wells are visited over the N_S period. If over the entire steady-state period the motion is confined to the well opposite the initial well, the motion is switching. If instead over the entire steady-state period the motion is confined to the initial well, the motion is described as reverting. If none of the above conditions are met, the motion is described as vacillating and the Lyapunov exponents are calculated to distinguish the motion between periodic and aperiodic [47]. If the largest Lyapunov exponent is less than zero, the motion is periodic. If instead the largest Lyapunov exponent is greater than zero, the system’s motion is chaotic (aperiodic).

We also track the energy of the system to test the strength of the condition for interwell motion described by the critical energy (Eq. (23)). Indeed, it is predicted that the system will have sufficient energy to move between wells when the total kinetic and potential energy exceeds the critical energy. For each trial, we tabulate the maximum total energy achieved in N_S period and compare it against the critical energy defined by Eq. (22). If the maximum combined kinetic and potential energy exceeds the critical energy, it is expected that the structure vacillates between wells. If the maximum combined kinetic and potential energy is less than the critical energy, the motion is expected to be confined to a well over the period N_S and either be intrawell, switching, or reverting motion. Since the critical energy is a necessary but not sufficient condition, we anticipate there to be cases with energy that exceeds the critical energy but does not escape the well.

We begin by exploring the design space of the one-dimensional model of an elastically deformed arch. Choosing an initial height of b = 25 and damping β = 1, the forcing frequency Ω/ω₁ and the magnitude of the force Q are varied from 0.2 to 1.2 and from 50 to 1500, respectively. For each pair of Ω and Q, the system is integrated in time for N_T and N_S cycles. Using the algorithm described in Fig. 4, the regimes of motion are determined, the Lyapunov exponents are calculated, values of the modal amplitude a₁ are sampled to generate bifurcation diagrams, and the maximum total energy is calculated over the N_S range.

The parameter space sampled points are colored according to the type of response (Fig. 5(a)). The majority of the plot is characterized by intrawell motion (plotted in yellow). As expected, low forcing magnitude result in intrawell motion. The critical force required to induce interwell motion varies with the forcing frequency as shown in the triangular structures descending from the top of the plot. These stalactite features representing interwell motion [49] are commonly called “Arnold’s” tongues and they have several practical implications. When large motion is favorable, such as in energy harvesting, the bistable structure should be operated within the bounds of these tongues to maximize motion [50]. An alternative use of the tongues is sensing critical frequencies for a constant force or critical forces for a constant frequency. The detection of the critical parameter could be aligned with the boundaries of the tongue, so interwell motion is achieved when the critical value is found [51].

Multiple analytical methods have been proposed for predicting the onset interwell motion [52,53,36]. To validate our findings, we consider the method proposed by Virgin et al. [48] which utilizes a harmonic balance method. Using a harmonic steady-state response, the first critical force that produces energy sufficient to overcome the energy barrier is determined. The corresponding force–energy relationship is then plotted in solid black (Fig. 5(a)). The analytical model presents very close agreement with the peak around the frequency Ω/ω₁ = 0.8. Since this model assumed that the structure’s motion has the same frequency as the forcing term, only the main tongue feature can be predicted.

As an additional method of visualizing the arch’s response, bifurcation diagrams are obtained using the method described by Parker and Chua [45]. For Q = 1000, we sweep over the same frequency range and sample the position of a₁ once every period of the forcing frequency is over the interval N_S. Plotting all values of a₁ for each frequency produces the bifurcation plot (Fig. 5(d)). With this method, period one motion appears as single points for a given value of Ω. Values of Ω that result in chaotic behavior appear on the bifurcation diagram as large sets of points.

Toward the right end of the plot in the region dominated by switching and reverting, there appears to be a region of period two motion from Ω/ω₁ > 1; however, this is misleading. In actuality, what is occurring is the motion changing from switching to reverting and vice versa for adjacent values of Ω. As this occurs, the bifurcation plot switches from having a cluster of points at a₁ ≈ 25 to having a cluster of points at a₁ ≈ −25. When this occurs for small changes in Ω, the bifurcation plot appears to be a period-2 motion. A similar phenomenon occurs around Ω/ω₁ ≈ 0.921.

Although only used to distinguish between periodic and aperiodic vacillation, the Lyapunov exponents were calculated for each pair of (Ω, Q) of the design space. The largest Lyapunov exponent λ₁ calculated for each trial is plotted within the parameter space (Fig. 5(b)) where the regions of periodic (λ₁ < 0) and aperiodic (λ₁ > 0) motion can be distinguished.

To visualize both Lyapunov exponents, a slice of the parameter space is taken for Q = 1000 while varying the value of Ω (Fig. 5(e)). It can be observed that the two Lyapunov exponents synchronously move with one another; as λ₁ increases, λ₂ decreases by approximately the same amount. In second-order non-autonomous systems, at least one Lyapunov exponent must be negative. Furthermore, for a chaotic attractor, since contraction must outweigh expansion when one Lyapunov exponent is positive, the other λ₂ < −λ₁ [45]. This relationship between the two Lyapunov exponents is seen as an inverse relationship between the two lines (Fig. 5(e)).

Over the final N_S period, the maximum of combined potential and kinetic energy was calculated in the parameter space (Fig. 5(c)). The maximum total energy has been compared with the critical energy: values below the critical energy predicted by Eq. (22) are plotted in blue and values which exceed the criteria are shaded yellow to red. Overall, our results show that critical energy is an excellent predictor of interwell motion. Although the energy criterion represents a necessary, but not sufficient condition, in only a few cases the maximum energy exceeds the critical energy but the motion remains intrawell.

A slice of the energy plot for Q = 1000 compares the critical energy plotted in a dashed line to the calculated energy (Fig. 5(f)). It is seen that all cases of interwell motion exceed this energy, and only a few cases of intrawell motion fall below the line. Reaffirming the necessary but not sufficient property of the critical energy, we can further examine why it is possible for the energy to exceed the critical energy, but motion not be interwell. Since the total energy is a combination of both kinetic and potential energy, the kinetic energy could be large enough to exceed the critical energy, but oscillating with very high frequency in a single well. For example, when the structure oscillates about a stable equilibrium with small amplitude motion but for a given set of frequencies, there might be enough kinetic energy to predict escape. The equivalent plots made for the two cases of the plastically deformed arch (i.e., high and low energy wells) are included in the appendix (Fig. 10) and they showcase a similar trend when compared to the elastic case.

Within the parameter space, the periodicity of the structure’s response can be further explored. For periodic motion, we can compare the ratio of the period of the structure to the period of the forcing frequency. The ratio between the forcing frequency and the structure’s frequency is plotted (Fig. 6). In most cases, the period of the structure is equal to the period of the forcing frequency, but higher period motion is also possible. The highest periodicity of response was found to be period-7 motion detected within the region dominated by the switching and reverting motion.

While a single-mode approximation of the arch presents strong connections with the Duffing oscillator, this model is not always sufficient to accurately describe the motion of the arch for large displacements. Experimental comparisons between analytical models and experimental results of the bistable arch suggested that a three-mode approximation can accurately describe the motion of the arch undergoing snap-through instability [23]. The analytical approach we introduced in this study can be easily applied to higher dimensional models of the bistable arch, and hence we report the investigation of the structure dynamical response for N = 3.

For N = 3, the parameter space is shown to be dominated by chaotic motion (Fig. 7(a)) and it significantly differs from the corresponding N = 1 parameter space (Fig. 5(a)). It is found that an approximate minimum value of Q ≈ 80 is required for the arch to escape the starting well. Above this threshold value of Q, the motion is almost entirely classified as interwell vacillation. Clearly defined tongues descending from the top of the plot observed in the SDOF approximation disappear in the three-DOF model of the arch. Very few cases of reverting or switching behavior are observed in the three-DOF design space.

The Lyapunov exponents have been determined within the parameter space and the largest exponent was used to verify if the interwell motion was periodic or chaotic (Fig. 7(b)). It is interesting to observe that for low forcing magnitudes, despite the motion being intrawell, the largest Lyapunov exponent is larger than zero in many cases. While the one-mode approximation saw very few points of intrawell motion with positive (or near zero) values of the largest Lyapunov exponent, the N = 3 model presents cases of chaotic intrawell motion. To visualize all six Lyapunov exponents, we consider the case of Q = 1000. Since modes are coupled via Eq. (6), when a single mode is chaotic, the chaos is transferred to the other two modes. Because of this coupling, it is expected that either all three modes will behave chaotically or all three modes will behave periodically. Indeed, this is observed in the bifurcation diagram for Q = 1000 (Fig. 7(d)) where either all three modes present a chaotic or a periodic behavior.

Lastly, we turn our attention to the critical energy criteria. As previously mentioned, the inclusion of the second mode introduces an unstable solution that lowers the energy barrier with respect to the N = 1 case. Thus, we expect the critical energy required to escape the wells to be much lower than the case of the single-mode approximation. This is consistent with the parameter space plot (Fig. 7(a)), where the majority of cases are characterized by interwell motion. Comparing the maximum combined energy from the parameter space with the critical energy shows that very few points are below the critical energy (Fig. 7(c)). Serving as a necessary but not sufficient criterion, the critical energy predicts if the structure is expected to escape an energy well, but due to such a small critical energy very few points do not escape. The energy plot for Q = 1000 (Fig. 7(f)) emphasizes the surplus of energy that the structure has with respect to the critical energy represented by the dashed horizontal line. Though higher dimensional systems may require both an excess of energy and a specific trajectory, the critical energy criteria can still provide powerful insight into the capacity of the system to demonstrate interwell motion.

In this work, we have analytically investigated the dynamics of elastically and plastically deformed bistable arches subjected to oscillatory loads. To characterize the response of the structure from intrawell to chaotic motion, we introduced a method based on the system’s location in the double-welled potential energy landscape. Thanks to this approach, a single scalar can conveniently encapsulate sufficient information to classify the motion for high-dimensional models of the system.

We first derived a closed-form solution for the natural frequencies of the structure and identified internal resonances. Tuning the natural frequencies based on the initial rise of the arches can represent a powerful tool in designing band gaps or frequency filters in arrays of arches. Moreover, we explored the nonlinear dynamic behavior of the system in which a softening behavior was observed for both the elastically and plastically deformed arches.

Finally, we thoroughly investigated the response of the structure within the parameter space by varying the force amplitude and excitation frequency. Using the energy wells to classify the motion, we identified five regimes of motion: intrawell, switching, reverting, interwell periodic, and aperiodic. Bifurcation diagrams, Lyapunov exponents, and critical energy plots have been introduced to identify the different regimes. These results were presented for both a single-mode and a three-mode approximation of the arch. In particular, the prediction of interwell motion and the onset of chaos are of significant relevance for engineering applications of bistable arches. To conclude, this study provides a deep understanding on the nonlinear dynamics of the elastically and plastically deformed shallow arches subjected to oscillatory loads. Thanks to its rich nonlinear behavior, the bistable arch is a promising building block for the design of next generation multistable mechanical metamaterials, passive control devices, and mechanical computers.

We would like to gratefully acknowledge Maryland High Power Computing at UMD.

• M.B. was supported by funding from the NDSEG Fellowship. E.T. was supported by funding from the University of Maryland Start-Up Package and the NSF CAREER Award (CMMI—2239841).

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

In Eq. (29), it was shown that the first natural frequency of the arch is a function of the initial rise of the arch. Higher natural frequencies are independent of the initial rise. The non-dimensional natural frequency is plotted against the non-dimensional rise of the arch (Fig. 8). From this plot, it can be seen that values for b can be chosen to induce integer ratios between the first and n-th frequencies. Integer ratios also exist between higher harmonics, but these harmonics are sparse and inherent to the structure of the arch.

The nonlinear phenomena of the plastically deformed arches were studied using pseudo-arclength continuation with the software auto. Choosing an initial rise of e = 25 for the plastically deformed arch and exciting the arch about its midpoint to activate only the odd-numbered modes, nonlinear phenomena can be seen for sufficiently large values of Q (Fig. 9).

The parameter space can be determined for plastically deformed arches with a SDOF approximation using the same method as described in the paper. For an initial rise e = 25, and performing time integration for N_T = 450 and N_S = 50, the method described in Fig. 4 is used to classify the motion into one of the five regimes (Fig. 10). The general structure of the stalactite-like tongues descending from the top of the graph remains consistent with what was seen for the elastically deformed arch.

The parameter space can be determined for plastically deformed arches with a three-DOF approximation using the same method as described in the paper. Choosing an initial rise e = 25, and performing time integration for N_T = 450 and N_S = 50, Fig. 11 was obtained using the algorithm shown in Fig. 4 to classify the motion into one of the five regimes. The parameter space is very similar to that seen in the elastically deformed cases, governed primarily by chaotic interwell motion.