Development of impact attenuator analysis tools in crash scenario using Euler method and finite element analysis

The problem considered in this work is the development of simulation method for simulating car crash which utilizes simple car—impact attenuator model developed in MATLAB. Usually, car crash simulation is done using full finite element simulation which could take hours or days depending on the model size. The purpose of proposed method is to achieve quick results on the car crash simulation. Past works which utilizes simple car—impact attenuator model to simulate car crash use continuous time model and the impact attenuator parameter is obtained from the experimental results. Different from the related works, this work uses discrete time model, and the impact attenuator parameter is obtained from finite element simulation. Therefore, the proposed simulation method is not only achieving quick simulation results but also minimizing the cost and time in obtaining the impact attenuator parameter. The proposed method is suitable for parametric study of impact attenuator.


Introduction
The road accident involving vehicles is a primary concern for cities all over the world. The accident happens because most drivers cannot avoid certain situations that result in a crash between a vehicle and other objects [1]. From the year 2000 to 2016, WHO (World Health Organization) data shows that the number of deaths per 100,000 population resulting from road accidents worldwide is increasing [2]. While the trend of road accidents may not be easy to alter due to increasing number of vehicles, mitigating the effects of an accident is possible, potentially reducing the number of fatalities.
Technologies to avoid or mitigate the effect of accidents or crashes have been developed. The works [3][4][5][6][7] are categorized as active technology which assists the driver to prevent an accident. The works [8][9][10][11] are categorized as passive technology where the structure of the vehicle is designed to protect the passengers and driver from fatal injury, minimizing fatalities. The ability of a vehicle to protect the passenger during a crash from fatal injury is called crashworthiness [12]. Crashworthiness is one of the most important aspects of vehicle design [13].
The two terms, impact attenuator (e.g., [8]) and crashbox (e.g., [11]) are technically similar since they are structures designed to absorb the energy during impact incidence. In this work, we will use the term impact attenuator. In automotive industry, the impact attenuator is classified as a passive safety system designed to reduce occupant injury during crashes [14]. An example of the study of impact attenuators in a formula race car can be found in the work of [15]. Another application of the impact attenuator is a rail guard in the public street and an external device carried by construction trucks [16]. In evaluating the effectiveness of the impact attenuator during the design stage, crash analysis is needed. Finite element analysis (FEA) using Finite Element suite such as ANSYS is a popular method in performing crash analysis. The works [15,17] use LS-DYNA software to simulate the impact attenuator for FSAE (Formula SAE, a student design competition in building a formula car) application and [13] uses LS-DYNA to simulate impact attenuator for racing application, [18,19] perform the crash simulation on FSAE model using RADIOSS Code software to evaluate the performance of the impact attenuator.
The crash analysis using FEA saves a lot of time and resources compared to an actual crash test on the actual vehicle in evaluating the impact attenuator, as suggested by [20]. However, FEA is usually done using high-performance computers to keep computational time minimum. Even with high-performance computer, crash simulation using FEA usually takes hours even days to complete the simulation. When high-performance computer is not available, static analysis can be done to estimate the crash analysis result, such as done by [1]. However, the static analysis does not consider the dynamic behavior of the model, which makes the static analysis results to be not accurate.
This paper proposes a numerical method to simulate a crash incidence on impact attenuator model by considering the model's dynamic behavior. Hybrid simulation using MATLAB and FEA software (ANSYS) approach is done therefore the simulation time can be minimized. Related works can be found in [21] where crash is modeled using two bodies and one isolator (spring and damper). In [21], the genetic algorithm is used to compute the model parameter (spring stiffness). The work [22] uses power law model to model the crash and experimental observation is done to obtain the impact attenuator parameters. The work [23] uses simple model of mass-spring to simulate the crash and the spring stiffness is obtained from full-scale experimental data analysis.
Similar to the works [21][22][23], the mathematical model of the impact attenuator in this work is built using the first principle as a continuous-time model, i.e., it comprises of mass-spring with an additional damper model. However, different from [21][22][23], we convert the continuous-time model into the discrete-time model using Euler method and we show the exact discrete-time model which is later coded into MATLAB's routine to simulate the problem. Additionally, different from the related works, the properties of the impact attenuator (the spring stiffness) are updated every time step by using the stiffness-strain data generated using quasi-static analysis in ANSYS Mechanical APDL which includes the elastic-plastic behavior of the impact attenuator. In other words, there is no need to perform an experiment to obtain the parameters. By using the proposed simulation method, crash simulation can be done relatively fast, including the time required to obtain the impact attenuator parameters. The proposed method is suitable to be used in the design stage, where a parametric study of the impact attenuator is done. The problem formulation and modeling of the impact attenuator were firstly demonstrated in this paper. The simulation results of the developed model are then comprehensively discussed.

