A CDM-based unified viscoplastic constitutive model of FLCs and FFLCs for boron steel under hot stamping conditions

Forming limit curves (FLCs), which are constructed using the limit strains at localised necking, are the most widely used tools for the evaluation of the formability of sheet metals. Fracture forming limit curves (FFLCs) are more recently developed, complementary tools for formability evaluation which are instead constructed using the limit strains at fracture. Since the formability depends strongly on forming conditions such as strain state, temperature and strain rate, models for predicting FLCs and FFLCs are essential for the optimisation and further application of hot forming processes in which these forming conditions vary significantly with both position and time. However, no model has so far been developed to predict FFLCs either alone or in conjunction with FLCs for sheet metals such as boron steel under hot stamping conditions. In this study, a set of unified viscoplastic constitutive equations for the prediction of both FLCs and FFLCs based on continuum damage mechanics (CDM) has been formulated from a set of recently developed constitutive equations for dislocation-based hardening, in combination with two novel coupled variables characterising the accumulated damage leading to localised necking and fracture. The novel variables take into account the effects of strain state, temperature and strain rate on the formability of sheet metals. The material constants in the CDM-based constitutive equations have been calibrated using experimental data comprising true stress-true strain curves and limit strains of a 22MnB5 boron steel obtained at a range of temperatures and strain rates. Investigation of the effect of varying selected parameters in the coupled damage variables on the resulting computed FLCs and FFLCs has demonstrated the flexibility of the model in enabling curves of different shapes and numerical values to be constructed. This indicates the potential of the CDM-based constitutive model for application to other materials for warm or hot stamping processes.


Introduction
Hot stamping of boron steels is commonly used for forming automotive panel components with complex shapes and ultra-high ultimate tensile strength (UTS) of up to 1500 MPa for a wide range of industrial applications (Karbasian and Tekkaya, 2010). For example, hot stamped steel components and assemblies were widespread in nearly all vehicle models which were to be launched onto the market in 2019 (Bachman, 2018). In a hot stamping process, boron steel is heated to above the austenitic temperature Ae 3 (around 900 C) to produce an austenitic microstructure and then transferred into water-cooled dies for simultaneous stamping and quenching . It is well known that the formability of sheet metals varies with strain path, temperature and strain rate (Bruni et al., 2010;Khan and Baig, 2011), all of which significantly vary with both position and time in industrial hot forming processes (Gao et al., 2017;Karbasian and Tekkaya, 2010). However, it is not possible to determine the formability of materials under arbitrary sets of conditions such as those encountered in hot stamping (Min et al., 2010;Turetta et al., 2006) using existing experimental procedures. Instead, formability is evaluated under a limited number of discrete, ideal conditions such as uniaxial, plane-strain and equi-biaxial tensions at constant temperatures and strain rates, and the resulting values are connected to construct forming limit diagrams (FLDs) (Shao et al., 2016). Therefore, developing a constitutive model which enables the formability of the materials under arbitrary forming conditions to be predicted is essential for the optimisation and further applications of practical hot stamping processes.
FLDs are the most popular tools used to represent the formability of sheet metals (Goodwin, 1968;Keeler, 1965). An FLD is a plot of the limit strains at the onset of localised necking in principal strain space. Forming limit curves (FLCs) are constructed by connecting limit strains in strain states from uniaxial to equi-biaxial obtained at a given temperature and strain rate. Due to their high sensitivity to strain paths, the limit strains at localised necking are usually obtained through loading with linear strain paths (Hu et al., 2018;Kuroda and Tvergaard, 2000). The fracture forming limit diagram (FFLD) is an extension of the FLD (Atkins, 1996), in which fracture forming limit curves (FFLCs) are constructed by connecting the limit strains at fracture in different strain states for a given temperature and strain rate. According to the experimental data ( Figure 3 in Section Determination of material constants from experimental data) on boron steel under hot stamping conditions reported by Zhang et al. (2022), there is a large difference between the limit strains at the onset of localised necking and those at fracture. From this experimental work, it has also been observed that with increasing deformation, necking progresses continuously, terminating at fracture. In some forming cases, however, it is very difficult to observe the onset of necking under hot stamping conditions. Therefore, in addition to FLCs, the use of FFLCs as complementary formability evaluation tools can improve the understanding of the extent of deformation before fracture in practical forming processes in ductile materials. The use of FFLCs to complement FLCs for the determination of formability has previously been reported in (Isik et al., 2014;Zahedi et al., 2019).
Considerable efforts have been devoted to the modelling of FLCs for metal forming applications (Bouktir et al., 2017;Brunet et al., 2004;Chan et al., 2005;Chow and Jie, 2009;Li et al., 2019;Thibaud et al., 2016). A brief review of the models for FLCs can be found in Zhang et al. (2018). Based on the theories adopted to describe failure mechanisms in the material, these models can be mainly divided into four categories: (i) models based on bifurcation theory, e.g. Hill's model (Hill, 1952), Swift's model (Swift, 1952), modified maximum force model (Hora et al., 2007(Hora et al., , 2013, St€ oren and Rice's (S-R) vertex theory (St€ oren and Rice, 1975), (ii) models based on geometrical imperfection theory, e.g. the Marciniak and Kuczy nski (M-K) model (Marciniak and Kuczy nski, 1967), modified M-K models (Hutchinson et al., 1978;Needleman and Triantafyllidis, 1978;Parmar et al., 1977), (iii) models based on continuum damage mechanics (CDM), e.g. plane-stress CDM-based models (Lin et al., 2014;Mohamed et al., 2015), and (iv) others, e.g. an empirical model proposed by the North American Deep Drawing Research Group (NADDRG) (Levy, 1996), and the Jones and Gillis (J-G) model (Jones and Gillis, 1984). Although most of the existing models for FLCs have been applied to cold forming processes, some have also been applied to high-temperature forming by taking into consideration the influence of temperature on deformation and failure in the material (Jie et al., 2011). For example, the S-R vertex theory was used to predict FLCs for boron steel (22MnB5) under hot stamping conditions by Min et al. (2010) and magnesium alloy at elevated temperatures by Min et al. (2013). The M-K model was adopted to predict FLCs at elevated temperatures for AA5083 by Zhang et al. (2009) and for AA3003 by Abedrabbo et al. (2006aAbedrabbo et al. ( , 2006b. By contrast, the plane-stress CDM-based model of Lin et al. (2014), in which a set of dislocation-based hardening constitutive equations and a damage variable was used, was developed purposely for predicting FLCs in sheet metals under hot stamping conditions and was successfully applied to the formability prediction for AA5754 in warm forming processes (Bai et al., 2016). A modification of this CDM-based model was adopted by Shao et al. (2017) for the prediction of FLCs in AA6082 under hot stamping conditions. Many approaches have been developed for modelling of FFLCs, but few of them have been applied to hot forming processes . In the existing models for FFLCs, various ductile fracture criteria which are not coupled with other parameters have been commonly used to model fracture behaviour in the material (Bao and Wierzbicki, 2004;Cao et al., 2018;Han and Kim, 2003;Li et al., 2020;Liu and Fu, 2014;Lou and Huh, 2013;Park et al., 2017;Tang et al., 1999;Wierzbicki et al., 2005). Since strain state has a significant effect on damage evaluation (Foster et al., 2011;Kaye et al., 2016), based on the McClintock fracture criterion (McClintock, 1968), Atkins (1996) constructed analytically a fracture locus (i.e. an FFLC) consisting of a line with a slope equal to À1 in the coordinate system of major strain e 1 against minor strain e 2 . In later modelling work, fracture loci with a similar monotonic linear dependence between e 1 and e 2 were developed (Lee et al., 2004;Mahalle et al., 2020;Mu et al., 2020); these were consistent with the experimental results obtained in studies on a variety of sheet materials (Gu¨ler and Efe, 2018;Isik et al., 2014;Jawale et al., 2018;Lee and Wierzbicki, 2005;Soeiro et al., 2015;Song et al., 2017). In other experimental studies, however, it was found that the FFLC had a similar shape to the FLC, with the FFLC approaching the FLC as the strain path approaches the equal-biaxial strain state (Jain et al., 1999;Tasan et al., 2009). The experimental studies cited above used a variety of methods for measuring fracture strains. Since the strain distribution around the fracture location is highly nonuniform at the point of fracture, the strain values determined are highly dependent on the method of measurement, which has not yet been standardised. We believe that this is the main reason for the observed discrepancies in the shape of experimentally determined FFLCs.
Furthermore, no single model has been developed that is capable of predicting both FLCs and FFLCs for sheet metals under hot stamping conditions. This is partly due to the lack of available experimental data on FLCs and FFLCs for the calibration of material models. Recently, both FLCs and FFLCs for 22MnB5 boron steel under hot stamping conditions were determined for the first time using a biaxial test method developed by Zhang et al. (2022). This provides fundamental experimental data on FLCs and FFLCs for the calibration of the model developed in the present study.
The aim of this study is to formulate a set of CDM-based unified viscoplastic constitutive equations for modelling both FLCs and FFLCs for boron steel for hot stamping applications. This model allows prediction of whether the material has passed the necking strain limit, and also estimation of how far the material is away from the fracture strain limit. For this purpose, two novel coupled damage variables are proposed to model the accumulated damage in the material leading to localised necking and fracture respectively. The model incorporates these new damage variables in combination with a set of constitutive equations for dislocation-based hardening. The material constants in the CDM-based constitutive equations are determined by curve fitting of the associated experimental data, i.e. true stress-true strain curves and limit strains for the boron steel under hot stamping conditions. Furthermore, the effects of varying selected parameters in the damage variables on the functional form and numerical values of both FLCs and FFLCs computed using the CDM-based equations are investigated and discussed.

CDM-Based unified viscoplastic constitutive equations for FLCs and FFLCs
This section presents a CDM-based unified viscoplastic constitutive model for FLCs and FFLCs. This model is formulated on the basis of the recently developed unified constitutive equations for dislocation-based hardening (Lin et al., 2014), coupled with two novel damage variables developed in the present study for modelling necking and fracture, respectively. Considering the relatively low thickness of the material used in hot stamping processes (e.g. 1.5 mm for the material in the present study), a plane-stress condition is assumed, i.e. r 3 ¼ 0, where r 3 is the stress in the through-thickness direction of the sheet.

Dislocation-based unified constitutive equations for hardening
Dislocation-based unified constitutive equations for hardening were originally developed for the uniaxial stress case (Lin and Liu, 2003), and then extended to multiaxial stress states (Lin et al., 2014). These equations are introduced briefly in the subsections below to enable ease of understanding of the subsequent implementation and calibration of the CDM-based constitutive model developed in the present study. Specifically, the constitutive equations for the uniaxial stress state will be calibrated using true stress-true strain curves in Subsection Dislocation-based hardening constitutive equations (Uniaxial), and those for the multiaxial stress states will be calibrated using limit strains at the onsets of necking and of fracture in Subsection CDM-based unified viscoplastic constitutive equations (Multiaxial).
Uniaxial relations. Under hot stamping conditions, the flow stress r can be divided into three parts: the stress r V due to viscoplasticity, the stress r R due to isotropic hardening, and the initial yield stress r k . Assuming a power-law relationship between the stress r V and the plastic strain rate _ e p , the latter can be formulated as : where r k , K and n 1 are temperature-dependent parameters. It is well known that the degree of isotropic hardening is directly related to the dislocation density. Based on this observation, an isotropic hardening equation was first developed as a function of average dislocation density by Sandstr€ om and Lagneborg (1975). In solving initial value problems, however, it is usually difficult to determine the initial value of the dislocation density which varies with material, chemical composition and processing route. In order to overcome this difficulty, a normalised dislocation density parameter q was originally proposed by Lin and Liu (2003); q ¼ ðq À q i Þ=q s , where q i is the initial dislocation density and q s is the saturation dislocation density. After the average dislocation density has been replaced with the normalised dislocation density, the modified dislocation-based isotropic hardening can be expressed as (Lin and Liu, 2003): where q is the normalised dislocation density which varies from 0 to 1 and B is a temperaturedependent parameter.
The dislocation density is a result of contributions from multiple mechanisms associated with high-temperature deformation. Equation (3) is an expression for dislocation density taking into account static and dynamic recovery (Garrett et al., 2005), in which the first term A _ e p represents the rate of increase of dislocation density due to plastic deformation, and the second term A q _ e p and the third term Cq n 2 describe the rate of decrease of dislocation density due to dynamic and static recovery, respectively: where A and n 2 are material constants, and C is a temperature-dependent parameter.
Hooke's law is used to describe the linear relationship between the flow stress and the elastic strain: where e is the total strain, e p the plastic strain and E the elastic modulus.
Multiaxial relations. At hot stamping temperatures, 22MnB5 boron steel sheet exhibits only very weak anisotropy, such that isotropic behaviour can be assumed for simplicity (Merklein and Lechler, 2009;Turetta et al., 2007). Therefore, the Von Mises yield criterion is used to describe the yield surface fðr ij Þ of the material in this study; this is given by: where r f is the flow stress which depends on the degree of plastic deformation, r ij are the stress components and S ij the deviatoric stresses which are given by: where d ij is the Kronecker delta, and r kk are the principal stresses which obey the summation convention.
Recalling the associated flow rule that describes the relationship between the yield surface fðr ij Þ and the increment of plastic strain component de p ij and its direction, the rate _ e p ij can be expressed as (Bland, 1957): where _ e p ij is the rate of change of the plastic strain component, dK a non-negative infinitesimal, t the time, and _ e p e the effective plastic strain rate. The multiaxial power-law viscoplastic equations are obtained by consideration of a dissipation potential function (Lin et al., 2014). With the initial yield stress r k , the energy dissipation can be expressed as: where r e is the Von Mises stress given by: Based on the associated flow rule (equation (7)), the strain rate equation then becomes: By introducing the isotropic hardening r R , the effective plastic strain rate _ e p e for the power-law viscoplastic material model can be written as (Lin et al., 2014): Note that equation (11) is applicable only when r e À r R À r k > 0. Otherwise _ e p e ¼ 0. The isotropic hardening under multiaxial stress conditions takes the same form as equation (2). The normalised dislocation density q is obtained by replacing the plastic strain rate _ e p with the effective plastic strain rate _ e p e . The generalised Hooke's law, which defines the general linear relation among all the components of the stress and strain tensors, is given by: where e ij are the components of the total strain and D ijkl the components of the fourth-order tensor of elastic moduli given by: where k and l are the Lam e constants expressed as: where is Poisson's ratio and G the shear modulus. In the present study, the value of Poisson's ratio is taken as ¼ 0.3 which is a widely accepted value for steels, e.g. (Gozzi et al., 2005;Shapiro, 2009;Song et al., 2017). Lin et al. (2014) were the first to develop a CDM-based constitutive model to predict FLCs for sheet metals under hot stamping conditions under the assumption of plane-stress conditions, coupled with a damage variable to characterise the extent of necking. Since an FLC is constructed using limit strains, Mohamed et al. (2015) modified this damage variable, replacing the stress components with strain components to model the effect of the stress or strain state on the damage evolution. The modified damage variable in multiaxial form is expressed as (Mohamed et al., 2015):

CDM-Based damage variables
where e 1 is the major strain, e 2 is the minor strain, e p e is the effective plastic strain, l 1 and g 2 are parameters which are dependent on temperature, and l 2 , /, g 1 and g 3 are material constants; c is a small constant introduced to avoid a denominator of zero, and takes the value of 0.0035 in the present study.
This modified damage variable x in equation (16) has previously been applied to predict FLCs for AA5754 under warm forming conditions (Mohamed et al., 2015), and for AA6082 under hot stamping conditions (Shao et al., 2017). It has also been shown that the damage variable has sufficient flexibility for the prediction of both the numerical values and the shapes of FLCs under a range of hot stamping conditions (Shao et al., 2017). However, attempts by some of the present authors to predict both FLCs and FFLCs for sheet metals under hot stamping conditions using critical values of the single damage variable in equation (16) have proved unsuccessful, demonstrating that a different approach is required. Previously, Othman et al. (1993) developed two coupled damage parameters to model the tertiary creep softening caused by multiplication of dislocations, and nucleation and growth of cavities near grain boundaries. These two damage parameters have been widely used for numerical analysis of the creep-rupture behaviour (Hayhurst et al., 2008).
Based on the above-mentioned studies, in the present work, two novel coupled damage variables x 1 and x 2 have been developed to model damage in the material leading to localised necking (i.e. FLCs) and to fracture (i.e. FFLCs), respectively. The novel damage variables x 1 and x 2 are expressed as: where D 1 , D 2 , / 1 and / 2 are temperature-dependent parameters.
Similarly to the damage variable of equation (16), the first two factors of equations (17) and (18) (square and round brackets) describe the effects of the strain state on the damage evolution in the material, with the coefficients l 1 and l 2 characterising the effects of the major and minor principal strains, respectively. The exponents / 1 and / 2 describe the intensity of the effects of the principal strains on the damage evolution. The third factors in each of equations (17) and (18) show the interaction between the two coupled damage variables. The fourth factor takes into account the effects of the temperature (via the temperature-dependent parameter g 2 ) and the strain rate on the failure evolution.
By introducing the two novel damage variables x 1 and x 2 , the effective plastic strain rate in equation (11) can be expressed as: The damage variable x 1 lies in the range 0 x 1 < 1; x 1 is zero at the start of deformation, and it is assumed that localised necking occurs when x 1 reaches 0.7. This critical value, in conjunction with suitable calibration of the material constants, has been used in the models of FLCs developed previously to avoid x 1 reaching unity (at which the model equations would become meaningless) and to save computing time (Lin et al., 2014;Mohamed et al., 2012). The damage variable x 2 has the same range as x 1 , and for convenience, it is also assumed that fracture occurs when x 2 reaches 0.7. When D 1 ¼ D 2 ¼ 0, the damage variables x 1 and x 2 become zero, and equation (19) reduces to equation (11).

CDM-based unified viscoplastic constitutive equations
In summary, a set of CDM-based unified viscoplastic constitutive equations has been developed for the modelling of both FLCs and FFLCs for sheet metals under hot stamping conditions, which is written as: In equation (20), the parameters r k , K, n 1 , B, C, E, D 1 , l 1 , / 1 , g 2 , D 2 and / 2 are temperaturedependent; r k , K, n 1 , B, E, D 1 , g 2 and D 2 have the form: where R is the universal gas constant, T the absolute temperature and X represents each of the parameters r k , K, n 1 , B, E, D 1 , g 2 and D 2 . The new material constants are introduced: r k0 , Q r k , K 0 , Q K , n 1 0 , Q n 1 , B 0 , Q B , E 0 , Q E , D 1 0 , Q D 1 , g 2 0 , Q g 2 , D 2 0 and Q D 2 . C, l 1 , / 1 and / 2 have the form: where Y represents each of the parameters C, l 1 , / 1 and / 2 ; the new material constants are introduced: C 0 , Q C , l 1 0 , Q l 1 , / 1 0 , Q / 1 , / 2 0 and Q / 2 .

Model implementation and discussion
Determination of material constants from experimental data Experimental data on 1.5 mm thick 22MnB5 boron steel sheet at hot stamping temperatures reported by some of the present authors are used for the determination of the material constants arising in equations (20) to (22). These data comprise true stress-true strain curves (Zhang et al., 2020), and limit strains at the onset of localised necking (i.e. FLCs) and at fracture (i.e. FFLCs) (Zhang et al., 2022) obtained under the conditions given in Table 1. These true stress-true strain curves of boron steel were obtained by carrying out uniaxial tensile tests under hot stamping conditions using a Gleeble thermal-mechanical simulator, together with the digital image correlation (DIC) technique for full-field strain measurement. To obtain experimental data for the construction of FLCs and FFLCs, cruciform specimens were deformed in different strain states and full-field strains within the gauge area were measured using the DIC technique, followed by the determination of the limit strain values using a spatio-temporal method originally developed by Zhang et al. (2021). The calibration process is divided into two steps. The first step is to determine the material constants within the dislocation-based hardening constitutive equations (equations (1) to (4)) using the true stress-true strain curves for the material. The damage variables can be neglected in this step since their values will be vanishingly small and have little influence on the viscoplastic deformation of the material in the initial stages of deformation. The second step is to determine the material constants in the equations for the two damage variables (equations (17) and (18)) using the experimentally determined limit strains at the onsets of localised necking and of fracture. The detailed procedures of the calibration are given in the following sections.
Dislocation-based hardening constitutive equations (uniaxial). In order to calibrate the material constants in the dislocation-based hardening constitutive equations, equations (1) to (4) were solved using a numerical integration method combined with the explicit Euler method. There are a relatively large number of material constants (14 in total), which causes difficulties in the curve fitting. Many studies have been carried out to develop mathematical fitting methods to overcome the difficulties (Abushawashi et al., 2013;Cao and Lin, 2008). In the present study, an optimisation method based on the trust-region-reflective least square algorithm, implemented by means of the Lsqnonlin solver in MATLAB (MathWorks, USA), was used for the curve fitting.
In the numerical computations, initial values were given to each of the material constants and upper and lower boundaries set, all according to experience. The time increment dt was set to 3.0E À 3, 3.0E À 4, and 3.0E À 5 for the strain rates of 0.02, 0.2 and 2 s À1 , respectively. The material constants determined are given in Table 2. Figure 1 shows a comparison of the true stress-strain curves computed using the calculated material constants and the corresponding experimental data at the indicated temperatures and strain rates. Good agreement can be seen.

CDM-Based unified viscoplastic constitutive equations (multiaxial).
In the second step, the CDM-based constitutive equations summarised in equation (20) were also numerically solved using the explicit Euler method. Figure 2 shows the associated flow chart for the numerical computations. At the start of each computation (j ¼ 1) under conditions of constant temperature T and effective strain rate _ e e , a non-zero small initial value (i.e. 1.0E À 08) was given to the variables such as time, strain and stress components (box A on the flow chart), to avoid the occurrence of a zero denominator. In the next time increment dt (j ¼ j þ 1), the rates of the effective plastic strain, the dislocation density and the isotropic hardening were computed (box B), and the rates of the plastic strain components and the two damage variables were obtained (box C). Using the explicit Euler method, the values of the associated variables, such as the plastic strain components and the damage variables, were then calculated for this time increment (box D). Given the effective strain rate _ e e and the strain ratio b, the principal strain components (box E) were computed using the  (25). Based on the assumption of the plane-stress condition, one has: where e 1 , e 2 and e 3 are the principal strain components. Under proportional straining conditions, the strain ratio b is given by: Furthermore, the effective strain e e , which can be computed as the product of the effective strain rate _ e e and the instantaneous time t, has the form:  The principal strain components were computed by solving equations (23) to (25) using the Newton-Raphson method, and the principal stress components were then obtained according to the generalised Hooke's law (box E). The above computations were repeated (j ¼ j þ 1) until the damage variable x 1 reached the critical value of 0.7, at which point the limit strains for necking were determined (box F), and then allowed to continue until x 2 reached 0.7, at which the limit strains for fracture were determined (box G). The material constants in the CDM-based constitutive equations, i.e. equations (20) to (22), were determined by fitting the experimental FLCs and FFLCs. It should be noted that the material constants determined in the first step (Table 2), for the uniaxial dislocation-based hardening constitutive equations in Section Dislocation-based hardening constitutive equations (Uniaxial), were retained. Thus, there are 15 material constants which need to be determined in this second step.
Due to the current limitations of the biaxial test rig, which is able to produce three distinct displacement ratios only, only a few experimental limit strain data points are available for constructing FLCs and FFLCs (Zhang et al., 2022). Therefore, the curve fitting was performed by trial and error, in combination with a thorough understanding of the physical meaning of each variable and associated material constants in the equations. Specifically, the following substeps were performed: • Substep 1: According to the form of equations (17) and (18), in the uniaxial stress state (b ¼ À0.5), the parts of the equation involving the temperature-dependent parameters l 1 , / 1 and / 2 and the material constant l 2 disappear, and these do not affect the computed limit strains. Experimental limit strain data obtained in uniaxial tension were therefore chosen for this first step of fitting; data from testing at 800 C and at 0.02, 0.1, 0.5 s À1 , were used to determine the values of the temperature-dependent parameters D 1 , D 2 and g 2 at 800 C and the material constant g 1 . • Substep 2: The values of the temperature-dependent parameters l 1 , / 1 and / 2 at 800 C and the material constants c and l 2 were then determined by fitting the experimental data obtained at 800 C and at 0.02, 0.1, 0.5 s À1 in the other strain states, i.e. the plane-strain (b ¼ 0) and equibiaxial (b ¼ 1) states. The data used for the fitting in both steps 1 and 2 were obtained at the single temperature of 800 C, so the effect of temperature on FLCs and FFLCs was not taken into account in this step. • Substep 3: The material constants in the temperature-dependent parameters (i.e. D 1 , D 2 , g 2 , l 1 , / 1 and / 2 ) in equations (21) to (22) were determined by fitting the experimental data obtained at the other temperatures (i.e. 750 and 850 C) and retaining the values of these parameters D 1 , D 2 , g 2 , l 1 , / 1 and / 2 determined at 800 C in steps 1 and 2.
In the above numerical computations of FLCs and FFLCs, the time increment dt was set to 4.0E À 4, 8.0E À 5 and 1.0E À 5 by trial and error, for fitting datasets with effective strain rates of 0.02, 0.1 and 0.5 s À1 , respectively. The material constants determined in this section are listed in Table 3. Figure 3 shows a comparison of the FLCs (solid lines) and FFLCs (dashed lines) computed using the material constant values given in Tables 2 and 3, with the associated experimental data (symbols) for the limit strains at the different temperatures (i.e. 750, 800 and 850 C, Figure 3(a)) and strain rates (i.e. 0.02, 0.1 and 0.5 s À1 , Figure 3(b)). Good agreement can be observed between the computed and experimental results. Figure 4 shows the evolution of the damage variables x 1 and x 2 in the computation of the FLC and the FFLC at 800 C and 0.1 s À1 in the three strain states: uniaxial (b ¼ À0.5), plane-strain (b ¼ 0) and biaxial (b ¼ 0.5). As can be seen, in all cases both x 1 and x 2 increase with increasing time, with x 1 always higher than x 2 . Furthermore, x 2 is relatively small when x 1 reaches the critical value of 0.7 that is here taken to represent the onset of localised necking; beyond this, x 2 increases sharply and reaches its critical value of 0.7 that is taken to correspond to fracture. The times taken to reach localised necking t LN and fracture t F in plane-strain stretching (Figure 4(b)) are shorter Table 3. Material constants determined in the CDM-based unified viscoplastic constitutive equations. D 10 (À) Q D1 (J/mol) l 10 (À) Q l 1 (J/mol) l 2 (À) / 10 (À) than those obtained either in uniaxial (Figure 4(a)) or biaxial stretching (Figure 4(c)). This agrees with the experimental observations in the formability tests in Zhang et al. (2022) used for the calibration of the CDM-based constitutive model.

Effects of novel damage variables on computed FLCs and FFLCs
This subsection aims to show the capability of the novel coupled damage variables to model FLCs and FFLCs of various shapes and numerical values. For this purpose, the effects of varying intentionally the values of selected parameters, i.e. l 1 , l 2 , / 1 , / 2 , D 1 and D 2 in the novel damage variables from equations (17) and (18), on both the shapes and values of the FLCs and FFLCs computed using the CDM-based model were investigated. The values chosen for each parameter, which are given in Table 4, include a 'basic' value, which is based on those obtained by calibration using the experimental data from the boron steel, and alternative values that are larger and smaller than this basic value. The other material constants determined in Table 2 and Table 3 are retained. An FLC and an FFLC were first computed using the basic values given in Table 4 (i.e.
The parameter l 1 was then changed to the lower alternative value (i.e. l 1 ¼ 0.3) in Table 4, maintaining the basic values for all the other parameters (i.e. l 2 ¼ 0.12, , and the computation was repeated. This process was repeated for the higher alternative value of l 1 (i.e. l 1 ¼ 0.5), and then for the lower and higher values of all the other parameters in Table 4. All computations were carried out for the same temperature of 800 C and the same effective strain rate of 0.1 s À1 , so as to avoid potential influence on any other parameters that are dependent on either temperature or strain rate.
Parameters l 1 and l 2 . Figure 5 shows the FLCs and FFLCs computed by varying the parameters l 1 and l 2 separately. In Figure 5(a), it can be seen that increasing l 1 from 0.3 to 0.5 leaves the limit major strain values in the uniaxial strain state (b ¼ À0.5) unaffected in both the FLCs and FFLCs. The limit major strain in the other strain states increases, with the increase becoming more significant at higher strain ratios, especially at b ¼ 1. Figure 5(b) shows that increasing l 2 from 0.1 to 0.15 has a similar effect to that of decreasing l 1 . However, referring to Figure 5, the computed FLCs and FFLCs are more sensitive to l 2 than l 1 in the ranges investigated with respect to the effect on the limit major strains.  Parameters / 1 and / 2 . Figure 6(a) presents the FLCs and FFLCs computed with / 1 values increasing from 4 to 8. It can be seen that the / 1 value has very little effect on the limit major strain values at the two extremes of the curves (strain ratios close to either À0.5 or 1) but in the centres, with strain ratios close to 0, increasing / 1 decreases the limit major strain. Figure 6(b) shows the effect of increasing / 2 from 1.5 to 4.5. The / 2 value affects the shapes of the FFLCs in a similar manner to / 1 . However, the variation of / 2 has almost no influence on the FLCs in the range investigated.
Parameters D 1 and D 2 . Figure 7 shows the FLCs and FFLCs computed using different values of the parameters D 1 and D 2 . In Figure 7(a), it can be seen that increasing D 1 from 3 to 5 decreases the limit major strains across the complete range of strain ratios b for both FLCs and FFLCs, i.e. both FLCs and FFLCs are shifted downwards but retain approximately the same shape. It can be seen from Figure 7(b) that increasing D 2 from 0.3 to 0.5 shifts the FFLCs downwards while keeping the curves approximately parallel, but has little effect on the computed FLCs. It is clear that the FFLCs are more sensitive to the D 2 value than to the D 1 value in the ranges investigated.

Discussion
The CDM-based unified constitutive equations are capable of modelling the limit strains at the onset of localised necking and at fracture for boron steel for hot stamping applications (Figure 3). According to Figures 5 to 7, the parameters (i.e. l 1 , l 2 , / 1 , / 2 , D 1 and D 2 ) appearing in the novel coupled damage variables provide sufficient flexibility to modify the shape and numerical value of the FLCs and FFLCs computed using the CDM-based constitutive model, such that the experimental limit strain data for the FLCs and FFLCs for the boron steel can be reproduced. It should be noted that thanks to its flexibility, the CDM-based model is also able to predict FFLCs with different shapes, according to the shape of the experimental data used for the model calibration. In the present study, the experimental data obtained by biaxial testing give a similar shape for the onset of fracture as for the onset of necking, and the fitted model reflects this. Importantly, the flexibility indicates that the CDM-based unified viscoplastic constitutive model developed in this study has the potential for application to other materials under various warm-or hot-stamping conditions. Further work is needed to verify this by applying the CDM-based model to other materials (e.g. aluminium alloys) for formability prediction under hot stamping conditions. In addition, the CDMbased constitutive model developed in the present study incorporates dependence on strain state, temperature and strain rate, so the model is able to take account of the history of strain state, temperature and strain rate rather than just the instantaneous strain state, temperature and strain rate during the forming processes. It is expected that this model will be able to predict the formability of sheet materials in practical sheet metal forming processes in which the forming conditions (e.g. strain path, temperature and strain rate) continuously vary as a function of location on the workpiece. However, this has not yet been validated in the present study due to the lack of availability of experimental data on sheet boron steel under conditions equivalent to industrial hot stamping. Therefore, additional work is needed to demonstrate the capability of the model for  flexible prediction of formability. It is interesting to compare the CDM-based constitutive model in this study with the modified Hosford-Coulomb model developed by Roth and Mohr (2014) in respect to the modelling of fracture initiation. In the latter, the combined effects of strain rate and temperature on damage evolution are taken into account on an empirical basis. In the CDM-based model, in contrast, hardening constitutive equations are expressed in more explicitly physically based terms of dislocation density, and the effects of strain path, temperature and strain rate on the damage evolution are formulated separately, but at a cost of introducing a larger number of material constants.
Since there are a relatively large number of material constants in the CDM-based constitutive equations, the determination of these material constants is complex and difficult. According to Figure 1, the optimisation method based on the trust-region-reflective least square algorithm enables the material constants in the dislocation-based hardening constitutive equations (uniaxial) to be determined such that experimental stress-strain curves are reproduced accurately. However, one of the difficulties of using this method is the need to define suitable initial values and lower and upper bounds for each of the material constants. Furthermore, the material constants in the CDM-based damage variables were successfully determined by trial and error ( Figure 3). However, the determination of material constants requires a deep understanding of the physical meaning of each of the variables and the associated material constants. Therefore, it is necessary to develop a simpler but still effective method for determining the material constants in the CDM-based constitutive equations, especially for industrial applications.

Conclusions
A set of CDM-based unified viscoplastic constitutive equations has been developed, for the first time, to model both FLCs and FFLCs for boron steel (22MnB5) sheet for hot stamping applications. These comprise a set of dislocation-based unified constitutive equations for hardening, and two novel coupled damage variables x 1 and x 2 to characterise the accumulated damage leading to localised necking and fracture, respectively. With increasing deformation, both x 1 and x 2 increase and reach critical values for the onsets of localised necking and fracture, respectively. The CDMbased unified viscoplastic constitutive equations have been shown to have the capability of modelling both FLCs and FFLCs for boron steel under hot stamping conditions. The two damage variables which are dependent on strain rate and temperature take into account the influence of the strain state on formability, thereby providing sufficient flexibility to precisely reproduce the various shapes and numerical values of experimental FLCs and FFLCs. This indicates that the CDM-based unified viscoplastic constitutive model has great potential for application to other materials for modelling the formability under warm or hot stamping conditions.
The material constants in the CDM-based constitutive equations have been calibrated from the curve fitting of associated experimental data from true stress-true strain curves and from FLCs and FFLCs obtained at different temperatures and strain rates. However, the calibration requires a thorough understanding of the equations and considerable practical experience. For industrial applications, a simpler but effective method of determining the material constants is necessary.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Engineering and Physical Sciences Research Council (EPSRC) [grant number EP/R001715/1] on "Lightform: Embedding Materials Engineering in Manufacturing with Light Alloys"; and the Chinese Scholarship Council (CSC) Imperial Scholarship [grant number 201700260069].