User friction model in Abaqus : how to code an hysteretic friction behaviour ?

Hello everyone,

I have to use an hysteretic friction model. I am trying to code this in Abaqus by using the FRIC subroutine. The language is Fortran 77.

I am testing my code on a very simple model : two plates sliding on eachother. I apply a 1 cycle sinus displacement at the end of one plate and a apply some boundaries at the other side.

The model works fine in loading, but it I want it change when unloading. I am trying to check for slip reversal at the end of loading but it does not seem to work. I have tried many other things but I cannot obtain a proper hysteresis loop...

I have attached a view of the model and the subroutine file.

Any help would be very appreciated ! Thanks.

Rémi

Here is the subroutine if you don't want to download the file : 

      SUBROUTINE FRIC(LM,TAU,DDTDDG,DDTDDP,DSLIP,SED,SFD,

     1     DDTDDT,PNEWDT,STATEV,DGAM,TAULM,PRESS,DPRESS,DDPDDH,

     2     SLIP,KSTEP,KINC,TIME,DTIME,NOEL,CINAME,SLNAME,

     3     MSNAME,NPT,NODE,NPATCH,COORDS,RCOORD,DROT,TEMP,

     4     PREDEF,NFDIR,MCRD,NPRED,NSTATV,CHRLNGTH,PROPS,NPROPS)

C                    

      INCLUDE 'ABA_PARAM.INC'

C                    

      CHARACTER*80 CINAME,SLNAME,MSNAME

      DIMENSION TAU(NFDIR),DDTDDG(NFDIR,NFDIR),DDTDDP(NFDIR),

     1     DSLIP(NFDIR),DDTDDT(NFDIR,2),STATEV(*),

     2     DGAM(NFDIR),TAULM(NFDIR),SLIP(NFDIR),

     3     COORDS(MCRD),RCOORD(MCRD),DROT(2,2),TEMP(2),

     4     PREDEF(2,*),TIME(2),PROPS(NPROPS)

C                    

      DIMENSION GAMMA(2),TAURV(2),GAMRV(2)

C                    

      PARAMETER(ZERO=0.0D0)

      REAL FINC,PM,DSLIPEQ,DGAMEQ,CPT

C  

C  

 A  = PROPS(1)

 B  = PROPS(2)

 C  = PROPS(3)

 D  = PROPS(4)

 E  = PROPS(5)

C                      

      IF (LM .EQ. 2) THEN 

       RETURN

      ENDIF

C  

      IF ((KSTEP .EQ. 1) .AND. (KINC .EQ. 1)) THEN 

        CPT = ZERO

      ENDIF

C

      LM = 0

C  

      IF (PRESS .LT. ZERO) THEN 

          RETURN 

      ELSE

C  

      GAMMA(1) = STATEV(1) + DGAM(1)

      GAMMA(2) = STATEV(2) + DGAM(2)      

C

      DGAMEQ = (DGAM(1) + DGAM(2))/2

      DSLIPEQ = (DSLIP(1)+ DSLIP(2))/2

C  

     IF (DGAMEQ*DSLIPEQ .GT. ZERO) THEN     

 CPT = 1

 PM = PRESS

 TAURV(1) = TAU(1)

 TAURV(2) = TAU(2)

 GAMRV(1) = GAMMA(1)

 GAMRV(2) = GAMMA(2)

  C

 DDTDDG(1,1) = ZERO

 DDTDDG(1,2) = ZERO

 DDTDDG(2,1) = ZERO  

 DDTDDG(2,2) = ZERO

 DDTDDP(1)   = TAURV(1)

 DDTDDP(2)   = TAURV(2)

C  

 TAU(1) = TAURV(1)*PRESS

 TAU(2) = TAURV(2)*PRESS

C  

 DSLIP(1)  = DGAM(1)

 DSLIP(2)  = DGAM(2)

C 

 STATEV(1) = GAMMA(1)

 STATEV(2) = GAMMA(2)

C  

          ELSE

     IF (CPT .EQ. ZERO) THEN

C   

 DDTDDG(1,1) = A*B*EXP(-B*GAMMA(1))*PRESS

 DDTDDG(1,2) = ZERO

 DDTDDG(2,1) = ZERO  

 DDTDDG(2,2) = A*B*EXP(-B*GAMMA(2))*PRESS

 DDTDDP(1)   = A*(1-EXP(-B*GAMMA(1)))

 DDTDDP(2)   = A*(1-EXP(-B*GAMMA(2)))

C

 TAU(1) = (A*(1-EXP(-B*GAMMA(1))))*PRESS

 TAU(2) = (A*(1-EXP(-B*GAMMA(2))))*PRESS

C

 DSLIP(1)  = DGAM(1)

 DSLIP(2)  = DGAM(2)

C

 STATEV(1) = GAMMA(1)

 STATEV(2) = GAMMA(2)

C

              ELSEIF (CPT .EQ. 1) THEN

C    

 DDTDDG(1,1) = C*PRESS

 DDTDDG(1,2) = ZERO

 DDTDDG(2,1) = ZERO  

 DDTDDG(2,2) = C*PRESS

 DDTDDP(1)   = (C*(GAMMA(1)-GAMRV(1))+TAURV(1))

 DDTDDP(2)   = (C*(GAMMA(2)-GAMRV(2))+TAURV(2))

C  

 TAU(1) = (C*(GAMMA(1)-GAMRV(1))+TAURV(1))*PRESS

 TAU(2) = (C*(GAMMA(2)-GAMRV(2))+TAURV(2))*PRESS

C

 DSLIP(1)  = DGAM(1)

 DSLIP(2)  = DGAM(2)

C 

 STATEV(1) = GAMMA(1)

 STATEV(2) = GAMMA(2)

              ENDIF  

          ENDIF

      ENDIF

      RETURN

      END