Springpot UMAT implementation (fractional order derivatives)

Hi everyone,

I'm currently trying to implement the constitutive model of a simple springpot element in abaqus UMAT My entire subroutine is based in this paper https://www.sciencedirect.com/science/article/abs/pii/S0168874X17309356?via%3Dihub .

However, it doesn't seem to work, the solution always blows-off after the first iteration (it doesn't really matter what the time step is).

I have mainly two doubts about my implementation; the first one is that in order to account for incompressibility I defined the poisson modulus as 0.499 and defined stress(2)=stress(3)= -poisson*stress(1), is this correct?

The second is that I've noticed that abaqus is passing very low values for stran, in the order of E-10. Is this normal? The simulated test is uniaxial traction with a ramp of 0.002 meters over the entire time step. The analysis is static.

Finally, if someone has an implementation of the springpot that could share with me, that would be awesome. I'm struggling a lot to find a single example of this online. 

Thanks !!

Below you'll find my routine

Best

Nelson Sanchez

SUBROUTINE SDVINI(STATEV,COORDS,NSTATV,NCRDS,NOEL,NPT,LAYER,KSPT)
C
INCLUDE 'ABA_PARAM.INC'
C
DIMENSION STATEV(NSTATV),COORDS(NCRDS)
INTEGER I
DO I=1,NSTATV
STATEV(I)=0.0
END DO
RETURN
END
C
SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDOT,DRPLDE,DRPLDT,
2 STRAN,DSTRAN,
3 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
4 NDI,NSHR,NTENS,NSTATEV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
5 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
C
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 CMNAME
DIMENSION STRESS(NTENS),STATEV(NSTATEV), DDSDDE(NTENS,NTENS),
1 DDSDDT(NTENS),DRPLDE(NTENS),STRAN(NTENS),DSTRAN(NTENS),
2 TIME(2),PREDEF(1),DPRED(1),PROPS(NPROPS),COORDS(3),DROT(3,3),
3 DFGRD0(3,3), DFGRD1(3,3),VOLSTRAN(1200,6),
4 DIRSTRANAVG(1200)

REAL*8 A, LAMBDA, DIRSTRANAVG, TERM, VISC, E, ALPHA,
1 TERM11, TERM22, TERM33, TERM12, TERM13, TERM23
2 GL11, GL22, GL33, GL12, GL13, GL23, POISSON, RSTEP
3 JTERM1, JTERM2, JTERM3
C ------------------------------------- * *----------------------------------------------------- C
C First attempt to implement the springpot as an UMAT based on the paper of Gioacchino Alotta, et. al
C 2018 , the finite element implementation of 3D fractional viscoelastic constitutive models
C Props(1) is the fractional viscosity (or whatever that coefficient may be called)
C Props(2) is the derivative order (alpha coefficient)
C--------------------------------------------------------------------------------------------------- C
C Get the parameters (easier to manipulate)
VISC = PROPS(1)
ALPHA=PROPS(2)
POISSON=0.499D0
C It seems necessary to initialize the array VOLSTRAN where each column collects the history of each stress
C component according to abaqus conventions
DO J=1,6
DO K=1,1200
VOLSTRAN(K,J) = 0.D0
ENDDO
ENDDO
C--------------------------------------------------------------------------------------------------- C
C Calculate the GL fractional derivative starts here
C Initialize everything
LAMBDA=0
TERM11 = 0
TERM22 = 0
TERM33 = 0
TERM12 = 0
TERM13 = 0
TERM33 = 0
ETERM = 0
GL11 = 0
GL22 = 0
GL33 = 0
GL12 = 0
GL13 = 0
GL33 = 0
GLE = 0
RSTEP= 0
JTERM1=0
JTERM2=0
JTERM3=0
C Incompressibility is considered by setting the poisson modulus to 0.499 and stran(2) = stran(3) = -poisson*stran(1)
C Might not be correct tho
 STRAN(2)=-POISSON*STRAN(1)
 STRAN(3)=-POISSON*STRAN(1)
