Strojniški vestnik - Journal of Mechanical Engineering 54(2008)4, 288-296	Paper received: 8.4.2008
UDC 621.778	Paper accepted: 15.5.2008
Impact of Young's Modulus Degradation on Springback Calculation in Steel Sheet Drawing
Marko Vrh1,2 - Miroslav Halilovič1 - Boris Stok1* •University of Ljubljana, Faculty of Mechanical Engineering, Slovenia 2Kovinoplastika Lož, Slovenia
The aim of the paper is to demonstrate, based on the results of respective numerical simulations of a drawing process, how springback of a steel sheet drawn part is affected by stiffness degradation, which results from the material damage evolved during the forming process. For the purpose of the presented paper the GTN (Gurson-Tvergaard-Needleman) damage constitutive model coupled with Mori-Tanaka approach has been implemented in ABAQUS/Explicit via VUMAT user material subroutine, while integration of the constitutive model is performed with a numerical scheme developed recently by the authors. Parameters of the constitutive model are identified, for authenticity of the paper, from the measurement data obtained from classical tensile tests ofcold rolled stainless steel EN1.4301, as well asfrom the measured stiffness degradation, as a function ofthe plastic strain. From a comparative analysis ofnumerical results ofthe springback behaviour, investigated on the Demeri springback cup test under assumption of different constitutive models, it can be concluded, that stiffness degradation should be incorporated in simulations of sheet metal forming. Its influence on a predicted final shape of a steel sheet drawn part is namely too significant to be neglected. © 2008 Journal of Mechanical Engineering. All rights reserved.
Keywords: drawing processes, steel sheet, springback, material damage, elastic properties, stiffness degradation
0 INTRODUCTION
A great concern in sheet metal forming design is devoted to aa appropriate prediction of the final shape of a formed part, or actually to springback, which plays a key role in a formation of this shape. Springback, a phenomenon that is related to the elastic strain recovery after removal of forming loads, is physically governed by the stress state achieved at the end of the forming process and by the elastic response of a formed part. Researchers dealing with springback prediction are mostly focused on the accurate prediction of the final stress state after forming. In this context many improvements were proposed from a numerical point of view, resulting in a development of finite elements that are capable to cover stress behaviour in a sheet more precisely [1] to [5]. Another wide attention and intense research, related to springback, was paid to the mathematical modelling of constitutive laws of respective materials, in which a special attention was paid on Bauschinger effect [6] and [7], strain dependent hardening [8] and plastic anisotropy [9].
Despite all that and though it was reported several times that effective Young's modulus drops with the increase of plastic strain [10] to [14], the elastic response of the formed part associated with the forming tools removal is in the most of computational analyses still performed according to Hooke's law using initial elastic properties. In fact, Young's modulus degradation can be noticed with several experimental techniques, like unloading of the prestrained specimen [10], ultrasonic testing [11], acoustic technique [12] and nano-indentation testing [13]. Experimental evidence, as well as theoretical results obtained from consideration of corresponding mathematical models, shows that the mechanism causing the degradation of stiffness of material is appearance of cavities and cracks in materials [12] and [14].
1 EXPERIMENTAL OBSERVATIONS AND MEASUREMENTS
Though in this paper the importance of considering stiffness degradation in the calculation of the final shape of a drawn part will be
*Corr. Author's Address: University of Ljubljana, Faculty of Mechanical Engineering, Aškerčeva 6, SI-1000 288	Ljubljana, Slovenia, boris.stok@fs.uni-lj.si
demonstrated comparatively by numerical means, i.e. by simulating springback under assumption of different constitutive models, the research will be performed by respecting macroscopic material behaviour and associated material data to its greatest extent. Thus, for the investigated stainless steel EN 1.4031 the yield curve and Young's modulus degradation have been measured on a sheet specimen. In addition, also microstructure of the stretched sheet has been observed in order to better understand physical background of the stiffness degradation.
Before proceeding we find necessary to explain the term "effective" which will be used regularly in a description of some quantities associated with material response. In the context followed in this research "effective" means that respective response is considered macroscopically, possibly hiding the actual reasons for evidenced behaviour. Thus, in a homogeneous material stiffness is directly correlated to Young's modulus of the material, and it is actually from the stiffness measurements that magnitude of the elastic modulus can be uniquely evaluated. The situation is not so straightforward when we consider a porous or damaged material. Namely, experimentally we can only evaluate stiffness of the considered material structure as a whole, i.e. including all possible inclusions and voids. In this case we will refer to effective Young's modulus which will definitely depend not only on the intrinsic material properties of the matrix, but in particular on the properties specifying departure from a homogeneous solid material.
1.1 Yield Curve
The effective yield curve of stainless steel EN 1.4031 has been obtained upon measurements
1200
200
o J-,-,-,-,-
00	0.1	0.2	0.3	0.4
< ["]
Fig. 1. Effective yield curve
from a standard tensile test that was performed on the Tira 2300 tensile test machine. The experimental results were fitted according to Ludwig's law which was found to be appropriate to fit experimental data (Fig. 1).
1.2 Effective Young's Modulus Degradation
To measure effective Young's modulus degradation as a function of equivalent plastic strain, which is for uniaxial stress case sj!q = ln (L / L0), standard specimens were first plastically prestrained in the tensile test machine to a certain degree of equivalent plastic strain. After that the thickness and width of each specimen were precisely measured in order to evaluate the respective cross-sectional area. Each specimen was then clamped again in the tensile test machine, with a dynamometer measuring force range up to ±10 kN and accuracy class being 1 (IS0376-EN10002-3). The guaranteed accuracy class of the strain transducer, which was mounted on the specimen as shown in Figure 2, is 0.1, its nominal displacement range being ±2.5 mm.
From the measured force-displacement relationship registered by elastic loading and unloading of the plastically prestrained specimen the effective elastic modulus was calculated considering Hooke's law, using interpolation of the measurement data of length, cross-sectional area and force. In order to retrieve Young's modulus degradation as a function of equivalent plastic strain the described procedure is repeated for different plastic prestrains. As it can be clearly seen from the plotted graph in
Fig. 2. Measurement of the elastic elongation
Fig. 3. Effective Young's modulus degradation
Figure 3, it is beyond all question that the evidenced degradation of the effective Young's modulus is directly correlated to the degree of plastic prestrain. The fact that in our experiment plastic prestraining was achieved under condition of uniaxial stress state certainly does not affect the general statement.
1.3 Observation of Microstructure
In order to gain a better insight into a cause of the proved stiffness degradation the plastically stretched material was inspected by electronic microscope JEOL KSM-5610. Small specimens of dimensions 10x10mm were cut out with plate shears from the prestrained steel specimens for that purpose. In the sequel we give some comments on findings of the observed microstructure.
First, we observe a specimen that was stretched until rupture, which happened at the magnitude of 0.55 of the equivalent plastic strain.
Fig. 5. Revealed voids at the edge of specimen at
£peq = 0.55 (microscope view) Figure 4 shows the detail of a fracture area and the respective material appearance along of it. The front (rolled) surface of the specimen does not give us much information about any damage in the material, even in the immediate vicinity of fracture. In the fracture area, however, this is observed to be surprisingly different. From inspection of this region it can be concluded that the main mechanism of occurred ductile damage in the investigated material is voids appearance. Since microstructure of a homogeneous stretched sheet metal must exhibit a certain level of continuation through its domain, one can assume, that microstructure in the inside of the deformed sheet metal should be similar to the microstructure in the fracture area.
Figure 5 shows the photography of the material internal structure at the cut edge of a same specimen as shown in Figure 4. By cutting a small piece of material was nipped off the surface because of a small defect in plate shears blade, thus
Fig. 4. Fracture area (microscope view)
Fig. 6. Revealed voids at the edge of specimen at £Pq = 0.2 (microscope view)
revealing fortunately the internal microstructure. From that surface, which is distinctively visible in Figure 5, damage in the form of voids may be easily perceived. Based on this evidence it can be concluded, that voids in material appear throughout the volume.
Let us now observe what is happening in the sheet metal before rupture. Again the microstructure at the cut edge of a specimen, this time prestrained to the magnitude of 0.20 of the equivalent plastic strain, was observed, where cutting again nips off a small piece of material and thus reveals the internal microstructure. Appearance of voids in microstructure can be seen in considerable smaller amount. From comparison of Figures 5 and 6 it can be concluded, that voids evolution in material depends on stretching. Their presence is obviously the main mechanism of ductile damage and the true reason of progressive stiffness degradation.
2 CONSTITUTIVE MODELLING AND NUMERICAL IMPLEMENTATION
For the purpose of this research we have adopted the isotropic GTN model which establishes the respective constitutive laws for the evolution of ductile damage in porous materials [15]. This model considers void nucleation due to plastic deformation and void growth due to hydrostatic stress, i.e. two essential elements of damage evolution besides void coalescence. In order to perform our investigation we have coupled the GTN model with a respective model considering stiffness degradation.
2.1 Plastic Potential
Plastic potential associated with the GTN model reads:

(a
-2 fq1cosh
3 % K 2 a,,
+ 93 f2) (1).
Above, oeq is Mises equivalent stress, aM is yield stress of the matrix material, aH = okk /3 is hydrostatic stress, q1, q2 and q3 are parameters of the damage model and f is void volume fraction or porosity in the material. Here, it should be reminded that due to porosity in material the effective yield stress oeS, which is obtained experimentally, does not equal the yield stress of the matrix material.
2.2 Evolution of Porosity
The law governing the porosity evolution considers two mechanisms, void growth and void nucleation, respectively:
df dfgrowth + dfnucleation
(2).
The first term on the right hand side can be formulated by considering mass conservation:
df^h = (1 - f) depkk
(3),
whereas the nucleation of voids due to microcracking and decohesion of particle-matrix interface is related to plastic deformation of the matrix material:
, = An dêp f
'dfnucleation
An =-
„V2n
exp
£F-£
(4).
Above, A follows a normal distribution about the mean nucleation strain £ with a
n
standard deviation sn. Parameter fn is maximum nucleated void volume fraction and £p is equivalent plastic strain of matrix material, obtained from the following equivalent plastic work expression:
(1 - f )oMd£p =Giyd£P
(5).
In this study decrease of strength of material due to extensive void coalescence is omitted.
2.3 Stiffness Degradation
For the characterization of the stiffness degradation due to damage Eshelby's equivalence principle and his solution of the elastic field of an ellipsoidal inclusion in an infinite elastic medium [16] can be used. As it was verified several times for different materials [17] to [19] Eshelby's principle is best combined with Mori-Tanaka's concept of average stress in a matrix [20] and [21]. Combination of the GTN model with Eshelby and Mori-Tanaka approach builds an appropriate constitutive model for a simulation of the measured material response. Here, linear and isotropic elastic law aiy = CIjkl£ekl will be used, with degradation of stiffness taken into account. According to Mori-Tanaka approach effective Young's modulus E and effective Poisson's ratio V are related to porosity for material which contains spherical voids:
s
E =
v =
(7),
2E (1 - f )(7 - 5V0) 2 (7 - 5v0 )+ f (13 - 2v0 -15 v02 ) (6)
2Vp (7 - 5V0)+ f (3 - 2V0 - 5v02 )
2 (7 - 5v0 )+ f (13 - 2v0 -15 V02)
where E0 and V0 are Young's modulus and effective Poisson's ratio of undeformed material.
2.4 Numerical Implementation Procedure
The above presented constitutive model has been implemented in a general purpose finite element code ABAQUS via VUMAT subroutine. In the implementation a new explicit integration scheme, developed recently by the authors, is used, its task being to find an appropriate increment of the plastic multiplier AA from given total strain increments AetJ. This is achieved by expanding the consistency condition, which is practically a condition of fulfilled yield criterion, O = 0, that must be respected through all the integration process during the evolution of plastic strains, into Taylor series, where higher order differentials are neglected. The numerical scheme is thus based, provided the values of state variables are known at the beginning of a considered increment, on imposing:
® + d® = 0
(8)
to be fulfilled in a considered increment. When considering (1) this leads to:
d® d® d® ® + ^-Act„ + ^-Act„ + Af = 0 (9).
With regard to the forward-Euler approach, which uses the differential form of the consistency condition, i.e. d® = 0, our approach considers the additional term ®. Though this term should be zero, because it represents a function whose value should obey the consistency condition ® = 0, numerically this is usually not true. This small difference between the two explicit schemes is, according to our experience, the key reason for a considerable improvement of the stability of numerical integration.
The remaining equations of the GTN model are the evolution equations, here expressed in the incremental form:
Aay = C1JM ( -As*) d®
Aep =oaa
a.
d®
Aep = -
J dJ aa
Act,, =
(1 - f)aM
da
del
M Ae p = daM
a.
d® (10).
daa
Af =
del (1 - f)a d®
-aa
(1 - f 4
aJ daj
da. n (1 - f )ak
AA
Considering the above evolution equations in (9) yields the increment of plastic multiplier:
a®
AA =-

a®
aa^.-
a® a® a® daM J ¿a„	a®
dOJ a "a(1 - f)aM	~df
a®
(1-f +4. f
(ii),
where the corresponding derivatives are defined as:
9zatt
d^A ( -att8j) + i^sinh
dOJ (Om ) J	aM
2a.
d® QK) - fg1?2 ak sinh
da
3®
df
- = -2-
M	(aM) (aM)
qzat, 2a
(12).
= 2 q1 cosh
/ qattA
2aM
- 2 9 f
Let us remind that according to explicit approach all the state variables appearing in equations (10) to (12) are written at the beginning of the increment. When the increment of plastic multiplier AA is calculated (11), the respective increments of the other state variables can be readily calculated using (10).
3 INVERSE IDENTIFICATION 3.1 Identification of the GTN Model Parameters
Since yield condition ® = 0 is fulfilled during plastic loading, the ratio a/aM in the case of uniaxial loading remains almost independent of the magnitude of the yield stress of matrix material. Therefore, the GTN model parameters can be identified separately from the yield curve.
Parameters q1, q2 and q, are essentially an improvement of the basic Gurson's constitutive model upgrading thus the original plastic potential. Values of those parameters are similar for all metals [22], so we adopt them accordingly, taking q1 = 1.5 , q2 = 1 and q3 = q2 = 2.25 . The remaining three parameters describing void nucleation, fn, £n and sn, can be estimated from the observed porosity evolution. The latter may be thus deduced by using equation (6):	_
E^
2 (7 - 5v0 )
f = -
1 - —
(13),
2 (7 - 5Vo ) + E (13 - 2Vo - 15Vo2 )
where E is measured effective Young's modulus at different loading stages and v0 = 0.3.
From observation of the measurements it can be concluded first, that the porosity evolution does not follow a normal distribution (Fig. 7), and second, that a linear trend of the porosity increase is reasonable to be assumed. In that case only one parameter is needed for the description of the evolution law instead of three:
df
, = A dem
(14).
Inverse identification of the parameter A is carried out using Gauss minimization method, in which the cost function, defined as a sum of squares of the differences between calculated and measured values of porosity, is iteratively minimized [23]. The inverse identification procedure gives as the final result A = 0.1458 .
3.2 Identification of the Yield Curve of Matrix Material
From the measured force F, elongation AL and initial cross-sectional area A0 effective (true)
stress offf, total logarithmic strain £ and equivalent plastic strain £^q were calculated, assuming that total strain can be additively decomposed into elastic and plastic strain £jj = £k + £?. For the identification of the yield curve of matrix material a rather simple method, based on a numerical simulation of the performed experiment, was used. The identification method consists of the following steps. In a computer simulation of the tensile test two geometrically equal specimens, say A and B as depicted in Figure 8, but having different material properties or obeying different constitutive laws, are considered. The specimens are completely separated in the analysis, but in order to provide for both the conditions of the experiment they are exposed to the same loading conditions, i.e. to a prescribed edge displacement in each time increment.
Specimen A follows the Mises constitutive model with yield curve being expressed by effective yield stress off, which corresponds to a classical continuum mechanics approach. For specimen B, on the contrary, a damage mechanics approach is applied with the GTN model prescribed. But to use the latter model yield curve of matrix material aM has to be known. Considering that from the experiment we have only effective yield curve the procedure of identifying yield curve of matrix material follows this iterative path. Assuming at the beginning that aM =ofS we expose both specimens to the prescribed loading conditions and trace the respective response, i.e. resultant axial force F(AL)and effective stress offf, during loading. Because of the porosity, which evolves in accordance with the assumed GTN model, specimen B exhibits smaller effective stress and resultant force than specimen A. Also, since the effective yield curve of the model A is calculated
Fig. 7. Growth of void volume fraction
Fig. 8. Schematic representation of models for identification of yield curve of matrix material
directly from the measured dependence F(AL), the latter should intrinsically perfectly fit the response of specimen A. The task is now to equalize curve F (AL) of specimen B with the one of specimen A, the latter being actually identical to the measured one. This is done so that the yield curve of specimen B is correspondingly scaled for the simulation in the next iteration. The scaling should however take into account that difference between the two curves may depend on loading, therefore, in each increment scaling must be done by considering the actual forces. Accordingly, for a specific increment the yield curve of specimen B to be used in the next iteration should be scaled up for the same factor as it is evidenced between the calculated forces of the two specimens:
a(i+1) =a
M °
p (i) 1 B
(15).
Since the yield curve of specimen B considered in the next iteration is higher, the resultant force should be closer to the resultant force of specimen A. The described procedure, which can be repeated until the difference between curves Fa (AL) and FB (AL) becomes small enough, gives the yield curve of matrix material, which we were looking for. The procedure has proved to be very effective since its convergence is fast. The final
Fig. 9. Effective yield stress and yield stress of matrix material
result, which is depicted in Figure 9 and is obtained in only three iterations, is characterized by the maximum relative difference between curves Fa (AL) and FB (AL) being smaller than 4-10-4 in any point of the measured force.
4 SIMULATION OF DEMERI SPRINGBACK CUP TEST (ASTM, WK8010)
For the research of how Young's modulus degradation affects calculation of springback in steel sheet drawing the Demeri springback cup test may be considered [24]. The test consists of a ring sample taken from the sidewall of a cylindrical deep drawn cup. The stress state in the ring is characterized by residual hoop stresses that evolved during deep drawing process in the wall of the cup. When the ring is split, large ring opening displacement appears due to release of those stresses (Fig. 10). The evidenced displacement can be considered as a measure of the actual springback.
In the numerical simulation a steel circular blank of 200 mm diameter and nominally 0.88 mm thick is considered. Geometry of the tools is as follows: diameter of the die is 110 mm and diameter of the punch is 100 mm with 10 mm radii on the die and the punch. Nominal depth of the drawn cup is 56 mm. The forming rate in simulation is 1m/s and the holder force is constant, its magnitude being 600 kN. To obtain a ring the cup wall is cut 11 mm and 26 mm below the cup's upper edge. Two material models have been used in simulations to study the influence of stiffness degradation and damage in material, the classical isotropic Mises model with yield curve being effective stress o fff and the GTN damage model with the material data identified as described in the previous sections. In addition, also effect of degradation of elastic properties is considered.
In the FEM simulations 2776 four-node linear shell elements with reduced integration were used for the sheet blank, while the tool is assumed
Fig. 10. Numerical simulation of Demeri springback cup test. a) circular cup, b) splitted ring
Table 1. Ring opening - comparison of the models
Used model	Ring opening [mm]
a) Mises model + initial elastic properties	156.2
b) GTN + initial elastic properties	165.0
c) GTN + degradation of elastic properties	177.2
to be rigid. The results of springback, which is defined with a ring opening displacement, are tabulated in Table 1. It is interesting to observe that even if in the considered Demeri cup test the maximum porosity in the ring reaches about 4%, the difference between the computed results is around 14%.
5 CONCLUDING REMARKS
From the comparison of numerical results of the springback behaviour, performed upon corresponding computer simulations of the considered sheet forming case and subsequent Demeri springback cup test, it can be concluded, that stiffness degradation in material has a significant influence on the final result of elastic strain recovery, and consequently on the final shape of the formed part. Therefore, the stiffness degradation in terms of Young's modulus degradation and damage occurrence cannot be neglected in calculations of springback. One way to take those effects into account is the approach, presented in this paper.
Acknowledgement
Authors would like to express their gratitude to Kovinoplastika Lož Ltd. and to prof. L. Kosec from Faculty of Natural Sciences and Technology, University of Ljubljana, who enabled experimental research, presented in the paper.
6 REFERENCES
[1]	Wagoner, R.H., Li, M. Simulation of springback, Through-thickness integration. International Journal of Plasticity, 2007, 23/ 3, p. 345-360.
[2]	Dong, W.P., Soo, I.O. A four-node shell element with enhanced bending performance for springback analysis. Comput. Methods Appl. Mech. Eng., 2004, 193, p. 2105-2138.
[3]	Guo, Y.Q., Gati, W., Naceur, H., Batoz, J.L. An efficient DKT rotation free shell element
for springback simulation in sheet metal forming. Computers and Structures, 2002, 80, p. 2299-2312.
[4]	Tabourot, L., Vacher, P., Coudert, T., Toussaint, F., Arrieux, R. Numerical determination of strain localisation during finite element simulation of deep-drawing operations. Journal of Materials Processing Technology, 2005, 159/2, p. 152-158.
[5]	Chun, B. K., Kim, H. Y., Lee, J. K. Modeling the Bauschinger effect for sheet metals, part II: applications. International Journal of Plasticity, 2002, 18/5-6, p. 597-616.
[6]	Yoshida, F., Uemori, T. A model of large-strain cyclic plasticity and its application to springback simulation. International Journal of Mechanical Sciences, 2003, 45/10, p. 1687-1702.
[7]	Gau, J.T., Kinzel, G.L. A new model for springback prediction in which the Bauschinger effect is considered. International Journal of Mechanical Sciences,
2001,	43/8, p. 1813-1832.
[8]	Geng, L., Shen, Y., Wagoner, R. H. Anisotropic hardening equations derived from reversebend testing. International Journal of Plasticity, 2002, 18/5-6, p. 743-767.
[9]	Geng, L., Wagoner, R. H. Role of plastic anisotropy and its evolution on springback. International Journal of Mechanical Sciences,
2002,	44/1, p. 123-148.
[10]	Morestin, F., Boivin, M. On the necessity of taking into account the variation in the Young modulus with plastic strain in elastic-plastic software. Nuclear Engineering and Design, 1996, 162(1), p. 107-116.
[11]	Cleveland, R. M., Ghosh, A. K. Inelastic effects on springback in metals. International Journal of Plasticity, 2002, 18(5-6), p. 769785.
[12]	Yeh, H.Y., Cheng, J. H. NDE ofmetal damage: ultrasonics with a damage mechanics model. International Journal of Solids and Structures,
2003,	40(26), p. 7285-7298.
[13]	Yang, M., Akiyama, Y., Sasaki, T. Evaluation of change in material properties due to plastic
deformation. Journal of Materials Processing Technology, 2004, 151(1-3), p. 232-236.
[14]	Bonora, N., Gentile, D., Pirondi, A., Newaz, G. Ductile damage evolution under triaxial state of stress: theory and experiments. International Journal of Plasticity, 2005, 21(5), p. 981-1007.
[15]	Gurson, A.L. Continuum Theory of Ductile Rupture by Void Nucleation and Growth: Part I—Yield Criteria and Flow Rules for Porous Ductile Materials. Journal of Engineering Materials and Technology, 1977, 99, p. 2-15.
[16]	Eshelby, J.D. The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc. Roy. Soc. London, 1957, A241, p. 376-396.
[17]	Carvalho, F.C.S., Labuz, J.F. Experiments on effective elastic modulus of two-dimensional solids with cracks and holes. International Journal of Solids and Structures, 1996, 33(28), p. 4119-4130.
[18]	Pundale, S.H., Rogers, R.J., Nadkarni, G.R. Finite element modeling of elastic modulus in ductile irons: Effect of graphite morphology. AFS Transactions, 1998, 106, p. 99-105.
[19]	Bohm, H. J., Eckschlager, A., Han, W. Multi-inclusion unit cell models for metal matrix composites with randomly oriented discontinuous reinforcements. Computational Materials Science, 2002, 25(1-2), p. 42-53.
[20]	Mori, T., Tanaka, K. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metall., 1973, 21, p. 571-574.
[21]	Zhao, Y. H., Tandon, G.P., Weng, G.J. Elastic Moduli for a Class of Porous Materials. Acta Mehanica, 1989, 76, p. 105-130.
[22]	Tvergaard, V. Influence of Voids on Shear Band Instabilities under Plane Strain Condition. International Journal of Fracture Mechanics, 1981, 17, p. 389-407.
[23]	Koc, P., Stok, B. Computer-aided identification of the yield curve of a sheet metal after onset of necking. Computational Materials Science, 2004, 31(1-2), p.155-168.
[24]	Foecke, T., Gnaeupel-Herold, T. Robustness of the sheet metal springback cup test. Metallurgical and Materials Transactions A, 2006, 37A, p. 3503-3510.