Dynamic buckling of an inextensible elastic ring: Linear and nonlinear analyses

TL;DR

This paper investigates the dynamic buckling behavior of elastic rings under pressure, combining experiments, numerical simulations, and asymptotic analysis to reveal how inertia influences shape evolution beyond static buckling predictions.

Contribution

It provides a comprehensive theoretical framework for understanding dynamic buckling of elastic rings, incorporating inertia effects through post-bifurcation analysis and numerical methods.

Findings

Inertia causes higher mode shapes during dynamic buckling.

Standard linear stability analysis fails to predict observed behaviors.

Theoretical results align with experimental observations.

Abstract

Slender elastic objects such as a column tend to buckle under loads. While static buckling is well understood as a bifurcation problem, the evolution of shapes during dynamic buckling is much harder to study. Elastic rings under normal pressure have emerged as a theoretical and experimental paradigm for the study of dynamic buckling with controlled loads. Experimentally, an elastic ring is placed within a soap film. When the film outside the ring is removed, surface tension pulls the ring inward, mimicking an external pressurization. Here we present a theoretical analysis of this process by performing a post-bifurcation analysis of an elastic ring under pressure. This analysis allows us to understand how inertia, material properties, and loading affect the observed shape. In particular, we combine direct numerical solutions with a post-bifurcation asymptotic analysis to show that inertia drives the system towards higher modes that cannot be selected in static buckling. Our theoretical results explain experimental observations that cannot be captured by a standard linear stability analysis.

Equations

$$\rho h\frac{{\partial }^{2}w}{\partial t^2}+B\frac{{\partial }^{4}w}{\partial x^4}+\mathcal{F}\frac{{\partial }^{2}w}{\partial x^2}=0,$$

$$w(-L/2)=w^{\prime \prime }(-L/2)=w(L/2)=w^{\prime \prime }(L/2)=0,$$

$$\rho h{\sigma }^2=\frac{(2n+1{)}^2{\pi }^2}{L^2}\left[\mathcal{F}-B\frac{(2n+1{)}^2{\pi }^2}{L^2}\right].$$

$$n\approx \frac{1}{2^{3/2}\pi }{\left(\frac{\mathcal{F}{L}^2}{B}\right)}^{1/2}-\frac{1}{2},$$

$$\frac{\partial \mathbf{F}}{\partial s}+p\mathbf{n}=\rho A_{\times }\frac{{\partial }^2\mathbf{r}}{\partial t^2},$$

$$\frac{\partial \mathbf{M}}{\partial s}+\mathbf{t}\times \mathbf{F}=0,$$

$$\frac{\partial x}{\partial s}=\cos \theta ,$$

$$\frac{\partial y}{\partial s}=\sin \theta ,$$

$$\frac{{\partial }^2\theta }{\partial s^2}=F_x\sin \theta -F_y\cos \theta ,$$

$$\frac{\partial F_x}{\partial s}-P\sin \theta =\frac{{\partial }^{2}x}{\partial t^2},$$

$$\frac{\partial F_y}{\partial s}+P\cos \theta =\frac{{\partial }^{2}y}{\partial t^2}.$$

$$P=\frac{p{a}^3}{B},$$

$$\theta (2\pi ,t)=\theta (0,t)+2\pi ,{\partial }_s\theta (2\pi ,t)={\partial }_s\theta (0,t),x(2\pi ,t)=x(0,t),$$

$$y(2\pi ,t)=y(0,t),F_x(2\pi ,t)=F_x(0,t),F_y(2\pi ,t)=F_y(0,t),$$

$$\theta (s_i,0)=s_i,x(s_i,0)=(1+{\mathcal{R}}_i)\sin (s_i),y(s_i,0)=-(1+{\mathcal{R}}_i)\cos (s_i).$$

$$\mathbf{M}={\mathbf{r}}^{\prime }\times {\mathbf{r}}^{\prime \prime }.$$

$${\mathbf{r}}^{\prime }\times ({\mathbf{r}}^{\prime \prime \prime }+\mathbf{F})=0,$$

$${\mathbf{r}}^{\prime \prime \prime }+\mathbf{F}=\lambda {\mathbf{r}}^{\prime }+\beta \mathbf{k},$$

$$\mathbf{F}=-{\mathbf{r}}^{\prime \prime \prime }+\lambda {\mathbf{r}}^{\prime }.$$

$$\ddot{\mathbf{r}}=-{\mathbf{r}}^{(4)}+(\lambda {\mathbf{r}}^{\prime }{)}^{\prime }+P\mathbf{n},{\mathbf{r}}^{\prime 2}=1,$$

$${\kappa }_0=1,$$

$${\mathbf{t}}_0=\cos s{\mathbf{e}}_x+\sin s{\mathbf{e}}_y,$$

$${\mathbf{n}}_0=-\sin s{\mathbf{e}}_x+\cos s{\mathbf{e}}_y.$$

$$\mathbf{r}={\mathbf{r}}_0+u{\mathbf{n}}_0+v{\mathbf{t}}_0.$$

$${\mathbf{r}}^{\prime }=\mathbf{t}=(\cos \varphi ){\mathbf{t}}_0+(\sin \varphi ){\mathbf{n}}_0=\left[1-u+v^{\prime }\right]{\mathbf{t}}_0+\left[u^{\prime }+v\right]{\mathbf{n}}_0.$$

$$\cos \varphi =1-u+v^{\prime },\sin \varphi =u^{\prime }+v.$$

$$\kappa =1+{\varphi }^{\prime }.$$

$$\ddot{\mathbf{t}}=-{\mathbf{t}}^{(4)}+(\lambda \mathbf{t}{)}^{\prime \prime }-P\kappa \mathbf{t}.$$

$$\cos \varphi =1-u+\frac{\partial v}{\partial s},$$

$$\sin \varphi =\frac{\partial u}{\partial s}+v,$$

$$\kappa =1+\frac{\partial \varphi }{\partial s},$$

$$\frac{{\partial }^2\varphi }{\partial t^2}=-\left(\frac{{\partial }^3\kappa }{\partial s^3}-6{\kappa }^2\frac{\partial \kappa }{\partial s}\right)+2\frac{\partial \lambda }{\partial s}\kappa +\lambda \frac{\partial \kappa }{\partial s}$$

