Identification and validation of constitutive model and fracture criterion for AlMgSi alloy with application to sheet forming

Materials and Design 30 (2009) 3005–3019 Contents lists available at ScienceDirect Materials and Design journal homepage: www.elsevier.com/locate/matdes Identification and validation of constitutive model and fracture criterion for AlMgSi alloy with application to sheet forming O.-G. Lademo a,c,*, O. Engler b,c, S. Keller b,c, T. Berstad a,c, K.O. Pedersen a,c, O.S. Hopperstad c a SINTEF Materials and Chemistry, Applied Mechanics and Corrosion, Richard Birkelandsvei 2, 7465 Trondheim, Norway Hydro Aluminium Deutschland GmbH, R&D Center Bonn, P.O. Box 2468, D–53104 Bonn, Germany c Structural Impact Laboratory (SIMLab), Centre for Research-based Innovation, Department of Structural Engineering, NTNU, Trondheim, Norway b a r t i c l e i n f o Article history: Received 17 October 2008 Accepted 23 December 2008 Available online 30 December 2008 Keywords: Instability Shear fracture Ductile fracture FEM FLC Polycrystalline analysis a b s t r a c t The present work addresses a modelling strategy for computer-aided formability evaluation of sheet metals and sheet metal forming operations using shell element analysis. The mechanical behaviour of the rolled AlMgSi alloy AA6016 in temper T4 is characterized by a series of different tests, including uniaxial tension tests, shear tests, plane-strain tension tests, plane-strain bending tests, bulge tests, square-cup drawing tests, and Nakajima and Marciniak–Kuczynski (M–K) formability tests. In addition, texturebased polycrystal yield surface (PCYS) calculations are carried out for identification of the anisotropic yield surface of the constitutive model. The uniaxial tension tests, shear tests and PCYS calculations are used to identify the material constants of the constitutive equation and failure criterion, while the remaining tests are left for validation. Finite element models of the tests are established and simulations are performed and compared with the experimental data. Ó 2009 Elsevier Ltd. All rights reserved. 1. Introduction In a variety of applications there is a demand for optimized components made out of sheet materials, often requiring exploitation of the material to the verge of strain localization and material failure. In particular, robust design and production of light but crashworthy structural components in aluminium for the automotive industry are challenging tasks, involving development of alloys and manufacturing processes, structural design and crashworthiness analysis. Finite element tools are routinely used to assess the feasibility of often complex industrial forming operations. The widely used industrial strategy is to perform the FE analyses without representing potential material failure, to post-process the analysis results and map the final principal strains into a forming limit diagram (FLD), and to compare the predicted strains with an experimentally determined forming limit curve (FLC) – considered being the limit beyond which the material will fail. The experimental FLC is determined through the use of a formability test set-up, e.g. the ones due to Marciniak and Kuczynski [35] (M–K) or Nakajima et al. [36], where the material is subjected to a set of more or less proportional deformation paths under membrane loading. Even * Corresponding author. Address: SINTEF Materials and Chemistry, Applied Mechanics and Corrosion, Richard Birkelandsvei 2, 7465 Trondheim, Norway. Tel.: +47 920 52 901. E-mail address: odd.g.lademo@sintef.no (O.-G. Lademo). 0261-3069/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.matdes.2008.12.020 though this method is proven powerful and has found widespread industrial utilization it has several deficiencies as discussed below. Another highly attractable strategy is to aim for direct prediction of material failure in the FE tool used to perform the forming analyses. Several distinct failure-related phenomena may contribute to the failure process of metal sheets. Teirlinck et al. [46] distinguished four different mechanisms for the failure under uniaxial tension: (i) plastic failure, (ii) cleavage and brittle intergranular fracture, (iii) ductile fracture and (iv) shear fracture. Direct prediction of failure by means of an FE tool generally requires appropriate representation of each relevant phenomenon. This is also reflected in the modelling approach proposed by Hooputra et al. [20]. Metals used in crash relevant automotive parts are usually ductile, and thus failure mechanism (ii) can be disregarded. In the present work the term ‘plastic failure’ of Teirlinck et al. [46] is generalized to account for necking under any loading situation. At this stage, it is important to notice that the instability modes in tension (i.e. plastic failure) and compression (buckling and folding) are here assumed to be represented by the finite element formulation and the elasto-(visco)plastic constitutive law. This necessitates proper attention to the constitutive model and may require the introduction of material or geometrical perturbations in the finite element mesh to trigger the instability modes and thereby getting realistic results. Furthermore, a relatively dense mesh with characteristic element size in the order of the sheet thickness is required to represent plastic failure. The reason is that 3006 O.-G. Lademo et al. / Materials and Design 30 (2009) 3005–3019 the width of the local neck is of the same size order as the sheet thickness [23]. Ductile fracture involves nucleation, growth and coalescence of voids. In the more advanced models, the damage evolution is coupled with the elasto-(visco)plastic constitutive laws, leading to material softening at some stage during the deformation process. The models based on porous plasticity [17] and continuum damage mechanics [29] belong to this class. The uncoupled damage and failure models are simpler. Damage is assumed to evolve as a function of the stress state and the plastic straining of the material, and at a certain critical damage the material fails. However, there is no influence of damage on the underlying elasto-(visco)plastic constitutive law before failure occurs and the material looses its loadcarrying capacity all together. The models proposed by Cockcroft and Latham [9] and Oyane et al. [37] fit into this class. Here the use of the phenomenological fracture criterion of Cockcroft and Latham [9] is explored. Teirlinck et al. [46] explain the shear fracture mechanism as follows: ‘‘. . . Under the right conditions, voids which nucleate in a slip band reduce the load-bearing area of the band so much that flow localizes there. Further shear increases the area of the void in the shear band, until separation occurs in the plane of the band.” During this process the initial voids are elongated in the plane of the localized shear planes, whereby the fracture mechanism has also been termed ‘void sheeting’. For polycrystalline materials, the mechanism may be influenced by texture and grain shape, particles, strain-rate softening and damage [44,13,22,46]. The shear fracture mechanism can not be spatially discretized in shell analyses, where plane stress conditions are assumed, so that it must be represented by a fracture criterion. Plastic failure may be the dominating failure mechanism for many sheet metals under membrane loading, like in most formability tests used to determine the FLC, but it does not occur in areas of the sheet that do not experience thinning, e.g. in zones of the sheet subjected to in-plane shear or sheet bending. The ductility in such areas is rather limited by ductile fracture or, potentially, shear fracture. In such a case, the current industrial strategy, as discussed above, could be overly conservative, for instance leading to the conclusion that a particular component can not be shaped from a given material. It is further well known that the experimental FLC is strongly strain path dependent, see, e.g. [15,16], which may lead to improper decisions about forming processes involving non-proportional strain paths. On the other hand, using direct predictions by means of the FE tool, physical inhomogeneity of any kind and general constitutive models, loading modes (bending vs. membrane) and loading paths can straightforwardly be represented. We consider the experimental FLC under proportional loading to be the realization of important, but idealized load cases, and use it for validation purposes. A recent summary on research related to the strain path dependence on FLCs and forming limit stress curves (FLSCs) is given by Yoshida et al. [49] and Yoshida and Kuwabara [48]. The latter work demonstrates that transients in the stress–strain response of the material due to strain path changes may be important for proper modelling of instability under general deformation paths. Such effects are not included in the present work where only isotropic hardening is considered. Rolled, heat treatable Al alloys are increasingly used for automotive skin sheets. These alloys obtain their preference due to the combination of good formability in solid solution and high strength in artificial aged conditions. Recent development has concentrated on the age hardening 6xxx series (AlMgSi–(Cu)) alloys. The present work focuses on the mechanical behaviour of the alloy AA6016 in temper T4, which is the preferred alloy in these applications within the European automotive industry [18,19,10]. The objective is to characterize the material with respect to its constitutive response and fracture characteristics and then to validate the predictive capability of a shell-based modelling approach towards a series of generic forming and formability tests. It is generally accepted that a proper constitutive model for aluminium alloys must rely upon a realistic yield criterion, in particular if plastic failure shall be predicted. The theoretical background is provided by Sowerby and Duncan [43] and Barlat [2]. Some examples demonstrating the practical consequence of using a reasonable anisotropic yield criterion, as opposed to the traditional engineering choice of the von Mises criterion, is presented by Lademo et al. [25,27]. When an anisotropic yield criterion is introduced in the constitutive model, a validated approach for identification of the anisotropy parameters is also requested. These parameters can be identified based upon experimental data, but parameter identification through polycrystalline analyses offers an attractive alternative. In the present work, both approaches are investigated. The experimental programme covers uniaxial tension tests, shear tests, plane-strain tension tests, plane-strain bending tests, bulge tests, Nakajima and Marciniak–Kuczynski formability tests and square-cup drawing tests. In addition, Polycrystal Yield Surface calculations (PCYS) are performed on the basis of the Visco-Plastic Self Consistent (VPSC) model and characterization by X-ray diffraction. The uniaxial tension tests, shear tests and PCYS calculations are used to identify the material constants of the constitutive equation and failure criterion, while the remaining tests are left for validation. Finite element models of the tests are established and simulations are performed and compared with the experimental data. The novelty of this work is related to the way existing models are combined, how parameters are identified, and to the broad measures taken to evaluate the resulting modelling concept. 2. Constitutive model and fracture criterion The full set of equations are presented elsewhere, see, e.g. [27,38], so that only a brief summary is given here. The ingredients of the constitutive model are an anisotropic yield criterion, the associated flow rule and a nonlinear isotropic hardening rule. The model is formulated for small elastic strains while plastic strains and rotations may be finite. To fulfil the principle of material frame indifference, a corotational stress formulation is adopted (see, e.g. [5]). The yield function f, which defines the elastic domain in stress space, is expressed in the form f ð^ r; RÞ ¼ f ð^ rÞ  ðY 0 þ RÞ  0; ð1Þ where Y0 is the reference yield stress and R is the isotropic strain  ¼ f ð^ rÞ is the effective stress, and hardening variable. In (1), r rY = Y0 + R is the flow stress (representing the strength of the material). The function f is then defined as [1] 2f m ¼ jr01 jm þ jr02 jm þ jr001  r002 jm ; ð2Þ where sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   ^ 11  a3 r ^ 22 2 a2 r ^ 212 ; þ a24 r 2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi    ^ 11  a6 r ^ 22 2 r001 r^ 11 þ r^ 22 a5 r ^ 212 : ¼ þ a27 r  2 2 r002 ^ 11 þ a1 r ^ 22 r01 a8 r  ¼ 2 r02  ð3Þ ð4Þ ^ 11 ; r ^ 22 ; r ^ 12 Þ are the stress components in- plane stress, Here ðr a1, . . . , a8 are material parameters that define the plastic anisotropy of the material, and m is a yield surface shape parameter. The Cartesian coordinate system is defined by the principal axes of anisotropy of the material, i.e. the x1-axis is parallel with the rolling O.-G. Lademo et al. / Materials and Design 30 (2009) 3005–3019 direction, the x2-axis is parallel to the in-plane transverse direction, and the x3-axis is parallel to the thickness direction of the sheet. An often assumed nonlinear isotropic strain hardening rule reads [30] 2 X RðeÞ ¼ i¼1 Q Ri ð1  expðC Ri eÞÞ; ð5Þ where e is the accumulated plastic strain and QRi and CRi are isotropic strain hardening constants. The rate form of the isotropic hardening rule is R_ ¼ HR e_ ; HR ¼ 2 X i¼1 C Ri Q Ri expðC Ri eÞ: ð6Þ It is seen from (5) that the strain hardening saturates, since P R ? QRi when e ! 1. This is not always confirmed by experiments, and thus the strain hardening rule is modified to read R_ ¼ maxðHR ; HR;min Þe_ ; ð7Þ where HR,min is identified from experiments at large plastic strains. It follows that (5) and (6) are only valid until an accumulated plastic strain e implicitly defined by HR;min ¼ 2 X i¼1 C Ri Q Ri expðC Ri e Þ: ð8Þ The utilized work hardening model, as represented by Eqs. (6)–(8), is an extension of the familiar Voce hardening rule, using two exponential terms instead of one, and with the possibility of a constant work-hardening rate at large plastic strains. The fracture criterion proposed by Cockcroft and Latham [9] is adopted in this study, i.e. W¼ Z 0 e ^ 1 ; 0Þde P W C ) r ^ ¼ 0; maxðr ð9Þ ^ 1 is the maximum principal stress and WC is a material conwhere r stant. Fracture occurs under conditions where the major principal stress is positive, and thus the strain to fracture depends on stress state (in particular the stress triaxiality), which is in agreement with experimental evidence. When the fracture criterion is fulfilled, the stress tensor in the actual integration point is set to zero. In a shell element, there are several through-thickness integration points. The user may define the number of through-thickness integration points that have to fail before the shell element entirely looses its load-carrying capacity and is eroded, i.e. the stress tensor is set to zero in all integration points. Usually this is assumed to happen when fracture occurs in the integration point located in the mid plane of the shell. In most simulations, under-integrated shell elements with one-byone quadrature in the plane of the shell are used. However, in some simulations it was necessary to use fully integrated shell elements, which have four-by-four in-plane quadrature. For these elements, it was decided to erode layer by layer of the shell element based upon the average in-plane value of W instead of dealing with erosion of each in-plane integration point. To avoid spurious strain localization and fracture when using shell elements with side lengths smaller than the width of the neck, which is in the order of the sheet thickness, the incremental thinning of the shell elements was calculated by use of an average or non-local value of the incremental plastic thickness strain. The non-local approach was originally proposed by Pijaudier-Cabot and Bazant [39] in order to solve the mesh dependence problem in materials with damage softening. In non-local damage mechanics, the damage evolution depends on the behaviour of the material within a radius of influence surrounding the integration point. Using non-local damage evolution the mesh sensitivity of fracture 3007 predictions is greatly reduced, leading to results that converge to a unique solution as the mesh is refined. In LS-DYNA [32], a non-local treatment of history variables is defined by the keyword * MAT_NONLOCAL. This option can be used with both solid and shell elements. In applying the non-local equations to shell elements, integration points lying in the same layer within the radius of influence, which is determined by a user-defined characteristic length L, are considered in the averaging procedure. In this study, a non-local value of the plastic thickness strain rate e_ p33 is used to regularize thinning instabilities. This will be denoted non-local plastic thinning. It should be noted that using the non-local plastic strain rate to compute the shell thinning violates plastic incompressibility locally, but it is convenient in the sense that it is uncoupled with the plane-stress constitutive equations. The method has been used successfully in prediction of strain localization in the heat-affected zone of welded aluminium connections by Wang et al. [47]. Another possibility would be to use a non-local value of the accumulated plastic strain e [8]. The constitutive model and fracture criterion have been implemented in LS-DYNA by Berstad et al. [6], and is available as material model No. 135 in the commercial version of this code, though with additional features than the ones which are presented here. 3. Identification method and model parameters Uniaxial tension tests, shear tests and PCYS calculations are used to identify the parameters of the constitutive equations and fracture criterion. The work hardening at small to medium plastic strain is defined from the uniaxial tensile tests in the rolling direction. The yield surface parameters are determined through additional uniaxial tensile tests in the 45° and 90° directions and PCYS calculations. Finally, the work-hardening rate at large plastic strain and the fracture parameter is estimated through numerical analyses of an in-plane shear tests oriented along the rolling direction. The nominal chemical composition of the investigated AA6016 alloy in wt% is Si: 1.0–1.5, Fe: 0.5, Cu: 0.2, Mn: 0.2, Mg: 0.25–0.6, Cr: 0.10, Zn: 0.20, Ti: 0.15, Al: balance. The aluminium alloy was produced as 1 mm thick sheet. The material was hot and cold rolled to final gauge, solution heat treated in a continuous annealing line and quenched. During room temperature storage between quenching and testing, the material will undergo natural ageing, which affects the strength and work hardening of the material. 3.1. Experimental work 3.1.1. Tensile tests Tensile tests were carried out according to the conditions described in EN 10002. The sheet was tested in three material directions (a = 0°, 45° and 90°) with respect to the rolling direction. Standard ISO tensile specimens with 80 mm original gauge length and 20 mm width were tested in a screw-driven testing machine Zwick Z100 at a nominal strain rate of 0.008 s–1. Plastic strains in the length and width directions were measured using two extensometers. Two or three duplicate tests were conducted in each direction. Since the duplicate tests showed insignificant scatter, only one representative test for each direction is presented below. The work-hardening curves (true stress vs. true strain) for the AA6016-T4 alloy are shown in Fig. 1. Only slight anisotropy in the flow stress is observed for this alloy. The flow stress ratio ra = ra/r0 as a function of accumulated plastic strain e was calculated, where ra is the true stress in a tensile test at an angle a from the rolling direction and r0 is the true stress in the rolling direction at the same level of accumulated plastic strain. The flow stress ratios were virtually constant for the AA6016-T4 alloy, and the 3008 O.-G. Lademo et al. / Materials and Design 30 (2009) 3005–3019 may be one order of magnitude larger than the strain at diffuse necking in the uniaxial tension test. 350 0 300 O 3.2. Polycrystal yield surface calculations O 45 σ [MPa] 250 90 O 200 150 100 50 0 0 0.05 0.1 0.15 0.2 0.25 ε p [-] Fig. 1. Experimental work-hardening curves in three directions given as true stress versus true strain (solid lines) and work-hardening curve in rolling direction fitted to the extended Voce law (dashed line). representative values were obtained as r0  1.00, r45 = 0.97 and r90 = 0.98. The length and width strains were measured during testing and used to calculate the plastic strain ratio Ra, which is the ratio of the plastic strain in the width direction to that in the thickness direction, both in a state of uniaxial tension. The plastic thickness strain was found from plastic incompressibility. The strain ratios in the three orientations were found to be R0 = 0.71, R45 = 0.45 and R90 = 0.79, which means that the alloy exhibits weak planar anisotropy. 3.1.2. Shear tests The specimen geometry of the in-plane shear tests is shown in Fig. 2. The specimen was clamped in hydraulic grips at each end and tested in a Dartec testing machine at a displacement rate of 1.2 mm/min. An extensometer with 75 mm gauge length was attached to the specimen to measure the elongation u of the centre region of the specimens. Duplicate tests were performed in the rolling direction of the sheet. The normalized force–elongation (s–u) curves are presented in Fig. 3. The normalized force s is the applied force divided by the nominal cross-sectional area of the shear zone. It is seen that duplicate tests give consistent results. Further, fracture is seen to occur shortly after reaching the maximum force. The fracture mode is illustrated in Fig. 4. It is important to realize that since the shear test is not prone to plastic instability, the strains in the gauge area 150 R15 R27 R1 12 30° 65.5 110 89.8 30 Fig. 2. Shear test specimen. Measures in mm. 80 R8 60° 20 Ø8 As already mentioned earlier in this paper, FE simulation of forming operations of anisotropic materials implies that the material input to the FE-code – the so-called constitutive equations – must incorporate the anisotropic properties of the sheet. This is commonly accomplished by using a suitable phenomenological yield function, e.g. the yield function due to Aretz [1]. Experimental determination of the required anisotropic materials parameters a1, . . . , a8 (Eqs. (3) and (4)) involves a number of mechanical tests – including uniaxial tensile tests in various in-plane directions, plane-strain tensile tests, shear tests, etc. – as described in Section 3.1. of this paper. As an alternative to the experimental determination of the anisotropy parameters, assuming that the anisotropy is entirely caused by the materials crystallographic texture, the yield surface and the associated flow rule can also be identified by texture-based polycrystal plasticity calculations directly from texture data. In some cases, the latter approach might be required, e.g. for materials with strong through-thickness inhomogeneity or for section shapes making it difficult to extract test coupons, while in general it offers an attractive alternative to an experimental test programme. Hybrid identification methods can also be chosen where available experimental data are used together with supplementary data obtained from the polycrystal plasticity calculations. For these reasons, it was found expedient to perform Polycrystal Yield Surface (PCYS) calculations and to evaluate the accuracy of this approach with respect to experimental data. The polycrystal plasticity calculations were performed using the Visco-Plastic Self Consistent (VPSC) formulation developed by Lebensohn and Tomé [28]. For analysis of the crystallographic texture of the sheet pole figures was measured by means of standard X-ray diffraction techniques and used to calculate the three-dimensional orientation distribution function (ODF) f(g) [40]. To be able to use texture data to predict material properties, an average texture – i.e. a texture integrated over all thickness layers – had to be determined. The ODF calculation was performed under the assumption of orthotropic sample symmetry – given by the rolling direction (RD), the transverse direction (TD) and the normal direction (ND) of the sheet – such that Euler angles are in the domain 0° 6 {u1, U, u2} 6 90°. The resulting ODF displays a quite sharp cube texture together with weak intensities of the R component (Fig. 5). Note that this alloy was developed for interior, not visible car body applications, such that ridging is not an issue [11]. For the subsequent polycrystal plasticity calculations, the ODF was discretized in an iterative procedure into sets of approximately 1000 individual orientations with equal weights. Note that for a consideration of the full polycrystal properties orthotropic sample symmetry should not be applied, i.e. 0° 6 u1 6 360°. The VPSC formulation developed by Lebensohn and Tomé [28] is a very versatile polycrystal plasticity approach which accounts for anisotropic properties of both the individual crystals and the aggregate as a whole. The deformation is simulated by imposing successive deformation increments; at each deformation step a set of boundary conditions (either strain rates or a combination of strain rates and stress components) is imposed to the (discretized) sheet texture, and the stresses and strain rates in each grain are calculated. The shear rates are used to determine changes in crystallographic orientation and to update grain shape and yield stresses in the individual grains. The overall (macroscopic) stresses and strains follow from averaging over the corresponding grain components. 3009 O.-G. Lademo et al. / Materials and Design 30 (2009) 3005–3019 a b 250 250 PCYS 200 Hybrid 150 s [MPa] s [MPa] 200 100 Experiments Analyses 50 150 100 Experiments HR, min = 0 MPa HR, min = 50 MPa HR, min = 100 MPa 50 0 0 0 0.5 1 1.5 2 2.5 0 0.5 u [mm] c 1 1.5 2 2.5 2 2.5 u [mm] d 250 250 Original mesh 200 s [MPa] s [MPa] 200 150 Refined mesh 100 150 100 Experiments Analyses Experiments WC=140 MPa 50 50 0 0 0 0.5 1 1.5 2 2.5 u [mm] 0 0.5 1 1.5 u [mm] Fig. 3. Shear test: experimental results and simulation with (a) two identification methods, (b) different minimum hardening rate, (c) coarse and refined mesh without nonlocal thinning and (d) non-local thinning and fracture criterion activated. The polycrystal yield surface was derived from the texture data with the VPSC model by probing the texture with equi-spaced strain rates contained in a given plane through the full five-dimensional stress space. The VPSC scheme offers several options to solve the equation of inclusion-matrix interaction. In the present application, best results were obtained by using a linearization scheme with a reduced effective inverse rate sensitivity, neff = 10, whereas the actual value of the inverse strain rate sensitivity n is usually larger. Two planes were probed in steps of 5°. Firstly, the conventional plane spanned by two in-plane stresses r11 and r22 is used, assuming that the through-thickness stress r33 is zero. Such yield locus computations were performed for the standard RD/TD plane (0°/90° section) and for a 45°-ND rotated frame (45°/135° section). Secondly, the yield locus was determined for a plane containing the shear stress r12. Here, the shear stress r12 is plotted versus pffiffiffi the value of ðr11 þ r22 Þ= 2. Finally, all data were normalized with the critical resolved shear stress. Calculated yield loci are depicted in Fig. 6. 3.3. Identification procedure and results The parameters of the extended Voce hardening rule were fitted to the experimental work-hardening curve in the rolling direction using the method of least squares. The parameters came out as: Y0 = 137 MPa, QR1 = 19.1 MPa, CR1 = 592, QR2 = 170 MPa and CR2 = 11.4. The fitted work-hardening curve in the rolling direction is plotted together with the experimental curve in Fig. 1. In the experiment, diffuse necking in the rolling direction occurred at a plastic strain equal to 17.5%, while the Considère criterion, dr/dep = r on the basis of the fitted curve predicts diffuse necking at 16.4%. The material Fig. 4. Shear test: experimental and predicted failure modes. 3010 O.-G. Lademo et al. / Materials and Design 30 (2009) 3005–3019 Fig. 5. Orientation distribution function identified by means of X-ray diffraction. parameter HR,min, which represents the minimum work-hardening rate, is determined from finite element analysis of the shear test, in which large plastic strains are obtained before fracture owing to the low stress triaxiality and the absence of plastic instabilities. From the uniaxial tensile tests three flow stress ratios ra and three plastic strain ratios Ra have been determined, while the yield criterion has eight free parameters a1, . . . , a8 in addition to the exponent m, requiring at least eight prescribed values (stress states at yielding or plastic flow directions) for proper parameter identification. Here, two approaches are used to determine the yield function. In the first approach, we determine a1, . . . , a8 entirely based upon the PCYS calculations. As shown in Fig. 6, the fitted yield locus with m = 6 accurately represented the PCYS data. The yield function was also fitted with m = 8, which is the default choice for aluminium alloys [33,2], but this led to less accurate representation of the PCYS data obtained using the VPSC model. In the second approach, the six experimental values are used together with two additional values extracted from the PCYS data. These were the equibiaxial flow stress ratio (i.e. the yield stress in equibiaxial tension divided by the yield stress in uniaxial tension in the rolling direction) and equibiaxial strain ratio (i.e. the ratio of the plastic strain in the rolling direction to the plastic strain in the transverse direction in equibiaxial tension). The exponent m was kept equal to 6. Fig. 6 compares the yield loci obtained with the two approaches. In the cross section shown in part (a) of the figure, only slight differences between the two approaches are observed, while a marked difference is observed in Fig. 6b. The hybrid approach gives lower shear strength. To evaluate the accuracy of the two fits of the yield function, variation of the flow stress ratio and the plastic strain ratio with tensile angle is shown in Fig. 7. The hybrid method captures the measured data, while the method entirely based upon PCYS calculations gives significant deviations from the experimental data. Table 1 gives the two sets of parameters for the anisotropic yield criterion. The shear test is used to determine the last two parameters of the constitutive model and fracture criterion, namely HR,min governing the minimum work-hardening rate and WC controlling fracture. Due to the complex stress state arising in the central shear zone, nonlinear finite element analyses are used to extract the parameters. Two different meshes were used to save computation time. A coarse mesh was used to determine HR,min, and then with the calibrated HR,min a finer mesh was used to find WC. The reason is that fracture is strongly localized and a fine mesh is needed to discretize the large strain gradients within the shear zone in the last stage before separation. The coarse mesh had 2538 fully integrated shell elements (type 16 of LS-DYNA) with a single integration point through the thickness, which gives a characteristic element size of 0.23 mm in the shear zone. The fine mesh had 8667 fully integrated shell elements in the shear zone with a characteristic element size of 0.07 mm. Simulations were first run with coarse mesh and the two set of identified parameters for the yield function without activating the fracture criterion. The minimum work-hardening rate HR,min was set to zero. The results are shown in Fig. 3a. It is seen that the calibration based on PCYS data overestimates the force, while the hybrid method gives a somewhat low estimate. Owing to the superior representation of the plastic anisotropy in uniaxial tension, it was decided to use the parameters obtained with the hybrid method in the remaining simulations. Fig. 3b shows simulations using the parameter set obtained with the hybrid method and HR,min equal to 0, 50 and 100 MPa. It is seen that HR,min = 50 MPa gives the best fit to the force–elongation curve for large strains, while the force 1.5 1.5 1 0.75 0.5 σy σ0 σxy σ0 0 0 -0.5 -0.75 -1 -1.5 -1.5 -1.5 -1 -0.5 0 σx σ0 0.5 1 1.5 -1.5 -0.75 0 0.75 1.5 √2 (σx+σy) 2σ0 Fig. 6. Yield locus obtained from PCYS calculations (circle symbols), yield locus fitted to PCYS data (dashed line) and yield locus determined by the hybrid method (solid line). Note: The yield surface shape parameter, m, is set to 6. 3011 O.-G. Lademo et al. / Materials and Design 30 (2009) 3005–3019 1.1 1.2 1.05 0.8 σα σ0 Rα 1 0.4 0.95 0.9 0 0 30 60 90 0 30 60 90 α [deg] α [deg] Fig. 7. Directional variation of the flow stress ratio and the plastic strain ratio with yield function determined from PCYS data (solid line) and by the hybrid method (dashed line). Experimental data are given by asterisks. level is somewhat underestimated. It was decided to use this value of HR,min in all subsequent simulations. Second, the fine mesh was used to determine WC. The fracture criterion is still not activated but the value of W is monitored. It is assumed that the maximum value of W within the gauge area at maximum force represents the critical value WC. It turned out that the fine mesh gave spurious strain localization leading to premature failure of the shear specimen by thinning of some elements. This is demonstrated in Fig. 3c, where it is seen that instability occurs at 1.25 mm elongation. To amend this problem, non-local plastic thinning was used with a radius of influence, L, equal to the sheet thickness using a uniform weight function. The result is shown in Fig. 3c together with the evolution of W, which indicates a critical value of W equal to WC = 140 MPa. Finally, Fig. 3d compares the experimental results with the simulations using HR,min = 50 MPa and WC = 140 MPa. To further evaluate the quality of the simulations, the experimental and predicted fracture modes are compared in Fig. 4. Both the predicted normalized force–elongation (s–u) characteristics and failure modes are seen to correlate well with the experiments. the model accurately represents the uniaxial plastic anisotropy and the strain hardening of the material (see [25,26]. The shear test characteristics are well represented through the chosen parameter identification routine, as demonstrated above. 4.1. Plane-strain tension tests The geometry of the plane-strain tension specimen tests is provided in Fig. 8. The specimens were clamped in hydraulic grips at each end, and tested in an Instron Tension/Torsion testing machine at a displacement rate of 1 mm/min. An extensometer with 27.5 mm gauge length was attached to measure the elongation u of the centre region of the specimen. Duplicate tests were successfully performed in the rolling direction of the sheet. The normalized tensile force vs. elongation (s–u) curve is presented in Fig. 9. The normalized force s is the applied force divided by the nominal cross-sectional area at minimum width. It is seen that the duplicate tests give consistent results. The finite element model of the plane-strain tension specimen consists of 3624 Belytschko-Tsay shell elements. In the gauge area, 36 elements were used across the minimum width, which gives a 4. Validation R15 12.4 60 26 R4 Plane-strain tension tests, plane-strain bending tests, bulge tests, square-cup drawing tests and Nakajima and M–K formability tests are performed for validation of the calibrated model. The results obtained from the Nakajima and M–K formability tests are evaluated towards shell-based FLC predictions. For the remaining validations tests, FE models are established, after which simulations are performed and compared with the experimental data. All simulations were run using the explicit solver of LS-DYNA, and mass or time scaling was used to get reasonable time step size and run times. In each case, the involved energies were controlled to ascertain quasi-static loading conditions, i.e. it was checked that the kinetic energy was only a small fraction of the internal energy of the system. The force-deformation characteristics and failure modes of the uniaxial tensile tests are expected to be well represented since 15 150 Fig. 8. Plane-strain tension test specimen. Measures in mm. Table 1 Parameters of anisotropic yield function. ID method m a1 a2 a3 a4 a5 a6 a7 a8 PCYS Hybrid 6 6 1.008 0.984 0.999 1.032 1.034 1.075 1.055 1.186 0.986 0.947 1.021 1.015 0.859 0.920 1.028 1.051 3012 O.-G. Lademo et al. / Materials and Design 30 (2009) 3005–3019 300 ri = 0.5t ri = t s [MPa] 200 Experiments Analyses 100 0 0 1 2 3 4 u [mm] Fig. 9. Plane-strain tension test: experimental and numerical normalized force vs. elongation curves. characteristic element size of 0.83 mm. Two simulations were run using non-local plastic thinning, one with a radius of influence L equal to the specimen thickness (as in the shear test analyses) and one with a radius of influence equal to half the specimen thickness. The predicted results are presented along with the experimental data in Fig. 9. The correlation between experiments and simulations is excellent, especially with L equal to the specimen thickness. 4.2. Plane-strain bending tests The plane-strain bending tests were carried out in accordance with the guidelines proposed by Daimler-Chrysler [4] and DIN 50 111 with a test set-up as shown in Fig. 10. Samples with size 60 mm  60 mm were machined from the sheet and tested in two directions, with the bending line being oriented parallel (90°) and perpendicular (0°) to the rolling direction. The tests were performed on a Zwick Z100 screw-driven testing machine, with five duplicates in each case. A constant roll gap of 2.04 mm was used. The starting conditions of the test were established by lowering the punch with a speed of 10 mm/min until a pre-force of 30 N was reached. The actual bending tests were performed at a speed of 20 mm/min. During the test, the applied force F and the punch displacement u were recorded. The test was automatically stopped when the force dropped by more than 15 N with respect to the maximum force, since this indicates failure initiation. The experimental scatter was in all cases very low and a representative experimental force–displacement curve for each case is presented in Fig. 11. The differences between the experiments in the 0° and 90° directions are small, though the 90° direction shows slightly higher force levels towards the end of the test. The test plates experienced fracture at a punch displacement of approximately 12 mm for both specimen orientations. Owing to symmetry of the test, a quarter part of the test set-up was modelled. The model consists of the punch and the roller support discretized with rigid shell elements and the test piece discretized with Belyschko-Tsay shell elements with one in-plane and nine through-thickness integration points. The element length along the longitudinal direction of the specimen in the mid of the gauge area, i.e. close to the bending zone, was 0.6 mm. The element length was gradually increased towards the supports. The punch was assigned a smooth velocity increase up to a constant speed within a solution time of 24 ms. Fig. 11 presents experimental and predicted force vs. displacement curves. The simulations show some oscillations, which are due to inertia effects in the explicit numerical solutions. Compared to the experiments, the simulations overestimate the force at large displacement, while the ductility is underestimated. This may also be a signal of significant inertia effect towards the end of the analyses. We note that explicit analyses of the problem at hand are prone to exhibit dynamic effects because of the small plastic zone compared with the relative large zones of elastic material that undergoes large translations and rotations. A number of analyses were performed with increased solution time, local mass scaling, alternative representation of the ‘wings’ of the specimen, without getting any improved prediction. For the presented analyses it was verified that the ratio between the global kinetic and internal energy was small and in this respect the analyses are considered being 2000 Incipient fracture 1600 F [N] 1200 800 6016-T4, 0o 6016-T4, 90o 400 0 0 4 8 12 16 u [mm] Fig. 10. Plane-strain bending test set-up [4], where sample length and width lp = bpb = 60 mm, punch radius rst = 0.1 mm, roll diameter droll = 30 mm, roll gap lab = 2a and specimen thickness a  1.0 mm. Fig. 11. Plane-strain bending test: experimental and numerical force vs. displacement curves in the rolling and transverse directions. 3013 O.-G. Lademo et al. / Materials and Design 30 (2009) 3005–3019 quasi-static and considered to be in acceptable accord with the experiments. There is hardly any difference between the predicted responses in the 0° and 90° orientations, which is in good accordance with the experiments. 16 12 Bulge tests were carried out on circular blanks of 100 mm diameter by forming a 50 mm diameter bulge as shown in Fig. 12. The blanks were locked with a blank-holder and a draw-bead. The room between the blank and the punch was filled with hydraulic oil. A cylindrical punch with 50 mm diameter was moved upwards with a constant velocity of 1.5 mm/s. During the tests, the force F necessary to move the punch and the resulting punch displacement u were recorded. After the tests, the bulge heights were measured manually. Three samples were loaded to fracture. The average bulge height and force at fracture were 17.0 mm and 29.6 kN, respectively. An additional test was performed for more detailed comparison with the finite element simulations. In this test, the blank was furnished with a grid consisting of a regular dot pattern. The test was stopped shortly before the maximum bulge height was reached, at a bulge height of 15.7 mm. After the test, the strain distribution on the formed bulge (in terms of major and minor strains as well as through-thickness strain) was determined using the ARGUS optical strain measurement system. The pressure–displacement curve is presented in Fig. 13, where the pressure is obtained by dividing the measured force by the punch area. A quarter part of the test set-up was modelled due to the initial two-fold symmetry. Hence, unsymmetrical failure modes are excluded from the solution. The blank-holder was not represented in the model. Instead, the outer periphery of the blank was constrained with respect to in-plane displacements. The die was modelled with rigid shell elements. The blank was discretized with 2700 Belyschko-Tsay shell elements with one in-plane and five through-thickness integration points. The characteristic element size was 0.5 mm. Non-local plastic thinning was utilized with a radius of influence equal to 0.75 mm. The reason for the slightly reduced radius of influence as compared to the above analyses is to reduce the number of ‘neighbours’ in the non-local calculations, in order to speed up the calculations. Loading is introduced as a smoothly increasing pressure to the outer face of the blank during a time period of 1 s, with a global mass scaling factor of 104. In the simulation, failure was caused by the local fracture criterion before onset of plastic failure, and thus the quality of the prediction relies entirely upon the identified fracture parameter WC. The obtained agreement is satisfactory. As can be seen from Fig. 13, the pressure vs. displacement curves from test and simulation compare well. The deviation in the initial phase is due to some transient dynamics in the numerical analysis. Furthermore, the experimental and numerical maximum bulge heights without Fig. 12. Bulge test set-up. p [MPa] 4.3. Bulge tests 8 Experimental Numerical 4 0 0 4 8 12 16 u [mm] Fig. 13. Bulge test: experimental and numerical pressure vs. displacement curves. failure were 17 mm and 16 mm, respectively. The predicted maximum punch pressure is somewhat less than in the experiment. 4.4. Square-cup drawing tests The square-cup drawing test has been established by Hydro R&D to investigate material damage caused by shear fracture [14]. Drawing of square cups was performed with a similar setup as described for the bulge tests. A sketch of the test set-up is presented in Fig. 14. Circular blanks with diameters D up to 100 mm were formed with a square punch with an edge length Lp = 40 mm. The rolling direction was parallel to one of the punch edges. The sheet was lubricated with lanolin and covered with a drawing foil to protect the grid applied for ex situ strain measurement. During the test, the blank was held down with a blankholder force of 7 kN. The punch was moved with a constant velocity of 8 mm/s. The force necessary to move the punch and the resulting punch displacement were recorded. Some samples were first tested to fracture qffiffiffiffiffiffiffiffiffiffiffiffiffiin order to determine the limiting drawing ratio bmax ¼ D= 4L2p =p. An additional test was performed with a sample having a grid and a slightly smaller blank diameter so as to avoid fracture. The smallest blank diameter leading to fracture initiation was 80 mm, while a reduced blank diameter of 79 mm did not cause fracture. The maximum drawing ratio was bmax = 2.0, and the maximum force at fracture was 21.9 kN. The force vs. displacement (F– u) curve from the experiment without fracture is presented in Fig. 15. On the basis of the applied grid, the strain distribution (major strain, minor strain, thickness strain) was determined for comparison with the FE simulations. Some results are presented in Fig. 16, which shows the major and minor principal strains along intersections with angles 0°, 45° and 90° with respect to the rolling direction. The finite element model of the square-cup drawing test set-up takes advantage of the initial symmetry to speed up the calculations. Accordingly, possible unsymmetrical failure modes can not be captured by the model. The model consists of the die, blankholder and punch discretized with rigid shell elements, while the blank was discretized with 7976 Belytschko-Tsay shell elements with one in-plane and seven through-thickness integration points. The characteristic element size was about 0.25 mm. The analyses were run using non-local plastic thinning with a radius of influence equal to 0.5 mm. Due to the experimental conditions, a coefficient 3014 O.-G. Lademo et al. / Materials and Design 30 (2009) 3005–3019 Fig. 14. (a) Square-cup drawing test set-up. Measures in mm. (b) Fractured sample. of friction of 0.02 was assumed between the various parts. As in the tests, the blank-holder was first moved down onto the blank until a contact force of 7 kN was achieved. Next, the punch was moved onto the blank within a quasi-static velocity regime during a period of 5 s, with a global mass scaling factor of 104. Blank diameters equal to 79 and 81 mm were simulated. Indeed, in accordance with the experiments, the simulations predicted that the blank with diameter 79 mm could be formed without fracture, while the larger blank was predicted to fracture. The force–displacement curves from test and simulation of the blank with diameter 79 mm are shown in Fig. 15 which also demonstrates the good accuracy of the numerical solution. Experimental and predicted principal strains along intersections in the rolling direction, 45° direction and transverse direction are shown in Fig. 16. Accurate predictions are obtained along intersections in the 0° and 90° directions, while some deviations are observed for the 45° direction. 30 Experimental Numerical 4.5. Formability tests and forming limit curves 20 F [kN] 4.5.1. Experimental work The Forming Limit Curve (FLC) of the material is traditionally established using either Nakajima or M–K formability tests. Here, both methods are explored. The Nakajima formability test set-up and the related extraction of an experimental FLC were determined on the basis of the guidelines of Liebertz et al. [31]. The blank is clamped between the die and the blank-holder under the action of a blank-holder force, after which a hemispherical punch with a diameter of 100 mm moves onto the plate with a velocity of about 90 mm/min. Fractured samples are presented in Fig. 17. It is seen that in most cases the local neck and fracture occur perpendicular to the direction of the major principal strain. Exceptions are specimens with negative strain ratio and the specimen in equibiaxial tension. In all the tests, the fracture occurred near the apex of the specimen. After testing the 10 0 0 10 20 30 u [mm] Fig. 15. Square-cup drawing test, blank diameter D = 79 mm: experimental and numerical force vs. displacement curves. a 0.4 c 0.4 b 0.8 Experiments Analyses 0.4 0 Experiments Analyses 0.2 Major strain ε1 , ε 2 Major strain ε1 , ε 2 ε1 , ε 2 0.2 Experiments Analyses 0 Major strain 0 Minor strain Minor strain -0.2 Minor strain -0.4 -0.4 -0.2 -0.8 0 10 20 30 lu 40 50 -0.4 0 10 20 30 lu 40 50 0 10 20 30 40 50 lu Fig. 16. Square-cup drawing test, blank diameter D = 79 mm: experimental and predicted principal strains along various intersection planes: (a) 0°, (b) 45° and (c) 90° intersections. 3015 O.-G. Lademo et al. / Materials and Design 30 (2009) 3005–3019 Fig. 17. Fractured Nakajima formability test samples. strain field was assessed by means of a speckle pattern and optical readings, again according to the guidelines of Liebertz et al. [31], and each of the specimens provides one point on the FLC. The experimental data is presented in Fig. 18. It is seen that minimum critical strain is obtained for a strain ratio of about 0.02, i.e. not under plane-strain conditions as found by classical M–K formability analysis. The reason for this deviation may be due to non-proportional strain paths during loading or to the curvature imposed to the Nakajima formability test specimens. The experimental FLC was also determined by means of doubleplate formability tests [35]. This test produces membrane straining of a flat test region, i.e. there is no superimposed bending of the test piece as with the Nakajima test. Four different plate geometries were tested, as shown in Fig. 19. All specimens were oriented along the rolling direction of the material. A minimum of two duplicate tests were performed for each of the geometries. The experimental procedure is presented elsewhere [38,27]. The present experiments are performed in an up-scaled experimental setup as compared with these past works. The samples were either marked with a 2 mm square grid by an etching technique or by a random speckle obtained by spray painting. In the former case the final strains were obtained using the automatic strain measurement equipment ASAME. The method has an accuracy of ±2% strain and is suited for measuring strains on curved surfaces and 0.6 Thinning instability Nominal strains at fracture Local fracture limit Nakajima test data MK test data 0.4 εx over a larger area. However, precaution needs to be taken when etching the grid on the surface to avoid grooves that could act as sites for strain localization. This was achieved by using short etching times. In the experiments performed using the alternative speckle pattern, strains were determined using Digital Image Correlation (DIC) analysis. The strains used to determine the formability diagrams were measured at a distance of 2.5 mm from the local instability. The experimental data is presented in Fig. 18, together with the results from the Nakajima tests. With this method, the minimum critical strain is achieved under plane-strain conditions. The numerical values of the minimum critical strain as obtained by the two methods correlate well, but in general the methods produce different experimental formability limits. Both in the second quadrant of the forming limit diagram and for equibiaxial straining the Nakajima test indicates better formability than the M–K test. Apparently, the M–K test data exhibit more scatter than the Nakajima test data. This is, however, partly due to the methods used to extract the experimental data points. Using the M–K test, a number of strain states outside the local neck are plotted for each test and some of the scatter merely reflects a somewhat inhomogeneous strain field at failure. For the Nakajima tests, a single critical strain state was extracted through the interpolation procedure proposed by Liebertz et al. [31]. For equibiaxial loading, four M– K tests were performed that actually showed strongly varying ductility. A wider scatter band under equibiaxial loading may generally be expected as it is predicted by plastic instability analysis 0.2 205 40 0 -0.2 -0.2 0 0.2 0.4 0.6 205, 165, 160, 155 εy Fig. 18. Forming limit strains from experiments and simulations. Fig. 19. Sample geometries tested in M–K double-plate formability test set-up [35]. Measures in mm. 3016 O.-G. Lademo et al. / Materials and Design 30 (2009) 3005–3019 [12]. Still, as will be seen, for the present material the ductility under equibiaxial loading does not seem to be limited by localized necking. In that case the scatter needs to be interpreted by the actual material fracture mechanisms, as discussed below. 4.5.2. Shell-based FLC predictions In previous studies [25,41,12], a shell-based Marciniak–Kuczynski analysis [35] has been implemented in LS-DYNA and validated against experiments. Predicted FLCs for the AA6016-T4 alloy are directly compared with the experimental results obtained through the Nakajima and M–K formability test set-ups, i.e. interpreting the experimental FLC as a material characteristic that is unaffected by the test set-up. A square patch of material is discretized into a grid of nel  nel square shell (or membrane) elements and plane stress states are assumed since the material is thin. The initial width of the patch is w0 and the initial thickness is t0. The patch is subjected to proportional straining defined by u1 in the x1-direction and u2 in the x2-direction. The strains are given by  e11 ¼ e cos h ¼ ln 1 þ where e ¼  u1 ; w0  e22 ¼ e sin h ¼ ln 1 þ  u2 ; w0 ð10Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e211 þ e222 and h is the angle between the x1-axis and the direction of straining. In terms of e and h, the displacements u1 and u2 read u1 ¼ w0 ðexpðe cos hÞ  1Þ; u2 ¼ w0 ðexpðe sin hÞ  1Þ: ð11Þ Since the simulations are performed by the explicit solver in the FEcode LS-DYNA, it is important to use a smooth strain–time function, viz.  e ¼ etot 1  cos  pt 2t tot  ; ð12Þ where t is the time, ttot is the termination time and etot is the strain e at termination. The displacements as function of time are found by inserting Eq. (12) into Eq. (11). It is assumed that the thickness at the element nodes is independent and normally distributed with mean value t0 and a userdefined coefficient of variation CoVt. During the biaxial stretching of the patch, the strain gradients will grow and eventually produce large localized strains in weak regions corresponding to the formation of a local neck. In order to predict the onset of localization, a criterion is required. The negative strain increment in the thickness direction is much larger within the neck than outside it when strain localization occurs. The average strain increment in the thickness direction for the whole patch, DeX 33 , can be found numerically by: DeX33 ¼ N 1X Dei N i¼1 33 ð13Þ where Dei33 is the strain increment in the thickness direction for element number i and N ¼ n2el is the total number of elements in the domain X, which is taken as the whole patch of elements. It has been assumed that the area of each element is approximately equal before localization occurs. The ratio bi ¼ Dei33 =DeX 33 between the element strain increment Dei33 for element i and the average strain increment DeX 33 in the thickness direction is used to monitor onset of strain localization. Plastic failure is assumed to occur at time tcr if any bi exceeds a critical value bC at tcr and during nit consecutive steps. By inserting the time at failure tcr into Eq. (12) and making use of Eq. (10), the forming limit strains of the patch may be found for the actual path. If the procedure is repeated for several different h, an FLC can be constructed. The calculations were performed on a patch of 50 by 50 elements and a coefficient of variation of the thickness, CoV(t), equal to 0.005. The characteristic element size was 1 mm, which means that the patch size was 50 by 50 mm2. The critical strains are defined by plastic failure (bi P bC) or fracture (Wi P WC), whatever appears first. In the simulations, bC = 2.0 and WC = 140 MPa. Nonlocal thinning was not used in these simulations. The results are shown in Fig. 18 together with the results from the Nakajima and M–K formability tests. In the figure, plastic failure refers to the predicted limit for incipient localized necking (bi P bC). The local fracture limit represents the local strains in the critical shell element at fracture (Wi P WC), while the nominal strains at fracture are the corresponding overall strains in the patch defined by Eq. (10). The local and nominal strains at fracture differ, since some strain localization occurs before fracture for all strain paths. In the major part of the FLD, the critical strains are defined by plastic failure, but in the neighbourhood of equibiaxial straining ductile fracture is critical. It is seen that the simulations are somewhat conservative for negative strain ratios. For equibiaxial straining the M–K formability data fall within the predicted fracture limit, i.e. the predictions are somewhat non-conservative, but this is not the case when the Nakajima and bulge-test data are considered. For strain ratios in the intermediate range, the simulations are quite accurate. 5. Discussion and conclusions An engineering modelling strategy for direct prediction of fracture in shell-based forming analyses is explored. A priori it is assumed that three distinct phenomena (plastic failure, ductile fracture and shear fracture) may contribute to the failure process of the AA6016-T4 material. A proper constitutive model is considered to be of general importance for accurate forming analyses [21] and a pre-requisite for realistic prediction of the phenomenon of plastic failure under arbitrary loading, Lademo et al. [25,27]. The selected constitutive model is based upon an anisotropic yield criterion for plane stress states, the associated flow rule and isotropic Fig. 20. Necking and fracture region of (a) plane-strain tensile test and (b) equibiaxial M–K formability test. O.-G. Lademo et al. / Materials and Design 30 (2009) 3005–3019 3017 Fig. 21. Close-up pictures of plane-strain tensile test (a) outside necking region and (b) fracture zone. work hardening. The phenomenological fracture criterion of Cockcroft and Latham [9] is used to represent ductile fracture. Shear fracture is considered to be a special case of ductile fracture and, at present, no additional fracture criterion is introduced to represent the phenomenon. The yield surface parameters were determined through a hybrid identification procedure based upon available uniaxial tensile test data and selected PCYS data. The isotropic work hardening parameters at low to intermediate strains were identified by a simple curve fit to the uniaxial tensile test data. It is known that the maximum effective strain in forming operations greatly exceeds the uniform strain in uniaxial tensile tests. A minimum work-hardening rate is introduced to achieve a more realistic representation of the work hardening at large plastic strains. This parameter was estimated through numerical analyses of the inplane shear test. Other authors report the use of the bulge-test data to determine the work-hardening curve at large effective strains [3]. The numerical model of the shear test was also used to identify the ductile fracture parameter, taken as the maximum value of the integral in Eq. (9) at a displacement correlating with the maximum force in the experimental test. Ideally, more work should have been performed to rule out influence of geometrical instability events in the performed shear test. This would require in-situ strain field measurements, preferably also supported by brick element analyses. At present, interested readers are referred to Tarigopula et al. [45], where such investigations have been undertaken for a dual-phase steel. Finally, to avoid unphysical strain localization and fracture in shell models a non-local value of the plastic thickness strain rate e_ p33 is used to regularize thinning instabilities. A number of different tests were performed to validate the model concepts and parameter identification procedure. The plane-strain tension tests produce membrane stretching under a biaxial stress state. The elongation is expected to be limited by plastic failure and to correlate with the minimum critical strain in the FLD, i.e. under plane-strain condition. The correlation was found to be excellent both in terms of the predicted force–displacement characteristics and ductility prediction. The radius of influence for the non-local calculation of plastic thinning had some influence on the post-necking response of the specimen. Fig. 20a shows the necking and fracture region of one of the test specimens. The picture bears evidence of a local neck with an extent close to the sheet thickness. The actual fracture appearing within this local neck seems to be caused by shear fracturing, but the limiting phenomenon for the specimen elongation is plastic failure. Plane-strain bending: Pure bending of a sheet does not cause any thinning and suppresses plastic failure. The specimen elongation is thus expected to be limited by ductile fracture of the material. Due to the mutual constraint between the tensile and compressive side within such a specimen, the stress state correlates to plane-strain conditions. The present experiments do not produce pure bending, but the specimen elongation is still thought to be limited by ductile fracture. The analyses and experiments correlated well at moderate deformations and in the sense that the 0° and 90° specimens showed the same response. The ductility of the specimens was underestimated. In the bulge test, an equibiaxial stress state and pre-dominant membrane stretching are experienced. The punch depth at fracture was predicted with very good accuracy, while the maximum pressure was somewhat underestimated. The limiting phenomenon, as predicted by the numerical model, is ductile fracture since no tendency to strain localization was seen in the model prior to satisfaction of the ductile fracture criterion. In the square-cup drawing test, a mixture of stress states and deformation conditions arise and evolve during the deformation process. The model predicted the critical blank diameter with good accuracy. The predicted thinning profiles and force-deformation characteristic for the largest blank diameter that did not fracture were also found to be in good accordance with the experiments. Two different formability test set-ups were explored to determine the FLC of the material. The M–K test methodology is considered to be most consistent with the FE-based formability predictions. The correlation between the predictions and these test data were found to be rather good, though non-conservative under equibiaxial stretching. As with the bulge tests, material fracture phenomena limit the specimen ductility under equibiaxial straining. In view of its simplicity, the modelling procedure and established parameters give a remarkably accurate representation of the 6016-T4 material. In uniaxial tension, plane-strain tension and several of the formability tests, plastic failure is seen to be the critical phenomenon causing strain localization and subsequent fracture. Under equibiaxial stretching, plane-strain bending, and in the critical zones of the square-cup drawing tests, the model indicates that ductile fracture is the limiting phenomenon. The FE models analyzed in the present work rely upon dense meshes, but somewhat coarser meshes could probably also give acceptable results. The modelling procedure is thus considered applicable for direct shell-based prediction of failure, and relevant for industrial process and product development. 3018 O.-G. Lademo et al. / Materials and Design 30 (2009) 3005–3019 Fig. 22. Close-up pictures of equibiaxial M–K formability test (a) outside necking region and (b) fracture zone. A pertinent question is: Does the model represent the actual fracture mechanism of the material? When inspecting the fracture region of one of the equibiaxially stretched M–K test specimens, see Fig. 20b, a local neck is also seen for these specimens. This strain localization is believed to occur during the fracture process of the specimen, i.e. not being the cause of the failure. Fig. 20 shows that the fracture process of the 6016-T4 material involves strong shear localization, i.e. according to Teirlinck et al. [46] it could be classified as a shear fracture phenomenon or void sheeting. This is probably the reason for the scatter in the experimental formability data for these test specimens. The utilized M–K formability test provokes fracture in a wide planar region subject to pure membrane stretching. Such conditions, with no through-thickness gradient in the damage evolution and stress state, are ‘ideal’ for a through-thickness shear instability/ fracture to arise. In the bulge test and the Nakazima formability tests, superimposed bending action probably stabilises this instability. Sarkar et al. [42] studied the ductility and bendability of another prominent age-hardenable Al–Mg–Si–Cu alloy, AA6111, as a function of iron content, ranging from 0.06 to 0.68 wt%. As for the present material, they found that both the low and high iron-containing alloys failed in a shear mode due to void sheeting. Figs. 21 and 22 show close-up pictures of the plane-strain tensile tests and equibiaxial formability test. The results so far indicate that the ductility of the former specimen is limited by plastic failure and the latter by shear localization connected to ductile damage. Part (a) of the figures shows close-up pictures outside the necking region, while part (b) of the figures shows the fracture zones. The equibiaxial formability tests have been deformed to higher effective strains under higher stress triaxiality than the plane-strain tensile tests. The ductile damage is expected to be larger in the former test specimen than in the latter. This is also supported by the pictures, unveiling a more deformed grain structure and more voids adjacent to the particles in the equibiaxial formability tests compared to the plane-strain tensile tests. It can safely be assumed that the material outside the local fracture zones of these membrane-stretched plates unloads elastically when either plastic failure or ductile fracture is attained. Hence, the damage exhibited in Fig. 22a should correspond to the critical damage for shear localization to occur. Even though the shear fracture phenomenon has received much attention in the last decades [24,7,46], it seems like full insight into this phenomenon requires multiscale FE modelling approaches. The present phenomenological and shell-based analyses give lim- ited insight into the fracture process and for this purpose more sophisticated modelling procedures are required. For instance, Madej et al. [34] explored the interplay between the development of microshear bands (micro level), shear bands (meso level), and shear localization (macro level) by means of a ‘Cellular Automata Finite Element’ model. They claim that the results achieved from this model qualitatively reproduced a realistic material behaviour and revealed limitations in conventional FE approaches. A more detailed investigation of the shear failure mechanisms of the AA6016 material is currently undertaken. Acknowledgement The present work has received funding from the COMPFORM project in the Norlight programme, supported by the Norwegian Research Council and Hydro Aluminium, and the Centre for Research-Based Innovation (CRI) SIMLab, hosted by Department of Structural Engineering, Norwegian University of Science and Technology. Helpful discussions with Dr. H. Aretz are gratefully acknowledged. References [1] Aretz H. A non-quadratic plane stress yield function for orthotropic sheet metals. J Mater Process Technol 2004;168:1–9. [2] Barlat F. Crystallographic texture, anisotropic yield surfaces and forming limits of sheet metals. Mater Sci Eng 1987;91:55–72. [3] Barlat F, Glazov MV, Brem JC, Lege DJ. A simple model for dislocation behavior, strain and strain rate hardening evolution in deforming aluminium alloys. Int J Plast 2002;18:919–39. [4] Behr H. 2000. Biegeversuch für Al-Blechwerkstoffe, Prüfanweisung DaimlerChrysler VWM, Sindelfingen, draft from 12.09.2000. [5] Belytschko T, Liu WK, Moran M. Nonlinear finite elements for continua and structures. Wiley; 2000. [6] Berstad T, Lademo O-G, Hopperstad OS. Weak and Strong Texture Models in LS-DYNA (WTM-2D/3D and STM-2D). SINTEF report STF80MK F05180 (proprietary), 2005. [7] Bressan JD, Williams JA. The use of a shear instability criterion to predict local necking in sheet metal deformation. Int J Mech Sci 1983;25:155–68. [8] Brunet M, Morestin F, Walter-Leberre H. Failure analysis of anisotropic sheetmetals using a non-local plastic damage model. J Mater Process Technol 2005;170:457–70. [9] Cockcroft MG, Latham DJ. Ductility and the workability of metals. J Inst Metals 1968. [10] Engler O, Brünger E. Microstructure and texture of aluminium alloys for autobody applications. Matériaux et Techniques 2002;90(5–6):71–8. [11] Engler O, Hirsch J. Texture control by thermomechanical processing of AA6xxx Al–Mg–Si sheet alloys for automotive applications – a review. Mater Sci Eng 2002;A336:249–62. [12] Fyllingen Ø, Hopperstad OS, Lademo OG, Langseth M. Estimation of forming limit diagrams by the use of the finite element method and Monte Carlo simulation. Comput Struct 2009;87:128–39. O.-G. Lademo et al. / Materials and Design 30 (2009) 3005–3019 [13] Gänser H-P, Werner EA, Fischer FD. Forming limit diagrams: a micromechanical approach. Int J Mech Sci 2000;42:2041–54. [14] Gese H, Keller S, Dell H. Verbesserte Plastizitäts- und Versagensmodelle in der Umformsimulation. Tagung Numerische Simulation Verarbeitungsprozesse und prozessgerechte Bauteilgestaltung, Bareuth, 2004. [15] Graf A, Hosford W. Effect of changing strain paths on forming limit diagrams of Al 2008-T4. Met Trans A 1993;24:2503–12. [16] Graf A, Hosford W. Influence of strain-path changes on forming limit diagrams of Al 6111 T4. Int J Mech Sci 1994;36(10):897–8910. [17] Gurson AL. Continuum theory of ductile rupture by void nucleation and growth: Part 1 – Yield criteria and flow rules for porous ductile media. J Eng Mater Technol 1977;99:2–15. [18] Hirsch J. Aluminium alloys for automotive applications. Mater Sci Forum 1997;242:33–50. [19] Hirth SM, Marshall GJ, Court SA, Lloyd DJ. Effects of Si on the aging behaviour and formability of aluminium alloys based on AA6016. Mater Sci Eng A 2001;319–321:452–6. [20] Hooputra H, Gese H, Dell H, Werner H. A comprehensive failure model for crashworthiness simulation of aluminium extrusions. Int J Crashworthiness 2004;9(5):449–63. [21] Hopperstad OS, Berstad T, Ilstad H, Lademo OG, Langseth M. Effects of the yield criterion on local deformations in numerical simulation of profile forming. J Mater Process Technol 1998;80–81:551–5. [22] Hopperstad OS, Børvik T, Berstad T, Lademo O-G, Benallal A. A numerical study on the influence of the Portevin–Le Chatelier effect on necking in an aluminium alloy. Modelling Simul Mater Sci Eng 2007;15:747–72. [23] Hosford WF, Caddel RM. Metal Forming. Mechanics and Metallurgy. PrenticeHall; 1993. [24] Hutchinson JW, Hatherly M, Malin AS, Anand L, Asaro RJ, Needleman A, et al. VIEWPOINT SET NO. 6. Scripta Met. 1984;18:421-458. [25] Lademo OG, Hopperstad OS, Berstad T, Langseth M. Prediction of plastic instability in extruded aluminium alloys using shell analysis and a coupled model of elasto-plasticity and damage. J Mater Process Technol 2004;166:247–55. [26] Lademo O-G, Hopperstad OS, Berstad T, Langseth M. Prediction of plastic instability and fracture in extruded aluminium alloys in shell element analyses. In: Presented at the ESAFORM 2005 conference, Cluj-Napoca, Romania, April 2005. p. 27–9. [27] Lademo O-G, Pedersen KO, Berstad T, Furu T, Hopperstad OS. An experimental and numerical study on the formability of textured AlZnMg alloys. Eur J Mech A/Solids 2008;27:116–40. [28] Lebensohn RA, Tomé CN. A self-consistent anisotropic approach for the simulation of plastic deformation and texture development of polycrystals: application to zirconium alloys. Acta Metall Mater 1993;41:2611–24. [29] Lemaitre J. A continuous damage mechanics model for ductile fracture. J Eng Mater Technol 1985;107:83–9. [30] Lemaitre J, Chaboche J-L. Mechanics of solid materials. Cambridge University Press: Cambridge; 1990. 3019 [31] Liebertz H, Duwel A, Illig R, Hotz W, Keller S, Koehler A, et al. Guideline for the determination of forming limit curves. Sindelfingen, Germany: IDDRG – International Deep Drawing Research Group; 2004. [32] Livermore Software Technology Corporation. LS-DYNA User’s Manuals, Version 971, 2006. [33] Logan RW, Hosford WF. Upper bound anisotropic yield locus calculations assuming h1 1 1i-pencil glide. Int J Mech Sci 1980;22:419. [34] Madej L, Hodgson PD, Pietrzyk M. The validation of a multiscale rheological model of discontinuous phenomena during metal rolling. Comput Mater Sci 2007;41:236–41. [35] Marciniak Z, Kuczynski K. Limit strains in the processes of stretch-forming sheet metal. Int J Mech Sci 1967;9:609–20. [36] Nakajima K, Kikuma T, Hasuka K. Study on the formability of steel sheets. Yawata Technical Report No. 264, 8517-8530, 1968. [37] Oyane M, Sato T, Okimoto K, Shima S. Criteria for ductile fracture and their applications. J Mech Work Technol 1980;4:65–81. [38] Pedersen KO, Lademo O-G, Berstad T, Furu T, Hopperstad OS. On the influence of texture and grain structure on strain localization and formability of AlMgSi alloys. J Mater Proc Technol 2008;200:77–93. [39] Pijaudier-Cabot G, Bazant ZP. Nonlocal damage theory. J Eng Mech ASCE 1987;113:1512–33. [40] Randle V, Engler O. Introduction to texture analysis: macrotexture, microtexture and orientation mapping. Amsterdam: Gordon and Breach Science Publ.; 2000. [41] Reyes A, Hopperstad OS, Lademo O-G, Langseth M. Modeling of textured aluminum alloys used in a bumper system: material tests and characterization. Comput Mater Sci 2006;37(3):246–68. [42] Sarkar J, Kutty TRG, Wilkinson DS, Embury JD, Lloyd DJ. Tensile properties and bendability of T4 treated AA6111 aluminum alloys. Mat Sci Eng A 2004;369(1– 2):258–66. [43] Sowerby R, Duncan JL. Failure in sheet metal in biaxial tension. Int J Mech Sci 1971;13:217–29. [44] Spencer K, Corbin SF, Lloyd DJ. The influence of iron content on the plane strain fracture behaviour of AA 5754 Al–Mg sheet alloys. Mater Sci Eng 2002;A325:394–404. [45] Tarigopula V, Hopperstad OS, Langseth M, Clausen AH, Hild F, Lademo O-G, et al. A study of large plastic deformations in dual phase steel using digital image correlation and FE analysis. Exp Mech 2008;48:181–96. [46] Teirlinck D, Zok F, Embury JD, Ashby MF. Fracture mechanism maps in stress space. Acta Metall 1988;36(5):1213–28. [47] Wang T, Hopperstad OS, Lademo O-G. Larsen PK. Finite element analysis of welded beam-to-column joints in aluminium alloy EN AW 6082 T6. Finite Elements in Analysis and Design, 441-16, 2007. [48] Yoshida K, Kuwabara T. Effect of strain hardening behavior on forming limit stresses of steel tube subjected to nonproportional loading paths. Int J Plast 2007;23:1260–84. [49] Yoshida K, Kuwabara T, Kuroda M. Path-dependence of the forming limit stresses in a sheet metal. Int J Plast 2007;23:361–84.