Predicting delayed instabilities in viscoelastic solids

See allHide authors and affiliations

Science Advances  04 Sep 2020:
Vol. 6, no. 36, eabb2948
DOI: 10.1126/sciadv.abb2948


Determining the stability of a viscoelastic structure is a difficult task. Seemingly stable conformations of viscoelastic structures may gradually creep until their stability is lost, while a discernible creeping in viscoelastic solids does not necessarily lead to instability. In lieu of theoretical predictive tools for viscoelastic instabilities, we are presently limited to numerical simulation to predict future stability. In this work, we describe viscoelastic solids through a temporally evolving instantaneous reference metric with respect to which elastic strains are measured. We show that for incompressible viscoelastic solids, this transparent and intuitive description allows to reduce the question of future stability to static calculations. We demonstrate the predictive power of the approach by elucidating the subtle mechanism of delayed instability in thin elastomeric shells, showing quantitative agreement with experiments.


The snapping of the Venus fly-trap leaf, one of the fastest motions in plant kingdom, is preceded by a relatively slow creeping motion (1). A similar creep is observed before the snap through of thin elastomeric shells known as jumping poppers (2). While the snap through itself lasts only a fraction of a second (2, 3), the slow creeping motion, during which the shells seems to be elastically stable, may last orders of magnitude longer. On a much larger scale, similar behavior is, in some cases, attributed to Earth’s crust before an earthquake aftershock (4, 5). The exact role of viscoelasticity in aftershocks is not fully understood, partially because of the absence of predictive theoretical framework of future stability in these systems. In all of the abovementioned systems, the slow viscoelastic flow in the material leads the system to an instability that abruptly releases the internally stored elastic energy. While in most cases we are able to explicitly write down the equations governing the viscoelastic behavior, the mechanism of delayed instabilities in viscoelastic solids remains poorly understood.

Viscoelastic solids display reversible elastic behavior when loaded at a fast rate. However, they also show a slow creeping behavior under a constant load and exhibit stress relaxation when held at constant displacement. Commonly, these solids are modeled by a constitutive law relating the stress rate to the stress, strain, and strain rate or equivalently by expressing the stress as a function of the history of strain rate and a material-dependent memory kernel (6). While these approaches capture the material response well and yield accurate results through simulations of viscoelastic structures (7, 8), they are rarely explicitly solvable and provide little intuition regarding the state of the viscoelastic material and, in particular, have little in the way of clarifying the governing processes in viscoelastic instabilities. Recent attempts to address viscoelastic instabilities by modeling the viscoelastic response as an elastic medium with temporally evolving stiffness (79) show only a qualitative agreement with experiment and have limited applicability as, in particular, they cannot capture creeping at zero load.

In this work, we quantitatively address the phenomenon of viscoelastic instability through a metric description. In this metric approach, we describe the materials’ behavior as a fast elastic response with respect to temporally evolving rest lengths that change because of the slow viscoelastic flow. The microscopic response in the material that leads to stress relaxation is interpreted as the evolution of the rest lengths of the system toward the lengths assumed in its present state (see Fig. 1). The instantaneous rest lengths in the system serve as a state variable accounting for the slow and inelastic evolution in the material. This approach is somewhat reminiscent of the additive decomposition of strains (10, 11); however, using the instantaneous reference metric rather than strain as the basic ingredient of the theory allows for a much more transparent treatment that, in particular, allows to predict future stability of unconstrained viscoelastic structures.

Fig. 1 The viscoelastic reference length evolution.

At the resting state, all three length measures on the body, its measured length g (marked red), its instantaneous reference length g¯ (marked gray), and its rest reference length g¯0 (marked black) are all equal. When subjected to a constant displacement extension, the instantaneous reference length evolves away from the rest length and toward the presently assumed length, thus resulting in stress relaxation. It asymptotically approaches the stationary state g¯stat=βg+(1β)g¯0, in which the initial stress is reduced by a factor of 1 − β. When released, the unconstrained system immediately adopts its favored instantaneous reference length, which, in turn, gradually creeps toward the rest lengths.

The key to our analysis of stability in viscoelastic materials is the discovery of a new stationarity property; if an incompressible viscoelastic system starting at rest is brought abruptly to another locally stable state, then despite the continuous evolution of the internal stresses in the structure, it will display no motion and, in particular, will never lose stability. This result both elucidates the subtle nature of delayed instabilities in these systems and paves the path for their quantitative understanding.


Linear viscoelastic continua