$${\left(\frac{\partial \varphi }{\partial t}\right)}^2=\left[{\kappa }^4-3{\left(\frac{\partial \kappa }{\partial s}\right)}^2-4\kappa \frac{{\partial }^2\kappa }{\partial s^2}\right]-\frac{{\partial }^2\lambda }{\partial s^2}+\lambda {\kappa }^2+P\kappa .$$

$$\frac{\partial v_1}{\partial s}=u_1,$$

$${\varphi }_1=v_1+\frac{\partial u_1}{\partial s},$$

$${\kappa }_1=\frac{\partial {\varphi }_1}{\partial s},$$

$$\frac{{\partial }^2{\varphi }_1}{\partial t^2}=(5-P)\frac{\partial {\kappa }_1}{\partial s}+2\frac{\partial {\lambda }_1}{\partial s}-\frac{{\partial }^3{\kappa }_1}{\partial s^3},$$

$$0=(P-2){\kappa }_1+{\lambda }_1+4\frac{{\partial }^2{\kappa }_1}{\partial s^2}+\frac{{\partial }^2{\lambda }_1}{\partial s^2}.$$

$${\sigma }^2=\frac{n^2(n^2-1)}{n^2+1}\left[P-(n^2-1)\right].$$

Full text

I Introduction

Buckling is a ubiquitous instability in elastic materials Timoshenko and Gere (1961). The static buckling problem consists of determining the value of the load at which an instability takes place as well as the shape observed close to the instability. It can be recast as a standard bifurcation problem for which there exists an extensive literature Antman (2005); Euler (2008); Levien (2008). The classic example of this is the Euler buckling of a column, which can be readily observed by attempting to compress a thin piece of paper along its axis Levien (2008): different buckling modes (eigenfunctions) exist, each with a different buckling load (the eigenvalue). In most situations it is the eigenfunction with the smallest eigenvalue that is observed experimentally; indeed, the higher modes are usually dynamically unstable, rapidly transitioning to the lowest mode Pandey et al. (2014).

Many other elastic structures exhibit a similar phenomenology: For instance, the static buckling of an elastic ring under external pressure is a physical problem that has received much attention Levy (1884); Carrier (1947); Tadjbakhsh and Odeh (1967); Antman and Warner (1965); Flaherty and Keller (1973); Giomi and Mahadevan (2012); Biria and Fried (2014); Chen and Fried (2014); Biria and Fried (2015). These seminal works provide an in-depth analysis of the stability of the ring, showing that the lowest mode is the ‘figure-of-eight’ shape, determining higher modes and also offering different analytical and numerical methods through which to find the equilibrium shapes after bifurcation. These post-bifurcation equilibria are relevant to biological problems such as fluid flow through blood vessels Rubinow and Keller (1972).

In contrast to the many studies on the static stability of elastic rings, the problem of dynamic buckling has received much less attention. In dynamic buckling, the problem is to understand both the onset of instability under possibly dynamic loading and the evolution of the object’s shape after the onset of instability. These problems are considerably harder as they involve both space and time Wah (1970); Caflish and Maddocks (1984) and simple energy arguments cannot be used to obtain the post-buckling dynamics. For instance, twisted elastic rings going through a Michell instability Goriely (2006) fold into a ‘figure-of-eight’ shape without ever following equilibrium solutions. Dynamic loading can also be used to excite higher unstable modes that cannot be reached through static methods Goriely and Tabor (1997); Audoly and Neukirch (2005); Gladden et al. (2005). As such, dynamic buckling is a much richer phenomenon that offers many possibilities for new behaviors and, potentially, for the development of new active devices.

In this paper we study the dynamic evolution of an elastic ring subject to a sudden jump in external pressure. When an elastic ring is subject to a slowly increasing external pressure, it buckles into its first unstable mode, a flattened, ellipse-like shape referred to as mode-2. Higher modes with more lobes can be excited in dynamic loading and the problem here is to understand how these modes are selected and how they evolve after the onset of the instability. In particular, we shall show that the presence of inertia implies that higher modes are selected, and we identify the mode number as a function of the driving pressure. Experimentally, these modes can be obtained by following the dynamics of an elastic ring embedded in a soap film, as presented in a companion paper Box et al. (2020). Initially, the elastic ring is maintained in equilibrium by the balance of surface tension between an internal and external soap film. Once the external soap film is removed (by puncturing it at a point and allowing it to retract), the force from the internal soap film is suddenly unbalanced and pulls the ring inwards. This is equivalent to an external pressure applied instantaneously to the ring, which is illustrated schematically in Fig. 1. The experiments described in ref. Box et al. (2020) demonstrate that a circular elastic ring buckling under surface tension in this way can exhibit a range of modes with higher modes selected dynamically as the importance of surface tension increases. An illustrative experimental time series of these experiments is shown in Fig. 2.

To develop some intuition for the problem of mode selection in dynamic buckling, we begin by considering the one-dimensional toy problem of an elastic beam (rather than a circular ring) of thickness, $h$, under an in-plane compressive force ${\cal F}$ per unit width. In linear beam theory the governing equation for the vertical deflection $w(x,t)$ of a pinned-pinned beam of length $L$ is Lange and Newell (1971):

[TABLE]

where $\rho$ is the density and $B$ is the bending stiffness of the beam. This equation exhibits growing buckled solutions of the form $w(x,t)=K\mathrm{exp}({\sigma t})\cos[(2n+1)\pi x/L]$ if

[TABLE]

By examining the maximum value of the growth rate, $\sigma$, as a function of $n$ in (3), we see that the $n^{\mathrm{th}}$ mode is the fastest growing mode (and hence expected to be selected) if the compressive force is suddenly raised to ${\cal F}_{n}=2(2n+1)^{2}\pi^{2}B/L^{2}$. Alternatively, if a particular force ${\cal F}$ is imposed then we might expect that the observed mode would be

[TABLE]

since it will be close to the fastest growing mode. This selection of a mode other than the lowest Euler mode (which corresponds to $n=0$) is made possible by the presence of the inertia term.

