Forums
Finite Element Mode...
User Material Subro...
Elastic and Plastic...
 
Notifications
Clear all

Elastic and Plastic deformation gradients in AB Viscoplastic model

 
Page 2 / 2
User Material Subroutines
Last Post by Leslie 3 months ago
15 Posts
3 Users
3 Likes
1,446 Views
RSS
Posts: 3990
 Jorgen
Nov 22, 2019 9:59 am
(@jorgen)
Member
Joined: 4 years ago

That is very impressive! ? 

Writing a functioning user-material model is not easy.

-Jorgen

Reply
Melly reacted
1 Reply
 Melly
(@melly)
Joined: 4 years ago

Eminent Member
Posts: 30
Nov 22, 2019 10:43 am
Reply toJorgen

@jorgen

Thank you for your compliment. 

All this is possible courtesy of your book and your answers to my questions. 

With thanks, 

Melly

Reply
Posts: 2
 Leslie
Jul 11, 2023 10:03 am
(@leslie)
New Member
Joined: 3 months ago

Hello @jorgen ,

i am new on this forum and have been reading a few of your interactions with other users. I have been trying to write a VUMAT subroutine based on a modified chaboche model proposed by F. Hu (see document attached). While testing my subroutine on a single element, it works but once i applied it to a normal tensile test specimen, abaqus displays the error "The elements contained in element set ErrElemExcessDistortion-Step1 have distorted excessively". I have reviewed my subroutine countless times and do not see an error. Please can you help me or any other person reading this?

Below is my VUMAT subroutine (only for the plateau region) :

      subroutine vumat(
C Read only -
     1  nblock, ndir, nshr, nstatev, nfieldv, nprops, jInfoArray,
     2  stepTime, totalTime, dtArray, cmname, coordMp, charLength,
     3  props, density, strainInc, relSpinInc,
     4  tempOld, stretchOld, defgradOld, fieldOld,
     3  stressOld, stateOld, enerInternOld, enerInelasOld,
     6  tempNew, stretchNew, defgradNew, fieldNew,
C Write only -
     7  stressNew, stateNew, enerInternNew, enerInelasNew )
C
      include 'vaba_param.inc'
C
C The state variables are stored as:
C      state(*,1) = elastic strain 11
C      state(*,2) = elastic strain 22
C      state(*,3) = elastic strain 33
C      state(*,4) = elastic strain 12
C      state(*,5) = elastic strain 13
C      state(*,6) = elastic strain 23
C      state(*,7) = plastic strain 11
C      state(*,8) = plastic strain 22
C      state(*,9) = plastic strain 33
C      state(*,10) = plastic strain 12
C      state(*,11) = plastic strain 13
C      state(*,12) = plastic strain 23
C      state(*,13) = equivalent plastic strain
C      state(*,14) = radius of the memory surface
C      state(*,15) = center of the memory surface 11
C      state(*,16) = center of the memory surface 22
C      state(*,17) = center of the memory surface 33
C      state(*,18) = center of the memory surface 12
C      state(*,19) = center of the memory surface 13
C      state(*,20) = center of the memory surface 23
C      state(*,21) = sback stress component 11  ! short-range backstress
C      state(*,22) = sback stress component 22
C      state(*,23) = sback stress component 33
C      state(*,24) = sback stress component 12
C      state(*,25) = sback stress component 13
C      state(*,26) = sback stress component 23
C      state(*,27) = sback stress component 11  ! 2nd short-range backstress
C      state(*,28) = sback stress component 22
C      state(*,29) = sback stress component 33
C      state(*,30) = sback stress component 12
C      state(*,31) = sback stress component 13
C      state(*,32) = sback stress component 23
C      state(*,33) = lback stress component 11  ! long-range backstress
C      state(*,34) = lback stress component 22
C      state(*,35) = lback stress component 33
C      state(*,36) = lback stress component 12
C      state(*,37) = lback stress component 13
C      state(*,38) = lback stress component 23
C      state(*,39) = lback stress component 11  ! 2nd long-range backstress
C      state(*,40) = lback stress component 22
C      state(*,41) = lback stress component 33
C      state(*,42) = lback stress component 12
C      state(*,43) = lback stress component 13
C      state(*,44) = lback stress component 23
C      state(*,45) = yield index  ! predefined variable
C
       dimension props(nprops), density(nblock),
     1  coordMp(nblock,*),
     2  charLength(*), dtArray(*), strainInc(nblock,ndir+nshr),
     3  relSpinInc(*), tempOld(*),
     4  stretchOld(*), defgradOld(*),
     5  fieldOld(*),  stressOld(nblock,ndir+nshr),
     6  stateOld(nblock,nstatev), enerInternOld(nblock),
     7  enerInelasOld(nblock), tempNew(*),
     8  stretchNew(*), defgradNew(*), fieldNew(*),
     9  stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
     2  enerInternNew(nblock), enerInelasNew(nblock), jInfoArray(*)
