Plastic Deformation Behavior of Metal Materials: A Review of Constitutive Models

: The deformation behavior of metal materials in plastic forming is intimately related to deformation conditions, which are greatly affected by deformation rate, forming temperature, and plastic variables. Macroscopic mechanical properties research is an important basis and technical means to analyze the process parameters and deformation process of metal plastic forming. Therefore, to reveal the inﬂuence mechanism of macroscopic mechanical properties of metal materials, and establish material constitutive models under different deformation conditions, it is of great signiﬁcance to choose reasonable forming parameters and prevent forming defects. There are substantial variances in the macroscopic mechanical characteristics of different materials in the deformation process. In order to accurately predict its deformation behavior, the phenomenological constitutive model, the microscopic constitutive model reﬂecting the microscopic deformation mechanism, and the artiﬁcial neural network constitutive model based on the neural network were constructed respectively on the basis of macroscopic mechanical tests and microscopic microstructure tests. On the basis of the existing research results, the advantages and disadvantages of phenomenological constitutive model, microscopic constitutive model, and neural network constitutive model are compared and analyzed, respectively. The research results of this paper will provide support for the selection of constitutive models for reasonably predicting the deformation behavior of metal materials.


Introduction
In recent years, with the rapid development of aviation, space flight, automobile and other fields, higher requirements have been put forward for the forming quality and service life of metal parts.The metal plastic forming process is the main forming technology for thin-walled and key stressed parts.Therefore, the formulation of suitable process parameters based on the plastic deformation behavior of materials is the key to the quality control of metal plastic forming [1,2].
In the process of metal plastic forming, the flow behavior of material is highly affected by deformation temperature, deformation velocity, and deformation amount.Especially in the thermoplastic forming process, the plastic deformation behavior of the material is the result of the combined action of work hardening and temperature softening, and is greatly affected by strain, strain rate, and forming temperature.On the one hand, under the conditions of given thermodynamic parameters, the deformation of materials is accompanied by microstructure evolution; on the other hand, the change of microstructure during plastic deformation also causes the change of flow stress of materials.Therefore, revealing the plastic deformation mechanism of metal materials under complicated deformation conditions is of great significance for the formulation and optimization of reasonable plastic forming process parameters [3][4][5][6][7].
In order to accurately predict the flow behavior of metal materials under deformation, a large number of studies have been carried out in recent years, and constitutive models have been established on this basis; the constitutive models can describe the effects of strain rate, forming temperature, strain and strain rate on the deformation behavior of materials during deformation.However, due to the obvious nonlinear characteristics of the plastic deformation behavior of the material during the plastic deformation process, it is difficult for a single constitutive model to fully describe the influence of all deformation parameters [8].
Therefore, it is the focus of research on material deformation behavior to create new constitutive models or alter existing constitutive models with the use of macroscopic and microscopic test data.Nonlinear flow behavior of different metals or alloys under many impact factors consider the particular influence of deformation process parameters.
In the plastic deformation process, the change of microstructure is the main reason for the change of macroscopic mechanical properties under different deformation conditions.When the deformation reaches a certain degree, with the increase of deformation, the microscopic damage of the material continues to accumulate and expand [8], which makes the macroscopic mechanical properties of the material change significantly.
The accumulation of damage during plastic deformation is the main cause of forming failure, which will also lead to premature fatigue failure of parts.Therefore, the prediction of plastic deformation damage accumulation is also an important factor to be considered in the design of forming process parameters.The damage evolution of metal materials during plastic deformation can be separated into two components: one is the growth of the original hole of the material, and the other is the nucleation and growth of new holes during deformation.Based on continuum damage theory, the accumulation of stress-strain is the main factor of plastic deformation damage.With the development of the research on the plastic deformation damage theory, the research on the deformation behavior of materials gradually changes from macro to micro.In the establishment of the constitutive model, in addition to considering the macroscopic deformation characteristics, the microscopic parameters are gradually introduced to characterize the microscopic deformation mechanism of materials.
The most commonly used constitutive models can be divided into three categories: phenomenological model, microscopic model, and artificial neural network model [9,10].Because of few parameters and it being easy to be established, the phenomenological model is the most important model in the engineering application, such as the Arrhenius constitutive equation, the Johnson-Cook (J-C) model, Fields-Backofen (FB) Model, Hansel-Spittel Model, and Molinari-Ravichandran Model.However, the phenomenological model cannot directly reflect the influence of material microscopic parameters on macroscopic mechanical properties.The microscopic model is based on the microscopic deformation mechanism of materials, such as the Zerilli-Armstrong Model (ZA Model), Preston-Tonks-Wallace (PTW) Model, Rusinek-Klepaczko (RK) Model, and the Cellular Automaton (CA) Model.Although the microscopic model directly reflects the influence of the microscopic parameters on the macroscopic performance, the microscopic model has more parameters compared with the phenomenological model, and it is more difficult to determine the model parameters.The artificial neural network model can be used to predict the flow stress of materials under complex deformation conditions, which avoids the difficult problem of material parameter determination of the traditional constitutive model, and is suitable for finite element simulation of the material deformation process.

Phenomenological Model
Based on the experimental data, the phenomenological model can determine the influence of stress, strain, temperature, and strain rate on the material deformation behavior through data fitting and other technical means, so as to construct a constitutive relationship model that can reflect the macroscopic deformation behavior of materials.Although the parameters of the phenomenological model are relatively few, it is easy to obtain by experimental means, and the phenomenological model does not consider the influence of material composition, grain size, dislocation density, and other factors on material properties.Therefore, it cannot reveal the micro-deformation mechanism of materials, which limits the application of the model.Usually, the prediction accuracy of the phenomenological model is high in the limited strain rate and temperature range, and the model cannot accurately reflect the deformation behavior of the material after exceeding the range.