All linear viscoelastic material constitutive relations and, in particular, the relations for every spring dashpot model, can be expressed through a convolution of the strain rate with a stress relaxation function, σ(t)=tG(ts)ε.(s)ds (6, 12). For one-dimensional systems, this relation can be easily brought into the form (13)σ(t)=C(ε(t)+β tϕ.(ts)ε(s)ds )(1)

Here, C is the material stiffness and ϕ(s) is the normalized memory kernel that satisfies ϕ(0) = 1, ϕ(∞) = 0, and as we expect, the memory to only decay in time also satisfies ϕ.(s)<0. 0 < β < 1 is a dimensionless constant of the viscoelastic material that measures the extent of viscoelasticity. For a general three-dimensional isotropic and incompressible material, the tensorial generalization of Eq. 1 reads (14)Sij(t)=Cijkl(εkl(t)+β tϕ.(ts)εkl(s)ds)(2)

Here, Cijkl is the isotropic stiffness tensor, which, in three dimensions, depends only on two constants, often interpreted as the Young’s modulus, Y, and the Poisson’s ratio, ν. As the material is assumed to be incompressible, ν = 1/2 for both the elastic and viscoelastic response tensors, rendering them proportional to each other. In this case, the viscoelastic evolution is fully captured by a single scalar function. Sij(t) denotes the second Piola-Kirchhoff stress tensor, which, through the equation above, depends on the full history of the the strain tensorεij(t)=12(gij(t)g¯ij0)(3)

The metric gij measures lengths in the material and uniquely describes the configuration of the elastic body. We have also introduced the rest reference metric, g¯ij0, the metric at which the system is locally stress free and stationary. For a more detailed description of covariant elasticity using metric tensors, the reader is referred to (14, 15).

Equation 2 predicts that instantaneous incremental deformations, Δgij, lead to linear stress increments, ΔSij=Cijkl12Δgkl, as one would predict for a purely elastic response. Inspired by this instantaneous elastic-like response, we define the time-dependent instantaneous reference metric g¯ij of the body such as to satisfySij(t)=Cijkl12(gkl(t)g¯kl(t))(4)

The temporal evolution of the instantaneous reference metric can be deduced from Eqs. 2 to 4 and readsg¯ij(t)=(1β) g¯ij0β tϕ.(ts) gij(s) ds(5)

It is important to stress that Eqs. 4 and 5 are completely equivalent to Eq. 2. However, the notion of an instantaneous reference metric provides a more transparent description of the viscoelastic system and allows a more intuitive understanding of its dynamics. Note that g¯ij(t) remains unchanged by instantaneous variations of gij. We may thus consider it as the slow-state variable describing the viscoelastic evolution of the material. At each moment in time, we may consider the system as an elastic system with respect to the instantaneous reference metric g¯ij (16, 17). The dimensionless factor β can be shown by Eq. 2 to quantify the fraction of stress asymptotically relaxed in a constant displacement experiment starting from rest. β is thus a measure of the degree of viscoelasticity in the system and can be directly deduced from simple measurements (see the Supplementary Materials). For a prescribed configuration given by the metric gij, the instantaneous reference metric g¯ij will asymptotically approach the stationary solution of Eq. 5g¯ijstat=βgij+(1β)g¯ij0(6)

Therefore, β also measures the material’s propensity to change its instantaneous reference metric toward the realized metric in its present configuration. β = 0 describes a purely elastic material, where the instantaneous reference metric remains g¯ij=g¯ij0 for all times (or if starting away from it, relaxes to g¯ij0 independently of the behavior of the body gij). β = 1 describes a viscoelastic fluid where the rest reference metric g¯ij0 has no meaning and the material approaches arbitrarily close to gij, diminishing the stress to zero in relaxation.

The quasi static approximation

Viscoelastic systems are dissipative; thus, the notion of an elastic free energy is ill defined. Nonetheless, the virtual work of displacements performed over a short period Δt → 0, coincides with the instantaneous elastic energy functional (16)E[εijel]=12Cijklεijelεklelg¯0 d3x(7)where the elastic strain is εijel=12(gijg¯ij). Typically, the elastic response time scale in elastomers is much smaller than the viscoelastic relaxation time [sometimes referred to as the large Deborah number limit (18)]. In these cases, we can eliminate inertia from the system and approximate the motion of the material by the quasi-static evolution between elastic equilibrium states. That is, the configuration at every instance in time, given by the metric gij(t), minimizes the instantaneous elastic energy functional Eq. 7. This yields the covariant elastic Euler-Lagrange equationsjSij+ΓjkiSjk+(Γ¯0)jkkSij=jSij+((Γ¯0)jkkΓjkj)Sij=0(8)where Γ and Γ¯0 denote the Christoffel symbols with respect to the metrics g and g¯0, respectively. This is reminiscent of the obtained formula for incompatible elasticity (17) with the exception that the volume form (and the thus the relevant Christoffel symbol in the variation) is inherited from g¯0 rather than g¯.

