Research ArticleMATERIALS SCIENCE

Strain rate–dependent mechanical metamaterials

See allHide authors and affiliations

Science Advances  17 Jun 2020:
Vol. 6, no. 25, eaba0616
DOI: 10.1126/sciadv.aba0616

Abstract

Mechanical metamaterials are usually designed to exhibit novel properties and functionalities that are rare or even unprecedented. What is common among most previous designs is the quasi-static nature of their mechanical behavior. Here, we introduce a previously unidentified class of strain rate-dependent mechanical metamaterials. The principal idea is to laterally attach two beams with very different levels of strain rate-dependencies to make them act as a single bi-beam. We use an analytical model and multiple computational models to explore the instability modes of such a bi-beam construct, demonstrating how different combinations of hyperelastic and viscoelastic properties of both beams, as well as purposefully introduced geometric imperfections, could be used to create robust and highly predictable strain rate-dependent behaviors of bi-beams. We then use the bi-beams to design and experimentally realize lattice structures with unique strain rate-dependent properties including switching between auxetic and conventional behaviors and negative viscoelasticity.

INTRODUCTION

The story of mechanical metamaterials started with researchers trying to achieve unusual properties through the rational design of architected materials. The quest that started with a relatively limited scope has, however, rapidly pullulated and now targets not only such well-established behaviors as auxeticity (1–3), negative thermal expansion (4, 5), snap-through instability (6, 7), and ultrahigh properties (8, 9) but also more novel functionalities including shape-morphing (10, 11), nonreciprocity in transferring motions (12), and even logic gate–like programmability (13, 14). In a way, the boundary between materials and machines is becoming more and more blurred, with some researchers even introducing concepts like “machine matter” (15, 16).

In this context, flexible metamaterials (2, 17–19) show the most promise. On the one hand, the large, reversible, and nonlinear deformations exhibited by these materials provide a fertile ground for cultivating unique functionalities. On the other hand, flexible mechanical metamaterials share multiple organic intersections with other current areas of research including soft robotics (9, 20, 21) and flexible electronics (22–24). They have, therefore, received the most attention during past few years.

Most previous studies have, however, focused on the quasi-static response of mechanical metamaterials. In the cases where the dynamic response of mechanical metamaterials is studied, the focus is usually on their acoustic and wave propagation properties [e.g., see (25–30)]. However, the viscoelastic and strain rate–dependent behaviors of mechanical metamaterials, particularly those exhibiting substantial flexibility, could present many novel opportunities. This is particularly important given the fact that most of the present designs of flexible mechanical metamaterials either are made from polymers that exhibit vastly strain rate–dependent behavior (31) or include other energy-dissipating mechanisms that often lead to strain rate dependency. For example, there is substantial friction and contact in origami-based flexible metamaterials (10, 24, 32) that causes energy dissipation and strain rate dependency. It is only natural that the strain rate–dependent behaviors of flexible mechanical metamaterials are taken into the account, or even better exploited to further expand the space of achievable properties and functionalities. However, very few [e.g., see (31, 33)] major attempts in that direction can be found in the literature.

Here, we used the concept of bi-beams to design, analyze, and fabricate metamaterials whose mechanical behavior either switches (e.g., from conventional to auxetic) depending on the applied strain rate or exhibits new types of viscoelastic behavior (e.g., negative viscoelasticity). Bi-beams are made by laterally attaching two beams made from two different materials of which one is both highly deformable and highly strain rate–dependent (i.e., visco-hyperelastic), while the other is highly deformable but largely strain rate–independent (i.e., hyperelastic) (Fig. 1A). Under quasi-static conditions, the stiffness of the hyperelastic beam is higher than that of the visco-hyperelastic beam. When subjected to high enough compressive forces, such a bi-beam construct will predictably buckle to either the left or the right side, depending on the rate of the applied force (strain) (Fig. 1, A and B). To study the behavior of a single bi-beam, let us assume that the hyperelastic layer is always positioned at the left side of a bi-beam that is clamped at its both ends (Fig. 1, A and B).

Fig. 1 Programming the buckling behavior of hyperelastic/visco-hyperelastic bi-beams.

(A) A schematic drawing of a bi-beam made of two identical beams with limited viscoelasticity (the hyperelastic beam) and substantial viscoelasticity (the visco-hyperelastic beam). Because of the higher stiffness of the hyperelastic beam, the bi-beam buckles to the left when loaded slowly in compression. Under high strain rates, the apparent stiffness of the visco-hyperelastic beam may exceed that of the hyperelastic beam, causing the bi-beam to buckle to the right. (B) Visualization of the stress distributions corresponding to the three potential states of sufficiently compressed bi-beams (i.e., left buckling, no buckling, or right buckling). (C) Maps showing the computationally predicted buckling direction for a bi-beam with an aspect ratio of 8 and different values of the mechanical properties of the hyperelastic and visco-hyperelastic beams. The results are presented for different values of the applied nondimensional strain rate, and the color code indicates the regions where right buckling occurs. (D) Evolution of the stresses and strains of a visco-hyperelastic material as a function of the applied nondimensional strain rate for four different sets of material properties. Analytical relationships are used for drawing the stress-strain curves and compared with the corresponding hyperelastic materials.

RESULTS AND DISCUSSION

Analytical and computational models of perfect bi-beams

Neo-Hookean models were used for representing the nonlinear elastic behaviors of the hyperelastic (Neo-Hookean constant, C10H) and visco-hyperelastic (long-term Neo-Hookean constant, C10V∞) beams. A single-term Prony series represented the viscous behavior of the visco-hyperelastic beam (dimensionless coefficient of the Prony series, g1; relaxation time, τ1). We define the nondimensional time, t*, as t* = t/τ1 and the nondimensional strain rate, ε′, as ε′=dε/dt*=τ1ε̇. There are three material parameters (i.e., C10V∞,C10H, and g1) the balance of which determines the threshold of the nondimensional strain rate above which the bi-beam buckles to the right. We developed an analytical model based on the Euler’s buckling theory for beams on elastic foundations and applied it to a linearized version of the governing equations of the beams (see Materials and Methods for the derivations). According to this model, the condition for buckling to the right is that the instantaneous elastic modulus of the visco-hyperelastic beam, EV, should exceed the elastic modulus of the hyperelastic beam, EH (i.e., EV > EH). Equation 17 (see Materials and Methods) can then be used to determine whether the bi-beam construct buckles to the left or to the right. The envelope of the possible range of parameters for which the bi-beam buckles to the right is obtained when the nondimensional strain rate is very large. In such a case, Eq. 17 can be simplified as C10V∞/C10H+g1>1 (Eq. 18).

