The Extreme Mechanics of Viscoelastic Metamaterials

Mechanical metamaterials made of flexible building blocks can exhibit a plethora of extreme mechanical responses, such as negative elastic constants, shape-changes, programmability and memory. To date, dissipation has largely remained overlooked for such flexible metamaterials. As a matter of fact, extensive care has often been devoted in the constitutive materials' choice to avoid strong dissipative effects. However, in an increasing number of scenarios, where metamaterials are loaded dynamically, dissipation can not be ignored. In this review, we show that the interplay between mechanical instabilities and viscoelasticity can be crucial and can be harnessed to obtain new functionalities. We first show that this interplay is key to understanding the dynamical behaviour of flexible dissipative metamaterials that use buckling and snapping as functional mechanisms. We further discuss the new opportunities that spatial patterning of viscoelastic properties offer for the design of mechanical metamaterials with properties that depend on loading rate.

In many of these dynamical responses, dissipation plays a key role.Specifically, dissipation plays a key role in vibration transmission.As such, dissipative vibration transmissions has been studied extensively in acoustic metamaterials featuring wave propagation and band gap structures [21,22,23,24,25,26,27,28,29,30] as well vibration and shock control [31,32,33,34,35,36].However, although the effect of the constitutive materials' dissipation on mechanical instabilities has been well studied [37,38,39,40,41,42,43,44], the role of dissipation in the instabilities that underpin the extreme mechanics of mechanical metamaterials remains poorly understood (Fig 1).Yet, there is growing evidence that dissipation can play a key role in the response of mechanical metamaterials that are based on buckling [45,46].In addition, a growing body of work has recently demonstrated that a wide variety of additional functionalities can be obtained by leveraging dissipation in mechanical metamaterials [47,46,48,49].
In this short review, we discuss the interplay between mechanical instabilities and viscoelasticity and how this interplay can be harnessed to offer new functionalities that depend on loading rate.

