## 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).

## 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, *C*_{10H}) and visco-hyperelastic (long-term Neo-Hookean constant,

) beams. A single-term Prony series represented the viscous behavior of the visco-hyperelastic beam (dimensionless coefficient of the Prony series, *g*_{1}; relaxation time, τ_{1}). We define the nondimensional time, *t**, as *t** = *t*/τ_{1} and the nondimensional strain rate, ε′, as

. There are three material parameters (i.e.,

${C}_{10V}^{\infty},{C}_{10H}$, and *g*_{1}) 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, *E _{V}*, should exceed the elastic modulus of the hyperelastic beam,

*E*(i.e.,

_{H}*E*>

_{V}*E*). 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

_{H}(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

${C}_{10V}^{\infty}/{C}_{10H}$ and *g*_{1} 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

or *g*_{1} 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

or *g*_{1} 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

and *g*_{1} 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

and *g*_{1} or when the sum of both is large enough (i.e.,

), 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 *g*_{1} and

(i.e.,

${g}_{1}+{C}_{10V}^{\infty}/{C}_{10H}\gg 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, *F _{el}*, is proportional to the lateral deformation of the beam it supports (

*F*=

_{el}*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

(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.

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,

$\stackrel{\u0307}{\mathrm{\epsilon}}$, 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

${C}_{10V}^{\infty}/{C}_{10H}=0.2,{g}_{1}=0.2,{\mathrm{\tau}}_{1}=1001/\mathrm{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.

### 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).

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 (

${C}_{10V}^{\infty}/{C}_{10H}=0.7,{g}_{1}=0.7,{\mathrm{\tau}}_{1}=0.011/\mathrm{s}$;

${C}_{10V}^{\infty}/{C}_{10H}=0.9,{g}_{1}=0.6,{\mathrm{\tau}}_{1}=11/\mathrm{s}$; and

${C}_{10V}^{\infty}/{C}_{10H}=0.5,{g}_{1}=0.8,{\mathrm{\tau}}_{1}=11/\mathrm{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*).

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., 10^{3} 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.

### 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.

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.

## 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

$\stackrel{\u0307}{\mathrm{\epsilon}}={10}^{-4}$and

$\stackrel{\u0307}{\mathrm{\epsilon}}={10}^{3}$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 as

$${C}_{10}^{R}(t)={C}_{10}^{0}(1-{\displaystyle {\sum}_{i=1}^{n}}{g}_{i}(1-{e}^{\raisebox{1ex}{$-t$}\!\left/ \!\raisebox{-1ex}{${\mathrm{\tau}}_{i}$}\right.}\left)\right)$$(1)in which the instantaneous modulus is determined from

$${C}_{10}^{\infty}={C}_{10}^{0}(1-{\displaystyle {\sum}_{i=1}^{n}}{g}_{i})$$(2)where *n*, *t*, τ* _{i}*,

*g*,

_{i}, and

${C}_{10}^{\infty}$ 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, *W _{V}*, and hyperelastic,

*W*, beams are, respectively, given as (

_{H}*39*)

(3)where *I*_{1} and *I*_{3} are the first and third invariants of the left Cauchy-Green deformation tensor, **B**, respectively. The terms *p _{V}* and

*p*represent the arbitrary hydrostatic pressures that are defined such that the imposed incompressibility conditions are satisfied. The Cauchy stress,

_{H}**σ**, could be derived from the strain energy function as

(4)

Given that

$\frac{\mathrm{\partial}{I}_{1}}{\mathrm{\partial}\mathbf{B}}=\mathbf{I}$and

$\frac{\mathrm{\partial}{I}_{3}}{\mathrm{\partial}\mathbf{B}}={I}_{3}{\mathbf{B}}^{-1}$, the Cauchy stresses of the visco-hyperelastic and hyperelastic beams can be calculated as

$$\begin{array}{c}{\mathbf{\sigma}}_{V}=\frac{2}{\sqrt{{I}_{3}}}\left\{{C}_{10V}^{0}\right(1-{\displaystyle \sum _{i=1}^{n}}{g}_{i}(1-{e}^{\raisebox{1ex}{$-t$}\!\left/ \!\raisebox{-1ex}{${\mathrm{\tau}}_{i}$}\right.}))\mathbf{B}+{p}_{V}{I}_{3}\mathbf{I}\}\\ {\mathbf{\sigma}}_{H}=\frac{2}{\sqrt{{I}_{3}}}\{{C}_{10H}\mathbf{B}+{p}_{H}{I}_{3}\mathbf{I}\}\end{array}$$(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.,

*I*

_{3}= λ

_{1}λ

_{2}λ

_{3}= 1) and considering uniaxial deformation in the direction 1, the left Cauchy-Green deformation tensor is obtained as

(6)

The incompressibility condition and uniaxial loading conditions (i.e., σ_{22} = σ_{33} = 0) could then be used to determine the arbitrary pressures *p _{V}* and

*p*. The stress in the direction 1 can be determined by combining Eqs. 5 and 6 as

_{H}(7)

Assuming a single-term Prony series, the instantaneous elastic modulus in the direction 1 is given by

$$\begin{array}{c}{E}_{V,11}=\frac{\mathrm{\partial}{\mathrm{\sigma}}_{V,11}}{\mathrm{\partial}\mathrm{\epsilon}}=\frac{\mathrm{\partial}{\mathrm{\sigma}}_{V,11}}{\mathrm{\partial}\mathrm{\lambda}}=2{C}_{10V}^{0}(1-{g}_{1}(1-{e}^{\raisebox{1ex}{$-t$}\!\left/ \!\raisebox{-1ex}{${\mathrm{\tau}}_{1}$}\right.}\left)\right)(2\mathrm{\lambda}+\frac{1}{{\mathrm{\lambda}}^{2}})\\ {E}_{H,11}=\frac{\mathrm{\partial}{\mathrm{\sigma}}_{H,11}}{\mathrm{\partial}\mathrm{\epsilon}}=\frac{\mathrm{\partial}{\mathrm{\sigma}}_{H,11}}{\mathrm{\partial}\mathrm{\lambda}}=2{C}_{10H}(2\mathrm{\lambda}+\frac{1}{{\mathrm{\lambda}}^{2}})\end{array}$$(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 as

$\begin{array}{c}{E}_{V}=6{C}_{10V}^{0}(1-{g}_{1}(1-{e}^{\raisebox{1ex}{$-t$}\!\left/ \!\raisebox{-1ex}{${\mathrm{\tau}}_{1}$}\right.}\left)\right)\\ {E}_{H}=6{C}_{10H}\end{array}$(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, *P _{cr}*, is higher than the Euler’s buckling load for an equivalent unsupported beam,

*P*, and can be obtained from the solution to the following implicit equation (

_{e}*34*)

(10)where

$x=\frac{{P}_{\mathit{cr}}}{{P}_{e}}$,

$y=\sqrt{\frac{k{l}^{4}}{\mathit{EI}}}$, *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

(or

$x=\frac{2y}{{\mathrm{\pi}}^{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}

(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 have

(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} to

(13)

Because the beams are geometrically identical and given that

$y=\sqrt{\frac{k{l}^{4}}{\mathit{EI}}}$, Eq. 13 can be rewritten as

$$\sqrt{\frac{{E}_{V}{l}^{4}}{{E}_{H}I}}>\sqrt{\frac{{E}_{H}{l}^{4}}{{E}_{V}I}}$$(14)or simply

$${E}_{V}>{E}_{H}$$(15)

Combining Eqs. 9 and 15 gives the condition for buckling to the right as

$$\frac{{C}_{10V}^{0}}{{C}_{10H}}>\frac{1}{1-{g}_{1}(1-{e}^{\raisebox{1ex}{$-t$}\!\left/ \!\raisebox{-1ex}{${\mathrm{\tau}}_{1}$}\right.})}$$(16)

Using Eq. 2, Eq. 16 could be rewritten as

$$\frac{{C}_{10V}^{\infty}}{{C}_{10H}}>\frac{1-{g}_{1}}{1-{g}_{1}(1-{e}^{\raisebox{1ex}{$-t$}\!\left/ \!\raisebox{-1ex}{${\mathrm{\tau}}_{1}$}\right.})}$$(17)

For extremely high strain rates (i.e., *t*/τ_{1} → 0), the condition for buckling to the right can be simplified as

(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., *E _{V}* >

*E*).

_{H}In our computational analyses, we varied the ratio of the Neo-Hookean coefficient of the visco-hyperelastic layer (

${C}_{10V}^{\infty}$) to the hyperelastic layer (*C*_{10H}) 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 (*g*_{1}) 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, *C*_{10H}, 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,

, and the Prony series parameters (τ* _{i}* and

*g*) were obtained on the basis of the relaxation tests performed for the visco-hyperelastic material using a previously described methodology (see the Supplementary Materials) (

_{i}*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.

**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.