Energy Absorption Performance of Bio-inspired Honeycombs: Numerical and Theoretical Analysis

Energy absorption performance has been a long-pursued research topic in designing desired materials and structures subject to external dynamic loading. Inspired by natural bio-structures, herein, we develop both numerical and theoretical models to analyze the energy absorption behaviors of Weaire, Floret, and Kagome-shaped thin-walled structures. We demonstrate that these bio-inspired structures possess superior energy absorption capabilities compared to the traditional thin-walled structures, with the specific energy absorption about 44% higher than the traditional honeycomb. The developed mechanical model captures the fundamental characteristics of the bio-inspired honeycomb, and the mean crushing force in all three structures is accurately predicted. Results indicate that although the basic energy absorption and deformation mode remain the same, varied geometry design and the corresponding material distribution can further boost the energy absorption of the structure, providing a much broader design space for the next-generation impact energy absorption structures and systems.


Introduction
Originated from biological evolution, natural materials offer a source of inspiration for incredibly efficient designs of structural applications [1]. Accordingly, advanced materials with high specific strength, stiffness, and energy absorption have been designed for aerospace and automotive engineering [2,3].
Plants or animals can barely survive if the desired biological structures with optimal geometric parameters are not achieved [4]. If a certain bio-structure has a function of impact resistance for protection, we can intuitively expect a good energy absorption capability of the structure toward the common scenarios. It is not uncommon, for example, for cellular structures in plants and animals that are known to support great weights, like bone [5], to contain a strong outer wall and a low density or entirely vacant central chamber such as grass [6]. Although these evolutions are not the same, they perform similar purposes, and as such much can be gleaned from their relationships.
Consequently, quite a few pioneering studies have been conducted. For example, nature has finetuned a variety of biological composites with remarkable mechanical properties at different scales and inspired researchers to design bio-inspired composites. Nacre-inspired composites can exhibit superior impact resistance [7] as well as ballistic performance [8]. Various helical structures were proposed to mimic a helicoidal arrangement of the mineralized area on the endocuticle of the dactyl club [9,10], while similar helix shapes can also be found in DNA-inspired structures for impact mitigation [11]. Meanwhile, primary efforts have been placed on developing novel structures with distinct orders at different scales [12,13] to mimic the hierarchical organization of biological materials (e.g., coconut shell [14], collagen [15], and conch shell [16]). Notably, each hierarchical order may work synergistically, and thus, result in an exceptional impact resistance for the overall structures. Furthermore, based on the traditional honeycombs, numerous bio-honeycombs have been designed for lightweight energy absorbers, such as pomelo peel-inspired honeycomb [17], wavy honeycombs inspired by woodpeckers' beaks [18], and spiderweb honeycombs [19].
Naturally, if a primary relationship and a combination of similarities between bio-inspired honeycomb structures are fully explored, we may obtain design inspiration to cope with the energy absorption problem using pre-existing models and assist further optimization of energy-absorbing structures. In this study, we investigate the out-of-plane energy absorption capabilities of three types of honeycombs that are inspired by biological constraints. In Sect. 2, the factors for structural design and the finite element (FE) models are introduced. Section 3 presents the typical results from different structures and the crashworthiness indicators for energy absorption capacity. In Sect. 4, a theoretical model is developed to analytically describe the energy absorption performances of the proposed structures. In addition, a parametric study is conducted to characterize the effect of geometric parameters on the dynamic behaviors of the structures and uncover the fundamental mechanisms on the bio-inspired structural design.