Delayed Snap-through Buckling by Dissipation
Mechanical metamaterials often harness mechanical instabilities, including snap-through instabilities [50].In particular, multistable buckling paths can be harnessed to switch between different responses [50].As such, we need to understand how dynamics and dissipation affect mechanical (snap-through) instabilities and multistability before we can develop new types of dissipative metamaterials.
On the one hand, while quasi-static buckling analyses do not consider any timescale, even purely elastic snap-through instabilities show a very distinct time dependence, transitioning from a very slow state of storing elastic energy to a fast snap-through motion [51].This snap-through motion is what allows for instance the Venus flytrap or the hummingbird to catch their prey [52,53].Such snap-through motion can further be tuned by controlling the loading rate [54].
However, most solid materials are not purely elastic.Most solid materials are viscoelastic instead, and exhibit both viscous and elastic characteristics when undergoing deformation.Viscous materials, like liquids, start flowing when a stress is applied.Purely elastic materials return to their original state once the stress is removed without dissipating any energy.Viscoelastic materials in turn exhibit dissipation and a time-dependent stress-strain response.Specifically, viscoelastic materials stiffen with increasing strain rates.Basic models used for viscoelastic solids are discussed in the Materials & Methods.
Adding dissipation complicates the dynamics of instabilities.Viscoelastic dissipation in instabilities has already been explored extensively since the pioneering works of Kempner in 1954 for beam buckling [37] and Nachbar and Huang in 1967 for snap-through instabilities [38].As one might expect, they showed that the buckling force increases with increasing loading rate, in line with what one might expect with linear viscoelastic materials.
From a qualitative perspective, an equivalent quasi-static analysis can offer significant insights, such as the Von Mises Truss set-up explored by Santer in Fig. 2A [39].Here, a rigid Von Mises truss is suspended by a vertical spring A and a horizontal spring B on the top and right, with spring stiffness K A and K B respectively.When pulling down the top of the Von Mises truss, it can be argued that increasing the ratio of K A /K B is equivalent to increasing the strain rate in a viscoelastic analysis: spring A experiences a much larger deflection rate than spring B when the Von Mises truss is deflected downward, and since both springs are viscoelastic, the instantaneous spring constant or spring A is effectively larger than that of spring B. The relation between force P, deflection delta and K A /K B is then captured in Fig. 2B [39].Two major effects can be identified in Fig. 2B.First of all, as K A /K B increases, not only does the buckling force of the first limit point increase, but so does its deflection.Second, as K A /K B increases, the second limit points moves upward and stops crossing 0 above a critical K A > K crit .At this point, the structure loses its bistability and becomes monostable.In other words, a viscoelastic Von Mises Truss can be bistable if it is loaded slowly but monostable when it is loaded fast.
The opposite effect can also occur.A viscoelastic structure can appear stable in its snapped state after it is loaded.But after viscoelastic relaxation occurs, the snapped state will lose its stability and the structure will snap back into the initial state.This behaviour is in fact very common in viscoelastic snap-through instabilities and has been referred to as temporary bistability [39], pseudo-bistability [40,41,42], acquired bistability [43] or metastability [40].For instance, viscoelastic jumping poppers (Fig. 2C) can appear bistable initially but may still snap back to the original state due to viscoelasticity [42].
However, while a quasi-static analysis can show some agreement with viscoelastic analyses, it can also differ substantially.To compare, we model the same Von Mises Truss (A) A single rigid truss supported at both ends by elastic springs (Reproduced from [39] with permission from Elsevier).Changing the ratio of the spring stiffnesses can be considered equivalent to increasing the loading rate in a viscoelastic analysis.The equilibrium paths in (B) (Reprinted from [39] with permission from Elsevier) show that both monostable and bistable configurations are possible.(C) A viscoelastic jumping popper shows pseudo-bistability (Reproduced from [42] under CC BY 4.0).(D) A viscoelastic version of the Von Mises Truss of (A), featuring Generalised Maxwell springs with 5 spring-dashpots (Modified from [39] with permission from Elsevier).(E) Force-Displacement curves due to moving at constant speed, with the colour indicating the total loading time in seconds.Both very high and very low speeds induce bistability, while intermediate speeds can generate monostability.
of Fig. 2A but by fully modelling the viscoelasticity of the springs (Fig. 2D) so that they faithfully describe the behaviour of a 3D printed rubber, which has also been explored in depth analytically by Gomez et al. [42].We do so with the Viscoelastic Solver in Abaqus 2021 (See Sec-tion5), using a three element truss model (T2D2).Specifically, we replace springs A and B in Fig. 2A by viscoelastic springs using a viscoelastic 5-term Prony series, which is equivalent to a Generalised Maxwell-Model with 5 springdashpots [56] (see Materials & Methods).For the strength and time scales of the Prony series, we assume that the springs are made from Stratasys Agilus30 material, a rubberlike material that is used in polyjet additive manufacturing and is commonly used as a prototypical example of a highly viscoelastic solid [45,57,58,49,59,60].The strength and time scales of the Prony series have been adapted from [45] S lo w F a s t  [55] with permission).(B) A laterally pre-confined viscoelastic mechanical metamaterial for shock absorption experiences a snap-through instability only when compressed slowly (Reproduced from [45] with permission from ASME). and can be found in the Materials & Methods.We employ w/h = 10 and K A /K B = 0.05.We then move the truss down with a constant rate and track the force-displacement response at various loading rates in Fig. 3E.
When going from small loading rates in yellow to intermediate loading rates in green and cyan, we observe a change from bistability (force goes through zero three times) to monostability (force crosses zero only once).Also, the reaction force P increases over the entire loading regime when increasing the loading rate.Both of these phenomena were also captured by the quasi-static analysis in Fig. 2B.However, when we further increase the loading rate, we again identify bistability, which is not captured by a simple quasistatic analysis.Relaxation from this state could induce another state of intermediate monostability instead, so the same structure could also capture pseudo-bistability.This shows that it is not self-evident to predict in what way viscoelasticity affects the pseudo-stability of a buckling structure.Moreover, while the quasi-static analysis in Fig. 2B suggests that the first limit point will shift to the right when increasing the equivalent loading rate, the viscoelastic analysis in Fig. 2E shows the opposite effect.This shows that it is not trivial to identify a priori whether viscoelasticity will increase or decrease the buckling deflection either.In both of these cases, increased viscoelasticity, as identified by a higher loading rate, has led to an increase in buckling force.
Yet, the opposite can also be true.Stein-Moldavo et al. recently showed that the buckling pressure of spherical shells could also drop due to viscoelastic effects [44].When a negative pressure was kept stable just below the elastic buckling limit, a viscoelastic shell would still buckle after relaxation.This is because the relaxation induces additional geometric imperfections, which in turn leads to a reduction in the buckling pressure.However, depending on the direction of the imperfections, the result could also act opposing and make the system more bistable [61].
In summary, the interplay between instabilities, imperfections and viscoelasticity all contribute to shaping how snap-through buckling is affected by viscoelasticity.
In turn, this knowledge of viscoelasticity-induced pseudo-stability can be used to understand how the behaviour of mechanical metamaterials is affected by viscoelasticity.First of all, Li et al. have shown that viscoelastic dissipation can delay buckling in viscoelastic metamaterials, which in turn can delay a shift from non-auxetic to auxetic in mechanical metamaterials [62].Moreover, Berwind et al. showed in Fig. 3A that pseudo-stability can occur in mechanical metamaterials just as it could for simpler visco-elastic snap-through instabilities [55].Pseudo-bistability can even be harnessed to program a morphing pattern of viscoelastic snap-through shells [63].On the other hand, Dykstra et al. showed in Fig. 3B that viscoelasticity can prevent a mechanical metamaterial featuring snap-through instabilities from snapping altogether under large strain rates [45].From a shock damping perspective, this means that a higher strain rate could actually induce poorer shock damping performance when no snap-through occurred.The strain at which this transition occurred was then tailored by shapechanging.This showed the very careful interrelation between geometry and dissipation in viscoelastic metamaterials.
In summary, delayed snap-through and modification of stability by viscoelasticity can occur because different parts of the structure or the metamaterial can experience different strain rates.As a result, these different parts experience dif-ferent levels of stiffening induced by viscoelasticity and this modifies the energy landscape.In the second section, we will show that suitable spatial patterning of the viscoelastic properties allows to further control how the energy landscape depends on strain rate.