The computationally predicted map of the buckling direction (for a beam with an aspect ratio of 8) is in agreement with the predictions of the analytical model and shows that when C10V∞/C10H and g1 are both very small, the bi-beam always buckles to the left regardless of the applied nondimensional strain rate (blank area, buckling to left; color codes, buckling to right; Fig. 1C). When C10V∞/C10H or g1 or the sum of both is large enough, the bi-beam buckles to the right for sufficiently large nondimensional strain rates (Fig. 1C). Regardless of how large the strain rate is, it is not possible to overcome the effects of low C10V∞/C10H or g1 values, as the boundaries of the right-buckling region remain largely unchanged despite multiple orders of magnitude of increase in the nondimensional strain rate (e.g., from 10−6 to 10) (Fig. 1C). The stress-strain curves calculated for the corresponding hyperelastic and visco-hyperelastic materials using uniaxial loading condition could capture the essence of the observations made regarding the effects of the different material parameters on the buckling behavior of the bi-beams (Fig. 1D). When C10V∞/C10H and g1 are both small (i.e., <0.8), the stresses generated in the visco-hyperelastic beam are far below those of the hyperelastic beam regardless of the value of the nondimensional strain rate (e.g., set I, Fig. 1D). For large enough values of C10V∞/C10H and g1 or when the sum of both is large enough (i.e., g1+C10V∞/C10H≳1), high enough nondimensional strain rates may cause the stresses generated in the visco-hyperelastic beam to saturate to (e.g., set II, Fig. 1D) or exceed (e.g., sets II and III, Fig. 1D) those of the hyperelastic beam. It is important to note that because buckling happens at relatively low strains, the direction of buckling can switch to the right side even when the visco-hyperelastic beam is softer than the hyperelastic beam at higher levels of strain (Fig. 1D, inset). A combination of high values of g1 and C10V∞/C10H (i.e., g1+C10V∞/C10H≫1) can result in the visco-hyperelastic beam being stiffer than the hyperelastic beam even for relatively low (i.e., 10−2) values of the applied nondimensional strain rate (e.g., set IV, Fig. 1D).

Effects of bi-beam geometry

To study the effects of bi-beam geometry on its buckling behavior, we built computational models corresponding to bi-beams with smaller and larger aspect ratios (Fig. 2A). For aspect ratios larger than 8, the buckling behavior of the bi-beams was found to be less predictable. Moreover, the higher modes of buckling appeared in some cases (depending on the balance between the material parameters and the nondimensional strain rate). According to the Euler’s buckling theory, as the length of the beams, l, increases, the critical strain decreases quadratically. Because the lateral expansion of the beams is related to the axial strain through the incompressibility constraint, a quadratically smaller value of the critical strain also means a quadratically smaller value of the lateral deformation, w. The force exerted by an elastic foundation, Fel, is proportional to the lateral deformation of the beam it supports (Fel = kw). A much smaller value of the lateral deformation, therefore, decreases the force exerted by the elastic foundation. Consequently, the relative importance of other minor forces (e.g., the slight forces caused by asymmetric meshes or small asymmetries in the numerical solutions) will increase, making the prediction of the buckling direction more challenging. In the theory of beams on elastic foundations, the parameter y=kl4/EI (where k is the elastic modulus of the elastic foundation and E is the elastic modulus of the beam itself) could be used to predict the occurrence of the higher modes of buckling (34). For sufficiently large values of the beam length, y increases enough to cause the appearance of the higher modes of buckling according to the theory of beams on elastic foundations (34). In addition to the beam length, the elastic moduli of both beams (i.e., k and E) are also important in determining whether the first or higher modes of buckling will appear. It is, however, important to realize that the theory of beams on elastic foundations does not take into account the global buckling of the bi-beam. The elastic modulus of the elastic foundations considered in the theory of beams on elastic foundation is not corresponding to the elastic modulus of the other beam when both beams buckle simultaneously. It is, therefore, unlikely that the higher modes of buckling predicted by the theory of beams on elastic foundation are observed experimentally, as they may be preceded by the global buckling of the bi-beam.

Fig. 2 Effects of aspect ratio and geometric imperfections on the buckling direction of bi-beams.

(A) Computationally predicted maps showing the buckling direction for bi-beams with four different aspect ratios. The color code indicates the regions where right buckling occurs as a consequence of applying the corresponding nondimensional strain rate. The gray regions are the regions with conflicting results, while the orange regions highlight the occurrence of some higher modes of buckling. (B) Sensitivity of the computational results to geometric imperfections (imp.). The geometric imperfections originating from the first buckling mode (calculated using a linear buckling analysis) are applied to the geometry of bi-beams either to the left side (negative values of the buckling mode) or to the right side (positive values of the buckling mode). The predictions of buckling to the right direction are presented for different magnitudes of the applied imperfection. (C) Modification of the bi-beams design to minimize their sensitivity to geometric imperfections.

Our computational models showed the diminished predictability of the buckling direction and the appearance of the higher modes of buckling for some higher values of the aspect ratio (i.e., 12 and 16, Fig. 2A). Regardless of whether global buckling or the higher modes of local buckling occur, we decided to discard aspect ratios above 8 for further analysis and the design of cellular structures. For a smaller value of the aspect ratio (i.e., 4), the sensitivity of the buckling direction to the applied nondimensional strain rate increased, meaning that the region within which right buckling occurs enlarges more prominently as a consequence of increasing the applied nondimensional strain rate (Fig. 2A). This is a positive development in terms of affording greater flexibility in the design of strain rate–dependent metamaterials. However, bi-beams with smaller aspect ratios show higher critical forces that may not be desirable, as we ultimately aim to use the bi-beams for designing cellular structures. If the forces required for the local buckling of the bi-beams are too high, the global buckling of the entire cellular solid starts to compete with the local buckling modes of individual bi-beams, rendering the buckling behavior of the cellular structure less predictable. An intermediate aspect ratio (e.g., 8 in the case of our simulations) balances the two abovementioned aspects to make the buckling behavior of the bi-beams as predictable as possible.

Imperfect bi-beams and purposefully introduced imperfections

Geometric imperfections and equivalent material imperfections are the other important parameters influencing the buckling direction of the bi-beams. In the buckling analyses presented so far, we assumed perfect beam geometries. To analyze the effects of imperfection size on the buckling direction predicted by our computational models, we performed a sensitivity analysis in which a gradually increasing percentage of the eigenvector corresponding to the first buckling mode was introduced as a geometric imperfection to the shape of the bi-beam. The imperfect geometries were then used to perform nonlinear buckling analyses. Only for very small values of the imperfection sizes (i.e., ≤0.001 mm or ≤0.02% of the width of the bi-beams), the buckling direction predicted by the perfect models remained fully valid (Fig. 2B). For an imperfection size of 0.01 mm (=0.2% of the width of the bi-beams), the buckling direction was already affected for some combinations of the material properties and nondimensional strain rates, while the entire switching effect was lost for an imperfection size as small as 0.1 mm (Fig. 2B). These results show the extremely high degree of sensitivity of the buckling direction to geometric (or equivalent material) imperfections. We used two strategies to mitigate the effects of imperfection dependency. As we aim to compare our computational results with experimental observations, we will use the dimensional version of the strain rate, ε̇, in the remainder of the paper.

The first strategy was to change the design of our bi-beams. When studying the evolution of the stresses induced in both hyperelastic and visco-hyperelastic beams, we found that the nonuniform distribution of the stresses caused by geometric imperfections results in stress concentrations at the interface of both beams (Fig. 2C). For example, when the direction of the geometric imperfection is to the left (e.g., in a design with C10V∞/C10H=0.2,g1=0.2,τ1=100 1/s), the maximum stress values are close to the interface in the stiffer beam (i.e., hyperelastic layer) and the direction of buckling is to the left (Fig. 2C). Switching the direction of the introduced imperfection results in higher values of stresses developing at the left side of the bi-beam, meaning that it buckles to the right. To minimize the sensitivity of the bi-beams to geometric imperfections, we modified the geometry of our bi-beams by cutting out square-shaped areas with side lengths equal to 25% of the width of the bi-beams at the either side of their clamped corners (Fig. 2C). This allowed both beams to relax the stresses induced at their interfaces and enhanced the overall uniformity of the stress distributions in both hyperelastic and visco-hyperelastic beams (Fig. 2C). As a consequence of this design modification, the bi-beam buckled to the left regardless of the direction of the applied imperfection even when the imperfection size was 0.1 mm (Fig. 2C).

