UMAT for muscle contraction

[Spring63]

Hi,

I am preparing a code for muscle contraction.I am working on a umat. However when I call the subroutine from ABAQUS, it gave me an error:"Problem during compilation". I appreciate if anyone can technically check my code and let me know where I have been wrong.

Thanks,
Spring


SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,
2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS,
3 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,
4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)

C IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 MATERL
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),
4 DFGRD0(3,3),DFGRD1(3,3)
C
DIMENSION PROPS1(100),PROPS2(100)

DIMENSION GC(3,3),CBAR(3,3),BBAR(3,3),AN0(3),ANN0(3,3),FN(3),
1 AN1(3),ANN1(3,3),UNIT2(3,3),ST(3,3),UNIT4(3,3,3,3),H4(3,3,3,3),
2 H2(9,9),STI(3,3),STFPE(3,3),STFSE(3,3),STJ(3,3)


PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0)
C DATA NEWTON,TOLER/10,1.D-6/
C
C -----------------------------------------------------------
C UMAT
C
C CAN NOT BE USED FOR PLANE STRESS
C -----------------------------------------------------------
C -----------------------------------------------------------

IF (NDI.NE.3) THEN
WRITE(6,1)
1 FORMAT(//,30X,'***ERROR - THIS UMAT MAY ONLY BE USED FOR ',
1 'ELEMENTS WITH THREE DIRECT STRESS COMPONENTS')
ENDIF


EPS=1.0E-6


C PROPS(x) ÊäÈëÄ£ÐͲÎÊý

c martins model
b= PROPS(1) ! 1.79
c= PROPS(2) ! 0.379e-3 MPa
TM0 = PROPS(8) ! 0.22mpa modified according to kojic model
D = PROPS(11) ! 0.001 (MPa)-1

c kojic model
AK= PROPS(3) ! 0.3
ALFA= PROPS(4) ! 10
BETA= PROPS(5) ! 1e-3
DDRSM0=PROPS(6) !2.0 ! strain rate? (m/s) should be (1/s)
AA= PROPS(7) ! 0.07
AB1= PROPS(9) ! 5
AD1= PROPS(10) ! 1.5


c as follow it is from the muscle model of Martins
FX0 = PROPS(12)
FY0 = PROPS(13)
FZ0 = PROPS(14)
c INITIAL VALUES
STATEV(1)=PROPS(15)
STATEV(2)=PROPS(16)
STATEV(3)=PROPS(17)
STATEV(4)=PROPS(18)
c number of the columns for s-lambda and activation fun
NVALUE1=PROPS(19)
NVALUE2=PROPS(20)

DIS = sqrt(FX0**2+FY0**2+FZ0**2)
FX0 = FX0/DIS
FY0 = FY0/DIS
FZ0 = FZ0/DIS

AN0(1)=FX0
AN0(2)=FY0
AN0(3)=FZ0

C21 Jacobian of deformation gradient

DET = DFGRD1(1,1) * DFGRD1(2,2) * DFGRD1(3,3)
2 + DFGRD1(1,2) * DFGRD1(2,3) * DFGRD1(3,1)
3 + DFGRD1(1,3) * DFGRD1(3,2) * DFGRD1(2,1)
4 - DFGRD1(1,2) * DFGRD1(2,1) * DFGRD1(3,3)
5 - DFGRD1(1,3) * DFGRD1(3,1) * DFGRD1(2,2)
6 - DFGRD1(2,3) * DFGRD1(3,2) * DFGRD1(1,1)

SCALE1 = DET**(-ONE /THREE)
SCALE2 = DET**(-TWO /THREE)

C22__the right Cauchy-Green strain tensor--GC(3,3)

DO I = 1,3
DO J = 1,3
GC(I,J)= 0.0
DO K=1,3
GC(I,J)= GC(I,J)+DFGRD1(K,I)*DFGRD1(K,J)
END DO
END DO
END DO

C23_right Cauchy-Green strain tensor with the volume change eliminated CBAR(3,3)

DO I=1,3
DO J=1,3
CBAR(I,J) = GC(I,J)* SCALE2
END DO
END DO

C24_the left Cauchy-Green strain tensor with volume change eliminated BBAR(3,3)

DO I = 1,3
DO J = 1,3
BBAR(I,J)= 0.0
DO K=1,3
BBAR(I,J)= BBAR(I,J)+DFGRD1(I,K)*DFGRD1(J,K)
END DO
END DO
END DO

DO I=1,3
DO J=1,3
BBAR(I,J) =BBAR(I,J)*SCALE2
END DO
END DO

C25_HIBARC1- the first invariant right Cauchy-Green strain tensor
C with the volume change eliminated

