# Inverse Strategies for Identifying the Parameters of Constitutive Laws of Metal Sheets

Authors: P. A. Prates; A. F. G. Pereira; N. A. Sakharova; M. C. Oliveira; J. V. Fernandes
Journal: Advances in Materials Science and Engineering
 (2016)

---

## Abstract

This article is a review regarding recently developed inverse strategies coupled with finite element simulations for the identification of the parameters of constitutive laws that describe the plastic behaviour of metal sheets. It highlights that the identification procedure is dictated by the loading conditions, the geometry of the sample, the type of experimental results selected for the analysis, the cost function, and optimization algorithm used. Also, the type of constitutive law (isotropic and/or kinematic hardening laws and/or anisotropic yield criterion), whose parameters are intended to be identified, affects the whole identification procedure.


---

## Body

## 1. Introduction

Finite Element Analysis (FEA) is now a well-established computational tool in industry for the optimization of sheet metal forming processes. The accurate modelling of these processes is a complex task due to the nonlinearities involved, such as those associated with (i) the kinematics of large deformations, (ii) the contact between the sheet and the tools, and (iii) the plastic behaviour of the metal sheet.The description of the plastic behaviour of metal sheets is usually performed using phenomenological constitutive models. In this context, the emergence of new steels and aluminium, magnesium, and other alloys, as well as their increasingly widespread use in the automotive and aeronautical industries, has encouraged the development of more reliable models, with increasing flexibility associated with a larger number of parameters to identify [1–14]. In fact, the accuracy of the numerical simulation results of sheet metal forming processes depends on the flexibility of a constitutive material model but also on the procedure adopted to identify its parameters. The complex nature of the plastic behaviour of metal sheets makes their characterization dependent upon factors such as (i) the constitutive model; (ii) the experimental tests performed, comprising the sample geometry, the testing conditions, and the analysis methodologies; (iii) the strategy for identifying the constitutive parameters.The strategy for identifying the model parameters is generally seen as an optimization problem, where the purpose is to minimise the difference between computed and experimental results of one or more experiments. Two main types of strategies for the identification of the constitutive parameters can be recognised in literature: classical and inverse strategies. The classical identification strategies for the constitutive parameters make use of a large number of standardised mechanical tests, with well-defined geometry and loading conditions, such that homogeneous stress and strain distribution develop in the region of interest (e.g., [15, 16]); nonstandard mechanical tests can also be performed to properly describe other biaxial stress states in the sheet plane (e.g., [17, 18]). However, sheet metal forming processes are carried out under strongly nonhomogeneous stresses and strains fields. Therefore, limiting the characterization of the mechanical behaviour of metal sheets to a restricted number of tests with linear strain paths and homogeneous deformation can lead to a somewhat incomplete characterization of the overall plastic behaviour of the material [19].Recent developments and accessibility of optical full-field measurement techniques, such as digital image correlation (DIC) technique coupled with FEA, make the inverse identification strategies a common current place. The full-field measurements allow the acquisition of enriched information from mechanical tests, such as displacement and strain fields; an overview on this topic can be found in [20]. This allows attenuating the constraints on the geometry and loading conditions of the mechanical tests used for the identification of materials parameters, so that nonhomogeneous stress and strain distributions can be developed in the region of interest (e.g., [21–29]). In this sense, the identification of constitutive parameters from nonhomogeneous strain fields and complex loading conditions provides a more reliable description of the material behaviour during real sheet metal forming processes [21]. In such complex mechanical tests, it is no longer possible to identify the constitutive parameters based on simple assumptions on the stress and/or strain states, as in the classical identification strategies. Instead, a finite element model of the mechanical test is established and cost functions are defined to minimise the gap between numerical and experimental results of the mechanical test, which demands efficient optimization algorithms. However, the efficiency of any inverse identification strategy directly depends on the information contained by the objective function. In the context of constitutive parameters identification, this is related to the type of experimental results included (e.g., loads, displacements, and strains) but also to the strain paths and levels of deformation attained by the experimental test. It turns out that there is no consensus about the experimental tests (sample geometry and loading conditions), the cost functions, and the optimization procedure that will lead to accurate constitutive parameters identification. Also, a major obstacle to the widespread use of advanced constitutive models in industrial simulations seems to result from the lack of an efficient strategy for parameters identification. In this sense, the developed strategy must be simple, from an experimental point of view, and allow evaluating to what extent the selected constitutive model allows perfectly describing the behaviour of a given material.The present paper describes recent inverse strategies coupled with FE simulations for the identification of the parameters of constitutive laws that describe the plastic behaviour of metal sheets, resorting to mechanical tests leading to nonuniform strain and stress states. Following this introduction, the paper addresses general concepts for the constitutive modelling and the optimization problem. Afterwards, an overview of inverse identification strategies for the constitutive parameters is presented, with emphasis on inverse identification strategies resorting to FE simulations.


## 2. Constitutive Modelling

Constitutive models have been developed to predict the onset and evolution of the plastic deformation of a deformable body undergoing a general state of stress. A phenomenological constitutive model is typically a combination of the following components:(i)
Yield criterion that describes the yield surface of the material in a multidimensional stress space: The metal sheets are usually assumed to be orthotropic, with invariant anisotropy during plastic deformation. With high incidence in the last decades, the emergence of anisotropic yield criteria with an increasing number of material parameters has been witnessed. They provide the flexibility required for accurately modelling the plastic behaviour of advanced metallic alloys, which are frequently used in automotive and aeronautical industries. Several approaches have been used for deriving yield criteria, based on
(1)
high-order polynomial functions (e.g., [1, 2]);
(2)
the generalization to anisotropy of the second and third invariants of the deviatoric stress tensor,J
2 and J
3, respectively (e.g., [3]);
(3)
one or more isotropic yield functions, using the linear Isotropic Plasticity Equivalent (IPE) stress space concept (e.g., [3–9]);
(4)
the construction of weighted sums of anisotropic yield criteria (e.g., [7]);
(5)
the capability to model the tension-compression asymmetry, particularly devoted to specific magnesium and titanium alloys (e.g., [3, 6, 10]);
(6)
the capability to model kinematic hardening [11];
(7)
the interpolation of second-order Bézier curves [12].
(ii)
Hardening laws that express the evolution of the yield surface during plastic deformation, as schematized in Figure1: The isotropic hardening law refers to the homothetic expansion of the yield surface (see Figure 1(a)) while the kinematic hardening law describes its translation in the stress space (see Figure 1(b)). Kinematic hardening laws are recommended for describing plastic deformation under strain path changes, mainly strain path reversal, in materials that exhibit Bauschinger effect (e.g., [14]). The combination of isotropic and kinematic hardening laws provides a flexible model, for simultaneously describing the change in size and the position of the centre of the yield surface, during plastic deformation. Isotropic hardening laws described by power laws (e.g., [32–37]), saturation laws (e.g., [38, 39]), and weighted combinations of isotropic hardening laws (e.g., [40, 41]) have been proposed. Linear (e.g., [42, 43]) and nonlinear (e.g., [13, 14, 44–47]) kinematic hardening laws were proposed, with the latter being more appropriate to describe the Bauschinger effect.
(iii)
Flow rule, to establish a relationship between the stress state and the plastic strain increment: Typically, an associated flow rule is adopted, that is, using the yield function as plastic potential, although some exceptions can be found in literature (see, e.g., [48]).The general representation of a constitutive model can be described through a function F
=
F
σ
′
-
X
′
,
ε
-
p
,
α
,
β:(1)
F
σ
′
-
X
′
,
ε
-
p
,
α
,
β
=
σ
-
σ
′
-
X
′
,
α
-
Y
ε
-
p
,
β
,where σ
-
(
σ
′
-
X
′
,
α
) is the equivalent stress defined by a given yield criterion and Y
(
ε
-
p
,
β
) is the hardening law that represents the evolution of the yield stress during the deformation. The equivalent stress, σ
-
(
σ
′
-
X
′
,
α
)
=
σ
-, is a function of the effective stress tensor, (
σ
′
-
X
′
), that includes the parameters of the yield criterion, α, for describing the anisotropy (σ
′ and X
′ are the deviatoric Cauchy stress and the deviatoric backstress tensors, resp.) and Y
(
ε
-
p
,
β
)
=
Y is a function of the equivalent plastic strain, ε
-
p, in which the parameters are represented by β. The yielding is defined based on the function F of (1) and can be written as follows:(2)
F
=
σ
-
-
Y
=
0
⟹
σ
-
=
Y
.If σ
-
<
Y, the stress state of the material remains inside the yield surface and only elastic deformation occurs. When plastic deformation occurs, the associated flow rule states that the increment of the plastic strain tensor is normal to the yield surface, for a stress state such that σ
-
=
Y. The normality condition, defined by the associated flow rule, assumes that the increment of the plastic strain tensor is normal to the yield surface and is expressed by(3)
d
ε
p
=
d
λ
∂
σ
-
σ
′
-
X
′
∂
σ
′
-
X
′
,where d
ε
p is the increment of the plastic strain tensor, d
λ is a scalar multiplier, and σ
-
σ
′
-
X
′
=
σ
- is the equivalent stress function, representing the plastic potential.Figure 1
Representation of the plastic behaviour of materials in tension-compression. The left side shows generic yield surfaces in the plane (σ
1; σ
2) and the right side shows the corresponding hardening law: (a) isotropic hardening and (b) kinematic hardening. See text for details. Adapted from [30].
(a)
(b)Even though a number of advanced constitutive models are available in literature, sheet metal forming simulations are still mostly performed in industry not taking into account kinematic hardening and with the well-known Hill’48 yield criterion [49], whose parameters identification can be easily assessed by uniaxial tensile tests. Mattiasson and Sigvant [50] mentioned some plausible explanations, still valid today, for this reality:(i)
The relative simplicity of the Hill’48 model that makes it attractive to use
(ii)
The unavailability of industry analysts for understanding to what extent the modelling of the material influences the simulation results
(iii)
The lack of knowledge, time, and money for performing the multiaxial tests required to identify reliable hardening curves and parameters of advanced yield criteria
(iv)
The additional cost in terms of CPU time for using more advanced constitutive models which is considered to be an effort that is not worth it.Nevertheless, in our view, the major obstacle to the widespread use of advanced constitutive models in industrial simulations comes from the large number of linear strain path tests, including multiaxial tests, required for the parameters identification. To overcome this barrier, a potential approach is to look for new constitutive parameters identification strategies that are alternative to the classical ones. In this sense, an accurate description of the material plastic behaviour could be attained from (i) a minimum number of mechanical tests and experimental data; (ii) flexible and user-friendly constitutive models; and (iii) an accessible identification procedure for the constitutive parameters, coupled with robust optimization algorithms. Therefore, Section 4 will discuss some identification procedures for the constitutive parameters based on inverse analysis, as an alternative to the classical approaches.


## 3. The Optimization Problem

The inverse identification of constitutive model parameters is generally seen as an optimization problem. The purpose is to minimise the difference between computed and experimental results of one or more experiments. This difference is expressed by a cost function and its minimisation is performed using optimization algorithms, which automatically operate on the values of the constitutive parameters.


### 3.1. Cost Function

A wide number of cost function formulations for the identification of constitutive parameters have been proposed in literature (e.g., [57, 58]). According to Cao and Lin [57], the cost function should operate as an “efficient guide” of the optimization procedure, in order to search for the best fit to the experimental results; therefore, the ideal cost function should comprise the following conditions:(i)
All measured points of a given experiment should be part of the optimization procedure and have equal opportunity to be optimized, provided that experimental errors are eliminated.
(ii)
All experiments should have equal opportunity to be equally optimized, and so the optimization should not depend on the number of points considered in each experiment.
(iii)
Different units of measure in the cost function should not affect the performance of the optimization.
(iv)
The identification procedure should not be dependent of the user, and so the values of the weighting factors should be optimized to achieve the abovementioned conditions.Cost functions are typically formulated under the concept of weighted least-squares, as follows: (4)
F
A
=
1
m
1
n
∑
i
=
1
m
w
i
∑
j
=
1
n
w
j
r
i
j
A
2
,where F
(
A
) is the cost function to minimise; A is the vector of constitutive parameters to optimize; m is the total number of experiments and n is the total number of points, considered in each experiment i; r
i
j
(
A
) is the residual between the numerically predicted results and those of the experiment i at point j; w
i and w
j are the weighting factors for each experiment i and for each point j, respectively. Within the context of inverse parameter identification, r
i
j
(
A
) can contain variables such as loads, pressures, angular moments, or those arising from full-field measurements (displacements or strains), as will be seen in detail later.The residuals can be expressed in terms of relative differences,(5)
r
i
j
A
=
u
i
j
N
u
m
A
-
u
i
j
E
x
p
u
i
j
E
x
p
,
w
i
t
h
i
=
1
,
…
,
m
,
j
=
1
,
…
,
n
,or in terms of absolute differences, (6)
r
i
j
A
=
u
i
j
N
u
m
A
-
u
i
j
E
x
p
,where u
i
j
N
u
m and u
i
j
E
x
p are, respectively, the numerically predicted and the experimental results at point j of experiment i. Residuals are often expressed using relative differences, which allows the use, in the same cost function, of several kinds of quantities exhibiting various orders of magnitude and units of measure [59]. When u
i
j
E
x
p admits values close to or equal to zero, the residuals should be expressed using absolute differences.


### 3.2. Optimization Algorithms

