User login

You are here

VUMAT of Drucker Prager, help needed please..

Gayan Aravinda's picture

Hi Everyone,

                    I have written Drucker Prager in VUMAT version as shown below..but still there are some deviation with the original in built version of ABAQUS. I have all the notes how i derived these equation.could you please show me where the code has gone wrong or have i missed a critical point in here..i tried my best to spot the error but could not..may be you all experts could find it. Thank you very much. 

 

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)

C

 character*80 cmname

C

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

     &   Three = 3.d0, op5 = 1.5d0,third = one/Three) 

C Assume purely elastic material at the beginning of the analysis

C      

      E = props(1)

      xnu = props(2)

      beta=props(3)

      sigmaT=props(4)

      K=props(5)

  yields1= props(6)

  eqpl1 = props(7)

               yields2= props(8)

  eqpl2 = props(9)      

      twomu = E /(one + xnu)

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

      thremu = op5 * twomu

      term1=one/K

      nvalue = (nprops/three)-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(6), nvalue)      

c   

    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     

C

            smean = (one/Three)*(s11+s22+s33)

C

            s11 = s11 - smean

            s22 = s22 - smean

            s33 = s33 - smean

C                

           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     

                          t=vmises

                          d=sigmaT*(term1+(third*tan(beta)))

                          DPYield=t-(-smean*tan(beta))-d

                          a=term1+(third*tan(beta))                     

c                                                                                                                                   

  sigdif=DPYield-yieldOld     

      facyld=zero   

  if (sigdif.le.zero)then

  factor =one  

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

               deqps=(sqrt(two)/three)*sqrt(((dps11-dps22)**two)+

     &         ((dps22-dps33)**two)+((dps11-dps33)**two)+six*

     &       ((dps12**two)+(dps23**two)+(dps13**two)))

                 yieldNew = yieldOld + hard * deqps          

    factor = yieldNew/(yieldNew + thremu * deqps) 

  else

c   facyld=one         

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

c  equivalant plastic strain expression 

               deqps=(sqrt(two)/three)*sqrt(((dps11-dps22)**two)+

     &         ((dps22-dps33)**two)+((dps11-dps33)**two)+six*

     &       ((dps12**two)+(dps23**two)+(dps13**two)))    

c plastic strain components for Drucker Prager criterion ( derived these equations from the flow rule shown in abaqus section 4.4.2)             

       dps11=((deqps/(a*two*vmises))*((two*s11)-s33-s22)

     &               - third*tan(beta))                   

               dps22=((deqps/(a*two*vmises))*((two*s22)-s11-s33)

     &               - third*tan(beta))

               dps33=((deqps/(a*two*vmises))*((two*s33)-s11-s22)

     &               - third*tan(beta))

               dps12=(deqps/(a*two*vmises))*six*s12

               dps13=(deqps/(a*two*vmises))*six*s13 

               dps23=(deqps/(a*two*vmises))*six*s23                   

               yieldNew = yieldOld + hard * deqps          

  factor = yieldNew/(yieldNew + thremu * deqps)

               end if                                                

C update the stress               

  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                   

C Update the state variables

              stateNew(k,1) = stateOld(k,1) + deqps

C Update the specific internal energy

          if (nshr.eq.1) then

          stressPower = half * (

     &        ( stressOld(k,1) + stressNew(k,1) ) * strainInc(k,1) +

     &        ( stressOld(k,2) + stressNew(k,2) ) * strainInc(k,2) +

     &        ( stressOld(k,3) + stressNew(k,3) ) * strainInc(k,3))+

     &        ( stressOld(k,4) + stressNew(k,4) ) * strainInc(k,4)   

          else  

          stressPower = half * (

     &        ( stressOld(k,1) + stressNew(k,1) ) * strainInc(k,1) +

     &        ( stressOld(k,2) + stressNew(k,2) ) * strainInc(k,2) +

     &        ( stressOld(k,3) + stressNew(k,3) ) * strainInc(k,3))+

     &        ( stressOld(k,4) + stressNew(k,4) ) * strainInc(k,4) + 

     &        ( stressOld(k,5) + stressNew(k,5) ) * strainInc(k,5) +    

     &        ( stressOld(k,6) + stressNew(k,6) ) * strainInc(k,6)    

          end if

         enerInternNew(k) = enerInternOld(k) + stressPower/density(k)

c

C Update the dissipated inelastic specific energy   

C

               plasticWorkInc = half*(yieldOld + yieldNew)*deqps

c

               enerInelasNew(k) = enerInelasOld(k)

     &         + plasticWorkInc / density(k)

 end do

end if

       return

end

c

c

subroutine vuhard(syield, hard, eqplas, table, nvalue)

        include 'vaba_param.inc'   

        dimension table(2,nvalue)

        parameter(zero=0.d0)   

C set yield stress to last value of table, hardening to zero

C

                  syield=table(1, nvalue)  

                  hard=zero            

C

C if more than one entry, search table  

    if (nvalue.gt.1)then

      do k1=1,nvalue-1

     eqpl1 = table(2,k1+1)  

if(eqplas.lt.eqpl1)then

 eqpl0 = table(2,k1)                          

C yield stress and hardening

                     deqpl=eqpl1-eqpl0                           

syield0=table(1,k1)

syield1=table(1,k1+1) 

dsyiel=syield1-syield0

hard=dsyiel/deqpl

syield=syield0+(eqplas-eqpl0)*hard

goto 40

            end if 

          end do

   40 continue

       end if

       

                return

   end

Comments

aslan mohammadpour's picture

Dear user

you had a post before and I answered  your post but you didn't reply that post, Plase let other know you have solved your problem...

Sometimes  "small" deviation in your answere and other Fe codes are quiet normal, I will tell you my experience ...

I wrote a full FE code in Fortran two years ago(Drucker prager model), and I  compared my code by ansys and I found some deviation...
My code was correct but the deviation was because of the number 1 in below. there may be  many  reason for deviation:

1- Ansys is using different shape function for elements(extra shape functions)  ( I used classical shape function)

2- Ansys may use different number of integration point

3- Ansys may use B_bar method for calculation of volumetric constrains...

After I changed the setting of the ansys to have same option as my code, my result became accurate up to four floating point.

In your case may be the same happens in abaqus, may be when you use Vumat the abaqus disables some option for elements...

The first step I suggest is to set initial yield condition vey high in order that plastic flow wont happen... then check your elastic answer with the one in original built-in abaqus... this helps you to make sure the deviation is not because of element options ....

 

 

Gayan Aravinda's picture

Hi, Aslan...thank you very much...i may forgotten soory about it...i solved the VUMAT for perfect plastic case.....thank you very much for the suggestions ..i will have a look at it...

ciprian.zub's picture

Hi Gayan Aravinda,

I am new with VUMAT. 

I kindly ask you for the VUMAT in the corrected version if possible.

I really would appreciate it.

ciprian.zub@student.upt.ro

Subscribe to Comments for "VUMAT of Drucker Prager, help needed please.."

Recent comments

More comments

Syndicate

Subscribe to Syndicate