C
      character*80 cmname
  pointer (ptrjElemNum, jElemNum)
       dimension jElemNum(nblock)
C
      parameter ( zero = 0.d0, one = 1.d0, two = 2.d0, three = 3.d0,
     1  third = one/three, half = .5d0, twoThirds = two/three,
     2  threeHalfs = three/two, six = 6.d0)
C
C Elastic constants
 
      e = props(1)          ! E-Mod
      xnu = props(2)        ! Poisson's ratio
      yield0 = props(3)     ! initial yield stress
      epst = props(4)       ! plastic strain at the end of the yield plateau in monotonic curve
      epst_amp = props(5)   ! threshold amplitude of plastic strain
      cskin1 = props(6)     ! short-range kinematic hardening parameter, c_s
      cskin2 = props(7)
      clkin1 = props(8)      ! long-range kinematic hardening parameter, c_l
      clkin2 = props(9)
      sgamma1 = props(10)    ! short-range kinematic hardening parameter, gamma_s
      sgamma2 = props(11)
      lgamma1 = props(12)    ! long-range kinematic hardening parameter, gamma_l
      lgamma2 = props(13)
      qishs = props(14)    ! short-range saturated value of the isotropic hardening component, Q
      qishl = props(15)    ! long-range
      bsats = props(16)    ! short-range rate in which the saturation is achieved, b_s
      bsatl = props(17)    ! long-range rate, b_l
      csp = props(18)    ! short-range scalar parameter, c_s
      clp = props(19)    ! long-range, c_l
 
C
      xmu = e / ( two * ( one + xnu ) )  ! mu (lamé constant) = G (shear modulus)
      twomu = e / ( one + xnu )  ! 2G = 2mu
  sixmu = three * twomu
      thremu =  three * e / ( two * ( one + xnu ) )  ! 3G = 3mu
      alamda = twomu * ( e - twomu ) / ( sixmu - two * e )  ! lamé constant
C      xkap = e / ( three * ( one - two*xnu ) )
      con = xnu / e
      coef1 = sqrt(twoThirds)
      coef2 = sqrt(threeHalfs)
  
  do i = 1, nblock
C Deviatoric strain increment
 
      trace = strainInc(i,1) + strainInc(i,2) + strainInc(i,3)
      de1 = strainInc(i,1) - (third*trace)
      de2 = strainInc(i,2) - (third*trace)
      de3 = strainInc(i,3) - (third*trace)
      de4 = strainInc(i,4)
      de5 = strainInc(i,5)
  de6 = strainInc(i,6)
  
C Trial elastic deviatoric stress increment
 
      s1_inc = twomu * de1
      s2_inc = twomu * de2
      s3_inc = twomu * de3
      s4_inc = twomu * de4
      s5_inc = twomu * de5
      s6_inc = twomu * de6
 
C Predictor stress
 
      sig1 = stressOld(i,1) + alamda*trace + twomu*strainInc(i,1)
      sig2 = stressOld(i,2) + alamda*trace + twomu*strainInc(i,2)
      sig3 = stressOld(i,3) + alamda*trace + twomu*strainInc(i,3)
      sig4 = stressOld(i,4) + twomu*strainInc(i,4)
      sig5 = stressOld(i,5) + twomu*strainInc(i,5)
      sig6 = stressOld(i,6) + twomu*strainInc(i,6)
 
C Trace of stress tensor
 
      smean = third*(sig1 + sig2 + sig3)
      
C Deviatoric stress from the backstress
      s1 = sig1 - smean - stateOld(i,21) - stateOld(i,27)
      s2 = sig2 - smean - stateOld(i,22) - stateOld(i,28)
      s3 = sig3 - smean - stateOld(i,23) - stateOld(i,29)
      s4 = sig4 - stateOld(i,24) - stateOld(i,30)
      s5 = sig5 - stateOld(i,25) - stateOld(i,31)
      s6 = sig6 - stateOld(i,26) - stateOld(i,32)
  