The minimization of the least-squares cost function, presented in (4), requires efficient and robust optimization algorithms, due to the strongly nonlinear nature of the least-squares cost function [60]. For this purpose, several optimization algorithms are described in the literature, which are commonly divided into two categories: gradient-free algorithms and gradient-based algorithms. Hybrid optimization strategies using both gradient-free and gradient-based algorithms are also proposed (e.g., [15, 60]).Gradient-free algorithms, such as evolutionary and SIMPLEX algorithms, have a great probability of achieving a global minimum due to their random search capability. They require a large number of cost function evaluations (i.e., iterations) and therefore the convergence can be very time-consuming. Because of this, gradient-free algorithms are not recommended within the context of inverse identification strategies, since they require a large number of finite element simulations and analyses [61].Gradient-based algorithms are most popular within inverse identification strategies, as they require far less cost function evaluations than gradient-free algorithms. As local optimizers, these algorithms use information of the gradient to update the vector of constitutive parameters in an adequate search direction [62]. Therefore, there is no guarantee that these algorithms converge to the global minimum, with the possibility of converging to undesirable local minima. This makes the optimization procedure dependent on the initial estimate for the parameters, and therefore the choice of convenient initial estimates for the constitutive parameters can be essential.Examples of gradient-based algorithms commonly used within the context of inverse identification strategies are the Gauss-Newton and Levenberg-Marquardt algorithms. The Gauss-Newton algorithm is described as follows:(7)
A
s
+
1
=
A
s
-
J
T
W
J
-
1
J
T
W
r
Α
s
,where s is the iteration step, A is the vector of constitutive parameters, W is the vector of weighting factors, J is the Jacobian matrix that expresses the sensitivity of the computed results to the constitutive parameters, and r
(
Α
s
) is the vector of residuals, which can be expressed in terms of relative or absolute differences (see (5) and (6), resp.). The dimension of the vector of residuals depends on the total number of experiments m and the total number of points n, in each experiment, that is, the dimension n
r
=
m
×
n. Considering that the total number of constitutive parameters to be identified is n
p, with n
r
≥
n
p, the Jacobian containing the partial derivatives of the residuals with respect to the constitutive parameters is defined: (8)
J
l
,
p
=
∂
r
l
A
p
∂
A
p
,
w
i
t
h
l
=
1
,
…
,
n
r
,
p
=
1
,
…
,
n
p
.An efficient method to compute the Jacobian matrix is finite differentiation. In order to improve convergence, the Jacobian matrix must be updated at each iteration step s. However, the calculation in each step requires high computational cost (at least one numerical simulation per constitutive parameter). To overcome this inconvenient, Endelt et al. [63, 64] and Cooreman [62] highlighted the possibility of computing the sensitivity matrix analytically, with the latter author having concluded the inability of this approach for computing the sensitivities of strain fields to the material parameters, in mechanical tests involving complex and/or heterogeneous deformation.In some cases, the Gauss-Newton algorithm can become unstable in the neighbourhood of the minimum, and so a stabilisation procedure is required. The Levenberg-Marquardt algorithm is similar to Gauss-Newton one, but includes a stabilising term, as follows [65]: (9)
A
s
+
1
=
A
s
-
J
T
W
J
+
λ
d
i
a
g
J
T
W
J
-
1
J
T
W
r
Α
s
,where λ is the stabilising parameter that is updated in each iteration according to the convergence rate [65]. When the Levenberg-Marquardt method shows stability, small values for λ are recommended for fast convergence; otherwise, large values of λ usually allow stable convergence, although slower, towards the minimum. Note that for λ
=
0 the Levenberg-Marquardt algorithm is equal to the Gauss-Newton one.A different type of optimization technique that has been recently used in the identification of constitutive parameters is the Response Surface Methodology (RSM) (e.g., [54, 66]). RSM is an optimization technique for generating smooth approximations of complex functions in a multidimensional design space. In the context of parameter identification, the design space contains all possible combinations for the constitutive parameters and related values of the cost function. The prohibitive size of the full design space requires a Design of Experiments (DoE), to efficiently construct an approximated design space from a few number of representative points (i.e., sets of constitutive parameters). The responses of the representative points (i.e., the values of the cost function) are used to fit a response surface, which is typically obtained from second-order polynomial regression, for the sake of simplicity. Finally, the minimum of the response surface is calculated using a gradient-based optimization procedure, which leads to an estimate of the optimal set of constitutive parameters. In brief, the RSM technique can be summarised as follows:(1)
An initial guess for the design space for the material parameters is selected.
(2)
Numerical simulations are performed with the different sets of parameters, representing the experimental design points needed for filling the design space.
(3)
For each simulation, the predicted results are compared with the experimental ones, and the cost function values are calculated according to (4).
(4)
A response surface is constructed to approximate the values of the cost function; typically, least-squares approximations are used to determine second-order polynomials.
(5)
An optimization algorithm is applied to determine the minimum point of the response surface (i.e., whereF
(
A
) is minimum), providing the optimal set of material parameters.If a converged solution is not found, the process starts all over again, adding a new region of interest to the design space.


## 3.1. Cost Function

A wide number of cost function formulations for the identification of constitutive parameters have been proposed in literature (e.g., [57, 58]). According to Cao and Lin [57], the cost function should operate as an “efficient guide” of the optimization procedure, in order to search for the best fit to the experimental results; therefore, the ideal cost function should comprise the following conditions:(i)
All measured points of a given experiment should be part of the optimization procedure and have equal opportunity to be optimized, provided that experimental errors are eliminated.
(ii)
All experiments should have equal opportunity to be equally optimized, and so the optimization should not depend on the number of points considered in each experiment.
(iii)
Different units of measure in the cost function should not affect the performance of the optimization.
(iv)
The identification procedure should not be dependent of the user, and so the values of the weighting factors should be optimized to achieve the abovementioned conditions.Cost functions are typically formulated under the concept of weighted least-squares, as follows: (4)
F
A
=
1
m
1
n
∑
i
=
1
m
w
i
∑
j
=
1
n
w
j
r
i
j
A
2
,where F
(
A
) is the cost function to minimise; A is the vector of constitutive parameters to optimize; m is the total number of experiments and n is the total number of points, considered in each experiment i; r
i
j
(
A
) is the residual between the numerically predicted results and those of the experiment i at point j; w
i and w
j are the weighting factors for each experiment i and for each point j, respectively. Within the context of inverse parameter identification, r
i
j
(
A
) can contain variables such as loads, pressures, angular moments, or those arising from full-field measurements (displacements or strains), as will be seen in detail later.The residuals can be expressed in terms of relative differences,(5)
r
i
j
A
=
u
i
j
N
u
m
A
-
u
i
j
E
x
p
u
i
j
E
x
p
,
w
i
t
h
i
=
1
,
…
,
m
,
j
=
1
,
…
,
n
,or in terms of absolute differences, (6)
r
i
j
A
=
u
i
j
N
u
m
A
-
u
i
j
E
x
p
,where u
i
j
N
u
m and u
i
j
E
x
p are, respectively, the numerically predicted and the experimental results at point j of experiment i. Residuals are often expressed using relative differences, which allows the use, in the same cost function, of several kinds of quantities exhibiting various orders of magnitude and units of measure [59]. When u
i
j
E
x
p admits values close to or equal to zero, the residuals should be expressed using absolute differences.


## 3.2. Optimization Algorithms

The minimization of the least-squares cost function, presented in (4), requires efficient and robust optimization algorithms, due to the strongly nonlinear nature of the least-squares cost function [60]. For this purpose, several optimization algorithms are described in the literature, which are commonly divided into two categories: gradient-free algorithms and gradient-based algorithms. Hybrid optimization strategies using both gradient-free and gradient-based algorithms are also proposed (e.g., [15, 60]).Gradient-free algorithms, such as evolutionary and SIMPLEX algorithms, have a great probability of achieving a global minimum due to their random search capability. They require a large number of cost function evaluations (i.e., iterations) and therefore the convergence can be very time-consuming. Because of this, gradient-free algorithms are not recommended within the context of inverse identification strategies, since they require a large number of finite element simulations and analyses [61].Gradient-based algorithms are most popular within inverse identification strategies, as they require far less cost function evaluations than gradient-free algorithms. As local optimizers, these algorithms use information of the gradient to update the vector of constitutive parameters in an adequate search direction [62]. Therefore, there is no guarantee that these algorithms converge to the global minimum, with the possibility of converging to undesirable local minima. This makes the optimization procedure dependent on the initial estimate for the parameters, and therefore the choice of convenient initial estimates for the constitutive parameters can be essential.Examples of gradient-based algorithms commonly used within the context of inverse identification strategies are the Gauss-Newton and Levenberg-Marquardt algorithms. The Gauss-Newton algorithm is described as follows:(7)
A
s
+
1
=
A
s
-
J
T
W
J
-
1
J
T
W
r
Α
s
,where s is the iteration step, A is the vector of constitutive parameters, W is the vector of weighting factors, J is the Jacobian matrix that expresses the sensitivity of the computed results to the constitutive parameters, and r
(
Α
s
) is the vector of residuals, which can be expressed in terms of relative or absolute differences (see (5) and (6), resp.). The dimension of the vector of residuals depends on the total number of experiments m and the total number of points n, in each experiment, that is, the dimension n
r
=
m
×
n. Considering that the total number of constitutive parameters to be identified is n
p, with n
r
≥
n
p, the Jacobian containing the partial derivatives of the residuals with respect to the constitutive parameters is defined: (8)
J
l
,
p
=
∂
r
l
A
p
∂
A
p
,
w
i
t
h
l
=
1
,
…
,
n
r
,
p
=
1
,
…
,
n
p
.An efficient method to compute the Jacobian matrix is finite differentiation. In order to improve convergence, the Jacobian matrix must be updated at each iteration step s. However, the calculation in each step requires high computational cost (at least one numerical simulation per constitutive parameter). To overcome this inconvenient, Endelt et al. [63, 64] and Cooreman [62] highlighted the possibility of computing the sensitivity matrix analytically, with the latter author having concluded the inability of this approach for computing the sensitivities of strain fields to the material parameters, in mechanical tests involving complex and/or heterogeneous deformation.In some cases, the Gauss-Newton algorithm can become unstable in the neighbourhood of the minimum, and so a stabilisation procedure is required. The Levenberg-Marquardt algorithm is similar to Gauss-Newton one, but includes a stabilising term, as follows [65]: (9)
A
s
+
1
=
A
s
-
J
T
W
J
+
λ
d
i
a
g
J
T
W
J
-
1
J
T
W
r
Α
s
,where λ is the stabilising parameter that is updated in each iteration according to the convergence rate [65]. When the Levenberg-Marquardt method shows stability, small values for λ are recommended for fast convergence; otherwise, large values of λ usually allow stable convergence, although slower, towards the minimum. Note that for λ
=
0 the Levenberg-Marquardt algorithm is equal to the Gauss-Newton one.A different type of optimization technique that has been recently used in the identification of constitutive parameters is the Response Surface Methodology (RSM) (e.g., [54, 66]). RSM is an optimization technique for generating smooth approximations of complex functions in a multidimensional design space. In the context of parameter identification, the design space contains all possible combinations for the constitutive parameters and related values of the cost function. The prohibitive size of the full design space requires a Design of Experiments (DoE), to efficiently construct an approximated design space from a few number of representative points (i.e., sets of constitutive parameters). The responses of the representative points (i.e., the values of the cost function) are used to fit a response surface, which is typically obtained from second-order polynomial regression, for the sake of simplicity. Finally, the minimum of the response surface is calculated using a gradient-based optimization procedure, which leads to an estimate of the optimal set of constitutive parameters. In brief, the RSM technique can be summarised as follows:(1)
An initial guess for the design space for the material parameters is selected.
(2)
Numerical simulations are performed with the different sets of parameters, representing the experimental design points needed for filling the design space.
(3)
For each simulation, the predicted results are compared with the experimental ones, and the cost function values are calculated according to (4).
(4)
A response surface is constructed to approximate the values of the cost function; typically, least-squares approximations are used to determine second-order polynomials.
(5)
An optimization algorithm is applied to determine the minimum point of the response surface (i.e., whereF
(
A
) is minimum), providing the optimal set of material parameters.If a converged solution is not found, the process starts all over again, adding a new region of interest to the design space.


## 4. Inverse Identification Strategies

While classical strategies make use of global measurements from experiments to infer the values of the constitutive parameters, using simple analytical relations to estimate the material response under the assumption of homogeneous stress and strain fields in the region of interest, the inverse identification strategies are much more flexible [20]. They make use of experiments allowing heterogeneous deformation and/or strain path changes, as close as possible of the conditions usually found during real sheet metal forming processes. In this perspective, some authors even proposed tests involving contact with friction, such as the punch stretch test (e.g., [67]) and the cylindrical cup test (e.g., [68]), for performing the inverse parameter identification. In these latter cases, the adequate description of the local contact with friction is of paramount importance because it can affect the final results of the parameter identification (e.g., see [69]).The inverse identification strategies make use of global measurements, such as tool loads and tool displacements, which are usually coupled with local measurements, represented as full-field states of displacements and/or strains on the surface of the sample. Then, a numerical analysis of the mechanical test is performed, assuming a constitutive model chosena priori and an initial estimate for its parameters. Finally, the experimental results of the mechanical test are iteratively compared with numerical results by acting on the values of the constitutive parameters until there is an adequate agreement between experimental and numerical results. Figure 2 shows a schematic representation of inverse identification strategies based on the comparison between experimental and numerical measurements using FE simulations.Figure 2
Schematic representation of inverse identification strategies. Adapted from [31].The advantages of this identification approaches include substantial amount of reliable data extracted from a single mechanical test, using full-field measurements, which enables the accurate identification of large sets of constitutive parameters taking into account a wide range of strain levels and strain paths; therefore, it does not require uniform stress and strain distributions, in the region of interest, and no particular restrictions to the sample geometry and/or loading conditions are imposed.Nevertheless, due to the design of the sample geometry, loading conditions, and induced strain paths, the inverse identification requires proper computational strategies [20]. The most common strategy uses Finite Element Model Updating (FEMU) and consists of performing successive finite element (FE) simulations of the physical experiment; the set of parameters are obtained by minimising the difference between the experimental and the numerical measurements. This difference is expressed by a cost function and its minimisation is performed using optimization algorithms, which automatically operate on the values of the constitutive parameters. Usually, cost functions compare the experimental and simulated loads and full-field measurements (e.g., [23, 59, 70–72]); less frequently, some authors propose to use only the load (e.g., [54, 66]), or full-field displacements (e.g., [24]) or strains (e.g., [21, 27–29]), at a given moment of loading.A promising alternative to the use of FEMU is the Virtual Fields Method (VFM), which is based on the principle of virtual work. This approach does not require using time-consuming FE analysis and therefore avoids potential drawbacks related to the accuracy of FE models, namely, the representation of the geometry and boundary conditions [73]. The VFM was successfully used in the identification of parameters of constitutive laws describing the plastic behaviour of metal sheets [73–76]. However, the accuracy of the parameter identification using VFM depends on the adequate choice of the virtual field, which is currently a challenge for problems involving large heterogeneity of deformation of anisotropic materials, as well as large plastic deformations. In fact, in this type of problems, the optimal virtual field has to be evaluated for each time increment, which makes it less attractive than that for linear problems.


### 4.1. Overview of FEMU Strategies

This subsection provides an overview of the literature on inverse methodologies for the identification of parameters of constitutive laws based on inverse strategies coupled with FE simulations. These cases highlight that the identification procedure is dictated by the loading conditions, the geometry of the sample, the type of experimental results selected for the analysis, the cost functions, and the optimization algorithm used. In this context, inverse strategies were developed for the identification of the anisotropy and hardening behaviour of metal sheets, simultaneously or separately. Table1 shows a comparative outline of these strategies, whose key details are highlighted in the following subsections.Table 1
Comparative overview of inverse identification strategies coupled with FE simulations.
Author
Type of test
Strain paths
Strain path reversal
Constitutive model
Optimization steps
Number of cost functions
Analysed results
Measuring instants
Cooreman et al. [21]
Cruciform biaxial tensile test
Uniaxial tension to biaxial stretching
No
Hill’48 criterion+Swift isotropic law
One
One
Strain field (2D)
Various(7)
Schmaltz and Willner [24]
Cruciform biaxial tensile test
Simple shear to biaxial stretching (but depending on the test geometry)
No
Hill’48 criterion+Hockett-Sherby isotropic law(only 2 of 4 parameters identified)
Two steps usingthe same cost function
One
Displacement field (2D)
One
Prates et al. [25]
Cruciform biaxial tensile test
Uniaxial tension to biaxial stretching
No
Hill’48 criterion+Swift isotropic law
Five steps usingfour cost functions
Four
Load+equivalent strain field (1D)
During the testOne
Prates et al. [26]
Cruciform biaxial tensile test
Uniaxial tension to biaxial stretching
No
Various criteria (4)+Swift and Voce isotropic laws
Two (Hill’48 criterion) or three (other criteria)
Two or three
Load+equivalent strain+strain path fields (1D)
During the testOne
Zhang et al. [28]
Cruciform biaxial tensile test
Uniaxial tension to biaxial stretching
No
Bron and Besson criterion
One
One
Strain field (1D)
One
Liu et al. [29]
Cruciform biaxial tensile test
Uniaxial tension to biaxial stretching
No
Voce-type isotropic law
One
One
Strain field (2D)
During the test
Güner et al. [22]
Uniaxial tensile test on a notched sample
Near uniaxial tension
No
Yld2000 (2D) criterion
One
One
Strain field (2D)+load
During the test
Pottier et al. [23]
Out-of-plane test
Simple shear to biaxial stretching
No
Hill’48 Criterion+Ludwick isotropic law
One
One
Load+displacement field (3D)
Various(6)
Chamekh et al. [51]
Circular and elliptical bulge tests
Biaxial stretching
No
Ludwick isotropic law+
r
0, r
45, and r
90
Two
One
Pressure
During the test
Bambach [52]
Circular bulge test
Biaxial stretching
No
Voce isotropic law
Two
Four
Pressure, pole strain, and pole thickness
During the test
Reis et al. [53]
Circular bulge test
Biaxial stretching
No
Swift isotropic law
Four
One
Pressure
During the test
Eggertsen and Mattiasson [54]
Three-point bending test
Tension-compression
Yes
Various kinematic laws (4)
One
One
Load
During the test
Yin et al. [55]
Twin bridge shear test
Pure shear to uniaxial tension
Yes
Armstrong-Frederick kinematic law+Voce isotropic law (only 2 of 3 parameters identified)
One
One
Angular moment
During the test
Pereira et al. [56]
Shear test
Shear with reversal
Yes
Lemaître-Chaboche kinematic law+Swift isotropic law
One
One
Load
During the test
Note: 1D: in one direction in the sheet plane; 2D: on the surface of the sheet plane; 3D: three-dimensional.


