Hi, this is my first time posting here.
I am trying to write a crystal plasticity UMAT that I will later modify as part of my research. I have been debugging it in Visual Studio and as far as I can tell each subroutine is calculating correctly (getting nonzero stress and Jacobian at the end of the routine, intermediate variables didn't go to zero or infinity). However when I try to run it in Abaqus with a single-element cube, it keeps throwing errors either that the minimum increment time is too big or that too many attempts were made per increment. Looking at the warnings, it keeps throwing "There is zero FORCE everywhere in the model based on the default criterion. please check the value of the average FORCE during the current iteration to verify that the FORCE is small enough to be treated as zero. if not, please use the solution controls to reset the criterion for zero FORCE." The .msg file also keeps saying that "THE FORCE IS ZERO EVERYWHERE BUT THE FORCE RESIDUAL OR THE DISP. CORRECTION IS NON-ZERO". Does anyone with more experience with Abaqus know what sort of issues cause these errors? I have provided my input file and UMAT subroutine below. Thanks in advance for any help!
Input file
*Heading
** Job name: Job-2 Model name: Job-1-UMAT_DD3
** Generated by: Abaqus/CAE 2023
*Preprint, echo=NO, model=NO, history=NO, contact=NO
**
** PARTS
**
*Part, name=Part-1
*Node
1, 5.99999985e-05, 5.99999985e-05, 5.99999985e-05
2, 5.99999985e-05, 1.62040791e-16, 5.99999985e-05
3, 5.99999985e-05, 5.99999985e-05, 0.
4, 5.99999985e-05, 1.62040791e-16, 0.
5, 1.62040791e-16, 5.99999985e-05, 5.99999985e-05
6, 1.62040791e-16, 1.62040791e-16, 5.99999985e-05
7, 1.62040791e-16, 5.99999985e-05, 0.
8, 1.62040791e-16, 1.62040791e-16, 0.
*Element, type=C3D8R
1, 5, 6, 8, 7, 1, 2, 4, 3
*Nset, nset=Set-1, generate
1, 8, 1
*Elset, elset=Set-1
1,
** Section: Section-1-SET-1
*Solid Section, elset=Set-1, controls=EC-1, material=MATERIAL-1
,
*End Part
**
**
** ASSEMBLY
**
*Assembly, name=Assembly
**
*Instance, name=Part-1-1, part=Part-1
*End Instance
**
*Nset, nset=Set-1, instance=Part-1-1, generate
1, 4, 1
*Elset, elset=Set-1, instance=Part-1-1
1,
*Nset, nset=Set-2, instance=Part-1-1
1, 2, 5, 6
*Elset, elset=Set-2, instance=Part-1-1
1,
*Nset, nset=Set-3, instance=Part-1-1, generate
2, 8, 2
*Elset, elset=Set-3, instance=Part-1-1
1,
*Nset, nset=Set-4, instance=Part-1-1
3, 4, 7, 8
*Elset, elset=Set-4, instance=Part-1-1
1,
*Elset, elset=_Surf-1_S5, internal, instance=Part-1-1
1,
*Surface, type=ELEMENT, name=Surf-1
_Surf-1_S5, S5
*End Assembly
**
** ELEMENT CONTROLS
**
*Section Controls, name=EC-1, hourglass=ENHANCED
1., 1., 1.
**
** MATERIALS
**
**
**
*Material, name=MATERIAL-1
*Depvar
246,
*User Material, constants=52, unsymm
1.08e+11, 6.13e+10, 2.85e+10, 2.5e+10, 0., 0., 0., 0.
0., 0., 0., 0., 0., 0., 0., 0.
0., 0., 0., 0., 0., 0., 0., 0.
1., 1., 1., 1., 1., 1., 0., 0.
4.16e+10, 4.16e+10, 8.33e+10, 8.33e+10, 0., 0., 0., 0.
1.1, 1.1, 0.141, 0.141, 1.86e-08, 9.3e-08, 3e-19, 3e-19
2e+06, 2e+06, 298., 2.863e-10
**
** BOUNDARY CONDITIONS
**
** Name: BC-1 Type: Displacement/Rotation
*Boundary
Set-1, 1, 1
** Name: BC-2 Type: Displacement/Rotation
*Boundary
Set-2, 3, 3
** Name: BC-3 Type: Displacement/Rotation
*Boundary
Set-3, 2, 2
** ----------------------------------------------------------------
**
** STEP: Step-1
**
*Step, name=Step-1, nlgeom=NO, inc=1000
*Static
1e-06, 0.1, 1e-20, 0.01
**
** LOADS
**
** Name: Load-1 Type: Pressure
*Dsload
Surf-1, P, -5000.
**
** CONTROLS
**
*Controls, reset
*Controls, parameters=time incrementation
, , , , , , , 15, , ,
**
** OUTPUT REQUESTS
**
*Restart, write, frequency=0
**
** FIELD OUTPUT: F-Output-2
**
*Output, field
*Element Output, directions=YES
LE, PE, PEEQ, PEMAG, S, SDV
**
** FIELD OUTPUT: F-Output-1
**
*Node Output
U,
**
** HISTORY OUTPUT: H-Output-1
**
*Output, history, variable=PRESELECT
*End Step
UMAT
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 DISLOCATION-DENSITY-BASED CRYSTAL PLASTICITY MODEL
C THIS MODEL IS FOR FCC SINGLE CRYSTAL
C DEFINE EXTERNAL VARIABLES
C STRESS(NTENS) = THIS ARRAY IS PASSED IN AS THE STRESS TENSOR AT
C THE BEGINNING OF THE INCREMENT AND MUST BE UPDATED IN THIS
C ROUTINE TO BE THE STRESS TENSOR AT THE END OF THE INCREMENT
C NTENS = SIZE OF THE STRESS OR STRAIN COMPONENT ARRAY (NDI + NSHR)
C NDI = NUMBER OF DIRECT STRESS COMPONENTS AT THIS POINT
C NSHR = NUMBER OF ENGINEERING SHEAR STRESS COMPONENTS AT
C THIS POINT
C STATEV(NSTATV) = ARRAY CONTAINING THE SOLUTION-DEPENDENT
C STATE VARIABLES
C DDSDDE(NTENS,NTENS) = JACOBIAN MATRIX OF THE CONSTITUTIVE MODEL.
C UNLESS YOU INVOKE THE UNSYMMETRIC EQUATION SOLUTION CAPABILITY
C FOR THE USER-DEFINED MATERIAL, ABAQUS/STANDARD WILL USE ONLY THE
C SYMMETRIC PART OF DDSDDE
C DDSDDT(NTENS) = VARIATION OF THE STRESS INCREMENTS WITH RESPECT
C TO THE TEMPERATURE
C STRAN(NTENS) = ARRAY CONTAINING THE TOTAL STRAINS AT THE
C BEGINNING OF THE INCREMENT (MECHANICAL STRAINS ONLY)
C DSTRAN(NTENS) = ARRAY OF STRAIN INCREMENTS
C (MECHANICAL STRAINS ONLY)
C TIME(1) = VALUE OF STEP TIME AT THE BEGINNING OF THE INCREMENT
C TIME(2) = VALUE OF TOTAL TIME AT THE BEGINNING OF THE INCREMENT
C DTIME = TIME INCREMENT
C TEMP = TEMPERATURE AT THE BEGINNING OF THE INCREMENT
C DTEMP = INCREMENT OF TEMPERATURE
C PREDEF = ARRAY OF INTERPOLATED VALUES OF PREDEFINED FIELD
C VARIABLES AT THIS POINT AT THE START OF THE INCREMENT, BASED ON
C VALUES READ IN AT THE NODES
C DPRED = ARRAY OF INCREMENTS OF PREDEFINED FIELD VARIABLES
C CMNAME = USER-DEFINED MATERIAL NAME, LEFT JUSTIFIED
C (DO NOT USE "ABQ_" AS THE LEADING STRING)
C NSTATV = NUMBER OF SOLUTION-DEPENDENT STATE VARIABLES
C PROPS(NPROPS) = USER-SPECIFIED ARRAY OF MATERIAL CONSTANTS
C NPROPS = USER-DEFINED NUMBER OF MATERIAL CONSTANTS
C COORDS = ARRAY CONTAINING COORDINATES OF THE POINT, CURRENT
C COORDINATES IF GEOMETRIC NONLINEARITY IS ACCOUNTED FOR
C DROT(3,3) = ROTATION INCREMENT MATRIX
C CELENT = CHARACTERISTIC ELEMENT LENGTH
C DFGRD0(3,3) = ARRAY CONTAINING DEFORMATION GRADIENT AT
C BEGINNING OF INCREMENT
C DFGRD1(3,3) = ARRAY CONTAINING DEFORMATION GRADIENT AT END
C OF INCREMENT
C NOEL = ELEMENT NUMBER
C NPT = INTEGRATION POINT NUMBER
C LAYER = LAYER NUMBER (FOR COMPOSITE SHELLS/LAYERED SOLIDS)
C KSPT = SECTION POINT NUMBER WITHIN CURRENT LAYER
C JSTEP(1) = STEP NUMBER
C JSTEP(2) = PROCEDURE TYPE KEY
C JSTEP(3) = 1 IF NLGEOM = YES FOR CURRENT STEP, OTHERWISE 0
C JSTEP(4) = 1 IF CURRENT STEP IS LINEAR PERTURBATION
C PROCEDURE, OTHERWISE 0
C KINC = INCREMENT NUMBER
C DEFINE EXTERNAL VARIABLES
IMPLICIT NONE
INTEGER::NTENS,NSTATV,NPROPS,NDI
INTEGER::NSLIP,I,J
REAL*8::STRESS(NTENS),STATEV(NSTATV),DDSDDE(NTENS,NTENS),
1 DDSDDT(NTENS),DRPLDE(NTENS),STRAN(NTENS),DSTRAN(NTENS),
2 TIME(2),PREDEF(1),DPRED(1),PROPS(NPROPS),COORDS(3),
3 DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),DTIME
C DEFINE UNUSED VARIABLES PASSED IN BY UMAT
INTEGER::KSTEP,KINC,NOEL,NPT,LAYER,KSPT,NSHR
REAL*8::SSE,SPD,SCD,RPL,DRPLDT,TEMP,DTEMP,PNEWDT,CELENT
CHARACTER*8::CMNAME
C DEFINE INTERNAL VARIABLES
REAL*8::SLIPSYSNORM(36),SLIPSYSDIR(36),SLIPSYSTAN(36),
1 F_OLD(NDI,NDI),F(NDI,NDI),FP_OLD(NDI,NDI),
2 C(6,6),GAMMA_DOT(12),
3 FP_INC(NDI,NDI),FP(NDI,NDI),FPINV(NDI,NDI),FE(NDI,NDI),
4 S(NDI,NDI),SIGMA(NDI,NDI),TAU(12)
C NUMBER OF SLIP SYSTEMS
C DEFINE PROPS
C PROPS(1-24) == ELASTIC PROPERTIES
C PROPS(1) = C11 (Pa)
C PROPS(2) = C12
C PROPS(3) = C44
C PROPS(4) = SHEAR MODULUS MU (Pa)
C PROPS(25-32) == INITIAL DISLOCATION MOBILITY PROPERTIES
C PROPS(25) = INITIAL EDGE DISLOCATION MOBILITY (m/s)
C PROPS(26) = INITIAL SCREW DISLOCATION MOBILITY (m/s)
C PROPS(27) = INITIAL EDGE DIS. MOBILITY, POSITIVE DIR (m/s)
C PROPS(28) = INITIAL EDGE DIS. MOBILITY, NEGATIVE DIR (m/s)
C PROPS(29) = INITIAL SCREW DIS. MOBILITY, RIGHT-HAND (m/s)
C PROPS(30) = INITIAL SCREW DIS. MOBILITY, LEFT-HAND (m/s)
C PROPS(33-40) == INITIAL DISLOCATION DENSITY PROPERTIES
C PROPS(33) = INITIAL EDGE DIS. DENSITY, POSITIVE DIR (m**-2)
C PROPS(34) = INITIAL EDGE DIS. DENSITY, NEGATIVE DIR (m**-2)
C PROPS(35) = INITIAL SCREW DIS. DENSITY, RIGHT-HAND (m**-2)
C PROPS(36) = INITIAL SCREW DIS. DENSITY, LEFT-HAND (m**-2)
C PROPS(41-48) == ADDITIONAL DISLOCATION DENSITY CONSTANTS
C PROPS(41) = EDGE DISLOCATION EXPONENT PARAMETER P_E
C PROPS(42) = SCREW DISLOCATION EXPONENT PARAMETER P_S
C PROPS(43) = EDGE DISLOCATION EXPONENT PARAMETER Q_E
C PROPS(44) = SCREW DISLOCATION EXPONENT PARAMETER Q_S
C PROPS(45) = EDGE DISLOCATION CAPTURE RADIUS R_E (m)
C PROPS(46) = SCREW DISLOCATION CAPTURE RADIUS R_S (m)
C PROPS(47) = EDGE DISLOCATION ACTIVATION ENERGY F_E (J/ATOM)
C PROPS(48) = SCREW DISLOCATION ACTIVATION ENERGY F_S (J/ATOM)
C PROPS(49-56) == ADDITIONAL DISLOCATION DENSITY CONSTANTS
C PROPS(49) = INTRINSIC LATTICE RESISTANCE S_P (Pa)
C PROPS(50) = RESISTANCE DUE TO FOREST DISLOCATION
C INTERACTIONS S_D (Pa)
C PROPS(51) = ABSOLUTE TEMPERATURE TEMP (K)
C PROPS(52) = MAGNITUDE OF A BURGERS VECTOR BMAG (m)
C DEFINE STATEV
C STATEV(1-12) == RESOLVED SHEAR STRESS PER SLIP SYSTEM
C STATEV(13-24) == SHEAR STRAIN PER SLIP SYSTEM
C STATEV(25-60) == COMPONENTS OF SLIP PLANE NORMAL VECTOR PER
C SLIP SYSTEM (123,456,...)
C STATEV(61-96) == COMPONENTS OF SLIP DIRECTION VECTOR PER SLIP
C SYSTEM (123,456,...)
C STATEV(97-132) == COMPONENTS OF DISLOCATION TANGENT VECTOR PER
C SLIP SYSTEM
C STATEV(133-144) == EDGE DIS. DENSITY, POSITIVE DIR PER SLIP SYS
C STATEV(145-156) == EDGE DIS. DENSITY, NEGATIVE DIR PER SLIP SYS
C STATEV(157-168) == SCREW DIS. DENSITY, RIGHT-HAND PER SLIP SYS
C STATEV(169-180) == SCREW DIS. DENSITY, LEFT-HAND PER SLIP SYS
C STATEV(181-192) == EDGE DIS. MOBILITY, POSITIVE DIR PER SLIP SYS
C STATEV(193-204) == EDGE DIS. MOBILITY, NEGATIVE DIR PER SLIP SYS
C STATEV(205-216) == SCREW DIS. MOBILITY, RIGHT-HAND PER SLIP SYS
C STATEV(217-228) == SCREW DIS. MOBILITY, LEFT-HAND PER SLIP SYS
C STATEV(229-237) == PLASTIC DEFORMATION GRADIENT
C STATEV(238-246) == ELASTIC DEFORMATION GRADIENT
WRITE(6,*) 'STRESS'
WRITE(6,*) STRESS
WRITE(6,*) 'DDSDDE'
WRITE(6,*) DDSDDE
WRITE(6,*) 'TIME'
WRITE(6,*) TIME
WRITE(6,*) 'DTIME'
WRITE(6,*) DTIME
WRITE(6,*) 'DFGRD0'
WRITE(6,*) DFGRD0
WRITE(6,*) 'DFGRD1'
WRITE(6,*) DFGRD1
C INITIAL TIME STEP PROCEDURE
NSLIP = 12
C INITIAL SLIP SYSTEM VECTORS
SLIPSYSNORM=(/1,1,1,
1 1,1,1,
2 1,1,1,
3 -1,1,1,
4 -1,1,1,
5 -1,1,1,
6 1,-1,1,
7 1,-1,1,
8 1,-1,1,
9 1,1,-1,
1 1,1,-1,
2 1,1,-1/)
SLIPSYSDIR=(/0,-1,1,
1 1,0,-1,
2 -1,1,0,
3 1,0,1,
4 1,1,0,
5 0,-1,1,
6 0,1,1,
7 1,1,0,
8 1,0,-1,
9 0,1,1,
1 1,0,1,
2 -1,1,0/)
SLIPSYSTAN=(/-2,1,1,
1 1,-2,1,
2 1,1,-2,
3 1,2,-1,
4 1,-1,2,
5 2,1,1,
6 2,1,-1,
7 -1,1,2,
8 1,2,1,
9 2,-1,1,
1 -1,2,1,
2 1,1,2/)
C INITIALIZE STATEV ARRAY
IF (TIME(2).EQ.0.) THEN
DO I=1,NSTATV
C INITIAL RESOLVED SHEAR STRESS
IF (I.LE.12) THEN
STATEV(I) = 0.
C INITIAL SHEAR STRAIN
ELSE IF (I.LE.24) THEN
STATEV(I) = 0.
C INITIAL SLIP SYSTEM NORMAL VECTORS
ELSE IF (I.LE.60) THEN
STATEV(I) = SLIPSYSNORM(I-24)
C INITIAL SLIP SYSTEM DIRECTION VECTORS
ELSE IF (I.LE.96) THEN
STATEV(I) = SLIPSYSDIR(I-60)
C INITIAL SLIP SYSTEM TANGENT VECTORS
ELSE IF (I.LE.132) THEN
STATEV(I) = SLIPSYSTAN(I-96)
C INITIAL EDGE DISLOCATION DENSITY PER SLIP SYSTEM
ELSE IF (I.LE.156) THEN
STATEV(I) = PROPS(33)
C INITIAL SCREW DISLOCATION DENSITY PER SLIP SYSTEM
ELSE IF (I.LE.180) THEN
STATEV(I) = PROPS(35)
C INITIAL EDGE DISLOCATION MOBILITY PER SLIP SYSTEM
ELSE IF (I.LE.204) THEN
STATEV(I) = PROPS(27)
C INITIAL SCREW DISLOCATION MOBILITY PER SLIP SYSTEM
ELSE IF (I.LE.228) THEN
STATEV(I) = PROPS(28)
END IF
END DO
C INITIAL DEFORMATION GRADIENT
STATEV(229) = 1.
STATEV(230) = 0.
STATEV(231) = 0.
STATEV(232) = 0.
STATEV(233) = 1.
STATEV(234) = 0.
STATEV(235) = 0.
STATEV(236) = 0.
STATEV(237) = 1.
STATEV(238) = STATEV(229)
STATEV(239) = STATEV(230)
STATEV(240) = STATEV(231)
STATEV(241) = STATEV(232)
STATEV(242) = STATEV(233)
STATEV(243) = STATEV(234)
STATEV(244) = STATEV(235)
STATEV(245) = STATEV(236)
STATEV(246) = STATEV(237)
END IF
C GENERAL PROCEDURE
C READ PREVIOUS TIME STEP AND CURRENT TIME STEP
C DEFORMATION GRADIENTS
F_OLD = DFGRD0
F = DFGRD1
C READ PREVIOUS TIME STEP PLASTIC DEFORMATION
C GRADIENT TO BE UPDATED
FP_OLD = 0.
DO I=1,NDI
DO J=1,NDI
FP_OLD(I,J) = STATEV(228+((I-1)*NDI)+J)
END DO
END DO
C DEFINE STIFFNESS MATRIX C
CALL STIFFNESSMATRIX (PROPS,NPROPS,C)
C CALCULATE RESOLVED SHEAR STRESS
CALL RSS (STRESS,SLIPSYSDIR,SLIPSYSNORM,NSLIP,TAU)
C CALCULATE SLIP RATE
CALL SLIPRATE (STATEV,NSTATV,PROPS(52),NSLIP,TAU,GAMMA_DOT)
C CALCULATE PLASTIC DEFORMATION GRADIENT INCREMENT
CALL FPSTEP (FP_OLD,SLIPSYSDIR,SLIPSYSNORM,DTIME,NSLIP,NDI,
1 GAMMA_DOT,FP_INC)
C UPDATE PLASTIC DEFORMATION GRADIENT
FP = FP_OLD + FP_INC
C CALCULATE NEW ELASTIC DEFORMATION GRADIENT
C FIND INVERSE OF FP
CALL INVERSEMATRIX (FP, NDI, FPINV)
FE = F*FPINV
C FIND RIGHT CAUCHY-GREEN DEFORMATION TENSOR (DO I NEED THIS HERE?)
C CALL RCGDEFTENSOR (FE, NDI, C)
C FIND CAUCHY STRESS
CALL CAUCHYSTRESS (FE,NDI,C,SIGMA)
C CALCULATE SLIP RATE GAMMA_DOT USING DISLOCATION
C DENSITY FORMULATION
CALL DISLOCATIONDENSITY (STATEV,NSTATV,PROPS,NPROPS,
1 DTIME,GAMMA_DOT)
C WRITE TO STRESS WITH CALCULATED CAUCHY STRESS
STRESS(1) = SIGMA(1,1)
STRESS(2) = SIGMA(2,2)
STRESS(3) = SIGMA(3,3)
STRESS(4) = SIGMA(1,2)
STRESS(5) = SIGMA(1,3)
STRESS(6) = SIGMA(2,3)
C KEEP JACOBIAN CONSTANT FOR RIGHT NOW
CALL STIFFNESSMATRIX (PROPS,NPROPS,DDSDDE)
C WRITE TO STATEV
C NEW PLASTIC DEFORMATION GRADIENT
DO I=1,NDI
DO J=1,NDI
STATEV(228+((I-1)*NDI)+J) = FP(I,J)
END DO
END DO
WRITE(6,*) 'NEW STRESS'
WRITE(6,*) STRESS
WRITE(6,*) 'NEW DDSDDE'
WRITE(6,*) DDSDDE
C READ(*,*)
C END UMAT
RETURN
END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C RESOLVED SHEAR STRESS SUBROUTINE
SUBROUTINE RSS (STRESS,M_LIST,N_LIST,NSLIP,TAU)
IMPLICIT NONE
INTEGER::NSLIP,ALPHA,I,J
REAL*8::STRESS(6),M_LIST(NSLIP*3),N_LIST(NSLIP*3),TAU(NSLIP),
1 SIGMA(3,3),M(3),N(3),PTENS(3,3)
TAU = 0.
C CONVERT STRESS ARRAY TO CAUCHY STRESS TENSOR
SIGMA(1,1) = STRESS(1)
SIGMA(2,2) = STRESS(2)
SIGMA(3,3) = STRESS(3)
SIGMA(1,2) = STRESS(4)
SIGMA(2,1) = STRESS(4)
SIGMA(1,3) = STRESS(5)
SIGMA(3,1) = STRESS(5)
SIGMA(2,3) = STRESS(6)
SIGMA(3,2) = STRESS(6)
C RESOLVED SHEAR STRESS PER SLIP SYSTEM
DO ALPHA=1,NSLIP
C SLIP PLANE DIRECTION AND NORMAL VECTORS FOR THE SLIP SYSTEM
M = M_LIST(((ALPHA-1)*3)+1 : ALPHA*3)
N = N_LIST(((ALPHA-1)*3)+1 : ALPHA*3)
C CALCULATE SCHMID TENSOR FOR THE SLIP SYSTEM
CALL SCHMID_TENSOR(M,N,PTENS)
C PERFORM TENSOR DOT PRODUCT OF CAUCHY STRESS WITH SCHMID TENSOR
C TO GET RESOLVED SHEAR STRESS FOR THE SLIP SYSTEM
DO I=1,3
DO J=1,3
TAU(ALPHA) = TAU(ALPHA) + SIGMA(I,J)*PTENS(I,J)
END DO
END DO
END DO
RETURN
END
C SLIP RATE SUBROUTINE
SUBROUTINE SLIPRATE (STATEV,NSTATV,BMAG,NSLIP,TAU,GAMMA_DOT)
INTEGER::NSTATV,ALPHA,NSLIP
REAL*8::STATEV(NSTATV),RHO_EP(NSLIP),RHO_EN(NSLIP),RHO_SP(NSLIP),
1 RHO_SN(NSLIP),V_EP(NSLIP),V_EN(NSLIP),V_SP(NSLIP),
2 V_SN(NSLIP),BMAG,TAU(NSLIP),GAMMA_DOT(NSLIP)
C READ DENSITY AND MOBILITY VALUES FROM STATEV
RHO_EP = STATEV(133:144)
RHO_EN = STATEV(145:156)
RHO_SP = STATEV(157:168)
RHO_SN = STATEV(169:180)
V_EP = STATEV(181:192)
V_EN = STATEV(193:204)
V_SP = STATEV(205:216)
V_SN = STATEV(217:228)
C CALCULATE SLIP RATE GAMMA_DOT PER SLIP SYSTEM
DO ALPHA=1,NSLIP
GAMMA_DOT(ALPHA) = (RHO_EP(ALPHA)*V_EP(ALPHA)
1 + RHO_EN(ALPHA)*V_EN(ALPHA) + RHO_SP(ALPHA)*V_SP(ALPHA)
2 + RHO_SN(ALPHA)*V_SN(ALPHA))*BMAG*SIGN(1.,TAU(ALPHA))
END DO
RETURN
END SUBROUTINE SLIPRATE
C PLASTIC DEFORMATION GRADIENT INCREMENT SUBROUTINE
SUBROUTINE FPSTEP (FP_OLD,M_LIST,N_LIST,DTIME,NSLIP,NDI,
1 GAMMA_DOT,FP_INC)
C DEFINE VARIABLES
INTEGER::I,NSLIP
REAL*8::FP_OLD(NDI,NDI),M_LIST(NSLIP*3),N_LIST(NSLIP*3),
1 LP(NDI,NDI),FP_INC(NDI,NDI), M(3), N(3), PTENS(3,3),
2 DTIME,GAMMA_DOT(NSLIP)
C CALCULATE PLASTIC VELOCITY GRADIENT
LP = 0.
DO I=1,NSLIP
C SLIP PLANE DIRECTION AND NORMAL VECTORS FOR THE SLIP SYSTEM
M = M_LIST(((I-1)*3)+1 : I*3)
N = N_LIST(((I-1)*3)+1 : I*3)
C CALCULATE SCHMID TENSOR FOR THE SLIP SYSTEM
CALL SCHMID_TENSOR(M,N,PTENS)
LP = LP + GAMMA_DOT(I)*PTENS
END DO
C CALCULATE INCREMENT
FP_INC = DTIME*LP*FP_OLD
RETURN
END SUBROUTINE FPSTEP
C SLIP SYSTEM SCHMID TENSOR SUBROUTINE
SUBROUTINE SCHMID_TENSOR(M,N,PTENS)
C DEFINE VARIABLES
REAL*8::M(3),N(3),PTENS(3,3)
C TENSOR OUTER PRODUCT OF 3X1 VECTOR WITH 1X3 VECTOR IS 3X3 MATRIX
DO I=1,3
DO J=1,3
PTENS(I,J) = M(I)*N(J)
END DO
END DO
RETURN
END SUBROUTINE SCHMID_TENSOR
C RIGHT CAUCHY-GREEN DEFORMATION TENSOR SUBROUTINE
SUBROUTINE RCGDEFTENSOR (FE, NDI, C)
INTEGER::NDI
REAL*8::FE(NDI,NDI),FE_TRANS(NDI,NDI),C(NDI,NDI)
C FIND TRANSPOSE OF DEFORMATION GRADIENT F
FE_TRANS = TRANSPOSE(FE)
C FIND CAUCHY-GREEN DEFORMATION TENSOR
C = FE_TRANS*FE
RETURN
END SUBROUTINE RCGDEFTENSOR
C CAUCHY STRESS TENSOR SUBROUTINE
SUBROUTINE CAUCHYSTRESS (FE,NDI,C,SIGMA)
INTEGER::NDI
REAL*8::FE(NDI,NDI),C(2*NDI,2*NDI),FE_TRANS(NDI,NDI),S(NDI,NDI),
1 SIGMA(NDI,NDI),FE_DET
C FIND 2ND PIOLA-KIRCHOFF STRESS TENSOR
CALL SECOND_PIOLA_KIRCHOFF (FE,C,NDI,S)
C FIND TRANSPOSE OF DEFORMATION GRADIENT FE
FE_TRANS = TRANSPOSE(FE)
C FIND DETERMINANT OF FE
CALL DETERMINANTMATRIX(FE,NDI,FE_DET)
C CALCULATE CAUCHY STRESS
SIGMA = (1/FE_DET)*FE*S*FE_TRANS
RETURN
END SUBROUTINE CAUCHYSTRESS
C 2ND PIOLA-KIRCHOFF STRESS TENSOR SUBROUTINE
SUBROUTINE SECOND_PIOLA_KIRCHOFF(FE,C,NDI,S)
INTEGER::NDI,I,J,K,L
REAL*8::FE(NDI,NDI),C(2*NDI,2*NDI),EE(NDI,NDI),
1 BIGC(NDI,NDI,NDI,NDI),S(NDI,NDI)
C FIND ELASTIC GREEN-LAGRANGE STRAIN TENSOR
CALL EGLEPSILON (FE,NDI,EE)
C FIND 4TH ORDER ELASTICITY TENSOR
CALL BIGELASTICITY (C,NDI,BIGC)
C CALCULATE 2ND PIOLA-KIRCHOFF STRESS TENSOR VIA CONTRACTION
S = 0.
DO I=1,NDI
DO J=1,NDI
DO K=1,NDI
DO L=1,NDI
S(I,J) = BIGC(I,J,K,L)*EE(K,L)
END DO
END DO
END DO
END DO
RETURN
END SUBROUTINE SECOND_PIOLA_KIRCHOFF
C ELASTIC GREEN-LAGRANGE STRAIN TENSOR SUBROUTINE
SUBROUTINE EGLEPSILON (FE, NDI, EE)
INTEGER::NDI,I
REAL*8::FE(NDI,NDI),FE_TRANS(NDI,NDI),IDENTITY(NDI,NDI),
1 EE(NDI,NDI)
C CREATE IDENTITY MATRIX
IDENTITY = 0.
DO I=1,NDI
IDENTITY(I,I) = 1.
END DO
C FIND TRANSPOSE OF ELASTIC DEFORMATION GRADIENT FE
FE_TRANS = TRANSPOSE(FE)
C CALCULATE ELASTIC GREEN-LAGRANGE STRAIN TENSOR
EE = 0.5*(FE_TRANS*FE - IDENTITY)
RETURN
END SUBROUTINE EGLEPSILON
C 4TH ORDER ELASTICITY TENSOR SUBROUTINE
SUBROUTINE BIGELASTICITY (C,NDI,BIGC)
INTEGER::NDI,I,J,K,L
REAL*8::C(2*NDI,2*NDI),BIGC(NDI,NDI,NDI,NDI)
C INITIALIZE 4TH ORDER TENSOR
BIGC = 0.
C MAP 2ND ORDER TENSOR C TO 4TH ORDER TENSOR BIGC
BIGC(1,1,1,1) = C(1,1)
BIGC(2,2,2,2) = C(2,2)
BIGC(3,3,3,3) = C(3,3)
BIGC(1,1,2,2) = C(1,2)
BIGC(1,1,3,3) = C(1,3)
BIGC(2,2,3,3) = C(2,3)
BIGC(2,2,1,1) = C(1,2)
BIGC(3,3,1,1) = C(1,3)
BIGC(3,3,2,2) = C(2,3)
BIGC(2,3,2,3) = C(4,4)
BIGC(3,2,3,2) = C(4,4)
BIGC(3,2,2,3) = C(4,4)
BIGC(2,3,3,2) = C(4,4)
BIGC(1,3,1,3) = C(5,5)
BIGC(3,1,3,1) = C(5,5)
BIGC(3,1,1,3) = C(5,5)
BIGC(1,3,3,1) = C(5,5)
BIGC(1,2,1,2) = C(6,6)
BIGC(2,1,2,1) = C(6,6)
BIGC(2,1,1,2) = C(6,6)
BIGC(1,2,2,1) = C(6,6)
C USE SYMMETRY TO POPULATE REST OF 4TH ORDER TENSOR
DO I=1,NDI
DO J=1,NDI
DO K=1,NDI
DO L=1,NDI
BIGC(J,I,K,L) = BIGC(I,J,K,L)
BIGC(I,J,L,K) = BIGC(I,J,K,L)
BIGC(J,I,L,K) = BIGC(I,J,K,L)
BIGC(K,L,I,J) = BIGC(I,J,K,L)
BIGC(L,K,I,J) = BIGC(I,J,K,L)
BIGC(K,L,J,I) = BIGC(I,J,K,L)
BIGC(L,K,J,I) = BIGC(I,J,K,L)
END DO
END DO
END DO
END DO
RETURN
END SUBROUTINE BIGELASTICITY
C DETERMINANT OF A MATRIX SUBROUTINE
SUBROUTINE DETERMINANTMATRIX (A, DIM1, A_DET)
INTEGER::DIM1,INFO,P(DIM1),I
REAL*8::A(DIM1,DIM1),L(DIM1,DIM1),U(DIM1,DIM1),A_DET
C CHECK IF A IS A SQUARE MATRIX
IF (SIZE(A,1).NE.SIZE(A,2)) THEN
PRINT*,"ERROR IN INVERSEMATRIX: MATRIX A MUST BE SQUARE"
RETURN
END IF
C INITIALIZE A_DET
A_DET = 1.
C IF A IS A 1X1 MATRIX
IF (DIM1.EQ.1) THEN
A_DET = A(DIM1,DIM1)
C RETURN
C IF A IS A 2X2 MATRIX:
ELSE IF (DIM1.EQ.2) THEN
A_DET = A(1,1)*A(2,2) - A(1,2)*A(2,1)
C RETURN
C OTHERWISE:
ELSE
C FIND THE UPPER TRIANGULAR MATRIX (REDUCING A TO ECHELON FORM)
C THE PRODUCT OF THE ENTRIES ON THE MAIN DIAGONAL IS A_DET
CALL LUDCMP (A, DIM1, L, U, P, INFO)
DO I=1,DIM1
A_DET = A_DET*U(I,I)
END DO
C RETURN
END IF
RETURN
END SUBROUTINE DETERMINANTMATRIX
C INVERSE MATRIX SUBROUTINE
SUBROUTINE INVERSEMATRIX (A, DIM1, A_INV)
INTEGER:: DIM1,INFO,P(DIM1),J,I
REAL*8::A(DIM1,DIM1),A_INV(DIM1,DIM1),L(DIM1,DIM1),U(DIM1,DIM1),
1 X(DIM1),B(DIM1),Y(DIM1)
C CHECK IF A IS A SQUARE MATRIX
IF (SIZE(A,1).NE.SIZE(A,2)) THEN
PRINT*,"ERROR IN INVERSEMATRIX: MATRIX A MUST BE SQUARE"
RETURN
END IF
C DETERMINE SIZE OF MATRIX A
DIM1 = SIZE(A,1)
C PERFORM LU DECOMPOSITION
CALL LUDCMP (A, DIM1, L, U, P, INFO)
C IF SINGULAR MATRIX FLAG IS RAISED
IF (INFO.NE.0) THEN
RETURN
END IF
C CALCULATE INVERSE OF MATRIX A
A_INV = 0.
C INVERT MATRIX COLUMN BY COLUMN
DO J=1,DIM1
C INITIALIZE COLUMN J OF THE IDENTITY MATRIX
B = 0.
B(J) = 1.
C SOLVE LY=PB FOR Y USING FORWARD SUBSTITUTION
CALL FORWARD_SUB (L, P, B, Y, DIM1)
C SOLVE UX=Y FOR X USING BACKWARD SUBSTITUTION
CALL BACKWARD_SUB (U, Y, X, DIM1)
C PLACE SOLUTION IN COLUMN J OF A_INV
DO I=1,DIM1
A_INV(I,J) = X(I)
END DO
END DO
RETURN
END SUBROUTINE INVERSEMATRIX
C FORWARD SUBSTITUTION SUBROUTINE
SUBROUTINE FORWARD_SUB (L, P, B, Y, DIM1)
INTEGER::DIM1,P(DIM1),I,J
REAL*8::L(DIM1,DIM1),B(DIM1),Y(DIM1),SUM
C SOLVE LY=PB
DO I=1,DIM1
SUM=0.
DO J=1,I-1
SUM = SUM + L(I,J)*Y(J)
END DO
Y(I) = B(P(I)) - SUM
END DO
RETURN
END SUBROUTINE FORWARD_SUB
C BACKWARD SUBSTITUTION SUBROUTINE
SUBROUTINE BACKWARD_SUB (U, Y, X, DIM1)
INTEGER::DIM1,I,J
REAL*8::U(DIM1,DIM1),Y(DIM1),X(DIM1),SUM
C SOLVE UX=Y
C GO IN REVERSE ORDER
DO I=DIM1,1,-1
SUM = 0.
DO J=I+1,DIM1
SUM = SUM + U(I,J)*X(J)
END DO
X(I) = (Y(I) - SUM) / U(I,I)
END DO
RETURN
END SUBROUTINE BACKWARD_SUB
C LOWER-UPPER DECOMPOSITION SUBROUTINE
SUBROUTINE LUDCMP (A, DIM1, L, U, P, INFO)
INTEGER::DIM1,INFO,K,I,J,MAX_I,P(DIM1)
REAL*8::A(DIM1,DIM1),A_MOD(DIM1,DIM1),L(DIM1,DIM1),U(DIM1,DIM1),
1 MAX_A
C INITIALIZE L AND U MATRICES, P ARRAY
L = 0.
U = 0.
C A_MOD IS COPY OF A THAT WILL BE REARRANGED AS NECESSARY
A_MOD = A
DO I=1,DIM1
L(I,I) = 1.
P(I) = I
END DO
C PERFORM LU DECOMPOSITION WITH PARTIAL PIVOTING TO REDUCE ERROR
INFO = 0
C COLUMN LOOP
DO K=1,DIM1
C PIVOTING
MAX_A = 0.
MAX_I = K
C ROW LOOP
DO I=K,DIM1
C CHECK IF VALUE WHICH ROW IS MAX VALUE IN COLUMN
IF (ABS(A_MOD(I,K)) > MAX_A) THEN
MAX_A = ABS(A_MOD(I,K))
MAX_I = I
END IF
END DO
C SINGULAR MATRIX FLAG
IF (MAX_A == 0.) THEN
INFO = 1
RETURN
END IF
C SWAP ROWS K AND MAX_I IN A AND P IF PIVOT VALUE IS NOT MAX
IF (MAX_I.NE.K) THEN
CALL PIVOT_SWAP(A_MOD, P, K, MAX_I, DIM1)
END IF
C POPULATE L AND U MATRICES
DO I=K,DIM1
U(K,I) = A_MOD(K,I)
END DO
DO I=K+1,DIM1
L(I,K) = A_MOD(I,K) / U(K,K)
DO J=K+1,DIM1
A_MOD(I,J) = A_MOD(I,J) - L(I,K)*U(K,J)
END DO
END DO
END DO
DO I=1,DIM1
DO J=I,DIM1
U(I,J) = A_MOD(I,J)
END DO
END DO
RETURN
END SUBROUTINE LUDCMP
C PIVOT ROW SWAP SUBROUTINE
SUBROUTINE PIVOT_SWAP (A, P, K, MAX_I, DIM1)
INTEGER::DIM1,K,J,TEMP_P,P(DIM1),MAX_I
REAL*8::A(DIM1,DIM1),TEMP
C SWAP ROWS K AND MAX_I IN MATRIX A, ARRAY LABEL IN ARRAY P
DO J=1,DIM1
C TEMPORARY STORAGE OF VALUE FROM ROW K
TEMP = A(K,J)
TEMP_P = P(K)
C REPLACE VALUE FROM ROW K WITH VALUE FROM ROW MAX_I
A(K,J) = A(MAX_I,J)
P(K) = p(MAX_I)
C REPLACE VALUE FROM ROW MAX_I WITH VALUE FROM ROW K
A(MAX_I,J) = TEMP
P(MAX_I) = TEMP_P
END DO
RETURN
END SUBROUTINE PIVOT_SWAP
C STIFFNESS MATRIX SUBROUTINE
SUBROUTINE STIFFNESSMATRIX (PROPS,NPROPS,C)
IMPLICIT NONE
INTEGER :: NPROPS, I, J
REAL*8 :: PROPS(NPROPS), C(6,6), C11, C12, C44
C11 = PROPS(1)
C12 = PROPS(2)
C44 = PROPS(3)
C INITIALIZE C
C = 0.
DO I=1,3
DO J=1,3
IF (I==J) THEN
C(I,J) = C11
ELSE
C(I,J) = C12
END IF
END DO
C(I+3,I+3) = C44
END DO
RETURN
END SUBROUTINE STIFFNESSMATRIX
SUBROUTINE DISLOCATIONDENSITY (STATEV,NSTATV,PROPS,NPROPS,
1 DTIME,GAMMA_DOT)
REAL*8::FCC_SLIP(12,9), MU, V_E0, V_S0, V_SP0, V_EP0, V_EN0,V_SN0
REAL*8 :: P, Q, R_E, R_S, S_EP, S_SP, F_E, BK, TEMP
REAL*8 :: TAU(12), RHO_EP(12), RHO_EN(12), RHO_SP(12), RHO_SN(12)
REAL*8 :: V_EP(12), V_EN(12), V_SP(12), V_SN(12)
REAL*8 :: N(12,3), B(12,3), T(12,3)
REAL*8 :: GBB, HBB, GAMMA_DOT(12), S_ED(12), S_SD(12)
REAL*8 :: L_SP(12), L_SN(12), L_EP(12), L_EN(12)
REAL*8 :: V_E(12), V_S(12), RHO_DOT_EGEN(12), RHO_DOT_SGEN(12)
REAL*8 :: RHO_DOT_GEN(12), RHO_DOT_EANN(12), RHO_DOT_SANN(12)
REAL*8 :: RHO_DOT_E(12), RHO_DOT_S(12), RHO_DOT_ANN(12)
REAL*8 :: RHO_DOT_FLUX(12), RHO_DOT(12)
REAL*8, dimension(12,12) :: G, H
REAL*8, parameter :: PI_8 = 4 * atan (1.0_8)
REAL*8 :: G0,G1,G2,G3,G4,G5,H0,H1,H2,H3,H4,H5
REAL*8 :: LG2,LG3,LG4,LG5,LH2,LH3,LH4,LH5,GHTERM
INTEGER :: ALPHA, BETA,I,J, NSTATV, NPROPS, NSLIP
REAL*8 :: STATEV(NSTATV), PROPS(NPROPS), DTIME
c Assign values (adjust as needed)
V_E0 = PROPS(25)
V_S0 = PROPS(26)
V_EP0 = PROPS(27)
V_EN0 = PROPS(28)
V_SP0 = PROPS(29)
V_SN0 = PROPS(30)
P_E = PROPS(41)
P_S = PROPS(42)
Q_E = PROPS(43)
Q_S = PROPS(44)
R_E = PROPS(45)
R_S = PROPS(46)
S_EP = PROPS(49)
S_SP = PROPS(50)
F_E = PROPS(47)
F_S = PROPS(48)
BK = 1.380649e-23
TEMP = PROPS(51)
BMAG = PROPS(52)
MU = PROPS(4)
C feed in n,b,t from statev
I = 1
DO ALPHA = 1, NSLIP
DO J = 1, 3
N(ALPHA,J) = STATEV(24+I)
B(ALPHA,J) = STATEV(60+I)
T(ALPHA,J) = STATEV(96+I)
IF (J .LT. 3) THEN
I = I + 1
END IF
END DO
I = I + 1
END DO
TAU = STATEV(1:12)
RHO_EP = STATEV(133:144)
RHO_EN = STATEV(145:156)
RHO_SP = STATEV(157:168)
RHO_SN = STATEV(169:180)
V_EP = STATEV(181:192)
V_EN = STATEV(193:204)
V_SP = STATEV(205:216)
V_SN = STATEV(217:228)
c Compute t
T = T/SQRT(6.)
c Compute G and Hnew matrices
do ALPHA = 1, 12
GBB = 0.
HBB = 0.
do BETA = 1, 12
c Define G and Hnew matrices elements (adjust as needed)
G0 = 0.10
G1 = 0.22
LG2 = 0.30
LG3 = 0.38
LG4 = 0.16
LG5 = 0.45
H0 = 0.0
H1 = 0.0
LH2 = 0.05
LH3 = 0.12
LH4 = 0.03
LH5 = 0.25
GHTERM = abs(dot_product(n(alpha,1:3),t(beta,1:3)))
C write (*,'(f0.4)') n(alpha,1:3)
G2 = LG2*GHTERM
G3 = LG3*GHTERM
G4 = LG4*GHTERM
G5 = LG5*GHTERM
H2 = LH2*GHTERM
H3 = LH3*GHTERM
H4 = LH4*GHTERM
H5 = LH5*GHTERM
C Populate G, H matrices
C G, H matrices are symmetric
G(1,:) = (/G0,G1,G1,G4,G5,G3,G4,G3,G5,G2,G3,G3/)
G(2,:) = (/G1,G0,G1,G5,G4,G3,G3,G2,G3,G3,G4,G5/)
G(3,:) = (/G1,G1,G0,G3,G3,G2,G5,G3,G4,G3,G5,G4/)
G(4,:) = (/G4,G5,G3,G0,G1,G1,G2,G3,G3,G4,G3,G5/)
G(5,:) = (/G5,G4,G3,G1,G0,G1,G3,G4,G5,G3,G2,G3/)
G(6,:) = (/G3,G3,G2,G1,G1,G0,G3,G5,G4,G5,G3,G4/)
G(7,:) = (/G4,G3,G5,G2,G3,G3,G0,G1,G1,G4,G5,G3/)
G(8,:) = (/G3,G2,G3,G3,G4,G5,G1,G0,G1,G5,G4,G3/)
G(9,:) = (/G5,G3,G4,G3,G5,G4,G1,G1,G0,G3,G3,G2/)
G(10,:) = (/G2,G3,G3,G4,G3,G5,G4,G5,G3,G0,G1,G1/)
G(11,:) = (/G3,G4,G5,G3,G2,G3,G5,G4,G3,G1,G0,G1/)
G(12,:) = (/G3,G5,G4,G5,G3,G4,G3,G3,G2,G1,G1,G0/)
H(1,:) = (/H0,H1,H1,H4,H5,H3,H4,H3,H5,H2,H3,H3/)
H(2,:) = (/H1,H0,H1,H5,H4,H3,H3,H2,H3,H3,H4,H5/)
H(3,:) = (/H1,H1,H0,H3,H3,H2,H5,H3,H4,H3,H5,H4/)
H(4,:) = (/H4,H5,H3,H0,H1,H1,H2,H3,H3,H4,H3,H5/)
H(5,:) = (/H5,H4,H3,H1,H0,H1,H3,H4,H5,H3,H2,H3/)
H(6,:) = (/H3,H3,H2,H1,H1,H0,H3,H5,H4,H5,H3,H4/)
H(7,:) = (/H4,H3,H5,H2,H3,H3,H0,H1,H1,H4,H5,H3/)
H(8,:) = (/H3,H2,H3,H3,H4,H5,H1,H0,H1,H5,H4,H3/)
H(9,:) = (/H5,H3,H4,H3,H5,H4,H1,H1,H0,H3,H3,H2/)
H(10,:) = (/H2,H3,H3,H4,H3,H5,H4,H5,H3,H0,H1,H1/)
H(11,:) = (/H3,H4,H5,H3,H2,H3,H5,H4,H3,H1,H0,H1/)
H(12,:) = (/H3,H5,H4,H5,H3,H4,H3,H3,H2,H1,H1,H0/)
Gbb = Gbb + G(alpha,beta)*(rho_ep(beta)+rho_en(beta))
1 + G(alpha,beta)*(rho_sp(beta)+rho_sn(beta))
Hbb = Hbb + H(alpha,beta)*(rho_ep(beta)+rho_en(beta))
1 + H(alpha,beta)*(rho_sp(beta)+rho_sn(beta))
end do
c Compute changes
gamma_dot(alpha) = (rho_ep(alpha)*v_ep(alpha)
1 + rho_en(alpha)*v_en(alpha)
2 + rho_sp(alpha)*v_sp(alpha) + rho_sn(alpha)*v_sn(alpha))
3 *bmag*sign(1.,tau(alpha))
s_ed(alpha) = mu*bmag*sqrt(Gbb)
s_sd(alpha) = s_ed(alpha)
l_sp(alpha) = 1 / sqrt(Hbb)
l_sn(alpha) = l_sp(alpha)
l_ep(alpha) = 1 / sqrt(Hbb)
l_en(alpha) = l_ep(alpha)
v_e(alpha) = v_e0*exp(-F_E/(bk*Temp)*
1 (1 - (abs(tau(alpha))/(s_ep+s_ed(alpha)))**P_E)**Q_E )
v_s(alpha) = v_s0*exp(-F_S/(bk*Temp)*
1 (1 - (abs(tau(alpha))/(s_sp+s_sd(alpha)))**P_S)**Q_S )
v_ep(alpha) = v_e(alpha)
v_en(alpha) = v_e(alpha)
v_sp(alpha) = v_s(alpha)
v_sn(alpha) = v_s(alpha)
rho_dot_egen(alpha) = 2*(rho_sp(alpha)*abs(v_sp(alpha))
1 /l_sp(alpha) + rho_sn(alpha)*v_sn(alpha)/l_sn(alpha))
rho_dot_sgen(alpha) = 2*(rho_ep(alpha)*abs(v_ep(alpha))
1 /l_ep(alpha) + rho_en(alpha)*v_en(alpha)/l_en(alpha))
rho_dot_gen(alpha) = rho_dot_egen(alpha) + rho_dot_sgen(alpha)
rho_dot_eann(alpha) = -2*rho_ep(alpha)*rho_en(alpha)
1 *R_e*(v_ep(alpha) + v_en(alpha))
rho_dot_sann(alpha) = -2*rho_sp(alpha)*rho_sn(alpha)
1 *R_s*(v_sp(alpha) + v_sn(alpha))
rho_dot_e(alpha) = rho_dot_egen(alpha) + rho_dot_eann(alpha)
rho_dot_s(alpha) = rho_dot_sgen(alpha) + rho_dot_sann(alpha)
rho_dot_ann(alpha) = rho_dot_eann(alpha) + rho_dot_sann(alpha)
rho_dot_flux(alpha) = 0.
rho_dot(alpha) = rho_dot_gen(alpha) + rho_dot_ann(alpha)
1 + rho_dot_flux(alpha)
rho_ep(alpha) = rho_ep(alpha) + DTIME*rho_dot_e(alpha)/2
rho_en(alpha) = rho_en(alpha) + DTIME*rho_dot_e(alpha)/2
rho_sp(alpha) = rho_sp(alpha) + DTIME*rho_dot_s(alpha)/2
rho_sn(alpha) = rho_sn(alpha) + DTIME*rho_dot_s(alpha)/2
end do
DO I = 1, NSLIP
STATEV(132+I) = RHO_EP(I)
STATEV(144+I) = RHO_EN(I)
STATEV(156+I) = RHO_SP(I)
STATEV(168+I) = RHO_SN(I)
STATEV(180+I) = V_EP(I)
STATEV(192+I) = V_EN(I)
STATEV(204+I) = V_SP(I)
STATEV(216+I) = V_SN(I)
END DO
RETURN
END
