User login

You are here

cohesive model of Needleman was input into abaqus, but wrong???

 I have written  vumat codes of a cohesive model of Needleman (whose paper entilted "Numerical simulations of fast crack growth in brittle solids"), and I have implemented it into ABAQUS/Explicit, but the predicted result is wrong. I don't know why. The corresponding fortran file is as follows:

      subroutine vumat(

C Read only (unmodifiable)variables -

     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 (modifiable) variables -

     7  stressNew, stateNew, enerInternNew, enerInelasNew )

C

      include 'vaba_param.inc'

C

      dimension props(nprops), density(nblock), coordMp(nblock,*),

     1  charLength(nblock), strainInc(nblock,ndir+nshr),

     2  relSpinInc(nblock,nshr), tempOld(nblock),

     3  stretchOld(nblock,ndir+nshr),

     4  defgradOld(nblock,ndir+nshr+nshr),

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

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

     7  enerInelasOld(nblock), tempNew(nblock),

     8  stretchNew(nblock,ndir+nshr),

     8  defgradNew(nblock,ndir+nshr+nshr),

     9  fieldNew(nblock,nfieldv),

     1  stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),

     2  enerInternNew(nblock), enerInelasNew(nblock)

C

      character*80 cmname

C

C  six parameters

      Parameter(EX=2.71828)

      snmax=props(1)  ! the maxinum normal stress

      deltn=props(2)  ! the maxinum normal displacement

      stmax=props(3)   ! the maxinum shear stress

      delts=props(4)  ! the maxinum shear displacement

      t=props(5)   !depth of cohesive elements

      r=props(6)   

 

C

C)

C the maxinum energy

      egyn=EX*snmax*deltn

      egyt=SQRT(EX/2)*stmax*delts

      q=egyt/egyn

C

C)

C

      do 100 i=1, nblock

C

C)))) calculate strain

        sta1=stateOld(i,1)+strainInc(i,1)

        sta2=stateOld(i,2)+strainInc(i,2)

        sta3=stateOld(i,3)+strainInc(i,3)

        sta4=stateOld(i,4)+strainInc(i,4)

calculate  displacement

        detn=sta1*t

        dets=sta2*t

        

        tn=detn/deltn

        ts=dets/delts

 

C)calculate stress and energy

 

        sig1=(egyn/deltn)*(EX**(-tn))*(tn*(EX**(-ts*ts))+((1-q)/(r-1))

     1   *(1-(EX**(-ts*ts)))*(r-tn))

        sig2=(egyn/deltn)*2*(deltn/delts)*ts*(q+((r-q)/(r-1))*tn)

     2   *(EX**(-tn))*(EX**(-ts*ts))

        sig3=0

        sig4=0

        

        ENGY=egyn+egyn*(EX**(-tn))*((1-r+tn)*((1-q)/(r-1))

     1   -(q+((r-q)/(r-1))*tn)*(EX**(-ts*ts)))

 

C calculate the pure normal and shear energy

        ENt=egyn+egyn*((1-r)*((1-q)/(r-1))

     1   -q*(EX**(-ts*ts)))

        ENn=egyn+egyn*(EX**(-tn))*((1-r+tn)*((1-q)/(r-1))

     1   -(q+((r-q)/(r-1))*tn))

 

C

C................upate the stress ##################################,,

C

        stressNew(i,1)=sig1

        stressNew(i,2)=sig2

        stressNew(i,3)=sig3

        stressNew(i,4)=sig4

        

 

C) calculate the energy difference between the predicted energy and the maxinum energy

        detENt=egyt-ENt

        detENn=egyn-ENn

 

C.....,update the state variables

        stateNew(i,1)=stateOld(i,1)+strainInc(i,1)

        stateNew(i,2)=stateOld(i,2)+strainInc(i,2)

        stateNew(i,3)=stateOld(i,3)+strainInc(i,3)

        stateNew(i,4)=stateOld(i,4)+strainInc(i,4)

        

C   update the state variables

        stateNew(i,5)=2*sta1

        stateNew(i,6)=dets

        stateNew(i,7)=detn

        stateNew(i,8)=ENGY

        stateNew(i,9)=ENt

        stateNew(i,10)=ENn

        stateNew(i,11)=detENt

        stateNew(i,12)=detENn

        

C..........predict the fracture........................,,

        

        if(stateNew(i,11)>0.01) then

          stateNew(i,13)=1

        else

          stateNew(i,13)=0

        endif

 

C...........predict the fracture....................,,

        if(stateNew(i,12)>0.01) then

          stateNew(i,14)=1

        else

          stateNew(i,14)=0

        endif 

 

C..........predict the fracture......,,

        

        if( (stateNew(i,13)==0) .or. (stateNew(i,14)==0) ) then

          stateNew(i,15)=1

        else

          stateNew(i,15)=0

        endif 

        

C

100    continue

C

      return

      end

the abaqus input file is at the attachment

AttachmentSize
abaqus input file17.72 KB
Subscribe to Comments for "cohesive model of Needleman was input into abaqus, but wrong???"

Recent comments

More comments

Syndicate

Subscribe to Syndicate