#### 4.1.1. Identification of Isotropic Hardening and/or Yield Criterion Parameters

(1) Strategies Using Cruciform Specimens. There has been a steadily growing interest in developing inverse identification strategies supported by the use of the biaxial tensile test of cruciform specimens, coupled with full-field displacement or strain measurements (e.g., [21, 24–29]). In general, this test allows (i) strain paths ranging from uniaxial tension (in the arms region of the specimen) to balanced biaxial tension (in the centre section of the specimen), (ii) high strain gradients from the centre region of the specimen to the end of the arms region, and (iii) no contact between surfaces and therefore no friction. Also, by changing the load and/or displacement ratio over the two normal loading axes, it is possible to achieve several biaxial stress states in the central region of the specimen. However, this kind of test only permits attaining low values of equivalent plastic strain (close to those obtained in uniaxial tension) before instability occurs and no occurrence of out-of-plane shear stresses is observed (i.e., the test is insensitive to the material parameters associated with these stresses). Figure 3 shows a set of cruciform geometries proposed in the literature for identifying material parameters using FEMU strategies.Figure 3
Cruciform specimen geometries proposed in the literature for FEMU strategies: (a) Schmaltz and Willner [24]; (b) Cooreman et al. [21]; (c) Prates et al. [25, 26]; (d) Zhang et al. [27, 28]; and (e) Liu et al. [29].
(a)
(b)
(c)
(d)
(e)In this context, Cooreman et al. [21] proposed the use of the biaxial tensile test on a perforated cruciform specimen (see Figure 3(b)), to simultaneously identify the material parameters of Swift hardening law [34] and Hill’48 yield criterion that describe the plastic behaviour of a 0.8 mm thick DC06 sheet steel. The identification makes use of strain field data measured at the central region of the specimen (see dashed area in Figure 3(b)) taken at 7 distinct load steps and iteratively compared with their numerical counterparts using Gauss-Newton algorithm, using the following cost function: (10)
F
A
=
1
3
n
∑
i
=
1
n
ε
x
x
E
x
p
i
-
ε
x
x
N
u
m
A
i
2
+
ε
y
y
E
x
p
i
-
ε
y
y
N
u
m
A
i
2
+
ε
x
y
E
x
p
i
-
ε
x
y
N
u
m
A
i
2
1
/
2
,where ε
E
x
p
i and ε
N
u
m
A
i are the experimentally determined and numerical values of the strain components ε
x
x, ε
y
y, and ε
x
y at point i, respectively, and n is the total number of measuring points. According to the authors, the results from the inverse strategy are similar to those from classical strategies, except for the parameter ε
0 of the Swift law, which leads to clearly different yield stress values. The authors attribute the discrepancy of ε
0 results to the use of strain field results from loading steps that neglect the beginning of the test. They suggest performing the inverse identification using additional strain fields from loading steps located near the onset of plastic deformation, which in our opinion can lead to high relative errors. A simpler alternative, from the experimental point of view, would be to include in the identification a cost function considering the load-displacement curves for both axes of the cruciform sample, as adopted in other identification strategies later described.Schmaltz and Willner [24] explored the usability of the biaxial tensile test on three cruciform specimen geometries, promoting different types of heterogeneous strain fields (see Figure 3(a)), in order to identify the plastic behaviour of a 2.0 mm thick DC04 sheet steel modelled via Hill’48 yield criterion and Hockett-Sherby hardening law [39]. The experimental and numerical displacement fields, in the regions in red (see Figure 3(a)) were compared and their difference was minimised using the Levenberg-Marquardt algorithm, with the following cost function:(11)
F
A
=
∑
l
=
1
2
∑
i
=
1
n
u
l
E
x
p
i
-
u
l
N
u
m
A
i
2
,where u
l
E
x
p
i and u
l
N
u
m
A
i are the experimentally determined and the numerically predicted values of the displacements at point i, respectively, in the 0
x and 0
y directions of the sheet plane (l
=
1
,
2), at a given moment of the test.The optimization procedure starts with three different initial sets of parameters that lead to quite similar optimized values, for each test geometry, suggesting that the global optimum is reached. Moreover, the identification is split in two sequential steps (using the same cost function): (i) the first step identifies two of the Hockett-Sherby hardening law parameters (the yield stress and the saturation stress), under the assumption of isotropic (von Mises) material; the two remaining parameters of this law (contained in the exponent) are obtained by fitting the experimental results determined from the biaxial tensile test and are kept fixed during the identification procedure; (ii) in the second step, the Hill’48 parameters are identified. The cruciform geometry with the central hole (in the middle of Figure3(a)) allows achieving the best results in terms of convergence of the iterative procedure and accuracy of the identified material parameters, which was assigned to the strong heterogeneity of the kinematic field.Prates et al. [25, 26] designed a cruciform sample and developed two strategies for simultaneously identifying the parameters of the anisotropic yield criteria and isotropic hardening law of sheet metals. Both strategies use the results of the load evolution during the test and of the major and minor principal strains distributions, along the axes of the sample, at a given moment of the test, preceding the maximum load. Both strategies were numerically tested. The optimization of the design of the cruciform sample (see Figure 3(c)) was performed by means of a numerical study with the purpose of maximising the sensitivity of the test results to the values of the constitutive parameters and for allowing a wide range of strain paths, from uniaxial tension, at the arms of the sample, to near equibiaxial tension, at the centre of the sample.The work was initially addressed for the identification of the parameters of the Hill’48 yield criterion and the isotropic Swift hardening law [25]. An inverse identification was performed without resorting to the traditional optimization algorithms (e.g., gradient-based algorithms, or others); that is, a specific algorithm was built for this purpose. The inverse analysis algorithm consists of a sequence of five optimization steps, using the results of load evolution with the sample boundaries displacement during the test, for the axes 0
x and 0
y (see Figure 3(c)) and the distributions of von Mises equivalent strain along both axes of the sample, for an instant preceding and close to the maximum load. At each step, one or more parameters or a relationship between them is identified. The strategy was tested using numerically generated results of fictitious materials, which proved to be competitive, when compared with classical strategies. This allowed understanding that a sequential optimization, since properly elaborated, is clearly advantageous when compared to most commonly inverse identifications, consisting of using a unique cost function including different types of results.The inverse analysis strategy mentioned above [25] enabled a good understanding of the issues involved, namely, concerning the delineation of the sequential algorithm leading to upper accuracy. This allowed extending the strategy to more complex constitutive models (yield criteria and isotropic hardening laws). Therefore, a general inverse identification strategy that sequentially uses three distinct cost functions was developed [26]. It resorts to the Levenberg-Marquardt algorithm for sequential optimization of the parameters of the yield criteria and isotropic hardening laws. More importantly, this strategy allows the identification of parameters of several yield criteria and hardening laws. It can be used directly for a given criterion or, sequentially, starting from the Hill’48 yield criterion and then using the Hill’48 criterion solution as an initial estimate for identifying the parameters of other criteria, on the condition that can be converted into the Hill’48 yield criterion for particular values of the parameters. In the last case, this strategy is detached in two stages and has the advantage of enabling the assessment of the adequacy of a number of constitutive models to describe the experimental results, starting from a simple anisotropic criterion, the Hill’48. The first stage consists of the simultaneous identification of the hardening law, Swift and/or Voce [38], in the case, and Hill’48 yield parameters, using the results of the load evolution in function of the displacement of the grips and the equivalent strain distribution at a given moment of the test, along the axes of the sample. The hardening parameters must be separately identified for the Swift and Voce laws and the one (Swift or Voce) that better describes the results of the cruciform test (if it is possible to distinguish) is selected for further optimization. The first stage involves the sequential minimisation of the following cost functions: (12)
F
1
A
=
1
Q
1
∑
i
=
1
Q
1
P
i
N
u
m
A
-
P
i
E
x
p
P
i
E
x
p
0
x
2
+
1
Q
2
∑
i
=
1
Q
2
P
i
N
u
m
A
-
P
i
E
x
p
P
i
E
x
p
0
y
2
,
F
2
B
=
1
R
1
∑
i
=
1
R
1
ε
-
i
N
u
m
B
-
ε
-
i
E
x
p
ε
-
i
E
x
p
0
x
2
+
1
R
2
∑
i
=
1
R
2
ε
-
i
N
u
m
B
-
ε
-
i
E
x
p
ε
-
i
E
x
p
0
y
2
,where P
i
N
u
m
(
A
) and P
i
E
x
p are the experimentally determined and numerical values of the load, respectively, and ε
-
i
N
u
m
(
B
) and ε
-
i
E
x
p are the experimentally determined and numerical values of the equivalent strain, respectively; Q
1, Q
2, R
1, and R
2 are the total number of measuring points for axes 0
x and 0
y of the sample; A and B are vectors of isotropic hardening law and yield criteria parameters, respectively.The second stage allows extending the parameters identification procedure to more complex yield functions, such as Barlat’91 [8], Karafillis and Boyce [9], and Drucker+L [3], the cases studied in the work. This second stage should be performed whenever the identification carried out during the first stage is found to be not enough satisfactory to capture the experimental strain paths results, along the axes of the sample, which are not considered in minimisation during the first stage. The second stage of this inverse strategy involves the minimisation of the following cost function: (13)
F
3
C
=
1
S
1
∑
i
=
1
S
1
ρ
i
N
u
m
C
-
ρ
i
E
x
p
0
x
2
+
1
S
2
∑
i
=
1
S
2
ρ
i
N
u
m
C
-
ρ
i
E
x
p
0
y
2
,where ρ
i
N
u
m
(
C
) and ρ
i
E
x
p are the experimentally determined and numerical values of the strain path, respectively, S
1 and S
2 are the total number of measuring points for axes 0
x and 0
y of the sample, and C is the vector of yield criteria parameters. This sequential optimization procedure is a successful alternative to the parameter identification by minimising a single cost function comprising all material parameters and results of different types, as commonly found in the literature. Namely, it is concluded that this last approach can deteriorate the description of the material behaviour, concerning the load versus displacement results, and therefore the parameters and the choice of the hardening law, without apparent improvement of the description of the results of equivalent strain and strain path distributions.Zhang et al. [28] identified the parameters of Bron and Besson yield criterion [4] for both AA5086 aluminium alloy and DP980 dual-phase steel sheets, using a cruciform sample previously designed by the authors [27] and shown in Figure 3(d). The inverse identification strategy consists of minimising the gap between the experimental and numerical distributions of the major and minor strains along the diagonal direction of the sample central area, at an instant immediately before rupture, using a SIMPLEX optimization algorithm. The cost function used is defined as follows:(14)
F
A
=
∑
l
=
1
2
∑
i
=
1
n
ε
l
E
x
p
i
-
ε
l
N
u
m
A
i
2
∑
i
=
1
n
ε
l
E
x
p
i
2
,where (
ε
l
E
x
p
)
i and (
ε
l
N
u
m
(
A
)
)
i are the experimentally determined and numerical values of the principal strain components (l
=
1
,
2) at point i, respectively, n is the total number of measuring points along the diagonal path, and A is the vector of 13 parameters to be identified: four isotropy and eight anisotropy parameters of the Bron and Besson yield model and the yield stress value. The hardening of the material is modelled with isotropic hardening described by the Voce law, in case of AA5086 aluminium, and an equation based on Swift and Voce laws, for DP980 steel. In both cases, the hardening parameters are directly identified from results of the tensile test in the rolling direction and were kept fixed during the inverse identification procedure (except the yield stress).Liu et al. [29] designed a cruciform sample with a thickness-reduced central zone and four slots at each arm (see Figure 3(e)), to perform the parameters identification of a modified Voce law, describing the hardening behaviour of AA5086 aluminium sheet. The inverse analysis procedure makes use of the SIMPLEX algorithm, to minimise the difference between experimental and numerical strains measured at the centre region of the sample during the test, expressed by the following cost function:(15)
F
A
=
∑
l
=
1
2
∑
i
=
1
n
ε
l
E
x
p
i
-
ε
l
N
u
m
A
i
2
∑
i
=
1
n
ε
l
E
x
p
i
2
1
/
2
,where (
ε
l
E
x
p
)
i and (
ε
l
N
u
m
(
A
)
)
i are the experimentally determined and numerical values of the principal strain components (l
=
1
,
2) at point i, respectively, at the central point of the sample, n is the total number of the time points of simulation, and A is the vector of three hardening parameters of the modified Voce law. The experimental force evolution along the two arms of the cruciform specimen was applied to the FE model for the numerical simulations, taking into account the lack of synchronization observed in the two tensile forces on each axis of the sample. In this inverse identification of the modified Voce law parameters, three yield functions were considered: von Mises, Hill’48, and the more advanced Bron and Besson criterion, whose parameters were previously identified [27]. The identified biaxial flow stress curves were then compared with an experimental curve obtained from a uniaxial tensile test, showing that a good agreement can be achieved if an adequate yield function is used in the FE model.(2) Strategies Using Bulge Test. Most of the strain paths observed in deep-drawing components are in the range between simple tension and balanced biaxial tension, which justifies the widespread use of the cruciform specimen in the framework of the methodologies of inverse analysis. Nevertheless, this test has a strong drawback related to low deformation levels achieved, particularly in the central region of the sample, in which the strain and stress paths can be close to equibiaxial, although being dependent on the applied load or displacement ratios along the two axes of the specimen. In contrast, the bulge test allows obtaining relatively high strain values before necking, and so the flow stress curves can be assessed up to large strain values, for several biaxial strain (or stress) paths depending on the geometry of the die (circular or elliptical). In this context, a few cases of inverse analysis methodologies were developed, with the aim of identifying the parameters of work-hardening laws.Chamekh et al. [51] describe an inverse approach, based on Artificial Neural Networks (ANN), to identify the material parameters of a stainless steel (AISI 304). They use the results of the evolutions of pressure with the pole height, which are transferred to the neural network. The ANN is trained by means of curves of pressure versus displacement of the central point of the cap, generated by finite element simulations of the circular bulge test. During the training process, the network computes the weight connections, minimising the total mean squared error between the actual output and the desired output. The neural network generates an approximated function for the material parameters depending on the profile of the evolution of pressure with the pole height curve. Then, it was exploited for the identification of material parameters from experimental results. The Ludwick hardening law [33] and the Hill’48 yield criterion were selected. Therefore, the set of parameters to be identified also comprise the Lankford coefficients, r
0, r
45, and r
90. The material parameters are identified according to the following two steps: (i) the first step, using the circular bulge test, is to find the Ludwick hardening law parameters (assuming the knowledge of Lankford’s coefficients determined from the tensile tests); (ii) the second step, using the elliptical bulge test for an off axis angle of 0°, is to recalculate the Lankford’s coefficients. An elliptical die for an off axis angle of 45° is used for the validation of the parameters identification. The authors conclude the following: (i) the ANN methodology can predict acceptable combination of the values of the material parameters; (ii) once the ANN was trained, output results for a given set of input data are available almost instantaneously. Despite these conclusions, it should be noted that the values of the experimental (from the tensile test) and identified hardening coefficients are far away (the experimental and identified values are 0.67 and 0.4, resp.). The remaining identified values of the parameters differ between 20 and 30% when compared with the tensile test except for the yield stress, whose values are approximately equal.Bambach [52] explored the usability of the circular bulge test to identify the parameters of the Voce law of a fictitious material, which is considered isotropic. Initially, the membrane theory is applied to the results as in experimental cases, in order to obtain a set of parameters of the Voce, by fitting the stress versus strain results. The inverse analysis strategy proposed resorts to a gradient-based optimization algorithm, which is known for being sensitive to the initial solution. Thus, by using an initial solution, the one previously obtained with the membrane theory, it is expected to avoid convergence problems. The work gives special focus on the choice of the cost function to be minimised, making use, separately or simultaneously, of results of pressure versus pole height, pole strain versus pole height, and pole thickness versus pole height, and formulated as follows:(16)
F
p
A
=
∫
0
h
m
a
x
p
E
x
p
h
-
p
N
u
m
h
,
A
2
d
h
,
F
ε
A
=
∫
0
h
m
a
x
ε
E
x
p
h
-
ε
N
u
m
h
,
A
2
d
h
,
F
t
A
=
∫
0
h
m
a
x
t
E
x
p
h
-
t
N
u
m
h
,
A
2
d
h
,where F
p
(
A
), F
ε
(
A
), and F
t
(
A
) are cost functions defined by the pressure, p, strain, ε, and thickness, t, with pole height, h, respectively. The author concluded that the reidentification procedure (so called by the author) is significantly improved when combining the first two types of results. This significantly improves the reidentification, since it will contribute to reducing the search area where the minimum value of the objective function is located. It should be noted that this proposal is a reidentification, which has its starting point in the parameters previously obtained from the use of the membrane theory, such as in the traditional procedure recommended by the ISO standard [77]. Furthermore, it needs to resort to strain results in the pole of the cap during the test, which does not simplify the experimental procedure.Reis et al. [53] proposed an inverse methodology for determining the hardening law of metal sheets, from the results of pressure versus pole height obtained in the bulge test, involving the identification of the parameters of the Swift law. The starting point of this analysis was to realize that it is possible to achieve a unified description (i.e., overlapping) of the evolution of the pressure with the pole height, for a given value of the hardening parameter of the Swift law, regardless of the yield stress and anisotropy of the material and sheet thickness. To achieve the overlapping of such curves, appropriate multiplying factors must be used for the values of pressure and pole height, depending on the yield stresses and thicknesses ratios of the sheets and also on their anisotropy. Thereafter, an inverse analysis methodology was developed, which consists in the search for the best coincidence between pressure versus pole height experimental and reference curves, with the latter being obtained by numerical simulation assuming isotropic material behaviour with various values of the hardening parameter in the range of the material under study. This methodology, when compared with the classical strategy, proves to be an efficient alternative avoiding the use of complex devices for measuring the radius of curvature and strain at the pole of the cap, during the bulge test. Moreover, the authors claim that it is easy to implement and it is more efficient than classical approach, since (i) a unique set of numerical curves can be used within a relatively wide range of hardening coefficients, that is, covering the values usually found within one or several class of materials, without having to remake the simulations every time an identification is performed; (ii) it is not exposed to experimental errors related to the evaluation of the strain at the pole of the bulge and the use of membrane theory approach for assessment of the stress from the radius of curvature, which is usually the major source of errors.(3) Other Specimens. Güner et al. [22] proposed an inverse analysis procedure for the identification of the Yld2000-2D yield criterion parameters [7]. This study uses a notched specimen submitted to a tensile test, in order to obtain an inhomogeneous deformation field. Moreover, a layer compression test is used in order to supply additional information, that is, the equibiaxial yield stress. The required data for the inverse identification are the major and minor principal strains in the sheet plane and the load and the equibiaxial yield stress. The cost function is minimised using the Levenberg-Marquardt algorithm and is written as follows:(17)
F
A
=
∑
i
=
1
n
i
n
c
∑
j
=
1
n
e
l
e
m
ε
1
,
i
j
N
u
m
A
-
ε
1
,
i
j
E
x
p
2
+
ε
2
,
i
j
N
u
m
A
-
ε
2
,
i
j
E
x
p
2
+
C
1
∑
i
=
1
n
i
n
c
P
i
N
u
m
A
-
P
i
E
x
p
2
+
C
2
σ
b
A
-
σ
b
E
x
p
2
,where ε
1 and ε
2 are the principal strains measured at each tool displacement increment i, at each element of the optical measurements, j; P represents the load and σ
b the equibiaxial yield stress; C
1 and C
2 are scale factors. The value of σ
b
(
A
) is analytically calculated at each iteration with the values of A predicted by the optimization algorithm. Different alternative orientations of the specimen with the rolling direction (0°, 45°, and 90°) and configurations of the cost function (setting one of the three terms of the cost function equal to zero) were considered to test the inverse procedure. The authors highlight the importance of including strain information and the equibiaxial stress on the cost function, in order to improve the characterization of the anisotropy coefficients.Pottier et al. [23] developed an out-of-plane testing procedure for the simultaneous identification of Hill’48 yield criterion and Ludwick hardening law parameters of a rolled titanium sheet. Figure 4 illustrates the experimental setup and the geometry of the sample developed by the authors. A hemispherical punch applies a prescribed displacement normal to the sheet plane, at the centre of the surface of the sample, using a simple uniaxial tensile test machine. Two cameras are located on the opposite side of the sample and the components of the displacement fields along the 0
x, 0
y, and 0
z axes are captured during the test using stereo digital image correlation. The sample was designed in order to exhibit multiaxial stress states, including shear, tension, and biaxial stretching. The numerical displacements fields along the 0
x, 0
y, and 0
z axes and the global load are obtained from a numerical model of the test and compared to the experimental ones using a single cost function, formulated as follows:(18)
F
A
=
∑
i
,
j
=
1
N
t
,
n
u
x
,
i
j
N
u
m
A
-
u
x
,
i
j
E
x
p
m
a
x
i
,
j
u
x
E
x
p
2
+
∑
i
,
j
=
1
N
t
,
n
u
y
,
i
j
N
u
m
A
-
u
y
,
i
j
E
x
p
m
a
x
i
,
j
u
y
E
x
p
2
+
∑
i
,
j
=
1
N
t
,
n
u
z
,
i
j
N
u
m
A
-
u
z
,
i
j
E
x
p
m
a
x
i
,
j
u
z
E
x
p
2
+
n
P
i
N
u
m
A
-
P
i
E
x
p
m
a
x
i
P
E
x
p
2
1
/
2
,where N
t is the number of time steps considered and n is the number of measured points; u
x, u
y, and u
z are the displacements along the 0
x, 0
y, and 0
z axes, respectively, P is the global load and A is the vector of parameters to identify; the subscripts Num and Exp refer to numerical and experimental results and i and j refer to the number of time steps and measured points, respectively. The minimisation of the cost function is performed with the Levenberg-Marquardt algorithm. To assess the quality of the identified set of constitutive parameters, the authors performed deep-drawing tests of a circular cup. Moreover, additional identifications of the constitutive parameters of the material were performed, following two different strategies: a classic strategy, based on three tensile tests cut along three different directions in the sheet plane, and an inverse identification strategy using heterogeneous planar shear-like tests, previously proposed [59]. The experimental results of the earing profile of the circular cup were then compared with the numerically predicted results from the different parameter identification strategies. The authors concluded that the use of the nonplanar sample allows a more accurate prediction of the earing profile than the planar shear-like tests and the three tensile tests, since the nonplanar sample test covers a wider range of strain paths. However, in our opinion, the chosen mechanical test shows a relative degree of complexity and has the inconvenience of presenting contact with friction between the punch and the metal sheet, which always raises questions concerning the impact of contact with friction on the identification of the constitutive parameters.Figure 4
Schematic representation of the out-of-plane test developed by Pottier et al. [23].


