Snapping of elastic strips with controlled ends

TL;DR

This paper analyzes the snapping behavior of elastic strips with controlled boundary conditions, introducing a universal snap surface that predicts critical conditions for snapping, validated by experiments and simulations, with applications in energy harvesting and robotics.

Contribution

It develops an analytical framework using Euler's elastica to identify stable configurations and defines a universal snap surface for elastic strips with controlled ends.

Findings

Universal snap surface accurately predicts snapping conditions.

Experimental and simulation results confirm theoretical predictions.

Insights into elastic energy release useful for impulsive motion optimization.

Abstract

Snapping mechanisms are investigated for an elastic strip with ends imposed to move and rotate in time. Attacking the problem analytically via Euler's elastica and the second variation of the total potential energy, the number of stable equilibrium configurations is disclosed by varying the kinematics of the strip ends. This result leads to the definition of a `universal snap surface', collecting the sets of critical boundary conditions for which the system snaps. The elastic energy release at snapping is also investigated, providing useful insights for the optimization of impulsive motion. The theoretical predictions are finally validated through comparisons with experimental results and finite element simulations, both fully confirming the reliability of the introduced universal surface. The presented analysis may find applications in a wide range of technological fields, as for instance energy harvesting and jumping robots.

Equations

$$x(s)={\int }_0^s\cos \theta (\varsigma )\text{d}\varsigma ,y(s)={\int }_0^s\sin \theta (\varsigma )\text{d}\varsigma ,$$

$${\int }_0^l\sin \theta (s)\text{d}s=0.$$

$$x(l)=d,\theta (0)={\theta }_0,\theta (l)={\theta }_l,$$

$$\begin{split}\mathcal{V}\big{(}\theta(s)\big{)}=&\,\mathcal{E}+R_{x}\left[-d+\int_{0}^{l}\cos{\theta(s)}\,\textup{d}s\right]+R_{y}\left[\int_{0}^{l}\sin{\theta(s)}\,\textup{d}s\right]+M_{l}\bigg{[}\theta_{l}-\theta(l)\bigg{]}-M_{0}\bigg{[}\theta_{0}-\theta(0)\bigg{]},\end{split}$$

$$\mathcal{E}=\frac{B}{2}{\int }_0^l{{\theta }^{\prime }(s)}^2\text{d}s,$$

$${\psi }^{\prime \prime }(s)+{\gamma }^2\sin \psi (s)=0\forall s\in [0,l].$$

$${\psi }^{\prime }(s)=\pm \gamma \sqrt{2}\sqrt{\cos \psi (s)+{\rm Y}},$$

$${\rm Y}=\frac{2-{\xi }^2}{{\xi }^2},\xi \in \left[0,\sqrt{\frac{2}{1-\underset{s\in (0,l)}{\min }\left\{\cos \psi (s)\right\}}}\right],$$

$$\begin{array}{r}{\chi }_0=\frac{{\theta }_0-\beta }{2},{\chi }_l=\frac{{\theta }_l-\beta }{2}.\end{array}$$

$${\rm Y}=-\cos {\hat{\psi }}_1.$$

$$\theta ({\hat{s}}_j)=-(-1{)}^j\theta ({\hat{s}}_1)+\left[1+(-1{)}^j\right]\beta .$$

