© A. Petit et al., Published by EDP Sciences, 2023

1 Introduction

Structured materials, known for their multifunctional application [1], are remarkable for their particular stress-strain curves: after an elastic part where the curve is quasi-linear, a plateau stress interval is reached where the stress is almost constant and finally the curves rise abruptly at densification strain [2]. During the past few decades [3,4], a new class of structured materials have emerged called ≪ auxetic ≫ (from Greek ≪ auxetos ≫: that may be increased). This kind of structure is characterized by a negative Poisson’s ratio which gives it, compared to conventional, cellular structure a higher shear modulus, better indentation resistance and energy absorption ([5–7]).

The present work focuses on one of these improved properties: the energy absorption capacity of auxetic structures (AS). Among the large range of structures having a negative Poisson’s ratio, the re-entrant honeycomb is one of the most investigated geometry. Experimental samples are relatively easy to develop using 3D printing technologies [8]. In addition, this geometry can be implemented in FE models without significant difficulties. That is why it has been chosen for this study.

In-plane auxetic crushing test is the benchmark analysis to determine quasi-static and dynamic crushing strength of AS. This strategy is adopted by Ruan et al. [9], Hu et al. [10,11], and Hou et al. [12] to develop their models. In 1993, Reid et al. [13], on the basis that wood behave like cellular structure, proposed a shock model for wood using the rigid, perfectly-plastic, locking (R-P-P-L) simplification. This model of the dynamic crushing of cellular structure is validated in 2003 by Ruan et al. [9] thanks to numerical in-plane dynamic crushing of honeycombs.

Based on energy conservation principle and shockwave analogy, a theoretical (analytical) formulation of stress will be established. Finite Element simulations have confirmed the analytical model.

The objective is to propose a general methodology to assess the strength of an AS loaded by the impact of a plate.

2 Characteristics of the auxetic structure under consideration

2.1 Geometry of the structure

The complete auxetic structure, which will be investigated in compression along the Y-axis in this article, is an auxetic “re-entrant honeycomb” represented in Figure 1.

The structure is composed of cells, as the one shown in Figure 2 with the notations used for the geometrical data. Note that S₀ is the total projected area on the X–Z plane, as shown in Figure 1.

A unit-cell geometry is characterised by its neutral axis (black dotted line in Fig. 2). The horizontal edge is denoted H and the bevel edge is denoted L. The angle formed by edges H and L is denoted θ. The cell-wall thickness is denoted t_thick and b is the out-of-plane thickness.

The edges dimensions h and ℓ are defined as follow :

$$\ell =L-\frac{t_{\text{thick}}}{2\text{\hspace{0.17em}sin\hspace{0.17em}}\theta }$$(1)

$$h=H-\frac{\left(1+\text{cos\hspace{0.17em}}\theta \right)}{\text{sin\hspace{0.17em}}\theta }{t}_{\text{thick}}$$(2)

The initial cell height is defined by:

$$d_{\text{0cell}}=2Lsin\theta$$(3)

At the initial time (t = 0), the total height in the Y-direction for n cell rows (layer) is:

$$D_0=n{d}_{0\text{cell}}+t_{\text{thick}}$$(4)

Two characteristic ratios are defined:

$$\tau =\frac{t_{\text{thick}}}{L}$$(5)

$$\eta =\frac{H}{L}.$$(6)

Fig. 1 View of the complete auxetic structure (“AS”).

Fig. 2 Re-entrant auxetic cell configuration.

2.2 Relative density and densification strain

The relative density of a cellular structure is defined as the ratio of the cellular structure mass density ρ_c to the mass density ρ_s of the material of which it is constituted. For the re-entrant honeycomb, it is defined as follow :

$$\frac{{\rho }_c}{{\rho }_s}=\frac{V_s}{V_c}$$(7)

where V_c and V_s are the cellular structure and material volumes, respectively. The relative density can be written in terms of the non-dimensional ratio τ and η :

$$\frac{{\rho }_c}{{\rho }_s}=1-\frac{\text{sin\hspace{0.17em}}\theta -\frac{\tau }{2}}{\text{sin\hspace{0.17em}}\theta \left(\eta -\text{cos\hspace{0.17em}}\theta \right)}\left[\eta -\text{cos\hspace{0.17em}}\theta -\frac{\tau }{\text{sin\hspace{0.17em}}\theta }\left(\frac{1+\text{cos\hspace{0.17em}}\theta }{2}\right)\right].$$(8)