Although the above calculation gives us an intuitive feel for the problem of dynamic buckling, it is expected to be only quantitatively accurate for an elastic beam subject to a controlled compressive force ${\cal F}$. Achieving such a condition experimentally is difficult and is perhaps closest to being achieved by an impactor with a known weight Gladden et al. (2005), though then the imposed force becomes dependent on the imposed boundary conditions. The case of an elastic ring embedded in a soap film comes closer to a truly controlled force, but the question then becomes how the calculation that leads to (4) is modified by the ring’s curvature $\kappa=1/a$, with $a$ the ring radius. In a related study Box et al. (2020), quantitative agreement with experiment was found by replacing the force ${\cal F}$ in (4) with the compressive force induced by the pressure difference $p$ combined with the ring curvature $1/a$ (Laplace’s law), which gives ${\cal F}=pa$. In this paper we go beyond this simple physical argument to account for the effects of curvature completely, both in terms of the onset of instability and also in the post-buckling of the ring.

The paper begins in §II with a formulation of the governing equations from first principles, together with our numerical technique and a presentation of the typical numerical results. To make analytical progress, §III first re-formulates the governing equations of §II in a form that is more convenient for analytical work. This allows us to present a detailed linear stability analysis that goes beyond the simple-minded prediction of (4). While the linear stability analysis is able to explain some features of our numerical results, other qualitative features of the numerics presented here and experiments of Box et al. (2020) are not explained by the solution of the linear problem. We therefore turn next to a weakly nonlinear analysis through which we derive an amplitude equation for the motion, which can be solved analytically. We illustrate the utility of this approach by applying the predictions of the weakly nonlinear analysis to explain various features of our numerical results, before summarizing our results in §IV.

II Theoretical formulation

II.1 Model problem

We consider a circular elastic ring of radius $a$ subject to an externally applied pressure $p$. We model the elastic ring as a rod of thickness $h$ and length $L=2\pi a$, whose shape is parametrized by its arc-length $s\in[0,L]$. Previous work has considered the vibrations of arches, accounting for the effects of extensibility Nayfeh et al. (1995); Neukirch et al. (2012). These works show that, under a controlled applied force, the effect of extensibility on the frequency of vibration is negligible provided that $\eta=h^{2}/(12L^{2})\ll 1$. The experiments of the companion to this paper Box et al. (2020) have $\eta\sim 10^{-4}$ and are performed with controlled applied force. We shall therefore model the ring as being inextensible and unshearable to simplify the analysis somewhat. The ring is composed of a material of density $\rho$, Young’s modulus $E$, and its circular cross section has cross sectional area $A_{\times}$ and second moment of area $I=A_{\times}^{2}/(4\pi)$; which gives a bending stiffness $B=EI=EA_{\times}^{2}/(4\pi)$.

The position of a material point on the ring centreline at a time $t$ is denoted by the vector $\bm{r}(s,t)=x\bm{e}_{x}+y\bm{e}_{y}+z\bm{e}_{z}\in\mathbb{R}^{3}$ (Fig. 1). Previous work has shown that the static equilibrium of a ring with an internal soap film may take a saddle shape Giomi and Mahadevan (2012); Chen and Fried (2014). However, in our problem the ring remains in the plane and we restrict our attention to planar deformations. With the assumption of planarity, we denote by $\bm{t}$, $\bm{n}$ and $\bm{k}$ the unit tangent, normal, and binormal vectors attached to the curve describing the centreline. They satisfy the property $\partial_{s}\bm{r}(s,t)=\bm{t}$ and $\partial_{s}\bm{t}=\kappa\bm{n}$, where $|\kappa|$ is the curvature.

The motion of the elastic ring is governed by Kirchhoff’s equations Love (1959); Audoly and Pomeau (2010). Let $p\bm{n}$ be the external body force per unit length due to the applied pressure (integrated over the ring thickness) and denote by $\bm{F}$, $\bm{M}$ the resultant force and moment, respectively, acting on the centreline. The balances of linear and angular momenta then lead to Antman (2005); Goriely (2017):

[TABLE]

where, in eqn (6), we have ignored the rotary inertia terms (Goriely, 2017, p. 115).

The governing equations are closed by the constitutive relation $\bm{M}=B\kappa(s,t)\bm{k}$ and the constraint that the ring is unshearable and inextensible. We note that, alternatively, an energy formulation can be used to establish the governing equations; for example, the static buckling of twisted elastic rings within a soap film is studied in Biria and Fried (2014, 2015).

We denote by $\theta$ the angle between the $x$-axis and the tangent $\bm{t}$, so that $\bm{t}=(\cos\theta,\sin\theta)$. We rescale all lengths by the radius of the initial unstressed ring $a$, pressure by $p_{*}=B/a^{3}$, and time by the inertial time scale $t_{*}=a^{2}\left({\rho A_{\times}}/{B}\right)^{1/2}$. In the new dimensionless variables, the governing equations projected on $\bm{e}_{x}$ and $\bm{e}_{y}$, are

[TABLE]

This non-dimensionalization introduces the dimensionless parameter

[TABLE]

which is the control parameter in the problem. It measures the importance of the work done by the external pressure compared to the bending energy of the ring, as described by Chen & Fried Chen and Fried (2014). We have the boundary conditions

[TABLE]

and initially the ring is stationary. The initial shape of the ring is close to that of a circle of unit radius, though a small amount of noise is added to initiate instability in our numerical simulations.

II.2 Numerical analysis

Before looking for asymptotic solutions, we present numerical solutions for the evolution of this elastic ring under imposed pressure. We solve (7–11) numerically using a finite difference scheme in which we discretize the system in space and time. This discretization produces a series of nonlinear equations, which are solved by using Newton’s method at each time step (see Appendix A for details). Note that we do not use a time-splitting projection method to enforce the inextensibility constraint; rather, our solution of the nonlinear equations obtained by discretization, ensures that the geometrical relationships (7) and (8) are satisfied, and hence that inextensibility is enforced at the same time. The initial velocity is zero, and the initial shape of the ring is circular with a uniformly distributed random number $\mathcal{R}_{i}\in[-\varepsilon,\varepsilon]$, with $\varepsilon\ll 1$, added to the local radius of curvature at each point; in particular, the discretized initial state at arc-length $s_{i}=(i-1)\Delta s$ is

