**1. Introduction**

Re-entrant honeycomb structures that display negative Poisson's ratios (NPR) are known to be one class of auxetic structures and have been used in many fields, such as aerospace and automotive industries. The multifunctionality of anisotropic re-entrant honeycomb has been widely studied for its static mechanical behavior [1–4], dynamic performance [5,6], thermal conductivity and heat transfer properties [7].

Hierarchy [8] is one of the most readily observed topological features in natural structures and now has been introduced to honeycomb and chiral lattice structures in pursuing ultralight materials with improving elastic properties and damage tolerance. Specifically, considering hierarchical sub-structures to honeycombs and designing novel metamaterials with tailorable multi-functional properties, have attracted increasing attention in recent years.

Extreme values of hierarchical metamaterial properties such as specific stiffness, toughness, strength, negative or complex Poisson's ratio, zero or negative thermal expansion, phononic band gaps as well as impact energy absorption have been reported in hierarchical architectures across multiple length scales [9–15]. Sun et al. [16] analytically studied the in-plane elastic moduli and thermal conductivity of a multifunctional hierarchical honeycomb (MHH), which is formed by replacing the solid cell walls of an original regular hexagonal honeycomb (ORHH) with three different isotropic honeycomb sub-structures possessing hexagonal, triangular or kagome lattices. Then the anisotropic multifunctional hierarchical honeycomb (AMHH) with triangular or kagome honeycomb substructures (OAHH) was proposed and the in-plane stiffness of these two kinds of AMHH was analytically studied with the help of Euler beam theory [17]. Taylor et al. [18] investigated the in-plane elastic

properties and structural hierarchy in honeycombs and explored the effects of adding hierarchy into a range of honeycombs, with hexagonal, triangular or square geometry super and sub-structure cells by using finite element simulation. Key parameters describing these geometries included the relative lengths of the sub- and super-structures, the fraction of mass shared between the sub- and super-structures, the co-ordination number of the honeycomb cells, the form and extent of functional grading, and the Poisson's ratio of the sub-structure. Mousanezhad et al. [19] studied the effects of chirality and hierarchy on elastic response of honeycombs, derived the closed-form expressions for elastic moduli of several chiral, anti-chiral and hierarchical honeycombs with hexagon and square based networks, and finally validated the analytical estimates of the elastic moduli by using finite element method. Gatt et al. [20] proposed a new class of hierarchical auxetics based on the rotating rigid unit mechanism. These systems retain the enhanced properties from having a negative Poisson's ratio with the added benefits of being a hierarchical system. Through design, one can control the extent of auxeticity, degree of aperture and size of the different pores in the system, which makes the system more versatile than similar non-hierarchical ones. Chen et al. [21] reported a group of hierarchically architected metamaterials constructed by replacing cell walls of regular honeycombs with hexagonal, kagome, and triangular lattices, respectively. The numerical and analytical studies indicate that the introduction of structural hierarchy in regular honeycombs results in improved heat resistance and thermal anisotropy. Then, Yin et al. [22] studied the in-plane crashworthiness of the hierarchical honeycomb group above, using the nonlinear finite element code LS-DYNA. The numerical simulation results indicate that the triangular hierarchical honeycomb provides the best performance compared to the other two hierarchical honeycombs and features more than twice the energy absorbed by the regular honeycomb under similar loading conditions. More recently, Wu et al. [23,24] proposed an innovative hierarchical anti-tetrachiral structure as well as a hierarchical anti-tetrachiral stent with circular and elliptical nodes, based on the auxetic deformation behaviors of anti-tetrachiral unit cell at different structural hierarchical levels. It was found that the mechanical behaviors of hierarchical anti-tetrachiral structure can be tailored through adjusting the levels of hierarchical structures and unit cell design, and the proposed hierarchical anti-tetrachiral stents exhibit remarkable radial expanding abilities while maintaining axial stability. Besides the above-mentioned hierarchical metamaterials, various types of alternative hierarchical structures have been proposed through the modification of the node or cell wall structural levels for generating enhanced and tunable mechanical properties through structural hierarchy approaches.

