Isotropic Hardening VUMAT Problem

Hi All,

      I have a problem of an Expression of the VUMAT for Isotropic hardening code. part of the code is as follows.

subroutine vumat(

c Read only -

     1  nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,

     2  stepTime, totalTime, dt, cmname, coordMp, charLength,

     3  props, density, strainInc, relSpinInc,

     4  tempOld, stretchOld, defgradOld, fieldOld,

     5  stressOld, stateOld, enerInternOld, enerInelasOld,

     6  tempNew, stretchNew, defgradNew, fieldNew,

c Write only -

     7  stressNew, stateNew, enerInternNew, enerInelasNew)

c

        include 'vaba_param.inc'

C

 dimension props(nprops), density(nblock),

     1  coordMp(nblock,*), nElement(nblock), nMatPoint(nblock),

     2  charLength(*), strainInc(nblock,ndir+nshr),

     3  relSpinInc(*), tempOld(*),nLayer(nblock),

     4  stretchOld(*), defgradOld(*),nSecPoint(nblock),

     5  fieldOld(*), stressOld(nblock,ndir+nshr),

     6  stateOld(nblock,nstatev), enerInternOld(nblock),

     7  enerInelasOld(nblock), tempNew(*),

     8  stretchNew(*), defgradNew(*), 

     9  fieldNew(*), stressNew(nblock,ndir+nshr), 

     1  stateNew(nblock,nstatev), 

     2  enerInternNew(nblock), enerInelasNew(nblock)

 character*80 cmname

        parameter(zero = 0.d0, one = 1.d0, two = 2.d0, half = .5d0,

     &   Three = 3.d0, op5 = 1.5d0) 

C Assume purely elastic material at the beginning of the analysis  

      E = props(1)

      xnu = props(2)

  yields1= props(3)

  eqpl1 = props(4)

               yields2= props(5)

  eqpl2 = props(6)      

      twomu = E /(one + xnu)

        alambda = xnu*twomu/(one-two*xnu)

      thremu = op5 * twomu

      nvalue = (nprops/two)-1       

       if (stepTime.eq.zero) then  

        do k = 1,nblock

              trace=strainInc(k,1)+strainInc(k,2)+strainInc(k,3)

   stressNew(k,1)=stressOld(k,1)

     &             + twomu*strainInc(k,1)+ alambda*trace

        stressNew(k,2)=stressOld(k,2)

     &             + twomu*strainInc(k,2)+ alambda*trace

               stressNew(k,3)=stressOld(k,3)

     &             + twomu*strainInc(k,3)+ alambda*trace

        stressNew(k,4)=stressOld(k,4)+twomu*strainInc(k,4)

        if (nshr.gt.1) then

      stressNew(k,5)=stressOld(k,5)+twomu*strainInc(k,5)

      stressNew(k,6)=stressOld(k,6)+twomu*strainInc(k,6)

      end if       

end do

else

do k=1,nblock   

c    

         peeqOld = stateOld(k,1)

  call vuhard(yieldOld, hard, peeqOld, props(3), nvalue)

 trace=strainInc(k,1)+strainInc(k,2)+strainInc(k,3)      

C Trial Stress

   s11=stressOld(k,1)+twomu*strainInc(k,1) 

     &  + alambda*trace

            s22=stressOld(k,2)+twomu*strainInc(k,2) 

     &  + alambda*trace

            s33=stressOld(k,3)+twomu*strainInc(k,3) 

     &  + alambda*trace

            s12=stressOld(k,4)+twomu*strainInc(k,4)

            if (nshr.gt.1) then 

   s13=stressOld(k,5)+twomu*strainInc(k,5) 

s23=stressOld(k,6)+twomu*strainInc(k,6)

end if      

            smean = (1/three)*(s11+s22+s33)

            s11 = s11 - smean

            s22 = s22 - smean

            s33 = s33 - smean               

           if (nshr.eq.1) then

            vmises = sqrt(op5*(s11*s11+s22*s22+s33*s33

     &      +two*s12*s12))  

           else

           vmises = sqrt(op5*(s11*s11+s22*s22+s33*s33

     &      +two*s12*s12+two*s13*s13+two*s23*s23))  

           end if                            

  sigdif=vmises-yieldOld     

      facyld=zero   

  if (sigdif.gt.zero)facyld=one  

C update the stress              

c   deqps = facyld * sigdif /(thremu + hard)

  yieldNew = yieldOld + hard * deqps

C update the stress

  factor = yieldNew/( yieldNew + thremu * deqps)             

  stressNew(k,1) = s11 * factor + smean

      stressNew(k,2) = s22 * factor + smean

      stressNew(k,3) = s33 * factor + smean

  stressNew(k,4) = s12 * factor

  if (nshr.gt.1)then

      stressNew(k,5) = s13 * factor

      stressNew(k,6) = s23 * factor 

      end if  

 This is only the part of the code..and the full code is working...my problem is, in this code..do anyone know how the expression called 'factor' is derived??

I have been looking.. how this expression has been derived..but could not find any clue for it...i mean..many people just mentioned this expression but did not mentioned how it has been derived.....if any body know it..could you please let me know..how it has been derived...any clue will be much apprciated..is there any book which this derivation is mentioned??? Thanks a lot...

regards