Dear all,
I am struggling with the VUMAT. I want to implement a strain rate dependence model regardless of the physical sense first. I suppose the strain rate can be calculate by STRAININC/DT. Its 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 wasnt working if using straininc(np,4)/dt in the if() then. I mean if the shear strain rate is calculated from straininc(np,4)/dt. But I found there is another variable named strainrate. I dont know what it is, maybe its an internal variable. If I use the variable named strainrate directly, the condition is working. I feel confused about the phenomenon. I want to track the straininc(np,4) and strainrate and check what happens. But I dont know how to state the two variables. I want to output the variables for all time steps. Could you 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