Dissipation to Seed Dynamic Imperfections
Flexible metamaterials based on buckling and snapthrough instabilities exhibit complex energy potentialswith mountains and valleys.We argue in this section that viscoelastic dissipation can be used to introduce a new stable bifurcation state [64].This in turn is an efficient way to control how metamaterials navigate such landscapes, based on the loading rate (Fig. 4A) [47].As such, the deformation pathways and ensuing shape-changes and effective properties of the metamaterials can depend dramatically on whether the metamaterials are loaded fast or slow.
The control over such deformation pathways can be achieved by spatially patterning the viscoelasticity of the flexible building blocks.The simplest example is perhaps that of Euler beam buckling.In such a case, the patterning of the viscoelastic relaxation strength allows to induce an imperfection that depends on strain rate [48] (Fig. 4B).In practice, take a column split in two parts: the left (right) half of the column is made of a material with a small (large) viscoelastic relaxation strength.When the column is loaded along its axis, it will buckle to the left (right) when loaded slower (faster) than the typical relaxation timescale of the more viscoelastic material.This happens due to a shift in the neutral axis.When loaded slowly (fast), the viscoelastic material is much softer (stiffer), so the neutral axis shifts left (right) of the point of load application and the beam therefore buckles to the left (right).When sufficiently compressed, nonlinear deformation prevents snap-back, implying a permanently stable bifurcation.This idea can be applied to construct flexible structures that exhibit switching between auxetic and non-auxetic behaviour (Fig. 4C) or show an apparent negative viscoelasticity-where the apparent stiffness is lower for larger loading rate and larger for lower loading rates [48].
The section above shows that the direction of buckling of a beam or a hinge can be controlled by spatially patterning its viscoelastic properties.In this case there is a single deformation mode (the out-of-plane bending) and one controls its direction by inducing a strain-rate dependent imperfection.
In the following, we also discuss a different approach, which consists of using structures with more than one deformation mode and where the spatial patterning allows to selectively actuate the mode of choice based on loading rate.Architectures that have been used for such an approach are oligomodal metamaterials (Fig. 4D) [49] or origami with multiple folding branches (Fig. 4E) [47].The idea is the following: a metamaterial has a distinct number of global deformation modes (or folding branches for origami), such as the metamaterial visualised in Fig. 4D which has two modes.One mode only actuates viscoelastic hinges, all of which act soft during slow loading of the structure.The other mode also actuates elastic hinges, which are much softer than the viscoelastic hinges during fast loading of the structure.As a result, fast loading will energetically favour deforming the mode featuring elastic components.Where there are more modes available, it is then also possible to obtain additional deformation modes at intermediate loading rates (Fig. 4E) [47].
However, using only viscoelasticity is often not enough to navigate the energy landscape effectively and to allow for a switchable actuation of the mode of choice [49].In some cases, some modes tend to have a lower energy barrier and as such can be preferentially actuated regardless of the rate of actuation or of the mismatch of viscoelasticity.Imperfections-traditionally considered a hazard in most mechanical systems-tailored to the preferred mode, can be used to tweak the energy landscape and hence to alleviate this problem.In quasi-static metamaterial designs they have allowed branching into specific chosen bifurcation paths [57,65], while they can enable switching between modes when combined with viscoelasticity [49].This means that imperfections can also be programmed specifically to change the loading rate at which mode switching occurs.
In summary, we have discussed how spatial patterning of viscoelastic relaxation strength allows to dynamically tune buckling imperfections or relative mode stiffness.We have discussed a few examples where this approach has been used to create metamaterials with shape-changes that depend on loading rate.In the next section, we will show that stored energy combined with viscoelasticity can also induce extreme mechanical properties.

