Large recoverable elastic energy in chiral metamaterials via twist buckling

Finite element analysis

All FEAs were performed using ANSYS, considering large nonlinear geometrical deformations. Typically, one end is fixed while a displacement boundary condition is applied to the other end, allowing us to obtain the von Mises stress, strain energy and reaction force.

Bending buckling mode

We consider the compressive buckling of a doubly clamped rod. Its bending deflection is assumed to be w(x) = a(1 − cos bx), b = 2nπ/L, where L is the rod length, n denotes the order of bending buckling and x is the position along the rod (Supplementary Note 1). The maximal deflection is w_max = 2a. The rod in the chiral metacell bends in the n = 1/2 mode. A slender primitive rod in non-chiral lattices typically follows the first-order buckling mode n = 1 under suitable constraints, but without stiff constraints, some global buckling mode may lead to n = 1/2 for the oblique rods (Extended Data Fig. 9b).

We note that the assumption for deflection, w(x) = a(1 − cos bx), is valid for moderate deformation ε < 0.2. When a rod is deeply buckled, it exhibits more complex deformation patterns. In FEA, these transformations in the bending deflection are adaptively accounted for through geometrical nonlinearity. This transformation is not considered in the analytical theory, which explains the increasing deviation from FEA results in Fig. 2 for ε > 0.2.

Compression buckling of non-chiral rods

Here we briefly describe the results of our analytical model for the compression buckling of a rod. The complete modelling process and mathematical details are provided in Supplementary Note 1. When an axial compression displacement (Δ) is applied to a doubly clamped rod of radius r, length L₀ and Young’s modulus E_s, the axial reaction force is F_1rod. At first, the stress σ_cpr is induced by pure compression. Buckling occurs at a critical force

where n denotes the order of buckling mode. The critical axial stress for buckling is

When F_1rod < F_c, the rod undergoes uniform strain ε = Δ/L₀ and stress σ_cpr = E_sε, storing energy

where V_s = πr²L₀ denotes rod volume. In this case, von Mises stress is σ_v = σ_cpr. Moreover, the critical energy for buckling is

This equation shows that the critical energy is proportional to n⁴. If the buckling order changes from n = 1 to 1/2, the critical energy greatly reduces to 1/16.

When buckling occurs, the stress induced by pure compression (σ_cpr) remains nearly constant. In the general case with n = 1, bending induces maximal tensile or compressive stresses on the surface near the clamped ends

where 2a denotes the bending deflection. Consequently, maximal von Mises stress σ_v = σ_bend + σ_cpr is found there. In this case, the strain energy stored in a compressive buckling rod is

This equation indicates that σ_cpr gets an eight times higher increment ratio of energy than σ_bend when increasing stress. Thus, for a specified moderate stress, such as σ_v = 0.1E_s, a thicker rod obtains higher U_1rod/U_s until buckling disappears at this stress level. This equation also indicates that combining bending and compression increases the increment ratio of U_1rod/U_s as strain ε = Δ/h₀ increases.

We note that when analysing the angled rod and chiral rod, the vertical compression displacement is defined as s. In the above equations, Δ denotes the axial compression displacement. Δ = s for a vertical rod with angle θ = 90°.

Parameter generalization

Based on the theory above, we adopt \(F_\rms=E_\rmsI/L_^2\) to normalize the compressive force. Moreover, specifying ξ = σ_cpr/σ_bend, the energy inside a buckled rod can be simplified as

Here, V_s = πr²L₀ is the rod volume and λ = σ_v/E_s is the normalized strength. Then, the enthalpy of a single rod is

The factor f can be approximated as f ≈ 2.5/1,000 for ξ → 0 and λ < 0.2. For convenience, it is reasonable to use ϕ_s = E_s/1,000 and U_s = V_sϕ_s to normalize the enthalpy and total energy, respectively.

Analytical model of chiral buckling

The process of deriving the analytical model for chiral buckling is extensive and involves rigorous mathematical calculations. In Supplementary Information section ‘Chiral twist buckling theory’, we systematically elaborate on the geometry, deformation, force, energy, and stress induced by in-plane bending, out-of-plane bending, in-rod twisting, helix and compression modes (Supplementary Notes 2 and 3). Definitions of all variables are also listed there. Here we summery the analytical theory (also in Supplementary Note 3.8), with key parameters labelled in Extended Data Fig. 1.

When a vertical compression displacement s is applied to half a chiral metacell, a variable rotation angle θ forms between the two tori. The equivalent strain is ε = s/h₀. The reaction force and elastic energy contributed by a single rod are F_1rod, U_1rod, respectively. Point B is fixed on torus O₁. Point A₀ is the original point on torus O₂, which moves to point A under the compression displacement s. Length is given as