#### 4.1.2. Identification of Kinematic Hardening

The identification of the kinematic hardening plays a significant role when phenomena such as the Bauschinger effect and permanent softening, due to reverse of strain path and other strain path changes, are relevant for the subsequent plastic deformation behaviour. Depending on the material and the strain path changes occurring in the deep-drawing process, these phenomena can be more or less noticeable and relevant. The Bauschinger effect is also associated with the springback, due to premature yielding after reversal of strain path [14]. Springback is the strain recovery when forming loads are removed, and its magnitude, which depends on the flow stress value [11] and Young’s modulus, can also be influenced by the Bauschinger effect [78]. Therefore, the proper modelling and identification of the kinematic hardening are also important in order to efficiently predict the springback.In this sense, Eggertsen and Mattiasson [54] were focused on the identification of the kinematic hardening law with the main concern of the accurate prediction of springback. The goal was to select the model able to accurately describe kinematic hardening features, such as the early reyielding, transient behaviour, work-hardening stagnation, and permanent softening, taking also into account the complexity on the evaluation of its parameters. Three-point bending tests (see Figure 5) on four typical materials from car manufacturing industry were performed: two dual-phase steels (TKS-DP600HF and SSAB-DP600), from different suppliers and with different thicknesses, a mild steel (Voest-DX56D), and an interstitial-free steel (TKS-220IF). Five different hardening models were considered: (i) a pure isotropic hardening law, used as comparative reference; (ii) a mixed isotropic-kinematic law [79, 80]; (iii) the Armstrong-Frederick model [44]; (iv) the Geng-Wagoner hardening law [47]; and (v) the Yoshida-Uemori hardening law [14]. The hardening parameters of all models were determined by inverse analysis, where the difference between the experimentally and numerically generated load-displacement curves of the three-point bending test is minimised. The inverse identification of the material parameters was performed resorting to Response Surface Methodology, using the following cost function:(19)
F
A
=
1
n
∑
i
=
1
n
w
i
P
i
N
u
m
A
-
P
i
E
x
p
s
i
2
,where P
i
N
u
m
(
A
) and P
i
E
x
p represent the calculated and measured values of the punch load as a function of the vector of hardening parameters A, respectively; s
i is the residual scale factor; and w
i is the weight applied to each component of the cost function. Both s
i and w
i were set equal to 1. The authors conclude that Yoshida-Uemori model provides the best result for all materials, while the isotropic hardening model gives the worst result. However, taking into account the accuracy and the complexity of the hardening model, the authors point out that Geng-Wagoner law corresponds to a better compromise. In fact, they state that about 30 simulations are needed to optimize the parameters of the mixed isotropic-kinematic hardening law, while up to 170 simulations are required to optimize the parameters of Yoshida-Uemori hardening law.Figure 5
Device for three-point bending test [54]: (a) experimental setup and (b) schematic representation of horizontal and vertical views.
(a)
(b)Pereira et al. [56] outline an inverse analysis methodology for simultaneously identifying the parameters of the isotropic Swift law and the Lemaître-Chaboche kinematic law [13] of metal sheets, using a reverse shear test. The outlined strategy uses a modified shear sample with a cylindrical notch along the axis of the sample, in order to confine the plastic deformation within the entire gauge section, which is not the case for the classical shear samples with constant thickness (see Figure 6). The geometry of the cylindrical notch was defined in order to ensure that all strain values, between the maximum (in the centre of the notch) and the minimum (zero, along the edge of the notch), are as much as possible represented at the moment of the strain path reversing. The geometry of the sample allows that the boundary conditions of experimental tests are accurately reproduced numerically and avoids the errors in the experimental determination of the stress versus strain curves, used in traditional methodologies, whose accuracy requires homogeneity of the stress and strain fields in the sample. The inverse analysis methodology consists of minimising the gap between experimental and numerical load versus displacement curves by making variations of the constitutive parameters, using the Levenberg-Marquardt algorithm. The following cost function is used:(20)
F
A
=
1
q
∑
i
=
1
q
P
i
E
x
p
-
P
i
N
u
m
A
P
i
E
x
p
2
+
1
p
∑
i
=
1
p
P
i
E
x
p
-
P
i
N
u
m
A
P
i
E
x
p
2
,where P
i
E
x
pand P
i
N
u
m
(
A
) are, respectively, the experimental and numerical values of load for the same tool displacement; q and p are the total number of points in the forward and reverse paths, respectively; and A is the vector of the parameters to be identified. The parameters of the Hill’48 yield criterion that best describe the anisotropy of the fictitious materials used in this work (described by Drucker+L yield criterion) were identified following the methodology mentioned above, proposed by Prates et al. [25]. This methodology also allows identifying the parameters of the isotropic Swift law, which were used as first estimate in this inverse analysis. If no identification of the Swift law parameters was previously performed, the first estimate of isotropic Swift hardening parameters can be obtained adopting typical values for the material under study. It is appropriate to experimentally test this methodology.Figure 6
Geometry and dimensions of the classical and proposed notched shear specimens used in [56].Yin et al. [55] proposed the use of the twin bridge cyclic shear test proposed by Brosius et al. [81] to evaluate the Bauschinger effect on three classes of steel sheets, DC06, DP600, and TRIP700. The twin bridge shear specimen has two gauge areas that are simultaneously deformed when a moment is applied. Due to the moment application, instead of load, no unwanted reaction moment is created, in contrast with the ASTM shear specimen. The inverse identification strategy involves minimising the difference between experimental and numerical angular moment results, resorting to a Trust Region Reflective Method [82], to identify the parameters of isotropic Voce hardening law and Armstrong-Frederick kinematic hardening law. The cost function is described as follows:(21)
F
A
=
∑
i
=
1
q
M
i
E
x
p
-
M
i
N
u
m
A
2
+
∑
i
=
1
p
M
i
E
x
p
-
M
i
N
u
m
A
2
,where M
i
E
x
p and M
i
N
u
m
(
A
) are, respectively, the experimental and numerical values of moment at the same rotation angles; q and p are the total number of points in the forward and backward paths, respectively; and A is the vector of kinematic hardening parameters of the Lemaître-Chaboche law and two of the three parameters of the Voce law. The remaining Voce law parameter, the yield stress, is obtained from uniaxial tensile tests in the rolling direction and is kept fixed during the identification procedure.


## 4.1. Overview of FEMU Strategies