Problem statement
In this work, a simple discrete-time model of impact attenuator crash will be developed. The physical phenomenon considered for the modeling is given in Fig. 1.
The car is modeled as a rigid body with mass M, the impact attenuator is attached to the car, and the impact attenuator has spring constant k and damping factor c. During the crash, the impact attenuator hits the wall and deformed. From the simulation tools that will be developed, the output variables are the deformation of the impact attenuator and the acceleration response of the car. If the deformation of the impact attenuator is sufficiently large, then it implies that the impact attenuator can efficiently absorb the impact energy. Hence, the acceleration of the car during the crash will be hugely reduced. This will lead to minimum injury of the passenger inside the car. In other words, the acceleration response of the car during the crash can be used to evaluate the effectiveness of the impact attenuator. In the model, the crash is assumed to happen in one direction.
The impact attenuator, when colliding with the wall, is expected to have large deformation. Thus, both elastic and plastic deformation of the impact attenuator must be considered. To this end, the spring constant k in the model is not only capturing the elastic behavior but also the plastic behavior of the impact attenuator. Meanwhile, the damping factor c is assumed to be constant. Damping is a function of energy dissipation property of the material (in this case the impact attenuator). In other words, we assume that the energy dissipation rate of the impact attenuator is constant throughout. To make such that the stiffness k capture the elastic-plastic behavior of the impact attenuator, the value of stiffness k must be updated according to the deformation level of the impact attenuator.
The crash incidence considered in this paper is the crash without rebound. That is, when the impact attenuator deforms plastically, we assume that the impact attenuator retains its deformation when there is no more force acting on the impact attenuator (the car loses all its kinetic energy).

Methods
In general, the proposed numerical method consists of two processes: 1. Mathematical modeling. This process includes the derivation of the discrete-time model. 2. Obtaining the stiffness-strain curve of the impact attenuator. In this process, a quasi-static compression analysis will be done to the impact attenuator and the analysis is done using FEA software (ANSYS). The stiffness-strain curve is used in the discrete-time model simulation in MATLAB. The stiffness obtained is the effective stiffness or the change in force divided by the change in length. 3. Simulation of the discrete-time model using MATLAB. The simulation in MATLAB is done to obtain the acceleration response of the car body for different impact attenuator designs when the impact attenuator crashes into the wall.

Mathematical model
In this work, the system in Fig. 1 is modeled as two degrees of freedom system (two masses), where the first body is the car's body with mass M and the second body is the shaker with mass M s , where M s ≫ M. The displacement of the car is denoted as y and x denotes the displacement of the shaker. A force F is given to the shaker such that € x or the acceleration of the shaker follows the half-sine pulse with peak P and pulse duration T (see Eq. (1)). Below is the formula used to construct the force F which is based on the work [24].
The impact attenuator is modeled as a spring and damper, and these are attached to the car's body. The crash model and the acceleration pulse can be seen in Fig. 2. Note that the direction (vertical or horizontal) does not matter for the simulation. The impact modeling in this work models the impact as a half-sine pulse acceleration at the end of the impact attenuator (the point in contact with the wall in the real impact model). Since we are interested only in the response of the car's body and the shaker's mass will be set to be very high in the simulation, the model in Fig. 2a is considered as base excitation problem. Base excitation problem as an approach to study the behavior of system under impact is commonly used such as found in the works of [25,26].
From Fig. 2a, assuming x > y, the free-body diagram can be constructed, as shown in Fig. 3.The equation of motion in matrix form can be expressed as where ϵ is the strain of the impact attenuator and the stiffness k(ϵ) is expressed as a function of the impact attenuator's strain. The state space equation can be expressed aṡ where y is the displacement of the car's body, and u(t) = F(t) where F(t) is given by (1). Meanwhile, matrices A(ϵ) and B will be defined shortly.
To simulate the crash, there are two conditions considered in the model: (i) The first condition is when the car is still compressing the impact attenuator, and the impact attenuator is still in the elastic region. (ii) The second condition is when the car compresses the impact attenuator until the impact attenuator reaches the plastic region. In this region, we assume that the impact attenuator has permanent deformation.
The simulation considers that condition (i) is continued by condition (ii) and then the car stops. Therefore, the matrices A(ϵ) and B in the state space equation are: The matrix A is a function of ϵ since it will be updated according to the strain of the impact attenuator.
To simulate the problem in MATLAB, the continuous-time model (3), is converted into the discrete-time model: where k is the time step, h is the sampling time, and I is the identity matrix. To obtain the time t in seconds from the discrete time simulation, we multiply k and h or t = kh.