$$\psi^{\prime}(s)=\left\{\begin{array}[]{lll}&(-1)^{p}\gamma\sqrt{2}\sqrt{\cos{\psi(s)}-\cos{\widehat{\psi}_{1}}},&\forall s\in[0,\,\widehat{s}_{1}],\\[8.53581pt] &(-1)^{p}(-1)^{j}\gamma\sqrt{2}\sqrt{\cos{\psi(s)}-\cos{\widehat{\psi}_{1}}},&\forall s\in[\widehat{s}_{j},\,\widehat{s}_{j+1}]\qquad\forall\,\,j\in[1,m-1],\\[8.53581pt] &(-1)^{p}(-1)^{m}\gamma\sqrt{2}\sqrt{\cos{\psi(s)}-\cos{\widehat{\psi}_{1}}},&\forall s\in[\widehat{s}_{m},\,l],\end{array}\right.$$

$$\eta =\sin \frac{{\hat{\psi }}_1}{2},\eta \sin \omega (s)=\sin \frac{\psi (s)}{2}.$$

$${\omega }_0=\arcsin \left(\frac{1}{\eta }\sin \frac{{\theta }_0-\beta }{2}\right),{\omega }_l=(-1{)}^m\arcsin \left(\frac{1}{\eta }\sin \frac{{\theta }_l-\beta }{2}\right)+(-1{)}^{p}m\pi .$$

$$\gamma l=\left\{\begin{array}[]{lll}\xi\left|K(\chi_{l},\,\xi)-K(\chi_{0},\,\xi)\right|,&m=0,\\[8.53581pt] \left|K(\omega_{l},\,\eta)-K(\omega_{0},\,\eta)\right|,&m\neq 0,\end{array}\right.$$

$$\theta(s)=\left\{\begin{array}[]{lll}\beta+2\,\textup{am}\left(\dfrac{s}{l}\left(K(\chi_{l},\,\xi)-K(\chi_{0},\,\xi)\right)+K(\chi_{0},\,\xi),\,\xi\right),&m=0,\\[8.53581pt] \beta+2\,\arcsin{\left[\eta\,\textup{sn}\left(\dfrac{s}{l}\left(K(\omega_{l},\,\eta)-K(\omega_{0},\,\eta)\right)+K(\omega_{0},\,\eta),\,\eta\right)\right]},&m\neq 0.\end{array}\right.$$

$$K\left(\sigma,\,\varphi\right)=\int_{0}^{\sigma}\frac{\textup{d}\phi}{\sqrt{1-\varphi^{2}\sin^{2}{\phi}}},\qquad\sigma=\textup{am}\bigg{(}K\left(\sigma,\,\varphi\right),\,\varphi\bigg{)},\qquad\textup{sn}(u,\,\varphi)=\sin{\left(\textup{am}(u,\,\varphi)\right)}.$$

$$\text{sn}\left(\sigma \phi ,\frac{1}{\phi }\right)=\phi \text{sn}\left(\sigma ,\phi \right),$$

$$\xi =1/\eta ,\sin \omega (s)=\xi \sin \chi (s),K\left(\omega (s),\eta \right)=\xi K\left(\chi (s),\xi \right).$$

$$\mathcal{E}=\left\{\begin{array}[]{lll}\displaystyle\frac{2B}{l}\big{[}K(\chi_{l},\,\xi)-K(\chi_{0},\,\xi)\big{]}\big{[}E(\chi_{l},\,\xi)-E(\chi_{0},\,\xi)\big{]},&m=0,\\[8.53581pt] \displaystyle\frac{2B}{l}\big{[}K(\omega_{l},\,\eta)-K(\omega_{0},\,\eta)\big{]}\bigg{\{}E(\omega_{l},\,\eta)-E(\omega_{0},\,\eta)-(1-\eta^{2})[K(\omega_{l},\,\eta)-K(\omega_{0},\,\eta)\big{]}\bigg{\}},&m\neq 0,\end{array}\right.$$

$$x(s)=\big{[}\cos{\beta}\,\mathcal{A}_{m}(s)+\sin{\beta}\,\mathcal{B}_{m}(s)\big{]}l,\qquad y(s)=\big{[}\sin{\beta}\,\mathcal{A}_{m}(s)-\cos{\beta}\,\mathcal{B}_{m}(s)\big{]}l$$

$$\displaystyle\mathcal{A}_{m}(s)=\left\{\begin{array}[]{lll}\displaystyle\frac{2}{\xi^{2}}\frac{\mathscr{E}\left(\dfrac{s}{l}\left(K(\chi_{l},\,\xi)-K(\chi_{0},\,\xi)\right)+K(\chi_{0},\,\xi),\,\xi\right)+\mathscr{E}\left(K(\chi_{0},\,\xi),\,\xi\right)}{K(\chi_{l},\,\xi)-K(\chi_{0},\,\xi)}-\frac{2-\xi^{2}}{\xi^{2}}\frac{s}{l},&m=0,\\[14.22636pt] \displaystyle 2\frac{\mathscr{E}\left(\dfrac{s}{l}\left(K(\omega_{l},\,\eta)-K(\omega_{0},\,\eta)\right)+K(\omega_{0},\,\eta),\,\eta\right)-\mathscr{E}\left(K(\omega_{0},\,\eta),\,\eta\right)}{K(\omega_{l},\,\eta)-K(\omega_{0},\,\eta)}-\frac{s}{l},&m\neq 0,\end{array}\right.$$

$$\displaystyle\mathcal{B}_{m}(s)=\left\{\begin{array}[]{lll}\displaystyle\frac{2}{\xi^{2}}\frac{\textup{dn}\left(\dfrac{s}{l}\left(K(\chi_{l},\,\xi)-K(\chi_{0},\,\xi)\right)+K(\chi_{0},\,\xi),\,\xi\right)-\textup{dn}\left(K(\chi_{0},\,\xi),\,\xi\right)}{K(\chi_{l},\,\xi)-K(\chi_{0},\,\xi)},&m=0,\\[14.22636pt] \displaystyle 2\,\eta\,\frac{\textup{cn}\left(\dfrac{s}{l}\left(K(\omega_{l},\,\eta)-K(\omega_{0},\,\eta)\right)+K(\omega_{0},\,\eta),\,\eta\right)-\textup{cn}\left(K(\omega_{0},\,\eta),\,\eta\right)}{K(\omega_{l},\,\eta)-K(\omega_{0},\,\eta)},&m\neq 0.\end{array}\right.$$

$$\text{cn}(u,\phi )=\cos \left(\text{am}(u,\phi )\right),\mathcal{E}\left(\sigma ,\phi \right)=E(\text{am}(\sigma ,\phi ),\phi ),\text{dn}(u,\phi )=\sqrt{1-{\phi }^2{\text{sn}}^2(u,\phi )},$$

$$E\left(\sigma ,\phi \right)={\int }_0^{\sigma }\sqrt{1-{\phi }^2{\sin }^2\varphi }\text{d}\varphi .$$

$$\{\begin{array}{rl}& d\cos \beta ={\mathcal{A}}_m(l)l,\\ & d\sin \beta ={\mathcal{B}}_m(l)l,\end{array}$$

$$\mathcal{A}_{m}(l)=\left\{\begin{array}[]{lll}\displaystyle 1-\frac{2}{\xi^{2}}\left(1-\frac{E(\chi_{l},\,\xi)-E(\chi_{0},\,\xi)}{K(\chi_{l},\,\xi)-K(\chi_{0},\,\xi)}\right),\\[14.22636pt] \displaystyle 2\frac{E(\omega_{l},\,\eta)-E(\omega_{0},\,\eta)}{K(\omega_{l},\,\eta)-K(\omega_{0},\,\eta)}-1,\end{array}\right.\,\mathcal{B}_{m}(l)=\left\{\begin{array}[]{lll}\displaystyle\frac{2}{\xi^{2}}\frac{\sqrt{1-\xi^{2}\sin^{2}{\chi_{l}}}-\sqrt{1-\xi^{2}\sin^{2}{\chi_{0}}}}{K(\chi_{l},\,\xi)-K(\chi_{0},\,\xi)},&m=0,\\[14.22636pt] \displaystyle\frac{2\eta\left(\cos{\omega_{l}}-\cos{\omega_{0}}\right)}{K(\omega_{l},\,\eta)-K(\omega_{0},\,\eta)},&m\neq 0.\end{array}\right.$$

$$\delta^{2}\mathcal{V}=-B\int_{0}^{l}{\bigg{[}\theta_{var}^{\prime\prime}(s)+\gamma^{2}\,\theta_{var}(s)\,\cos{\big{(}\theta(s)-\beta\big{)}}\bigg{]}\,\theta_{var}(s)}\,\textup{d}s,$$

$${\int }_0^l{\theta }_{var}(s)\sin \theta (s)\text{d}s=0,{\int }_0^l{\theta }_{var}(s)\cos \theta (s)\text{d}s=0.$$

$${\varphi }_n(0)=0,{\varphi }_n(l)=0,{\int }_0^l{\varphi }_n(s)\sin \theta (s)\text{d}s=0,{\int }_0^l{\varphi }_n(s)\cos \theta (s)\text{d}s=0,$$

$${\varphi }_n^{\prime \prime }(s)+{\zeta }_{n}w(s){\varphi }_n(s)=C_3\sin \theta (s)+C_4\cos \theta (s),$$

Full text

11footnotetext: Corresponding author: Francesco Dal Corso fax: +39 0461 282599; tel.: +39 0461 282522; web-site: http://www.ing.unitn.it/$\sim$dalcorsf/; e-mail: francesco.dalcorso@unitn.it

Keywords: Elastica, snap-back, bifurcation, stability.

1 Introduction

Snap instability is a well-known phenomenon in mechanics for which a structural system suffers a sudden and dramatic change in the deformed configuration triggered by a small variation in the loading conditions. This behaviour is explained as the consequence of the stability loss for the deformed configuration, so that the structure dynamically moves towards a non-adjacent configuration through a partial release of its elastic energy. Classical examples of snap mechanisms can be found in shell structures, also in everyday life, as for instance when squeezing an empty can of soda. Other examples may be found in nature, as in the case of click-beetles (elateridae) [15], insects able to turn on their side when initially lying on their back by means of a jump realized by a snap mechanism.

Following the new paradigm of exploiting (instead of avoiding) instabilities in the structural design for advanced applications [18, 30], in the last years many researchers have investigated the snap mechanisms towards the realization of bistable or multistable devices [3, 10, 11, 12, 13, 14, 31, 33], metamaterials [16, 26, 29], locomotion [25, 34, 35, 36], and energy harvesting [17, 19]. Because the typical approach adopted in these investigations is to focus on specific evolutive mechanical problems and specific structural properties, the evaluation of the whole set of critical snap conditions is missed from a general perspective.

The aim of the present study is to provide a general criterion to be exploited in the nonlinear design of structures to obtain (or to avoid) snap mechanisms. With reference to the model of an inextensible (weightless) linearly elastic strip, the number of stable equilibrium configurations is disclosed by varying the parameters defining the kinematics of the strip ends, which are reduced to a normalized distance between the two ends and the two rotations of the ends. The stability of the equilibrium configurations, expressed in closed form through elliptic integrals, is assessed analyzing the sign of the second variation of the total potential energy at fixed ends conditions. This analysis allows us to define the correspondence between the ends’ kinematics and the presence of multiple stable configurations, so that the set of boundary conditions for which one of these configurations loses its stability follows.

The obtained results are used to define universal snap surfaces, collecting the whole sets of critical boundary conditions for which the system snaps. The energy release at snapping is also estimated and investigated by means of a dimensionless analysis for varying the snap conditions. The theoretical predictions, obtained analytically within a quasi-static framework, are experimentally validated through comparisons with data available for specific symmetric boundary conditions and observations on a physical model proposed for investigating non-symmetric boundary conditions. Finally, finite element simulations performed with ABAQUS show the reliability of the presented universal surface in the case of evolutions with moderate velocity and its limits in the case of very fast ends evolutions.

2 Governing equations and the planar elastica

An inextensible elastic strip of length $l$ is considered to be deformed within a plane orthogonal to the axis of minimum momenta of inertia of the strip’s cross section. The strip has a uniform cross section and is initially flat, so that its centerline is described by a straight line in the undeformed configuration. Neglecting the effects of self weight and disregarding rigid body motions, the mechanical fields along the generic curvilinear coordinate $s\in[0,l]$ can be represented in the local reference system $x-y$ (Fig. 1), with origin at one strip’s end ($s=0$), and $x$-axis passing by the other one ($s=l$), so that $x(0)=y(0)=y(l)=0$, and pointing from the initial to the final curvilinear coordinate. The primary kinematic field is the rotation angle $\theta(s)$, which measures the rotation of the structure’s centerline with respect to the $x$-axis, from which, considering the inextensibility constraint, the position fields can be obtained as

[TABLE]

so that the condition $y(l)=0$ implies the following isoperimetric constraint on the rotation field $\theta(s)$

[TABLE]

The strip is considered subject to kinematic boundary conditions in terms of position and rotation at both ends, which slowly move in time along quasi-static evolutions. Considering the $x-y$ reference system allows for the rational interpretation of the boundary conditions imposed at the two rod’s ends. Indeed, the six kinematic boundary conditions (two positions and one rotation at each end) affect the strip configuration by means of the following three primary kinematical quantities

[TABLE]

being $d$ the distance between the two clamps, $d\in[0,l]$ where the lower bound is given by the definition of the x-axis direction while the upper bound is related to the inextensibility assumption (although extensibility may affect the mechanical response even in the proximity of the limit condition, $d\simeq l$).

The planar behaviour of the considered strip is modelled as the Euler elastica so that the bending moment $M$ is given by $M(s)=B\theta^{\prime}(s)$ where the symbol ′ stands for the derivative with respect to the curvilinear coordinate $s$ and $B$ is the bending stiffness, constant because the strip’s cross section is uniform. Considering that the deformed configuration is described by imposing the distance $d$ and the two rotations $\theta_{0}$ and $\theta_{l}$, the total potential energy $\mathcal{V}$ under quasi-static conditions is given by

[TABLE]

where $\mathcal{E}$ is the elastic energy stored within the strip

[TABLE]

while the quantities $R_{x}$ and $R_{y}$ are the Lagrangian multipliers representing the reaction forces components at both ends along the $x$ and $y$ directions and, similarly, $M_{l}$ and $M_{0}$ are those referred to the rotational degrees of freedom of the strip’s coordinates $s=l$ and $s=0$.

The annihilation of the first variation of the functional $\mathcal{V}$ leads to the differential equation of the Euler elastica [7]

[TABLE]

being $\psi(s)=\theta(s)-\beta$ the auxiliary angle, $\gamma^{2}=R/B$ the normalized load parameter, $R=\sqrt{R_{x}^{2}+R_{y}^{2}}$ the total reaction force and $\beta$ its inclination with respect to the $x$-axis (Fig 1).

A first integration of the equation (6) leads to

[TABLE]

where $\Upsilon$ is the constant of integration to be evaluated with reference to the boundary conditions and defining the elastica as a part of a non-inflectional mother curve ($\Upsilon\in[1,\infty]$) or an inflectional mother curve ($\Upsilon\in[-1,1]$) [23].

In the following paragraphs, the analytical description of the deformed configurations is presented with reference to the number $m\in\mathbb{N}^{0}$ of ‘inflection points’, corresponding to the curvilinear coordinates within the set $s\in(0,l)$ where the curvature vanishes. It is remarked that points with null curvature located at both ends are not considered in the definition of $m$. For convenience, the constant of integration $\Upsilon$ is defined distinguishing the fundamental cases of the absence ($m=0$) and the presence ($m\neq 0$) of inflection points and the deformed configurations are described, respectively, in terms of the equations related to the so-called non-inflectional and inflectional elastica. Therefore it is worth to remark that, differently from Love [23], the deformed configurations with $m=0$ associated with an inflectional mother curve are here described for simplicity through the expressions used for the non-inflectional mother curve (for which the ideal unlimited elastica, $-\infty<s<\infty$, has no inflection points) but restricted to the physical range of the curvilinear coordinate, $s\in[0,l]$ (as shown in Sect. 3).

Absence of inflection points along the strip ($m=0$).

In this case $\psi^{\prime}(s)\neq 0$ for $s\in(0,l)$ and the integration constant $\Upsilon$ can be defined as a function of an unknown parameter $\xi$ as

[TABLE]

where the upper bound for the parameter $\xi$ restricts the solution at equilibrium to non-null curvature values within the set $s\in(0,l)$. It is noted that, differently from Love [23], the range of $\xi$ is extended to values higher than 1, which are related to inflectional mother curves.

Furthermore, an auxiliary rotation field is introduced as $\chi(s)=\psi(s)/2$ whose values at the two ends are defined as $\chi(0)=\chi_{0}$ and $\chi(l)=\chi_{l}$, given by

[TABLE]

Presence of $m$ inflection points along the strip ($m\neq 0$).

In this case $\psi^{\prime}(\widehat{s}_{j})=0$ at the curvilinear coordinates $\widehat{s}_{j}\in(0,l)$ with $j\in[1,m]$ ordering the inflection points with respect to their curvilinear coordinate, $\widehat{s}_{j}<\widehat{s}_{j+1}$ . Eqn (7) implies that the integration constant $\Upsilon$ is a function of the angle $\psi(\widehat{s}_{1})=\widehat{\psi}_{1}$ measured at the inflection point along the strip and closest to the origin as

[TABLE]

Being $\cos\widehat{\psi}_{j}=\cos{\widehat{\psi}_{1}}$ with $\psi(\widehat{s}_{j})=\widehat{\psi}_{j}$ ($j=1,...,m$), the rotation angle at the $j$-th inflection point is given by

[TABLE]

Because the curvature changes sign at each inflection point, the solution (7) for the curvature can be rewritten as

[TABLE]

where $p$ is a boolean parameter defining the curvature sign at the left end $s=0$, namely $p=0$ ($p=1$) when the curvature is positive (negative) at the left end, $\theta^{\prime}(s=0)>0$ ($\theta^{\prime}(s=0)<0$).111In the case of null curvature at the initial coordinate, $\theta^{\prime}(s=0)=0$, the boolean $p$ is defined by the sign of the curvature at positive infinitesimal values of the coordinate, $\theta^{\prime}(s=0^{+})$.

Towards the solution achievement it is instrumental to introduce the parameter $\eta\in\,[0,1]$ and the auxiliary field $\omega(s)$ defined as

[TABLE]

The values of the function $\omega(s)$ attained at the two ends, $\omega(s=0)=\omega_{0}$ and $\omega(s=l)=\omega_{l}$, represent the fundamental parameters involved in the problem resolution and depend on the parameters $\beta$ and $\eta$ as follows:

[TABLE]

3 The elasticae joining two constrained ends

A further integration of the differential equation (7) leads to the following relations for the normalized load $\gamma$

[TABLE]

and the rotation field

[TABLE]

In equations (15) and (16), $K$ is the Jacobi’s incomplete elliptic integral of the first kind, am is the Jacobi’s amplitude function, sn is the Jacobi’s sine amplitude function,

[TABLE]

It is remarked that, when the deformed configuration is associated with an inflectional mother curve, $\Upsilon\in[-1,1]$, but there are no inflection points along the strip, $m=0$, the analytical representation of the normalized load $\gamma$, eqn (15), and of the rotation field $\theta(s)$, eqn (16), for the case $m=0$ coincides with that for the case $m\neq 0$.222In addition to the general property for the reciprocal modulus transformation for the Jacobi’s sine amplitude function (Byrd [9], his eqn 162.01, pag 38)

$\textup{sn}\left(\sigma\varphi,\frac{1}{\varphi}\right)=\varphi\,\,\textup{sn}\left(\sigma,\varphi\right),$

(18) under the circumstance $\Upsilon\in[-1,1]$ and $m=0$ the following properties hold:

$\xi=1/\eta,\qquad\sin\omega(s)=\xi\sin\chi(s),\qquad K\left(\omega(s),\,\eta\right)=\xi K\left(\chi(s),\,\xi\right).$

(19)

From the rotation field (16), the elastic energy $\mathcal{E}$, eqn(5), is given by

[TABLE]

and the deformed shape at equilibrium can be evaluated from the position field (1) and expressed as

[TABLE]

with the functions $\mathcal{A}_{m}(s)$ and $\mathcal{B}_{m}(s)$ given by

[TABLE]

In eqns (22) the function cn is the Jacobi’s cosine amplitude function, $\mathscr{E}$ is the Jacobi’s epsilon function, and dn is the Jacobi’s elliptic function, defined as

[TABLE]

and $E$ is the Jacobi’s incomplete elliptic integral of the second kind

[TABLE]

Imposing specific boundary conditions at the strip’s ends, the elastica joining the two constrained ends can be found. The two isoperimetric constraints, $x(l)=d$ and $y(l)=0$, are enforced through the following nonlinear system:

[TABLE]

where, from equation (22), the quantities $\mathcal{A}_{m}(l)$ and $\mathcal{B}_{m}(l)$ are given by

[TABLE]

The nonlinearity of the system (25) leads to possible non-uniqueness for the equilibrium configuration, providing interesting conditions of bifurcation or snap instability during a quasi-static variation of the boundary conditions as shown in the next Section.

As a final remark, equations (8)-(22) can be directly exploited to attack various boundary value problems by properly imposing the nonlinear system to be solved, so that the present formulation can also be considered to treat structural systems as an alternative to previously adopted procedures [37].

3.1 Stability of the elastica with constrained ends

Extending previously presented procedures [4],[8], [20], the stability of the achieved equilibrium configurations is judged with reference to small and compatible perturbations in the rotation field $\theta(s)$. In particular, the stability is connected to the sign of the second variation of the total potential energy (4), which referring to a weak perturbation $\theta_{var}(s)$ can be written as

[TABLE]

so that an equilibrium configuration is stable when $\delta^{2}\mathcal{V}$ is positive definite ($\delta^{2}\mathcal{V}>0$ for every compatible perturbation $\theta_{var}(s)$), while unstable when $\delta^{2}\mathcal{V}$ is indefinite or negative definite. Differently, in the case that $\delta^{2}\mathcal{V}$ is semi-positive definite, higher order variations have to be considered (as described at the end of this subsection) to judge the stability.

The compatibility restricts the perturbation field $\theta_{var}(s)$ to have null values at the two ends, $\theta_{var}(0)=\theta_{var}(l)=0$, and to satisfy the following isoperimetric constraints (derived from the first variation of the displacement boundary conditions $x(l)=d$ and $y(l)=0$, eqns (3)1 and (2))

[TABLE]

Introducing the eigenfunction $\phi_{n}(s)$, subject to the boundary conditions

[TABLE]

the positive definiteness of the functional $\delta^{2}\mathcal{V}$, eqn (27), can be analyzed investigating the non-trivial solutions of the following Sturm-Liouville problem:

[TABLE]

where the weight function is defined as

[TABLE]

and where the resultant inclination $\beta$, the normalized load $\gamma$ and the rotation field $\theta(s)$ are known at this stage from the definition of equilibrium configuration, while $C_{3}$ and $C_{4}$ are Lagrangian multipliers, and $\zeta_{n}$ is the eigenvalue related to the eigenfunction $\phi_{n}(s)$. The following property holds:

[TABLE]

so that any compatible perturbation $\theta_{var}(s)$ can be expressed as the Fourier series expansion of the eigenfunctions $\phi_{n}(s)$

[TABLE]

where $c_{n}$ are the Fourier coefficients. Considering the properties (32) and the Fourier series expansion for the perturbation field $\theta_{var}(s)$, the second variation of the total potential energy (27) can be rewritten as

[TABLE]

so that the equilibrium configuration is stable whenever $\zeta_{n}\notin\left[0,1\right]$ for every $n$, unstable when $\zeta_{n}\in\left(0,1\right)$ for at least one value of $n$, and to be investigated through higher-order variations otherwise. It follows that the condition of smallest eigenvalue greater than one implies that the configuration is stable. Following Kuznetsov and Levyakov [22], the eigenfunction $\phi_{n}(s)$ can be represented as the linear combination of the functions $\varphi_{n}^{(j)}(s)$ ($j=1,...,4$)

[TABLE]

where the four functions are defined respectively as the solutions of the following four second order differential problems

[TABLE]

Considering the boundary condition (29)1, and the boundary conditions defined for each function $\varphi_{n}^{(j)}(s)$ ($j=1,...,4$) in the differential problems (36), it follows that $C_{1}=0$ so that the evaluation of the function $\varphi_{n}^{(1)}(s)$ can be disregarded. Imposing the remaining boundary conditions (29)2, (29)3, and (29)4, the homogenous linear problem $\textbf{A}(\zeta_{n})\textbf{C}=\textbf{0}$ is obtained for the vector $\textbf{C}=\left\{C_{2},C_{3},C_{4}\right\}$, where

[TABLE]

so that the eigenvalues $\zeta_{n}$ can be finally evaluated from the condition of vanishing determinant, $\det{\textbf{A}(\zeta_{n})}=0$. From the operational point of view, with reference to every equilibrium rotation field $\theta(s)$, the differential system (36) can be numerically solved as a function of $\zeta_{n}\in[0,1]$ and the configuration judged stable when $\det{\textbf{A}(\zeta_{n})}\neq 0$ for $\zeta_{n}\in[0,1]$, while unstable if this condition is not fulfilled. It is noted that in the limiting case when $\det{\textbf{A}(\zeta_{n})}$ vanishes only for $\zeta_{n}=1$, the related equilibrium configuration is unstable if the third variation is non-null for the related eigenfunction $\phi_{n}(s)$, while if this variation vanishes the analysis should consider higher-order variations. In particular, if the next even variation is positive (negative) for $\theta_{var}(s)=\phi_{n}(s)$, the configuration is stable (unstable). Otherwise, if the next even variation is null for $\theta_{var}(s)=\phi_{n}(s)$, the next odd variation should be considered and if not null the configuration is unstable. In general, the $k$-th variation of the total potential energy can be written for $k\geq 3$ as

[TABLE]

4 Number of stable solutions, bifurcations and universal snap conditions

For any triad of parameters $\left\{d,\theta_{0},\theta_{l}\right\}$, the nonlinear system (25) can be solved to numerically evaluate the existing pairs of parameters ($\beta$ and $\xi$ for $m=0$, or $\beta$ and $\eta$ for $m\neq 0$) associated with the possible stable equilibrium configurations for such boundary conditions. Reference is made solely to configurations with $m\in[0,2]$ because those with $m>2$ are numerically found to be always unstable (although a general analytical proof seems awkward, [23]). The number of stable solutions has been observed to vary from 1 to 3 within the kinematical parameter space defined by $\left\{d,\theta_{0},\theta_{l}\right\}$.333It is remarked that the end rotations $\theta_{0}$ and $\theta_{l}$ have no restriction on their value, being unlimited their difference $\theta_{l}-\theta_{0}=\int_{0}^{l}\theta^{\prime}(s)\textup{d}s$. If the angles $\theta_{0}$ and $\theta_{l}$ were referred to the end inclinations (instead of being referred to the end rotations, as in the present analysis), their respective sets would be limited due to angular periodicity, for example to $\theta_{0}\in(-\pi,\pi)$ and $\theta_{l}\in(-\pi,\pi)$, and more than 3 stable configurations could be found as solution to the same end inclinations problem [2].

Furthermore, being the developed model referred to the case of strips with loads applied only at the two ends, configurations of self-intersecting elastica are found through the presented numerical evaluation for a set of boundary conditions, while the respective self-contact configurations between different parts of the strip are excluded (see Section 4.4 for further details).

A typical map in the plane $\left[\theta_{0},\,\theta_{l}\right]$ representing the number of stable solutions is shown in Fig. 2 (left) for a distance $d=0.6l$, where the existence of one, two, and three stable configurations for a specific triad $\left\{d,\theta_{0},\theta_{l}\right\}$ is identified by the regions I, II, and III, respectively. The figure shows that the regions within this plane are characterized by a periodicity vector $[2\pi,2\pi]$ related to the shifting in the solution of $2j\pi$ ($j\in\mathbb{Z}$) when the imposed rotations are modified by $2j\pi$ at both ends,

[TABLE]

With reference to such periodicity property, it is instrumental to introduce the angles $\theta_{A}$ and $\theta_{S}$, respectively defined as the antisymmetric and symmetric parts of the imposed rotations,

[TABLE]

which are reported in Fig. 2 (left) through grey axes inclined at an angle $\pi/4$ with respect to the axes $\theta_{0}-\theta_{l}$. Considering the shifting of the rotation field expressed by eqn (39) it follows that

[TABLE]

highlighting the mentioned periodicity property of the equilibrium configurations which is given only in the $\theta_{A}$ variable within the $\theta_{A}-\theta_{S}$ reference system.444 The periodicity for the solution in the rotation field for the elastica is similar to that observed in the kinematic description of the physical pendulum, which is insensitive to an increase of an angle $2j\pi$ ($j\in\mathbb{Z}$).

Another interesting property is now pointed out. A change in sign for each of the two parameters $\theta_{A}$ and $\theta_{S}$ is related to a specific ‘mirror’ in the boundary conditions for the parameters $\theta_{0}$ and $\theta_{l}$, namely

[TABLE]

It follows that the respective deformed configurations of the cases $\mathrm{b}$, $\mathrm{c}$ and $\mathrm{d}$ can be obtained through the mirroring of the reference deformed configurations of the case $\mathrm{a}$, Fig. 2 (right, upper part). The mirroring properties can be summarized as:

• •

a change in sign for the parameter $\theta_{A}$ defines a configuration obtained as the mirroring with respect to the line orthogonal to that joining the two ends and passing at its center;

• •

a change in sign for the parameter $\theta_{S}$ defines a configuration obtained as two mirrorings, one with respect to the line joining the two ends and the other with respect to the orthogonal line at its mid point;

• •

a change in sign for both the parameters $\theta_{A}$ and $\theta_{S}$ defines a configuration obtained as the mirroring with respect to the line joining the two ends.

These mirroring properties allow a simplified representation of the number of stable solutions within the plane $|\theta_{A}|-|\theta_{S}|$, so that only the first quadrant is drawn and restricted to the condition $|\theta_{A}|\leq\pi$ because of the periodicity vector. In this way, the map of solutions number reported in Fig. 2 (left) for a distance $d=0.6l$ is represented in Fig. 3 (right), where the inflection points number $m$ related to each stable solution is specified as listed in the subscript within parentheses. As a further example, the representation of the number of stable solutions is also reported for a distance $d=0.3l$ in Fig. 3 (left). Lines crossing and lines bounding the regions I, II, and III are reported in Fig. 3, in particular:

• •

the grey thin lines crossing the regions define the transition for which the number $m$ of inflection points along the strip changes, while the number of stable solutions is kept constant. These lines can be defined imposing null curvature at one of the two ends (so that the strip has one hinged end while the other is a rotating clamp);

• •

the thick lines bounding the regions define the transition for which the number of stable solutions changes, so that they represent the critical condition of snap for one of the stable solutions. The lines are reported as thick orange and thick brown in order to identify the two possible different snap-back conditions (described in the following).

Because of the simple reference (clamped-hinged) structure to which are related, the grey lines can be found from straightforward computations. Indeed, for a given distance $d$, the relation $\theta_{S}=\theta_{S}(\theta_{A})$ tracing the grey line can be computed from the nonlinear system (25) by imposing the condition of null moment at one of the two strip’s ends, for example considering $\omega_{0}=(-1)^{p}\pi/2$ if the hinge is located at the coordinate $s=0$. Differently, drawing the orange and brown thick lines is a more complex task because they are related to the issue of stability loss for one of the possible equilibrium configurations. Such analysis requires a further representation for the solution domains which is now introduced.

4.1 Equilibrium paths for a fixed distance $d$

In the attempt to simplify the visualization of the solution domains for a fixed distance $d$, reference is made to the quasi-static evolutions of the rotations $\theta_{A}$ and $\theta_{S}$ starting from a stable configuration for null angles at both clamps, $\theta_{0}=\theta_{l}=0$ so that $\theta_{A}=\theta_{S}=0$ (namely, a shortening of the distance between the two clamps is imposed starting from the straight configuration). Only two stable solutions can be obtained for these ‘initial’ boundary conditions, both characterized by the same number of inflection points ($m=2$) but differing in the value of $p$, and representative of the two possible buckled configurations for a strip clamped at both ends and subject to a shortening $l-d$. Considering this, the solution maps of Fig. 2 and of Fig. 3 (right) for a fixed distance $d=0.6l$ can be split in two separate representations, each of these restricted to a specific value of $p$ (identifying the positiveness, $p=0$, or negativeness, $p=1$, of the curvature at the coordinate $s=0$) as shown in Fig. 4 and therefore related to one of the two initial configurations possible for $\theta_{A}=\theta_{S}=0$. In this way, the generic evolution for the boundary conditions can be represented with a continuous curve in the plane $\theta_{A}-\theta_{S}$ starting from its origin. Despite the non-linearity of the problem, each point along this curve is related to a (when existing) unique deformed configuration whenever snap conditions or bifurcation points are not encountered along the considered evolution. Consequently, such a representation provides a fundamental tool in the definition of the solution domains and of the snap-back conditions for every possible evolution of the imposed rotations at a fixed distance $d$, showing the critical and bifurcation conditions depending on the value of $p$.

The moment-rotation response curve associated with each evolution of the boundary conditions reveals the presence of snap-back instabilities (related to the annihilation of the second variation, $\delta^{2}\mathcal{V}=0$) when a point with vertical tangent in the response curve is reached and no smooth stable555In general, the equilibrium path smoothly continues after the snap point with an opposite change of the rotation value. However, this path is not considered here due to its unstable nature. evolutions of the deformed configuration can be obtained for a further monotonic variation in the rotation. The set of boundary conditions $\left\{d,\theta_{A},\theta_{S}\right\}$ corresponding to the vertical tangency in the moment-rotation response are identified by means of a standard bisection algorithm applied to the analyzed equilibrium path.

The above described analysis provides the domains and lines as reported in Fig. 4, where also specific equilibrium configurations are displayed for some critical pairs of boundary conditions. The solution domains are reported with different colors, identifying different properties as follows:

• •

blue/green/red regions – the unique stable equilibrium configuration corresponding to the considered $p$ has zero, one, and two inflection points along the strip for the blue, green, and red regions, respectively;

• •

white regions – no stable equilibrium configuration is possible for that specific value of $p$. However, for the mirroring property, the stable equilibrium configuration exists for the other value of $p$;

• •

grey regions – the equilibrium configuration related to a pair $\overline{\theta}_{A},\overline{\theta}_{S}$ belonging to this region can be represented by that related to the pair within the colored or white regions (namely, outside the grey regions) considering the shifting of the solution as highlighted by eqn (41). The corresponding dual values can be evaluated from $\overline{\theta}_{A}$ and $\overline{\theta}_{S}$ as

[TABLE]

where the symbol $\lfloor\cdot\rfloor$ represents the floor function, which evaluates the greatest integer that is less than or equal to the relevant argument. Note that, from eqn (43) follows that $\theta_{A}\in[-\pi,\pi]$ for every $\overline{\theta}_{A}$.

Similarly, the lines separating the different domains are drawn with various styles, representing the type of transition in the equilibrium configuration occurring when the evolution of the boundary conditions passes from one region to another.

With reference to any continuous evolution in the set of applied rotations $\left\{\theta_{A}(\tau),\,\theta_{S}(\tau)\right\}$ from an initial ($\tau=\tau_{i}$) to a final instant ($\tau=\tau_{f}$), where $\tau$ is a time-like parameter, the following cases are possible:

• •

if a dashed grey line is crossed, the equilibrium configuration smoothly varies changing the number $m$ of the inflection points, while the value of $p$ is kept constant;

• •

if a continuous grey line is crossed, the equilibrium configuration smoothly varies changing both the number $m$ of the inflection points and the value of $p$, in particular $p(\tau_{f})=1-p(\tau_{i})$. Therefore, the final configuration is related to a solution map dual to that of the initial value of $p(\tau_{i})$. As it can be noted in Fig. 4, the continuous grey lines are present only at the borders between one of the colored regions with a white region;

• •

if a dashdotted grey line is crossed, the equilibrium configuration smoothly varies keeping fixed both the number $m$ of the inflection points and the value of $p$. As it can be noted in Fig. 4, these lines appears only at $\theta_{S}=\pm\pi$ and are connected to the shifting property of the solution, eqn (41), so that, considering eqn (43), the final equilibrium configuration has to be referred to the values $\left\{\theta_{A}(\tau_{f})-2\pi\lfloor 1/2+\theta_{A}(\tau_{f})/(2\pi)\rfloor\right.,$ $\left.\theta_{S}(\tau_{f})\right\}$;

• •

if a thick orange line is crossed, a critical configuration of snap-back type 1 instability is encountered so that a small variation in the boundary conditions realizes a large variation in the equilibrium configuration. Therefore, the final equilibrium configuration is achieved by means of a dynamic motion from the initial one. The initial and final configurations are described by a different value of $p$, namely, $p(\tau_{f})=1-p(\tau_{i})$. From the present analysis it is also observed that before and after the snap-back type 1 mechanism there are always two inflection points, $m(\tau_{f})=m(\tau_{i})=2$;

• •

if a thick brown line is crossed, the critical condition of snap-back type 2 instability is encountered. However, differently from crossing the thick orange line, in this case the boundary conditions of the final configuration $\left\{\theta_{A}(\tau_{f})/(2\pi),\theta_{S}(\tau_{f})\right\}$ lie within the grey region, so that the interpretation of the considered solution map requires a further effort. Indeed, from the solution shifting principle, the final configuration should be referred to the values $\left\{\theta_{A}(\tau_{f})-2\pi\lfloor 1/2+\theta_{A}(\tau_{f})/(2\pi)\rfloor,\theta_{S}(\tau_{f})\right\}$. If these values correspond to a white region for the solution map with $p(\tau_{i})$, the final configuration is the one associated with the same values but related to the dual map for which $p(\tau_{f})=1-p(\tau_{i})$;

• •

if no line is crossed, the initial and final configurations have the same values for $p$ and $m$, and therefore are represented within the same solution map corresponding to $p(\tau_{i})=p(\tau_{f})$.

It is worth to remark that crossing a colored thick line does not always provide a snap-back instability. Indeed, this phenomenon is strictly related to the equilibrium configuration taken by the structural system before crossing this condition. More specifically, referring to the $\theta_{A}-\theta_{S}$ plane and a fixed distance $d$, the snap mechanism is realized whenever the scalar product between the incremental vector connecting the initial to the final boundary conditions and the normal (defined as the derivative of the tangent) to the snap-back curve is non-negative.

Moreover, the snap-back curves display the typical shape of catastrophic cusps [5], revealing how the magnitude of the so-called control parameter at a critical point decreases with respect to the perfect case (maximum critical rotation at the cusp) due to the presence of increasing imperfections. In particular, the role of such control parameter is played by the symmetric part of the rotation $\theta_{S}$ for snap-back curves of type 1 and by the antisymmetric part $\theta_{A}$ for snap-back curves of type 2, while the role of the imperfection is respectively played by $\theta_{A}$ and $\theta_{S}$.

It is finally noted that all the regions and lines reported in Fig. 4 satisfy the mirroring properties shown in eqns (42) and highlighted in Fig. 2 (right). Therefore, when $p$ is switched from 0 to 1, the domains of the case $p=1$ (Fig. 4, right) can be obtained from those of the case $p=0$ (Fig. 4, left) through a mirroring with respect to the axis $\theta_{S}$ for the configurations with $m=1$ and through a double mirroring (with respect to both axes, $\theta_{A}$ and $\theta_{S}$) for those with $m=0$ and $m=2$.

Similarly to Fig. 4, the domains are reported in Figs. 5 and 6 for different values of the end’s distance, respectively for $d=\left\{0.05,0.1,0.2,0.3\right\}l$ and $d=\left\{0.4,0.5,0.7,0.8\right\}l$. The maps are reported only for $p=0$, since those are related to the case $p=1$ through mirroring.

A last comment is made about the two limit cases of maximum and minimum distance between the ends, respectively $d\simeq l$ and $d=0$:

• •

In the former case ($d\simeq l$), although the analysis should be improved by considering stretching energy, the present model (based on the inextensibility assumption) shows that the curve of snap-back type 1 reduces to the point with coordinates $\theta_{A}=\theta_{S}=0$ while the curve of snap-back type 2 reduces to the four line segments defined by $|\theta_{S}|=2\pi-|\theta_{A}|$ with $|\theta_{A}|\in[\pi,2\pi]$;

• •

Differently, in the case when the two ends have the same position ($d=0$), the mechanical system shows the independence of the parameter $\theta_{A}$, because in the case of null distance the angle $\theta_{A}$ merely expresses a rigid rotation of the entire structure and the elastic energy stored within the structure is only a function of the angle $\theta_{S}$. In this case, the curve of snap-back type 1 becomes the segment lines given by $\theta_{S}=0$ and $|\theta_{A}|\lesssim 0.726\pi$ while the curve of snap-back type 2 becomes the segment lines given by $\theta_{S}=0$ and $|\theta_{A}|\in[1,1.726]\pi$. It is also worth to mention that in the very special case of $\theta_{S}=0$, an infinite set of stable and equivalent (namely, corresponding to the same elastic energy) solutions exists for the strip with fixed end rotations, provided by a ‘8-shaped’ configuration [21, 32].

4.2 Equilibrium paths for a variable distance $d$

The families of snap curves reported for fixed values of $d$ in Figs. 4, 5 and 6 suggest the existence of snap surfaces within the space $d/l-\theta_{A}-\theta_{S}$. These snap surfaces can be disclosed by interpolating the snap curves within specific planes. The snap curves evaluated for twenty planes, taken for computational convenience at constant $\theta_{A}$ and $\theta_{S}$ for snap types 1 and 2, respectively, have been exploited to generate the snap-back surfaces by means of the function Interpolation in Mathematica (v.10). Due to mirroring properties, the snap surfaces are entirely described through their representations within one octant of the space $d/l-\theta_{A}-\theta_{S}$, displayed in Fig. 7 for snap-back type 1 (left) and type 2 (right) surfaces. Such snap surfaces represent a fundamental tool in the investigation about snap mechanisms and the definition of the correspondent critical boundary conditions during a quasi-static evolution of strips with controlled ends. With the purpose to facilitate the use of this tool, the numerically evaluated dataset of about 2000 critical conditions for snap type 1 is made available as Supplementary material. Before exploiting the concept of universal snap surfaces in some applicative examples presented in Sect. 5, it is worth to discuss the possibility of bifurcations during a loading process, the amount of elastic energy release, and the conditions of self-intersecting elastica.

4.3 Bifurcations at snap

The boundary conditions corresponding to (almost all) the intersections of the orange and brown snap-back curves with the lines $\theta_{A}=\left\{-\pi,0,\pi\right\}$ and with the line $\theta_{S}=0$ (corresponding to the cusps and end points of the snap curves) are marked with colored spots in Figs. 4, 5, and 6. For the boundary conditions corresponding to these points, the semi-positive definiteness of the second variation $\delta^{2}\mathcal{V}$ is complemented by the annihilation of the third variation $\delta^{3}\mathcal{V}$ for the related eigenfunction (while $\delta^{3}\mathcal{V}$ is different from zero for the other points along the snap curve) and a positive value is found for the corresponding fourth variation $\delta^{4}\mathcal{V}$. These special points (some of them also observed in [24]) are marked as red and yellow spots, respectively corresponding to pitchfork (with no snapping) and unstable-symmetric bifurcation points (with snapping), and defined as follows.

Red spots.

These points are located at the intersections of the orange thick curve with the condition $\theta_{S}=0$ and of the brown thick curve with the condition $\theta_{A}=\pm\pi$. Snap never occurs when the snap curve is crossed through these points; however when the snap curve is crossed from outside to inside, a bifurcation may be encountered. More specifically, a pitchfork bifurcation is found when crossing the orange (brown) curve from outside to inside at the red point with a null increment in the angle $\theta_{S}$ ($\theta_{A}$), namely a purely antisymmetric (symmetric) variation in the boundary conditions at the critical point. These bifurcative behaviours are displayed in Fig. 8 and Fig. 9, where specific loading processes are considered for a fixed distance $d=0.4l$ and crossing the red spots located at the edges of the snap-back type 1 and type 2 curves.

The antisymmetric case, $\theta_{S}=0$, at increasing $\theta_{A}$ is considered in Fig. 8, where in the first line the loading path is reported on the left and the moment $M_{0}$ at the left clamp as a function of the angle parameter $\theta_{A}$ on the right. The evolution of the deformed shape is displayed in the second, third, and fourth lines and is reported for the six specific stages of the boundary conditions, highlighted both in the solution maps and moment-angle response through the alphabetic letters $a$, $b$, $c$, $d$, $e$, $f$, $g$, $h$, and $i$. During the evolution, the bifurcation occurs when the stage $c$ is attained, namely, when the snap curve is crossed from outside to inside, and corresponds to the condition of null moment at both clamps. Just after the stage $c$ is passed, the structure may equally evolve through two different loading paths, namely, the structure may equally reach configuration $d1$ or $d2$. Once that one of these two branches is undertaken, the evolution continues on that specific branch. However, at increasing the rotation parameter, both branches finally join together when the stage $g$ is reached and, after this stage (from $g$ to $i$), the structure follows the only possible evolution. It is also important to highlight that the configurations $c$ and $g$ are stable solutions having $m=1$ internal inflection point and null curvature at both ends (similar configurations are also present for all the bifurcation conditions, red spots, on snap-back curve type 1 for every distance $d$, and for some bifurcation conditions, yellow spots, on snap-back curve type 2 as discussed below).

A similar behaviour is also displayed in Fig. 9 with reference to $\theta_{A}=-\pi$ and decreasing $\theta_{S}$. Differently from Fig. 8, here the bifurcation occurs for a symmetric configuration and is associated with a non-null moment value at both clamps. The pitchfork bifurcation at the point $c$ reveals a symmetry-breaking behaviour. In fact, the initial symmetric configuration becomes unstable after the bifurcation point is reached, and the system may only follow two other stable paths described by the mirrored configurations. A similar behaviour has also been detected in the case of a ring pinched by two radial loads [22], while the existence of a central unstable and symmetric solution has also been reported by [7] and [22], where the symmetric configuration for the double-clamped rod with null rotations at its ends is proven to snap towards the S-shaped configurations $e1$ or $e2$.

Yellow spots.

These points are located at the intersections of the orange thick curve with the condition $\theta_{A}=0$ and, only for $d\leq 0.369\,l$, of the brown thick curve with the condition $\theta_{S}=0$.666It is remarked that the yellow points given by the intersection of the brown snap-back curve with the axis $\theta_{S}=0$ exist only for $d\leq 0.369\,l$ and correspond to stable configurations with $m=1$ inflection point and null curvature at both ends. This behaviour is not observed for $d>0.369\,l$, where the third variation $\delta^{3}\mathcal{V}$ is not null and snap mechanism occurs without any bifurcation as soon as the snap curve is crossed at any other point. It follows that no yellow point appears in Fig. 6 for snap-back curves of type 2, being the considered distances $d>0.369\,l$. Snap occurs when the snap curve is crossed through these points from inside to outside and a bifurcation may be even encountered. More specifically, an unstable-symmetric bifurcation is found at the snap when crossing the orange (brown) curve from inside to outside at the yellow point with a null increment in the angle $\theta_{A}$ ($\theta_{S}$), namely a purely symmetric (antisymmetric) variation in the boundary conditions at the critical point. The moment rotation response at increasing value of $\theta_{S}$ crossing the snap type 1 curve is displayed in Fig. 10 for a fixed distance $d=0.4l$ and at very small fixed values of $\theta_{A}=\{0,10^{-4},10^{-3},10^{-2}\}$. An example of unstable-symmetric behaviour can be envisaged in the purely symmetric case, $\theta_{A}=0$, where the symmetric response (for which the moments at the two clamps have the same value, $M_{0}=M_{l}$) intersects the two other unstable paths. Differently, if a small constant value is assumed for the antisymmetric part of rotation $\theta_{A}$, the moment rotation response reaches a critical condition of snap-back, for which the tangent of the response curve is vertical.

4.4 Self-intersecting elastica

The snap surfaces reported in Fig. 7 have been obtained assuming that the strip is only loaded at its ends. It follows that the obtained critical conditions hold whenever the development of self-contact points along the strip is excluded.777Analysis of elasticae with self-contact points requires the resolution of two or more elasticae subject only to end loadings [27]. This circumstance is trivially realized when the deformed configuration is not self-intersecting, but also when the self-intersection is made possible by the out-of plane geometry. The latter case is realized with strips shaped along the out-of-plane direction in such a way that during the planar intersection two external halves of the strip contain a central strip, namely a Y-shaped strip with specific out-of-plane variations [8, 22].

In order to detect when self-intersection does occur before snapping and, equivalently, when it does not, it is of practical interest to define the boundary conditions for which a contact point is first realized. This information can be consequently used to determine which portions of snap surfaces are attained only after an evolution involving the self-intersection. The boundary conditions of first self-contact can be found imposing that for one and only one pair of curvilinear coordinates $s^{(1)}$ and $s^{(2)}$ have the same position,

[TABLE]

Evaluating the conditions of first self-contact, it is observed that:

• •

the conditions of snap-back type 1 (orange surface, Fig. 7 left) are always reached without developing a self-intersecting elastica if $d\gtrsim l/4$. For $d\lesssim l/4$, the snap-back most likely occurs after developing a self-intersecting configuration, however the exact limit distance for which the self-intersection is realized depends on the values of $\theta_{A}$ and $\theta_{S}$, Fig. 11;

• •

the conditions of snap-back type 2 (brown surface, Fig. 7 right) are always reached after developing a self-intersecting elastica.

It is finally observed that self-contact may occur during snapping although the strip has no self-contact at the critical snap condition. Indeed, the dynamic transition from the pre and post snap configurations may evolve requiring self-intersecting shapes, which could be not geometrically feasible even for Y-shaped strips.

4.5 Energy release at snapping

The design of snapping devices can be optimized through the maximization of the energy release during mechanisms. In a first approximation, the energy release $\Delta\mathcal{E}$ can be estimated as the difference in the energy amounts $\mathcal{E}$, eqn (5), between the stages before and after the snap-back and evaluated under the quasi-static assumption through eqn (20).

Restricting attention only to snap-back type 1, the elastic energy difference $\Delta\mathcal{E}$ is reported in Fig. 12. Two nondimensionalizations are considered, division by $B/l$ (Fig. 12, upper part) and division by the elastic energy before snapping $\mathcal{E}_{0}$ (Fig. 12, lower part). On the left part of Fig. 12, the elastic energy difference $\Delta\mathcal{E}$ is shown for different distances $d/l=\left\{\right.0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.7,$ $0.8,$ $0.9$ $\left.\right\}$ for varying the angle for which the snap-back type 1 occurs. The angle measure is expressed as the modulus of $\theta_{A}^{sb1}$, the antisymmetric critical angle at the snap-back, normalized through division by $\theta_{A}^{sb1,bif}$, the antisymmetric angle at bifurcation for snap-back (the red spots on the $\theta_{A}$ axis in Fig. 4, 5 and 6 and for which no snap occurs, $\theta_{A}^{sb1,bif}(d)=\max_{\theta_{S}}\left\{\theta_{A}^{sb1}(\theta_{S},d)\right\}=\theta_{A}^{sb1}(\theta_{S}=0,d)$). For all the reported cases it can be concluded that the maximum elastic energy release for a fixed distance $d$ is always attained under symmetric conditions, $\theta_{A}^{sb1}=0$. For such symmetric condition, the elastic energy difference $\Delta\mathcal{E}$ is shown on the right part of Fig. 12 as a function of the distance $d$ between the strips’ ends. The plot on the upper right part shows that the energy release $\Delta\mathcal{E}$ has a maximum value of about $76.67B/l$ at $d\approx 0.189l$ and has null values for both the limit-cases $d=l$ and $d=0$. The plot on the lower right part shows that the relative energy release $\Delta\mathcal{E}/\mathcal{E}_{0}$ monotonically increases with the increase of the distance $d$ and attains its maximum ratio of about $0.889$ in the limit condition of $d=l$. With reference to this last case, it is also worth to highlight that the relative energy release is approximately constant for $d>l/4$, more specifically it varies from $\Delta\mathcal{E}/\mathcal{E}_{0}(d=l/4)=0.884$ to $\Delta\mathcal{E}/\mathcal{E}_{0}(d\rightarrow l)=0.889$.

5 Validation of the analytical predictions

The universal critical surface for snap type 1 is validated through comparison with experimental observations and results from numerical simulations. The available experimental data [6],[28], restricted to symmetric boundary conditions ($\theta_{A}=0$), are complemented by testing a physical model developed to cover non-symmetric paths ($\theta_{A}\neq 0$). Finally, the influence of dynamical effects on the system response is assessed through numerical simulation of evolutive problems performed in ABAQUS. Reliability of the quasi-static predictions obtained through the universal surface is shown in the case of evolutions with moderate velocity.

5.1 Experimental results

A physical model (Fig. 13, left) is developed in order to experimentally investigate the snap conditions of the structural system under non-symmetric paths. Two forks in steel are exploited to control in practice the clamps, constraining the position and the rotation at both ends of an elastic strip. The strip is obtained from cutting a transparency film by Folex and has cross section 12 mm width times 0.1 mm height and length $l=200$ mm. Restricting the kinematics of the two forks, specific non-symmetric paths are covered with the developed device. More specifically, the fork constraining the left end has a fixed position while the fork constraining the right end imposes null inclination, $\theta_{l}=0$, and may move along the $x$-axis. From these boundary conditions it follows that the evolution expressed in terms of the three main kinematic quantities is given by

[TABLE]

so that the non-symmetric paths characterized by $\theta_{A}(\tau)=-\theta_{S}(\tau)=\theta_{0}(\tau)/2$ can be investigated. Varying only one kinematical parameter, the two following types of experiments are performed:

1. Exp. A

– keeping a fixed distance $d(\tau)=\overline{d}$ between the two forks, the rotation at the left end $\theta_{0}(\tau)$ changes in time; 2. Exp. B

– keeping a fixed rotation at the left end $\theta_{0}(\tau)=\overline{\theta}_{0}$, the distance $d(\tau)$ between the two forks changes in time.

During each experiment, the rotation $\theta_{0}(\tau)$ or the distance $d(\tau)$ is slowly varied by hand. The variation in the kinematical parameter is stopped as soon as the strip snaps and the respective critical value is measured for the rotation $\theta_{0,cr}(\overline{d})$ in Exp. A or for the distance $d_{cr}(\overline{\theta_{0}})$ in Exp. B with the goniometer or the ruler mounted on the device, respectively. The critical conditions experimentally collected from Exp. A and Exp. B are respectively reported as dots and crosses within the plane $\theta_{0}-d/l$ in Fig. 13 (right) together with the theoretical critical curve, namely the intersection of the universal snap surface type 1 (Fig. 7, left) with the plane $\theta_{A}=-\theta_{S}$.

These experimental measures are also reported together with those measured in the case of purely symmetric rotation $\theta_{A}=0$ by Beharic et al. [6] (for $\overline{d}/l\backsimeq\left\{0.833,0.862,0.893,0.926,0.962,0.980,0.990\right\}$) and by Plaut and Virgin [28] (for $\overline{d}/l\backsimeq\left\{0.413,0.637,0.827,0.955\right\}$) in Fig. 14, where the universal snap surface type 1 and its intersection with planes at constant values of $d/l$, $|\theta_{A}|$, and $|\theta_{S}|$ are represented, fully confirming the reliability of the predictions from the present model.

5.2 Numerical simulations and dynamic effects

In order to evaluate the possible influence of inertia on the snapping conditions of the system, the presented quasi-static predictions are finally compared with the response obtained from the numerical simulations performed in ABAQUS (v.6.13) for two different evolutive problems.

1. Sim. I

– Clamps rotating cyclically and with opposite velocity. Considering a fixed distance $\overline{d}$, the kinematics for the strip’s ends is described by

[TABLE]

or equivalently, through eqn (40), in terms of symmetric and antisymmetric parts of the angles as

[TABLE]

from which it is evident that the function $\vartheta(\tau)$ represents the evolution of the symmetric part of the two angles in the time-like parameter $\tau$ while $\overline{\theta}_{A}$ represents a constant anti-symmetric angle during the evolution. Referring to the physical time $t=T\tau$, where $T=\sqrt{\rho l^{4}/B}$ is the characteristic time for the system having a linear mass density $\rho$, results from Sim. I are reported for a fixed antisymmetric rotation $\overline{\theta}_{A}=0.16\pi$ and two fixed distances $\overline{d}/l=0.4$ (Fig. 15, upper part, left) and $\overline{d}/l=0.8$ (Fig. 15, lower part, left). In both cases, the cyclic evolution in the boundary conditions is realized through the succession of the increase and decrease in the symmetric rotation within the range $\theta_{S}\in[-\pi/2,\pi/2]$ keeping a constant modulus in the velocity,

[TABLE]

where a superimposed dot corresponds to the derivative in the physical time $t$ and $\Omega$ is the dimensionless (angular) velocity. The cyclic path prescribed in Sim. I theoretically encounters two snap conditions given by the same modulus of $\theta_{S}$. 2. Sim. II

– Monotonic variation in the clamps distance. Constant rotations $\overline{\theta}_{0}$ and $\overline{\theta}_{l}$ are assumed for the two ends, so that the kinematical evolution of the strip’s ends is given by

[TABLE]

where the function $d(\tau)$ defines a monotonic increase or decrease in the distance between the two clamps. Results from Sim. II are reported for fixed rotations $\overline{\theta}_{0}$ and $\overline{\theta}_{l}$ such that $\overline{\theta}_{S}=\pi/8$ (Fig. 15, upper part, right) and $\overline{\theta}_{S}=\pi/4$ (Fig. 15, lower part, right), and in both cases $\overline{\theta}_{A}=0.16\pi$. In both cases, the monotonic variation in the distance $d$ is considered from the value $d(\tau=0)/l=0.336$ with a constant velocity,

[TABLE]

where $\Delta$ is the dimensionless velocity. During both the monotonic shortening and the monotonic lengthening in the clamps distance prescribed in Sim. II, a snap mechanism is theoretically predicted for each evolution.

A linear viscous Rayleigh damping acting on the mechanical system, modeled as 100 planar beam elements with linear elastic constitutive behaviour, is considered in all the simulations through the mass-proportional and the stiffness-proportional damping coefficients respectively as $A_{d}=8.25\times 10^{-3}/T$ and $B_{d}=6.06\times 10^{-3}T$. The inherent extensibility of the strip modeled in ABAQUS has been considered through the axial stiffness $EA=10^{6}\times B/l^{2}$. All the presented analyses are performed using the nonlinear geometry option and started from the undeformed straight configuration with null rotations, $d=l$ and $\theta_{0}=\theta_{l}=0$. All the simulations share the first two static steps, while are different in the last dynamic step as follows:

1. Step 1

– Static: An end’s distance $d(0)$ is imposed and a transversal load is applied in order to achieve the buckled configuration; 2. Step 2

– Static: the transversal load is removed and the clamp rotations are imposed in order to set the initial values of $\overline{\theta}_{A}$ and $\theta_{S}(0)$; 3. Step 3

– Dynamic implicit: inertial effects are analyzed during the evolution in the boundary conditions at velocity with constant modulus from $t=0$ to $t=t_{f}$ as

1. Sim. I

– initial and final configurations have $\theta_{S}(0)=\vartheta(0)=\vartheta(t_{f})=0$. Introducing the reference time $T_{r}=\pi T/(2\Omega)$, the duration of the evolution is given by $t_{f}=4T_{r}$. The rotation velocity is assumed $\dot{\vartheta}=\Omega/T$ for $t\in[0,T_{r}]$ and $t\in[3T_{r},4T_{r}]$, while is assumed $\dot{\vartheta}=-\Omega/T$ as for $t\in[T_{r},3T_{r}]$, so that the velocity changes sign whenever the modulus of rotation reaches $|\theta_{S}|=\pi/2$; 2. Sim. II

– the initial and final configurations have the ends distance given by $d(0)$ and $d(t_{f})$. The transition between these two conditions occurs with constant velocity and has a duration $t_{f}=T_{r}$, where $T_{r}=0.637T/\Delta$. The constant velocity is taken positive or negative in order to investigate the snap mechanism during lengthening or shortening.

The results of the simulations of the two evolutive problems are reported in Fig. 15 and compared with the respective quasi-static behaviour predictions using the inextensible elastica (and highlighting the snap angles and snap distances, defined by the snap curve $\theta_{A}=0.16\pi$ in Fig.14 right, lower part). In particular, results for the Sim. I are reported in Fig. 15 (left) in terms of the moment (at the right end) at cyclically varying (symmetric) rotation $\theta_{S}$, while those for the second evolutive problem are reported in Fig. 15 (right) in terms of the horizontal reaction (at the right end) for monotonic increase and decrease in the clamps distance. The respective numerical results are reported for three values of dimensionless velocities, $\Omega=\left\{0.01,\,0.1,\,1\right\}$ and $\Delta=\left\{0.001,\,0.01,\,0.1\right\}$, showing that the quasi-static model accurately describes the mechanical behaviour of the structural system until approaching the snap-back conditions, identified as the intersection of the loading path with presented snap-back surface type 1. Due to the presence of dissipative effects, the post-snap quasi-static path is reached after a transient time from the snap for which the dynamical effects are decayed. More specifically, non-negligible dynamic effects lead to a delay in the occurrence of snap with respect to the quasi-static prediction for high velocities. Oppositely, the dynamic response becomes almost completely superimposed (except for a small transient) to the quasi-static curve in the case of a velocity $\Omega=0.01$ and $\Delta=0.001$, values defining the velocity orders below which the present model, although obtained within the quasi-static framework, fully represents a reliable model.

6 Conclusions

Within a quasi-static framework, the number of stable equilibrium configurations has been disclosed for an elastic strip for varying the parameters controlling the kinematics of its ends. This analysis has led to the definition of universal snap surfaces, collecting the critical values of ends distance and rotations for which the strip shows a snap mechanism. Available experimental data and experimental results from testing a developed physical model fully validate the presented theoretical universal snap surface. Finally, finite element simulations show the influence of inertia on the snapping mechanisms and, in the case when the controlled ends move moderately, confirm the theoretical predictions based on the present quasi-static model. These results are complemented by the dimensionless analysis of elastic energy release at snapping, towards the optimal design of impulsive motion. In addition to the relevant contribution to the stability of structures, the present results may find application in a wide range of technological fields, ranging from energy harvesting to jumping robots.

Acknowledgments The authors gratefully acknowledge financial support from the ERC advanced grant ‘Instabilities and nonlocal multiscale modelling of materials’ ERC-2013-AdG-340561-INSTABILITIES (2014-2019).