Viscoelastic instabilities through the metric description

A given instantaneous reference metric g¯ij can yield multiple elastically stable configurations in the instantaneous elastic energy functional Eq. 7. As g¯ij evolves viscoelastically according to Eq. 5, it could acquire new stable configurations, merge existing stable points, or cause stable elastic configurations to lose their stability. The latter phenomenon manifests as a slow viscoelastic evolution followed by a rapid snap to a near-stable configuration and forms the main difficulty in predicting the stability of viscoelastic structures. This phenomenon was termed temporary bistability (7), pseudo bistability (8), or creep buckling (19, 20). Each of these terms captures some of the essence of the phenomenon, but not all. A locally stable state arrived to instantaneously from rest will remain stable for all times. We therefore argue that for an incompressible linearly viscoelastic solids to creep into instability, two distinct processes need to take place: First, an elastically unstable state needs to acquire stability through viscoelastic relaxation under some external load exerted for a finite time. The second process commences with the removal of the load as the body assumes the newly acquired stable state and consists of a viscoelastic creep that inevitably results in instability. An acquired stability state cannot remain stable, indefinitely rendering this stability transient. We next use the metric description of viscoelasticity to prove these results, which, when joined together, elucidate the mechanism governing the stability viscoelastic structures.

Locally stable configuration stationarity. The first result, termed locally stable configuration stationarity, states that if a system at rest is brought instantaneously to a locally stable configuration, then not only it will never lose stability but also the configuration will remain stationary; i.e., the configuration will display no temporal variation despite the continuous evolution of the instantaneous reference metric and the relaxation of the corresponding stress. To prove this result, consider a system at rest, for which g=g¯=g¯0, that is brought instantaneously at time t = 0 to another configuration g=g*g¯, which is elastically stable, i.e., satisfies Eq. 8. It is easy to see by direct substitution in Eq. 5 that the constant solution g(t) = g* causes the instantaneous reference metric to evolve according tog¯ij(t)=(1β)g¯ij0βg¯ij0 0ϕ.(ts) dsβgij* 0tϕ.(ts) ds=(1β)g¯ij0+βgij*+βϕ(t) (g¯ij0gij*)=g¯*stat+β(g¯ij0gij*)ϕ(t)

This evolution of the reference metric, in turn, causes the stress to diminish by a uniform time dependent factor according toSij(t)=Cijkl12((1β)+βϕ(t)) (gij*g¯ij0)=((1β)+βϕ(t))Sij(0)(9)

It is thus straightforward to see that g* remains the mechanical equilibrium despite the evolution of ḡ, and no external force is required to fix g = g*. As S(0) satisfies mechanical equilibrium condition Eq. 8, so will S(t), and the state of the system will remain in equilibrium for all t. A simple and graphical interpretation of the metric collinearity which leads to the locally stable configuration stationarity is provided in Fig. 2.

Fig. 2 Schematic representation of the metrics collinearity.

The minimization of the metric g (marked by a full black circle) is constrained and performed with respect to the subset of metrics that correspond to realizable configurations (thick black line). These metrics are in, particular, orientation preserving and Euclidean. Given an instantaneous reference metric, g¯ (marked by a full gray circle), the realized metric will correspond to the closest point from the set of admissible metrics to g¯ according to the distance function given by the instantaneous elastic energy Eq. 7. Starting from rest, g¯ evolves from g¯0 (marked by a full red circle) toward the g, which remains the closest admissible metric to g¯ due to collinearity of the three metrics. As g remains stationary, the evolution of g¯ will preserve the collinearity, asymptotically approaching g¯stat (marked by an open circle), which is also collinear. We stress that throughout this evolution, g remains unchanged; thus, no variation of the configuration will be observed despite the stress relaxation.

Conversely, if a system converges to a fixed stable state (wherein the configuration and the instantaneous reference metric do not evolve), then the corresponding configuration must also be extremal with respect to the rest reference metric. This again could be shown by substituting Eq. 6 into Eq. 8. We note that these results are not limited to a particular memory kernel. Intuitively, both claims are due to collinearity of gij,ḡij0 and ḡij under the appropriate initial conditions, as described in Fig. 2.