The strain will be defined by the ratio of the height variation D₀ − D to the initial height D₀, D being the height at time t:

$$\epsilon =\frac{D_0-D}{D_0}.$$(9)

At densification the strain is ε_DS and the height is D_s, while the travelled distance of the upper wall is ε_DSD₀.

At this characteristic strain, the cells are all collapsed so that all cell-walls are in contact and the compressive stress is almost the compressive strength of the material itself. So, the height is D_s = (4 n + 1) t_thick and the densification strain ε_DS is therefore expressed as:

$${\epsilon }_{\text{DS}}=1-\frac{\left(4n+1\right)}{D_0}{t}_{\text{thick}}.$$(10)

3 Analytic prediction of stress

3.1 Stress enhancement through the shock front

Now, an impact on the AS is considered. It corresponds to impact the top surface (y = D₀, Fig. 1) with a plate of mass M, while the other face (y = 0) is fixed. The velocity of the plate is υ₀ when it impacts the AS at time t₀ = 0, and υ for t > 0.

When the AS is dynamically crushed, the compaction starts from the impacted face and propagates like a shock-wave through the structure, involving a strong discontinuity (or “front”) between the undeformed and deformed regions. The main steps of the shock-wave analogy applied to cellular material have been presented by [13]. As a similar approach is adopted in this study, these steps are briefly recalled in Figure 3.

The AS is idealized with the R-P-P-L model. Therefore, the stress is raised instantaneously to the static plateau stress σ_stat at the initial time. The σ_stat expression can be obtained by considering the plastic collapse of the cells network [14] :

$${\sigma }_{\text{stat}}={\sigma }_y\frac{{\tau }^2}{2\left(\eta -\text{cos\hspace{0.17em}}\theta \right)\text{\hspace{0.17em}cos\hspace{0.17em}}\theta }$$(11)

with σ_y the yield stress of the material.

Figures 3c and d show the deformed states at times t and t + dt. Behind the front, the stress is raised to the dynamic value σ and the structure is compressed up to densification strain ε_DS. The structure is supposed to keep the same area S₀ and to deform only in the uniaxial direction.

In Figure 3, the shock wave propagates from proximal to distal end with the wave speed υ_f. From kinematic analysis, this velocity in lagrangian coordinates is expressed as:

$${\upsilon }_f=-\frac{\text{d}D}{\text{d}t}=\frac{\text{d}u}{\text{d}t}+\frac{\text{dλ}}{\text{d}t}=\frac{\upsilon }{{\epsilon }_{\text{DS}}}$$(12)

where u is the displacement of the plate and λ is the deformed length of the crushed AS. As the current velocity is υ = du/dt, the front velocity relative to the impactor is :

$$\frac{\text{d}\lambda }{\text{d}t}=\frac{1-{\epsilon }_{\text{DS}}}{{\epsilon }_{\text{DS}}}\upsilon$$(13)

The conservation of mass gives :

$$\frac{{\rho }_c}{{\rho }_{\text{DS}}}=1-{\epsilon }_{\text{DS}}=\frac{\text{λ}}{\text{λ}+u}$$(14)

where ρ_DS is the AS density when the structure is crushed up to densification.

During an infinetely small duration dt, the element of length dD (hashed in Fig. 3c) passes from an undeformed to a crushed state (hashed in Fig. 3d) and its lengths is dλ. The mass included within this volume hashed in Figure 3d, S₀ dλ, is expressed in terms of the cellular density ρ_c :

$$\text{d}m=\frac{{\rho }_c{S}_0}{1-{\epsilon }_{\text{DS}}}\text{d}\lambda .$$(15)

During this time interval, the volume element S₀ dλ considered above is submitted to the stress σ_stat on its right face and to the stress σ on its left face. The velocity of the element increases from 0 to υ during this duration. Thus the conservation of momentum of this element reads:

$$\upsilon \text{\hspace{0.17em}d}m=\left(\sigma -{\sigma }_{\text{stat}}\right)S_0\text{\hspace{0.17em}d}t.$$(16)

