# Delayed buckling of spherical shells due to viscoelastic knockdown of the
critical load

Lucia Stein-Montalvo<EMAIL_ADDRESS>Department of Mechanical Engineering,
Boston University, Boston, MA, 02215.  Douglas P. Holmes<EMAIL_ADDRESS>Department of Mechanical Engineering, Boston University, Boston, MA, 02215.
Gwennou Coupier<EMAIL_ADDRESS>Université Grenoble
Alpes, CNRS, LIPhy, F-38000 Grenoble, France.

###### Abstract

We performed dynamic pressure buckling experiments on defect-seeded spherical
shells made of a common silicone elastomer. Unlike in quasi-static
experiments, shells buckled at ostensibly subcritical pressures (i.e. below
the experimentally-determined load at which buckling occurs elastically),
often following a significant time delay. While emphasizing the close
connections to elastic shell buckling, we rely on viscoelasticity to explain
our observations. In particular, we demonstrate that the lower critical load
may be determined from the material properties, which is rationalized by a
simple analogy to elastic spherical shell buckling. We then introduce a model
centered on empirical quantities to show that viscoelastic creep deformation
lowers the critical load in the same predictable, quantifiable way that a
growing defect would in an elastic shell. This allows us to capture how both
the critical deflection and the delay time depend on the applied pressure,
material properties, and defect geometry. These quantities are straightforward
to measure in experiments. Thus, our work not only provides intuition for
viscoelastic behavior from an elastic shell buckling perspective, but also
offers an accessible pathway to introduce tunable, time-controlled actuation
to existing mechanical actuators, e.g. pneumatic grippers.

##  I Introduction

Shell structures are lightweight and flexible. Largely owing to their
curvature, they offer considerable strength with little material. As a result,
shells are abundant in nature (e.g. eggshells and blood vessels) and design
(e.g. fuel tanks and soda cans). However, slenderness also brings
susceptibility to abrupt and often catastrophic deformations. Clearly,
understanding how a thin, curved structure will lose stability – and in
particular at what load value this will occur – is crucial.

We restrict our attention to spherical shells in the present work. The first-
known quantitative prediction for the critical load $P_{c}$ in a perfect,
spherical, elastic shell subjected to uniform pressure was produced by Zoelly
in 1915 via linear eigenvalue analysis, and is given as:

$P_{c}=\frac{2E}{\sqrt{3(1-\nu^{2})}}\eta^{-2}$ (1)

for a shell with Young’s modulus $E$, Poisson’s ratio $\nu$, and radius ($R$)
to thickness ($h$) ratio $\eta\equiv R/h$.

Although this result is still widely accepted today, it severely overpredicts
the buckling load observed in experiments. Recognizing this disrepancy, which
is due to the extreme sensitivity to imperfections inherent to thin shells,
scientists at the space agency NASA and collaborators introduced the knockdown
factor ($k_{d}$) in 1930. The quantity is defined as the ratio of the observed
critical load $P_{c}^{e}$, to that predicted by the theory, i.e. $k_{d}\equiv
P_{c}^{e}/P_{c}$. Based on surveyed experimental results Lee2016 ; Wagner2020
, engineers at NASA settled for the extremely conservative design code of
$k_{d}\approx 0.2$ for spherical shell structures NASASP8032 .

Over the decades that followed, extensive work was dedicated to correcting the
persistent overprediction of the critical pressure. This involved studies of
the post-buckling behavior VonKarman1939 ; Tsien1942 , and the imperfection
sensitivity Koiter1945 ; Koiter1969 ; Hutchinson1967 ; Bushnell1967 ;
Krenzke1965 ; Koga1969  of thin spherical shells. Yet, marked success arrived
only recently, after Lee et al. developed a new fabrication technique for
polymeric spherical caps Lee2016a . The authors used this method to create
dominant dimple-like defects (larger than those naturally occurring in the
shell) with systematic size variations Lee2016 . These experiments, validated
with finite element modeling (FEM) and numerical analysis, quantitatively
showed for the first time how the imperfection depth lowers the critical load.

Other contributions followed, including studies on spheres with similar dimple
defects and sinusoidal equatorial undulations Hutchinson2016 , large-amplitude
dimples Jimenez2017 , through-thickness defects Yan2020  (a notable
predecessor is Ref. Paulose2013 ), dent defects Gerasimidis2020 , and probing
force imperfections Gerasimidis2018 , which collectively clarify the effect of
the type of defect on the knockdown factor for spherical shells. More broadly,
this long-awaited breakthrough provoked a new surge of progress in spherical
shell theory (see e.g. Refs. Hutchinson2016 ; Hutchinson2017 ; Audoly2020 ;
HutchinsonEMLoverview ). These developments afford engineers the opportunity
to design sturdy structures with more specific – and permissive Wagner2020  –
lower bounds on the load carrying capacity. This was the initial goal. In more
recent years, though, a community of researchers has adopted, in a sense, the
opposite goal: to design structures that buckle and snap on command Holmes2019
; Reis2015  for functions like colloidal self-assembly Sacanna2010 ,
encapsulation Shim2012 , inflatable snapping actuation Keplinger2012 ;
Gorissen2020 , and artificial muscle actuation Yang2016 .

The new understanding of defect sensitivity in a different light demonstrates
the tunability of spherical shell buckling, and thus serves this
“buckliphilic” Reis2015  community equally. In particular, the geometry and
placement of a dominant defect prescribe, respectively, the buckling strength
and the spot where an instability localizes. Recent extensions of this concept
couple geometric defects with differentially swelling Lee2019  or magneto-
responsive Yan2020b  materials to modify the knockdown factor over time.
Relatedly, a more general study showed how a homogenous natural curvature –
which can be a proxy for nonmechanical stimuli like thermal expansion, changes
in pH, or differential growth – acts to raise or lower the knockdown factor in
spherical shells Holmes2020 .

Besides control over when instability occurs – which depends on geometry,
loading, and mechanics – mechanical actuators generally rely on reversibility.
Repeatable actuation calls for robust, elastic materials, and silicone rubbers
like polydimethylsiloxane (PDMS) and vinyl polysiloxane (VPS) have answered
this call in mechanics research Reis2015b . In addition to their elastic
behavior, these elastomers are readily accessible and allow for fast, easy
fabrication Lee2016a .

Recently, Djellouli, et al. combined these ingredients to produce a mechanical
swimmer Djellouli2017 . Quasi-static pressure cycles drove the device, a
defect-seeded spherical shell made of the elastomer Dragon SkinTM 30, to
propel forward through a viscous fluid by buckling and unbuckling
controllably. The authors propose that maintaining dimensionless quantities
constant would allow for miniturization, with implications for drug delivery.
A natural extension of this work, and the motivation for the present study, is
to seek control over the speed of swimming by adjusting the frequency and/or
amplitude of pressure cycles. Largely because shell and fluid motion are
highly coupled in swimming, and because of possible resonance with postbucking
oscillations Mokbel2021 , we expect this phase space to be complex. Thus, we
set out to first isolate the shell buckling response, independent of fluid
motion, to dynamic loading at pressures in the vicinity of the critical load.

We fabricated imperfect spherical shells like those in Ref. Djellouli2017 ,
and fixed them in place surrounded by air. A small nozzle allows for internal
pressure control, through which we step-load the shells – that is, we abruptly
apply, and then maintain, a pressure load. Explicitly, we reduce the pressure
inside the shell cavity, creating a negative inside-outside differential
pressure. For simplicity, we will refer to this pressure difference in terms
of its magnitude. These straightforward experiments produced surprisingly rich
results. Even for loads below the experimentally measured elastic critical
pressure (which we loosely call “subcritical” herein), we consistently observe
buckling. Further, this buckling at subcritical loads occurs abruptly, often
after an extended period of very slow deformation perhaps mistakeable for
stability. We observe singular thresholds and a delay time which increases
monotonically as pressure decreases – this contrasts the findings of a recent
numerical study on dynamic step loading of spherical shells which are much
thinner than our own Sieber2019 .

As discussed, geometric imperfections can lead to buckling at lower-than-
expected loads. In this case, though, the shell defects are already accounted
for. Devoid of any plausible geometric explanation for this strange buckling
behavior, our results require a closer examination of the materials. Although
silicone elastomers are selected precisely because they behave elastically in
most settings, they are in fact prone to time-dependent molecular
rearrangement, and hence are viscoelastic Lin2009 ; Yuk2017 ; Urbach2020 .
Thus, they behave differently depending on how fast they are loaded, and
exhibit both stress-relaxation (softening when subjected to a constant strain)
and creep (deformation over time under constant stress).

In the present work, we rely on this viscoelasticity to account for our
observations, while emphasizing the close connections to elastic spherical
shell buckling. The structure is as follows: First, we introduce key aspects
of our experiments in Sect. II. Specifically, we report the viscoelastic
material properties, describe the geometry of our imperfect spherical shells,
and briefly introduce our experimental setup. In Sect. III, we present an
overview of our findings. We address the critical pressure conundrum in Sect.
IV by defining two pressure thresholds: the elastic critical pressure, and the
lower viscoelastic critical pressure, which can be related through the
limiting material properties. These thresholds separate three regimes:
immediate buckling, delayed buckling, and no buckling (stable). In Sect. V, we
introduce an analogy wherein viscoelastic creep deformation lowers the
critical load in the same way that a growing dimple-like defect would in an
elastic shell. This provides insight about the pre-buckling deformation
(Sect.V.1), and reveals how the delay time preceding buckling depends on the
imposed pressure, shell geometry and material properties (Sect.V.2). Finally,
we offer concluding remarks in Sect. VI.

