### Experimental validation

The analytical model is validated using the results from our compressive tests on the bone-like cellular structure, which are presented in Fig. 3^{16}. The goal is to verify that buckling is the primary collapse mechanism of the structure as opposed to plastic deformation in the material. The following parameters are used for validating the analytical model, which are identical to the tested cellular structure (see Fig. 3): (W = H = 20,{text{mm}}); (alpha = 63.43^circ , beta = 25^circ), (gamma = 41.14^circ); (l_{ut} = 10,{text{mm}}), (l_{lt} = 5,{text{mm}}); (t_{w} = 2,{text{mm}}); (t = 30,{text{mm}}); (E_{S} = 320,{text{MPa}}) (Markforged Tough Nylon, 3D printing material). A force of 1,000 N is chosen arbitrarily within the linear elastic region of the force–displacement curve (Fig. 3), and divided over the number of half unit cells (six); the force per half unit cell is thereby (166.67,{text{N}}). A vertical displacement of (v_{Gf} = 0.84,{text{mm}}) at node F (left( {v_{Gf} } right)) was predicted from the analytical model by solving Eq. (2). The vertical displacement at node F (v_{Gf}) is used to calculate the strain and Young’s modulus from Eq. (3). The critical buckling load was then calculated using Eq. (9) for each column and presented in Table 3.

It can be inferred from Table 2 that the predicted Young’s modulus of the half unit cell is in agreement with the experiment within 7%. The experimental critical buckling load is also within the range of the predicted results. Although column AB provided the closest buckling load (Table 3) with the experiment, it is difficult to predict that this column will first buckle, as the condition of instability depends on imperfections that are inherent in the 3D printing process, as well as the load transfer through the structure. Regardless, the analytically predicted Young’s modulus and critical buckling load are conservative for design purposes. Importantly, the predicted minimum critical buckling stress (in column BC) is well within the linear elastic region of the material curve of the Markforged Tough Nylon 3D printing material^{16}, which verifies that buckling is the primary collapse mechanism of the structure as opposed to plastic collapse. The effective lengths and (k) factors are calculated from the Euler buckling load and presented in Table 3, which are slightly larger than 1 as expected, noting that (k = 1) represents a sway rigid-to-rigid column with one point of inflexion or a braced pinned-to-pinned column. The critical buckling loads listed in Table 3 signify that buckling can occur in either column AB, BC, CD or DE. Column EF is ignored because the critical buckling stress exceeds the linear elastic limit of the material (see^{16}). Hence, for the optimisation studies conducted hereafter, the critical buckling load is calculated for columns AB to EF for each unit cell configuration and the minimum buckling load is used to compute the strain energy density.

### Optimal design parameters

A computer algorithm was developed to execute many parallel simulations to obtain the design parameters ((alpha ,beta ,gamma ,l_{ut} ,l_{lt})) that produce the biomimetic unit cell with the optimal energy absorption. To this effect, the sub-cell angles and tie lengths were modified as follows: (alpha ,beta ,gamma{:} , left[ {5,175^circ } right]) in increments of (1^circ); and (l_{ut} ,l_{lt}{:} left[ {1,19,{text{mm}}} right]) in increments of 1 mm. This equates to 1.7 billion cases, which were run on a supercomputer using 170 threads (total runtime of 2.3 h). The results are reported in Table 4 and analysed herein. Case 1 ((alpha = beta = gamma = 90^circ , l_{ut} = l_{lt} = 1,{text{mm}})) results in vertical columns and produces the stiffest structure with the lowest energy absorption. This is expected because the columns are more significantly rigid under axial deformation as opposed to flexure. Case 2 ((alpha = 79^circ , beta = 30^circ , gamma = 34^circ , l_{ut} = l_{lt} = 9,{text{mm}})) results in a structure with a significantly low stiffness and a high buckling load compared to Case 1. This effectively increases the energy absorption relative to Case 1 according to Eq. (10). The design parameters are analysed individually in the following sub-section to gain further insight into this behaviour. Case 3 ((alpha = 122^circ , beta = 20^circ ,gamma = 29^circ , l_{ut} = 14,{text{mm}}, l_{lt} = 18,{text{mm}})) produces the most optimal structure with the lowest stiffness and highest energy absorption. Although the Euler buckling load is nearly identical to Case 1, the stiffness is significantly lower. According to Eq. (9), the buckling load depends on the ratio of the length of the tie restraint (L_{t}) to the length of the restrained column. The mechanical properties of the re-entrant auxetic structure, which is extensively investigated in the literature for protective applications, are also presented in Table 4 for comparison. It is evident that the bioinspired structure shows a significantly higher energy absorption capacity with a slight increase in mass (around 4%) under quasi-static loading. This promising result can motivate further numerical investigations to obtain the optimal bioinspired structure under extreme loadings in future work. The design parameters are also analysed individually in the following sub-section to gain further insight into this behaviour.

### Influence of individual design parameters

In this set of studies, the influence of each individual design parameter of the biomimetic unit cell (see Fig. 2) on its energy absorption, stiffness and Euler buckling stress is analysed. To this effect, the parameters of the baseline unit cell are modified individually from the reference baseline configuration ((alpha = 63.43^circ , beta = 25^circ), (gamma = 41.14^circ), (l_{ut} = 10,{text{mm}}) and (l_{lt} = 5,{text{mm}})). It is evident in Fig. 4 that the sub-cell angle (beta) has the most significant influence on the stiffness of the unit cell, followed by the length of the lower tie (l_{lt}). It is expected that (beta = 90^circ) results in the stiffest unit cell because the lower sub-cell columns become parallel such that they transmit the load via axial deformation predominantly, as opposed to flexure.