[TABLE]

Note that this initial condition is only consistent with the geometrical constraint (7)–(8) to $O(\epsilon/(\Delta s))$. However, after the first time step, our numerical scheme ensures that the inextensibility constraint is satisfied subsequently (see Appendix A for details). Snapshots of the evolution of the ring profile for different values of the imposed pressure, $P$, are shown in Fig. 3. These show two trends: First, the selected mode number of the instability (defined as the number of lobes) increases with the imposed pressure $P$. Second, the instability develops notably faster for larger values of $P$. These initial observations will motivate our asymptotic analyses of the problem.

III Analysis

To facilitate the subsequent analysis, we first express the governing equations in the moving frame of the ring $(\bm{t},\bm{n},\bm{k})$, which is a right-handed orthonormal frame attached to the centreline of the ring at arc-length position $s$. We start by eliminating the force $\bm{F}$ from equation (5) to obtain a single equation for the vector position $\bm{r}$ together with the inextensibility constraint. The (dimensionless) bending moment is given by the constitutive equation $\bm{M}=\kappa(s,t)\bm{k}$. Here, and henceforth, we use a prime to denote derivatives with respect to the arc-length $s$ and an overdot for derivatives with respect to time $t$. Since $\kappa\bm{k}=\bm{r}^{\prime}\times\bm{r}^{\prime\prime}$, the resultant moment can be written

[TABLE]

Substituting (16) into the dimensionless version of (6) then gives

[TABLE]

which immediately implies that

[TABLE]

for some unknown functions $\lambda(s,t)$ and $\beta(s,t)$. However, since we restrict attention to planar deformations, we have $\beta=0$ and

[TABLE]

We note that, if redimensionalized, the term involving $\lambda$ in (19) would have the units of force and contributes to the tangential component of the internal force Singh and Hanna (2018). Using the expression for the force (19) in (5) gives

[TABLE]

where the last equation is the local inextensibility constraint. The two equations (20) form a system for the two unknowns $\bm{r}(s,t)$ and $\lambda(s,t)$.

Since we are interested in the evolution of the ring away from its initially circular configuration, we take advantage of the circular geometry to consider displacements from this shape with respect to the polar vectors associated with the circle. The curvature, $\kappa_{0}$, of the ring in its initial configuration, together with the initial tangential and normal vectors (denoted $\bm{t}_{0}$ and $\bm{n}_{0}$, respectively), are:

[TABLE]

Note that $\bm{t}^{\prime}_{0}=\bm{n}_{0}$, $\bm{n}^{\prime}_{0}=-\bm{t}_{0}$ and the undeformed ring is parametrized by $\bm{r}_{0}=-\bm{n}_{0}$.

We denote the displacements in the normal and tangential directions from the initial shape by $u(s,t)$ and $v(s,t)$, respectively. The ring shape may then be written as

[TABLE]

We introduce the angle $\phi(s,t)$ between the tangent of the deformed state and the tangent of the undeformed state by $\bm{t}=(\cos\phi)\,\bm{t}_{0}+(\sin\phi)\,\bm{n}_{0}$. Taking the derivative of (24) gives

[TABLE]

Hence, we have:

[TABLE]

To relate the curvature to the angle $\phi$ we use $\bm{t}^{\prime}=\kappa\bm{n}$, where $\bm{n}=-(\sin\phi)\bm{t}_{0}+(\cos\phi)\bm{n}_{0}$ so that

[TABLE]

We see that the curvature has two components: a contribution, $\kappa_{0}=1$, coming from the fact that the ring is initially curved, and a contribution from the deformation of the ring relative to its initial circular shape. Now taking the derivative of (20) with respect to arc-length gives

[TABLE]

In the absence of pressure, namely $P=0$, we recover the equation obtained by Burchard and Thomas (2003) for a closed loop elastica. (We note that in this particular case, local existence and uniqueness of the elastica solutions has been established for initial data in suitable Sobolev spaces, but that no global existence result exists Burchard and Thomas (2003). Similarly, Caflish and Maddocks (1984) proved global existence for particular initial data in a similar system that includes rotary inertial terms, but again with $P=0$.) We can resolve equation (28) into tangential and normal components giving a complete system of five equations for the five variables $(u,v,\lambda,\kappa,\phi)$ for the shape of an elastic ring subject to a normal pressure $P$.

[TABLE]

We now make use of the formulation presented in eqns (29–33) to perform a linear stability analysis of the problem; this is followed by a multiple scales analysis that allows us to examine the development of the instability beyond infinitesimal deformations (linear stability).

III.1 Linear stability analysis

We begin by noting that (29–33) admit as a solution the undeformed circle: $v_{0}=u_{0}=\phi_{0}=0$, $\kappa_{0}=1$, provided that $\lambda_{0}=-P-1$. Physically, this corresponds to the ring having an internal compressive force, $\lambda_{0}=-P-1$, which, combined with its initial curvature, balances the external pressure (via Laplace’s law). Our interest lies in determining the stability of this equilibrium state. Letting $u=u_{0}+\varepsilon u_{1}$, $v=v_{0}+\varepsilon v_{1}$, $\phi=\phi_{0}+\varepsilon\phi_{1}$, $\kappa=\kappa_{0}+\varepsilon\kappa_{1}$, $\lambda=\lambda_{0}+\varepsilon\lambda_{1}$, with $\varepsilon\ll 1$ arbitrary, we expand (29–33) to first order in $\varepsilon$ to obtain:

[TABLE]

Seeking solutions of the form $u_{1}=a_{1}\exp{(\sigma t+\text{i}n\theta)}$, $v_{1}=b_{1}\exp{(\sigma t+\text{i}n\theta)}$, $\kappa_{1}=c_{1}\exp{(\sigma t+\text{i}n\theta)}$ , $\phi_{1}=d_{1}\exp{(\sigma t+\text{i}n\theta)}$, $\lambda_{1}=e_{1}\exp{(\sigma t+\text{i}n\theta)}$, we obtain a homogeneous linear system that has a non-trivial solution if, and only if, the associated determinant vanishes. This happens when

[TABLE]

