FORMING OF ALUMINUM ALLOYS AT ELEVATED TEMPERATURES By Nader Elias Abedrabbo
A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mechanical Engineering 2005
ABSTRACT FORMING OF ALUMINUM ALLOYS AT ELEVATED TEMPERATURES By Nader Elias Abedrabbo Weight reduction has long been identified as a key priority for improving automotive fuel economy, and many studies often suggest substituting lightweight materials for typical steel applications. However, replacing steel in the structure and body of an automobile with lighter materials such as aluminum, magnesium, plastics, composites, etc., can be costly and is not simple or straightforward. Aluminum sheet, in particular, has lower formability at room temperature than typical sheet steel. Numerical analysis is critically important to understanding the complex deformation mechanics that occur during sheet forming processes. Finite element analysis (FEA) and simulations are used in automotive design and formability processes to predict deformation behavior accurately during stamping operations. Although available commercial FEA codes offer a library of material models applicable to a variety of applications, they often do not offer highly specialized material models developed for a specific material and process. Also, very few available material models are capable of handling complex forming simulations that incorporate the temperature-dependence of materials, especially for anisotropic materials such as aluminum sheets.
The goal of this research was to develop a temperature-dependent anisotropic material model for use in a coupled thermo-mechanical finite element analysis of the forming of aluminum sheets. The anisotropic properties of the aluminum alloy sheet AA3003-H111 were characterized for a range of temperatures 25°C–260°C (77°F–500°F) and for different strain rates. After performing standard ASTM uniaxial tensile tests, the anisotropy coefficients for two accurate material models: Barlat YLD96 and Barlat YLD2000-2d, both in the plane stress condition, were calculated for several elevated temperatures. The developed temperature-dependent anisotropic material models (YLD96 and YLD2000-2d) were then implemented in the commercial FEM code LS-Dyna as a user material subroutine (UMAT) using the cutting-plane algorithm proposed by Simo et al. (1985) for the integration of a general class of elastoplastic constitutive models. Exact implementation of the yield functions and their derivatives are presented. The temperature-dependent material models were then used to numerically simulate the thermo-coupled finite element model for pure stretch conditions in order to compare the accuracy of the UMAT’s ability to predict both forming behavior and failure locations with experimental results of the pure stretch forming process for AA3003-H111 under several elevated temperatures. The favorable comparison found between the numerical and experimental data shows that a promising future exists for the development of more accurate temperature-dependent yield functions to apply to thermo-hydroforming process.
Copyright By Nader Elias Abedrabbo 2005
To My Family
v
ACKNOWLEDGMENTS
I would like to convey my sincere thanks to Dr. Farhang Pourboghrat for his advice and his continuous support during the course of my studies at Michigan State University. His leadership and encouragement were instrumental in achieving this goal. I would like also to thank my Ph.D. committee members: Dr. Ron Averill, Dr. Frédéric Barlat from ALCOA, Dr. Charles MacCluer and Dr. Thomas Pence, their constructive advice and suggestions helped me throughout the work. Special thanks are due to Dr. John Carsley from General Motors. Without his advice and help in running experiments at GM it would not have been possible to complete this work. Sincere thanks also for Drs. Paul Krajewski and Anil Sachdev from the GM Research and Development Center for their assistance and helpful discussions in support of this research. This research was made possible through funding from General Motors. Their continuous support for our research activities during the years is greatly appreciated. Most of all I would like to thank my family for their support during these years. Without their support and encouragement during the years I would not have achieved this.
vi
TABLE OF CONTENTS
LIST OF TABLES .................................................................................................ix LIST OF FIGURES ...............................................................................................xi LIST OF ABBREVIATIONS ................................................................................xxi CHAPTER 1 INTRODUCTION .................................................................................................. 1 1.1 Research Objectives ................................................................................... 4 1.2 Breakdown of Research Objectives ............................................................ 6 CHAPTER 2 LITERATURE REVIEW ........................................................................................ 8 2.1 Warm Forming Literature Review ............................................................... 8 2.2 Sheet Metal Hydroforming Literature Review............................................ 11 2.3 Wrinkling Literature Review ...................................................................... 16 CHAPTER 3 THEORETICAL ANALYSIS ................................................................................ 20 3.1. Anisotropic Constitutive Models ............................................................... 20 3.1.1 Barlat YLD96 ...................................................................................... 20 3.1.2 Barlat YLD2000-2d ............................................................................. 23 3.2. Plastic Anisotropy Parameters ................................................................. 27 3.3. Anisotropy Coefficients Calculation.......................................................... 28 3.3.1 Barlat YLD96 ...................................................................................... 28 3.3.2 Barlat YLD200-2d ............................................................................... 34 3.4. Constitutive Equations (Flow Stress) ....................................................... 37 CHAPTER 4 MATERIAL CHARACTERIZATION .................................................................... 40 4.1. Materials .................................................................................................. 40 4.2. Experimental Procedure........................................................................... 41 4.3. Hardening Model...................................................................................... 45 4.4. Barlat’s Yield 96 Anisotropy Coefficients Calculation ............................... 54 4.5. Barlat’s YLD2000-2d Anisotropy Coefficients Calculation........................ 71
vii
CHAPTER 5 STAMPING EXPERIMENTS .............................................................................. 87 5.1. Experimental Apparatus........................................................................... 87 5.2. Experimental Work................................................................................... 92 CHAPTER 6 NUMERICAL ANALYSIS setup .......................................................................... 99 6.1. Preliminary Results ................................................................................ 100 6.2. Coupled Thermal Structural Finite Element Model................................. 105 6.3. Failure Criteria ....................................................................................... 109 CHAPTER 7 STRESS INTEGRATION FOR ELASTO-PLASTICITY USING ANISOTROPIC YIELD FUNCTIONS ......................................................................................... 113 7.2. Explicit derivation of YLD96 and its derivative ....................................... 128 7.3. Explicit derivation of YLD2000-2d and its derivative .............................. 133 CHAPTER 8 DISCUSSION OF RESULTS: NUMERICAL Vs. EXPERIMENTAL ................. 136 CHAPTER 9 CONCLUSIONS ............................................................................................... 150 BIBLIOGRAPHY............................................................................................... 154
viii
LIST OF TABLES
Table 3.1
Experimental data needed to calculate yield function coefficients for YLD96 in the plane stress state. .....................
28
Experimental data needed to calculate yield function coefficients for YLD2000-2d. ...................................................
35
Table 4.1
Chemical composition of AA3003-H111 (wt. %)......................
40
Table 4.2
As-received mechanical properties of AA3003-H111..............
40
Table 4.3
Material properties of AA3003-H111 at elevated temperatures. The values shown below represent an average of three tensile tests. (The hardening values are for the rolling direction).................................................................
49
Summary of equations used to fit hardening parameters for the power law flow rule. Temperatures in °C...........................
53
Barlat’s YLD96 model anisotropy coefficients calculated at several elevated temperatures. ...............................................
54
Barlat YLD96 material model anisotropy coefficients for AA3003-H111 as a function of temperature using two different order polynomial fit functions. Temperatures in °C....
65
Barlat’s YLD2000-2d model anisotropy L′ coefficients calculated at several elevated temperatures. ..........................
71
Barlat’s YLD2000-2d model anisotropy L′′ coefficients calculated at several elevated temperatures. ..........................
72
Table 3.2
Table 4.4
Table 4.5
Table 4.6
Table 4.7
Table 4.8
ix
Table 4.9
Barlat YLD2000-2d material model anisotropy coefficients for AA3003-H111 as a function of temperature using a 4th order polynomial fit function. Temperatures in °C. ..................
80
Table 4.10 Barlat YLD2000-2d material model anisotropy coefficients for AA3003-H111 as a function of temperature using a 4th order polynomial fit function. Temperatures in °C. .................
81
Table 6.1
Thermal properties of material used in numerical analysis. .... 107
Table 7.1
Stress update algorithem based on incremental theory of plasticity. ................................................................................. 123
Table 8.1
Measured temperatures of dies, punch and blank during experiments............................................................................. 137
Table 8.2
Numerical prediction of failure depth compared to exprimental results. ................................................................. 148
x
LIST OF FIGURES
Figure 3.1
Determination of the anisotropy parameter 2β from Mohr’s circle.....................................................................................
22
Characterization of required tests to be used in calculation of the anisotropic values in relation to the yield locus. .........
29
Figure 3.3
Flow stresses at equal amounts of plastic work (Wb=Wu). ..
30
Figure 4.1
Rectangular tension test specimen, ASTM-E8. All dimensions inmm. ................................................................
41
Figure 4.2
Tensile specimen location with respect to sheet. .................
42
Figure 4.3
Stress values at equal values of plastic work (20 MPa/unit-vol) as a function of temperature. Data shown for tests at 0°, 45° and 90° directions at multiple elevated temperatures and for balanced biaxial tension (bulge) data at room temperature. ...........................
44
True stress strain data of AA3003-H111 at room temperature for uniaxial tests at 0°, 45° and 90° directions and for balanced biaxial tension (bulge) data.......................
45
Engineering stress strain curves of AA3003-H111 at several elevated temperatures for the rolling direction. ........
46
True stress strain data of AA3003-H111 from uniaxial tests at 204°C (400°F) at several strain rates for the rolling direction. ....................................................................
47
Strength hardening coefficient (K) of AA3003-H111 along the rolling direction as a function of temperature at 0.5 mm/mm strain rate. Solid line is a linear curve fit. ...............
48
Figure 3.2
Figure 4.4
Figure 4.5
Figure 4.6
Figure 4.7
xi
Figure 4.8
Figure 4.9
Figure 4.10
Figure 4.11
Figure 4.12
Figure 4.13
Figure 4.14
Figure 4.15
Figure 4.16
Strength hardening coefficient (K) of AA3003-H111 along the rolling direction as a function of temperature at 0.5 mm/mm strain rate. Solid line is a linear curve fit. ...............
50
Strain-hardening exponent (n) of AA3003-H111 along the rolling direction as a function of temperature at 0.5 /min strain rate. Solid line is a linear curve fit. .............................
51
Strain-rate sensitivity index (m) of AA3003-H111 along the rolling direction as a function of temperature. Solid line is an exponential curve fit.....................................................
52
Plastic anisotropy parameters (Rθ) of AA3003-H111 as a function of temperature, calculated at 0.5/min strain rate (ASTM-E517). ......................................................................
53
2D-Plot of von Mises yield function vs. Barlat’s YLD96 yield function at 25°C (77°F) for AA3003-H111....................
55
2D-Plot of Barlat’s YLD96 yield function for AA3003-H111 at several elevated temperatures. Normalized stress values...................................................................................
56
2D-Plot of Barlat’s YLD96 yield function for AA3003-H111 at several elevated temperatures. Plot of actual stresses to show the change in the size of the yield surface. .............
57
Effect of αx, αy and αz0 on the yield function (Barlat et al., 1997). ...................................................................................
58
Plot of experimental values of plastic anisotropy parameter (R0) of AA3003-H111 compared with predictions from the Barlat YLD96 material model. Calculations were done first using variable values of the yield function with each temperature and then by using a fixed value for all temperatures. ...........................................
60
xii
Figure 4.17
Figure 4.18
Figure 4.19
Figure 4.20
Figure 4.21
Figure 4.22
Figure 4.23
Figure 4.24
Figure 4.25
Plot of (c1) as a function of temperature for AA3003-H111. 3rd and 5th order polynomial curve fits are shown also. ........
61
Plot of (c2) as a function of temperature for AA3003-H111. 3rd and 5th order polynomial curve fits are shown also. ........
62
Plot of (c3) as a function of temperature for AA3003-H111. 3rd and 5th order polynomial curve fits are shown also. ........
62
Plot of (c6) as a function of temperature for AA3003-H111. 3rd and 5th order polynomial curve fits are shown also. ........
63
Plot of (αx) as a function of temperature for AA3003-H111. 3rd and 5th order polynomial curve fits are shown also. ........
63
Plot of (αy) as a function of temperature for AA3003-H111. 3rd and 5th order polynomial curve fits are shown also. ........
64
Plot of (αz1) as a function of temperature for AA3003H111. 3rd and 5th order polynomial curve fits are shown also.......................................................................................
64
Plot of experimental values of plastic anisotropy parameter (R0) for AA3003-H111 compared with predictions from Barlat’s YLD96 material model where the yield function parameters are a function of temperature. Calculations were done using 3rd and 5th order curve fits for the yield function anisotropy coefficients. ........................
66
Plot of experimental values of plastic anisotropy parameter (R45) for AA3003-H111 compared with predictions from Barlat’s YLD96 material model where the yield function parameters are a function of temperature. Calculations were done using 3rd and 5th order curve fits for the yield function anisotropy coefficients. ........................
67
xiii
Figure 4.26
Figure 4.27
Figure 4.28
Figure 4.29
Figure 4.30
Figure 4.31
Figure 4.32
Figure 4.33
Plot of experimental values of plastic anisotropy parameter (R90) for AA3003-H111 compared with predictions of Barlat’s YLD96 material model where the yield function parameters are a function of temperature. Calculations were done using 3rd and 5th order curve fits for the yield function anisotropy coefficients. ........................
68
Plot of Barlat’s YLD96 yield function for AA3003-H111 using anisotropy coefficients extracted from experimental data compared to anisotropy coefficients from both the 3rd and 5th order curve fits shown in Table 4.6. Temperature was 25°C (77°F). ..................................................................
69
Plot of Barlat’s YLD96 yield function for AA3003-H111 using anisotropy coefficients extracted from experimental data compared to anisotropy coefficients from both the 3rd and 5th order curve fits shown in Table 4.6. Temperature was 204°C (400°F). ..............................................................
70
2D-Plot of von Mises yield function vs. Barlat’s YLD20002d yield function at 25°C (77°F) for AA3003-H111...............
73
2D-Plot of Barlat’s YLD000-2d yield function for AA3003H111 at several elevated temperatures. Normalized stress values. .......................................................................
74
Plot of ( L′ 11) as a function of temperature for AA3003H111 Experimental data is fitted using a 4th order polynomial curve. .................................................................
75
Plot of ( L′ 12) as a function of temperature for AA3003H111 Experimental data is fitted using a 4th order polynomial curve. .................................................................
76
Plot of ( L′ 21) as a function of temperature for AA3003H111 Experimental data is fitted using a 4th order polynomial curve. .................................................................
76
xiv
Figure 4.34
Figure 4.35
Figure 4.36
Figure 4.37
Figure 4.38
Figure 4.39
Figure 4.40 Figure 4.41
Plot of ( L′ 22) as a function of temperature for AA3003H111 Experimental data is fitted using a 4th order polynomial curve. .................................................................
77
Plot of ( L′ 66) as a function of temperature for AA3003H111 Experimental data is fitted using a 4th order polynomial curve. .................................................................
77
Plot of ( L′′ 11) as a function of temperature for AA3003H111 Experimental data is fitted using a 4th order polynomial curve. .................................................................
78
Plot of ( L′′ 12) as a function of temperature for AA3003H111 Experimental data is fitted using a 4th order polynomial curve. .................................................................
78
Plot of ( L′′ 21) as a function of temperature for AA3003H111 Experimental data is fitted using a 4th order polynomial curve. .................................................................
79
Plot of ( L′′ 22) as a function of temperature for AA3003H111 Experimental data is fitted using a 4th order polynomial curve. .................................................................
79
Plot of ( L′′ 66) as a function of temperature for AA3003H111 Experimental data is fitted using a 4th order polynomial curve. ................................................................. Plot of experimental values of plastic anisotropy parameter (R0) for AA3003-H111 compared with predictions from Barlat’s YLD2000-2d material model where the yield function parameters are a function of temperature. Calculations were done using the 4th order curve fit for the yield function anisotropy coefficients from Tables 4.9 and 4.10. ............................................................
xv
80
82
Figure 4.42
Figure 4.43
Figure 4.44
Figure 4.45
Figure 5.1
Figure 5.2
Figure 5.3
Figure 5.4
Plot of experimental values of plastic anisotropy parameter (R45) for AA3003-H111 compared with predictions from Barlat’s YLD2000-2d material model where the yield function parameters are a function of temperature. Calculations were done using the 4th order curve fit for the yield function anisotropy coefficients from Tables 4.9 and 4.10. ............................................................
83
Plot of experimental values of plastic anisotropy parameter (R90) for AA3003-H111 compared with predictions from Barlat’s YLD2000-2d material model where the yield function parameters are a function of temperature. Calculations were done using the 4th order curve fit for the yield function anisotropy coefficients from Tables 4.9 and 4.10. ............................................................
84
Plot of Barlat’s YLD2000-2d yield function for AA3003H111 using anisotropy coefficients extracted from experimental data compared to anisotropy coefficients from both the 4th order curve fit shown in Tables 4.9 and 4.10. Temperature was 25°C (77°F). ...................................
85
Plot of Barlat’s YLD2000-2d yield function for AA3003H111 using anisotropy coefficients extracted from experimental data compared to anisotropy coefficients from both the 4th order curve fit shown in Tables 4.9 and 4.10. Temperature was 260°C (500°F). ...............................
86
The modified Interlaken servo press 75 used for the forming of aluminum sheets. ................................................
88
Schematics of a warm-forming press using a hemispherical punch with a blank holding support. ..............
90
Experimental stamping setup showing heating elements and controller. ......................................................................
91
Pure stretch experimental results using the 101.6mm (4in) hemispherical punch at 25°C (77°F). ...................................
93
xvi
Figure 5.5
Figure 5.6
Figure 5.7
Figure 5.8
Figure 5.9
Figure 5.10
Pure stretch experimental results using the 101.6mm (4in) hemispherical punch at 93°C (200°F). .................................
93
Pure stretch experimental results using the 101.6mm (4in) hemispherical punch at 149°C (300°F). ...............................
94
Pure stretch experimental results using the 101.6mm (4in) hemispherical punch at 177°C (350°F). ...............................
94
Pure stretch experimental results using the 101.6mm (4in) hemispherical punch at 204°C (400°F). ...............................
95
Pure stretch experimental results using the 101.6mm (4in) hemispherical punch at several elevated temperatures as indicated along with failure punch depth. The failure locations are also shown. .....................................................
96
Punch force vs. punch depth for the hemispherical punch experiments at several elevated temperatures.....................
97
Figure 6.1
Initial results of punch force vs. punch depth for the hemispherical punch compared to numerical predictions at 204°C (400°F) using only structural FEM analysis. .......... 101
Figure 6.2
A schematic showing the effect of temperature drop away from the band heaters. The dies which are contacting the band heaters are at a higher temperature than the punch (which is not heated directly), and the center of the sheet. .. 102
Figure 6.3
A blank divided into multiple sections with each section having different temperature. Temperatures used in the analysis are shown on the graph.......................................... 103
Figure 6.4
Punch force vs. punch depth for the hemispherical punch compared to numerical predictions for the sectioned blank in Figure 6.3. ........................................................................ 104
xvii
Figure 6.5
LS-Dyna Full 3-D model created for stamp warm forming process with a hemispherical punch, using square blank..... 106
Figure 6.6
Punch velocity used for the finite element simulation. .......... 108
Figure 6.7
Forming limit diagrams for AA3003-H111 at 25°C (77°F). Calculations based on the M-K model and Barlat’s YLD96 anisotropic material model using two types of hardening law: Voce and Power law. .................................................... 111
Figure 6.8
Forming limit diagrams (FLD’s) for AA3003-H111 based on the M-K model, Barlat’s YLD96 anisotropic yield function, and Voce hardening law at several elevated temperatures. ....................................................................... 112
Figure 7.1
Computed values of elastic, plastic and thermal strain increments for an element at different times. The results shown correspond to thermoforming analysis at the elevated temperature of 204°C (400°F). .............................. 114
Figure 7.2
Geometrical interpretation of the cutting plane algorithm. ) The trial stress state σ n(trial +1 is returned iteratively to the yield surface. ........................................................................ 124
Figure 7.3
Accuracy of the developed UMAT with respect to step size compared to true stress strain results of a uniaxial test performed at 25 °C. ....................................................... 126
Figure 7.4
Determination of the anisotropy parameter 2β from Mohr’s circle..................................................................................... 130
Figure 8.1
Finite element results from a fully coupled thermomechanical simulation of thermo-forming at 25°C (77°F). Failure location is shown. The sheet failed at a punch depth of 28 mm (1.1 in). ....................................................... 138
xviii
Figure 8.2
Distribution of the major and minor strains on the FLD. The sheet failed where the strains crossed the FLD curve... 139
Figure 8.3
Results of numerical simulation of thermo-coupled FEM analysis at 93°C (200°F). Failure location shown on the graph. Punch depth = 32mm (1.26in)................................... 140
Figure 8.4
Results of numerical simulation of thermo-coupled FEM analysis at 93°C (200°F). The graph shows the distribution of the major and minor strains on the FLD. The sheet failed where the strains crossed the FLD curve... 141
Figure 8.5
Results of numerical simulation of thermo-coupled FEM analysis at 149°C (300°F). Failure location shown on the graph. Punch depth = 35mm (1.37in)................................... 142
Figure 8.6
Results of numerical simulation of thermo-coupled FEM analysis at 149°C (300°F). The graph shows the distribution of the major and minor strains on the FLD. The sheet failed where the strains crossed the FLD curve... 143
Figure 8.7
Results of numerical simulation of thermo-coupled FEM analysis at 204°C (400°F). Failure location shown on the graph. Punch depth = 38mm (1.5in)..................................... 144
Figure 8.8
Results of numerical simulation of thermo-coupled FEM analysis at 204°C (400°F). The graph shows the distribution of the major and minor strains on the FLD. The sheet failed where the strains crossed the FLD curve... 145
Figure 8.9
Numerical simulation results of pure stretch using the 101.6mm (4in) hemispherical punch at several elevated temperatures. Punch depth at failure and failure locations based on FLD prediction are also shown. ............................ 146
xix
Figure 8.10
Experimental results of pure stretch using the 101.6mm (4in) hemispherical punch at several elevated temperatures. Punch depth at failure and failure locations are shown............................................................................. 147
Figure 8.11
Punch force vs. punch depth for the hemispherical punch stretch experiments at several elevated temperatures. Numerical results closely match the experimental results. ... 148
xx
LIST OF ABBREVIATIONS
σ - Flow stress K- Strength hardening coefficient n- Strain-hardening exponent m- Strain-rate sensitivity index
ε p - Effective plastic strain
ε - Strain rate. ε0 - A constant
ε sr 0 - A base strain rate (constant) T- Temperature in °C
εn +1 - Total strain increment
e - Elastic strain increment εn+ 1
C - Fourth order elasticity tensor
σ - Cauchy stress tensor
ε p - Plastic strain tensor
Δε p - Equivalent plastic strain increment (trial ) - Trial stress states σ n+ 1
σ y = Yield Stress E- Young’s modules Φ - Yield Function
xxi
S1, S2 & S3- Principal Values of the Stress Deviator S – Fourth Order Tensor L – Fourth Order Linear Operator
c1, c2, c3, c6, αx, αy, αz0 and αz1: Coefficients that describe the anisotropy of the material dεw - Width Strain Increment d ε t - Thickness Strain Increment
εxx - Strain Rate in xx Direction Rθ - Plastic Anisotropy. θ = 0, 45 and 90
xxii
CHAPTER 1
INTRODUCTION
Weight reduction has long been identified as a key priority for improving automotive fuel economy (Greene et al., 2000), and many studies often suggest substituting lightweight materials for typical steel applications. However, replacing steel in the structure and body of an automobile with lighter materials such as aluminum, magnesium, plastics, composites, etc., can be costly and is not simple or straightforward. Aluminum sheet, in particular, has lower formability at room temperature than typical sheet steel (Ayres et al., 1979). High temperature forming methods based on superplastic behavior of Al-Mg alloys have been used to produce automotive closure panels that far exceed the conventional stamping formability of steel (Schroth, 2004). Superplastic formability however, requires fine-grained microstructures and slower strain rates which can affect production cost. The formability of typical ‘off-the-shelf’ automotive aluminum sheet alloys (5182-O and 5754-O) can be greatly improved by warm forming (Li et al., 2003) without the increased costs of refining the microstructure. The elevated temperature corresponds with decreased flow stress and increased ductility in the sheet, which can lead to deeper drawing and more stretching to form panels without design modifications to the stamped steel product. Furthermore, the serrated flow behavior of Al-Mg sheet alloys (dynamic strain aging / PLC effect) (Robinson et al., 1994) and corresponding Lüder’s line
1
surface defects that result from the interactions of solutes with mobile dislocations can be avoided by deforming the sheet metal above a critical temperature. Numerical analysis is critically important to understanding the complex deformation mechanics that occur during sheet forming processes. Finite element analysis (FEA) and simulations are used in automotive design and formability processes to predict deformation behavior accurately during stamping operations. Confidence in the numerical analysis of formability depends on the accuracy of the constitutive model describing the behavior of the material (Chung et al. 1992). This is especially important when the material exhibits anisotropic characteristics, as do most cold rolled sheet metals. Previous research (Zampaloni et al., 2003; Abedrabbo et al., 2005c) demonstrated the importance of using appropriate material models with respect to wrinkling and ironing during the sheet hydroforming process of an Al-Mg-Si alloy (AA6111-T4). The material model used was the anisotropic yield function proposed by Barlat et al. (1997a) (YLD96), which is based on a phenomenological description of the material. Although available commercial FEA codes offer a library of material models applicable to a variety of applications, they often do not offer highly specialized material models developed for a specific material and process. Also, very few available material models are capable of handling complex forming simulations that incorporate the temperature-dependence of materials. The process becomes increasingly complicated when materials exhibit anisotropic behavior. Currently available commercial codes do not offer material models that are appropriate for
2
simulating the thermo-mechanical forming processes of anisotropic materials such as aluminum sheet alloys. The use of anisotropic material models in FEA is further complicated by the material dependence on the anisotropy coefficients that require thorough material characterization under multiple loading conditions (Chung et al., 1996). This difficulty is further exacerbated when temperature effects are introduced. Because material hardening behavior and material response to loading conditions change at elevated temperatures, the anisotropy coefficients must be determined as a function of temperature so that the material model can account for these changes. In this research, material characterization of the anisotropy coefficients of two accurate material models: YLD96 (Barlat et al., 1997) model and YLD2000-2d (Barlat et al., 2003), and the modified power-law flow stress were evaluated as a function of temperature. Also, the effect of temperature on the evolution of the yield surface and on the plastic anisotropy parameters R0, R45 and R90 is studied. Several fitting functions were compared for both the anisotropy coefficients of the yield function and for the flow stress. The temperature dependent anisotropic materials models are then implemented as a user material subroutine (UMAT) into the explicit finite element code LS-Dyna. The UMAT was then used in a coupled thermo-mechanical finite element analysis of stamping of the aluminum sheet with a hemispherical punch, under pure stretch boundary condition, and deformation behavior and failure locations were compared against experimental results. Next, a detailed description of the research objectives is presented.
3
1.1 Research Objectives
1. Material Characterization of AA3003-H111
To extract the anisotropy coefficients required for numerical analysis, a number of uniaxial and bulge experimental tests are carried out to characterize the mechanical behavior of the aluminum sheet alloy 3003-H111 at elevated temperatures. The parameters of interest to be measured using the standard ASTM-8 tests are: 1. Material tensile strength, σ α (α = 0 , 45 o ,90 o ) , along the rolling (0°), 45°, and transverse (90°) directions, respectively.
)
(
2. Anisotropy parameters, Rα α = 0,45o ,90o , along the rolling (0°), 45°, and transverse (90°) directions, respectively. 3. Strain-rate sensitivity (m) at elevated temperatures.
2. Experimental Heating Capability
The existing experimental setup will be fitted with a new heating capability that would enable performing stamp forming at elevated temperatures. A control process is to be done through a closed loop control system with thermocouples.
3. Warm-Forming Experiments
Once the heating capability is added; warm-forming of pure stretch experiments using a hemispherical punch with 101.6mm (4in) diameter will be 4
performed. These experiments will be conducted at temperatures ranging from 25°C–260°C (77°F–500°F).
4. User Material Subroutine (UMAT) of the Constitutive Model
To develop a temperature-dependant constitutive model, based on two accurate anisotropic yield functions: Barlat’s YLD96 and Barlat’s YLD2000-2d, to represent the behavior of the aluminum alloy sheet at elevated temperatures and to incorporate the developed constitutive models into the commercial finite element code LS-Dyna, using a User-defined Material subroutine (UMAT).
5. Numerical Analysis
An important goal in manufacturing research is determining the optimum method of production of efficient products with less cost. The optimization criterion varies, depending on the products, but having a thorough understanding of the manufacturing processes is an essential step. Sheet metal forming design requires the understanding of the fundamentals of deformation mechanics involved in the processes. Without proper understanding of the effect of different variables such as material properties, friction and geometry the design process would be difficult, time consuming, and expensive. Also it would not be possible to predict and prevent defects form occurring until it is too late. Numerical analysis of the warm-forming process will be carried out to verify the experimental results.
5
1.2 Breakdown of Research Objectives
1. Experiments:
o Characterization of the mechanical properties of an automotive aluminum sheet alloys (AA3003-H111) at elevated temperatures 25°C–260°C (77°F–500°F), using uniaxial tensile test. o Adding heating capability to current experimental setup. o Warm forming of aluminum alloy sheets with a hemispherical punch at elevated temperatures (e.g., up to a maximum of 400 °F).
2. Constitutive Modeling:
o Develop constitutive models, based on two accurate anisotropic yield functions, to characterize the behavior of the aluminum alloy sheet at elevated temperatures. o Incorporate the constitutive model into LS-Dyna finite element code, using a User-defined Material subroutine (UMAT).
3. Numerical Simulation:
o Simulate the pure stretch warm-forming of aluminum alloy sheets with a hemispherical punch at various elevated temperatures, and compare the results with experiments performed at MSU. o Simulate thermo-hydroforming of a complex automotive aluminum sheet metal part using the developed model. 6
Literature review of the forming process is presented in Chapter 2. Theoretical analysis is presented in Chapter 3.
In Chapter 4 material
characterization is presented with coefficients of the two yield functions, i.e. YLD96 and YLD2000-2d. Experimental stamping experiments and results using the warm forming process are presented in Chapter 5. The setup of the numerical analysis and the stress integration algorithm used to setup the FEM UMAT is presented in Chapter 6 and Chapter 7, respectively. In Chapter 8 numerical results using the developed material models are presented and compared to experimental ones. Conclusions of the research presented in Chapter 9.
7
CHAPTER 2
LITERATURE REVIEW
Numerous studies have been conducted on traditional sheet metal stamping methods such as mechanical stretch forming and deep drawing where both methods require a male and a female die for the proper forming of a finished part. In this section a review of the literature focusing on the most recent experimental and numerical studies conducted on sheet metal stamping, sheet hydroforming processes, wrinkling and warm forming of aluminum alloys will be presented.
2.1 Warm Forming Literature Review
Ayres (1979), Painter et al. (1980), Taleff et al. (1998 and 2001), Takata et al. (2000), Naka et al. (2001) and Li et al. (2003; 2004) studied the effects of elevated temperatures between 25°C – 410°C (77°F – 770°F) on the formability of aluminum alloys. Serrated flow behavior (dynamic strain aging) characteristic of AA5182 at room temperature (RT) was shown to diminish with elevated temperature deformation. Also, a larger total elongation in the post-uniform strain region was revealed along with a reduction in tensile strength. This indicates a shifting in the hardening behavior with temperature to a greater emphasis on strain-rate hardening relative to strain hardening.
8
Several research papers characterized the tensile behavior of aluminum alloys at elevated temperatures (Li et al., 2003; Ayres, 1979; Painter et al., 1980; Takata et al., 2000; Naka et al. 2001; Boogaard et al., 2001); the majority however, only studied the effect of elevated temperature on the evolution of the flow (hardening) stress. The effect of elevated temperature on the anisotropy of the material and how the yield surface of the aluminum alloy evolved as a function of temperature was not fully explored. In most cases, either Hill’s 1948 model (Hill, 1948) or the von Mises isotropic yield functions were used. Boogaard et al. (2001) characterized the behavior of AA5754-O in which two types of functions representing flow stress was used: the modified power law model and the Bergström model. The coefficients of the power law model were curve-fit exponentially as a function of temperature. The predictions of the material model, however, underestimated the values of the punch load in both models (PowerLaw and Bergström models). Keum et al. (2001; 2002) studied the effect of temperature on the anisotropic coefficients of the Barlat strain rate potential model (Barlat et al., 1993) for AA5052-H32. In this study, a third order polynomial function was used to fit the anisotropy coefficients, and a fifth order polynomial function was used to fit the coefficients for the power law model of flow stress. The material model, however, was assumed to reduce to the von Mises isotropic model as temperature increases. Furthermore, the plastic anisotropy parameters, R0, R45 and R90 almost remained constant with increasing temperature, which is opposite of what was found in the current research study. Other studies by Cleveland et al. (2003)
9
and Masuda et al. (2003) also describe temperature and strain rate effects on flow stress and deformation behavior. Naka et al. (2003) conducted experimental tests to determine the yield locus of AL5083 at elevated temperatures. It was shown that the yield locus changes as a function of temperature and a higher order yield function (e.g. Barlat YLD96) is best for representing such a material. Vial-Edwards (1997) show methods for determining the yield loci for FCC and BCC sheet metals. In warm forming methods, the physics of dislocation movement suggests a thermally activated process. Zener et al. (1944) studied the effects of strain rate and suggested a Zener-Hollomon parameter in which the relation between strain rate and temperature can be derived from statistical mechanics. However, as explained in Boogaard (2002), the Zener-Hollomon approach can only be used for small strain rates and temperature variations. When the strain rate itself is a function of temperature these types of models that incorporate the ZenerHollomon parameter are inappropriate for the simulation of warm forming of aluminum. Gronostajski (2000) provides a list of different types of deformation dependent flow stresses for FEM analysis. Håkansson et al. (2005) made a comparison of isotropic hardening and kinematic hardening in thermoplasticity. Prior research available for simulation of warm forming processes focuses only on the effect of elevated temperature on the evolution of the flow (hardening) stress. These include Li et al. (2003), Ayres and Wenner (1979), Painter et al. (1980), Takata et al. (2000), Naka et al. (2001) and Boogaard et al. (2001). The evolution of the yield surface of aluminum alloys as a function of
10
temperature and the effect on the anisotropy coefficients were not fully explored. In most cases, either Hill’s 1948 model (Hill, 1948) or the von Mises isotropic yield functions was used. Boogaard et al. (2001) characterized the behavior of AA5754-O for which two types of functions representing the flow stress were used: the modified power law model and the Bergström model. The yield surface used in this case was assumed to remain constant with respect to changing temperatures. Only the coefficients of the power law model were curve-fit exponentially as a function of temperature. The predictions of the material model, however, underestimated the values of the punch load in both models (PowerLaw and Bergström models). Čanađija et al. (2004) presented an associative coupled thermoplasticity model for J2 plasticity model to represent internal heat generated due to plastic deformation. In it, temperature-dependent material parameters developed were used.
2.2 Sheet Metal Hydroforming Literature Review
McClintock (1968) and Rice et al. (1969) conducted studies on sheet metals demonstrating a rapid decrease in fracture ductility as a hydrostatic pressure, applied across the material, was increased. Clift et al. (1990) and Hartley et al. (1992) demonstrated that for sheet metals, the use of a hydrostatic pressure prevented the initiation and spreading of micro-cracks within the metallic material. Yossifon and Tirosh (1977–1988) published a series of articles dealing with simple analysis of the hydroforming deep drawing process as applied to the
11
formation of cups from metallic materials such as copper, aluminum, steel and stainless steel. The goal of the studies was to establish a hydroforming fluid pressure path, relative to the punch stroke, that would prevent part failure due to rupture or wrinkling. Their earlier studies demonstrated the effect that excessive and insufficient fluid pressures have on the premature failure of hydroformed parts. The purpose of the later investigations was to determine a predetermined path that could be followed in order to produce parts that are free from these types of defects. In order to minimize wrinkling instabilities the fluid pressure was held to the possible minimum. The pressure relationship calculated by Yossifon and Tirosh (1977–1988), based on equating the bending energy of the buckled plate and the work done against lateral load (spring-type blankholder or fluid pressure) to the work done by the in-plane compressive membrane forces, included the governing parameters of friction coefficient and anisotropy. Through their work they were able to show that rupture instabilities occur when the fluid pressure being used for the hydroforming process was too high.
The fluid pressure
constrained the motion of the part and assisted the punch to deform the material. The fluid pressure to prevent rupture was evaluated in terms of average friction coefficient, material properties, and geometrical considerations. Using these two fluid pressure values a range was determined that allowed for the manufacture of parts without the occurrence of wrinkling or rupturing. This theory was tested experimentally and the results were very favorable with the predicted outcomes.
12
Lo, Hsu and Wilson (1993) expanded upon the earlier work of Yossifon and Tirosh by applying the deep drawing hydroforming theory to the analysis of the hemispherical punch hydroforming process. The purpose of this work was to determine a theoretical method of predicting failure due to wrinkling (buckling) or rupture (tensile instability) during the punch hydroforming of hemispherical cups. This work was basically an extension of the work done by Yossifon and Tirosh by incorporating a general friction-force expression into the analysis and expanding to more complicated geometries. In order to predict failure, the part was split into three regions based on the geometric characteristics of this operation. First there was a region where the part was free from contact with the die, a second region that consisted of the unsupported area termed the “lip area”, and the third region that was the area of the part that had already come into contact with the surface of the punch. Along with the determination of the failure areas, the study also attempted to identify an upper and lower bound for manufacturing, a region termed the “work zone”. Lo et al. (1993) proposed that if processes were run within these limits then there should be limited potential for failure.
They were able to conclude that the
working zone could be expanded by low friction forces, high strain hardening exponents, small drawing ratios, thick workpieces, and through the use of orthotropic materials. Hsu and Hsieh (1996) attempted to verify the theory developed by Lo et al. through a series of experimental procedures. The purpose was the validation and verification of the failure prediction method for wrinkling and rupture
13
instabilities during the punch hydroforming of sheet metal hemispherical cups. Various hydroforming pressure paths were tested during the process to validate the theory. They determined conclusively that a path that intersected the lower boundary of the working zone would lead to premature material failure due to wrinkling in every case. The same result was found for the pressure paths that intersected the upper boundary of the working zone. Through a series of varying parameter
experiments
the
results
achieved
experimentally
were
very
comparable to the theoretical predicted results. Gelin et al. (1994) experimentally and numerically studied the effects of process parameters during the aquadraw deep drawing process. The purpose of the study was to determine the main parameters that influence the aquadraw deep drawing process, specifically, the determination of the pressure in the cavity and under the blankholder as functions of process geometry, material parameters, and fluid parameters. Aquadraw deep drawing differs from hydroforming due to the use of a thin layer of water, subjected to fluid flow that replaces the thin rubber diaphragm between the material and the die cavity. The investigation, limited to axisymmetric sheet metal materials, proposed a cavity pressure modeling technique based on the optimal parameters of the process instead of being modeled by the Reynolds equation. A relationship to determine the cavity pressure was also derived based on the material behavior, the material thickness, the die entrance radius, and the drawing ratio. The paper evaluated the influence of each of these parameters on the overall cavity pressure. To demonstrate the effectiveness of these
14
parameters on the determination of the cavity pressure, the study referenced other researcher’s experiments.
Numerical analysis predictions of the deep
drawing process were in good agreement with the experimental behavior of the parts analyzed. Gelin et al. (1998) and Baida et al. (1999) both expanded upon the numerical work conducted in the Gelin et al. work dealing with the aquadraw deep drawing process. These two investigations expanded upon the numerical work by adding the process parameters monitoring, identification tools and general sensitivity analyses to the numerical method used as a predictor of the die cavity pressure during the deep drawing process. Overall their respective results showed very good correlation between the numerical and experimental behavior of the material. Shang et al. (1997) spent time on the evaluation of the copper spherical shell hydroforming process by studying the effects of intermittent draw-in during the operation. The purpose of this investigation was to examine, experimentally and numerically, the effects these intermittent changes would have on the formability of the blank material. During the processing of the cups there were two main formability factors that were investigated; the radius of the die shoulder and the blank holding force. Reducing the die shoulder radius increased formability but the use of a small radius had the potential of causing premature tearing of the blank along the die shoulder. Reducing the blank holding load encouraged drawin, inward flow of the flange material, thereby increasing the average thickness of the product and delayed the onset of material failure.
15
Since the radius of the die shoulder is normally fixed or limited by the product specifications then the logical approach to increasing formability would be to vary the blank holding load. During this study the copper material was formed into a nearly spherical shell using four different approaches. The first approach was a single-stage hydroforming process using two different deformation paths, one that allowed for the draw-in of the flange, and one that did not allow the draw-in to occur. The second approach evaluated the effect of a double-stage hydroforming process also using two different flow paths. The first path allowed for the draw-in during the first stage, and restricted it in the second. The second path was just the opposite, draw-in was not allowed during the first stage yet was permitted during the second stage. The results showed that during the singlestage hydroforming process, the formability of the material was greatly improved. For the double-stage hydroforming operation, the best results were achieved during the path that did not allow for the draw-in of the flange during the first stage, but did during the second stage.
2.3 Wrinkling Literature Review
Wrinkling in sheet metal forming, with tearing, is one of the most important instabilities that occur in parts formed using stamp forming and deep drawing processes. This phenomenon limits the type of parts and geometries that could be formed using these techniques. Simulation of wrinkling behavior using the finite element method (FEM) in sheet metal stamping is an important predictive tool. An accurate finite element model that could accurately predict the formation
16
of wrinkling could also be used at the tooling design stage of parts of various shapes. Many studies have been made to study the wrinkling behavior in sheet metal forming. These could be traced back to Geckler (1924), Baldwin et al. (1947), Senior (1981), Yoshida et al. (1981) and others. Triantafyllidis and Needlmen (1980) studied the effect of compressive bifurcation instabilities on the onset of flange wrinkling using the Swift cup test. Using numerical analysis linked with previous experimental work, they established the limiting drawing ratios (LDR), defined as the largest drawing ratio from which a cup can be drawn without fracture. The critical conditions governing the onset of wrinkling were studied. Many researchers have studied simulation of wrinkling behavior in sheet metals using the FEM (Finite Element Method) method. Doege et al. (1995) studied the necking and wrinkling behavior using several techniques. Necking, caused by tensile instability, was studied utilizing the results of the Continuum Damage Mechanics (CDM) and using the Gurson constitutive model. For studying the wrinkling behavior, they start by studying the buckling of onedimensional long column. The FEA method was used, where the problem was solved both implicitly and explicitly. Numerical results were then compared against experimental results of deep drawing of a 50mm cup. Necking behavior requires taking into account of the material microstructure, in particular microscopic processes that precede rupture. Therefore, the Gurson model, which accounts for the microscopic processes, was used.
17
Wrinkles that form during the sheet forming are due to internal compressive instabilities. Two types of wrinkles occur: (i)
wrinkles of first order in the flange; and
(ii)
wrinkles of second order in the free-forming (unsupported) zone between the punch radius and the die radius.
Since wrinkling is a problem of equilibrium state, the prediction of wrinkles is more difficult for implicit codes than for explicit codes. Introducing statistical geometrical imperfections in the blank is necessary in order to be able to simulate the wrinkling behavior using the implicit method. Boyce and Cao (1997), and Wang et al. (2000) studied the problem of wrinkling simulation using the implicit and explicit methods. Using the ABAQUSimplicit code the problem of forming a square cup was studied. Also Ls-Dyna explicit code was used. In order to solve the problem using the implicit code, initial imperfections had to be introduced into the blank mesh to take into account the instabilities. No imperfections were introduced in the explicit part of the analysis. Using these results they showed that the implicit FEM model with shell elements may overwhelmingly over-predict the failure heights, and the predictions of the explicit FEM models are sensitive to the selected critical wrinkling heights, the mesh density, the punch velocity, etc. Since wrinkling is a geometrical and a material dependent phenomenon, an anisotropic constitutive model was used. The anisotropic constitutive model introduced by Barlat et al. (1997) was used in this study. Prior to this work, finite
18
element analyses assumed material isotropy in their simulations, although the sheet metals usually showed anisotropic properties. Using experimental and numerical studies with ABAQUS of the cup drawing test using AA2008-T4 material, they showed the importance of using the correct constitutive model in representing material behavior in sheet metal forming processes using the explicit method. And the accuracy that could be achieved using this model versus experimental data was established. Several other studies have been also made, such as Kim et al. (2000, 20001), Kawaka et al. (2001). In these papers, verifying numerical data against experimental ones established the importance of using the correct constitutive model in the FEM analysis to represent anisotropic behavior of the material.
19
CHAPTER 3
THEORETICAL ANALYSIS
Accuracy of numerical analysis depends on the use of a constitutive model that precisely describes the behavior of the material. In previous research by the author (Abedrabbo et al., 2005c), the anisotropic behavior of AA6111-T4 aluminum alloy sheets was studied and the importance of using an appropriate material model that captures the anisotropic behavior was shown. In this research, two accurate anisotropic material models were used: Barlat’s YLD96 and Barlat’s YLD2000-2d. Description of both models is presented next.
3.1. Anisotropic Constitutive Models
3.1.1 Barlat YLD96
The YLD96 anisotropic yield function presented by Barlat et al. (1997a) is one of the most accurate yield functions for aluminum alloy sheets (Barlat et al. 2003), because it simultaneously accounts for yield stress and r-value directionalities. This yield function is based on a phenomenological description of the material, as an improvement to the YLD91 model (Barlat et al., 1991). To represent the behavior of these materials mathematically, a yield function was generalized to include the most general stress tensor with six components as shown in Equation (3.1)
20
Φ = α1 S2 − S3
a
a
+ α 2 S3 − S1 + α 3 S1 − S2
a
= 2σ a
(3.1)
with a = 6 and a = 8 for BCC and FCC materials, respectively. S1, S2 and S3 are the principal values of the stress tensor Sij, which is defined later. σ is the flow stress. The isotropic plasticity equivalent (IPE) stress is: (3.2)
S = Lσ
where, for orthotropic symmetry, L is defined below with anisotropy coefficients ck. ⎡ (c2 + c3 ) ⎢ 3 ⎢ ⎢ −c3 ⎢ 3 ⎢ Lij = ⎢ −c2 ⎢ 3 ⎢ 0 ⎢ ⎢ 0 ⎢ 0 ⎣
−c3 3 (c3 + c1) 3 −c1 3 0 0 0
−c2 0 3 −c1 0 3 (c1 + c2 ) 0 3 0 c4 0 0 0 0
0 0 0 0 c5 0
⎤ 0⎥ ⎥ ⎥ 0⎥ ⎥ 0⎥ ⎥ 0 ⎥⎥ 0⎥ ⎥ c6 ⎦
(3.3)
For the plane stress case (σz = σyz = σzx = 0), Equation (3.2) reduces to: ⎡ (c2 + c3 ) ⎢ 3 ⎡ Sx ⎤ ⎢ ⎢ S ⎥ ⎢ −c3 y ⎥ Sij = ⎢⎢ =⎢ 3 Sz ⎥ ⎢ ⎢ ⎥ ⎢ −c2 ⎢⎣Sxy ⎥⎦ ⎢ 3 ⎢ 0 ⎣
−c3 3 (c3 + c1) 3 −c1 3 0
−c2 ⎤ 0⎥ 3 ⎥⎡ σx ⎤ ⎢ ⎥ −c1 0 ⎥⎥ ⎢ σ y ⎥ 3 ⎥⎢ 0 ⎥ ⎢ ⎥ (c1 + c2 ) 0 ⎥ ⎢σ ⎥ ⎥ ⎣ xy ⎦ 3 ⎥ c6 ⎦ 0
(3.4)
The principal values of Sij, as defined in Equation (3.1), are found as follows:
21
2
S x + Sy ⎛ S x − Sy ⎞ 2 S1,2 = ± ⎜⎜ ⎟⎟ + Sxy 2 2 ⎝ ⎠ and S3 = -(S1+S2)
(3.5)
because of the deviatoric nature of Sij. The anisotropy
coefficients αi of the yield function (3.1) are defined as:
α1 = α x cos2 β + α y sin2 β α 2 = α x sin2 β + α y cos2 β
(3.6)
α3 = α z0 cos2 2β + α z1 sin2 2β In the above equations, c1, c2, c3, c6, αx, αy, αz0 and αz1 are coefficients that describe the anisotropy of the material. The parameter 2β represents the angle between the line OA and the axis of the principal IPE stress, S1, as shown in Figure 3.1, and is equal to: ⎛ 2Sxy 2β = tan−1 ⎜ ⎜ S x − Sy ⎝
τ
2β
⎞ ⎟ ⎟ ⎠
A
⎫ ⎬ Sxy ⎭
Sy P S2
(3.7)
O
S1
σ
Sx
Figure 3.4 Determination of the anisotropy parameter 2β from Mohr’s circle.
22
With this yield function it is possible to increase the yield stress at pure shear without increasing the other plane strain yield stresses. It should be noted that by setting
the
anisotropy
coefficients
c1=c2=c3=c6=α1=α2=α3=1.0
and
a=2
(quadratic), the von Mises isotropic yield function is recovered.
3.1.2 Barlat YLD2000-2d
Although YLD96 (Barlat et al., 1997) is one of the most accurate anisotropic yield functions for aluminum and its alloys (Lademo, 1999; Lademo et al., 1999) where it takes seven parameters into account in the plane stress condition. There are three problems associated with Yld96 (Barlat et al., 2003) with respect to FE simulations:
There is no proof of convexity, which is an important requirement in numerical simulations to ensure the uniqueness of a solution.
The derivatives are difficult to obtain analytically which is inconvenient again for FE simulations.
Numerical problems solved using the full stress state formulation of YLD96, which might be difficult to solve because of the relative complexity of the Yld96 formulation, have been encountered (Becker, 1998; Szabo, 2001).
Therefore, a better incompressible anisotropic plasticity formulation that can guarantee convexity, make FE implementation and application simpler, and take
σ0, σ45, σ90, r0, r45, r90 and σb into account for plane stress was developed by Barlat et al. (2003).
23
Starting with the function proposed by Hershey (1954) and Hosford (1972) Φ1 = s1 − s2
a
+ s2 − s3
a
a
+ s3 − s1 = 2σa
(3.8)
an isotropic yield function that reduces to Equation (3.8) can be simply written as φ = φ′ + φ′′ = 2σa
(3.9)
where a φ′ = s1 − s2
(3.10)
a a φ′′ = 2s2 + s1 + 2s1 + s2
(3.11)
and
Because a plane stress state can be described by two principal values only, φ′ and φ′′ are two isotropic functions since it is possible to permute the (in-plane)
indices 1 and 2 in each function. For the anisotropic case, a linear transformation reduces to ⎡ X′ ⎤ ⎡ C′ ′ C12 0 ⎤ ⎡sxx ⎤ ⎢ xx ⎥ ⎢ 11 ⎢ ⎥ 0 ⎥⎥ ⎢syy ⎥ ⎢ X′yy ⎥ = ⎢C′21 C′22 ⎢ ⎥ ⎢ ⎢ ⎥ 0 C′66 ⎥⎦ ⎢s ⎥ ⎢⎣ X′xy ⎥⎦ ⎣ 0 ⎣ xy ⎦
(3.12)
⎡ X′′ ⎤ ⎡ C′′ ′′ C12 0 ⎤ ⎡sxx ⎤ ⎥ ⎢ xx ⎥ ⎢ 11 ⎢ 0 ⎥⎥ ⎢syy ⎥ ⎢ X′′yy ⎥ = ⎢C′′21 C′′22 ⎥ ⎢ ⎥ ⎢ ⎢ 0 C′′66 ⎥⎦ ⎢s ⎥ ⎢⎣ X′′xy ⎥⎦ ⎣ 0 ⎣ xy ⎦
(3.13)
and
or, using
24
⎡ 2 3 −1 3 0 ⎤ T = ⎢⎢ −1 3 2 3 0 ⎥⎥ ⎢⎣ 0 1⎥⎦ 0
(3.14)
X′ = C′.s = C′.T.σ = L′.σ X′′ = C′′.s = C′′.T.σ = L′′.σ
(3.15)
they become
The anisotropic yield function is given by Equation (3.9), where a φ′ = X1′ − X′2
(3.16)
a a φ′′ = 2X′′2 + X1′′ + 2X1′′ + X′′2
(3.17)
and
reduces to the isotropic expression when the matrices C′ and C′′ are both taken as the identity matrix so that X′ = X′′ = s . Because φ′ depends on X1′ − X′2 , only three coefficients are independent in C′ . For convenience, the coefficients of L′ and L′′ can be expressed as follows ′ ⎤ ⎡23 0 ⎡ L11 ⎢ L′ ⎥ ⎢ 0 ⎢ 12 ⎥ ⎢ − 1 3 ⎢ L′21 ⎥ = ⎢ 0 −1 3 ⎢ ⎥ ⎢ 23 ⎢L′22 ⎥ ⎢ 0 ⎢⎣L′66 ⎥⎦ ⎢⎣ 0 0
0⎤ 0 ⎥⎥ ⎡ α1 ⎤ ⎢ ⎥ 0 ⎥ ⎢α2 ⎥ ⎥ 0 ⎥ ⎢⎣α7 ⎥⎦ 1⎥⎦
(3.18)
and ′′ ⎤ ⎡ L11 ⎡ −2 2 8 −2 ⎢L′′ ⎥ ⎢ ⎢ 12 ⎥ 1 ⎢ 1 −4 −4 4 ⎢ L′′21 ⎥ = ⎢ 4 −4 −4 1 ⎢ ⎥ 9⎢ 2 −2 ⎢L′′22 ⎥ ⎢ −2 8 ⎣⎢ 0 0 0 0 ⎣⎢L′′66 ⎦⎥
25
0 ⎤ ⎡ α3 ⎤ ⎢ ⎥ 0 ⎥⎥ ⎢α 4 ⎥ 0 ⎥ ⎢ α5 ⎥ ⎥⎢ ⎥ 0 ⎥ ⎢α6 ⎥ 9 ⎦⎥ ⎣⎢ α8 ⎥⎦
(3.19)
where all the independent coefficients αk (for k from 1 to 8) reduce to 1 in the isotropic case. Only seven coefficients are needed to account for the seven input data mentioned above. There are several possibilities to deal with the eighth ′′ = C′′21 or L12 ′′ = L′′21, or use another input coefficient, for instance, assuming C12
data. This additional input data can be the ratio rb = ε yy ε xx , which characterizes the slope of the yield surface in balanced biaxial tension ( σ yy = σ xx ). This parameter, which is denoted rb by analogy with the r value obtained in uniaxial tension, can be determined with three different methods: Experimentally measured, calculated with another yield function, for instance Yld96, or computed from a polycrystal model if the crystallographic texture of the material is known. The principal values of X′ and X′′ are ⎛ X1 = 1 ⎜ X xx + X yy + 2⎜ ⎝
( Xxx − Xyy )
⎛ X2 = 1 ⎜ X xx + X yy − 2⎜ ⎝
( Xxx − Xyy )
2
⎞ + 4X2xy ⎟ ⎟ ⎠
(3.20)
⎞ + 4X2xy ⎟ ⎟ ⎠
(3.21)
and 2
with the appropriate indices (prime and double prime) for each stress. Assuming the associated flow rule, the normal direction to the yield surface, which is needed to calculate the strain rates (or strain increments), is given by
∂φ ∂φ′ ∂X′ ∂φ′′ ∂X′′ ∂φ′ ∂φ′′ • • • L′ + • L′′ = + = ∂σ ∂X′ ∂σ ∂X′′ ∂σ ∂X′ ∂X′′
26
(3.22)
In the calculation of the derivatives, there are two singular cases, namely when X1′ = X′2 and X1′′ = X′′2 . However, the normal directions to the yield surface can still be obtained for these two special cases.
3.2. Plastic Anisotropy Parameters
The plastic anisotropy parameter Rθ is defined as the ratio of width-tothickness increments: dε Rθ = w dε t
(3.23)
The thickness strain, however, is difficult to measure accurately in a thin sheet. Thickness strains are instead calculated from measurements of the longitudinal and width strains assuming constancy of volume as follows: d ε t = −(d ε l + d εw )
(3.24)
Therefore, two strain rate measurements are required, namely the longitudinal (dεl) and the width strains (dεw), in order to calculate the plastic anisotropy parameters. For isotropic materials, R-values are equal to 1.0 for any direction θ. R-values not equal to 1.0 indicate that plastic anisotropy exists in the material. A high Rvalue suggests that the material has a high resistance to thinning and thickening, which implies better formability of the material. If R-values depend on θ, then the material is planar anisotropic, otherwise it is planar isotropic.
27
3.3. Anisotropy Coefficients Calculation
3.3.1 Barlat YLD96
In order for the constitutive model to accurately represent the aluminum alloy sheet at large strains, the anisotropic coefficients describing the behavior of the material, c1, c2, c3, c6, αx, αy, and αz1 must be calculated (αz0 can be set equal to 1.0) using seven test results. αz0 can be set to 1.0 because increasing αz0 will decrease αx and αy simultaneously as these parameters are related. Four stress states, uniaxial tension along the rolling, 45°, and transverse directions and balanced biaxial stress (Green et al., 2004), provide the required seven data points (Barlat et al. 1997b). A description of the data points from experiments with respect to the yield function is shown in Figure 3.2. A summary of the experimental data needed to calculate the coefficients is given in Table 3.1.
Table 3.1 Experimental data needed to calculate yield function coefficients for YLD96 in the plane stress state. Test
Orientation
Balanced Biaxial
N/A
Uniaxial
0°
45°
90°
Flow Stress
σb
σ0
σ45
σ90
R-Value
N/A
R0
R45
R90
28
σ2
Uniaxial tension
Plane strain εx=0 Balanced biaxial tension
+Y
45°
−σ1
Plane strain εy=0
. -X
+X
σ1 Uniaxial tension
Pure shear -Y
−σ2
Figure 3.5 Characterization of required tests to be used in calculation of the anisotropic values in relation to the yield locus. The yield stress could be used as input data instead of the flow stress. However, it may be difficult to accurately measure yield stress in the bulge test as well as in uniaxial tension because the slope of the stress-strain is steep and yielding is not always a discrete event in many aluminum alloys. Also, the yield stress is associated with very small plastic strain and might not reflect the anisotropy of the material over a larger strain range. For these reasons, flow stresses at equal amount of plastic work (Wb=Wu) were selected as input data rather than yield stress, as shown in Figure 3.3.
29
In this research, a value of Wb= Wu= 20 MPa per unit volume was used to extract flow stresses σ 0 , σ 45 and σ 90 , as the material exhibited a sufficient amount of plastic strain (~18%) at that level of plastic work. The normalized values σ 0 σ 0 , σ 45 σ 0 and σ 90 σ 0 were then used as input to calculate the anisotropy coefficients of the yield function. Flow stresses extracted at different values of W would not significantly change those stress ratios, and therefore consistent results can be obtained as long as sufficiently large values of W is used.
Uniaxial Tension Test
σb
Longitudinal Stress
Membrane Stress
Bulge Test
σu
Longitudinal plastic strain
Thickness plastic strain
Figure 3.6 Flow stresses at equal amounts of plastic work (Wb=Wu). The anisotropy coefficients were calculated as follows. For uniaxial tension in the rolling direction (0°), the principal IPE stresses become:
30
(c2 + c3 ) σ0 3 c = − 3 σ0 3 c = − 2 σ0 3
S1 = S2 S3
(3.25)
The yield function (3.1) can then be written as: F1 = α x c2 − c3
a
+ α y 2c2 + c3
a
a
+ α z0 c2 + 2c3
a
⎛ 3σ ⎞ − 2⎜ ⎟ =0 ⎝ σ0 ⎠
(3.26)
Similar equations could be written for uniaxial tension in the 45°-direction, the transverse direction (90°) and the balanced biaxial tension as follows: F2 = α x 2c1 + c3
a
⎛ 3σ ⎞ + α y c3 − c1 + α z0 2c3 + c1 − 2 ⎜ ⎟ =0 ⎝ σ 90 ⎠
a
a
a
a
2
⎛c −c ⎞ F3 = α1 ( c1 + c2 ) − ⎜ 2 1 ⎟ + (2c6 )2 ⎝ 3 ⎠
+ a
2
⎛c −c ⎞ α 2 (c1 + c2 ) + ⎜ 2 1 ⎟ + (2c6 )2 ⎝ 3 ⎠
α3
a
2
⎛ c2 − c1 ⎞ 2 ⎜ 3 ⎟ + (2c6 ) ⎝ ⎠
(3.27)
+
(3.28)
a
⎛ 2σ ⎞ − 2⎜ ⎟ =0 ⎝ σ 45 ⎠
where α1, α2 and α3 are As defined by Equation (3.6). F4 = α x 2c1 + c2
a
a
⎛ 3σ ⎞ + α y 2c2 + c1 + α z0 c2 − c1 − 2 ⎜ ⎟ =0 σ b ⎝ ⎠ a
31
a
(3.29)
The other necessary sets of equations are derived from the plastic anisotropy definitions for R0, R45 and R90, as follows: Rθ = −
d εw (d ε w + d ε l )
(3.30)
From the normality rule:
εxx = λ εyy = λ
∂Φ ∂σ xx ∂Φ
(3.31)
∂σ yy
εxy + εyx = λ
∂Φ ∂σ xy
where λ is a proportionality factor. After substitution, the plastic anisotropy is found as follows. For uniaxial tension performed along an arbitrary orientation θ measured counterclockwise from the x-axis, the ratio of the strain rates Rθ is defined by
∂Φ Rθ =
∂Φ / ∂σ 90 +θ ∂σ 90 +θ =− ∂Φ / ∂σ zz ⎛ ∂Φ ∂Φ ⎞ + ⎜ ⎟ ⎝ ∂σ 90 +θ ∂σθ ⎠
(3.32)
where
∂Φ ∂Φ ∂σ xx ∂Φ ∂σ yy ∂Φ ∂σ xy = + +2 ∂σθ ∂σ xx ∂σθ ∂σ yy ∂σθ ∂σ xy ∂σθ
(3.33)
∂Φ ∂Φ ∂Φ ∂Φ = cos 2 θ + sin 2 θ + 2 cos θ sin θ ∂σθ ∂σ xx ∂σ yy ∂σ xy
(3.34)
or
32
and
∂Φ ∂σ 90 +θ
=
∂Φ ∂Φ ∂Φ sin 2 θ + cos 2 θ − 2 cos θ sin θ ∂σ xx ∂σ yy ∂σ xy
(3.35)
which leads to
Rθ = −
∂Φ ∂Φ ∂Φ sin 2 θ + cos 2 θ − 2 cos θ sin θ ∂σ xx ∂σ yy ∂σ xy ⎛ ∂Φ ∂Φ + ⎜ ⎜ ∂σ xx ∂σ yy ⎝
⎞ ⎟ ⎟ ⎠
(3.36)
where θ = 0°, 45° and 90°, and Φ is the yield function. These systems of equations were solved for the values of the anisotropy coefficients c1, c2, c3, c6, αx, αy, and αz1 using a non-linear solver (e.g. NewtonRaphson) with initial values corresponding to the isotropic situation (c1 = c2 = c3 = c6 = 1.0).
33
3.3.2 Barlat YLD200-2d
Coefficients α1 to α6
Three stress states, namely uniaxial tension along the rolling and the transverse directions, and the balanced biaxial stress state, provide six data points, σ0, σ90, σb, r0, r90, and rb. rb defines the slope of the yield surface at the balanced biaxial stress state ( rb = ε yy ε xx ). This ratio can be evaluated by performing compression of circular disks in the sheet normal direction and measuring the aspect ratio of the specimen after deformation. rb can also be estimated by calculations using either a polycrystal model or the yield function Yld96. The loading for each stress state can be characterized by the two deviatoric components, sx = γσ and s y = δσ . There are two equations to solve per stress state, one for the yield stress and the other for the r value. a
F = φ − 2 (σ σ) = 0
(3.37)
which satisfies the yield stress, and
G = qx
∂φ ∂φ − qy =0 ∂sxx ∂syy
(3.38)
which satisfies the r value. Function ф can be rewritten as a
a
a
a
φ = α1γ − α 2δ + α3 γ + 2α 4δ + 2α5 γ + α 6δ − 2 ( σ σ ) = 0
34
(3.39)
where γ, δ, σ, qx and qy for the tests mentioned above are given in the following in Table 3.2.
Table 3.2 Experimental data needed to calculate yield function coefficients for YLD20002d. Test
γ
δ
0° tension
2/3
-1/3
90° tension
-1/3
Balanced biaxial tension
-1/3
σ
qx
qy
0
1-r0
2+r0
2/3
90
2+r0
1-r90
-1/3
b
1+2rb
2+rb
The six coefficients αk can be computed by solving the two Equations (3.37) and (3.38). These equations were obtained using the linear transformations with the Ckl (prime and double prime) coefficients, which are related to the αk coefficients in the following way: ′ α1 = C11 α2 = C′22 ′′ α3 = 2C′′21 + C11 2α 4 = 2C′′22 + C1′′ 2 ′′ + C′′21 2α5 = 2C11
(3.40)
′′ + C′′22 α6 = 2C12 The αk coefficients were used because the yield function reduces to its isotropic form when all these coefficients are 1. combinations of the αk .
35
The Ckl and Lkl coefficients are linear
Coefficients α7 and α8
Uniaxial tension tests loaded at 45° to the rolling direction give two data points, σ = σ45 and r = r45 . The stress state is on the yield surface if the following equation is satisfied:
F=
k′22 + 4α72 2
a
3k1′′ − k ′′22 + 4α82 + 4
3k1′′ + k ′′22 + 4α82 4
a
+
(3.41)
a a
− 2 ( σ σ45 ) = 0
where α − α2 k′2 = 1 3 2α5 + α6 + α3 + 2α 4 k1′′ = 9 2α5 + α6 − α3 − 2α 4 k′′2 = 3
(3.42)
The associated flow rule can be expressed as: G=
∂φ ∂φ 2aσa + − =0 ∂σ xx ∂σ yy σ (1 + r45 )
(3.43)
The Newton-Raphson numerical procedure is used to solve for the eight αk coefficients simultaneously. The two matrices L′ and L′′ are completely defined with these eight coefficients.
36
3.4. Constitutive Equations (Flow Stress)
Flow stress ( σ ) represents the radius of the yield function during deformation. Metals undergoing plastic deformation at high temperatures and different strain rates should be modeled according to the physical behavior of the material (Gronostajski, 2000). An appropriate constitutive equation describing changes in the flow stress of the material depends on deformation conditions such as temperature and strain rate. Wagoner et al. (1988) proposed a new flow rule that includes the strain-rate sensitivity: ⎛ ε ⎞ σ (ε p ,ε ) = K (ε p + ε 0 )n ⎜ ⎟ ⎝ ε sr 0 ⎠
m
(3.44)
where K (strength hardening coefficient), n (strain-hardening exponent) and m (strain-rate sensitivity index) are material constants. ε p is the effective plastic strain and ε is the strain rate. ε0 is a constant representing the elastic strain to yield and ε sr 0 is a base strain rate (a constant). This model was primarily selected over other types of hardening laws (e.g. Voce) because it represented the hardening behavior of the current material accurately, and it incorporates strain rate effects. Gronostajski (2000) describes other types of hardening laws that could be used (e.g. Backofen, Grosman) to represent other materials, including those with hardening saturation behavior. The coefficients K, n and ε0 of Equation (3.44) were solved for simultaneously using Levenberg-Marquardt algorithm (Levenberg, 1944; Marquardt, 1963). The
37
strain-rate sensitivity index m was determined by the following equation using experimental data m=
ln(σ 2 /σ1) ln(ε2 / ε1)
(3.45)
From the data acquired for ε0, it was found that its value is very small (as will be seen later). Another method was also used to solve for ε0, which represents the elastic strain to yield by solving for the intersection of the linearly elastic loading equation with the strain hardening equation:
σ = Eε σ = Kε n
(3.46)
This provides the elastic strain at yield as: ⎡ 1 ⎤
⎥ ⎛ E ⎞⎢ ε 0 = ⎜ ⎟ ⎣ n −1⎦ ⎝K ⎠
(3.47)
where E is Young’s modulus. Comparing the values of ε0 calculated using Equation (3.47), to the ones determined using the Levenberg-Marquardt algorithm, it was found that the differences between them are very small. Also the values of ε0 are generally too small to have a great effect on resulting stresses in forming processes. Furthermore, it was observed that ε0 does not show significant temperature dependence. In this research ε0 calculated from Equation (3.47) was therefore used. Material properties change with increasing temperature. To include temperature effects in the flow rule, the material constants K, n and m are expressed as a function of temperature, and therefore the flow stress equation becomes: 38
⎛ ε ⎞ σ (ε p , ε,T ) = K (T )(ε p + ε 0 )n(T ) ⎜ ⎟ ⎝ ε sr 0 ⎠
m(T )
(3.48)
The hardening model described by Equation (3.48) was used assuming isotropic hardening behavior.
39
CHAPTER 4
MATERIAL CHARACTERIZATION
4.1. Materials
Chemical composition and initial mechanical properties of the aluminum alloy (AA3003-H111) used in this study are shown in Tables 4.1 and 4.2. All samples were taken from a single lot of material. Because AA3003 is very soft in the annealed condition, one cold roll pass (H111) was applied to these sheets in order to provide some rigidity during handling.
Table 4.1 Chemical composition of AA3003-H111 (wt. %). Al Mg Cu Mn Fe Si bal
0.02
0.07
1.10
0.50
0.21
Cr
Zn
Ti
Ni
0.005
0.01
0.02
0.005
Table 4.2 As-received mechanical properties of AA3003-H111. Thickness UTS Yield Strength (mm) (MPa) (MPa) 0.96
114.0
63.0
40
Total Elongation (%) 32.0
4.2. Experimental Procedure
Uniaxial tests with standard ASTM-E8 rectangular dog-bone shaped samples shown in Figure 4.1, were performed on an Instron Model 1127, screw-driven frame with a 4.5 kN load cell, 25 mm axial extensometer (50% strain max) and 12.7 mm transverse extensometer (30% strain max).
50
200 50
12.5
20
Figure 4.46 Rectangular tension test specimen, ASTM-E8. All dimensions in mm. The tensile samples were prepared from the aluminum sheet metal at 0° (rolling direction, RD), 45° and 90° (transverse direction, TD) from the rolling direction of the sheet, as shown in Figure 4.2.
41
90º to Rolling Direction
Rolling Direction (0º)
Figure 4.47 Tensile specimen location with respect to sheet. For the measurement of plastic anisotropy parameters, ASTM-E517 specifies that the test be performed at a 0.5/min (0.0083 s-1) straining rate. Sample temperature was controlled with an Instron Model 3119 oven with a convection heating system. Tests were performed for several elevated temperatures in the range of 25°C – 260°C (77°F – 500°F), with the results of three tests averaged for each temperature. It should be mentioned that only a small variation in the stress-strain behavior was noticed between the three tests for each temperature. Tests in the experiments were performed at a fixed ambient temperature. To study the strain-rate sensitivity of the material, uniaxial tests under several strain rates (0.001 s-1 – 0.08 s-1) were performed at each temperature. Bulge testing was performed on sheet samples from the same lot of material in order to determine σb at room temperature. At the moment, experimental calculation of the balanced biaxial data at elevated temperature is under 42
development. Based on experimental observations from uniaxial tests at 0°, 45° and 90°, see Figure 4.3, it was noticed that at room temperature the stress values at 45° and 90° coincide with the stress value from bulge test. Since stress values at 45° and 90° continued to coincide at elevated temperatures, it was assumed that a good estimate for the bugling stress at elevated temperature would be the average of the stress values at 45° and 90°. This assumption is unproven and will be studied and reported separately in a future publication. It should be noted that other methods could be used to predict σb. For example, bulge stress σb could be calculated from Hill’s 1948 model (Hill,
1948) or YLD91 model (Barlat et al., 1991) using only the available uniaxial tension tests data. In the current research however, it was not deemed appropriate to present this approach until after experimental bulge tests are completed, at which time a comparison of all these methods will be reported.
43
160 140
True Stress (MPa)
120 100 80 60 40 20 σ0
σ45
σ90
σb
0 0
50
100
150
200
250
300
Temperature (°C)
Figure 4.48 Stress values at equal values of plastic work (20 MPa/unit-vol) as a function of temperature. Data shown for tests at 0°, 45° and 90° directions at multiple elevated temperatures and for balanced biaxial tension (bulge) data at room temperature.
44
4.3. Hardening Model
Figure 4.4 shows the true stress-true strain behavior at room temperature for rolling direction (RD), 45° and TD uniaxial tests as well as the balanced biaxial (bulge) test. 200 180 160
Balanced Biaxial
True Stress (MPa)
140 120 100 80 60 40
Rolling Direction
20
45°-Direction Transverse Direction
0 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
True strain
Figure 4.49 True stress strain data of AA3003-H111 at room temperature for uniaxial tests at 0°, 45° and 90° directions and for balanced biaxial tension (bulge) data. As seen from the bulge test, much larger strains can be obtained compared to uniaxial strains. Figure 4.5 shows the engineering stress-strain behavior at several temperatures for the rolling direction. As temperature increased, flow stress of the material decreased with a corresponding increase in the elongation to failure.
45
120 25°C 100 66°C 121°C Tensile Stress ( MPa)
80 177°C 204°C
60
232°C 40
20
0 0
10
20
30
40
50
Tensile Strain (%)
Figure 4.50 Engineering stress strain curves of AA3003-H111 at several elevated temperatures for the rolling direction. Two methods were used to measure strain-rate sensitivity of the material. In the first, several tensile tests were performed on identical samples at different strain-rates (from 0.001 s-1 to 0.08 s-1) and at several elevated temperatures. For example, Figure 4.6 shows the true stress-true strain curves at 204°C (400°F) in the rolling direction.
46
-1
0.08s
100
-1
0.05s
-1
0.01s
True Stress (MPa)
80
0.001s
-1
60
40
20
0 0.00
0.05
0.10
0.15
0.20
0.25
0.30
0.35
True Strain
Figure 4.51 True stress strain data of AA3003-H111 from uniaxial tests at 204°C (400°F) at several strain rates for the rolling direction. Considering that the sample-to-sample variation (thickness, width, etc…) may cause scatter in the results, it was desirable to perform the strain-rate sensitivity analysis using a single sample with the “Jump-rate Test” method (Wagoner et al., 1996). In this method, the crosshead speed is increased to produce a “jump” in the strain rate at some predetermined level of strain. Figure 4.7 shows results of this method at several elevated temperatures in which crosshead speed jumps were 10-50-150 mm/mm/min. As seen in Figures 4.6 and 4.7, AA3003-H111 exhibits very small strain-rate sensitivity at room temperature, but with increasing temperature, the material becomes more strain rate sensitive.
47
160
150mm/mm/min
140
25°C 93.3°C
50mm/mm/min
120 True Stress (MPa)
10mm/mm/min 100 204.4°C 80 260°C 60 40 20 0 0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
True Strain
Figure 4.52 True stress strain data of AA3003-H111 from uniaxial tests using the jump-rate strain method. The jumps in crosshead speed were 10-50-150 mm/mm/min. (at temperature of 260 °C, at cross head speed of 150 mm/mm/min the specimen failed, therefore there is no information available). From the results of uniaxial tension tests performed at different constant strain rates, and at several elevated temperatures, values for the Holloman hardening rule (K, n, ε0 and m) and the plastic anisotropy parameters (Rθ) were calculated as a function of temperature. For example, Table 4.3 shows values of K, n, ε0 obtained with Levenberg-Marquardt algorithm, and m for the rolling direction. Similar results were also found for 45° and the transverse directions. For the finite element simulation, values of K, n and m only from the rolling direction were used to represent isotropic hardening. (The two other results, found for 45° and the transverse directions, could be used on the condition that
48
the material axes definition of the anisotropic material in the FEA coincides with the desired test direction).
Table 4.3 Material properties of AA3003-H111 at elevated temperatures. The values shown below represent an average of three tensile tests. (The hardening values are for the rolling direction). Temperature K n ε0 m R0 R45 R90 °C MPa 25
199.82
0.215
8.30E-04
0.003
0.827
1.126
0.773
38
186.41
0.200
5.00E-04
0.004
1.074
1.424
0.892
66
175.78
0.187
3.16E-04
0.004
1.092
1.457
0.978
93
168.41
0.179
6.64E-04
0.005
1.247
1.690
0.879
121
146.98
0.175
3.20E-04
0.010
1.262
1.742
0.982
149
139.18
0.163
7.40E-04
0.015
1.273
1.787
0.943
177
119.65
0.157
6.20E-04
0.030
1.413
1.975
1.007
204
106.32
0.137
2.21E-04
0.045
1.614
1.977
1.105
232
93.82
0.132
3.71E-04
0.065
1.688
2.281
1.236
260
77.32
0.116
5.05E-04
0.080
2.035
2.485
1.328
Figures 4.8 & 4.9 show that the values of K and n decrease linearly with temperature.
49
250
K-Values (MPa)
200
150 K(T) = -0.5058*T + 210.4 100
50
Experimental Data
Curve Fit
0 0
50
100
150
200
250
300
Temperature (°C)
Figure 4.53 Strength hardening coefficient (K) of AA3003-H111 along the rolling direction as a function of temperature at 0.5 mm/mm strain rate. Solid line is a linear curve fit.
50
0.25
0.20
n(T) = - 0.0004*T + 0.2185 n-Value
0.15
0.10
0.05
Experimental Data
Curve Fit
0.00 0
50
100
150
200
250
300
Temperature (°C)
Figure 4.54 Strain-hardening exponent (n) of AA3003-H111 along the rolling direction as a function of temperature at 0.5 /min strain rate. Solid line is a linear curve fit. This is in line with the decreasing flow stress and flattening of the hardening curve shown in Figure 4.5. Figure 4.10 shows the variation of the strain-rate sensitivity index, m, as a function of temperature for AA3003-H111. At lower temperatures, m values are very small, indicating the material is strain-rate insensitive; but at higher temperatures, the material exhibits a significant sensitivity to strain-rate. An exponential function, as shown on the graph, was used to represent the behavior of m as a function of temperature.
51
0.09 0.08 0.07 0.06 0.0147*T
m-value
m(T) = 0.0018*e 0.05 0.04 0.03 0.02 0.01
M-value Expon. (M-value)
0.00 0
50
100
150 Temperature (°C)
200
250
300
Figure 4.55 Strain-rate sensitivity index (m) of AA3003-H111 along the rolling direction as a function of temperature. Solid line is an exponential curve fit.
Figure 4.11 shows the variation of the plastic anisotropy parameters R0, R45 and R90, measured at a straining rate of 0.5/min, with respect to temperature. As mentioned earlier, R-values higher than 1.0 indicate good formability and resistance to thinning. As can be seen from Figure 4.11, the values of R0, R45 and R90 increase with temperature, which suggests that the formability of the aluminum sheet also enhances at elevated temperature. Table 4.4 shows a summary of the equations used to fit the hardening parameters for the flow rule, along the rolling direction, as a function of temperature. Curve fitting at 45° and the transverse directions are not required, since only the values of hardening in the rolling direction are used in the finite element simulation.
52
3.0
2.5
R - Value
2.0
1.5
1.0
0.5
R0
R45
R90
0.0 0
50
100
150
200
250
300
Temperature °C
Figure 4.56 Plastic anisotropy parameters (Rθ) of AA3003-H111 as a function of temperature, calculated at 0.5/min strain rate (ASTM-E517).
Table 4.4 Summary of equations used to fit hardening parameters for the power law flow rule. Temperatures in °C. Hardening Rolling Direction, 0° Unit Parameter K(T)
= -0.5058*T + 210.40
n(T)
= -0.0004*T + 0.2185
m(T)
= 0.0018*exp(0.0147*T)
53
MPa
4.4. Barlat’s Yield 96 Anisotropy Coefficients Calculation
Results from the uniaxial and bulge tests were implemented into Equations (3.26―3.36), and Barlat’s YLD96 model coefficients were calculated at several temperatures with the results shown in Table 4.5. It should be noted that in Barlat’s yield function αz0 = 1.0, and the non-quadratic parameter for aluminum alloys (FCC) is a = 8.
Table 4.5 Barlat’s YLD96 model anisotropy coefficients calculated at several elevated temperatures. Temperature c2 c3 c6 αx αy αz1 c1 °C 25
1.1169
0.9545
1.0030
1.0429
0.9130
1.3960
1.250
38
1.1622
1.0424
0.9877
1.0985
0.6516
0.7866
1.040
66
1.2299
1.0232
0.9968
1.0888
0.5028
0.8513
1.156
93
1.2871
1.0478
0.9953
1.0721
0.4810
0.7034
1.508
121
1.0847
1.0171
1.0116
0.9639
0.8060
0.7864
2.110
149
1.1804
1.0318
1.0049
1.0185
0.6182
0.7380
1.806
177
1.1280
1.0619
1.0003
0.9782
0.6650
0.5930
2.146
204
1.1742
1.0681
1.0058
0.9932
0.5382
0.5294
1.901
232
1.2526
1.2095
0.9564
0.9939
0.3160
0.2427
1.969
260
1.2329
1.0994
0.9939
0.9904
0.3920
0.4515
2.123
Figure 4.12 is a plot of the yield function for AA3003-H111, comparing Barlat’s YLD96 model at 25°C (77°F) with the isotropic von Mises yield function (normalized stress values with respect to σ ).
54
σ / ⎯σ x
Figure 4.57 2D-Plot of von Mises yield function vs. Barlat’s YLD96 yield function at 25°C (77°F) for AA3003-H111. Figure 4.13 is a similar plot comparing Barlat’s YLD96 function at several elevated temperatures (normalized stress values with respect to σ ), indicating that temperature has an effect on the shape of the yield surface (change in the local slopes at different temperatures). This change in the shape is particularly noticed near the balanced biaxial region of the yield surface. This is to be expected, since any change in the plastic anisotropy parameters (R-values) would affect the shape of the yield surface shape (Barlat et al., 1997a, b).
55
1.5 1
y
σ /⎯σ
0.5 0
-0.5 25°C 93°C 149°C 204°C 260°C
-1 -1.5 -1.5
-1
-0.5
0
0.5
1
1.5
σ /⎯σ x
Figure 4.58 2D-Plot of Barlat’s YLD96 yield function for AA3003-H111 at several elevated temperatures. Normalized stress values. Figure 4.14 shows the change in the size of the yield surface. It is evident from the figure that temperature has a distinct effect on the yield surface’s shape and size, which must be accounted for during thermoforming analysis.
56
60 40
0
y
σ (MPa)
20
-20 25°C 93°C 149°C 204°C 260°C
-40 -60 -60
-40
-20
0 20 σ (MPa)
40
60
x
Figure 4.59 2D-Plot of Barlat’s YLD96 yield function for AA3003-H111 at several elevated temperatures. Plot of actual stresses to show the change in the size of the yield surface.
Figure 4.15 shows the effect of αx, αy and αz0 on the yield surface shape (Barlat et al., 1997). The αz1 parameter primarily affects the yield surface for the stress state with σxy and therefore its effect cannot be seen in a 2D plot (Figure 4.15) where a constant shear stress is assumed.
57
1.5
α
x
1
α
x
σy/ ⎯σ
0.5
α
α
y
y
0
α
z0
-0.5
α
z0
-1 -1.5 -1.5
-1
-0.5
0
0.5
1
1.5
σ / ⎯σ x
Figure 4.60 Effect of αx, αy and αz0 on the yield function (Barlat et al., 1997).
Due to the significant effect that temperature has on material properties, the application of a temperature-dependent constitutive model for accurate analysis of warm forming becomes imperative. It is also necessary that this temperaturedependent constitutive model be used in a coupled thermo-mechanical finite element analysis of warm forming process where the thermal analysis provides temperature as input to the mechanical model (Čanađija et al., 2004). It is from such a coupled analysis that deformation stresses corresponding to both thermal and mechanical deformation could be accurately calculated. This issue is discussed in more detail in a later chapter. 58
At a first glance, it seems that the coefficients of Barlat YLD96 in Table 4.5 are fluctuating randomly without any apparent trend. This is attributed to the fact that the experimental data needed to calculate the anisotropic coefficients of YLD96 yield function (Table 3.1) do not linearly increase with temperature. Looking at the stress values in Figure 4.3 it is noticed that they decrease with temperature, while plastic anisotropy parameters (R0, R45 and R90) in Figure 4.11 increase with temperature. This opposite behavior causes the anisotropy coefficients of the yield function to fluctuate. Therefore, it is important to use higher order fitting functions to capture these variations. To show the importance of varying the anisotropy coefficients with temperature, the plastic anisotropy parameters R0, R45 and R90 were first calculated directly from the yield function using Equation (3.36) and the values of the yield function coefficients at the discrete temperatures listed in Table 4.5. Then, the plastic anisotropy parameters (R0, R45 and R90) were calculated at those same temperatures, but this time values of the yield function anisotropy coefficients at 25°C were used for all temperatures. In both cases, however, the temperature dependent hardening values from Table 4.4 were used, as flow stress in general only affects the radius of the yield function and not the plastic anisotropy parameters. Figure 4.16 shows the result of the two tests for the R0 case.
59
2.5
2.0
R0
1.5
1.0
0.5
Experimental-R0 Variable anisotropy coefficients Constant values @25°C
0.0 0
50
100
150
200
250
300
Temperature (°C)
Figure 4.61 Plot of experimental values of plastic anisotropy parameter (R0) of AA3003-H111 compared with predictions from the Barlat YLD96 material model. Calculations were done first using variable values of the yield function with each temperature and then by using a fixed value for all temperatures. From Figure 4.16, it could be seen that the plastic anisotropy parameter R0 closely followed the experimental data, when the values of anisotropy coefficients at each temperature from Table 4.5 were used. However, the resulting plastic anisotropy parameters remained constant when the yield function coefficients were assumed to remain constant for all temperatures. This indicates the importance of using temperature dependent anisotropy coefficients, when accurate modeling results are expected. In order to develop an anisotropic material model for use in a coupled thermomechanical finite element analysis of sheet metal forming processes, the anisotropy coefficients of a yield function must be represented as a function of
60
temperature. Thus, curve-fitting methods were used to fit the anisotropy coefficients shown in Table 4.5. Figures 4.17 through 4.23 show a plot of the YLD96 coefficients, c1, c2, c3, c6, αx, αy, αz0 and αz1 as a function of temperature. In these plots, two levels of polynomial curve fitting, 3rd and 5th order were used to represent variation with respect to temperature. The fit functions used for all seven anisotropy coefficients are shown in Table 4.6. It should be noted that no extrapolation to temperatures outside the experimental temperature range, i.e. 25°C – 260°C, using the supplied polynomial equations can be performed since no experimental data is available.
1.2
1.0
C1(T)
0.8
0.6
0.4
0.2 Experimental
3th Order Fit
5th Order Fit
0.0 0
50
100
150
200
250
300
Temperature (°C)
Figure 4.62 Plot of (c1) as a function of temperature for AA3003-H111. 3rd and 5th order polynomial curve fits are shown also.
61
1.2
1.0
C2(T)
0.8
0.6
0.4
0.2 Experimental
3th Order Fit
5th Order Fit
0.0 0
50
100
150
200
250
300
Temperature (°C)
Figure 4.63 Plot of (c2) as a function of temperature for AA3003-H111. 3rd and 5th order polynomial curve fits are shown also. 1.2
1.0
C3(T)
0.8
0.6
0.4
0.2 Experimental
3th Order Fit
5th Order Fit
0.0 0
50
100
150
200
250
300
Temperature (°C)
Figure 4.64 Plot of (c3) as a function of temperature for AA3003-H111. 3rd and 5th order polynomial curve fits are shown also. 62
1.2
1.0
C6(T)
0.8
0.6
0.4
0.2 Experimental
3th Order Fit
5th Order Fit
0.0 0
50
100
150
200
250
300
Temperature (°C)
Figure 4.65 Plot of (c6) as a function of temperature for AA3003-H111. 3rd and 5th order polynomial curve fits are shown also. 1.4
1.2
1.0
αx(T)
0.8
0.6
0.4
0.2 Experimental
3th Order Fit
5th Order Fit
0.0 0
50
100
150
200
250
300
Temperature (°C)
Figure 4.66 Plot of (αx) as a function of temperature for AA3003-H111. 3rd and 5th order polynomial curve fits are shown also. 63
1.8 1.6 1.4
αy(T)
1.2 1.0 0.8 0.6 0.4 0.2 Experimental
3th Order Fit
5th Order Fit
0.0 0
50
100
150 Temperature (°C)
200
250
300
Figure 4.67 Plot of (αy) as a function of temperature for AA3003-H111. 3rd and 5th order polynomial curve fits are shown also. 2.5
2.0
αz1(T)
1.5
1.0
0.5
Experimental
3th Order Fit
5th Order Fit
0.0 0
50
100
150 Temperature (°C)
200
250
300
Figure 4.68 Plot of (αz1) as a function of temperature for AA3003-H111. 3rd and 5th order polynomial curve fits are shown also. 64
Table 4.6 Barlat YLD96 material model anisotropy coefficients for AA3003-H111 as a function of temperature using two different order polynomial fit functions. Temperatures in °C. YLD96 5th Order Fit 3rd Order Fit Coefficient = 1.0446 + 0.0045834T - = 0.90902 + 0.0097045T 3.9615E-5T2 + 9.7409E-8*T3 6.2517E-5T2 - 3.3687E-7T3 + c1 3.7316E-9T4 - 7.6631E-12T5 = 0.9867 + 0.000373T + = 0.90308 + 0.0023977T + 1.1717E-7T2 + 3.0052E-9T3 3.5398E-5T2 - 8.9808E-7T3+ c2 5.4066E-9T4 - 9.8763E-12T5 =0.9860 + 0.0003334T - = 0.98895 + 0.0010519T 1.7857E-6T2 + 1.593E-9T3 3.7432E-5T2 + 4.7418E-7T3 c3 2.3519E-9T4 + 3.9325E-12T5 = 1.0592 + 0.000855T - = 0.76161 + 0.017737T 1.3728E-5T2 + 3.696E-8T3 0.0003099T2 + 2.2357E-6T3 c6 7.2652E-9T4 + 8.8146E-12T5 = 1.0146 - 0.010636T + = 1.8838 - 0.05496T + 8.3876E-5T2 - 2.0664E-7T3 0.00073166T2 - 3.8433E-6T3 + αx 7.5607E-9T4 - 3.0894E-12T5 =1.4253 - 0.012009T + = 2.683 - 0.079188T + 6.4874E-5T2 - 1.3432E-7T3 0.0011348T2 - 7.0955E-6T3 + αy 1.9271E-8T4 - 1.8287E-11T5 = 0.86396 + 0.007906T + = 2.7451 - 0.096914T + 3.0245E-6T2 - 6.2854E-8T3 0.0017925T2 - 1.2899E-5T3 + αz1 4.0712E-8T4 - 4.7071E-11T5 To determine the validity and accuracy of the 3rd and 5th order polynomial curve fit functions, the plastic anisotropy parameters R0, R45 and R90 were calculated directly from the yield function using Equation (3.36) and compared to experimentally measured values. All material parameters used in these calculations, i.e. hardening parameters and Barlat’s YLD96 function anisotropy coefficients, were set as temperature-dependent equations from Tables 4.4 and 4.6. The results at discrete points are shown in Figures 4.24 through 4.26.
65
2.5
2.0
R0
1.5
1.0
0.5
Experimental-R0
3rd-Order-Fit
5th-Order-Fit
0.0 0
50
100
150 Temperature (°C)
200
250
300
Figure 4.69 Plot of experimental values of plastic anisotropy parameter (R0) for AA3003-H111 compared with predictions from Barlat’s YLD96 material model where the yield function parameters are a function of temperature. Calculations were done using 3rd and 5th order curve fits for the yield function anisotropy coefficients.
66
3.0
2.5
R45
2.0
1.5
1.0
0.5 Experimental-R45
3rd-Order-Fit
5th-Order-Fit
0.0 0
50
100
150 Temperature (°C)
200
250
300
Figure 4.70 Plot of experimental values of plastic anisotropy parameter (R45) for AA3003-H111 compared with predictions from Barlat’s YLD96 material model where the yield function parameters are a function of temperature. Calculations were done using 3rd and 5th order curve fits for the yield function anisotropy coefficients.
67
1.4
1.2
1.0
R90
0.8
0.6
0.4
0.2
Experimental-R90
3rd-Order-Fit
5th-Order-Fit
0.0 0
50
100
150
200
250
300
Temperature (°C)
Figure 4.71 Plot of experimental values of plastic anisotropy parameter (R90) for AA3003-H111 compared with predictions of Barlat’s YLD96 material model where the yield function parameters are a function of temperature. Calculations were done using 3rd and 5th order curve fits for the yield function anisotropy coefficients. While both results of 3rd and 5th order polynomial curve fit functions of the plastic anisotropy parameters R0 and R45 closely agree with the experimental data, the 5th order fit function appears to match the data slightly better than the 3rd order function. On the other hand, in the case of plastic anisotropy in the transverse direction, R90, neither of the two functions agrees very well with experimental data at temperatures lower than 150 oC. This difference could be attributed to inaccuracies associated with both experimental tests as well as with the fitting of the yield function anisotropy coefficients.
68
The experimentally determined anisotropy coefficients were compared to those calculated from the 3rd and 5th order polynomial curve fit functions from Table 4.6 for a temperature of 25°C (77°F) in the plot of Barlat’s YLD96 yield function shown in Figure 4.27. There is very little difference between these two yield surface curves.
σ / ⎯σ x
Figure 4.72 Plot of Barlat’s YLD96 yield function for AA3003-H111 using anisotropy coefficients extracted from experimental data compared to anisotropy coefficients from both the 3rd and 5th order curve fits shown in Table 4.6. Temperature was 25°C (77°F).
69
Figure 4.28 shows the yield function plot for a temperature of 204°C (400°F). Again in both cases, the yield functions coincide with very small differences. (All yield function plots show normalized values with respect to σ ).
σ / ⎯σ x
Figure 4.73 Plot of Barlat’s YLD96 yield function for AA3003-H111 using anisotropy coefficients extracted from experimental data compared to anisotropy coefficients from both the 3rd and 5th order curve fits shown in Table 4.6. Temperature was 204°C (400°F).
70
4.5. Barlat’s YLD2000-2d Anisotropy Coefficients Calculation
Results from the uniaxial and bulge tests were implemented into Equations (3.37―3.43), and Barlat’s YLD2000-2d model independent coefficients αk (for k from 1 to 8) were calculated at several temperatures. Implementing these values into Equations (3.18) and (3.19), the coefficients of L′ and L′′ were calculated. The results are shown in Table 4.7 for L′ and Table 4.8 for L′′ . The nonquadratic parameter for aluminum alloys (FCC) is a = 8.
Table 4.7 Barlat’s YLD2000-2d model anisotropy L′ coefficients calculated at several elevated temperatures. Temperature L′ 12 L′ 21 L′ 66 L′ 11 L′ 22 °C 25
0.620326
-0.310163
-0.365965
0.731930
1.069533
38
0.644069
-0.322034
-0.361830
0.723660
1.103172
66
0.624512
-0.312256
-0.381120
0.762240
1.108886
93
0.621924
-0.310962
-0.390826
0.781651
1.128550
121
0.660493
-0.330246
-0.352849
0.705697
1.058723
149
0.643380
-0.321690
-0.370560
0.741120
1.096767
177
0.665820
-0.332910
-0.354979
0.709958
1.076084
204
0.660504
-0.330252
-0.366006
0.732013
1.076314
232
0.656426
-0.328213
-0.372631
0.745263
1.088619
260
0.662875
-0.331437
-0.374668
0.749336
1.081612
71
Table 4.8 Barlat’s YLD2000-2d model anisotropy L′′ coefficients calculated at several elevated temperatures. Temperature L′′ 12 L′′ 21 L′′ 22 L′′ 66 L′′ 11 °C 25
0.681653
-0.331290
-0.350798
0.699980
1.060969
38
0.656793
-0.311223
-0.319738
0.678789
1.053135
66
0.659869
-0.301271
-0.325601
0.679240
1.041682
93
0.652165
-0.300168
-0.318332
0.693067
1.007063
121
0.656470
-0.328693
-0.327657
0.680490
0.938710
149
0.653195
-0.314113
-0.321840
0.684933
0.970763
177
0.644020
-0.317444
-0.312396
0.674525
0.923349
204
0.638760
-0.310674
-0.310110
0.672910
0.918824
232
0.631735
-0.293619
-0.300007
0.659711
0.890074
260
0.617323
-0.284566
-0.285783
0.652661
0.847258
Figure 4.29 is a plot of the yield function for AA3003-H111, comparing Barlat’s YLD200-2d model at 25°C (77°F) with the isotropic von Mises yield function (normalized stress values with respect to σ ). Figure 4.30 is a similar plot comparing Barlat’s YLD2000-2d function at several elevated temperatures (normalized stress values with respect to σ ), indicating that temperature has an effect on the shape of the yield surface. This change in the shape is particularly noticed near the balanced biaxial region of the yield surface. This is to be expected, since any change in the plastic anisotropy parameters (R-values) would affect the shape of the yield surface.
72
σ / ⎯σ x
Figure 4.74 2D-Plot of von Mises yield function vs. Barlat’s YLD2000-2d yield function at 25°C (77°F) for AA3003-H111.
73
σ / ⎯σ x
Figure 4.75 2D-Plot of Barlat’s YLD000-2d yield function for AA3003-H111 at several elevated temperatures. Normalized stress values. In order to develop an anisotropic material model for use in a coupled thermomechanical finite element analysis of sheet metal forming processes, the anisotropy coefficients of a yield function must be represented as a function of temperature. Curve-fitting methods were used to fit the anisotropy coefficients shown in Tables 4.7 and 4.8. Figures 4.31 through 4.40 show a plot of the YLD2000-2d coefficients, L′ 11 - L′ 66 and L′′ 11 - L′′ 66, as a function of temperature. In these plots, a 4th level of polynomial curve fitting was used to
74
represent variation with respect to temperature. The fit functions used for all anisotropy coefficients are shown in Table 4.9 for L′ and Table 4.10 for L′′ using a 4th order polynomial curve fit.
0.70
0.60
L'11(T)
0.50
0.40
0.30
0.20
0.10 Experimental
4th Order Fit
0.00 0
50
100
150
200
250
300
Temperature (°C)
Figure 4.76 Plot of ( L′ 11) as a function of temperature for AA3003-H111 Experimental data is fitted using a 4th order polynomial curve.
75
0.00 0
50
100
150
200
250
300
-0.05 Experimental
-0.10
4th Order Fit
L'12(T)
-0.15 -0.20 -0.25 -0.30 -0.35 -0.40 Temperature (°C)
Figure 4.77 Plot of ( L′ 12) as a function of temperature for AA3003-H111 Experimental data is fitted using a 4th order polynomial curve. 0.00 0
50
100
150
200
250
300
-0.05 Experimental
4th Order Fit
-0.10
L'21(T)
-0.15 -0.20 -0.25 -0.30 -0.35 -0.40 -0.45 Temperature (°C)
Figure 4.78 Plot of ( L′ 21) as a function of temperature for AA3003-H111 Experimental data is fitted using a 4th order polynomial curve. 76
0.9 0.8 0.7
L'22(T)
0.6 0.5 0.4 0.3 0.2 0.1 Experimental
4th Order Fit
0.0 0
50
100
150
200
250
300
Temperature (°C)
Figure 4.79 Plot of ( L′ 22) as a function of temperature for AA3003-H111 Experimental data is fitted using a 4th order polynomial curve. 1.2
1.0
L'66(T)
0.8
0.6
0.4
0.2 Experimental
4th Order Fit
0.0 0
50
100
150
200
250
300
Temperature (°C)
Figure 4.80 Plot of ( L′ 66) as a function of temperature for AA3003-H111 Experimental data is fitted using a 4th order polynomial curve. 77
0.8 0.7 0.6
L''11(T)
0.5 0.4 0.3 0.2 0.1 Experimental
4th Order Fit
0.0 0
50
100
150
200
250
300
Temperature (°C)
Figure 4.81 Plot of ( L′′ 11) as a function of temperature for AA3003-H111 Experimental data is fitted using a 4th order polynomial curve. 0.00 0
50
100
150
200
250
300
-0.05 Experimental
4th Order Fit
-0.10
L''12(T)
-0.15 -0.20 -0.25 -0.30 -0.35 -0.40 Temperature (°C)
Figure 4.82 Plot of ( L′′ 12) as a function of temperature for AA3003-H111 Experimental data is fitted using a 4th order polynomial curve. 78
0.00 0
50
100
150
200
250
300
-0.05 Experimental
4th Order Fit
-0.10
L''21(T)
-0.15 -0.20 -0.25 -0.30 -0.35 -0.40 Temperature (°C)
Figure 4.83 Plot of ( L′′ 21) as a function of temperature for AA3003-H111 Experimental data is fitted using a 4th order polynomial curve. 0.8 0.7 0.6
L''22(T)
0.5 0.4 0.3 0.2 0.1 Experimental
4th Order Fit
0.0 0
50
100
150
200
250
300
Temperature (°C)
Figure 4.84 Plot of ( L′′ 22) as a function of temperature for AA3003-H111 Experimental data is fitted using a 4th order polynomial curve. 79
1.2
1.0
L''66(T)
0.8
0.6
0.4
0.2 Experimental
4th Order Fit
0.0 0
50
100
150
200
250
300
Temperature (°C)
Figure 4.85 Plot of ( L′′ 66) as a function of temperature for AA3003-H111 Experimental data is fitted using a 4th order polynomial curve.
Table 4.9 Barlat YLD2000-2d material model anisotropy coefficients for AA3003-H111 as a function of temperature using a 4th order polynomial fit function. Temperatures in °C. YLD2000-2d 4th Order Fit Coefficient L′ 11
= 0.64823 - 0.0010529*T + 1.5487e-05*T2 - 7.0602e-08*T3 + 1.0518e-10*T^4
L′ 12
= -0.32412 + 0.00052647*T - 7.7433e-06*T2 + 3.5301e-08*T3 5.2588e-11*T4
L′ 21
= -0.31409 - 0.0025156*T + 3.106e-05*T2 - 1.4139e-07*T3 + 2.1408e-10*T4
L′ 22
= 0.62818 + 0.0050312*T - 6.2119e-05*T2 + 2.8279e-07*T3 4.2815e-10*T4
L′ 66
= 0.97613 + 0.0053338*T - 6.7055e-05*T2 + 3.1471e-07*T3 4.9904e-10*T4
80
Table 4.10 Barlat YLD2000-2d material model anisotropy coefficients for AA3003-H111 as a function of temperature using a 4th order polynomial fit function. Temperatures in °C. YLD2000-2d 4th Order Fit Coefficient L′′ 11
= 0.71012 - 0.0018074*T + 2.0358e-05*T2 - 9.3684e-08*T3 + 1.4202e-10*T4
L′′ 12
= -0.39394 + 0.0035533*T - 4.4893e-05*T2 + 2.1215e-07*T3 3.3022e-10*T4
L′′ 21
= -0.3855 + 0.0022626*T - 2.6669e-05*T2 + 1.2449e-07*T3 1.916e-10*T4
L′′ 22
= 0.72272 - 0.0015865*T + 2.0637e-05*T2 - 1.0375e-07*T3 + 1.6886e-10*T4
L′′ 66
= 1.039 + 0.0017293*T - 4.0784e-05*T2 + 2.3049e-07*T3 4.2396e-10*T4
To determine the validity and accuracy of the curve fit function for YLD20002d coefficients, the plastic anisotropy parameters R0, R45 and R90 were calculated directly from the yield function using Equation (3.36) and compared to experimentally measured values. All material parameters used in these calculations, i.e. hardening parameters and Barlat’s YLD2000-2d function anisotropy coefficients, were set as temperature-dependent equations from Tables 4.4, 4.9 and 4.10, respectively. The results at discrete points are shown in Figures 4.41 through 4.43.
81
2.5
2.0
R0
1.5
1.0
0.5
Experimental-R0
4rd-Order-Fit
0.0 0
50
100
150 Temperature (°C)
200
250
300
Figure 4.86 Plot of experimental values of plastic anisotropy parameter (R0) for AA3003-H111 compared with predictions from Barlat’s YLD2000-2d material model where the yield function parameters are a function of temperature. Calculations were done using the 4th order curve fit for the yield function anisotropy coefficients from Tables 4.9 and 4.10.
82
3.0
2.5
R45
2.0
1.5
1.0
0.5 Experimental-R45
4rd-Order-Fit
0.0 0
50
100
150 Temperature (°C)
200
250
300
Figure 4.87 Plot of experimental values of plastic anisotropy parameter (R45) for AA3003-H111 compared with predictions from Barlat’s YLD2000-2d material model where the yield function parameters are a function of temperature. Calculations were done using the 4th order curve fit for the yield function anisotropy coefficients from Tables 4.9 and 4.10.
83
1.6 1.4 1.2
R90
1.0 0.8 0.6 0.4 0.2 Experimental-R90
4rd-Order-Fit
0.0 0
50
100
150 Temperature (°C)
200
250
300
Figure 4.88 Plot of experimental values of plastic anisotropy parameter (R90) for AA3003-H111 compared with predictions from Barlat’s YLD2000-2d material model where the yield function parameters are a function of temperature. Calculations were done using the 4th order curve fit for the yield function anisotropy coefficients from Tables 4.9 and 4.10. The experimentally determined anisotropy coefficients were compared to those calculated from the 4th order polynomial curve fit functions from Tables 4.9 and 4.10 for a temperature of 25°C (77°F) in the plot of Barlat’s YLD2000-2d yield function shown in Figure 4.44. There is very little difference between these two yield surface curves.
84
σ / ⎯σ x
Figure 4.89 Plot of Barlat’s YLD2000-2d yield function for AA3003-H111 using anisotropy coefficients extracted from experimental data compared to anisotropy coefficients from both the 4th order curve fit shown in Tables 4.9 and 4.10. Temperature was 25°C (77°F).
Figure 4.45 shows the yield function plot for a temperature of 260°C (500°F). Again in both cases, the yield functions coincide with very small differences. (All yield function plots show normalized values with respect to σ ).
85
σ / ⎯σ x
Figure 4.90 Plot of Barlat’s YLD2000-2d yield function for AA3003-H111 using anisotropy coefficients extracted from experimental data compared to anisotropy coefficients from both the 4th order curve fit shown in Tables 4.9 and 4.10. Temperature was 260°C (500°F).
86
CHAPTER 5
STAMPING EXPERIMENTS
In this section, the experimental apparatus for sheet stamping will be described, followed by a discussion of the warm forming experimental results obtained using the 101.6mm (4 in) diameter hemispherical punch.
5.1. Experimental Apparatus
The apparatus used in these experiments was built around an Interlaken 75 double action servo press, Figure 5.1, manufactured by Interlaken Technology Corporation, Eden Prairie, Minnesota. The double action refers to the clamping mechanism moving independently of the punch mechanism. This allows for the boundaries of the sheet blank to be clamped while the punch pushes the sheet into the die cavity filled with supporting fluid. The ability to independently control both the clamp and the punch affords the opportunity for various modifications of the experimental procedure. The apparatus uses the LDH (Limiting Dome Height) setup used in industry for the evaluation of lubricants in the sheet metal. Some modifications were made to change the setup to research requirements. The LDH die is essentially a pair of cylinders that are clamped together after placing a draw blank between them. The punch moves through one chamber, meets the material and stretches it into the second cylinder. The clamping
87
mechanism typically contains a draw bead that would cause the sheet metal to get stretched only.
Figure 5.11 The modified Interlaken servo press 75 used for the forming of aluminum sheets. Sheet metal stamping uses a male (punch), a female and a blank holding die, to plastically (permanently) deform a blank sheet of metal into a desired shape. This technique is widely used in order to fabricate thousands of sheet metal structures per day in several industries, e.g. automotive, aerospace, beverage industry etc… Depending upon several factors such as geometry, volume and material type, deep drawing or stretch forming is used as different methods to form sheet metals (Cao et al., 1997). In sheet-forming processes however, several types of failures could occur, such as rupturing, necking, wrinkling and springback (Kawaka et al., 2001), that are undesirable. Also, there are significant expenses associated with the necessary tooling and the success of the process largely depends on the ingenuity of skilled machinists to bring the economic costs down to suitable levels. Sheet metal stamp hydroforming is currently being 88
considered as a viable alternative to conventional sheet metal stamping. Stamp hydroforming offers many advantages over conventional stamping when fabricating difficult-to-form parts. These advantages include improved formability due to the applied fluid pressure, low wear rate of tools, better distribution of plastic deformation, significant economic savings due to using less number of tooling, and finally the potential for reduced amount of finishing work required (Youssef et al., 1998). Warm forming of aluminum parts is a method for improving formability of aluminum sheets compared to cold (room temperature) forming. Increasing the temperature of the sheet during the forming process reduces the tensile strength of the material and increases the elongation. This would cause deeper draws in the material and the formed products would not differ significantly from parts formed at room temperature. The sheet metal stamp warm-forming method, shown in Figure 5.2, is a process in which a part is formed by a punch (in this case a hemispherical punch). The workpiece is placed on the clamping mechanism, as shown in Figure 5.2(1). Figure 5.2(2) shows the upper chamber being lowered and the workpiece clamped securely between the two die halves. Band heaters positioned on the outer body of the upper, middle and lower dies heat the chamber and the fluid inside the chamber is raised to the desired temperature. A set of thermocouples connected to a controller measures and maintains the temperature level to achieve the isothermal condition.
89
Figure 5.12 Schematics of a warm-forming press using a hemispherical punch with a blank holding support. After achieving a constant temperature level, the punch moves up, and the sheet metal begins to deform and take on the shape of the punch, in this case a hemispherical shape, as shown in Figure 5.2(3). Once the punch reaches the prescribed depth, the upper chamber is raised, as shown in Figure 5.2(4), to remove the formed part from the die. The experimental setup was used to form 101.6 mm (4 in) diameter hemispherical cups from 1 mm (0.04 in) thick, 177.8 mm (7 in) diameter round blanks of AA3003-H111under pure stretch condition. The blank was placed over a draw bead and clamped with a blank holding force (BHF) of approximately 267
90
kN (60000 lbf) to insure no material will draw in during the pure stretch experiments. Heating elements with an active control device were added to the LDH machine in order to raise the temperature of the dies and the sheet to the desired elevated temperature. The active control was achieved by using two thermocouples linked to the die and blank system. Other thermocouples were installed to measure the temperature directly of both the blank and the punch during the forming process. These are critical measurements needed in order to be able to perform accurate numerical analysis of the experiment. The new setup is shown in Figure 5.3.
Figure 5.13 Experimental stamping setup showing heating elements and controller.
91
5.2. Experimental Work
The procedure for performing the experiment at a specific elevated temperature is as follows. The blank was clamped between the dies with three round heating elements positioned around the perimeter of the dies. The system was then insolated to minimize heat loss to the environment. Using the temperature control device, the desired temperature was set and maintained for a period of about 20 min or until a constant and isothermal condition was achieved. Temperature was monitored using several thermocouples within the system. Initially, the temperature of the punch was not controlled independently. A mechanism to control the temperature of the punch is currently under development. For the current research, the punch temperature was found to be cooler than the blank, see Table 8.1. With an isothermal blank condition, the punch was then actuated to stretch the blank. During the process, the punch force-displacement curve was recorded. This process was repeated several times for each temperature to establish repeatability Pure stretch experiments were performed at several elevated temperatures in the range of 25°C – 204°C (77°F – 400°F). Figures 5.4 to 5.8 show the experimental results at different temperatures with failure location and punch depth at failure shown.
92
Failure Depth: 27mm (1.1in)
Figure 5.14 Pure stretch experimental results using the 101.6mm (4in) hemispherical punch at 25°C (77°F).
Failure Depth: 31mm (1.2in)
Figure 5.15 Pure stretch experimental results using the 101.6mm (4in) hemispherical punch at 93°C (200°F).
93
Failure Depth: 34mm (1.3in)
Figure 5.16 Pure stretch experimental results using the 101.6mm (4in) hemispherical punch at 149°C (300°F).
Failure Depth: 36mm (1.4in)
Figure 5.17 Pure stretch experimental results using the 101.6mm (4in) hemispherical punch at 177°C (350°F).
94
Failure Depth: 37mm (1.5in)
Figure 5.18 Pure stretch experimental results using the 101.6mm (4in) hemispherical punch at 204°C (400°F).
Figure 5.9 shows the pure stretch results at these temperatures side by side for comparison reasons. In the figure, punch depth at failure and failure locations are shown. Figure 5.10 shows the punch load vs. punch depth at several elevated temperatures.
95
25°C, 27 mm
93°C, 31 mm
149°C, 34 mm
Failure 177°C, 36 mm
204°C, 37 mm
Figure 5.19 Pure stretch experimental results using the 101.6mm (4in) hemispherical punch at several elevated temperatures as indicated along with failure punch depth. The failure locations are also shown.
96
18 38 °C
16
93 °C
149 °C 177 °C
25 °C 14 204 °C
Punch Load (kN)
12 10 8 6 4 2 0 0
5
10
15
20 25 Punch Depth ( mm)
30
35
40
Figure 5.20 Punch force vs. punch depth for the hemispherical punch experiments at several elevated temperatures.
From these figures (Figures 5.4 and 5.10), it is noticed that as the test temperature increased the following is noticed: i.
Forming depth increased from 27 mm (1.06 in) at room temperature to 37 mm (1.46 in) at 204°C (400°F). This represents a 37% increase in forming depth.
ii.
The location of failure (where material ruptures) changed dramatically during the forming process (Figure 5.9). At room temperature, failure occurred at the contact edge with the punch, but as temperature increased the failure point migrated from the punch contact region
97
toward the blank holder contact edge. At temperatures of 177°C (350°F) and 204°C (400°F), failure occurred at the die contact edge. iii. The punch load required to form the part reduced as temperature increased. This was expected since the tensile strength of the material reduced with temperature as shown in Chapter 4 Figure 4.5.
The greater forming depths achieved at elevated temperatures are attributed to the temperature gradient between the blank and the punch. In the experimental setup, the punch itself was not heated directly and was found to be at a lower temperature than the blank (see Table 8.1). This condition was advantageous because areas of the blank in contact with the punch were lower in temperature and thus stronger than the unsupported areas of the blank. As the punch traveled, the unsupported regions of the blank (regions not in contact with the punch) stretched more due to their lower tensile strength. At the highest possible temperature for this setup, 204°C (400°F), the unsupported regions continued to stretch until critical thinning occurred at the edge of the sheet and blank holder.
98
CHAPTER 6
NUMERICAL ANALYSIS SETUP
An important goal in manufacturing research is determining the optimum method of production of efficient products with less cost. The optimization criterion varies, depending on the products, but having a thorough understanding of the manufacturing processes is an essential step. Sheet metal forming design requires the understanding of the fundamentals of deformation mechanics involved in the processes. Without proper understanding of the effect of different variables such as material properties, friction and geometry the design process would be difficult, time consuming, and expensive. Also it would not be possible to predict and prevent defects form occurring until it is too late. Failure modes such as necking, wrinkling and springback may occur in sheet metal forming process. The automotive industry in recent years has seen more use of very thin high strength materials in which defects like folding and wrinkling occur more often. The finite element method (FEM) gives an advantage in predicting such defects, before the real stamping operation takes place (Kawaka et al. 2001). For the stamp warm forming process the numerical study was performed using the explicit finite element code, LS-Dyna 3D. In the following a general description of the FEA code built for this analysis is discussed.
99
6.1. Preliminary Results
In the early stages of the research, a structural only finite element analysis (FEA) was performed of the warm forming process, where all parts were assumed to be in an isothermal condition during the forming process, i.e. every part in the process had the same temperature and there was no temperature gradient. Preliminary results from this model did not compare well with the experiments. That is, both failure location and punch load vs. punch depth curve were different from the experiments. Figure 6.1 shows the predicted punch load vs. punch depth curve as compared with the experimental curve at the temperature of 204°C (400°F). In the numerical analysis the blank was assumed to have an isothermal condition, i.e. the temperature of the blank did not change as it contacted the dies and the punch. In the structural FEM analysis, the anisotropy coefficients for the material model were set to be constant and equal to those at 204°C (400°F)
100
16 Experimental Result @ 204°C
14
Numerical Result @ 204°C
Punch Load (kN)
12 10 8 6 4 2 0 0
5
10
15
20
25
30
35
40
Punch Depth (mm)
Figure 6.9 Initial results of punch force vs. punch depth for the hemispherical punch compared to numerical predictions at 204°C (400°F) using only structural FEM analysis. Upon careful observation of the experimental process and using multiple thermocouples to record the temperature of the blank, punch, and dies, it was noted that the punch is at a lower temperature than the blank and the dies. This was due to the fact that the punch itself was not directly heated; rather its temperature was indirectly raised through heat transfer. Therefore, as the punch moved into the die cavity and contacted the blank, parts of the blank contacting the punch lost some of its heat to the punch. Figure 6.2 shows a schematic of the effect of having the band heaters placed on the outside and temperature drop toward the center of the die cavity.
101
Band Heaters
Hotter
Colder
Hotter
Figure 6.10 A schematic showing the effect of temperature drop away from the band heaters. The dies which are contacting the band heaters are at a higher temperature than the punch (which is not heated directly), and the center of the sheet. Next, to better understand the effect of temperature gradient on the deformation, the blank in the numerical analysis was divided into four sections, with each section having constant material properties at different temperatures. The center of the blank was assumed to be at the lowest temperature, while the edge of the blank was assumed to be at the maximum temperature. This model crudely approximated the measured temperature gradient in the sheet in our experimental setup when band heaters were placed on the outside of the die. Figure 6.3 shows the setup used for the blank and the temperatures assigned to each section.
102
204°C 177°C 163°C 149°C
Figure 6.11 A blank divided into multiple sections with each section having different temperature. Temperatures used in the analysis are shown on the graph. A structural only analysis was then performed using constant material properties appropriate for that section’s temperature. Figure 6.4 shows a comparison of the predicted punch load vs. punch depth curve with the experimental curve at 204°C (400°F). As could be seen from the graph, there is a significant improvement, as compared to Figure 6.1, in the match between the experimental and numerical curves. Although the punch force-punch depth curve nicely matched the experimental results, the model’s prediction of the punch depth at the time of failure did not compare well with the experiment.
103
16 Experimental Result @ 204°C
14
Numerical Result @ 204°C Using Sections
Punch Load (kN)
12 10 8 6 4 2 0 0
5
10
15
20
25
30
35
40
Punch Depth (mm)
Figure 6.12 Punch force vs. punch depth for the hemispherical punch compared to numerical predictions for the sectioned blank in Figure 6.3. To obtain accurate predictions of both force-displacement as well as failure of the sheet a fully coupled thermo-mechanical finite element analysis of warm forming process is required.
104
6.2. Coupled Thermal Structural Finite Element Model
Finite element analysis (FEA) was performed using the commercial explicit finite element code LS-Dyna (Hallquist, 1999) to understand the deformation behavior of the aluminum sheet during the thermoforming process. The UMAT option was used to build the user material subroutine in FORTRAN (COMPAQ VISUAL FORTRAN PROFESSIONAL EDITION 6.6B®), which was then linked with the library files supplied by LSTC. The finite element model used in the simulations was first created using Unigraphics® and imported as IGES (Initial Graphics Exchange Specification) files. Hypermesh® was used to create the finite element mesh, assign the boundary conditions and to build the LS-Dyna input deck for the analysis. The full size finite element model, Figure 6.5, used approximately 48000 four- and three-node shell elements. The punch, die, and the blank-holder were created using rigid materials (Material 20 in LS-Dyna). First trials with the adaptivity option in LS-Dyna to reduce calculation time revealed errors and problems in convergence of the thermal analysis. After contacting LSTC about this problem, they confirmed that this is indeed a bug in the LS-Dyna code that occurs in the adaptive algorithm when used in a coupled thermal mechanical analysis, the bug will be solved in the next version of LSDyna. Therefore, for the current analysis, the blank was modeled with a fine mesh of about 30,500 four-node shell elements without using adaptive meshing schemes.
105
Figure 6.13 LS-Dyna Full 3-D model created for stamp warm forming process with a hemispherical punch, using square blank. The thermal analysis was performed first, during which the temperature of each element was calculated and supplied as input to the UMAT. Using the temperature value for each element, the temperature dependent anisotropic material model coefficients were calculated. Before every structural iteration step, two thermal analysis steps were performed with a controlled time step to insure that the temperature update was adequate In this research, a linear fully implicit transient thermal analysis was performed with the diagonal scaled conjugate gradient iterative solver type in LSDyna. The die and blank materials were assumed to behave with isotropic
106
thermal properties. Table 6.1 shows the thermal properties defined in the analysis for the die and blank.
Table 6.1 Thermal properties of material used in numerical analysis.
Material
Density kg/mm3
Heat Capacity kN.mm/kg.K
Thermal conductivity kN/ms.K
Rigid Dies (FE)
7.85 E-6
450.0
7.00 E-5
Blank (AL)
2.71 E-6
904.0
2.22 E-4
The lower die, blank holder and punch were assigned a constant temperature boundary condition, while the blank was given an initial temperature boundary condition equal to the upper and lower dies. The temperature of the punch was set to the lower temperature based on experimental data (See Chapter 8). Thermal properties were assigned to contact surfaces to enable heat transfer at appropriate areas of contact between the blank and tooling during the analysis. Subsequently, areas of the blank that made contact with the punch lost heat to the punch while the unsupported regions of the blank remained at higher temperatures In the FEA simulations, the punch was given a trapezoidal velocity profile to fit the curve shown in Figure 6.6.
107
Punch Velocity 2.5
Velocity (mm/ms)
2.0
1.5
1.0
0.5
0.0 0
2
4
6
8
10
12
14
16
18
20
22
24
Time (ms)
Figure 6.14 Punch velocity used for the finite element simulation.
108
26
28
30
6.3. Failure Criteria
In sheet-forming operations, the deformation is characterized by biaxial stretching (Hosford et al., 1983). Failure in stretching operations normally occurs by the development of a sharp localized neck on the surface. By measuring minor and major strains in a specimen during deformation and plotting them, a Forming Limit Diagram (FLD) can be constructed. In sheet forming, the value of the measured strain near the necked region of the sheet is considered as ”failed strain”, while strain away from the necked region is considered as “safe strain”. Failure criteria used in the analysis are based on forming limit diagrams (FLD). FLD’s for AA3003-H111 at multiple temperatures were calculated with the M-K model (Marciniak et al, 1967) using Barlat’s YLD96 anisotropic yield function and appropriate coefficients for each temperature. Yoa et al. (2002) describe methods for extracting FLD for prediction of forming limit curves using an anisotropic yield function. In the current process, it was assumed that the loading path is sufficiently close to being linear that the use of a strain-based FLD to assess accurately the failure of the sheet is justified. For a general forming process in which the loading path may not be linear, it would be necessary to either integrate the M-K model into the FEM analysis to assess each element separately according to its loading path, or to use a stress-based FLD, which is less sensitive to strain path (Stoughton 2001). A review of different types of FLD and their use in FEA can be found in Stoughton et al. (2004).
109
The FLD curves for the current material were calculated using two types of hardening laws: Voce hardening law and Hollomon’s power law, as represented in Equation (6.1) and Equation (6.2) ⎛ ε ⎞ σ (ε ,ε ) = K (ε p + ε 0 )n ⎜ ⎟ ⎝ ε sr 0 ⎠ σ (ε p ) = A − B exp( −Cε p )
m
(6.1)
(6.2) where K (strength hardening coefficient), n (strain-hardening exponent) and m (strain-rate sensitivity index), A , B and C are material constants. ε p is the effective plastic strain and ε is the strain rate. ε0 is a constant representing the elastic strain to yield and ε sr 0 is a base strain rate (a constant). It should be noted that strain rate sensitivity (m) when calculating the FLD using Voce hardening law was incorporated into the algorithm using a multiplicative method. This insures that the Voce code has the strain rate sensitivity at elevated temperatures. It was found that there is a noticeable difference between the predictions of the two models. As seen in Figure 6.7, Voce hardening law predicts a lower failure limit curve as compared to the power law. A recent paper by Aghaie et al. (2004) also reports a similar difference between the two models and also shows that a FLD based on the Voce hardening law better predicts the experimental data. In general, the forming limits predicted by the Voce hardening law offer a more conservative prediction of the failure as compared to the power law.
110
0.45 0.40 0.35
Major Strain
0.30 0.25 0.20 0.15 0.10 0.05 0.00 -0.30
Power Law Voce -0.20
-0.10
0.00
0.10
0.20
0.30
0.40
Minor Strain
Figure 6.15 Forming limit diagrams for AA3003-H111 at 25°C (77°F). Calculations based on the M-K model and Barlat’s YLD96 anisotropic material model using two types of hardening law: Voce and Power law. In this paper, FLD curves based on the Voce hardening law (with strain rate sensitivity included) were used. Figure 6.8 shows the FLD curves for AA3003H111 at different temperatures. As seen from this figure, the forming limit curves increase with temperature, suggesting that AA3003-H111 can be stretched to higher levels before failure occurs.
111
0.90 0.80 0.70
Major Strain
0.60 0.50 0.40
25°C 38°C 93°C 149°C 177°C 204°C 232°C 260°C
0.30 0.20 0.10 0.00 -0.60
-0.40
-0.20
0.00
0.20
0.40
0.60
Minor Strain
Figure 6.16 Forming limit diagrams (FLD’s) for AA3003-H111 based on the M-K model, Barlat’s YLD96 anisotropic yield function, and Voce hardening law at several elevated temperatures.
112
CHAPTER 7
STRESS INTEGRATION FOR ELASTO-PLASTICITY USING ANISOTROPIC YIELD FUNCTIONS
Generally, heat generation due to dissipated mechanical work during plastic deformation leads to temperature rise in the specimen (Wriggers et al., 1992; Armero et al., 1993). This temperature rise is a local phenomenon however and dependent on the forming speed. During the characterization process of the material, the anisotropy coefficients were evaluated at the same temperature and the same forming speed, and therefore any thermal strain effects present were included in the material model. In the current thermoforming analysis, since the forming dies were maintained at a constant elevated temperature and the speed at which the test was performed was slow, the magnitude of thermal strains turned out to be of the same order as elastic strains and negligible compared to plastic strains. Therefore, to save on computational time during thermoforming analysis, the effect of thermal strain in the integration of the elastoplastic constitutive model was neglected. This is particularly important when performing thermoforming analysis of large industrial components. Figure 7.1 shows the plot of thermal, elastic and plastic strain increments calculated for the case of thermoforming at the elevated temperature of 204°C (400°F). It could be seen that thermal strains are between 2-3 orders of magnitude smaller than plastic
113
strains and the assumption to neglect them in the stress integration algorithm was well justified. Time (ms) 0 0.00
5.00
10.00
15.00
20.00
25.00
-1
Strain Values (Log Scale)
-2 -3
Elastic Strain Increment Plastic Strain Increment Thermal Strain Increment
-4 -5 -6 -7 -8 -9
Figure 7.5 Computed values of elastic, plastic and thermal strain increments for an element at different times. The results shown correspond to thermoforming analysis at the elevated temperature of 204°C (400°F). Stress integration of the elasto-plastic yield functions is explored in numerous publications (Armero and Simo, 1993; Auricchio and Taylor, 1999; Tuğcu and Neale, 1999; Hashiguchi, 2005). The temperature dependent YLD96 and YLD2000-2d model developed in the previous section was implemented within the framework of rate-independent plasticity and plane stress conditions using an efficient integration algorithm proposed by Simo et al. (1985) and further analyzed by Ortiz and Simo (1986). The current analysis is based on incremental theory of plasticity (Chung et al., 1993; Yoon et al., 1999; Han et al., 2003).
114
These algorithms, which fall within the class of cutting-plane methods of constrained optimization, was proposed to bypass the need for computing the gradients of the yield function and the flow rule as required by the closest point projection iterative methods (Simo and Hughes, 1998). The general closest point projection procedure usually leads to systems of nonlinear equations, the solution of which by the Newton-Raphson method requires evaluation of the gradients of system equations. While the previous approach might be applicable to simple plasticity models (e.g. von Mises), its application to complex yield functions such as Barlat’s YLD96 YLD2000-2d is not only exceedingly laborious, but also computationally extensive and makes the FEM code run slower for industrial applications. In the following, the stress integration algorithm for the YLD96 YLD2000-2d yield functions will be presented and the implementation of it as a UMAT into the explicit finite element code LS-Dyna will be described. First, the yield functions represented by Equation (3.1) for YLD96 and Equation (3.8) for YLD2000-2d can be written as:
Φ(σ , ε p , ε,T ) = σ (σ ,T ) − H (ε p , ε,T ) = 0
(7.1)
where H is the hardening rule defined by: ⎛ ε ⎞ H (ε p , ε,T ) = K (T )(ε p + ε 0 )n(T ) ⎜ ⎟ ⎝ ε sr 0 ⎠
m(T )
(7.2)
where T is the temperature calculated during the thermal analysis step and supplied as input to UMAT during the structural analysis part. In a displacement finite element formulation, the nature of the FEM code is strain driven. The cutting plane algorithm falls within the operator splitting
115
methodology (Ortiz, 1981) in which the strain is decomposed into two parts: elastic and plastic. The method proposed by Simo et al. (1985) and Ortiz and Simo (1986), however, is based on the total deformation theory. In this method, the history of total strain and total plastic strain are saved as history variables for the next step. This adds an unnecessary step, and in some cases where loading and unloading occurs, it might affect the accuracy of the code. Using the incremental theory of plasticity eliminates this step. The incremental theory of plasticity (Chung et al., 1993; and Yoon et al., 1999) was applied to the elasto-plastic formulation based on the materially embedded coordinate system. Under this scheme, the strain increments in the flow formulation are the discrete true (or logarithmic) strain increments, and the material rotates by the incremental angle obtained from the polar decomposition at each discrete step (Yoon et al., 2004). In addition, a multiplicative decomposition theory can be also utilized, especially when material deformation follows minimum plastic work path (or logarithmic strain path). Multiplicative theory formulation coincides with the current additive decomposition theory based on the incremental theory (Han et al., 2003). In the general commercial codes, e.g. LS-Dyna and Abaqus, the strain increment ( εn +1 ), the previous stress state value ( σ n ) and any history variables saved at the previous stress update step are provided at the beginning of each time step. The new strain increment is then assumed to be elastic, and an elastic predictor stress state “trial stress” is calculated through the customary elasticity
116
relations. Using the cutting plane algorithm, the actual stress state is then restored (plastic corrector) and other plastic variables are calculated. The basic steps in the numerical procedure for iteratively integrating the elastoplastic constitutive equations for rate independent plasticity with associated flow rule are:
εne+1 = εn +1
(7.3)
σ = C : εne+1
(7.4)
Associative flow rule:
ε p = λ
∂Φ ∂σ
(7.5)
Yield function:
Φ≤0
(7.6)
Normality parameter:
λ ≥ 0
(7.7)
Kuhn-Tucker condition:
λΦ = 0
(7.8)
Consistency condition:
=0 λΦ
(7.9)
where σ , εe and ε p are the stress, elastic strain and plastic strain, respectively.
C is the fourth order elastic tensor which is assumed to be a constant. The associated flow rule is expressed by Equation (7.5) in which λ is the plastic multiplier and Φ is the yield function as defined by Equation (7.1). The yielding criterion check and the loading-unloading conditions are expressed in the standard Kuhn-Tucker form (Simo, 1998) in which the constraints in Equations (7.6)-(7.9) are satisfied. The exact implementation of the iterative procedure for integrating the elastoplastic constitutive equations for YLD96 and YLD2000-2d is explained
117
next. In the finite element analysis using isoparametric elements, stress is updated at the Gauss (or Lobatto) points and the incremental deformation is given (Ortiz, 1986). The function of the user material subroutine (UMAT) is therefore to update the state variables (stress, plastic strain) from a converged configuration Bn to their corresponding values on the updated configuration Bn+1. To obtain the associated plastic strain increment the normality rule is utilized. From the associative flow rule:
ε p = λ
∂σ (σ ) ∂σ
(7.10)
The numerical procedure in updating the stress state involves finding the unknown λ (normality parameter). Using λ , all kinematics and stresses are updated at the end of the iteration procedure. In incremental theory, it should be noted that Δε p = λ as follows:
∂σ (σ ) σ : λ (σ ) λσ ∂σ Δε p = = = = λ σ (σ ) σ (σ ) σ (σ )
σ : ε p
(7.11)
where σ (σ ) is a first order homogenous function, i.e. σ (σ ) = σ
∂σ (σ ) p , and Δε ∂σ
is the equivalent plastic strain increment. In order to calculate Δε p , the calculation of σ (σ ) and
∂σ (σ ) is required. The explicit forms for these terms for ∂σ
both YLD96 and YLD2000-2d will be presented in the following section. At the beginning of each time step t, the strain increment ( εn +1 ), the previous total stress state value ( σ n ) and any history variables saved at the previous
118
stress update step are provided by the FEM code as input. In the elastic predictor step, the strain increment is assumed to be elastic and combined with the previous converged values of the state variables, a trial elastic stress state is calculated. If the new stress state lies outside the yield surface, this trial state is taken as the initial condition for the plastic corrector part. First, a trial stress state is calculated using the elasticity tensor ) σ n(trial = σ n + C : εn +1 +1
(7.12)
Using this trial stress, the yield function and its first derivative are calculated (see Section 7.2 and Section 7.3 for explicit derivation for YLD96 and YLD20002d, respectively) ) (trial ) σ n(trial +1 = σ (σ n +1 )
) ∂σ n(trial +1
∂σ
=
(7.13)
) ∂σ n(trial +1
∂σ
(σ (trial ) ) n +1
The size of the yield locus H (ε np+1, ε,T ) is calculated using the hardening flow rule as presented in Equation (7.2). Next, a check is performed to test whether the calculated trial stress state lies outside the yield surface, i.e. plastic state, by rewriting Equation (7.1) as ) p (trial ) p Φ(σ n(trial +1 , ε n +1, ε ,T ) = σ (σ n +1 ) − H (ε n +1, ε ,T ) ≤ 0
(7.14)
If the condition is met, then the trial stress state is elastic and the trial stress is the actual stress state that is returned to the FEM code. If the condition is not met, then the material has yielded and the stress state is plastic. The Newton-
119
Raphson method is then used to iteratively return the trial stress state to the yield surface by calculating the normality parameter λ using sub-steps i. After calculating the normality parameter λ , the sequential update procedure for the next stress state is done as follows:
σ n( i++11) = C : ⎡ε n +1 − ε np+1⎤
(7.15)
ε np+1 = ε np + Δε p
(7.16)
⎢⎣
⎥⎦
and since,
then
σ n( i++11) = C : ⎡ε n +1 − ε np ⎤ − C : Δε p
⎣⎢
⎦⎥
(7.17)
al ) , this becomes Combined with Equation (7.10) and setting σ n( i+) 1 = σ n(tri + 1
∂σ n( i+) 1(σ ) ( i +1) (i ) σ n +1 = σ n +1 − λC : ∂σ
(7.18)
With this new stress state, the yield function and the hardening rule are calculated at this new step and the yielding check is preformed again Φ( i +1) (σ n( i++11), ε np+( i1+1), ε,T ) = σ (σ n( i++11) ) − H (ε np+( i1+1), ε,T ) ≤ 0
(7.19)
The iterative procedure is repeated for a number of sub-steps until the plastic consistency is restored within a prescribed tolerance, i.e. Φ( i +1) (σ n( i++11), ε np+(i1+1), ε,T ) ≤ δ
where δ is a small number, e.g. δ = 10−7 .
120
(7.20)
The tangent modulus involved in the Newton-Raphson iteration to solve nonlinear equilibrium equations was obtained by consistently linearlizing the response function obtained from the integration algorithm (Simo et al., 1985; Ortiz et al., 1986). To solve for the normality parameter λ , the system equation as shown in Equation (7.14) is linearized using Taylor’s expansions:
)
(
(
∂Φ(i ) (i +1) ∂Φ(i ) p(i +1) p( i ) (i ) (i ) i ( ) 0 = Φ (σ n +1, ε n +1 ) + σ n +1 − σ n +1 + ε n +1 − ε np+( i1) p ∂σ ∂ε
)
(7.21)
From Equation (7.18)
∂σ n(i+) 1(σ ) ( i +1) (i ) σ n +1 − σ n +1 = −λC : ∂σ
(7.22)
ε np+(i1+1) − ε np(i ) = λ
(7.23)
and
Therefore, Equation (7.21) becomes (i ) ⎛ ⎞ ∂Φ(i ) ⎜ ∂σ n +1(σ ) ⎟ ∂Φ(i ) (i ) p( i ) ( i ) +λ 0 = Φ (σ n +1, ε n +1 ) + : −λC : ⎟ ∂σ ⎜ ∂σ ∂ε p
⎝
(7.24)
⎠
The normality parameter λ is then found from
λ =
Φ(i ) (σ n(i+) 1, ε np+(i1) )
∂σ n(i+) 1(σ ) ∂σ n(i+) 1(σ ) ∂H (i ) :C : − ∂σ ∂σ ∂ε p
(7.25)
In finite element commercial codes, e.g. LS-Dyna and Abaqus, for the plane stress case, the thickness strain needs to be updated and reported to the FEM code. This is done with the following equation:
121
p p p p εxx + εyy + ν ( εxx + εyy ) − 2ν ( εxx + εyy ) ) ( εzz =
ν −1
(7.26)
where the plastic strain is calculated using Equation (7.10) as follows:
εnp+1 = λ
∂σ (σ ) ∂σ
(7.27)
The cutting plane algorithm used to update the stress state according to previous equations is summarized in Table 7.1.
122
Table 7.1 Stress update algorithm based on incremental theory of plasticity. (i) Geometric update: (given by FEM code and user history variables)
εn +1, σ n , ε np
(ii) Elastic predictor:
σ n(0) +1 = σ n + C : ε n +1
p ε np+(0) 1 = εn (0) σ n(0) +1 = σ (σ n +1,T ) H (0) = σ (ε np+(0) 1 , ε ,T ) Φ(0) = σ n(0) − H (0) n +1 +1
(iii) Check for yielding: Φ(0) n +1 ≥ 0
NO: Material is elastic. Set trial state to be final state:
ε np+1 = ε np+(0) 1
σ n +1 = σ n(0) +1 ε np+1 = ε np+(0) 1 YES: Material is plastic. Set i=0 (iv) Plastic corrector:
Φ( i )
λ =
∂σ (i ) (σ ) ∂σ (i ) (σ ) ∂H (i ) :C : − ∂σ ∂σ ∂ε p ⎡ ∂σ ( i ) (σ ) ⎤ σ n(i++11) = σ n(i+) 1 − λ ⎢C : ⎥ ∂σ ⎥ ⎢⎣ ⎦ ∂σ (i +1) (σ ) εnp+(i1+1) = λ ∂σ ε np+(i1+1) = ε np+(i1) + λ
σ n(i++11) = σ (σ n( i++11) )
H (i +1) = H (ε np+(i1+1), ε,T ) ~ Continued
123
Table 7.1 (~ Continued) Stress update algorithm based on incremental theory of plasticity. (v) Convergence check:
(σ n(i++11) − H(i +1) ) ≤ δ
Where δ is a small number, e.g. 10-7.
NO: set i ← i + 1and GO TO (iv) YES: set states to converged values and exit
σ n +1 = σ n(i++11)
p εn +1 = εnp+(i1+1) p ε n +1 = ε np+( i1+1)
The cutting-plane algorithm described above can be interpreted geometrically as shown in Figure 7.2.
σ ntrial +1 σ n(1)+1 σ n(2) +1
σ n +1 σn Figure 7.6 Geometrical interpretation of the cutting plane algorithm. The trial ) stress state σ n(trial +1 is returned iteratively to the yield surface.
124
The linearized yield function as shown in Equation (7.21) defines a tangent “cut” at each new iteration of the yield function, and then the new variables are projected to define the next iteration. It was found that this implementation of the code, is relatively fast and a converged solution is reached within 1-3 iterations, which is considered efficient for large-scale finite element problems. Since this implementation of the constitutive equations is used in the explicit part of LSDyna, there is no need to calculate the consistent tangent modulus since it is only required in the implicit formulations. Accuracy of the developed user material subroutine (UMAT) was tested by comparing the calculated values of the plastic anisotropy parameters (R0, R45 and R90) using the developed code with the ones extracted from the experimental uniaxial tests. Results of this comparison were presented in Chapter 4 and were shown to be satisfactory (e.g. Figure 4.23 and Figure 4.40). The algorithm used to iteratively return the trial stress state to the yield surface is found to be sensitive to the number of (strain) increments used to calculate the stress state; this is shown in Figure 7.3. As step size increases, the strain increment value decreases and vice versa. As can be seen, with a small step size of 30 (strain increment = 0.01), the predicted stress values do not match the experimental results very well. As step size is increased, the calculated stress values start to converge toward the actual stress values. The stress values calculated using 3000 steps (strain increment = 0.0001), match the experimental stress values perfectly. Therefore, an important factor affecting the accuracy of the result is the number of steps or the size of the strain increment
125
used in the stress integration algorithm. In an implicit code, where larger strain increments can be used, this integration algorithm will not perform satisfactorily. However, in an explicit finite element analysis, where small step sizes are used, to ensure convergence and numerical stability, good accuracy could be expected.
200
True Stress (MPa)
150
100
50
0 -0.05
Tensile Test UMAT- 30 Steps UMAT- 300 Steps UMAT- 3000 Steps
0
0.05
0.1 0.15 0.2 True Strain
0.25
0.3
0.35
Figure 7.7 Accuracy of the developed UMAT with respect to step size compared to true stress strain results of a uniaxial test performed at 25 °C.
126
It should be noted that planar anisotropy was incorporated into the formulation for sheet forming simulations using the plane stress version of Barlat's YLD96 and YLD2000-2d models. When the deformation of the workpiece is not limited to the plane of the sheet, it is important to impose the requirement that for the plane-stress application, the in-plane material axes have to remain in the plane of the sheet during the deformation (Tuğcu et al., 1999). In the current application using the LS-Dyna FEA code, the initial anisotropy of the material is introduced by defining two local vectors in the plane of the material. Then, in the LS-Dyna code, all transformations into the element local system are performed prior to entering into this user material subroutine. Transformations back to the global system are performed after exiting the user material subroutine (Hallquist, 1999).
127
7.2. Explicit derivation of YLD96 and its derivative
The YLD96 yield function is written as Φ = α1 S2 − S3
a
a
+ α 2 S3 − S1 + α 3 S1 − S2
a
= 2σ a
(7.28)
First, an expression for σ (σ ) must be obtained. The plane stress can be expressed as ⎡σ xx ⎤ ⎢ ⎥ ⎢σ yy ⎥ σk = ⎢ ⎥ ⎢0 ⎥ ⎢σ xy ⎥ ⎣ ⎦
(7.29)
The deviatoric stress is defined as c ⎡ (c2 + c3 ) ⎤ σ xx − 3 σ yy ⎥ ⎢ 3 3 ⎥ ⎡Sxx ⎤ ⎢ (c + c ) c ⎢ 3 1 σ 3σ ⎥ ⎢ ⎥ yy − xx ⎥ ⎢ S 3 3 ⎢ yy ⎥ ⎡ξk ⎤ = ⎢ ⎥ ⎥ = ⎢ ⎛c ⎣ ⎦ c1 ⎞ S ⎢ 2 zz ⎢ ⎥ − ⎜ σ xx + σ yy ⎟ ⎥ ⎢ 3 ⎢Sxy ⎥ ⎝ 3 ⎠ ⎥ ⎥ ⎣ ⎦ ⎢ c6σ xy ⎢ ⎥ ⎢⎣ ⎥⎦
(7.30)
Although that σ 3 = 0 , Szz ≠ 0 and therefore must be included in the calculations
(
of
the
principal
stresses.
It
should
)
Szz = − Sxx + Syy .
The principal values of the deviatoric stress are
128
be
noted
also
that
⎡ ⎢ Sxx ⎢ ⎢ ⎡ S1 ⎤ ⎢ ⎢S ⎡ηk ⎤ = ⎢S2 ⎥ = ⎢ xx ⎣ ⎦ ⎢ ⎥ ⎢⎣S3 ⎥⎦ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣
⎤ 2 ⎛ Sxx − Syy ⎞ ⎥ 2 + ⎜⎜ ⎟⎟ + Sxy ⎥ 2 2 ⎝ ⎠ ⎥ ⎥ 2 + Syy ⎛ Sxx − Syy ⎞ 2 ⎥ − ⎜⎜ ⎟⎟ + Sxy ⎥ 2 2 ⎝ ⎠ ⎥ ⎥ − ( S1 + S2 ) ⎥ ⎥ ⎥ ⎦
+ Syy
(7.31)
where S3 = -(S1+S2) because of the deviatoric nature of Sij. The principal values are ordered such that S1 ≥ S2 ≥ S3. The anisotropy coefficients αi of the yield function (7.28) are defined as
⎡ α cos2 β + α sin2 β ⎤ x y ⎥ ⎡ α1 ⎤ ⎢ ⎢ ⎢ ⎥ 2 2 [α k ] = ⎢α 2 ⎥ = ⎢ α x sin β + α y cos β ⎥⎥ ⎥ ⎢⎣α 3 ⎥⎦ ⎢ 2 2 ⎢⎣α z0 cos 2β + α z1 sin 2β ⎥⎦
(7.32)
In the above equations, c1, c2, c3, c6, αx, αy, αz0 and αz1 are coefficients that describe the anisotropy of the material. The parameter 2β represents the angle between the line OA and the axis of the principal IPE stress, S1, as shown in Figure 7.4, and is equal to:
⎛ 2Sxy 2β = tan−1 ⎜ ⎜ Sxx − Syy ⎝
129
⎞ ⎟ ⎟ ⎠
(7.33)
τ
2β
A
⎫ ⎬ Sxy ⎭
Sy P S2
O
S1
σ
Sx Figure 7.8 Determination of the anisotropy parameter 2β from Mohr’s circle. Therefore, σ (σ ) defined in Equation (7.28) can be written as: 1
(
)
1
⎧1 ⎫ ⎧1 a a a ⎫ σ (σ ) = ⎨ Ψ ⎬ a = ⎨ α1 S2 − S3 + α 2 S3 − S1 + α3 S1 − S2 ⎬ a (7.34) ⎩2 ⎭ ⎩2 ⎭ ∂σ (σ ) is obtained by applying the chain ∂σ
The derivative of the yield function rule
⎧ −1 ⎛ 1 ⎞ ⎫ ∂σ (σ ) ⎪⎪ 2 a ⎜⎝ a −1⎟⎠ ⎪⎪ 3 4 ⎛ ∂σ (σ ) ∂ηi ∂ξ j ∂σ (σ ) ∂α i ∂β ∂ξ j ⎞ σ + ⎟⎟ (7.35) =⎨ ⎬ ∑∑ ⎜⎜ a ∂σ ∂ η ∂ ξ ∂ σ ∂ α ∂ β ∂ ξ ∂ σ i j k i j k ⎪ ⎪ i j ⎝ ⎠ ⎪⎩ ⎪⎭
where k=1~3. The components of the previous equation are as follows
{α (S − S ) S − S {α (S − S ) S − S {α (S − S ) S − S
⎡a ⎢ ⎢ ∂σ (σ ) ⎢ = a ⎢ ∂ηi ⎢ ⎢a ⎣
3
1
2
1
2
1
2
3
2
3
2
3
1
3
1
}⎥⎥ }⎥⎥ ⎥ }⎥⎦
a −2
− α 2 ( S3 − S1 ) S3 − S1
a −2 ⎤
a −2
− α 3 ( S1 − S2 ) S1 − S2
a −2
a −2
− α1 ( S2 − S3 ) S2 − S3
a −2
130
(7.36)
By defining
(
2 + S −S ϕ = 4Sxy xx yy
)
2
(7.37)
then,
∂ηi ∂ξ j
⎡ 1 ⎛ Sxx − Syy ⎢ ⎜⎜ 1 + ϕ ⎢2 ⎝ ⎢ 1 ⎛ Sxx − Syy = ⎢ ⎜⎜ 1 − ⎢2 ϕ ⎢ ⎝ ⎢ −1 ⎢ ⎢⎣
⎞ ⎟⎟ ⎠ ⎞ ⎟⎟ ⎠
1 ⎛ Sxx − Syy ⎜1− ϕ 2 ⎜⎝
⎞ ⎟⎟ 0 ⎠ 1 ⎛ Sxx − Syy ⎞ ⎜1+ ⎟⎟ 0 ϕ 2 ⎜⎝ ⎠ 0 −1
⎡ (c2 + c3 ) ⎢ 3 ⎢ c ∂ξ j ⎢ − 3 ⎢ 3 = ∂σ k ⎢ ⎢ − c2 ⎢ 3 ⎢ 0 ⎣
c ⎤ − 3 0⎥ 3 ⎥ (c1 + c3 ) 0 ⎥⎥ 3 ⎥ c1 − 0⎥ ⎥ 3 c6 ⎥⎦ 0
2Sxy ⎤ ⎥ ϕ ⎥ ⎥ −2Sxy ⎥ ϕ ⎥ ⎥ 0 ⎥ ⎥ ⎥⎦
(7.38)
(7.39)
The derivative with respect to αi is
⎡S −S a⎤ 3 ⎥ ⎢ 2 ∂σ (σ ) ⎢ a⎥ = ⎢ S3 − S1 ⎥ ∂α i ⎢ a⎥ ⎢ S1 − S2 ⎥ ⎣ ⎦
∂α i ∂β
( (
) )
⎡ α y − α x sin 2β ⎤ ⎢ ⎥ ⎢ ⎥ = α x − α y sin 2β ⎢ ⎥ ⎢2 (α − α ) sin 4 β ⎥ z1 z0 ⎢⎣ ⎥⎦
131
(7.40)
(7.41)
Finally,
∂β ∂ξ j
⎡ −Sxy ⎢ ϕ2 ⎢ ⎢ S xy ⎢ ⎢ = ϕ2 ⎢ 0 ⎢ ⎢S − S yy ⎢ xx ⎢ ϕ2 ⎣
132
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
(7.42)
7.3. Explicit derivation of YLD2000-2d and its derivative
The YLD2000-2d yield function is written as Φ1 = s1 − s2
a
+ s2 − s3
a
a
+ s3 − s1 = 2σa
(7.43)
First, an expression for σ (σ ) must be obtained. The plane stress can be expressed as ⎡σ xx ⎤ ⎢ ⎥ ⎢σ yy ⎥ σk = ⎢ ⎥ ⎢0 ⎥ ⎢σ xy ⎥ ⎣ ⎦
(7.44)
The symmetric matrices are defined as ′ σ xx + L12 ′ σ yy ⎤ ⎡ X xx ′ ⎤ ⎡L11 ⎥ ⎢ ⎥ ⎢ ′ ⎥ = ⎢L21 ′ σ xx + L22 ′ σ yy ⎥ X k′ = ⎢ X yy ⎥ ⎢ X ′xy ⎥ ⎢L′ σ ⎣ ⎦ ⎢⎣ 66 xy ⎥⎦
(7.45)
′′ σ xx + L12 ′′ σ yy ⎤ ⎡ X xx ′′ ⎤ ⎡L11 ⎥ ⎢ ⎢ ⎥ ′′ ⎥ = ⎢L21 ′′ σ xx + L22 ′′ σ yy ⎥ X k′′ = ⎢ X yy ⎥ ⎢ X ′′xy ⎥ ⎢L′′ σ ⎣ ⎦ ⎢⎣ 66 xy ⎥⎦
(7.46)
and
The principal values of X ′ and X ′′ are ⎡ ⎤ 2 ⎢ X1′ + X 2′ + ⎛ X1′ − X 2′ ⎞ + X ′2 ⎥ ⎜ ⎟ 3 ⎥ 2 2 ⎡ X1′ ⎤ ⎢ ⎝ ⎠ ⎡ηk′ ⎤ = ⎢ ⎥ = ⎢ ⎥ ⎣ ⎦ 2 ⎥ ⎣ X 2′ ⎦ ⎢ X ′ + X ′ 2 − ⎛ X1′ − X 2′ ⎞ + X ′2 ⎥ ⎢ 1 ⎜ ⎟ 3 ⎥ ⎢⎣ 2 2 ⎝ ⎠ ⎦
133
(7.47)
and ⎡ ⎤ 2 ⎢ X1′′ + X 2′′ + ⎛ X1′′ − X 2′′ ⎞ + X ′′2 ⎥ ⎜ ⎟ 3 ⎥ 2 2 ⎡ X ′′ ⎤ ⎢ ⎝ ⎠ ⎡ηk′′ ⎤ = ⎢ 1 ⎥ = ⎢ ⎥ ⎣ ⎦ ′′ X 2 ⎢ ⎥ ⎣ 2⎦ ′′ ′′ ′′ ′′ ⎢ X1 + X 2 − ⎛⎜ X1 − X 2 ⎞⎟ + X ′′2 ⎥ 3 ⎥ ⎢⎣ 2 2 ⎝ ⎠ ⎦
(7.48)
Therefore, σ (σ ) defined in Equation (7.43) can be written as 1
(
)
1
⎧1 ⎫ ⎧1 a a a ⎫ σ (σ ) = ⎨ Ψ ⎬ a = ⎨ X1′ − X 2′ + 2 X 2′′ + X1′′ + 2 X1′′ + X 2′′ ⎬ a ⎩2 ⎭ ⎩2 ⎭ The derivative of the yield function
(7.49)
∂σ (σ ) is obtained by applying the chain ∂σ
rule
⎧ −1 ⎛ 1 ⎞ ⎫ ∂σ (σ ) ⎪⎪ 2 a ⎜⎝ a −1⎟⎠ ⎪⎪ 2 3 ⎛ ∂σ (σ ) ∂ηi′ ∂X ′j ∂σ (σ ) ∂ηi′′ ∂X ′′j ⎞ + ⎟ σ =⎨ ⎬ ∑∑ ⎜⎜ ′ ′ ′′ ′′ X σ η σ η a X ∂σ ∂ ∂ ∂ ∂ ∂ ∂ i j k i j k⎟ ⎪ ⎪ i j ⎝ ⎠ ⎪⎩ ⎪⎭
(7.50)
where k=1~3. The components of the previous equation are as follows
{
} ⎤⎥⎥ }⎥⎥⎦
⎡a X ′ − X ′ X ′ − X ′ a − 2 ( 1 2) 1 2 ∂σ (σ ) ⎢ =⎢ ∂ηi′ ⎢ −a ( X ′ − X ′ ) X ′ − X ′ a − 2 1 2 1 2 ⎢⎣
{
{ {
⎡a 2 X ′′ + X ′′ 2 X ′′ + X ′′ a − 2 + 2 2 X ′′ + X ′′ 2 X ′′ + X ′′ a − 2 ( 2 1) 2 1 ( 1 2) 1 2 ⎢ ∂σ (σ ) =⎢ ∂ηi′′ ⎢a 2 ( 2 X ′′ + X ′′ ) 2 X ′′ + X ′′ a − 2 + ( 2 X ′′ + X ′′ ) 2 X ′′ + X ′′ a − 2 2 1 2 1 1 2 1 2 ⎣⎢
By defining
134
(7.51)
}⎥⎥⎤ (7.52) }⎥⎦⎥
2
⎛ X ′ − X 2′ ⎞ ϕ′ = ⎜ 1 + X 3′2 ⎟ 2 ⎝ ⎠ (7.53) 2
⎛ X ′′ − X 2′′ ⎞ 2 ϕ ′′ = ⎜ 1 ⎟ + X 3′′ 2 ⎝ ⎠ then,
∂ηi′ ∂X ′j
⎡ 1 ⎛ ( X1′ − X 2′ ) ⎞ ⎢ ⎜1+ ⎟ 2ϕ ′ ⎢2 ⎝ ⎠ =⎢ ⎢ 1 ⎛ 1 − ( X1′ − X 2′ ) ⎞ ⎟ ⎢ ⎜ 2ϕ ′ ⎠ ⎣2 ⎝
1 ⎛ ( X1′ − X 2′ ) ⎞ ⎜1− ⎟ 2⎝ 2ϕ ′ ⎠
X 3′ ⎤ ⎥ ϕ′ ⎥ ⎥ X 3′ ⎥ 1 ⎛ ( X1′ − X 2′ ) ⎞ ⎜1+ ⎟ − ϕ ′ ⎦⎥ 2⎝ 2ϕ ′ ⎠
∂ηi′′ ∂X ′′j
⎡ 1 ⎛ ( X1′′ − X 2′′ ) ⎞ ⎢ ⎜1+ ⎟ 2ϕ ′′ ⎢2 ⎝ ⎠ =⎢ ⎢ 1 ⎛ 1 − ( X1′′ − X 2′′ ) ⎞ ⎟ ⎢ ⎜ 2ϕ ′′ ⎠ ⎣2 ⎝
1 ⎛ ( X1′′ − X 2′′ ) ⎞ ⎜1− ⎟ 2⎝ 2ϕ ′′ ⎠
(7.54)
and
X 3′′ ⎤ ⎥ ϕ ′′ ⎥ ⎥ X 3′′ ⎥ 1 ⎛ ( X1′′ − X 2′′ ) ⎞ ⎜1+ ⎟ − ϕ ′′ ⎦⎥ 2⎝ 2ϕ ′′ ⎠
(7.55)
Finally, L′ 0 ⎤ ⎡ L′ ∂X ′j ⎢ 11 12 = L21 ′ ′ L22 0 ⎥⎥ ⎢ ∂σ k ⎢⎣ 0 ′ ⎥⎦ 0 L66
(7.56)
0 ⎤ ⎡ L′′ L′′ ∂X ′′j ⎢ 11 12 = L21 ′′ L22 ′′ 0 ⎥⎥ ∂σ k ⎢ ⎢⎣ 0 ′′ ⎥⎦ 0 L66
(7.57)
and
135
CHAPTER 8
DISCUSSION OF RESULTS: NUMERICAL Vs. EXPERIMENTAL
Finite element analysis of pure stretch experiments was performed with the thermal-structural
finite
element
model
described
previously
using
the
temperature dependant user material subroutine (UMAT) of both Barlat’s YLD96 and Barlat’s YLD2000-2d anisotropic material models. The purposes of the numerical analysis were first to check the validity of the assumption that thermal strains are negligible in thermoforming process, and then to verify the accuracy of both the FEM model and the developed user material subroutine to predict failure in the aluminum sheet at elevated temperatures. From the simulation results, there was no difference in predicted deformation behavior or failure location between the two material models. YLD2000-2d model, however, proved to be more efficient than YLD96 with respect to solution time. This is due to the simplicity of the model. The simulation results shown here are from both models. Experimental test results for the various temperatures are shown in Figure 5.9. Temperatures of the dies, punch and blank were measured with thermocouples as listed in Table 8.1. Since the material model is temperature dependant, these temperature values were input to the numerical model to insure accurate analysis corresponding to the experimental tests.
136
Table 8.2 Measured temperatures of dies, punch and blank during experiments. Test Die Punch Blank Temperature Temperature Temperature Temperature °C °C °C °C 25
25
25
25
38
38
38
38
93
93
90
92
149
149
120
142
177
176
144
170
204
204
170
200
Figure 8.1 shows the results of the coupled thermo-mechanical pure stretch simulation at 25°C (77°F). The major and minor strains of each element were projected on the FLD calculated at the same temperature in Figure 8.2. The punch depth at which the sheet failed in this analysis was 28 mm (1.1 in). This compares well with the experimental data shown in Figure 5.4 of about 27 mm (1.06 in) for the test at 25°C (77°F).
137
Failure
Figure 8.9 Finite element results from a fully coupled thermo-mechanical simulation of thermo-forming at 25°C (77°F). Failure location is shown. The sheet failed at a punch depth of 28 mm (1.1 in).
138
Figure 8.10 Distribution of the major and minor strains on the FLD. The sheet failed where the strains crossed the FLD curve.
Figures 8.3 and 8.4 show results of the coupled thermo-mechanical simulation at 93°C (200°F), Figures 8.5 and 8.6 show results of the coupled thermo-mechanical simulation at 149°C (300°F), and Figures 8.7 and 8.8 the results at a temperature of 204°C (400°F). All simulations used temperatures for dies and sheet as specified in Table 8.1. Failure locations and strain distribution are also shown. The punch depth at which the sheet failed in the 93°C (200°F) case was 32mm (1.26in), for 149°C (300°F) case was 35mm (1.37in) and for 204°C (400°F) case was 38mm (1.5in). This compares well with the experimental results of Figure 5.9. (31mm, 34mm and 37mm, respectively).
139
Failure
Figure 8.11 Results of numerical simulation of thermo-coupled FEM analysis at 93°C (200°F). Failure location shown on the graph. Punch depth = 32mm (1.26in).
140
Figure 8.12 Results of numerical simulation of thermo-coupled FEM analysis at 93°C (200°F). The graph shows the distribution of the major and minor strains on the FLD. The sheet failed where the strains crossed the FLD curve.
141
Failure
Figure 8.13 Results of numerical simulation of thermo-coupled FEM analysis at 149°C (300°F). Failure location shown on the graph. Punch depth = 35mm (1.37in).
142
Figure 8.14 Results of numerical simulation of thermo-coupled FEM analysis at 149°C (300°F). The graph shows the distribution of the major and minor strains on the FLD. The sheet failed where the strains crossed the FLD curve.
143
Failure
Figure 8.15 Results of numerical simulation of thermo-coupled FEM analysis at 204°C (400°F). Failure location shown on the graph. Punch depth = 38mm (1.5in).
144
Figure 8.16 Results of numerical simulation of thermo-coupled FEM analysis at 204°C (400°F). The graph shows the distribution of the major and minor strains on the FLD. The sheet failed where the strains crossed the FLD curve. In Figure 8.9, numerical simulation results with failure location and punch depth at failure are shown at several elevated temperatures. Figure 8.10 shows the corresponding experimental results at the same temperatures. As seen from these figures, the fully coupled thermo-mechanical finite element analysis model was able to predict accurately both the failure depth and location in the sheet at various temperatures. Table 8.2 compares experimental and numerical punch depths at failure for all measured temperatures. Figure 8.11 shows a comparison between experimental and numerical results of the punch load vs. punch depth curve at several elevated temperatures. As
145
could be seen from this plot, the developed fully coupled thermo-mechanical model is capable of accurately predicting the punch load curves.
25°C, 28mm
149°C, 35mm
Failure
93°C, 32mm
Failure
204°C, 38mm
Failure
Failure
Figure 8.17 Numerical simulation results of pure stretch using the 101.6mm (4in) hemispherical punch at several elevated temperatures. Punch depth at failure and failure locations based on FLD prediction are also shown.
146
25°C, 27mm
93°C, 31mm
149°C, 34mm
204°C, 37mm
Figure 8.18 Experimental results of pure stretch using the 101.6mm (4in) hemispherical punch at several elevated temperatures. Punch depth at failure and failure locations are shown.
147
16 14
Punch Load (kN)
12 10 8 6 25°C -Experimental 25°C -Numerical 149°C-Experimental 149°C-Numerical 204°C-Experimental 204°C-Numerical
4 2 0 0
5
10
15
20
25
30
35
40
Punch Depth (mm)
Figure 8.19 Punch force vs. punch depth for the hemispherical punch stretch experiments at several elevated temperatures. Numerical results closely match the experimental results. Table 8.3 Numerical prediction of failure depth compared to experimental results.
Temperature °C
Experimental punch depth at failure, mm
Numerical punch depth at failure, mm
25
27
28
25
29
30
93
31
32
149
34
35
177
36
36
204
37
38
148
As shown from the previous results, the numerical model and the UMAT developed in this research were able to predict accurately both the behavior and failure locations of an aluminum sheet during warm forming process at several elevated temperatures.
149
CHAPTER 9
CONCLUSIONS
In this research, a temperature-dependent anisotropic material model to be used in the finite element analysis model has been developed for the aluminum alloy sheet AA3003-H111. Using experimental data from uniaxial and bulge tests in different directions, the anisotropy coefficients for two accurate material models: Barlat’s YLD96 and YLD2000-2d for several elevated temperatures in the range of 25°C- 260°C (77°F–500°F) have been calculated. Using several polynomial curve fitting functions (3rd, 4th, and 5th order), the anisotropy coefficients of the yield functions were determined as a function of temperature. From experimental tests, a strain-rate dependent hardening flow rule was also determined as a function of temperature. The curve fitting functions were verified by a comparison of the plastic anisotropy parameters (R0, R45 and R90) that were extracted from experimental tests to the values predicted from the yield function in which the anisotropy coefficients vary with temperature. Also, the normalized yield functions at several temperatures compared favorably to experimental results. After development and verification of the accuracy of the temperaturedependent anisotropic material model for the aluminum alloy 3003-H111 was successfully implemented as a user material subroutine (UMAT) in the explicit finite element code LS-Dyna. A fully coupled thermo-mechanical finite element
150
model was developed for the analysis of the warm forming of the aluminum sheet with a hemispherical punch under pure stretch condition. The thermal part of the numerical analysis was used to solve the heat transfer behavior of the aluminum sheet based on the assigned boundary conditions. Temperatures from the thermal analysis were then applied to the UMAT, and the anisotropic behavior of AA3003-H111 was then determined in the structural analysis portion of the solution. The adaptive meshing capability in LS-Dyna used in the coupled thermomechanical analysis of thermoforming produced errors at contact regions. Temperatures of new elements created using adaptive meshing became negative. As mentioned earlier, this was idinetifed by LSTC to be a bug in LSDyna version 970 which will be corrected in version 971. The modified power law was used as flow stress in the thermoforming analysis, while the FLD used for AA3003-H111 was calculated based on the Voce hardening law. This was done for following two reasons; (a) at elevated temperatures
the
aluminum
sheet
exhibited
strain
rate
characteristics
(Abedrabbo et al., 2005a), while the conventional Voce model does not account for strain rate, (b) failure locations predicted by FLD curves based on power law were incorrect, as also reported recently by Aghaie et al. (2004). Voce hardening law does not include strain rate sensitivity directly in the equation; a simple multiplication method was used to add the strain rate sensitivity in order to calculate the FLD. Although this was ok for calculating the FLD, the modified power law was preferred when using the flow stress.
151
Finite element analysis with the developed thermo-mechanical constitutive model accurately predicted both the deformation behavior and the failure location in the blank and compared favorably to the experimental results. The current research shows the importance of using both thermal analysis and an accurate anisotropic temperature-dependant material model in a fully coupled mode in order to model accurately the warm forming process. Although the current thermoforming analysis was only verified for biaxial stretching, its application to more complex parts is expected also to yield accurate results. This is because the accuracy of both YLD96 and YLD2000-2d yield functions at room temperature already has been thoroughly verified by many researchers, therefore, to expect a similar performance at elevated temperature is not unreasonable. However, to verify the accuracy of this model’s prediction further, numerical analysis of thermoforming of deep-drawn automotive parts should be performed and compared to experimental results.
152
BIBLIOGRAPHY
BIBLIOGRAPHY
Abedrabbo, N., Pourboghrat, F., Carsley, J., In Press 2005a. Forming of aluminum alloys at elevated temperatures –Part 1: Material characterization. Int. J. of Plasticity. Abedrabbo, N., Pourboghrat, F., Carsley, J., In Press 2005b. Forming of Aluminum Alloys at Elevated Temperatures – Part 2: Numerical Modeling and Experimental Verification. Int. J. Plasticity. Abedrabbo, N., Zampaloni, M., Pourboghrat, F., 2005c. Wrinkling control in aluminum sheets using stamp hydroforming. Int. J. Mech. Sci. 47 (3), 333358. Aghaie, M., Mahmudi, R., 2004. Predicting of plastic instability and forming limit diagrams. Int. J. Mech. Sci. 46, pp. 1289-1306. Armero, F., Simo, J. C., 1993. A priori stability estimates and unconditionally stable product formula algorithms for nonlinear coupled thermoplasticity. Int. J. of Plasticity 9(6) 749-782. Auricchio, F., Taylor, R. L., 1999. A return-map algorithm for general associative isotropic elasto-plastic materials in large deformation regimes. Int. J. of Plasticity, 15 (12), 1359-1378. Ayres, R. A., Wenner, M. L., 1979. Strain and strain-rate hardening effects in punch stretching of 5182-O aluminum at elevated temperatures. Met. Trans. A 10, pp. 41-46. Ayres, R. A., 1979. Alloying aluminum with magnesium for ductility at warm temperatures (25 to 250°C). Met. Trans. A 10, pp. 849-854. Baida, M., Gelin, J. C., Ghouati, O., 1999. Modeling the hydroforming of thin metallic components. In: Proceedings of the Seventh International Symposium on Plasticity and Its Current Applications (PLASTICITY ’99), Edited by Khan A., Cancun, Mexico, January 5-13, pp. 293-296.
154
Baldwin W. M Jr., Howald T. S., 1947. Folding in the cupping operation. In: Transactions of A.S.M. 38, pp. 757-788. Barlat, F., Lege, D. J., Brem, J. C., 1991a. A six-component yield function for anisotropic metals. Int. J. Plasticity 7, p. 693-712. Barlat, F., Chung, K., 1993. Anisotropic potentials for plastically deforming metals. Modeling and Simulation in Materials Science and Engineering 1, pp. 403–416. Barlat, F., Maeda, Y., Chung, K., Yanagawa, M., Brem, J. C., Hayashida, Y., Lege, D. J., Matsui, K., Murtha, S. J., Hattori, S., Becker, R. C., Makosey, S., 1997a. Yield function development for aluminum alloy sheets. J. Mech. Phys. Solids 45 (11/12), pp. 1727–1763. Barlat, F., Becker, R. C., Hayashida, Y., Maeda, Y., Yanagawa, M., Chung, K., Brem, J.C., Lege, D. J., Matsui, K., Murtha, S. J., Hattori, S., 1997b. Yielding description of solution strengthened aluminum alloys. Int. J. of Plasticity 13, 385-401. Barlat, F., Brem, J. C., Yoon, J. W., Chung, K., Dick, R. E., Lege, D. J., Pourboghrat, F., Choi, S,-H., Chu, E., 2003. Plane stress yield function for aluminum alloy sheets-part 1: theory. Int. J. of Plasticity 19, pp. 1297-1319. Becker, R. C., 1998. Private information, Alcoa Technical Center, Pennsylvannia, December 1998. Boogaard, A. H. van den, Bolt P.J., Werkhoven, R.J., 2001. Aluminum sheet forming at elevated temperatures. In: Simulation of Materials Processing: Theory, Methods and Applications, Swets & Zeltinger (Eds.), Lisse, pp. 819824. Boogaard, A. H. van den, 2002. Thermally enhanced forming of aluminum sheet: modeling and experiments. The Netherlands, Ponson & Looijen. Čanađija, M., Brnić, J., 2004. Associative coupled thermoplasticity at finite strain with temperature-dependent material parameters. Int. J. of Plasticity 20 (10), 1851-1874.
155
Cao J., Boyce M. C., 1997. A Predictive Tool for Delaying Wrinkling and Tearing Failures in Sheet Metal Forming. J. of Eng. Mat. Technology 119, pp. 354365. Chung K., Shah K., 1992. Finite element simulation of sheet metal forming for planer anisotropic metals. Int. J. of Plasticity 8, pp. 453-476. Chung, K., Richmond O., 1993. A deformation theory of plasticity based on minimum work paths. Int. J. of Plasticity 9(8), 907-920. Chung, K., Lee, S. Y., Barlat, F., Keum, Y. T., Park, J. M., 1996. Finite element simulation of sheet forming based on a planar anisotropic strain-rate potential. Int. J. of Plasticity 12(1), 93-115. Cleveland, R. M., Ghosh, A. K., Bradley, J. R., 2002. Comparison of superplastic behavior in two 5083 aluminum alloys. Materials Science and Engineering A 351, pp. 228-236. Clift, S. E., Hartley, P., Sturgess, C. E. N., Rowe, G. W., 1990. Fracture prediction in plastic deformation process. Int. J. of Mech. Sci. 32 (1), pp. 117. Doege E., El-Dsoki T., Seibert, D., 1995. Prediction of necking and wrinkling in sheet-metal forming”, J. of Mat. Proc. Tech. 50, pp. 197-206. Garrett, R.P., Lin, J., Dean, T. A., 2005. An investigation of the effects of solution heat treatment on mechanical properties for AA6xxx alloys: experimentation and modeling. Int. J. of Plasticity 21 (8), 1640-1657. Green, D. E., Neale, K. W., MacEwen, S. R., Makinde, A., Perrin, R., 2004. Experimental investigation of the biaxial behaviour of an aluminum sheet. Int. J. of Plasticity 20 (8-9), 1677-1706. Geckler J. W., 1924. Plastisches knicken der wandung von hohlzylnder und einige faltungsverscheinungen an schalen und blechen. Ziet. Andew. Math. Mech. 8, p. 341-352.
156
Gelin, J. C., Delassus, P., Fontaine, J. F., 1994. Experimental and numerical modeling of the effects of process parameters in the aquadraw deep drawing. J. of Mat. Proc. Tech. 45, pp. 329-334. Gelin, J. C., Ghouati, O., Paquier, P., 19988. Modeling and control of hydroforming processes for flanges forming. In: CIRP Annals - Manufacturing Technology , Hallwag Publ. Ltd, Berne, Switzerland 47 (1), pp. 213-216. Greene, D. L., DiCicco J., 2000. Engineering-economic analyses of automotive fuel economy potential in the United States. ORNL/TM-2000/26, Oak Ridge National Laboratory, Oak Ridge, TN. Gronostajski, Z., 2000. The constitutive equations for FEM analysis. J. Mat. Proc. Tech. 106, pp. 40-44. Håkansson, P., Wallin, M., Ristinmaa, M., 2005. Comparison of isotropic hardening and kinematic hardening in thermoplasticity. Int. J. of Plasticity 21 (7), 1435-1460. Han, C. S., Chung, K., Wagoner, R. H., Oh, S. I., 2003. A multiplicative finite elasto-plastic formulation with anisotropic yield functions. Int. J. of Plasticity 19(2), 197-211. Hallquist, J. O., 1999. Ls-Dyna User’s manual, Livermore Software Technology Corporation, CA. Hartley, P., Pillinger, I., Sturgess, C., 1992. Numerical modeling of material deformation processes research development and applications. SpringerVerlag. Hashiguchi, K., 2005. Generalized plastic flow rule. Int. J. of Plasticity 21 (2), 321-35. Hershey, A. V., 1954. The plasticity of an isotropic aggregate of anisotropic face centred cubic crystals. J. Appl. Mech. Trans. ASME 21, 241-249. Hill, R., 1948. A theory of the yielding and plastic flow of anisotropic metals. In: Proc. Royal Soc. London Series A 193, pp. 281-297.
157
Hosford, W. F., 1972. A generalized isotropic yield criterion. J. Appl. Mech. Trans. ASME 39, 607-609. Hosford, W. F., Caddell, R. M., 1983. Metal forming – mechanics and metallurgy”, Prentice-Hall, Inc. Hsu, T. C., Hsieh, S. J., 1996. Theoretical and experimental analysis of failure for the hemisphere punch hydroforming processes. J. of Man. Sci. and Eng. 118, pp. 434-438. Kawaka, M., Olejnik, L., Rosochowski A., Sunaga H., Makinouchi A., 2001. Simulation of wrinkling in sheet metal forming. J. of Mat. Processing Tech. 109, pp. 283-289. Keum, Y. T., Ghoo, B. Y., 2001. 3-Dimensional finite element analysis of nonisothermal forming processes for non-ferrous sheets. In: Proc. of the 7th Int. Conf. on Num. Meth. in Industrial Forming Processes. Swets & Zeltinger (Eds.), Lisse, 813-818. Keum, Y. T., Ghoo, B. Y., 2002. Anisotropy at high temperatures. J. Ceramic Proc. Res. 3 (3), pp. 178-181. Kim, J. B., Yoon, J. W., Yang, D. Y., Barlat, F., 2001. Investigation into wrinkling behavior in the elliptical cup deep drawing process by finite element analysis using bifurcation theory. J. of Mat. Proc. Tech. 111, pp. 170-174. Kim, J. B., Yang, D. Y., Yoon, J. W., Barlat, F., 2000. The effect of plastic anisotropy on compressive instability in sheet metal forming. Int. J. of Plasticity 16, pp. 649-676. Lademo, O. G., 1999. Engineering models of elastoplasticity and fracture for aluminum alloys. PhD. Thesis, Norvegian Institute of Science and Technology, Trondheim, Norway, 1999. Lademo, O. G., Hopperstad, O. S., Langseth, M., 1999. An evaluation of yield criteria and flow rules for aluminum alloys. Int. J. Plasticity 15, 191-208. Levenberg, K., 1944. A Method for the Solution of Certain Problems in Least Squares. Quart. Appl. Math. 2, 164-168.
158
Li, D., Ghosh, A., 2004. Biaxial warm forming behavior of aluminum sheet alloys. J. Mat. Proc. Tech. 145, pp. 281-293. Li, D., Ghosh, A., 2003. Tensile deformation behavior of aluminum alloys at warm forming temperatures. Mat. Sci. Eng. A. 352, pp. 279-286. Lo, S. W., Hsu, T. C., Wilson, W. R. D., 1993. An analysis of the hemisphericalpunch hydroforming process. J. of Mat. Proc. Tech. 37, pp. 225-239. Marciniak, Z., Kuczynski, K., 1967. Limit strains in the processes of stretchforming sheet metal. Int. J. Mech. Sci 9 (9), pp. 609-612. Marquardt, D., 1963. An Algorithm for Least-Squares Estimation of Nonlinear Parameters. SIAM J. Appl. Math. 11, 431-441. Masuda, T., Kobayashi, T., Wang, L., Toda, H., 2003. Effects of strain rate on deformation behavior of A6061-T6 aluminum alloy. Mat. Sci. Forum (426432), pp. 285-290. McClintock, F. A., 1968. A criterion for ductile fracture by the growth of holes. J. of App. Mech. 35, pp. 363-371. Naka, T., Torikai, G., Hino, G., Yoshida, F., 2001. The effects of temperature and forming speed on the forming limit diagram for type 5083 aluminummagnesium alloy sheet. J. Mat. Proc. Tech. 113, pp. 648-653. Naka, T., Nakayma, Y., Uemori, T., Hino, R., Yoshida, F. 2003. Effects of temperature on yield locus for 5083 aluminum alloy sheet. J. of Mat. Processing Tech. 140, pp. 494-499. Ortiz, M., 1981. Topics in constitutive theory for inelastic solids. Ph.D. Dissertation. Dept. of Civil Engineering, Univ. of California, Berkeley, CA. Ortiz, M., Simo, J. C., 1986. An analysis of a new class of integration algorithms for elastoplastic constitutive relations. Int. J. Num. Meth. Eng. 23, pp. 353366.
159
Painter, M. J., Pearce, R., 1980. The elevated temperature behavior of some AlMg alloys. Les Mâemoires Scientifiques De La Revue Mâetallurgie 77, pp. 617-634. Rice, J. R., Tracey, D. M., 1969. On the ductile enlargement of voids on triaxial stress fields. J. of Mech. Phy. and Solids 17, pp. 201-217. Robinson, J. M., Shaw, M. P., 1994. Microstructural and mechanical influences on dynamic strain aging phenomena. Int. Mat. Reviews 39, 113-122. Schroth, J. G., 2004. General Motors’ quick plastic forming process. In: Taleff, E., Friedman, P.A., Krajewski, P.E., Mishra, R.S., Schroth, J.G., (Eds.), TMS, Warrendale, Advances in Superplasticity and Superplastic Forming, PA, 920. Senior B .W., 1981. Flange wrinkling in deep-drawing operations. J. of Mech. Phy. Solids 4, pp. 235-246. Shang, H. M., Qin, S., Tay, C. J., 1997. Hydroforming sheet metal with intermittent changes in the draw-in condition of the flange. J. of mat. Processing Tech. 63, pp. 72-76. Simo, J. C., Ortiz, M., 1985. A unified approach to finite deformation elastoplastic analysis based on the use of hyperelastic constitutive equations. Comp. Meth. Appl. Mech. Eng. 49, pp. 221-245. Simo, J. C., Hughes, T. J. R., 1998. Computational Inelasticity. Springer, New York, pp. 143-149. Stoughton, T. B., 2001. Stress-based Forming Limits in Sheet Metal Forming. J. Eng. Mat. Tech. 123, 417- 422. Stoughton, T., Zhu, X., 2004. Review of theoretical models of the strain-based FLD and their relevance to the stress-based FLD. Int. J. of Plasticity 20 (8-9), 1463-1486. Szabo, 2001. Private communication at Euromech Colloquium 430, Formulations and Constitutive Laws for Very Large Strains, Prague, Czech Republic, Oct. 2001.
160
Takata, K., Ohwue, T., Saga, M., Kikuchi, M., 2000. Formability of Al-Mg alloy at warm temperature. Mat. Sci. Forum (331-337), pp. 631-636. Taleff, E. M., Nevland, P. J., Krajewski, P.E., 2001. Tensile ductility of several commercial aluminum alloys at elevated temperatures. Met. Trans. 32A, pp. 1119-1130. Taleff, E. M., Henshall, G. A., Nieh, T. G., Lesuer, D. R., Wadsworth, J., 1998. Warm-temperature tensile ductility in Al-Mg alloys. Met. Trans. 29A, pp. 10811091. Tirosh, J., Yossifon, S., Eshel, R., Betzer, A., 1977. Hydroforming process of uniform wall thickness products. ASME J. of Eng. for Industry 99, pp. 685691. Triantafyllidis, N., Needlmen, A., 1980. An analysis of wrinkling in the swift cup test. J. of Eng. Mat. Tech. 102, pp. 241-248. Tuğcu, P., Neale, K. W., 1999. On the implementation of anisotropic yield functions into finite strain problems of sheet metal forming. Int. J. of Plasticity 15(10), 1021-1040. Vial-Edwards, C., 1997. Yield loci of FCC and BCC sheet metals. Int. J. of Plasticity 13 (5), 521-531. Wagoner, R. H., Nakamachi, E., Germain, Y., 1988. Analysis of sheet forming operations using the finite element method. In: Proceedings of IDDRG working groups, Toronto, p.1. Wagoner, R. H., Chenot, J-L, 1996. In: Fundamentals of Metal Forming, John Wiley & Sons, New York, pp. 10-31. Wang, X., Cao, J., 2000. On the predication of side-wall wrinkling in sheet metal forming processes. Int. J. of Mech. Sc. 42, pp. 2369-2394. Wilson, D., 1990. Representation of material behavior in finite element methods of modeling sheet forming processes. In: International Materials Review, Rev. 35, pp.329.
161
Wriggers, P., Miehe, C., Kleiber, M., Simo, J. C., 1992. On The Coupled Thermomechanical Treatment of Necking Problems via Finite-Element Methods. Int. J. Num. Methods Eng. 33 (4), 869-883. Yoa, H., Cao, J., 2002. Prediction of forming limit curves using an anisotropic yield function with prestrain induced backstress. Int. J. Plasticity 18, 1013– 1038. Yoon, J. W., Yang, D. Y., Chung, K., Barlat, F., 1999. A general elasto-plastic finite element formulation based on incremental deformation theory for planar anisotropy and its application to sheet metal forming. Int. J. of Plasticity 15(1), 35-67. Yoon, J. W., Barlat, F., Dick, R., Chung, R., Kang, T. J., 2004. Plane stress yield function for aluminum alloy sheets—part II: FE formulation and its implementation. Int. J. of Plasticity 20 (3), 495-522. Yoshida, K., Hayashi, H., Hirata, M., Hira, T., Ujihara, S., 1981. Yoshida buckling test”, In: IDDRG Paper, DDR/WG III/81, Japan. Yossifon, S., Tirosh, J., 1988. On the permissible fluid-pressure path in hydroforming deep drawing processes - analysis of failures and experiments. J. of Eng. for Industry 110, pp. 146-152. Yossifon, S., Tirosh, J., 1985. Buckling prevention by lateral fluid pressure in deep drawing. Int. J. of Mech. Sci. 27, pp. 177-185. Yossifon, S., Tirosh, J., 1985. Rupture instability in hydroforming deep-drawing process. Int. J. of Mech. Sci. 27, pp. 559-570. Yossifon, S., Tirosh, J., Kochavi, E., 1984. On suppression of plastic buckling in hydroforming processes. Int. J. of Mech. Sci. 26, pp. 389-402. Youssef, Y., Denault, J., 1998. Thermoformed glass fiber reinforced polypropylene: microstructure, mechanical properties and residual stresses. Polymer Composites 19 (3), pp. 301-309.
162
Zampaloni, M., Abedrabbo, N., Pourboghrat, F., 2003. Experimental and numerical study of stamp hydroforming of sheet metals. Int. J. Mech. Sci. 45 (11), pp. 1815-1848. Zener, C., Hollomon, J. H., 1944. Effect of strain rate upon plastic flow of steel. J. of Appl. Phys. 15, p.22.
163