Johnson-Cook (J-C) Model
The Johnson-Cook model (J-C) is a typical representative of the phenomenological model, which is widely used in industrial production and material property prediction research [11].As shown in Equation ( 1), the J-C model assumes that the material is isotropic and decouples the work hardening effect, strain rate effect, and temperature effect to simplify the constitutive equation, and the parameters of J-C model can be calibrated by uniaxial tensile or compressive tests of materials: where σ e is the Von-Mises equivalent stress, A is the yield stress (MPa) under the reference strain rate, B is the hardening coefficient (MPa), C is the strain rate sensitivity coefficient, n is the hardening exponent, m is the temperature sensitivity coefficient, ε p is the equivalent plastic strain, .
ε is the strain rate (s −1 ), ε p is the reference strain rate (s −1 ), T is the actual temperature, T 0 is the ambient temperature, and T m is the melting point temperature of the material.
When the deformation temperature T of the material is basically the same as the ambient temperature T 0 , and the material is loaded under quasi-static conditions, the flow stress of the material only needs to consider the effect of strain hardening on the stress.Thus, the J-C model can be simplified as The relationship between the equivalent stress σ e and the equivalent plastic strain ε p in Equation ( 2) can be adapted to all stress states, and the three parameters A, B, and n are coupled with each other.
The material parameters in the J-C model can be determined by fitting with less experimental data; however, this will not take into account the influence of complex stress states on material properties.When the traditional model parameter calibration method is used to determine the model parameters, there will be a large error between the numerical simulation results and the actual deformation of the material under complex stress conditions.In order to reduce this difference, Ru [12] took 6005A-T6 aluminum alloy as the research object and obtained the optimal solution of J-C model parameters based on the test data under different stress states with the help of the inversion calibration method based on a genetic algorithm.
The J-C model assumes that thermal softening, strain rate hardening, and strain hardening are three independent phenomena, which can be isolated from each other without considering the effect of the interaction among the three.That is to say, the original J-C model does not include any historical effect of thermal or strain rate.Although the J-C model can easily realize the calibration of parameters, it can only describe the relationship at the macro level such as strain rate, temperature, flow stress, and strain.It lacks the essence of material deformation and cannot describe the flow stress of material from the micro deformation mechanism of material.It can only describe the relationship of strain rate, temperature, flow stress, and strain at the macroscopic level, but cannot represent the microscopic deformation mechanism of the material.
Zhang [13] modified the original J-C model by considering the effect of forming temperature on strain hardening behavior of IC10 alloy.The modified J-C model can be expressed as where B(T * ) is a temperature effect function, which can be defined as where m 1 , P 1 , P 2 , and P 3 are material constants.σ br is the fracture stress at room temperature.ε br is the fracture strain at room temperature.
The flow stress of IC10 alloy was fitted by the modified J-C model, and the fitting results are shown in Figure 1.The results show that the flow stress predicted by the modified J-C model is in good agreement with the experimental results, but the prediction accuracy is poor when the temperature is 800 • C, and the strain rate is 0.0001 s −1 .
Lin and Chen [14] studied the deformation behavior of 42CrMo alloy steel under hot compression based on Johnson-Cook and Zerili-Armstrong (JC-ZA) coupled model, and established the relationship between flow stress, strain rate, and forming temperature: where A, B, n, C 3 , and C 4 are material parameters.σ is equivalent stress.ε is equivalent strain.T is the absolute temperature. .
ε 0 is dimensionless strain rate.Since about 3-5% plastic deformation energy is retained in the material as storage energy, this will increase the temperature of the sample; thus, Equation ( 6) can be modified to The flow stress of 42CrMo steel obtained by JC-ZA is shown in Figure 2. Since the JC-ZA model not only considers the strain hardening phenomena, but also considers the coupling effect of temperature and strain rate on the flow behavior, the accuracy of the prediction results of the model is improved.
( )  ( ) ̇< ̇ Under the condition of high strain rate, since the thermal effect of plastic deformation is significantly correlated with the strain rate, the original J-C model cannot reflect the strain rate effect under the high strain rate.Therefore, the prediction accuracy of J-C model will be reduced.To solve this problem, Vural and Caro [15] modified the strain rate sensitivity parameter C in J-C model to improve the prediction accuracy of the model.The modified J-C model of Vural and Caro makes temperature directly couple with strain hardening.
where T is room temperature.T r is the reference temperature.Based on the tensile tests and compressive tests, Vural and Caro [15] determined the parameters of the modified J-C model of 2139 Aluminum Alloys.By enhancing the coupling effect between temperature and strain hardening, the flow stress changes related to rate and temperature were better captured, as shown in Figure 3. Comparing the tensile test results with the compression test results, the strain rate sensitivity under dynamic state is enhanced through smooth and continuous transition.
1 ()  2 () σ  In 1992, Khan and Huang [16] proposed a viscoplastic constitutive model (KH model) to simulate the behavior of coarse-grained Al1100 alloy in a wide strain rate range: where J 2 is the second invariant of the stress deviator, D p 2 is the second invariant of plastic deformation rate, and ε 2 is the equivalent strain.
In a one-dimensional case, Equation (11) can be rewritten as where .
] n 1 (14) where g 1 (ε) represents the correlation between stress and plastic strain at the reference strain rate.g 2 .
ε is a function of the strain rate.σ represents von-Mises equivalent stress, ε is the cumulative plastic strain, and .ε is the equivalent plastic strain rate.D p 0 is the upper limit of arbitrary strain rate.n, E ∞ , σ 0 , and α are material constants.
Yu et al. [17] used the KH model to describe the flow behavior of DP600 steel, as shown in Figure 4.The results show that the KH model can almost perfectly describe the correlation between the stress and plastic strain under the strain rate 10 −4 s −1 .However, the KH model neglects the coupling effect of strain and strain rate on work hardening, and does not consider the effect of temperature, so there is a large error between the predicted results of KH model and the measured results of mechanical properties of DP600 steel at high strain rate (600 s −1 ).
In view of the above shortcomings, Equation ( 12) is modified to better describe the effect of strain rate on flow stress: where σ is flow stress, ε p is plastic strain,    With the help of the modified model (Equations ( 15) and ( 16)), the rheological properties of DP600 steel were predicted at strain rates of 10 −3 s −1 , 10 −2 s −1 , 500 s −1 , 1100 s −1 , and 1600 s −1 , respectively, and the predict results are shown in Figure 6.Compared with the experimental results, the modified KH model can better describe the deformation behavior of DP600 steel at a high strain rate.Khan and Liang [18] modified the relationship between J 2 and D P 2 in the original KH model by considering the influence of temperature effect on the mechanical properties of materials, and established a new KHL constitutive model, as shown in Equation (18): where T is the absolute temperature.ε is equivalent plastic strain.f 1 ε 2 , D p 2 is a function characterizing the coupling effect of strain and strain rate on work hardening.T r is the reference temperature.T m is the melting temperature of the material, A, B, n, n 0 , and m are fitting parameters based on test data.
In 2000, Khan et al. [19] introduced the Hall-Petch relation into the KHL model, and replaced constant A (yield stress of material under reference strain rate) with the Hall-Petch relation.The modified KHL model comprehensively characterized the effects of grain size, work hardening, rate, and temperature on the plastic deformation behavior of materials.The modified KHL model was as follows: where σ y is yield stress, d is the average grain size of polycrystalline, and a and k are material constants.Khan et al. [20] used the modified KHL constitutive model (Equation ( 20)) compared with J-C model, as shown in Figure 7.The results show that the prediction accuracy of the modified KHL model for the rheologiacal characteristics of Ti-6AL-4V alloy is significantly higher than that of the J-C model.
Figure 7.The modified KHL model compared with J-C model [19].
However, when the deformation temperature of the constitutive Equation ( 20) is lower than the reference temperature, the temperature term (1 − T * m ) will return to a negative value, which will make the model unable to predict the viscoplastic response of the material below room temperature.Therefore, in order to characterize the mechanical properties of metal materials at low temperature, Khan [20] revised Equation (20), and the modified flow stress model is shown as follows: where T is the absolute temperature, T m is the melting temperature, T r is the reference temperature,a and k are material constants, and can be determined from experiments at a strain rate of .ε = 1.d* is defined as the reference grain size when the yield stress begins to increase significantly.a 1 , a 2 , k 1 , and k 2 are constants based on the yield stress and grain size (d −0.5 ).
As the grain size of the material decreases to the nanometer range, the strain hardening behavior also increases.The modified constitutive Equation ( 22) is very close to the experimental results of nanocrystalline aluminum, as shown in Figure 8.Based on the KHL model, Farrokh and Khan [21] established a new viscoplastic constitutive equation Khan-Liang-Farrokh (KLF) model (Equation ( 24)) related to grain size and temperature: where ε p is plastic strain.
. ε p is the current strain rate.σ is the flow stress.T m represents the melting temperature of the material.T is ambient temperature.T r is the reference temperature.D p 0 is the optional upper limit strain rate.
. With the help of the above theoretical models, the uniaxial compressive deformation behavior of mechanically milled ultrafine grain and nanocrystalline Cu and Al under quasi-static and dynamic strain rates was investigated [21], as shown in Figures 9 and 10.The KLF model can well capture the variation of yield stress and work hardening behavior with grain size (to large strain), and can respond to different grain sizes (from coarse grain to nanocrystalline materials), different strain rates (quasi-static to dynamic), and different temperatures (low temperature and high temperature).

