Shock wave speed and stress-strain relation of aluminium honeycombs under dynamic compression

The propagation of layer-wise crushing bands in cellular materials under dynamic impact can be described by the plastic shock wave model. A cell-based finite element model of irregular aluminum honeycomb is constructed to carry out several constant-velocity compression tests. The shock wave speed is obtained by the one-dimensional stress distribution in the specimen along the loading direction. The relation between the shock wave speed and impact velocity is obtained and analyzed. It is found that the relation tends to be linear with the increase of the impact velocity. But the shock wave speed tends to be a constant value with the decrease of the impact velocity. A piecewise model is proposed to describe the dynamic stress-strain relation of aluminum honeycombs based on a piecewise hypothesis of the relation between the shock wave speed and the impact velocity together with the one-dimensional shock wave theory. Different stress-strain relations corresponding to different impact velocity regions and different deformation modes are obtained.


Introduction
Cellular materials have a typical characteristic of deformation localization under dynamic crushing due to the complexity of microstructures. Reid et al. put forward the concept of "structural shock wave" in the impact test of the one-dimensional ring system [1], and confirmed the shock wave propagation characteristics of the crushing deformation bands in the dynamic impact test of wood [2]. Liu et al. [3] observed the layer-wise propagation of deformation crushing bands in the dynamic impact of honeycomb and proposed three kinds of deformation modes, namely homogeneous mode, transitional mode and shock mode. Zou el al. [4] believed that the shock front in the honeycomb has a width of about cell size under high velocity impact. Liao el al. [5] proposed a local strain field calculation method to characterize the propagation of one-dimensional shock wave in the aluminum honeycomb and obtained shock wave speed. The researchers have proposed a series of one-dimensional shock models to describe the shock wave propagation in cellular materials, such as R-PP-L (rate independent, rigid -perfectly plastic -locking) model [2,6,7], E-PP-R (elastic-perfectly plastic-rigid) model [8][9] and R-PH (rate independent, rigid-plastic hardening) model [10,11]. But most of the shock models are obtained based on quasi-static nominal stress-strain curves. It is difficult to experimentally get dynamic stress-strain relationship directly. In fact, the nominal stress-strain curve loses the physical meaning due to the deformation localization and the stress enhancement effect under dynamic impact [3]. Zheng et al. [11] obtained the stress-strain state points behind the shock front by analysing the three-dimensional finite element model of aluminum foam, and put forward a D-R-PH (dynamic, rate independent-rigid-plastic hardening) model to characterize the dynamic stress-strain behaviour. The dynamic stress-strain states of cellular materials were obtained experimentally and numerically by Barnes et al. [12] and Sun et al. [13] respectively. But these stress-strain states are obtained under high impact velocity (shock mode), the stress-strain state under moderate velocity (transitional mode) is still not very clear. In addition, the stress behind the shock in the literature was mostly obtained by the nominal stress on the impact end, and it lacks direct representation of the local stress information inside the material. Thus in this paper, the stress distribution of honeycomb under constant-velocity compression is studied numerically and the behaviour of shock wave propagation under moderate-velocity impact is analysed.

Cell-based finite element model
An irregular honeycomb specimen is constructed by the Voronoï technique from Ref. [14], see Fig. 1. The irregularity is 0.3, the number of cells is 500. The relative density ρ 0 /ρ s is 0.1, where ρ 0 is initial density of the specimen，ρ s is the density of the base material. The length and width of the specimen are 50 mm and 50 mm. The mean cell size of the specimen is defined as the diameter of the circle equal to the average cell area, which is about 2.5 mm.
The finite element code ABAQUS/Explicit is employed to perform virtual tests. Cell walls are modelled with S4R (a 4-node doubly curved, reduced integration, hourglass control, finite membrane strains) shell elements. The material of cell walls is assumed to be elastic, perfectly plastic with Young's modulus E = 66 GPa, Poisson's ratio ν = 0.3, yield stress σ ys = 175 MPa, and density ρ s = 2700 kg/m 3 . The in-plane size of the element is about 0.2 mm determined through a mesh sensitivity analysis. All nodes are constrained in the out-of-plane direction to simulate a plane strain situation.
The specimen is placed between two rigid plates, one is moving with a constant velocity V along the X direction and the other is fixed.