##  II Materials, geometry and methods

We have performed dynamic, step-loaded pressure buckling experiments on soft,
viscoelastic spherical shells. Here, we briefly summarize the material
properties, shell fabrication and geometry, and dynamic loading methods.

###  II.1 Material characterization

The shells used in our experiments are made of the elastomer Dragon SkinTM 30
(manufactured by Smooth-On; Poisson’s ratio $\nu=0.5$ Djellouli2017 ). We
assume that our viscoelastic material can be described by the Standard Linear
Solid (SLS) model, the simplest linear model that captures both stress-
relaxation (the decreasing stress response over time for a structure subjected
to a constant strain) and creep (deformation under a prolonged constant
stress) Zener1948 . The SLS model describes limited creep behavior, i.e. creep
deformation does not progress indefinitely, nor does the modulus eventually go
to zero.

According to the SLS model, the modulus relaxes over time according to:

$E(t)\equiv\frac{\sigma(t)}{\varepsilon}=E_{\infty}+E_{1}e^{-t/\tau_{\sigma}}$
(2)

where $\sigma(t)$ is the time-varying stress, $\varepsilon$ the constant
strain, and $\tau_{\sigma}$ the relaxation time. The parameter $E_{1}\equiv
E_{0}-E_{\infty}$ quantifies the total stiffness lost as the elastic modulus
$E_{0}$ decreases to the long-term (equilibrium) modulus $E_{\infty}$, where
$E_{0}\geq E_{\infty}$. The function (2) is known as the relaxation modulus.
It is related by Laplace transform to the creep compliance function, which
describes the temporally increasing strain $\varepsilon(t)$ of the SLS element
under imposed constant stress $\sigma$ Lakes2009 :

$J(t)\equiv\frac{\varepsilon(t)}{\sigma}=J_{0}+J_{1}\big{(}1-e^{-t/\tau_{\varepsilon}}\big{)}$
(3)

where $J_{0}=E_{0}^{-1}$, $J_{1}\equiv J_{\infty}-J_{0}$ with
$J_{\infty}=E_{\infty}^{-1}\geq J_{0}$, and $\tau_{\varepsilon}$ is the
retardation time.

We performed uniaxial tension and stress-relaxation tests using the tensile
testing machine Instron 5943 to identify the parameters in Eq. (2) (details
are provided in Sect. VII.1.) We found the relaxation time to be
$\tau_{\sigma}=0.78\pm 0.49$ s. As for the moduli, we determined
$E_{0}=0.59\pm 0.04$ MPa and $E_{\infty}=0.54\pm 0.08$ MPa. Reported errors
throughout the text correspond to one standard deviation unless otherwise
noted111For quantities with few measurements, the approximate standard
deviation is reported as one-fourth of the error range.. The resulting ratio
of the mean long-term modulus to the mean instantaneous one is $\bar{E}\equiv
E_{\infty}/E_{0}=0.91$, which is central to the analysis beginning in Sect.
IV.

Since our materials have a relatively low relaxation strength, defined as
$\Delta=E_{1}/E_{\infty}$ Lakes2009 , Eq. (2) and the inverse of Eq. (3)
differ negligibly, i.e. $J^{-1}(t)\approx E(t)$ (see Sect. VII.1, Fig. 7).
While creep is the relevant process in our experiments, our primary aim is to
draw connections to elastic shell theory, which relies on the modulus $E$.
Thus, we will use this convenient fact to interchange the representation of
these two mechanisms in our analysis.

###  II.2 Shell geometry

Spherical shells were fabricated following the bi-molding method from Ref.
Djellouli2017 , which is detailed in Sect. VII.2. The shells, made of two
hemispheres seamed with a thin layer of diluted polymer, all have an outer
radius of $R_{o}=R+h/2=25$ mm. The thickness $h\in[1$ mm, $5$ mm$]$, such that
$\eta\in\\{24.5,12.0,6.6,4.5\\}$.

Figure 1: Schematic of the clamped shell with an inlet for pressure control.
The flexible tube connects at its other end to a vacuum tank so that when the
line is open, $P_{\text{in}}<P_{\text{atm}}$ (or
$P\equiv|P_{\text{in}}-P_{\text{atm}}|>0$). The relevant parameters are
labeled: the nominal thickness $h$ and midline radius $R$, imperfection depth
$\delta$, and half-angular width of the imperfection $\beta_{0}$.

To control the location and direction of buckling, each shell is seeded with
an imperfection, where the thickness is reduced by an amount
$\delta\in[0.40,0.81]$ mm, such that
$\bar{\delta}\equiv\delta/h\in[0.08,0.76]$ in a circular region spanning a
half-angle of $\beta_{0}\approx\pi/24$ radians. The shell parameters are shown
schematically in Fig. 1, and the effect of this defect is discussed further in
Sect. VII.3.

Following e.g. Ref. Lee2016 , we also introduce the parameter $\lambda$,
defined as

$\lambda=\big{(}12(1-\nu^{2})\big{)}^{1/4}\eta^{1/2}\beta_{0},$ (4)

which describes the defect geometry in the context of spherical caps
Kaplan1954 ; Taffetani2018 . For our shells in order of increasing thickness,
$\lambda\in\\{1.038,0.716,0.612,0.492\\}$.

###  II.3 Step loading

A flexible tube was connected on one end to the inside of the shell, and on
the other to a vacuum tank. An electrovalve interrupting this channel allowed
us to abruptly remove air from the inner volume of the shell, creating a
pressure difference of magnitude $P\in[0.3$ kPa, $46.5$ kPa$]$ (see Fig. 1),
which we monitor with a pressure sensor. The pressure load is maintained for
either the time it takes the shell to buckle, or $t_{hold}\in[5s,360s]$. For
details on the experimental setup, see Sect. VII.4.

##  III Three regimes

For each of our shells, the immediate response to any non-negligible pressure
load was qualitatively the same: The shell compresses as soon as the pressure
is felt, and deformation quickly localizes at the unclamped pole (in the
vicinity of the imperfection), forming a dimple-like depression with a
deflection depth $w$ (see Fig. 2 a.) Beyond this early behavior, which occurs
in approximately the first 0.05 seconds, three regimes were evident from our
experiments.

Figure 2: (a.) Selected high-speed camera images and (b.) corresponding plots
of the pole deflection normalized by the shell thickness, i.e. $\bar{w}=w/h$,
show the typical response of a shell ($\eta=6.6$) to step loading at
relatively i. high (supercritical), ii. moderate (subcritical), and iii. low
pressures. The pole deflection $w$ is measured from the initial state, marked
in a,i. by the upper dashed gray line. At or above the elastic buckling
pressure (e.g. $\bar{P}\equiv P/P_{c}^{e}=1.02$), the shell quickly buckles
($t_{c}\approx 0.07$ s). The corresponding critical pole deflection $w_{c}$ is
indicated in a,i. Even at subcritical pressures, e.g. $\bar{P}=0.93$, the
shell eventually buckles ($t_{c}\approx 1.26$ s). This collapse follows a
deceleration in the pole deformation at $t\approx 0.05$ s, and a subsequent
period of slow, constant-rate deformation. For pressures below a second
threshold, e.g. $\bar{P}=0.90$, the shell does not buckle.

If the imposed pressure is high, the initial fast rate of pole deformation is
maintained, and the shell quickly buckles – that is, the pole inverts, driving
global collapse (Fig. 2, i). In postbuckling, which we do not study in detail
here, the pole region is completely inverted and oscillations occur before
stability is reached again when $w\approx 2R$. We refer to this regime,
wherein the shell behaves elastically throughout deformation, as the immediate
buckling regime. The lowest pressure at which we observe this buckling
behavior defines the experimental elastic critical load $P_{c}^{e}$. Due to
the seeded defects and the non-negligible thickness of our shells, the
experimentally-determined value $P_{c}^{e}$ differs from the theoretical
critical pressure $P_{c}$ (Eq. (1)) for thin, perfect elastic shells. For
details, see Sect. VII.3.

At slightly lower pressures, the deformation rate slows considerably following
the fast response to loading, at a transition time ($t\approx 0.05$ s). The
dimple slowly deepens (i.e. $w$ increases; See Fig. 2, ii), before an abrupt
acceleration after some time $t_{c}\in[0.09$ s, $17.09$ s$]$ signifies
buckling. We define this intermediate regime as the delayed buckling regime.
At still-lower pressures, slow pole motion eventually stops, and the shell
settles into indefinite stability for as long as the load is maintained. We
call this third regime the stable regime (Fig. 2, iii).

The value of the pressure which separates the delayed buckling and stable
regimes, and hence marks the boundary of whether collapse will occur, is
clearly of interest. This lower pressure threshold was more or less constant
for all of our shells when normalized by the elastic load. In other words, the
reduced critical pressure is independent of geometry. With this nudge toward
the materials, we proceed to rationalize these findings.

##  IV Pressure thresholds via modulus ratio