Structural Design Concept
This study reviews two biological parameters of structural design and their relation to energy absorption capabilities. The first is how biological cells naturally grow to best fit into their environment via space-saving geometries, while the second is how geometries of known bio-inspired honeycombs relate to each other. This involves regarding one structure as the graph of another and exploring how shapes evolve over continued iterations.
One shape highly optimized for space-filling geometries is the Kelvin cell ( Fig. 1a, b), which is meant to fit three-dimensional (3D) spaces well, while their cross section can be an inspiration for honeycomb tiling as well (Fig. 1c). Moreover, for the Weaire-Phelan structure (shortened to Weaire in the following sections), the cross section is comprised of three different geometries that depend on the angle (Fig. 1d), making it difficult to find a representative unit cell. Therefore, an approximate cross section is taken from one of its geometries (Fig. 1e), yielding the unit cell not directly built from a Weaire cell but rather an inspiration to help gauge the effects of this geometric interpretation.
The relationships between previously known geometries are further investigated by the graph strategy in Fig. 1f-k. While the process for determining these structures may be conducted over complex 3D geometries, for this study, this method is performed over 2D honeycomb cross sections extended into the third dimension. Surprisingly the most basic honeycombs are graph pairs to each other to a minimal extent, after which a more stable graph pairing arises. Kagome devolves into triangles, which in turn devolves into hexagonal and self-stable ones for all future graph expansions. Similarly, a square is its graph pair where the following graph is 45 degrees offset from its previous iteration. Furthermore, certain more complicated geometries are paired with each other in similar manners until hierarchy overtakes. When transformed into graphs, spiderweb [19] and pomelo-inspired honeycombs [17], for example, will yield similar surfaces if the rhombuses in the pomelo-inspired honeycomb are neglected. The center hexagon would be filled with triangles if the shapes were truly graph pairs in one of the geometries. Despite this, the fact that similar structures are created independently but share a great similarity in their graphs suggests a more primal relationship between these structures. For these graph-created structures, hierarchy plays a factor that can only be noted in the more complex geometry. As the spiderweb and pomelo devolve, the size of the cells decreases, and the hierarchal  [20] (a); model of a 3D Kelvin cell [21] (b); and Kelvin-inspired honeycomb geometry (c). Derivation of Weaire-inspired honeycomb from segments of Weaire Phelan structures [22] (d) and interpretations onto a two-dimensional (2D) plane that allows it to be connected with similar geometries in a Weaire-inspired unit cell (e). Relationships between graph pairs, geometries whose intersections of nodes match pre-existing geometries of known honeycombs: Triangular, hexagonal, and Kagome (f, g, and h) are all paired to triangular cells, while Pomelo and spiderweb are paired to spiderweb [17,19] (i and j). K shows the derivation of floret from a hierarchical honeycomb [23,24] detail decreases until the structure stabilizes into a triangular honeycomb. From a single hierarchy of hexagonal honeycomb with triangular walls [23,24], this method may be used to produce a floret honeycomb (Fig. 1k).

Finite Element Modeling
Although all unit cells were designed to be fitted in a repetitive array, each honeycomb was modeled as a singular unit cell with a characteristic wall length of λ = 10 mm. Accordingly, FE models of different unit cells were developed in ABAQUS/Explicit with a unit height of 50 mm. The base material was 6061-T6 aluminum with Young's modulus E = 68 GPa, yield stress σ y = 240 MPa, ultimate stress σ u = 290 MPa, Poisson's ratio υ = 0.35, and density ρ = 2700 kg/m 3 . Each cell was modeled by shell elements (S4R) with thickness t = 0.5 mm. Based on a mesh convergence study, a mesh size of 0.5 mm was selected to balance the computational cost and accuracy. A rigid plate with a uniform loading rate of 2 m/s was applied on the top for crushing the honeycombs, while the bottom edges of the honeycombs were fixed. The loading rate was chosen to aid in computational speed. Impact velocities from 0.5 to 2 m/s were measured and no notable change in force-displacement curve was observed. The measured internal energy also did not imply inertial effects became prevalent. This impact velocity is within a range of velocities used for quasi-static compression in other studies [17,23]. Surface-to-surface contact with a friction coefficient of 0.35 was used as it is proportional to existing friction conditions in literature [25], and a general contact constraint was used to restrict the model from clipping into itself. The material was modeled as linear hardening in the plastic stage.

Model Validation
To ensure accuracy of the FE models, samples were machined of each cell in question by wire EDM under quasi-static compression. By comparing the experimental compression, which was performed at 5 mm/min, the numerical analysis of 2 m/s displacement used in the model was confirmed to be representative of the quasi-static regime. The experiment also showed elastic, plateau, and densification stages, which are consistent with literature and the FE model, confirming that the numerical analysis was performed reasonably and is a viable tool to perform parametric analysis upon.

