Forming limit curve determination of AA6061-T6 aluminum alloy sheet

TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 20, SỐ K2-2017 51  Abstract—The forming limit curve (FLC) is used in sheet metal forming analysis to determine the critical strain or stress values at which the sheet metal is failing when it is under the plastic deformation process, e.g. deep drawing process. In this paper, the FLC of the AA6061-T6 aluminum alloy sheet is predicted by using a micro-mechanistic constitutive model. The proposed constitutive model is implemented via a vectorized user-defined material subroutine (VUMAT) and integrated with finite element code in ABAQUS/Explicit software. The mechanical behavior of AA6061-T6 sheet is determined by the tensile tests. The material parameters of damage model are identified based on semi-experience method. To archive the various strain states, the numerical simulation is conducted for the Nakajima test and then the inverse parabolic fit technique that based on ISO 124004-2:2008 standrad is used to extracted the limit strain values. The numerical results are compared with the those of M- K, Hill and Swift analytical models. Index Terms—forming limit curve, void growth, Nakajima drawing, Dung model. 1 INTRODUCTION ver many years, the aluminum alloy sheets was widely applied in automotive and civil industries because of their outstanding advantages in high strength and light weight. Therefore, it is necessary to accurately describe their forming behaviors at large strains. Manuscript Received on July 13th, 2016. Manuscript Revised December 06th, 2016. This research is funded by Vietnam National University Ho Chi Minh City (VNU-HCM) under grant number C2017-20-05. Nguyen Huu Hao author is with the Engineering Mechanics Department, Ho Chi Minh City University of Technology, VNU-HCM, Vietnam (e-mail: nguyenhuuhao@tdnu.edu.vn). Nguyen Ngoc Trung author is with School of Mechanical Engineering, 585 Purdue Mall, ME3011, Purdue University, West Lafayette, IN 47907, USA. (e-mail: trungnguyen@perdue.edu). Vu Cong Hoa is with the Engineering Mechanics Department, Ho Chi Minh City University of Technology, VNU-HCM, Vietnam (e-mail: vuconghoa@hcmut.edu.vn). The FLC curve is usually predicted by the Marciniak-Kuczynski (M-K) theory model [1] that based on an inconsistency in sheet. Beside the FLC theory prediction, the Nakajima deep drawing model is also applied widely in experiment and numerical simulation to determine the forming limit curve. Accordingly, the Nakajima test is usually conducted for the several specimens to find the various strain paths that presents forming response of material from uniaxial to biaxial stretched loading state. In this method, the limit strains are determined by an inverse parabolic fit [2, 3] or time-dependent technique [2, 4] at or after the onset of necking. The ductile fracture mechanism of metallic materials and their alloys has been proved to be due to the micro-void nucleation, growth and coalescence in matrix material [5, 6]. A cylindrical micro-void growth in rigid-plastic material based ductile fracture criterion was proposed by McClintock [5]. Dung [7] has modified the McClintock model for the ellipsoidal and cylindrical void growth in hardening matrix material under the remoted stress field and has proposed a constitutive model for porous ductile material. Employing a ductile fracture model to predict FLC is widely applied because it is considered as an effective remedy for saving more time than that of the experiment [3, 8]. In this study, we use a Dung’s porous ductile material model [7], conjugated with the Hill’48 quadratic yield function to predict the FLC of AA6061-T6 aluminum alloy. The ductile fracture model is implemented by a vectorized user-defined material subroutine (VUMAT) in ABAQUS/Explicit software package. The seven specimens with various waist width would be used to numerical simulation and then the limit strains were attained by the inverse parabolic fit in accordance with ISO 12004-2:2008 standard. The present results are compared with the those of theory FLC models. Forming limit curve determination of AA6061-T6 aluminum alloy sheet Nguyen Huu Hao, Nguyen Ngoc Trung, and Vu Cong Hoa O 52          SCIENCE & TECHNOLOGY DEVELOPMENT, Vol 20, No.K2- 2017 2 CONSTITUTIVE MODEL  The  sheet  metal  is  usually  showing  an  anisotropy  so  that  the  von  Mises  equivalent  stress  function  in  the  original  Dung’s  model  [7]  is  replaced by the Hill48 quadratic criterion [9].        2 2 22 33 33 11 1/22 2 2 2 11 22 23 31 122 2 2 e F G H L M N           = -  -   -       (1)  Here  ij    , 1,2,3i j =   are  components  of  Cauchy  stress  tensor,  , , , , ,F G H L M N are  anisotropic coefficients.         0 0 90 0 0 0 0 90 45 90 0 1 1 1 1 1 2 3 2 1 2 R R F , G , H , R R R R R R R N , L M R R = = =      = = =    (2)  The  Lankford’s  coefficients  0R ,  45R   and  90R   are determined by uniaxial tensile tests at  0o ,  45o   and  90o  in rolling direction.  Noting  that  for  isotropic  material,  the  Lankford’s  coefficients  0 45 90 1R R R= = = ,  stress  equivalent  Hill’48  becomes  stress  equivalent  von  Mises [10].  The hardening rule of matrix material,   pf f =   (3)  Here  p is  equivalent  plastic  strain  of  matrix  material.  Gurson  [11]  has  been  introduced  a  yield  function  based  on  mechanism  of  void  nucleation,  growth and coalescences in matrix material. Based  on McClintock’s void growth model  [5], Dung [7]  proposed  not  only  a  yield  function  that  similar  to  Gurson-Tvergaard-Needleman  (GTN)  model  [12]  but  also  addition  of  a  explicitly  hardening  parameter n to consider hardening effects of matrix  material under deformation as follow:      2 2* * 1 22 3 1 2 cosh 1 3 ij ije ff n f q q f     -  =  - -       (4)  Where,  the parameters  1q ,  2q   are proposed by  Tvergaard  and  Needleman  [12],  n   is  hardening  exponent  of  matrix  material,  e is  Hill’48  equivalent  stress,  *f   is  function  of  void  volume  fraction (VVF),  ij  is delta Kronecker.    *                                    if      if  c F c c c c F u c f f f f f f f f f f f f f f   = -  -   -  (5)  Here  cf and  Ff  are critical and onset of fracture  void  volume  fraction,  respectively,  1 2/uf q q=   is  ultimate void volume fraction.  The  evolution  of  void  volume  fraction  is  computed as follow:  growth nucleationdf df df=    (6)  Here,  the  void  volume  fraction  growth  of  the  presence voids in matrix material:   *1 pgrowth ij ijdf f d = -   (7)  Here  pijd  is plastic strain rate tensor.  The evolution of nucleated void volume fraction  during matrix material under deformation:  p nucleationdf Ad=   (8)  The number of nucleated voids A  is  a  function  of equivalent plastic strain of matrix material.  2 1 exp 22 p N N NN f A ss      -  = -         (9)  Where,  Nf ,  Ns ,  N are the parameters relative  to the void nucleation during matrix material under  deformation.  3 NUMERICAL IMPLEMENTATION   Based  on  the  numerical  algorithm  of  Aravas  [13],  the  Dung’s  porous  ductile  model  is  implemented  by  a  vectorized  user-defined  subroutine  (VUMAT)  and  conjugated  with  finite  element  code  of  ABAQUS/Explicit  software.  The  implemented procedure for Dung’s model has been  completed by Hao et al. [14].  4 EXPERIMENTAL WORKS  The experimental works adopted in this section  to identify the mechanical behavior of AA6061-T6  aluminum alloy. The specimens to be designed and  tested according to the ASTM-E8 standard [15].   Tensile  tests  were  accomplished  with  a  thin  sheet  that  its  nominal  thickness  of  2.0  mm.  To  identify  Lankford’s  coefficients  (R0,  R45,  R90),  having  least  three  dog-bone  specimens  on  each  direction  of  the  rolling,  transverse  and  45  degrees  to rolling direction have used. The initial  length of  the gage marks is 50 mm for all tests. The geometry  and  dimension of dog-bone  specimen are  given  in  Figure 1.   TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 20, SỐ K2-2017 53  Figure 1 Dog-bone specimen  The  tensile  tests  help  to  obtain  the  mechanical  properties of AA6061-T6 aluminum alloy as shown  in TABLE 1. The difference of engineering stress- strain behavior in three directions of 0o, 45o, 90o to  rolling direction  is presented as Figure 2. The  true  strain-stress  curve  that  used  to  fit  Swift  hardening  rule is given in Figure 3.  Assuming  that  the  isotropic  hardening  rule  obeys  Swift  model  [16],  fitting  true  strain-stress  curve,  the  hardening  parameters  ( K ,  0 ,  n )  is  obtained  as  TABLE  .  The  Lankford’s  coefficients  are calculated by the eq. (10).      02 3 0 0 ln / ln / f f f w w R l w l w    = =   (10)  Where  2  and  3 are the transverse and normal  strains,  respectively.  l0,  lf,  w0,  wf  (0  and  f  indexes  imply  initial  and  final  values)  are  the  gage  length  and width of dog-bone specimen,  0 ,45 ,90o o o = .  Figure 2. Experimental load behavior in various directions  Figure 3. True stress-strain curve  5 PARAMETER CALIBRATION  To  apply  the  porous  plastic  material  model  to  prediction  of  ductile  fracture,  8  parameters   1 2 0, , , , , , ,F C N N nq q f f f s f   must  be  identified.  In general, any identification procedure that used to  identify  all  these  parameters  would  be  still  requirement  of  the  computational  time  cost.  In  addition, for each material type, may be have more  one  set  of  material  parameter  (non-uniqueness  of  the  solution)  [17-19].    A  literature  review  of  material parameter identification for porous ductile  model is necessary in  this work. On that basis,  the  material  parameters  can be  selected  and  calibrated  for Dung’s model.  Two  parameters  1 1.5q =   and 2 2 1 2.25q q= =   that proposed by Tvergaard [20] to correct result of  numerical calculation and original Gurson model.   The  initial  VVF  parameter  f0  is  determined  by  observation  of  micrograph  of  virgin  material  [21,  22]  or  calibration  [23].For  AA6061  aluminum  alloy,  value  of  initial  VVF  is  provided  by  several  researchers  such  as  Agarwal  et  al.  [22]    (f0    =  0.0014),  Xu  et  al.  [21]  (f0  =  0.0025),  Shen  et  al.  [23]  (f0  =  0.0005).  Therefore,  a  suitable  range  for  the  value  of  f0  VVF  of  AA6061-T6  can  be  lie  in  (0.0005-0.0025).  The  parameter  sN  can  be  explained  through  a  little  metrology  significance  of  nucleated  strain  measurements. The distribution of nucleated  strain  values εN  is assumed  to obey a normal distribution  with  a  standard  deviation  sN.  Qualitatively,  a  low  standard  deviation  shows  that  the  values  of  nucleated  strain  εN  tend  to  be  close  to  the  mean  (also  called  the  expected  value)  of  the  data  set,  while  a  high  standard  deviation  indicates  that  the  values of nucleated  strain εN  are  spread out over  a  wider range of values. In this work, a good quality  of  nucleated  strain  measurements  is  assumed  to  obtain so that value of standard deviation sN of 0.05  is selected.  Two  parameters  εN  and  fN  are  usually  used  as  TABLE 1  MECHANICAL PROPERTIES OF AA6061-T6 ALUMINUM  ALLOY  Young’s  modulus   ( E )  Yield  stress   ( 0 )  Poisson’s  ratio   ( )  74.6 GPa  244 MPa  0.314  TABLE 2  MATERIAL PARAMETERS OF AA6061-T6 ALUMINUM  ALLOY  K (MPa)  0   n  Lankford’s  coefficients  R0  R45  R90  489.74  0.02  0.179  0.55  0.52  0.53  54          SCIENCE & TECHNOLOGY DEVELOPMENT, Vol 20, No.K2- 2017 the fitting parameters. In practice, it is difficulty to  recognize exactly the moment at void nucleation so  that  the  value  of  nucleated  strain  εN  is  relatively  selected  based  on  onset  of  material  damage  [19].  Accordingly,  a  comparison  of  force  vs.  displacement  curve  between  experiment  and  finite  element  method  (FEM)  result  of  pure  Hill48  plasticity  theory  (no  damage)  is  performed  to  estimate the value of nucleated strain εN.  The  mesh  size of 0.5  mm  x  0.5  mm at  critical  zone  is  used  to  mesh  for  dog-bone  specimen.  The  displacement controlled load is applied to top edge  of  specimen.    The  element  type  of  3D,  reduced  integration,  8-nodes  (C3D8R)  used  for  dog-bone  specimen.  Figure 4. Graphics of determination of nucleated strain εN  Values of nucleated, critical and fractured VVF  (fN,  fC,  fF)  are  calibrated  by  matching  load- displacement curve of dog-bone specimen between  experiment and FEM.  The  FEM  simulations  are  performed  for  nine  values of fN = 0.01, 0.015, 0.02, 0.025, 0.03, 0.035,  0.04,  0.045,  0.05.  A  best  matched  result  of  load- displacement  curve  between  FEM  and  experiment  is selected to fit the values of fC and fF in next step.  The  evolution  stage  of  VVF  from  fC  to  fF  increased  more  rapid  than  that  of  previous  period  due to the coalescence of micro-voids lead to quick  losing of loading carrying of matrix material. There  are  25  possible  combinations  of  fC  and  fF  from  TABLE  .  However,  because  of  the  constrain  C Ff f  so that have only 24 runs in ABAQUS is  possible to obtain a best combination of (fC, fF) pair  that matches the experimental curve.  Finally, the best fit parameters for predicting of  ductile fracture are given in TABLE .  The displacement – load curve corresponding to  the best  fitted material parameters are presented  in  Figure 5.  Figure 5. The displacement – load curve after calibration  6 FORMING LIMIT CURVE  6.1 Nakajima test  The Nakajima’s type deep drawing is conducted  for  the  seven  specimens  with  waist  width  w =  30,  55,  70,  90,  120,  145  and  the  circular  shape  as  Figure 6a. The  setup of deep drawing  is presented  in Figure 6b. The blank used mesh  type of 3D, 8- nodes,  reduced  integration  (C3D8R)  whereas  the  punch,  holder  and  die  are  assumed  absolute  hard  with 3D analytical rigid type. The initial mesh size  at  analysis  zone  is  1.0  mm  x  1.0  mm.  Three  element  layers  through  the  thickness  of  blank  are  used. The blank holding force  450holdF   kN is used  to avoid any sliding phenomenon and early damage  at the blank holding region. The friction coefficient  between  the  blank  and  punch  surfaces  is  0.03  whereas  the  friction coefficient value of 0.1 on all  remain contact surfaces is adopted.   Figure 6 (a) Blank and (b) deep drawing setup (unit: mm)  TABLE 3  THE VALUES OF CRITICAL AND FRACTURE VVF FOR CALIBRATION  fC  0.015  0.035  0.06  0.08  0.15  fF  0.08  0.15  0.17  0.2  0.25  TABLE 4  BEST FIT VALUES OF MATERIAL PARAMETERS FOR DUNG MODEL  q1  q2  fF  fC  f0  εN  sN  fN  1.5  2.25  0.15  0.035  0.0018  0.09  0.05  0.03  TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 20, SỐ K2-2017 55  After the blank is clamped and the die is fixed,  the blank is stretched by moving the punch in  vertical direction until its fracture occurs.  (a)                         (b)  Figure 7 The extracted path along cross section of W30 and  W120 specimens  The limit strains are then determined based on  cross-section method that its basic concept is the  analysis of the measured strain data along  predefined cross sections at onset necking time.  The detail procedure of this method is given in ISO  12004-2:2008 – part 2 standard. Accordingly, the  principal strain average value of three extracted  paths along cross section of each specimen (Figure  7) are taken to fit an inverse parabola.  The best fit inverse parabola is limited by the fit  boundaries that presented through the inner (purple  dot  square  line)  and  outer  (green  solid  line)  fit  window limits as Figure 8.   The size of inner fit window (L0) is determined  by the highest peaks of the second derivative of the  second  order  parabola  that  regressed  by  three  consecutive points of principal strain data within a  range of 6 mm.   (a)  (b)  Figure 8 The curve fit of the principal strain data and the limit  strain determination . (a) W30 and (b) W120 specimens  The size of the outer limit window should have  at least 5 points and calculated as follows:  L R FW =W =W /2   (11)  Where  WL  left  fit  window  width,  WR  right  fit  window width   F 2 1W 10 1 / =       (12)  With   2 2, 2,1 / 2 BL BR  =     (13)   1 1, 1,1 / 2 BL BR  =    (14)  The  subscripts  “BL”  and  “BR”  are  used  for  ε1  and  ε2  of  the  left  and  right  inner  boundaries,  respectively.  After  determination  of  fit  boundaries,  the  inverse best  fit parabola  is  fitted by all data points  within  fit  window  (WL  and  WR).  The  resulting  value in the crack position  is  the wanted limits for  principal strains ε1 and ε2.  56          SCIENCE & TECHNOLOGY DEVELOPMENT, Vol 20, No.K2- 2017 6.2 M-K model  The  Marciniak-Kuczynski  (M-K)  model  is  probably  the  most  well-known and  widely  used  to  predict  analytical  FLC  curve  [24].  Marciniak  and  Kuczynski  introduced  imperfections  into  sheets  to  describe  necking  condition.  This  theory  based  on  the  material  inhomogeneity  assumption,  i.e.,  there  is  groove  which  is  perpendicular  to  the  axial  of  maximum principal stress on the sheet surface (see  Figure  9).  This  initial  inhomogeneity  grows  continuously  and  eventually  becoming  a  localized  necking.  From  the  Fig.9,  the  zone  (b)  is  groove  zone,  it  is  assumed  the  zone  (a)  is  homogeneous  zone and obey uniform proportional loading states.  The  x,  y,  z  axes  correspond  to  rolling,  transverse  and normal directions of the sheet, whereas 1 and 2  represent the principal stress and strain directions in  the  homogeneous  region.  Meanwhile,  the  set  of  axes aligned to  the groove is represented by n, t, z  axes,  where  t  is  the  longitudinal  one.  In  the  sheet  metal forming process,  the material  is firstly under  plastic  deformation  with  constant  incremental  stretching  until  maximum  force  happen.  The  M-K  model  assumes  the  flow  localization  occurs  in  the  groove  when  a  critical  strain  is  reached  in  the  homogeneous  region.  Then,  the  values  of  strain  increment  in  two  regions  are  compared  with  specific criterion  (e.g., dε1b 10dε1a)  and  finally  the  material major and minor strain limits are obtained  on the forming limit diagrams.  Because of M-K model based on an assumption  of  plane  stress  state  so  that  Hill48  yield  criterion  can be written as follow      1/2 0 902 2 0 1 2 1 2 90 0 0 1 2 1 1 i R R R R R R        =  -       (15)      1/2 0 90 2 0 1 90 0 0 1 2 1 1 1 i R R R R R R        = =  -         (16)  Figure 9 Marciniak-Kuczynski (M-K) model  Because of M-K model based on an assumption  of  plane  stress  state  so  that  Hill48  yield  criterion  can be written as follow      1/2 0 902 2 0 1 2 1 2 90 0 0 1 2 1 1 i R R R R R R        =  -       (17)      1/2 0 90 2 0 1 90 0 0 1 2 1 1 1 i R R R R R R        = =  -        (18)  The behavior of material  can be  represented  in  the form of power law  n m i i iK  =     (19)  Where  n  hardening  exponent,  m  strain  rate  exponent.  The  ratio  of  the  principal  stress  and  strain  are  defined as follows:  2 2 2 1 1 1 , d d         = = =    (20)  The associated flow rule is expressed by  1 1 id d      =   and  2 2 id d      =    (21)  The yield criterion can be rewritten as follows        1 2 90 1 0 90 1 2 0 2 0 90 1 2 3 90 1 0 2 90 01 i i d d R R R R R R d d R R R R              =  - - - - - = =     (22)  Thus, the strain rate can be written as follows:      0 90 0 902 1 90 0 0 90 1 1 R R R Rd d R R R R      - = =  -    (23)  The ratio of strain rate can be calculated  3 /d dt t =   (24)      90 0 90 0 3 1 2 90 0 90 0 0 901 1 R R R R d d d R R R R R R           = - = -  - - -   (25)  Introducing  a  new  parameter  β  and  using  eq.  (22)      90 0 1 90 0 1 1 1 i R Rd d R R      = =  -     (26)  The ratio of initial thickness between (b) and (a)  zones     0 0 0 b mk a t f t =   (27)  Because of thickness strain   3 0ln /t t =   so that  the current thickness of sheet can be calculated as   0 3expt t = . The present thickness ration is  determined as follow   0 3 3 0 expb b b a a a t t t t  = -   (28)  TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 20, SỐ K2-2017 57  or     3 30 expmk b amkf f  = -   (29)  The equilibrium condition requires that the applied  load remains constant between (a) and (b) zone,  therefore  1 1a bF F=   (30)  If the sheet width is a unity then  1 1a a b bt t =   (31)  or  1 1a mk bf =   (32)  From eq. (18)  a ia mk b ibf   =   (33)  From eq. (19)      n nm m a ia ia ia mk b ib ib ibd f d        =     (34)  From eq. (29)         3 30 exp n nm m a ia ia ia b a b ib ib ibmk d f d          = -     (35)  From eq. (23), (30) and (31), the strain relation  between the (a) and (b) zones is given as follow         3 30 exp m n a a ia ia a m n b b a b ib ibmk b d f d                      = -        (36)  In general, the equilibrium equation (36)  can be  solved  numerically  by  using  the  supplementary  equations (18), (23) and (26). Given a stress ratio in  (a) zone (αa) and a finite increment of strain is also  imposed  in  (a)  zone  ( εa  =  0.001).  The  values  of  hardening  exponent  n  =  0.179  and  of  strain  rate  exponent  m  =  0  are  chosen.  Ratio  of  initial  thickness between (b) and (a) zones   0 0.996mkf = .  Then,  the  numerical  computation  is  performed  by  using  a  computational  program,  e.g.  MatLab  language,  to  determine  the  limit  strain  of  each  strain path in the FLC. The limit strains in (a) zone  (ε1a, ε2a) are determined once condition (dε1b/dε1a >  10) is satisfied.  6.3 Hill model  Hill [25] proposed a model to describe the curve  on  the  left  side  of  the  FLC  (ε2  <  0)  based  on  the  local  necking  condition.  Principal  strains  are  calculated as follows,  1 1 n   =    and  2 1 n   =    (37)  Where  2 1/d d  =   denotes  strain  ratio.  According to eq. (37), the FLC calculated based on  Hill’s  model  is  only  dependent  on  the  hardening  coefficient of n = 0.179 and strain ratio    that lie  in range from -0.5 to zero.  6.4 Swift model  Swift  [16]  introduced  a  criterion  for  predicting  FLC  based  on  the  onset  of  diffuse  necking  criterion.       2 1 2 2 1 1 2 2 n         =  -     (38)       2 1 2 2 1 1 2 2 n         =  -    (39)  It is important to remark that, for plane strain (β  = 0) and equibiaxial tension (β = 1). Similar to Hill  model,  given  hardening  coefficient  of  n  =  0.179  and strain ratio    that lie in range from zero to 1.0,  the right side of FLC curve is plotted in Fig. 10.  Figure 10. FLC curve of AA6061-T6 aluminum alloy sheet  Finally,  the  FLC curve  of  AA6061-T6  sheet  is  obtained by the Nakajima deep drawing simulation  using  a  porous  ductile  model  and  the  analytical  model as shown in Figure 0. The results show that  the  FLC  of  three  models  are  consistent  with  each  other at plane strain state and the FLC curve shape  of Dung-Hill48 model is agree with that of Hill and  Swift models. While that the M-K model displays a  big shift large compare to two remaining models.  7 CONCLUSION  In  this  study,  we  present  a  FLC  determination  of  AA6061-T6  aluminum  alloy  sheet.  Material  properties  and  anisotropy  coefficients  were  obtained from tensile test. Applying Dung’s porous  ductile  model  to  determined  FLC  through  numerical  simulation  of  Nakajima  deep  drawing.  The  inverse  parabolic  fit  technique  that  based  on  58          SCIENCE & TECHNOLOGY DEVELOPMENT, Vol 20, No.K2- 2017 ISO 12004-2:2008-part standard is used to achieve  the limit strain values in forming process. Using the  famous  theory  models  of  the  FLC  calculation  by  M-K,  Swift  and  Hill,  the  analytical  FLC  curve  is  proposed.  The  analytical  FLC  curve  shape  of  Hill  and Swift models agrees with that of the numerical  data  whereas  the  predicted  FLC  curve  in  biaxial  loading  states  using  M-K  model  is  fairly  large  deviation from that of remaining models.   REFERENCES  [1]  Z.  Marziniak  and  K.  Kuczynski,  "Limit  strain  in  the  process  of  stretch-forming  sheet  metals,"  Int.  J.  Mech.  Sci, vol. 9, pp. 609-620, 1967.  [2]  S.  Dicecco,  C.  Butcher,  M.  Worswick,  E.  Boettcher,  E.  Chu,  and  C.  Shi,  "Determination  of  forming  limit  diagrams  of  AA6013-T6  aluminum  alloy  sheet  using  a  time and position dependent  localized necking criterion,"  in  IOP  Conference  Series:  Materials  Science  and  Engineering,  2016,  vol.  159,  no.  1,  p.  012009:  IOP  Publishing.  [3]  A. Kami, B. M. Dariani, A. S. Vanini, D. S. Comsa, and  D.  Banabic,  "Numerical  determination  of  the  forming  limit  curves  of  anisotropic  sheet  metals  using  GTN  damage  model,"  Journal  of  Materials  Processing  Technology, vol. 216, pp. 472-483, 2015.  [4]  W.  Hotz,  M.  Merklein,  A.  Kuppert,  H.  Friebe,  and  M.  Klein,  "Time  dependent  FLC  determination  comparison  of  different  algorithms  to  detect  the  onset  of  unstable  necking  before  fracture,"  in  Key  Engineering  Materials,  2013, vol. 549, pp. 397-404: Trans Tech Publ.  [5]  F.  A.  McClintock,  "A  Criterion  for  Ductile  Fracture  by  the Growth of Holes," Journal of Applied Mechanics, vol.  35, no. 2, pp. 363-371, 1968.  [6]  J. R. Rice and D. M. Tracey, "On the ductile enlargement  of  voids  in  triaxial  stress  fields,"  Journal  of  the  Mechanics and Physics of Solids, vol. 17, no. 3, pp. 201- 217, 6// 1969.  [7]  N. L. Dung, "Three Dimensional Void Growth in Plastic  Materials," Mechanics Research Comunications, vol. 19,  no. 3, p. 227, 1992.  [8]  K.  Kolasangiani,  M.  Shariati,  and  K.  Farhangdoost,  "Prediction  of  forming  limit  curves  (FLD,  MSFLD  and  FLSD)  and  necking  time  for  SS304L  sheet  using  finite  element method and ductile fracture criteria," Journal Of  Computational  And  Applied  Research  In  Mechanical  Engineering, vol. 4, 2015.  [9]  R.  Hill,  "A  theory  of  the  yielding  and  plastic  flow  of  anisotropic metals,"  in  Proceedings  of  the  Royal Society  of  London  A:  Mathematical,  Physical  and  Engineering  Sciences,  1948,  vol.  193,  no.  1033,  pp.  281-297:  The  Royal Society.  [10]  R.  v.  Mises,  "Mechanik  der  plastischen  Formänderung  von  Kristallen,"  ZAMM‐Journal  of  Applied  Mathematics  and  Mechanics/Zeitschrift  für  Angewandte  Mathematik  und Mechanik, vol. 8, no. 3, pp. 161-185, 1928.  [11]  A.  L.  Gurson,  "Continuum  theory  of  ductile  rupture  by  void  nucleation  and  growth:  Part  I—Yield  criteria  and  flow  rules  for  porous  ductile  media,"  Journal  of  engineering materials and technology, vol. 99, no. 1, pp.  2-15, 1977.  [12]  V.  Tvergaard  and  A.  Needleman,  "Analysis  of  the  cup- cone  fracture  in  a  round  tensile  bar," Acta Metallurgica,  vol. 32, no. 1, pp. 157-169, 1// 1984.  [13]  N.  Aravas,  "On  the  numerical  integration  of  a  class  of  pressure-dependent  plasticity  models,"  International  Journal  for  Numerical  Methods  in  Engineering,  vol.  24,  no. 7, pp. 1395-1416, 1987.  [14]  H.  H.  Nguyen,  T.  N.  Nguyen,  and  H.  C.  Vu,  "Implementation  and  Application  of  Dung’s  Model  to  Analyze Ductile Fracture of Metallic Material,"  in AETA  2015:  Recent  Advances  in  Electrical  Engineering  and  Related Sciences: Springer, 2016, pp. 903-913.  [15]  E.  ASTM,  "Standard  test  methods  for  tension  testing  of  metallic  materials,"  Annual  book  of  ASTM  standards.  ASTM, 2001.  [16]  H.  W.  Swift,  "Plastic  instability  under  plane  stress,"  Journal  of  the  Mechanics  and  Physics  of  Solids,  vol.  1,  no. 1, pp. 1-18, 10// 1952.  [17]  Z.  Li, B. Bilby,  and  I.  Howard,  "A  study of  the  internal  parameters  of  ductile  damage  theory,"  Fatigue  &  Fracture of Engineering Materials & Structures, vol. 17,  no. 9, pp. 1075-1087, 1994.  [18]  Z.  Zhang,  "A  sensitivity  analysis  of  material  parameters  for  the  Gurson  constitutive  model,"  Fatigue  &  Fracture  of Engineering Materials & Structures, vol. 19, no. 5, pp.  561-570, 1996.  [19]  R. Kiran and K. Khandelwal,  "Gurson model parameters  for  ductile  fracture  simulation  in  ASTM  A992  steels,"  Fatigue  &  Fracture  of  Engineering  Materials  &  Structures, vol. 37, no. 2, pp. 171-183, 2014.  [20]  V.  Tvergaard,  "On  localization  in  ductile  materials  containing  spherical  voids,"  (in  English),  International  Journal  of  Fracture,  vol.  18,  no.  4,  pp.  237-252,  1982/04/01 1982.  [21]  F. Xu, S. Zhao, and X. Han,  "Use of a modified Gurson  model  for  the  failure  behaviour  of  the  clinched  joint  on  Al6061  sheet,"  Fatigue  &  Fracture  of  Engineering  Materials & Structures, vol. 37, no. 3, pp. 335-348, 2014.  [22]  H.  Agarwal,  A.  Gokhale,  S.  Graham,  and  M.  Horstemeyer,  "Void  growth  in  6061-aluminum  alloy  under  triaxial  stress  state,"  Materials  Science  and  Engineering: A, vol. 341, no. 1, pp. 35-42, 2003.  [23]  Y.  Shen,  T.  F.  Morgeneyer,  J.  Garnier,  L.  Allais,  L.  Helfen, and J. Crépin, "Three-dimensional quantitative in  situ study of crack initiation and propagation in AA6061  aluminum alloy sheets via synchrotron laminography and  finite-element  simulations," Acta Materialia, vol. 61, no.  7, pp. 2571-2582, 2013.  [24]  Z.  Marciniak  and  K.  Kuczyński,  "Limit  strains  in  the  processes  of  stretch-forming  sheet  metal,"  International  journal  of  mechanical  sciences,  vol.  9,  no.  9,  pp.  609IN1613-612IN2620, 1967.  [25]  R.  Hill,  "On  discontinuous  plastic  states,  with  special  reference to  localized necking in thin sheets," Journal of  the Mechanics and Physics of Solids, vol. 1, no. 1, pp. 19- 30, 1952/10/01 1952.  TẠP CHÍ PHÁT TRIỂN KH&CN, TẬP 20, SỐ K2-2017 59  Nguyen Huu Hao received  the  B.E.  (2004),  Mechanical  Engineering,  Military  Technical  Academy  (MTA),  Vietnam,  M.E.  (2012)  degree  in  engineering  mechanics  from  Ho  Chi  Minh  City  University  of Technology - VNU-HCM.  He  is  a  Ph.D.  Candidate  at  Department  of  Engineering  Mechanics,  Ho  Chi  Minh  City  University  of  Technology  -  VNU-HCM.  His  current  interests  include  the  modelling  and  predicting of ductile fracture of metallic materials.  Nguyen Ngoc Trung  B.E.  (2001,  Aeronautical  Engineering,  Ho  Chi  Minh  City University of Technology  - VNU-HCM, Vietnam), M.Sc.  (2003,  Mechanics  of  Construction,  Liège  University,  Belgium),  Ph.D.  (2008, Aerospace Engineering, Konkuk University,  South Korea)   He is current working as a postdoctoral research  associate  at  Purdue  University,  IN,  USA.  His  research  interests  include  constitutive  modeling  of  materials,  materials  design,  and  mechanics  of  advanced manufacturing processes.  Vu Cong Hoa received the  B.E. (1995), M.E. (2000,), and  Dr.Eng. (2006) in Mechanical  Design from Chonbuk  National University, South  Korea.  He is an Associate Professor,  at Department of Engineering  Mechanics, Ho Chi Minh City University of  Technology - VNU-HCM. His current interests  include estimating of microscopic ductile fracture  of metallic materials and micro-mechanical  behavior of composite materials.  60          SCIENCE & TECHNOLOGY DEVELOPMENT, Vol 20, No.K2- 2017 Tóm tắt - Đường cong giới hạn gia công được sử dụng trong phân tích gia công kim loại dạng tấm nhằm xác định các giá trị ứng suất hoặc biến dạng tới hạn mà tại các giá trị tới hạn này vật liệu sẽ bị hư hòng khi chịu biến dạng dẻo, ví dụ như quá trình dập kim loại. Bài báo này nhằm dự đoán giới hạn gia công của tấm hợp kim nhôm AA6061-T6 dựa trên mô hình nứt dẻo vi mô. Mô hình cơ sở được lập trình dưới dạng chương trình vật liệu người dùng kết hợp với mã phần tử hữu hạn trong phần mềm ABAQUS/Explicit. Các thí nghiệm kéo đơn trục được thực hiện để xác định ứng xử cơ tính của vật liệu. Các tham số đầu vào của mô hình cơ sở được xác định dựa trên phương pháp bán kinh nghiệm. Để đạt được các trạng thái biến dạng khác nhau, các mẫu dập sâu Nakajima được sử dụng để mô phỏng và kỹ thuật hồi quy parapol ngược theo chuẩn ISO 124004-2:2008 được áp dụng để tính các giá trị biến dạng giới hạn. Các kết quả đạt được thông qua mô phỏng số sẽ được so sánh với các mô hình giải tích như M-K, Hill và Swift. Từ khóa - đường cong giới hạn gia công, tăng trưởng lỗ hổng vi mô, dập Nakajima, mô hình N. L. Dung  Xác định đường cong giới hạn gia công   tấm hợp kim nhôm AA6061-T6  Nguyễn Hữu Hào, Nguyễn Ngọc Trung, Vũ Công Hoà