In the current work, a novel auxetic hierarchical metamaterial was designed, which consisted of a re-entrant representative unit cell as well as a star subordinate cell with zero Poisson's ratio; both the main and sub cells were planar symmetric. Full dimensional models of the new hierarchical metamaterial to describe the effective elasticity as well as loading-bearing capability in plane were simulated by finite element method (FEM). Comprehensive parametric studies for both parent RUC and subordinate cell were performed to evaluate in-plane non-dimensional moduli, effective Poisson's ratio and relative density of the new hierarchical structure, corresponding optimum structure-property relations were explored for the designed metamaterial.

### **2. Materials and Methods**

The geometry of the novel hierarchical star-re-entrant metamaterial is presented in Figure 1. The representative unit cell (RUC), composed of star subordinate cell can be described as a re-entrant hexagonal structure, and apparently both the cells are horizontal and vertical symmetry. Analogy to a zero Poisson's ratio cellular structure [25], the star subordinate cell is represented by length *L*0, inclined angle θ as well as in-plane thickness *t*, respectively (Figure 1a). Then dimensions of the new hierarchical star-re-entrant RUC is described by length *Lx*, *Ly* and symmetric inclined angle 45◦ (Figure 1b), where *Li* = *Ni* × *L*0 (*i* = *x*, *y*), *Ni* is the number of star subordinate cells along *x* or *y* direction. The detail of star joint in re-entrant structure is presented in Figure 1c, the length of jointing star is expressed as *W* = ( √2*t*)/(4 sin θ cos(45◦ − <sup>θ</sup>)). To avoiding the overlapping and contact of cell walls, geometric constraints 0◦ <θ< 45◦ and 0 < *t* < *<sup>L</sup>*0(cos θ − sin θ) need to be satisfied.

**Figure 1.** Geometry of the novel star-re-entrant hierarchical metamaterial: (**a**) star subordinate cell, (**b**) re-entrant representative unit cell, (**c**) the jointing between two neighboring stars.

The relative density is an important parameter for cellular structures, and it is defined by:

$$\frac{\rho}{\rho\_c} = \frac{A}{A\_c} \tag{1}$$

where *A* and *Ac* are respectively the cross-sections perpendicular to the out-plane thickness direction and the load bearing area. Here, for this new hierarchical RUC, relative density results in:

$$\frac{\rho}{\rho\_c} = \frac{36L\_0 t \sin 3\theta - 36L\_0 t \sin \theta - 10t^2 \sin 2\theta - 8t^2 \cos 2\theta - 8t^2}{9L\_0^2 \sin 4\theta} \tag{2}$$

Finite element analysis was performed with Abaqus/CAE 6.13–4 commercial package standard for one RUC as well as the whole hierarchical metamaterial structure. For all the simulations in this paper, geometric dimensions of the new RUC are defined as *L*0 = 20, θ = 30◦, *t* = 2 and *Lx* = *Ly* = 6*L*0 with symmetric inclined angle 45◦, respectively. Acrylonitrile butadiene styrene (ABS) plastic with a rapid prototyping Fusion Deposition Molding (FDM) Stratasys machine was used to manufacture all the experimental samples and the elastic mechanical properties of the core material for finite element simulation were set as *Ec* = 2265MPa and Poisson's ratio ν*c* = 0.25 [25,26]. An elastic shell element with reduced integration (S4R) and element size of 0.8 for convergence were chosen for the simulation of in-plane effective moduli, shown in Figure 2a. Accounting for the symmetry, the moduli of elasticity and the Poisson's ratio were determined by one quarter of the metamaterial structure. Taking full-size representative volumes with 6 × 6 cells for example, boundary conditions for in-plane tensile Young's moduli and Poisson's ratio were established based on References [23,27], where nodes on the left and bottom edge were constrained from out-plane rotation and translation normal to the edge direction, respectively, and displacements in the *x*-direction (*y*-direction) were applied to the ligament nodes on the right (top) edge, which was also constrained from in-plane rotation, shown in Figure 2b. In the case of the in-plane shear simulation, biaxial loading was introduced as close as possible a pure shear deformation field [27,28], corresponding boundary conditions above as well as displacements in both *x* and *y* directions were applied and are presented in Figure 2c.

Nominal strain and stress of the hierarchical representative unit cell in *i*( = *x* or *y*) direction were calculated from:

$$
\sigma\_i = \frac{\delta\_i}{L\_i}, \sigma\_i = \frac{F\_i}{A\_j} \tag{3}
$$

where δ*i* is the applied displacement, *Fi* is the sum of the nodal reaction forces on the edge to which displacement was applied, *Li* and *Aj* are the initial length and cross-sectional area of the hierarchical structure in the *i* and *j* (= *y* or *x*) directions, respectively. According to Equation (3), effective Young's modulus as well as Poisson's ratio are calculated by:

$$E\_i = \frac{\sigma\_i}{\varepsilon\_i}, \nu\_{ij} = -\frac{\delta\_j L\_i}{\delta\_i L\_j} \tag{4}$$

where, *i* is the loading direction.

**Figure 2.** Numerical model description: (**a**) Mesh of one representative unit cell, (**b**) boundary conditions of axial tension along x direction, (**c**) boundary conditions of biaxial shear test.

Then, effective shear modulus under biaxial loading in plane was obtained as following [27,28]:

$$G\_{xy} = \frac{\tau}{\gamma} = \frac{\sigma\_x - \sigma\_y}{2(\varepsilon\_x - \varepsilon\_y)} = \frac{R\_x L\_x - R\_y L\_y}{2h(\delta\_x L\_y - \delta\_y L\_x)}\tag{5}$$

In Equation (5), *Rx* and *Ry* are reaction force along *x* and *y* direction, *h* is the out-plane thickness of hierarchical structure.

In order to highlight the influence of the cell numbers on the convergence of the results, computations were undertaken starting by a number of 2 × 2 cells to a maximum of 40 × 40 cells. The convergence was found to be achieved at the number of 40 cells; corresponding dependence of mechanical property on the computations number of cells as well as the results of effective non-dimensional moduli and Poisson's ratio are demonstrated in Figure 3, respectively.