The second strategy was the purposeful introduction of geometric imperfections into the shape of the bi-beams to predispose them such that they tend to buckle to the left (i.e., the side of the hyperelastic beam) under low strain rates. This predisposition ensures that the buckling direction of a bi-beam is robustly predictable under low strain rates, regardless of any random imperfections that may result from the applied fabrication process. However, the purposefully introduced imperfections change the shape of the regions within which switching to right buckling could take place (Fig. 3). There are two changes in the shape of the switching zones. First, the size of the switching zones decreases as the imperfection size increases from 0.1 to 0.5 mm (corresponding to 2 to 10% of the width of the bi-beam) (Fig. 3). A smaller size of the switching zone is caused by the fact that the elevated stiffness of the visco-hyperelastic beam should overcome not only the higher stiffness of the hyperelastic beam but also the forces created because of the presence of the purposefully introduced geometric imperfections. Second, the boundary between both zones of left and right buckling, which was initially diagonal, starts to rotate and becomes fully horizontal for the larger sizes of the geometric imperfection (e.g., a left imperfection size of 0.5 mm) (Fig. 3). This indicates the increased relative importance of the properties of the visco-hyperelastic beam in determining the switching behavior of the bi-beam as compared to those of the hyperelastic beam (Fig. 3). In addition to the effects of the imperfection size, it is important to note that higher values of τ1 minimize the size of the switching zone, because, at lower strain rates, the incidence of right buckling for high τ1 values is high.

Fig. 3 Computationally predicted boundaries that define the regions within which a switching in the buckling direction happens with strain rate for bi-beams with purposefully introduced geometric imperfections.

The switching is from left (at lower strain rate) to right (at higher strain rates). These regions are determined as the intersections of the regions corresponding to right buckling at higher strain rate with the regions defining left buckling for lower strain rates.

A combined experimental and computational study of imperfect bi-beams

To experimentally evaluate our predictions, we fabricated four bi-beams from a hyperelastic silicone rubber and a 3D (three-dimensional)–printed visco-hyperelastic rubber and introduced geometric imperfections to predispose the bi-beams to left buckling (Fig. 4A). The dimensions of the purposefully introduced imperfections were measured to be <5% of the width of the bi-beams (visible in Fig. 4A). For low rates of the applied deformation (e.g., 100 mm/min, strain rate = 1.7 1/s), all four specimens reliably buckled to the left (Fig. 4A). For a higher rate of the applied deformation (e.g., 500 mm/min, strain rate = 8.3 1/s), all four specimens reliably buckled to the right (Fig. 4A). The experiments were repeated more than 10 times, and the buckling direction was consistent between all experiments and all four beams. The effects of the applied strain rate are demonstrated in movie S1. The measured stiffness values of the bi-beams clearly increased with the applied strain rate (Fig. 4A).

Fig. 4 Strain rate–based programming of the buckling direction in realistic bi-beams.

(A) Stress-strain curves of actual bi-beams compressed at different loading rates. The warm colors refer to left buckling, while the cold colors stand for right buckling. The bi-beams were predisposed to left buckling. Switching in the direction of buckling (from left to right) happens by increasing the loading speed (from 100 to 500 mm/min). (B) Comparison between the buckling direction and the stress-strain curves predicted by our computational models and the experimental observations of the same quantities for two different aspect ratios of the bi-beams (i.e., 6 and 8) and for bi-beams with a variation in the thickness of their hyperelastic beam. The lengths of the beams with aspect ratios of 8 and 6 (short bi-beam) were, respectively, 120 and 90 mm (beam thickness = 21 mm for all types of beams). The length of the thick bi-beams was 120 mm, while the width of the hyperelastic beam was 12.5 mm, and the width of the visco-hyperelastic beam was 7.5 mm. Two types of computational models were used with either single-term or seven-term Prony series representing the viscous behavior. Photo credit: Shahram Janbaz, Delft University of Technology.

To assess the validity of our computational predictions, we compared the buckling directions predicted by our computational models, where single-term Prony series were used for capturing the viscoelastic behavior of the bi-beams, with those predicted by more complex computational models (i.e., with seven-term Prony series) and experimental observations (Fig. 4B). For all considered design variations (i.e., different aspect ratios and a variation in the geometrical design of bi-beams), the strain rate–dependent buckling directions predicted by both types of computational models agreed with each other and with experimental observations (Fig. 4B). Moreover, the stress-strain curves predicted by the different types of computational models largely overlapped and were in reasonable agreement with the experimentally measured stress-strain curves (Fig. 4B). These observations confirm that the results presented in Figs. 1 to 3, which were predicted using a single-term Prony series, are predictive of the actual buckling direction of the bi-beams. Despite the general agreement between the computational and experimental results, there were some differences between the stress-strain curves predicted by our numerical models and those obtained experimentally. It is important to note that the estimated and measured stress-strain curves and their corresponding stiffness values highly agreed with each other at the beginning of the loading (i.e., for small strains) (Fig. 4B). For large strains, however, the experiments consistently showed a softer response than what was predicted by our computational models (Fig. 4B). This softer response could be attributed to imperfections, such as air bubbles and imperfect curing and some omissions in the computational models, such as the effects of the adhesive used for attaching the two beams on the mechanical response of the materials from which the beams are made.

Application of bi-beams for the design of strain rate–dependent metamaterials

The reliable and highly predictable switching of the buckling direction opens up many opportunities for the design of strain rate–dependent mechanical metamaterials whose structural elements are bi-beams. We designed and created two such types of cellular structures with the following properties: (i) strain rate–dependent switching between auxetic and conventional (i.e., nonauxetic) behaviors and (ii) negative viscoelasticity. We also explored the possibility of applying bi-beams for the design of soft mechanisms.

