Using Finite Element Approach for Crashworthiness Assessment of a Polymeric Auxetic Structure Subjected to the Axial Loading

Polyurethane foams are one of the most common auxetic structures regarding energy absorption enhancement. This present study evaluates the result reliability of two different numerical approaches, the H-method and the P-method, to obtain the best convergence solution. A polymeric re-entrant cell is created with a beam element and the results of the two different methods are compared. Additionally, the numerical results compare well with the analytical solution. The results show that there is a good agreement between converged FE models and the analytical solution. Regarding the computational cost, the P-method is more efficient for simulating the re-entrant structure subjected to axial loading. During the second part of this study, the re-entrant cell is used for generating a polymeric auxetic cellular tube. The mesh convergence study is performed on the cellular structures using the H- and P- methods. The cellular tube is subjected to tensional and compressive loading, the module of elasticity and Poisson’s ration to calculate different aspect ratios. A nonlinear analysis is performed to compare the dynamic response of a cellular tube versus a solid tube. The crashworthiness indicators are addressed and the results are compared with equivalent solid tubes. The results show that the auxetic cellular tubes have better responses against compressive loading. The primary outcome of this research is to assess a reliable FE approach for re-entrant structures under axial loading.


Introduction
Modern technology requires new physical material properties. Study of the material using negative Poisson's ratio (NPR) is taken into account as an option for material property enhancement [1]. Utilizing re-entrant structures with NPR provides the ability to improve mechanical properties in different applications such as medical, automotive, textile engineering etc. [2,3]. Generally, changing the chemical properties and cell shape are two methods for improving mechanical behaviour [4,5]. Most of the researches in this area focus on changing the chemical elements [6][7][8]. Recently, an attempt is being made to enhance the mechanical properties by changing the cell shape of the auxetic structures. Concerning the elastic region, the mechanical properties of materials are mainly influenced by four elastic constants that are elastic moduli or Young's moduli (E), shear moduli (G), bulk moduli (K) and Poisson's ratio (ν). Poisson's ratio is defined as the ratio of transverse contraction strain to longitudinal extension strain concerning the direction of stretching force applied. The formula of Poisson's ratio contains a minus sign. Therefore, general materials possess a positive Poisson's ratio [9]. Based on elasticity theorem, the mentioned constants for isotropic materials are dependent on the following set of equations (Wei Yang et al., 2004) [4]: To enhance the G and K of the structure, considering that E is not a variable, Poisson's ratio is the option that can be improved to obtain the desired G and K [9,10]. Factually, most of the materials have a positive Poisson's ratio and the materials with a negative Poisson's ratio are made by foam and porous materials such as polyurethane or aluminium foams [11][12][13].
Some research has been conducted to convert the material with positive Poisson's ratios to the NPR [14]. However, the NPR structure can be made by isotropic material in the macro scale, thus, in inter-atomic bonds the re-entrant strucre with NPR can be observed. Several studies have been conducted to examine the global stiffness of auxetic structures. Generally, cellular auxetic structures are designed by repeating a unit cell, and their effective stiffness can be determined through a unit cell [15,16]. To obtain a negative Poisson's ratio, many different designs such as re-entrant cells as well as rotating rectangles and triangles, arrow-heads, and star-shaped configurations are proposed [17][18][19]. One of the most interesting auxetic structures is the re-entrant cell with scale-independent properties [20,21]. The deformation behaviour of re-entrant structures as a cellular material can be implemented at every scale range from nano scale [22] to a macro level.
The effect of cell-wall alignment on the dynamic responses of a re-entrant honeycomb structure was studied by Zhang et al., [23]. It was recognized that increasing the impact speed, cell angle and relative density leads to increasing the crashworthiness capability of the auxetic structure. A similar investigation on the hexagonal structure has been carried out by Hu and his colleagues [24], where it was realized that a honeycomb with a cell angle of 45 • had better energy absorption performance when subjected to impact loading. Three different deformation patterns under impact loading of the cellular structures were discussed by Zou et al., [25]. Based on the collapsing mechanism of hexagonal honeycomb structures, an analytical formulation for the energy absorption ability was derived by Hu et al., [26]. They later developed an analytical model validated by simulations to anticipate the crashworthiness of hexagonal honeycombs under low impact loading [27]. To exhibit their auxetic property, auxetic materials should possess substantial porosity in their microstructures [1]. Therefore, geometric complexity and porosity in these structures under numerical study make their analysis cumbersome. Regarding terms of finite element work, previous studies in this area show that the numerical simulation of energy absorption by auxetic materials, especially in three-dimensional re-entrant structures, is still limited and sparse, thus needing further development [28].
Impact resistance and energy absorption of auxetic structures is also an interesting topic to which researchers have been paying much attention [29][30][31]. Reid and Peng [32] developed the one-dimensional shock theory estimating the characteristics of a crushing front through the wood subjected to uniaxial impact. Ruan et al., [33] performed finite element analysis (FEA) to study the effectiveness of impact speed and the wall thickness of a cell on the localized deformation state and plateau stress. The FE result highly depends on the element size and degree of the polynomial solution. The solution convergence in the FE problem is the most crucial matter for obtaining a single correct solution [34]. Different methods are presented to demonstrate the numerical convergence. The different ways for increasing the model's degree of freedom can be categorized into two main branches: element refinement (H-method) and increased polynomial degree with the highest accuracy (P-method) [35]. The H-method result will be more accurate by increasing the number of elements. Using a finer mesh, in other words, brings more accuracy to the model. However, choosing the number of fine mesh needs expertise and increases the computational cost as well. The P-method utilizes the complex shape function, keeping constant the number of elements. The initial iteration uses the first-order polynomial shape function and, in the following run, the order of the shape function can be increased [34]. While the difference between the two last iterations remains in a specific tolerance, the simulation is running. Recently, a combination of the H-and P-methods has been used to take advantage of the desired element size as well as the exponential rate of convergence.
Several studies have been carried out on cellular structures in experimental, numerical, and analytical approaches [36][37][38]. However, when the structures are small in scale, especially in micro or nano-scale, the mathematical models are not precise enough to be validated with an experiment. Therefore, a reliable and fast response model is crucial to simulate the dynamic response of the cellular structure. This present study has several aims: (1) To formulate the total strain energy of a 3D re-entrant cell analytically; (2) To evaluate the effectiveness of the different numerical convergence methods for the re-entrant structure subjected to axial loading using the FE model with the beam element generated; (3) Examining the strain energy of the H-version compared to the P-version regarding CPU time and number of elements; (4) Development of a unit cell that generates a cellular tube with NPR behaviour, as well as, (5) The crashworthiness indicators of the proposed structures compared to the conventional tube with different aspect ratios. The FE formulation can be proposed as an efficient and reliable model when a complex cellular structure is subjected to axial loading.