Extreme Mechanics due to Negative Stiffness
Traditionally, the mechanical properties of a composite material typically cannot surpass those of its constituents [66,67].For example, for the Young's modulus, upper and lower bounds for laminates are provided by the Voigt [68] and Reuss [69] bounds, while upper and lower bounds for isotropic composite materials are provided by the Hashin-Shtrikman bounds [70] (Fig. 5A).A similar compromise is present for damping.As a result, a combination of high damping and high stiffness seems to be mutually exclusive [67].
However, these bounds all assume that the constituent materials have positive definite elastic moduli, also known as positive stiffness [66,71,72,73,67].In other words, they assume that the constituents do no not contain stored energy.However, while a freestanding material made out of a single material must display positive stiffness, confined systems can store energy.For example, post-buckled systems can display negative stiffness by releasing mechanical energy [67].Viscoelasticity in turn stabilises this state of negative stiffness [74,75,76,71,77,72,67].In turn, the negative stiffness in these post-buckled systems can be harnessed to generate extreme damping [78,79,67].
Furthermore, confined material unit cells can also display negative stiffness.In a notable example, Jaglinski et al. showed in Fig. 5B  sition from positive to negative stiffness when it is cooled without allowing deformation [71].When the normalised temperature T N becomes negative, the Landau energy [80] as function of strain ε becomes non-positive definite, implying negative stiffness (Fig. 5B).While such a negative stiffness unit cell would be unstable in an unconfined configuration, Jaglinski et al. embedded such unit cells in a tin matrix, creating a mechanical metamaterial featuring negative stiffness inclusions [71].As a result, they were able to generate extreme damping in cyclic loading (Fig. 5C) [71].
In turn, it even allowed them to generate a viscoelastic modulus greater than diamond, the hardest natural material known to humankind (Fig. 5C) [71].
Other authors have shown similar results, both theoretically and experimentally: negative stiffness stabilised by viscoelasticity can generate extreme mechanical properties, including extreme mechanical damping, extreme viscoelastic modulus and extreme thermal expansion [74,75,76,71,77,72,67].This shows the great potential of using negative stiffness in viscoelastic metamaterials.GPa; and tin, 50 GPa (Reproduced from [71] with permission from AAAS).

