Dynamic Crushing Analysis of a Three-Dimensional Re-Entrant Auxetic Cellular Structure

Dynamic behaviors of the three-dimensional re-entrant auxetic cellular structure have been investigated by performing beam-based crushing simulation. Detailed deformation process subjected to various crushing velocities has been described, where three specific crushing modes have been identified with respect to the crushing velocity and the relative density. The crushing strength of the 3D re-entrant auxetic structure reveals to increase with increasing crushing velocity and relative density. Moreover, an analytical formula of dynamic plateau stress has been deduced, which has been validated to present theoretical predictions agreeing well with simulation results. By establishing an analytical model, the role of relative density on the energy absorption capacity of the 3D re-entrant auxetic structure has been further studied. The results indicate that the specific plastic energy dissipation is increased by increasing the relative density, while the normalized plastic energy dissipation has an opposite sensitivity to the relative density when the crushing velocity exceeds the critical transition velocity. The study in this work can provide insights into the dynamic property of the 3D re-entrant auxetic structure and provides an extensive reference for the crushing resistance design of the auxetic structure.


Introduction
Auxetic cellular structures have attracted great attention since, unlike conventional honeycomb, they expand laterally when stretched and contract laterally when compressed, featuring a negative Poisson's ratio (NPR). The termed 'auxetic effect' brings comprehensive improvements of mechanical properties, like increased stiffness, high shear modulus, indentation resistance and energy absorbing performance, to the cellular structure. These fascinating properties all contribute to augmenting their applications in automotive, packaging and structural protection [1][2][3]. In recent years, increasing interests have been devoted in this field and the research results can be summarized to cover the topology design, property investigation, rapid prototyping and utilization of the auxetic structures [4][5][6][7][8].
Due to superior performance in dissipating energy, auxetic cellular structures can work as a better alternative in the design of sandwich structure and vehicle crash box, which illustrates the significance of dynamic crushing analysis of auxetic structures [9,10]. To investigate the dynamic properties of auxetic cellular structures, a few studies of theoretical models and experiments have been conducted and some remarkable conclusions have been demonstrated. Grujicic et al. [11] have investigated the dynamic response of the auxetic sandwich structures by using a multiphysics computational The materials of the struts may be varied to achieve different properties. As a preferred embodiment, they are assumed to have the same material in this study. Then, as shown in Figure 1, the relative density of the 3D re-entrant auxetic structure, , that is, the ratio between the effective density of cellular structure and the density of strut material, can be obtained by dividing the volume occupied by all the struts and the total volume. The volume occupied by all the struts is specified as = 2 2 ⁄ + 4 2 ⁄ + 8 . The total volume is calculated as = (2 sin + 2 ) • (2 − 2 cos ). Combining the independent parameters defined above, the relative density is therefore described as:

Geometry and Modeling
2.1. Geometry Figure 1a depicts the employed 3D re-entrant auxetic cellular structure, which is arrayed by repeating the representative volume element (RVE) shown in Figure 1b in all three principle directions. As can be seen from Figure 1b, the RVE can be described with a pair of orthogonally arranged re-entrant hexagons, which consist of several vertical and inclined struts. By referring to the previous publications, the struts are assumed to have the identical square-shape cross-sections. Figure 1c illustrates the whole view of the geometrical parameters for characterizing the RVE. The lengths of the inclined struts and vertical struts are defined as L and H, respectively. The thickness of the strut cross-section is defined as t. The re-entrant angle between the vertical and inclined struts is defined as θ. Therefore, the geometry of the RVE can be determined independently with four parameters namely L, H, t and θ. Noting that to ensure re-entrant configuration, a geometric constraint H > 2L cos θ should be satisfied in initial design. Considering the symmetry of the 3D structure, it is apparent that the mechanical properties along the x and y directions are of the same.
The materials of the struts may be varied to achieve different properties. As a preferred embodiment, they are assumed to have the same material in this study. Then, as shown in Figure 1, the relative density of the 3D re-entrant auxetic structure, ρ R , that is, the ratio between ρ 0 the effective density of cellular structure and ρ s the density of strut material, can be obtained by dividing the volume occupied by all the struts and the total volume. The volume occupied by all the struts is specified as V s = 2Ht 2 /2 + 4Ht 2 /2 + 8Lt 2 . The total volume is calculated as V 0 = (2L sin θ + 2t) 2 ·(2H − 2L cos θ). Combining the independent parameters defined above, the relative density ρ R is therefore described as:

Description of Finite Element Modeling
In this section, a full size finite element (FE) model of 3D re-entrant auxetic cellular structure for dynamic crushing simulation has been built using Abaqus ® , as shown in Figure 2a. The 3D re-entrant auxetic structure adopted for simulation have been given with 12 × 9 × 6 representative cells (twelve unit cells in x-direction, nine unit cells in y-direction and six unit cells in z-direction). The geometrical parameters of the 3D specimens in crushing analysis are set as L = 20 mm, H = 30 mm, θ = 60 • and t varies from 1 mm to 4 mm, featuring a variable relative density closely between 1% and 10%. The dimensions of whole model are then given as 416 mm by 360 mm by 208 mm in Figure 2a.
In meshing of the model, by referring to the publication [22], the two-node beam element BEAM188 has been used to mesh the struts of the structure as this element type can ensure high accuracy at a low computational cost. As demonstrated in the paper [22], this element type is suitable to simulate the linear small-deflection and non-linear large-deformation of the strut, which both occur in dynamic crushing. Figure 2a also shows an example of the FE beam model of a half unit cell. As can be observed from the figure, the mesh model is constructed based on the design scheme shown in Figure 1c (red dash), as the schematic lines of beam element are directly assigned along the central axes of the struts. And all the corners of struts are defined as default nodes before meshing. As a preferred material for additive manufacturing, PA6 has been adopted as the matrix material of strut and been assumed to be a rate-independent elastic-perfectly plastic model with mass density ρ s = 1.14 g/cm 3 , Young's modulus E s = 2.32 GPa, yield stress σ y = 115 MPa and Poisson's ratio ν s = 0.34 [23].
The detailed boundary conditions for crushing simulations are depicted in Figure 2b. To simulate the dynamic crushing process, a rigid-wall has been modeled to impact the upper cells of the specimen with a constant velocity, v c , along the y-direction, as shown in Figure 2b. The crushing velocity v c varies from 2 m/s to 50 m/s in this study. All the nodes at the bottom boundary have been fixed along the y-direction as u y = 0. To fix the model, a central node at the bottom boundary has been fixed by constraining all degrees of freedom. Nodes at the lateral boundaries have been set with the coupled boundaries so as to avoid the influence of boundary effect [24]. An automatic single contact algorithm in Abaqus ® has been adopted to model the contact between surfaces with a friction coefficient of 0.35. contact algorithm in Abaqus ® has been adopted to model the contact between surfaces with a friction coefficient of 0.35. A mesh sensitivity investigation has been carried out to determine the proper mesh arrangement. The basic element size in meshing the strut varies from 4 mm to 10 mm for different model complexities. The different models are investigated under crushing at a constant velocity of 30 m/s. The convergence of the simulations is depicted using crushing reaction force in Figure 3a and hourglass error in Figure 3b. Abaqus explicit analysis performed in the present study uses hourglass control mode. Once the artificial strain energy is greater than approximately 5% of the total energy, that is, hourglass error exceeds 5%, hourglassing could be a problem and should be further analyzed by refining mesh [25]. It can be concluded from Figure 3 that the hourglass error of meshing with basic element size of 6 mm can be limited in 5% (in Figure 3b) and is able to achieve a similar result compared with the other A mesh sensitivity investigation has been carried out to determine the proper mesh arrangement. The basic element size in meshing the strut varies from 4 mm to 10 mm for different model complexities. The different models are investigated under crushing at a constant velocity v c of 30 m/s. The convergence of the simulations is depicted using crushing reaction force in Figure 3a and hourglass error in Figure 3b. contact algorithm in Abaqus ® has been adopted to model the contact between surfaces with a friction coefficient of 0.35. A mesh sensitivity investigation has been carried out to determine the proper mesh arrangement. The basic element size in meshing the strut varies from 4 mm to 10 mm for different model complexities. The different models are investigated under crushing at a constant velocity of 30 m/s. The convergence of the simulations is depicted using crushing reaction force in Figure 3a and hourglass error in Figure 3b. Abaqus explicit analysis performed in the present study uses hourglass control mode. Once the artificial strain energy is greater than approximately 5% of the total energy, that is, hourglass error exceeds 5%, hourglassing could be a problem and should be further analyzed by refining mesh [25]. It can be concluded from Figure 3 that the hourglass error of meshing with basic element size of 6 mm can be limited in 5% (in Figure 3b) and is able to achieve a similar result compared with the other Abaqus explicit analysis performed in the present study uses hourglass control mode. Once the artificial strain energy is greater than approximately 5% of the total energy, that is, hourglass error exceeds 5%, hourglassing could be a problem and should be further analyzed by refining mesh [25]. It can be concluded from Figure 3 that the hourglass error of meshing with basic element size of 6 mm can be limited in 5% (in Figure 3b) and is able to achieve a similar result compared with the other refined models (in Figure 3a). Thus, combining the accuracy and computational cost, the basic element size in meshing the model has been set as 6 mm for all simulation in the present work.