TRCBAR=CBAR(1,1)+CBAR(2,2)+CBAR(3,3)
HIBARC1=TRCBAR

C26_ANMDF--fiber stretch ratio in the direction N of undeformed fiber
DO I=1,3
DO J=1,3
ANN0(I,J)=AN0(I)*AN0(J)
END DO
END DO

ANMDF = 0.0
DO I=1,3
DO J=1,3
ANMDF =ANMDF + CBAR(I,J)*ANN0(J,I)
END DO
END DO

ANMDF = SQRT(ANMDF)

C27_AN1(I)--the current muscle fiber direction

DO I=1,3
FN(I) =0.0
DO K=1,3
FN(I) =FN(I)+DFGRD1(I,K)*AN0(K)
END DO
AN1(I) = SCALE1*FN(I)/ANMDF
END DO


C
C the relationship between time and active function
C

DO K=1,2*NVALUE1
PROPS1(K) = PROPS(21+K)
END DO
C
C the relationship between rate of stretch and tension stress of maximum
C

DO K=1,2*NVALUE2
PROPS2(K) = PROPS(51+K)
END DO
C
C giving the initial value of some variables
C
RSS = STATEV(1)
RSM = STATEV(2)
SIGM0 = STATEV(3)
DRSM = STATEV(4)

AAA0=0.0
DO I=1,4
AAA0=AAA0+ABS(STATEV(I))
END DO

IF (AAA0.EQ.0.) THEN
RSS = 1.0
RSP = 1.0 ! the total rate of stretch
RSM = (1.0+AK)*RSP-AK ! the rate of stretch of contractile element
DRSM=0.0
CALL GETCURVE(SIGM0,RSP,PROPS2,NVALUE2)
ENDIF

C
C CALL ACTIVE FUNCTION at time t+dt
C
c ATIME=time(1)+dtime
c ATIME=time(2)+dtime

! if (KSTEP.eq.1) then
! ALFAA=0.0
! end if

! if (KSTEP.eq.2) then
ATIME=time(1)+dtime
CALL GETCURVE(ALFAA,ATIME,PROPS1,NVALUE1)
! end if

c
c giving the stretch of the material point in current time
c
RSP=ANMDF
C
C hill muscle with algorithm presented by kojic
C
IF ((DRSM+1.0).GT.1.001) THEN
AB = AB1
AD = AD1
ELSE
AB = 1.0
AD=0.0
END IF
cccccccccccccccccccccccccccccctest
C only kojic Hill's model
c AB = 1.0
c AD=0.0

C computing the stress of contractile element at time t

EXPARS=EXP(ALFA*( RSS-1.0))
SIGMS= BETA*( EXPARS-1.0)

C
C (1)computing to be solved cofficient A1
C

A1=(1.+AK)*RSP-RSM-AK*RSS
C
C (2)computing to be solved cofficient A2
C
DRSM0=DTIME*DDRSM0
AC= SIGM0/AA
A2=(SIGMS+BETA)*(1.0-AC*AC*A1/DRSM0)
C
C (3)computing to be solved cofficient A3
C
A3=(SIGMS+BETA)*AB*AK*AC/DRSM0
C
C (4)computing to be solved cofficient A4
C
A4=(AB*BETA*AC+SIGM0*ALFAA*(AD*AB*AC+AD-1.0))*AK/DRSM0
C
C (5)computing to be solved cofficient A5
C

A5= BETA+ SIGM0*ALFAA+
1 (SIGM0*ALFAA*(1.0-AD-AD*AB*AC)-AB*BETA*AC)*A1/DRSM0
C
C solve the unknown RSS by standard Newton's method
C
DRSS=0.0
CALL DNEWT(DRSS,EPS,L,A2,A3,A4,A5,ALFA)
C KNUMBER=60-L ! iter number
c
c update of quantities for the next step
c
c update RSS
RSS= RSS+DRSS

c update RSM and computing the velocity in CE

DRSM=A1-AK*DRSS

V=DRSM/DTIME

RSM=RSM+DRSM

c update sigms and sigmm

EXPARS=EXP(ALFA*( RSS-1.0))
SIGMS= BETA*(EXPARS-1.0)

SIGMM= SIGM0* ALFAA*(AD-(AD-1.0)*(1.0+ DRSM/ DRSM0)
1 /(1.0-AB*AC* DRSM/DRSM0))
C
C
C28_
IF(ANMDF.GT.1.0) THEN
FPE = 4.0*(ANMDF-1.0)**2
ELSE
FPE = 0.0
END IF

C2_10_

U1PE=TM0*FPE

C2_11

U1SE = SIGMS ! changed to the active part in kojic model

C2_12
U1F=U1PE+U1SE