One-dimensional stress distribution
For the honeycomb subjected to constant-velocity compression, the one-dimensional approximation is applicable. So the local stress information can be characterized by cross-sectional engineering stress. The method is detailed in Ref. [15]. The stress distribution at the impact velocity of 200 m/s is shown in Fig. 2. At any nominal strain, the stress will jump at a certain position which moves to the right of the specimen with the increase of the nominal strain. This shows that the shock wave propagation does exist in the specimen under high velocity impact, and the position of the stress jump is the shock front position. The stress enhancement phenomenon can be observed obviously from the figure. Across the shock front, the stress increases rapidly from ahead of the shock front to behind the shock front and almost keeps invariable.

Shock wave speed
The stress distribution gradient can be used to determine the position of the shock wave front, as shown in Fig. 3. The stress gradient reflects the changing degree of one-dimensional stress in the loading direction. The position of the stress changes fastest when the absolute value of the stress gradient reaches the maximum value. This position is the shock front position at the current time. The relationship between the shock front position and the impact time under different impact velocities is shown in Fig. 4. Under constant-velocity compression, the shock front position is linear with the impact time. Therefore, the shock wave speed can be predicted by the slope of the relation from linear fitting. The greater the impact velocity is, the greater the slope is and the greater the shock wave speed is.

.1 Piecewise hypothesis of shock wave speed
In the previous section, we have directly obtained the relationship between the shock wave speed and the impact velocity under the high velocity impact. But under a moderate impact velocity, the shock wave propagation phenomenon can still be observed [3]. The stress gradient distribution is still used to determine the shock front position and the relationship between the shock wave speed and the impact velocity is shown in Fig. 5. It is found that when the impact velocity is very high, the shock wave speed has a linear relationship with the impact velocity. However, with the decrease of the impact velocity, the shock wave speed enters a transitional stage and tends to a constant. Therefore, we can divide the relationship between the shock wave speed and impact velocity into three regions, namely, uniform-velocity region, transitional region and linear region. It is conceivable that these three regions correspond to the three types of deformation modes of cellular materials, namely, homogeneous mode, transitional mode and shock mode.
In the stress wave theory [16], the mass and momentum conservation equations across the shock front are respectively ( ) and ( ) From the combination of the two equations, the following equation can be obtained, i.e., where V s (t) is shock wave speed, {σ B , ε B , v B } are physical quantities behind the shock front，{σ A , ε A , v A } are physical quantities ahead of the shock front ， respectively are stress, strain and velocity. Under constant-velocity compression scenario, the conditions of ε A = 0, v A = 0, v B = V and σ A = σ 0 are satisfied based on the assumption of one-dimensional shock theory. Thus Eq. (3) can be rewrite as Therefore, if we can get the relation between the shock wave speed V s and the impact velocity V, we can get the relation between the strain behind the shock front and the impact velocity and the shock wave speed based on the mass conservation relation (1). By eliminating the impact velocity term in the formula (4), the relation between the stress and strain can be obtained. It is difficult to directly use a continuous function to characterize the relationship between the shock wave speed and the impact velocity in Fig. 5, thus we consider a piecewise hypothesis of the relation based on the characteristic of the relation between the shock wave speed and impact velocity, shown in Fig. 6. The function g is the first derivative of the shock wave speed. The function g is determined by the following equation, i.e., where a and b represent the two critical impact velocities between the deformation modes. Then we integrate the function g and consider the condition that the shock wave speed approaches a constant V 0 in low-velocity uniform region. The integral constant is determined by the continuous condition and boundary condition. The piecewise function of shock wave speed and impact velocity is expressed as

