UMAT subroutine - Element Deletion
UMAT subroutine - Element Deletion
(OP)
Hello guys,
i am writing a UMAT user subroutine in abaqus to remove element that has reached damage criteria. I have been trying different ways but no success. Not only do the elements remain but they also still have high stress values. So I am wondering if there would be a way/routine to set state or status of failed elements in next increment to be deleted/deactivated. I am working with Abaqus/Standard. U guys are my last hope.
Thank you very much.
Here is my code
SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1RPL,DDSDDT,DRPLDE,DRPLDT,
2STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
3NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
4CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
IMPLICIT NONE
! EXPLICIT DECLARATIONS OF ABAQUS UMAT INTERFACE VARIABLES
INTEGER :: LAYER, NPT, NOEL, KSPT, NPROPS,KSTEP,
1NSTATV, NTENS, NSHR, NDI, KINC
REAL*8 :: SSE, SPD, SCD, RPL, DRPLDT, DTIME, TEMP, DTEMP, PNEWDT,
2CELENT, PREDEF(1), DPRED(1)
REAL*8 :: STRESS(NTENS),STATEV(NSTATV),
1DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
2STRAN(NTENS),DSTRAN(NTENS),TIME(2),PROPS(NPROPS),COORDS(3),
4DROT(3,3),DFGRD0(3,3),DFGRD1(3,3), JSTEP(4)
CHARACTER*80 CMNAME
!!!
REAL*8 :: KAPPA,KAPPANP1,KAPPA_NEU,S,DF,DFNEU,FD,KAPPA_D,PSIH,ETA
INTEGER :: I,J,K,L,itw
! INITIALISE TENSOR-VALUED 3X3 COUNTERPARTS TO ABAQUS 1X6 VOIGT VECTOR
REAL*8 :: SIG(3,3),EPS(3,3),EPSE(3,3),SIGMA(3,3),TENSOR(3,3,3,3),
1DSIGDEPS(3,3,3,3),E,NU
REAL*8 :: ENU,OMN,OPN,OMZN,EOMN,HILF1,HILF2,HILF3,FLAG
! REAL*8 :: ELASTI(3,3,3,3) !,LAM,MUE,MUE2,KRON(3,3)
!!
itw=6
! if ((KSTEP.eq.1).and.(KINC.eq.0)) THEN
! if (TIME(1).LE.1.0d-08) then
! der integrationspunkt ist vorhanden.
! FLAG=1.0d+00
! else
! FLAG=STATEV(3)
! FLAG=1.1d+01
! end if
IF (STATEV(3).LT.0.1D+00) THEN
RETURN
END IF
!NTENS = 6
!DSIGDEPS = TENSOR
! KAPPA = 0.0D0
ETA = 0.15D0
KAPPA_D = 0.05D0
!!!
!E = 200
!NU = 0.3
!DO I=1,3
! DO J=1,3
! SIG(I,J)=0.0
! ENDDO
!ENDDO
!DO I = 1,NTENS
! STRESS(I) = 0.0
!ENDDO
CALL TENSOR3X3FROMABAQUSVOIGT(STRESS,SIG)
!!!
!DO S = 0.0,0.05,0.0005
! DO I = 1,NTENS
! DSTRAN(I) = 0.0
! ENDDO
! DSTRAN(1) = S
! !!
! STRAN(1) = 0.0D0
! STRAN(2) = 0.01D0
! STRAN(3) = 0.01D0
! STRAN(4) = 0.01D0
! STRAN(5) = 0.01D0
! STRAN(6) = 0.01D0
!!!!
EPS(1,1) = STRAN(1)
EPS(2,2) = STRAN(2)
EPS(3,3) = STRAN(3)
EPS(1,2) = 0.5D0*STRAN(4)
EPS(1,3) = 0.5D0*STRAN(5)
EPS(2,3) = 0.5D0*STRAN(6)
EPS(2,1) = EPS(1,2)
EPS(3,1) = EPS(1,3)
EPS(3,2) = EPS(2,3)
!!!!
EPSE(1,1) = STRAN(1) + DSTRAN(1)
EPSE(2,2) = STRAN(2) + DSTRAN(2)
EPSE(3,3) = STRAN(3) + DSTRAN(3)
EPSE(1,2) = 0.5D0*(STRAN(4) + DSTRAN(4))
EPSE(1,3) = 0.5D0*(STRAN(5) + DSTRAN(5))
EPSE(2,3) = 0.5D0*(STRAN(6) + DSTRAN(6))
EPSE(2,1) = EPSE(1,2)
EPSE(3,1) = EPSE(1,3)
EPSE(3,2) = EPSE(2,3)
!DO I = 1,NTENS
! DO J = 1,NTENS
! EPS(I,J) = STRAN(I)+DSTRAN(I)
! ENDDO
!ENDDO
!!!!
E = PROPS(1)
NU = PROPS(2)
ENU=E*NU
OMN=1.0D+00-NU
OPN=1.0D+00+NU
OMZN=1.0D+00-2.0D+00*NU
EOMN=E*OMN
HILF1=EOMN/(OPN*OMZN)
HILF2=E/(2.0D0*OPN)
HILF3=ENU/(OPN*OMZN)
!!!
DO I = 1,NTENS
DO J = 1,NTENS
DDSDDE(I,J) = 0.0D0
ENDDO
ENDDO
DDSDDE(1,1) = HILF1
DDSDDE(2,2) = HILF1
DDSDDE(3,3) = HILF1
DDSDDE(4,4) = HILF2
DDSDDE(5,5) = HILF2
DDSDDE(6,6) = HILF2
DDSDDE(1,2) = HILF3
DDSDDE(1,3) = HILF3
DDSDDE(2,1) = HILF3
DDSDDE(2,3) = HILF3
DDSDDE(3,1) = HILF3
DDSDDE(3,2) = HILF3
!!!!
!CALL STRAIN_TO_TENSOR1(STRAN,DSTRAN,EPS)
CALL T4FROMVOIGT(DDSDDE,TENSOR)
! DO I=1,3,1
! DO J=1,3,1
! DO K=1,3,1
! DO L=1,3,1
! TENSOR(I,J,K,L)=ELASTI(I,J,K,L)
! END DO
! END DO
! END DO
! END DO
!WRITE(ITW,*)PSIH(EPS,TENSOR)
DF=STATEV(2)
KAPPA = STATEV(1)
KAPPANP1 = KAPPA_NEU(EPS,TENSOR,KAPPA,itw)
! KAPPANP1 = KAPPA_NEU(EPSE,TENSOR,KAPPA,ETA,KAPPA_D,DF,itw)
!PRINT*,'KAPPANP1 = '
!WRITE(ITW,*)KAPPANP1
DFNEU = FD(KAPPANP1,ETA,KAPPA_D)
! IF (DFNEU.LE.0.5d+00) THEN
! der integrationspunkt ist defekt und soll geloescht werden.
! FLAG=0.0d+00
! ENDIF
IF (DFNEU.LT.0.2D+00) THEN
STATEV(3)=0.0D+00
do i=1,ntens,1
stress(i)=0.0d+00
end do
do i=1,ntens,1
do j=1,ntens,1
ddsdde(i,j)=1.0d+03
end do
end do
return
END IF
! DFNEU = 1.0D0
!PRINT*,'DF ='
!WRITE(ITW,*)DF
!!
! WRITE(*,'(1X,A)') 'TENSOR, ELASTI'
! DO I=1,3,1
! DO J=1,3,1
! DO K=1,3,1
! WRITE(*,'(1X,3I1,3(2E14.7,2X))') I,J,K,(TENSOR(I,J,K,L),
! A ELASTI(I,J,K,L),L=1,3,1)
! END DO
! END DO
! END DO
! WRITE(*,'(1X,A)') 'EPSE:'
! DO I=1,3,1
! WRITE(*,'(1X,I1,3E14.7)') I,(EPSE(I,J),J=1,3,1)
! END DO
! if (FLAG.gt.0.99d+00) then
DO I=1,NDI
DO J=1,NDI
SIGMA(I,J)=0.0D0
DO K=1,NDI
DO L=1,NDI
SIGMA(I,J) = SIGMA(I,J) + (TENSOR(I,J,K,L)*EPSE(K,L))
END DO
END DO
SIGMA(I,J)=DFNEU*SIGMA(I,J)
END DO
END DO
! else
! DO I=1,NDI
! DO J=1,NDI
! SIGMA(I,J)=1.0D-02
! END DO
! END DO
! end if
! WRITE(*,'(1X,A)') 'SIGMA:'
! DO I=1,NDI,1
! WRITE(*,'(1X,I1,3E14.7)') I,(SIGMA(I,J),J=1,3,1)
! END DO
!!!!!!
!!!!!
! DO I=1,NDI
! DO J=1,NDI
! DO K=1,NDI
! DO L=1,NDI
! SIGMA(K,L)=0.0D0
! SIGMA(K,L) = SIGMA(K,L) + (TENSOR(I,J,K,L)*EPS(I,J))
! END DO
! END DO
! END DO
! END DO
!!!!!
!!!!!
DO I=1,NDI
DO J=1,NDI
DO K=1,NDI
DO L=1,NDI
DSIGDEPS(I,J,K,L) = DFNEU*TENSOR(I,J,K,L)
ENDDO
ENDDO
ENDDO
ENDDO
! WRITE(*,'(1X,A)') 'TENSOR, DSIGDEPS'
! DO I=1,3,1
! DO J=1,3,1
! DO K=1,3,1
! WRITE(*,'(1X,3I1,3(2E14.7,2X))') I,J,K,(TENSOR(I,J,K,L),
! A DSIGDEPS(I,J,K,L),L=1,3,1)
! END DO
! END DO
! END DO
!!!!
! ROLAND SCHROEDER
! DO I=1,NDI
! DO J=1,NDI
! SIG(I,J)=0.0D0
! DO K=1,NDI
! DO L=1,NDI
! SIG(I,J) = SIG(I,J) + (DSIGDEPS(I,J,K,L)*EPS(K,L))
! ENDDO
! ENDDO
! ENDDO
! ENDDO
! ROLAND SCHROEDER
! WRITE(*,'(1X,A)') 'SIG:'
! DO I=1,NDI,1
! WRITE(*,'(1X,I1,3E14.7)') I,(SIG(I,J),J=1,3,1)
! END DO
! CALL TENSOR3X3TOABAQUSVOIGT(SIG,STRESS,ITW)
CALL TENSOR3X3TOABAQUSVOIGT(SIGMA,STRESS,itw)
!!!!
!DO I=1,NDI
! DO J=1,NDI
! DO K=1,NDI
! DO L=1,NDI
! WRITE(ITW,*)'DSIGDEPS(',I,J,K,L,')',DSIGDEPS(I,J,K,L)
! ENDDO
! ENDDO
! ENDDO
!ENDDO
! AENDERUNG ROLAND SCHROEDER
!CALL T4TOVOIGT(DSIGDEPS, DDSDDE,ITW)
!!!!!
!PRINT*,'SIG(I,J) = '
!DO I=1,NDI
! WRITE(ITW,*)(SIG(I,J),J=1,NDI)
!ENDDO
! AENDERUNG ROLAND SCHROEDER
!CALL TENSOR3X3TOABAQUSVOIGT(SIG,STRESS,ITW)
! AENDERUNG ROLAND SCHROEDER
! ROLAND SCHROEDER
! CALL T4TOVOIGT(TENSOR, DDSDDE,ITW)
CALL T4TOVOIGT(DSIGDEPS, DDSDDE,itw)
STATEV(1) = KAPPANP1
STATEV(2) = DFNEU
! STATEV(3) = FLAG
! WRITE(ITW,'(A,2D14.7)') 'KAPPANP1, DF:',KAPPANP1, DF
!!
!READ(*,*)
!ENDDO
!STOP
RETURN
END SUBROUTINE UMAT
NB : I used a damage model so changing eta und kappa either speeds the damage or slows it.
i am writing a UMAT user subroutine in abaqus to remove element that has reached damage criteria. I have been trying different ways but no success. Not only do the elements remain but they also still have high stress values. So I am wondering if there would be a way/routine to set state or status of failed elements in next increment to be deleted/deactivated. I am working with Abaqus/Standard. U guys are my last hope.
Thank you very much.
Here is my code
SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1RPL,DDSDDT,DRPLDE,DRPLDT,
2STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
3NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
4CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
IMPLICIT NONE
! EXPLICIT DECLARATIONS OF ABAQUS UMAT INTERFACE VARIABLES
INTEGER :: LAYER, NPT, NOEL, KSPT, NPROPS,KSTEP,
1NSTATV, NTENS, NSHR, NDI, KINC
REAL*8 :: SSE, SPD, SCD, RPL, DRPLDT, DTIME, TEMP, DTEMP, PNEWDT,
2CELENT, PREDEF(1), DPRED(1)
REAL*8 :: STRESS(NTENS),STATEV(NSTATV),
1DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
2STRAN(NTENS),DSTRAN(NTENS),TIME(2),PROPS(NPROPS),COORDS(3),
4DROT(3,3),DFGRD0(3,3),DFGRD1(3,3), JSTEP(4)
CHARACTER*80 CMNAME
!!!
REAL*8 :: KAPPA,KAPPANP1,KAPPA_NEU,S,DF,DFNEU,FD,KAPPA_D,PSIH,ETA
INTEGER :: I,J,K,L,itw
! INITIALISE TENSOR-VALUED 3X3 COUNTERPARTS TO ABAQUS 1X6 VOIGT VECTOR
REAL*8 :: SIG(3,3),EPS(3,3),EPSE(3,3),SIGMA(3,3),TENSOR(3,3,3,3),
1DSIGDEPS(3,3,3,3),E,NU
REAL*8 :: ENU,OMN,OPN,OMZN,EOMN,HILF1,HILF2,HILF3,FLAG
! REAL*8 :: ELASTI(3,3,3,3) !,LAM,MUE,MUE2,KRON(3,3)
!!
itw=6
! if ((KSTEP.eq.1).and.(KINC.eq.0)) THEN
! if (TIME(1).LE.1.0d-08) then
! der integrationspunkt ist vorhanden.
! FLAG=1.0d+00
! else
! FLAG=STATEV(3)
! FLAG=1.1d+01
! end if
IF (STATEV(3).LT.0.1D+00) THEN
RETURN
END IF
!NTENS = 6
!DSIGDEPS = TENSOR
! KAPPA = 0.0D0
ETA = 0.15D0
KAPPA_D = 0.05D0
!!!
!E = 200
!NU = 0.3
!DO I=1,3
! DO J=1,3
! SIG(I,J)=0.0
! ENDDO
!ENDDO
!DO I = 1,NTENS
! STRESS(I) = 0.0
!ENDDO
CALL TENSOR3X3FROMABAQUSVOIGT(STRESS,SIG)
!!!
!DO S = 0.0,0.05,0.0005
! DO I = 1,NTENS
! DSTRAN(I) = 0.0
! ENDDO
! DSTRAN(1) = S
! !!
! STRAN(1) = 0.0D0
! STRAN(2) = 0.01D0
! STRAN(3) = 0.01D0
! STRAN(4) = 0.01D0
! STRAN(5) = 0.01D0
! STRAN(6) = 0.01D0
!!!!
EPS(1,1) = STRAN(1)
EPS(2,2) = STRAN(2)
EPS(3,3) = STRAN(3)
EPS(1,2) = 0.5D0*STRAN(4)
EPS(1,3) = 0.5D0*STRAN(5)
EPS(2,3) = 0.5D0*STRAN(6)
EPS(2,1) = EPS(1,2)
EPS(3,1) = EPS(1,3)
EPS(3,2) = EPS(2,3)
!!!!
EPSE(1,1) = STRAN(1) + DSTRAN(1)
EPSE(2,2) = STRAN(2) + DSTRAN(2)
EPSE(3,3) = STRAN(3) + DSTRAN(3)
EPSE(1,2) = 0.5D0*(STRAN(4) + DSTRAN(4))
EPSE(1,3) = 0.5D0*(STRAN(5) + DSTRAN(5))
EPSE(2,3) = 0.5D0*(STRAN(6) + DSTRAN(6))
EPSE(2,1) = EPSE(1,2)
EPSE(3,1) = EPSE(1,3)
EPSE(3,2) = EPSE(2,3)
!DO I = 1,NTENS
! DO J = 1,NTENS
! EPS(I,J) = STRAN(I)+DSTRAN(I)
! ENDDO
!ENDDO
!!!!
E = PROPS(1)
NU = PROPS(2)
ENU=E*NU
OMN=1.0D+00-NU
OPN=1.0D+00+NU
OMZN=1.0D+00-2.0D+00*NU
EOMN=E*OMN
HILF1=EOMN/(OPN*OMZN)
HILF2=E/(2.0D0*OPN)
HILF3=ENU/(OPN*OMZN)
!!!
DO I = 1,NTENS
DO J = 1,NTENS
DDSDDE(I,J) = 0.0D0
ENDDO
ENDDO
DDSDDE(1,1) = HILF1
DDSDDE(2,2) = HILF1
DDSDDE(3,3) = HILF1
DDSDDE(4,4) = HILF2
DDSDDE(5,5) = HILF2
DDSDDE(6,6) = HILF2
DDSDDE(1,2) = HILF3
DDSDDE(1,3) = HILF3
DDSDDE(2,1) = HILF3
DDSDDE(2,3) = HILF3
DDSDDE(3,1) = HILF3
DDSDDE(3,2) = HILF3
!!!!
!CALL STRAIN_TO_TENSOR1(STRAN,DSTRAN,EPS)
CALL T4FROMVOIGT(DDSDDE,TENSOR)
! DO I=1,3,1
! DO J=1,3,1
! DO K=1,3,1
! DO L=1,3,1
! TENSOR(I,J,K,L)=ELASTI(I,J,K,L)
! END DO
! END DO
! END DO
! END DO
!WRITE(ITW,*)PSIH(EPS,TENSOR)
DF=STATEV(2)
KAPPA = STATEV(1)
KAPPANP1 = KAPPA_NEU(EPS,TENSOR,KAPPA,itw)
! KAPPANP1 = KAPPA_NEU(EPSE,TENSOR,KAPPA,ETA,KAPPA_D,DF,itw)
!PRINT*,'KAPPANP1 = '
!WRITE(ITW,*)KAPPANP1
DFNEU = FD(KAPPANP1,ETA,KAPPA_D)
! IF (DFNEU.LE.0.5d+00) THEN
! der integrationspunkt ist defekt und soll geloescht werden.
! FLAG=0.0d+00
! ENDIF
IF (DFNEU.LT.0.2D+00) THEN
STATEV(3)=0.0D+00
do i=1,ntens,1
stress(i)=0.0d+00
end do
do i=1,ntens,1
do j=1,ntens,1
ddsdde(i,j)=1.0d+03
end do
end do
return
END IF
! DFNEU = 1.0D0
!PRINT*,'DF ='
!WRITE(ITW,*)DF
!!
! WRITE(*,'(1X,A)') 'TENSOR, ELASTI'
! DO I=1,3,1
! DO J=1,3,1
! DO K=1,3,1
! WRITE(*,'(1X,3I1,3(2E14.7,2X))') I,J,K,(TENSOR(I,J,K,L),
! A ELASTI(I,J,K,L),L=1,3,1)
! END DO
! END DO
! END DO
! WRITE(*,'(1X,A)') 'EPSE:'
! DO I=1,3,1
! WRITE(*,'(1X,I1,3E14.7)') I,(EPSE(I,J),J=1,3,1)
! END DO
! if (FLAG.gt.0.99d+00) then
DO I=1,NDI
DO J=1,NDI
SIGMA(I,J)=0.0D0
DO K=1,NDI
DO L=1,NDI
SIGMA(I,J) = SIGMA(I,J) + (TENSOR(I,J,K,L)*EPSE(K,L))
END DO
END DO
SIGMA(I,J)=DFNEU*SIGMA(I,J)
END DO
END DO
! else
! DO I=1,NDI
! DO J=1,NDI
! SIGMA(I,J)=1.0D-02
! END DO
! END DO
! end if
! WRITE(*,'(1X,A)') 'SIGMA:'
! DO I=1,NDI,1
! WRITE(*,'(1X,I1,3E14.7)') I,(SIGMA(I,J),J=1,3,1)
! END DO
!!!!!!
!!!!!
! DO I=1,NDI
! DO J=1,NDI
! DO K=1,NDI
! DO L=1,NDI
! SIGMA(K,L)=0.0D0
! SIGMA(K,L) = SIGMA(K,L) + (TENSOR(I,J,K,L)*EPS(I,J))
! END DO
! END DO
! END DO
! END DO
!!!!!
!!!!!
DO I=1,NDI
DO J=1,NDI
DO K=1,NDI
DO L=1,NDI
DSIGDEPS(I,J,K,L) = DFNEU*TENSOR(I,J,K,L)
ENDDO
ENDDO
ENDDO
ENDDO
! WRITE(*,'(1X,A)') 'TENSOR, DSIGDEPS'
! DO I=1,3,1
! DO J=1,3,1
! DO K=1,3,1
! WRITE(*,'(1X,3I1,3(2E14.7,2X))') I,J,K,(TENSOR(I,J,K,L),
! A DSIGDEPS(I,J,K,L),L=1,3,1)
! END DO
! END DO
! END DO
!!!!
! ROLAND SCHROEDER
! DO I=1,NDI
! DO J=1,NDI
! SIG(I,J)=0.0D0
! DO K=1,NDI
! DO L=1,NDI
! SIG(I,J) = SIG(I,J) + (DSIGDEPS(I,J,K,L)*EPS(K,L))
! ENDDO
! ENDDO
! ENDDO
! ENDDO
! ROLAND SCHROEDER
! WRITE(*,'(1X,A)') 'SIG:'
! DO I=1,NDI,1
! WRITE(*,'(1X,I1,3E14.7)') I,(SIG(I,J),J=1,3,1)
! END DO
! CALL TENSOR3X3TOABAQUSVOIGT(SIG,STRESS,ITW)
CALL TENSOR3X3TOABAQUSVOIGT(SIGMA,STRESS,itw)
!!!!
!DO I=1,NDI
! DO J=1,NDI
! DO K=1,NDI
! DO L=1,NDI
! WRITE(ITW,*)'DSIGDEPS(',I,J,K,L,')',DSIGDEPS(I,J,K,L)
! ENDDO
! ENDDO
! ENDDO
!ENDDO
! AENDERUNG ROLAND SCHROEDER
!CALL T4TOVOIGT(DSIGDEPS, DDSDDE,ITW)
!!!!!
!PRINT*,'SIG(I,J) = '
!DO I=1,NDI
! WRITE(ITW,*)(SIG(I,J),J=1,NDI)
!ENDDO
! AENDERUNG ROLAND SCHROEDER
!CALL TENSOR3X3TOABAQUSVOIGT(SIG,STRESS,ITW)
! AENDERUNG ROLAND SCHROEDER
! ROLAND SCHROEDER
! CALL T4TOVOIGT(TENSOR, DDSDDE,ITW)
CALL T4TOVOIGT(DSIGDEPS, DDSDDE,itw)
STATEV(1) = KAPPANP1
STATEV(2) = DFNEU
! STATEV(3) = FLAG
! WRITE(ITW,'(A,2D14.7)') 'KAPPANP1, DF:',KAPPANP1, DF
!!
!READ(*,*)
!ENDDO
!STOP
RETURN
END SUBROUTINE UMAT
NB : I used a damage model so changing eta und kappa either speeds the damage or slows it.