Three dimensional forging simulation with finite element control

Journalof

ELSEVIER

J. Mater. Process. Tcchnol. 45 (1994) 75-80

Materials Processing Technology

T h r e e d i m e n s i o n a l forging simulation with finite e l e m e n t control Jan Terziyski~, Tatsuhiko Aizawa~ and Junji Kihara ~ ~Department of Metallurgy, The University of Tokyo, Hongo 7-3-1, Bunkyo-ku, 113 Tokyo, Japan

The three dimensional rigid-plastic code NS2EC (Nonlinear Solid Analysis by Element Control) with new ideas has been developed and utilized to forging simulation. New robust finite element scheme has been used to describe the materials incompressiblebehaviors and to trace large geometrical changes. The mixed type variational principle has been employed with penalty function method. The finite element control by direct tessellation algorithm has been used for generation of the initial model and has been proved promising for adaptive remeshing and Arbitrary Lagrangian-Eulerian geometry control.

1. I N T R O D U C T I O N Two analytical approaches are widely used in metal forming simulations: Eulerian and LagrangJan. The Eulerian formulation, where the the mesh or grid frame is fixed in space and the material particles flow through elements, is well adaptive to analyses of liquid-like mechanical behaviors and rapidly deforming materials. In the Lagrangian approach, the mesh model varies with the progress of simulation, so that the geometric configuration is changed through deformation. For engineering solids, the total Lagrangian method is preferable, where the current stress state can be obtained by means of integration from the initial configuration. When considering strain-path dependent materials, the incremental procedure in terms of the updated Lagrangian technique can precisely describe the visco/elastoplastic and visco/rigid-plastic mechanical behaviors. Among the recent advances in the finite element methodology, new approaches are progressing for essential innovation of metal forming analyses. As pointed out in Ref[1], the mesh adaptivity with mapping/remapping features is indispensable for the practical simulations without element distortions; the Eulerian approach incorporating the Lagrangian step[2] can provide successful simulation up to failure or fracture; the Ar-

bitrary Lagrangian-Eulerian(ALE) method[3] is needed to trace the varying tool-work boundary in order to consider directly the tool-work interaction. In this paper, the three dimensional solid analysis code NS/EC is introduced to discuss its capabilities to deal with forging simulations. Firstly, the three dimensional composite finite element is proposed to trace large deformations and to consider the material incompressibility. Next, the mixed type variational principle is formulated by using the penalty function method, and the updated Lagrangian approach is used for time marching scheme. The finite element control by direct tessellation method is shown to be an indispensable tool for automatic mesh generation and adaptation during forging simulations; the ALE formulation is derived in the frame of the updated Lagrangian method to describe the singular velocity field in the vicinity of die corners. Finally, several numerical examples are taken to verify the validity of our developing code and the robustness of the composite finite element.

2. R I G I D - P L A S T I C A N A L Y S I S 2.1. C o m p o s i t e t e t r a h e d r a l e l e m e n t When considering actual processes together with material incompressibility, a reliable finite

0924-0136/94/$07.00 © 1994 - Elsevier Science B.V. All rights reserved. SSD1 0924-0136(94)00163-4

76

element scheme must be used to prevent severe mesh locklng[4,5] and to keep numerical stability in calculations. In order to consider the incompressibility condition, both velocity and pressure fields must satisfy the adequate requirements. Since the pressure might be in the class of 12, the element piece-wise constant approximation is chosen to trace its change. Through the theoretical patch test[5], the three dimensional composite element(Fig.l), where the velocity is approximated by piece-wise continuous function and the pressure is element-wise constant, has been found to have sufficient capability to consider incompressibility by enough degrees of freedom to satisfy the LBB condition[4,5].