From equation (16) it can be deduced that:

$$\sigma ={\sigma }_{\text{stat}}+\frac{{\rho }_c}{{\epsilon }_{\text{DS}}}{\upsilon }^2.$$(17)

This formulation is the same as the one proposed by Reid [13] to assess the stress inside the crushed part (shaded part in Fig. 3). After the impact, the crushing velocity υ and the reaction forces of the AS vary constantly with time due to the dynamic coupling between the plate and the AS at the interface. It is proposed to evaluate this evolution with respect to time by using the energy conservation approach.

Fig. 3 Schematics and parameters defining a shock propagation into cellular bar after a rigid mass impact.

3.2 Analytical determination of displacement, velocity and stress

The conservation of energy of the whole system during the period t to t + dt (Figs. 3c and 3d) reads:

$$\begin{array}{l}\frac{1}{2}\left[M+{\rho }_{\text{DS}}{S}_0\left(\text{λ}+\text{dλ}\right)\right]{\left(\upsilon +\text{d}\upsilon \right)}^2-\frac{1}{2}\left[M+{\rho }_{\text{DS}}{S}_0\text{λ}\right]{\upsilon }^2\hfill \\ =\frac{1}{2}\left({\sigma }_{\text{stat}}+\sigma \right){\epsilon }_{\text{DS}}{S}_0\text{\hspace{0.17em}d}D.\hfill \end{array}$$(18)

The left hand side of equation (18) includes a part of the structure which is set in motion by the impactor. This part shares a fraction of the system kinetic energy with the impactor. A fraction of the initial kinetic energy is transferred to the structure in the form of internal (strain) energy.

Equation (18) is expanded with the help of equations (13), (12), (14) for dλ, dD and ρ_DS respectively. In addition the second order term ρ_c S₀ du dυ/ɛ_DS is neglected in the left hand side of equation (18). Finally, the non-linear differential equation describing the plate decelleration is:

$$\frac{\text{d}\upsilon }{\text{d}t}=\frac{{\sigma }_{\text{stat}}+C{\upsilon }^2}{B+Cu}.$$(19)

Where the constants B and C are defined by:

$$B=\frac{M}{S_0}$$(20)

$$C=\frac{{\rho }_c}{{\epsilon }_{\text{DS}}}.$$(21)

Equation (19) can be written in term of the displacement u only by using the second and first time derivatives of u:

$$B\frac{{\text{d}}^{2}u}{\text{d}{t}^2}+Cu\frac{{\text{d}}^{2}u}{\text{d}{t}^2}+C{\left(\frac{\text{d}u}{\text{d}t}\right)}^2+{\sigma }_{\text{stat}}=0.$$(22)

Note that equation (22) holds only if u, du/ dt are all both non constant functions of time.

Solving equation (22) with the initial conditions du/ dt = υ₀ and u = 0 at t = 0 gives :

$$u\left(t\right)=-\frac{B}{C}+\frac{\sqrt{C^2{\sigma }_{\text{stat}}^2\left(B^2+2BC{\upsilon }_{0}t-C{\sigma }_{\text{stat}}{t}^2\right)}}{C^2{\sigma }_{\text{stat}}}$$(23)

The time derivative of of u(t) gives the impactor velocity:

$$\upsilon \left(t\right)=+\frac{{\sigma }_{\text{stat}}\left(BC{\upsilon }_0-C{\sigma }_{\text{stat}}t\right)}{\sqrt{C^2{\sigma }_{\text{stat}}^2\left(B^2+2BC{\upsilon }_{0}t-C{\sigma }_{\text{stat}}{t}^2\right)}}$$(24)

Substituting equation (24) into equation (17) gives the formulation of dynamic stress with respect to time and impact conditions.

When u(t) and υ(t) are known the stress σ is directly calculated by equation (17). In addition, the following physical quantities of interest can also be determined :

– Stop time for the rigid plate. According to equation (24), if the time t equals

$$t_{stop}=\frac{B}{{\sigma }_{\text{stat}}}{\upsilon }_0$$(25)

the velocity of the plate is theoretically zero.

– Time for densification. At the particular time value t_DS the AS begins to densified and the corresponding strain is ε_DS given by equation (10). By equating the crushed length to ε_DSD₀, the time for densification is the the appropriate (physical) time root t_DS of the following equation:

$$u\left(t_{\text{DS}}\right)-{\epsilon }_{\text{DS}}{D}_0=0.$$(26)

4 Finite element model

A finite element model was developed to simulate the impact of the plate on the AS represented in Figure 1 and to evaluate the validity of the previous formulations (Sect. 3.2). The commercial explicit non-linear FE software RADIOSS^® (Altair HyperWorks) [15] is used to implement this model. The model is composed of a top falling plate, an auxetic structure supported by a fixed plate at y = 0. The auxetic is made of six cells in height Y-direction and five cells in horizontal X-direction. The characteristics dimensions for the AS and the plates are given in Tables 1 and 2.

An elastic-perfectly plastic aluminium is used for the auxetic material. A perfectly-elastic material is used for the two plates. The material properties are given in Tables 3. In the simulation, the crushing plate is supposed elastic with an elastic Young modulus equal to that of steel. The mass density is adjusted so that the crushing plate (with its imposed dimensions) has a mass equal to 4 g. This gives the plate a density of ρs = 1.11 × 10⁴ kg/m³.

In this simulation, Poisson’s ratio is not a critical parameter, as the behaviour of the plates is assumed to be elastic. As its value can vary between 0.28 and 0.34, we have chosen 0.33 to represent an average value for the steel. Poisson’s ratio also has no influence on the perfectly plastic elastic behaviour of the auxetic structure material. Thus, this material property is not a critical parameter in the simulation. For the sake of simplicity, we have set the same value of 0.33 for plates and aluminium.

The mesh size sensitivity has been analyzed with a coarse mesh and four gradually refined meshes (from 2 to 5 elements in the thickness, corresponding to sizes from t_thick/2 to t_thick/5). Variations of the impactor velocity results υ for these meshes are presented in Figure 4. It is found that the order of the deviation is of 1% between the one and five elements in thickness meshes. If the element sizes are equal or lower than t_thick/2, the results differ of less than 0.5%. Therefore, a mesh with two elements in the thickness (0.15 mm size) was retained.

The mesh of the three parts is composed of under-integrated 8-node hexaedral solid elements with physical hourglass stabilization. The plates are meshed with elements of 1 mm length. The auxetic structure is meshed with three elements in Z-direction. The fixed plate is constraint in all directions, the top plate is able to move only in Y-translation direction and, to avoid out-of-plane displacement, the auxetic is fixed in Z-translation direction. A penalty method general contact is applied to model all contact interactions at the plate/auxetic and auxetic/auxetic interfaces. The contact friction is not taken into account.

Fig. 4 Finite element impactor velocity results υ for several mesh sizes (from one to five elements in the thickness tthick) at initial velocity υ0 = 200 ms−1.

5 Results and dicussion

We consider the crushing of the AS by a falling plate of initial velocity υ₀ in Y-direction. The simulations are done for a 4 g mass plate with initial velocity υ₀: 35, 50 and 200 ms⁻¹. The dynamic response of the AS obtained through the FE simulation is expressed in terms of stress σ (left scale of the graph) and the instantaneous plate velocity υ (right scale of the graph) υs. time are plotted in Figures 5b, 6b and 7b, respectively. Figures 5a, 6a and 7a give respectively the AS deformation at selected times.

The plate is supposed always in contact with the AS. The stress σ_FE is calculated by the ratio of the contact force (crushing plate – AS) to the initial transverse (X-Z plane) area S₀. The dynamic stress σ is calculated with equation (17) by using the instantaneous velocity equation (24). The constant plateau stress σ_stat, equation (11), is also plotted. The times t_stop and t_DS are calculated with equation (25) and equation (26), respectively.

The calculated stress, equation (17), leads to accurate estimation of the AS average stress values for the different impact conditions.

In Figures 5b, 6b and 7b an initial peak stress can be observed for a short duration which causes an important absorption of energy. Hence, during this same time lapse, the FE velocity decreases sharply. This initial peak is not taken into account in the theory developed in Section 3.2. That is why FE and analytic velocity differ from each other at the beginning.