This subsection provides an overview of the literature on inverse methodologies for the identification of parameters of constitutive laws based on inverse strategies coupled with FE simulations. These cases highlight that the identification procedure is dictated by the loading conditions, the geometry of the sample, the type of experimental results selected for the analysis, the cost functions, and the optimization algorithm used. In this context, inverse strategies were developed for the identification of the anisotropy and hardening behaviour of metal sheets, simultaneously or separately. Table1 shows a comparative outline of these strategies, whose key details are highlighted in the following subsections.Table 1
Comparative overview of inverse identification strategies coupled with FE simulations.
Author
Type of test
Strain paths
Strain path reversal
Constitutive model
Optimization steps
Number of cost functions
Analysed results
Measuring instants
Cooreman et al. [21]
Cruciform biaxial tensile test
Uniaxial tension to biaxial stretching
No
Hill’48 criterion+Swift isotropic law
One
One
Strain field (2D)
Various(7)
Schmaltz and Willner [24]
Cruciform biaxial tensile test
Simple shear to biaxial stretching (but depending on the test geometry)
No
Hill’48 criterion+Hockett-Sherby isotropic law(only 2 of 4 parameters identified)
Two steps usingthe same cost function
One
Displacement field (2D)
One
Prates et al. [25]
Cruciform biaxial tensile test
Uniaxial tension to biaxial stretching
No
Hill’48 criterion+Swift isotropic law
Five steps usingfour cost functions
Four
Load+equivalent strain field (1D)
During the testOne
Prates et al. [26]
Cruciform biaxial tensile test
Uniaxial tension to biaxial stretching
No
Various criteria (4)+Swift and Voce isotropic laws
Two (Hill’48 criterion) or three (other criteria)
Two or three
Load+equivalent strain+strain path fields (1D)
During the testOne
Zhang et al. [28]
Cruciform biaxial tensile test
Uniaxial tension to biaxial stretching
No
Bron and Besson criterion
One
One
Strain field (1D)
One
Liu et al. [29]
Cruciform biaxial tensile test
Uniaxial tension to biaxial stretching
No
Voce-type isotropic law
One
One
Strain field (2D)
During the test
Güner et al. [22]
Uniaxial tensile test on a notched sample
Near uniaxial tension
No
Yld2000 (2D) criterion
One
One
Strain field (2D)+load
During the test
Pottier et al. [23]
Out-of-plane test
Simple shear to biaxial stretching
No
Hill’48 Criterion+Ludwick isotropic law
One
One
Load+displacement field (3D)
Various(6)
Chamekh et al. [51]
Circular and elliptical bulge tests
Biaxial stretching
No
Ludwick isotropic law+
r
0, r
45, and r
90
Two
One
Pressure
During the test
Bambach [52]
Circular bulge test
Biaxial stretching
No
Voce isotropic law
Two
Four
Pressure, pole strain, and pole thickness
During the test
Reis et al. [53]
Circular bulge test
Biaxial stretching
No
Swift isotropic law
Four
One
Pressure
During the test
Eggertsen and Mattiasson [54]
Three-point bending test
Tension-compression
Yes
Various kinematic laws (4)
One
One
Load
During the test
Yin et al. [55]
Twin bridge shear test
Pure shear to uniaxial tension
Yes
Armstrong-Frederick kinematic law+Voce isotropic law (only 2 of 3 parameters identified)
One
One
Angular moment
During the test
Pereira et al. [56]
Shear test
Shear with reversal
Yes
Lemaître-Chaboche kinematic law+Swift isotropic law
One
One
Load
During the test
Note: 1D: in one direction in the sheet plane; 2D: on the surface of the sheet plane; 3D: three-dimensional.


### 4.1.1. Identification of Isotropic Hardening and/or Yield Criterion Parameters

(1) Strategies Using Cruciform Specimens. There has been a steadily growing interest in developing inverse identification strategies supported by the use of the biaxial tensile test of cruciform specimens, coupled with full-field displacement or strain measurements (e.g., [21, 24–29]). In general, this test allows (i) strain paths ranging from uniaxial tension (in the arms region of the specimen) to balanced biaxial tension (in the centre section of the specimen), (ii) high strain gradients from the centre region of the specimen to the end of the arms region, and (iii) no contact between surfaces and therefore no friction. Also, by changing the load and/or displacement ratio over the two normal loading axes, it is possible to achieve several biaxial stress states in the central region of the specimen. However, this kind of test only permits attaining low values of equivalent plastic strain (close to those obtained in uniaxial tension) before instability occurs and no occurrence of out-of-plane shear stresses is observed (i.e., the test is insensitive to the material parameters associated with these stresses). Figure 3 shows a set of cruciform geometries proposed in the literature for identifying material parameters using FEMU strategies.Figure 3
Cruciform specimen geometries proposed in the literature for FEMU strategies: (a) Schmaltz and Willner [24]; (b) Cooreman et al. [21]; (c) Prates et al. [25, 26]; (d) Zhang et al. [27, 28]; and (e) Liu et al. [29].
(a)
(b)
(c)
(d)
(e)In this context, Cooreman et al. [21] proposed the use of the biaxial tensile test on a perforated cruciform specimen (see Figure 3(b)), to simultaneously identify the material parameters of Swift hardening law [34] and Hill’48 yield criterion that describe the plastic behaviour of a 0.8 mm thick DC06 sheet steel. The identification makes use of strain field data measured at the central region of the specimen (see dashed area in Figure 3(b)) taken at 7 distinct load steps and iteratively compared with their numerical counterparts using Gauss-Newton algorithm, using the following cost function: (10)
F
A
=
1
3
n
∑
i
=
1
n
ε
x
x
E
x
p
i
-
ε
x
x
N
u
m
A
i
2
+
ε
y
y
E
x
p
i
-
ε
y
y
N
u
m
A
i
2
+
ε
x
y
E
x
p
i
-
ε
x
y
N
u
m
A
i
2
1
/
2
,where ε
E
x
p
i and ε
N
u
m
A
i are the experimentally determined and numerical values of the strain components ε
x
x, ε
y
y, and ε
x
y at point i, respectively, and n is the total number of measuring points. According to the authors, the results from the inverse strategy are similar to those from classical strategies, except for the parameter ε
0 of the Swift law, which leads to clearly different yield stress values. The authors attribute the discrepancy of ε
0 results to the use of strain field results from loading steps that neglect the beginning of the test. They suggest performing the inverse identification using additional strain fields from loading steps located near the onset of plastic deformation, which in our opinion can lead to high relative errors. A simpler alternative, from the experimental point of view, would be to include in the identification a cost function considering the load-displacement curves for both axes of the cruciform sample, as adopted in other identification strategies later described.Schmaltz and Willner [24] explored the usability of the biaxial tensile test on three cruciform specimen geometries, promoting different types of heterogeneous strain fields (see Figure 3(a)), in order to identify the plastic behaviour of a 2.0 mm thick DC04 sheet steel modelled via Hill’48 yield criterion and Hockett-Sherby hardening law [39]. The experimental and numerical displacement fields, in the regions in red (see Figure 3(a)) were compared and their difference was minimised using the Levenberg-Marquardt algorithm, with the following cost function:(11)
F
A
=
∑
l
=
1
2
∑
i
=
1
n
u
l
E
x
p
i
-
u
l
N
u
m
A
i
2
,where u
l
E
x
p
i and u
l
N
u
m
A
i are the experimentally determined and the numerically predicted values of the displacements at point i, respectively, in the 0
x and 0
y directions of the sheet plane (l
=
1
,
2), at a given moment of the test.The optimization procedure starts with three different initial sets of parameters that lead to quite similar optimized values, for each test geometry, suggesting that the global optimum is reached. Moreover, the identification is split in two sequential steps (using the same cost function): (i) the first step identifies two of the Hockett-Sherby hardening law parameters (the yield stress and the saturation stress), under the assumption of isotropic (von Mises) material; the two remaining parameters of this law (contained in the exponent) are obtained by fitting the experimental results determined from the biaxial tensile test and are kept fixed during the identification procedure; (ii) in the second step, the Hill’48 parameters are identified. The cruciform geometry with the central hole (in the middle of Figure3(a)) allows achieving the best results in terms of convergence of the iterative procedure and accuracy of the identified material parameters, which was assigned to the strong heterogeneity of the kinematic field.Prates et al. [25, 26] designed a cruciform sample and developed two strategies for simultaneously identifying the parameters of the anisotropic yield criteria and isotropic hardening law of sheet metals. Both strategies use the results of the load evolution during the test and of the major and minor principal strains distributions, along the axes of the sample, at a given moment of the test, preceding the maximum load. Both strategies were numerically tested. The optimization of the design of the cruciform sample (see Figure 3(c)) was performed by means of a numerical study with the purpose of maximising the sensitivity of the test results to the values of the constitutive parameters and for allowing a wide range of strain paths, from uniaxial tension, at the arms of the sample, to near equibiaxial tension, at the centre of the sample.The work was initially addressed for the identification of the parameters of the Hill’48 yield criterion and the isotropic Swift hardening law [25]. An inverse identification was performed without resorting to the traditional optimization algorithms (e.g., gradient-based algorithms, or others); that is, a specific algorithm was built for this purpose. The inverse analysis algorithm consists of a sequence of five optimization steps, using the results of load evolution with the sample boundaries displacement during the test, for the axes 0
x and 0
y (see Figure 3(c)) and the distributions of von Mises equivalent strain along both axes of the sample, for an instant preceding and close to the maximum load. At each step, one or more parameters or a relationship between them is identified. The strategy was tested using numerically generated results of fictitious materials, which proved to be competitive, when compared with classical strategies. This allowed understanding that a sequential optimization, since properly elaborated, is clearly advantageous when compared to most commonly inverse identifications, consisting of using a unique cost function including different types of results.The inverse analysis strategy mentioned above [25] enabled a good understanding of the issues involved, namely, concerning the delineation of the sequential algorithm leading to upper accuracy. This allowed extending the strategy to more complex constitutive models (yield criteria and isotropic hardening laws). Therefore, a general inverse identification strategy that sequentially uses three distinct cost functions was developed [26]. It resorts to the Levenberg-Marquardt algorithm for sequential optimization of the parameters of the yield criteria and isotropic hardening laws. More importantly, this strategy allows the identification of parameters of several yield criteria and hardening laws. It can be used directly for a given criterion or, sequentially, starting from the Hill’48 yield criterion and then using the Hill’48 criterion solution as an initial estimate for identifying the parameters of other criteria, on the condition that can be converted into the Hill’48 yield criterion for particular values of the parameters. In the last case, this strategy is detached in two stages and has the advantage of enabling the assessment of the adequacy of a number of constitutive models to describe the experimental results, starting from a simple anisotropic criterion, the Hill’48. The first stage consists of the simultaneous identification of the hardening law, Swift and/or Voce [38], in the case, and Hill’48 yield parameters, using the results of the load evolution in function of the displacement of the grips and the equivalent strain distribution at a given moment of the test, along the axes of the sample. The hardening parameters must be separately identified for the Swift and Voce laws and the one (Swift or Voce) that better describes the results of the cruciform test (if it is possible to distinguish) is selected for further optimization. The first stage involves the sequential minimisation of the following cost functions: (12)
F
1
A
=
1
Q
1
∑
i
=
1
Q
1
P
i
N
u
m
A
-
P
i
E
x
p
P
i
E
x
p
0
x
2
+
1
Q
2
∑
i
=
1
Q
2
P
i
N
u
m
A
-
P
i
E
x
p
P
i
E
x
p
0
y
2
,
F
2
B
=
1
R
1
∑
i
=
1
R
1
ε
-
i
N
u
m
B
-
ε
-
i
E
x
p
ε
-
i
E
x
p
0
x
2
+
1
R
2
∑
i
=
1
R
2
ε
-
i
N
u
m
B
-
ε
-
i
E
x
p
ε
-
i
E
x
p
0
y
2
,where P
i
N
u
m
(
A
) and P
i
E
x
p are the experimentally determined and numerical values of the load, respectively, and ε
-
i
N
u
m
(
B
) and ε
-
i
E
x
p are the experimentally determined and numerical values of the equivalent strain, respectively; Q
1, Q
2, R
1, and R
2 are the total number of measuring points for axes 0
x and 0
y of the sample; A and B are vectors of isotropic hardening law and yield criteria parameters, respectively.The second stage allows extending the parameters identification procedure to more complex yield functions, such as Barlat’91 [8], Karafillis and Boyce [9], and Drucker+L [3], the cases studied in the work. This second stage should be performed whenever the identification carried out during the first stage is found to be not enough satisfactory to capture the experimental strain paths results, along the axes of the sample, which are not considered in minimisation during the first stage. The second stage of this inverse strategy involves the minimisation of the following cost function: (13)
F
3
C
=
1
S
1
∑
i
=
1
S
1
ρ
i
N
u
m
C
-
ρ
i
E
x
p
0
x
2
+
1
S
2
∑
i
=
1
S
2
ρ
i
N
u
m
C
-
ρ
i
E
x
p
0
y
2
,where ρ
i
N
u
m
(
C
) and ρ
i
E
x
p are the experimentally determined and numerical values of the strain path, respectively, S
1 and S
2 are the total number of measuring points for axes 0
x and 0
y of the sample, and C is the vector of yield criteria parameters. This sequential optimization procedure is a successful alternative to the parameter identification by minimising a single cost function comprising all material parameters and results of different types, as commonly found in the literature. Namely, it is concluded that this last approach can deteriorate the description of the material behaviour, concerning the load versus displacement results, and therefore the parameters and the choice of the hardening law, without apparent improvement of the description of the results of equivalent strain and strain path distributions.Zhang et al. [28] identified the parameters of Bron and Besson yield criterion [4] for both AA5086 aluminium alloy and DP980 dual-phase steel sheets, using a cruciform sample previously designed by the authors [27] and shown in Figure 3(d). The inverse identification strategy consists of minimising the gap between the experimental and numerical distributions of the major and minor strains along the diagonal direction of the sample central area, at an instant immediately before rupture, using a SIMPLEX optimization algorithm. The cost function used is defined as follows:(14)
F
A
=
∑
l
=
1
2
∑
i
=
1
n
ε
l
E
x
p
i
-
ε
l
N
u
m
A
i
2
∑
i
=
1
n
ε
l
E
x
p
i
2
,where (
ε
l
E
x
p
)
i and (
ε
l
N
u
m
(
A
)
)
i are the experimentally determined and numerical values of the principal strain components (l
=
1
,
2) at point i, respectively, n is the total number of measuring points along the diagonal path, and A is the vector of 13 parameters to be identified: four isotropy and eight anisotropy parameters of the Bron and Besson yield model and the yield stress value. The hardening of the material is modelled with isotropic hardening described by the Voce law, in case of AA5086 aluminium, and an equation based on Swift and Voce laws, for DP980 steel. In both cases, the hardening parameters are directly identified from results of the tensile test in the rolling direction and were kept fixed during the inverse identification procedure (except the yield stress).Liu et al. [29] designed a cruciform sample with a thickness-reduced central zone and four slots at each arm (see Figure 3(e)), to perform the parameters identification of a modified Voce law, describing the hardening behaviour of AA5086 aluminium sheet. The inverse analysis procedure makes use of the SIMPLEX algorithm, to minimise the difference between experimental and numerical strains measured at the centre region of the sample during the test, expressed by the following cost function:(15)
F
A
=
∑
l
=
1
2
∑
i
=
1
n
ε
l
E
x
p
i
-
ε
l
N
u
m
A
i
2
∑
i
=
1
n
ε
l
E
x
p
i
2
1
/
2
,where (
ε
l
E
x
p
)
i and (
ε
l
N
u
m
(
A
)
)
i are the experimentally determined and numerical values of the principal strain components (l
=
1
,
2) at point i, respectively, at the central point of the sample, n is the total number of the time points of simulation, and A is the vector of three hardening parameters of the modified Voce law. The experimental force evolution along the two arms of the cruciform specimen was applied to the FE model for the numerical simulations, taking into account the lack of synchronization observed in the two tensile forces on each axis of the sample. In this inverse identification of the modified Voce law parameters, three yield functions were considered: von Mises, Hill’48, and the more advanced Bron and Besson criterion, whose parameters were previously identified [27]. The identified biaxial flow stress curves were then compared with an experimental curve obtained from a uniaxial tensile test, showing that a good agreement can be achieved if an adequate yield function is used in the FE model.(2) Strategies Using Bulge Test. Most of the strain paths observed in deep-drawing components are in the range between simple tension and balanced biaxial tension, which justifies the widespread use of the cruciform specimen in the framework of the methodologies of inverse analysis. Nevertheless, this test has a strong drawback related to low deformation levels achieved, particularly in the central region of the sample, in which the strain and stress paths can be close to equibiaxial, although being dependent on the applied load or displacement ratios along the two axes of the specimen. In contrast, the bulge test allows obtaining relatively high strain values before necking, and so the flow stress curves can be assessed up to large strain values, for several biaxial strain (or stress) paths depending on the geometry of the die (circular or elliptical). In this context, a few cases of inverse analysis methodologies were developed, with the aim of identifying the parameters of work-hardening laws.Chamekh et al. [51] describe an inverse approach, based on Artificial Neural Networks (ANN), to identify the material parameters of a stainless steel (AISI 304). They use the results of the evolutions of pressure with the pole height, which are transferred to the neural network. The ANN is trained by means of curves of pressure versus displacement of the central point of the cap, generated by finite element simulations of the circular bulge test. During the training process, the network computes the weight connections, minimising the total mean squared error between the actual output and the desired output. The neural network generates an approximated function for the material parameters depending on the profile of the evolution of pressure with the pole height curve. Then, it was exploited for the identification of material parameters from experimental results. The Ludwick hardening law [33] and the Hill’48 yield criterion were selected. Therefore, the set of parameters to be identified also comprise the Lankford coefficients, r
0, r
45, and r
90. The material parameters are identified according to the following two steps: (i) the first step, using the circular bulge test, is to find the Ludwick hardening law parameters (assuming the knowledge of Lankford’s coefficients determined from the tensile tests); (ii) the second step, using the elliptical bulge test for an off axis angle of 0°, is to recalculate the Lankford’s coefficients. An elliptical die for an off axis angle of 45° is used for the validation of the parameters identification. The authors conclude the following: (i) the ANN methodology can predict acceptable combination of the values of the material parameters; (ii) once the ANN was trained, output results for a given set of input data are available almost instantaneously. Despite these conclusions, it should be noted that the values of the experimental (from the tensile test) and identified hardening coefficients are far away (the experimental and identified values are 0.67 and 0.4, resp.). The remaining identified values of the parameters differ between 20 and 30% when compared with the tensile test except for the yield stress, whose values are approximately equal.Bambach [52] explored the usability of the circular bulge test to identify the parameters of the Voce law of a fictitious material, which is considered isotropic. Initially, the membrane theory is applied to the results as in experimental cases, in order to obtain a set of parameters of the Voce, by fitting the stress versus strain results. The inverse analysis strategy proposed resorts to a gradient-based optimization algorithm, which is known for being sensitive to the initial solution. Thus, by using an initial solution, the one previously obtained with the membrane theory, it is expected to avoid convergence problems. The work gives special focus on the choice of the cost function to be minimised, making use, separately or simultaneously, of results of pressure versus pole height, pole strain versus pole height, and pole thickness versus pole height, and formulated as follows:(16)
F
p
A
=
∫
0
h
m
a
x
p
E
x
p
h
-
p
N
u
m
h
,
A
2
d
h
,
F
ε
A
=
∫
0
h
m
a
x
ε
E
x
p
h
-
ε
N
u
m
h
,
A
2
d
h
,
F
t
A
=
∫
0
h
m
a
x
t
E
x
p
h
-
t
N
u
m
h
,
A
2
d
h
,where F
p
(
A
), F
ε
(
A
), and F
t
(
A
) are cost functions defined by the pressure, p, strain, ε, and thickness, t, with pole height, h, respectively. The author concluded that the reidentification procedure (so called by the author) is significantly improved when combining the first two types of results. This significantly improves the reidentification, since it will contribute to reducing the search area where the minimum value of the objective function is located. It should be noted that this proposal is a reidentification, which has its starting point in the parameters previously obtained from the use of the membrane theory, such as in the traditional procedure recommended by the ISO standard [77]. Furthermore, it needs to resort to strain results in the pole of the cap during the test, which does not simplify the experimental procedure.Reis et al. [53] proposed an inverse methodology for determining the hardening law of metal sheets, from the results of pressure versus pole height obtained in the bulge test, involving the identification of the parameters of the Swift law. The starting point of this analysis was to realize that it is possible to achieve a unified description (i.e., overlapping) of the evolution of the pressure with the pole height, for a given value of the hardening parameter of the Swift law, regardless of the yield stress and anisotropy of the material and sheet thickness. To achieve the overlapping of such curves, appropriate multiplying factors must be used for the values of pressure and pole height, depending on the yield stresses and thicknesses ratios of the sheets and also on their anisotropy. Thereafter, an inverse analysis methodology was developed, which consists in the search for the best coincidence between pressure versus pole height experimental and reference curves, with the latter being obtained by numerical simulation assuming isotropic material behaviour with various values of the hardening parameter in the range of the material under study. This methodology, when compared with the classical strategy, proves to be an efficient alternative avoiding the use of complex devices for measuring the radius of curvature and strain at the pole of the cap, during the bulge test. Moreover, the authors claim that it is easy to implement and it is more efficient than classical approach, since (i) a unique set of numerical curves can be used within a relatively wide range of hardening coefficients, that is, covering the values usually found within one or several class of materials, without having to remake the simulations every time an identification is performed; (ii) it is not exposed to experimental errors related to the evaluation of the strain at the pole of the bulge and the use of membrane theory approach for assessment of the stress from the radius of curvature, which is usually the major source of errors.(3) Other Specimens. Güner et al. [22] proposed an inverse analysis procedure for the identification of the Yld2000-2D yield criterion parameters [7]. This study uses a notched specimen submitted to a tensile test, in order to obtain an inhomogeneous deformation field. Moreover, a layer compression test is used in order to supply additional information, that is, the equibiaxial yield stress. The required data for the inverse identification are the major and minor principal strains in the sheet plane and the load and the equibiaxial yield stress. The cost function is minimised using the Levenberg-Marquardt algorithm and is written as follows:(17)
F
A
=
∑
i
=
1
n
i
n
c
∑
j
=
1
n
e
l
e
m
ε
1
,
i
j
N
u
m
A
-
ε
1
,
i
j
E
x
p
2
+
ε
2
,
i
j
N
u
m
A
-
ε
2
,
i
j
E
x
p
2
+
C
1
∑
i
=
1
n
i
n
c
P
i
N
u
m
A
-
P
i
E
x
p
2
+
C
2
σ
b
A
-
σ
b
E
x
p
2
,where ε
1 and ε
2 are the principal strains measured at each tool displacement increment i, at each element of the optical measurements, j; P represents the load and σ
b the equibiaxial yield stress; C
1 and C
2 are scale factors. The value of σ
b
(
A
) is analytically calculated at each iteration with the values of A predicted by the optimization algorithm. Different alternative orientations of the specimen with the rolling direction (0°, 45°, and 90°) and configurations of the cost function (setting one of the three terms of the cost function equal to zero) were considered to test the inverse procedure. The authors highlight the importance of including strain information and the equibiaxial stress on the cost function, in order to improve the characterization of the anisotropy coefficients.Pottier et al. [23] developed an out-of-plane testing procedure for the simultaneous identification of Hill’48 yield criterion and Ludwick hardening law parameters of a rolled titanium sheet. Figure 4 illustrates the experimental setup and the geometry of the sample developed by the authors. A hemispherical punch applies a prescribed displacement normal to the sheet plane, at the centre of the surface of the sample, using a simple uniaxial tensile test machine. Two cameras are located on the opposite side of the sample and the components of the displacement fields along the 0
x, 0
y, and 0
z axes are captured during the test using stereo digital image correlation. The sample was designed in order to exhibit multiaxial stress states, including shear, tension, and biaxial stretching. The numerical displacements fields along the 0
x, 0
y, and 0
z axes and the global load are obtained from a numerical model of the test and compared to the experimental ones using a single cost function, formulated as follows:(18)
F
A
=
∑
i
,
j
=
1
N
t
,
n
u
x
,
i
j
N
u
m
A
-
u
x
,
i
j
E
x
p
m
a
x
i
,
j
u
x
E
x
p
2
+
∑
i
,
j
=
1
N
t
,
n
u
y
,
i
j
N
u
m
A
-
u
y
,
i
j
E
x
p
m
a
x
i
,
j
u
y
E
x
p
2
+
∑
i
,
j
=
1
N
t
,
n
u
z
,
i
j
N
u
m
A
-
u
z
,
i
j
E
x
p
m
a
x
i
,
j
u
z
E
x
p
2
+
n
P
i
N
u
m
A
-
P
i
E
x
p
m
a
x
i
P
E
x
p
2
1
/
2
,where N
t is the number of time steps considered and n is the number of measured points; u
x, u
y, and u
z are the displacements along the 0
x, 0
y, and 0
z axes, respectively, P is the global load and A is the vector of parameters to identify; the subscripts Num and Exp refer to numerical and experimental results and i and j refer to the number of time steps and measured points, respectively. The minimisation of the cost function is performed with the Levenberg-Marquardt algorithm. To assess the quality of the identified set of constitutive parameters, the authors performed deep-drawing tests of a circular cup. Moreover, additional identifications of the constitutive parameters of the material were performed, following two different strategies: a classic strategy, based on three tensile tests cut along three different directions in the sheet plane, and an inverse identification strategy using heterogeneous planar shear-like tests, previously proposed [59]. The experimental results of the earing profile of the circular cup were then compared with the numerically predicted results from the different parameter identification strategies. The authors concluded that the use of the nonplanar sample allows a more accurate prediction of the earing profile than the planar shear-like tests and the three tensile tests, since the nonplanar sample test covers a wider range of strain paths. However, in our opinion, the chosen mechanical test shows a relative degree of complexity and has the inconvenience of presenting contact with friction between the punch and the metal sheet, which always raises questions concerning the impact of contact with friction on the identification of the constitutive parameters.Figure 4
Schematic representation of the out-of-plane test developed by Pottier et al. [23].


