subroutine vumat( 1 nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, 2 stepTime, totalTime, dt, cmname, coordMp, charLength, 3 props, density, strainInc, relSpinInc, 4 tempOld, stretchOld, defgradOld, fieldOld, 3 stressOld, stateOld, enerInternOld, enerInelasOld, 6 tempNew, stretchNew, defgradNew, fieldNew, 5 stressNew, stateNew, enerInternNew, enerInelasNew ) C INCLUDE 'VABA_PARAM.INC' C C C All arrays dimensioned by (*) are not used in this algorithm dimension props(nprops), density(nblock), 1 coordMp(nblock,*), 2 charLength(*), strainInc(nblock,ndir+nshr), 3 relSpinInc(*), tempOld(*), 4 stretchOld(*), defgradOld(*), 5 fieldOld(*), stressOld(nblock,ndir+nshr), 6 stateOld(nblock,nstatev), enerInternOld(nblock), 7 enerInelasOld(nblock), tempNew(*), 8 stretchNew(*), defgradNew(*), fieldNew(*), 9 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev), 1 enerInternNew(nblock), enerInelasNew(nblock) C character*80 cmname C C *********************************************************************************** ** VUMAT,PROBABILISTIC DAMAGE MODEL FOR ABAQUS/STANDARD INCORPORATING ** ** ELASTIC PROPERTIES AND WEIBULL PARAMETERS ** ** ** *********************************************************************************** *********************************************************************************** ** ** ** *USER SUBROUTINE PARAMETER (M=3,N=3,ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0, + SIX=6.D0, NINE=9.D0,XALPHA=0.31,S=4.18) C DIMENSION XIDEN(M,N),XC(6,6),XS(6,6), + STR(M,N),XDAM(nblock,M),AJAC(6,6),STRN(M,N), + XPRIN(nblock,M),XDEF(nblock,M),XPRINN(nblock,M), + XFAIL(nblock),XPRINR(nblock,M) dimension eigVal(nblock,M), eigVec(nblock,M) real V,B,RKO,RKT,RKTH,RKF,RK integer seed seed=7654321 C C XDAM=DAMAGE VALUE CORRESPONDING TO EACH PRICIPAL DIRECTION,XPRIN=PRINCIPAL STRESS C XDEF=MODIFIED GROWTH DEFECT DENSITY CORRESPONDING TO EACH PRINCIPAL STRESS C XPRINN-PRINCIPALSTRESS NEW,XPRIN-PRICIPAL STRESS OLD C XPRINR-PRINCIPAL STRESS RATE C C-------------------------------------------------------------------- C C SPECIFY MATERIAL PROPERTIES C ALL DIMENSIONS ARE IN meter C E=props(1) XNUE=props(2) XDENSITY=props(3) XPOROSITY=props(4) XM=props(5) XMEANFS=props(6) XEFFVOL=props(7) ALAMBDA=E*XNUE/(ONE+XNUE)/(ONE-TWO*XNUE) AMU=E/(ONE+XNUE)/2 V=SQRT(E/XDENSITY) DO I=1,M DO J=1,N STR(I,J)=ZERO END DO END DO NTENS=ndir+nshr DO I=1,NTENS DO J=1,NTENS AJAC(I,J)=ZERO ENDDO ENDDO C STATE VARIABLES:1 TO 3 ARE DAMAGE VALUE VARYING FROM 0 TO 1 C STATE VARIABLES:4-FAILURE STRESS C STATE VARIABLES:5,6,7-PRINCIPAL STRESSES C STATE VARIABLES:8,9,10-DEFECT DENSITY CORRESPONDING TO EACH PRINCIPAL STRESS if (stepTime.eq.zero) then do 101 i = 1,nblock B=((XEFFVOL)*((XMEANFS)**XM)) PK=ran(seed) Q=ONE/(ONE-PK) XFAIL(i)=(XMEANFS)*((LOG(Q))**(ONE/XM)) AJAC(1,1)= (ALAMBDA+TWO*AMU) AJAC(2,2)= (ALAMBDA+TWO*AMU) AJAC(3,3)= (ALAMBDA+TWO*AMU) AJAC(4,4)= AMU AJAC(5,5)= AMU AJAC(6,6)= AMU AJAC(1,2)= ALAMBDA AJAC(1,3)= ALAMBDA AJAC(2,1)= ALAMBDA AJAC(2,3)= ALAMBDA AJAC(3,1)= ALAMBDA AJAC(3,2)= ALAMBDA print*, 'Failure stress' print*, XFAIL(i) DO j=1,NTENS stressNew(i,j)=stressOld(i,j)+ 2 (DOT_PRODUCT(AJAC(j,:), strainInc(i,:))) ENDDO stateNew(i,4)=XFAIL(i) seed=seed+1 101 continue endif do 100 i = 1,nblock C C WRITE STRESSES IN MATRIX FORM C DO j = 1,3 STR(j,j) = stressOld(i,j) END DO STR(1,2) = stressOld(i,4) STR(2,1) = stressOld(i,4) STR(1,3) = stressOld(i,5) STR(3,1) = stressOld(i,5) STR(2,3) = stressOld(i,6) STR(3,2) = stressOld(i,6) c FIND THE EIGEN VALUES OF THE STRESSOLD MATRIX CALL VSPRINC( NBLOCK, STR, eigVal, ndir, nshr ) C EVR IS THE PRICIPAL STRESSES DO j=1,3 XPRIN(i,j)=eigVal(i,j) ENDDO print*, 'old principal stress' print*, XPRIN(i,:) C C WRITE STRESSES IN MATRIX FORM C DO j = 1,3 STRN(j,j) = stressNew(i,j) END DO STRN(1,2) = stressNew(i,4) STRN(2,1) = stressNew(i,4) STRN(1,3) = stressNew(i,5) STRN(3,1) = stressNew(i,5) STRN(2,3) = stressNew(i,6) STRN(3,2) = stressNew(i,6) c FIND THE EIGEN VALUES OF THE STRESSNEW MATRIX CALL VSPRINC( NBLOCK, STRN, eigVal, ndir, nshr ) C EVR IS THE PRICIPAL STRESSES DO j=1,3 XPRINN(i,j)=eigVal(i,j) ENDDO print*, 'new principal stress' print*, XPRINN(i,:) C FINDING THE PRINCIPAL STRESS RATE DO j=1,3 XPRINR(i,j)=(XPRINN(i,j)-XPRIN(i,j))/dt ENDDO print*, 'principal stress rate' print*, XPRINR(i,:) C FINDING THE DEFECT DENSITY CORRESPONDING TO EACH PRICIPAL STRESS C FINDING THE DAMAGE VALUE CORRESPONDING TO EACH PRINCIPAL STRESS C REFER EQUATION 83 DO j=1,3 IF (XPRINR(i,j).LE.ZERO) THEN stateNew(i,j)=stateOld(i,j) ELSE IF ((XPRIN(i,j).GT.ZERO).AND.(XPRINR(i,j).GT.ZERO))THEN IF (XPRIN(i,j).LE.XFAIL(i)) THEN stateNew(i,j+7)=ZERO stateNew(i,j)=ZERO ELSE stateNew(i,j+7)=K*((XPRIN(i,j))**XM) C XBETA IS THE VALUE OF THE RIGHT HAND SIDE IN THE EQUATION 81 XBETA=6*S*((0.33*V)**3)*stateNew(i,j+7) RKO=((1-XDAM(i,j))*(XBETA*(totalTime**2)))/2 RKT=((1-(XDAM(i,j)+((dt*RKO)/2))) 1 *(XBETA*((totalTime+dt/2)**2))/2) RKTH=(1-(XDAM(i,j)+(dt*RKT)/2)) 1 *(XBETA*((totalTime+dt/2)**2))/2 RKF=(1-(XDAM(i,j)+(dt*RKTH))) 1 *(XBETA*((totalTime+dt)**2))/2 RK=RKO+2*RKT+2*RKTH+RKF stateNew(i,j)=stateOld(i,j)+((dt/6)*RK) ENDIF ELSE stateNew(i,j+7)=ZERO stateNew(i,j)=ZERO ENDIF ENDIF ENDDO print*, 'Damage values' print*, stateNew(i,1) print*, stateNew(i,2) print*, stateNew(i,3) C STATE VARIABLES:1 TO 3 ARE DAMAGE VALUE VARYING FROM 0 TO 1 C STATE VARIABLES:4-FAILURE STRESS stateNew(i,4)=XFAIL(i) stateNew(i,5)=XPRIN(i,1) stateNew(i,6)=XPRIN(i,2) stateNew(i,7)=XPRIN(i,3) C UPDATING THE STRESS MATRIX:TO DO THAT UPDATE THE COMPLIANCE TENSOR C MATRIX AND INVERSE THE COMPLIANCE TENSOR TO OBTAIN THE CONSTITUTIVE TENSOR C XC=COMPLIANCE TENSSOR,XS=INVERSE OF COMPLIANCE TENSOR XC(1,1)=ONE/E/(1-stateNew(i,1)) XC(1,2)=(-XNUE)/E XC(1,3)=(-XNUE)/E XC(1,4)=0.0 XC(1,5)=0.0 XC(1,6)=0.0 XC(2,1)=(-XNUE)/E XC(2,2)=ONE/E/(1-stateNew(i,2)) XC(2,3)=(-XNUE)/E XC(2,4)=0.0 XC(2,5)=0.0 XC(2,6)=0.0 XC(3,1)=(-XNUE)/E XC(3,2)=(-XNUE)/E XC(3,3)=ONE/E/(1-stateNew(i,3)) XC(3,4)=0.0 XC(3,5)=0.0 XC(3,6)=0.0 XC(4,1)=0.0 XC(4,2)=0.0 XC(4,3)=0.0 XC(4,4)=(ONE+XNUE)/((ONE-stateNew(i,2))**(XALPHA)) 3 /((ONE-stateNew(i,3))**(XALPHA))/E XC(4,5)=0.0 XC(4,6)=0.0 XC(5,1)=0.0 XC(5,2)=0.0 XC(5,3)=0.0 XC(5,4)=0.0 XC(5,5)=(ONE+XNUE)/((ONE-stateNew(i,3))**(XALPHA)) 4 /((ONE-stateNew(i,1))**(XALPHA))/E XC(5,6)=0.0 XC(6,1)=0.0 XC(6,2)=0.0 XC(6,3)=0.0 XC(6,4)=0.0 XC(6,5)=0.0 XC(6,6)=(ONE+XNUE)/((ONE-stateNew(i,1))**XALPHA) 5 /((ONE-stateNew(i,2))**XALPHA)/E C FINDING XS=INVERSE OF COMPLIANCE TENSOR BY GUASS JORDAN METHOD CALL INVERSEA(XC,XS) print*, 'inverse of compliance tensor' print*, XS(:,:) DO j=1,NTENS stressNew(i,j)=stressOld(i,j)+DOT_PRODUCT(XS(j,:), strainInc(i,:)) ENDDO 100 continue return end ********************END OF VUMAT SUBROUTINE****************************** ************************************************************************* C ** INVERSE OF A 6X6 MATRIX BY GAUSS JORDAN METHOD ************************************************************************* C ** INVERSE OF A 6X6 MATRIX BY GAUSS JORDAN METHOD SUBROUTINE INVERSEA(A,C) C INCLUDE 'VABA_PARAM.INC' C PARAMETER (M=6,N=6) DIMENSION A(M,N),B(6,12),C(M,N) INTEGER I,J,N C ** BUILDING AUGUMENTED MATRIX B DO 40 J=1,N DO 30 I=1,N B(I,J)=A(I,J) 30 CONTINUE 40 CONTINUE DO I=1,N DO J=N+1,2*N IF ((J-I).EQ.N) THEN B(I,J)=1.0 ELSE B(I,J)=0.0 ENDIF ENDDO ENDDO DO 160 I=1,N B(I,:)=B(I,:)/B(I,I) DO 150 J=1,N IF (I.NE.J) THEN B(J,:)=B(J,:)-B(J,I)*B(I,:) ENDIF 150 CONTINUE 160 CONTINUE DO I=1,N DO J=N+1,2*N C(I,J-6)=B(I,J) ENDDO ENDDO PRINT*,'C MATRIX ROW 1' PRINT*, C(1,:) RETURN END