Since we know our materials are viscoelastic, we can presume that the slow
pole deformation under constant pressure is an exhibition of creep. This would
situate our observations in the terrain of creep buckling. Creep buckling was
introduced in the literature in 1951 Rosenthal1951 . The bulk of the work in
this field was developed in the thirty or so years that followed, and was
aimed at understanding creep collapse that occured on timescales of hours or
even days in metallic, mono-resin materials, and reinforced concrete. However,
general theories emerged, which are illuminating when applied to our elastomer
shells. In particular, Hayman Hayman1978 ; Hayman1981  and others Hoff1973 ;
Minahen1993  proposed that a viscoelastic structure that buckles due to creep
may be treated as an equivalent elastic structure with a lower critical load.
The main result is that the lower threshold – which we will henceforth refer
to as the viscoelastic critical pressure $P_{c}^{v}$ – is directly related to
the long-term modulus. For spherical shells, the lower critical pressure
$P_{c}^{v}$ may be found by simply replacing Young’s elastic modulus $E$ with
$E_{\infty}$ in calculating the elastic critical load (Eq. (1)). Denoting
normalization by the experimentally measured elastic critical pressure
$P_{c}^{e}$ with an overbar, i.e. $\bar{P}\equiv P/P_{c}^{e}$ and
$\bar{P}_{c}^{v}\equiv P_{c}^{v}/P_{c}^{e}$, that is:

Figure 3:  Phase plot depicting the three regimes, separated by two pressure
thresholds. The elastic critical pressure (thin dotted line, $\bar{P}=1$) was
experimentally determined as the minimum pressure at which a slowdown of
dynamics does not preceed buckling. We expect “immediate” (elastic) buckling
for an imposed pressure $\bar{P}\geq 1$, which corresponds to the green
region. The theoretical viscoelastic critical pressure
$\bar{P}_{c}^{v}=\bar{E}$ (thick dashed line) is determined from (5) using the
material parameters measured independently in stress-relaxation tests (see
Sect. II.1). Between $\bar{P}=1$ and $\bar{P}_{c}^{v}$ (yellow region),
viscoelastic creep can explain delayed instability at subcritical pressures.
Below $\bar{P}=\bar{E}$ (red region), limited creep for our material is
insufficient to cause buckling, so we expect indefinite stability. Green
triangles, yellow circles, and blue squares represent, respectively,
experiments which buckled elastically, buckled after a time delay, and did not
buckle. Error bars represent one standard deviation.

$\bar{P}_{c}^{v}=\bar{E}.$ (5)

In Fig. 3, we show that (5) agrees with our experimental data very well,
solidifying the notion that viscoelasticity is indeed the cause for the
subcritical buckling we observe. The dashed line marks the theoretical lower
limit where creep bucking may be observed, $\bar{P}_{c}^{v}=\bar{E}$, which
for our materials (see Sect.II.1 & Sect. VII.1) is $0.91$.

Explicitly, Eq. (5) ignores the actual mechanism that leads to instability,
creep deformation, in favor of a straightforward way to determine the minimum
buckling pressure for creep-limited materials, which is perhaps the most
crucial information for any design goal. Given the elastic (viscoelastic)
critical pressure and the experimentally-determined instantaneous and long-
term moduli, the ratio in Eq. (5) readily predicts the viscoelastic (elastic)
critical pressure.

This view is effective and completely general with respect to geometry and
material properties. However, it provides no information about deformation or
the time delay that preceeds buckling. Besides that these features are of
fundamental interest, understanding this delay mechanism will offer a route to
including controllable delays in elastomer device design. We address these
open questions in the following section.

##  V Creep deformation as an evolving defect

As we have seen, the efficient, modulus-based approach in Sect. IV connects
the limiting critical pressure of a viscoelastic shell to that of the
equivalent elastic shell. It leaves questions, however, about the time it
takes a subcritically-loaded shell to buckle, and the underlying pre-buckling
deformation. Traditional analytical approaches to capture the critical time
and/or deflection for creep buckling involve incorporating calculated
quantities for stress and strain into the constitutive model (in our case Eq.
3). Instability may be identified by solving the eigenvalue problem of the
governing differential equations, or by the quasi-static “critical strain
approach” Gerard1958  wherein the critical strain must be known or assumed a
priori, and the corresponding time is directly solved for Miyazaki2015 . These
methods require precise representations of the stresses and strains throughout
deformation, and have met moderate success in capturing experimental behavior
for simple structures like columns Hoff1954 ; Hoff1956 , trusses and arches
Huang1967 , plates and even cylinders Gerard1958 . (See Ref. Miyazaki2015  for
a review of the relatively recent work on creep buckling of shell structures,
or Ref. Hoff1973  for an earlier review on creep buckling of plates and
shells.)

A clear problem with these approaches is that it is generally assumed that a
shell undergoing creep will lose stability at the same strain as its elastic
counterpart Gerard1956 ; Hoff1975 . However, it has been noted that this
assumption often leads to underprediction of the critical displacement and
time Obrecht1977 . Indeed, although the immediate and delayed buckling regimes
appear qualitatively very similar in terms of deformation in our experiments
(see Fig. 2, i. & ii.), we observe that shells which creep for longer sustain
more deformation before buckling.

Complex geometries like imperfect spherical shells and their nonlinear
deformations introduce significant analytical difficulties of their own. Creep
buckling in spherical caps and complete spherical shells has primarily been
studied with numerical analyses Huang1965 ; Shi1970 ; Jones1976 ; Miyazaki1977
; Xirouchakis1980 , which do not produce closed-form solutions for when
instability occurs. Few experiments exist for comparison to these results, and
attempts to replicate limited experimental creep buckling behavior for
spherical shells have largely been unsuccessful Shi1970 ; Leckie1976 .

A more enlightening approach relies on the observation that the pre-buckling
deformation approximately amplifies the initial defect (see Fig. 4a). This
suggests that we may be able to draw an analogy between creep deformation and
a growing imperfection in an elastic shell. This concept was proposed by
Hayman in 1981, but the author conjectured that while his “locus of critical
points” approach offered intuition, it lacked predictive power for all but
simple, statically determinate structures Hayman1981 . To our knowledge, the
approach has not been implemented besides in the original work, when it was
validated against a small number of experiments on concrete three-pin arches.

One key challenge was that the effect of the defect must be quantifiable.
Indeed, as discussed in Sect. I, the imperfection sensitivity of shell
structures has only recently become experimentally tractable. In what follows,
we rely on the work of Lee and collaborators Lee2016 , and center our analysis
upon quantities that are readily determined through experiments, to derive a
practical and predictive model for creep buckling in our shells.

Knowing that an imperfection “knocks down” the critical load in an elastic
structure, in this view quasi-static creep deformation has the same effect. In
other words, creep deformation behaves like an evolving imperfection by
progressively knocking down the critical load. It follows that creep collapse
occurs when the critical load associated with the “imperfect” creep-deformed
structure falls to the value of the applied load. Conversely, the imperfection
size associated with a given applied load should correspond to the critical
deformation at which creep buckling occurs. It is important to note that this
was buried but implied in the approach of Sect. IV; the two arguments are
complementary.

To make this analogy quantitative, we turn to the recent literature on
geometric imperfections in elastic spherical shells. The true defects in our
shells are characterized by a local reduction in thickness as in Ref. Yan2020
(see Sect. VII.3). The pole deflection generated during creep, however, is
qualitatively more similar to a dimple-like imperfection where the thickness
does not change, but the curvature of the shell midline does. Because the
midline curvature of our shells is unaffected by the thickness reduction
except at the discontinuity at either edge of the defect profile, we consider
our shells initially “perfect” in the dimple sense. For the present analysis,
we rely on the findings of Lee et al. Lee2016 , which are in agreement with
Refs. Hutchinson1967 ; Hutchinson2016 ; Jimenez2017 .

Figure 4: Pre-buckling pole deformation ($w\leq w_{c}$) is qualitatively
similar to a geometric defect of increasing depth $\delta$. (a.) i. Schematic
of dimple-like defect in an elastic spherical shell, after e.g. Ref. Lee2016 .
ii. Edge contours from high-speed images at $t=0,0.07,0.14,0.21,0.28$ s and
$t_{c}=0.35$ s (darkening blue corresponds to increasing time) for $\eta=4.5$
and $\bar{P}\approx 0.97$. After initial compression, deformation localizes to
the pole and progresses. (b.) Our fitting parameter $a\approx
f(\lambda)=1.08-0.44\lambda$.

The key finding of their work, for our purposes, is that the knockdown factor
for a given shell $k_{d}=P_{c}^{e}/P_{c}$ is a function of
$\bar{\delta}=\delta/h$ which initially decreases for increasing
$\bar{\delta}$, then reaches a plateau. The authors present an empirically-
determined function describing the lower bounding envelope over the range of
$\lambda$ they study, which takes the form $k_{d}=a+b/(c+\bar{\delta}$). We
find that this functional form describes their individual curves sufficiently
well.

Primarily because $\lambda$ for our thick shells is below the range studied in
Ref. Lee2016 , we cannot directly extract the results relevant to our work.
Instead, we expect that these general trends will hold. We assume that we can
simply replace the imperfection depth $\delta$ with the pole deflection $w$.
Because we only consider the additional knockdown due to creep deformation,
and not that due to the inital defect, we take the reference pressure as the
experimental critical pressure $P_{c}^{e}$ that corresponds to the initially
imperfect shell. Then we define a general viscoelastic knockdown function

$k_{d}^{v}(w)=\bar{P}_{c}(\bar{w})=a+\frac{b}{c+\bar{w}},$ (6)

