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
