2006:068 CIV
M A S T E R ' S T H E S I S
Finite Element Simulation of Punching
Magnus Söderberg
Luleå University of Technology MSc Programmes in Engineering
Mechanical Engineering
Abstract
This work is a study of a punching process using the ABAQUS/Explicit FE-code. A damage model based on equivalent plastic strain has been used to describe crack initiation and propagation in the blank.
An axisymmetric model created in ABAQUS/CAE was used throughout the study. The model consists of a punch, blank holder, die and blank. In this work the punch, blank holder and die are modelled as rigid bodies while the blank is considered elastic-plastic.
The materials used in simulations are DC04, Docol 350 YP, Docol 800DP and DC1400M.
This selection of steel qualities has a wide span in mechanical properties and fields of application.
Material characterisation for simulation purposes has been done by conducting uniaxial tensile tests. The data from these tests have been extrapolated beyond the necking strain by double Voce extrapolation.
Punching experiments have been performed in an excenter driven press equipped with position and force measurement devices. The experiments result in load curves for the different steel qualities and gives information about the characteristic zones of the sheet edges.
The results from simulations have been evaluated by comparisons with experiments concerning load curves, characteristic zones of the sheet edge, work done during the operation, maximum forces and stroke to failure.
A good compliance with experiments was found concerning maximum forces, characteristic zones and work done during the operation.
Table of Contents
1. Introduction...4
2. Influence of simulation parameters...5
2.1 FE-Model ...5
2.1.1 Geometry...5
2.1.2 Element type ...6
2.1.3 Material data and model...6
2.1.4 Contact and boundary conditions ...7
2.1.5 Failure criterion...8
2.2 Fe-Simulations ...9
2.2.1 Mesh density ...9
2.2.2 Adaptive mesh ...14
2.2.3 Failure strain value...16
2.3 Comparison with results from DEFORM ...21
2.3.1 Failure criterion in DEFORM...21
2.3.2 Load curves...21
2.3.3 Simulated edge geometry...22
2.4 Elastic tools...23
2.5 Comments ...27
3. Material properties ...28
3.1 Materials ...28
3.2 Tensile tests...28
3.2.1 Mechanical properties...28
3.2.2 Extrapolation of test data ...29
4. Experiments ...30
4.1 Experimental setup ...30
4.2 Data from experiments...31
5. Simulation model...32
5.1 Geometry and mesh density...32
5.2 Material model ...33
5.2.1 Strain rate dependency...34
5.3 Failure criterion...34
5.4 Contact and boundary conditions ...35
5.4.1 Contact conditions ...35
5.4.2 Boundary conditions ...37
6. Comparison with experiments ...38
6.1 Method of evaluating simulations...38
6.1.1 Evaluation of simulated characteristic edge zones ...38
6.1.2 Evaluation of work...39
6.1.3 Evaluation of stroke to failure...39
6.1.4 Evaluation of maximum punch force...40
6.2 Load curves and edge profiles ...40
6.2.1 DC04...40
6.2.2 Docol 350 YP...41
6.2.3 Docol 800 DP...42
6.2.4 Docol 1400M ...43
6.3.3 Fracture zone...46
6.4 Maximum punch force...47
6.5 Work ...48
6.6 Stroke to failure...49
7. Discussion...50
8. Conclusions...51
9. Future work...52
10. Acknowledgements...53
11. References...54
Appendix 1 Mechanical properties ...55
Appendix 2 Data from punching experiments ...57
1. Introduction
Punching is a sheet metal forming procedure were material shearing mechanisms are used to produce the desired geometry. The development today is focused against using
advanced high strength steel qualities to make lighter products with sustained or increased product durability. These steel qualities puts higher demands on the tooling used to cut or form the sheet steel, which makes the choice of process parameters even more important than with ordinary mild steels, because of the narrow process window when using high strength steels. Identification of the process parameters is today mostly done by large series of experiments, which is costly and time-consuming [1]. In earlier work performed by D.
Thunvik [2] and others [3], [4], the DEFORM FE-code has been used to simulate the punching process with results close to experimental. Simulations with the
ABAQUS/Standard FE-code and damage implemented by means of a user subroutine have also shown results in compliance with experiments [5], [6]. This shows the possibility of using FE-simulations to reduce the number of experiments that has to be conducted and give a increased understanding concerning the influence of process parameters.
This work utilises the ABAQUS/Explicit FE-code and a failure criterion that is readily implemented in the code. The objective has been to resolve the capacity of this code concerning simulations of punching and establishment of a modelling procedure for this type of problems. Simulations will be compared to experimental data concerning geometry of the edges, load curves, work done by the punch during the process, maximum forces involved and stroke to failure of the blank.
The first section is devoted to finding a suitable modelling approach in ABAQUS concerning mesh densities, failure criterion and features in the software such as adaptive mesh. The results from this section will be compared to simulations made in DEFORM [2].
2. Influence of simulation parameters
In this section a FE model is created in ABAQUS to resolve the influence of the variations in the mesh density, value of the failure criterion and the adaptive mesh function that is implemented in the software. A single case of punching in Docol 800 DP will be simulated using the dimensions and material model as earlier in the DEFORM [2] simulations, to make a comparison between the results from the two different software packages possible.
Simulations with elastic tools have been made to study if it is possible to evaluate the stress state in the tools during punching and the elastic response in the tools using this modelling approach.
2.1 FE-Model
All FE-models evaluated in this study have been created using the ABAQUS/CAE pre- processor. The solution is obtained with the ABAQUS/Explict solver, which uses a direct integration scheme on the structural dynamic equations. This approach is suitable for highly dynamic problems with non-linear behaviour, which is the case when simulating fracture of metals.
2.1.1 Geometry
The tools and blank have been modelled using an axisymmetrical model that can be seen in Figure 2.1, values of the dimensions can be found in Table 2.1.
Figure 2.1: The axisymmetric FE-model in ABAQUS.
Table 2.1: Dimensions corresponding to the model in Figure 2.1
r1 [mm] r2 [mm] r3 [mm] r4 [mm] c1 [mm] c2 [mm] t [mm]
2.44 0.01 0 0.01 0.5 0.06 1
2.1.2 Element type
The elements in the blank consist of four node bilinear axisymmetric quadrilateral elements with reduced integration (CAX4R) and three node axisymmetrical triangles (CAX3). The triangles are used to coarsen the mesh for increased computational
efficiency, see Figure 2.4. Both elements belongs to the family of solid elements and are of the first order, which means that the strain is computed as an average over the element volume instead of the first order gauss point. The feature of reduced integration used in the CAX4R element causes the integration order to be lower than full integration, in this case only one integration point in the centre of the element is used. By using reduced integration the numbers of constraints which are introduced by the elements are decreased, this
prevents “locking” in the elements causing a stiff response. The drawback of this technique is that for certain modes of deformation no energy is registered in the element integration point. These modes are usually referred to as “hourglass modes”. This problem is
addressed in ABAQUS using a “hourglass control” algorithm [9].
2.1.3 Material data and model
The blank material consists of Docol 800 DP in all simulations in this section. Yield data were taken from [2] in an effort to make the simulations as comparable as possible, see Figure 2.2. The yield data is extrapolated using the double Voce method, which is a way of producing reliable data at high strains. This procedure is explained in more detail in section 3.2.2, which covers the topic of extrapolation of yield curves.
Double Voce extrapolated yield data
0 200 400 600 800 1000 1200 1400
0 0.2 0.4 0.6 0.8 1
True plastic strain
True stress [MPa]
Yield stress Docol 800 DP
Figure 2.2: Yield data for Docol 800DP extrapolated with the double Voce method.
The material model used is the isotropic von Mises hardening model. In this model the material is assumed to have similar properties in all directions. As no reversed loading occur during the simulations the hardening was described as isotropic i.e. the Bauchinger effect is not modelled. The yield criterion for this model can be expressed by means of the principal stresses as:
] ) (
) (
) 2[(
1 2
3 2 2 2 1 2 2
1 σ σ σ σ σ
σ
σy = − + − + − (1)
be interpreted by saying that yielding occurs at a certain level of deviatoric strain energy [10].
2.1.4 Contact and boundary conditions
Contact between the tools and the blank is enforced by a kinematic contact condition, using pure master-slave surface pairs established in the first step of the solution. The master and slave designation must be chosen so the rigid tool forms the master surface and the surfaces defined on the blank act as slave [11].
The formulation of this contact condition is of predictor/corrector type were the model is first advanced to a kinematic state without consideration to the contact condition. The slave nodes witch penetrates the master surface is then determined and the force to move the slave nodes on to the master surface is calculated. These forces are calculated based on the depth of penetration, nodal masses and time increment. The acceleration correction for the master and slave surfaces is then calculated based on the forces needed to oppose
penetration and the inertia of the contacting bodies [11]. The surfaces which form the contact pairs in this model are summarised in Figure 2.3 and Table 2.2.
Figure 2.3: The surfaces used in the contact pairs. 1) is the surface of the punch, 2) is the top nodes of the blank, 3) is an internal nodal based surface, 4) is the surface of the blank holder, 5) is the bottom nodes of the blank and 6) is the surface of the die
Table 2.2: A summarisation of the contact conditions numbered as in Figure 2.3 Master surface Slave surface/surfaces
1 2, 3
4 2
6 3, 5
Friction between the surfaces is implemented with a Coulomb model defined as:
n
f µσ
τ = (2)
Where τf is the friction shear stress, µ is the friction coefficient and σn is the normal pressure. A friction coefficient of 0.1 has been used in all simulations.
The punch displacement is applied as a prescribed velocity of 40 mm/s and the holder force is set to be 1350 N which is about 15 percent of the maximum punch force. The nodes at the symmetry line are locked in the radial direction by a displacement boundary condition. When an axisymmetrical definition is adapted, no material flow in the angular direction is presumed.
2.1.5 Failure criterion
To estimate the start and propagation of fracture, a local fracture criterion is used in the simulations. The criterion used in this set of simulations is the shear failure model in ABAQUS. This model is based on a value of equivalent plastic strain at element
integration points and failure is assumed to occur when the damage parameter exceeds 1.
The damage parameter, ω, is defined as:
pl f
pl pl
ε ε
ω =ε0 +¦∆ (3)
Where ε0pl is the initial value of the equivalent plastic strain, ∆εpl is an increment of plastic strain and εfpl is the strain at failure [11].
2.2 Fe-Simulations
2.2.1 Mesh density
Because of the high strain gradients in the cut region, a sufficiently dense mesh must be applied in this area. The mesh density will also affect the geometry of the cut edges as a result of fracture being simulated with element deletion.
To resolve how the mesh density will affect the results four models with different mesh densities have been analysed. The number of elements was 656, 1488, 3901 and 7693 and the plastic failure strain was set to 1.5. Mesh densities were chosen by starting out with 16 through thickness (t in Figure 2.1) elements in the sheared zone and then dividing them by 2 until 128 through thickness elements were reached. The maximum number of elements through thickness will hereafter be used when referring to a mesh density i.e. 16, 32, 64 and 128. For additional data about the different mesh sizes see Table 2.3. The different mesh sizes are shown in Figure 2.4.
Table 2.3: Data about the different mesh densities Maximum number of through
thickness (t) elements Number of elements Smallest characteristic length of elements [mm].
16 656 0.0625
32 1488 0.03125
64 3901 0.01563
128 7693 0.007813
Figure 2.4: a) 656 elements, 16 through thickness elements, b) 1488 elements, 32 through thickness elements, c) 3901 elements, 64 through thickness elements and d) 7693 elements with 128 through thickness elements. The number of through thickness elements stated is the maximum number in the refined zone.
It can be seen from the load curves plotted in Figure 2.5 that the mesh density influences the rupture behaviour and stroke to failure to a greater extent than the maximum force.
This is believed to be an effect of the more localised plastic deformation when the mesh density is increased, which results in higher strains at an earlier stage of the simulation thus triggering the failure criterion with following element deletion, see Figure 2.6.
Mesh density Docol 800 DP
0 2000 4000 6000 8000 10000
0 0.2 0.4 0.6 0.8
Stroke [mm]
Force [N] 64 elements
32 elements 16 elements 128 elements
Figure 2.5: Force versus stroke curves using the above mentioned mesh densities. The plastic failure strain was set to 1.5 for all simulations.
Figure 2.6: Field plots of the equivalent plastic strain with a) 64 elements through thickness and b) 128 elements through thickness.
When coarser meshes are used the final rupture is delayed because of the cracks in the blank that originates from the punch and die miss each other, as shown in Figure 2.7. This behaviour is not observed when using 128 elements through thickness because of the possibility for cracks to initiate closer to the edges of the punch and die. Another effect of the increased mesh is that it seems to facilitate for the cracks to change columns in the mesh as they propagate, see Figure 2.8 for the simulated edge geometry using different mesh densities. The ability for the cracks to change column is of great importance if the crack tips shall be able to meet each other and produce rupture behaviour similar to the one observed in experiments.
Figure 2.7: Crack tips missing each other during fracture of the blank. The mesh density is 64 elements through thickness.
Figure 2.8: The simulated sheet edges with a) 16 through elements, b) 32 through thickness elements, c) 64 elements and d) 128 elements. The number of elements refers to the maximum number of elements through thickness as mentioned above A plastic failure stain of 1.5 was used in all simulations.
2.2.2 Adaptive mesh
Adaptive meshing in ABAQUS is carried out in two steps: creating a new mesh and remapping of the solution variables from the old mesh to the new mesh by an advection sweep. The adaptive meshing task is triggered at user specified intervals with a default value of 10 increments. During an adaptive mesh increment, the new mesh is created by one or many mesh sweeps, which moves nodes to reduce element distortion. The number of mesh sweeps can be specified by user input and has a default value of one. This value has to be increased if the deformation rate is high because the numerical stability of the advection sweep is maintained only if the difference between the new and old mesh is small [11]. For the simulations in this section, the default values will be used.
In the section above the crack tips miss each other when using coarse meshes, which causes a section of elements between the two cracks to experience large distortion before final rupture. This results in an increased stiffness causing a longer stroke to fracture. To investigate if the adaptive mesh function in ABAQUS is able to improve the solution, the four mesh densities used in the previous section were simulated using this function on the elements in the zone between the tool edges. The results of the simulations can be seen in Figure 2.9. The location of crack initiation is improved with adaptivity i.e. the cracks is starting nearer to the tool edges which is an effect of the mesh ability to follow the tool geometry closer. For 16 and 32 elements the cracks propagates in the same element column during the whole rupture process which gives a realistic load drop but fails to represent the edge geometry by producing a very short sheared zone, see Figure 2.10. With 64 elements the simulations produces results that are improved over the ones without adaptivity but still suffers from the fact that a element column is “trapped” between the two crack tips which results in a deviation from the simulation with 128 elements at a stroke of 0.16 mm. This leads to the conclusion that between 64 and 128 through thickness elements are needed for the simulations. For capturing the edge geometry, 128 elements will probably be required but 64 elements can be used to study the influence of parameter values on the force stroke curve to reduce the cpu time required to perform the simulations.
Mesh density Docol 800 DP
0 2000 4000 6000 8000 10000
0.00 0.10 0.20 0.30 0.40
Stroke [mm]
Force [N]
16 elements 32 elements 64 elements 128 elements
Figure 2.9: Force versus stroke curves for different mesh densities using adaptive meshing. The plastic failure strain was set to 1.5 for all simulations.
Figure 2.10: The simulated sheet edges with a) 16 elements, b) 32 elements, c) 64 elements and d) 128 elements in the thickness direction, using adaptive mesh.
2.2.3 Failure strain value
To investigate the influence of the failure strain value five models with the same mesh density, 64 through thickness elements with adaptive meshing has been simulated. The values of failure strain have been set to 1.5, 2.0, 3.0, 3.5 and 4. The results from the simulations showed that the failure strain value mainly affected the stroke to fracture and had a small influence on the maximum force, which can be seen in Figure 2.11.
Shear failure value Docol 800 DP
0 2000 4000 6000 8000 10000
0.00 0.10 0.20 0.30 0.40 0.50
Stroke [mm]
Force [N]
Shear 1.5 Shear 2.0 Shear 3.0 Shear 3.5 Shear 4.0 Experiment
Figure 2.11: Force versus stroke curves for different values of plastic failure strain. A mesh density of 64 through thickness elements with adaptive mesh was used in all simulations. The simulated curves are denoted with “Shear” and a corresponding value of the failure strain.
The influence of failure strain value on the characteristic zones have been investigated using a results from simulations in DEFORM[2] as a benchmark for a realistic edge geometry. This was done because no experimental data for the edge geometry could be found for this set of process parameters.
Measurments of the characteristic zones were taken by using the distance query tool in ABAQUS/CAE, which calculate the distance between nodes picked by the user. Defintion of the zones have been done by taking the rollover as the distance between the top of the blank and the first node at the transition to shear zone. The shear zone is measured between the transition point until initiation of fracture and the remaining part down to the lower edge of the blank is considered to be the fracture zone. Burr is measured from the lower edge of the blank. See Figure 2.12 for details concerning how the zones are measured.
Figure 2.12: a) Rollover is measured from the top of the sheet edge to the first node of the shear zone. b) The shear zone is measured from the end of rollover until the start of fracture. c) The fracture zone is measured from the start of fracture until the lower edge of the sheet. d) Burr is measured from the lower edge of the sheet.
The characteristic zones were measured using the procedure outlined above and the results are presented together with results from DEFORM for comparison, see Figure 2.13. The measured edge profiles are presented in Figure 2.14.
Characteristic zones
0 0.2 0.4 0.6 0.8 1 1.2
1 2 3 4 5 6
Length of characteristic zones [mm]
Rollover Shear zone Fracture zone Burr height
Figure 2.13: The characteristic zones measured with different values on the failure strain. Bar 1-5 are measurements from simulations done in ABAQUS with failure strain: 1) 1.5, 2) 2, 3) 3, 4) 3.5, 5) 4. Bar 6 is measurements from a reference case taken from DEFORM. All ABAQUS simulations
Figure 2.14: Simulated sheet edges with value of strain at failure of a) 1.5, b) 2.0, c) 3.0, d) 3.5 and e) 4.0. The edges were simulated with 64 through thickness elements and adaptive meshing.
values of 3.5 and 4 were used. The model with failure strain of 3.5 were simulated both with and without adaptive mesh but as can bee seen from Figure.2.15 the results with adaptive mesh were not as close to experimental results as the ones with a standard Lagragian mesh definition. Because of the erroneous results produced with adaptive meshing, it was discontinued for the simulations with a failure strain of 4, which produced a stroke to failure similar to experiments. The deterioration of the results with adaptive mesh can probably be explained by the large number of mappings of the solutions to a new mesh that has to be done when the failure strain value is set to high values, yielding an iteration count of about 5 million increments. Also some kind of incompatibility with adaptivity and failure combined is suspected. The simulated edge geometry for the models with 128 elements through thickness can be seen in Figure 2.16. Measurements of the characteristic zones are presented in Figure 2.17.
Docol 800 128 elements through thickness
0 2000 4000 6000 8000 10000 12000 14000 16000
0 0.1 0.2 0.3 0.4
Stroke [mm]
Force [N]
Failure strain = 3.5 Failure strain = 4.0 Experimental
Failure strain = 3.5 adap
Figure.2.15: Force versus stroke curves with 128 through thickness elements. “adap” indicates that adaptive mesh have been used.
Figure 2.16: The simulated sheet edges with 128 through thickness elements with a failure strain value of a) 3.5, b) 4 and c) 3.5 with adaptive mesh.
Characteristic zones
0 0.2 0.4 0.6 0.8 1 1.2
1 2 3 4
Length of characteristic zones [mm] Rollover
Shear zone Fracture zone Burr height
Figure 2.17: Measurements of the characteristic zones with different values of failure strain. Bar 1- 2 are measurements from ABAQUS using a failure strain of 1) 3.5 and 2) 4. Bar 3 is simulated in ABAQUS using a failure strain of 3.5 and adaptive mesh. Bar 4 is a reference case from DEFORM.
All simulations in ABAQUS used 128 elements through thickness.
The difference in the initial response in simulations and experiments can probably be explained by the fact that the experimental punching setup has elastic properties when loaded while the tools in the simulations are modelled as rigid. Another possible source of error is that the material model used in the simulations is unable to predict the actual behaviour of the material when subjected to this kind of deformation.
2.3 Comparison with results from DEFORM
DEFORM is an implicit code that is mainly used to simulate different kind of metal forming operations such as forging, machining and rolling. To facilitate these kind of simulations it has a more advanced mesh updating function were an adaptive refinement procedure is implemented in the code, which allows severe case of deformation without distorting the mesh. The code also has more models for fracture initiation and propagation than is implemented in ABAQUS today. A drawback with such a specialised code is the reduced capability to simulate more general problems were for example non-metal
materials is included. Another issue is that DEFORM is not as commonly used in industry as ABAQUS. Therefore a comparison between the results produced with the two different codes is of interest to see if the process can be simulated with no mesh updating and the less advanced failure criterion used in ABAQUS.
2.3.1 Failure criterion in DEFORM
The Cockroft and Latham failure criterion were chosen to predict failure initiation and propagation in the work used for comparison [2]. Other authors have also used this failure criterion with acceptable results [4]. The formula for the Cockroft and Latham failure criterion is given by:
C d
f
³σσ ε =
ε *
(4)
Whereσ is the maximum principal stress,σ is the effective stress, ε is the effective strain,* ε is the effective strain at fracture and C is the damage value which is a materialf
parameter. When the damage value C reaches a critical value in an element, given as user input, the element is deleted from the mesh. The interpretation of this criterion is that a high degree of triaxiality invokes fracture at a lower strain level compared to the uniaxial case.
2.3.2 Load curves
Two important results from the simulations are the maximum force and the stroke to failure. ABAQUS and DEFORM shows good agreement in predicted maximum force and stroke to failure compared to experiments. The noise produced in the ABAQUS curve is a result of mass scaling techniques being used to reduce the cpu time necessary to perform the simulations.
ABAQUS vs DEFORM and experiment
-2000 0 2000 4000 6000 8000 10000
0 0.1 0.2 0.3 0.4 0.5
Stroke [mm]
Force [N]
Experiment DEFORM C=1.5 ABAQUS Failure strain=4
Figure 2.18: Comparison between load curves from ABAQUS, DEFORM and experiments.
2.3.3 Simulated edge geometry
The dimensions of the different characteristic edge zones have been compared earlier and can be found in Figure 2.13 and Figure 2.17. The measured zones from the two codes differ somewhat but are in reasonably good agreement, keeping in mind that the
measurements are dependent on how the author chooses to define the limits of the different zones. The appearance of the edges is presented in Figure 2.19. The edges simulated in ABAQUS appear to be smoother than in DEFORM because of the higher degree of discretization necessary to capture the geometry and fracture behaviour when a adaptive mesh refinement method are not available.
Figure 2.19: Simulated sheet edges in DEFORM (left) and ABAQUS (right). The ABAQUS simulation has 128 elements through thickness with a plastic failure strain of 4,Adaptive meshing was not used.
2.4 Elastic tools
In order to resolve the difference between using elastic and rigid tools, a simulation model using elastic tools were created see Figure 2.20. The boundary conditions is applied in the same manner as the simulations using rigid tools except for the contact between the tools and the blank, which is modelled as a pure master-slave contact pair with element-node contact surfaces. This contact condition prescribes that nodes on slave surfaces cannot penetrate the master surfaces but nodes from master surfaces can penetrate slave surfaces.
In this case, the blank is modelled as a nodal based slave surface and the tools are element based master surfaces. The motivation for the choice of master and slave surfaces is that ABAQUS requires parts that can fail to be modelled with nodal surfaces and nodal surfaces is always slave using this contact definition [11].
Figure 2.20: The FE-model with elastic tools.
Figure 2.21 shows the difference between rigid and elastic tools using 128 through
thickness elements in the blank and a failure strain of 3.5. The elastic response in the tools can be seen as a slight displacement of the curve to the right in the beginning of the simulation and faster rupture of the blank. The faster rupture is assumed to be an effect of the elastic springback in the tools during load drop. A greater elastic response would probably be achieved by using tools that are longer in the 2 direction in Figure 2.20.
Elastic vs Rigid tools
-500 1500 3500 5500 7500 9500 11500
0 0.1 0.2 0.3
Stroke [mm]
Force [N]
Shear 3.5 128 el elast Shear 3.5 128 el
Figure 2.21: Force stroke curves with rigid and elastic tools. Both simulations had a plastic failure strain of 3.5 and 128 elements through thickness.
The stress in the tools is shown for different levels of stroke in Figure 2.22-19. The Von Mises stress, x-stress component (S11) and y-stress component (S22) is presented in the mentioned order, the orientation of the directions can be seen in figure 2.22 (1 is the x direction and 2 is the y direction). It can bee seen that the blank mesh is too coarse to capture the sharp corner radius of the punch. This leads to that only a few nodes from the blank are in contact with the punch creating very high local stress levels (over 6500 MPa in Figure 2.23) and sharp gradients in the stress field. Severe penetration of the elements in the slave surface by the punch master surface is also observed.
Figure 2.22: Effective stress and components in the punch and blank at a stroke of 0.05 mm. The stresses are given in MPa.
Figure 2.23: Effective stress and components in the punch and blank at a stroke of 0.15 mm. The
To resolve these problems the elements in the refined area was divided by 2 in the x direction. This resulted in a better interaction between the punch and blank and
smoothened the gradients in the stress field see Figure 2.24. The effective stress level was about 1000 MPa evenly distributed around the radius instead of 2600 MPa directly under the node contact points, which were the case in Figure 2.23. With these results in mind, conclusions can be made that a very high level of mesh refinement is needed to resolve the state of stress in the punch during the cutting operation when using this type of contact condition.
Figure 2.24: Effective stress and stress components in the punch and tool at a stroke of 0.05mm.
The stresses are given in MPa.
Figure 2.25: Effective stress and components in the punch and blank at a stroke of 0.15mm. The stress is given in MPa.
2.5 Comments
It has been shown in this section that simulations with ABAQUS using an axisymmetrical model can predict the maximum punch force and stroke to failure with good accuracy, using a plastic failure strain of 4. Calibrations against experiments have to be done to determine the value of the plastic failure strain.
Comparisons with simulations made with the DEFORM FE-code showed good agreement in terms of the characteristic edge zones and stroke to failure when using a value on the failure parameter that produces a stroke as measured in experiments.
Simulations using elastic tools showed that the elastic response had a small effect on the load curve (Figure 2.21), causing a delay in the rise of the punch force and a shorter stroke to failure. A greater difference is expected if simulations with tool dimensions similar to the ones used in experiments would be made, this has not been done because of the long cpu time such a simulation would require. The tool stress could not be resolved using the mesh density from the simulations with rigid tools because of the single sided contact condition allowed penetration by the punch radius into the mesh of the blank. This was improved when the mesh density in the blank was increased in the radial direction,
showing more reasonable stress levels, but further studies has to be made if estimations of the accuracy shall be possible.
3. Material properties
3.1 Materials
The materials that will be used in the simulations are DC04, Docol 350 YP, Docol 800 DP and Docol 1400M. The Docol steels are cold reduced high and ultra high strength qualities.
The corresponding numerical value to each Docol grade corresponds to the lowest yield strength for the YP qualities and lowest tensile strength for the DP and M qualities. DC04 is a mild steel quality suitable for deep drawing.
3.2 Tensile tests
To generate the necessary input data for the simulations tensile tests have been performed for the above mentioned materials. Specimens for each material have been created in the rolling direction (0 degrees), 45 degrees and transverse to the rolling direction (90
degrees). The tests were made and documented by SSAB Tunnplåt. True stress versus true strain curves generated during the tests can be found in Appendix I.
3.2.1 Mechanical properties
The mechanical properties for the different materials in the directions mentioned earlier are given in Table 3.1. The parameter r is the plastic strain ratio that is used to assess the anisotropy of blank materials [7]. It can be expressed in terms of strains as:
t
r w
ε
=ε (5)
Were εw is the strain in the width direction and ε is the strain in the thickness direction.t Another parameter that is useful for determining the properties of sheet metal is rwhich is measure of the normal anisotropy i.e. the average anisotropy in the plane of the sheet. A high r value shows that the preferred flow direction is in the plane of the sheet (r>1) while a low value (r <1) indicates a preferred flow direction in the thickness direction. The normal anisotropy parameter can be calculated from the plastic strain ratios with the
following formula [8]:
4 2 45 90
0 r r
r = r + + (6)
Were the indices denote the orientation of the tensile direction given as the angle to the rolling direction.
Table 3.1: Mechanical properties of the different sheet materials.ris the normal anisotropy.
Material designation
Direction Rp0.2 [MPa] Rm [MPa] r r n
Docol 350 YP 0° 363 423 0.67 0.11
45° 369 416 0.97 0.11
90° 403 432 1.05 0.92 0.11
Docol 800 DP 0° 653 873 0.70 0.11
45° 645 862 0.98 0.11
90° 659 887 0.87 0.88 0.11
Docol 1400 M 0° 1277 1486 0.75 0.11
45° 1258 1459 0.61 0.11
90° 1278 1496 0.62 0.65 0.11
DC04 0° 184 308 2.26 0.21
45° 187 315 1.55 0.20
90° 176 299 2.86 2.06 0.20
3.2.2 Extrapolation of test data
Data produced with uniaxial tensile tests is only valid until the point of necking. Necking causes a triaxial stress state in the specimen that alters the strain hardening behaviour. The strain developed in the blank during a cutting operation exceeds the necking strain in uniaxial tension by far, which calls for other procedures to characterise the material response under high strain. One way to overcome this problem is extrapolating the yield data from the tensile test, starting with the point at the onset of necking. This has been done using a method called double Voce extrapolation were the yield data is fitted to the expression given by:
) 1
( ) 1
( 2 3 4
1 0
P
P C
C C e
e
C ε ε
σ
σ = + − − + − − (7)
Were σ is the initial yield strength of the material and 0 ε is the plastic strain. A methodp for determining the constants in Equation (7) has been developed by Dr Mats Sigvant [13].
This method is based on using the Ag value and two points before to resolve the constants K and n in a power law extrapolation such as:
n
K(εp)
σ = (8)
Thereafter the stress at a fourth point is computed by adding a strain of four percent to the strain that is found at the Ag and inserting this value into Equation (8). This set of points is then used to compute the constants in Equation (7). The resulting curve increases almost like a power law after the Ag value but as the strain increases the exponential terms
vanishes and leaves only the initial yield stress and the constants C1 and C3 thus forming a plateau at higher strains. See Figure 3.1 for extrapolated yield curves concerning the different materials used in this study.
Double Voce extrapolated yield curves
0.00 200.00 400.00 600.00 800.00 1000.00 1200.00 1400.00 1600.00 1800.00
0.0 0.5 1.0 1.5 2.0 2.5 3.0 Plastic strain
Yield stress [MPa]
Docol 800 DP Docol 1400 M Docol 350 YP DC04
Figure 3.1: Double Voce extrapolated yield data for the different materials.
4. Experiments
To generate data for validation of the simulations, punching tests with the before mentioned materials have been performed by Uddeholm Tooling AB. The experiments give information about the characteristic zones of the cut edges, forces used to penetrate the blank and the stroke required to complete the punching operation. For a summarisation of experimental data see Appendix 2.
4.1 Experimental setup
The machine used to produce the experimental data is an excenter driven press with sensors for measuring position and forces added to the main tripod. A piezoelectric load cell is used for registering forces while an inductive position sensor is measuring the position with a sampling frequency of 10kHz. The layout of the machine and position of the sensors can be seen in Figure 4.1-2.
Figure 4.2: The main tripod with the position measuring sensor marked with a red arrow.
4.2 Data from experiments
In order to be able to compare the load curves from the experimental punching machine with simulations, some postprocessing had to be done. Because of the limited space in the machine, the load cell used to register the forces had to be placed in such a way that all forces used in the operation is registered see Figure 4.3 for an example of an experimental load curve.
Punch forces DC350YP
6000 8000 10000 12000 14000 16000 18000 20000 22000
-0.2 0.8 1.8 2.8 3.8
Distance betw een edges [mm]
Total force [N]
Figure 4.3: Experimental load curve
To extract the data of interest from the experimental load curve a straight line was fitted to the part of the force curve were the punch had not yet made contact with the blank and only the force of the blank holder is being registered, see Figure 4.4. The equation of this line gives the distance dependency of the holder force, which then can be withdrawn from the total force. The stroke is also corrected so the zero level is when the punch makes first contact with the blank.
Blank holder force vs stroke
0 2000 4000 6000 8000 10000
0 0.5 1 1.5 2 2.5 3 3.5 4
Distance between edges [mm]
Force [N]
Figure 4.4: Straight line fitted to the blank holder force
The above corrections have been made to all experimental load curves and the result is presented in Figure 4.5
Experimental load curves
-5000 0 5000 10000 15000 20000 25000 30000 35000
0 0.2 0.4 0.6 0.8 1
Stroke [mm]
Force [N]
DC04 Experiment DC350YP Experiment DC800DP Experiment DC1400M Experiment
Figure 4.5: Corrected experimental load curves.
5. Simulation model
The experimental data presented in Section 4.2 have been produced with different tool geometry than were used during the parameter study in Section 2. A different material model was also to be used for DC04. Therefore a new simulation model had to be made to simulate the punching operation.
5.1 Geometry and mesh density
An axisymmetric FE-model was used in all simulations with tool and blank dimensions matching those that were used in the experiments, see Figure 5.1 and Table 5.1. The mesh density was 128 elements through thickness in the sheared zone and were chosen based on
Figure 5.1: The axisymmetric FE-model.
Table 5.1: Tool and blank dimensions used with the different materials.
Material r1 [mm] r2 [mm] r3 [mm] r4 [mm] c1 [mm] c2 [mm] t [mm]
DC04 4.94 0.02 0.01 0.02 0.06 0.26 0.98
DC350YP 4.94 0.02 0.01 0.02 0.06 0.26 1
DC800DP 4.90 0.02 0.01 0.02 0.10 0.3 0.96
DC1400M 4.90 0.02 0.01 0.02 0.10 0.3 0.97
5.2 Material model
The isotropic von Mises material model described in Section 2.1.3 have been used for all materials except DC04 for which the model proposed by Hill in 1948 have been evaluated.
Usage of this model is motivated by the high values of the plastic strain ratio shown by DC04. The yield criterion for this material model can be expressed in terms of stress components as:
2 12 2
31 2
23 2
22 11 2
11 33 2
33
22 ) ( ) ( ) 2 2 2
( )
(σ F σ σ G σ σ H σ σ Lσ Mσ Nσ
f = − + − + − + + + (9)
Were the constants F, G, H, L, M and N are obtained by tests of the material in different directions. The definition of these constants is:
2 , 3 2 ,
3 2 ,
3
1 ), 1 ( 1
2 1
1 ), 1 ( 1
2 1
1 ), 1 ( 1
2 1
2 12 2 13 2 23
2 33 2 22 2 11
2 22 2 11 2 33
2 11 2 33 2 22
N R M R L R
R R H R
R R G R
R R F R
=
=
=
− +
=
− +
=
− +
=
(10)
The parameters R11, R22, R33, R12, R13 and R23 are anisotropic yield stress ratios defined as follows:
0 23 0 13 0 12 0 33 0 22 0
11 , , , , ,
τ σ τ σ τ σ σ σ σ σ σ
σ (11)
The type of anisotropy that will be modelled here is the normal anisotropy, which describes the material as isotropic in the plane with increased yield strength in the
thickness direction [11]. This kind of behaviour is modelled by setting (with the directions 1 and 3 in the plane of the blank and 2 in the thickness direction):
R11=R33=1 (12)
2 1
22
= r+
R (13)
Were ris the normal anisotropy parameter given by Equation (6). The remaining parameters are set to a value of unity.
5.2.1 Strain rate dependency
The material is modelled with a strain rate dependent function to account for the increased yield strength during high strain rate deformation. This is implemented in ABAQUS by specifying a yield stress ratio as a tabular function of the equivalent plastic strain rate:
0
) ( σ
ε σ pl
R
= (14)
Were )σ is the yield stress at a specified plastic strain rate and (εpl σ is the reference0 yield stress at a zero strain rate. The yield stress ratios used for the different materials and plastic strain rates can be found in Table 5.2.
Table 5.2: Yield stress ratios for the materials used in the simulations at different levels of equivalent plastic strain rate
Yield stress ratio Equivalent plastic strain rate
DC04 DC350YP DC800DP DC1400M
1 1 1 1 0
1.064 1.069 1.04 1.054 10
1.145 1.126 1.072 1.078 100
1.28 1.191 1.108 1.106 1000
1.28 1.191 1.145 1.127 10000
1.28 1.191 1.193 1.160 100000
5.3 Failure criterion
The shear failure model that is described in Section 2.1.5 has been used in all simulations to predict the rupture behaviour.
5.4 Contact and boundary conditions
5.4.1 Contact conditions
At first attempts were made using the pure master-slave kinematic contact definition as in Section 2. This proved to be insufficient when simulating materials that showed large plastic deformation before rupture. When the mesh became distorted the distance between the nodes in the blank became too great which caused penetration at the punch edge as can bee seen in Figure 5.2.
Figure 5.2: Penetration of the punch edge into the blank. The part of the punch that has penetrated the blank is marked with red.
This effect altered the state of deformation around the punch edge, allowing high strains to be generated in the elements thus triggering the failure criterion at an early stage of the simulation. To prevent this phenomenon a new approach was made by defining a balanced penalty contact condition between the edge and side of the punch and the top elements of the blank as shown in Figure 5.3. A similar definition was also added between the blank and die.
Figure 5.3: The penalty contact pair consisting of a) the discrete rigid punch and b) the top element surface of the blank.
The penalty contact is a less stringent enforcement of contact than the kinematic contact
slave surfaces. With the penalty contact condition the algorithm searches for slave node penetrations and if penetration is found a force is applied to the slave nodes to oppose the penetration while an force that is equal acts on the master surface on the penetration point.
This is often referred to in the literature as adding a spring between the surfaces to prevent penetration [14]. When the balanced penalty contact algorithm is used the forces are computed in two steps were the master and slave designation is switched between the surfaces. An average of the two forces is then applied [11]. When using this approach of contact modelling the penetration of the punch edge into the blank is eliminated, see
Figure 5.4: Deformation state around the punch edge with balanced penalty contact.
When the first element is removed from the top surface of the blank the penalty contact is no longer active. The contact between the blank and punch is then enforced by a kinematic contact condition. The motivation for two different contact algorithms to handle the a single contact between two surfaces is that ABAQUS does not allow using the same contact algorithm to enforce contact with both the top elements of the blank and the internal nodal surface that becomes exposed during the rupture.
5.4.2 Boundary conditions
As the model is axisymmetric the only boundary condition that needs to be applied directly to the blank is a locking of the nodes at the symmetry line in the radial direction. To hold the blank in place during the punching operation a force of 10kN is applied in the
downward direction on the blank holder. This force is applied during the first step of the solution with a ramp function during 0.001s and is then held at this level throughout the simulation.
Friction between the blank and tools are of Coulomb type, which is described by Equation (2). A friction coefficient of 0.2 has been used throughout the simulations.
The punch is controlled by a prescribed velocity boundary condition were it is given a velocity of 70 mm/s. This velocity has been determined by a linear regression on experimental data, see Figure 5.5. The velocity is then found by differentiating the equation of the regression line with respect to time.
Punch stroke vs time y = 70.144x + 0.0261
0 0.2 0.4 0.6 0.8 1 1.2
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 Time [s]
Stroke [mm]
Figure 5.5: A straight line fit to experimental data concerning the stroke as a function of time. The zero level for both time and stroke is set at the position were the punch makes first contact with the blank.
