The stability of circular rings is a classic and important topic that has attracted considerable interest among scholars and researchers. The problem may affect many fields with endless practical applications: structural, nuclear and offshore engineering, robotics, electromechanics, biomechanics, biomedical and biotechnological sciences, just to mention a few. The ways of approaching the stability of a ring can be very different depending on the perspective from which one decides to face the problem: Buckliphobia or Buckliphilia [1, 2].

From a Buckliphobia standpoint, the buckling of a ring is a failure mechanism to be avoided. Consider, for example, the loss of stability of collapsible tubes or pliable shells conveying pressurized and high-velocity fluid, such as gas, fuel, and blood [3–5]. The interaction between the internal flow and the tube, under certain critical conditions (e.g., negative transmural pressure of critical magnitude), can induce self-excited oscillations with finite amplitudes and, consequently, catastrophic collapse of the tube/shell [6]. There are a number of examples of such instability phenomena in everyday life: the aeroelastic stability (i.e., fluttering) of curved aircraft wings, for instance, is a topic of great interest in engineering since it induces painful limitations in the operational condition of the aircraft [7]; in the same area, one can mention the buckling and the postbuckling of chemical plant piping, offshore pipelines, or undersea vessels induced by hydroelastic interaction [8, 9] as they can affect the normal operation of plants or vessels.

More complex dynamical phenomena can be encountered in the physiological and medical fields due to the high flexibility of the biological tubes and the low Reynolds numbers of the biological flows: the high-frequency Korotokoff sound, i.e., a puzzling phenomenon induced by the resonant motion of arterial blood vessels, which trigger a dynamic transition from a buckled state to an expanding state (see Figure 1 and [10]); crackling sounds generated when a collapsed airway opens, e.g., the snoring produced by the oscillations of the soft palate, pharyngeal walls, epiglottis, and tongue during sleep when the pharyngeal muscles relax and block the airway [11]; the flow limitation during forced expiration induced by the interaction between airways and air [12]; the micturition induced by the interaction between urethras and urine [13]; and the pulse wave propagation in the cardiovascular system that can induce several cardiovascular pathologies, such as heart failure, stenosis [14], flow-induced deformation, and ultimate rupture of aneurysms; all come to mind as just a few of the more outstanding examples. A comparatively large share of the scientific and technological interest focuses on the instability of elastic rings subjected to external confinement [15, 16]. Curved film-substrate structures, in particular, have drawn special attention in the biomedical field because of the wide variety of folding phenomena that are observable in the human body: the formation of creases in the colon due to the growth of mucous membrane or the wrinkling patterns that can be observed on the surface of the bronchus lumen (Figure 2A), crumpled gut surfaces, the sulci and gyri of the brain, or the morphological evolution of tumor growth (Figure 2B), which is guided by the buckling mechanism of layered structures [17, 18].

On the contrary, from a Buckliphilia standpoint, the nonlinear ring dynamics can be regarded as a potential resource to be exploited. Biological phenomena are often based on the control ad tuning of buckling. Take, for instance, the wavy or wrinkled morphology of petals (Figure 2C), leaves, stems, or fungi that confers them unique physical, mechanical, and optical properties, such as the capability to actively control morphological transitions [19].

Mimicking such an aptitude, a new generation of tunable/adaptive micro-/macro-electronic devices, actuators, and controlled fabric folding can take advantage of the control of ring dynamics to attain a desirable feature. The issue is scientifically fascinating and technically important, especially in the field of soft matter mechanics and structures [20]. Soft materials undergo large and reversible deformations with a variety of instability modes, such as wrinkling and creasing modes. Such a behavior can be usefully exploited in a variety of functional smart applications by controlling, switching, and tuning the formation and disappearance of wrinkles and creases in such a way as to create custom-made topography, control chemical patterning, adhesion properties, shape morphing, and configuration, and so on [21, 22]. The range of innovative applications in the biomedical engineering field is equally wide. Venous valves used in vascular surgery and annuloplasty devices (Figure 2D), for instance, may benefit from the control of ring dynamic behavior providing a major socio-economic impact and/or potential for investment in addition to an improvement of their performances, efficiency, and reliability. In the optic field, intensive investigations are going on tunable resonators, photonic modulators, wavelength division multiplexing filters, and switches [23–26]. Taking advantage of the resonance characteristics of micro-rings, such devices offer many desirable characteristics and unique performance metrics such as ultra-compact physical size, ultra-high modulation speed, low power consumption, ultra-low modulation power, high energy efficiency, and so on. However, to date, the high susceptibility of micro-ring resonators/modulators to thermal fluctuations and manufacturing tolerances and their resonant nature make their practical applications still challenging. Another interesting field of application is voice control. During phonation, in fact, the position, the shape, and the internal surface tension of the glottis actively control the frequency and amplitude of the vocal cords oscillations [27].