Dynamic Crushing Modes
As indicated in previous publications [26][27][28], the crushing modes for regular hexagonal honeycombs with same relative density are highly dependent on the crushing velocity. Inspired by this work, keeping a constant ρ R = 1%, dynamic crushing simulation of the 3D re-entrant auxetic cellular structure has been performed under a wide range of crushing velocity v c varying from 2 m/s to 50 m/s. By defining the y-directional strain rate as the ratio of v c to the height of the structure h y , the strain rate varies approximately from 5 s −1 to 120 s −1 in this study. Figures 4-7 depict the deformation processes of the auxetic structures subjected to dynamic crushing along the y-direction. As seen from the simulation results, similar to the studies for regular hexagonal honeycombs [26], three specific crushing modes have also been observed from the numerical simulations for auxetic structures, namely quasi-static, transition and dynamic mode. But due to the auxetic effect, these specific crushing modes all exhibit different deformation shapes compared to the crushing patterns for regular hexagonal honeycombs reported as 'X,' 'V' and 'I' types [26]. As shown in Figure 4, when the auxetic structure is crushed under a low velocity (v c = 5 m/s), localized deformation occurs as the cells collapse layer by layer and obvious contraction is observed along the lateral plane. The first collapse band is initiated randomly at the 5th layer. As the crushing continues, the cells near the proximal end gradually pile up towards the collapse band, forming a densification core. Until all the cells above the collapse band densified, the densification core keeps spreading towards the distal end. No 'X' or 'V' type modes as for regular hexagonal honeycombs are observed in crushing case for the auxetic structure. Since the lateral contraction around the densification core is larger than those at the proximal and distal ends, the compressed structure shows a necking shape, which is seen as the diagnostic feature of quasi-static mode. Moreover, simulations for low-velocity cases show an instability in the growing of densification core, as shown in Figure 4. The possible cause for this phenomenon can be attributed to the reduction of cross-section area in densification core, leading to a dislocation between the densified layers. refined models (in Figure 3a). Thus, combining the accuracy and computational cost, the basic element size in meshing the model has been set as 6 mm for all simulation in the present work.