To achieve strain rate–dependent switching between a positive (i.e., conventional behavior) and a negative (i.e., auxetic behavior) value of the Poisson’s ratio, we connected two parallel but mirrored bi-beams (with purposefully introduced geometric imperfections) to each other through two highly stiff connectors to create the basic unit cell Fig. 5A, i). The unit cell was then repeated in different directions to create a lattice structure (Fig. 5A, iii). We created nonlinear finite-element models of the lattice structures with three different sets of material properties (C10V∞/C10H=0.7,g1=0.7,τ1=0.01 1/s; C10V∞/C10H=0.9,g1=0.6,τ1=1 1/s; and C10V∞/C10H=0.5,g1=0.8,τ1=1 1/s) that were chosen to demonstrate some of the possible choices of the material properties that lead to a strain rate–dependent switching behavior (Fig. 5A, ii and iii). The evolution of the Poisson’s ratio, ν, with the applied strain was calculated for all the three combinations of the material properties and for some different values of the applied strain rate (Fig. 5A, ii). In all cases, the Poisson’s ratio of the lattice structure was calculated to be initially close to zero, which is expected given the near-orthogonal angles of the struts constituting the unit cell. As the strain increased, however, all three lattice structures exhibited either a positive or negative value of the Poisson’s ratio depending on the applied strain rate (Fig. 5A, ii). The switching from an auxetic to a conventional behavior with the applied strain rate was also observed in our experiments on lattice structures [Fig. 5A (iv) and movie S2]. In practice, a fivefold increase in the deformation rate (from 100 to 500 mm/s) was sufficient to trigger the desired switching from an auxetic behavior to a conventional one (Fig. 5A, iv). If an opposite switching behavior (i.e., from conventional to auxetic) with the applied strain rate is desired, the bi-beams that constitute the repeating unit cell simply need to be flipped. There are many potential applications for such a switching behavior among which protection against impact or other types of sudden movements is an important area. In particular, the Poisson’s ratio of soft wearable devices may be chosen such that, under small strains, they morph the contours of the human body [e.g., similar to (35)], while all of the constituting unit cells should turn auxetic for high enough strain rates (e.g., in the event of a patient’s fall) to increase the energy absorption capacity of the wearable device and protect the patient against low-energy bony fractures that are a major complication associated with osteoporosis (36).

Fig. 5 Application of bi-beams for designing metamaterials and mechanisms with novel properties and functionalities.

(A) (i) A near-orthogonal unit cell made from predisposed bi-beams was used to form a lattice structure whose Poisson’s ratio is highly dependent on the applied strain rate and could switch between an auxetic (i.e., negative Poisson’s ratio) and a conventional (i.e., a positive Poisson’s ratio) behavior. (ii) Three different sets of material properties were chosen to evaluate the evolution of the Poisson’s ratio with strain under different rates of the applied strain. (iii) Switching from an auxetic to a conventional behavior. Some unit cells from the middle part of the cellular lattice have been magnified for better visualization of the strain distributions. (iv) Experimental observation of strain rate–based switching from an auxetic to a conventional behavior. (B) A circular arrangement of bi-beams was used to transform axial compression to either clockwise or counterclockwise rotation, depending on the applied strain rate. Moreover, the stiffness of the rotational mechanism is highly dependent on the applied strain rate. Photo credit: Shahram Janbaz, Delft University of Technology.

A circular arrangement of the bi-beams was then used to transform axial (compressive) loads to either clockwise or counterclockwise rotation, depending on the applied strain rate (Fig. 5B and movie S3). Together with the direction of rotation, the stiffness of the mechanism changed as well (Fig. 5B), meaning that the characteristic curve of such a mechanism would be also dependent on the applied strain rate. Such possibilities in transforming axial deformation to rotation, in switching between clockwise and counterclockwise rotation depending on the applied strain rate, and in adjusting the characteristic curve of a soft mechanism could be instrumental in the development of dedicated powertrains for soft machines.

A third application of bi-beams is in creating mechanical metamaterials with negative viscoelastic behavior. We defined negative viscoelasticity as the property of a material that shows lower instantaneous stiffness under higher strain rates. Solid materials usually only show positive viscoelasticity, meaning that their stiffness increases with the strain rate. In some ways, negative viscoelasticity is similar to “thixotropy” (37) in non-Newtonian fluids where the viscosity and, thus, the stiffness of the viscoelastic fluid decrease upon the application of a high–strain rate load. Negative elasticity is, however, distinct from thixotropy, as it is defined for permanently solid materials and is not associated with large changes in the viscosity of the material. To achieve negative viscoelasticity, we changed the structural element of our cellular structure from one single bi-beam to two bi-beams that were positioned with a small distance, a, between them (Fig. 6A). The applied strain rate determines the buckling direction of each of the constituting bi-beams, creating three types of contact situations between them (Fig. 6A). The direction of the bi-beams and the material properties can be selected such that for very small strain rates (e.g., 10−4 1/s), both constituting bi-beams buckle toward each other, thereby giving rise to rapidly increasing contact forces that increase the stiffness of the cellular structure (Fig. 6A). For intermediate strain rates (e.g., 1 1/s), both bi-beams initially buckle toward each other followed by a deformation to the same direction after a second instability (Fig. 6A). For high strain rates (e.g., 103 1/s), the bi-beams buckle toward opposite directions and away from each other, avoiding any contacts and the associated increase in the stiffness (Fig. 6A). The stress-strain curves and instantaneous stiffness values calculated for all the three ranges of the applied strain rates showed that, for large enough strains (e.g., >0.15), the stiffness associated with a high strain rate is many times smaller than that of a low–strain rate case (Fig. 6A). When such double bi-beams were arranged in a cellular structure using highly stiff connecting elements, the same type of negative viscoelasticity was observed (Fig. 4B and movie S4). We also performed experiments with double bi-beams to explore their response to different rates of the applied strain (Fig. 6C). For a deformation rate of 500 mm/s, the bi-beams buckled away from each other, similar to what was predicted by the computational models (Fig. 6C). For a deformation rate of 100 mm/s, the expected stiffening behavior takes place where the bi-beams buckle toward each other (Fig. 6C). With an increasing level of the applied strain, a second instability mode is activated, where one of the bi-beams pushes the other bi-beam to straighten and buckle to the other side, ultimately resulting in the separation of the bi-beams (Fig. 6C). For high enough deformations (e.g., >7 mm), the instantaneous stiffness and the maximum force were much lower for a deformation rate of 500 mm/s as compared to 100 mm/s (Fig. 6C). Of course, the negative viscoelasticity described here is a post-buckling phenomenon. For infinitesimal strains, the stiffness of a double bi-beam is slightly higher in the case of higher strain rates, as a higher strain rate initially increases the stiffness of the visco-hyperelastic beam present in the construction of each of the bi-beams (Fig. 6C). Because such post-buckling behaviors are highly dependent on the contact between the bi-beams, geometric parameters such as the distance between both bi-beams, a, could be used to adjust the observed properties and their strain rate dependency to befit the specific application at hand.

Fig. 6 Demonstration of negative viscoelasticity in a unit cell made by parallel arrangements of the bi-beams.

(A) At lower strain rates, the parallel bi-beams buckle toward each other, resulting in contact forces that increase the stiffness of the cellular structure, once the threshold of bi-beam instability has been reached. At higher strain rates, the bi-beams buckle away from each other. The stress-strain curves and the evolution of the elastic modulus with strain are compared for different loading rates. (B) The arrangement of the unit cells in a lattice structure and the response of the lattice structure to different rates of the applied strain. (C) Experiments demonstrating the same behavior as predicted by our computational models. Initially (i.e., before the onset of instability), the unit cells exhibit a stiffer behavior under higher strain rates. This is reversed after the onset of instability beyond which point the stiffness corresponding to a higher strain rate is lower than the one corresponding to a lower strain rate. Photo credit: Shahram Janbaz, Delft University of Technology.

Practical consideration: Imperfect loading conditions and scalability