Representative Results
Representative force-displacement curves together with the deformation modes of the three structures are presented in Fig. 2, where three stages that are analogous to traditional honeycombs can be observed, including an elastic stage followed by an initial peak, a plateau region, and a densification stage. Meanwhile, a progressive buckling deformation mode accounts for the fluctuations in the plateau regime.
According to the force-displacement curves, a series of crashworthiness criteria is defined to characterize the energy absorption performance of the honeycombs, including specific energy absorption (SEA), mean crushing force (MCF), and crushing force efficiency (CFE). Herein, SEA is a critical criterion to evaluate the absorbed energy per mass during the crushing process, which is given by where m is the total mass of the honeycomb. EA is the total energy absorption until densification and is defined as Another key factor to regard for the validity of the structures, especially in automotive safety environments, is the crushing force efficiency, CFE, which is given by Here, PCF is the peak crushing force, which can be directly obtained from the curve, typically before the onset of the plateau region. MCF is defined according to Eq. (4) as As the energy absorption of the samples is taken from the beginning of the load to the point of densification, these bounds must be determined. The onset of the load is taken at the initialization of with R being the ratio of EA to the maximum force. For its definition, d occurs when the slope of the curve changes from positive to negative. This peak defines the rapid increase in force consistent with the densification region. We compare the crashworthiness criteria of the proposed honeycombs and the traditional honeycomb (i.e., hexagonal honeycomb) with consistent wall length and thickness in Table 1. Given the same geometric parameters, the proposed honeycombs have superior energy absorption performances over the hexagonal honeycomb.

Theoretical Analysis
A theoretical model is developed in this section to achieve a reasonable approximation of the crushing resistance of specified unit cells [27]. The absorbed energy during the crushing process consists of bending energy W B and membrane energy W m . For a full folding wavelength 2H, the bending energy can be determined by calculating the total wall length of each honeycomb Nλ. Thus, where λ is a predefined length, and M is the fully plastic bending moment being defined as where σ o is the flow stress where n is the strain hardening coefficient measured as 0.06 [28]. The numbers N of Kelvin, Weaire, and floret structures are 14, 50.8, and 15, respectively. The membrane energy can be calculated based on the constitutive elements of different honeycombs (Fig. 3). According to [29], the membrane energy of the corner element (type I) in a full wavelength is given by where θ represents the internal angle of the element (0 ≤ θ ≤ 180), and r (θ) refers to a rolling that occurs in the intersection of panels in membrane elements. Typically, the rolling of the corners is neglected, yet it is considered here due to the prevalence of the corners in a single cell column, which is expressed as where a 1 , a 2 , and a 3 are parameters intrinsic to the geometry, and B refers to a half of the wall length at the corners. Zhang and Zhang [29] found that these values match well with [0.163, 0.6, 0.06] for a corner element of 120 degrees. The Kelvin cell requires a primary angle of 135 degrees; however, these values must be resolved in this case through experimental means. As the change in corner angle for 120 to 135 degrees is relatively small, resulting in a negligible change in energy absorption for the corner elements [30], these values were left unchanged as they are related to geometry but not wall length or thickness. From this, a 1 was solved as specified by Zhang and Zhang [29] with the known predetermined a 2 and a 3 values for a sole corner element and found to be 0.173.
The membrane energy of the T intersection (type II) is determined by [31] The membrane energy of the 6-panel intersection element (type III) is given by [32] Similarly, the membrane energy of the 4-panel intersection (type IV) may use [30] The perfect 3-panel intersection of 120 degrees, denoted as type V, uses Eq. (14) for an approximation [31] and the more specialized 3-panel intersection, denoted as type VI, uses Eq. (15) for a primary equation [31] The total membrane energy is given by a summation of membrane energies from different constitutive elements where T refers to the number of membranes in the unit cell of concern, n and p refer to the structure in question and the membrane of that structure, respectively. The equilibrium equation is As the folds are not perfect due to geometric effects and buckling, this value is best approximated from 0.7 to 0.75 [31], and this study uses 0.73 as its value. Based on the stationary condition It can be substituted back into Eq. (15), and MCF can be resolved. The general equation for MCF of Kelvin and floret can be defined as where A and B are specific to each geometry, being 145.32 and 72.56 for Kelvin honeycomb and 181.9 and 92.28 for floret honeycomb, respectively. Similarly, the MCF formula for Weaire honeycomb is with A being 229.1, B being 384, and C being 371.7. Based on the comparison in Fig. 4, the MCFs obtained from the theoretical model are close to what the FE model yielded.
More specifically, the Weaire structure's theoretical approximation is 42.8 kN, which is at a 8.9% difference from the 46.6 kN gathered from the numerical simulation. Similarly, the differences for the Kelvin and floret cases are 2.5% and 1.7%, respectively, indicating the validity of such approximations to determine the MCF of the proposed honeycombs.