C Von Mises equivalent stress
 
      smises = s1**2 + s2**2 + s3**2 + two*(s4**2 + s5**2 + s6**2)
      smises = sqrt(threeHalfs*smises)
 
C Equivalent plastic strain and radius of memory surface
 
C      eqps0 = stateOld(i,13)
C      rad_ms0 = stateOld(i,14)
 
C Trial elastic and plastic strain
 
      eltr1 = stateOld(i,1) + strainInc(i,1)   
  eltr2 = stateOld(i,2) + strainInc(i,2)
  eltr3 = stateOld(i,3) + strainInc(i,3)
  eltr4 = stateOld(i,4) + strainInc(i,4)
  eltr5 = stateOld(i,5) + strainInc(i,5)
  eltr6 = stateOld(i,6) + strainInc(i,6)
  
  pls1 = stateOld(i,7)
  pls2 = stateOld(i,8)
  pls3 = stateOld(i,9)
  pls4 = stateOld(i,10)
  pls5 = stateOld(i,11)
  pls6 = stateOld(i,12)
  
C Radius of yield surface      
C   radpr0 = yield0 + qishs*( one - exp(-bsats*eqps0) )
C      dradpr = bsats*(yield0 + qishs - radpr0)*eqps0
 
      radpr0 = yield0 + qishs*( one - exp(-bsats*stateOld(i,13)) )
      dradpr = bsats*(yield0 + qishs - radpr0)*stateOld(i,13)   
  
  yindex = zero     ! yield index
C   stateOld(i,45) = yindex
  
      if ( (smises - yield0) .ge. zero ) then   
  
  yindex = one
C   stateNew(i,45) = yindex
  smises = smises + (one - yindex)
  
C On the bounding surface
 
      delta = yindex*(smises - yield0)/thremu    ! equivalent plastic strain increment
 
      flow1 = s1/(yield0 + (thremu*delta))       ! n_i+1
      flow2 = s2/(yield0 + (thremu*delta))   
      flow3 = s3/(yield0 + (thremu*delta))
  flow4 = s4/(yield0 + (thremu*delta))
  flow5 = s5/(yield0 + (thremu*delta))
  flow6 = s6/(yield0 + (thremu*delta))
 
      ps_inc1 = threeHalfs*flow1*delta      ! updated plastic strain increment
  ps_inc2 = threeHalfs*flow2*delta
      ps_inc3 = threeHalfs*flow3*delta
  ps_inc4 = threeHalfs*flow4*delta
  ps_inc5 = threeHalfs*flow5*delta
  ps_inc6 = threeHalfs*flow6*delta
  
  stateNew(i,7) = stateOld(i,7) + ps_inc1    ! updated plastic strain
      stateNew(i,8) = stateOld(i,8) + ps_inc2
      stateNew(i,9) = stateOld(i,9) + ps_inc3 
  stateNew(i,10) = stateOld(i,10) + ps_inc4
  stateNew(i,11) = stateOld(i,11) + ps_inc5
  stateNew(i,12) = stateOld(i,12) + ps_inc6 
 
      ds1 = s1 - (twomu*ps_inc1)       ! update deviatoric stress
      ds2 = s2 - (twomu*ps_inc2)
  ds3 = s3 - (twomu*ps_inc3)
  ds4 = s4 - (twomu*ps_inc4)
  ds5 = s5 - (twomu*ps_inc5)
  ds6 = s6 - (twomu*ps_inc6)
 
      epxi1 = stateNew(i,7) - stateOld(i,15)
      epxi2 = stateNew(i,8) - stateOld(i,16)
      epxi3 = stateNew(i,9) - stateOld(i,17)
      epxi4 = stateNew(i,10) - stateOld(i,18)
      epxi5 = stateNew(i,11) - stateOld(i,19)
      epxi6 = stateNew(i,12) - stateOld(i,20)
 
      term0 = epxi1**2 + epxi2**2 + epxi3**2 + two*(epxi4**2 + epxi5**2     
     1      + epxi6**2)
  glamda = sqrt(twoThirds*term0) - stateOld(i,14)
  
  if (glamda .gt. zero) then
       gdelta = glamda
  else
       gdelta = zero
      end if
 
      radpr1 = (radpr0 + ((yield0 + qishs)*bsats*gdelta ))     ! update - R_i+1 
     1       /(one + (bsats*gdelta))
 
  factor0 = one - (radpr1/yield0)
      alfa1 = factor0*ds1    ! updating the backstress at the bounding surface level
  alfa2 = factor0*ds2
  alfa3 = factor0*ds3
  alfa4 = factor0*ds4
  alfa5 = factor0*ds5
  alfa6 = factor0*ds6
  