C2_13

U1I=b*c*EXP(b*(HIBARC1-3.0))
C2_14
U1J=2.*(DET-1.0)/D

C computer Cauchy stress

DO I=1,3
DO J=1,3
ANN1(I,J)=AN1(I)*AN1(J)
END DO
END DO

DO I=1,3
DO J=1,3
UNIT2(I,J)=0.0
END DO
UNIT2(I,I)=1.0
END DO

DO I=1,3
DO J=1,3
STI(I,J) = (U1I*(2.0*BBAR(I,J)-2.0*HIBARC1*UNIT2(I,J)/3.0))/DET
END DO
END DO

DO I=1,3
DO J=1,3
STFPE(I,J)=(U1PE*(ANMDF*ANN1(I,J)-1.0*ANMDF*UNIT2(I,J)/3.0))/DET
END DO
END DO

DO I=1,3
DO J=1,3
STFSE(I,J)=(U1SE*(ANMDF*ANN1(I,J)-1.0*ANMDF*UNIT2(I,J)/3.0))/DET
END DO
END DO

DO I=1,3
DO J=1,3
STJ(I,J)=U1J*UNIT2(I,J)
END DO
END DO

DO I=1,3
DO J=1,3
ST(I,J) = STI(I,J)+STFPE(I,J)+STFSE(I,J)+STJ(I,J)
END DO
END DO



STRESS(1) = ST(1,1)
STRESS(2) = ST(2,2)
STRESS(3) = ST(3,3)
STRESS(4) = ST(1,2)
STRESS(5) = ST(1,3)
STRESS(6) = ST(2,3)

C
C CALCULATE THE STIFFNESS
C
C31_UNIT4(I,J,K,L)--symmetric 4th order unit tensior

DO I=1,3
DO J=1,3
DO K=1,3
DO L=1,3
UNIT4(I,J,K,L)=0.5*(UNIT2(I,K)*UNIT2(J,L)+UNIT2(I,L)*UNIT2(J,K))
END DO
END DO
END DO
END DO


C32_
U2I=(b**2)*c*exp(b*(HIBARC1-3.0))


C33_
U2J=2.0/D
C34_
IF(ANMDF.GT.1.0) THEN
U2PE=TM0*8.0*(ANMDF-1.0)
ELSE
U2PE=0.0
END IF

C35_

U2SE=((1.+AK)/AK)*ALFA*BETA*EXP(ALFA*(RSS-1.)) ! Changed to the kojic
C36_
U2F=U2PE+U2SE

C_Stiffness

DO I=1,3
DO J=1,3
DO K=1,3
DO L=1,3

H4(I,J,K,L)=4.0*U2I*BBAR(I,J)* BBAR(K,L)/DET-
1 (4./3.)*(U1I+HIBARC1*U2I)*
2 (BBAR(I,J)*UNIT2(K,L)+UNIT2(I,J)*BBAR(K,L))/DET+
3 ((4./9.)*HIBARC1*(U1I+HIBARC1*U2I)+
4 (1./9.)*ANMDF*(U1F+ANMDF*U2F)+
5 DET*(U1J+DET*U2J))*
6 (UNIT2(I,J)*UNIT2(K,L))/DET+
7 (ANMDF**2)*(U2F-U1F/ANMDF)*
8 (AN1(I)*AN1(J)*AN1(K)*AN1(L))/DET-
9 (1./3.)*((ANMDF)**2)*(U2F+U1F/ANMDF)*
10 (AN1(I)*AN1(J)*UNIT2(K,L)+UNIT2(I,J)*AN1(K)*AN1(L))/DET+
11 (2./DET)*((2./3)*HIBARC1*U1I+(1./3.)*ANMDF*U1F-DET*U1J)*
12 UNIT4(I,J,K,L)

END DO
END DO
END DO
END DO

CALL T42(H4,H2)

C DDSDDE
DO I=1,6
DO J=1,3
DDSDDE(I,J)=H2(I,J)
END DO
END DO

DO I=1,6
DO J=4,6
DDSDDE(I,J)=0.5*(H2(I,J)+H2(I,J+3))
END DO
END DO

c
C to get the tension from Isometric tension-length curve of muscle in time t
c
CALL GETCURVE(SIGM0,RSP,PROPS2,NVALUE2)
c
c pass to main program
c
STATEV(1) = RSS ! rate of stretch in SEE
STATEV(2) = RSM ! rate of stretch in CE
STATEV(3) = SIGM0 ! maximum stress of CE at time t+dt
STATEV(4) = DRSM ! increment rate of stretch of CE
C Only as output as follow
STATEV(5) = SIGMS ! stress of SEE
STATEV(6) = SIGMM ! stress of CE
STATEV(7) = V ! the contractional Veolcity of CE
STATEV(8) = ALFAA
c Only as output as follow
statev(9) = STI(3,3)
statev(10) = STFPE(3,3)
statev(11) = STFSE(3,3)
statev(12) = STJ(3,3)