Modelling Viscoelastic Metamaterials
Another issue that comes with designing nonlinear viscoelastic metamaterials, is that of computational effort.While modelling quasi-static or linear systems require relatively little computational effort, modelling nonlinear dissipative systems often requires computationally expensive nonlinear dynamics.In most cases, the response can fully and accurately be modelled using nonlinear dynamic finite element methods.For relatively smooth dynamic analysis, implicit finite element methods can be used.Implicit dynamic analyses are analyses where relatively large time steps can be taken which each have to converge individually within numerical accuracy to the exact solution.We have provided the codes for two of such cases on Zenodo [81], corresponding to the analyses of Fig. 2DE and Fig. 4BC, which have been solved using the viscoelastic solver in Abaqus (Simulia).For more complicated dynamic analyses, explicit analyses may be required.Explicit dynamic analyses are analyses, which use very small time steps that do not need to converge to the exact solution individually.Instead, such analyses may jump around the exact solution from one step to another.While implicit dynamic analyses can already be computationally expensive for a detailed mesh, explicit dynamics can be orders of magnitude more expensive when modelling long model time durations.Often, using such dynamic finite element solvers is not computationally viable for more complicated systems, parametric studies or optimisation.As such, reducing this computational effort is another key point to advance the field of viscoelastic shape-changing metamaterials.
Several strategies have already been identified.For one, many dissipative shape changes without inertial effects can be approximated with fully elastic quasi-static analyses, by scaling local time-dependent material stiffnesses with a fixed value.Such analyses can qualitatively match the key shapechange, which has been done among other for the results in Fig. 2AB [39] and Fig. 4D [49].However, these analyses only give a qualitative description and can fall short in obtaining the requested response [49].On the other hand, a simplified model employing full nonlinear dynamics, such as that of Fig. 2DE, can also offer qualitative agreement but may still fall short quantitatively [45].In such cases, the challenge remains in finding a minimum viable model without quantitative sacrifice.For systems that are overdamped or quasi-static, a viscoelastic finite element solver which neglects the effects of mass can prove viable (Fig. 7 [48]).Such solvers are available in both an implicit and explicit variant.For these overdamped or quasi-static analyses, neglecting the effects of mass can lead to multiple orders of magnitude improvement in computational effort, without quantitative sacrifice.In cases where mass is not negligible or where many unit cells have to be analysed simultaneously, another option would be to model mechanical metamaterials as a continuum.Statically, a conformal mapping approach has showed excellent quantitative agreement in modelling mechanism based metamaterials [82].More-over, Glaesener et al. have even shown that viscoelastic truss-based metamaterials can be modelled as time-dependent continua [83].In doing so, they managed to capture the full nonlinear dynamics of the metamaterials accurately.
These works have shown that various strategies can be taken to model viscoelastic metamaterials more efficiently.These strategies, and more efficient efficient strategies still, will allow us to predict and design a new generation of vis-coelastic metamaterials.

Outlook
We have reviewed here how dissipative effects can be used to control the extreme mechanics of metamaterials in the context of viscoelasticity.Further progress could be achieved in this direction by further improvements in additive manufacturing [84], material science [85] and computational design methods.In particular, better control of multimaterial printing of materials with strong contrast in their viscoelastic properties could help in taking the ideas presented in this article further-one could imagine more complex geometries or different ways to embed internal energy.Printing at the microscale could help in designing a wider range of shape-changing metamaterials with multiple functionalities.One could also harness viscosity to turn fluids subjected to acoustic vibrations into mechanical metamaterials [86].
One can also look beyond the case of viscoelasticity treated here to manage the energy landscape of adaptive and interactive materials and materials systems [87].For instance, Che et al. [46] showed that temperature can be used to expand the design space of the previously discussed pseudo-bistability in mechanical metamaterials [40,41,42].In particular, they showed that by changing the temperature, they could control the amount of time it took to snapback.This is because the viscoelastic properties of rubbers are highly temperature dependent.Further down the line, one could envision that recent ideas on odd viscoelasticity [88,89,90,91,92,93]-where the viscoelastic tensors are not symmetrical and which have so far required active elements-could be realised via suitable combinations of purely passive viscoelastic elements.
Alternatively, dissipation can be harnessed from different mechanisms than viscoelasticity.For instance, internal adhesive contacts can be used to generate extreme dissipation and recovery after loading of nanofoams [94,95].Active dissipation on the other hand can be used to control waves [96,88] or achieve robotic functionalities [89] in mechanical metamaterials.Otherwise, when combining elastomers with Low Melting Point Alloys, plasticity and heat be used to generate reversible plastic deformation [97].This generates a response similar to what can be obtained using shape memory materials.One can also think of a plethora of alternative ways to harness dissipation, including for instance piezoelectric [98], magnetic [99,100], photo-responsive [101] or superelastic [102] materials or fluidic devices [103].In an even broader sense, nearly any type of time-dependent mechanical properties can be exploited for 4D printing: 3D printing with a fourth dimension in time [104,105].4D printing can be used among others to generate morphing or reconfigurable structures, where viscoelastic creep and stress relaxation can play a vital role [105].
With so many possibilities to harness dissipation in the design of mechanical metamaterials, the next step is to go from generating new types of mechanical behaviour to new functionalities to actual applications.Obvious applications include using extreme vibration and shock damping due to negative stiffness in e.g.aerospace, automotive or high-tech machinery [67].However, viscoelastic metamaterials also show potential in synthetic soft tissues and implants [106].Dissipative materials can also be beneficial to enrich the design of machine-like devices.For example, the locomotion of soft robots can be influenced by embedding dissipative components in their bodies [89].For instance, the control of motion can be reduced to the control of flow rate in pneumatic soft actuators [107,103].For example, at larger scales, dissipative mechanisms in wearable devices and soft exosuits [108] can promote the design of such soft systems by reducing the duty of active elements and annexing memory.