This dispersion relation may be recovered from that obtained previously by Wah Wah (1970) if one neglects rotational inertia in their formulation (corresponding to the limit in which their parameter $\mu\propto a^{2}/A_{\times}\to\infty$). We also note that if $P=0$ then $\sigma^{2}$ is negative so that $\sigma$ is imaginary and all the modes oscillate with a frequency given by the imaginary part of $\sigma$, reproducing a result first obtained by Hoppe in 1871 Hoppe (1871).

Our main interest lies in the instability that may be observed with positive external pressure, $P>0$. In this regard, note that $\sigma=0$ when $P=p_{n}=n^{2}-1$, giving another perspective on the $n^{\mathrm{th}}$ static buckling load of a ring under normal pressure Tadjbakhsh and Odeh (1967); Flaherty et al. (1972). For a given $P>0$, all modes with $n\leq(P+1)^{1/2}$ are unstable. However, the mode that is expected to be most relevant in the development of instability from arbitrary initial conditions is the most unstable mode (the mode with the largest value of $\sigma$) for a given $P$. Since the mode number $n$ is an integer, the calculation of the most unstable mode number is obtained for a pressure interval as shown in Fig. 4 (see red dots). This value agrees well with the mode number observed in numerical simulations. Note that the numerically determined mode number is determined from the simulations at the point at which the instability is first visible in the shape. For large $n$ this appears to introduce an error in $n$ of 1, perhaps because one lobe disappears before instability has become numerically observable. An alternative, if approximate, approach is to treat $n$ as a continuous variable and obtain the value of $n$ that maximizes $\sigma$ in (39) using standard methods. This calculation reveals that the most unstable mode, $n_{c}$, is a solution of

[TABLE]

For $P\gg 1$, we find $n_{c}\sim(P/2)^{1/2}$, which in practice gives a very good account of the numerical results, and the more detailed analytical calculation, as shown by the dashed line in Fig. 4. Note that this result is precisely that obtained by letting $L=2\pi a$ and ${\cal F}=pa$ in (4): for large mode numbers (large dimensionless pressures) the most important role of the ring’s curvature is to convert the externally applied pressure $P$ into an in-plane compression $\lambda$ within the ring (corresponding to ${\cal F}$ in the introductory beam analysis).

III.2 Weakly nonlinear analysis

The predictions of the linear stability analysis hold only for early times, $t\ll 1$. As the amplitude of the perturbation to the circular state grows, various nonlinear terms compete with the exponentially growing terms found by solving the linear system; eventually, these nonlinearities can no longer be ignored. As a simple illustration of the importance of nonlinear effects, note that the linear theory predicts that the area enclosed by the ring is constant to leading order: $A(t)\approx\tfrac{1}{2}\int_{0}^{2\pi}(1+u)^{2}\leavevmode\nobreak\ \mathrm{d}s=\pi+\mathcal{O}(\varepsilon^{2})$. Moreover, if one naively calculates the correction at $\mathcal{O}(\varepsilon^{2})$ using the linear expansion, one finds that the area enclosed should increase with time, in clear contradiction to the numerical results of Fig. 3 and experimental results Box et al. (2020). This prediction is a result of the inconsistency of using a result determined at $\mathcal{O}(\varepsilon)$ to make predictions about a correction at $\mathcal{O}(\varepsilon^{2})$.

To go beyond the linear theory we perform a weakly nonlinear analysis. From (39), we see that as the external pressure increases, the ring remains stable as long as $\sigma^{2}\leq 0$ for all $n\geq 0$. Therefore the ring is stable only when $P\leq p_{n}$ for all $n\geq 2$ (with $p_{2}=3$ the critical pressure above which the circular solution of the ring becomes unstable Tadjbakhsh and Odeh (1967)). We are interested in the dynamic evolution of a critical mode $n$ when the imposed pressure slightly exceeds the critical pressure $p_{n}$. To find this evolution, we introduce a new small parameter, $\epsilon=\sqrt{P-p_{n}}\ll 1$ (distinct from the arbitrary $\varepsilon$ used previously), to measure how far the pressure is above the critical pressure for a given mode. Therefore, we introduce $P=n^{2}-1+\epsilon^{2}$ in the equations that now depend explicitly on $\epsilon$. Next, we expand all variables to third order in $\epsilon$. For example, we write $u=u_{0}+\epsilon u_{1}+\epsilon^{2}u_{2}+\epsilon^{3}u_{3}+\mathcal{O}(\epsilon^{4})$ with analogous expansions for the other variables $v,\kappa,\phi,\lambda$. These expansions are substituted into (29-33) and the resulting hierarchy of linear systems solved. Further details are given in Appendix B.

At $\mathcal{O}(\epsilon)$, the solution is given by:

[TABLE]

where we have an arbitrary amplitude $\alpha_{1}$, that evolves on the long time scale $t=O(1/\epsilon)$.

At $\mathcal{O}(\epsilon^{2})$, the solution is given by:

[TABLE]

Note that in the above expressions the $\mathcal{O}(\epsilon^{2})$ radial displacement, $u_{2}$, is uniform, independent of arc-length $s$. This is contrary to the higher order radial displacement, $u_{1}$, which is oscillatory and, hence, does not contribute to a net radial displacement (once integrated over the whole ring). From this finding, we conclude that the clear visual impression that the ring shrinks radially (as seen in both the experiments Box et al. (2020) and numerical simulations of Figs 2 and 3, respectively) is a second-order effect — this striking feature of the problem can only be properly understood by going beyond the standard linear stability analysis.

The value of the amplitude $\alpha_{1}$ as a function of time is determined by a multiple-time expansion to derive the so-called amplitude equation. This equation is obtained as a compatibility condition for the existence of bounded solutions by using the Fredholm alternative Lange and Newell (1971); Goriely et al. (2001) for the linear system to third order (see Appendix B for details). This condition reads

[TABLE]

with

[TABLE]

If we neglect the non-linear term, $\mu\alpha^{3}$, we recover the results of the linear stability theory, as expected, since (48) gives the same linear growth rate $\sigma$ as (39).

