Hi all, let me start by saying that I wish I'd found this forum a month or 3 earlier.
For my graduation at the Erasmus MC (Netherlands) I'm researching the role of fibre dispersion and orientation in the damage initiation and propagation of bio-engineered tissue (before moving to actual tissue one day). My predecessor made a UMAT which incorporates the HGO model. From experiments in our group I was able to find the local fibre orientation and dispersion factor.
Now to calculate/represent the damage propagation in the modelled tissue, I added a damage factor to the model, which said that the element should be deleted once damage reached 1. However, as soon as the first element is deleted, the simulation diverges, gets negative eigenvalues, and stops working. I guess that is because when the element gets damaged, it is unable to take stresses, leading to no stiffness, leading to the infamous and dreaded buckling/ negative eigenvalues. From what I've read, my conclusion now is that only an Explicit simulation with a VUMAT will allow for the deletion of elements during a simulation.
Now the options I consider are the following:
1. Use partitions in my CAE, and give local material properties to certain tiles (including fibre dispersion and orientation) using the built in HGO model and add the damage to the model using VUSDFLD. However the damage is based on the strain energy density, which I don't think can be found as an output using built in material models.
2. Convert my UMAT into VUMAT. This seems the better option as it leads to no black boxes.
My question (finally) is:
1. Could anyone help out with converting my UMAT to a VUMAT?
2. On a more helicopter view, do you think this method of simulating damage propagation in fibrous tissue seems logical?Â
I've added the code below, so hopefully we don't have to reinvent the wheel here together.
C 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),DDSDDT(NTENS),DRPLDE(NTENS), 2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1), 3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3) C C ---------------------------------------------------------------- C LOCAL ARRAYS C ---------------------------------------------------------------- C DIMENSION XKIRCH(3,3), DUMSTRSS(3,3), XAMAT(3,2), XI4(2), R(3,3), 1 U(3,3), RT(3,3), XADUM(3,1), XADUM1(3,1) C C PARAMETER(ZERO=0.D0, ONE=1.D0, TWO=2.D0, THREE=3.D0, FOUR=4.D0, 1 SIX=6.D0) C C ---------------------------------------------------------------- C NUMBER OF FIBRE FAMILIES C ---------------------------------------------------------------- C NANISO = 1 C C ---------------------------------------------------------------- C MATERIAL PROPERTIES C ---------------------------------------------------------------- C10 = PROPS(1) D1 = PROPS(2) xk1 = PROPS(3) xk2 = PROPS(4) xkap = PROPS(5) THETAD1 = PROPS(6) THETAD2 = PROPS(7) T1 = PROPS(8) T2 = PROPS(9) C Convert degrees to radians XPI = 3.14159265359 THETAR1 = THETAD1*XPI/180.0 THETAR2 = THETAD2*XPI/180.0 C C Fibre vectors in the reference configuration XAMAT(1,1)=COS(thetar1); XAMAT(2,1)= SIN(thetar1); XAMAT(3,1)=0. XAMAT(1,2)=COS(thetar2); XAMAT(2,2)= SIN(thetar2); XAMAT(3,2)=0. C C This is related to the calculation of the numberical Jacobian C 0 indicates that the stress calculation is not part of the C stress calculation. iter=0 C C--------------------------------------------------------------------------- C POLAR DECOMPOSITION C--------------------------------------------------------------------------- call kpolarDecomp(DFGRD1, U, R) C Calculate the transpose of the rotation matrix RT = transpose(R) call kprinter(R,3,3) C C----------------------------------------------------------------- C CALCULATE LOCAL FIBRE VECTOR IN THE REFERENCE CONFIG. C----------------------------------------------------------------- C DO IAN = 1,NANISO C Pick out fibre vector XADUM(:,1) = XAMAT(:,IAN) C C Update to local basis system XADUM1 = MATMUL(RT,XADUM) C C Store result XAMAT(:,IAN) = XADUM1(:,1) END DO C C----------------------------------------------------------------- C ZERO THE TANGENT MATRIX C----------------------------------------------------------------- C DO I=1,NTENS DO J=1,NTENS DDSDDE(I,J) = ZERO END DO END DO C C----------------------------------------------------------------- C CALCULATE THE KIRCHHOFF STRESS C----------------------------------------------------------------- C call kstress_calc(DFGRD1, C10, D1, xkap, xk1, xk2, T1,T2, XKIRCH, XJ, 1 XAMAT, XI1, XI4, NPT, NANISO, iter, XWE, XWK, RFACTOR,TK, EA, D2) C STATEV(1) = XI1 STATEV(2) = XJ STATEV(3) = EA STATEV(4) = XI4(1) STATEV(5) = XI4(2) STATEV(6) = XWE STATEV(7) = XWK STATEV(8) = RFACTOR STATEV(9) = TK STATEV(10) = D2 C C----------------------------------------------------------------- C CONVERT KIRCHHOFF STRESS TO CAUCHY STRESS C----------------------------------------------------------------- DUMSTRSS = XKIRCH / XJ C C Convert the stress tensor to a Voigt vector call kmatrix2vector(DUMSTRSS, STRESS, nshr) C C----------------------------------------------------------------- C CALCULATE THE TANGENT MATRIX C----------------------------------------------------------------- C call kCTM(XKIRCH,DFGRD1,NTENS,PROPS,XAMAT,NPT,ITER, + NANISO,DDSDDE,NPROPS,XJ,NSHR) C RETURN C C ENDSUBROUTINE UMAT C C----------------------------------------------------------------------------- C SUBROUTINES C----------------------------------------------------------------------------- C SUBROUTINE kstress_calc(DGRAD, C10, D1, xkap, xk1, xk2, T1,T2, 1 XKIRCH, XJ, XAMAT, XI1, XI4, NPT, NANISO, iter, XWE, XWK, RFACTOR,TK,EA,D2) C INCLUDE 'ABA_PARAM.INC' C INTENT(IN) :: DGRAD, C10, D1, xkap, xk1,xk2,XAMAT,NANISO,iter,T1,T2 INTENT(OUT):: XKIRCH, XJ, XI1, XI4, XWE, XWK, RFACTOR, TK, EA, D2 C DIMENSION DGRAD(3,3), BMAT(3,3), XKIRCH(3,3), XAMAT(3,2), 1 B2MAT(3,3), XA(3), a_curnt(3), aoa(3,3), XANISOK(3,3), 2 XANISOTOT(3,3), XI4(2) C PARAMETER(ZERO=0.D0, ONE=1.D0, TWO=2.D0, THREE=3.D0) C C LOCAL ARRAYS C ---------------------------------------------------------------- C C Determinant of the deformation gradient XJ=DGRAD(1, 1)*DGRAD(2, 2)*DGRAD(3, 3) 1 -DGRAD(1, 2)*DGRAD(2, 1)*DGRAD(3, 3) 2 +DGRAD(1, 2)*DGRAD(2, 3)*DGRAD(3, 1) 3 +DGRAD(1, 3)*DGRAD(3, 2)*DGRAD(2, 1) 4 -DGRAD(1, 3)*DGRAD(3, 1)*DGRAD(2, 2) 5 -DGRAD(2, 3)*DGRAD(3, 2)*DGRAD(1, 1) C C C Calculate left Cauchy-Green deformation tensor BMAT = MATMUL(DGRAD,TRANSPOSE(DGRAD)) C C----------------------------------------------------------------------- C CALCULCATE THE INVARIANTS C----------------------------------------------------------------------- C C I1 invariant XI1 = BMAT(1,1)+BMAT(2,2)+BMAT(3,3) C C Square of the B matrix B2MAT = MATMUL(BMAT,BMAT) C C-------------------------------------------------------------------- C CALCULATE THE ANISOTROPIC PORTION OF THE KIRCH STRESS C-------------------------------------------------------------------- C C Zero the matrix DO I=1,3 DO J=1,3 XANISOTOT(I,J) = ZERO END DO END DO C C-------------------------------------------------------------------- C Loop over the number of fibre families C-------------------------------------------------------------------- DO IANISO = 1,NANISO C C Fibre vector in reference configuration XA = XAMAT(1:3,IANISO); C C Calculate fibre vector in current configuration: a = F*A a_curnt = MATMUL(DGRAD,XA) C C Calculate structural tensor: a \otimes a DO i=1,3 DO j=1,3 aoa(i,j) = a_curnt(i) * a_curnt(j) ENDDO ENDDO C C Calculate the fibre invariant XI4(IANISO) = aoa(1,1) + aoa(2,2) + aoa(3,3) C EA = xkap*(XI1 - 3.D0) + (1.D0 - 3.D0*xkap)*(XI4(IANISO)-1.D0) EAA = xk2*(EA)**2 C XWE = C10*(XI1 - 3.D0) 1 + (xk1/2.D0*xk2)*((exp(EAA)) - 1.D0) 2 + (1.D0/D1)*((XJ**TWO - 1.D0)/2.D0 - LOG(XJ)) C TK = SQRT(TWO * XWE) C IF (TK .LT. T1) THEN DAM = ZERO ELSE IF (TK .LE. T2) THEN DAM = ((TK - T1) / (T2 - T1))**TWO ELSE DAM = ONE END IF C XWK = (1.D0 - DAM)*(C10*(XI1 - 3.D0) 1 + (xk1/2.D0*xk2)*((exp(EAA)) - 1.D0)) 2 + (1.D0/D1)*((XJ**TWO - 1.D0)/2.D0 - LOG(XJ)) C RFACTOR = ONE - DAM C If damage is higher than 0.2, delete element. C IF (DAM .LE. 0.2) THEN D2= 1 ELSE D2=0 END IF C CALCULATE THE ISOTROPIC PORTION OF THE KIRCHHOFF STRESS C C PART 1 COEFF1 = TWO*(C10 + exp(EAA)*xkap*EA*xk1)/(XJ**(TWO/THREE)) C C KIRCHHOFF STRESS PART 2 C TRBMAT= COEFF1*(BMAT(1,1)+BMAT(2,2)+BMAT(3,3))/THREE C C KIRCHHOFF STRESS PART 1 C XKIRCH = COEFF1 * BMAT C C SUBTRACT THE PART 2 C DO I = 1,3 XKIRCH(I,I) =XKIRCH(I,I) - TRBMAT END DO XKIRCH = XKIRCH * RFACTOR C C FORM VOLUMETRIC PART, 3 C COEFF3 = 2*(XJ-ONE)*XJ/D1 C C ADD TO THE PREVIOUS PARTS C DO I = 1,3 XKIRCH(I,I) = XKIRCH(I,I) + COEFF3 END DO C C CALCULATE THE ANISOTROPIC PORTION OF THE KIRCHHOFF STRESS C C Calculate the various parts of the equation for kirch stress p1 = 2.0*xk1; C p2 = EA*(1.D0 - 3.D0*xkap) C p4 = exp(EAA) C XANISOK = RFACTOR*p1*p2*p4*aoa C C Sum the stress over the fibre families XANISOTOT = XANISOTOT + XANISOK C END DO C C Add the isotropic and anisotropic parts XKIRCH = XKIRCH + XANISOTOT C RETURN END SUBROUTINE kstress_calc C SUBROUTINE kpolarDecomp(G,S,R) c c declarations INCLUDE 'ABA_PARAM.INC' dimension G(3,3),R(3,3),S(3,3),COG(3,3),XP(3,3),ADP(3,3), 1 XPBI(3,3) c c constants EPS=1.0E-5 c c cofactors and determinant of g COG(1,1)=G(2,2)*G(3,3)-G(2,3)*G(3,2) COG(2,1)=G(1,3)*G(3,2)-G(1,2)*G(3,3) COG(3,1)=G(1,2)*G(2,3)-G(1,3)*G(2,2) COG(1,2)=G(2,3)*G(3,1)-G(2,1)*G(3,3) COG(2,2)=G(1,1)*G(3,3)-G(1,3)*G(3,1) COG(3,2)=G(1,3)*G(2,1)-G(1,1)*G(2,3) COG(1,3)=G(2,1)*G(3,2)-G(2,2)*G(3,1) COG(2,3)=G(1,2)*G(3,1)-G(1,1)*G(3,2) COG(3,3)=G(1,1)*G(2,2)-G(1,2)*G(2,1) G3=G(1,1)*COG(1,1)+G(2,1)*COG(2,1)+G(3,1)*COG(3,1) c c P = trans(G) * G = S * S DO 10000 I=1,3 XP(I,1)=G(1,I)*G(1,1)+G(2,I)*G(2,1)+G(3,I)*G(3,1) XP(I,2)=G(1,I)*G(1,2)+G(2,I)*G(2,2)+G(3,I)*G(3,2) XP(I,3)=G(1,I)*G(1,3)+G(2,I)*G(2,3)+G(3,I)*G(3,3) 10000 CONTINUE c c adjoint of P ADP(1,1)=XP(2,2)*XP(3,3)-XP(2,3)*XP(3,2) ADP(2,2)=XP(1,1)*XP(3,3)-XP(1,3)*XP(3,1) ADP(3,3)=XP(1,1)*XP(2,2)-XP(1,2)*XP(2,1) c c G invariants G1SQ=XP(1,1)+XP(2,2)+XP(3,3) G1=SQRT(G1SQ) G2SQ=ADP(1,1)+ADP(2,2)+ADP(3,3) G2=SQRT(G2SQ) c c initialize iteration H1=G2/G1SQ H2=G3*G1/G2SQ X=1.0 Y=1.0 c c iteration loop 10001 CONTINUE DEN=2.0*(X*Y-H1*H2) RES1=1.0-X*X+2.0*H1*Y RES2=1.0-Y*Y+2.0*H2*X DX=(Y*RES1+H1*RES2)/DEN DY=(H2*RES1+X*RES2)/DEN X=X+DX Y=Y+DY IF(ABS(DX/X).GT.EPS.OR.ABS(DY/Y).GT.EPS)GO TO 10001 c c BETA invariants BETA1=X*G1 BETA2=Y*G2 c c invert ( trans(G) * G + BETA2 * identity ) XP(1,1)=XP(1,1)+BETA2 XP(2,2)=XP(2,2)+BETA2 XP(3,3)=XP(3,3)+BETA2 XPBI(1,1)=XP(2,2)*XP(3,3)-XP(2,3)*XP(3,2) XPBI(1,2)=XP(1,3)*XP(3,2)-XP(1,2)*XP(3,3) XPBI(1,3)=XP(1,2)*XP(2,3)-XP(1,3)*XP(2,2) XPBI(2,1)=XP(2,3)*XP(3,1)-XP(2,1)*XP(3,3) XPBI(2,2)=XP(1,1)*XP(3,3)-XP(1,3)*XP(3,1) XPBI(2,3)=XP(1,3)*XP(2,1)-XP(1,1)*XP(2,3) XPBI(3,1)=XP(2,1)*XP(3,2)-XP(2,2)*XP(3,1) XPBI(3,2)=XP(1,2)*XP(3,1)-XP(1,1)*XP(3,2) XPBI(3,3)=XP(1,1)*XP(2,2)-XP(1,2)*XP(2,1) DETPBI=XP(1,1)*XPBI(1,1)+XP(2,1)*XPBI(1,2)+XP(3,1)*XPBI(1,3) c c R = (cofac(G)+BETA1*G) * inv(trans(G)*G+BETA2*identity) DO 10002 I=1,3 R(I,1)=((COG(I,1)+BETA1*G(I,1))*XPBI(1,1) 1 +(COG(I,2)+BETA1*G(I,2))*XPBI(2,1) 2 +(COG(I,3)+BETA1*G(I,3))*XPBI(3,1))/DETPBI R(I,2)=((COG(I,1)+BETA1*G(I,1))*XPBI(1,2) 1 +(COG(I,2)+BETA1*G(I,2))*XPBI(2,2) 2 +(COG(I,3)+BETA1*G(I,3))*XPBI(3,2))/DETPBI R(I,3)=((COG(I,1)+BETA1*G(I,1))*XPBI(1,3) 1 +(COG(I,2)+BETA1*G(I,2))*XPBI(2,3) 2 +(COG(I,3)+BETA1*G(I,3))*XPBI(3,3))/DETPBI 10002 CONTINUE c c S = trans(R) * G DO 10003 I=1,3 S(I,1)=R(1,I)*G(1,1)+R(2,I)*G(2,1)+R(3,I)*G(3,1) S(I,2)=R(1,I)*G(1,2)+R(2,I)*G(2,2)+R(3,I)*G(3,2) S(I,3)=R(1,I)*G(1,3)+R(2,I)*G(2,3)+R(3,I)*G(3,3) 10003 CONTINUE c c done RETURN END C C Subroutine to calculate the tangent matrix C Subroutine to calculate the consistent tangent matrix using the C perturbation method outlined by Miehe and Sun C C subroutine kCTM(XKIRCH1,DFGRD1,NTENS,PROPS,XAMAT,NPT,ITER, + NANISO,DDSDDE,NPROPS,XJ1,NSHR) C INCLUDE 'ABA_PARAM.INC' C DIMENSION XKIRCH1(3,3), DFP(3,3), DFGRD_PERT(3,3),XKIRCH_PERT(3,3) 1 , CMJ(3,3), CMJVEC(NTENS), ilist(6), jlist(6), XAMAT(3,2), XI4(2) 2 , DFGRD1(3,3), DDSDDE(NTENS,NTENS), PROPS(NPROPS) C C C LOCAL ARRAYS C ---------------------------------------------------------------- C XKIRCH1 - True Kirchhoff stress (i.e. not based on perturbed def. grad.) C DFGRD1 - True deformation gradient for this increment C DFP - Increment of the perturbed deformation gradient C DFGRD_PERT - Perturbed def. grad. C XKIRCH_PERT - Perturbed Kirchhoff stress C CMJ - (:,:,I,J) Components of material jacobian C CMJVEC - Above in vector form C ILIST, JLIST - Set of index components to be perturbed upon C DUMSTRSS - DUMMY STRESS TENSOR C XAMAT - array containing reference configuration fibre vectors C XI4 - array storing the fibre invariant (will not be used in this) C DDSDDE - Voigt notation matrix to store the tangent matrix C PROPS - Properties read in from Abaqus .inp file C ---------------------------------------------------------------- C C Perturbation parameter eps = 1.0e-08 C C Material properties C10 = PROPS(1) D1 = PROPS(2) xk1 = PROPS(3) xk2 = PROPS(4) xkap = PROPS(5) THETAD1 = PROPS(6) THETAD2 = PROPS(7) T1 = PROPS(8) T2 = PROPS(9) C C The Kirchhoff stress is perturbed six times over the C i, j'th components of the deformation gradient. This C array lists each of the componets to be perturbed over. C i | j ilist(1) = 1; jlist(1) = 1; ilist(2) = 2; jlist(2) = 2; ilist(3) = 3; jlist(3) = 3; ilist(4) = 1; jlist(4) = 2; ilist(5) = 1; jlist(5) = 3; ilist(6) = 2; jlist(6) = 3; C C Loop over each of the six components of the deformation gradient C which are to be perturbed, each time calculating a new column of C the tangent matrix. Perturbation: do iter = 1,NTENS C C----------------------------------------------------------------- C CREATE THE PERTURBED DEFORMATION GRADIENT C----------------------------------------------------------------- C Pick out the (i,j) component of the deformation gradient C to be perturbed ii = ilist(iter) jj = jlist(iter) C C Calculate the increment of the perturbed deformation gradient call kdelF(ii, jj, DFGRD1, eps, DFP) C C Add to the "true" deformation gradient to create the perturbed C deformation gradient DFGRD_PERT = DFGRD1 + DFP C C----------------------------------------------------------------- C CALCULATE KIRCHHOFF STRESS BASED ON THE PERTURBED DEFORMATION GRADIENT C----------------------------------------------------------------- C N.B. call kstress_calc(DFGRD_PERT, C10, D1, xkap, xk1, xk2, T1,T2, 1 XKIRCH_PERT, XJP, XAMAT, XI1, XI4, NPT, NANISO, 2 iter, XWE, XWK, RFACTOR,TK, EA, D2) C C----------------------------------------------------------------- C C Difference between the perturbed(i,j) and unpert. stress CMJ = XKIRCH_PERT - XKIRCH1 C C Normalize by perturbation parameter and return to Cauchy stress CMJ = CMJ/XJ1/eps C C C----------------------------------------------------------------- C FORM THE TANGENT MATRIX C----------------------------------------------------------------- C C Convert tensor to Voigt vector call kmatrix2vector(CMJ, CMJVEC, NSHR) C C Insert this vector into the iter'th column of DDSDDE do insert = 1,NTENS C DDSDDE(insert,iter) = CMJVEC(insert) C end do C C end do Perturbation end subroutine kCTM C C Helper subroutines C * KDELF - Calculate the increment of the deformation gradient for C a given perturbation in (i,j), with epsilon C C * KPRINTER - Print out a matrix of any size C C * KMTMS - Multiply two matrices C C * KMATRIX2VECTOR - Convert a 3x3 matrix to a 6x1 vector C C * KDOTPROD - Dot product of two vectors C C------------------------------------------------------------ subroutine kdotprod(A, B, dotp, n) INCLUDE 'ABA_PARAM.INC' intent(in) :: A, B, n intent(out):: dotp dimension A(n), B(n) dotp = 0.0 do i = 1,n dotp = dotp + A(i)*B(i) end do end subroutine kdotprod C C------------------------------------------------------- C C------------------------------------------------------- C subroutine kmatrix2vector(XMAT, VEC, NSHR) INCLUDE 'ABA_PARAM.INC' intent(in) :: XMAT, NSHR intent(out):: VEC dimension xmat(3,3), vec(6) do i=1,3 vec(i) = xmat(i,i); end do vec(4) = xmat(1,2); IF (NSHR==3) then vec(5) = xmat(1,3); vec(6) = xmat(2,3); END IF end subroutine kmatrix2vector C C------------------------------------------------------- C C------------------------------------------------------- C SUBROUTINE KMTMS (M, N, L, A, KA, B, KB, C, KC) INCLUDE 'ABA_PARAM.INC' C intent(in) :: M, N, L, A, KA, B, KB, KC intent(out):: C C C C PRODUCT OF REAL MATRICES C DIMENSION A(KA,N), B(KB,L), C(KC,L) DOUBLE PRECISION W C C DO 30 J = 1,L DO 20 I = 1,M W = 0.D0 DO 10 K = 1,N W = W + A(I,K) * B(K,J) 10 CONTINUE C(I,J) = W 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE KMTMS C C------------------------------------------------------- C C------------------------------------------------------- C subroutine kprinter(tens, m, n) INCLUDE 'ABA_PARAM.INC' intent(in):: tens, m, n dimension tens(m,n) write(6,*) do i = 1,m do j = 1,n write(6,'(e19.9)',advance='no'),tens(i,j) end do write(6,*) end do write(6,*) return end subroutine kprinter C C------------------------------------------------------- C C------------------------------------------------------- C subroutine kdelF(m, n, DGRAD, eps, DF) INCLUDE 'ABA_PARAM.INC' intent (in) :: DGRAD, eps, m, n intent (out):: DF C Input: the index's i & j; The current deformation gradient (DGRAD). The perturbation C increment (eps) C C Output: The perturbed increment DF C dimension dyad1(3,3), dyad2(3,3), DGRAD(3,3), DF(3,3), DFp1(3,3) c Zero the dyad matrices c do i = 1,3 do j = 1,3 dyad1(i,j) = zero dyad2(i,j) = zero end do end do c Place the 1's in the correct location dyad1(m,n) = 1.0; c dyad2(n,m) = 1.0; c c KMTMS (M, N, L, A, KA, B, KB, C, KC) call KMTMS(3, 3, 3, dyad1, 3, DGRAD, 3, DFp1, 3) DF = DFp1 call KMTMS(3, 3, 3, dyad2, 3, DGRAD, 3, DFp1, 3) DF = DF + DFp1 DF = 0.5*DF*eps end subroutine kdelF C