GTN 模型参考vumat
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
!
!
