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) double precision E,XNUE,XDENSITY,XM,XMEANFS,XEFFVOL real V,B,RKO,RKT,RKTH,RKF,RK,ALAMBDA,AMU integer seed,mflag seed=654321 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 STRN(I,J)=ZERO END DO END DO NTENS=ndir+nshr DO I=1,NTENS DO J=1,NTENS AJAC(I,J)=ZERO XC(I,J)=ZERO XS(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-defect density corresponding to each prin. dir if (stepTime.gt.0.0) then print*, 'Failure stress stateNew, stateOld' print*, stateNew(i,4), stateold(i,4) endif C Intialisation of state variables if (stateOld(i,8).lt.0) then stateOld(i,8)=0 endif do j=1,M XDAM(i,j)=zero XPRIN(i,j)=zero XDEF(i,j)=zero enddo mflag=stateOld(i,8) do 100 i = 1,nblock if (mflag.gt.5) then XDAM(i,1)=stateOld(i,1) XDAM(i,2)=stateOld(i,2) XDAM(i,3)=stateOld(i,3) XFAIL(i)=stateOld(i,4) XDEF(i,1)=stateOld(i,5) XDEF(i,2)=stateOld(i,6) XDEF(i,3)=stateOld(i,7) endif if (stepTime.eq.zero) then B=((XEFFVOL)*((XMEANFS)**XM)) PK=ran(seed) Q=ONE/(ONE-PK) XFAIL(i)=(XMEANFS)*((LOG(Q))**(ONE/XM)) seed=seed+1 mflag=10 endif 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 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-XDAM(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-XDAM(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-XDAM(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-XDAM(i,2))**(XALPHA)) 1 /((ONE-XDAM(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-XDAM(i,3))**(XALPHA)) 2 /((ONE-XDAM(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-XDAM(i,1))**XALPHA) 3 /((ONE-XDAM(i,2))**XALPHA)/E C FINDING XS=INVERSE OF COMPLIANCE TENSOR BY GUASS JORDAN METHOD CALL INVERSEA(XC,XS) DO j=1,NTENS stressNew(i,j)=stressOld(i,j)+DOT_PRODUCT(XS(j,:),strainInc(i,:)) ENDDO 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 C FINDING THE PRINCIPAL STRESS RATE DO j=1,3 XPRINR(i,j)=(XPRINN(i,j)-XPRIN(i,j))/dt ENDDO 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 XDAM(i,j)=XDAM(i,j)+zero ELSE IF ((XPRIN(i,j).GT.ZERO).AND.(XPRINR(i,j).GT.ZERO))THEN IF (XPRIN(i,j).LE.XFAIL(i)) THEN XDEF(i,j)=ZERO XDAM(i,j)=ZERO mflag=10 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 XDAM(i,j)=XDAM(i,j)+((dt/6)*RK) mflag=10 ENDIF ELSE XDEF(i,j)=ZERO XDAM(i,j)=ZERO mflag=10 ENDIF ENDIF ENDDO print*, 'Damage values' print*, XDAM(i,1) print*, XDAM(i,2) print*, XDAM(i,3) if (XDAM(i,1).gt.1.0) then XDAM(i,1)=0.99 endif if (XDAM(i,2).gt.1.0) then XDAM(i,2)=0.99 endif if (XDAM(i,3).gt.1.0) then XDAM(i,3)=0.99 endif C STATE VARIABLES:1 TO 3 ARE DAMAGE VALUE VARYING FROM 0 TO 1 C STATE VARIABLES:4-FAILURE STRESS stateNew(i,1)=XDAM(i,1) stateNew(i,2)=XDAM(i,2) stateNew(i,3)=XDAM(i,3) stateNew(i,4)=XFAIL(i) stateNew(i,5)=XDEF(i,1) stateNew(i,6)=XDEF(i,2) stateNew(i,7)=XDEF(i,3) stateNew(i,8)=mflag print*, 'step time,total time=t' print*, stepTime,totalTime print*, 'Failure stress stateNew,Xfail at step t=t' print*, stateNew(i,4),XFAIL(i) if (stateNew(i,1).lt.10.0**-6) then stateNew(i,1)=zero endif if (stateNew(i,2).lt.10.0**-6) then stateNew(i,2)=zero endif if (stateNew(i,3).lt.10.0**-6) then stateNew(i,3)=zero endif 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 RETURN END