Analytical Solution for 3D Re-Entrant Cellular Structure
The analytical model used in the present study was developed by Shokri Rad et al., [9], which is based on Castigliano's theorem. According to this model, for a linear elastic solid, the energy theorem mentioned below presents the relationship among the potential strain energy saved in a beam element, the load, and the bending moment distribution across the same component: where: δ is the corresponding deflection by the bending moment magnitude, M. P, is the applied axial force and ∅ is the beam slope provided by the bending moment. V is the unit volume of structure and, also, σ ii and ε ii are the elements of stress and strain tensors, respectively. The strain energy for a one-dimensional (1D) beam element can be calculated along the length of the beam. Therefore, if a beam is subjected to a bending moment composed by an axial load, then the energy can be calculated as: Referring to Figure 1a, m is located at the centre of the rod connected to the neighbouring cell and the structure is symmetric. Also, it is assumed that the connecting rod, mn, is rigid. Hence, the point n is fixed as well. As a result, considering the symmetric property of the cellular structure and point n being fixed, Figure 1b can be considered to render the analysis easier. Regarding this figure, P q and m q are the axial force and bending moment at the centre of the beam bc. (point q in Figure 1a), respectively. Modula of elasticity, second moment of inertia and cross section area are denoted by E, I, respectively. Table 1 shows the nomenclatures used in this paper.  Due to the symmetric characteristic of the model, it can be assumed that the beam has no horizontal deflection and slope at the center (Point ). Therefore, the total strain energy in Figure 1b can be written as: To simplify the above equation, we define: = −( sin 2 )( + )  Due to the symmetric characteristic of the model, it can be assumed that the beam bc has no horizontal deflection and slope at the center (Point q). Therefore, the total strain energy in Figure 1b can be written as: To simplify the above equation, we define: To define the virtual force and moment as a function of applied force it can be written as follows: Let us define: To obtain: Therefore, by substituting Equation (14) into Equation (5), the total strain energy for a re-entrant cell can be calculated.