Dynamic Crushing Modes
As indicated in previous publications [26][27][28], the crushing modes for regular hexagonal honeycombs with same relative density are highly dependent on the crushing velocity. Inspired by this work, keeping a constant = 1%, dynamic crushing simulation of the 3D re-entrant auxetic cellular structure has been performed under a wide range of crushing velocity varying from 2 m/s to 50 m/s. By defining the y-directional strain rate as the ratio of to the height of the structure ℎ , the strain rate varies approximately from 5 s −1 to 120 s −1 in this study. Figures 4-7 depict the deformation processes of the auxetic structures subjected to dynamic crushing along the y-direction. As seen from the simulation results, similar to the studies for regular hexagonal honeycombs [26], three specific crushing modes have also been observed from the numerical simulations for auxetic structures, namely quasi-static, transition and dynamic mode. But due to the auxetic effect, these specific crushing modes all exhibit different deformation shapes compared to the crushing patterns for regular hexagonal honeycombs reported as 'X,' 'V' and 'I' types [26]. As shown in Figure 4, when the auxetic structure is crushed under a low velocity ( = 5 m/s), localized deformation occurs as the cells collapse layer by layer and obvious contraction is observed along the lateral plane. The first collapse band is initiated randomly at the 5th layer. As the crushing continues, the cells near the proximal end gradually pile up towards the collapse band, forming a densification core. Until all the cells above the collapse band densified, the densification core keeps spreading towards the distal end. No 'X' or 'V' type modes as for regular hexagonal honeycombs are observed in crushing case for the auxetic structure. Since the lateral contraction around the densification core is larger than those at the proximal and distal ends, the compressed structure shows a necking shape, which is seen as the diagnostic feature of quasi-static mode. Moreover, simulations for low-velocity cases show an instability in the growing of densification core, as shown in Figure 4. The possible cause for this phenomenon can be attributed to the reduction of cross-section area in densification core, leading to a dislocation between the densified layers. When the crushing velocity is increased ( = 15 m/s in Figure 5), the auxetic structure deforms in the same way as in the low-velocity crushing cases through accumulation of collapse layers into the densification core. But as can be seen from the figure, there is an obvious difference that the initial collapse band is generated in the layer close to the proximal end. Additionally, the lateral contractions at the proximal and distal ends almost occur simultaneously. As the crushing continues, the growing When the crushing velocity is increased (v c = 15 m/s in Figure 5), the auxetic structure deforms in the same way as in the low-velocity crushing cases through accumulation of collapse layers into the densification core. But as can be seen from the figure, there is an obvious difference that the initial collapse band is generated in the layer close to the proximal end. Additionally, the lateral contractions at the proximal and distal ends almost occur simultaneously. As the crushing continues, the growing of the densification core at the upper structure and the contraction of the both ends are the major deformation pattern for the auxetic structure, featuring a spindle shape for the transition mode.
When the auxetic structure is subjected to a high-velocity crushing (v c = 50 m/s in Figure 6), a shockwave like crushing front is generated at the proximal end only and then spreads towards the distal end, which is similar to the deformation pattern for regular hexagonal honeycombs in dynamic mode. Upon crushed, the cells near the proximal end rapidly collapse into a densification core and a localized contraction is also expected around the densification core, exhibiting a top-shrinking shape as a whole. In this crushing mode, no obvious collapse and contraction are observed near the distal end. This might be because that most of the crushing energy is absorbed through plastic collapse of cells in the densification core. Hence, the stress waves travelled to the distal ends are not intense enough to collapse the cells at this area [29]. This phenomenon can be further investigated by integrating with the crushing strength analysis given below.
The detailed crushing processes for auxetic structures with respect to the analyzed range of crushing velocities are shown in Figure 7, which clearly illustrates the gradual transition among three specific crushing modes. As described above, the main differences among the three crushing modes lie in the location of initial collapse band (marked in Figure 7 with red box) and the growth process of densification core. In addition to the crushing velocity, the crushing modes of auxetic structures are highly affected by the relative density. Therefore, auxetic structures with variable relative densities have been modeled by changing the thickness of strut and used to perform crushing simulations with respect to different crushing velocities. of the densification core at the upper structure and the contraction of the both ends are the major deformation pattern for the auxetic structure, featuring a spindle shape for the transition mode. When the auxetic structure is subjected to a high-velocity crushing ( = 50 m/s in Figure 6), a shockwave like crushing front is generated at the proximal end only and then spreads towards the distal end, which is similar to the deformation pattern for regular hexagonal honeycombs in dynamic mode. Upon crushed, the cells near the proximal end rapidly collapse into a densification core and a localized contraction is also expected around the densification core, exhibiting a top-shrinking shape as a whole. In this crushing mode, no obvious collapse and contraction are observed near the distal end. This might be because that most of the crushing energy is absorbed through plastic collapse of cells in the densification core. Hence, the stress waves travelled to the distal ends are not intense enough to collapse the cells at this area [29]. This phenomenon can be further investigated by integrating with the crushing strength analysis given below. of the densification core at the upper structure and the contraction of the both ends are the major deformation pattern for the auxetic structure, featuring a spindle shape for the transition mode. When the auxetic structure is subjected to a high-velocity crushing ( = 50 m/s in Figure 6), a shockwave like crushing front is generated at the proximal end only and then spreads towards the distal end, which is similar to the deformation pattern for regular hexagonal honeycombs in dynamic mode. Upon crushed, the cells near the proximal end rapidly collapse into a densification core and a localized contraction is also expected around the densification core, exhibiting a top-shrinking shape as a whole. In this crushing mode, no obvious collapse and contraction are observed near the distal end. This might be because that most of the crushing energy is absorbed through plastic collapse of cells in the densification core. Hence, the stress waves travelled to the distal ends are not intense enough to collapse the cells at this area [29]. This phenomenon can be further investigated by integrating with the crushing strength analysis given below.  Figure 7, which clearly illustrates the gradual transition among three specific crushing modes. As described above, the main differences among the three crushing modes lie in the location of initial collapse band (marked in Figure 7 with red box) and the growth process In Figure 8, the corresponding crushing mode for each simulation case has been identified with a solid marker, forming a mode classification map for the 3D re-entrant auxetic structure. As depicted in the figure with the dashed lines, an approximately linear relativity can be recorded between the critical transition velocity, and the relative density, showing an increase of with increasing . This can be understood by referring to the shock propagation research of cellular system, where there is a widely accepted one-dimensional shock theory using rigid, perfectly plastic, locking (R-PP-L) idealization model, revealing the dependence of critical velocity on the density of the structure [30][31][32]. By employing this theory, a similar conclusion for the regular hexagonal honeycomb has been drawn in previous studies [26,27] and the transition velocity among different modes has been further expressed intuitively [28,31].  In Figure 8, the corresponding crushing mode for each simulation case has been identified with a solid marker, forming a mode classification map for the 3D re-entrant auxetic structure. As depicted in the figure with the dashed lines, an approximately linear relativity can be recorded between the critical transition velocity, v t and the relative density, showing an increase of v t with increasing ρ R . This can be understood by referring to the shock propagation research of cellular system, where there is a widely accepted one-dimensional shock theory using rigid, perfectly plastic, locking (R-PP-L) idealization model, revealing the dependence of critical velocity on the density of the structure [30][31][32]. By employing this theory, a similar conclusion for the regular hexagonal honeycomb has been drawn in previous studies [26,27] and the transition velocity among different modes has been further expressed intuitively [28,31]. In Figure 8, the corresponding crushing mode for each simulation case has been identified with a solid marker, forming a mode classification map for the 3D re-entrant auxetic structure. As depicted in the figure with the dashed lines, an approximately linear relativity can be recorded between the critical transition velocity, and the relative density, showing an increase of with increasing . This can be understood by referring to the shock propagation research of cellular system, where there is a widely accepted one-dimensional shock theory using rigid, perfectly plastic, locking (R-PP-L) idealization model, revealing the dependence of critical velocity on the density of the structure [30][31][32]. By employing this theory, a similar conclusion for the regular hexagonal honeycomb has been drawn in previous studies [26,27] and the transition velocity among different modes has been further expressed intuitively [28,31].

