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