where $\bar{w}=w/h$, and $a$, $b$, and $c$ are yet unknown.

Since $w$ increases according to the creep strain rate, an alternative form of
Eq. (6) specifies the time dependence. Approximating the magnitude of the
circumferential strain to first order as $\varepsilon\approx w/R$, we can say
$\bar{w}\approx\eta\varepsilon(t)$. From Eq. (3) this means
$\bar{w}\approx\eta\sigma J(t)$, which is approximately $\eta\sigma/E(t)$
since the creep compliance function $J(t)$ and the inverse of the relaxation
modulus $1/E(t)$ are nearly indistinguishable for our material. At early times
the shell behaves elastically, so we assume Hooke’s Law applies, i.e.
$\sigma\approx E_{0}\epsilon(t_{t})$ at a transition time $t_{t}$ between the
elastic and creep stages of deformation. Further, we assume
$\epsilon(t_{t})\approx\epsilon_{c}^{e}$, where $\epsilon_{c}^{e}\approx
w_{c}^{e}/R$ is the critical strain corresponding to $\bar{P}=1$, when the
shell buckles immediately following elastic deformation. Then by Eq. (2), the
normalized pole deflection increases with time following the relation

$\bar{w}(t)\approx\frac{E_{0}\>\bar{w}_{c}^{e}}{E_{\infty}+E_{1}e^{-t/\tau_{\sigma}}}.$
(7)

Note that due to our simplified representations of stresses and strains in Eq.
(7), all relations that follow are approximations, despite that we present
them as equalities for simplicity. Substituting Eq. (7) in Eq. (6), gives a
time-dependent version of the generalized viscoelastic knockdown function:

$k_{d}^{v}(t)=\bar{P}_{c}(t)=a+\frac{b}{c+\frac{E_{0}\>\bar{w}_{c}^{e}}{E_{\infty}+E_{1}e^{-t/\tau_{\sigma}}}}.$
(8)

It remains to determine the three unknown quantities, which should explain how
sensitive the critical pressure is to deformation (Eq. (6)) and how quickly
the critical pressure decreases (Eq. (8)) for each shell. To do so, we
constrain the functions, enforcing what we know about the limiting behavior.
The critical deflection required for buckling at the elastic limit where
$\bar{P}=1$ is an experimentally-determined value, $w_{c}^{e}$, which differs
for each shell based on geometry. Since for buckling to occur, $P_{c}(w)=P$,
from Eq. (6) this condition is stated as:

$1=a+\frac{b}{c+\bar{w}_{c}^{e}}.$ (9)

We also know from Eq. (5) that buckling will not occur below
$k_{d}=\bar{P}_{c}^{v}=\bar{E}$. Taking the limit of the right hand side of
Eq. (8) as $t$ approaches infinity gives a second constraint:

$\bar{E}=a+\frac{b}{c+\frac{\bar{w}_{c}^{e}}{\bar{E}}}.$ (10)

Solving Eqs. (9) & (10) simultaneously gives

$b=\frac{(a-1)(a-\bar{E})\bar{w}_{c}^{e}}{\bar{E}}$ (11)

and

$c=\frac{-a\bar{w}_{c}^{e}}{\bar{E}},$ (12)

which we insert into Eq. (6) to arrive at:

$\bar{P}_{c}(w)=a+\frac{(1-a)(\bar{E}-a)\bar{w}_{c}^{e}}{\bar{E}\bar{w}-a\bar{w}_{c}^{e}}$
(13)

which describes the critical pressure for a given degree of pole deflection.
If $\bar{w}$ is large enough that $\bar{P}_{c}$ is lowered to the imposed
dimensionless pressure $\bar{P}$, in theory buckling will occur.

Similarly, inserting the expressions for $b$ and $c$ in Eq. (8) gives

$\bar{P}_{c}(t)=a+(a-1)(a-\bar{E})\bigg{(}\frac{E_{0}\bar{E}}{E_{\infty}+E_{1}e^{-\frac{t}{\tau_{\sigma}}}}-a\bigg{)}^{-1}$
(14)

which specifies the time-dependence, according to the SLS model, of the
knockdown to the critical pressure that occurs as the pole deflection
progresses. Again, if $t$ is such that $\bar{P}_{c}=\bar{P}$, we expect
collapse to occur.

Note that from Eq. (14), it is clear that any explicit dependence on geometry
is contained in $a$. When $a$ is left as a fitting parameter for the curves
defined by Eq. (16) (see Fig. 6a), we find that it is a decreasing function of
the geometric defect parameter $\lambda$, which is in general agreement222We
fit a function of the form $k_{d}=a+b/(c+\bar{\delta})$ to the numerical data
over $\lambda\in[1,5]$ in Fig. 6b of Ref. Lee2016 , and found $a\approx
0.5-0.14\lambda$. with Ref. Lee2016 . In particular, we determine $a\approx
1.08-0.44\lambda$ (Fig. 4b) for our shells. We refer to this linear function
henceforth as $f(\lambda)$, and substitute it accordingly for $a$ to emphasize
the deduced connection between our fitting parameter and the defect geometry.

The deflection (Eq. (13)) and time (Eq. (14)) forms of the viscoelastic
knockdown relations lead to expectations for how the critical displacement and
critical time should depend on the material properties, geometry and the
imposed pressure. We assess both in the following subsections.

###  V.1 Critical deflection

Figure 5: The critical pole deflection $w_{c}$ depends on $P$, the material
properties, and the geometry through $a=f(\lambda)$ and $P_{c}^{e}$: when
creep occurs (delayed buckling, yellow region, to the left of the dotted line
at $1$ on the horizontal axis which marks $\bar{P}=1$), the shell must deform
more than $w_{c}^{e}=w_{c}(P_{c}^{e})$ to “knock down” the critical pressure
until it coincides with the imposed pressure $P$, and stability is lost. This
is represented by the diagonal line (Eq. (15)), which matches our data well
until the lower pressure limit for buckling is approached. This asymptote is
marked by the loosely dotted line with abscissa $1/\bar{E}$, which corresponds
to $\bar{P}=\bar{E}=P_{c}^{v}/P_{c}^{e}$ (Eq. (5)). Open symbols in the red
(stable) region represent the maximum (equilibrium) $w/w_{c}^{e}$ for
experiments that did not buckle. As expected, these points fall below the
theoretical line, which projects the deflection necessary for buckling. Our
model does not extend to the immediate buckling regime (green region). Error
bars are approximately one standard deviation.

Evaluating Eq. (13) at $P=P_{c}(\bar{w}_{c})$, replacing $a$ with
$f(\lambda)$, and solving for $\bar{w}_{c}/\bar{w}_{c}^{e}$ gives the
following expression for the dimensionless critical deflection:

$\frac{w_{c}}{w_{c}^{e}}=\frac{f(\lambda)(1+\bar{E}-\bar{P})-\bar{E}}{\bar{E}(f(\lambda)-\bar{P})}$
(15)

which is valid when creep occurs ($\bar{P}_{c}^{v}=\bar{E}<\bar{P}<1$, where
the first equality refers to Eq. (5)) and tells us the pole deflection we
should expect at (or, that is required for) buckling. Eq. (15) is plotted
against our data in Fig. 5, showing good agreement until
$\bar{P}\approx\bar{P}_{c}^{v}$. Near this asymptote, marked by the loosely
dotted line, the critical deflection observed in experiments exceeds the
predicted value. We discuss this deviation further in Sect. V.2.

The results in Fig. 5 are consistent with the trend observed in the creep
buckling literature, where it is suggested that a shell undergoing creep prior
to buckling sustains more deformation before collapse than its counterpart
that buckles elastically. The explanation is that as the imposed pressure
decreases, more deformation is required to sufficiently “knock down” the
critical pressure for buckling to occur.

Creep does not occur before buckling when $\bar{P}\geq 1$ (the green region in
Fig. 5), where the shell behaves elastically. In this region, where Eq. (15)
does not apply, we do not observe a clear trend in the critical deflection.
Further, for experimental points that do not buckle, the deformation was
insufficient to reduce the critical pressure to the value of the relatively
low applied load. Accordingly, the maximum pole deflection (open symbols in
Fig. 5) falls below theoretical curves when $\bar{P}<\bar{P}_{c}^{v}$ (red
region in Fig. 5). Taking the limit as $t$ approaches infinity in either Eq.
(3) or Eq. (2) gives $\varepsilon_{\text{max}}=\sigma/E_{\infty}$, so
$w_{\text{max}}\approx R\sigma/E_{\infty}$. Meanwhile, Hooke’s law provides an
estimate for the critical deflection if $\bar{P}=1$ is imposed:
$w_{c}^{e}=R\sigma_{e}/E_{0}$. Then in theory, if $\bar{P}<\bar{P}_{c}^{v}$,
$w_{\text{max}}/w_{c}^{e}\approx\sigma/(\bar{E}\sigma_{e})\sim\bar{P}/\bar{E}$.
We do not attempt to verify this scaling with our sparse data in this region,
besides noting that as $\bar{P}$ approaches $\bar{P}_{c}^{v}$, we would expect
$w_{\text{max}}/w_{c}^{e}$ to approach 1 (from below). Our data appears to
support this conjecture.

###  V.2 Critical time