Dynamic Crushing Strength
To further investigate the dynamic response of auxetic structure, the crushing stresses of auxetic structure (ρ R = 1%) produced at the proximal and distal ends under different crushing velocities are given in Figure 9a,b with respect to compression strain. All the nominal stress-strain curves show an overall piecewise characteristic during crushing, featuring with the ephemeral elasticity, long plateau and final densification periods. Another characteristic is that the stress-strain curves are observed with high oscillations. The possible cause can be attributed to the localization as the cells collapse layer by layer, resulting in a sudden drop of the stress when an adjacent layer of cells starts to collapse. Thus, the fluctuations in curves are therefore expected, especially in the dynamic crushing mode (v c = 50 m/s).

Dynamic Crushing Strength
To further investigate the dynamic response of auxetic structure, the crushing stresses of auxetic structure ( = 1%) produced at the proximal and distal ends under different crushing velocities are given in Figure 9a,b with respect to compression strain. All the nominal stress-strain curves show an overall piecewise characteristic during crushing, featuring with the ephemeral elasticity, long plateau and final densification periods. Another characteristic is that the stress-strain curves are observed with high oscillations. The possible cause can be attributed to the localization as the cells collapse layer by layer, resulting in a sudden drop of the stress when an adjacent layer of cells starts to collapse. Thus, the fluctuations in curves are therefore expected, especially in the dynamic crushing mode ( = 50 m/s).
As the crushing velocity increases, the crushing stress at the proximal end increases remarkably, while the stress at the distal end has less dependence on the crushing velocity. Additionally, the crushing stress at the proximal end is higher than that at the distal end. Since the rate-independent elastic-perfectly plastic model has been adopted and the effect of strain-rate is ignored in this study, the reason for the enhanced strength can be reasonably attributed to the change of deformation mode of struts as discussed above. The above two claims indicate that the analyzed re-entrant auxetic structure has great potential in protecting the contacted part at the distal end. Densification strain is defined as a strain at which densification begins. Another trend in Figure  9a,b is observed that the stress-strain curves at both ends show a lengthened plateau with increasing the crushing velocity, increasing the strain at which the final densification period starts, which indicates an crushing velocity dependent densification strain in auxetic cellular structures. By referring to the publication [33], the dynamic densification strain is defined as the intersection of the tangents to the plateau and densification regimes in this study and the value is then determined numerically from the obtained stress-strain curves within FE simulation. Figure 9c depicts the dynamic densification strains with respect to different crushing velocities. Obviously, the dynamic densification strain increases gradually with increasing the crushing velocity and converges to a peak value when the velocity reaches around 35 m/s, showing an increase of about 10% compared As the crushing velocity increases, the crushing stress at the proximal end increases remarkably, while the stress at the distal end has less dependence on the crushing velocity. Additionally, the crushing stress at the proximal end is higher than that at the distal end. Since the rate-independent elastic-perfectly plastic model has been adopted and the effect of strain-rate is ignored in this study, the reason for the enhanced strength can be reasonably attributed to the change of deformation mode of struts as discussed above. The above two claims indicate that the analyzed re-entrant auxetic structure has great potential in protecting the contacted part at the distal end.
Densification strain is defined as a strain at which densification begins. Another trend in Figure 9a,b is observed that the stress-strain curves at both ends show a lengthened plateau with increasing the crushing velocity, increasing the strain at which the final densification period starts, which indicates an crushing velocity dependent densification strain in auxetic cellular structures. By referring to the publication [33], the dynamic densification strain ε D is defined as the intersection of the tangents to the plateau and densification regimes in this study and the value is then determined numerically from the obtained stress-strain curves within FE simulation. Figure 9c depicts the dynamic densification strains ε D with respect to different crushing velocities. Obviously, the dynamic densification strain increases gradually with increasing the crushing velocity and converges to a peak value when the velocity reaches around 35 m/s, showing an increase of about 10% compared with that for the quasi-static case. In this study, the dynamic plateau stress is defined as the average value of compression stresses for a compressive strain of between the 0.5% strain and the dynamic densification strain ε D [30]. Since in the R-PP-L idealization model [31,32], the densification strain is adopted as a rate-independent value when formulating the explicit form of dynamic plateau stress. The trend revealed in Figure 9c indicates an amendment needed to determine the dynamic plateau stress with rate sensitivity being considered, which is further discussed in the next section.
Further effect of the relative density has been analyzed, as shown in Figure 9d. As the crushing strength with high strain rate is of high interest in crashworthiness design, the density related research at the proximal end is performed within the dynamic crushing mode, where the auxetic structures with different relative densities are crushed with the same velocity, v c = 50 m/s. From the shown stress-strain curves, the peak stress and plateau stress increase with the relative density. As the variation of relative density in this study is achieved by varying the strut thickness, the thickened strut can present a higher resistance to collapse loading, increasing the crushing strength of cellular structure. Although a long plateau with high stress level can guarantee a satisfactory energy absorption performance, the increased relative density simultaneously generates a higher peak stress, which is not appropriate in crashworthiness design. Thus, the function-oriented design for auxetic structure deserves a more detailed guideline, not merely increasing the relative density and this is not included in the present study.
As analyzed above, the theoretical prediction of dynamic plateau stress is highly related to the crushing velocity and relative density. Ideally, in the shock theory proposed by Reid, assuming the cellular structure as the R-PP-L model, the dynamic plateau stress σ D , at the proximal end is given in a general form as in Reference [30]: where σ 0 is defined as the quasi-static plateau stress, ε D is the densification strain, ρ 0 is the effective density of cellular structure and v c represents the crushing velocity. For the quasi-static plateau stress σ 0 and the densification strain ε D , Gibson and Ashby [34] formulated them as a function of the relative density ρ 0 /ρ s and the yield stress of strut material σ y as below: where ρ s is the density of strut material and C and γ are dimensionless constants. The marks shown in Figure 10 are the simulation results of quasi-static plateau stress σ 0 and densification strain ε D for the structures with different relative densities. By employing Equations (3) and (4), two fitting curves have been drawn on basis of the discrete marks, where the dimensionless constants of equations are derived as: C = 1.83 and γ = 0.85. The simulation results for dynamic plateau stress of auxetic structures with different combinations of crushing velocities and relative density are listed in Table 1, where the theoretical predictions using the Equations (2)-(4) are also listed as a contrast. As shown in the table, with increasing crushing velocity, an increasing overestimation amounting to about 13% has been observed for the theoretical results above the simulation results. The deviation can be reasonably attributed to the constant assumption of ε D in the R-PP-L model, which is actually a rate-dependent coefficient [35]. Thus, to improve the prediction accuracy of the theoretical model, a modification has been performed within the empirical formula of dynamic plateau stress. Materials 2019, 12, x FOR PEER REVIEW 10 of 15 Figure 10. The quasi-static plateau stress and the densification strain with respect to relative density. As revealed in Figure 9c, the dynamic densification strain can be related to the crushing velocity and thus a piecewise-defined function has been formulated herein by two sub-functions, a linear function fitting for all values of less than the wave trapping speed [30] and a constant function fitting for the values of greater than or equal to . The formulated piecewise-defined function is given as: (2), the dynamic plateau stress considering rate-dependent is given as: Figure 10. The quasi-static plateau stress σ 0 and the densification strain ε D with respect to relative density. As revealed in Figure 9c, the dynamic densification strain can be related to the crushing velocity and thus a piecewise-defined function has been formulated herein by two sub-functions, a linear function fitting for all values of v c less than the wave trapping speed v w [30] and a constant function fitting for the values of v c greater than or equal to v w . The formulated piecewise-defined function is given as:

Substituting Equations (3) and (5) into Equation
Substituting Equations (3) and (5) into Equation (2), the dynamic plateau stress σ D considering rate-dependent ε D is given as: Figure 11 shows the comparison between the simulation results and the theoretical predictions calculated following Equation (6) as for dynamic plateau stress σ D . As clearly shown in the figure, for the given range of relative densities and velocities, the theoretical model with Equation (6) agrees well with the FE model with a maximum deviation of 4.5%. Considering the slight error due to the approximate expression of ε D , the comparison has validated the correlations given in Equation (6). Therefore, the modified formula of dynamic plateau stress σ D provides a good guidance for predicting crushing strength of auxetic structure in crashworthiness design.
Materials 2019, 12, x FOR PEER REVIEW 11 of 15 Figure 11 shows the comparison between the simulation results and the theoretical predictions calculated following Equation (6) as for dynamic plateau stress . As clearly shown in the figure, for the given range of relative densities and velocities, the theoretical model with Equation (6) agrees well with the FE model with a maximum deviation of 4.5%. Considering the slight error due to the approximate expression of , the comparison has validated the correlations given in Equation (6). Therefore, the modified formula of dynamic plateau stress provides a good guidance for predicting crushing strength of auxetic structure in crashworthiness design. Figure 11. The Dynamic plateau stress from numerical simulation and theoretical prediction using Equation (6).

