Recent developments on the analysis and optimum design of sheet metal forming parts using a simplified inverse approach

Computers and Structures 78 (2000) 133±148

www.elsevier.com/locate/compstruc

Recent developments on the analysis and optimum design of sheet metal forming parts using a simpli®ed inverse approach Y.Q. Guo, J.L. Batoz *, H. Naceur, S. Bouabdallah, F. Mercier, O. Barlet Labo. LG2mS/UPRESA 6066 du CNRS, D epartement GSM, Universit e de Technologie de Compi egne, B.P. 20529, 60205 Compi egne Cedex, France Received 30 November 1998; accepted 17 July 1999

Abstract A simpli®ed ecient ®nite element method called the inverse approach (IA) has been developed to estimate the large elasto-plastic strains in thin metallic panels obtained by deep drawing. This paper deals with the main recent developments introduced by the authors on the IA to improve its eciency in the analysis and optimum design of blank contours of complicated industrial parts. The IA mainly exploits the knowledge of the 3D shape of the ®nal workpiece. An iterative scheme is used to ®nd the original position of each material point in the initial ¯at blank after which it is possible to estimate the strains and stresses in the ®nal workpiece. Important assumptions are adopted regarding the constitutive equations (the deformation theory of plasticity) and the action of the tools (the punch, die and blank holders). The IA implies only two degrees of freedom per node even if bending e€ects are considered. In this paper, we present several recent developments: (1) The bending e€ects are taken into account using a simple triangular shell element without increasing the number of dof per node. (2) Some analytical formulas are introduced to consider the restraining forces due to the drawbeads. (3) Some improvements of resolution algorithms such as the introduction of a relaxation coecient, a damping factor and a good initial solution are realized. (4) Shape optimization of blank contours is performed using a numerical procedure based on the coupling of the IA and a sequential quadratic programming method (SQP). In this work, all sensitivities are computed analytically using the adjoint variable method. The numerical results of the IA on two benchmark tests are compared with experimental and other numerical results. The optimization procedure is applied to the blank optimum design of the Renault/Twingo dashpot cup where the objective function is de®ned to minimize the maximum of the thickness variations. Ó 2000 Civil-Comp Ltd. and Elsevier Science Ltd. All rights reserved. Keywords: Inverse approach; Sheet metal forming; Optimization; Drawbead; Large elasto-plastic strain; Shell element; Resolution algorithms

1. Introduction The production of high quality formed panels in a short time and at a low cost is an ultimate goal for manufacturing companies. To reach this goal, continuous progress are made at the design and at the ¯oor shop

*

Corresponding author.

levels of forming tools. The stage of research and development in CAD/CAM/CAE in relation with the analysis and design of auto-body panels can be found in several papers published in proceedings of international meetings [1±9,38,44]. To avoid trial and error tryout procedures, the use of sheet metal forming simulation is in constant progress in the stamping industry to evaluate the deformation paths and the forming defects such as fracture and wrinkling. Many research groups are still developing and improving ®nite element codes for the analysis when an initial

0045-7949/00/$ - see front matter Ó 2000 Civil-Comp Ltd. and Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 0 ) 0 0 0 9 5 - X

134

Y.Q. Guo et al. / Computers and Structures 78 (2000) 133±148

design of the blank and the tools is done and when the forming conditions are de®ned. Several analysis tools are described in the proceedings [1±9,38,44]. They are based on membrane, shell or even solid elements, considering static or dynamic, implicit or explicit approaches. These analysis tools can be very precise if used by well-trained engineers, but they are time consuming and need expensive computer resources. Computer simulation can also be very attractive and helpful for the process and tooling engineers to de®ne the initial blank (thickness, contour and surface), some process parameters (boundary conditions, holding forces, lubrication conditions, drawbead types and positions, etc.) and the material properties (yield stress, hardening, anisotropy, etc.). This has been recognized by some industrial and academic research and development groups. As a result a number of simpli®ed methods have been developed in the last 10 years. They are mainly based on the fact that the shape of the desired part is known. A comparison is then made with the initial ¯at blank to estimate the deformation of the ®nal product taking into account simple constitutive equations and assumptions regarding the tool actions. These simpli®ed procedures have been called di€erently: ``Geometrical mapping method'' by Gerdeen et al. [10], Sowerby et al. [11], ``Single and multi-stage forming formulations'' by Chung and Lee [12], ``One step solution'' by Sklad [13], El Mouatassim et al. [14], Liu et al. [15,16], ``ideal forming theory'' by Chung and Richmond [17], ``inverse approach'' (IA) by Batoz, Guo et al. [18±23,25,26,35,36,40,42,43], ``simpli®ed approach'' by Chateau [24]. The IA developed by the senior authors and their collaborators at UTC is based on a discretization of the ®nal workpiece by simple triangular ¯at facet shell elements. Up to 1995, only membrane CST elements and conical membrane elements were considered. Since 1995, bending e€ects are taken into account using a simple discrete Kirchho€ shell element [31±36]. For a large number of industrial applications, the membrane e€ects are dominant, but sometimes it is necessary to consider bending e€ects [27±30]. Assumptions are made regarding the action of the tools (punch and die) at the end of the forming process. Logarithmic strains and total deformation theory of plasticity are considered. The equilibrium of the workpiece leads to a set of nonlinear equations with only two dof per node even if bending e€ects are taken into account. These nonlinear equations can be solved by di€erent techniques such as the Newton±Raphson static implicit approach, the dynamic relaxation method or the dynamic explicit algorithm [18,19,40]. We observe that convergence diculties can be encountered for practical situations involving deep drawn workpieces (with almost vertical walls) and a low plastic hardening law. We also developed a simpli®ed scheme in order to estimate the trimming part of the workpiece for a given blank shape, ¯at or curved

[18,19,22]. Our inverse approach has been continuously evaluated by comparing the numerical results with experimental and other numerical results obtained by incremental approaches. The procedure has been found very ecient and quite precise at the preliminary tool design stage. Our in house academic code is called FLECHE. Based on the ideas of the simpli®ed IA two industrial codes have been developed in France. One is ISOPUNCH from SOLLAC (USINOR-SACILOR group) [25,26]. The second code is SIMEX [14] from RENAULT and SIMTECH. Both codes are routinely used at the preliminary forming design stage of car panels and thin walled structural members. Other backward simulation codes are developed and used in the industry, such as FAST3D, from FTI of Canada [15,16] and ICEM stamp from control data of Germany [41]. The forward simulation codes are appropriate when an initial set of all forming parameters is available. However, they cannot be used to generate the state variables in an automatic optimization of process parameters, mainly for prohibitive computer costs. The present IA appears to be an attractive way to generate the state variables in the optimization loops of process parameters. Since 1995, our IA has been combined with mathematical programming techniques for the optimum design of axi-symmetrical [53,54] or general 3D parts [55,56]. To our knowledge, such optimum blank design procedures have not been described so far. In this paper, we ®rst review the detailed derivation of the formulation of the IA, including the bending effects. Then, we present the analytical calculation of the restraining forces due to drawbeads. We also discuss some improvements of resolution algorithms. Several numerical results will show the eciency and accuracy of the IA. We then present the coupling of the IA and an SQP algorithm with emphasis on the computation of the sensitivities using an adjoint variable method. The objective function for an acceptable forming part and the design variables and limitations are also discussed before presenting the numerical results we have obtained for the optimum design of the blank contour of an industrial problem (Renault/Twingo dashpot cup).