C Update the different variables
 
      stateNew(i,13) = stateOld(i,13) + delta    ! update - equivalent plastic strain
      stateNew(i,14) = stateOld(i,14) + (gdelta*csp)       ! update - radius of the memory surface
      
      xms1 = (stateNew(i,7) - stateOld(i,15))/(stateOld(i,14) + gdelta)
      xms2 = (stateNew(i,8) - stateOld(i,16))/(stateOld(i,14) + gdelta)
      xms3 = (stateNew(i,9) - stateOld(i,17))/(stateOld(i,14) + gdelta)
      xms4 = (stateNew(i,10) - stateOld(i,18))/(stateOld(i,14) + gdelta)
      xms5 = (stateNew(i,11) - stateOld(i,19))/(stateOld(i,14) + gdelta)
      xms6 = (stateNew(i,12) - stateOld(i,20))/(stateOld(i,14) + gdelta)
  
  stateNew(i,15) = stateOld(i,15) + (one - csp)*gdelta*xms1      ! update - center of the memory surface
  stateNew(i,16) = stateOld(i,16) + (one - csp)*gdelta*xms2
  stateNew(i,17) = stateOld(i,17) + (one - csp)*gdelta*xms3
  stateNew(i,18) = stateOld(i,18) + (one - csp)*gdelta*xms4
  stateNew(i,19) = stateOld(i,19) + (one - csp)*gdelta*xms5
  stateNew(i,20) = stateOld(i,20) + (one - csp)*gdelta*xms6
  
C Updated stress   
  trace1 = ps_inc1 + ps_inc2 + ps_inc3
C   term1 = third*trace1
C   hyd1 = ps_inc1 - term1
C   hyd2 = ps_inc2 - term1
C   hyd3 = ps_inc3 - term1
  
  stressNew(i,1) = sig1 - ((alamda*trace1) + (twomu*ps_inc1))
  stressNew(i,2) = sig2 - ((alamda*trace1) + (twomu*ps_inc2))
  stressNew(i,3) = sig3 - ((alamda*trace1) + (twomu*ps_inc3))
  stressNew(i,4) = sig4 - twomu*ps_inc4
  stressNew(i,5) = sig5 - twomu*ps_inc5
  stressNew(i,6) = sig6 - twomu*ps_inc6
  
  else
  
      yindex = one
C      stateNew(i,33) = yindex 
  tol = 1.d-7
  delta0 = zero
      delta = yindex*(smises - yield0)/thremu