Transient acquired elastic stability. The notion of locally stable state stationarity does not imply that viscoelastic systems cannot creep into instabilities, but rather constrains the settings in which these instabilities can occur. Starting from rest, a configuration may be elastically stable or elastically unstable. Neither case shows a slow creeping motion. Thus, if a system displays creeping, then it was not brought to its state from rest; it must have been held under external load for some time, allowing its instantaneous reference metric to evolve. This evolution of the instantaneous reference metric may lead unstable configurations to acquire stability. These acquired stability configurations are configurations that, if arrived to from rest, will be unstable (7, 8). However, after a long enough holding period under external load, they will become locally elastically stable. The instantaneous reference metric in this case would not remain stationary, and its evolution will manifest as the slow creeping of the configuration. It is straightforward to see that this creeping motion must lead to instability, rendering the acquired stability transient. Assume in contradiction that an acquired stability configuration becomes stationary, i.e., both g and g¯ cease to evolve in time, then much like in Eq. 9, it can be shown that the stress with respect to g¯ coincides, up to a constant factor, with the one measured with respect to g¯0. Thus, the same state must also be extremal with respect to g¯0 in contradiction with the assumption that the state was initially unstable.


The results presented above explain many of the qualitative phenomenology of viscoelastic instabilities. We next come to test the quantitative predictions of the theory by experimentally examining the response of silicone-rubber conical poppers. We cast silicon rubber poppers in the geometry of truncated conical shells (Fig. 3A; see also the Supplementary Materials). The conical popper geometry behaves similarly to the spherical poppers, yet this geometry allows simpler control over the thickness of the sheet. The conical shells have an apex angle of 45, inner and outer radius radii rmin,rmax, and thickness h. Sufficiently thin poppers show bistability with the inverted shape close to a mirror image of the rest state. As the thickness is increased, this bistability is lost, and when brought from rest to the inverted state, the popper immediately snaps back. If the thickness is large enough, then this instability will persist regardless of how long we hold the conical shell in its inverted state. For intermediate values of the thickness, we expect to observe unstable states that could acquire stability if held long enough in their inverted state. This acquired stability is expected to be lost in a finite time. These three phases are plotted in Fig. 3 (B and C).

Fig. 3 Experimental verification of the viscoelastic stability diagram.

(A) Straight and inverted conical poppers. Photo credit: Erez Y. Urbach, Weizmann Institute. (B) The two axes span the dimensionless geometrical properties of the truncated conical poppers. The background colors represent the theoretically predicted regions of each of the phases. Each marker corresponds to a different popper; different shaped (and colored) markers indicate the different phases observed in experiment. (C) Numerically calculated flipping time as a function of the normalized thickness of the conical popper for immediate release and long holding time. The different poppers were simulated by varying their thicknesses and constant radii rmin = 10 mm, rmax = 25 mm. The material properties taken were β = 0.1, and the memory kernel was assumed to be exponential with τ = 0.1 s, Young’s modulus E = 2.5 MPa, and Poisson’s ratio v = 0.47. Varying the kernel may lead to varying rate of divergence of the flipping time between the stable and acquired stability region, yet the location of this divergence will remain unchanged. The divergence of flipping times was addressed in (8), and more recently, the rate of divergence was studied in (18).

The conical shell popper geometry is completely determined by two dimensionless variables, h/rmin and rmax/rmin. The boundaries between the different stability phases of Fig. 3C depend on these geometric parameters and on the material constant β. We produced ∼50 different conical popper of different geometries and tested their phases. All popper were produced from the same material and share the same β. Starting from rest gij=g¯ij=g¯ij0, the poppers were abruptly brought to the inverted shape and held for time thold ∼ 1 min. The poppers were then released and their stability recorded.