Quasi-static compression test analysis of the impact attenuator
In this work, the impact attenuator is a hollow box made from Aluminum 6082-T6. The reason of using a hollow box as the impact attenuator model is so that we can perform a parametric study on its stiffness by changing the thickness of the shell, while the outer dimension of the box remains the same. To obtain the stiffness data to be used in the simulation, a compression test analysis is done to the hollow box and the compression is done with a strain rate of 0.1/s. The compression test analysis is done using ANSYS Mechanical APDL with the following steps: -The Aluminum 6082-T6 compression stress-strain curve is first obtained from the work of [27]. The curve is obtained from compression test (experiment) on standard compression test specimen.
-The stress-strain curve from the above step is keyed into the material properties of the finite element model of the hollow box.
-Static analysis in ANSYS is done by giving one end of the hollow box a constant velocity such that the strain rate becomes 0.1/s. The strain rate is a function of the velocity and the length of the hollow box.
-Stiffness-strain curve can be obtained easily using post-processing menu in ANSYS Mechanical APDL.

Simulation set-up
There are two types of impact attenuator that will be simulated in this work. Both are hollow box, one with thickness of 3 mm and another with thickness of 8 mm. The hollow box is first created in SolidWorks and then imported to ANSYS. Figure 4 shows the isometric view of the hollow box. The impact attenuator model is 2 m wide, 1.25 m high, and it has a length of 3 m. The damping factor, c, is assumed to be constant and it is computed from the loss material coefficient of Aluminum 6082-T6. We assume the loss material to be 0.05. Let η denote the loss material coefficient. Thus, the damping factor is as follows: where M is the mass of the car and k is the stiffness of the impact attenuator on the elastic region. The damping factor is varied to see the effect of damping to the acceleration response of the car during the impact. The shock duration or pulse duration T is set to be 0.1 s (see [28]) and the peak of the acceleration is set to be 300 m/s 2 . With this pulse duration and peak acceleration, the considered initial velocity before impact is 48 km/h. This number is within the range of recommended vehicle's crash test (see [29]). Meanwhile, the shaker mass is set to be M s = 1 × 10 6 kg and the car's mass is set to be M = 1045 kg. The compression test data keyed into material library of ANSYS for the purpose of compression test analysis is given by Fig. 5 [27].  The boundary condition of the compression analysis is given in Fig. 6, the right side of the box is pushed to the left with constant velocity while the left side of the box is constrained in all direction. To illustrate the deformation of the box, the contour plot of the displacement of the box in the direction of the crash is given in Fig. 7. Figure 7 shows the contour plot of the displacement when the free side of the box is displaced 0.8004 m.
To obtain the stress-strain curve of the box, the stress and strain is first calculated from displacement and force data on the rightmost node. To get the stress data, σ, the force is divided by the box cross section area or σ ¼ F A , where F is the force at the node and A is the cross-section area of the box. Meanwhile, to get the strain data, ϵ, the displacement of the same node is divided by the initial length of the box or ϵ ¼ disp L 0 , where disp is the displacement data and L 0 is the initial length of the box. To obtain the stiffness of the box, k, the force F is divided by the displacement of the same node or k ¼ F disp . The force and displacement data can be easily obtained from the TimeHist Postpro menu in ANSYS. The stress-strain curve and the stiffness-strain curve of the 3 mm box and 8 mm box are given in Figs. 8 and 9, respectively. We can clearly see that the 8 mm box is stiffer than the 3 mm box. On the other hand, from Figs. 8 and 9, we can see that the elastic region for both boxes is in between 0 and 0.02 strain.