2. IA formulation aspects including the bending e€ects 2.1. Kinematic relations In the IA, only two con®gurations are considered: the initial ¯at blank C 0 and the ®nal 3D workpiece C. Using the Kirchho€ assumption, the initial and ®nal position vectors of a material point can be expressed on C (Fig. 1):

Y.Q. Guo et al. / Computers and Structures 78 (2000) 133±148

2

ÿ1

‰F0 Š

ÿu;y 1 ÿ v;y ÿw;y

3 kx =k3 ky =k3 5; kz =k3

…9†

znx;y 1 ‡ zny;y znz;y

3 nx ny 5; nz

…10†

1 ÿ u;x ˆ 4 ÿv;x ÿw;x 2

1 ‡ znx;x ‰Fz Š ˆ 4 zny;x znz;x

135

where u, v, w are the displacements from p0 to p in the local coordinates. The inverse Cauchy±Green left tensor can be obtained by D En o   dx0q dx0q ˆ dxq ‰ F ŠÿT ‰ F Šÿ1 dxq

  ˆ dxq ‰ BŠÿ1 dxq ; …11† ‰BŠÿ1 ˆ ‰F ŠÿT ‰F Šÿ1 :

The thickness stretch is supposed to be constant through the thickness. So, the tensor ‰BŠÿ1 takes the following simple form: 2 3 a b 0 ÿ1 ‰BŠ ˆ 4 b c 0 5: …13† 0 0 kÿ2 3

Fig. 1. Shell geometry and kinematics.

x0q ˆ x0p ‡ z0 n0 ˆ xp ÿ up ‡ z0 n0 ;

…1†

xq ˆ xp ‡ zn;

…2†

where up is the displacement of the point p on the midsurface of the sheet, n0 and n are the normals of the midsurface at p0 and p. Let x ˆ hx y zi local orthogonal curvilinear coordinates. The deformation gradient tensors at the points q0 and q with respect to p are given by dx0q ˆ ‰F0 Šÿ1 dx;

…3†

dxq ˆ ‰Fz Š dx;

…4†

  ‰F0 Šÿ1 ˆ xp;x ÿ up;x ; xp;y ÿ up;y ; n0 =k3 ;

…5†

  ‰Fz Š ˆ xp;x ‡ zn;x ; xp;y ‡ znp;y ; n ;

…6†

where k3 ˆ z=z0 ˆ h=h0 is the thickness stretch (assumed to be constant through the thickness). The inverse deformation gradient tensor at q is obtained from Eqs. (3) and (4): dx0q

ˆ ‰F Š

ÿ1

dxq

…7†

with ‰F Šÿ1 ˆ ‰F0 Šÿ1 ‰Fz Šÿ1 : The tensors ‰F0 Šÿ1 and ‰Fz Šÿ1 take very simple forms in the local coordinate system de®ned by t1 ˆ xp;x ;

t2 ˆ xp;y ;

n ˆ t1  t2 =kt1  t2 k;

…12†

…8†

The eigenvalue calculation of ‰BŠÿ1 gives two principal plane stretches k1 ; k2 and their direction transformation matrix ‰MŠ. Then, the thickness stretch k3 is calculated by the incompressibility assumption. Finally, the logarithmic strains are obtained: 2 3 ex exy 0 ‰eŠ ˆ 4 exy ey 0 5 0 0 e3 2 3 ln k1 0 0 ˆ ‰MŠ4 0 0 5‰MŠT : …14† ln k2 0 0 ln k3 In the IA, the ®nal shape in C is fully given (except for the thickness h), so the vertical displacement and rotations at every point from C0 to C are known. The unknown quantities describing the deformation are the horizontal displacement components U and V (the thickness h is in function of U,V according to the incompressibility assumption). This means that we have only two dof per node even if shell elements, including bending e€ects, are considered. 2.2. Internal force vector of DKT12 shell element The known tri-dimensional workpiece is discretized by ¯at triangular shell elements of constant thickness having three corner nodes P1 , P2 , P3 and three mid-side nodes P4 , P5 , P6 (Fig. 2). The principle of virtual work (PVW) is considered with the ®nal con®guration C as reference:

136

Y.Q. Guo et al. / Computers and Structures 78 (2000) 133±148

Fig. 2. Shell element in the initial and ®nal con®gurations. (a) 2D local coordinates, (b) 3D global coordinates.

W ˆ

X Xÿ
e e ÿ Wext ˆ 0; We ˆ Wint elt

…15†

elt

where Wint is the internal virtual work, Wext is the external virtual work which results from the tool actions (Section 2.4). Neglecting the transverse shear e€ects of the thin sheet, the element internal virtual work is given by the following expression: 8 9 Z Z < rx =

 e ex ey cxy ry dv Wint …16† ˆ he ifrg dv ˆ : ; ve ve rxy  with

he i ˆ he i ‡ zhv i

ÿ

 h h 6z6 ; 2 2

…17†

where he i is composed of the virtual membrane and bending strains, {r} is the Cauchy plane stress vector, x and y are local coordinates in the element plane (Fig. 2), the thickness h is assumed to be constant per element. For the present ¯at triangular element, the virtual membrane strains are expressed in terms of the virtual displacement components u , v along the local coordinates x, y of the element (Fig. 2). Linear approximations are used for u and v (CST membrane element):

  he i ˆ u;x v;y u;y ‡ v;x ˆ ‰Bm Š un …18† with

 un ˆ h ui vi wi

2 y23 1 6 ‰Bm Š ˆ 4 0 2A x32

i ˆ 1; 2; 3 i;

…19†

0

0

y31

0

0

y12

0

x32 y23

0 0

0 x13

x13 y31

0 0

0 x21

x21 y12

y23 ˆ y2 ÿ y3 ; . . . ;

0

3

7 0 5; 0

where …x1 ; y1 †; …x2 ; y2 †; …x3 ; y3 † are the coordinates of the three nodes in the local reference. The virtual bending strains are expressed in a simple form considering the discrete Kirchho€ triangular plate element called DKT6 presented by Batoz and Dhatt [31]. This bending element leads to a constant strain operator [Bf ]. The DKT6 sti€ness matrix is similar to that of the hybrid stress model or the Morley incompatible displacement model based also on the Kirchho€ plate theory [31,34]. The DKT6 element was used in previous works for incremental analysis of large strain shell problems [32,33]. The DKT6 element has six dof (three transverse displacements w at the corner nodes and three mid-side rotations hs along the element sides, Fig. 2). To de®ne the virtual curvatures of the DKT6 element, the rotation components bx and by of the normal (rotations from z to x and y) are linearly expressed with semi-C 0 approximations in terms of the rotations at the mid-nodes: X bx ˆ Nk bxk ; …20† kˆ4;5;6

by ˆ

X

Nk byk ;

…21†

kˆ4;5;6