Parametric Study
In accordance with Eq. (20), t and λ are vital for the energy absorption performances of the proposed honeycombs, and as such, we carry out a parametric study to analyze the effects of these parameters in this section (Fig. 4). The range for wall length is limited to 5-20 mm since the model begins to exit the standard deformation mode for all cells beyond these values. The varying range of the thickness is between 0.3 and 1 mm. As t increases, all the crashworthiness criteria show an increasing trend in the proposed honeycombs, except for Kelvin (Fig. 4a) which approaches a maximum for SEA and CFE at 7.5 mm. When it comes to the variation of λ, an evident decrease in SEA and CEF is found for all the honeycombs, however, MCF is increasing as λ, indicating that the efficiency of energy absorption is closely dependent on the overall size of the structures. Meanwhile, for all unit cells, as either parameter increases the MCF increases, yet it occurs much faster with the increase in t than with the increase in λ. Furthermore, the MCFs obtained from the theoretical model are compared with the numerical results in Table 2, where large differences can be seen in unit cells near the extremes of each parameter (0.3 and 1mm for t and 0.5 and 20 mm for λ). This matches with the FE data beginning to transform the failure mechanisms (e.g., from plastic collapse to Euler buckling) at these extremes, which cannot be well described by the proposed theoretical model, and in turn, establishing the limits for the range of the parameters. As the parameters of thickness and wall length categorize the overall cross-sectional area of each cell, a third parameter of respective area, A/A o , or ratio of parameterized area to original area was taken [33]. As it is dependent on t and λ, wall lengths from the parameterized set were chosen and thicknesses were derived to minimize a change in resulting area from the baseline. This allows for an optimization of SEA and MCF based on the respective mass of the unit cell, further optimizing each cell for performance. This methodology can assist finding the most optimized structure for energy absorption characteristics; however, when characterizing differences across structures, the ratios of area should be considered between structures as well. Based on the above discussions, the comparison  between the geometries mentioned here based on A/A o cannot be conducted, nor can it be done for other known structures without further refining sizes of these cells. The benefits of the adapted structures are evident, but a closer inspection of the underlying mechanisms may yield more critical conclusions. The Kelvin cells are shown as 2D representations of a 3D Kelvin geometry, but can also be regarded as an overlapping of octagonal sections. In comparison with traditional honeycombs, it is shown that Kelvin under out-of-plane crushing has an improvement over octagonal honeycombs [25] but also yields limitations of these structures. The larger relative size of the octagonal honeycomb with the same side length of λ, as compared with hexagon, and the shape being not perfectly meshed with its respective cells, requiring square cells to fill gaps, may yield problems when attempting to provide more hierarchy to the cells.
The floret cells, while at first appearing unique, also can be viewed as an array of smaller hexagonal honeycombs of which a center hexagon is voided of sides and replaced with the 6-panel intersections. As it is known that smaller cell sizes, and therefore higher relative densities, have an increased effect on energy absorption capabilities, this may explain the beneficial properties of the floret design. Also, the Weaire structure is analogous to an array of hexagons of a side length λ with larger hexagons of a side length of √ 2λ. This relationship will only make sense to improve the energy absorption since the addition of the honeycombs also adds in cross-intersection membranes.

Conclusions
In this paper, we comprehensively investigated the dynamic behaviors of three types of bio-inspired honeycombs under out-of-plane crushing, including Kelvin, Weaire, and floret honeycombs. Compared to traditional honeycombs, superior energy absorption performances are achieved in these structures. Based on the deformation mechanisms, a theoretical model was developed to derive the analytical approximation of MCFs for each structure, which can accurately match the results from the finite element analysis. Parametric studies were also conducted to explore the effects of various geometric parameters on the energy absorption performances for these structures. Meanwhile, the limits of the theoretical models were analyzed and the underlying mechanisms for structural design were revealed, which may offer insights and broader design considerations for next-generation energy absorber in various areas.