Simulation result
In this section, simulation results of the car with impact attenuator crashing the wall will be presented. There are three sets of simulation that will be done as shown in Table 1. The three sets of simulation will be presented in the following subsections.
Result of loss material coefficient = 0.05 The acceleration response of the car when using η = 0.05 is given in Fig. 10, while the strain of the impact attenuator is given in Fig. 11. The strain responses in Fig. 12 show that the plastic region starts at 0.6 s of the response (refer to Figs. 9 and 10 that the plastic region starts after strain of 0.02). Meanwhile, from Fig. 10, we can see a spike around 0.6 s; and this is due to the change in material behavior from elastic to plastic. From Fig. 10, we can see that the acceleration response is way below the shaker acceleration. This means that both impact attenuator design (3 mm and 8 mm thickness)  can suppress more than twice of the input acceleration. Although the 3-mm-thick impact attenuator results in slightly lower acceleration response, the 8-mm-thick impact attenuator has lower strain; or that its deformation is smaller. The final acceleration value (zero) in Fig. 10 and final strain value (constant) in Fig. 11 indicate that the car finally stops after the impact. Figure 12 shows the graph of mass times acceleration versus strain for both impact attenuator. The area under the graph is the energy absorbed per unit length. The energy absorber per unit length for 3 mm box is 1:523 Â 10 4 J m , while the energy absorbed per unit length for 8 mm box is 1:512 Â 10 4 J m .
Result of loss material coefficient = 0.5 The acceleration response of the car when using η = 0.5 is given in Fig. 13, while the strain of the impact attenuator is given in Fig. 14. The strain responses in Fig. 14 show that the plastic region starts at 0.6 s of the response (refer to Figs. 8 and 9 that the plastic region starts after a strain of 0.02). Meanwhile, from Fig. 13, we can see a spike around 0.6 s; and this is due to the change in material behavior from elastic to plastic. The peak acceleration response for η = 0.5 is higher compared to the case when η = 0.05, but the acceleration becomes zero (the car stopped) at a shorter time. Compared to the case when η = 0.05, the acceleration response when η = 0.5 oscillates less. This is the result of a higher damping factor. From the strain response in Fig. 14, we can see that the deformation of impact attenuator when η = 0.5 is smaller compared to the impact attenuator when η = 0.05. The final strain is constant because the impact attenuator is in the plastic region. Figure 15 shows the mass times acceleration versus strain for both impact attenuators. For 3 mm

Result of loss material coefficient = 1
The acceleration response of the car when using η = 1 is given in Fig. 16. We can see that the acceleration response is higher compared to cases when η = 0.5 and η = 0.05. The higher acceleration response peak is due to higher damping factor of the impact attenuator. This situation is related to the force transmissibility in the base excitation problem. The force transmitted from the shaker to the car will be higher when the damping factor is higher, which results in higher acceleration response.
The strain response of the impact attenuator when using η = 1 is given in Fig. 17. The strain is also lower compared to the case when using η = 0.5, but the strain shows that the impact attenuator is in the plastic region. Meanwhile, Fig. 18 shows the mass times acceleration versus strain for both impact attenuators. For 3 mm box, the energy absorbed per unit length is 1:0298 Â 10 4 J m , while for 8 mm box is 1:0116 Â 10 4 J m . These also show that higher damping results in lower energy absorption capacity during the impact.
The summary of the results is given in Figs. 19, 20, and 21. From Fig. 19, we can see that the loss material coefficient increases the acceleration response peak. We can also see that the 3-mm and 8-mm-thick impact attenuator has no to little difference in attenuating the acceleration when η is small, but the difference is clearer when η is larger. From Fig. 20, we can see that when η is small, the difference of final strain between the 3 mm thick and 8 mm thick is clearer compared to larger η. The energy absorbed is the area under the curve of mass× acceleration versus strain. The total energy absorbed cannot be concluded from the peak acceleration and the final strain value since it depends on the time domain response of both acceleration and strain. Thus, the dynamics of the car-impact attenuator affect the energy absorbed. From Fig. 21, we can see that the energy absorbed between the 3 mm and 8 mm box when η is small has very little difference. Meanwhile, the gap is wider when η is large. A good impact attenuator must have high energy absorption during the crash. Thus, a high damping factor is not desirable for such purpose. On the other hand, lower stiffness is not necessarily yielding a better impact attenuator. The 8-mm-thick impact attenuator with elastic stiffness of 2:088 Â 10 6 N m has very little difference with 3-mm-thick impact attenuator with elastic stiffness of 1:64 Â 10 6 N m in terms of acceleration attenuation.

Conclusions
This paper presents a method of simulating car crashes with a simple car-impact attenuator model. The model is based on two degrees of freedom system, where the two masses represent the car and the shaker. The shaker gives the impact force to the model. To capture the dynamics of the impact, the stiffness of the impact attenuator model is updated by using the stiffness-strain data generated from compression test analysis done in finite element software, ANSYS. The study shows that the method can be used to perform a parametric study on the impact attenuator design, with acceleration response and impact attenuator strain as the two variables being monitored. Using the simulation method, we study two impact attenuator models with different thicknesses; and we show that a hollow box with 8 mm thickness attenuates the acceleration response quite well while retaining lower deformation.
The parameter that is not considered to be changing due to plastic deformation is damping. It will be an interesting future work to include the effect of plastic deformation on the damping value of the impact attenuator.