Viscoelastic Material Models
The two most commonly used viscoelastic material models for viscoelastic solids are the Kelvin-Voigt model (Fig. 6A) and the Generalised Maxwell model (Fig. 6B) [109].
The Kelvin-Voigt model (Fig. 6A) consists of a spring, with a dashpot in parallel and is typically used to describe solids featuring fluid or frictional dissipation.The Kelvin-Voigt model is governed by the following equation [109], where σ is stress, E is the Young's modulus, η is the dashpot stiffness, ε is strain and t is time: The Generalised Maxwell model (Fig. 6B) consists of a spring, with multiple spring-dashpots in parallel and is typically used to describe material dissipation.It is also referred to as the Maxwell-Wiechert model [56].The spring-dashpots are each known as Maxwell elements, each with a time scale τ i = η i /E i .When only one Maxwell elements is used, the Generalised Maxwell model reduces to the Standard Linear Solid [109].
When a material based on the Generalised Maxwell model is subjected to a stress relaxation test, the time response of the Young's modulus, E(t), is described by: with E 0 the peak Young's modulus under instantaneous load, β n the dimensionless relaxation strength, τ n the timescale of the individual Maxwell elements [56], N the number of Maxwell elements considered and E e = E 0 1 − ∑ N n=1 β n the fully relaxed Young's modulus.Eq. 2 is also known as a Prony series.Whether the response at short time scales can

Viscoelastic Properties of Agilus30
To determine the visco-hyperelastic properties of Ag-ilus30, we performed stress relaxation experiments, which have largely been recorded previously in [45].
An Agilus 30 dogbone section with a length L = 50 mm, depth d = 5 mm and width w = 10 mm has been 3D printed using a Stratasys Objet500 Connex3 printer.The sample is stretched quickly at a strain rate of ε = 0.4 s −1 , using a uniaxial testing device (Instron 5943) to a strain ε = 20%, after which the force is allowed to relax for one hour.The data is measured with a frequency of 1000 Hz with t = 0 s defined at the point of highest load.
The parameters of Eqn. ( 2) are then fitted using a leastsquares fit to the test data, which was interpolated on a logarithmic time scale from 0.01 s to 1 hour.The relaxation test results are presented in Fig. 6C, including fits with N=1, N = 3 and N = 5, with corresponding material properties in Tab. 1.

Switchable Bi-beam
Here, we provide an example of a model of the strainrate dependent bi-beams described in Fig. 4BC [48].The code is available on Zenodo [81].
Bi-beams are made by attaching two similar flexible beams made from two different materials.The left beam (Fig. 7A) is made of a stiff elastic material.The right beam is highly viscoelastic and softer than the elastic one in the relaxed state but stiffer at high strain rates.Such a bi-beam will predictably buckle to either the left (Fig. 7B) or the right (Fig. 7B) side, depending on the rate of the applied strain.
Nonlinear buckling analysis was performed using the non-linear visco solver in Abaqus (Abaqus 2021.Simulia, USA) [48]: an implicit dynamic solver (see Section 5), which ignores inertial effects.In order to discretise the geometry of the bi-beam, we used 8-node biquadratic hybrid elements (CPE8H) while 12 elements were seeded through the width of the bi-beam.The bi-beam is clamped at its both ends while the top clamped nodes are moved longitudinally in our simulations.A linear elastic material model (E = 1.5 MPa) was used to define the elastic properties of the left beam and the viscoelastic properties of the right beam is according to the five term Prony series materials parameters that is proposed for Agilus in Table 1.The bi-beam in our model has a length and width of 80 mm and 20 mm respectively.
The code shared [81] can be used to model buckling of bi-beam (to the left) at a strain rate ε = 1×10 −3 1/s (Fig. 7B) and (to the right) at a strain rate ε = 1 1/s (Fig. 7C).It can also serve as a general example on how to model viscoelastic metamaterials.

Data and Code Availability
The data and codes that support the figures within this paper are publicly available on a Zenodo repository [81].