We have seen that the pole deflection is analogous to a dimple-like defect.
Incorporating the viscoelastic material model, in turn, tells us how the
critical pressure is expected to decrease over time. This offers a means to
explain the critical buckling time, which increases monotonically for
decreasing pressure. To do so, we evaluate Eq. (14) at
$\bar{P}=\bar{P}_{c}(t_{c})$ and solve for the dimensionless critical time
$t_{c}/\tau_{\sigma}$:

$\frac{t_{c}}{\tau_{\sigma}}=\ln\bigg{(}\frac{(\bar{E}-1)(f(\lambda)(\bar{E}-\bar{P}+1)-\bar{E})}{\bar{E}(f(\lambda)-1)(\bar{E}-\bar{P})}\bigg{)}.$
(16)

Eq. (16) is plotted against our data in Fig. 6. Like for the critical
deflection (Fig. 5), the knockdown theory captures the critical time well in
the intermediate range. At or above the elastic limit $\bar{P}=1$, Eq. (16)
predicts nonphysical critical times of $t_{c}\leq 0$. This is because the SLS
model assumes that both loading and initial elastic deformation happens
instantaneously, so creep begins at $t=0$ s. Of course, in reality elastic
instability occurs on a timescale associated with the elastic wavespeed.
Accordingly, the inertial timescale $t_{*}$ begins to dominate the
viscoelastic one as $\bar{P}$ approaches $1$. Following e.g. Refs. Pandey2014
; Gomez2016 ; Sieber2019 , we expect that the elastic timescale
$t_{*}\sim(2R)^{2}/(ch)$, where $c=\sqrt{E_{0}/\rho}=23.37$ m/s is the speed
of sound within the material and $\rho=1080$ kg/m3 is the material density
according to the manufacturer. For our shells in order of increasing
thickness, this gives $t_{*}\sim\\{0.1027,0.0493,0.0264,0.0173\\}$ s. We have
indicated the elastic snap-through time for an arch
$\tau_{inertial}=2\sqrt{3}t_{*}$ Pandey2014 , with horizontal dashed lines in
Fig. 6a.

Figure 6: As the dimensionless applied pressure increases from
$\bar{P}_{c}^{v}=0.91$ to the elastic critical pressure $\bar{P}=1$, the
buckling time decreases monotonically according to Eq. (16), which specifies
how long it takes for the pole to deform enough to decrease the critical
pressure to the applied one. (a.) Experimental data (markers) and color-
corresponding curves from Eq. (16), which were fit to obtain
$f(\lambda)=0.62,0.76,0.81,0.85$ (plotted in Fig. 4b). for
$\eta=4.5,6.6,12.0,24.5$, respectively. The densely dotted vertical line marks
$\bar{P}=1$. Near this limit, the inertial (elastic) timescale, indicated by
horizontal dashed lines, sets $t_{c}$ for each shell. The loosely dotted
vertical line marks the asymptotic pressure set by
$\bar{E}=0.91=\bar{P}_{c}^{v}$ by Eq. (5). (b.) Eq. (16) collapses the same
experimental data plotted in (a.), and captures the critical time for
intermediate pressure values. When $\bar{P}$ decreases to about $0.94$, the
theory underpredicts the critical time. Error bars correspond to one standard
deviaton in $\bar{P}$. Horizontal dashed lines in a. & b. indicate the elastic
snap-through time $\tau_{inertial}$, which sets the minimum $t_{c}$.

When the imposed pressure nears the lower limit $\bar{P}_{c}^{v}$, the
asymptotic behavior is captured qualitatively (see Fig. 6a). However, the
theory underestimates the critical time in this region. This divergence occurs
around $\bar{P}\approx 0.94$, which corresponds to when the predicted buckling
time surpasses $\tau_{\sigma}$. Thus, deformation has slowed considerably
prior to buckling for these experiments and the inertia that was present at
early times is no longer available. We have identified elastic buckling in our
experiments as when the deformation rate exceeds that at very early times
(when the shell also behaves elastically). However, it is possible that
buckling initiated sooner in reality, but that the shell needs further
perturbation – that is, to deform more, which requires more time – before we
detect collapse. This is reminiscient of critical slowing down phenomena,
wherein dynamics slow considerably near instability Gomez2016 ; Gomez2018 .
Other possible explanations for the deviation from our model are the
simplified representations of stresses and strains, or the inability of the
SLS model to capture the material behavior exactly. Nonetheless, we conclude
that despite the notable simplicity of our assumptions, the knockdown theory
explains our observations quite well, and does so while emphasizing the close
connections between elastic and viscoelastic shell buckling.

##  VI Conclusion

In summary, we subjected thick, spherical, defect-seeded viscoelastic shells
to step pressure loading. We observed three regimes: When the pressure load
was at or above the experimentally determined elastic critical load, we
observed predictible elastic behavior, i.e. prompt buckling. At intermediate
loads – below the elastic critical pressure – the shell buckled, albeit after
a time delay during which deformation slowly progressed. At still-lower
pressures, the shell deformed but collapse never occured. Our aim in this work
was to rationalize our findings in a way that maintains close ties to elastic
shell buckling, and is readily useable for experiment or design goals. To this
end, we demonstrated that the load thresholds, critical deflection, and
critical time may all be captured by a framework that treats creep deformation
like an evolving defect in an elastic shell.

In particular, the ratio of the long-term modulus to the short-term (elastic)
one is the same as as the ratio of the two critical pressures. This result is
rooted in elastic shell theory, but practically, the material properties alone
can explain the two pressure thresholds. This finding was suggested in various
theoretical works on creep buckling Hayman1978 ; Hayman1981 ; Hoff1973  and is
independent of geometry, and hence is completely general.

We used this fact and existing work Lee2016  on defects in spherical shells to
discern an expression for how deformation due to creep acts to “knock down”
the critical pressure. In this view, the shell loses stability when creep
deformation, which localizes at pole and amplifies the initial imperfection,
progresses enough to reduce the critical pressure to the value of the imposed
one. This allowed us to capture the dependence of the critical deflection on
the imposed pressure (normalized by the experimentally-determined elastic
critical pressure), the modulus ratio, and the defect geometry (Eq. (15)).
This offers an explanation rooted in elastic shell behavior for a decidedly
viscoelastic phenomenon: in the delayed buckling regime, higher deflection is
required for instability as the pressure decreases.

Because deformation occurs on a timescale sufficiently well-described by our
chosen viscoelastic material model (SLS), a time-dependent form of the
viscoelastic knockdown function immediately follows. From this we devise an
expression for how the pre-buckling delay time depends on the same quantities:
the modulus ratio, the dimensionless pressure, and the defect geometry (Eq.
(16)). The buckling time increases monotonically but non-linearly as the
pressure decreases, which is generally captured by our model. While the
modulus ratio $\bar{E}$ was fixed in our experiments, we expect that our model
is valid for any material with relatively low relaxation strength (i.e. the
relaxation modulus and creep compliance functions do not differ
significantly). Further, we note that the success of our model does not depend
on the approximation $J^{-1}(t)\approx E(t)$ that we have employed. Rather, if
one were to obtain the creep compliance function in experiments, using this
directly in place of the relaxation modulus would likely improve the accuracy
of predictions.

While viscoelastic shells behave elastically during buckling, viscous effects
re-enter in later stages of postbuckling, as discussed in Ref. Coupier2019 .
Unbuckling was studied in detail in Ref. Djellouli2017 , as the non-reciprocal
nature of the buckling-unbuckling cycle is the source of motility. We did not
examine unbuckling in the present work. However, we could reasonably expect
that another viscoelastic delay phenomenon, termed “pseudo-bistability”
Brinkmeyer2012  (first introduced in Ref. Santer2010  as “temporary
bistability”), could be seen in our shells.

Pseudo-bistability refers to the delayed “snap-back” instability that occurs
in the unloaded state, following a loading-unloading sequence that induces
both stress-relaxation and creep Gomez2018 . Early works modeled
viscoelasticity via an evolving stiffness, acheiving qualitative agreement
with experiments Santer2010 ; Brinkmeyer2012 . A recently-proposed metric
framework introduces viscoelasticity as a temporally evolving (fictitious)
reference length instead, and among its merits is the ability to predict
delayed snap-back instability Urbach2020 ; Urbach2018 . Delayed snap-back
instability occurs in the unloaded state, and in our setting would require
maintaining the pressure load for a sufficiently long time while the shell is
fully buckled before unloading. As such, creep instability and pseudo-bistable
snap-back can in theory be induced in the same thick structure, with one or
the other surpressed at will.

Our findings highlight the importance of the load rate, and the sensitivity of
the shell to pressure variations in the vicinity of the elastic critical load
(when approached from below). These are important considerations for
structural designs using silicone rubber where either stability or elastic
behavior is desired. For instance, knowledge of the critical load thresholds
is clearly important when designing an efficient spherical swimmer, or indeed
any viscoelastic structure that should undergo oscillatory instability. If the
goal is fast motility, it would likely be desirable to minimize the time delay
before buckling by avoiding the delayed buckling (intermediate pressure)
regime altogether.