From a purely theoretical point of view, there is extremely rich and interesting literature dealing with elastic ring instability problems (buckling in the plastic range poses much more complex and still debated problems [3–5, 28]). Most of these studies are devoted to the buckling and the initial postbuckling of rings under external pressure, while only a few articles are devoted to the complex nonlinear dynamics of elastic rings. In the present study, a unified approach will be pursued; that is, a set of the possible equations of motion will be derived from the static formulation of ring buckling.

The first contributions in this respect to the topic were given in the second half of the 19th century by Bresse [29], Boussinesq [30], and Levy [31] and were concerned with the static in-plane buckling of inextensible rings subjected to uniform external pressure. Significant advancements were successively provided by several authors in the 20th century [32–40]. The general agreement is that the critical value of the pressure magnitude, by neglecting the effect of the axial and the shear deformations, is

where EI is the bending stiffness of the ring wall (E is the Young's modulus of the material and I is the inertia of the cross section, respectively), R is the mean ring radius, and k is the buckling coefficient, which is dependent on the boundary and pressure-type loading conditions. In particular, for the hydrostatic pressure it is [37, 41–46]

so that for a two-wave deformation, n = 2,

holds true (note that n = 1 corresponds to a rigid-body motion). However, despite the long history of studies on the subject, the accuracy of the results, especially numerically, may still be an issue sometimes. Indeed, considering different constitutive assumptions, such as linear elastic material and Neo-Hookean materials, either compressible or incompressible, the case of the static, axially symmetric, deflection of a ring under internal and external pressure is analyzed in [47], and the obtained results show that also widely used commercial finite element programs may give inaccurate results in the case of linearly elastic rings undergoing large strains. Similar findings have been found and discussed in [48].

One of the first contributions to the study of the dynamic stability of the circular ring dates back to 1871 by Hoppe [49], who derived the frequency equation for the in-plane vibrations of a thin circular ring. Later, Love [50] and other authors proposed different linear and nonlinear approaches (e.g., [51–58]). In the midst of the 20th century, several authors analyzed the ring under uniform radial pressure periodically varying with time (e.g., [59–63]). More specifically, they generally analyzed the nonlinear parametric stability problem in the hypothesis of small perturbations from the initial motion, ring midline inextensibility, and negligible shear deformation by linearizing the problem in some way. For thick rings, some authors took into account a sort of shear deformability in the vibration analysis (e.g., [64–69]). With respect to the inextensibility, Cooley and Parkey made interesting considerations regarding the vibration of high-speed rotating elastic rings [70]. Other authors dealt with the stability of rings subjected to impulsive or step loads (e.g., [56, 71–77]). They essentially showed that the dynamic buckling of elastic rings is triggered by flexural modes governed by a Mathieu equation. The solution of this type of equation is expressed in terms of Mathieu functions, which are periodic in their argument only for certain so-called characteristic values of the ratio between the system and the loading frequencies. The solution may be stable, i.e., bounded, or unstable, i.e., diverging. The transition curves between stable and unstable domains in the plane of parameters are related to the characteristic values [78]. The dynamical post-buckling behavior of circular rings under external pressure shows an imperfection sensitivity for some types of loading different from the hydrostatic one [9, 39, 79–84]. In this regard, it is noteworthy that the omission of some nonlinear terms led Rehfield [79, 85] to the misleading prediction of imperfection sensitivity also for rings under hydrostatic pressure. Inconsistent results can also be found in [86] again because of erroneously neglected cubic and quartic terms in the energy of the loads [79].

As stated earlier, in the next sections, the present study will move from the basic formulation of the static buckling problem (Section 2), highlighting some underlying subtleties in mathematical modeling and will provide in this framework a set of very simple equations that nevertheless allow sketching the very complex and fascinating nonlinear behavior of the ring dynamics (Sections 3 and 4).

For the classical static buckling analysis, it is assumed that the ring wall is undeformable from the axial and the shear standpoint, and it can only experience flexural deformation. Here, it will be shown that the way in which axial inextensibility is imposed may have significant implications for the nonlinear dynamics of the problem.

With reference to the geometry of a buckled circular ring (see Figure 3), the axial strain ε_θθ = ε_θθ(θ, t) at the ring centerline, the cross-sectional rotation ϕ = ϕ(θ, t), and the change of curvature κ = κ(θ, t) of the ring, respectively, are [37, 38]

where R is the mean ring radius, v = v(θ, t) and w = w(θ, t) are the tangential and normal displacement to the ring wall midline (see Figure 4), and κ_ℓ = κ_ℓ(θ, t) is the linear change of curvature,

In all formulae, the subscript comma indicates differentiation with respect to the subsequent variables (in the case of the previous equations, the polar coordinate θ). No shear deformation is allowed for the ring wall.

The total potential energy of the ring is

where the strain energy V is the sum of the stretching and bending energies as follows

It is worth noticing that, to obtain the buckling pressure given by Equation (1), it is sufficient to assume κ ≈ κ_ℓ.

The load potential due to a uniform external load per unit length q = q(t) is

where ΔA describes the difference between the area enclosed in the ring centerline before and after the deformation (Figure 3), which can be expressed as follows [38]

The kinetic energy is

where μ stands for the mass per unit length of the ring wall.

The Lagrangian function is

and, by virtue of the Euler–Lagrange equation, the partial differential equations (PDEs) of motion can be written as follows

and

where F = F(θ, t) and G = G(θ, t) are dimensionless functions defined as follows

Substituting Equation (16) in Equation (4) and imposing the inextensibility of the ring (i.e., ε_θθ = 0) gives

which leads to four solutions in the form

In the following, we will make reference to Equations (18) that are independent of the radius R. It is interesting to recall that to evaluate the instability pressure given by Equation (1), classical static solutions (see, e.g., [37, 38]) assume that the inextensibility of the centerline is obtained by simply requiring that

which makes the linear part of Equation (4) null in virtue of the first of Equations (18). Indeed, the displacements corresponding to the first classic buckling pressure are given in Equation (1), i.e.,

which do not satisfy the equality, ε_θθ = 0 in toto.

Interestingly, should one use Equation (11) in the case of a hydrostatic load, the coefficient k of the buckling pressure would result 36/13, instead of the classic value given in Equation (3) [48]. To overcome this drawback, Equations (18) are both taken into consideration only to reduce Equation (11) [48] to

On the basis of the formulation introduced in Section 3.1 for the nonlinear dynamics of an elastic ring, results from several numerical simulations are collected in Figures 7, 8 in terms of Poincaré maps. To understand the dynamic behavior, several sets of initial conditions are considered. In particular, Equations (34) and (35) are integrated assuming that all the initial conditions, but φ_k(0) and φ.k(0), vanish, and we call the mode k the dominant degree of freedom. Therefore, it is set as follows:

It can be observed, and numerically verified, that such a choice implies that all the coupling terms in Equation (35) vanish. The results reported below are obtained by choosing the second mode (k = 2) as the dominant one. With reference to Figure 6, the motion is assumed to start at a configuration that is perturbed from the initial perfect round one (Figure 6A), and this allows the circular ring to deform dynamically (Figure 6B) also for maximum values of the forcing load, which would result below the statical critical pressure.

Provided that the forcing load is set as

the equation of motion of the parametrically excited dominant shape results from the Mathieu–Duffing type [78]. In the performed calculations, it is set Ω = 4.0, and the intensity δ ranges from 0.5 to 5.0 (δ = 1.0 being the static buckling load).

To depict the Poincaré sections with a sufficiently defined pattern, a sufficient number of points must be collected. Therefore, the simulation time spans over the interval 0 ≤ τ ≤ 2cπ/Ω, with c being a sufficiently large number of cycles of Q. In the following, c = 5 × 10³ is set if not otherwise stated.

Figures 7, 8 appear as layered tori, each torus being an invariant manifold fixed by the initial conditions. These stable tori occupy a precise, simply or multiply connected, region of the phase space. The outside regions are unstable with rapidly divergent solutions. For instance, in the case of a simply connected stable region, the unstable one can be found by choosing initial conditions such that Φ ≈ ±5/9, and Φ.=0. This is interesting since, according to the definitions given in Equation (3.1), it is 5/9=1/A2(2). Focusing on the stable solution of which the vast majority is quasiperiodic, it can be noticed that some periodic solutions can be found, too.

At δ = 0.5, a resonant solution is found for Φ = 0.4185636 and Φ.=0. Since the trajectory in the phase space goes through the Poincaré plane 10 times per cycle of the forcing load, this is a subharmonic of order 10, i.e., a solution of period 10 with respect to Q. Such a periodic solution is shown in Figure 5A and can contribute to explaining the behavior of a resonant torus. In the corresponding power spectrum, (Figure 5B), the first pulse appears at a frequency of exactly Ω/(10 × 2π). The 10 points are located inside islands, which can be detected with the proper refinement of the grid of the initial conditions, as can be seen in Figure 9A. The islands correspond to tori showing multiply connected sections (i.e., multiperiodic tori). In the case under consideration, since the cross sections of the tori are made of 10 islands, we have 10-periodic tori. The points corresponding to the periodic solutions are called elliptic points and are alternated with hyperbolic ones. These latter are actually saddle points (see Figure 9B). The lines connecting the hyperbolic points represent the boundaries of the islands and separate regions of tori with different periodicities.

At δ = 0.75, the resonant solution of period 10 can still be found closer to the unstable region, but the small islands are still not visible at the resolution of the graphics in Figure 7. More interestingly, a new stratum of islands appears. Differently from the case of the 10-periodic tori, here two families of 3-periodic tori (i.e., their cross sections are made of three islands), develop approximately two 3-periodic solutions. In the case of the 10-periodic tori, initial conditions that are symmetric with respect to the origin of the axes lead to the same trajectory. By contrast, in the case of the two families of 3-periodic tori, two sets of symmetric conditions lead to two different trajectories, one on a torus of the first family and the other on a torus of the second. The islands belonging to a family of tori are alternated with those belonging to the other.

At δ = 1.00, the islands of the external 10-periodic tori and the two internal 3-periodic tori are larger. At both strata, elliptic and hyperbolic points are found.

The section at δ = 1.25 shows that the islands of the 3-periodic tori are still increasing and, at the same time, radially pushed toward the unstable region by that of the internal (simply connected or single) tori. These islands progressively enlarge as the load intensity δ increases. As a consequence, the 10-periodic tori enter the unstable region, and their islands progressively lose their stability.

Such a phenomenon is confirmed at δ = 1.50 and it is even more evident at δ = 1.75 where 3-periodic tori are detached from the stable region and show distorted shapes. More interestingly, other small islands, apparently corresponding to 18-periodic tori, develop inside the 3-periodic tori islands.

As δ increases, the onset of new periodic solutions can be observed (see δ = 2 and δ = 2.25).

About δ = 2.25, the cross section of internal tori starts to shrink along the φ.2−axis. At δ = 2.5, a saddle at the origin of axes separates an external stratum of single tori from an internal family of 2-periodic tori. These latter tori will survive for a while (δ = 2.75 and more), even after the destruction of the external single tori that vanish in the unstable region.

However, for some values of δ above 3.25, where the islands are already distorted, no stable solutions can be found, at least in the region of parameters that have been considered (0.5 ≤ δ ≤ 5).

Overall, the variety of dynamic behaviors of the elastic rings yielded by the solution of Equations (34) and (35) summarizes some results of the Kolmogorov–Arnold–Moser (KAM) theorem and of the related theory. The KAM theorem proved for the first time that, in Hamiltonian systems, quasiperiodic motions are generic and, thus, it radically changed the previously accepted idea that their dynamics were close to being integrable [88]. Without entering details, if a Hamiltonian (as in the present case, where no source of dissipation is considered) dynamical system is subject to a weakly nonlinear perturbation, some of the invariant tori are deformed while others are destroyed. In the present study, the process has been driven by changing initial conditions since neighboring initial conditions can be thought of as mutual perturbations. The criterion by which this occurs is a condition of “quasiresonance” on the frequencies of the motions (commensurability), and the KAM theorem quantifies the conditions on the perturbations for this to happen. Indeed, the deformed tori have an even number of points, appearing in pairs of elliptic and hyperbolic fixed points, in common with the undeformed tori. At the elliptic fixed points, the dynamics is the same as the main system, giving rise to a self-similar structure. Hyperbolic points, being fixed saddle points, present trajectories entering or moving away from the fixed point, and they are both invariant sets, i.e., stable and unstable manifolds, respectively. Noteworthy, the considered system may exhibit chaotic behavior originating from hyperbolic points. The region of phase space where the chaos develops must be confined and that happens in the case an erratic motion arises between invariant tori, which play the role of barriers. Since often the erratic motions that have been found in the present study develop in places very close to the unstable region, additional investigations would be necessary to determine if any of the motions shown in Figures 7, 8 may be classified as chaotic, or they are just unstable solutions that will diverge after a sufficiently long time.

In the present study, starting from a brief resumè of the possible implications in various fields of physics, engineering, and biology of the buckling instability and nonlinear dynamics of circular elastic rings, attention has been focused on the fundamental hypotheses at the basis of most classic treatments of the problem.

In this respect, a simple model for analyzing the nonlinear dynamic behavior of elastic rings under external pulsating loads has been proposed, and a depiction of the possible evolution of the dynamical system has been discussed in detail.

In the literature, the behavior of rings under external pressure load has been recognized to be governed by the Mathieu equation, which can actually be used as a mathematical model in a number of problems in engineering and physics. To cite but a few, these are the vibrations in elliptic membranes, the motion of an inverted pendulum with the suspension point undergoing vertical vibration, the stability of floating bodies, the ion trap, and the elastic oscillations of a ferromagnetic material. The solution of the Mathieu equation is expressed in terms of Mathieu functions, which are periodic in their argument only for certain so-called characteristic values of the ratio between the system and the loading frequencies. The solution may be stable, i.e., bounded, or unstable, i.e., diverging, and the transition curves between stable and unstable domains in the plane of parameters are related to the characteristic values.

The simple modeling proposed in this study leads to a nonlinear version of the Mathieu equation called the Mathieu–Duffing equation due to the presence of cubic nonlinearity. Indeed, the Mathieu–Duffing equation is a relevant combination of the Mathieu and Duffing equations, the latter being perhaps the most well-known and celebrated oscillator with cubic nonlinearity.

The analysis of the obtained results has shown that even the simple problem under consideration can give origin to very complex dynamics in the phase space. Since the considered equation of motion does not contain any dissipative term, the resulting dynamics is Hamiltonian, with a wide prevalence of quasiperiodic responses. This means that invariant tori are characterized by two not commensurate frequencies and curves over these tori are not closed and fill the tori, which are non-resonant. In the case of commensurate frequencies, the responses are periodic and close to themselves.

The implications for many of the practical problems of elasticity described in Section 1 can be countless, depending on the specific application. It goes without saying that adaptations of the governing equations may be needed with respect to the type of loading, the material model, and so on.

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

IM and FG conceived the study. EB, IM, and FG contributed to the design of the study. IM developed the literature review related to the state of the art. EB developed the code and performed the numerical computations. EB and IM prepared the figures, edited the visualization, and drafted the manuscript. All authors contributed to the revision of the manuscript and read and approved the submitted version.