Hi,
I am trying to implement probabilistic damage model of the dynamic fragmentation of brittle materials,Hild,Advanaces in applied mechanics,2010
I have got excessive distortion error.Please anyone suggest me how to rectify the error.Please see my subroutine if any error correct me.my model includes anisotropic damage(D1,D2,D3 in principal directions which are found in each time step,alpha is known). Please check the attached cae,vumat code files
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
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)*XDEF(i,j)
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
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)
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
RETURN
END
Did you try it with a single-element simulation model and what are your boudary and load conditions?
The problem is impact of bullet on ceramic plate.The quarter of a plate is the model. I attach a jpeg file of the model boundary conditions and velocity applied.
As Morasch already suggested, you should first check your VUMAT with single element simulations. That is the easiest and probably the only way to debug these scripts. In single element simulation keep track of variables which blow up and cause excessive distortions. This should help in debugging your script.
- Badges
- Science topic