**Figure 3.** Finite element simulations of effective mechanical property for the new hierarchical metamaterial: (**a**) *<sup>E</sup>*<sup>∗</sup>1/*Ec* and <sup>ν</sup><sup>∗</sup>12, (**b**) *<sup>E</sup>*<sup>∗</sup>2/*Ec* and <sup>ν</sup><sup>∗</sup>21, (**c**) *<sup>G</sup>*<sup>∗</sup>12/*Ec*.

### **3. Results and Discussion**

To understand how the geometrical parameters of the star-re-entrant RUC influence the e ffective mechanical properties of the new hierarchical metamaterial designed, parametric studies were conducted by using finite element models described in Section 3, and the numerical results are presented and discussed as follows.

### *3.1. The Geometry E*ff*ects of Star Subordinate Cell*

Figures 4–6 demonstrate the FE homogenization of the non-dimensional in-plane elastic moduli and corresponding Poisson's ratio versus various parameters of lengths *L*0, thickness *t* and cell inclined angle θ. In general, all the non-dimensional elastic moduli decrease with increasing *L*0, shown in Figure 4 and increase with increasing thickness *t*, seen in Figure 5, when the other geometrical parameters keep constant. The variations of non-dimensional e ffective moduli with cell inclined angle θ are presented in Figure 6. For the increase of θ, *E*∗ 1/*Ec* exhibits an up-down-up trend, while the other two non-dimensional moduli display a first descent and then ascent with di fferent gradients, for which *E*∗ 2/*Ec* and *<sup>G</sup>*<sup>∗</sup>12/*Ec* increased by 9.89% and 35.9%, respectively. The impact of geometric parameters on e ffective Poisson's ratio of the new hierarchical metamaterial are discussed as following. Poisson's ratio <sup>ν</sup><sup>∗</sup>12 only increases with increasing *L*0 (Figure 4a) and declines with thickness *t* and cell inclined angle θ, presented in Figures 5a and 6a, respectively. While the variations of Poisson's ratio <sup>ν</sup><sup>∗</sup>21 with rising *L*0, θ and *t* exhibit as constant (Figure 4b), up-down (Figure 5b) as well as ascending (Figure 6b), separately.

**Figure 4.** The variations of effective mechanical properties verse length *L*0: (**a**) *E*∗ 1/*Ec* and <sup>ν</sup><sup>∗</sup>12, (**b**) *E*∗ 2/*Ec* and <sup>ν</sup><sup>∗</sup>21, (**c**) *<sup>G</sup>*<sup>∗</sup>12/*Ec*.

**Figure 5.** The variations of effective mechanical properties verse thickness *t*: (**a**) *<sup>E</sup>*<sup>∗</sup>1/*Ec* and <sup>ν</sup><sup>∗</sup>12, (**b**) *<sup>E</sup>*<sup>∗</sup>2/*Ec* and <sup>ν</sup><sup>∗</sup>21, (**c**) *<sup>G</sup>*<sup>∗</sup>12/*Ec*.

**Figure 6.** The variations of effective mechanical properties verse inclined angle θ: (**a**) *<sup>E</sup>*<sup>∗</sup>1/*Ec* and <sup>ν</sup><sup>∗</sup>12, (**b**) *<sup>E</sup>*<sup>∗</sup>2/*Ec* and <sup>ν</sup><sup>∗</sup>21, (**c**) *<sup>G</sup>*<sup>∗</sup>12/*Ec*.

In view of increasing cell thickness *t* mainly enhancing the structural weight and stiffening mechanical behaviors, therefore, the effects of cell inclined angle θ for different *L*0 on in-plane

mechanical property were studied in detail and are demonstrated in Figure 7. It may be observed that variation and magnitude of the effective mechanical properties for *L*0 = 20 with enhancive θ are completely different from those of *L*0 = 40 − 120. With the increasing θ, non-dimensional modulus *<sup>E</sup>*<sup>∗</sup>1/*Ec* displays an increasing and decreasing trend for *L*0 = 40, 60, and exhibits a gradual decrease for *L*0 = 80 − 120, presented in Figure 7a; *<sup>E</sup>*<sup>∗</sup>2/*Ec* and *<sup>G</sup>*<sup>∗</sup>12/*Ec* are both observed declines with slower slopes, for *L*0 = 40 to *L*0 = 120, shown in Figure 7b,c, respectively. For the study of Poisson's ratio, the increasing cell angle θ makes <sup>ν</sup><sup>∗</sup>12 decrease and <sup>ν</sup><sup>∗</sup>21 increase inversely, seen in Figure 7d,e. Additionally, parameter *L*0 makes no apparent effect on Poisson's ratio, except the scenarios of *L*0 = 20, where <sup>ν</sup><sup>∗</sup>12 remains constant in the range of θ = 35◦ to θ = 40◦ as well as <sup>ν</sup><sup>∗</sup>21 behaviors a relatively grea<sup>t</sup> variation from −3.456 with θ = 10◦ to −1.271 with θ = 40◦.

**Figure 7.** Effects of parameters *L*0 and θ on the effective mechanical property for *t* = 2: (**a**) *<sup>E</sup>*<sup>∗</sup>1/*Ec*, (**b**) *<sup>E</sup>*<sup>∗</sup>2/*Ec*, (**c**) *<sup>G</sup>*<sup>∗</sup>12/*Ec*, (**d**) <sup>ν</sup><sup>∗</sup>12, (**e**) <sup>ν</sup><sup>∗</sup>21.

### *3.2. The E*ff*ects of Subordinate Cell Amount*

The amount of star subordinate cells along *x* and *y* directions are defined and presented in Figure 8, where *Nx* is the number of half star subordinate cell along *x* direction and *Ny* is the number of entire star subordinate cell along *y* direction. The effects of subordinate cell amount on effective mechanical behavior were then studied. It may be clearly observed from Figure 9a that a growing amount of *Nx* increases *<sup>E</sup>*<sup>∗</sup>1/*Ec* and <sup>ν</sup><sup>∗</sup>12, when *Nx* ≥ 7, <sup>ν</sup><sup>∗</sup>12 turns positive. However, increasing *Nx* makes both *<sup>E</sup>*<sup>∗</sup>2/*Ec* and *<sup>G</sup>*<sup>∗</sup>12/*Ec* decline and remains <sup>ν</sup><sup>∗</sup>21 to be constant, presented in Figure 9b,c, respectively. Figure 10 shows how the effective properties of the new metamaterials vary with subordinate cell amount along y direction. Non-dimensional effective modulus *<sup>E</sup>*<sup>∗</sup>1/*Ec* and *<sup>G</sup>*<sup>∗</sup>12/*Ec* decrease with the increasing number of *Ny*, Poisson's ratio <sup>ν</sup><sup>∗</sup>12 increases and stays auxetic, whereas, *<sup>E</sup>*<sup>∗</sup>2/*Ec* and <sup>ν</sup><sup>∗</sup>21 exhibit the opposite variations under an increasing number of *Ny*. Consequently, it is found that a small amount of *Nx* and appropriate number of *Ny* can satisfy the new hierarchical metamaterial with wholly auxetic behavior and strong stiffness in-plane simultaneously.

**Figure 8.** The number variations of star subordinate cell: (**a**) *Nx* along *x* direction, (**b**) *Ny* along *y* direction.

**Figure 9.** Variation of subordinate cell amount *Nx* verse: (**a**) *<sup>E</sup>*<sup>∗</sup>1/*Ec* and Poisson's ratio <sup>ν</sup><sup>∗</sup>12, (**b**) *<sup>E</sup>*<sup>∗</sup>2/*Ec* and Poisson's ratio <sup>ν</sup><sup>∗</sup>21, (**c**) *<sup>G</sup>*<sup>∗</sup>12/*Ec*.

### *3.3. Relative Density Study*

The formula of the relative density for one representative unit cell is given by Equation (2); the variation of relative density ρ/ρ*c* with parameters *L*0,θ and *t* were obtained and are shown in Figure 11a–c. For Figure 11d, parameter θ was valid in a range of 2.5◦ − 42.5◦, corresponding ρ/ρ*c* exhibits a non-monotonic going up and down variation. Therefore, the values of θ in the range of 35◦ to 40◦ made a different influence on the results of relative density as well as previous effective mechanical properties shown in Figure 7.

**Figure 10.** Variation of subordinate cell amount *Ny* verse: (**a**) *<sup>E</sup>*<sup>∗</sup>1/*Ec* and Poisson's ratio <sup>ν</sup><sup>∗</sup>12, (**b**) *<sup>E</sup>*<sup>∗</sup>2/*Ec* and Poisson's ratio <sup>ν</sup><sup>∗</sup>21, (**c**) *<sup>G</sup>*<sup>∗</sup>12/*Ec*.

**Figure 11.** *Cont*.

**Figure 11.** Geometric parameters effects on relative density of the hierarchical RUC: (**a**) *t* = 2, (**b**) θ = 30◦, (**c**) *L*0 = 20, (**d**) *L*0 = 20 and *t* = 2.

The relationship between effective mechanical properties and relative density ρ/ρ*c* simulated by FEM finite element method were investigated as following. Variations of specific stiffness (*E*<sup>∗</sup>i/*Ec*)/(ρ/ρ*c*) with different parameters are represented in Figures 12 and 13, respectively. In Figure 12a, specific stiffness all decline with increasing *L*0, extremum values of (*E*<sup>∗</sup>2/*Ec*)/(ρ/ρ*c*) = 0.05963 and (*<sup>G</sup>*<sup>∗</sup>12/*Ec*)/(ρ/ρ*c*) = 0.00563 were achieved with *L*0 = 20, θ= 30◦ and *t* = 2. Similarly, for the increase of cell thickness *t*, specific stiffnesses all exhibit a growing variation with different gradient, seen in Figure 12b. The impact of parameter θ on effective mechanical behavior of the proposed new metamaterial were investigated and are represented in Figure 12c. For the increase of θ, it may be observed that (*E*<sup>∗</sup>2/*Ec*)/(ρ/ρ*c*) and (*<sup>G</sup>*<sup>∗</sup>12/*Ec*)/(ρ/ρ*c*) have a resembling variation of first decline and then ascent, while (*E*<sup>∗</sup>1/*Ec*)/(ρ/ρ*c*) exhibit an up-down-up variation. In the range of θ= 35◦ to θ= 40◦, specific stiffnesses (*E*<sup>∗</sup>i/*Ec*)/(ρ/ρ*c*) all increase monotonically.

**Figure 12.** Variation of specific stiffness (*E*<sup>∗</sup>i/*Ec*)/(ρ/ρ*c*) with parameters: (**a**) *L*0, (**b**) *t*, (**c**) θ.

Figure 13 presents the variation of in-plane specific stiffness with the cell angle θ for various parameter *L*0, while *t* = 2. The specific stiffness along 1-direction exhibits three different variational trends with a rising θ for different *L*0, respectively, shown in Figure 13a: (1) for *L*0 = 20, (*E*<sup>∗</sup>1/*Ec*)/(ρ/ρ*c*) varies as an up-down-up curve and reaches the maximum value of 0.04829 with θ = 40◦, which is 25 times greater than the one for *L*0= 120; (2) for *L*0 = 40, 60,(*E*<sup>∗</sup>1/*Ec*)/(ρ/ρ*c*) presents a first increasing and then decreasing variation; (3) for *L*0 ≥ 80, (*E*<sup>∗</sup>1/*Ec*)/(ρ/ρ*c*) decreases monotonically. From Figure 13b,c, it can be observed that (*E*<sup>∗</sup>2/*Ec*)/(ρ/ρ*c*) and (*<sup>G</sup>*<sup>∗</sup>12/*Ec*)/(ρ/ρ*c*) decrease clearly and then ascend with the variation of cell angle θ from 35◦ to 40◦ when *L*0 = 20, a comparison of the specific stiffness in this range shows that (*E*<sup>∗</sup>2/*Ec*)/(ρ/ρ*c*) varies slightly from 0.03817 to 0.04422, however, (*<sup>G</sup>*<sup>∗</sup>12/*Ec*)/(ρ/ρ*c*) increases significantly from 0.0049 to 0.00704. When *L*0 ≥ 40, (*E*<sup>∗</sup>2/*Ec*)/(ρ/ρ*c*) and (*<sup>G</sup>*<sup>∗</sup>12/*Ec*)/(ρ/ρ*c*) both monotonically decline with increasing θ and *L*0. Figure 13c illustrates that an increase of more than 143% of the specific shear stiffness show up when the parameter *L*0 varies from 20 to 120 with θ = 40◦, which makes varying *L*0 also a good design method for (*<sup>G</sup>*<sup>∗</sup>12/*Ec*)/(ρ/ρ*c*). Therefore, it can be determined that all the high specific stiffness in plane can be achieved simultaneously by choosing proper parameter θ and *L*0.

**Figure 13.** Variations of specific stiffness (*E*<sup>∗</sup>i/*Ec*)/(ρ/ρ*c*) verse parameters *L*0 and θ, when *t* = 2: (**a**) (*E*<sup>∗</sup>1/*Ec*)/(ρ/ρ*c*), (**b**) (*E*<sup>∗</sup>2/*Ec*)/(ρ/ρ*c*), (**c**) (*G*<sup>∗</sup>12/*Ec*)/(ρ/ρ*c*).

Finally, the ratio between effective Poisson's ratio and relative density verse parameters *L*0 and θ were investigated. In Figure 14a, it is seen that <sup>υ</sup><sup>∗</sup>12/(ρ/ρ*c*) performs from positive to remarkable auxetic behavior and it declines significantly with both increasing *L*0 and θ, the maximal descending slope is achieved with *L*0= 120. Comparing with Figure 13a, the optimum values of both high specific stiffness and large auxetic deformation in 1-direction can be selected widely for the special curve shape of (*E*<sup>∗</sup>1/*Ec*)/(ρ/ρ*c*) with *L*0 = 20. Figure 14b reveals that the auxetic <sup>υ</sup><sup>∗</sup>21/(ρ/ρ*c*) increases with greater θ and decreases with increasing *L*0, oppositely. Contrast to Figure 13b, the optimal values of both high specific stiffness and large auxetic deformation in 2-direction is acquired for θ = 10◦ and *L*0 = 20, where (*E*<sup>∗</sup>2/*Ec*)/(ρ/ρ*c*)= 0.57228 and <sup>υ</sup><sup>∗</sup>21= −3.45649.

**Figure 14.** Variations of the ratio between effective Poisson's ratio and relative density verse parameters *L*0 and θ, when *t* = 2: (**a**) <sup>υ</sup><sup>∗</sup>12/(ρ/ρ*c*), (**b**) <sup>υ</sup><sup>∗</sup>21/(ρ/ρ*c*).