At 35 ms⁻¹, the FE simulation shows that the plate is completely halted 1.85 ms after impact. The time t_stop, calculated with equation (25), is 2.16ms, which indicates that the method presented in Section 3.2 overestimates the stop time for low initial velocity. After the period of initial peak stress, values of plate FE velocity and analytical velocity, determined with equation (24), are very close. However, once the velocity is lower than about 20 m s⁻¹ , the two values come apart from each other.

It can be seen, in Figure 5a, that the deformed region differs from that of the shock-wave model used to established the formulation of σ. In fact, the deformed region observed in the FE results corresponds to deformation modes, which depend upon υ₀ and the AS topology and geometry, (see for example Ruan et al. [9] for constant crushing velocity).

At 50 m s⁻¹ and 200 m s⁻¹ , the AS was not able to stop the plate before densification. When densification begins the plate velocity is not zero. The plate will continue to crush the structure while the stress increases up to the peak and rebound of the plate will occurs. The prediction of t_DS with equation (26), is reasonably accurate regarding FE results. For both, the evaluation of analytic plate velocity are very close to the FE velocity.

The formulations established in Section 3.2 seem to be more accurate regarding FE results for medium ( 50 m s⁻¹ ) and high initial impact velocity (200 ms⁻¹).

The alloy properties chosen in this article are close to those of the AL3024-T351 alloy. Khan et al. [16] have shown that if the strain rate does not exceed 1500 s⁻¹ , the behaviour of the material is not very sensitive to strain rate. As we have verified that in the simulations of the dynamic compression of the structure, the dε/dt rate remains below 100 s⁻¹, we have chosen not to take the strain rate into account in the behaviour of the material making up the AS.

6 Conclusion

In this paper, closed-form formulas have been proposed for predicting the crushing strength of architectured auxetic lattices loaded by the impact of a plate. The approach is based on energy conservation principle and shock-waves propagation analogy in a rigid, perfectly plastic, locking material model. It enables the analysis of periodic collapse of the structure. We have extended formulation in order to include more design parameters (specific geometrical ratios and/or material properties). Our analysis enables to theoretically predict the dynamic crushing strength. The formulation depends on the geometric and the material characteristics of the auxetic but also on the impact velocity.

The results of analytical calculations are consistent with the results of finite-element simulations. It is then possible to adapt the geometrical and material parameters in order to optimize the energy absorption ability for facing specific load cases. In-depth studies are in progress to understand the influence of the deformation modes on energy absorption.

Fig. 5 AS crushing simulations results for impact conditions: M = 4 g and υ0 = 35 ms−1. (a) AS deformation at t = 0.77 ms. (b) Stress and plate velocity υs. time.

Fig. 6 AS crushing simulations results for impact conditions: M = 4 g and υ0 = 50 ms−1. (a) AS deformation at t = 0.45 ms. (b) Stress and plate velocity υs. time.

Fig. 7 AS crushing simulations results for impact conditions: M = 4 g and υ0 = 200 ms−1. (a) AS deformation at t = 0.17 ms. (b) Stress and plate velocity υs. time.

Acknowledgement

The authors acknowledge the French Region Centre-Val de Loire for having provided the funds for the study under the project AxSur.

References

All Tables

All Figures

	Fig. 1 View of the complete auxetic structure (“AS”).
In the text

	Fig. 2 Re-entrant auxetic cell configuration.
In the text

	Fig. 3 Schematics and parameters defining a shock propagation into cellular bar after a rigid mass impact.
In the text

	Fig. 4 Finite element impactor velocity results υ for several mesh sizes (from one to five elements in the thickness tthick) at initial velocity υ0 = 200 ms−1.
In the text

	Fig. 5 AS crushing simulations results for impact conditions: M = 4 g and υ0 = 35 ms−1. (a) AS deformation at t = 0.77 ms. (b) Stress and plate velocity υs. time.
In the text

	Fig. 6 AS crushing simulations results for impact conditions: M = 4 g and υ0 = 50 ms−1. (a) AS deformation at t = 0.45 ms. (b) Stress and plate velocity υs. time.
In the text

	Fig. 7 AS crushing simulations results for impact conditions: M = 4 g and υ0 = 200 ms−1. (a) AS deformation at t = 0.17 ms. (b) Stress and plate velocity υs. time.
In the text