where N4 ˆ 1 ÿ 2g, N5 ˆ ÿ1 ‡ 2n ‡ 2g, N6 ˆ 1 ÿ 2n and n, g are area coordinates. Then, a Kirchho€ constraint is introduced along the element sides to express the rotation bsk in terms of the normal displacement wi and wj at the corner nodes, i.e.: Z 0

Lk

csk ds ˆ

Z

Lk 0



w;s ‡ bs



ds ˆ 0

…22†

Y.Q. Guo et al. / Computers and Structures 78 (2000) 133±148

so that

bsk ˆ

  wi ÿ wj Lk

:

…23†

The virtual rotations bnk ˆ hsk are zero since the rotations at mid-nodes are known. We ®nally obtain the matrix [Bf ] as (see Ref. [31], Section 6.4.3 for details):

 hv i ˆ bx;x



Bf



   bx;y ‡ by;x ˆ Bf un ;

by;y

2 0 0 14 0 0 ˆ A 0 0

S4 C4 ÿ S6 C6 ÿS4 C4 ‡ S6 C6 S42 ÿ C42 ÿ S62 ‡ C62

…24†

with

 Un ˆ h Ui 2

‰T Š ˆ 4

‰tŠ

Vi  3

‰t Š

‰tŠ

5;

i ˆ 1; 2; 3 i; 2

t1X ‰tŠ ˆ 4 t2X nX

2.3. Constitutive equations for a planar isotropic elastoplastic sheet In the present IA, only the initial and ®nal con®gurations (C 0 and C) are considered. The elasto-plastic deformation is assumed to be independent of the loading path. As in our previous work [18±23], the so-called Hencky±Ilyushin deformation theory of plasticity is adopted.

0 0 S5 C5 ÿ S4 C4 0 0 ÿS5 C5 ‡ S4 C4 0 0 S52 ÿ C52 ÿ S42 ‡ C42

where A is the area of the triangular element, L2k ˆ x2ji ‡ yji2 , Ck ˆ xji =Lk , Sk ˆ yji =Lk …k ˆ 4; 5; 6 for ji ˆ 21; 32; 13†. The matrices [Bm ] and [Bf ] are constant per element and unchanged during the iteration solution process. Using Eqs. (17), (18) and (24), we obtain ÿ  
  fe g ˆ ‰Bm Š ‡ z Bf un : …26†   We then transform the local un in terms of the global  Un (considering hsk ˆ 0):    un ˆ ‰T Š Un …27† …28† t1Y t2Y 5: nY

3 0 0 S6 C6 ÿ S5 C5 0 0 ÿS6 C6 ‡ S5 C5 5; 0 0 S62 ÿ C62 ÿ S52 ‡ C52

…29†

So, the element internal force vector in the global coordinates is obtained:

 e e Wint …30† ˆ Un Fint    e   R T with Fint ˆ ‰T ŠT V e ‰Bm ŠT ‡ z Bf frg dz dA. For a given value of z, the strains and stresses are constant on each element, so that    e  T Fint ˆ ‰T ŠT ‰Bm ŠT f N g ‡ Bf f M g A; …31† where the {N} and {M} are the membrane forces and bending moments. In the above equation, only {N} and {M} are changing during the iteration process. These resultant forces can be obtained by numerical integration through the thickness at the element centroõd using three or ®ve Lobatto integration points.

…25†

If we consider planar isotropic sheets, then yielding can be described by the Hill criterion with the plane stress conditions [37]: / ˆ hri‰ P Šfrg ÿ r2 ˆ 0

…32†

with hri ˆ h rx ry rxy i, where r is the equivalent yield stress. The matrix [P] is in function of the mean planar isotropic coecient r obtained from the three anisotropic coecients: 2 3 r 1 ÿ 1‡ 0 r 6 7 ÿ r 1 0 7; ‰P Š ˆ 6 4 1‡r 5 2…1‡2 r† 0 0 1‡ r  r ˆ 14…r0 ‡ 2r45 ‡ r90 †;

3

137

with

ra ˆ e_P

…a‡90 † e_PZ

…33†

.

The 0° direction is de®ned as the rolling direction. The proportional (or radial) loading assumption allows to obtain the total plastic strains in terms of the total stresses: feP g ˆ

eP ‰ P Šfrg r

…34†

with eP ˆ …heP i‰ P Šÿ1 feP g†1=2 the equivalent plastic strain. If we assume the above anisotropy (planar isotropy) and incompressibility to be valid also for the (small) elastic deformation, then the Poisson coecient is related to the mean anisotropy coecient  r: mˆ

r : 1 ‡ r

…35†

So, the total constitutive equation takes the following simple form: ÿ1 frg ˆ ES ‰ P Š feg;

…36†

138

Y.Q. Guo et al. / Computers and Structures 78 (2000) 133±148

where ES is the secant modulus of the uniaxial stress± strain curve:  1=2 e ˆ hei‰ P Šÿ1 feg :

r ES ˆ ; e

…37†

With an estimation of the total strains feg, we can compute e and ES , then estimate the stresses. This operation is performed for each numerical integration point through the thickness. 2.4. External nodal forces A number of assumptions can be made regarding the actions of the punch, die and blank holder. Some details are given in previous papers [19,20]. In the case of no friction the normal pressure force is assumed as F ˆ F n;

…38†

where n is the nodal unit normal vector to the sheet, F is the unknown force intensity. The consideration of friction between the tools and the sheet consists in modifying the direction of the resultant force in the following simpli®ed manner [19]: F ˆ F nf

…39†

p with nf ˆ …1= 1 ‡ l2 †…n ÿ lt† where l is the friction coecient and t is a unit vector de®ned by the relative movement between the tool and the sheet. An equivalent external load vector due to the friction under the blank holder can be de®ned as follows: f ˆ ÿ2lqn

u ; kuk

…40†

where f is the tangential force per unit area, qn is the normal blank holder pressure acting on the two sheet faces and u the tangential displacement (varying during the iteration process [19]). An equivalent load nodal vector can be de®ned for each element as Z ff n g ˆ

Ae

‰ N ŠT f f g dA;

…41†

where [N] is composed of the linear shape function of the triangular elements. The e€ect of drawbead can be replaced by equivalent restraining forces. For a given drawbead geometry and sheet material properties, the simpli®ed formulas of Stoughton [39] can give a satisfactory estimation. The nodal forces resulting from the tool action are assembled in the global load vector {Fext }.

3. Estimation of the restraining forces due to drawbeads In a sheet forming process, the drawbeads play an important role for the control of the material ¯ow. A good drawbead design allows us to avoid the necking, wrinkling and other imperfections. The entire simulation with the drawbead geometry involves too many small elements and local contact problems. So, in most of incremental simulation codes, the actions of the drawbeads are represented by some restraining forces and strains obtained by a separated 2D simulation [48]. Many researchers have been working on this subject to ®nd empirical or analytical formulas for drawbead restraining forces [39,49±52]. Nine [50] has done numerous experiments to study the in¯uence of the drawbead geometry, the sheet thickness and material, the holding force and friction coecient on the restraining force. Stoughton [39] has proposed a set of analytical formulas to determine the drawbead restraining force and drawbead hold down force. These formulas are based on the principle of virtual work: the work required to pull (or to restrain) a sheet through a drawbead is equal to the work required to overcome the sheet±tool friction and to bend±unbend the sheet. At each radius changing site i (Fig. 3), the restraining force Fi per unit of the width due to the bending or unbending is given by ! Z ‡ti =2 Z e‡i Fi ˆ rde dz …i ˆ 1; 2; . . . ; 6†; …42† ÿti =2

