“这么厉害的嘛！ ”
前言
今天有个北理工的博士，发邮件给我，问我：“许老师，我在看你的B-W 模型……”
正文
开源了解一下？
C
SUBROUTINE VUMAT
C Read only -
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 -
7 STRESSNEW, STATENEW, ENERINTERNNEW, ENERINELASNEW)
C
INCLUDE 'vaba_param_dp.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),DEFGRADOLD(NBLOCK,NDIR+NSHR+NSHR),
4 FIELDOLD(NBLOCK, NFIELDV), STRESSOLD(NBLOCK, NDIR+NSHR),
5 STATEOLD(NBLOCK, NSTATEV), ENERINTERNOLD(NBLOCK),
6 ENERINELASOLD(NBLOCK), TEMPNEW(NBLOCK),
7 STRETCHNEW(NBLOCK, NDIR+NSHR),DEFGRADNEW(NBLOCK,NDIR+NSHR+NSHR),
8 FIELDNEW(NBLOCK, NFIELDV), STRESSNEW(NBLOCK,NDIR+NSHR),
9 STATENEW(NBLOCK, NSTATEV), ENERINTERNNEW(NBLOCK),
1 ENERINELASNEW(NBLOCK)
CHARACTER*80 CMNAME
DOUBLE PRECISION DEQPLAS,HARD,SYIELD,SYIELD0,
+ S1, S2, S3, S4, S5, S6,dsmag,dsmag0, EQPLAS,
+ EQUSTRESS,FLOW1, FLOW2, FLOW3, FLOW4, FLOW5, FLOW6,
+ ETA, XI, BETA, THETA, GAMMA,
+ p_cp, p_eta0, p_cths, p_cthc, p_m,factor,
+ fmag, ext_p_th, ext_p_th1, fsum, lodetol
c
DIMENSION EELAS(6),EPLAS(6),FLOW(6),olds(6)
C
PARAMETER(ZERO=0.0D0, ONE=1.D0, TWO=2.D0, THREE=3.D0, SIX=6.D0,
1 ENUMAX=.4999D0, NEWTON=150,HALF=0.5D0,
2 con_gamma=6.4641, UTHETA=6.44641D0,PI=3.14159D0 )
C ----------------------------------------------------------------
C PROPS(1) - E
C PROPS(2) - NU
C PROPS(3..15) - SYIELD AN HARDENING DATA
C CALLS UHARD FOR CURVE OF YIELD STRESS VS. PLASTIC STRAIN
C ----------------------------------------------------------------
PI6 = PI/6.0
FLAG = ZERO
DEQPLAS = ZERO
eqplas = zero
ETATOL = zero
ext_p_th = 1.0
fsum = 1.0
x0 = zero
factor = zero
dgamma = zero
const = sqrt(THREE/TWO)
c
EMOD = PROPS(1)
ENU = MIN(PROPS(2), ENUMAX)
p_cp = PROPS(10)
p_eta0 = PROPS(11)
p_cths = PROPS(12)
p_cthc = PROPS(13)
p_m = PROPS(14)
EBULK3 = EMOD/(ONE-TWO*ENU)
EG2 = EMOD/(ONE+ENU)
EG = EG2/TWO
EG3 = THREE*EG
ELAM = (EBULK3-EG2)/THREE
c
c
IF ( STEPTIME .EQ. ZERO ) THEN
DO K = 1, NBLOCK
TRACE = STRAININC(K,1) + STRAININC(K,2) + STRAININC(K,3)
STRESSNEW(K,1) = STRESSOLD(K,1)
* + EG2 * STRAININC(K,1) + ELAM * TRACE
STRESSNEW(K,2) = STRESSOLD(K,2)
* + EG2 * STRAININC(K,2) + ELAM * TRACE
STRESSNEW(K,3) = STRESSOLD(K,3)
* + EG2 * STRAININC(K,3) + ELAM * TRACE
STRESSNEW(K,4)=STRESSOLD(K,4) + EG2 * STRAININC(K,4)
IF ( NSHR .GT. 1 ) THEN
STRESSNEW(K,5)=STRESSOLD(K,5) + EG2 * STRAININC(K,5)
STRESSNEW(K,6)=STRESSOLD(K,6) + EG2 * STRAININC(K,6)
END IF
END DO
ELSE
DO K = 1, NBLOCK
TRACE = STRAININC(K,1) + STRAININC(K,2) + STRAININC(K,3)
S1 = STRESSOLD(K,1)+ EG2 * STRAININC(K,1) + ELAM * TRACE
S2 = STRESSOLD(K,2)+ EG2 * STRAININC(K,2) + ELAM * TRACE
S3 = STRESSOLD(K,3)+ EG2 * STRAININC(K,3) + ELAM * TRACE
S4 = STRESSOLD(K,4)+ EG2 * STRAININC(K,4)
S5 = ZERO
S6 = ZERO
IF ( NSHR .GT. 1 ) THEN
S5 = STRESSOLD(K,5) + EG2 * STRAININC(K,5)
S6 = STRESSOLD(K,6) + EG2 * STRAININC(K,6)
END IF
* Deviatoric part of trial stress
smean = ( s1 + s2 + s3 ) / three
ds1 = s1 - smean
ds2 = s2 - smean
ds3 = s3 - smean
* Magnitude of the deviatoric trial simses stress
dsmag0 = const*sqrt ( ds1*ds1 + ds2*ds2 + ds3*ds3
1 + two*s4*s4 + two*s5*s5 + two*s6*s6 )
* read for the equivalent plastic strain
EQPLAS = STATEOLD(K,1+2*(NDIR+NSHR))
ETATOL = STATEOLD(K,15)
LODETOL= STATEOLD(K,16)
* read the old stress state in time t for calculate flow direction (explicit)
do i=1,6
olds(i) = zero
end do
c
do i=1,ndir+nshr
olds(i) = STRESSOLD(k,i)
end do
* CALCULATE THE CURRENT YIELD STRESS
CALL CAL_SYIELD( olds(1), olds(2), olds(3), olds(4),olds(5),
+ olds(6), EQPLAS, PROPS(3), HARD, SYIELD0,
+ ext_p_th, p_cp, p_eta0, p_cths, p_cthc, p_m )
c
c write(99,*) dsmag0,ext_p_th*SYIELD0,eqplas
if(dsmag0.gt.ext_p_th*SYIELD0) then
c plastic yielding actually £¬determine the flow direction £¨explicitly£©
SYIELD = SYIELD0
call CAL_DFLOW( olds(1), olds(2), olds(3), olds(4),
+ olds(5), olds(6), eqplas,SYIELD,
+ FLOW1, FLOW2, FLOW3, FLOW4, FLOW5, FLOW6,
+ ext_p_th, p_cp, p_eta0, p_cths, p_cthc, p_m)
fsum = sqrt( FLOW1**2 +FLOW2**2 + FLOW3**2
+ 2.*FLOW4**2 +2.* FLOW5**2 +2.*FLOW6**2 )
c
c determine dgamma and plastic correction
call cal_pp_explicit (S1, S2, S3, S4, S5, S6,
+ EG2, EXT_P_TH, SYIELD,dSmag0,
+ flow1, flow2, flow3, flow4, flow5, flow6, dgamma)
else
dgamma = zero
flow1 = zero
flow2 = zero
flow3 = zero
flow4 = zero
flow5 = zero
flow6 = zero
end if
if (dgamma.lt.zero) dgamma = zero
deqplas = fsum*dgamma/const
factor = EG2 * dgamma
c write(100,*) n_iter,factor
stressNew(k,1) = s1 - factor * FLOW1
stressNew(k,2) = s2 - factor * FLOW2
stressNew(k,3) = s3 - factor * FLOW3
stressNew(k,4) = s4 - factor * FLOW4
if(nshr.gt.1) then
stressNew(k,5) = s5 - factor * FLOW5
stressNew(k,6) = s6 - factor * FLOW6
end if
c calculate the equstress
SMISES =(STRESSNEW(K,1)-STRESSNEW(K,2))**2 +
1 (STRESSNEW(K,2)-STRESSNEW(K,3))**2 +
2 (STRESSNEW(K,3)-STRESSNEW(K,1))**2
DO K1=NDIR+1,NDIR+NSHR
SMISES=SMISES+SIX*STRESSNEW(K,K1)**2
ENDDO
SMISES = SQRT(SMISES/TWO)
C Update equivalent plastic strain and other state variable
STATENEW(K,1+2*(NDIR+NSHR)) = STATEOLD(K,1+2*(NDIR+NSHR))
+ + deqplas
C------------------------------------------------------------
C calculate the pressure, stress triaxiality,lode angle
SMEAN = ( STRESSNEW(K,1)+ STRESSNEW(K,2) +
+ STRESSNEW(K,3))/ 3.0
ST5 = 0.0D0
ST6 = 0.0D0
IF ( NSHR .GT. 1 ) THEN
ST5 = STRESSNEW(K,5)
ST6 = STRESSNEW(K,6)
END IF
C stress triaxiality: eta
IF (SMISES .GT. ZERO) THEN
ETA = SMEAN/ SMISES
ETATOL = ETATOL + ETA*deqplas
ELSE
ETA = ZERO
END IF
STATENew(K,14) = ETA
STATENEW(K,15) = ETATOL
C lode parameter: gamma and beta
CALL CAL_DSJ3(stressNew(K,1), stressNew(K,2),stressNew(K,3),
1 stressNew(K,4), ST5, ST6, DSJ3)
IF (SMISES .GT. ZERO) THEN
XI = 13.5*DSJ3/(SMISES**3)
ELSE
XI = ZERO
END IF
BETA = XI
c
IF (XI .GT. ONE) XI= ONE
IF (XI .LT. -ONE) XI=-ONE
THETA=ACOS(XI)/THREE
STATENew(K,17) = THETA
STATENew(K,18) = 1 - THETA/PI6
GAMMA=CON_GAMMA*(ONE/COS(THETA-PI6)-ONE)
LODETOL = LODETOL + STATENew(K,18)*deqplas
STATENEW(K,16) = LODETOL
STATENEW(K,20) = GAMMA
C
c CALCULATE THE AVG ETA AND NORMALIZED LODE ANGEL
if (STATENEW(K,1+2*(NDIR+NSHR)).gt. zero) then
STATENEW(K,21) = ETATOL/STATENEW(K,1+2*(NDIR+NSHR))
STATENEW(K,22) = LODETOL/STATENEW(K,1+2*(NDIR+NSHR))
end if
C ------------------------------------------
C CALCULATE PLASTIC DISSIPATION
SPD=DEQPLAS*(SYIELD0+SYIELD)/TWO
c
IF ( NSHR .EQ. 1 ) THEN
STRESSPOWER = HALF * (
* ( STRESSOLD(K,1) + STRESSNEW(K,1) ) * STRAININC(K,1) +
* ( STRESSOLD(K,2) + STRESSNEW(K,2) ) * STRAININC(K,2) +
* ( STRESSOLD(K,3) + STRESSNEW(K,3) ) * STRAININC(K,3))+
* ( STRESSOLD(K,4) + STRESSNEW(K,4) ) * STRAININC(K,4)
ELSE
STRESSPOWER = HALF * (
* ( STRESSOLD(K,1) + STRESSNEW(K,1) ) * STRAININC(K,1) +
* ( STRESSOLD(K,2) + STRESSNEW(K,2) ) * STRAININC(K,2) +
* ( STRESSOLD(K,3) + STRESSNEW(K,3) ) * STRAININC(K,3))+
* ( STRESSOLD(K,4) + STRESSNEW(K,4) ) * STRAININC(K,4) +
* ( STRESSOLD(K,5) + STRESSNEW(K,5) ) * STRAININC(K,5) +
* ( STRESSOLD(K,6) + STRESSNEW(K,6) ) * STRAININC(K,6)
END IF
ENERINTERNNEW(K) = ENERINTERNOLD(K) +
& STRESSPOWER / DENSITY(K)
C UPDATE THE DISSIPATED INELASTIC SPECIFIC ENERGY -
C
PLASTICWORKINC = HALF * ( SYIELD0 + SYIELD ) * DEQPLAS
ENERINELASNEW(K) = ENERINELASOLD(K) +
& PLASTICWORKINC / DENSITY(K)
END DO
END IF
C
RETURN
END
c this subroutine is calculate the product of plastic (explicit)
SUBROUTINE cal_pp_explicit (S1, S2, S3, S4, S5, S6,
+ EG2, EXT_P_TH, SYIELD,dSmag0,
+ flow1, flow2, flow3, flow4, flow5, flow6, dgamma)
c
INCLUDE 'vaba_param_dp.inc'
c--------------------------------------
CHARACTER*80 CMNAME
DOUBLE PRECISION S1, S2, S3, S4, S5, S6,
+ flow1, flow2, flow3, flow4, flow5, flow6,
+ EG2, EXT_P_TH, SYIELD,dSmag0, fmag,
+ temp1,temp2, a, b, c, dgamma,CONST
c
CONST=SQRT(2.0/3.0)
C
call CAL_EQUSTRESS(FLOW1, FLOW2,FLOW3,FLOW4, FLOW5,
+ FLOW6, fmag)
a = 2.*fmag**2*(EG2**2)
temp1 = (s1-s2)*(flow1-flow2)+(s2-s3)*(flow2-flow3)
+ + (s3-s1)*(flow3-flow1)
temp2 = 6.*(s4*flow4 + s5*flow5+ s6*flow6)
b = -2.*EG2*(temp1 + temp2)
c = 2.*dSmag0**2 - 2.* (SYIELD**2)*(EXT_P_TH**2)
c slove the quadratic equation with one unknown for calulate dgamma
dgamma = -b/2./a - sqrt(b**2-4.*a*c)/2./a
c write(97,*) a, b, c ,dgamma
return
end
c SUBROUTINE for calculate the harding parameter and yield stress
SUBROUTINE CAL_SYIELD( S1, S2, S3, S4, S5, S6,
+ EP, TABLE, HARD,SYIELD1,
+ EXT_P_TH, p_cp, p_eta0, p_cths, p_cthc, p_m )
c
C NVALUE,
INCLUDE 'vaba_param_dp.inc'
c--------------------------------------
CHARACTER*80 CMNAME
DOUBLE PRECISION A0,HARD,B,EN,A1,EN1,A2,EN2,SYIELD1,
+ S1, S2, S3, S4, S5, S6, equStress,EP,TABLE(7),
+ p_cp, p_eta0, p_cths, p_cthc, p_cthax, p_m,
+ DSJ3, THETA, GAMMA, ETA, con_gamma, pi6, XI,
+ ext_p, ext_th, ext_p_th
C
PARAMETER(ZERO=0.D0, con_gamma = 6.4641, pi6 = 0.523599)
c------------------------------------------------------
A0 = TABLE(1)
B = TABLE(2)
EN = TABLE(3)
A1 = TABLE(4)
EN1 = TABLE(5)
A2 = TABLE(6)
EN2 = TABLE(7)
c
ext_p_th = 1.0d0
HARD = A1*EN1 + A2*EN2
C -------------------------------------------------------
SM=(S1+ S2 + S3)/3.0
CALL CAL_EQUSTRESS(S1, S2, S3, S4, S5, S6, equStress)
if (equStress .gt. 0.0) then
ETA=SM/equStress
else
ETA=0.0
end if
C -------------------------------------------------------
CALL CAL_DSJ3(S1, S2, S3, S4, S5, S6, DSJ3)
if (equStress .gt. 0.0) then
XI = 13.5*DSJ3/(equStress**3)
else
XI = 0.0
end if
if (XI .gt. 1.0) XI=1.0
if (XI .lt. -1.0) XI=-1.0
THETA=ACOS(XI)/3.0
GAMMA=con_gamma*(1.0/cos(THETA-pi6)-1.0)
if (XI .ge. 0.0) then
p_cthax=1.0
else
p_cthax=p_cthc
end if
C -----------
ext_p = 1.0 - p_cp*(ETA-p_eta0)
ext_th= p_cths + (p_cthax-p_cths)
+ *(1.0+1.0/p_m)*(GAMMA-GAMMA**(p_m+1.0)/(p_m+1.0))
c PROTECT CONDITION !!!!!!!!
if (ext_p .ge. 2.0) then
ext_p = 2.0
else if (ext_p .le. 0.2) then
ext_p = 0.2
end if
if (ext_th .ge. 2.0) then
ext_th = 2.0
else if (ext_th .le. 0.2) then
ext_th = 0.2
end if
ext_p_th = ext_p * ext_th
c
if (ext_p_th .ge. 2.0) then
ext_p_th = 2.0
else if (ext_p_th .le. 0.2) then
ext_p_th = 0.2
end if
C CURRENT YIELD STRESS AND HARDENING
IF(EP.EQ.0.0) THEN
SYIELD1 = A0
c ext_p_th = 1.0
ELSE
HARD= EN*B*EP**(EN-1) +
& A1*EN1*DEXP(-EN1*EP) +
& A2*EN2*DEXP(-EN2*EP)
SYIELD1= A0 + B*EP**EN +
& A1*(1-DEXP(-EN1*EP)) +
& A2*(1-DEXP(-EN2*EP))
ENDIF
c write(101,*) syield1, ep
RETURN
END
C------------------------------------------------------------------
C SUBROUTINE for calculate the third invariant
SUBROUTINE CAL_DSJ3(S1, S2, S3, S4, S5, S6, DSJ3)
C --
C -- THIS SUBROUTINE CALCULATES THE THIRD STRESS PARAMETER: J3 (DSJ3)
C -- OF A 3 BY 3 STRESS MATRIX [SIG].
C
C VARIABLES
C
DOUBLE PRECISION S1, S2, S3, S4, S5, S6, SM, DSIG(3,3), DSJ3
SM=(S1+S2+S3)/3.0
DSIG(1,1)=S1-SM
DSIG(2,2)=S2-SM
DSIG(3,3)=S3-SM
DSIG(1,2)=S4
DSIG(2,3)=S5
DSIG(3,1)=S6
C
DSIG(2,1)=S4
DSIG(3,2)=S5
DSIG(1,3)=S6
CALL MDET(DSIG, DSJ3)
RETURN
END
C
C ------------------------------------------------------------------
C
SUBROUTINE MDET(A,DET)
C --
C -- THIS SUBROUTINE CALCULATES THE DETERMINANT
C -- OF A 3 BY 3 MATRIX [A].
C
C VARIABLES
C
DOUBLE PRECISION A(3,3), DET
C
C COMPUTATION
C
DET = A(1,1)*A(2,2)*A(3,3)
+ + A(1,2)*A(2,3)*A(3,1)
+ + A(1,3)*A(2,1)*A(3,2)
+ - A(3,1)*A(2,2)*A(1,3)
+ - A(3,2)*A(2,3)*A(1,1)
+ - A(3,3)*A(2,1)*A(1,2)
RETURN
END
C ------------------------------------------------------------------
SUBROUTINE CAL_EQUSTRESS(S1, S2, S3, S4, S5, S6, equStress)
C --
DOUBLE PRECISION S1, S2, S3, S4, S5, S6, SM, equStress
SM = (S1+S2+S3)/3.0
DS1 = S1-SM
DS2 = S2-SM
DS3 = S3-SM
equStress = sqrt(1.5 *( DS1**2+ DS2**2+ DS3**2
+ + 2.0 * S4**2 + 2.0* S5**2 + 2.0 * S6**2))
RETURN
END
C
C ------------------------------------------------------------------
SUBROUTINE CAL_DFLOW(S1, S2, S3, S4, S5, S6, EP,SYIELD1,
+ FLOW1, FLOW2, FLOW3, FLOW4, FLOW5, FLOW6,
+ ext_p_th, p_cp, p_eta0, p_cths, p_cthc, p_m)
C
C This subroutine is to calculate the flow direction
C
C
DOUBLE PRECISION S1, S2, S3, S4, S5, S6, SM, DS(3,3), equStress,
+ FLOW1, FLOW2, FLOW3, FLOW4, FLOW5, FLOW6, EP,SYIELD1,
+ p_cp, p_eta0, p_cths, p_cthc, p_cthax, p_m, fmag,
+ DSJ3, COS3TH, SIN3TH, THETA, GAMMA, ETA, con_gamma, pi6, pi3
+ C1, C2, C3, A1, A2, A3, B1, B2, F(3,3), SS, QMAG, ext_p_th,
+ DELTA, XI
con_gamma = 6.4641
pi6 = 0.523599
pi3 = 1.047198
C
SM=(S1+S2+S3)/3.0
DS(1,1)=S1-SM
DS(2,2)=S2-SM
DS(3,3)=S3-SM
DS(1,2)=S4
DS(2,3)=S5
DS(3,1)=S6
DS(2,1)=S4
DS(3,2)=S5
DS(1,3)=S6
C
CALL CAL_EQUSTRESS(S1, S2, S3, S4, S5, S6, equStress)
if (equStress .le. 0.0) then
equStress=1.0
ETA = 0.0
else
ETA=SM/equStress
end if
C ------------
CALL CAL_DSJ3(S1, S2, S3, S4, S5, S6, DSJ3)
if (equStress .gt. 0.0) then
XI = 13.5*DSJ3/(equStress**3)
else
XI = 0.0
end if
if (XI .gt. 1.0) XI=1.0
if (XI .lt. -1.0) XI=-1.0
COS3TH=XI
THETA=ACOS(XI)/3.0
if ((THETA .le. 0.0) .or. (THETA .ge. pi3)) then
SIN3TH=0.00000001
else
SIN3TH=SIN(3.0*THETA)
end if
GAMMA=con_gamma*(1.0/cos(THETA-pi6)-1.0)
if (XI .ge. 0.0) then
p_cthax=1.0
else
p_cthax=p_cthc
end if
c ------------
C1=3.0/2.0/equStress
C2=SYIELD * p_cp *(1.0+1.0/p_m)*(p_cths +(p_cthax-p_cths)
+ *(GAMMA-(GAMMA**(p_m+1.0))/(p_m+1.0)))
C3=SYIELD * (1.0 - p_cp*(ETA-p_eta0)) * ((p_cthax-p_cths)
+ *(1.0-GAMMA**p_m)) *(1.0+1.0/p_m)* con_gamma * tan(THETA-pi6)
+ /cos(THETA-pi6) * 3.0/(equStress*SIN3TH)
A1=1.0/3.0
A2=COS3TH/2.0/equStress
A3=3.0/2.0/(equStress**2)
B1=1.0/3.0/equStress
B2=3.0*SM/2.0/(equStress**3)
C
C ------------
DO 300 II=1, 3
DO 300 JJ=II, 3
C
IF (II .eq. JJ) THEN
DELTA=1.0
ELSE
DELTA=0.0
END IF
C
SS=0.0
DO 310 KK=1,3
SS=SS+DS(II,KK)*DS(KK,JJ)
310 CONTINUE
C
F(II,JJ)=C1*DS(II,JJ) -C2*B2*DS(II,JJ)
+ -C3*(A1*DELTA + A2*DS(II,JJ) - A3*SS)
300 CONTINUE
C
C -------------
CALL CAL_EQUSTRESS(F(1,1), F(2,2), F(3,3), F(1,2),
+ F(2,3), F(1,3), fMAG)
if (fMAG .eq. 0.0) fMAG=1.0
FLOW1 = F(1,1) / fMAG
FLOW2 = F(2,2) / fMAG
FLOW3 = F(3,3) / fMAG
FLOW4 = F(1,2) / fMAG
FLOW5 = F(2,3) / fMAG
FLOW6 = F(1,3) / fMAG
C -------------
C
RETURN
END