It can be observed in Fig. 5 that both sub-cell angles (beta) and (gamma), and the tie lengths (l_{ut}) and (l_{lt}) produce the unit cell with the highest Euler buckling stress. According to Eq. (9), the Euler buckling load is governed by the ratio of the length of the tie restraint ((L_{t})) to the length of the restrained column ((L)). The ratio (L_{t} {/}L) and the critical column are listed in Table 5 for each individual design parameter that produced the maximum buckling stress. It can be deduced that the plateau region e.g. for (45^circ < alpha < 103^circ) gives the same (L_{t} {/}L) ratio and thereby produces an identical Euler buckling stress.

It is evident in Fig. 6 that the sub-cell angle (gamma =)(31^circ) and upper cell wall angle of (l_{ut} = 7.8,{text{mm}}) result in unit cells with the highest strain energy density. According to Eq. (10), the optimal energy absorption can be obtained by minimising the Young’s modulus E whilst maximising the critical buckling stress (sigma_{CR}). According to the five-parameter optimisation in Table 4, modifying each design parameter individually will significantly underestimate the unit cell with the maximum strain energy density ((U_{max} = 18.8 ,{text{mJ/mm}}^{3})). It can be deduced from Figs. 4 and 5 that minimising the sub-cell angle (gamma) produces a relatively constant low stiffness (sim 13,{text{MPa}}) and a maximum buckling stress of (sim 19,{text{MPa}}), respectively. The sub-cell angle (beta) and lower tie length (l_{lt}) can overcome this restriction to reduce the stiffness further and thereby obtain the unit cell with the optimal energy absorption (see Fig. 4).

### Influence of combined design parameters

In this set of studies, the combined effects of two design parameters of the biomimetic unit cell (see Fig. 2) on its energy absorption, stiffness and Euler buckling stress are analysed. The parameters of the baseline unit cell are modified from the reference configuration (alpha = 63.43^circ ,{ }beta = 25^circ), (gamma = 41.14^circ); (l_{ut} = 10,{text{mm}}), (l_{lt} = 5,{text{mm}}). The corresponding unit cells with the maximum values are also illustrated near each surface plot. It is evident from Fig. 7a,c that the sub-cell angle ({upbeta }) (see Fig. 2) is a key design parameter that dictates stiffness, whereby the lower columns of the unit cell are relatively vertical. This trend is consistent with that observed in Fig. 4, which showed that modifying (beta) alone produces the stiffest unit cell by far. The unit cell with the highest stiffness ((E_{max} = 50.3,{text{MPa}})) is obtained when (beta = 99^circ) and (gamma = 69^circ) (Fig. 7c). This is an expected result, given that the columns of both sub-cells are relatively vertical, which promotes axial deformation over flexure. A relatively stiff unit cell ((E = 44.7,{text{MPa}})) can also be obtained by other combinations of sub-cell angles e.g. (alpha = 76^circ) and (beta = 88^circ) (Fig. 7a), which further reinforces that the stiffness is controlled by (beta) as a governing parameter. In contrast, Fig. 7b,d show that the unit cells with the highest stiffness have more oblique members such that flexural deformation becomes more prominent, which consequently reduces the stiffness of the structure.

The highest critical buckling stress (({upsigma }_{{{text{CR}}}} = 19.7,{text{MPa}})) is obtained from ({upalpha } = 80^circ) and ({upbeta } = 33^circ) (Fig. 8a). This is followed by ({upsigma }_{{{text{CR}}}} = 19.3,{text{MPa}}), which results from the tie lengths ({text{l}}_{{{text{ut}}}} = {text{l}}_{{{text{lt}}}} = 9,{text{mm}}) (Fig. 8d). From the discussion preceding Fig. 5 and Table 5, and referring to Eq. (9), the critical buckling load depends on the ratio of the restraining tie length to the length of the restrained column ({text{L}}_{{text{t}}} {text{/L}}). It is evident that the unit cell in Fig. 8a has the shortest tie (attached to longer columns CD and DE), which effectively yields a small ({text{L}}_{{text{t}}} {text{/L}}) ratio and thereby attracts a higher buckling stress. Longer plateau regions can be observed in the plots of Fig. 8b,c because the unit cells that produce the highest buckling stress have a more uniform structure, whereby the relative lengths between the struts and ties are similar.

The unit cell with the highest strain energy density (({text{U}}_{{max}} = 14.4,{text{MPa}})) is obtained when ({upalpha } = 88^circ) and ({upgamma } = 32^circ) (Fig. 9b). Recalling from Eq. (10) and the discussion preceding Fig. 6 that the peak strain energy density is obtained when the critical buckling stress is maximised and the stiffness is minimised, it is evident that the unit cell illustrated Fig. 9b is designed to balance axial and flexural deformation, whereby the oblique orientations of its columns reduce the stiffness. The unit cell also has relatively similar column and tie lengths, which thereby produces a relatively high critical buckling load ((sigma_{CR} = 18.9,{text{MPa}})). Note that the critical buckling load does not vary within a large interval compared to the stiffness, given that the (L_{t} {/}L) is constrained by the size of the unit cell. Importantly, other design parameters can be modified to produce unit cells with high energy absorption e.g. ({text{U}} = 14.2,{text{MPa}}) is obtained from (beta = gamma = 153^circ) (Fig. 9c). This reinforces that (gamma) is a governing parameter that has a significant influence on the strain energy density of the bone-like unit cell. The unit cells in Fig. 9a,d that produce the peak strain energy density have relatively short tie lengths compared to the struts, which has the effect of increasing the critical buckling stress. The oblique struts of these unit cells produce more prominent flexural deformations, which consequently reduces the Young’s modulus and thereby increases the strain energy density according to Eq. (10).