Energy Absorption Capacity
To explore the effects of crushing velocity and relative density on the energy absorption capacity of the re-entrant auxetic structure ( = 8%), the plastic energy dissipated during the crushing process has been monitored and plotted against the overall compression strain, as shown in Figure  12. The frictional dissipated energy due to friction between the contacted struts has also been plotted in the figure.
As can be seen, the frictional dissipated energy accounts for a small proportion of total internal energy and the major energy absorption mechanism is dominated by the plastic bending and axial compression of struts [36]. The plotted curves show that all the plastic energy dissipations under different crushing velocities feature an almost linearly increasing trend with the increase of compression strain. In addition, the plastic energy dissipation in the same compression strain increases remarkably as the crushing velocity increases, especially for higher velocities. This might be because of the increased crushing strength of structure under high crushing velocity.

Energy Absorption Capacity
To explore the effects of crushing velocity and relative density on the energy absorption capacity of the re-entrant auxetic structure (ρ R = 8%), the plastic energy dissipated U p during the crushing process has been monitored and plotted against the overall compression strain, as shown in Figure 12. The frictional dissipated energy U f due to friction between the contacted struts has also been plotted in the figure. for the given range of relative densities and velocities, the theoretical model with Equation (6) agrees well with the FE model with a maximum deviation of 4.5%. Considering the slight error due to the approximate expression of , the comparison has validated the correlations given in Equation (6). Therefore, the modified formula of dynamic plateau stress provides a good guidance for predicting crushing strength of auxetic structure in crashworthiness design.

Energy Absorption Capacity
To explore the effects of crushing velocity and relative density on the energy absorption capacity of the re-entrant auxetic structure ( = 8%), the plastic energy dissipated during the crushing process has been monitored and plotted against the overall compression strain, as shown in Figure  12. The frictional dissipated energy due to friction between the contacted struts has also been plotted in the figure.
As can be seen, the frictional dissipated energy accounts for a small proportion of total internal energy and the major energy absorption mechanism is dominated by the plastic bending and axial compression of struts [36]. The plotted curves show that all the plastic energy dissipations under different crushing velocities feature an almost linearly increasing trend with the increase of compression strain. In addition, the plastic energy dissipation in the same compression strain increases remarkably as the crushing velocity increases, especially for higher velocities. This might be because of the increased crushing strength of structure under high crushing velocity.  As can be seen, the frictional dissipated energy accounts for a small proportion of total internal energy and the major energy absorption mechanism is dominated by the plastic bending and axial compression of struts [36]. The plotted curves show that all the plastic energy dissipations under different crushing velocities feature an almost linearly increasing trend with the increase of compression strain. In addition, the plastic energy dissipation in the same compression strain increases remarkably as the crushing velocity increases, especially for higher velocities. This might be because of the increased crushing strength of structure under high crushing velocity.
To evaluate the energy absorption efficiency from a standardized perspective, two indices have been defined by involving the quality of structure and the quasi-static dissipation energy, giving the specific plastic energy dissipation U s p and the normalized plastic energy dissipation U p [37], respectively. They are given herein as and respectively, where m = ρ 0 V 0 = ρ 0 Ah z is the mass of cellular structure, A is the effective cross-section area, h z is the effective height along crushing direction and σ yc is the effective yield strength of cellular structure and can be calculated as σ yc = 0.5ρ 2 R σ y . According to the shock wave theory, for the continuum-based R-PP-L model, the plastic energy dissipation can be explicitly written as Applying the rate-dependent ε D function in Equations (6)- (9) gives Substituting Equation (10) into Equations (7) and (8), the theoretical predictions of U s p and U p can thus be obtained.
By adopting the above model, the theoretical predictions of the specific plastic energy dissipation U s p and the normalized plastic energy dissipation U p at densification crushing are presented in Figure 13 with respect to the normalized crushing velocity parameter V = v c /v w , respectively. To compare with the theoretical predictions, the two energy absorption indices calculated directly from the crushing simulations within the analyzed velocity and density range have also been presented in Figure 13. The comparison between the predictions and the FE results shows a satisfactory agreement, validating the effectiveness of theoretical model. To evaluate the energy absorption efficiency from a standardized perspective, two indices have been defined by involving the quality of structure and the quasi-static dissipation energy, giving the specific plastic energy dissipation and the normalized plastic energy dissipation [37], respectively. They are given herein as = (7) and = ℎ (8) respectively, where = = ℎ is the mass of cellular structure, is the effective crosssection area, ℎ is the effective height along crushing direction and is the effective yield strength of cellular structure and can be calculated as = 0.5 . According to the shock wave theory, for the continuum-based R-PP-L model, the plastic energy dissipation can be explicitly written as Applying the rate-dependent function in Equation (6) to Equation (9) gives Substituting Equation (10) into Equations (7) and (8), the theoretical predictions of and can thus be obtained.
By adopting the above model, the theoretical predictions of the specific plastic energy dissipation and the normalized plastic energy dissipation at densification crushing are presented in Figure 13 with respect to the normalized crushing velocity parameter V = / , respectively. To compare with the theoretical predictions, the two energy absorption indices calculated directly from the crushing simulations within the analyzed velocity and density range have also been presented in Figure 13. The comparison between the predictions and the FE results shows a satisfactory agreement, validating the effectiveness of theoretical model. It is also found from Figure 13a that, for the analyzed relative density range, a consistent trend can be summarized from the curves as the specific plastic energy dissipation of higher relative density is obviously larger than that of lower relative density. Given the lightweight as design consideration, this conclusion can be an important reference in the engineering application of auxetic cellular structures. Differing from the specific plastic energy dissipation, when crushing at low velocities, the normalized plastic energy dissipation increases as the relative density is increased (Figure 13b). While It is also found from Figure 13a that, for the analyzed relative density range, a consistent trend can be summarized from the curves as the specific plastic energy dissipation of higher relative density is obviously larger than that of lower relative density. Given the lightweight as design consideration, this conclusion can be an important reference in the engineering application of auxetic cellular structures. Differing from the specific plastic energy dissipation, when crushing at low velocities, the normalized plastic energy dissipation increases as the relative density is increased (Figure 13b). While at high crushing velocities, the structure relative density possesses negative effect on the capacity of energy absorption, demonstrating a shifting point at the intersection of the curves. So under high-velocity crushing, lower relative density results in higher normalized plastic energy dissipation during the crushing process and exhibits a higher rate of increase in energy dissipation capacity with an increasing crushing velocity. The shifting point for the proportional relation between the crushing velocity and the relative density is observed to be around V = 1, which actually represents the critical transition velocity between the transition mode and the dynamic mode. The results suggest that the crush hardening effect is more prominent for cellular structure with lower relative density under high-velocity crushing. To further validate the theoretical model, extensive experiments are suggested, which will be conducted in future work.

Conclusions
In the present work, a three-dimensional auxetic cellular structure consisting of re-entrant hexagonal cells has been proposed. A beam-based finite element method has been utilized to investigate the dynamic behaviors of the 3D re-entrant auxetic structure under direct crushing. Insights into the influence of crushing velocity and relative density on the crushing mode, crushing strength and energy absorption have been presented.
The detailed deformation process of the 3D re-entrant auxetic structure under the crushing velocity ranging from 2 m/s to 50 m/s has been described with the features of local densification and lateral contraction due to auxetic effect and three specific crushing modes have been observed, namely quasi-static, transition and dynamic mode. Further, the dependence of the crushing velocity and the relative density on the crushing mode has been identified, forming a mode classification map for the 3D re-entrant auxetic structure. The crushing strength of the 3D re-entrant auxetic structure reveals to increase with the increase in crushing velocity and relative density. Moreover, an analytical formula of dynamic plateau stress has been deduced by introducing the rate-dependent densification strain, which proves to present theoretical predictions agreeing well with simulation results. By establishing an analytical model, we have also studied the role of relative density on the energy absorption capacity of the 3D re-entrant auxetic structure with respect to crushing velocity. The theoretical and simulation results both show that the specific plastic energy dissipation of the structure increases consistently with the increase in the relative density, while the normalized plastic energy dissipation has an opposite sensitivity to the relative density when the crushing velocity exceeds the critical transition velocity.
The results in this work can be a preliminary complement to the dynamic analysis of the 3D re-entrant auxetic structure and provides an extensive reference for the crushing resistance design of the auxetic structure. Future work will focus on the experimental validation of the results and considering more factors in mechanical analysis such as the size effect and the matrix material properties.