the velocity field. When a tends to zero, the current problem described by equations (1) and (2) converges to an exact incompressible rigid-plastic problem. The pressure in each composite element can be expressed by: P = -JV[

,_m

= ~ ~' ~ ,

2.2. Mixed type variational equation Consider a solid medium of the volume V in rigid-plastic deformation. When the internal forces are absent and the mechanical and geometrical boundary conditions are prescribed on the boundaries Sa and S=, respectively, the mixed type variational equations can be written by:

= 0,

(4)

where t~, n and m are experimentally determined material constants and the corresponding multiaxial relation is deduced by the Levy-Mises condition. By substituting equations (2), (3) and (4) into Eq(1), the mixed type variational formulation with penalty function for the pressure field for the rigid-plastic solid is obtained:

~vdV

6ivdV =

Ti6uidSo

(5)

¢r

Fig.l: Three dimensional composite element

v(~V - ap)@dV

(3)

where Ve is the element volume. The rigid-plastic material constitutive equation is given by:

ate 0 Pressure node(free) 0 Velocity node(free) 0 Velocity node(constrained)

/vSiJ~i,jdV + j[vp~,,dV = fs Ti~u, dS ,

~vdV,

(1)

(2)

where Sij is the deviatoric stress tensor, tij the strain rate tensor, iv the volumetric strain rate, p the hydrostatic pressure, a a small positive penalty number, T~ the prescribed traction and u~

The iterative Newton algorithm with linesearch step should determine the convergent solution velocity field for the linearized version of Eq(5). The pressure, equivalent strain rate and stress tensor are readily obtained from the Eq(3), the Levy-Mises equation and the constitutive Eq(4). To deal with rigid inclusions and dead zones, where no material displacements can be observed, small tolerance for the equivalent plastic strain rate by ~ = 10-s is allowed.

3. F I N I T E E L E M E N T C O N T R O L In actual processing simulations, two main issues in the finite element modeling must be taken into account: 1) automatic mesh generation algorithm with capabilities for remeshing with mapping/remapping material history, and 2) ALE features to account for the tool-work interaction directly. However, many of the present codes for

77

forging analysis have left those tasks to semiautomatic procedures, which inevitably produces errors and irregularities due to human participation. 3.1. D i r e c t tessellation m e t h o d The full-automatic direct tessellation method with Delaunay algorithm has been developed and employed in metal forming analyses[4,6]. In Delaunay algorithm, the tetrahedral elements are connected from nodes ~i sharing common edges and intermediate faces. The mesh represented by such tetrahedrons is dual to the Voronoi polyhedron representation of the analytical region[6]. The initial Delaunay tessellation (DT) [~i, ~j, ~k, ~l] is created by series of minimizations for the edge distance [~i, ~j], the surface [~i, ~j, ~k] and the tetrahedron [~i, ~j, ~k, ~l]. The successive DT's can be obtained directly from the first one by sequence of eliminating loops for the control intermediate faces[~j, ]k, ~l], sharing edge in common, the control edges [~i, ~j], sharing node in common, and control nodes ~i, as shown in Fig.2. ml

? ,,~,, ! m

J

I

/ l~

! ~ z

~.

Con~'oJnode: #i

the convergent characteristics of the solution and to reduce numerical errors due to element distortions. 3.2. A d a p t i v e r e m e s h i n g b y D T m e t h o d In the frame of the direct tessellation method, the global and the local mesh refinement are available, where the new nodes generation should be based on the numerical error distribution over the whole domain or in the subdomain, respectively. The numerical error measure can be defined by 1) the geometrical changes in aspect the ratio or volume/surface ratio for the elements, and 2) purely physical criteria such as the velocity norm, equivalent strain deviation is elements and so forth. For the case of solid analysis, several error estimates must be examined to determine exactly the magnitude of numerical errors and to perform adequate regeneration of nodes[l]. The direct tessellation method can automatically reconnect the newly generated nodes preserving the domain contiguity and with little geometrical contaminations on the boundary between the initial mesh and the new one. Old tessellation

Old Voronol polyhedron PJ

Control edge: [#i, #j] k Control face: [#j, #k, #1] Initial DT: [#i, #j, #k, #1]

Bestcandidate

for connection:

#m eshlr~

Fig.2: Determination of new DT In order to preserve the contiguity of the domain without gaps and overlaps between elements and to speed up the mesh subdivision, the candidate nodes for each one DT must be selected among the set of approximate neighbors of intermediate face nodes ~j, ~k and ~1. It is mathematically proved that the elements generated by the direct tessellation method are semi-regular tetrahedrons and, in the case of orthogonal net modeling in the region, the aspect ratios, defined by the radii of circumscribed and inscribed spheres, are often located between 4.0 and 6.0 and seldom above 10.016]. This regularity is of primary importance for the finite element analysis to improve

