C 4X4 JACOBIAN THEN CONDENSATION WITH REGULARIZATION
c (energy due to viscous regularization is calculated)
SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,
2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 CMNAME
DIMENSION STRESS(NTENS),STATEV(NSTATV),
1 DDSDDE(NTENS,NTENS),
2 DDSDDT(NTENS),DRPLDE(NTENS),
3 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
4 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)
DIMENSION STRANT(6),TSTRANT(4)
DIMENSION CFULL(6,6),CDFULL(6,6)
DIMENSION DDFDE(6), DDMDE(6), DCDDF(6,6), DCDDM(6,6)
DIMENSION ATEMP1(6), ATEMP2(6), TDDSDDE(6,6)
DIMENSION OLD_STRESS(6)
DIMENSION DOLD_STRESS(6),D_STRESS(6)
PARAMETER (ZERO = 0.D0,ONE = 1.D0,TWO = 2.D0, HALF = 0.5D0)
C****************************
C STRANT..... STRAIN AT THE END OF THE INCREMENT
C TSTRANT.....TEMPORARY ARRAY TO HOLD THE STRAIN FOR PLANE STRESS PROBLEM
C CFULL.......FULL 6X6 ELASTICITY MATRIX
C CDFULL......FULL 6X6 DAMAGED ELASTICITY MATRIX
C DDFDE....... D DF/D E
C DDMDE....... D DM/D E
C DCDDF....... D C/ D DF THE DERIVATIVE OF THE FULL MATRIX OVER DF
C DCDDM........D C/ D DM THE DERIVATIVE OF THE FULL MATRIX OVER DM
C ATEMP1,ATEMP2...TEMPORARY ARRAY USED IN JACOBIAN CALCULATION
C TDDSDDE.....UNCONDENSED JACOBIAN MATRIX FOR PLANE STRESS PROBLEM
C OLD_STRESS...STRESS AT THE BEGINNING OF THE INCREMENT, SAVED FOR THE ENERGY
C COMPUTATION
C DOLD_STRESS...STRESS AT THE BEGINNING OF THE INCREMENT,
C IF THERE'S NO VISCOUS REGULARIZATION
C D_STRESS...STRESS IF THERE'S NO VISCOUS REGULARIZATION, THE ABOVE IS CALCULATED
C TO CALCULATE THE SCD, ENERGY CAUSED BY VISCOUS REGULARIZATION
C STATEV(1) damage variable df
C STATEV(2) damage variable dm
C STATEV(3) regularized damage variable dfv
C STATEV(4) regularizaed damage variable dmv
C STATEV(5:10) TEMPORARY ARRAYS TO SAVE DOLD_STRESS
C************
C
C GET THE MATERIAL PROPERTIES---ENGINEERING CONSTANTS
C
TENL = PROPS(1) !YOUNG'S MODULUS IN DIRECTION 1 (L)
TENT = PROPS(2) !YOUNG'S MODULUS IN DIRECTION 2 (T)
SHRLT = PROPS(3) !SHEAR MODULUS IN 12 PLANE
SHRTT = PROPS(4) !SHEAR MODULUS IN 23 PLANE
XNULT = PROPS(5) !POISON'S RATIO POI_12
XNUTT = PROPS(6) !POISON'S RATIO POI_23
XNUTL = XNULT / TENL * TENT !POI_21
C
C GET THE FAILURE PROPERTIES
C
SIGTL = PROPS(7) !FAILURE STRESS IN 1 DIRECTION IN TENSION
SIGCL = PROPS(8) !FAILURE STRESS IN 1 DIRECTION IN COMPRESSION
SIGTT = PROPS(9) !FAILURE STRESS IN 2 DIRECTION IN TENSION
SIGCT = PROPS(10) !FAILURE STRESS IN 2 DIRECTION IN COMPRESSION
SIGSLT = PROPS(11) !FAILURE STRESS IN SHEAR IN 1-2 PLANE
GFMAT = PROPS(12) !FRACTURE ENERGY IN MATRIX
GFFIB = PROPS(13) !FRACTURE ENERGY IN FIBER
ETA = PROPS(14) ! VISCOSITY FOR REGULARIZATION
C
C CALCULATE THE STRAIN AT THE END OF THE INCREMENT
C
DO I = 1, NTENS
STRANT(I) = STRAN(I) + DSTRAN(I)
END DO
c
C FILL THE 6X6 FULL STIFFNESS MATRIX
DO I = 1, 6
DO J = 1, 6
CFULL(I,J)=ZERO
END DO
END DO
ATEMP = ONE - TWO * XNULT * XNUTL - XNUTT ** TWO
1 - TWO * XNULT * XNUTL * XNUTT
CFULL(1,1) = TENL * (ONE - XNUTT ** TWO) / ATEMP
CFULL(2,2) = TENT * (ONE - XNULT * XNUTL) / ATEMP
CFULL(3,3) = CFULL(2,2)
CFULL(1,2) = TENT * (XNULT + XNULT * XNUTT) / ATEMP
CFULL(1,3) = CFULL(1,2)
CFULL(2,3) = TENT * (XNUTT + XNULT * XNUTL) / ATEMP
CFULL(4,4) = SHRLT
CFULL(5,5) = SHRLT
CFULL(6,6) = SHRTT
DO I = 2, 6
DO J = 1, I-1
CFULL(I,J) = CFULL(J,I)
END DO
END DO
c calculate the failure strain by failure stress
EPITL = SIGTL / cfull(1,1) !FAILURE STRAIN 1 DIRECTION IN TENSION
EPICL = SIGCL / cfull(1,1) !FAILURE STRAIN 1 DIRECTION IN COMPRESSION
EPITT = SIGTT / cfull(2,2) !TENSILE FAILURE STRAIN 2 DIRECTION
EPICT = SIGCT / cfull(2,2) !COMPRESSIVE FAILURE STRAIN 2 DIRECTION
EPISLT = SIGSLT/ SHRLT ! FAILURE SHEAR STRAIN ...ENGINEERING STRAIN
C
C CHECK THE FAILURE INITIATION CONDITION
c
DFOLD = STATEV(1)
DMOLD = STATEV(2)
DFVOLD = STATEV(3)
DMVOLD = STATEV(4)
CALL CheckFailureIni(EPITL,EPICL,EPITT,EPICT,EPISLT,STRANT,
1 GFMAT,GFFIB, CELENT, CFULL, DF, DM, DDFDE, DDMDE, NTENS,
2 DFOLD, DMOLD,NDI)
C
C ! USE VISCOUS REGULARIZATION
C
DFV = ETA / (ETA + DTIME) * DFVOLD + DTIME / (ETA + DTIME) * DF
DMV = ETA / (ETA + DTIME) * DMVOLD + DTIME / (ETA + DTIME) * DM
C SAVE THE OLD STRESS TO OLD_STRESS
DO I = 1, NTENS
OLD_STRESS(I) = STRESS(I)
END DO
C CALL ROUTINE TO CALCULATE THE STRESS
C CALCULATE THE STRESS IF THERE'S NO VISCOUS REGULARIZATION
CALL GetStress(CFULL,CDFULL,DF,DM,D_STRESS,STRANT,NDI,NTENS)
C CALCULATE THE STRESS IF THERE'S VISCOUS REGULARIZATION
CALL GetStress(CFULL,CDFULL,DFV,DMV,STRESS,STRANT,NDI,NTENS)
C GET THE OLD STRESS IF THERE'S NO VISCOUS REGULARIZATION
DO I=1,NTENS
DOLD_STRESS(I)=STATEV(I+4)
END DO
C SAVE THE CURRENT STRESS IF THERE'S NO VISCOUS REGULARIZATION
DO I=1,NTENS
STATEV(I+4)=D_STRESS(I)
END DO
C
C CALCULATE THE DERIVATIVE MATRIX DC/DDM, DC/DDF OF THE DAMAGED MATRIX
C
CALL ElasticDerivative(CFULL,DMV,DFV, DCDDM,DCDDF)
C
C UPDATE THE JACOBIAN
C
C FULL 3D CASE
IF (NDI .EQ. 3) THEN
DO I = 1, NTENS
ATEMP1(I) = ZERO
DO J = 1, NTENS
ATEMP1(I) = ATEMP1(I) + DCDDM(I,J) * STRANT(J)
END DO
END DO
DO I = 1, NTENS
ATEMP2(I) = ZERO
DO J = 1, NTENS
ATEMP2(I) = ATEMP2(I) + DCDDF(I,J) * STRANT(J)
END DO
END DO
DO I = 1, NTENS
DO J = 1, NTENS
DDSDDE(I,J)=CDFULL(I,J) + ( ATEMP1(I) * DDMDE(J)
1 + ATEMP2(I) * DDFDE(J) ) * DTIME / (DTIME + ETA)
END DO
END DO
C
C ! PLANE STRESS CASE
C
ELSE IF (NDI .EQ.2) THEN
TSTRANT(1) = STRANT(1)
TSTRANT(2) = STRANT(2)
TSTRANT(3) = -CDFULL(1,3) / CDFULL(3,3) * STRANT(1)
1 - CDFULL(2,3) / CDFULL(3,3) * STRANT(2)
TSTRANT(4) = STRANT(3)
DO I = 1, 4
ATEMP1(I) = ZERO
DO J = 1, 4
ATEMP1(I) = ATEMP1(I) + DCDDM(I,J) * TSTRANT(J)
END DO
END DO
DO I = 1, 4
ATEMP2(I) = ZERO
DO J = 1, 4
ATEMP2(I) = ATEMP2(I) + DCDDF(I,J) * TSTRANT(J)
END DO
END DO
DO I = 1,6
DO J = 1,6
TDDSDDE(I,J) = ZERO
END DO
END DO
C TO GET THE UNCONDENSED JACOBIAN FOR PLANE STRESS CASE
DO I = 1, NTENS
DO J = 1, NTENS
DDSDDE(I,J) = ZERO
END DO
END DO
DO I = 1, 4
DO J = 1, 4
TDDSDDE(I,J)=CDFULL(I,J) + ( ATEMP1(I) * DDMDE(J)
1 + ATEMP2(I) * DDFDE(J) ) * DTIME / (DTIME + ETA)
END DO
END DO
C
C CONDENSE THE JACOBIAN MATRIX FOR PLANE STRESS PROBLEM
C
CALL MatrixCondense(TDDSDDE,DDSDDE)
END IF
C
C TO UPDATE THE STATE VARIABLE
C
STATEV(1) = DF
STATEV(2) = DM
STATEV(3) = DFV
STATEV(4) = DMV
C
C TO COMPUTE THE ENERGY
C
DO I = 1, NDI
SSE = SSE + HALF * (STRESS(I) + OLD_STRESS(I)) * DSTRAN(I)
END DO
DO I = NDI+1, NTENS
SSE = SSE + (STRESS(I) + OLD_STRESS(I)) * DSTRAN(I)
END DO
C TO COMPUTE THE INTERNAL ENERGY WITHOUT VISCOUS REGULARIZATION
DO I = 1, NDI
SCD = SCD + HALF * (STRESS(I) + OLD_STRESS(I)
1 -D_STRESS(I)-DOLD_STRESS(I)) * DSTRAN(I)
END DO
DO I = NDI+1, NTENS
SCD = SCD + (STRESS(I) + OLD_STRESS(I)
1 -D_STRESS(I)-DOLD_STRESS(I)) * DSTRAN(I)
END DO
RETURN
END
C******************************************************************************
C CALCULATE THE STRESS BASED ON THE DAMAGE VARAIBLES***************************
C******************************************************************************
SUBROUTINE GetStress(CFULL,CDFULL,DFV,DMV,STRESS,STRANT,NDI,NTENS)
INCLUDE 'ABA_PARAM.INC'
DIMENSION CFULL(6,6),CDFULL(6,6),STRESS(NTENS),
1 STRANT(6),CDTHREE(3,3)
PARAMETER (ZERO=0.D0, ONE=1.D0)
C CDTHREE.....DAMAGED CONDENSED-ELASTICITY MATRIX FOR PLANE STRESS PROBLEM
DO I = 1, 6
DO J = 1, 6
CDFULL(I,J)=CFULL(I,J)
END DO
END DO
IF ( (DFV .NE. ZERO) .OR. (DMV .NE. ZERO)) THEN
CDFULL(1,1) = (ONE - DFV) * CFULL(1,1)
CDFULL(1,2) = (ONE - DFV) * (ONE - DMV) * CFULL(1,2)
CDFULL(2,1) = CDFULL(1,2)
CDFULL(2,2) = (ONE - DMV) * CFULL(2,2)
CDFULL(1,3) = (ONE - DFV) * CFULL(1,3)
CDFULL(3,1) = CDFULL(1,3)
CDFULL(2,3) = (ONE- DMV) * CFULL(2,3)
CDFULL(3,2) = CDFULL(2,3)
CDFULL(4,4) = (ONE - DMV) * (ONE - DFV) * CFULL(4,4)
END IF
C UPDATE THE STRESS STATE IF 3D CASE
C
IF (NDI .EQ. 3) THEN
DO I = 1, NTENS
STRESS(I)=ZERO
DO J = 1, NTENS
STRESS(I)=STRESS(I)+CDFULL(I,J) * STRANT(J)
END DO
END DO
C
C INITIALIZE THE 3X3 CONDENSED STIFFNESS MATRIX IF PLANE STRESS CASE
C
ELSE IF ( NDI .EQ. 2) THEN
DO I = 1, NTENS
DO J = 1, NTENS
CDTHREE(I,J)=ZERO
END DO
END DO
C
C
C CONDENSE THE UNDAMAGED STIFFNESS MATRIX
C
CALL MatrixCondense(CDFULL,CDTHREE)
C
C UPDATE THE STRESS
C
DO I = 1, NTENS
STRESS(I)=ZERO
DO J = 1, NTENS
STRESS(I)=STRESS(I)+CDTHREE(I,J) * STRANT(J)
END DO
END DO
END IF
RETURN
END
C******************************************************************************
C TO CHECK THE FAILURE INITIATION AND THE CORRESPONDING DERIVATIVE*********
C******************************************************************************
SUBROUTINE CheckFailureIni(EPITL,EPICL,EPITT,EPICT,EPISLT,STRANT,
1 GFMAT,GFFIB, CELENT, CFULL, DF, DM, DDFDE, DDMDE, NTENS,
2 DFOLD,DMOLD, NDI )
INCLUDE 'ABA_PARAM.INC'
DIMENSION DDFDE(6), DDMDE(6), STRANT(6), CFULL(6,6)
DIMENSION DFMNDE(6), DFFNDE(6)
PARAMETER (ZERO = 0.D0, ONE = 1.D0, TWO = 2.D0, HALF = 0.5D0)
C
C CHECK THE INITIATION CONDITION FOR MATRIX
C FMN=FM/EPITT > 1 THEN EVALUATE THE DAMAGE VARIABLE AND DERIVATIVE
C
TERM1 = STRANT(2)**TWO / EPICT / EPITT
TERM2 = (EPICT - EPITT) / EPICT / EPITT * STRANT(2)
IF (NDI .EQ. 3) THEN
TERM3 = (STRANT(4))**TWO / EPISLT**TWO
ELSE IF (NDI .EQ. 2) THEN
TERM3 = (STRANT(3))**TWO / EPISLT**TWO
END IF
TERM = TERM1 + TERM2 + TERM3
IF (TERM .GT. ZERO) THEN
FMN = SQRT(TERM)
ELSE
FMN = ZERO
END IF
C
C INITIALIZE THE ARRAY AND VARIABLE
C
DM = ZERO
DDMDFMN = ZERO
DO I = 1, 6
DFMNDE(I) = ZERO
DDMDE(I) = ZERO
END DO
IF (FMN .GT. ONE) THEN
C CALCULATE DM, DDMDFMN
CALL DamageEvaluation( CFULL(2,2), FMN, GFMAT, CELENT,
1 EPITT, DM, DDMDFMN)
C CALCULATE DFMNDE
IF (DM .GT. DMOLD) THEN
DFMNDE(2) = HALF / FMN * (TWO * STRANT(2) + EPICT - EPITT)
1 / EPICT / EPITT
IF (NDI .EQ. 3) THEN
DFMNDE(4) = ONE / FMN * STRANT(4) / EPISLT**TWO
ELSE IF (NDI .EQ. 2) THEN
DFMNDE(4) = ONE / FMN * STRANT(3) / EPISLT**TWO
END IF
DO I = 1, 6
DDMDE(I) = DFMNDE(I) * DDMDFMN
END DO
END IF
END IF
DM = MAX (DM, DMOLD)
C
C CHECK THE INITIATION CONDITION FOR FIBER
C FFN=FF/EPITL>1 THEN CALCULATE THE DAMAGE VARIABLE AND DERIVATIVE
C
TERM1 = STRANT(1)**TWO / EPICL / EPITL
TERM2 = (EPICL - EPITL) / EPICL / EPITL *STRANT(1)
TERM = TERM1 + TERM2
IF (TERM .GT. ZERO) THEN
FFN = SQRT(TERM)
ELSE
FFN = ZERO
END IF
DF = ZERO
DDFDFFN = ZERO
DO I = 1, 6
DFFNDE(I) = ZERO
DDFDE(I) = ZERO
END DO
IF (FFN .GT. ONE) THEN
C CALCULATE DF, DDFDFFN
CALL DamageEvaluation( CFULL(1,1), FFN, GFFIB, CELENT,
1 EPITL, DF, DDFDFFN)
C CALCULATE DFFNDE
IF (DF .GT. DFOLD) THEN
DFFNDE(1) = HALF / FFN * (TWO * STRANT(1) + EPICL - EPITL)
1 / EPICL / EPITL
DDFDE(1) = DFFNDE(1) * DDFDFFN
END IF
END IF
DF = MAX (DF, DFOLD)
RETURN
END
C******************************************************************************
C SUBROUTINE TO EVALUATE THE DAMAGE AND THE
c DERIVATIVE************************
C******************************************************************************
SUBROUTINE DamageEvaluation(STIFF, FN, GF, CELENT, EPIT, D,
1 DDDFN)
C CALCULATE DAMAGE VARIABLE
INCLUDE 'ABA_PARAM.INC'
PARAMETER (ONE = 1.D0,tol=1d-3, zero = 0.d0)
TERM1 = STIFF * EPIT**2 * CELENT / GF
TERM2 = (ONE - FN) * TERM1
D = ONE - EXP(TERM2) / FN
C CALCULATE THE DERIVATIVE OF DAMAGE VARIABLE WITH RESPECT TO FAILURE
C RITERION
DDDFN = (ONE / FN + TERM1) * (ONE - D)
RETURN
END
C******************************************************************************
C SUBROUTINE TO CONDENSE THE 4X4 MATRIX INTO 3X3 MATRIX********************
C******************************************************************************
SUBROUTINE MatrixCondense(CFULL,CTHREE)
INCLUDE 'ABA_PARAM.INC'
DIMENSION CFULL(6,6),CTHREE(3,3)
C
CTHREE(1,1) = CFULL(1,1) - CFULL(1,3) * CFULL(3,1) / CFULL(3,3)
CTHREE(1,2) = CFULL(1,2) - CFULL(1,3) * CFULL(3,2) / CFULL(3,3)
CTHREE(2,1) = CFULL(2,1) - CFULL(2,3) * CFULL(3,1) / CFULL(3,3)
CTHREE(2,2) = CFULL(2,2) - CFULL(2,3) * CFULL(3,2) / CFULL(3,3)
CTHREE(3,3) = CFULL(4,4)
RETURN
END
C*******************************************************************************
C SUBROUTINE TO GET THE DERIVATIVE MATRIX OF CONDENSE DAMAGED MATRIX OVER
C**** THE DAMAGE VARIABLE******************************************************
C*******************************************************************************
SUBROUTINE ElasticDerivative(CFULL,DMV,DFV, DCDDM,DCDDF)
INCLUDE 'ABA_PARAM.INC'
DIMENSION CFULL(6,6), DCDDM(6,6),
1 DCDDF(6,6)
PARAMETER (ZERO = 0.D0, ONE = 1.D0, HALF = 0.5D0)
C initialize the data to zero
DO I = 1, 6
DO J = 1, 6
DCDDM(I,J) = ZERO
DCDDF(I,J) = ZERO
END DO
END DO
C
C CALCULATE DC/DDF
C
DCDDF(1,1) = - CFULL(1,1)
DCDDF(1,2) = - (ONE - DMV) * CFULL(1,2)
DCDDF(2,1) = DCDDF(1,2)
DCDDF(1,3) = -CFULL(1,3)
DCDDF(3,1) = DCDDF(1,3)
DCDDF(4,4) = -(ONE - DMV) * CFULL(4,4)
C
C CALCULATE DC/DDM
C
DCDDM(1,2) = - (ONE - DFV) * CFULL(1,2)
DCDDM(2,1) = DCDDM(1,2)
DCDDM(2,2) = -CFULL(2,2)
DCDDM(2,3) = -CFULL(2,3)
DCDDM(3,2) = DCDDM(2,3)
DCDDM(4,4) = -(ONE - DFV) * CFULL(4,4)
RETURN
END