### 4.1.2. Identification of Kinematic Hardening

The identification of the kinematic hardening plays a significant role when phenomena such as the Bauschinger effect and permanent softening, due to reverse of strain path and other strain path changes, are relevant for the subsequent plastic deformation behaviour. Depending on the material and the strain path changes occurring in the deep-drawing process, these phenomena can be more or less noticeable and relevant. The Bauschinger effect is also associated with the springback, due to premature yielding after reversal of strain path [14]. Springback is the strain recovery when forming loads are removed, and its magnitude, which depends on the flow stress value [11] and Young’s modulus, can also be influenced by the Bauschinger effect [78]. Therefore, the proper modelling and identification of the kinematic hardening are also important in order to efficiently predict the springback.In this sense, Eggertsen and Mattiasson [54] were focused on the identification of the kinematic hardening law with the main concern of the accurate prediction of springback. The goal was to select the model able to accurately describe kinematic hardening features, such as the early reyielding, transient behaviour, work-hardening stagnation, and permanent softening, taking also into account the complexity on the evaluation of its parameters. Three-point bending tests (see Figure 5) on four typical materials from car manufacturing industry were performed: two dual-phase steels (TKS-DP600HF and SSAB-DP600), from different suppliers and with different thicknesses, a mild steel (Voest-DX56D), and an interstitial-free steel (TKS-220IF). Five different hardening models were considered: (i) a pure isotropic hardening law, used as comparative reference; (ii) a mixed isotropic-kinematic law [79, 80]; (iii) the Armstrong-Frederick model [44]; (iv) the Geng-Wagoner hardening law [47]; and (v) the Yoshida-Uemori hardening law [14]. The hardening parameters of all models were determined by inverse analysis, where the difference between the experimentally and numerically generated load-displacement curves of the three-point bending test is minimised. The inverse identification of the material parameters was performed resorting to Response Surface Methodology, using the following cost function:(19)
F
A
=
1
n
∑
i
=
1
n
w
i
P
i
N
u
m
A
-
P
i
E
x
p
s
i
2
,where P
i
N
u
m
(
A
) and P
i
E
x
p represent the calculated and measured values of the punch load as a function of the vector of hardening parameters A, respectively; s
i is the residual scale factor; and w
i is the weight applied to each component of the cost function. Both s
i and w
i were set equal to 1. The authors conclude that Yoshida-Uemori model provides the best result for all materials, while the isotropic hardening model gives the worst result. However, taking into account the accuracy and the complexity of the hardening model, the authors point out that Geng-Wagoner law corresponds to a better compromise. In fact, they state that about 30 simulations are needed to optimize the parameters of the mixed isotropic-kinematic hardening law, while up to 170 simulations are required to optimize the parameters of Yoshida-Uemori hardening law.Figure 5
Device for three-point bending test [54]: (a) experimental setup and (b) schematic representation of horizontal and vertical views.
(a)
(b)Pereira et al. [56] outline an inverse analysis methodology for simultaneously identifying the parameters of the isotropic Swift law and the Lemaître-Chaboche kinematic law [13] of metal sheets, using a reverse shear test. The outlined strategy uses a modified shear sample with a cylindrical notch along the axis of the sample, in order to confine the plastic deformation within the entire gauge section, which is not the case for the classical shear samples with constant thickness (see Figure 6). The geometry of the cylindrical notch was defined in order to ensure that all strain values, between the maximum (in the centre of the notch) and the minimum (zero, along the edge of the notch), are as much as possible represented at the moment of the strain path reversing. The geometry of the sample allows that the boundary conditions of experimental tests are accurately reproduced numerically and avoids the errors in the experimental determination of the stress versus strain curves, used in traditional methodologies, whose accuracy requires homogeneity of the stress and strain fields in the sample. The inverse analysis methodology consists of minimising the gap between experimental and numerical load versus displacement curves by making variations of the constitutive parameters, using the Levenberg-Marquardt algorithm. The following cost function is used:(20)
F
A
=
1
q
∑
i
=
1
q
P
i
E
x
p
-
P
i
N
u
m
A
P
i
E
x
p
2
+
1
p
∑
i
=
1
p
P
i
E
x
p
-
P
i
N
u
m
A
P
i
E
x
p
2
,where P
i
E
x
pand P
i
N
u
m
(
A
) are, respectively, the experimental and numerical values of load for the same tool displacement; q and p are the total number of points in the forward and reverse paths, respectively; and A is the vector of the parameters to be identified. The parameters of the Hill’48 yield criterion that best describe the anisotropy of the fictitious materials used in this work (described by Drucker+L yield criterion) were identified following the methodology mentioned above, proposed by Prates et al. [25]. This methodology also allows identifying the parameters of the isotropic Swift law, which were used as first estimate in this inverse analysis. If no identification of the Swift law parameters was previously performed, the first estimate of isotropic Swift hardening parameters can be obtained adopting typical values for the material under study. It is appropriate to experimentally test this methodology.Figure 6
Geometry and dimensions of the classical and proposed notched shear specimens used in [56].Yin et al. [55] proposed the use of the twin bridge cyclic shear test proposed by Brosius et al. [81] to evaluate the Bauschinger effect on three classes of steel sheets, DC06, DP600, and TRIP700. The twin bridge shear specimen has two gauge areas that are simultaneously deformed when a moment is applied. Due to the moment application, instead of load, no unwanted reaction moment is created, in contrast with the ASTM shear specimen. The inverse identification strategy involves minimising the difference between experimental and numerical angular moment results, resorting to a Trust Region Reflective Method [82], to identify the parameters of isotropic Voce hardening law and Armstrong-Frederick kinematic hardening law. The cost function is described as follows:(21)
F
A
=
∑
i
=
1
q
M
i
E
x
p
-
M
i
N
u
m
A
2
+
∑
i
=
1
p
M
i
E
x
p
-
M
i
N
u
m
A
2
,where M
i
E
x
p and M
i
N
u
m
(
A
) are, respectively, the experimental and numerical values of moment at the same rotation angles; q and p are the total number of points in the forward and backward paths, respectively; and A is the vector of kinematic hardening parameters of the Lemaître-Chaboche law and two of the three parameters of the Voce law. The remaining Voce law parameter, the yield stress, is obtained from uniaxial tensile tests in the rolling direction and is kept fixed during the identification procedure.


