Numerical Product Design: Springback Prediction, Compensation
and Optimization
T. Meinders
a
,I.A. Burchitz
a,b
, M.H.A. Bonte
a
, R.A. Lingbeek
a,b,c
a
University of Twente, P.O.Box 217, 7500AE Enschede, The Netherlands
b
Netherlands Institute for Metals Research, P.O.Box 5008, 2600GA Delft, The Netherlands
c
Inpro Innovationsgesellschaft fur fortgeschrittene Produktionsysteme in der Fahrzeugindustrie, Hallerstrasse 1, D-10587
Berlin, Germany
Abstract
Numerical simulations are being deployed widely for product design. However, the accuracy of the numerical to ols is
not yet always sufficiently accurate and reliable. This article focuses on the current state and recent developments in
different stages of p roduct design: springback prediction, springback compensation and optimization by Finite Element
(FE) Analysis. To improve the springback prediction by FE Analyis, guidelines regarding the m esh discretization are
provided and a new t hrough-thickness integration scheme for shell elements is launch ed. In the next stage of virtual
product design the product is compensated for springback. Currently, deformations due to springback are manually
compensated in the industry. Here, a procedure to automatically compensate the tool geometry, including the CAD
description, is presented and it is successfully applied to an industrial automotive part. The last stage in virtual
product design comprises optimization. This article presents an optimization scheme which is capable of designing
optimal and robust metal forming processes efficiently.
Key words: springback, compensation, optimization, metal forming
1. Introduction
Sheet meta l forming is a widely used production
process, e.g. in the automotive industry (o uter pan-
els, inner panels, s tiffeners ), the packaging industry
(pet food co ntainers, beverage cans) and the house-
hold appliances industry (housings, razor caps). Al-
though most sheet metal forming processes have
been successfully applied for decades, the entire pro-
cess is still not fully understood. As a result, a costly
and time consuming trial and error process must
be started to determine the proper process design,
leading to the desired product. The lack of under-
standing is becoming even more pronounced, since
Email address: v.t.meinders@ctw.utwente.nl (T.
Meinders).
industries tend to favor light construction principles,
leading to the usage of new materials and new pro-
duction processes like hydroforming and incremen-
tal forming.
Presently, Finite Element (FE) Analyses are used
to enhance the understanding of sheet metal form-
ing processe s. Althoug h FE Programs are quite so-
phisticated nowadays, their accuracy and reliabil-
ity do not yet satisfy the industrial requirements.
Amongst others, so me of the main reas ons are the in-
accuracy of the used numerical algorithms (e.g. nu-
merical integration, e lement formulation) and lack
of experience in using the FE Code. As a result, the
FE Method is not sufficiently capable of s imulat-
ing forming processes and subsequently unable to
accurately predict the springback behavior. There-
fore, r e search in these fields is necessary to improve
Preprint submitted to Elsevier 30 July 2007
the usability of numerical simulations in sheet metal
forming. This article focuses on numerics in relation
to springback prediction, as explained in Section
2. Various simplifications, introduced for making a
simulation of forming more efficient, may have a sig-
nificant influence on the accuracy of springback pre-
diction and are reanalyzed in this article. Besides, a
new integration algorithm is developed to improve
the accurac y of a springback analysis at minimal
costs.
Springback is a major problem for process-
planning engineers. In many cases the shape de-
viation of the spr ung back part and the desired
product is so large that springback co mpensation
is needed. The tools o f the sheet forming process
must be changed such tha t the product becomes
geometrically accurate after spring back. In indus-
trial practice, deformations due to springback are
compensated manually by doing ex tensive mea-
surements on prototype parts, and altering the tool
geometry by hand, which is a time cons uming pro-
cess. The FE Method can be used fo r springback
compensation. In this article a procedure is de-
scribed to automatically compensate the CAD tool
shape to obtain the desired product shape (Sec-
tion 3). The potential of this method is successfully
demonstrated by an industrial automotive part.
The final step in virtual pr oduct design is op-
timization. Cost saving and product improvement
have always been important goals in the metal form-
ing industry. To achieve these goals, metal form-
ing proc e sses need to be optimized. Until recently,
trial-and-error methods were used in factories for
process optimization. Nowadays, numerical simu-
lations, and the possibility of coupling these nu-
merical simulations to mathematical optimization
algorithms, are offering a promising alternative to
design optimal metal forming processes instead of
only feasible ones, as explained in Section 4. An
overview of possible optimization a lgorithms that
can be applied to optimize metal forming processes
using time-consuming FE simulations will be given.
A promising optimization strategy tha t assists a n
engineer to efficiently model an optimization pro-
cedure is proposed. It includes an efficient prob-
lem solving algorithm and addresse s a future trend
in metal forming simulation: optimization of robust
metal fo rming processes.
2. Springback Prediction
Springback may be defined as an e lastically-
driven cha nge of shape of a product which o ccurs
when external loads are removed. In many cases the
shape deviation of the sprung back part and the
desired product is so large that springback compen-
sation is needed to obtain the desired product. In
virtual product design, this compens ation will be
based on the springback prediction provided by a n
FE Analysis. An efficient springback comp e nsation
therefore requires an accurate spring back predic-
tion. Various experimental procedures used to study
springback revealed tha t it is a complex physical
phenomenon which involves small scale plasticity
effects and thus, depends on a deformatio n path,
crystallographic texture and its evolution [1]. There-
fore, to accurately model the pheno menon in FE
Analysis it is preferable to use physically-based ma-
terial models, which are fully capable o f describing
complex material b e havior such as the B auschinger
effect and inelastic effects during unloading.
Accuracy of the prediction of springback is also
influenced by factors that are responsible for quality
of simulation of a forming step. In the past decades,
various assumptions were introduced to make form-
ing simula tio ns more efficient at the cost of ac c u-
racy. Applicability of most of the assumptions to the
springback analysis should be reanalyzed. For exam-
ple, some studies indicate that simplification of load-
ing and unloading conditions [2,3] and assumptions
underlying a shell elements theory [4] may be con-
trary to reality and are not applicable to the spr ing-
back analysis. Spatial discretization introduces an-
other appr oximation and it is recommended that
finer meshes a nd more integration points in thick-
ness direction are needed for an accurate springback
simulation [4,5].
Some results of a numerical study, which was per-
formed to develop guidelines on a level of blank dis-
cretization in a springback analysis o f real industrial
components, are discussed in Section 2.1. An adap-
tive through-thickness integration strategy for shell
elements is presented in Section 2.2, which may help
to improve springback prediction accuracy at mini-
mal costs.
2.1. Mesh Density in Springback Analysis
The theory of FE Analysis suggests that to re-
duce an e rror due to discretization, sufficiently fine
2
meshes must be used in places of high stress gra-
dients. In sheet metal forming high stress gradients
appear in regions of abrupt changes of geometry, for
example at the tool radius. It is often recommended
in literature, that accurate springback analysis re-
quires a blank to be discretized so that an element
which is in contact with the tool r adius covers 5
◦
-10
◦
of turning angle. An angle of 5
◦
per element places
high C PU power demands and is no t desirable in
simulations of industrial compo nents. To define the
level of blank discre tization that can help to estab-
lish a balance between efficiency and accurac y of
springback analysis, a numerical study is performed
which is based on a test of bending of a beam under
tension.
In the test a strip spec imen is bent over a tool
radius (R = 10mm) and dur ing bending a constant
value of tension is applied. A number of simulations
is performed with the inhouse FE Code DiekA with
a blank discretized using Discrete Kirchhoff trian-
gular elements with 7 integration points through the
thickness. A dual phase steel DP965 with a Young’s
modulus of 205 GPa, a Poisson ratio of 0.3 a nd a
yield stress of 650 MPa is us e d in the analysis.
A value of internal bending moment is used to
quantify springback. The bending moment obtained
from a simulation with the finest mesh is chosen
to be a refere nce . A relative difference be tween the
reference bending moment and a moment calculated
from the simulations with other mesh densities gives
a relative moment e rror. Hence, the springback error
induced by the discre tization error is proportional
to this relative moment error.
Figure 1 shows results of simulations in which the
relative moment error is plotted versus the number
of elements that are in contact with the tool radius.
The tension is equal to 0.1 times the tension needed
to initiate material yielding.
0
5
10
15
20
0 5 10 15 20 25 30 35 40
Relative moment error, [%]
Number of elements
real error
trendline
Fig. 1. Variation of the relative moment error depending on
mesh density.
The figure shows that to have an accuracy of
about 1% in this springback analysis it is necessary
to minimize the dis c retization error by using about
10 elements over the tool radius, an angle of 9
◦
per
element. Simulations with other materials under
the same proce ss conditions lead to similar results
and hence not presented here.
A top-hat section test is used to check the infer-
ence, mentioned above. The test resembles the NU-
MISHEET’93 benchmark [6] except for the tools’
radii. A number of simulations is performed in which
blanks with various mesh densities are used. Results
of the simulations w ith aluminum alloy (Young’s
modulus of 71 GPa, a Poisson ratio of 0.34 and an
initial yield stress of 125 MPa) are rep orted in Fig-
ure 2, where a shap e of the blank after unloading is
shown. Starting from 9
◦
per element a round the too l
radius (the c urve with 10 elements) the final shape
does not differ significantly, verifying the rule that
a discre tization density of 10 elements over a tool
radius gives accurate results.
-90
-80
-70
-60
-50
-40
-30
-20
-10
0
0 20 40 60 80 100 120 140 160
Current z coordinate, [mm]
Current x coordinate, [mm]
3 elements
5 elements
7 elements
15 elements
20 elements
10 elements
Fig. 2. Top-hat section test. Shape after springback.
However, note that the c onclusions presented here
are applicable to situations when Kirchhoff elements
are used for the blank discretization. It is likely that
fewer higher-order elements are re quired to achieve
similar accuracy of the springback prediction. In
practice, it is beneficial to employ an adaptive mesh
refinement procedure that includes a type of poste-
riori error e stimate which will trigger mesh revision
in the regions with high stress gradients.
2.2. Advanced Integration Scheme for Shell
Elements
In FE Analysis characteristic matrices of ele-
ments are c omplex volume integrals that are s olved
3
by numerical integration. In shell elements this in-
tegration is split into an in-plane integration and
a through-thickness integration. In this article the
fo c us is on the latter.
Standard integration rules are available to perform
the through- thickness integration, e.g. Trapez oidal
integration. To describe bending effects, usually
more than one integration point in the thickness
direction is needed. When a material undergoes
plastic deformations there appea r points of discon-
tinuity in through-thickness stress profile. None
of the standard rules can guar antee accurate re-
sults in this case and the number of the integration
points required to represent the non-linear stress
profile increases [7]. It was shown that depending
on a material, process parameters and the inte-
gration rule, 10 − 50 integration points are needed
through the thickness to minimize influence of nu-
merical integration error on springback prediction
[8]. Very likely, the reason for this is the inability
of the integration schemes to accurately cope with
discontinuities. Such a broad range of the required
integration points shows inefficiency of the stan-
dard integration rules. Furthermore, using more
than 20 integration points places high demands on
computational power and is very undesirable.
To overco me this problem one may use a stra tegy
that is based on adaptive integration. An integration
rule of the strategy may change its definition and
placement of the integ ration points to adapt itself
to a varying stress profile. Ideally, integration points
will be placed in the locations where the disconti-
nuities appear, which assure that these discontinu-
ities do not influence the accuracy. Thus a number
of the integration points needed to guarantee a cer-
tain level of accuracy of spr ingback prediction may
be reduced significantly.
A flexible framework for c onstructing algorithms
for adaptive integration is suggested by Rice [9]. Us-
ing the framework, an adaptive noniterative integra-
tion strategy for shell elements is defined. It consists
of two gro ups of c omponents: an interval manager
and an interval processor. The components are de-
scribed in Figur e 3.
The first task of the interval ma nager is to lo -
cate points of discontinuity in the integrand profile.
Several points of discontinuity may be present in a
stress profile when a material undergoes cyclic bend-
ing and unbending while passing, for exa mple, a die
radius. Figure 4 shows a fictive through-thickness
stress profile which occurs in a beam after bending
and reverse bending without tension. Elastic, per -
Interval processor
- applies integration rule
Interval manager
- locates discontinuities
- defines subintervals
- adapts integration points
- prepares integration
- calculates internal variables
Fig. 3. Schematic of adaptive through-thickness integration
algorithm for shell elements.
fectly plastic material is chosen and, due to a sym-
metry, only half of the stress profiles is shown.
0
t/2
0
Thickness
Stress
z
2
z
1
σ
f
σ
t
σ
y0
σ
y0
POD
1
POD
2
1
2
Bending
Reverse bending
Fig. 4. Through-thickness stress profiles for bending/reverse
bending.
While bending, as soon as yielding of the material
is initiated, an elastic-plastic boundary or a point
of discontinuity (POD1 in Figure 4 ) app e ars in the
stress profile. σ
y0
is the initial yield stress and σ
t
is
the trial stre ss. Yielding of the material during load-
ing in the reverse direction produces a new elastic-
plastic boundary and a new point of discontinuity
(POD2). σ
f
is a fictitious elastic stress that would
exist if the material did not b e c ome plastic during
the reverse loading. It can be shown that a slope of
the fictitious elas tic stress (line 1− 2 in Figure 4) in-
tersects the neutral axis at the location of the yield
stress during the bending part of the cycle. After
some work, the position of the transition points z
1
and z
2
can be calculated [10].
Following similar considerations and assuming
linear stress variation in the region of plastic de-
formations, it is possible to define equations for
locating the points of discontinuity after bending
and reverse bending for an elastic-plastic material
with hardening. Generalization from the cons idered
uniaxial case to a multiaxial problem is done in [10].
Having found the location of discontinuities in the
stress profile the interval manager divides the com-
plete integration interval [−
t
2
,
t
2
] into several parts.
After the subdivision the discontinuities coincide
4
with endpoints of the subintervals and the integrand
on eve ry subinterval is a smoo th function. Integra-
tion points are rearr anged so that there ar e several
points inside every subinterval and two of the points
are lying on its limits. If required, the location and
number of the integration points on each subinter-
val can be adapted dep e nding on used numerical in-
tegration scheme and smoothness of the integrand.
Applying an integration rule to every subinterval
and adding results gives a numerical value of the in-
tegral.
For high flexibility, the interval proces sor may em-
ploy numerical s chemes that can perform integration
using unequally distributed points [11]. Implemen-
tation of the adaptive through-thickness integration
strategy for real material behavior requires ex tra at-
tention. As soon as an integration point is relocated
or newly introduced, there is no history information
available for this point. Therefo re, it is necessary
to calculate new values of internal variables, such
as stress and strain vectors and hardening param-
eters. These values can be calculated by interpola-
tion using the information of old integration points.
If unloading occurs, the adapted integration points
will guarantee more accurate stress resultants and,
therefore, more ac curate change of shape.
To test the perfo rmance of the adaptive integra-
tion strategy, a moment M resulting from bend-
ing o f a beam to a radius under in-plane tension T
is calcula ted. A model of bending under tension is
schematically represented in Figure 5. Simplicity of
the model allows finding a closed form solution for
the bending moment. Details of the analytical cal-
culation of the bending moment are presented in lit-
erature, for example in [8].
(a)
R
unloading
forming
C
(b)
r
A
M
T
zz
D
B
t/2
00
t
central line
neutral line
Stress
Strain
Fig. 5. (a): Deformation of a beam under combined influence
of moment and tension. (b): Typical through-thickness strain
and stress profiles.
Calculations are performed for a range of values
of the in-plane tension. Prior to integration the in-
terval manager of the adaptive integration strategy
evaluates the integrand and performs the following
tasks:
(i) identifies location of two points of discontinu-
ity (POD1) in the stress profile (symmetry);
(ii) divides integration interval [−
t
2
,
t
2
] into three
subintervals AB, BC and CD, shown in Fig-
ure 5;
(iii) adapts location of the integration points on
every subinterval.
For integration on every subinterval, the interval
processor employs a r ule that uses a natural cubic
spline to approximate the integrand.
To show the per formance of the traditional inte-
gration schemes the trapez oidal rule is used to cal-
culate the bending moment. A value of the in-pla ne
tension that causes a neutra l line shift of −0.4mm
is chosen. A set of calculations is performed with
a varying number of integration points. An error
due to applying the numerical integration is quan-
tified by finding a relative difference between values
of the bending moment c alculated analytically and
numerically. Figure 6 shows the relative moment er-
ror plotted versus the number of integration points.
Presence of points of discontinuity in the stress pro-
file leads to a considerably high error when using
less than 10 integration points.
0
20
40
60
80
100
0 10 20 30 40 50 60
Relative moment error, [%]
Number of integration points
real error
trendline
Fig. 6. Relative moment error due to i ntegration with trape-
zoidal rule.
Results of numerical integration with a fixed num-
ber of the integration po ints are shown in Figure 7,
where the relative moment erro r varies as a func-
tion of the in-plane tension. In this figure, the in-
plane tension is represented by the normalized shift
of the neutral line. The numerical integration is per-
formed using the trapezoidal rule with 50 integra-
tion points and the adaptive spline integration with
only 11 points.
The integration error, obtained when using the
trapezoidal rule, oscillates. It occur s because for dif-
ferent in-plane tension the location of the points of
discontinuity in the stress profile vary to the lo c a-
tion of the integration points; the closer the loc a-
tions are, the smaller the error [8]. This fact makes
5
