strain rate in VUMAT

Dear all,



I am struggling with the VUMAT subroutine. I want to implement a strain rate dependence model regardless of the physical sense first. I suppose that the strain rate can be calculate by STRAININC/DT. It's supposed that the shear stiffness component depends on the shear strain rate (it equals straininc(np,4)/dt) if the strain rate is higher than 10 s-1. But now I found it wasn't working if using "straininc(np,4)/dt" in the condition "if() then". I mean if the shear strain rate is calculated from straininc(np,4)/dt, the VUMAT code wasn's working. But I found there is another variable named strainrate. I don't know what it is, maybe it's an internal variable. If I use the variable named "strainrate" directly as the consition, the below VUMAT code is working. I feel confused about the issue. I want to track the "straininc(np,4)/dt" and "strainrate" and check what happens. But I don't know how to state the two variables. In fact, I have tired "write(*,*)", but I cannot find the output. Maybe I made a mistatke when I use the "write". Now I want to output the variables for all time steps and to trach the "SR" in the below code. Could anybody help me? Or do you know anything about the two variables? Thanks a million.



Haibin





I copy my code as follows.

C************************************************* *********

C VUMAT FOR STRAIN RATE DEPENDENT STIFFNESS MODEL, ORTHOTROPIC MATERIELS

C WRITTEN BY HAIBIN ZHU

C************************************************* *********

SUBROUTINE VUMAT(

C READ ONLY (UNMODIFIABLE)VARIABLES -

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,

5 STRESSOLD, STATEOLD, ENERINTERNOLD, ENERINELASOLD,

6 TEMPNEW, STRETCHNEW, DEFGRADNEW, FIELDNEW,

C WRITE ONLY (MODIFIABLE) VARIABLES -

7 STRESSNEW, STATENEW, ENERINTERNNEW, ENERINELASNEW )

C

INCLUDE 'VABA_PARAM.INC'

C

DIMENSION PROPS(NPROPS), DENSITY(NBLOCK), COORDMP(NBLOCK,*),

1 CHARLENGTH(NBLOCK), STRAININC(NBLOCK,NDIR+NSHR),

2 RELSPININC(NBLOCK,NSHR), TEMPOLD(NBLOCK),

3 STRETCHOLD(NBLOCK,NDIR+NSHR),

4 DEFGRADOLD(NBLOCK,NDIR+NSHR+NSHR),

5 FIELDOLD(NBLOCK,NFIELDV), STRESSOLD(NBLOCK,NDIR+NSHR),

6 STATEOLD(NBLOCK,NSTATEV), ENERINTERNOLD(NBLOCK),

7 ENERINELASOLD(NBLOCK), TEMPNEW(NBLOCK),

8 STRETCHNEW(NBLOCK,NDIR+NSHR),

8 DEFGRADNEW(NBLOCK,NDIR+NSHR+NSHR),

9 FIELDNEW(NBLOCK,NFIELDV),

1 STRESSNEW(NBLOCK,NDIR+NSHR), STATENEW(NBLOCK,NSTATEV),

2 ENERINTERNNEW(NBLOCK), ENERINELASNEW(NBLOCK)

C

CHARACTER*80 CMNAME

C

PARAMETER( ZERO = 0.D0, ONE = 1.D0, TWO = 2.D0, THREE = 3.D0,

1 THIRD = ONE/THREE, HALF = .5, TWOTHIRDS = TWO/THREE,

2 THREEHALFS = 1.5D0 )

Q11 = PROPS(1)

Q12 = PROPS(2)

Q22 = PROPS(3)

Q13 = PROPS(4)

Q23 = PROPS(5)

Q33 = PROPS(6)

A = 7E9

B = 7E4

C

C

DO NP = 1, NBLOCK

C Copy the values of strains into temp variables



ES1 = strainInc(NP,1)

ES2 = strainInc(NP,2)

ES3 = strainInc(NP,3)

ES4 = strainInc(NP,4)





C COMPUTE THE STRAIN RATE DEPEDENT STIFFNESS AND LAM CONSTANTS

SR = strainrate

C SR = strainInc(NP,4)/DT

C IF THE STRAINRATE IS NOT 0, Q22 IS STRAINRATE, OTHERWISE, Q22 IS 10E9

IF (SR.GT.10) THEN

Q66 = A+B*(LOG(SR))

C stressNew(NP,1) = stressOld(NP,1)+Q11*ES1+Q12*ES2+Q13*ES3

C stressNew(NP,2) = stressOld(NP,2)+Q12*ES1+Q22*ES2+Q23*ES3

C stressNew(NP,3) = stressOld(NP,3)+Q13*ES1+Q23*ES2+Q33*ES3

C stressNew(NP,4) = stressOld(NP,4)+Q66*ES4*2

ELSE

Q66 = 7E9

END IF

C

stressNew(NP,1) = stressOld(NP,1)+Q11*ES1+Q12*ES2+Q13*ES3

stressNew(NP,2) = stressOld(NP,2)+Q12*ES1+Q22*ES2+Q23*ES3

stressNew(NP,3) = stressOld(NP,3)+Q13*ES1+Q23*ES2+Q33*ES3

stressNew(NP,4) = stressOld(NP,4)+Q66*ES4*2

C



C

C-----UPDATE THE SPECIFIC INTERNAL ENERGY--------

C

STRESSPOWER = HALF *

1 ( ( STRESSOLD(NP,1)+STRESSNEW(NP,1) )*STRAININC(NP,1)

2 + ( STRESSOLD(NP,2)+STRESSNEW(NP,2) )*STRAININC(NP,2)

3 + ( STRESSOLD(NP,3)+STRESSNEW(NP,3) )*STRAININC(NP,3)

4 + ( STRESSOLD(NP,4)+STRESSNEW(NP,4) )*STRAININC(NP,4)*TWO )

C

ENERINTERNNEW(NP) = ENERINTERNOLD(NP)

1 + STRESSPOWER / DENSITY(NP)

C

C-----UPDATE THE DISSIPATED INELASTIC SPECIFIC ENERGY--------

C

PLASTICWORKINC = 0

C

ENERINELASNEW(NP) = ENERINELASOLD(NP)

1 + PLASTICWORKINC / DENSITY(NP)

C

END DO

C

RETURN

END

C

C