In the last step of the study, we examined a few aspects that are expected to be of importance for the practical applications of bi-beams. First, the real-world applications of bi-beams are likely to be associated with imperfections in the way the loads are applied even when the design of the metamaterials is based on fully orthogonal loading conditions. Moreover, the design of strain rate–dependent mechanical metamaterials or the devices that incorporate them may necessitate loading conditions that are not fully axial. We, therefore, examined the effects of two different types of deviations from the previously discussed loading conditions on the strain rate–dependent buckling behavior of bi-beams. In the first experiment, we applied 20-mm offsets to the left and to the right (Fig. 7A and movie S5). Compressive loads with different rates (corresponding to the same approximate strain rates as applied to the bi-beams tested in the previous experiments) were then applied to the bi-beams. The bi-beams continued to exhibit the desired strain rate–dependent behavior despite the applied offset (Fig. 7A and movie S5). In the second experiment, we applied ≈45° rotation to the bi-beams before applying the compressive loads (Fig. 7A and movie S5). Once more, the bi-beams continued to exhibit the desired strain rate–dependent behavior despite the applied rotation (Fig. 7A and movie S5). These experiments show the robustness of the strain rate–dependent buckling behavior of the bi-beams that remain unchanged despite such large deviations from the axial loading conditions and underscore the potential of the bi-beam concept for practical applications.

Fig. 7 Robustness of the strain rate–dependent behavior and the scalability of the bi-beams.

(A) The bi-beams (thick bi-beam, aspect ratio of 6, similar to those used in Fig. 4B) exhibited the same strain rate–dependent buckling behavior despite imperfect loading conditions, including 20-mm offsets to the left and right and ≈45° rotation in the bi-beams. (B) Two miniaturized bi-beams (left: length = 22 mm, aspect ratio = 3.7, and thickness = 10 mm; right: length = 32 mm, aspect ratio = 5.3, and thickness = 10 mm) were manufactured and tested under low and high loading speeds that were scaled with the length of the bi-beams so as to maintain the same approximate strain rates as was used in the previous experiments (i.e., Figs. 4 to 6). Scaling down the length of the bi-beams by up to ≈5.3 times did not change the strain rate–dependent buckling behavior of the bi-beams. Photo credit: Shahram Janbaz, Delft University of Technology.

The other important aspect is the size of the bi-beams and how that may affect their buckling behavior. To examine the effects of specimen size, we manufactured two smaller versions of bi-beams (i.e., 5.5 and 3.8 times smaller) with different aspect ratios (i.e., 3.7 and 5.3) (Fig. 7B and movie S6). Both smaller versions of the bi-beams exhibited the expected strain rate dependency in their buckling behavior (Fig. 7B and movie S6). This confirms the scalability of the presented approach at least in the range of the sizes considered here. However, the concept of bi-beams and their strain rate–dependent behavior could be of interest for a wide range of practical applications. While structural application in which the length scales of the bi-beams is larger and in the range of those demonstrated here are relatively straightforward to realize, many other applications particularly in microelectromechanical systems, (steerable) medical instruments, and microrobots will require manufacturing of (a large number of) bi-beams at much smaller length scales. While the mechanical principles that are used in the design of bi-beams and the demonstrated strain rate–dependent behavior are highly scalable, the realization of such bi-beam constructs at smaller scales could pose some challenges. In terms of manufacturing, the microscale fabrication techniques required for such applications are already available. In principle, (multi-material) additive manufacturing techniques such as material jetting (9) and direct laser writing (38) can be used to manufacture similar constructs at the microscale and submicrometer scale, respectively. The main challenge, however, is in finding photopolymerizable materials that can be processed using the currently available techniques and exhibit the differential viscoelastic properties required for creating bi-beams. Moreover, the adhesion between both types of polymers (i.e., visco-hyperelastic and hyperelastic) should be very strong so as to support the transfer of forces between them. These are very specific requirements that, according to our preliminary screening of the currently available materials, are difficult to simultaneously satisfy. Alternatively, the differential strain-rate response of the beams of a bi-beam construct could be created in other ways of which inertia forces are a prime example. The use of inertia forces is much easier to realize in practice, and such types of bi-beams would be highly scalable. However, there is a limit (i.e., the difference between the range of the density of processable materials and voids) that determines the maximum amount of the differential behavior achievable using inertia forces. The application of additive manufacturing techniques would also mean that it is, at least in principle, relatively easy to fabricate 3D version of bi-beam arrays. This is likely to constitute a worthy avenue for further research, as the mechanical behavior of such bi-beam arrays is expected to be quite rich.

CONCLUSIONS

In summary, we demonstrated how strain rate–dependent mechanical metamaterials can be realized using bi-beams, which consist of two laterally attached beams of which one is hyperelastic and the other is visco-hyperelastic. Our computational models show the importance of geometrical design including the shape of the bi-beams (e.g., cut-out patterns), aspect ratio, and purposefully introduced imperfections in the reliability of the strain rate–dependent control of their buckling direction. We also showed how cellular structures and soft mechanisms exhibiting novel strain rate–dependent behaviors could be realized using (double) bi-beams as their structural elements. The strain rate–dependent switching between auxetic and conventional behaviors, negative viscoelasticity, and strain rate–dependent transformation of axial deformation to (counter) clockwise rotation are expected to find potential applications in various areas, including soft robotics, the development of wearable medical devices, and flexible electronics.

MATERIALS AND METHODS

Nonlinear buckling analysis was performed using an implicit nonlinear solver (Abaqus 2016. Simulia, USA) to study the effects of material parameters on the buckling behavior of bi-beams. Quadrilateral elements (type CPE4H) were used for the discretization of the bi-beam geometries (length, 40 mm; width, 5 mm; thickness, 25 mm). On the basis of the results of a mesh convergence study, at least 12 elements were seeded through the width of the bi-beams. The aspect ratio of the beams was between 4 and 16 (varied in steps of 4). The bi-beams were clamped at both their ends, and a reference point was created to apply the loads to the top side of the bi-beams. A series of simulations were performed for applied strain rates ranging between ε̇=10−4 and ε̇=103 1/s. The Neo-Hookean material model was used to define the nonlinear elastic behavior of both beams. Prony series were used to describe the viscous behavior of the visco-hyperelastic beam. The strain rate–dependent constant that defines the Neo-Hookean material model can be expressed asC10R(t)=C100(1−∑i=1ngi(1−e−tτi))(1)in which the instantaneous modulus is determined fromC10∞=C100(1−∑i=1ngi)(2)where n, t, τi, gi, C100, and C10∞ are the number of the series terms, real time, relaxation time, dimensionless coefficients of the Prony series, instantaneous modulus, and long-term modulus, respectively. The strain energy functions of the visco-hyperelastic, WV, and hyperelastic, WH, beams are, respectively, given as (39)WV=C10V0(I1−3)(1−∑i=1ngi(1−e−tτi))+pV(I3−1)WH=C10H(I1−3)+pH(I3−1)(3)where I1 and I3 are the first and third invariants of the left Cauchy-Green deformation tensor, B, respectively. The terms pV and pH represent the arbitrary hydrostatic pressures that are defined such that the imposed incompressibility conditions are satisfied. The Cauchy stress, σ, could be derived from the strain energy function asσ=2I3B∂W∂B(4)