Fields-Backofen (FB) Model
In 1957, Fields and Backofen [22] proposed the FB model for most metal materials: where k is the strength coefficient, n is the strain hardening index, and m is the strain rate sensitivity index.The FB model is widely used to describe the stress-strain relationship, and can be used to express the work hardening phenomenon in the process of plastic deformation through the strain hardening index (n) and strain rate sensitive index (m).
Based on the FB model, Cheng et al. [23] studied the uniaxial tensile deformation behavior of an AZ31 magnesium alloy sheet under strain rates of 0.1-0.0001s −1 and deformation temperatures of 423-573 K, as shown in Figure 11.By comparing the predicted results of FB model with the experimental results, the predicted results are in good agreement with the experimental results at low temperature and high strain rate.However, the rheological deformation of AZ31 magnesium alloy sheet at higher temperature and lower strain rate shows obvious recrystallization softening characteristics; the FB model cannot characterize the softening behavior of the material.In view of the inaccuracy of the softening description of the FB model, Zhang [24] and Chen [25] introduced a softening term into the FB model to characterize the softening behavior of materials in the process of plastic deformation, and modified the FB model on this basis.The modified FB model established is shown below: where K is the strength coefficient, n is the strain hardening index, m is the strain rate sensitivity index, b is material constant, and T is the absolute temperature.s represents the softening ratio of magnesium alloy with increasing of the strain.The modified magnesium model was used to predict the deformation behavior of AZ31 magnesium alloy [23], as shown in Figure 12.The results showed that the modified FB model could better characterize the softening behavior of the material during plastic deformation.
 

Hansel-Spittel Model
In 1978, Hansel and Spittel [26] combined the classical strain rate power-law model with the Hollomon and Swift model, treating the corresponding material constants as functions of strain and temperature.Based on this, a new equation of thermal deformation flow stress is proposed: where A is the stress intensity coefficient, m 1 is the temperature-dependent coefficient, m 2 is the strain strengthening coefficient, m 3 is the strengthening coefficient of strain rate, m 4 is the strain softening coefficient, m 5 is the temperature-dependent strain strengthening coefficient, m 6 is strain correlation coefficient, m 7 is the strengthening coefficient of strain rate related to temperature, and m 8 is the temperature coefficient.Due to the numerous parameters of the Hansel-Spittel model, it is very difficult for regression analysis.Therefore, Spigarelli and El Mehtedi [27] modified Equation (28) to reduce the regression coefficient and make it have better prediction ability: Under thermal deformation conditions, the strain rate sensitivity of material flow stress decreases with the increase of strain.However, the Hansel-Spittel model is not sufficient to describe the change of strain rate sensitivity under given temperature and strain conditions; Brown proposed a modified Hansel-Spittel model (Equations ( 30)-( 32)) to describe the strain rate sensitivity more accurately [28].In the Brown modified model, the stress component is assumed to be temperature-independent under a given strain (or steady state) condition; thus, the parameters m 7 and m 8 were neglected to simplify the constitutive model: . .
where σ is the flow stress, T is the absolute temperature, A, A 2 , and β are the material constant, n is the strain rate sensitive index, R is the molar constant of gas R = 8.31 J/mol•K, α is the stress level parameter, and Q is the activation energy of the material.
The results show that the modified Hansel-Spittel model can better characterize the rheological properties of CoCrFeMnNi alloy than the original Hansel-Spittel model, as shown in Figure 13, where the real lines represent simulated data, and the virtual lines represent experimental data. ---

Arrhenius Model
The Arrhenius equation is widely used to describe the relationship between flow stress and strain rate and deformation temperature under high temperature deformation conditions.In the exponential equation, the Zener-Hollomon parameter is used to represent the influence of temperature and strain rate on the deformation behavior.Moreover, a more Metals 2022, 12, 2077 13 of 33 accurate approximation between the Zener-Hollomon parameter and the flow stress is obtained by hyperbolic function.The Arrhenius constitutive equation is as follows [29]: where A 1 , A 2 , A, n 1 , n, α and β are material constants, Q is the activation energy of deformation (K J/mol), R is a general gas constant (8.31 J/(mol • K)), σ is flow stress.The relationship between α, β and n 1 can be expressed as α = β/n 1 .
According to the definition of hyperbolic function, the flow stress σ can be written as a function of Z parameter: In recent years, many scholars [30,31] have determined the flow stress of different metals and alloys during thermal deformation by the Arrhenius model.Since the Arrhenius model does not consider the influence of strain on flow stress, the prediction accuracy of this model is not high under large deformation conditions.By constructing the polynomial relationship between the model parameters Q, A, N, and α and the strain, respectively.The modified Zener-Hollomon model coupled with the influence of temperature, strain rate, and deformation was established [32], as shown in Equations ( 35) and ( 36): ln

Molinari-Ravichandran Model
In 2005, Molinari and Ravichandran proposed a phenomenological constitutive model (MR model) that can describe the flow behavior of metals under various loading conditions [33,34].The MR model is based on a single internal variable δ, which is related to the characteristic length of the microstructure formed during metal deformation.The flow stress σ is expressed as a function of the inherent resistance σ0 and strain rate .ε of the material: where .ε 0 is the reference strain rate, d is the average grain size, and m represents the parameter of instantaneous strain rate sensitivity dependent on temperature.The reduction Metals 2022, 12, 2077 14 of 33 of the internal characteristic length δ with an increase of increasing equivalent plastic strain ε, which can be expressed by the following phenomenological evolution equation: where δ r is a dimensionless parameter characterizing the refinement rate of microstructure.δ s is the saturation length at relatively large strains, which both can be a function of strain rate and temperature.Assume δ r and δ s are constants, which is the case for constant strain rate and temperature.Equation ( 41) can be expressed as follows: where δ 0 is the initial internal characteristic length.However, when the saturation of the "microstructural length" reaches the limiting value δ s , the flow stress will also be induced to reach saturation values.
The laws proposed for δ r and δ s can be based on the theory of thermally activated processes or on empirical relations depending on temperature and strain rate: where δ s and δ s0 are the saturation values of effective microstructure length δ and the reference values at zero strain rate, respectively.The material constants α s , ξ s , and v s control the strain rate and the dependence of δ s on temperature; similarly, δ r depends on δ r0 , α r , ξ r , and v r .
In order to avoid negative values of the MR model at relatively low strain rates, Molinari and Durrenberger [35] proposed a new phenomenological model (DMR model).The flow stress is defined as the sum of the internal stress representing the long-term interaction and the thermally activated effective stress: In the isothermal condition, Durrenberger [35] pointed out that Equation ( 45) can be simplified as where σ ref is a reference stress, T 0 is the reference temperature, .ε 0 is the reference strain-rate, ν* is the temperature sensitivity, a and b are material constants, and m characterizes the instantaneous rate sensitivity.
The DMR model is considered to be able to accurately describe the thermal and viscoplastic response behavior of metals under harsh loading conditions [36].Based on the DMR model, the strain rate sensitivity simulation and experimental data of four kinds of sheets for automobile forming are compared as shown in Figure 14 (Circle points represent experimental data, real lines represent simulation data) [36].It can be seen that the DMR model can effectively describe the strain hardening effect of materials under different strain plastic rates.

Comparison and Analysis of Phenomenological Models
In addition to the above widely used models, there are many phenomenological constitutive models.Due to the limitation of the length of the paper, other phenomenological constitutive models are not introduced in detail in this paper.
In order to facilitate readers to understand the features and applicable conditions of the phenomenological constitutive model introduced in the article, the main features of each model are summarized in Table 1.
Table 1.The characteristics of the phenomenological constitutive model.

Model Name Characteristics and Application Conditions
Johnson-Cook (J-C) Model

1.
The model has few parameters and is easy to solve 2.
The influence of deformation temperature and strain rate is taken into account 3.
Easy model correction, modified J-C model Khan-Huang (KH) Model

1.
It consists of elastic component and plastic component 2.
Temperature effects are not taken into account Khan-Huang-Liang (KHL) Model 1.
The modified model of the KH model 2.
Consider the temperature effect of deformation
The modified model of the KHL model 2.
Both temperature effect and strain rate effect are taken into account 3.
The grain size effect of the material is considered Fields-Backofen (FB) Model 1.
Both temperature effect and strain rate effect are taken into account 2.
The model has few parameters 3.
The parameter fitting process is simple Combined the classical strain rate power-law model with the Hollomon and Swift model 2.
The effects of both temperature and strain rate are considered 3.
The model can be solved by regression analysis Arrhenius Model 1.
The influence of strain rate and temperature at high temperature is considered 2.
It can predict the peak stress of the material 3.
The modification of model parameters can be suitable for large deformation Molinari-Ravichandran Model 1.
The effects of temperature, strain rate, and grain size are considered comprehensively 2.
Based on the microscopic characteristics of the material 3.
The DMR model can characterize the thermal and viscoplastic behavior

Constitutive Model Based on Microscopic Characteristics
The phenomenological model can obtain the corresponding model parameters through the fitting and the regression analysis of experimental data, which makes the process of solving the model easy.However, under the extreme deformation conditions, such as large deformation and high temperature, the phenomenological model cannot characterize the effect of the change of microstructure on the macroscopic mechanical properties of materials [37].
In the thermal deformation process, the plastic deformation behavior of the material can be divided into four stages, as shown in Figure 15.In the first deformation stage (stage 1), the accumulation and multiplication of dislocations are the main reasons for material hardening.However, in the second stage of deformation (stage 2), the temperature effect and strain hardening simultaneously affect the flow stress of the material.In this stage, the strain hardening effect plays a leading role, and the strain hardening rate of the material gradually decreases due to the effect of temperature effect.In the third stage of deformation (stage 3), with the increase of deformation amount, dislocations accumulate continuously and grain boundary energy storage increases, which enhances dynamic recrystallization (DRV) under temperature and reduces strain hardening effect.In the fourth stage (stage 4), strain hardening, dynamic recovery, and dynamic re-crystallization are relatively balanced, and the flow stress of the material is in the steady state deformation stage.Based on the above micro deformation theory, He [38] believed that the total strain rate .ε in isothermal deformation was composed of elastic strain rate .ε e and inelastic part .ε in , and proposed a multi-scale constitutive model based on dislocation density to describe the thermal deformation behavior of HEAs materials.As shown in Equations ( 47) and ( 48):
where ρ cm is the dislocation density, b is magnitude of the Burgers vector, λ is the average distance, v a denotes the attempt frequency, m is the Taylor orientation factor, ∆G is the activation energy, k is the Boltzmann constant, T is the absolute temperature, σ eq is the equivalent stress, and S is the deviatoric part of local stress tensor.

Zerilli-Armstrong Model (ZA Model)
The ZA model takes into account the effects of strain hardening, strain rate hardening, and temperature softening on the plastic deformation behavior of metals and alloys, and is one of the most widely used dislocation density models.The ZA model divides the flow stress into thermal component and athermal component, as shown in Equation ( 49) [39]: The athermal component σ a considering the effect of grain size is expressed as where l is the average grain size, k is the Hall-Petch constant, C 0 is the fitting parameter.The thermal stress component σ th can be calculated by Equations ( 51) and ( 52): where M is the direction factor, ∆G 0 is the free energy of thermal activation at an absolute temperature of zero, A is the thermal activation area at at an absolute temperature of zero, b is the Burgers vector, β is a parameter related to strain rate.For BCC materials, A is a constant, that is, the deformation activation energy of BCC material is not affected by plastic strain.For FCC materials, the relationship between A and plastic strain can be expressed as ε 1/2 .Therefore, according to the Equations ( 50)-( 52), the constitutive equation of Z-A model can be expressed as ε , for FCC materials (54) where C 0 , C 1 , C 2 , C 3 , C 4 , C 5 , and n are fitting parameters.It can be seen from Equations ( 53) and ( 54) that the strain hardening index is affected by strain rate and temperature for FCC materials, while for BCC materials, the strain hardening index is constant.The ZA constitutive model does not fully realize the coupling effect of temperature, strain, and strain rate.Therefore, some improved ZA models have been proposed in recent years.Zhang et al. [40] proposed an improved flow stress model of IC10 alloy considering the effects of temperature, strain rate, and strain on the flow behavior of IC10 alloy under a wide range of temperature and strain rate deformation conditions, as shown in Equations ( 55) and (56).Moreover, the flow behavior of IC10 under different deformation conditions is predicted by the modified model, and the prediction results are shown in Figure 16.The results show that the modified ZA model can predict the rheological behavior of metal materials more accurately: f (T) , for FCC materials ( 56) where ρ 0 is the initial dislocation density of materials.M is the material constant related to the increase rate of dislocation density.C 0 , C 1 , C 2 , C" 3 , C' 4 , and n are fitting parameters.
̇   In order to better characterize the flow stress of D9 alloy under the high temperature, based on the J-C model, Samantaray [41] modified the ZA model (Equation ( 54)) by considering the coupling effect of temperature, strain, strain rate, and temperature on flow stress.The modified model proposed by Samantaray is shown in Equation ( 58), and the predicted results of D9 alloy are shown in Figure 17:

Preston-Tonks-Wallace (PTW) Model
In 2003, Preston et al. [42] proposed a numerical simulation model of metal plastic flow suitable for explosive load and high-speed impact.At low strain rate, the dependence of plastic strain rate on applied stress is in the form of Arrhenius, but the activation energy is singular at zero stress, so the deformation rate disappears at this limit.The work hardening effect is modeled as a generalized Voce law.The constitutive model of the low strain rate region can be expressed as where ψ is the plastic strain, p is a dimensionless material parameter, and θ is determined by the slope of work hardening data.τ is the normalized flow stress.τs is the work hardening saturated stress under thermal activation condition, τy is the yield stress, where .
Ψ is the plastic strain rate.s 0 and s ∞ are the values of τs when temperature is zero and infinite, respectively.y 0 and y ∞ are the values of τy when temperature is zero and infinite, respectively.k and γ are dimensionless material constants.T = T/T m (T is Metals 2022, 12, 2077 20 of 33 deformation temperature, T m is melting temperature).
. ξ is the equivalent scaling factor, and it can be defined as where ρ is material density, M is atomic mass, G is the shear modulus of material, which is the function of density and temperature, G = G 0 1 − α T , where G 0 is the shear modulus at absolute temperature, and α is the temperature dependence coefficient.
Under ultra-high strain rate loading under strong shock wave, the relationship between plastic strain rate and critical shear stress satisfies an exponential function, and the yield stress is equal to the saturation stress.Therefore, the model can be expressed as follows: (63) where A and β are the material constants.
In order to maintain the continuity of the model, the saturation stress and yield stress in the transition zone from the low-strain rate region to the high-strain rate region are expressed as In order to explain the strain hardening at 0-10 6 s −1 strain rate, Kim [43] used the Voce equation to modify the strain hardening term in the original PTW model, and the modified constitutive equation is shown in Equations ( 66) and (67).Compared with the other constitutive models, the Kim's modified PTW model can better describe the deformation behavior of tantalum in a wider strain, strain rate, and temperature range, as shown in Figure 18. ξ̇= where τ is the initial yield stress, and it can be expressed as where A, B, and C are material constants.A represents the maximum strain hardening at unit strain rate, B represents the saturation velocity of control hardening, and C represents the dependence of maximum hardening on strain rate.

Rusinek-Klepaczko (RK) Model
Considering the influence of strain hardening, strain rate sensitivity, temperature effect and other factors on the material deformation process, Rusinek and Klepaczko [44] proposed a viscoplastic constitutive model (RK model), which decomposed the flow stress σ into the sum of internal stress σ u and effective stress σ*: where E 0 , T m , and θ* represent the elastic modulus, melting point temperature, and characteristic assimilation temperature of the material at absolute zero, respectively.The internal stress σ u can be expressed as where B 0 and D 2 are material parameters.v is the temperature sensitive coefficient.n 0 is the strain hardening index at absolute zero.
ε max are the minimum strain and maximum strain rate accepted by material, respectively.
The effective stress can be obtained according to the Arrhenius equation, as shown in Equation ( 73): In the effective stress equation, σ* is the effective stress at absolute zero.D 1 is the material parameters.m* is a constant allowing the definition of strain rate and temperature dependence.
Aiming at the problem that the original RK model cannot accurately describe the viscoplastic behavior of aluminum alloy, Rusinek et al. [45] modified the RK model.They defined the macroscopic negative strain rate sensitivity (NSRS) and viscous resistance, and added a negative strain rate sensitivity term on the basis of the original RK model, as shown in Equation ( 74): where σ ns .ε p , T is the stress component of NSRS, which depends on the strain rate and temperature: In the above formula, σ ns 0 and D 3 are material parameters, which describe the stress drop caused by the interaction of dynamic strain aging (DSA) with strain rate and tempera- .ε trans is the turning point from positive to negative strain rate sensitivity, which can be obtained by fitting experimental results.
Based on the Equations ( 74) and ( 75), the viscoplastic behavior of AA5083 aluminum alloy was researched, as shown in Figures 19 and 20.It can be seen that the analytical prediction results are consistent with the experimental results under different strain rates, and the results also prove that the effectiveness of the model is describing correctly the DSA effect [45]. ( ) Another extension of the RK model is by adding a viscosity component related to the strain rate on the basis of the the equivalent Huber-Misses stress σ of the original model.In this way, the viscous resistance effect of face-centered cubic metal at high deformation rate is characterized.As shown in Equation ( 76), where σ ath .ε p is the stress component accounting for the viscous resistance which is just dependent on strain rate.σ ath .ε p can be specifically expressed as where M is the Taylor factor, B is the resistance coefficient, ρ m is the mobile dislocation density, and b is the magnitude of the Burgers vector, χ is a material constant, α represents an effective damping coefficient affecting the dislocation motion, and τ y denotes the effective damping coefficient that affects the dislocation motion.
In 2010, Rusinek et al. [10] proposed a modified RK model based on physical concepts.Compared with the original RK model, the internal stress σ u is defined as a constant Y in the new model, and the effect of plastic strain is considered in the effective stress component σ*: where E(T) is an expression for the temperature-dependent Young's modulus, as shown in Equation ( 80).E 0 , T m , and θ * are Young's modulus, melting point temperature, and characteristic assimilation temperature of the material at absolute zero, respectively.σ µ the internal stress of the material, in Equation ( 79), and the internal stress σ µ is equal to the flow stress Y of the materials without deformation, namely σ µ = Y.σ * is effective stress of the material, and the expression is shown in Equation (81).σ vs is the viscous stress component, and its expression as shown in Equation ( 82): where T m and θ * are melting point temperature and characteristic assimilation temperature of the material at absolute zero, respectively: where ξ 1 and ξ 2 are material constants describing temperature and rate sensitivities of the material, respectively.
. ε max the maximum strain-rate level for material:

Cellular Automaton (CA) Model
In the dynamic re-crystallization (DRX) process, the accuracy of nucleation parameters is the key factor to simulate the micro-structure evolution.Goetz et al. [46] proposed the cellular automaton model (CA model) to describe the microstructure evolution during dynamic recrystallization.With plastic deformation, dislocation density evolution can be expressed by the following equation: where k 1 and k 2 are model parameters representing strain hardening and dynamic recovery, respectively.When the dislocation density exceeds the critical value ρ c of DRX nucleation, a new atomic nucleus is formed on the grain boundary.For two-dimensional cellular automaton model, I(ρ>ρ c ) is defined as the number of nuclei formed per unit time.For boundary cells with dislocation density exceeding the critical value, the nucleation probability P N in the time step ∆t is where N CA is the number of cells along the unit boundary length.L CA is cell length.
Transformation probability P G can be expressed as In order to ensure that P N and P G are less than 1, an upper limit of ∆t must be set.Therefore, the upper limit of ∆t can be defined as the ratio of cell size L CA to the maximum growth rate v max , namely, The study of dynamic re-crystallization behavior by CA method is 1 based on the physical metallurgical principles.However, it is difficult to accurately estimate the nucleation rate due to the complexity of DRX and limitation of experiments.The method proposed in the literature [47-51] can be used to identify the value of model parameters and explore their correlation with thermomechanical process parameters, which will contribute to the development of thermomechanical process analysis model.With the help of the new model parameter prediction method, Jin et al. [47] studied the microstructure evolution and flow stress behavior of oxygen-free high-conductivity (OFHC) copper during DRX at the temperature of 775 K and the strain rate of 0.002 s −1 ; the results are shown in Figure 21.Although the influence of strain rate has been taken into account for most constitutive models, the constitutive models mentioned above cannot accurately reflect the deformation behavior of materials under high strain rate conditions, such as impact and explosion.In order to accurately describe the deformation characteristics of materials under high strain rate and high-pressure coupling.Steinberg et al. [52] proposed the SG constitutive model in 1980.In the SG model, the shear modulus G and yield stress Y of the material are expressed as functions of pressure P and temperature T: The SG model needs to satisfy Equations ( 90) and (91): where G 0 and Y 0 are shear modulus and yield stress at temperature T = 300 K, presser P = 0 and strain ε = 0, respectively.η is the volume compression ratio of the material during deformation, η=V 0 /V.G' P is the first partial derivative of shear modulus G with respect to pressure P, G' T is the first partial derivative of shear modulus G with respect to temperature T, Y' P is the first partial derivative of yield stress Y with respect to pressure P. β and n are strain hardening parameters.Y max is the maximum allowable value of material flow stress Y during strain hardening.ε i is the initial strain, usually ε i = 0. ε i is the plastic strain of the material.Under the condition of high-speed impact, Steinberg believed that the softening effect caused by the increase of temperature counteracts the hardening effect caused by the strain rate.When the strain rate is greater than 10 5 s −1 , the flow stress of the material reaches the maximum value.Therefore, the expression of SG model has nothing to do with strain rate and cannot characterize the mechanical properties of materials at low strain rate.To overcome this shortcoming, Steinberg and Lund [53] is the shear modulus associated with temperature T and pressure P: .
where Y P is the Peierls-Nabarro resistance to plastic deformation of metal.ρ is the dislocation density.b is the Burgers vector.w is the width between the dislocation rings.v is the Debye frequency.L is the length of the dislocation.D is the resistance coefficient.K is Boltzmann's constant.U K is the energy required to form dislocation rings on a dislocation segment of length L.
Compared with the SG model, the SL model combines the microscopic deformation mechanism with the macroscopic mechanical properties, and it can be used to fit the plastic deformation properties of metal materials with strain rates of 10 −4 and 10 6 s −1 .
However, the SG model or SL model can comprehensively reflect the dependence of the shear modulus and yield strength of the material on pressure and temperature.However, the physical mechanism of material deformation under strong impact conditions is very complex, so the applicability of the model has always been the focus of scholars' research and exploration [54,55].In recent years, in order to better describe the dependence of shear modulus and yield strength on temperature and pressure, the revision of SG model based on physical deformation characteristics of materials has become one of the focus of scholars' attention and research [56,57].

Comparison and Analysis of Microscopic Constitutive Models
Based on dislocation dynamics and thermal activation theory, the microscopic constitutive model combines strain hardening, temperature, and strain rate effects with dislocation motion and plastic deformation.However, the deformation mechanism of materials under different deformation conditions is very different.Moreover, the applicable conditions and characteristics of the microscopic constitutive model are different.In order to facilitate the application of the model, the characteristics and application conditions of common microscopic constitutive models described in this paper are shown in Table 2.
Table 2.The characteristics of microscopic constitutive model.

Model Name Characteristics and Application Conditions
Zerilli-Armstrong Molde (ZA Model) 1.
Based on the dislocation dynamics theory 2.
The influence of deformation temperature, strain rate, and the average size of grains are taken into account 3.
The influence of crystal structure is distinguished 4.
It can be used to characterize the plastic deformation behavior of the thermal activated region Preston-Tonks-Wallace (PTW) Model

1.
The effect of strain rate mechanism is considered, dislocation slip region, transition region, and ultra-high strain rate region 2.
Applicable to a wide range of strain rates 3.
It can characterize the nonlinear dislocation effect under impact load Rusinek-Klepaczko Model (RK Model) 1.
The flow stress is decomposed into internal stress and effective stress 2.
The influence of temperature on Young's modulus is considered 3.
The temperature effect and strain hardening are considered simultaneously 4.
The modified RK model can describe the viscoplastic behavior of FCC material Cellular Automaton Model (CA Model) 1.
Based on physical metallurgical principles 2.
It can describe the microstructure evolution during dynamic recrystallization 3.
The effects of temperature and strain rate on dynamic recrystallization are considered Steinberg-Guinan Model (SG Model) 1.
The effects of temperature and pressure are considered 2.
The effect of strain rate is not considered 3.
It is suitable for high strain rate (>10 5 s −1 ) and high-pressure conditions Steinberg-Lund Model (SL Model) 1.
The effects of temperature and pressure are considered 2.
The flow stress is the sum of the thermal and non-thermal stress components 3.
It can be used over a relatively wide range of strain rates (10 −4 -10 6 s −1 )

Artificial Neural Network Model (ANN Model)
In general, the thermal deformation behavior of materials is usually described by the phenomenological model or dislocation density model mentioned above.However, the deformation behavior of materials at high temperature and strain rate is highly nonlinear, and many factors affecting flow stress are also nonlinear, which results in low accuracy of flow stress predicted by the regression method and limited scope of application.
Artificial neural network (ANN) is the most suitable method for solving the problem which is difficult to solve by the traditional calculation method.The main advantage of this method is that there is no need to assume a mathematical model.At present, the artificial neural network model is widely used in industrial production and academic research [58,59].

Back-Propagation (BP) Neural Network
The BP (back-propagation) algorithm is one of the most commonly used feedforward neural network learning algorithms in materials modeling.The BP model has strong capacity for learning, generalization, fault tolerance, and nonlinear mapping.The typical structure of BP model is schematically illustrated in Figure 22, which consists of one input layer, one output layer, and one or more hidden layers [59][60][61][62][63].
In the BP model, the number of nodes in the input layer depends on the number of parameters affecting the actual deformation of the material, and the output layer parameter is the true stress of the material.The selection of the number of hidden layer neurons is a complex and important problem, which will have a great impact on the performance of the BP neural network.If the number of hidden layer neurons is too small, the BP network cannot converge during model training.On the contrary, if there are too many neurons in the hidden layer, the training data will be over-fitted, which will prolong the learning time of the network and reduce the generalization ability of the network.There is no rigid rule to find the number of neurons in the hidden layer, but it can be determined according to empirical Equation (96) [59] or Equation (97) [60]: where M is the number of hidden layer neurons, n is the number of neurons in the input layer, m is the number of neurons in the output layer, a is a positive integer ranging from 1 to 10, and k is the number of training set.In addition to the number of hidden layer neurons, other structural parameters of BP neural network, such as the number of hidden layers, transfer function, and training function, will also have a great impact on the convergence speed and accuracy of the neural network model [61,62].The Sigmoid function is usually used as the excitation function of hidden layer neurons, and the Purelin function as the excitation function of output layer neurons.The output O pj of the p training sample unit j from the input layer to the hidden layer can be expressed as where N j is the input of jth neuron in the hidden layer, w ij is the connection weight between the input layer and the hidden layer, X j is the output of jth neuron in the input layer to the hidden layer, and θ j is the threshold of the input layer to the hidden layer.Through modifying the connection weight among different layers and the threshold among different neurons, the high precision modeling with complex nonlinear characteristics can be completed.However, in order to realize the prediction of material flow stress, it is necessary to train the BP model with the help of experimental data.
Since the experimental data were measured in different units and value ranges, such as temperature, strain, strain rate, and stress, if the experimental data are directly used for model training, the convergence speed and accuracy of the model will decrease, or even not converge [63,64].Therefore, the input and output of training samples and test samples based on experimental data must be normalized into dimensionless units before training.Through normalization processing, all input and output data are converted into numbers between 0 and 1, so as to eliminate the difference of orders of magnitude between experimental data of different dimensions.
In addition, the choice of the training function and parameters of the neural network model also has a great influence on the convergence speed and prediction accuracy.When the tansig and purelin functions are selected as network inter-layer transfer functions, the trainlm function of the Levenberg-Marquardt algorithm is used as a training function to optimize the BP algorithm; this will effectively reduce the number of iterations of the model and improve the convergence speed and accuracy [61][62][63].

Applications of (BP) Neural Network
At present, the ANN model has been widely used to predict the microstructure evolution and mechanical properties of metal materials in plastic deformation.Sheikh et al. [65] used the artificial neural network model to predict the flow stress of aluminum alloy (AA5083), trained the neural network with the feedforward and back propagation algorithm, and calculated the stress-strain curves in different regions, which were compared with the experimental data to verify the reliability of the prediction, as shown in Figure 23.Yan et al. [66] established the Arrhenius constitutive equation and artificial Neural network (ANN) model to study the flow behavior of Al-6.2Zn-0.70Mg-0.30Mn-0.17Zraluminum alloy; the experimental and computational results are shown in Figure 24.The results show that the well-trained neural network model can describe the real stresstrue strain curve well.Sabokpa et al. [67] have successfully predicted the mechanical properties of AZ81 magnesium alloy under different temperature conditions by using the BP algorithm, and pointed out that the BP network model can reflect the complex nonlinear relationship between the ultimate properties of alloys and the deformation temperature.In addition, Li et al. [68] and Haghdadi et al. [69] also compared the prediction ability of the traditional constitutive model and the BP model for different metal materials, and their research results showed that the prediction ability of the BP model was more accurate than that of the traditional constitutive model.
However, the BP neural network has some inherent flaws, such as the need to learn a large number of sample data to ensure the prediction accuracy, the slow learning speed, and tendency of falling into local minimum value and network instability [70][71][72][73].Therefore, it is necessary to optimize the BP neural network to obtain a more accurate and reliable network model and better prediction performance.Ding et al. [74] used the genetic algorithm (GA) and particle swarm optimization (PSO) algorithm to optimize BP neural ical properties under complex deformation conditions.These characteristics determine that the artificial neural network model is suitable for finite element simulation of the material deformation process.In modern scientific research, the finite element numerical simulation plays an increasingly important role in many fields [76,77].However, the accuracy of numerical simulation results is closely related to geometric boundaries such as material models.
In the numerical simulation model of aluminum alloy deformation behavior, Stoffel [77] used the ANN constitutive model instead of the traditional physical constitutive model to accurately predict the dynamic response of the aluminum alloy sheet sample under impact load.Bobbili [78] applied the ANN constitutive model to ABAQUS numerical simulation software and successfully predicted the dynamic impact deformation behavior of Ti-13NB-13Zr titanium alloy.
However, the prediction accuracy of the ANN model is greatly affected by the training database, so it is necessary to have a large enough test results database to train the model, which will lead to the increase of the test cost.Therefore, based on the limited test data, a numerical simulation method is developed which can quickly build a database for neural network model training.This will be an effective way to improve the prediction accuracy and application of a neural network model.For example, Gao [75] used the crystal plastic finite element model to build a training database, which is an effective means.
The ANN model has obvious advantages in finite element simulation.However, the ANN model cannot reflect the action mechanism of various internal influencing factors during plastic deformation of materials, and also cannot replace the physical constitutive model.Therefore, the physical constitutive model and ANN model will have their respective advantages in different application fields.

Conclusions
In this paper, the constitutive models of metal material are reviewed.According to their different mechanisms and characteristics, they can be divided into three types, namely the phenomenological model, microscopic model, and artificial neural network model.For each model, typical examples are given and their applicable conditions, advantages, and disadvantages are analyzed: (1) The phenomenological model can be used to determine the constitutive relationship model of materials by means of regression analysis of experimental data.The solution process is relatively simple, but it cannot reflect the microscopic deformation mechanism of materials.Therefore, it is necessary to continuously modify to couple the effects of strain, strain rate, and temperature on the deformation behavior of materials.
(2) The microscopic constitutive model based on thermodynamics and dynamics theory can accurately characterize the microscopic deformation mechanism of materials under complex deformation conditions, and it also can accurately predict the flow stress of materials.However, compared with the macroscopic constitutive model, the microscopic constitutive model has more material parameters, which requires more complex experiments and processing methods.Therefore, the engineering application of the model is limited to a certain extent.
(3) As a new method, the artificial neural network (ANN) models are increasingly used to predict the flow stress of materials under complex deformation conditions, and it is suitable for finite element simulation of the material deformation process.However, the application of the ANN models largely depends on the training database, and cannot reflect the action mechanism of various internal influencing factors during plastic deformation of materials.
Author Contributions: X.J., methodology, investigation, review and editing, and supervision.K.H., original draft preparation and visualization.Z.L., visualization and review.Z.F., visualization and review.All authors have read and agreed to the published version of the manuscript.

Figure 3 .
Figure 3. Vural and Caro modified J-C model [15], (a) temperature dependence of flow stress at 10 −4 s −1 reference strain rate; (b) effects of strain rate and temperature under dynamic loading.

Figure 4 .Figure 5 .
Figure 4. Comparison of KH model predictions and experimental data at different strain rates [17].

Figure 6 .
Figure 6.Comparison of the new constitutive model prediction data and experimental data [17].

Figure 8 .
Figure 8. Prediction results of mechanical properties of nanocrystalline aluminum based on the KHL model [20].

ε p *
is the reference strain rate.n * is the material following the Hall-Petch relationship.d is the average size of grains.d 0 is the grain size of coarse grains.a, k, B, n o , n 1 , n 2 , c, and m are the material constants of the KHL model.

Figure 12 .
Figure 12.The stress prediction results of AZ31 magnesium alloy based on the modified FB model [23].

Figure 17 .
Figure 17.The predicted flow stresses of D9 alloy based on the Samantaray modified model [41].
the Kim's modified PTW model can better describe the defor

3. 5 .
Steinberg-Guinan (SG) Model and Steinberg-Lund (SL) Model introduced the calculation method of strain rate in the HM model into the SG model and established the Steinberg-Lund model (SL model), which decomposed the yield stress Y into thermal stress component Y T and non-thermal stress component Y A : Y = Y T .ε p , T + Y A f ε p G(P, T) G 0 (92) where Y T .ε p , T is the thermally activated part of the yield strength and is a function of .ε p and T. Y A f ε p is the non-thermal part of the flow stress, f ε

Figure 22 .
Figure 22.Schematic diagram of the BP neural network model.

Figure 23 .
Figure 23.Comparison of prediction results between experimental data and non-sampling results [65].