To gain analytical insight into the behavior of the solution of the amplitude equation (47) we consider the particular problem of introducing a finite perturbation of size $\alpha(0)$ from rest, so that $\dot{\alpha}(0)=0$. The solution of the linearized version of (47) would then be simply $\alpha(t)=\alpha(0)\cosh(\sigma t)$. This linear solution, together with the form of (47), suggests that we first introduce the rescaled variables $\tilde{t}=\sigma t$, $\tilde{\alpha}=\alpha/\alpha(0)$, which transform (47) to

[TABLE]

where

[TABLE]

is the sole remaining dimensionless parameter and measures the strength of the nonlinearity in the amplitude equation.

Eqn (49) may readily be solved subject to the initial conditions $\tilde{\alpha}=1$ and $\mathrm{d}\tilde{\alpha}/\mathrm{d}\tilde{t}=0$ at $\tilde{t}=0$ to give

[TABLE]

where $I_{E}(\phi,m)=\int_{0}^{\phi}(1-m\sin^{2}\theta)^{-1/2}\mathrm{\leavevmode\nobreak\ d}\theta$ is the Elliptic Integral of the first kind Olver (2010). The solution in (51) holds strictly only while $\tilde{\alpha}<[(1-\tilde{\mu})/\mu]^{1/2}$, which in turn requires $\tilde{t}<\tilde{T}/2$ where the period of the motion

[TABLE]

This analytical solution allows us to discuss briefly the effect of the nonlinearity. We can see in Fig. 5 the solution for three values of the parameter $\tilde{\mu}$. As might be expected from the rescaled amplitude equation (49), the results are closer to the linear solution, $\tilde{\alpha}(\tilde{t})=\cosh\tilde{t}$, when $\tilde{\mu}$ is small. More quantitatively, the nonlinear term may be neglected while $\tilde{\alpha}\ll\tilde{\mu}^{-1/2}$, i.e. while $\tilde{t}\ll\cosh^{-1}(\tilde{\mu}^{-1/2})\sim-(\log\tilde{\mu})/2$. In particular, since for large $n$, $\sigma\sim n(P-p_{n})^{1/2}$ and $\mu\sim 3n^{6}/8$ so that $\tilde{\mu}\sim 3n^{4}[\alpha(0)]^{2}/[16(P-p_{n})]$, we expect that the role of nonlinearities will become more important more quickly for larger imposed pressures. (Note that this sense of ‘faster’ onset of nonlinearity is measured in the natural time variable $\tilde{t}=\sigma t$, and so goes beyond the prediction of linear stability analysis that the linear growth rate of instability also increases with pressure.)

III.3 Applications of the weakly nonlinear analysis

Since the amplitude $\alpha$ may readily be computed as a function of time (either using the above expressions or by direct numerical solution of the amplitude equation), we are now in a position to compare directly the numerical solutions of the full system with the post-bifurcation solution to second-order that is provided by our weakly nonlinear analysis. We study several aspects of the problem that are accessible in the full numerical solutions to illustrate the application of our weakly nonlinear analysis. Before doing so, we note that the shape of the ring is related to the radial and tangential displacements ($u$ and $v$, respectively) via (24), which may be written in component form as

[TABLE]

As a result, the displacement field can easily be inferred from the shape, and vice versa.

To facilitate the comparison between the prediction of the weakly nonlinear analysis and our numerical solution of the problem, we ensure that we use the same initial conditions in both problems; in particular, we solve (7)–(11) numerically with an initial condition seeded by the weakly nonlinear analysis with a particular mode number, i.e. we use (53)–(54) with $u$, $v$ being the perturbation solution provided by (41)–(46) at $t=0$,

[TABLE]

rather than the random initial condition (15) used for computations with fixed pressure, but $n$ not known a priori. Here, we choose $\alpha(0)=0.01$.

Shape

We compare the shape obtained by numerically solving the governing equations with that predicted by the weakly nonlinear analysis. The predicted shape is reconstructed from the amplitude $\alpha(t)$ and using the predictions (41) and (44). In detail, the profile predicted by the weakly nonlinear analysis may be constructed as follows: using (41) and (44), the displacement field $(u,v)$ is given by $u=u_{0}+\epsilon u_{1}+\epsilon^{2}u_{2}$ and $v=v_{0}+\epsilon v_{1}+\epsilon^{2}v_{2}$, which may then be directly substituted into (53)–(54). The results of this reconstruction are shown in Fig. 6 together with the numerical results for $P=3.5$ and $P=8.7$, which correspond to the $n=2$ and $n=3$ modes, respectively. The results of this comparison are very favourable.

Amplitude

We compare the amplitude $\alpha(t)$ with that obtained from the numerical simulations. To do so, we compute numerically the function $u(s,t)$ and extract the amplitude of its first Fourier component:

[TABLE]

We see from Fig. 7 that the amplitude equation captures quantitatively the evolution of the $n=2$ mode but only qualitatively captures the evolution of the $n=3$ mode.

Area

From the post-bifurcation solution to second order, we can compute the area of the ring including a perturbation with mode number $n$. We find that this area is

[TABLE]

The change in area, $\Delta A(t)=A(0)-A_{n}(t)=\pi-A_{n}(t)$ and the relative change in area is

[TABLE]

The plots in Fig. 8a show that the prediction of (59) is in good agreement with the detailed numerical simulations of the full problem for early times, while the plots in Fig. 8b show that at early times the relative change in area grows according to $\Delta A/A(0)\propto\exp(2\sigma t)$ with $\sigma$ the growth rate of the linear instability. Finally, Fig. 8c shows that the error in the area change occurs at higher order in $\alpha$ than $\mathcal{O}(\alpha^{2})$, confirming that the weakly nonlinear analysis presented here is correct to this order.

Compressive force $\lambda$

The amplitude equation can also be used to study the evolution of the compressive force $\lambda$ within the ring: the predicted behavior of $\lambda(s,t)$ can readily be computed once the amplitude equation for $\alpha(t)$ (47) is solved numerically.