ew te

per =-~-Zc (#ki)

Fig.3: Mapping/remapping procedure Once the remeshing is completed, the material history variables must be mapped to the new spatial frame in order to preserve the continuity of the solution in time and to provide an initial approximation for the incremental procedure. The material variables are transferred from the old frame to the new one in three steps(Fig.3): 1) the time pass dependent variables such as the equivalent plastic strain rate are averaged on the nodes of the old model; 2) the node-associated

78

variables must be transferred from the old nodes to the newly generated ones; at this stage, the I N / O U T check can determine the old Voronoi polyhedron encompassing the new node, and the mapping from an old node to a new node is oneto-one; 3) the history variables are interpolated from the new nodes to the new elements.

fv s~JbeijdV+ ~-~'~ffv evdV /y 6~jdV+

3.3. A L E g e o m e t r y control When selecting an analytical method for three dimensional simulation with consideration of toolwork interaction even in the cases of simply shaped dies, the nodal positions in space must be controlled to describe accurately the changing tool-work boundary and the singularities in velocity field in the vicinity of die corners[3,7]. The ALE method must be a powerful tool to control element nodes in order that the grid points should move in spatial(laboratory) coordinate system with a velocity different from this of the material particles. Consider a continuum of the volume V, where the coordinates of material particles are xi at time to. Under the external tractions, the particles Px move to new positions Px at time t = to + At, and we denote the corresponding material domain by Vx. For a reference coordinate system V×, moving in space with an arbitrary velocity, the particles P~ can be described by their referential coordinates Xi. The mapping J = [x~/xj[ is referred as the Jacobian transformation between the referential and the spatiM coordinates. For a physical quantity f(:e, t), we can define a corresponding mapping function ](x,t) such, that f(x, t) = f(~e,t). It can be shown[7], that

o~Se

O (Jf(x,t))

=

J[O (](x,t) + v ( / ( x , t ) ~ ) ]

(6)

where w is the relative material velocity with respect to the reference frame. Eq(6) provides the time derivative of physical quantities in the reference coordinates identical to that in the spatial coordinates. Taking together equations (1), (2), (3) and (6), after applying Gauss divergence theorem, we obtain the rigid-plastic ALE variational formulation with penalty function method:

~Q SijwlnjS~ijdS+ Q(¢)

JSQ

where SQ is the boundary on which ALE formulation is used and nj is the normal vector to this boundary. The Eq(7) can be viewed as the sum of the conventional Lagrangian Eq(5) and the energy of the transported material through the boundary SQ due to the referential frame movement. In practical applications, only the nodes in the vicinity of the contact surfaces can be treated by ALE description. The velocity field for such nodes should have two vectors: the material and the grid velocities. Since the velocities in cosine directions to the tool surface are prescribed for those nodes and the mesh movement is usually known, the convective terms for those degrees of freedom are not calculated. Hence, the surface integrals over SQ must be taken only in regions, where 1) no boundary conditions are applied and 2) the ALE description is used. It is to be noted that, to generate an undistorted reference mesh, the direct tessellation method should be available for automatic reconnection of ALE nodes at each step of reference movement.

4. N U M E R I C A L

RESULTS

4.1. Upsetting of cylindrical billet To verify our developing simulator, let us solve the incompressible problem depicted in Fig.4a. The velocity of the tool in contact with the upper surface is prescribed, so that the reduction in height should be 1% at each time step. No friction is prescribed on the tool-work boundary and the material constants in Eq(4) are chosen ~ = 10 MPa and n = 0, and m is varied from 0.1 to 1.0. The initial configuration and the deformed

79

shape at reduction in height of 50% are shown in Fig.4b and Fig.4c, respectively. The calculated results show good agreement with the exact solution, and convergent solutions for the the velocity field is obtained in two iterations for each time pass. The number of iterations is found to be indifferent to the values of the penalty number a and the volume change is less than 0.01% for a --: 10 -14.