Fig. 1 :
Fig. 1: Viscoelastic Metamaterials.Viscoelasticity affects the dynamics of buckling and can host new responses and functionalities.

Fig. 2 :
Fig. 2: Pseudo-bistability controlled by viscoelasticity.(A)A single rigid truss supported at both ends by elastic springs (Reproduced from[39] with permission from Elsevier).Changing the ratio of the spring stiffnesses can be considered equivalent to increasing the loading rate in a viscoelastic analysis.The equilibrium paths in (B) (Reprinted from[39] with permission from Elsevier) show that both monostable and bistable configurations are possible.(C) A viscoelastic jumping popper shows pseudo-bistability (Reproduced from[42] under CC BY 4.0).(D) A viscoelastic version of the Von Mises Truss of (A), featuring Generalised Maxwell springs with 5 spring-dashpots (Modified from[39] with permission from Elsevier).(E) Force-Displacement curves due to moving at constant speed, with the colour indicating the total loading time in seconds.Both very high and very low speeds induce bistability, while intermediate speeds can generate monostability.

Fig. 3 :
Fig. 3: Delayed buckling by dissipation in mechanical metamaterials.(A) A viscoelastic unit cell (a) after compression (b), showing initially bistable behaviour (c) before snapping back after relaxation (d) (Copyright Wiley-VCH GmbH.Reproduced from[55] with permission).(B) A laterally pre-confined viscoelastic mechanical metamaterial for shock absorption experiences a snap-through instability only when compressed slowly (Reproduced from[45] with permission from ASME).

Fig. 4 :
Fig.4: Dissipation to seed dynamic imperfections.(A) Viscoelasticity and loading rate can be used to navigate an energy landscape with multiple valleys (Reproduced from from[47] under CC BY 4.0).(B) Elastic-viscoelastic bi-beams can switch buckling direction when changing loading rate (Reproduced from[48] under CC BY-NC 4.0) and (C) can allow switching between auxetic and non-auxetic when turned into a mechanical metamaterial (Reproduced from[48] under CC BY-NC 4.0).(D) Loading rate can also be used to switch modes altogether in metamaterials with a distinct number of global modes and both elastic and viscoelastic hinges (Reproduced from[49] under CC BY-NC-ND 4.0).(E) An origami sheet featuring a higher number of modes and both elastic and viscoelastic hinges can also experiences alternative modes at intermediate loading rates (Reproduced from from[47] under CC BY 4.0).

Fig. 5 :
Fig. 5: Negative Stiffness Inducing Extreme Dynamic Properties.(A) Young's modulus bounds for composite materials made out of constituents with positive definite elasticity.The Young's modulus for laminates is bound between the Voigt [68] (blue dashed line and inset) and Reuss [69] (red solid line and inset) bounds.The Young's modulus for isotropic materials is bound between the upper (green dotted line) and lower (cyan dash-dotted line) Hashin-Shtrikman [70] bounds.The two constituents have a stiffness ratio of 10 and Poisson's ratio ν = 0.3.E n is the composite Young's modulus normalised with the Young's modulus of the stiff constituent.f is the volume fraction of the stiff constituent.(B) Landau energy function of strain ε and normalised temperature T N = (αγ/β 2 )(T -T 1 )-0.25, with unit cells of BaTiO 3 in cubic and tetragonal phases.α,γ, and β are constants that depend on the material (Reproduced from[71] with permission from AAAS).(C) Young's modulus |E * | and viscoelastic damping tanδ of tin-BaTiO 3 composite metamaterial, showing extremely high modulus over a range of temperatures.δ is the phase angle between stress and strain.Young's moduli of the constituents are BaTiO 3 , 100 GPa; and tin, 50 GPa (Reproduced from[71] with permission from AAAS).

Fig. 6 :
Fig. 6: Viscoelastic material models.(A)Representation of the Kelvin-Voigt model of viscoelasticity.(B) Representation of the Generalised Maxwell model.(C) Relaxation test for a dogbone sample of Agilus 30 3D printed material.Time-dependent Young's modulus vs. relaxation time.The thick blue line indicates the test results, the thin dashed orange, green and red curves present fitted data with N = 1, N = 3 and N = 5 respectively.

Fig. 7 :
Fig. 7: Bi-beam model.(A) Bi-beam geometry.(B) Buckling to the left at low strain rates.(C) Buckling to the right at high strain rates.