C      cnt = 0
      do while (abs(delta - delta0) .ge. tol)   
  
  delta0 = delta
  denom1 = one + (sgamma1*delta)
  denom2 = one + (sgamma2*delta)
  
  zt1 = s1 - ( (stateOld(i,21)/denom1) + (stateOld(i,27)/denom2) )        ! eqn(77) - zeta
  zt2 = s2 - ( (stateOld(i,22)/denom1) + (stateOld(i,28)/denom2) )
  zt3 = s3 - ( (stateOld(i,23)/denom1) + (stateOld(i,29)/denom2) )
  zt4 = s4 - ( (stateOld(i,24)/denom1) + (stateOld(i,30)/denom2) )
  zt5 = s5 - ( (stateOld(i,25)/denom1) + (stateOld(i,31)/denom2) )
  zt6 = s6 - ( (stateOld(i,26)/denom1) + (stateOld(i,32)/denom2) )
  
  zta = zt1**2 + zt2**2 + zt3**2 + zt4**2 + zt5**2 + zt6**2
  ztmag = sqrt(zta)
  sigt = sqrt(threeHalfs)*ztmag    ! eqn(78a) - sigma_tilde
  
  flo1 = zt1/sigt       ! eqn(78b) - n_i+1
  flo2 = zt2/sigt
  flo3 = zt3/sigt
  flo4 = zt4/sigt
  flo5 = zt5/sigt
  flo6 = zt6/sigt
  
  gam1 = sgamma1/(denom1**2)
  gam2 = sgamma2/(denom2**2)
  
  fac1 = stateOld(i,21)*gam1 + stateOld(i,27)*gam2       ! 2nd term in eqn(81)
  fac2 = stateOld(i,22)*gam1 + stateOld(i,28)*gam2
  fac3 = stateOld(i,23)*gam1 + stateOld(i,29)*gam2
  fac4 = stateOld(i,24)*gam1 + stateOld(i,30)*gam2
  fac5 = stateOld(i,25)*gam1 + stateOld(i,31)*gam2
  fac6 = stateOld(i,26)*gam1 + stateOld(i,32)*gam2
  
  ffm1 = threeHalfs*(flo1*fac1 + flo2*fac2 + flo3*fac3)
  ffm2 = threeHalfs*(flo4*fac4 + flo5*fac5 + flo6*fac6)
  ffmag = ffm1 + ffm2
  
  sumt = cskin1/(denom1**2) + cskin2/(denom2**2)      ! 2nd SUMMATION in eqn(81)
  sumt2 = (cskin1/denom1) + (cskin2/denom2)       ! 2nd to last term in eqn(79)
  
  dff1 = (fac1 - (ffmag*flo1))/sigt       ! 2nd term in eqn(84)
  dff2 = (fac2 - (ffmag*flo2))/sigt
  dff3 = (fac3 - (ffmag*flo3))/sigt
  dff4 = (fac4 - (ffmag*flo4))/sigt
  dff5 = (fac5 - (ffmag*flo5))/sigt
  dff6 = (fac6 - (ffmag*flo6))/sigt
  
  sflo1 = flo1 + (delta*dff1)     ! 2nd term in eqn(83)
  sflo2 = flo2 + (delta*dff2)
  sflo3 = flo3 + (delta*dff3)
  sflo4 = flo4 + (delta*dff4)
  sflo5 = flo5 + (delta*dff5)
  sflo6 = flo6 + (delta*dff6)
  
  ps_inc1 = (threeHalfs*flo1)*delta      ! eqn(57)
  ps_inc2 = (threeHalfs*flo2)*delta
  ps_inc3 = (threeHalfs*flo3)*delta
  ps_inc4 = (threeHalfs*flo4)*delta
  ps_inc5 = (threeHalfs*flo5)*delta
  ps_inc6 = (threeHalfs*flo6)*delta
  
  stateNew(i,7) = stateOld(i,7) + ps_inc1
      stateNew(i,8) = stateOld(i,8) + ps_inc2
      stateNew(i,9) = stateOld(i,9) + ps_inc3
      stateNew(i,10) = stateOld(i,10) + ps_inc4
      stateNew(i,11) = stateOld(i,11) + ps_inc5
      stateNew(i,12) = stateOld(i,12) + ps_inc6
  
C   ds1 = s1 - (twomu*ps_inc1)
C      ds2 = s2 - (twomu*ps_inc2)
C      ds3 = s3 - (twomu*ps_inc3)
C      ds4 = s4 - (twomu*ps_inc4)
C      ds5 = s5 - (twomu*ps_inc5)
C      ds6 = s6 - (twomu*ps_inc6)
  
  epxi1 = stateOld(i,7) - stateOld(i,15)
      epxi2 = stateOld(i,8) - stateOld(i,16)
      epxi3 = stateOld(i,9) - stateOld(i,17)
      epxi4 = stateOld(i,10) - stateOld(i,18)
      epxi5 = stateOld(i,11) - stateOld(i,19)
      epxi6 = stateOld(i,12) - stateOld(i,20)
  
  sfep1 = sflo1*epxi1 + sflo2*epxi2 + sflo3*epxi3
  sfep2 = sflo4*epxi4 + sflo5*epxi5 + sflo6*epxi6
  factor1 = sfep1 + sfep2
  
  factor2 = two*(epxi4**2 + epxi5**2 + epxi6**2)
  term2 = epxi1**2 + epxi2**2 + epxi3**2 + factor2 
  
  glamda = sqrt(twoThirds*term2) - stateOld(i,14)
  gamadel = (delta + factor1)/(stateOld(i,14) + glamda)
      
  if (glamda .gt. zero) then
   gdelta = gamadel
  else
       gdelta = zero
      end if
  
      term3 = yield0 + qishs - radpr0
  dRddel = (bsats*term3*gdelta)/(one + (bsats*glamda))**2        ! eqn(82)
      term4 = yield0 + qishs
  radpr1 = (radpr0 + (term4*bsats*glamda))/(one + (bsats*glamda))
  
  dhdel = ffmag - thremu - sumt - dRddel
  hdelta = sigt - ((thremu + sumt2)*delta) - radpr1
      
  delta = delta - (hdelta/dhdel)