The in-plane bending deflection is denoted as 2a_in and the out-of-plane bending deflection is denoted as 2a_out. Specifying β = ∠ABA₀, 2a_in = L_s sin β, L_x = L_s cos β. The term 2a_out is defined by length |CD| shown in Extended Data Fig. 1b.

The moments induced by in-plane bending and in-rod twisting (Supplementary Note 3) are

Here b = π/L_x, I = πr⁴/4 is the second moment of the rod, G is the shear modulus, θ_r = θ-γ. The bulge-out angle γ arises from the helical deformation of rod, as shown in Extended Data Fig. 1c.

The components of the two moments in three-dimensional spaces are

where n_in is the unit normal vector of the in-plane bending and n_twist is the unit normal vector in the twisting deformation. The concentrated force applied on a single rod denotes (F_x, F_y, F_z). We use the following equations to analyse chiral buckling deformation and calculate energy.

1. 1.

   Force equilibrium equation (Supplementary Fig. 7)

   $$\left\{\beginarrayc2M_\rmin,y-F_x(h_-s)-F_z(R_1-R_2\cos \alpha )=0,\\ 2M_\rmin,\rmz+F_xR_2\sin \alpha +F_y(R_1-R_2\cos \alpha )=0,\\ M_\rmin,\rmz+T_\rmtwist,z+F_xR_2\sin \alpha -F_yR_2\cos \alpha =0.\endarray\right.$$

   Here α = θ + α₀. Based on this equation, the compressive load applied on a chiral rod is F_1rod = −F_z.

2. 2.

   Compatibility equation of deformation

   $$L_\rmin=L_-\eta \varDelta _\rmout-\varDelta _\rmcpr,$$

   where \(\Delta _\rmo\rmu\rmt=\rm\pi ^4a_\rmo\rmu\rmt^2/4L_\rms\) denotes the shortening made by out-of-plane deformation; ∆_cpr is the pure compression displacement along the rod, and \(\Delta _\rmc\rmp\rmr < \rm\pi r^2/L_\); \(L_\rmi\rmn=\fracE_\rmk(\rm\pi ,-a_\rmi\rmn^2b^2)b\) denotes the length of the rod under in-plane deformation, and E_k(·) is the elliptical equation; η denotes a correction factor explained in Supplementary Note 3.4.

3. 3.

   Energy equation

   $$U_1\rmr\rmo\rmd=\fracE_\rmsI\rm\pi ^4a_\rmi\rmn^24L_x^3+\fracE_\rmsI\rm\pi ^4a_\rmo\rmu\rmt^2L_^3+\fracGI_\rmp\theta ^22L_+\fracK_\rmc\,\Delta _\rmc\rmp\rmr^22.$$

   Here K_c = E_sπr²/L₀ denotes the longitudinal stiffness of the rod. Based on the energy method, the compressive load is

   $$F_1\rmrod=\frac\partial U_1\rmrod\partial s.$$

4. 4.

   Equation of stresses

   The maximal stresses induced by in-plane bending, in-rod compression and twisting are

   $$\beginarrayccc\sigma _\rmi\rmn=\fracE_\rmsr\rm\pi ^2a_\rmi\rmnL_\rms^2, & \sigma _\rmc\rmp\rmr=\fracF_1\rmr\rmo\rmd\rm\pi r^2, & \tau =\fracG\theta _\rmrrL_.\endarray\,$$

   Here σ_in and σ_cpr are in the same direction, σ_norm = σ_in + σ_cpr. At the same point, the direction of twisting shear stress τ is perpendicular to σ_norm. The stress in the radial direction of the rod is negligible. According to the theory of equivalent von Mises stress,

   $$\sigma _\rmv=\sqrt\sigma _\rmnorm^2+3\tau ^2.$$

We adopt this equivalent stress to evaluate the strength of the rods in the chiral structure.

Stress evaluation

The shear modulus is G = E_s/2(1 + υ), where Poisson’s ratio is υ = 0.3. In the chiral metacell, both the in-plane and out-of-plane bending deflection follow the function 1 − cos(πx/L), with their deflection amplitudes denoted as 2a_in and 2a_out, respectively. In-rod twisting generates shear stress given by τ = E_sπr/2.6L₀ = σ_bend(L₀ − s)²/2.6L₀πa. A proper chiral structure sets L₀ ≈ 4R; under large deformation, s ≈ R and the in-plane bending deflection 2a ≈ 2R. Thus, we find that τ ≈ 0.3σ_bend. Under the same compressive displacement s, σ_norm in the chiral rod and the lattice rod (Fig. 1a) are approximately equal.

Performance evaluation

We use the following quantities to evaluate the performance.

Here, force F, energy U, volume V, mass m, the projected area in load-bearing direction A and the vertical compression displacement s, with the subscript ‘cell’ are measured on the unit cell. The equivalent elastic modulus E denotes the slope of the ε–σ_eq curve at a small strain ε < 0.02. The maximal stress on the buckling plateau, σ_bk = max(σ_eq), symbolizes the load-bearing strength of the entire metamaterial. These variables are normalized for evaluation as E/E_s, σ_bk/E_s, ϕ/ϕ_s and ρ/ρ_s, where ρ_s is the mass density of the rod material, ϕ_s = E_s/1,000 is defined in the section ‘Parameter generalization’.

For a chiral metamaterial, there are N rods around the circle with radius R. Here we adopt the integral N = Round(πR/2r) with enough space left to avoid contact between the deformed arms. Thus, V_cell = 4h₀(R + r)², A_cell = 4(R + r)², F_cell = NF_1rod, U_cell = NU_1rod, m_cell = Nm_rod, where m_rod is the mass of a rod. Here we adopt the maximal cubic volume and square area, instead of the smaller cylinder and circular sizes, to evaluate the performance of chiral metamaterials.

For octahedron lattices, we calculate only 1/8 metacell corresponding to one rod. The parameters are defined as follows: F_cell = F_1rod, U_cell = U_1rod, m_cell = m_rod, \(V_\rmcell=(L_0\rmr^3\cos ^2\theta \sin \theta )/2\), \(A_\rmcell=(L_0\rmr^2\cos ^2\theta )/2\) where θ denotes the oblique angle between the rod and the horizontal. L_0r denotes the distance between the neighbour connection nodes.

For prism lattices, we calculate only 1/4 metacell corresponding to one rod. F_cell = F_1rod, U_cell = U_1rod, m_cell = m_rod, \(V_\rmcell=l_\rmpL_0\rmr^2\cos ^2\theta \sin \theta \), \(A_\rmcell=l_\rmpL_0\rmr\cos \theta \), where l_p denotes the distance between neighbouring rhombuses, that is, the lattice constant along the x-axis in Fig. 4. Here we adopt l_p = 4r to form a dense lattice.

For octahedron and prism lattices, the useful length of a rod is L₀. In the above equations, L_0r denotes the distance between two connection nodes in lattices. If L_0r = L₀, we will obtain ρ/ρ_s > 1, when the oblique angle θ > 70°, which is impractical. In reality, the connection nodes also occupy space, as shown by the samples in Fig. 4. Therefore, the practical lattices have L_0r > L₀. For our calculations, we adopt a reasonable value L_0r = 1.1L₀ for the rod with r = 1.5 mm.

Samples

The rubber chiral metacell in Fig. 4a features N = 8 rods with r = 1.5 mm, whereas the titanium version has N = 20 and r = 0.6 mm. All samples have h₀ = 30 mm or L₀ = 30 mm. Both the rubber and titanium chiral samples in Fig. 4a,b have R = 7.5 mm and α₀ = 5°. The oblique angle in prism lattices is θ = 40°. The distance between the neighbour parallel rods is 3 mm for the TC4 lattice and 10 mm for the rubber lattice. The rods inside octahedron and tetradecahedron lattices have equal length. Based on Fig. 3, the in-rod maximal Mises stresses in chiral, prism and octahedron lattices are approximate when the same global strain ε is specified. Although tetradecahedron and tensegrity lattices can theoretically endure higher global deformation, their neighbouring rods make contact when ε > 0.4. Our performance evaluation takes these differences into account. Other parameters and results are listed in Extended Data Table 1.

In Fig. 4c, chiral-20° has R = 5.5, 2r = 1.8, α₀ = 20°, h₀ = 20 mm. Chiral-50° has R = 6, 2r = 1.7, α₀ = 50°, h₀ = 20 mm. E_s ≈ 5.5 MPa for rubber samples in Fig. 4e, whereas other rubber samples have E_s ≈ 15 MPa.

Experiments

Cyclic compression experiments were performed. The rods in newly fabricated metal samples are initially straight. However, after the first compressive cycle, the all-metal samples exhibit some residual plastic deformation. In this case, α₀ increases, and the rod becomes slightly bent. Moreover, the first cycle strengthens the material, improving its overall strength. Subsequent compression cycles are fully repeatable, that is, no further plastic deformation occurs. The enthalpy is calculated as \(\phi =\int \sigma _\rmeq\rmd\varepsilon \).

For the prism lattice, the one-layer sample shown in Extended Data Fig. 9a consists of two half metacells. We tested the deformation of two-layer and four-layer samples to confirm consistency. For multi-layer prism lattices, out-of-plane bulge-out deformation (corresponding to 1/2-order buckling) occurs when there are no lateral constraints (see Extended Data Fig. 9b and Supplementary Video 5). This deformation reduces the buckling strength σ_bk and enthalpy ϕ by a factor of 10 compared with the desired in-plane buckling mode shown in Fig. 1a (Extended Data Fig. 9h). Therefore, as shown in Extended Data Figs. 9c and 10e, we placed the prism lattices inside a box to constrain the lateral bulge-out deformation.