Given that ∂I1∂B=I and ∂I3∂B=I3B−1, the Cauchy stresses of the visco-hyperelastic and hyperelastic beams can be calculated asσV=2I3{C10V0(1−∑i=1ngi(1−e−tτi))B+pVI3I}σH=2I3{C10HB+pHI3I}(5)where I is the second-order identity tensor. The eigenvalues of the left and right Cauchy-Green deformation tensors are called the principal stretches, λi (i = 1, …,3). By imposing the incompressibility condition (i.e., I3 = λ1λ2λ3 = 1) and considering uniaxial deformation in the direction 1, the left Cauchy-Green deformation tensor is obtained asB=[λ2000λ−1000λ−1](6)

The incompressibility condition and uniaxial loading conditions (i.e., σ22 = σ33 = 0) could then be used to determine the arbitrary pressures pV and pH. The stress in the direction 1 can be determined by combining Eqs. 5 and 6 asσV,11=2C10V0(1−∑i=1ngi(1−e−tτi))(λ2−1λ)σH,11=2C10H(λ2−1λ)(7)

Assuming a single-term Prony series, the instantaneous elastic modulus in the direction 1 is given byEV,11=∂σV,11∂ε=∂σV,11∂λ=2C10V0(1−g1(1−e−tτ1))(2λ+1λ2)EH,11=∂σH,11∂ε=∂σH,11∂λ=2C10H(2λ+1λ2)(8)

For the small axial deformations preceding buckling, we have λ ≅ 1. The elastic moduli of the visco-hyperelastic and hyperelastic beams could then be written asEV=6C10V0(1−g1(1−e−tτ1))EH=6C10H(9)

To study the buckling behavior of the bi-beams, each of the beams (i.e., visco-hyperelastic and hyperelastic beams) can be modeled as a beam on an elastic foundation with each beam acting as the elastic foundation of the other beam. The buckling load for a clamped-clamped beam on elastic foundation, Pcr, is higher than the Euler’s buckling load for an equivalent unsupported beam, Pe, and can be obtained from the solution to the following implicit equation (34)sin12π2x+2ysin12π2x−2y=∓π2x+2yπ2x−2y(10)where x=PcrPe, y=kl4EI, k is the modulus of the elastic foundation, l is the length of the Euler beam, E is the elastic modulus of the material from which the beam itself is made, and I is the second moment of area of the cross section of the beam. This implicit equation can be numerically solved [e.g., see (34)]. The resulting solutions show relatively small deviations from the line y=π22(x−4) (or x=2yπ2+4). For the sake of the approximate analysis performed here, we use this line to estimate the buckling loads of the beams. In a bi-beam construct, both beams undergo the same strain but different forces. Therefore, to determine which beam buckles first, we need to do the analysis in terms of the critical strains εcr,H and εcr,Vεcr,H=Pcr,HAHEH=xPe,HAHEH=Pe,HAHEH(2yHπ2+4)εcr,V=Pcr,VAVEV=xPe,VAVEV=Pe,VAVEV(2yVπ2+4)(11)where A stands for the cross-sectional areas of the beams and subscripts H and V refer to the hyperelastic and visco-hyperelastic beams. Assuming that the visco-hyperelastic and hyperelastic beams are geometrically identical, we havePe,HAHEH=Pe,VAVEV(12)

The condition for buckling to the right side (i.e., the side of the hyperelastic beam) can then be simplified from εcr,H> εcr,V toyH>yV(13)

Because the beams are geometrically identical and given that y=kl4EI, Eq. 13 can be rewritten asEVl4EHI>EHl4EVI(14)or simplyEV>EH(15)

Combining Eqs. 9 and 15 gives the condition for buckling to the right asC10V0C10H>11−g1(1−e−tτ1)(16)

Using Eq. 2, Eq. 16 could be rewritten asC10V∞C10H>1−g11−g1(1−e−tτ1)(17)

For extremely high strain rates (i.e., t/τ1 → 0), the condition for buckling to the right can be simplified asC10V∞C10H+g1>1(18)

One of the limitations of the beams on elastic foundation theory as applied here is that it neglects the global buckling of the beams. When both beams buckle simultaneously, the modulus of the elastic foundation, k, is different from the elastic modulus of the other beam. One way to estimate the elastic modulus of the elastic foundation would be to impose the requirement that the critical load predicted by the beams on the elastic foundation theory equals that of the Euler’s theory for a bi-beam made from two identical beams (i.e., with identical geometries and material properties). Such an approach would suggest that the modulus of the elastic foundation is a function of not only the elastic modulus of the beams but also their geometry, particularly their slenderness. This additional geometry-dependent coefficient will, however, not change the main requirement for buckling to the right (i.e., EV > EH).

In our computational analyses, we varied the ratio of the Neo-Hookean coefficient of the visco-hyperelastic layer (C10V∞) to the hyperelastic layer (C10H) between 0.1 and 1. To be able to study the effects of the parameters of Prony series, we initially used only a single term. The predictions of a single-term model were later compared with those of a seven-term model. The dimensionless coefficient (g1) was varied in the range of 0.1 to 0.99, and the relaxation time (τ1) was 0.01, 0.1, 1, 10, or 100 s. A MATLAB (MathWorks, USA) code that computed the instantaneous stiffness of the visco-hyperelastic material as a function of the applied strain rate was used to explore the reason for the switching behavior in the buckling direction.

To evaluate the effects of geometric imperfections on the buckling direction of bi-beams, we used an eigenvalue buckling analysis to predict the first mode of buckling for a single-material beam. Having added positive and negative values of the buckling mode as a geometric imperfection to the geometry of our bi-beams, we performed nonlinear buckling analyses. In the context of these analyses, the imperfection size refers to the maximum size of the offset from the centerline of the bi-beam geometry presented (dimension, mm).

To manufacture the bi-beams used in the experiments, visco-hyperelastic beams were separately fabricated using an advanced multi-jet printer (Stratasys Objet350 Connex3, USA) by using a ultraviolet-curable polymer (Agilus, Stratasys, USA). The surface of the visco-hyperelastic beams was then fully covered with a thin layer of a rubber silicone adhesive (Sil-Poxy, Smooth-On) to ensure good adhesion between Agilus layer and silicone rubber. We used a silicone rubber (Dublisil 30, Dreve Dentamid GmbH, Germany) to mold the hyperelastic part of the bi-beams. To fabricate the other half of the bi-beams, we printed PLA (poly lactic acid) molds using a fused deposition modeling (FDM) 3D printer (Ultimaker 3, the Netherlands). We then positioned the visco-hyperelastic beam in the PLA molds and used the same rubber silicone adhesive to ensure that they remained attached to the walls of the molds during the curing process of the silicone rubber used for creating the hyperelastic beam. Subsequently, we filled the empty parts of the molds with silicone rubber. The bi-beams were kept in the molds for 8 hours at room temperature after which they were removed from the molds. We characterized the Neo-Hookean parameter, C10H, of the hyperelastic rubber based on the compression tests of rubber disks (D = 28.5 mm and h = 12.5 mm) under quasi-static conditions, while the Neo-Hookean parameter, C10V∞, and the Prony series parameters (τi and gi) were obtained on the basis of the relaxation tests performed for the visco-hyperelastic material using a previously described methodology (see the Supplementary Materials) (40). A Lloyd mechanical testing machine (LR5K) equipped with a 5-kN load cell was used for the compression and relaxation tests of the materials. In the case of the pretwisted bi-beams (Fig. 7), we used a plexiglass piece to maintain the pretwisting of the specimens throughout the duration of the mechanical tests.