eÿ i

‡ where ti is the thickness of the sheet, eÿ i and ei are the strains before and after the ith radius change, r is the normal stress along the pull direction. The total restraining force, including the friction, was obtained by Stoughton and presented in Ref. [39]. A simple drawbead (Fig. 3 with g ˆ 0, d ˆ 2R) is studied with di€erent parameters (E ˆ 210 GPa, ry ˆ 516e0:23 , average anisotropy  r ˆ 1:6, other data in Table 1) [50]. The comparison in Table 2 shows that the results obtained by the formulas of Stoughton and the simulation of OPTRISTM (2D numerical simulation based on a dynamic explicit incremental algorithm) are both in good agreement with the experimental results of Nine. We have improved and generalized the formulas of Stoughton to consider more complicated drawbeads

Fig. 3. Drawbead geometry.

Y.Q. Guo et al. / Computers and Structures 78 (2000) 133±148

139

Table 1 Data for ®ve drawbeads Specimens

Radius (mm)

Thickness (mm)

Friction on ¯at zones

Friction on curve zones

1 2 3 4 5

4.75 4.75 5.5 5.5 5.5

0.86 0.86 0.76 0.86 0.86

0.002 0.002 0.002 0.002 0.002

0.16 0.07 0.18 0.16 0.05

(with various shapes, di€erent radius and angles of contact). The following modi®ed restraining force, including friction, is obtained (Fig. 3): ÿÿ
F ˆ F1 elh1 ‡ lHe ‡ F2 ‡ F3 e2l…h2 ‡h3 †
lh ‡ lHe ‡ F4 ‡ F5 e 4 ‡ F6 ; …43† where l is the friction coecient, h1 , h2 ‡ h3 and h4 are the angles of contact in the zones 1±2, 3±4 and 5±6 (Fig. 3), He is the elastic part of the holding force required to close the drawbead. He can be obtained by the classical beam theory: He ˆ

16Ebdt3

…44† …R1 ‡ 2R2 ‡ R3 ‡ g1 ‡ g2 †3 ÿ ÿ

with d ˆ min d; …R1 ‡ 2R2 ‡ R3 ‡ g1 ‡ g2 † R2 ry =Et , ry is the yield point, R1 , R2 , R3 are the radii in the zones 1±2, 3±4 and 5±6. The benchmark test of NumisheetÕ96 called ``limiting dome height'' has a drawbead with a complicated geometry (Fig. 4). The properties for three materials are given in Table 3. The drawbead restraining force and the binder hold down force are calculated with the modi®ed formulas of Stoughton and the 2D simulation of OPTRISTM . We note a good agreement for IF and HS steel but some di€erences for aluminum (Table 4).

4. Improvements in resolution algorithms The assemblage of the internal and external force vectors leads to the global residual vector:

Fig. 4. Geometry of the LDH drawbead (NumisheetÕ96).

f R…U ; V ; F †g ˆ fFint …U ; V †g ÿ fFext …F †g

…45†

with three unknowns per node, i.e. the global displacements components U, V and the nodal force intensity F (Eq. (39)). The knowledge of the vertical displacement allows one to eliminate the increment DF at the element level so that we have only two unknowns DU , DV at each node. The inverse approach leads to highly nonlinear equation systems because of the large displacements, large strains and elasto-plastic materials. Sometimes, convergence diculties are encountered. Di€erent resolution algorithms have been developed and tested [40]: · N±R static implicit method, · static explicit method, · method of Levenberg±Marquardt, · dynamic explicit algorithm, · dynamic relaxation method, · dynamic viscous relaxation method, etc. Among numerous resolution methods, the N±R method is very attractive (eciency, simplicity, veri®cation of equilibrium). However, this method requires a good initial solution and a well-conditioned tangent matrix. Several improvements have been made for the N±R method used in our IA. The ®rst one is to use a good initial solution. The most simple solution guess was

Table 2 Comparison of the drawbead restraining forces (DBRF) and the blank holding forces (BHF) DBRF (N/mm) 1 2 3 4 5

BHF (kN)

Nine (experimental)

Stoughton

Optris

Nine (experimental)

Stoughton

Optris

138 106 112 128 92

139 98 103 121 80

118 85 102 104 67

98 82 84 96 82

122 94 87 106 77

150 130 130 140 115

140

Y.Q. Guo et al. / Computers and Structures 78 (2000) 133±148

Table 3 Properties of the three blank materials of the LDH test (NumisheetÕ96) q (kg/mm3 )

E (GPa) IF steel HS steel Alu

ÿ6

206 206 69

7:8  10 7:8  10ÿ6 2:7  10ÿ6

m

ry (MPa)

r ˆ Ken (MPa)

0.3 0.3 0.3

158 370 155

K ˆ 526, n ˆ 0:233 K ˆ 734, n ˆ 0:156 K ˆ 488, n ˆ 0:232

Table 4 Comparison of DBRF and BHF for the LDH drawbead Material of blank IF steel HS steel Alu

Restraining force (kN/mm)

Blank holding force (kN/mm)

Optris

Stoughton

Optris

Stoughton

0.26 0.32 0.12

0.250 0.314 0.168

0.24 0.31 0.125

0.222 0.291 0.149

obtained by vertical node projection (i.e. fU g ˆ fV g ˆ f0g), and this initial solution leads to very large strains and stresses for the elements belonging to the almost vertical walls of the workpiece. Therefore, very small steps are needed to avoid divergence problems. Recently, we have developed a geometric mapping method to de®ne a good initial guess in the following manner: (1) for each node D, of the element mesh, a vertical plane is de®ned by D and the centre of gravity C, (2) the arclength from C to D is developed to locate the ``guessed'' initial position of D in the ¯at blank. This initial guess considerably reduces the number of steps and the CPU time. The standard N±R method cannot always ensure the convergence for highly nonlinear systems; a relaxation factor is introduced to keep the solution within the radius of convergence:  i KT fDU i g ˆ fRi g; …46† 

U i‡1 ˆ fU i g ‡ xi fDU i g with …0 < xi 6 1†:

5. Numerical results for forming analysis 5.1. Benchmark test ``limiting dome height'' of Numisheet '96 A limiting dome height test was proposed for the NumisheetÕ96 international conference [38]. The LDH test tooling geometry is presented in Fig. 5. This example has a complicated drawbead and a very weak hardening material in the case of the high strength steel (HS). The data of the steel sheet are E ˆ 206 GPa, m ˆ 0.3, h0 ˆ 0.92 mm, l ˆ 0.11 (friction), r ˆ 734e0:15 MPa,  r ˆ 1:04 (average anisotropy). A punch displacement of 30 mm was considered. Using the modi®ed Stoughton formulas, we obtain the drawbead restraining force (DBRF ˆ 314 N/mm). This loading is introduced in an equivalent manner at the nodes situated near the drawbead line. The ®nite element mesh shown on Fig. 6 for a quarter of the workpiece involves 1623 elements and 875 nodes (around 1700 dof). For our simpli®ed IA, we had to build the deformed middle surface of the workpiece based on geometrical considerations. In par-