In other settings, viscoelastic behavior can enhance the functionality of
reversibly actuatable structures. This has been demonstrated recently in
designs that rely on pseudo-bistable snap-back, e.g. 3D-printed viscoelastic
metastructures whose time-dependent properties are tunable based on
temperature Che2019 , with even more flexibility afforded by using multiple
viscoelastic materials Che2018 . Another study examines the interplay between
viscous dissipation and geometric hysteresis, as a function of the strain
rate, for the design of optimal energy dissipating metamaterials Dykstra2019 .
Introducing tunable delays via creep buckling has not yet been explored, but
the potential is vast: A switch or capsule could be loaded subcritically, so
that the shell has time to move to a desired location before buckling occurs.
Over time as deformation reduces the critical load and thus the energy barrier
Hutchinson2018 , a much smaller probing force or other perturbation could
trigger buckling. This could be useful for pneumatic gripping. A mechanical
signal of fixed input frequency (and varying amplitude) could produce a varied
output frequency, which has implications for mechanical computing. Because our
analysis relies only on quantities that are straightforward to determine in
experiments, our findings are especially amenable to accessing such
tunability. As we have shown, these possibilities are achievable with the
nearly-elastic materials that are already common and cherished in mechanics
research.

In 1956, Gerard wrote about creep buckling that “there are almost more
theories than reliable test points which can be used to check the theories”
Gerard1956 . While a limited number of experimental contributions have come
about since, to our knowledge our experiments on full spheres are novel to the
existing creep buckling literature. Further, we have demonstrated that some
concepts central to general creep buckling theories, which previously were
mostly tested on metallic, mono-resin materials, and reinforced concrete
Miyazaki2015 , are indeed applicable to soft, rubbery materials. We expect
that the concepts we have studied can be extended to other geometries, other
loading methods, and similar polymers.

## Conflicts of interest

We declare that we have no competing interests.

## Funding

D.P.H. and L.S.M. are grateful for financial support from the NSF
CMMI-1824882. G.C. acknowledges financial support from the European
Community’s Seventh Framework Programme (FP7/2007-2013) ERC Grant Agreement
Bubbleboost no. 614655.

## Acknowledgements

We thank Catherine Quilliet and Philippe Marmottant for many helpful
discussions, Adel Djellouli for the shell fabrication protocols and equipment,
Patrice Ballet for his technical support, and Francisco Lopez-Jimenez and
Pedro Reis for providing data from their work in Ref. Lee2016

##  VII Appendix

###  VII.1 SLS parameters

To determine the material properties of Dragon Skin 30TM, we fabricated 7
dogbone specimens (ASTM D412 Type C). The molds, cut using the Epilog Helix
laser machine, were designed to give each sample a protruding defect – two
thin horizontal lines 500 microns thick and separated by 500 microns – at the
middle of the gauge section. These lines were tracked with a zoom lens
attached to a Nikon D610, allowing for accurate strain measurements taken
later in ImageJ. Samples were then tested after 16, 18, 24, 44, 65, and 94
hours at room temperature post-fabrication, as well as 1 hour after 25 minutes
of curing at an elevated temperature of 65$\degree$C. The manufacturer lists
the cure time as 16 hours.

Figure 7: Averaged values for (a.) the elastic and long-term moduli of Dragon
Skin 30TM and (b.) the ratio of their mean values, as determined from tensile
stress-relaxation tests. The x-axes indicate the time elapsed since
fabrication; The cure time is listed by the manufacturer as 16 hours post-
fabrication. (c.) Plot of the relaxation time from stress-relaxation tests.
The dashed lines in (b.) & (c.) show the averaged values used in our analysis,
$\bar{E}=E_{\infty}/E_{0}=0.91$ and $\tau_{\sigma}=0.78$ s, which were
determined from the samples whose ages were relevant to our experiments (24-48
hours cured at room temperature, and 0.4 hours cured at elevated temperature).
(d.) SLS relaxation modulus ($E(t)$) and inverse creep compliance
($J^{-1}(t)$) functions constructed from our averaged quantitites. Because the
material is nearly elastic, the $E(t)\approx J^{-1}(t)\>\forall\>t>0$. All
error bars represent one standard deviation.

The long-term (equilibrium) modulus was determined from tensile stress-
relaxation tests using the tensile testing machine Instron 5943.
Displacements, which resulted in strains $\varepsilon_{0}\in[2\%,20\%]$, were
imposed at rates of 150-500 mm/min, then maintained for 120 seconds while the
force, which decreases over a timescale $t\sim\tau_{\sigma}$ before
plateauing, was measured. The plateau force was used to calculate
$\sigma_{\infty}$ where we took
$E_{\infty}\equiv\sigma_{\infty}/\varepsilon_{0}$. The values for $E_{\infty}$
plotted in Fig. 7a are calculated as the average of 12 total measurements from
4 tests at varied strain levels and strain rates.

This relatively fast loading resulted in measurement uncertainty at early
times, so we did not fit a curve to the entire range of stress-relaxation
data, nor did we extract the elastic modulus from stress-relaxation tests.
Instead, tensile tests on the same samples were conducted at rates of 10 and
20 mm/min up to $\varepsilon_{\text{max}}\in[8\%,26\%]$. This resulted in
linear stress-strain curves, and the slopes were used to calculate $E_{0}$.
Each data point in Fig. 7a represents the average of 6 measurements from 2
tests for each sample.

Between 16 and 24 hours post-fabrication at room temperature (approximately
20°C), we observe an increase in the long-term modulus from $0.20\pm 0.05$ MPa
to $0.55\pm 0.09$ MPa. By 24 hours, the long-term modulus reaches a plateau.
The elastic modulus follows a similar trend: the material stiffens from
$E_{0}=0.48\pm 0.06$ MPa at 16 hours post-fabrication to $0.60\pm 0.02$ MPa at
24 hours, by which time the elastic modulus has plateaued. The plateau value
we measure is in agreement with the 100% modulus value listed by the
manufacturer, Smooth-On, of 0.59 MPa. More or less the same plateau values
result from curing the sample in the oven at 65°C for 25 minutes
($E_{\infty}=0.53\pm 0.08$ MPa, $E_{0}=0.58\pm 0.05$ MPa).

The ratio $\bar{E}=E_{\infty}/E_{0}$ is central to our analysis, and is
plotted in Fig. 7b. For the shell-buckling experiments discussed in the body
of this paper, the relevant times are 24 and 44 hours, and the 25 minute oven-
cured sample (see Sect. VII.2). Averaging this data gives the value we use
throughout our analysis (and the dashed line in Fig. 7b), $\bar{E}=0.91$.
Instead including all of the samples whose moduli appear to have saturated
(i.e. all except the 16 and 18 hour samples) does not change this averaged
value (it only slightly increases the standard deviation, to 0.3 MPa.)

With the functional form of the relaxation modulus in mind, we identified the
relaxation time $\tau_{\sigma}$ as the time when the modulus has reduced to
$E_{\infty}+E_{1}e^{-1}$. The averaged values are shown in Fig. 7c. The
overall average was $\tau_{\sigma}=0.78\pm 0.49$ s. This corresponds to a
retardation time $\tau_{\varepsilon}\approx 0.86$ s Lakes2009 . The relatively
large error range on the characteristic timescale results from initial
uncertainty in stress-relaxation curves.

In Sects. IV & V, we blur the lines between the creep compliance and the
relaxation modulus: While the active process is creep, we use not the creep
compliance, but the relaxation modulus – which describes stress-relaxation –
in our analysis. This is for two reasons: First, our primary aim is to
describe creep buckling in the language of elastic buckling – which refers to
Young’s Modulus $E$ – rather than to provide an exact description of creep in
this material. Second, we found that stress-relaxation (displacement-
controlled) tests provided much more reliable data than creep (force-
controlled) tests using Instron 5943.

The relaxation modulus (Eq. (2) in the main text) is related by Laplace
transform to the creep compliance (Eq. (3)). Thus, the limiting values
$E_{0}=\sigma(t=0)/\varepsilon$ and $J_{0}=\varepsilon(t=0)/\sigma$;
$E_{\infty}=\sigma(t=\infty)/\varepsilon$ and
$J_{\infty}=\varepsilon(t=\infty)/\sigma$ are exact inverses.

The relaxation time $\tau_{\sigma}$ is less than the retardation time
$\tau_{\varepsilon}$. The two are related through the relaxation strength,
defined as $\Delta=E_{1}/E_{\infty}$, according to
$\tau_{\varepsilon}=\tau_{\sigma}(1+\Delta)$, and thus $E(t)\neq
J^{-1}(t)\forall t$ Lakes2009 . However, the difference between the two
functions is very small, as shown in Fig. 7d. Due to this negligible
difference, we conclude that our approximation in Sect. V is also reasonable.

###  VII.2 Shell fabrication

Figure 8: Shell fabrication process. i. Mixed and degassed DragonSkin 30 is
poured into two hemispherical aluminum cavities of $R_{o}=25$ mm. Only one is
shown in i.-iv., and several layers of tape have been cut into a circle and
adhered to the center of one of the hollow spheres. ii. After degassing again,
the filled cavity is fitted with an alignment sleeve. iii. A half-sphere of
$R_{i}<R_{o}$ is inserted. iv. After curing for 25 minutes in the oven, a
spherical hemisphere has formed. v. Heptane-diluted polymer is deposited with
a syringe to glue two cured hemispheres. vi.-vii. The edges of the two
hemispheres are joined and aligned by the alignment sleeve. viii. After curing
at room temperature for at least $16$ hours, a sealed sphere is removed from
the mold.

The process for fabricating polymeric spherical shells, which was developed in
Djellouli2017 , is as follows: Custom aluminum molds consist of female and
male components along with alignment sleeves. The female mold is a cylinder of
equal radius and height ($30$ mm) with a hemispherical cavity, which sets the
outer radius of the shells to $R_{o}=25$ mm. The male component consists of a
shouldered half-sphere whose size determines the inner shell radius,
$R_{i}=\\{24,23,21.5,20\\}$ mm. The shoulder (chamfered to $5\degree$) is $10$
mm tall, and its maximum diameter matches that of the female mold as well as
the inside of the guiding sleeve. The latter is a $40$ mm tall cylindrical
tube, internally chamfered up to a depth of $10$mm to accommodate the male
mold.

To control the location of the onset of buckling, shells are seeded with a
circular imperfection. This is achieved by affixing 1-4 layers of adhesive
tape cut to $R_{\delta}\approx 6$ mm (resulting in imperfection depth of
$\delta\in\\{0.76,0.81,0.84,0.40\\}$ mm, in order for the thinnest to thickest
shell) to the center of one of the two female molds used to make each shell.

To make shells, the polymer is prepared according to package instructions,
degassed in a vacuum, and poured into the female molds. After degassing again,
each polymer-filled cavity is fitted inside an alignment sleeve, and the male
mold is inserted. The assembly is tightly clamped between the two plates of a
simple mechanical press and cured at 65°C for 25 minutes.

After curing, the alignment sleeve and the male component are removed,
revealing two hemispherical shells resting in the female molds. To join the
two halves, a glue is prepared. The viscosity of the liquid polymer is reduced
via dilution with heptane at a 2:1 ratio. This allows for the application of a
sufficiently thin layer of glue, deposited with a syringe around the equator
of each hemisphere. Again a sleeve is used to align the two halves, whose
contact is ensured by the mechanical press, and the shell is left to cure at
room temperature for 16 hours.

A drill press is used to create a 1 mm diameter hole, into which a small
nozzle connected to a tube allows for internal pressure control. Lastly, a
suction cup ($2-3$ cm diameter) is glued to the shell surface opposite the
buckling spot with cyanoacrylate (Loctite). During experiments, the shell is
fixed in place via a screw attached to the back of the suction cup.

###  VII.3 Knockdown of elastic critical pressure due to through-thickness
defects

Unsurprisingly, the classical prediction for the buckling pressure of a
perfect elastic spherical shell (Eq. (1)) does not capture the behavior of our
imperfect shells. Recently, predictions for the knockdown factor as a function
of the size of an axisymmetric imperfection have been presented for dimple-
like Lee2016 ; Jimenez2017  and through-thickness Yan2020  defects.

Figure 9: Comparison of the theoretical elastic buckling pressure $P_{c}^{e}$
for thin, shells with through-thickness imperfections with our experimental
value. The experimental value is taken as the minimum pressure at which no
slowdown of deformation occurs before buckling. The theoretical value
$k_{d}P_{c}$ is calculated using $k_{d}$ extracted from Fig. 11b in
Ref.Yan2020 , and $P_{c}$ from Eq (1). Error bars correspond to approximately
one standard deviation, and are smaller than the markers in most cases.

The reduced-thickness defects in our experiments (see Fig. 1) are like those
in the work of Yan et al. Yan2020 . The authors present data from experiments
and FEM simulations on the knockdown factor for varied depth, angular width,
and transition width of the imperfection. We use this data333Knockdown values
were extracted from Figure 11b of Ref. Yan2020 , and found to be
$k_{d}=\\{0.2719\pm 0.007,0.7423\pm 0.000,0.8664\pm 0.008,0.93+0.070\\}$ for
$\eta=\\{24.5,12.0,6.4,4.5\\}$, respectively. to calculate the “knocked-down”
theoretical elastic critical pressure for each of our shells. This theoretical
value is plotted against our experimental data in Fig. 9. We have taken
$P_{c}^{e}$ to be the minimum pressure at which a plot of volumetric or pole
deformation versus time shows no slope change until buckling (the rate of
deformation is much slower when creep occurs.)

Up to relatively thick shells, the theoretical knockdown predicts our
experimentally measured critical pressure relatively well
($\eta=24.5,12.0,6.6$) – perhaps surprsingly so, given that the data in Ref.
Yan2020  was collected in quasi-static experiments on shells of fixed
$\eta=100$. However, as $\eta$ decreases even further from the thin shell
limit, the thick shell withstands higher pressures than predicted. This
inapplicability of thin shell knockdown theories is significant for the
thickest shell we tested ($\eta=4.5$). For consistency, then, in all arguments
throughout the main text we rely on an experimentally determined value for the
elastic buckling pressure. Thus, we take the experimental value for
$P_{c}^{e}$ to be the minimum pressure difference where elastic buckling
occurs.

###  VII.4 Methods for dynamic pressure loading experiments

Before each experiment, a vacuum pump (Becker U 4.40) is used to reduce the
pressure inside a 50 liter tank. A flexible tube ($3$ mm inner diameter, $12$
cm length) connects the tank to a differential pressure sensor (Freescale
Semiconductor MPX5100DP, sensitivity 45mV/kPa). Pressure readings were
recorded using a microcontroller (Arduino UNO) every $0.02$ s with a
resolution of $0.1$ kPa.

As this pressure resolution is slightly coarse for the thinner shells, which
buckle at pressures on the order of $1$ kPa, we account for rounding (and
small pressure fluctuations) by reporting the mean and standard deviation of
the pressure values reported at intermediate444For shells subjected to large
$P$ which buckle elastically, we simply report the (non-fluctuating) imposed
pressure, recorded before the shell begins to deform, with a default maximum
error of $\pm 0.05$ kPa. times (after initial elastic deformation and before
buckling, should it occur).

Via a T-junction, a second tube ($6$ cm long) connects the tank to the inside
of the shell, by way of an electrovalve (Matrix Pneumatics Solenoid Valve,
MX891.901C224). The response time of the valve is $<1$ ms. Based on a back-of-
the-envelope calculation of the amount of gas that needs to travel from the
shell to the tank for pressure equilibrium and the corresponding mass flow
rate, we expect that shell-vacuum tank system equilibriates within
approximately $4$ ms for pressures imposed on the thickest shell (where higher
pressure gradients drive faster air flow), to about $40$ ms for the thinnest.

Arduino IDE software enables synchronization of the sudden opening of the
valve with the digital recording of the pressure difference (at a rate of 500
Hz), as well as with the triggering of the high-speed camera (Phantom Miro
310). Depending on the empirically determined delay time before buckling,
images were captured at a rate of 2000-9000 frames per second. Image
processing of TIFF stacks was completed using ImageJ and custom Python
scripts.