To ensure well-defined clamping conditions during compression tests, the bi-beams were placed in a compression test setup, which was 3D-printed from PLA using the same FDM printer. The speed of the crosshead of the testing machine was allowed to stabilize before compressing the bi-beams by introducing a gap (10 mm) between the moving crosshead and the test setup. Before placing the bi-beams in the test setup, their end parts were cut in such a way that a small geometric imperfection toward the hyperelastic side was materialized once they were clamped in the test setup. The lattice structures and parallel bi-beams were fabricated by printing the additional parts from PLA using the same FDM 3D printer (Figs. 5 and 6). The bi-beams and the cellular designs were compressed using the same Lloyd test bench equipped with a 100-N load cell. A digital camera (Sony A7R with a Sony FE 90 mm f/2.8 macro OSS lens) was used to record the deformation patterns of the bi-beams and the cellular designs during the compression tests.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/25/eaba0616/DC1

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.

REFERENCES AND NOTES

  1. ↵
    1. R. Lakes
    , Foam structures with a negative Poisson’s ratio. Science 235, 1038–1040 (1987).
    OpenUrlAbstract/FREE Full Text
  2. ↵
    1. S. Babaee,
    2. J. Shim,
    3. J. C. Weaver,
    4. E. R. Chen,
    5. N. Patel,
    6. K. Bertoldi
    , 3D soft metamaterials with negative Poisson’s ratio. Adv. Mater. 25, 5044–5049 (2013).
    OpenUrl
  3. ↵
    1. K. Bertoldi,
    2. P. M. Reis,
    3. S. Willshaw,
    4. T. Mullin
    , Negative Poisson’s ratio behavior induced by an elastic instability. Adv. Mater. 22, 361–366 (2010).
    OpenUrl
  4. ↵
    1. Q. Wang,
    2. J. A. Jackson,
    3. Q. Ge,
    4. J. B. Hopkins,
    5. C. M. Spadaccini,
    6. N. X. Fang
    , Lightweight mechanical metamaterials with tunable negative thermal expansion. Phys. Rev. Lett. 117, 175901 (2016).
    OpenUrl
  5. ↵
    1. J. Qu,
    2. M. Kadic,
    3. A. Naber,
    4. M. Wegener
    , Micro-structured two-component 3D metamaterials with negative thermal-expansion coefficient from positive constituents. Sci. Rep. 7, 40643 (2017).
    OpenUrl
  6. ↵
    1. A. Rafsanjani,
    2. A. Akbarzadeh,
    3. D. Pasini
    , Snapping mechanical metamaterials under tension. Adv. Mater. 27, 5931–5935 (2015).
    OpenUrl
  7. ↵
    1. A. Rafsanjani,
    2. L. Jin,
    3. B. Deng,
    4. K. Bertoldi
    , Propagation of pop ups in kirigami shells. Proc. Natl. Acad. Sci. U.S.A. 116, 8200–8205 (2019).
    OpenUrlAbstract/FREE Full Text
  8. ↵
    1. X. Zheng,
    2. H. Lee,
    3. T. H. Weisgraber,
    4. M. Shusteff,
    5. J. DeOtte,
    6. E. B. Duoss,
    7. J. D. Kuntz,
    8. M. M. Biener,
    9. Q. Ge,
    10. J. A. Jackson,
    11. S. O. Kucheyev,
    12. N. X. Fang,
    13. C. M. Spadaccini
    , Ultralight, ultrastiff mechanical metamaterials. Science 344, 1373–1377 (2014).
    OpenUrlAbstract/FREE Full Text
  9. ↵
    1. S. Janbaz,
    2. F. S. L. Bobbert,
    3. M. J. Mirzaali,
    4. A. A. Zadpoor
    , Ultra-programmable buckling-driven soft cellular mechanisms. Mater. Horiz. 6, 1138–1147 (2019).
    OpenUrl
  10. ↵
    1. J. T. Overvelde,
    2. J. C. Weaver,
    3. C. Hoberman,
    4. K. Bertoldi
    , Rational design of reconfigurable prismatic architected materials. Nature 541, 347–352 (2017).
    OpenUrlCrossRefPubMed
  11. ↵
    1. C. Yang,
    2. M. Boorugu,
    3. A. Dopp,
    4. J. Ren,
    5. R. Martin,
    6. D. Han,
    7. W. Choi,
    8. H. Lee
    , 4D printing reconfigurable, deployable and mechanically tunable metamaterials. Mater. Horiz. 6, 1244–1250 (2019).
    OpenUrl
  12. ↵
    1. C. Coulais,
    2. D. Sounas,
    3. A. Alu
    , Static non-reciprocity in mechanical metamaterials. Nature 542, 461–464 (2017).
    OpenUrlCrossRefPubMed
  13. ↵
    1. D. J. Preston,
    2. P. Rothemund,
    3. H. J. Jiang,
    4. M. P. Nemitz,
    5. J. Rawson,
    6. Z. Suo,
    7. G. M. Whitesides
    , Digital logic for soft devices. Proc. Natl. Acad. Sci. U.S.A. 116, 7750–7759 (2019).
    OpenUrlAbstract/FREE Full Text
  14. ↵
    1. B. Treml,
    2. A. Gillman,
    3. P. Buskohl,
    4. R. Vaia
    , Origami mechanologic. Proc. Natl. Acad. Sci. U.S.A. 115, 6916–6921 (2018).
    OpenUrlAbstract/FREE Full Text
  15. ↵
    1. C. Coulais,
    2. E. Teomy,
    3. K. de Reus,
    4. Y. Shokef,
    5. M. van Hecke
    , Combinatorial design of textured mechanical metamaterials. Nature 535, 529–532 (2016).
    OpenUrlCrossRefPubMed
  16. ↵
    1. D. Yang,
    2. B. Mosadegh,
    3. A. Ainla,
    4. B. Lee,
    5. F. Khashai,
    6. Z. Suo,
    7. K. Bertoldi,
    8. G. M. Whitesides
    , Buckling of elastomeric beams enables actuation of soft machines. Adv. Mater. 27, 6323–6327 (2015).
    OpenUrl
  17. ↵
    1. K. Bertoldi,
    2. V. Vitelli,
    3. J. Christensen,
    4. M. van Hecke
    , Flexible mechanical metamaterials. Nat. Rev. Mater. 2, 17066 (2017).
    OpenUrl
    1. S. Janbaz,
    2. M. McGuinness,
    3. A. A. Zadpoor
    , Multimaterial control of instability in soft mechanical metamaterials. Phys. Rev. Appl. 9, 064013 (2018).
    OpenUrl
  19. ↵
    1. T. Frenzel,
    2. M. Kadic,
    3. M. Wegener
    , Three-dimensional mechanical metamaterials with a twist. Science 358, 1072–1074 (2017).
    OpenUrlAbstract/FREE Full Text
  20. ↵
    1. A. Rafsanjani,
    2. K. Bertoldi,
    3. A. R. Studart
    , Programming soft robots with flexible mechanical metamaterials. Sci. Robot. 4, eaav7874 (2019).
    OpenUrl
  21. ↵
    1. A. Rafsanjani,
    2. Y. Zhang,
    3. B. Liu,
    4. S. M. Rubinstein,
    5. K. Bertoldi
    , Kirigami skins make a simple soft actuator crawl. Sci. Robot. 3, eaar7555 (2018).
    OpenUrl
  22. ↵
    1. Y. Jiang,
    2. X. Liu,
    3. N. Matsuhisa,
    4. D. Qi,
    5. W. R. Leow,
    6. H. Yang,
    7. J. Yu,
    8. G. Chen,
    9. Y. Liu,
    10. C. Wan,
    11. Z. Liu,
    12. X. Chen
    , Auxetic mechanical metamaterials to enhance sensitivity of stretchable strain sensors. Adv. Mater. 30, e1706589 (2018).
    OpenUrl
    1. J. A. Rogers,
    2. T. Someya,
    3. Y. Huang
    , Materials and mechanics for stretchable electronics. Science 327, 1603–1607 (2010).
    OpenUrlAbstract/FREE Full Text
  24. ↵
    1. H. Zhao,
    2. K. Li,
    3. M. Han,
    4. F. Zhu,
    5. A. Vázquez-Guardado,
    6. P. Guo,
    7. Z. Xie,
    8. Y. Park,
    9. L. Chen,
    10. X. Wang,
    11. H. Luan,
    12. Y. Yang,
    13. H. Wang,
    14. C. Liang,
    15. Y. Xue,
    16. R. D. Schaller,
    17. D. Chanda,
    18. Y. Huang,
    19. Y. Zhang,
    20. J. A. Rogers
    , Buckling and twisting of advanced materials into morphable 3D mesostructures. Proc. Natl. Acad. Sci. U.S.A. 116, 13239–13248 (2019).
    OpenUrlAbstract/FREE Full Text
  25. ↵
    1. S. Babaee,
    2. N. Viard,
    3. P. Wang,
    4. N. X. Fang,
    5. K. Bertoldi
    , Harnessing deformation to switch on and off the propagation of sound. Adv. Mater. 28, 1631–1635 (2016).
    OpenUrl
    1. P. Celli,
    2. B. Yousefzadeh,
    3. C. Daraio,
    4. S. Gonella
    , Exploring heterogeneity and disorder in tunable elastic metamaterials. J. Acoust. Soc. Am. 143, 1917–1917 (2018).
    OpenUrl
    1. P. Celli,
    2. B. Yousefzadeh,
    3. C. Daraio,
    4. S. Gonella
    , Bandgap widening by disorder in rainbow metamaterials. Appl. Phys. Lett. 114, 091903 (2019).
    OpenUrl
    1. S. A. Cummer,
    2. J. Christensen,
    3. A. Alù
    , Controlling sound with acoustic metamaterials. Nat. Rev. Mater. 1, 16001 (2016).
    OpenUrl
    1. T. Frenzel,
    2. J. Köpfler,
    3. E. Jung,
    4. M. Kadic,
    5. M. Wegener
    , Ultrasound experiments on acoustical activity in chiral mechanical metamaterials. Nat. Commun. 10, 3384 (2019).
    OpenUrl
  30. ↵
    1. M. Kadic,
    2. G. W. Milton,
    3. M. van Hecke,
    4. M. Wegener
    , 3D metamaterials. Nat. Rev. Phys. 1, 198–210 (2019).
    OpenUrl
  31. ↵
    1. D. M. J. Dykstra,
    2. J. Busink,
    3. B. Ennis,
    4. C. Coulais
    , Viscoelastic snapping metamaterials. J. Appl. Mech. 86, 111012 (2019).
    OpenUrl
  32. ↵
    1. S. Janbaz,
    2. N. Noordzij,
    3. D. S. Widyaratih,
    4. C. W. Hagen,
    5. L. E. Fratila-Apachitei,
    6. A. A. Zadpoor
    , Origami lattices with free-form surface ornaments. Sci. Adv. 3, eaao1595 (2017).
    OpenUrlFREE Full Text
  33. ↵
    1. K. Che,
    2. M. Rouleau,
    3. J. Meaud
    , Temperature-tunable time-dependent snapping of viscoelastic metastructures with snap-through instabilities. Extreme Mech. Lett. 32, 100528 (2019).
    OpenUrl
  34. ↵
    M. Hetényi, Beams on Elastic Foundation: Theory with Applications in the Fields of Civil and Mechanical Engineering (University of Michigan, 1971).
  35. ↵
    1. G. P. T. Choi,
    2. L. H. Dudte,
    3. L. Mahadevan
    , Programming shape using kirigami tessellations. Nat. Mater. 18, 999–1004 (2019).
    OpenUrl
  36. ↵
    1. N. Papadimitriou,
    2. K. K. Tsilidis,
    3. P. Orfanos,
    4. V. Benetou,
    5. E. E. Ntzani,
    6. I. Soerjomataram,
    7. A. Künn-Nelen,
    8. U. Pettersson-Kymmer,
    9. S. Eriksson,
    10. H. Brenner,
    11. B. Schöttker,
    12. K.-U. Saum,
    13. B. Holleczek,
    14. F. D. Grodstein,
    15. D. Feskanich,
    16. N. Orsini,
    17. A. Wolk,
    18. A. Bellavia,
    19. T. Wilsgaard,
    20. L. Jørgensen,
    21. P. Boffetta,
    22. D. Trichopoulos,
    23. A. Trichopoulou
    , Burden of hip fracture using disability-adjusted life-years: A pooled analysis of prospective cohorts in the CHANCES consortium. Lancet Public Health 2, E239–E246 (2017).
    OpenUrl
  37. ↵
    1. J. Mewis,
    2. N. J. Wagner
    , Thixotropy. Adv. Colloid Interface Sci. 147-148, 214–227 (2009).
    OpenUrlCrossRefPubMed
  38. ↵
    1. T. Bückmann,
    2. N. Stenger,
    3. M. Kadic,
    4. J. Kaschke,
    5. A. Frölich,
    6. T. Kennerknecht,
    7. C. Eberl,
    8. M. Thiel,
    9. M. Wegener
    , Tailored 3D mechanical metamaterials made by dip-in direct-laser-writing optical lithography. Adv. Mater. 24, 2710–2714 (2012).
    OpenUrl
  39. ↵
    1. A. Ghorbanoghli,
    2. K. Narooei
    , A new hyper-viscoelastic model for investigating rate dependent mechanical behavior of dual cross link self-healing hydrogel. Int. J. Mech. Sci. 159, 278–286 (2019).
    OpenUrl
  40. ↵
    1. K. Narooei,
    2. M. Arman
    , Generalization of exponential based hyperelastic to hyper-viscoelastic model for investigation of mechanical behavior of rate dependent materials. J. Mech. Behav. Biomed. Mater. 79, 104–113 (2018).
    OpenUrl
Acknowledgments: Funding: This work was not externally funded. Author contributions: S.J. and A.A.Z. conceived the study. K.N., A.A.Z., and S.J. developed the analytical model. S.J. built the computational models and performed the simulations. S.J. and K.N. performed the material characterization study. S.J. and T.v.M. manufactured the specimens and the test setup, performed the experiments, and conducted the subsequent data analysis. S.J. and A.A.Z. wrote the manuscript. All authors edited the manuscript and approved its final version. 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.
View this article with LENS

Article Tools

My saved folders

Stay Connected to Science Advances

Related Content

Similar Articles in:

Citing Articles in:

Navigate This Article