Development of piecewise model
Stress-strain relation in different velocity regions are obtained from the combination of Eq. (6) and the one-dimensional shock wave theory. The transitional region is derived from the high-velocity and low-velocity regions. When V < a, the constant shock wave speed V s = V 0 = 55.83 m/s is obtained by the average of the low-velocity region. In this case, the deformation mode of the honeycomb is quasi-static homogeneous mode, so it may be inaccurate to call the shock wave. In order to continue to use shock wave theory to solve the problem, it is supposed that a wave with the speed of V 0 can make stress and strain jump from initial collapse state under a small strain case. The formula V = V 0 ε can be obtained by the mass conservation relation (1) and then the stress-strain relationship in this low-velocity uniform region is obtained from the combination of Eq. (4), i.e., This is a linear hardening process. The critical impact velocity a is about 20 m/s, which is determined by the stress uniformity index [3]. The corresponding characteristic strain ε 1 = 0.358, the corresponding characteristic stress σ 1 = 1.05 MPa.
The stress behind the shock front can be obtained by the stress distribution in Section 3.1. The strain behind the shock front can be obtained by the strain distribution by the local strain calculaiton method [5]. Fitting the equation above leading to the parameters c = 1.037, d = 13.92 m/s. And then b = 61 m/s, k = 0.0254 (m/s) −1 is calculated by the combination of the critical impact velocity of a.
When a < V < b, i.e., in the moderate-velocity region, the relation of shock wave speed and the impact velocity is formulated as The stress-strain relation is obtained from Eq. (4), which is formulated as In order to ensure the hardening relation of stress to strain, the negative sign is considered and the characteristic stress σ 2 = 2.01 MPa.
Therefore, the corresponding dynamic stress-strain relation is respectively obtained in the low-velocity uniform region, moderate-velocity transitional region and high-velocity linear region. Here, we call this dynamic constitutive model D-R-P (dynamic, rate-independent, rigid-piecewise) model, which is obtained from the piecewise hypothesis of the relation between shock wave speed and impact velocity together with the one-dimensional shock theory. Dynamic stress-strain relation obtained by different models and the quasi-static stress-strain curve are compared in Fig. 7. It can be seen from the graph that the D-R-P model and D-R-PH model are all located on the right side of the quasi-static stress-strain curve when the strain is relatively large, which obviously shows dynamic plastic hardening feature. Under the shock mode with high impact velocity, the stress-strain relationship obtained by the two models is very close. This is because the two models are both based on the linear relationship between the shock wave speed and the impact velocity, but the parameters have slight differences. It is different that when the strain is relatively small, especially when it is less than the characteristic strain ε 1 , the stress-strain relationship obtained from the D-R-P model is linear hardening, which is located above the quasi-static stress-strain curve and the D-R-PH model. The locking strain ε L in the R-PP-L model is located between the two characteristic strains of the D-R-P model, and the plateau stress is located between the quasi-static initial crush stress and the characteristic stress σ 1 . Under this small-strain condition, the D-R-PH [11] model is conservative, and the predicting stress-strain state is relatively small. This is because the basic relationship of the model is based on the dynamic stress-strain state under high-velocity impact, lacking the research on the stress-strain state at low and moderate impact velocity. The D-R-P model proposed here considers the three cases of shock wave velocity based on the eyesight of shock wave propagation. A relatively complete piecewise model that can describe the dynamic mechanical properties of honeycomb is obtained.

Verification of D-R-P model
The relationship between the impact stress and time obtained by the D-R-P model and the numerical results at three different impact velocities is shown in Fig. 8. It can be observed from the graph that whether the impact velocity is low, moderate or high, the stress behind the wave front predicted by the D-R-P model can well coincide with the platform segment of stress on the impact end, which indicates that the D-R-P model can be used to characterize the mechanical properties of honeycomb material under dynamic compression.