C---------------------------------------------------------------------------------------------------
C Now calculate the terms associated with the "memory" effect of the fractional derivative
C---------------------------------------------------------------------------------------------------
DIRSTRANAVG(KINC)=(STRAN(1) + STRAN(2) + STRAN(3))/3
VOLSTRAN(KINC,1) = STRAN(1) - DIRSTRANAVG(KINC)
VOLSTRAN(KINC,2) = STRAN(2) - DIRSTRANAVG(KINC)
VOLSTRAN(KINC,3) = STRAN(3) - DIRSTRANAVG(KINC)
VOLSTRAN(KINC,4) = STRAN(4)
VOLSTRAN(KINC,5) = STRAN(5)
VOLSTRAN(KINC,6) = STRAN(6)
C Additive part
C Leading GL coefficient is always positive = 1 and the rest negative
LAMBDA=1
DO J=1,KINC
TERM11 = TERM11 + LAMBDA*VOLSTRAN(KINC-J+1,1)
TERM22 = TERM22 + LAMBDA*VOLSTRAN(KINC-J+1,2)
TERM33 = TERM33 + LAMBDA*VOLSTRAN(KINC-J+1,3)
TERM12 = TERM12 + LAMBDA*VOLSTRAN(KINC-J+1,4)
TERM13 = TERM13 + LAMBDA*VOLSTRAN(KINC-J+1,5)
TERM23 = TERM23 + LAMBDA*VOLSTRAN(KINC-J+1,6)
ETERM = ETERM + LAMBDA*DIRSTRANAVG(KINC-J+1)
C Update the GL coefficient
LAMBDA = (J-1.D0-ALPHA)*(1.D0/J)*LAMBDA
ENDDO
C Multiply the addition of all terms involved GL derivative
C with the prefactor Cp*DT^alpha as in Alotta 2018 eq.12
GL11 = VISC*(DTIME**-ALPHA)*TERM11
GL22 = VISC*(DTIME**-ALPHA)*TERM22
GL33 = VISC*(DTIME**-ALPHA)*TERM33
GL12 = VISC*(DTIME**-ALPHA)*TERM12
GL13 = VISC*(DTIME**-ALPHA)*TERM13
GL23 = VISC*(DTIME**-ALPHA)*TERM23
GLE = VISC*(DTIME**-ALPHA)*ETERM
C---------------------------------------------------------------------------------------------------
C Now the stress strain stuff, first the volumetric terms according to Alotta 2018, eq. 11a
STRESS(1) = 2*VISC*GL11 + 3*VISC*GLE
STRESS(2) = 2*VISC*GL22 + 3*VISC*GLE
STRESS(3) = 2*VISC*GL33 + 3*VISC*GLE
C Now the deviatoric parts as in Alotta 2018, eq. 11b
STRESS(4) = GL12
STRESS(5) = GL13
STRESS(6) = GL23
LAMBDA=1
C Initialize the jacobian
DO K1 = 1,NTENS
DO K2 = 1,NTENS
DDSDDE = 0.D0
ENDDO
ENDDO
C Now give it values, brother
C Define first the components to be allocated in the Jacobian, to make it a bit clearer
C Normal and shear stresses have the same power law order in this case
JTERM1 = VISC*(DTIME**-ALPHA) + (4.D0/3.D0)*(VISC*(DTIME**-ALPHA))
JTERM2 = VISC*(DTIME**-ALPHA) - (2.D0/3.D0)*(VISC*(DTIME**-ALPHA))
JTERM3 = VISC*(DTIME**-ALPHA)
C Allocate terms in their respective places, eqs. 14 a,b, and c
DO K1=1,NDI
DO K2=1,NDI
DDSDDE(K2,K1)= JTERM2
ENDDO
ENDDO
DO K1=1,NDI
DDSDDE(K1,K1) = JTERM1
ENDDO
C Now define the terms related to the deviatoric terms
DO K1=NDI+1,NTENS
DDSDDE(K1,K1)= JTERM3
ENDDO
END