Finite Element Approaches and Simulation
Here, the nonlinear finite element code ABAQUS/Implicit is employed to conduct the computer simulations. The geometry should be presented in micro or nanoscale, but to illustrate the effect of the FE result, the geometry is assumed in macro scale. Both unit cell and re-entrant structures are generated from the beam element (B3). To compare the effect of mesh refinement through the H-and P-versions, the elastic properties of aluminium alloy taken from a previous study [34] are considered for the unit cell and re-entrant structures. These properties are shown in Table 2.  26 36 To calculate the structure's module of elasticity and Poisson's ratio, the simulation is performed in the elastic region. Both the H-and P-versions' refinements are employed for simulations. To reduce the time of calculation, and to simplify the FE formulation, the material properties of aluminium are used rather than polymeric material. Actually, the material properties should be taken from a polymeric material, however, to study the overall structural behaviour, an elastic-plastic material like aluminium shows a clear response when a structure is subjected to a crushing load. Here, the region of interest is only the elastic region hence, the overall behaviour of the structure, regardless of material type, can be used as the guideline for using another kind of material such as auxetic polyurethane foam [14]. It can be said that the outcome results of this simulation are far from a real test, however, to perform a qualitative study this procedure is acceptable [1,9,22].

Unit Cell Modelling
Two models with different mesh refinement approaches were used. Regarding the H-version approach, the size of the mesh was refined gradually from 10 to 2 mm to achieve convergence and result stability. The one-dimensional beam element with a linear interpolation function was used in this case. Conversely, the P-version approach used only the one-dimensional element, so the size of the element remained 10 mm, however, the interpolation function was different by using first-, secondand third-order polynomial interpolation functions. The dimensions to generate the FE model were taken from Table 1. Figure 2 illustrates the applied loads and boundary conditions assigned to the rendered beam model. A couple of 10 N axial forces were applied at the centre of the structure along the Z direction, and all four neighbouring connecting rods were fixed in all degrees of freedom.
Polymers 2020, 12, x FOR PEER REVIEW 6 of 16 polyurethane foam [14]. It can be said that the outcome results of this simulation are far from a real test, however, to perform a qualitative study this procedure is acceptable [1,9,22].

Unit Cell Modelling
Two models with different mesh refinement approaches were used. Regarding the H-version approach, the size of the mesh was refined gradually from 10 to 2 mm to achieve convergence and result stability. The one-dimensional beam element with a linear interpolation function was used in this case. Conversely, the P-version approach used only the one-dimensional element, so the size of the element remained 10 mm, however, the interpolation function was different by using first-, second-and third-order polynomial interpolation functions. The dimensions to generate the FE model were taken from Table 1. Figure 2 illustrates the applied loads and boundary conditions assigned to the rendered beam model. A couple of 10 N axial forces were applied at the centre of the structure along the Z direction, and all four neighbouring connecting rods were fixed in all degrees of freedom.

Cellular Tube Modelling
To generate a polymeric auxetic cellular structure, a unit cell was repeated in a circular direction. Repeating the single layer along the axial direction generated the cellular tube. Figure 3 shows a unit cell, a single layer, and the final geometry of a cellular tube. Dividing different lengths of tubes over the average diameter of the tube (382 mm) made different aspect ratios (L/D). Different aspect ratios (L/D = 1-5) were considered to investigate the mechanical properties of the polymeric auxetic structure subjected to the tension and compressive loading. A reference point was connected to the lower set of nodes at the bottom centre of the cellular tube. The upper set of nodes of the tube were connected to the upper reference point to apply the axial loading. To verify the mesh convergent study, the element size from the unit cell was used for both H-version and P-version methods. Regarding this, the results for 2, 3 and 5 mm element sizes for the H-version were compared with a second-order polynomial function. As a result, the converged values of stresses, strain energy and CPU time were compared.