The norm of the residual forces kf R…U i ‡ xi DU i †gk is a very complicated function in terms of xi . A line search technique is used to ®nd xi minimizing kf Rgk. The large strains and weak material hardening often induce an illconditioned tangent matrix and convergence troubles. A diagonal damping matrix k‰ I Š…k > 0† is introduced to improve the conditioning of the sti€ness matrix: ÿ i 
KT ‡ ki ‰ I Š fDU i g ˆ fRi g: …47† The solution of the new system will converge towards the solution of the original system if ki is chosen within certain limits. However, the factor ki decreases the convergence speed. We can reduce the ki value when the conditioning of the tangent matrix is improved in course of iterations.

Fig. 5. LDH test tooling geometry (NumisheetÕ96).

Y.Q. Guo et al. / Computers and Structures 78 (2000) 133±148

141

and Brunet using a dynamic explicit incremental approach are also reported (see ``Benchmark problems and results'' in Ref. [38]). We can observe the quality of our results using the simpli®ed IA approach for this problem. Of course, less satisfactory results are obtained at the point D and its neighbourhood due to the assumption on the shape of the ®nal workpiece. The good quality of the results is mainly attributed to the validity of the total deformation theory of plasticity for this case and to a good estimation of the drawbead restraining forces. 5.2. Benchmark test ``oil pan'' of Baden±Baden

Fig. 6. Finite element mesh for 30 mm punch height.

ticular, the free surface between the punch area and the die radius is supposed to be a tangential conical one. However, the free edge CD is unknown (D is not on the conical surface) and will be a result of the iteration procedure (in correlation with the given initial ¯at blank of 90 ´ 50 mm2 ). The analysis was performed with the shell elements in 30 iterations and 5 min of CPU time on a DEC ALPHA 3000-900 workstation. Friction between the sheet and the tools was considered. Very good results were obtained (the results for IF steel were reported in the proceedings [38]). The punch force obtained by the IA (57.1 kN) is a little overestimated with respect to the values obtained experimentally (49 kN average) and numerically using an incremental code (48.7 kN). In Fig. 7, we report the results for the principal strains at the outer surface along the x direction (AB and DC). The average experimental results and the numerical results obtained by Sabourin

This benchmark test was proposed by the international conference of Baden±Baden on ``metal forming process simulation in industry'' [6]. The diculties for this example are characterized by its great depth with almost vertical walls, the weak hardening of material and the separation in two parts. The data of the test are E ˆ 210 GPa, m ˆ 0:33, h0 ˆ 0.63 mm, l ˆ 0:165 (friction), r0 ˆ 2:12, r45 ˆ 1:73, r90 ˆ 2:53 (anisotropy), F ˆ 400 kN (blank holding force). The elasto-plastic curve is given point by point. The element mesh is composed of 10 650 triangular elements and 5385 nodes (Fig. 8). It was impossible to get the convergence with the standard N±R method and the vertical projection as initial guess because of the highly geometric and material nonlinearities. However, good results are obtained using a good initial solution guess (Section 4) and the improvements of the conditioning of the tangent matrix (Eq. (46)). The IA calculation is carried out on an HP9000/780 Workstation in only 1 h 45 min of CPU time (in 10 equal steps for the punch travel and 60 equilibrium iterations). The thickness distributions along the AB and RS lines are presented and compared with the average experimental values and other

Fig. 7. Principal strains along AB and DC in the LDH-HS test.

142

Y.Q. Guo et al. / Computers and Structures 78 (2000) 133±148

workpiece since an increase of the thickness is often correlated with wrinkling and a reduction of thickness may lead to rupture (necking). For example, the thickness increase and decrease are often limited to 15% and )25% in practice. These considerations are relevant for the de®nition of the objective function of the problem of nonlinear optimization [58] de®ned in the standard form: min J …v†; v ˆ hv1 ; . . . ; vn i: with cj …v† ˆ 0; j ˆ 1; 2; . . . ; m; gk …v† 6 0; k ˆ 1; 2; . . . ; p; vil 6 vi 6 viu ; i ˆ 1; 2; . . . ; n; v 2 Rn ;

Fig. 8. Element mesh and positions of the lines AB and RS.

numerical simulation results (Fig. 9). The global length and width of the oil pan after drawing are compared in Table 5.

…48†

where J(v) is the objective function, v is the vector of design variables, cj (v) and gk (v) are, respectively, equality and inequality constraints, and vil , viu the lower and upper geometrical constraints (bounds) on every design variable. Our objective function depends on the thickness distribution and is expressed as

6. Shape optimization problem The shape of the blank has a great in¯uence upon the quality of the ®nal workpiece. When the blank is too large, the material ¯ow in the die cavity is limited and the risk of the excessive thinning (necking) is increased. On the other hand, if the blank size is too small, wrinkling can occur. The study of the in¯uence of the blank shape has been studied both experimentally and numerically by some authors [46,47,56,57]. In this study, we assume a close relationship between the variation of the thickness of the workpiece and the de®nition of a defect-free product. For practical applications, the thickness variation appears to be a good indicator of the formability of the

Jˆ

nelt X …he ÿ h0 †p

!1=p ;

…49†

e

where h0 is initial thickness, he the element thickness. This de®nition of J(v) allows one to minimize the thickness variation. The exponent p is introduced in order to emphasize the in¯uence of maximal thickness

Fig. 9. Thickness distribution along the lines AB and RS. Table 5 Comparison of the global length and width after drawing (L0  W0 ˆ 492  380 mm2 ) Length Width

Average experimental [6]

Autoform [6]

Icem stamp [6]

LS Dyna3D [6]

Abaqus [6]

Our IA

445 305

450 305

445 305

475 350

450 325

424 300

Y.Q. Guo et al. / Computers and Structures 78 (2000) 133±148

variation. Other objective functions have been considered in Refs. [53±56]. The design variables v are the positions of the control knots describing the 3D workpiece contour (Fig. 11). A parametrization with a B-spline curve is used to reduce the number of design variables. Some constraints on the design variables are introduced in order to satisfy technological limitations (bounds). At every optimization loop, the objective function, constraints and sensitivities are computed using strains and stresses obtained by a ®nite element simulation based on the simpli®ed IA. The optimization algorithm is a sequential quadratic programming method with a line search technique [59] to solve a general nonlinear optimization problem (Eq. (48)). This algorithm is quite robust and particularly ecient for optimization problems with highly nonlinear constraints. The optimization process starts with a preliminary blank contour, a nonlinear analysis is carried out using the IA in order to obtain thickness distribution, strains and stresses. Using these results, the objective function and constraints are evaluated; then, the gradients are computed with an explicit di€erentiation presented in Section 7. These results are used to update the design variables minimizing the objective function. If the new design is acceptable, the process stops. Otherwise, the iterative process continues, producing a succession of designs, until an optimal design is achieved.

Sensitivity calculations are based on the di€erentiation of the equilibrium equations [59] (the adjoint variable method). These nonlinear equations resulting from the IA are given by …50†