C   cnt = cnt + 1
  end do
      
C Updating the variables
 
      den1 = one + (sgamma1*delta)
      den2 = one + (sgamma2*delta)
  
      ztn1 = s1 - ( (stateOld(i,21)/den1) + (stateOld(i,27)/den2) )        ! eqn(77) - zeta
  ztn2 = s2 - ( (stateOld(i,22)/den1) + (stateOld(i,28)/den2) )
  ztn3 = s3 - ( (stateOld(i,23)/den1) + (stateOld(i,29)/den2) )
  ztn4 = s4 - ( (stateOld(i,24)/den1) + (stateOld(i,30)/den2) )
  ztn5 = s5 - ( (stateOld(i,25)/den1) + (stateOld(i,31)/den2) )
  ztn6 = s6 - ( (stateOld(i,26)/den1) + (stateOld(i,32)/den2) )
  
      ztm1 = ztn1**2 + ztn2**2 + ztn3**2
      ztm2 = ztn4**2 + ztn5**2 + ztn6**2
      ztnmag = ztm1 + ztm2
      signt = sqrt(threeHalfs*ztnmag)       ! eqn (78b)
  
      xflo1 = ztn1/signt         ! update n_i+1; eqn (78a)
      xflo2 = ztn2/signt
      xflo3 = ztn3/signt
      xflo4 = ztn4/signt
      xflo5 = ztn5/signt
      xflo6 = ztn6/signt
  
      ps_inc1 = (threeHalfs*delta)*xflo1      ! updating the plastic strain increment
      ps_inc2 = (threeHalfs*delta)*xflo2
      ps_inc3 = (threeHalfs*delta)*xflo3
      ps_inc4 = (threeHalfs*delta)*xflo4
      ps_inc5 = (threeHalfs*delta)*xflo5
      ps_inc6 = (threeHalfs*delta)*xflo6
  
      stateNew(i,7) = stateOld(i,7) + ps_inc1        ! updated plastic strain
      stateNew(i,8) = stateOld(i,8) + ps_inc2
      stateNew(i,9) = stateOld(i,9) + ps_inc3
      stateNew(i,10) = stateOld(i,10) + ps_inc4
      stateNew(i,11) = stateOld(i,11) + ps_inc5
      stateNew(i,12) = stateOld(i,12) + ps_inc6
 
C      ds1 = s1 - (twomu*ps_inc1)       ! update deviatoric stress
C      ds2 = s2 - (twomu*ps_inc2)
C   ds3 = s3 - (twomu*ps_inc3)
C   ds4 = s4 - (twomu*ps_inc4)
C   ds5 = s5 - (twomu*ps_inc5)
C   ds6 = s6 - (twomu*ps_inc6)
 
      pxi1 = stateNew(i,7) - stateOld(i,15)
      pxi2 = stateNew(i,8) - stateOld(i,16)
      pxi3 = stateNew(i,9) - stateOld(i,17)
      pxi4 = stateNew(i,10) - stateOld(i,18)
      pxi5 = stateNew(i,11) - stateOld(i,19)
      pxi6 = stateNew(i,12) - stateOld(i,20)
 
      term5 = pxi1**2 + pxi2**2 + pxi3**2 + two*(pxi4**2 + pxi5**2     
     1      + pxi6**2)
  glamda = sqrt(twoThirds*term5) - stateOld(i,14)
  
  if (glamda .gt. zero) then
       gdelta = glamda
  else
       gdelta = zero
      end if
 
      radpr1 = (radpr0 + ((yield0 + qishs)*bsats*gdelta))     ! update - R_i+1 
     1       /(one + (bsats*gdelta))
 
 
      stateNew(i,13) = stateOld(i,13) + delta      ! updated equivalent plastic strain
      stateNew(i,14) = stateOld(i,14) + (gdelta*csp)       ! update - radius of the memory surface
      
      rms1 = (stateNew(i,7) - stateOld(i,15))/(stateOld(i,14) + gdelta)
      rms2 = (stateNew(i,8) - stateOld(i,16))/(stateOld(i,14) + gdelta)
      rms3 = (stateNew(i,9) - stateOld(i,17))/(stateOld(i,14) + gdelta)
      rms4 = (stateNew(i,10) - stateOld(i,18))/(stateOld(i,14) + gdelta)
      rms5 = (stateNew(i,11) - stateOld(i,19))/(stateOld(i,14) + gdelta)
      rms6 = (stateNew(i,12) - stateOld(i,20))/(stateOld(i,14) + gdelta)
  
  stateNew(i,15) = stateOld(i,15) + (one - csp)*gdelta*rms1      ! update - center of the memory surface
  stateNew(i,16) = stateOld(i,16) + (one - csp)*gdelta*rms2
  stateNew(i,17) = stateOld(i,17) + (one - csp)*gdelta*rms3
  stateNew(i,18) = stateOld(i,18) + (one - csp)*gdelta*rms4
  stateNew(i,19) = stateOld(i,19) + (one - csp)*gdelta*rms5
  stateNew(i,20) = stateOld(i,20) + (one - csp)*gdelta*rms6
  