Cellular Tube Modelling
To generate a polymeric auxetic cellular structure, a unit cell was repeated in a circular direction. Repeating the single layer along the axial direction generated the cellular tube. Figure 3 shows a unit cell, a single layer, and the final geometry of a cellular tube. Dividing different lengths of tubes over the average diameter of the tube (382 mm) made different aspect ratios (L/D). Different aspect ratios (L/D = 1-5) were considered to investigate the mechanical properties of the polymeric auxetic structure subjected to the tension and compressive loading. A reference point was connected to the lower set of nodes at the bottom centre of the cellular tube. The upper set of nodes of the tube were connected to the upper reference point to apply the axial loading. To verify the mesh convergent study, the element size from the unit cell was used for both H-version and P-version methods. Regarding this, the results for 2, 3 and 5 mm element sizes for the H-version were compared with a second-order polynomial function. As a result, the converged values of stresses, strain energy and CPU time were compared.
When a compressive loading is applied to the structure it is possible to have a contact between beam links. To avoid penetration of different beam links together, a general contact algorithm was considered. The general contact method includes surface-to-surface self-contacting between different beams. The tangential behaviour followed the penalty method with a 0.09 friction coefficient, and the normal behaviour was considered as the hard contact algorithm [39]. The main objective of a mesh convergence study for FE simulation is to reach a stable value of stress or magnitude of energy. Therefore, showing strain energy and, especially, von Mises stress are two main criteria that prove the validity of the numerical simulation, whereas, using principal stresses in this study due to the non-linear behaviour of the material and geometry was not appropriate. Hence, the von Mises stress is a key value to justify the FE method results [40][41][42][43][44][45]. When a compressive loading is applied to the structure it is possible to have a contact between beam links. To avoid penetration of different beam links together, a general contact algorithm was considered. The general contact method includes surface-to-surface self-contacting between different beams. The tangential behaviour followed the penalty method with a 0.09 friction coefficient, and the normal behaviour was considered as the hard contact algorithm [39]. The main objective of a mesh convergence study for FE simulation is to reach a stable value of stress or magnitude of energy. Therefore, showing strain energy and, especially, von Mises stress are two main criteria that prove the validity of the numerical simulation, whereas, using principal stresses in this study due to the non-linear behaviour of the material and geometry was not appropriate. Hence, the von Mises stress is a key value to justify the FE method results [40][41][42][43][44][45].
During the next step, to study the energy absorption of the cellular structures with different aspect ratios under compressive loading, a non-linear plastic behaviour was used [34]. Considering each tube, the applied load was increased as ramp until the maximum deflection reached 70% of the initial length of the cellular tube. Occurring there, the energy absorption of the structure from the area under the load-displacement curve was estimated. To compare the impact resistance capability, an equivalent solid tube was designed. Both structures (cellular tubes and solid tubes) had the same average diameter and the same weight. Figure 4 shows the deformed and unreformed shape of cellular tubes and conventional tubes with different aspect ratios. During the next step, to study the energy absorption of the cellular structures with different aspect ratios under compressive loading, a non-linear plastic behaviour was used [34]. Considering each tube, the applied load was increased as ramp until the maximum deflection reached 70% of the initial length of the cellular tube. Occurring there, the energy absorption of the structure from the area under the load-displacement curve was estimated. To compare the impact resistance capability, an equivalent solid tube was designed. Both structures (cellular tubes and solid tubes) had the same average diameter and the same weight. Figure 4 shows the deformed and unreformed shape of cellular tubes and conventional tubes with different aspect ratios.

H-Method Results for Convergence Study of a Unit Cell
The key point in the FE solution problem is the rate of convergence and, subsequently, the