The comparison between the prediction of the weakly nonlinear analysis and full numerical simulations is shown in Fig. 9. As with other variables, the comparison between numerics and the weakly nonlinear analysis is very good, particularly at early times. However, this plot of $\lambda(s,t)$ reveals two features that are not so readily observed in other variables: the oscillations in $\lambda(s,t)$ for a fixed time $t$ are not up-down symmetric and, in particular, the crest of these oscillations splits in two, showing the importance of a higher frequency oscillation in the arc-length $s$. Both of these features can be understood by observing from (44) that the prefactor of $\lambda_{2}\sim n^{4}\alpha_{1}^{2}$ while (41) shows that the prefactor of $\lambda_{1}\sim n^{2}\alpha_{1}$. As a result, $\lambda_{2}$ is able to compete with $\lambda_{1}$ despite being strictly at higher order in $\alpha$: the $\cos 2ns$ term in the expansion of $\lambda(s,t)$ can compete with the $\sin ns$ term, breaking up-down symmetry and leading to the splitting of the crests. This is not so apparent in other variables, since the higher order terms remain sub-dominant for longer; for example, $\kappa_{2}\sim\alpha_{1}^{2}n^{2}$ while $\kappa_{1}\sim\alpha_{1}n^{2}$.

The appearance of second-order variables with frequency $2n$ might suggest this as a precursor of a secondary bifurcation with frequency doubling. Strictly speaking, this is not a secondary instability as these higher modes are present at the primary bifurcation but have very low amplitudes. Further, we note a competing effect that tends to prevent frequency doubling. Indeed, in the last panel of Fig. 3 (with $P=256$), we observe several events in which two crests merge into a single crest which implies that lobes progressively disappear with time. We hypothesize that this merger occurs via the ‘snap-through’ of the internal element, as observed in a simple elastic arch Pandey et al. (2014).

IV Conclusion

We have presented both linear and weakly nonlinear analyses of the dynamic buckling of an elastic ring subject to a suddenly imposed external pressure, $p$. Returning to our initial discussion of the dynamic buckling of a beam under a compressive force ${\cal F}$, we speculated that the same argument might hold for a ring under pressure $p={\cal F}/a$ and, further, that this is the most important effect of the ring’s curvature. This analogy is enough to reproduce the main results of the linear stability analysis for sufficiently large pressures; detailed calculations showed that the most unstable mode is

[TABLE]

which grows with dimensional growth rate

[TABLE]

These large pressure results are precisely the same as those for a beam with ${\cal F}=pa$. We note that while inertia is required to observe buckling at a high mode number, (60) shows that the mode number selected is independent of the inertia of the ring, $\rho A_{\times}$. Rather the inertia of the ring sets the time scale for the growth of instability, as shown in (61). The apparent contradiction that inertia is required to observe a high mode number but is not involved in the selection of a mode number is resolved Box et al. (2020) by varying the loading rate to be below that in (61).

While a linear stability analysis is sufficient to understand the early stages of the dynamics, it does not provide any insight into how the instability develops further. In particular, the linear stability analysis shows that the area enclosed by the loop does not change to leading order; moreover, if one attempts to use its results at higher order, one comes to the obviously erroneous conclusion that the enclosed area increases as instability proceeds. Therefore, we performed a weakly nonlinear analysis that allows us to show that the relative change in area $\Delta A/A(0)$ grows exponentially in time, with growth rate $2\sigma_{\mathrm{max}}$. The weakly nonlinear analysis is able to explain other features of our numerical simulations, including an apparent frequency doubling of the compressive force $\lambda$ within the ring. While this might be expected to lead to an increase in the mode number with time, we note that our numerical simulations show that crests tend to merge with time, presumably due to the increase in confinement that they experience. The time scale of this coarsening remains to be understood.

The dynamic buckling we have studied here may have analogues in other systems of interest. For example, similar dynamics may be relevant in the collapse of veins Chow and Mak (2006) for which it is known that inertial effects due to the fluid–structure interaction are important Jensen and Heil (2003); Grotberg and Jensen (2004); however, we also note that in these systems axial modes of deformations often dominate the radial modes considered here. Finally, we note that the ring geometry we have studied may be a more promising paradigm within which to understand dynamic buckling instabilities in tissue and soft materials. While the uniaxial compression of a rod has often been studied Gladden et al. (2005), it can be difficult to impose a known force experimentally and the resulting instability is sensitive to the imposed boundary conditions. The circular symmetry of the ring problem renders the latter issue irrelevant while our analysis has shown that, up to the instability, the ring curvature translates the normal force given by the applied pressure into a uniform in-plane compression.

Appendix A Details of the numerical scheme

In this Appendix we describe the numerical procedure used to determine the evolution of an elastic ring subject to a normal pressure; this motion is governed by the system (7)–(15) and so it is a numerical solution of these equations that we seek. We begin by discussing the discretization of the equations used numerically, before turning to an examination of the constraint associated with the ring’s inextensibility.

A.1 Discretization

Following Santillan (2007) we discretize (7)–(15) in time $t$ and in arc-length $s$; this results in a system of nonlinear equations that are solved using Newton’s method. In particular, we consider discrete time points $t_{j}=(j-1)\Delta t$, where typically $\Delta t=5\times 10^{-3}$ in the simulations reported in the main paper. The arc-length along the ring is discretized on $N+1$ grid-points at $s_{i}=(i-1)\Delta s$ where $\Delta s=2\pi/N$, and $i\in[1,N+1]$.

The system (7)–(15) does not involve any time derivative of the forces $F_{x}$ and $F_{y}$; we therefore cannot use a simple forward time stepping scheme. Instead, we need to solve for these forces and the rest of the equations together, in a single time step. To do this we use a one-sided time derivative Chung (2002) accurate to $\mathcal{O}(\Delta t^{2})$, namely:

[TABLE]

Derivatives with respect to arc-length are discretized using central differences. The discretized versions of (7)–(11) read:

[TABLE]

Note that here we use the notational convention that $f_{i+1/2}=(f_{i}+f_{i+1})/2$ for $f\in\{v_{x},v_{y},\theta,\kappa,F_{x},F_{y}\}$.

Given the values of the variables at times $t_{j-2}$ and $t_{j-1}$, the system (64) – (71) can be solved at time $t_{j}$ using a Newton solver; we use the MATLAB routine fsolve to find these solutions.

A.2 Inextensibility constraint

Here we discuss the constraint of inextensibility focussing, in particular, on whether the total length of the ring is preserved over time during the integration. Our approach is intended to ensure that the inextensibility constraints (7)–(8) are automatically satisfied; this is achieved by solving (64)–(65) together with the remaining equations (66)-(71) at each time step.

In the following, we check that the inextensibility constraint is satisfied for the simulations reported in Fig. 3. To do this, we numerically determine the arc-length of the ring as

[TABLE]

In Fig. 10, we plot the magnitude of the relative error in the total arc-length of the ring, $|(\ell(t_{j})-2\pi)/2\pi|$, as a function of time for different values of the driving pressure and for different spatio-temporal resolutions. Fig. 10a shows simulations for different values of pressures, $P=8,32,128$, at a fixed value of the time step and number of grid points ($\Delta t=10^{-3}$ and $N=256$). In contrast, Fig. 10b shows the relative error for simulations with a fixed pressure, $P=128$, but different spatial and temporal resolutions. We see here that the ring does indeed maintain its length over time as guaranteed by our numerical scheme — the magnitude of the relative error in the total length for the simulations reported in Fig. 3 remains less than $10^{-12}$ throughout.

Finally, we discuss two limitations of our numerical scheme. First, the initial conditions we used, i.e. (15), contains a random noise at each spatial grid point. This means that the initial condition is mesh dependent. However during the simulations reported in Fig. 3 we noticed that, at the first time step, the solutions initially converge to a circular ring shape, and then the fastest growing modes start emerging. This suggests that the observed instability is an intrinsic feature of the solutions of (7) -(11). Second, it is to be noted that since we have used a BDF2 (Backward Differentiation Formula) in (62) and (63), it follows that our scheme will be dissipative, and only accurate at $\mathcal{O}(\Delta t^{2})$. This limitation implies that the solutions reported in this paper cannot be expected to be accurate at very late time — a limitation that might be overcome by using a symplectic integrator or, in general, a geometric integrator Marsden et al. (1998).

Appendix B Details of the weakly nonlinear analysis

In this Appendix, we provide more details of the weakly nonlinear analysis that is performed to describe the behavior when the pressure is just above the critical pressure of buckling of the $n^{\text{th}}$ mode. We introduce $\bm{U}=[u,v,\phi,\kappa,\lambda]^{T}$ and linearize equations (29)–(33) around the static solution $\bm{U}_{0}=[0,0,0,1,-P-1]^{T}$ by computing $\bm{U}=\bm{U}_{0}+\epsilon\bm{U}_{1}$ to first order. We obtain a linear homogeneous system of the form

[TABLE]

where $\mathcal{L}_{P}$ is a linear differential operator in $s$ and $t$ with constant coefficients and is dependent on the control parameter $P$. The solution of this equation is given by (41) with $\alpha_{1}(t)=\exp(\sigma t)$. Here, the growth rate $\sigma$ is given by the dispersion relation depending on the mode $n$, and pressure $P$: $\sigma=\sigma(n,P)$. We define the critical pressure $p_{n}$ by the first positive root of $\sigma(n,p_{n})=0$, which leads to $p_{n}=n^{2}-1$. The leading-order problem therefore precisely reproduces the results of the linear stability analysis, as it must.

Next, we perform the weakly nonlinear analysis. We assume that the pressure $P=p_{n}+\epsilon^{2}$ is slightly above the critical pressure $p_{n}$. Doing so, we relate the distance to the bifurcation to a small parameter $\epsilon$. Furthermore, since we place ourselves at the bifurcation, the only dependence on time is through a slow time scale. An analysis of the dispersion relation shows that the relevant growth rate is $O(\epsilon)$ and hence we expect evolution to occur on a time scale $\tau=\epsilon t$. This motivates the expansion

[TABLE]

where $\bm{U}_{0}$ is as before and $\bm{U}_{j}=[u_{j}(s,\epsilon t),v_{j}(s,\epsilon t),\phi_{j}(s,\epsilon t),\kappa_{j}(s,\epsilon t),\lambda_{j}(s,\epsilon t)]^{T}$. Substituting this expansion into the governing equations (29)–(33) results in a hierarchy of systems for the $\bm{U}_{j}$. To first order, the system reads $\mathcal{L}_{p_{n}}\bm{U}_{1}=0$ and due to the particular choice of $P$, we obtain (41) where $P=p_{n}$ and $\alpha_{1}(\tau)$ is now arbitrary.

To second order, we obtain a system of the form

[TABLE]

whose unknown is $\bm{U}_{2}$ and where $\mathbf{G}_{2}(\bm{U}_{1})$ is a function of $\bm{U}_{1}$. The solution of this system is

[TABLE]

To third order, the system for $\bm{U}_{3}$ reads

[TABLE]

where the extra dependence on the second derivative of $\alpha_{1}(\tau)$, $\alpha_{1,\tau\tau}$, comes from the assumption of the long time scale. In general, this system does not have a solution as the inhomogenous term is not in the image of the operator $\mathcal{L}_{p_{n}}$. Therefore, a compatibility condition, the so-called Fredholm alternative Goriely et al. (2001); Neukirch et al. (2012), needs to be used. Let $\bm{\xi}$ be the solution of the homogeneous adjoint problem $\mathcal{L}_{p_{n}}^{\ast}\bm{\xi}=\mathbf{0}$. Then, this compatibility condition to third order simply reads

[TABLE]

This condition provides a differential equation for $\alpha_{1}$. Taking into account the different change of variable and defining $\alpha(t)=\epsilon\alpha_{1}(\epsilon t)$, we obtain the amplitude equation (47).

Once $\bm{U}$ is known to second order, we may use the relations (24) to express the position of the ring in Cartesian coordinates. Accurate to $\mathcal{O}(\epsilon^{2})$ we find that:

[TABLE]

We may use (27) to express $\theta(s,t)$ which is given by:

[TABLE]

Similarly, the expression for the internal forces may be expressed in Cartesian coordinates using (19) which is given component-wise as:

[TABLE]

Acknowledgements.

**Acknowledgments– ** The research leading to these results has received funding from the European Research Council under the European Union’s Horizon 2020 Programme/ERC Grant Agreement no. 637334 (DV) and the Engineering and Physical Sciences Research Council grant EP/R020205/1 (AG). We are grateful to Finn Box for comments on an earlier version of this article and to Peter Howell for discussions about dynamic buckling that led to this work.