be strongly dependent on the lubrication condition between the punch and the material and on the initial work diameter. The deformed shapes and the pressure distributions at the reductions in height of 20 % and 40% are shown in Fig.6a and Fig.6b, respectively.

a) a)

b)

b)

Fig.5: a) Impression die forging; b) Initial 1/8-modeh 595 nodes,288 elements

c) Fig.4: a) Experimental setup; b)Initial configuration: 225 nodes, 96 elements; c) Deformed shape at the reduction in height of 50% 4.2. I m p r e s s i o n - d i e forging The initial positions and the geometries of the tool, the dies and the material are shown in Fig.Sa. The axisymmetric punch is pushed down towards the material with a prescribed velocity. The outer die and the inner core constrain the material flow in radial direction. Assuming the field symmetry, 1/8-model of the original workpiece is discretized into the mesh grid as shown in Fig.5b. It was observed that, beyond 40% of reduction in height, the deformed material contacts with the inner die core, which yields severe strain localization in its central region and possible damage. The impression filling and the flash were found to

a)

b)

Fig.6: Deformed shapes and pressure distributions at the reductions in height of a) 20% and b) 40% 4.3. Forward e x t r u s i o n The capability of the three dimensional composite element to deal with incompressible forging analyses has been proved by forward extrusion simulation, the initial setup of which is illustrated in Fig.7a. On the upper work surface, the velocity of the tool is prescribed to achieve reduction in height of 1% at each time step. The original material shape is a parallelepiped(Fig.7b), and

80 the square die consists of several supporting side walls, as shown in the figure.

i a)

b)

Fig.7: a) Forward extrusion; b) Initial 1/4-modeh 343 nodes,162 elements

a)

5. C O N C L U S I O N The three dimensional rigid-plastic analysis code with element control NS~C has been developed and utilized to several forging analyses. Through the numerical examples, the composite tetrahedral finite element and the computational procedure has been proved to be robust against severely distorted geometries and highly nonlinear material behaviors. The pressure field is smoothed over the analytical domains, which certifies the composite tetrahedron for incompressible analyses. The direct tessellation method has been used for finite element control on the stages of initial modeling and has been found favorable for remeshing and ALE control. However, to improve the convergent characteristics of the numerical solution, the mesh adaptivity scheme by the direct tessellation method with mapping/remapping functions should be utilized to trace large deformations and to deal with strain localization zones. The ALE formulation for rigidplasticity together with element control can describe precisely the changeable tool-work boundary, and its application to metalworking analysis is to be reported.

b)

Fig.8: Deformed shapes and pressure distributions at the reductions in height of a) 9% and b) 27% The material flows inside the die and the pressure gradients at reductions in height of 9% and 27% are shown in Fig.8a and Fig.8b, respectively. Although the nodes C, located at the edge-corners of the sidewalls, were fixed in order to account for the changing tool-work boundary, there was no mesh locking induced by the incompressibility constraints[5]. The time step predictor-corrector scheme exhibited good approximation of the die shape and the numerical behavior was stable with only six to seven iterations for the velocity field solution in each time increment. However, the numerical solution could not converge beyond 36% of reduction in height due to severe element distortions in the model front tip.

REFERENCES 1. G.Lin et al.: Numerical Methods for Simulation of Industrial Metal Forming Processes, ASME, AMD-156(1992), 33. 2. D.Benson: Comp. Meth. Appl. Mech. Eng, 99(1992), 235. 3. T.Aizawaand J.Kihara: IUTAMSymposium, Hannover, 1991, 449. 4. T.Aizawa and J.Kihara: J. Mat. Proc. Tech., 27(1991), 55. 5. O.Zienkiewicz et al.: Int. J. Num. Meth. Eng., 23(1986), 1873. 6. T.Aizawa and Y.Yoshino: Proceedings of tile Postconference seminar, ICES, 1989, 29. 7. J.Donea: Arbitrary Lagrangian-Eulerian Methods, in T.Belytchko and T.Hughes: Computational Methods for Transient Analysis, 1983, Chapter 10.