Theoretically determining the phase boundaries of delayed stability can be achieved by performing only static calculations, which we solved numerically (see the Supplementary Materials). Determining the line separating the stable configurations from unstable configurations (top line separating the white and blue regions in Fig. 3) is carried out by stability analysis with respect to g¯=g¯0. This phase boundary is independent of the viscoelastic properties of the material. We next need to determine which of the unstable states acquired stability after being held for thold in the inverted state ginverted. This is performed by considering the stability with respect to the instantaneous reference metric at time thold given byg¯(thold)=g¯0βA(thold)(g¯0ginverted)where A(t) = (1 − ϕ(t)) starts at zero at t = 0 and asymptotically approaches unity. The value of A(thold) can be directly measured in a uniaxial experiment. Moreover, in most cases, for long enough holding time (or for fast enough relaxation mechanisms), A ∼ 1, rendering the degree of viscoelasticity, β, the most important parameter for predicting asymptotic stability. For a variety of polymeric elastomers, the value β assumes is roughly around 1/10. For the silicone rubber used in our experiment, we measured β = 0.083 ± 0.005 through a stress relaxation experiment (see the Supplementary Materials). The value of A(thold) used for the stability calculations can be directly measured as well. This phase boundary is sensitive to variations in the material properties and the aging protocol; for example, holding the thin shell for a longer period would shift the curve down. Figure 3B shows the experimentally obtained phases superimposed on the phases predicted from the theory, showing good agreement.


The metric description of viscoelasticity presented here is closely related to the additive decomposition of strains in elastoplasticity (11). One can further generalize the theory presented here to account for plasticity or growth within the material by allowing g¯ij0(t) to vary in time, resulting in a covariant viscoelastoplastic theory. The configurational response in these growth processes will show the evolution of g¯ij rather than directly that of g¯ij0. Thus, the viscoelastic response may cause a simple growth rule to appear complex or cause an abrupt and localized growth event to seem diffuse and gradual.

When implemented to isotropic incompressible viscoelastic solids, the metric theory provides the basic rules for viscoelastic instabilities. For a given structure to creep into instability, the creeping must have been preceded by an aging stage where the structures were held under external load for a finite time. When viscoelastic structures are held fixed in a nonextremal configuration, stress relaxation could lead to new acquired stable configurations, yet this acquired stability must be lost in time. The general metric theory of viscoelasticity applies to all linear-viscoelastic materials regardless of the form of their memory kernel. However, the result of locally stable configuration stationarity was obtained under the restricting assumptions that the Poisson’s ratio associated with the elastic and viscoelastic response are similar, and that the system is brought instantaneously from its rest state to the locally stable state. The former assumption is plausible for incompressible elastomers and gels. The latter assumption, however, is expected to be violated by every physical experiment. While one could not prove full stationarity if these assumptions are lifted, the material’s response is continuous in the deviations from the idealized conditions and small deviations will result in very little motion. In particular, in the experiments that we performed, no discernible motion was observed for locally stable states.

The theory proves particularly powerful when applied to describe the delayed instability in elastomeric thin shells. We focused on stability classification in the present work, yet the theory is capable of quantitatively describing the nontrivial dependence of jumping time on the holding time in the inverted state and, in particular, its divergence near the boundary between the unstable regimes and stable regime in Fig. 3C. In the vicinity of this boundary, very short holding times could lead to arbitrarily long delays before the inevitable instability.

While in some specific cases we have compelling evidence for the role of viscoelasticity in the delayed triggering of earthquake aftershocks (4), in the general case, its role as the instigator of delayed instabilities is still unknown. In lieu of a theoretical framework for viscoelastic instabilities, the evidence for the relevance of viscoelasticity to aftershocks mostly relies on numerical simulations relevant to a specific setting. The metric description proposed here aims to provide a theoretical framework for these delayed viscoelastic instabilities. Slip events are described as fast variations in g¯0. Consequently, the instantaneous reference metric g¯ creeps viscoelastically. However, instead of attempting to infer the future stability by examining the creeping motion, we focus on estimating the variation in the rest reference metric; determining asymptotic stability of a viscoelastic system requires understanding the stability landscape with respect to g¯0 and metrics in its vicinity. Thus, the dynamical question of future stability in a creeping system is substituted by examining the stability of a few appropriately chosen elastic static systems.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank D. Biron, G. Cohen, A. Grosberg, S. M. Rubinstein, E. Sharon, Y. Bar-Sinai, and D. Vella for helpful discussions and acknowledge the Aspen Center for Physics, which is supported by National Science Foundation grant PHY-1607761, for hosting many of these discussions. Funding: This work was funded by the Minerva foundation grant number 712273 and by the ISF grant number 1479/16. E.E. acknowledges the supported by the Ascher Foundation and by the Alon Fellowship. Author contributions: This work was initiated and designed by E.E. The theory and experiments were developed by E.E. and E.Y.U and performed by E.Y.U. Simulations were performed by E.Y.U. Both authors participated in the data analysis and contributed to the writing of this manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article