Prediction of the ductility limit of magnesium AZ31B alloy

Prediction of the ductility limit of magnesium AZ31B alloy Mohamed Yassine Jedidi, Mohamed Ben Bettaieb, Anas Bouguecha, Farid Abed-Meraim, Mohamed Taoufik Khabou, Mohamed Haddar To cite this version: Mohamed Yassine Jedidi, Mohamed Ben Bettaieb, Anas Bouguecha, Farid Abed-Meraim, Mohamed Taoufik Khabou, et al.. Prediction of the ductility limit of magnesium AZ31B alloy. Second International Conference on Advanced Materials, Mechanics and Manufacturing (A3M’2018), Dec 2018, Hammamet, Tunisia. pp.182-193. ฀hal-03048492฀ HAL Id: hal-03048492 https://hal.science/hal-03048492 Submitted on 9 Dec 2020 HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. PREDICTION OF THE DUCTILITY LIMIT OF MAGNESIUM AZ31B ALLOY Mohamed Yassine JEDIDI1,2, Mohamed BEN BETTAIEB2,3, Anas BOUGUECHA1, Farid ABED-MERAIM2,3, Mohamed Taoufik KHABOU1, Mohamed HADDAR1 Laboratoire de Mécanique, Modélisation et de Production (LA2MP), École Nationale d’Ingénieurs de Sfax (ENIS), Route Soukra Km 3.5 1 [email protected], [email protected], [email protected], [email protected] 2 Arts et Métiers ParisTech, Université de Lorraine, CNRS, LEM3, F-57000, France 3 DAMAS, Laboratory of Excellence on Design of Alloy Metals for low-mAss Structures, Université de Lorraine, France [email protected], [email protected] Abstract. In many engineering applications (automotive, computer and mobile device industries, etc.), magnesium alloys have been widely used owing to their interesting physical and mechanical parameters. However, magnesium alloys are identified by the low ductility at room temperature, due to their strong plastic anisotropy and the yielding asymmetry between tension and compression. In this work, the ductility limit of a rolled magnesium AZ31 sheet metal at room temperature is numerically investigated. This investigation is based on the coupling between a reduced-order crystal plasticity model and the Marciniak– Kuczyński localized necking approach. This reduced-order model is used to describe the anisotropic behavior of this material taking into account the strong plastic anisotropy (e.g., yielding asymmetry between tension and compression) due to the limited number of slip systems (i.e., twinning mode). To accurately describe the plastic anisotropy due to slip and twinning modes, a combination of two separate yield functions (according to Barlat and Cazacu) is used. The coupling between the adopted constitutive framework and the Marciniak–Kuczyński instability approach is numerically implemented via an implicit algorithm. Comparisons between experimental results from the literature and numerical results obtained by using our calculation tool are carried out to validate the choice of the reducedorder crystal plasticity model. Keywords: ductility limit; reduced-order crystal plasticity model; Marciniak–Kuczyński instability approach; plastic anisotropy; yielding asymmetry. 1 Introduction The magnesium alloy AZ31B has been used in a wide range of engineering applications, notably in the aeronautic and automotive industries (Yu et al. 2005). Compared to other materials with hexagonal structures, the magnesium alloy AZ31B is characterized by its lightweight and high specific strength. However, it is known to have low formability at room temperature, due to its strong plastic anisotropy and the yielding asymmetry between tension and compression. These characteristics are respectively caused by the limited number of slip systems and the activation of the Mohamed Yassine JEDIDI, Mohamed BEN BETTAIEB, Anas BOUGUECHA, Farid ABED-MERAIM, Mohamed Taoufik KHABOU, Mohamed HADDAR 2 twinning deformation mode. To accurately predict the plastic response of the magnesium alloy AZ31B, several researchers have developed constitutive models that take into account both deformation modes (namely slip and twinning). In this work, we have chosen a reduced-order crystal plasticity model. This model has been recently developed for magnesium alloys (Steglich et al. 2016; Madi et al. 2017; Kondori et al. 2018). It allows describing the strong plastic anisotropy and its evolution as well as the yielding asymmetry between tension and compression, which results from the competition between slip and twin mechanisms. Hence, the plastic anisotropy is described in this reduced-order model by two separate yield functions: the Cazacu yield function (Cazacu et al. 2006) to describe the plastic anisotropy due to the activation of twinning mode and the Barlat yield function (Barlat et al. 1991) to model slipinduced plastic anisotropy. Plastic deformation during sheet metal forming processes is limited by the occurrence of localized necking. To predict the ductility limit, the concept of forming limit diagram (FLD) is introduced and applied in the present investigation. The determination of the ductility limit by experimental methods is often complex and relatively expensive. To avoid these practical difficulties, several theoretical approaches for localized necking analysis have been developed by different researchers to predict FLDs, such as: the Swift model (Swift 1952), the Hill model (Hill 1952), the Marciniak–Kuczyński (M–K) imperfection model (Marciniak and Kuczyński 1967), the bifurcation theory (Rudnicki and Rice 1975), the stability analysis by linear perturbation theory (Dudzinski and Molinari 1991) and the modified maximum force criterion (Aretz 2004). The M–K formability limit analysis is one of the most widely-used ductility approaches in the literature. Several researchers have coupled this approach with various constitutive models to predict forming limit diagrams of HCP materials. In this field, one may quote (Wu et al. 2015) who have coupled this approach with an elaborate version of the Cazacu model (Plunkett et al. 2006) to predict the FLD for magnesium AZ31B alloy. It is well known that FLDs are strongly influenced by the mechanical properties of the used sheet metal, such as plastic anisotropy and hardening parameters. Hence, two hardening laws are applied with the combined yield function (Cazacu–Barlat). In the current study, a generic Mathematica script is developed to couple the M–K approach with the reduced-order crystal plasticity model. The numerical predictions are compared with experimental results of (Kondori et al. 2018) to validate the implementation of this reduced-order model. Then, the ductility limit of the AZ31B magnesium alloy is analyzed. Furthermore, the effect of the activated deformation mode on the ductility limit is also investigated. Notations, conventions and abbreviations Vectors and tensors are indicated by bold letters and symbols. Scalar variables and parameters are designated by thin letters or symbols. The used Void convention leads to express all second-order symmetric tensors in  vector form as: V = 11 ,  22 , 33 , 2 12 , 2  23 , 2 13  T .  PREDICTION OF THE DUCTILITY LIMIT OF MAGNESIUM AZ31B ALLOY 2 3 Reduced-order crystal plasticity model The constitutive model proposed in (Kondori et al. 2018) is adopted here to describe the material behavior of magnesium sheets. This model is formulated within the framework of rate-independent finite plasticity. As this model is used to predict the ductility limit of thin sheets, the plane-stress assumption is adopted (see, e.g., Hutchinson et al. 1978). To describe the anisotropic plasticity due to slip and twinning modes, two independent yield surfaces are introduced. 2.1 Twinning mode Twinning leads to yielding asymmetry between tension and compression. For this reason, the non-quadratic yield criterion developed by (Cazacu et al. 2006) is used. Within this criterion, the effective stress σ for the twinning mode is expressed as:  σ ,  ,  t t 1 t 2 t 3     t 1 k t 1   at t t 2 k t 2   at t 3 k t 3  at  1 at , (1) where 1 ,  2 and  3 are the eigenvalues of tensor  and k (comprised bet t t t tween -1 and 1) and a t (superior to 1) are material parameters describing the initial stress differential effect and the shape of the yield surface, respectively. Tensor  is obtained by a linear transformation from the deviatoric Cauchy stress t tensor s ;   L  s . The transformation tensor L is given in the following matrix form: t t  l LL  t  L LT  LtLS t L   0  0   0  t t t 0 L LT L LS t 0 0 t TT t TS 0 0 lSS t 0 0 0 0 0 0 t LT 0 0 0 0 0 l TS 0 0 0 0 l t L TS L l t  , (2) 0  lSL  t where Lt is expressed onto three orthotropic axes using three directions: longitudinal (L), transverse (T) and short-transverse (S). The activation of the twinning deformation mode is governed by the following Kuhn–Tucker inequality constraints:   σ Y 0 ; t t t ε t  0 ;  tεt  0 , (3) Mohamed Yassine JEDIDI, Mohamed BEN BETTAIEB, Anas BOUGUECHA, Farid ABED-MERAIM, Mohamed Taoufik KHABOU, Mohamed HADDAR  4 where Y is the hardening yield function and ε is the equivalent plastic strain t t t rate due to twinning. Y is expressed as a function of the cumulated plastic strain due to slip ε s and twinning ε t :  Y t  ε s , ε t   R 0t +Q1t 1  e  b1t ε t  +Q  e t qt εt   1  H t ε t  H st ε s . (4) The plastic strain rate ε due to twinning is given by the normality law: t ε =ε t t  t    . (5) The equivalence relationship of the plastic work rate allows us to determine the expression of ε : t ε t 2.2 = ε σ t . t Slip mode (6) Limited number of slip systems for HCP materials leads to a strong initial and evolving plastic anisotropy. For this reason, the yield criterion proposed by (Barlat et al. 1991) and based on a linear transformation technique is used in this work. Within this criterion, the effective stress σ for the slip mode is expressed as: s   a a  1 s s as s s s s s s  σ  ,  ,     1   2   2   3   3  1  2  s s 1 s 2 s 3 1 as , (7) where 1 ,  2 and  3 are the eigenvalues of tensor  and a s (superior to 1) is s s s s a material parameter describing the shape of the yield surface. Tensor  is obs tained by a linear transformation from the Cauchy stress tensor:   L  σ . The s s transformation tensor L is given in the following matrix form: s PREDICTION OF THE DUCTILITY LIMIT OF MAGNESIUM AZ31B ALLOY 1 s 1 s 1 s  s 0 0 0   LSS  L TT  3  L TT  LSS  3 3   1 1 s 1 s s   LsSS LSS  L LL  0 0 0   L LL  3 3 3   s 1 s 1 s . s L    1 Ls L TT  L LL  0 0 0   L LL  TT  3 3 3   s 0 0 0 L TL 0 0   s  0 0 0 0 L TS 0   s  0 0 0 0 0 LSL   5 (8) The activation of the slip deformation mode is governed by the following Kuhn– Tucker inequality constraints:   σ Y 0 ; s s s ε s   0 ;  s εs  0 , (9) where Y is the hardening yield function and ε is the equivalent plastic strain s s s rate due to slip. Y is expressed as a function of the cumulated plastic strain due to slip ε s and twinning ε t :  Y s  ε s , ε t   R s0 +Q1s 1  e  b1s εs  +Q  1  e s 2  bs2 εs H ε . ts t (10) The plastic strain rate ε is given by the normality law: s ε =ε s s  s    . (11) The equivalence relationship of the plastic work rate allows us to determine the expression of ε : s ε = s 2.3 ε σ s s . Combined yield surface (12) The plastic anisotropy of the studied HCP materials is defined by a combination of the Cazacu and Barlat yield functions as shown in Fig. 1. Mohamed Yassine JEDIDI, Mohamed BEN BETTAIEB, Anas BOUGUECHA, Farid ABED-MERAIM, Mohamed Taoufik KHABOU, Mohamed HADDAR 6 Fig. 1. Schematic illustration of the combined yield surface. The plastic strain rate given by Eq. (5) (resp. Eq. (11)) represents only the contribution to the total plastic strain rate ε when the twinning mode (resp. the slip mode) is activated. When both deformation modes are activated simultaneously, the stress state is located at the intersection of the two yield surfaces as illustrated p in Fig. 1, and the plastic strain rate ε is expressed as: p ε =ε p 3 t  t    ε s  s   . (13) M–K localized necking criterion The theoretical M–K approach assumes that the sheet is composed of two zones: the homogeneous zone (a) and the band zone (b). Fig. 2 illustrates this concept, where the band is characterized by its orientation  with respect to the major strain direction, and its relative thickness f with respect to the homogeneous zone ( f  t b / t a ). The relative thickness f is also called the imperfection factor or imperfection ratio, and its initial value classically ranges between 0.98 and 0.999 (Eyckens et al. 2011). PREDICTION OF THE DUCTILITY LIMIT OF MAGNESIUM AZ31B ALLOY 7 Fig. 2. Marciniak–Kuczyński (M–K) model. The current band orientation θ is related to its initial counterpart θ0 by the following relation: tang  θ   tang  θ0  e  ε11a  εa22  , (14) where 11 and a22 are the major and the minor strain in the plane of the sheet (in the homogenous zone). a In addition to the constitutive equations described in Section 2, the M–K approach is defined by the following equations:  The equation describing the force equilibrium between the band and the homogeneous zone: σa .n t a  σ b .n t b . (15)  The geometric compatibility condition linking the velocity gradient in the homogeneous zone La to its counterpart in the band zone Lb : allowing us to relate the strain rate fields in the two zones: Lb  La    n , (16) where  is the jump vector. To predict the occurrence of strain localization, a biaxial loading is applied to the sheet metal, which is defined by the following velocity gradient in the homogeneous zone:  La11 La    0 0  ,  La11  (17) where  is the strain-path ratio. Due to the presence of the geometrical imperfection, the plastic strain will become more and more concentrated within the band (as compared to the homogeneous zone). Strain localization is assumed to occur when the strain-rate ratio ε 33 / ε 33 reaches a critical value (this critical value is b a typically set to 10). When this strain localization condition is satisfied, the current strain component ε 11 is referred to as the critical strain ε 11 corresponding to the strain-path ratio  a * and to the initial band orientation θ0 . The prediction of the entire FLD using the M–K approach is based on the following incremental algorithm with two nested loops:  For each strain-path ratio ranging from ρ  0.5 to ρ  1 (we take intervals of 0.1). Mohamed Yassine JEDIDI, Mohamed BEN BETTAIEB, Anas BOUGUECHA, Farid ABED-MERAIM, Mohamed Taoufik KHABOU, Mohamed HADDAR 8  For each initial band orientation ranging from θ0  0 to θ0  90 (with typical intervals of 1°). Apply the implicit incremental algorithm based on the theoretical framework and constitutive equations presented above to calculate the critical * . strain ε 11 The smallest critical strain ε 11 over all possible initial band orientations θ0 de* fines the localization limit strain ε 11 corresponding to the strain-path ratio  . L This major limit strain is used to calculate the associated minor limit strain ε 22 usL ing the corresponding strain-path ratio  . 4 Results and discussions The numerical simulations in the current work are based on the parameters for the AZ31 magnesium alloy, as shown in Table 1. For elastic behavior, this material is characterized by Young’s modulus E  43 GPa and Poisson’s ratio υ  0.29 . Table 1. Isotropic hardening and plastic anisotropy parameters using slip and twin model for magnesium AZ31B alloy (Kondori et al. 2018). Slip hardening anisotropy R s0  206 MPa Q1s  126 MPa Hts  637 MPa Qs2  40.3 MPa bs2  230 b1s  23.6 LsLL  1 LsTT  1.09 LsSS  1.25 s LT L Twin hardening anisotropy  1.22 s TS L  1.69 s SL L  1.55   R  30 MPa Q  0 MPa b Qt  19.4 MPa q t  6.3 Ht  0 t lLL 1 t lTT  0.97 t lSS  0.33 t 0 l t LT t LT L  0.78  0.22 t 1 l t TS t TS L  0.61  0.76 t 1 t SL l t LS L as  4  0.60  0.73 Hst  386 MPa at  4 k  0.53 PREDICTION OF THE DUCTILITY LIMIT OF MAGNESIUM AZ31B ALLOY 9 In the current paper, the amount of deformation due to slip and twinning modes is predicted by numerical simulations. In the literature, this amount can be experimentally investigated by using the electron backscatter diffraction (EBSD) which is capable to observe the activation of different twin and slip systems. Based on EBSD, twin mode is identified from local intergranular misorientation maps, and active slip system is identified from long range intragranular misorientation maps (Tong et al. 2018). Another diffraction type denoted X ray diffraction (XRD) is used in the literature to calculate the amount of deformation, but for materials with high symmetry. Thus, XRD is not valid for HCP materials due to their yielding asymmetry between tension and compression. In a preliminary stage, some relevant comparisons are carried out in the current work, using the parameters identified in (Kondori et al. 2018), to assess and validate the numerical implementation of the constitutive framework. Then, predictions of forming limit diagrams (FLDs) are undertaken for magnesium AZ3B alloy. Experimental results were carried out based on uniaxial tension and compression tests for the AZ31B magnesium alloy to predict the engineering stress F S0 , where F is the applied force in the tensile and compression tests, and S 0 represents the nominal section of the sheet. Moreover, numerical simulations are carried out, taking into account the anisotropy parameters as well as the hardening function. Despite the complexity of the material behavior, Fig. 3 clearly shows that the reduced-order model is able to reproduce the experimental results with good agreement, owing to flow stresses and plastic anisotropy parameters. Fig. 3. Engineering stress versus axial strain. Accordingly, the constitutive modeling of the AZ31B magnesium alloy is validated compared with experiments in (Kondori et al. 2018). The aim of this investigation is to predict FLDs for this material based on the M–K approach using the reduced-order model. For this reason, the numerical simulation is carried out with a fixed imperfection factor f  0.98 . The onset of localized necking is determined by the numerical necking criterion I n (where I n represents the strain-rate ratio ε 33b / ε 33a ). More specifically, local- Mohamed Yassine JEDIDI, Mohamed BEN BETTAIEB, Anas BOUGUECHA, Farid ABED-MERAIM, Mohamed Taoufik KHABOU, Mohamed HADDAR 10 ized necking occurs when the strain-rate ratio I n becomes larger than 10 for a given strain path and initial band orientation. Fig. 4 shows the evolution of this localized necking indicator for four representative strain paths and an initial band orientation 0  0 : uniaxial tension    0.5 , plane-strain tension    0  , biaxial tension  e.g.;   0.5 , and equibiaxial expansion    1 . Fig. 4. Prediction of localized necking. The imperfection factor depends on thicknesses t a , t b and the initial imperfection factor f 0 . It is always smaller than f 0 during sheet deformation (see Fig. 5). Fig. 5. Imperfection factor evolution. The necking criterion of Section 3 allows us to predict the FLD of the AZ31B magnesium sheet metal. Fig. 6 compares the FLD predicted by the proposed reduced-order Cazacu–Barlat model to that obtained by the simpler Barlat–von Mises model. PREDICTION OF THE DUCTILITY LIMIT OF MAGNESIUM AZ31B ALLOY 11 Fig. 6. Influence of the adopted yield criterion on the predicted FLD. To recover the classical von Mises criterion from the Cazacu criterion, the following parameters should be taken: a t  2 , k   , and Lt  I , where I is the identity tensor. The FLD predicted with the reduced model is limited to the strain between -0.07 and 0.15, whereas the FLD based on von Mises-Barlat model is ranged from -0.10 to 0.20. This is due to the limits of validity of the proposed model. From these predicted FLDs, it is clearly seen that the twinning mode has a clear distinct influence on the forming limit diagram. As expected, the mechanical behavior of the material and its evolution have a strong influence on the predicted FLDs. Due to the twinning mechanism, the FLD predicted by the Barlat–Cazacu model is lower than that predicted by the more classical Barlat–von Mises model. 5 Conclusions A reduced-order crystal plasticity method has been used in this work to model the plastic anisotropy while taking into account the twinning and slip mechanisms. The resulting two-surface constitutive model has been coupled with the M–K imperfection approach to predict the ductility limit of the AZ31B magnesium alloy. The proposed numerical tool has been developed as a generic Mathematica script so that it can be easily applied to other HCP materials. The predicted material behavior has been compared with experimental results, which allowed us to validate the proposed reduced-order model. Also, it has been shown that the combined yield surface (Cazacu–Barlat) has a great influence on the predicted FLD. References Aretz, H., 2004. Numerical restrictions of the modified maximum force criterion for prediction of forming limits in sheet metal forming. Modelling and Simulation in Materials Science and Engineering, 12(4), 677-692. Barlat, F., Lege, D. J., & Brem, J. C. (1991). A six-component yield function for anisotropic materials. International Journal of Plasticity, 7(7), 693-712. Cazacu, O., Plunkett, B., & Barlat, F. (2006). Orthotropic yield criterion for hexagonal closed packed metals. International Journal of Plasticity, 22(7), 1171-1194. Mohamed Yassine JEDIDI, Mohamed BEN BETTAIEB, Anas BOUGUECHA, Farid ABED-MERAIM, Mohamed Taoufik KHABOU, Mohamed HADDAR 12 Dudzinski, D., & Molinari, A. (1991). Perturbation analysis of thermoviscoplastic instabilities in biaxial loading. International Journal of Solids and Structures, 27(5), 601-628. Eyckens, P., Van Bael, A., & Van Houtte, P. (2011). An extended Marciniak–Kuczynski model for anisotropic sheet subjected to monotonic strain paths with through-thickness shear. International Journal of Plasticity, 27(10), 1577-1597. Hill, R. (1952). On discontinuous plastic states, with special reference to localized necking in thin sheets. Journal of the Mechanics and Physics of Solids, 1(1), 19-30. Hutchinson, J. W., Neale, K. W., & Needleman, A. (1978). Sheet necking—I. Validity of plane stress assumptions of the long-wavelength approximation. In Mechanics of sheet metal forming (pp. 111126). Springer, Boston, MA. Kondori, B., Madi, Y., Besson, J., & Benzerga, A. A. (2018). Evolution of the 3D plastic anisotropy of HCP metals: Experiments and modeling. International Journal of Plasticity (in press). Marciniak, Z., & Kuczyński, K. (1967). Limit strains in the processes of stretch-forming sheet metal. International Journal of Mechanical Sciences, 9(9), 609IN1613-612IN2620. Madi, Y., Benzerga, A., & Besson, J. (2017). Modeling the 3D plastic anisotropy of magnesium AZ31B alloy. In Contributions to the Foundations of Multidisciplinary Research in Mechanics, Proceedings of the XXIV International Congress of Theoretical and Applied Mechanics (ICTAM). IUTAM (pp. 2730-2731). Plunkett, B., Lebensohn, R. A., Cazacu, O., & Barlat, F. (2006). Anisotropic yield function of hexagonal materials taking into account texture development and anisotropic hardening. Acta Materialia, 54(16), 4159-4169. Rudnicki, J. W., & Rice, J. R. (1975). Conditions for the localization of deformation in pressuresensitive dilatant materials. Journal of the Mechanics and Physics of Solids, 23(6), 371-394. Swift, H. (1952). Plastic instability under plane stress. Journal of the Mechanics and Physics of Solids, 1(1), 1-18. Steglich, D., Tian, X., & Besson, J. (2016). Mechanism-based modelling of plastic deformation in magnesium alloys. European Journal of Mechanics-A/Solids, 55, 289-303. Tong, V., Wielewski, E., & Britton, B. (2018). Characterisation of slip and twinning in high rate deformed zirconium with electron backscatter diffraction. arXiv preprint arXiv:1803.00236. Yu, K., Li, W. X., & Wang, R. C. (2005). Plastic deformation mechanism of magnesium alloys. Chinese Journal of Nonferrous Metals, 15(7), 1081.