H-Method Results for Convergence Study of a Unit Cell
The key point in the FE solution problem is the rate of convergence and, subsequently, the accuracy of the results. Figure 5 demonstrates the mesh convergence analysis using the H-method FE results. It can be seen that the initial course mesh was not appropriate for prediction of the structure behaviour and, by refining the element size, the results became stable and finally remained at a certain amount of strain energy. The converged strain energy and stress at the end of the H-method refinement were 1.7 kJ and 25.02 MPa, respectively.

H-Method Results for Convergence Study of a Unit Cell
The key point in the FE solution problem is the rate of convergence and, subsequently, the accuracy of the results. Figure 5 demonstrates the mesh convergence analysis using the H-method FE results. It can be seen that the initial course mesh was not appropriate for prediction of the structure behaviour and, by refining the element size, the results became stable and finally remained at a certain amount of strain energy. The converged strain energy and stress at the end of the H-method refinement were 1.7 kJ and 25.02 MPa, respectively. Additionally, Figure 6 a,b shows the von Mises stress and overall deformation of the re-entrant cell, respectively, for the H-method analysis. The degrees of freedom for the structure was varying Additionally, Figure 6a,b shows the von Mises stress and overall deformation of the re-entrant cell, respectively, for the H-method analysis. The degrees of freedom for the structure was varying from 300 to 1452, and the CPU time was 0.8 to 1.7 s as well. Decreasing the size of the element from 10 mm to 2 mm made it complex to attain the converged solution. To evaluate a structure with various numbers of cells, the computational cost would be increased significantly.
Polymers 2020, 12, x FOR PEER REVIEW 9 of 16 from 300 to 1452, and the CPU time was 0.8 to 1.7 s as well. Decreasing the size of the element from 10 mm to 2 mm made it complex to attain the converged solution. To evaluate a structure with various numbers of cells, the computational cost would be increased significantly.  Figure 7 depicts the convergence study for the P-method analysis with three different polynomial orders. The second-order polynomial error showed that increasing the higher order polynomial function was not necessary and the second-order polynomial function could be utilized for simulation reliability.  Figure 7 depicts the convergence study for the P-method analysis with three different polynomial orders. The second-order polynomial error showed that increasing the higher order polynomial function was not necessary and the second-order polynomial function could be utilized for simulation reliability.  Figure 7 depicts the convergence study for the P-method analysis with three different polynomial orders. The second-order polynomial error showed that increasing the higher order polynomial function was not necessary and the second-order polynomial function could be utilized for simulation reliability. The degree of freedom and CPU time for the second-order polynomial function were 612 and 0.8 s, respectively. Comparison between the P-method and H-method, regarding corresponding stress after convergence, represents a less than 1% error. The stress and deformation contour indicate that the re-entrant cells are capable of absorbing energy. Table 3 compares the result of different simulations between the H-method and P-method. The degree of freedom and CPU time for the second-order polynomial function were 612 and 0.8 s, respectively. Comparison between the P-method and H-method, regarding corresponding stress after convergence, represents a less than 1% error. The stress and deformation contour indicate that the re-entrant cells are capable of absorbing energy. Table 3 compares the result of different simulations between the H-method and P-method.

Unit Cell Strain Energy
The total strain energy from equation 5 was calculated and the result can be seen in Figure 8. The magnitude of the strain energy using the analytical method for the re-entrant cell was 1.703 kJ and, in comparison with the converged numerical solution, the error was less than 1%. However, the exact solution can be useful to highlight the effect of convergence in an FE simulation. Table 4 compares the results of different FE methods for a cellular tube with L/D = 1. The second-order polynomial P-version results presented greater efficiency in terms of strain energy and CPU time. It can be said that using a second-order P-version FE formulation increased the speed of calculation significantly. Comparing the results of the second-order FE formulation with converged H-version results (2 mm element size) shows that this method can accelerate the calculation time 287 s faster.

Unit Cell Strain Energy
The total strain energy from equation 5 was calculated and the result can be seen in Figure 8. The magnitude of the strain energy using the analytical method for the re-entrant cell was 1.703 kJ and, in comparison with the converged numerical solution, the error was less than 1%. However, the exact solution can be useful to highlight the effect of convergence in an FE simulation.  Table 4 compares the results of different FE methods for a cellular tube with L/D = 1. The secondorder polynomial P-version results presented greater efficiency in terms of strain energy and CPU time. It can be said that using a second-order P-version FE formulation increased the speed of calculation significantly. Comparing the results of the second-order FE formulation with converged H-version results (2 mm element size) shows that this method can accelerate the calculation time 287 s faster.

Structure Stiffness and Poisson's Ratio
To calculate the structure elasticity and Poisson's ration, the cellular tubes with different aspect ratios were subjected to axial loading. The load and displacement data for every 10% loading increment were recorded, and corresponding elastic modules were estimated. The simulation was stopped when the elongation reached 50% of the initial length. Moreover, the Poisson's ratio of each cellular tube was calculated over time. Figure 9 shows the load-displacement curve of the polymeric auxetic tube under tensile loading with different aspect ratios.
Having transferred the load-displacement to the stress-strain value, the module of elasticity and Poisson's ratio can be calculated from equation (1). Figure 10 shows the average module of elasticity and corresponding Poisson's ratio of structure. Considering this figure, it can be interpreted that the Poisson's ratio can be varied in the elastic region using negative values. After applying 50% of the yield stress on the cellular tubes, the structures still had a negative Poisson's ratio. The average module of elasticity for all cellular tubes was 55.39 GPa and the Poisson's ratio of the structure, after 50% elongation, reached −0.5. The error bars in this figure present the upper and lower bounds of different stress magnitudes related to the different aspect ratios of the cellular tubes. The red curve is the mean value of the stress-strain curve of the cellular tube. The black line also represents the slope of the stress-strain curve as the average module of elasticity of all cellular tubes.
To calculate the structure elasticity and Poisson's ration, the cellular tubes with different aspect ratios were subjected to axial loading. The load and displacement data for every 10% loading increment were recorded, and corresponding elastic modules were estimated. The simulation was stopped when the elongation reached 50% of the initial length. Moreover, the Poisson's ratio of each cellular tube was calculated over time. Figure 9 shows the load-displacement curve of the polymeric auxetic tube under tensile loading with different aspect ratios. Having transferred the load-displacement to the stress-strain value, the module of elasticity and Poisson's ratio can be calculated from equation (1). Figure 10 shows the average module of elasticity and corresponding Poisson's ratio of structure. Considering this figure, it can be interpreted that the Poisson's ratio can be varied in the elastic region using negative values. After applying 50% of the yield stress on the cellular tubes, the structures still had a negative Poisson's ratio. The average module of elasticity for all cellular tubes was 55.39 GPa and the Poisson's ratio of the structure, after 50% elongation, reached −0.5. The error bars in this figure present the upper and lower bounds of different stress magnitudes related to the different aspect ratios of the cellular tubes. The red curve is the mean value of the stress-strain curve of the cellular tube. The black line also represents the slope of the stress-strain curve as the average module of elasticity of all cellular tubes. Auxetic structures are intrinsically weak under bending and torsion. These kinds of structures are designed for applications which are subjected to axial loading, in other words. This is the reason, proven in several studies [1,2,14], that using auxetic polyurethane foam as a filler foam for structures under compression increases the impact resistance of structures significantly. The value of the shear modulus for the present auxetic structure can be calculated using Equation (1), however, due to the negative Poisson's ratio of the structure, the shear modulus has a negative or near to zero value. Therefore, in this study the shear modulus calculation is neglected. Figure 11 demonstrates the dynamic responses of cellular tubes and conventional tubes with different aspect ratios. It can be seen that the load-displacement curve from the cellular tube presented a better energy absorption capability. Concerning the conventional tubes, the first peak load was higher than the cellular tube, however, the applied load had a sharp drop over time. Conversely, the loading response from the cellular tube showed that the compressive behaviour of the structure had a good progressive collapse as opposed to the conventional tubes. This phenomenon can be considered as the significant characteristic of the auxetic cellular structure subjected to compressive loading. Table 5 presents the crashworthiness indicator for conventional and cellular tubes with different aspect ratios. Calculation of the area under load-displacement cure yielded the energy absorption of Auxetic structures are intrinsically weak under bending and torsion. These kinds of structures are designed for applications which are subjected to axial loading, in other words. This is the reason, proven in several studies [1,2,14], that using auxetic polyurethane foam as a filler foam for structures under compression increases the impact resistance of structures significantly. The value of the shear modulus for the present auxetic structure can be calculated using Equation (1), however, due to the negative Poisson's ratio of the structure, the shear modulus has a negative or near to zero value. Therefore, in this study the shear modulus calculation is neglected. Figure 11 demonstrates the dynamic responses of cellular tubes and conventional tubes with different aspect ratios. It can be seen that the load-displacement curve from the cellular tube presented a better energy absorption capability. Concerning the conventional tubes, the first peak load was higher than the cellular tube, however, the applied load had a sharp drop over time. Conversely, the loading response from the cellular tube showed that the compressive behaviour of the structure had a good progressive collapse as opposed to the conventional tubes. This phenomenon can be considered as the significant characteristic of the auxetic cellular structure subjected to compressive loading. Comparing the results of solid tubes and cellular tubes show that using a cellular structure can increase the crash resistance more than 30%. The interesting point regarding the cellular structure is the progressive collapse and increase of the mean crash load in comparison with conventional tubes. To design lightweight components, polymeric auxetic cellular structures could be a good option rather than foam-filled tubes when subjected to axial loading. Comparison of the results of cellular tubes with different L/D shows that increasing the value of the L/D will affect the crashworthiness indicators. Considering all crashworthiness indicators, an increasing trend alongside an increasing L/D can be observed. However, from Figure 10, it can be said that the stiffness of structure for all L/D were almost the same. This value can be seen by comparing the value of Pmax from Table 5. Figure 11. Load-displacement curve for solid and cellular tubes with different aspect ratios. Figure 11. Load-displacement curve for solid and cellular tubes with different aspect ratios. Table 5 presents the crashworthiness indicator for conventional and cellular tubes with different aspect ratios. Calculation of the area under load-displacement cure yielded the energy absorption of different samples. The other crashworthiness indicators such as maximum peak load (P max ), mean crush load (P ave ), crash force efficiency (CFE) and specific energy absorption (SEA) can be estimated from ref. [5].

