1. Introduction
Nowadays, magnesium alloys have been used in transportation, electronics, medical industries due to their good lightweight properties, machinability, corrosion resistance, shock absorption, dimensional stability and impact resistance. For example, in automotive applicants, low-temperature components like brackets, covers, cases of modern automotive are made of magnesium alloys [
1,
2,
3,
4,
5,
6,
7]. In biomedical applications, bioimplants, which are devices that replace the affected or damaged part of the human body and assist in the normal functioning of the human body with a high degree of physiological acceptance, were made by magnesium alloys [
8,
9,
10,
11]. For example, Amerinatanzi et al. performed a prediction of the biodegradation of magnesium alloy implants [
12]. With good application prospects and great potential, magnesium alloys have become one of the hot issues of new materials in the future [
13,
14].
Numerical simulation is an important procedure for the design and optimization of magnesium alloy structures. The reliability of numerical analysis depends on the accuracy of the material model. Developing a constitutive model of the magnesium alloys, however, is a rather challenging task compared to cubic metals because magnesium alloys have a hexagonal lattice structure, which affects the fundamental properties of these alloys. The hexagonal structure limits easy dislocation motion to the close-packed direction on the basal planes, and plastic deformation of the hexagonal lattice is more complicated than in cubic latticed metals like aluminum, copper and steel. A material model for magnesium alloys should be able to describe these anisotropy, asymmetry and temperature-dependent behaviors [
15,
16,
17,
18,
19,
20,
21,
22].
For quite some time, many researchers have been working on the simulation of magnesium behavior. Saniee et al. [
23] considered the flow behaviors of several magnesium alloys in tension and compression. On the other hand, Agnew et al. [
24] validated an elastoplastic self-consistent polycrystal model that was used to simulate the macroscopic flow curves and internal strain developments within the distinctly textured magnesium alloy samples. Staroselsky et al. [
25] also applied a constitutive model hexagonal close-packed (HCP) material deforming by slip and twinning to magnesium alloy AZ31B. In the application aspect, Zhu et al. [
26] put forward two optimization-based methodologies to calibrate material parameters for the application of AM60B magnesium alloy material model to structure component crush analysis. A three-dimensional finite element model (FEM) was established by Shen et al. [
27] to simulate the temperature distribution, flow activity, and deformation of the melt pool of selective laser melting (SLM) AZ91D magnesium alloy powder. Samuha et al. [
28] improved formability for the commercial magnesium AZ80 alloy through the application of the high-rate electromagnetic forming (EMF) technique. Ma et al. [
29] developed a temperature-dependent anisotropic material model combined with the Hill48 yield function for the warm-drawing of magnesium alloys. A new extrusion process for Mg–Al–Mn–Ca magnesium alloys using rapidly solidified powders produced very fine grains was also developed by Ma et al. [
30].
Kim et al. [
31] proposed a constitutive model that combined both Cazacu’06 yield criterion [
32] and Hill’48 yield criterion [
33] to describe the temperature-dependent asymmetric cyclic behavior of magnesium alloy with the concept of the dominant deformation mode twinning (T), untwinning (U) and slip (S) in 2013. Although effective, this theory was complex to be realized in a computation program.
Therefore, this paper develops a new material model facing practical industry application based on Kim’s theory [
31] in LS-DYNA. First, we performed a cyclic loading test and employed a one-shell element model to verify the accuracy of the new material model. Then, this new material model was applied to predict the deformation behavior of magnesium alloy made crash box. The simulation results showed good agreement with the experiment. Several factors were investigated, such as wall thickness, friction coefficient and initial deflection, to try to find which factor led to the deviation between the simulation results and the experiment results. It was found that the initial deflection, friction coefficient and thickness of the wall may be the considerable and reasonable factors that led to the error.
2. Materials and Methods
To describe the material properties of an HCP metal and alloys, Cazacu et al. proposed a yield function (Cazacu’06) [
32] that can capture the asymmetric yield surface of magnesium alloys. The asymmetry of the yield condition of magnesium alloy was caused by deformation based on twinning and untwinning. Schematic stress–strain curves in tension–compression and then compression–tension is shown in
Figure 1.
For the accurate simulation of the deformation of the HCP material, considering the effect of the twinning was important. Kim et al. [
31] introduced an idea that considered three modes as twinning, untwinning and slip, to describe the cyclic loading behavior of magnesium alloys by using the combination of the Hill’48 [
33] and Cazacu’06 [
32] yield functions. To conduct a finite elment analysis (FEA) simulation based on the constitutive model, which considered twinning, untwinning and slip modes, the introduced model was implemented to commercial FEA code LS-DYNA by using a user material subroutine. The new material model focused more on the practical application. Hence, some parts of Kim’s model were simplified to reach high-efficiency with acceptable accuracy. The Hill’48 [
33] yield function was used for describing the slip-dominant deformation. For simplification, plane stress conditions and only transverse anisotropy were assumed. Then, its effective stress was given by:
where
was the average Lankford value.
,
,
and
were normal components of the Cauchy stress in the loading direction. The slip deformation was dominant in the tensile-loading condition without twinning history. On the other hand, the twinning deformation was dominant in the compression loading condition. Cazacu’06 [
32] yield function was used for the twinning and untwinning dominant deformation. Assuming the plane stress condition and planar isotropy; then, effective stress is given by:
where
was the exponent of the yield function, and
was the parameter of stress ratio.
and
were the in-plane principal values of the Cauchy stress. In this study, exponent
a was set 2.0.
was the parameter of the stress ratio, which was defined by using compressive yield stress
and tensile yield stress
for the twinning and untwinning dominant deformation in Equation (3):
During the tensile deformation with twinning history, if the following condition Equation (4) was satisfied, the constitutive equation was computed as untwinning deformation was dominant.
where
C was the scale value,
was effective plastic strain by twinning deformation and
was effective plastic strain by untwinning deformation.
Then, how the new material model works in the plane stress state is shown in
Figure 2a. The thick gray line shows the normalized yield surface used in this research, which is the combination of the Hill’48 [
33] and the Cazacu‘06 [
32] yield functions. Numerical_TEST1 and Numerical_TEST2 were the verification points at different stress ratios.
Let
be the mean stress defined as below:
When the tension,
and
, yield stress was under slip control, the Hill’48 [
33] yield function was used for elastic-plastic calculation. When compression,
, yield stress was controlled by twinning, and the Cazacu’06 [
32] yield function was used for elastic-plastic calculation. When tension,
and
, yield stress was controlled by untwinning, so in the elastic-plastic calculation, the Cazacu’06 yield function could be used.
As an example of the validation of the new material model, tension–compression–tension loading tests using AZ91 and numerical calculations were carried out. Their results are shown in
Figure 2b. The strain was set as 0% → +3% → −3% → +3% and strain rate was 0.1. For the comparison, the calculation result of the fundamental elastic-plastic model (von Mises, V-M) is also included in
Figure 2b. The new material model considered two stress–strain curves; one was a stress–strain curve for slip deformation, the other one was a stress–strain curve for twinning deformation as input data. The yield stress during untwinning dominant deformation was described by using yield-stress based on twinning and the stress ratio parameter
β. This result was consistent with the result given by Kim’s model. In this study, parameter
β was treated as a constant value for the simplification. Specifically, parameter
β was set as 0.5 in the validation example; how this was determined is introduced in
Section 4. When the stress ratio
β could be defined as nonconstant data, the relation between twinning and untwinning dominant yield stress could be reproduced more flexibly.