C Updated back stresses and stress
  
  stateNew(i,21) = ( stateOld(i,21) + (cskin1*flon1*delta) )/den1     ! 1st short-range backstress
  stateNew(i,22) = ( stateOld(i,22) + (cskin1*flon2*delta) )/den1
  stateNew(i,23) = ( stateOld(i,23) + (cskin1*flon3*delta) )/den1
  stateNew(i,24) = ( stateOld(i,24) + (cskin1*flon4*delta) )/den1
  stateNew(i,25) = ( stateOld(i,25) + (cskin1*flon5*delta) )/den1
  stateNew(i,26) = ( stateOld(i,26) + (cskin1*flon6*delta) )/den1
  
  stateNew(i,27) = ( stateOld(i,27) + (cskin2*flon1*delta) )/den2     ! 2nd short-range backstress
  stateNew(i,28) = ( stateOld(i,28) + (cskin2*flon2*delta) )/den2
  stateNew(i,29) = ( stateOld(i,29) + (cskin2*flon3*delta) )/den2
  stateNew(i,30) = ( stateOld(i,30) + (cskin2*flon4*delta) )/den2
  stateNew(i,31) = ( stateOld(i,31) + (cskin2*flon5*delta) )/den2
  stateNew(i,32) = ( stateOld(i,32) + (cskin2*flon6*delta) )/den2
  
      trace2 = ps_inc1 + ps_inc2 + ps_inc3
C   factor3 = third*trace2
C   diff1 = ps_inc1 - factor3
C   diff2 = ps_inc2 - factor3
C   diff3 = ps_inc3 - factor3
  
  stressNew(i,1) = sig1 - ((alamda*trace2) + (twomu*ps_inc1))
  stressNew(i,2) = sig2 - ((alamda*trace2) + (twomu*ps_inc2))
  stressNew(i,3) = sig3 - ((alamda*trace2) + (twomu*ps_inc3))
  stressNew(i,4) = sig4 - (twomu*ps_inc4)
  stressNew(i,5) = sig5 - (twomu*ps_inc5)
  stressNew(i,6) = sig6 - (twomu*ps_inc6)
  
      end if
  
      end do
 
      return
      end     

 

Reply
Posts: 3990
 Jorgen
Jul 11, 2023 1:02 pm
(@jorgen)
Member
Joined: 4 years ago

Debugging a VUMAT code can be quite time consuming, and unfortunately I don't have time review your code at this time. Perhaps someone else can. 

Best of luck,
Jorgen

Reply
1 Reply
 Leslie
(@leslie)
Joined: 3 months ago

New Member
Posts: 2
Jul 11, 2023 7:25 pm
Reply toJorgen

Hi @jorgen. Thanks for your answer. The problem with the subroutine comes from everything after the "else"-function. The 1st part of the subroutine (if  (smises - yield0 .ge. zero)then ...) works without any problem. For the 2nd part, i'm just following what the author of the paper i attached in the previous message gave us. So as i said, i am really confused.

Reply
Page 2 / 2
Forum Jump:
  Previous Topic
Next Topic  
Share: