subroutine vumat(! Read only - . nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, . stepTime, totalTime, dt, cmname, coordMp, charLength, . props, density, strainInc, relSpinInc, . tempOld, stretchOld, defgradOld, fieldOld, . stressOld, stateOld, enerInternOld, enerInelasOld, . tempNew, stretchNew, defgradNew, fieldNew,! Write only - . stressNew, stateNew, enerInternNew, enerInelasNew) include 'vaba_param.inc' DIMENSION PROPS(NPROPS), DENSITY(NBLOCK), COORDMP(NBLOCK), . CHARLENGTH(NBLOCK), STRAININC(NBLOCK, NDIR+NSHR), . RELSPININC(NBLOCK, NSHR), TEMPOLD(NBLOCK), . STRETCHOLD(NBLOCK, NDIR+NSHR),DEFGRADOLD(NBLOCK,NDIR+NSHR+NSHR), . FIELDOLD(NBLOCK, NFIELDV), STRESSOLD(NBLOCK, NDIR+NSHR), . STATEOLD(NBLOCK, NSTATEV), ENERINTERNOLD(NBLOCK), . ENERINELASOLD(NBLOCK), TEMPNEW(NBLOCK), . STRETCHNEW(NBLOCK, NDIR+NSHR),DEFGRADNEW(NBLOCK,NDIR+NSHR+NSHR), . FIELDNEW(NBLOCK, NFIELDV), STRESSNEW(NBLOCK,NDIR+NSHR), . STATENEW(NBLOCK, NSTATEV), ENERINTERNNEW(NBLOCK), . ENERINELASNEW(NBLOCK)! CHARACTER*8 CMNAME!!-----------------------------------------------------------------------! Internal VUMAT variables!----------------------------------------------------------------------- REAL*8, PARAMETER :: PI = 3.1415927 parameter ( zero = 0.d0, one = 1.d0, two = 2.d0, . third = 1.d0 / 3.d0, half = 0.5d0, op5 = 1.5d0) REAL*8 SigDevOld(6) ! Old deviatoric stress REAL*8 SigDevNew(6) ! New deviatoric stress REAL*8 kappa, Eff_Void REAL*8 stresstr(6)!!!-----------------------------------------------------------------------! Read parameters from ABAQUS material card!----------------------------------------------------------------------- YOUNG = props(1) ! Elastic/Youngs Modulus (e.g: E=210MPa for steel) POISS = props(2) ! Poisson's ratio (often: \mu=0.33) SIGMA0 = props(3) ! Initial Yield stress (at EpsPEqv=0) BETA1 = props(4) ! Modified Gursons 1st constant: q_1 BETA2 = props(5) ! Modified Gursons 2nd constant: q_2 BETA3 = props(6) ! Modified Gursons 3rd constant: q_3 TOL = props(8) ! TOLerance on the plastic return mapping MXITER = props(9) ! Maximum number of iterations in the return map T_1 = props(10) ! Q_1 = props(11) ! 3 term Voce hardening rule... T_2 = props(12) ! T_i * (1 - EXP(-Q_i*EpsPEqv)) Q_2 = props(13) ! ... "" ... T_3 = props(14) Q_3 = props(15) FN = props(16) ! f_N, the volume fraction of void nucleating particles. SN = props(17) ! S_N, i.e the standard deviation in the nucleation distribution. EPSN = props(18) ! \epsilon_N, i.e mean plastic strain for the distribution. f_c = PROPS(19) ! Void volume fraction upon which coalescence occurs. f_f = PROPS(20) ! Critical void volume fraction upon which failure occurs. k_s = PROPS(21) ! Shear damage term, k_s=0 means no contribution! ! !-----------------------------------------------------------------------! Compute elasticity matrix!----------------------------------------------------------------------- G = 0.5*YOUNG/(1.0+POISS) K = YOUNG/(3.0*(1.0-2.0*POISS))! Elastic constants: twomu = YOUNG / ( one + POISS ) alamda = POISS * twomu / ( one - two * POISS ) thremu = op5 * twomu!!---- Calculate the preprocessor step -------------------------------------------- IF ((totalTime.eq.0.0).and.(stepTime.eq.0.0)) THEN DO km = 1, NBLOCK trace = strainInc(km,1) + strainInc(km,2) + strainInc(km,3) stressNew(km,1) = stressOld(km,1) . + twomu * strainInc(km,1) + alamda * trace stressNew(km,2) = stressOld(km,2) . + twomu * strainInc(km,2) + alamda * trace stressNew(km,3) = stressOld(km,3) . + twomu * strainInc(km,3) + alamda * trace stressNew(km,4)=stressOld(km,4) + twomu * strainInc(km,4) IF ( nshr .gt. 1 ) THEN stressNew(km,5)=stressOld(km,5) + twomu * strainInc(km,5) stressNew(km,6)=stressOld(km,6) + twomu * strainInc(km,6) END IF END DO ELSE! !--------------------------------------------------------------------------------------! Continue with the return mapping algorithm:!-------------------------------------------------------------------------------------- DO km = 1, nblock trace = strainInc(km,1) + strainInc(km,2) + strainInc(km,3) stresstr(1) = stressOld(km,1) . + twomu * strainInc(km,1) + alamda * trace stresstr(2) = stressOld(km,2) . + twomu * strainInc(km,2) + alamda * trace stresstr(3) = stressOld(km,3) . + twomu * strainInc(km,3) + alamda * trace stresstr(4)=stressOld(km,4) + twomu * strainInc(km,4) IF ( nshr .gt. 1 ) THEN stresstr(5)=stressOld(km,5) + twomu * strainInc(km,5) stresstr(6)=stressOld(km,6) + twomu * strainInc(km,6) END IF ! !-----------------------------------------------------------------------! Read the void content either from material card or memory:!----------------------------------------------------------------------- IF (totalTime.eq.dt) THEN !read the void content from material properties, as it is the first increment in the simulation VoidTr = PROPS(7) ELSE VoidTr = STATEOLD(km, 1) END IF! !-----------------------------------------------------------------------! Calculate pTr!----------------------------------------------------------------------- SigHNew = (stresstr(1)+stresstr(2)+stresstr(3))/3.0 pTr = -SigHNew ! !-----------------------------------------------------------------------! Calculate qTr!----------------------------------------------------------------------- SigDevNew(1) = stresstr(1)-SigHNew SigDevNew(2) = stresstr(2)-SigHNew SigDevNew(3) = stresstr(3)-SigHNew SigDevNew(4) = stresstr(4) IF ( nshr .gt. 1 ) THEN SigDevNew(5) = stresstr(5) SigDevNew(6) = stresstr(6) END IF IF ( nshr .eq. 1 ) THEN SigEqv = sqrt((3.0/2.0)*(SigDevNew(1)*SigDevNew(1) . +SigDevNew(2)*SigDevNew(2) . +SigDevNew(3)*SigDevNew(3)) . +3.0*(SigDevNew(4)*SigDevNew(4))) ELSE SigEqv = sqrt((3.0/2.0)*(SigDevNew(1)*SigDevNew(1) . +SigDevNew(2)*SigDevNew(2) . +SigDevNew(3)*SigDevNew(3)) . +3.0*(SigDevNew(4)*SigDevNew(4) . +SigDevNew(5)*SigDevNew(5) . +SigDevNew(6)*SigDevNew(6))) END IF qTr = SigEqv !!-----------------------------------------------------------------------! Third deviatoric stress invariant!----------------------------------------------------------------------- IF ( nshr .eq. 1 ) THEN J_3 = SigDevNew(1)*SigDevNew(2)*SigDevNew(3) . - SigDevNew(4)*SigDevNew(4)*SigDevNew(3) ELSE IF (nshr .gt. 1) THEN J_3 = SigDevNew(1)*SigDevNew(2)*SigDevNew(3) . - SigDevNew(1)*SigDevNew(5)*SigDevNew(5) . - SigDevNew(4)*SigDevNew(4)*SigDevNew(3) . + 2*SigDevNew(4)*SigDevNew(5)*SigDevNew(6) . - SigDevNew(2)*SigDevNew(6)*SigDevNew(6) END IF ! !-----------------------------------------------------------------------! State variable from previous time step!----------------------------------------------------------------------- EpsPEqvTr = STATEOLD(km,2)! !-----------------------------------------------------------------------! Update yield stress!----------------------------------------------------------------------- CALL VOCE(SigM, SIGMA0, EpsPEqvTr, T_1, Q_1, T_2, . Q_2, T_3, Q_3)! !-----------------------------------------------------------------------! Calculate the effective void volume fraction f^*!----------------------------------------------------------------------- IF ((VoidTr .ge. f_c).and.(f_f.ne.0)) THEN C_Factor = ((1.0/BETA1) - f_c) / (f_f - f_c) Eff_Void = f_c + C_Factor * (VoidTr - f_c) ELSE Eff_Void = VoidTr END IF !!---------- In order to avoid numerical errors ------------------------- IF (Eff_Void .ge. 0.99*(1.0/BETA1)) THEN Eff_Void = 0.99*(1.0/BETA1) END IF ! !-----------------------------------------------------------------------! Compute Bracket Term!----------------------------------------------------------------------- TH=(-3.0*BETA2*pTr)/(2.0*SigM)! !-----------------------------------------------------------------------! Compute Yield Function!----------------------------------------------------------------------- F = (qTr**2.0/SigM**2.0)+(2.0*Eff_Void*BETA1*cosh(TH)) . -1.0-BETA3*Eff_Void**2.0!!-----------------------------------------------------------------------! Check for plasticity!----------------------------------------------------------------------- IF(F.GT.0.0)THEN ! Plastic flow!-----------------------------------------------------------------------! Initialize Variables!----------------------------------------------------------------------- p = pTr q = qTr Void = VoidTr EpsPEqv = EpsPEqvTr DeltaEpsP = 0. DeltaEpsQ = 0. omega_shear = 1.0 - (27*J_3 / (2.0*qTr**3))**2 IF(q.LE.1.0E-20)THEN q=1.0E-20 ENDIF!!-----------------------------------------------------------------------! Start iteration!----------------------------------------------------------------------- DO ITER=1,MXITER IF(q.LE.1.0E-20)THEN q=1.0E-20 ENDIF!-----------------------------------------------------------------------! Calculate Derivatives!----------------------------------------------------------------------- SigMDerivX2 = T_1*Q_1*EXP(-Q_1*EpsPEqv) ! dSigM / deps_eq_p . + T_2*Q_2*EXP(-Q_2*EpsPEqv) . + T_3*Q_3*EXP(-Q_3*EpsPEqv) GDerivQ = 2.0*q/SigM**2.0 ! dG/dq GDerivP = -3.0*Eff_Void*BETA1*BETA2/SigM*SINH(TH) ! dG/dp GDerivX1 = 2.0*BETA1*COSH(TH)-2.0*BETA3*Eff_Void ! dG/df^* GDerivX2 = -2.0*q**2.0 * SigMDerivX2 /SigM**3.0 ! dG/deps_eq_p . + 3.0*Eff_Void*BETA1*BETA2*p*SigMDerivX2 . / SigM**2.0*SINH(TH) G2DerivPQ = 0.0 G2DerivP2 = (9.0*Eff_Void*BETA1*BETA2**2.0)/(2.0*SigM**2.0) . *COSH(TH) G2DerivQ2 = 2.0/SigM**2.0 G2DerivQX1 = 0.0 G2DerivQX2 = -4.0*q*SigMDerivX2/SigM**3.0 G2DerivPX1 = -3.0*BETA1*BETA2/SigM*SINH(TH) G2DerivPX2 = (3.0*SigMDerivX2*Eff_Void*BETA1*BETA2/SigM**2.0) . *(SINH(TH)+TH*COSH(TH)) h1DerivDeltaEpsP = 1.0-Void h1DerivP = 0.0 h2DerivDeltaEpsP = (-p)/((1-Void)*SigM) h2DerivP = -DeltaEpsP/((1-Void)*SigM) h1DerivDeltaEpsQ = k_s*Void*omega_shear ! 0 Initially h1DerivQ = 0.0 h2DerivDeltaEpsQ = q/((1-Void)*SigM) h2DerivQ = (DeltaEpsQ) /((1-Void)*SigM) h1DerivX1 = -DeltaEpsP IF((p.LE.0.0).and.(FN.ne.0))THEN h1DerivX2 = -DeltaEpsPEqv*FN/(SN*sqrt(2.0*PI)) . *(EpsPEqv-EPSN)/(SN**2.0)*exp(-0.5*((EpsPEqv-EPSN)/SN)**2.0) ELSE h1DerivX2 = 0.0 ENDIF h2DerivX1 = (-p*DeltaEpsP+q*DeltaEpsQ) . /(SigM*(1.0-Void)**2.0) h2DerivX2 = SigMDerivX2*(p*DeltaEpsP-q*DeltaEpsQ) . /((1.0-Void)*SigM**2.0) C11Inv = 1.0-h1DerivX1 C12Inv = h1DerivX2 C21Inv = h2DerivX1 C22Inv = 1.0-h2DerivX2 DetC = C11Inv*C22Inv-C12Inv*C21Inv C11 = C22Inv/DetC C12 = -C12Inv/DetC C21 = -C21Inv/DetC C22 = C11Inv/DetC X1DerivDeltaEpsP = . +C11*((h1DerivDeltaEpsP)+K*(h1DerivP)) . +C12*((h2DerivDeltaEpsP)+K*(h2DerivP)) X2DerivDeltaEpsP = . +C22*((h2DerivDeltaEpsP)+K*(h2DerivP)) . +C21*((h1DerivDeltaEpsP)+K*(h1DerivP)) X1DerivDeltaEpsQ = . +C11*((h1DerivDeltaEpsQ)-3.0*G*(h1DerivQ)) . +C12*((h2DerivDeltaEpsQ)-3.0*G*(h2DerivQ)) X2DerivDeltaEpsQ = . +C22*((h2DerivDeltaEpsQ)-3.0*G*(h2DerivQ)) . +C21*((h1DerivDeltaEpsQ)-3.0*G*(h1DerivQ))!-----------------------------------------------------------------------! Calculate Constants For Update!----------------------------------------------------------------------- A11 = GDerivQ+DeltaEpsP*(K*G2DerivPQ . +G2DerivQX1*X1DerivDeltaEpsP . +G2DerivQX2*X2DerivDeltaEpsP) . +DeltaEpsQ*(K*G2DerivP2 . +G2DerivPX1*X1DerivDeltaEpsP . +G2DerivPX2*X2DerivDeltaEpsP) A12 = GDerivP+DeltaEpsP*(-3.0*G*G2DerivQ2 . +G2DerivQX1*X1DerivDeltaEpsQ . +G2DerivQX2*X2DerivDeltaEpsQ) . +DeltaEpsQ*(-3.0*G*G2DerivPQ . +G2DerivPX1*X1DerivDeltaEpsQ . +G2DerivPX2*X2DerivDeltaEpsQ) A21 = K*GDerivP+GDerivX1*X1DerivDeltaEpsP . +GDerivX2*X2DerivDeltaEpsP A22 = -3.0*G*GDerivQ+GDerivX1*X1DerivDeltaEpsQ . +GDerivX2*X2DerivDeltaEpsQ b1 = -DeltaEpsP*GDerivQ-DeltaEpsQ*GDerivP b2 = -F T1 = b2-A21*b1/A11 T2 = A22-A21*A12/A11 Cq = T1/T2 Cp = (b1-A12*Cq)/A11!!-----------------------------------------------------------------------! Update Values!----------------------------------------------------------------------- DeltaEpsP = DeltaEpsP+Cp DeltaEpsQ = DeltaEpsQ+Cq!!---------- Update pressure and eqv stress ---------------------------------- p=pTr+K*DeltaEpsP q=qTr-3.0*G*DeltaEpsQ!!---------- Update equivalent plastic strain ---------------------------------- DeltaEpsPEqv = (-p*DeltaEpsP+q*DeltaEpsQ)/((1.0-Void)*SigM)! !---------- Update Void Growth term ---------------------------------- DeltaVoidG = DeltaEpsP*(1.0-Void)! !---------- Update Void Nucleation term ---------------------------------- IF ((p.LE.0.0).and.(FN.ne.0)) THEN DeltaVoidN = DeltaEpsPEqv*FN/(SN*sqrt(2.0*PI)) . * exp(-0.5*((EpsPEqv-EPSN)/SN)**2) ELSE DeltaVoidN = 0.0 END IF!!---------- Update Void Shear term ---------------------------------- IF (k_s .ne. 0) THEN DeltaVoidS = k_s*VoidTr*omega_shear*DeltaEpsQ ELSE DeltaVoidS = 0. END IF !!---------- Update the equivalent plastic strain --------------------------------------- EpsPEqv = EpsPEqvTr + DeltaEpsPEqv! !---------- Update yield strength ------------------------------------------------ CALL VOCE(SigM, SIGMA0, EpsPEqv, T_1, Q_1, T_2, Q_2, T_3, Q_3)!---------- Update bracket term ------------------------------------------------ TH = (-3.0*BETA2*p)/(2.0*SigM)!!----------- Sum all void contributions ------------------------------------------------ DeltaVoid = DeltaVoidG + DeltaVoidN + DeltaVoidS Void = VoidTr+DeltaVoid!! ! VOID COALESCENCE! !----------- Calculate the effective void volume fraction f^* -------------------------- IF ((Void .ge. f_c).and.(f_f.ne.0)) THEN C_Factor = ((1/BETA1) - f_c) / (f_f - f_c) Eff_Void = f_c + C_Factor * (Void - f_c) ELSE Eff_Void = Void END IF !!---------- In order to avoid numerical errors ------------------------- IF (Eff_Void .ge. 0.99*(1/BETA1)) THEN Eff_Void = 0.99*(1/BETA1) END IF ! !---------- Check the modified yield function -------------------------- F = (q**2.0/SigM**2.0)+ . (2.0*Eff_Void*BETA1*COSH(TH)) . -1.0-BETA3*Eff_Void**2.0 ! !-----------------------------------------------------------------------! Compute convergence criterion!----------------------------------------------------------------------- RESNOR = ABS(F/SigM)!-----------------------------------------------------------------------! Check for convergence!----------------------------------------------------------------------- IF(RESNOR.LE.TOL)THEN ! RMAP has converged!-----------------------------------------------------------------------! Now, introducing the \kappa parameter in order to speed up the void deformations.! This code is only executed if there are plastic deformations in the element, which is correct.! The adjusted update in the void volume fraction is local, but the rate is adjusted dependent on its ! neighbors plastics increments...!----------------------------------------------------------------------- ! Based on the plastic increment ratio NL/L, we will speed up / slow down the void growth:! kappa = 1. !If \kappa=1 is used, the local solution is obtained Void = VoidTr + kappa * DeltaVoid! !-------------- In order to avoid numerical errors ----------------------- IF (f_f.ne.0) THEN ! Then we have coalescence, and Void should be less than f_f IF (Void .ge. 0.99*f_f) THEN Void = 0.99*f_f END IF ELSE IF (Void .ge. 0.99*(1.0/BETA1)) THEN Void = 0.99*(1.0/BETA1) END IF END IF!!-----------------------------------------------------------------------! Update the stress tensor!----------------------------------------------------------------------- STRESSNEW(km,1) = (stresstr(1)+pTr)*(q/qTr) - p STRESSNEW(km,2) = (stresstr(2)+pTr)*(q/qTr) - p STRESSNEW(km,3) = (stresstr(3)+pTr)*(q/qTr) - p STRESSNEW(km,4) = stresstr(4)*(q/qTr) IF ( nshr .gt. 1 ) THEN STRESSNEW(km,5) = stresstr(5)*(q/qTr) STRESSNEW(km,6) = stresstr(6)*(q/qTr) END IF!-----------------------------------------------------------------------! Update the history variables!----------------------------------------------------------------------- IF(abs(Void).LE.1.0E-10)THEN Void=0.0 ENDIF STATENEW(km,1) = Void ! VoidFrac STATENEW(km,2) = EpsPEqv ! EpsPEqv STATENEW(km,3) = q ! SigEqv STATENEW(km,4) = p ! SigH STATENEW(km,5) = SigM ! SigM STATENEW(km,6) = F ! F STATENEW(km,8) = EpsPEqv ! Local plastic strain!---------------------------------------------------------------------------! Delete the element if the void content has reached the critical value: IF (STATEOLD(km,7) .eq. 1) THEN ! If the element hasn't been deleted yet IF (f_f .ne. 0) THEN IF (Void .ge. 0.98*f_f) THEN STATENEW(km,7) = 0 ELSE STATENEW(km,7) = 1 END IF ELSE IF (Void.ge.(0.98*(1/BETA1))) THEN STATENEW(km,7) = 0 ! If critical void content: delete the element ELSE STATENEW(km,7) = 1 END IF END IF ELSE STATENEW(km,7) = 0 ! The element was deleted in an earlier increment.. END IF !--------------------------------------------------------------------------- GOTO 90 ELSE ! RMAP has not converged yet IF(ITER.eq.(MXITER-1))THEN PRINT*, "RMAP has not converged!" PRINT*, "Convergence: ", RESNOR PRINT*, Void STOP ENDIF ENDIF ENDDO !The plastic convergence loop in the vumat ELSE!!-----------------------------------------------------------------------! If this case, the element has undergone purely elastic deformation: !-----------------------------------------------------------------------! Update the stress tensor!----------------------------------------------------------------------- STRESSNEW(km,1) = stresstr(1) STRESSNEW(km,2) = stresstr(2) STRESSNEW(km,3) = stresstr(3) STRESSNEW(km,4) = stresstr(4) IF ( nshr .gt. 1 ) THEN STRESSNEW(km,5) = stresstr(5) STRESSNEW(km,6) = stresstr(6) END IF!-----------------------------------------------------------------------! Update the history variables!----------------------------------------------------------------------- STATENEW(km,1) = VoidTr ! VoidFrac STATENEW(km,2) = EpsPEqvTr ! EpsPEqv STATENEW(km,3) = SigEqv ! SigEqv STATENEW(km,4) = pTr ! SigH STATENEW(km,5) = SigM ! SigM STATENEW(km,6) = F ! F STATENEW(km,7) = 1 ! Delete the element or not.. STATENEW(km,9) = EpsPEqvTr ! EpsPEqv ENDIF 90 CONTINUE!-----------------------------------------------------------------------! Store the specific internal energy -!----------------------------------------------------------------------- if ( nshr .eq. 1 ) then stressPower = half * ( . ( stressOld(km,1) + stressNew(km,1) ) * strainInc(km,1) + . ( stressOld(km,2) + stressNew(km,2) ) * strainInc(km,2) + . ( stressOld(km,3) + stressNew(km,3) ) * strainInc(km,3) ) + . ( stressOld(km,4) + stressNew(km,4) ) * strainInc(km,4) else stressPower = half * ( . ( stressOld(km,1) + stressNew(km,1) ) * strainInc(km,1) + . ( stressOld(km,2) + stressNew(km,2) ) * strainInc(km,2) + . ( stressOld(km,3) + stressNew(km,3) ) * strainInc(km,3) ) + . ( stressOld(km,4) + stressNew(km,4) ) * strainInc(km,4) + . ( stressOld(km,5) + stressNew(km,5) ) * strainInc(km,5) + . ( stressOld(km,6) + stressNew(km,6) ) * strainInc(km,6) end if enerInternNew(km) = enerInternOld(km) + stressPower / density(km) ! !-----------------------------------------------------------------------! END the material point loop and vumat branch:!----------------------------------------------------------------------- END DO ! End of the material point loop END IF ! End of the preprocessor / vumat branch!!-----------------------------------------------------------------------! END SUBROUTINE!----------------------------------------------------------------------- RETURN END! ! !!! 3 term Voce hardening rule: SUBROUTINE VOCE(sigma_new, sigma0, peeqOld, T_1, Q_1, T_2, .Q_2, T_3, Q_3) include 'vaba_param.inc'! term1 = T_1 * (1 - EXP(-Q_1 * peeqOld)) term2 = T_2 * (1 - EXP(-Q_2 * peeqOld)) term3 = T_3 * (1 - EXP(-Q_3 * peeqOld)) sigma_new = sigma0 + term1 + term2 + term3 RETURN END!!材料卡片来源：我的博士日记