## 4.1.1. Identification of Isotropic Hardening and/or Yield Criterion Parameters

(1) Strategies Using Cruciform Specimens. There has been a steadily growing interest in developing inverse identification strategies supported by the use of the biaxial tensile test of cruciform specimens, coupled with full-field displacement or strain measurements (e.g., [21, 24–29]). In general, this test allows (i) strain paths ranging from uniaxial tension (in the arms region of the specimen) to balanced biaxial tension (in the centre section of the specimen), (ii) high strain gradients from the centre region of the specimen to the end of the arms region, and (iii) no contact between surfaces and therefore no friction. Also, by changing the load and/or displacement ratio over the two normal loading axes, it is possible to achieve several biaxial stress states in the central region of the specimen. However, this kind of test only permits attaining low values of equivalent plastic strain (close to those obtained in uniaxial tension) before instability occurs and no occurrence of out-of-plane shear stresses is observed (i.e., the test is insensitive to the material parameters associated with these stresses). Figure 3 shows a set of cruciform geometries proposed in the literature for identifying material parameters using FEMU strategies.Figure 3
Cruciform specimen geometries proposed in the literature for FEMU strategies: (a) Schmaltz and Willner [24]; (b) Cooreman et al. [21]; (c) Prates et al. [25, 26]; (d) Zhang et al. [27, 28]; and (e) Liu et al. [29].
(a)
(b)
(c)
(d)
(e)In this context, Cooreman et al. [21] proposed the use of the biaxial tensile test on a perforated cruciform specimen (see Figure 3(b)), to simultaneously identify the material parameters of Swift hardening law [34] and Hill’48 yield criterion that describe the plastic behaviour of a 0.8 mm thick DC06 sheet steel. The identification makes use of strain field data measured at the central region of the specimen (see dashed area in Figure 3(b)) taken at 7 distinct load steps and iteratively compared with their numerical counterparts using Gauss-Newton algorithm, using the following cost function: (10)
F
A
=
1
3
n
∑
i
=
1
n
ε
x
x
E
x
p
i
-
ε
x
x
N
u
m
A
i
2
+
ε
y
y
E
x
p
i
-
ε
y
y
N
u
m
A
i
2
+
ε
x
y
E
x
p
i
-
ε
x
y
N
u
m
A
i
2
1
/
2
,where ε
E
x
p
i and ε
N
u
m
A
i are the experimentally determined and numerical values of the strain components ε
x
x, ε
y
y, and ε
x
y at point i, respectively, and n is the total number of measuring points. According to the authors, the results from the inverse strategy are similar to those from classical strategies, except for the parameter ε
0 of the Swift law, which leads to clearly different yield stress values. The authors attribute the discrepancy of ε
0 results to the use of strain field results from loading steps that neglect the beginning of the test. They suggest performing the inverse identification using additional strain fields from loading steps located near the onset of plastic deformation, which in our opinion can lead to high relative errors. A simpler alternative, from the experimental point of view, would be to include in the identification a cost function considering the load-displacement curves for both axes of the cruciform sample, as adopted in other identification strategies later described.Schmaltz and Willner [24] explored the usability of the biaxial tensile test on three cruciform specimen geometries, promoting different types of heterogeneous strain fields (see Figure 3(a)), in order to identify the plastic behaviour of a 2.0 mm thick DC04 sheet steel modelled via Hill’48 yield criterion and Hockett-Sherby hardening law [39]. The experimental and numerical displacement fields, in the regions in red (see Figure 3(a)) were compared and their difference was minimised using the Levenberg-Marquardt algorithm, with the following cost function:(11)
F
A
=
∑
l
=
1
2
∑
i
=
1
n
u
l
E
x
p
i
-
u
l
N
u
m
A
i
2
,where u
l
E
x
p
i and u
l
N
u
m
A
i are the experimentally determined and the numerically predicted values of the displacements at point i, respectively, in the 0
x and 0
y directions of the sheet plane (l
=
1
,
2), at a given moment of the test.The optimization procedure starts with three different initial sets of parameters that lead to quite similar optimized values, for each test geometry, suggesting that the global optimum is reached. Moreover, the identification is split in two sequential steps (using the same cost function): (i) the first step identifies two of the Hockett-Sherby hardening law parameters (the yield stress and the saturation stress), under the assumption of isotropic (von Mises) material; the two remaining parameters of this law (contained in the exponent) are obtained by fitting the experimental results determined from the biaxial tensile test and are kept fixed during the identification procedure; (ii) in the second step, the Hill’48 parameters are identified. The cruciform geometry with the central hole (in the middle of Figure3(a)) allows achieving the best results in terms of convergence of the iterative procedure and accuracy of the identified material parameters, which was assigned to the strong heterogeneity of the kinematic field.Prates et al. [25, 26] designed a cruciform sample and developed two strategies for simultaneously identifying the parameters of the anisotropic yield criteria and isotropic hardening law of sheet metals. Both strategies use the results of the load evolution during the test and of the major and minor principal strains distributions, along the axes of the sample, at a given moment of the test, preceding the maximum load. Both strategies were numerically tested. The optimization of the design of the cruciform sample (see Figure 3(c)) was performed by means of a numerical study with the purpose of maximising the sensitivity of the test results to the values of the constitutive parameters and for allowing a wide range of strain paths, from uniaxial tension, at the arms of the sample, to near equibiaxial tension, at the centre of the sample.The work was initially addressed for the identification of the parameters of the Hill’48 yield criterion and the isotropic Swift hardening law [25]. An inverse identification was performed without resorting to the traditional optimization algorithms (e.g., gradient-based algorithms, or others); that is, a specific algorithm was built for this purpose. The inverse analysis algorithm consists of a sequence of five optimization steps, using the results of load evolution with the sample boundaries displacement during the test, for the axes 0
x and 0
y (see Figure 3(c)) and the distributions of von Mises equivalent strain along both axes of the sample, for an instant preceding and close to the maximum load. At each step, one or more parameters or a relationship between them is identified. The strategy was tested using numerically generated results of fictitious materials, which proved to be competitive, when compared with classical strategies. This allowed understanding that a sequential optimization, since properly elaborated, is clearly advantageous when compared to most commonly inverse identifications, consisting of using a unique cost function including different types of results.The inverse analysis strategy mentioned above [25] enabled a good understanding of the issues involved, namely, concerning the delineation of the sequential algorithm leading to upper accuracy. This allowed extending the strategy to more complex constitutive models (yield criteria and isotropic hardening laws). Therefore, a general inverse identification strategy that sequentially uses three distinct cost functions was developed [26]. It resorts to the Levenberg-Marquardt algorithm for sequential optimization of the parameters of the yield criteria and isotropic hardening laws. More importantly, this strategy allows the identification of parameters of several yield criteria and hardening laws. It can be used directly for a given criterion or, sequentially, starting from the Hill’48 yield criterion and then using the Hill’48 criterion solution as an initial estimate for identifying the parameters of other criteria, on the condition that can be converted into the Hill’48 yield criterion for particular values of the parameters. In the last case, this strategy is detached in two stages and has the advantage of enabling the assessment of the adequacy of a number of constitutive models to describe the experimental results, starting from a simple anisotropic criterion, the Hill’48. The first stage consists of the simultaneous identification of the hardening law, Swift and/or Voce [38], in the case, and Hill’48 yield parameters, using the results of the load evolution in function of the displacement of the grips and the equivalent strain distribution at a given moment of the test, along the axes of the sample. The hardening parameters must be separately identified for the Swift and Voce laws and the one (Swift or Voce) that better describes the results of the cruciform test (if it is possible to distinguish) is selected for further optimization. The first stage involves the sequential minimisation of the following cost functions: (12)
F
1
A
=
1
Q
1
∑
i
=
1
Q
1
P
i
N
u
m
A
-
P
i
E
x
p
P
i
E
x
p
0
x
2
+
1
Q
2
∑
i
=
1
Q
2
P
i
N
u
m
A
-
P
i
E
x
p
P
i
E
x
p
0
y
2
,
F
2
B
=
1
R
1
∑
i
=
1
R
1
ε
-
i
N
u
m
B
-
ε
-
i
E
x
p
ε
-
i
E
x
p
0
x
2
+
1
R
2
∑
i
=
1
R
2
ε
-
i
N
u
m
B
-
ε
-
i
E
x
p
ε
-
i
E
x
p
0
y
2
,where P
i
N
u
m
(
A
) and P
i
E
x
p are the experimentally determined and numerical values of the load, respectively, and ε
-
i
N
u
m
(
B
) and ε
-
i
E
x
p are the experimentally determined and numerical values of the equivalent strain, respectively; Q
1, Q
2, R
1, and R
2 are the total number of measuring points for axes 0
x and 0
y of the sample; A and B are vectors of isotropic hardening law and yield criteria parameters, respectively.The second stage allows extending the parameters identification procedure to more complex yield functions, such as Barlat’91 [8], Karafillis and Boyce [9], and Drucker+L [3], the cases studied in the work. This second stage should be performed whenever the identification carried out during the first stage is found to be not enough satisfactory to capture the experimental strain paths results, along the axes of the sample, which are not considered in minimisation during the first stage. The second stage of this inverse strategy involves the minimisation of the following cost function: (13)
F
3
C
=
1
S
1
∑
i
=
1
S
1
ρ
i
N
u
m
C
-
ρ
i
E
x
p
0
x
2
+
1
S
2
∑
i
=
1
S
2
ρ
i
N
u
m
C
-
ρ
i
E
x
p
0
y
2
,where ρ
i
N
u
m
(
C
) and ρ
i
E
x
p are the experimentally determined and numerical values of the strain path, respectively, S
1 and S
2 are the total number of measuring points for axes 0
x and 0
y of the sample, and C is the vector of yield criteria parameters. This sequential optimization procedure is a successful alternative to the parameter identification by minimising a single cost function comprising all material parameters and results of different types, as commonly found in the literature. Namely, it is concluded that this last approach can deteriorate the description of the material behaviour, concerning the load versus displacement results, and therefore the parameters and the choice of the hardening law, without apparent improvement of the description of the results of equivalent strain and strain path distributions.Zhang et al. [28] identified the parameters of Bron and Besson yield criterion [4] for both AA5086 aluminium alloy and DP980 dual-phase steel sheets, using a cruciform sample previously designed by the authors [27] and shown in Figure 3(d). The inverse identification strategy consists of minimising the gap between the experimental and numerical distributions of the major and minor strains along the diagonal direction of the sample central area, at an instant immediately before rupture, using a SIMPLEX optimization algorithm. The cost function used is defined as follows:(14)
F
A
=
∑
l
=
1
2
∑
i
=
1
n
ε
l
E
x
p
i
-
ε
l
N
u
m
A
i
2
∑
i
=
1
n
ε
l
E
x
p
i
2
,where (
ε
l
E
x
p
)
i and (
ε
l
N
u
m
(
A
)
)
i are the experimentally determined and numerical values of the principal strain components (l
=
1
,
2) at point i, respectively, n is the total number of measuring points along the diagonal path, and A is the vector of 13 parameters to be identified: four isotropy and eight anisotropy parameters of the Bron and Besson yield model and the yield stress value. The hardening of the material is modelled with isotropic hardening described by the Voce law, in case of AA5086 aluminium, and an equation based on Swift and Voce laws, for DP980 steel. In both cases, the hardening parameters are directly identified from results of the tensile test in the rolling direction and were kept fixed during the inverse identification procedure (except the yield stress).Liu et al. [29] designed a cruciform sample with a thickness-reduced central zone and four slots at each arm (see Figure 3(e)), to perform the parameters identification of a modified Voce law, describing the hardening behaviour of AA5086 aluminium sheet. The inverse analysis procedure makes use of the SIMPLEX algorithm, to minimise the difference between experimental and numerical strains measured at the centre region of the sample during the test, expressed by the following cost function:(15)
F
A
=
∑
l
=
1
2
∑
i
=
1
n
ε
l
E
x
p
i
-
ε
l
N
u
m
A
i
2
∑
i
=
1
n
ε
l
E
x
p
i
2
1
/
2
,where (
ε
l
E
x
p
)
i and (
ε
l
N
u
m
(
A
)
)
i are the experimentally determined and numerical values of the principal strain components (l
=
1
,
2) at point i, respectively, at the central point of the sample, n is the total number of the time points of simulation, and A is the vector of three hardening parameters of the modified Voce law. The experimental force evolution along the two arms of the cruciform specimen was applied to the FE model for the numerical simulations, taking into account the lack of synchronization observed in the two tensile forces on each axis of the sample. In this inverse identification of the modified Voce law parameters, three yield functions were considered: von Mises, Hill’48, and the more advanced Bron and Besson criterion, whose parameters were previously identified [27]. The identified biaxial flow stress curves were then compared with an experimental curve obtained from a uniaxial tensile test, showing that a good agreement can be achieved if an adequate yield function is used in the FE model.(2) Strategies Using Bulge Test. Most of the strain paths observed in deep-drawing components are in the range between simple tension and balanced biaxial tension, which justifies the widespread use of the cruciform specimen in the framework of the methodologies of inverse analysis. Nevertheless, this test has a strong drawback related to low deformation levels achieved, particularly in the central region of the sample, in which the strain and stress paths can be close to equibiaxial, although being dependent on the applied load or displacement ratios along the two axes of the specimen. In contrast, the bulge test allows obtaining relatively high strain values before necking, and so the flow stress curves can be assessed up to large strain values, for several biaxial strain (or stress) paths depending on the geometry of the die (circular or elliptical). In this context, a few cases of inverse analysis methodologies were developed, with the aim of identifying the parameters of work-hardening laws.Chamekh et al. [51] describe an inverse approach, based on Artificial Neural Networks (ANN), to identify the material parameters of a stainless steel (AISI 304). They use the results of the evolutions of pressure with the pole height, which are transferred to the neural network. The ANN is trained by means of curves of pressure versus displacement of the central point of the cap, generated by finite element simulations of the circular bulge test. During the training process, the network computes the weight connections, minimising the total mean squared error between the actual output and the desired output. The neural network generates an approximated function for the material parameters depending on the profile of the evolution of pressure with the pole height curve. Then, it was exploited for the identification of material parameters from experimental results. The Ludwick hardening law [33] and the Hill’48 yield criterion were selected. Therefore, the set of parameters to be identified also comprise the Lankford coefficients, r
0, r
45, and r
90. The material parameters are identified according to the following two steps: (i) the first step, using the circular bulge test, is to find the Ludwick hardening law parameters (assuming the knowledge of Lankford’s coefficients determined from the tensile tests); (ii) the second step, using the elliptical bulge test for an off axis angle of 0°, is to recalculate the Lankford’s coefficients. An elliptical die for an off axis angle of 45° is used for the validation of the parameters identification. The authors conclude the following: (i) the ANN methodology can predict acceptable combination of the values of the material parameters; (ii) once the ANN was trained, output results for a given set of input data are available almost instantaneously. Despite these conclusions, it should be noted that the values of the experimental (from the tensile test) and identified hardening coefficients are far away (the experimental and identified values are 0.67 and 0.4, resp.). The remaining identified values of the parameters differ between 20 and 30% when compared with the tensile test except for the yield stress, whose values are approximately equal.Bambach [52] explored the usability of the circular bulge test to identify the parameters of the Voce law of a fictitious material, which is considered isotropic. Initially, the membrane theory is applied to the results as in experimental cases, in order to obtain a set of parameters of the Voce, by fitting the stress versus strain results. The inverse analysis strategy proposed resorts to a gradient-based optimization algorithm, which is known for being sensitive to the initial solution. Thus, by using an initial solution, the one previously obtained with the membrane theory, it is expected to avoid convergence problems. The work gives special focus on the choice of the cost function to be minimised, making use, separately or simultaneously, of results of pressure versus pole height, pole strain versus pole height, and pole thickness versus pole height, and formulated as follows:(16)
F
p
A
=
∫
0
h
m
a
x
p
E
x
p
h
-
p
N
u
m
h
,
A
2
d
h
,
F
ε
A
=
∫
0
h
m
a
x
ε
E
x
p
h
-
ε
N
u
m
h
,
A
2
d
h
,
F
t
A
=
∫
0
h
m
a
x
t
E
x
p
h
-
t
N
u
m
h
,
A
2
d
h
,where F
p
(
A
), F
ε
(
A
), and F
t
(
A
) are cost functions defined by the pressure, p, strain, ε, and thickness, t, with pole height, h, respectively. The author concluded that the reidentification procedure (so called by the author) is significantly improved when combining the first two types of results. This significantly improves the reidentification, since it will contribute to reducing the search area where the minimum value of the objective function is located. It should be noted that this proposal is a reidentification, which has its starting point in the parameters previously obtained from the use of the membrane theory, such as in the traditional procedure recommended by the ISO standard [77]. Furthermore, it needs to resort to strain results in the pole of the cap during the test, which does not simplify the experimental procedure.Reis et al. [53] proposed an inverse methodology for determining the hardening law of metal sheets, from the results of pressure versus pole height obtained in the bulge test, involving the identification of the parameters of the Swift law. The starting point of this analysis was to realize that it is possible to achieve a unified description (i.e., overlapping) of the evolution of the pressure with the pole height, for a given value of the hardening parameter of the Swift law, regardless of the yield stress and anisotropy of the material and sheet thickness. To achieve the overlapping of such curves, appropriate multiplying factors must be used for the values of pressure and pole height, depending on the yield stresses and thicknesses ratios of the sheets and also on their anisotropy. Thereafter, an inverse analysis methodology was developed, which consists in the search for the best coincidence between pressure versus pole height experimental and reference curves, with the latter being obtained by numerical simulation assuming isotropic material behaviour with various values of the hardening parameter in the range of the material under study. This methodology, when compared with the classical strategy, proves to be an efficient alternative avoiding the use of complex devices for measuring the radius of curvature and strain at the pole of the cap, during the bulge test. Moreover, the authors claim that it is easy to implement and it is more efficient than classical approach, since (i) a unique set of numerical curves can be used within a relatively wide range of hardening coefficients, that is, covering the values usually found within one or several class of materials, without having to remake the simulations every time an identification is performed; (ii) it is not exposed to experimental errors related to the evaluation of the strain at the pole of the bulge and the use of membrane theory approach for assessment of the stress from the radius of curvature, which is usually the major source of errors.(3) Other Specimens. Güner et al. [22] proposed an inverse analysis procedure for the identification of the Yld2000-2D yield criterion parameters [7]. This study uses a notched specimen submitted to a tensile test, in order to obtain an inhomogeneous deformation field. Moreover, a layer compression test is used in order to supply additional information, that is, the equibiaxial yield stress. The required data for the inverse identification are the major and minor principal strains in the sheet plane and the load and the equibiaxial yield stress. The cost function is minimised using the Levenberg-Marquardt algorithm and is written as follows:(17)
F
A
=
∑
i
=
1
n
i
n
c
∑
j
=
1
n
e
l
e
m
ε
1
,
i
j
N
u
m
A
-
ε
1
,
i
j
E
x
p
2
+
ε
2
,
i
j
N
u
m
A
-
ε
2
,
i
j
E
x
p
2
+
C
1
∑
i
=
1
n
i
n
c
P
i
N
u
m
A
-
P
i
E
x
p
2
+
C
2
σ
b
A
-
σ
b
E
x
p
2
,where ε
1 and ε
2 are the principal strains measured at each tool displacement increment i, at each element of the optical measurements, j; P represents the load and σ
b the equibiaxial yield stress; C
1 and C
2 are scale factors. The value of σ
b
(
A
) is analytically calculated at each iteration with the values of A predicted by the optimization algorithm. Different alternative orientations of the specimen with the rolling direction (0°, 45°, and 90°) and configurations of the cost function (setting one of the three terms of the cost function equal to zero) were considered to test the inverse procedure. The authors highlight the importance of including strain information and the equibiaxial stress on the cost function, in order to improve the characterization of the anisotropy coefficients.Pottier et al. [23] developed an out-of-plane testing procedure for the simultaneous identification of Hill’48 yield criterion and Ludwick hardening law parameters of a rolled titanium sheet. Figure 4 illustrates the experimental setup and the geometry of the sample developed by the authors. A hemispherical punch applies a prescribed displacement normal to the sheet plane, at the centre of the surface of the sample, using a simple uniaxial tensile test machine. Two cameras are located on the opposite side of the sample and the components of the displacement fields along the 0
x, 0
y, and 0
z axes are captured during the test using stereo digital image correlation. The sample was designed in order to exhibit multiaxial stress states, including shear, tension, and biaxial stretching. The numerical displacements fields along the 0
x, 0
y, and 0
z axes and the global load are obtained from a numerical model of the test and compared to the experimental ones using a single cost function, formulated as follows:(18)
F
A
=
∑
i
,
j
=
1
N
t
,
n
u
x
,
i
j
N
u
m
A
-
u
x
,
i
j
E
x
p
m
a
x
i
,
j
u
x
E
x
p
2
+
∑
i
,
j
=
1
N
t
,
n
u
y
,
i
j
N
u
m
A
-
u
y
,
i
j
E
x
p
m
a
x
i
,
j
u
y
E
x
p
2
+
∑
i
,
j
=
1
N
t
,
n
u
z
,
i
j
N
u
m
A
-
u
z
,
i
j
E
x
p
m
a
x
i
,
j
u
z
E
x
p
2
+
n
P
i
N
u
m
A
-
P
i
E
x
p
m
a
x
i
P
E
x
p
2
1
/
2
,where N
t is the number of time steps considered and n is the number of measured points; u
x, u
y, and u
z are the displacements along the 0
x, 0
y, and 0
z axes, respectively, P is the global load and A is the vector of parameters to identify; the subscripts Num and Exp refer to numerical and experimental results and i and j refer to the number of time steps and measured points, respectively. The minimisation of the cost function is performed with the Levenberg-Marquardt algorithm. To assess the quality of the identified set of constitutive parameters, the authors performed deep-drawing tests of a circular cup. Moreover, additional identifications of the constitutive parameters of the material were performed, following two different strategies: a classic strategy, based on three tensile tests cut along three different directions in the sheet plane, and an inverse identification strategy using heterogeneous planar shear-like tests, previously proposed [59]. The experimental results of the earing profile of the circular cup were then compared with the numerically predicted results from the different parameter identification strategies. The authors concluded that the use of the nonplanar sample allows a more accurate prediction of the earing profile than the planar shear-like tests and the three tensile tests, since the nonplanar sample test covers a wider range of strain paths. However, in our opinion, the chosen mechanical test shows a relative degree of complexity and has the inconvenience of presenting contact with friction between the punch and the metal sheet, which always raises questions concerning the impact of contact with friction on the identification of the constitutive parameters.Figure 4
Schematic representation of the out-of-plane test developed by Pottier et al. [23].