where f Rg is the residual vector, fFint g and fFext g are, respectively, the internal and external force vectors (31), (39) and (41). In our problem, the objective function is a highly non linear function in terms of the design variables and nodes displacements: J ˆ J …v; U…v††;



oJ oU



oU ovi

 ;

i ˆ 1; n:

…51†

where v is the vector of the n design variables, and U the nodal displacement vector. The shape of the blank contour (in 2D) is parametrized by using cubic B-spline functions. The coordinates of the k control knots ([Xj , Yj ] with j ˆ 1; k) of the Bspline curve are used as design variables (vector v). The sensitivity analysis consists in calculating the total derivative of J at (v, U):

…52†

The derivation of the IA equilibrium equations (50) leads to  ‰KT Š

oU ovi



 ˆÿ

 oR ; ovi

i ˆ 1; n;

…53†

where ‰KT Š is the nonsymmetric tangent sti€ness matrix. Using Eq. (53), Eq. (52) becomes   dJ oJ oR ˆ ‡ hP i ; i ˆ 1; n …54† dvi ovi ovi with f P g the solution of the following linear system:    T oJ KT f P g ˆ ÿ : oU

…55†

The total derivatives of all constraints, can be performed in the same manner, by replacing J with gj (jth constraint function). The value of the term foR=ovi g can be used not only for the sensitivities of the objective function, but also for the sensitivities of the constraints. In practice, for a given objective function J, its partial derivatives with respect to the design variable vi can be calculated as follows: k k oJ X oJ oxj X oJ oyj ˆ ‡ ; ovi ox oy ov j i j ovi jˆ1 jˆ1

7. Sensitivity analysis

f Rg ˆ fFint g ÿ fFext g ˆ 0;

dJ oJ ˆ ‡ dvi ovi

143

i ˆ 1; n;

…56†

where xj and yj are the coordinates of the nodes belonging to the workpiece contour. The coordinates x, y (piecewise cubic B-spline) are expressed in terms of the coordinates of the control knots (Xj , Yj ) by 8 x ˆ 1=6‰…1 ÿ s†3 Xj ‡ …4 ÿ 6s2 ‡ 3s3 †Xj‡1 > > < ‡…1 ‡ 3s ‡ 3s2 ÿ 3s3 †Xj‡2 ‡ s3 Xj‡3 Š; 3 > ÿ 6s2 ‡ 3s3 †Yj‡1 > : y ˆ 1=6‰…1 ÿ s† Yj ‡ …4 2 ‡…1 ‡ 3s ‡ 3s ÿ 3s3 †Yj‡2 ‡ s3 Yj‡3 Š

…57†

with 0 6 s 6 1: These coordinates (Xj , Yj ) constitute the design variable vector v, so for each control knot j (j ˆ 1; k), we can easily compute the values of oxj =ovi and oyj =ovi . The two other terms oJ =oxj and oJ =oyj in (Eq. (56)) can be computed by di€erentiation of the objective function. In our IA, the elementary internal vector is given by  e Fint ˆ ‰Bm ŠT frgkz h0 Ae ; …58† where ‰Bm Š is the constant membrane strain matrix, frg is the Cauchy stress vector, kz is the thickness stretch and Ae is the area of the triangular element.

144

Y.Q. Guo et al. / Computers and Structures 78 (2000) 133±148

The di€erentiation of Eq. (58) with respect to xj gives       e oFint oBm or ˆ r ‡ ‰Bm ŠT kz h0 Ae oxj oxj oxj  T   okz e oAe ‡ Bm frgh0 A ‡ kz : …59† oxj oxj   Terms oBm =oxj and oAe =oxj can be calculated easily for the triangular CST membrane element. In Eq. (59), xi represents the control node coordinates where  we are looking for the sensitivity. Terms okz =oxj and or=oxj are more complicated and can be obtained from di€erentiation of the IA equations. In the IA, stresses are given by (36) 

frg ˆ ES ‰ P Šÿ1 ‰C Šfeg;

…60†

where [C] is the transformation matrix from local to principal directions [45,53,55,56]. The di€erentiation of (Eq. (60)) with respect to xj gives      or oES oC ‰C Šfeg ‡ ES ˆ ‰ P Šÿ1 feg oxj oxj oxj   oe ‡ ES ‰C Š : …61† oxj   The computation of oC=oxj is straightforward. The term oES =oxj can be obtained by di€erentiation of the relation ES ˆ r=e, and using the explicit uniaxial strain± stress curve [45,53]. The last vector oe=oxj can be performed by di€erentiation of the eigenvalues of the left Cauchy±Green tensor [B] expressed in terms of the nodal displacements [45]. The use of the adjoint variable method allows one to solve Eq. (55) only one time at each loop of the optimization process, but before solving this linear system of equations, we have to store the global tangent matrix at the end of the IA analysis and we need to calculate the vector f oJ =oU g. For example, for the objective function given by Eq. (49), its di€erentiation with respect to the global displacement vector {U} is reduced to the computation of fokz =oU g. This can be done by di€erentiating the incompressibility relation: kz ˆ

1 : k1 k2

Fig. 10. Finite element mesh of the dashpot cup of Twingo (1634 elements).

surface is ¯at but inclined. The ®nite element mesh consists of 1634 triangular membrane elements (Fig. 10). The material characteristics are E ˆ 206 800 MPa, 0:176 MPa,  r ˆ 0:995, h0 ˆ 1:97 mm. r ˆ 624…0:007 ‡ ep † The friction e€ects are not considered in the optimization procedure. The contour of the workpiece is described by a Bspline curve. Eight design variables are the curve control knot positions. They are constrained to move along the normal directions of the contour of the useful part. These directions and the lower and upper bounds of the contours are illustrated in Fig. 11. One of the industrial objectives in optimization is also to minimize the amount of material. For the initial conditions of the optimization loop we consider the lower bound contour (where the volume is minimum) as the preliminary solution for the design variables fv0 g ˆ fvmin g. The actual TWINGO dashpot cup is produced with a rectangular blank (Fig. 10). The optimization results will

…62†

The technique exposed above, has been implemented in our code. Very good results are obtained in point of the gradient accuracy and computer time. 8. Numerical application of the optimization procedure We present in this section, the optimization results obtained for the dashpot cup of the Renault Twingo [14]. It presents a signi®cant sheet thinning, specially in the zone under the punch ()20%). The blank holder

Fig. 11. Parametrization of the contour.

Y.Q. Guo et al. / Computers and Structures 78 (2000) 133±148

145

Table 6 Optimization results

Rect. blank (reference) Lower bound Optimization results pˆ1 pˆ2 pˆ4 pˆ8

ÿDhmax …%†

‡Dhmax …%†

Volume (cm3 )

)20.27

+17.44

314

)13.28

+17.87

256

)18.85 )16.55 )13.03 )12.92

+13.86 +13.09 +13.56 +13.41

311 269 263 261

be compared with this reference solution. The average CPU time for the complete optimization process is about 2 h on an HP J280 (PA 8000) machine. The variations of thickness and the volume according to di€erent values of p are presented in Table 6. The evolution of the objective functions is shown in Fig. 12. Around 10±20 iterations are needed for the convergence of the optimization process. For p ˆ 1, we minimize the sum of the thickness variations on the whole workpiece for the most uniform distribution of the thickness. The exponent p ˆ 2; 4; 8; . . . is introduced in order to emphasize the reduction of the most important thickness changes. The results show that the values of the maximum and minimum of the thickness changes are minimized compared to the reference solution (rectangular blank). For p ˆ 8, a reduction of 16.6% of the material volume was obtained with respect to the reference solution. We