statev(13)=RSP
statev(14)=ANMDF

statev(15)=AN1(1)
statev(16)=AN1(2)
statev(17)=AN1(3)

RETURN
END


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!--------!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine getcurve(y,x,table,nvalue)
C IMPLICIT DOUBLE PRECISION (A-H,O-Z)
INCLUDE 'ABA_PARAM.INC'
dimension table(2,nvalue)
! set y to last of table, slope to zero
c write (7,*), "table1", table(1,
c write (7,*), "table2", table(2,
y=table(1,nvalue)
slope=0.0
C
do K1=1,nvalue-1
x2=table(2,K1+1)
if(x.LT.x2) then
x1=table(2,K1)
! computing current value y
dx=x2-x1
y1=table(1,K1)
y2=table(1,K1+1)
dy=y2-y1
slope=dy/dx
y=y1+(x-x1)*slope
goto 20
endif
enddo
20 continue
return
end

c the subroutine is used to solve equation numerically
c by standard Newton's method
C source program from FOR1,
C modified by zhangguang


C subroutine program:compute f(x)=F and df/dx=DY
SUBROUTINE FS(X,F,DY,A2,A3,A4,A5,ALFA)
C IMPLICIT DOUBLE PRECISION (A-H,O-Z)
INCLUDE 'ABA_PARAM.INC'
F = (A2+A3*X)*EXP(ALFA* X)-A4* X-A5
DY=(A3+ ALFA *(A2+A3* X))* EXP(ALFA* X)-A4
RETURN
END

C subroutine program: newton iteration method
SUBROUTINE DNEWT(X,EPS,L,A2,A3,A4,A5,ALFA)
INCLUDE 'ABA_PARAM.INC'
C IMPLICIT DOUBLE PRECISION (A-H,O-Z)
L=60 ! the most iterative number to be given
CALL FS(X,F,DY,A2,A3,A4,A5,ALFA)

30 IF (ABS(DY)+1.0.EQ.1.0) THEN
L=0
WRITE(*,2000)
RETURN
END IF

X1=X-F/DY
CALL FS(X1,F,DY,A2,A3,A4,A5,ALFA)
IF ((ABS(X1-X).GE.EPS).OR.(ABS(F).GE.EPS)) THEN
L=L-1
c WRITE(*,*) "L=",L,x1
X=X1
IF (L.EQ.0) RETURN
GOTO 30
END IF
X=X1
RETURN
2000 FORMAT(1X,' ERR')
END

C ËÄάÊý×éDR4(3x3x3x3)ת´æ¶þάÊý×éDR2(9x9)
C
SUBROUTINE T42(DR4,DR2)
INCLUDE 'ABA_PARAM.INC'
c IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C PASS IN
DIMENSION DR4(3,3,3,3)
C PASS OUT
DIMENSION DR2(9,9)
C
DO M=1,9
CM1
IF (M.EQ.1) THEN
I=1
J=1
END IF
CM2
IF (M.EQ.2) THEN
I=2
J=2
END IF
CM3
IF (M.EQ.3) THEN
I=3
J=3
END IF
CM4
IF (M.EQ.4) THEN
I=1
J=2
END IF
CM5
IF (M.EQ.5) THEN
I=1
J=3
END IF
CM6
IF (M.EQ.6) THEN
I=2
J=3
END IF
CM7
IF (M.EQ.7) THEN
I=2
J=1
END IF
CM8
IF (M.EQ.8) THEN
I=3
J=1
END IF
CM9
IF (M.EQ.9) THEN
I=3
J=2
END IF
CM END START N
DO N=1,9
CN1
IF (N.EQ.1) THEN
K=1
L=1
END IF
CN2
IF (N.EQ.2) THEN
K=2
L=2
END IF
CN3
IF (N.EQ.3) THEN
K=3
L=3
END IF
CN4
IF (N.EQ.4) THEN
K=1
L=2
END IF
CN5
IF (N.EQ.5) THEN
K=1
L=3
END IF
CN6
IF (N.EQ.6) THEN
K=2
L=3
END IF
CN7
IF (N.EQ.7) THEN
K=2
L=1
END IF
CN8
IF (N.EQ.8) THEN
K=3
L=1
END IF
CN9
IF (N.EQ.9) THEN
K=3
L=2
END IF
C
DR2(M,N)=DR4(I,J,K,L)
C
END DO
END DO

RETURN
END