## 4.1.2. Identification of Kinematic Hardening

The identification of the kinematic hardening plays a significant role when phenomena such as the Bauschinger effect and permanent softening, due to reverse of strain path and other strain path changes, are relevant for the subsequent plastic deformation behaviour. Depending on the material and the strain path changes occurring in the deep-drawing process, these phenomena can be more or less noticeable and relevant. The Bauschinger effect is also associated with the springback, due to premature yielding after reversal of strain path [14]. Springback is the strain recovery when forming loads are removed, and its magnitude, which depends on the flow stress value [11] and Young’s modulus, can also be influenced by the Bauschinger effect [78]. Therefore, the proper modelling and identification of the kinematic hardening are also important in order to efficiently predict the springback.In this sense, Eggertsen and Mattiasson [54] were focused on the identification of the kinematic hardening law with the main concern of the accurate prediction of springback. The goal was to select the model able to accurately describe kinematic hardening features, such as the early reyielding, transient behaviour, work-hardening stagnation, and permanent softening, taking also into account the complexity on the evaluation of its parameters. Three-point bending tests (see Figure 5) on four typical materials from car manufacturing industry were performed: two dual-phase steels (TKS-DP600HF and SSAB-DP600), from different suppliers and with different thicknesses, a mild steel (Voest-DX56D), and an interstitial-free steel (TKS-220IF). Five different hardening models were considered: (i) a pure isotropic hardening law, used as comparative reference; (ii) a mixed isotropic-kinematic law [79, 80]; (iii) the Armstrong-Frederick model [44]; (iv) the Geng-Wagoner hardening law [47]; and (v) the Yoshida-Uemori hardening law [14]. The hardening parameters of all models were determined by inverse analysis, where the difference between the experimentally and numerically generated load-displacement curves of the three-point bending test is minimised. The inverse identification of the material parameters was performed resorting to Response Surface Methodology, using the following cost function:(19)
F
A
=
1
n
∑
i
=
1
n
w
i
P
i
N
u
m
A
-
P
i
E
x
p
s
i
2
,where P
i
N
u
m
(
A
) and P
i
E
x
p represent the calculated and measured values of the punch load as a function of the vector of hardening parameters A, respectively; s
i is the residual scale factor; and w
i is the weight applied to each component of the cost function. Both s
i and w
i were set equal to 1. The authors conclude that Yoshida-Uemori model provides the best result for all materials, while the isotropic hardening model gives the worst result. However, taking into account the accuracy and the complexity of the hardening model, the authors point out that Geng-Wagoner law corresponds to a better compromise. In fact, they state that about 30 simulations are needed to optimize the parameters of the mixed isotropic-kinematic hardening law, while up to 170 simulations are required to optimize the parameters of Yoshida-Uemori hardening law.Figure 5
Device for three-point bending test [54]: (a) experimental setup and (b) schematic representation of horizontal and vertical views.
(a)
(b)Pereira et al. [56] outline an inverse analysis methodology for simultaneously identifying the parameters of the isotropic Swift law and the Lemaître-Chaboche kinematic law [13] of metal sheets, using a reverse shear test. The outlined strategy uses a modified shear sample with a cylindrical notch along the axis of the sample, in order to confine the plastic deformation within the entire gauge section, which is not the case for the classical shear samples with constant thickness (see Figure 6). The geometry of the cylindrical notch was defined in order to ensure that all strain values, between the maximum (in the centre of the notch) and the minimum (zero, along the edge of the notch), are as much as possible represented at the moment of the strain path reversing. The geometry of the sample allows that the boundary conditions of experimental tests are accurately reproduced numerically and avoids the errors in the experimental determination of the stress versus strain curves, used in traditional methodologies, whose accuracy requires homogeneity of the stress and strain fields in the sample. The inverse analysis methodology consists of minimising the gap between experimental and numerical load versus displacement curves by making variations of the constitutive parameters, using the Levenberg-Marquardt algorithm. The following cost function is used:(20)
F
A
=
1
q
∑
i
=
1
q
P
i
E
x
p
-
P
i
N
u
m
A
P
i
E
x
p
2
+
1
p
∑
i
=
1
p
P
i
E
x
p
-
P
i
N
u
m
A
P
i
E
x
p
2
,where P
i
E
x
pand P
i
N
u
m
(
A
) are, respectively, the experimental and numerical values of load for the same tool displacement; q and p are the total number of points in the forward and reverse paths, respectively; and A is the vector of the parameters to be identified. The parameters of the Hill’48 yield criterion that best describe the anisotropy of the fictitious materials used in this work (described by Drucker+L yield criterion) were identified following the methodology mentioned above, proposed by Prates et al. [25]. This methodology also allows identifying the parameters of the isotropic Swift law, which were used as first estimate in this inverse analysis. If no identification of the Swift law parameters was previously performed, the first estimate of isotropic Swift hardening parameters can be obtained adopting typical values for the material under study. It is appropriate to experimentally test this methodology.Figure 6
Geometry and dimensions of the classical and proposed notched shear specimens used in [56].Yin et al. [55] proposed the use of the twin bridge cyclic shear test proposed by Brosius et al. [81] to evaluate the Bauschinger effect on three classes of steel sheets, DC06, DP600, and TRIP700. The twin bridge shear specimen has two gauge areas that are simultaneously deformed when a moment is applied. Due to the moment application, instead of load, no unwanted reaction moment is created, in contrast with the ASTM shear specimen. The inverse identification strategy involves minimising the difference between experimental and numerical angular moment results, resorting to a Trust Region Reflective Method [82], to identify the parameters of isotropic Voce hardening law and Armstrong-Frederick kinematic hardening law. The cost function is described as follows:(21)
F
A
=
∑
i
=
1
q
M
i
E
x
p
-
M
i
N
u
m
A
2
+
∑
i
=
1
p
M
i
E
x
p
-
M
i
N
u
m
A
2
,where M
i
E
x
p and M
i
N
u
m
(
A
) are, respectively, the experimental and numerical values of moment at the same rotation angles; q and p are the total number of points in the forward and backward paths, respectively; and A is the vector of kinematic hardening parameters of the Lemaître-Chaboche law and two of the three parameters of the Voce law. The remaining Voce law parameter, the yield stress, is obtained from uniaxial tensile tests in the rolling direction and is kept fixed during the identification procedure.


## 5. Final Remarks

This review shows that a great investment has been made lately in the development of inverse strategies, namely, in FEMU strategies for the identification of parameters of constitutive laws describing the plastic behaviour of metal sheets. This includes the sample design, loading conditions, and optimization procedure and intends to make the identification of the material parameters easier and more reliable. Currently, some of these strategies allows determining simultaneously the parameters of isotropic hardening law and/or anisotropic yield criteria, using only a single test. Others only allow the evaluation of the parameters of kinematic hardening law. Also, most strategies are limited to specific types of work-hardening laws and plasticity criteria. In this context, the investment in such strategies should be directed towards the simultaneous identification of parameters of any constitutive law, including isotropic and kinematic hardening, and any anisotropic yield criterion. This must be accomplished using the results of the minimum number of mechanical tests and results, for example, from the biaxial test of the cruciform sample and the shear test, and building their own optimization procedure, eventually sequential. In fact, given that the optimization procedure for parameters identification can influence the solution, it is important to examine the possibility of resorting to sequential optimization procedures, especially when using different types of results.



---

*Source: 4152963-2016-09-20.xml*