## References

  * (1) Lee A, Jiménez FL, Marthelot J, Hutchinson JW, Reis PM. 2016 The geometric role of precisely engineered imperfections on the critical buckling load of spherical elastic shells. Journal of Applied Mechanics 83, 111005.
  * (2) Wagner H, Hühne C, Zhang J, Tang W. 2020 On the imperfection sensitivity and design of spherical domes under external pressure. International Journal of Pressure Vessels and Piping 179, 104015.
  * (3) NASA. 1969 Buckling of thin-walled doubly-curved shells. NASA space vehicle design criteria (structures). NASA SP-8032.
  * (4) Karman TV, Tsien HS. 1939 The buckling of spherical shells by external pressure. Journal of the Aeronautical Sciences 7, 43–50.
  * (5) Tsien HS. 1942 A theory for the buckling of thin shells. Journal of the Aeronautical Sciences 9, 373–384.
  * (6) Koiter WT. 1945 Over de stabiliteit van het elastisch evenwicht. PhD thesis.
  * (7) Koiter WT. 1969 Proc Kon Ned Ak WetB72, 40–123.
  * (8) Hutchinson JW. 1967 Imperfection sensitivity of externally pressurized spherical shells. Journal of Applied Mechanics 34, 49–55.
  * (9) Bushnell D. 1967 Nonlinear axisymmetric behavior of shells of revolution. AIAA Journal 5, 432–439.
  * (10) Krenzke M, Kiernan T. 1965 The effect of initial imperfections on the collapse strength of deep spherical shells. .
  * (11) Koga T, Hoff NJ. 1969 The axisymmetric buckling of initially imperfect complete spherical shells. International Journal of Solids and Structures 5, 679–697.
  * (12) Lee A, Brun PT, Marthelot J, Balestra G, Gallaire F, Reis PM. 2016 Fabrication of slender elastic shells by the coating of curved surfaces. Nature Communications 7.
  * (13) Hutchinson JW. 2016 Buckling of spherical shells revisited. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 472, 20160577.
  * (14) Jiménez FL, Marthelot J, Lee A, Hutchinson JW, Reis PM. 2017 Technical brief: Knockdown factor for the buckling of spherical shells containing large-amplitude geometric Defects. Journal of Applied Mechanics 84, 034501.
  * (15) Yan D, Pezzulla M, Reis PM. 2020 Buckling of pressurized spherical shells containing a through-thickness defect. Journal of the Mechanics and Physics of Solids 138, 103923.
  * (16) Paulose J, Nelson DR. 2013 Buckling pathways in spherical shells with soft spots. Soft Matter 9, 8227.
  * (17) Gerasimidis S, Hutchinson JW. 2020 Dent imperfections in shell buckling: The role of geometry, residual stress, and plasticity. Journal of Applied Mechanics 88.
  * (18) Gerasimidis S, Virot E, Hutchinson JW, Rubinstein SM. 2018 On Establishing Buckling Knockdowns for Imperfection-Sensitive Shell Structures. Journal of Applied Mechanics 85.
  * (19) Hutchinson JW, Thompson JMT. 2017 Nonlinear buckling behaviour of spherical shells: barriers and symmetry-breaking dimples. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 375, 20160154.
  * (20) Audoly B, Hutchinson JW. 2020 Localization in spherical shell buckling. Journal of the Mechanics and Physics of Solids 136, 103720. The Davide Bigoni 60th Anniversary Issue.
  * (21) Hutchinson JW. 2020 EML Webinar overview: New developments in shell stability. Extreme Mechanics Letters 39, 100805.
  * (22) Holmes DP. 2019 Elasticity and stability of shape-shifting structures. Current Opinion in Colloid & Interface Science 40, 118–137.
  * (23) Reis PM. 2015 A perspective on the revival of structural (in)stability with novel opportunities for function: From buckliphobia to buckliphilia. Journal of Applied Mechanics 82.
  * (24) Sacanna S, Irvine WTM, Chaikin PM, Pine DJ. 2010 Lock and key colloids. Nature 464, 575–578.
  * (25) Shim J, Perdigou C, Chen ER, Bertoldi K, Reis PM. 2012 Buckling-induced encapsulation of structured elastic shells under pressure. Proceedings of the National Academy of Sciences 109, 5978–5983.
  * (26) Keplinger C, Li T, Baumgartner R, Suo Z, Bauer S. 2012 Harnessing snap-through instability in soft dielectrics to achieve giant voltage-triggered deformation. Soft Matter 8, 285–288.
  * (27) Gorissen B, Melancon D, Vasios N, Torbati M, Bertoldi K. 2020 Inflatable soft jumper inspired by shell snapping. Science Robotics 5.
  * (28) Yang D, Verma MS, So JH, Mosadegh B, Keplinger C, Lee B, Khashai F, Lossner E, Suo Z, Whitesides GM. 2016 Buckling pneumatic linear actuators inspired by muscle. Advanced Materials Technologies 1, 1600055.
  * (29) Lee A, Yan D, Pezzulla M, Holmes DP, Reis PM. 2019 Evolution of critical buckling conditions in imperfect bilayer shells through residual swelling. Soft Matter 15, 6134–6144.
  * (30) Yan D, Pezzulla M, Cruveiller L, Abbasi A, Reis PM. 2020 Magneto-active elastic shells with tunable buckling strength. .
  * (31) Holmes DP, Lee JH, Park HS, Pezzulla M. 2020 Nonlinear buckling behavior of a complete spherical shell under uniform external pressure and homogenous natural curvature. Physical Review E 102.
  * (32) Reis PM, Jaeger HM, van Hecke M. 2015 Designer matter: A perspective. Extreme Mechanics Letters 5, 25–29.
  * (33) Djellouli A, Marmottant P, Djeridi H, Quilliet C, Coupier G. 2017 Buckling instability causes inertial thrust for spherical swimmers at all scales. Phys. Rev. Lett. 119, 224501.
  * (34) Mokbel M, Djellouli A, Quilliet C, Aland S, Coupier G. 2021 Post-buckling dynamics of spherical shells. Forthcoming.
  * (35) Sieber J, Hutchinson JW, Thompson JMT. 2019 Nonlinear dynamics of spherical shells buckling under step pressure. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 475, 20180884.
  * (36) Lin IK, Ou KS, Liao YM, Liu Y, Chen KS, Zhang X. 2009 Viscoelastic Characterization and Modeling of Polymer Transducers for Biological Applications. Journal of Microelectromechanical Systems 18, 1087–1099.
  * (37) Yuk H, Zhao X. 2017 A New 3D Printing Strategy by Harnessing Deformation, Instability, and Fracture of Viscoelastic Inks. Advanced Materials 30, 1704028.
  * (38) Urbach EY, Efrati E. 2020 Predicting delayed instabilities in viscoelastic solids. Science Advances 6, eabb2948.
  * (39) Zener C. Elasticity and anelasticity of metals.
  * (40) Lakes RS. 2009 Viscoelastic materials. Cambridge ; New York: Cambridge University Press. OCLC: ocn286432282.
  * (41) Kaplan A, Fung Y. 1954 A nonlinear theory of bending and buckling of thin elastic shallow spherical shells. .
  * (42) Taffetani M, Jiang X, Holmes DP, Vella D. 2018 Static bistability of spherical caps. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 474, 20170910.
  * (43) Rosenthal D, Baer H. 1951 An elementary theory of creep buckling of columns. In Journal of Applied Mechanics-Transactions of the ASME vol. 18 pp. 321–321. ASME.
  * (44) Hayman B. 1978 Aspects of creep buckling. I. The influence of post-buckling characteristics. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 364, 393–414.
  * (45) Hayman B. 1981 Creep buckling — A general view of the phenomena. In Creep in Structures pp. 289–307. Springer Berlin Heidelberg.
  * (46) Hoff NJ. 1973 Creep buckling of plates and shells. In Theoretical and Applied Mechanics pp. 124–140. Springer Berlin Heidelberg.
  * (47) Minahen T, Knauss W. 1993 Creep buckling of viscoelastic structures. International Journal of Solids and Structures 30, 1075–1092.
  * (48) Gerard G, Gilbert AC. 1958 A critical strain approach to creep buckling of plates and shells. Journal of the Aerospace Sciences 25, 429–434.
  * (49) Miyazaki N, Hagihara S. 2015 Creep buckling of shell structures. Mechanical Engineering Reviews 2, 14–00522–14–00522.
  * (50) Hoff NJ. 1954 Forty-first Wilbur Wright memorial lecture: Buckling and stability. The Journal of the Royal Aeronautical Society 58, 3–52.
  * (51) Hoff NJ. 1956 Creep buckling. Aeronautical Quarterly 7, 1–20.
  * (52) Huang NC. 1967 Nonlinear creep buckling of some simple structures. Journal of Applied Mechanics 34, 651–658.
  * (53) Gerard G. 1956 A creep buckling hypothesis. Journal of the Aeronautical Sciences 23, 879–882.
  * (54) Hoff NJ. 1975 The effect of geometric nonlinearities on the creep buckling time of axially compressed circular cylindrical shells. Journal of Applied Mechanics 42, 225–226.
  * (55) Obrecht H. 1977 Creep buckling and postbuckling of circular cylindrical shells under axial compression. International Journal of Solids and Structures 13, 337–355.
  * (56) Huang NC. 1965 Axisymmetrical creep buckling of clamped shallow spherical shells. Journal of Applied Mechanics 32, 323–330.
  * (57) Shi JJJ, Johnson CD, Bauld NR. 1970 Application of the variational theorem for creep of shallow spherical shells. AIAA Journal 8, 469–476.
  * (58) Jones N. 1976 Creep buckling of a complete spherical shell. Journal of Applied Mechanics 43, 450–454.
  * (59) Miyazaki N, Yagawa G, Ando Y. 1977 A parametric analysis of creep buckling of a shallow spherical shell by the finite element method. Nuclear Engineering and Design 41, 257–263.
  * (60) Xirouchakis P, Jones N. 1980 Axisymmetric and bifurcation creep buckling of externally pressurised spherical shells. International Journal of Solids and Structures 16, 131–148.
  * (61) Leckie FA, Hayman B. 1976 Creep instability of thick shell structures. In Buckling of Structures pp. 86–94. Springer Berlin Heidelberg.
  * (62) Pandey A, Moulton DE, Vella D, Holmes DP. 2014 Dynamics of snapping beams and jumping poppers. EPL (Europhysics Letters) 105, 24001.
  * (63) Gomez M, Moulton DE, Vella D. 2016 Critical slowing down in purely elastic ‘snap-through’ instabilities. Nature Physics 13, 142–145.
  * (64) Gomez M, Moulton D, Vella D. 2018 Dynamics of viscoelastic snap-through. Journal of the Mechanics and Physics of Solids 124, 781–813.
  * (65) Coupier G, Djellouli A, Quilliet C. 2019 Let’s deflate that beach ball. The European Physical Journal E 42.
  * (66) Brinkmeyer A, Santer M, Pirrera A, Weaver P. 2012 Pseudo-bistable self-actuated domes for morphing applications. International Journal of Solids and Structures 49, 1077 – 1087.
  * (67) Santer M. 2010 Self-actuated snap back of viscoelastic pulsing structures. International Journal of Solids and Structures 47, 3263 – 3271\.
  * (68) Urbach EY, Efrati E. 2018 The metric description of viscoelasticity and instabilities in viscoelastic solids. .
  * (69) Che K, Rouleau M, Meaud J. 2019 Temperature-tunable time-dependent snapping of viscoelastic metastructures with snap-through instabilities. Extreme Mechanics Letters 32, 100528.
  * (70) Che K, Yuan C, Qi HJ, Meaud J. 2018 Viscoelastic multistable architected materials with temperature-dependent snapping sequence. Soft Matter 14, 2492–2499.
  * (71) Dykstra DMJ, Busink J, Ennis B, Coulais C. 2019 Viscoelastic snapping metamaterials. Journal of Applied Mechanics 86.
  * (72) Hutchinson JW, Thompson JMT. 2018 Imperfections and energy barriers in shell buckling. International Journal of Solids and Structures 148-149, 157–168.