Fig. 13. Optimal workpiece contours obtained by using di€erent p values.

Fig. 14. Thickness distribution obtained with a rectangular blank (reference solution).

Fig. 12. Evolution of objective functions with di€erent p values.

Fig. 15. Thickness distribution obtained by optimization with p ˆ 8.

146

Y.Q. Guo et al. / Computers and Structures 78 (2000) 133±148

present in Fig. 13 the optimal contours of the ®nal workpiece obtained for di€erent p values. The optimal solution for p ˆ 8 is quite close to the lower bounds, but it gives smallest values of Dhmax …%†. The solutions for p ˆ 2; 4; 8 are very close. The thickness distributions for the reference solution and p ˆ 8 are given in Figs. 14 and 15. 9. Conclusions A simpli®ed inverse approach is developed to estimate the large elasto-plastic strains in thin metallic sheets obtained by deep drawing. The approach takes into account the a priori knowledge of the workpiece shape. It involves a discretization of the 3D surface by triangular facet shell elements and the computation of the inversedeformation gradient tensor to estimate the large logarithmic strains. The total deformation theory is considered to de®ne the constitutive equations. The nonlinear equilibrium equations result from the principle of virtual work for some assumed load vectors due to the tool actions on the sheet. The IA leads to only two dof per node, even if bending e€ects are taken into account. The analytical formulas of Stoughton are generalized to consider the drawbead restraining force and blank holding force. Some improvements on resolution algorithms make the IA code more ecient and robust. From our experience, the present formulation appears to be quite ecient in terms of CPU time and suciently precise for the preliminary design phase of blanks and tools. The eciency and convergence rate are found quite dependent upon the initial solution (the guess of the initial mesh of the blank for a given mesh of the workpiece), upon the nature of the stress±strain curve and also on the magnitude of plastic strains. The quality of the results compared to those obtained by incremental codes depends on the path dependent character of the deformation process (elasto-plasticity, multiple stage forming, friction). An automatic procedure combining the IA with an SQP algorithm has been developed to optimize the shape of the blank. Objective functions based on the thickness variations are considered. The parametrization of the contour limits the number of design variables. An analytical (explicit) calculation of the sensitivities leads to a€ordable computer times for the optimum design of blank contours. The numerical procedure has given satisfactory results for the industrial application (Twingo dashpot cup). The IA combined with a mathematical programming algorithm appears to be a valuable tool for the automatic design of 3D industrial metallic blanks subjected to deep drawing. Other process parameters (material

properties, draw bead geometry) can also be optimized using similar procedures.

References [1] Chenot JL, Wood RD, Zienkiewicz OC, editors. Numerical methods in industrial forming processes. Proceedings of the NUMIFORM 92 Conference, Rotterdam: Balkema, 1992. [2] Wang Y, Zhou S, Li C. International Conference on Advanced Technology and Machinery in Metal Forming. Huazhong, University of Science and Technology Press, 1992. [3] Teodosiu C, Raphanel JL, Sidoro€ F, editors. Large plastic deformations, fundamental aspects and applications to metal forming. Rotterdam: Balkema, 1993. [4] Makinouchi A, Nakamachi E, Onate E, Wagoner RH, editors. Numisheet'93. In: Second International Conference on Numerical Simulation of 3D Sheet Metal Forming Processes ± Veri®cation Simulation with Experimental, Isehara, Japan, 1993. [5] Wang ZR, He, Y. Advanced technology of plasticity. In: Proceedings of the Fourth International Conference on Technology of Plasticity. Beijing, China: International Academic Publishers, 1993. [6] Kroplin B, Luckey E, editors. Metal forming process simulation in industry. In: International Conference and Workshop, Baden±Baden, Germany, 1994. p. 28±30. [7] Altan T, editor. Sheet metal forming technology. J Mater Process Tech 1994;46(3±4). [8] Owen DRJ, Onate E, Hinton E. Computational plasticity ± fundamentals and applications, COMPLAS IV. Swansea, UK: Pineridge Press, 1995. [9] Shen SF, Dawson P, editors. Simulation of material processing: theory, methods and applications, NUMIFORM 95. Rotterdam, Brook®eld: Balkema, 1995. [10] Gerdeen JC, Chen P. Geometric mapping method of computer modeling of sheet metal forming, NUMIFORM 89, 1989. p. 437±44. [11] Sowerby R. Determination of large strains in metal forming. J Strain Anal 1982;7:95. [12] Chung K, Lee D. In: First International Conference on Technology of Plasticity. Computer-aided analysis of sheet material forming processes, vol. 1. Tokyo, Japan, 1984. p. 660±5. [13] Sklad MP, Yungblud BA. Analysis of multi-operation sheet forming processes, NUMIFORM 92, 1992. p. 543±7. [14] El Mouatassim M, Thomas B, Jameux JP, Di Pasquale E. An industrial ®nite element code for one step simulation of sheet metal forming, NUMIFORM 95. [9] 1995. p. 761±6. [15] Liu SD, Kolodziejski J, Assempoor A, Aboutour T, Cheng W. Development of a fast design and trouble-shooting FEM in sheet metal forming. 19th IDDRG Biennal Congress, 1996. p. 265±76. [16] Liu SD, Assempoor A. Development of FAST 3D ± a design-oriented one step FEM in sheet metal forming, COMPLAS IV, Part II, 1995. p. 1515±26. [17] Chung K, Richmond O. Sheet forming process design on ideal forming theory, NUMIFORM 92, 1992. p. 455±60.

Y.Q. Guo et al. / Computers and Structures 78 (2000) 133±148 [18] Batoz JL, Duroux P, Guo YQ, Detraux JM. An ecient algorithm to estimate the large strains in deep drawing, NUMIFORM 89, 1989. p. 383±8. [19] Guo YQ, Batoz JL, Detraux JM, Duroux P. Finite element procedures for strain estimations of sheet metal forming parts. Int J Num Meth Engng 1990;39:1385±401. [20] Batoz JL, Guo YQ, Detraux JM. An inverse ®nite element procedure to estimate the large plastic strain in sheet metal forming. In: Third International Conference on Technology of Plasticity. Kyoto, 1990. p. 1403±8. [21] Guo YQ, Batoz JL, El Mouatassim M, Detraux JM. On the estimation of thickness strains in thin car panels by the inverse approach, NUMIFORM 92. Rotterdam, Brook®eld: Balkema. 1992. p. 1403±8. [22] Guo YQ, Batoz JL, El Mouatassim M, Detraux JM. Determination of the trimming contour under the curved blank holders with the inverse approach. In: Proceedings of the International Conference, ATMMF'92. Wuhan, China, 22±24 October, 1992. [23] Batoz JL, Narainen R, Duroux P, Guo YQ. Comparison of the implicit, explicit and inverse approaches for the estimation of the strain distribution in axi-symmetrical thin sheets. In: Fourth Conference on the Technology of Plasticity. Beijing, China, 1993. p. 1695±700. [24] Chateau X. A simpli®ed approach for sheet forming processes design. Int J Mech Sci 1994;36(6):579±97. [25] Duroux P. ISOPUNCH, Un outil d'aide  a l'analyse de l'emboutissabilite des pieces de carrosserie,  a la disposition des bureaux d'etudes et des methodes. Congres International SIA, Lyon, 23±24 Mars 1994. p. 5. [26] Ravier P. Aide a la conception des pieces par des methodes de calcul simpli®e, CETIM (Calculs en emboutissage), 1994. p. 67±79. [27] Wang NM, Tang SC. Analysis of bending e€ects in sheet forming operations. Int J Num Meth Engng 1988;25(1): 253±67. [28] Peric D, Owen DRJ, Honnor ME. Simulation of thin sheet metal forming processes employing a thin shell element. FE Simulation of 3-D Sheet Metal Forming Processes in Automotive Industry. Z urich: VDI Verlag, 1991. p. 569±600. [29] Lee D. Recent innovations in sheet material forming. J Mater Process Technol 1994;46(3±4):333±49. [30] Brunet M, Sabourin F. A simpli®ed triangular shell element with necking criterion for 3D sheet forming analysis, NUMISHEET 93. [9] Rotterdam, Brook®eld: Balkema. 1993. p. 229±38. [31] Batoz JL, Dhatt G. Modelisation des structures par elements ®nis. vol. 3. Coques, Hermes Editeur, Paris, 1992. [32] Roelandt JM, Batoz JL. Shell ®nite element for deep drawing problems: computational aspects and results. In: Besdo D, Stein E, editors. Finite inelastic deformationstheory and applications, IUTAM Symposium. Berlin: Springer. 1992. p. 8. [33] Roelandt JM, Liu XJ, Batoz JL, Jameux JP. Axi-symmetric and general shell elements for large transformations (formulation and applications). Swanesa, UK: Pineridge Press, 1993. p. 457±63. [34] Onate E, Zarate F, Flores F. A simple triangular element for thick and thin plate and shell analysis. IJNME 1994;37(15):2569±82.

147

[35] Batoz JL, Guo YQ, Mercier F. Accounting for bending e€ects in sheet metal forming using the inverse approach, COMPLAS IV, Part I, 1995. p. 707±18. [36] Batoz JL, Guo YQ, Mercier F. The inverse approach including bending e€ects for the analysis and design of sheet metal forming parts, NUMIFORM 95, 1995. p. 661±7. [37] Hosford WF, Caddell RM. Metal forming. Mechanics and metallurgy. Englewood Cli€s, NJ: Prentice Hall. 1983. [38] Lee JK, Kinzel GL, Wagoner R. Numerical simulation of 3D sheet metal forming processes, veri®cation of simulations with experiments. In: Third International Conference on NUMISHEET 96. The Ohio State University, Dearborn, MI, September 29±October 3, 1996. [39] Stoughton TB. Model of drawbead forces in sheet metal forming. In: Proceedings of the 15th Bienal IDDRG Congress. Dearborn, 1988. p. 205±15. [40] Hu XB, Guo YQ, Batoz JL. Inverse approach for sheet metal forming using explicit, implicit static and DR algorithms. In: Fifth ICTP. Columbus, OH, USA, October 1996. [41] Roll K. Anforderungen der praxis an programme zur umforsimulation. J strain anal 1994;1:166±82. [42] Batoz JL, Guo YQ. Analysis and design of sheet forming parts using a simpli®ed inverse approach, COMPLAS V. Barcelona, Spain, 1997. p. 178±95. [43] Batoz JL, Guo YQ, Mercier F. The inverse approach with simple triangular shell elements for large strain predictions of sheet metal forming parts. J Engng Comput 1998;15(7): 864±92. [44] Altan T. Advance technology of plasticity. In: Fifth ICTP. Columbus, OH. October 1996. [45] Duroux P, de Smet G, Batoz JL. In¯uence of the blank contour on the thickness variation in sheet metal forming. In: Predeleanu M, editor. Computional methods for prediction material processing defects. Amsterdam: Elsevier, 1989. p. 103±12. [46] Toh CH, Kobayashi S. Deformation analysis and blank design in square cup drawing. Int J Machine Tool Design Res 1985;25:15. [47] Tatenami T, Nakamura Y, Ohata T. E€ect of pro®le of blankholder surface on drawbility in cylindrical deep drawing process. In: Third ICTP, Kyoto, 1990. p. 1237±42. [48] Carleer B. Equivalent drawbead model in ®nite element simulations, NUMISHEET' 96, 1996. p. 25±31. [49] Levy BS. Development of a predictive model for draw bead restraining force utilizing work of Nine and Wang. J Appl Metal Working 1983;3(1):38±44. [50] Nine HD. Drawbead forces in sheet metal forming. In: Koisinen DP, Wang NM, editors. Mechanics of sheet metal forming. New York: Plenum Press, 1978. p. 179±211. [51] Yellup JM. The prediction of strip shape and restraining force for shallow drawbead systems. J Appl Metal 1985; 4:30±8. [52] Wang NM. A mathematical model of drawbead forces in sheet metal forming. J Appl Metal working 1982;2:193±9. [53] Batoz JL, Naceur H, Barlet O, Knopf-Lenoir C. Optimum design of blank contours in axisymmetrical deep drawing process. In: Wan-Xie Z, Geng-Dong C, Xi-Kui Li, editors. Advances in computational mechanics. Beijing, China: International Academic Publisher, 1996.

148

Y.Q. Guo et al. / Computers and Structures 78 (2000) 133±148

[54] Naceur H, Batoz JL, Barlet O, Knopf-Lenoir C. Optimisation de forme du contour du ¯an dans l'emboutissage de t^ oles axisymetriques. Dans Actes du 13eme Congres Francßais de Mecanique, 1±5 Septembre, Poitiers, France, 1997. [55] Barlet O, Batoz JL, Guo YQ, Mercier F, Naceur H, Knopf-Lenoir, C. The inverse approach and mathematical programming techniques for optimum design of sheet forming parts. In: Third Biennial European Joint Conference on Engineering Systems Design and Analysis, ESDA'96, July 1±4. Montpellier, France, 1996. [56] Barlet O, Batoz JL, Guo YQ, Mercier F, Naceur H, Knopf-Lenoir C. Optimum design of blank contours using the inverse approach and mathematical programming

techniques. In: Third International Conference On Numerical Simulation of 3D Sheet Forming Process, NUMISHEETÕ96, September 29±October 3. Dearborn, MI, USA, 1996. [57] Kobayashi S, et al. Metal forming and the ®nite element method. Oxford: Oxford university Press, 1989. [58] Powell MJD. A fast algorithm for nonlinearly constrained optimization calculations. In: Watson GA, editor. Proceedings of the Biennial Conference Held at Dundee, June 1977. Lecture Notes in Mathematics, Numerical analysis, vol. 630. Berlin: Springer, 1978. [59] Ryu YS, Haririan M, Wu CC, Arora JS. Structural design sensitivity analysis of nonlinear response. Comput Struct 1985;21(1/2):245±55.