Energy Absorption Evaluation
Comparing the results of solid tubes and cellular tubes show that using a cellular structure can increase the crash resistance more than 30%. The interesting point regarding the cellular structure is the progressive collapse and increase of the mean crash load in comparison with conventional tubes. To design lightweight components, polymeric auxetic cellular structures could be a good option rather than foam-filled tubes when subjected to axial loading. Comparison of the results of cellular tubes with different L/D shows that increasing the value of the L/D will affect the crashworthiness indicators. Considering all crashworthiness indicators, an increasing trend alongside an increasing L/D can be observed. However, from Figure 10, it can be said that the stiffness of structure for all L/D were almost the same. This value can be seen by comparing the value of P max from Table 5.

Conclusions
The different approaches to FE simulation were proposed to evaluate the effectiveness of convergence solutions on the corresponding behaviour of a re-entrant cell, such as in polymeric auxetic materials. Regarding DOF, CPU time, and number of elements, the results show that the second-order P-method was the most efficient solution for a re-entrant cell subjected to axial loading. Moreover, there was a good agreement between the converged numerical results (H-and P-methods) and the analytical solution. This research presents that the P-method is the most reliable approach to simulate the re-entrant cell. During the second section of this study, the re-entrant cell was used for generating an auxetic cellular tube. The convergence study was repeated for a cellular tube and the results showed that the second-order P-version formulation was the more efficient method to model re-entrant structures. Polyurethane foams are one of the most common auxetic structures in the field of energy absorption enhancement. The cellular tube was subjected to tensional and compressive loading, and the module of elasticity and Poisson's ration for different aspect ratios were calculated. A nonlinear analysis was performed to compare the dynamic response of a cellular tube and solid tube. The crashworthiness indicators were addressed and the results were compared with equivalent solid tubes. The results show that the auxetic cellular tubes had better responses against compressive loading.

Conflicts of Interest:
The authors declare no conflicts of interest.