VUMAT for Hyperelastic Materials - Update to Bergstrom's Example

Hi Travis, that I great. Are you able to also share the source code here?

-Jorgen

Hello,

I have confirmed that I can post the code. For a detailed explanation of the statements, refer to the article above which summarizes the finite elasticity theory that is necessary to build a hyperelastic VUMAT. A discussion of implementation of the theory in Abaqus is given in Appendix C. The article was built largely from the more detailed accounts in the books by Holzapfel (2000) and Bergström (2015).

Some pertinent updates to the Neo-Hookean subroutine that is commonly shared online are:

- updated names for Abaqus variables.
- a gen-Yeoh SEF is used. Set:
----- (K2,K3) = (0,0) for Neo-Hookean SEF.
----- (m,p,q) = (1,2,3) for Yeoh SEF.
----- (K2,K3) = (0,2) for Davies-De-Thomas SEF (1994, Rubber Chem. & Tech. 67 (4), 716).
- includes linear elastic initialization step required per Abaqus VUMAT subroutine template.
- includes logic to handle indeterminate stress condition that can arise with the gen-Yeoh SEF.

The code only handles 3D elements. A 2D plane-strain or axisymmetric (but not plane-stress) version can be written by removing references to the 5th and 6th elements in the stress, strain, and strain increment vectors. Alternatively, the logic in Bergström's Neo-Hookean example can be included to handle 2D and 3D cases.

Regards,
Travis

C *******************************************************************
C  VUMAT (3D only) for gen-Yeoh Strain-Energy Function
C     Compatibility confirmed for Abaqus v6.14 & R2018
C  ---------------------------------------------------
C  Authors:  Travis Hohenberger & Richard Windslow
C  Date:     2019-07-11
C  E-mail:   twhohen@gmail.com
C  Source:   Hohenberger et al (2019). Rubber Chem. & Tech., 92 (4), 653
C *******************************************************************
C
C  Strain-energy function:
C
C  W = K1*(I1-3)^m + K2*(I1-3)^p + K3*(I1-3)^q + (1/D1)*(J-1)^2
C
C  where K1, K2, K3, m, p, q are distortional fitting parameters and D1
C  is a volumetric fitting parameter. I1 is the first invariant of the
C  modified stretch tensor. J is the volumetric ratio.
C
C **********************************************************************
C
      SUBROUTINE VUMAT(
     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,
     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),
     5          fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
     6          stateOld(nblock,nstatev), enerInternOld(nblock),
     7          enerInelasOld(nblock), tempNew(nblock),
     8          stretchNew(nblock,ndir+nshr),
     9          defgradNew(nblock,ndir+nshr),
     1          fieldNew(nblock,nfieldv),
     2          stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
     3          enerInternNew(nblock), enerInelasNew(nblock)
C
      CHARACTER*80 cmname
C
C     PARAMETERS
C     ----------
      REAL*8    oneThrd, half, twoThrd, one, two, three, thresh
      PARAMETER(oneThrd=1.d0/3.d0, half=0.5d0, twoThrd=2.d0/3.d0,
     $          one=1.d0, two=2.d0, three=3.d0, thresh=10.d0**-12.d0)
C
C     LOCAL VARIABLES
C     ---------------
      REAL*8 g0      , k0     , twoG   , lmda   , trace  , d1      ,
     $       k1      , k2     , k3     , em     , pe     , qu      ,
     $       Bxx     , Byy    , Bzz    , Bxy    , Bxz    , Byz     ,
     $       BbarXX  , BbarYY , BbarZZ , BbarXY , BbarXZ , BbarYZ  ,
     $       dBbarXX , dBbarYY, dBbarZZ, dBbarXY, dBbarXZ, dBbarYZ ,
     $       J       , J23    , duDi1  , duDi3  , I1     , g1      ,
     $       p0      , j1     , u1
C
      k1 = props(1)
      k2 = props(2)
      k3 = props(3)
      em = props(4)
      pe = props(5)
      qu = props(6)
      d1 = props(7)
C
      g0 = two * k1
      k0 = two / d1
C
C     ***************************************************************
C     ------------ INITIALIZE MATERIAL AS LINEARLY ELASTIC -------------
C     ***************************************************************
C
      twoG = two * g0
      lmda = k0 - twoG * oneThrd
C
      IF (totalTime.EQ.0.0) THEN
C
         DO k = 1,nblock
            trace = strainInc(k,1) + strainInc(k,2) + strainInc(k,3)
            stressNew(k,1) = stressOld(k,1) + twoG*strainInc(k,1) +
     $                       lmda*trace
            stressNew(k,2) = stressOld(k,2) + twoG*strainInc(k,2) +
     $                       lmda*trace
            stressNew(k,3) = stressOld(k,3) + twoG*strainInc(k,3) +
     $                       lmda*trace
            stressNew(k,4) = stressOld(k,4) + twoG*strainInc(k,4)
            stressNew(k,5) = stressOld(k,5) + twoG*strainInc(k,5)
            stressNew(k,6) = stressOld(k,6) + twoG*strainInc(k,6)
         END DO
C
         RETURN
C
      END IF
C
C     ***************************************************************
C     ----------- START LOOP FOR MATERIAL POINT CALCULATIONS -----------
C     ***************************************************************
C
      DO k = 1,nblock
C
C        CALCULATE LEFT CAUCHY-GREEN STRAIN TENSOR, B = U*U
C        --------------------------------------------------
         Bxx = stretchNew(k,1) * stretchNew(k,1) +
     $         stretchNew(k,4) * stretchNew(k,4) +
     $         stretchNew(k,6) * stretchNew(k,6)
         Byy = stretchNew(k,2) * stretchNew(k,2) +
     $         stretchNew(k,4) * stretchNew(k,4) +
     $         stretchNew(k,5) * stretchNew(k,5)
         Bzz = stretchNew(k,3) * stretchNew(k,3) +
     $         stretchNew(k,5) * stretchNew(k,5) +
     $         stretchNew(k,6) * stretchNew(k,6)
         Bxy = stretchNew(k,1) * stretchNew(k,4) +
     $         stretchNew(k,4) * stretchNew(k,2) +
     $         stretchNew(k,6) * stretchNew(k,5)
         Bxz = stretchNew(k,1) * stretchNew(k,6) +
     $         stretchNew(k,4) * stretchNew(k,5) +
     $         stretchNew(k,6) * stretchNew(k,3)
         Byz = stretchNew(k,4) * stretchNew(k,6) +
     $         stretchNew(k,2) * stretchNew(k,5) +
     $         stretchNew(k,5) * stretchNew(k,3)
C
C        CALCULATE J = |F| = |U|
C        -----------------------
C
         J =    stretchNew(k,1) *
     $        ( stretchNew(k,2) * stretchNew(k,3)   -
     $          stretchNew(k,5) * stretchNew(k,5) ) +
     $          stretchNew(k,4) *
     $        ( stretchNew(k,5) * stretchNew(k,6)   -
     $          stretchNew(k,3) * stretchNew(k,4) ) +
     $          stretchNew(k,6) *
     $        ( stretchNew(k,4) * stretchNew(k,5)   -
     $          stretchNew(k,2) * stretchNew(k,6) )
C
C        CALCULATE MODIFIED STRAIN TENSOR, Bbar = J^(-2/3)*B 
C        ---------------------------------------------------
         J23 = J**(-twoThrd)
C
         BbarXX = J23 * Bxx
         BbarYY = J23 * Byy
         BbarZZ = J23 * Bzz
         BbarXY = J23 * Bxy
         BbarXZ = J23 * Bxz
         BbarYZ = J23 * Byz
C
C        FIRST INVARIANT OF Bbar = tr(Bbar)
C        ----------------------------------
         I1 = BbarXX + BbarYY + BbarZZ
C
C        DEVIATORIC PART OF Bbar
C        -----------------------
         p0 = oneThrd * I1
C
         dBbarXX = BbarXX - p0
         dBbarYY = BbarYY - p0
         dBbarZZ = BbarZZ - p0
         dBbarXY = BbarXY
         dBbarXZ = BbarXZ
         dBbarYZ = BbarYZ
C
C        DERIVATIVES OF STRAIN-ENERGY FUNCTION
C        -------------------------------------
         j1 = I1 - three
C
         IF (j1.LT.thresh) THEN
            duDi1 = zero
         ELSE
            duDi1 = em * k1 * j1**(em-one) +
     $              pe * k2 * j1**(pe-one) +
     $              qu * k3 * j1**(qu-one)
         END IF
C
         duDi3 = two/d1 * (J - one)
C
C        COROTATIONAL CAUCHY (TRUE) STRESSES
C        -----------------------------------
         g1 = two/J * duDi1
C
         stressNew(k,1) = g1 * dBbarXX + duDi3
         stressNew(k,2) = g1 * dBbarYY + duDi3
         stressNew(k,3) = g1 * dBbarZZ + duDi3
         stressNew(k,4) = g1 * dBbarXY
         stressNew(k,5) = g1 * dBbarYZ
         stressNew(k,6) = g1 * dBbarXZ
C
C        UPDATE SPECIFIC INTERNAL ENERGY
C        -------------------------------
         u1 = 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) +
     $                  two * ( (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) ) )
C
         enerInternNew(k) = enerInternOld(k) + u1 / density(k)
C
      END DO
C
      RETURN
C
      END SUBROUTINE VUMAT

So it seems I cannot edit my post, so I apologize for the double post. I wanted to add the following general question:

In my VUMATs, I use matrix multiplications rather than computing the components individually, for example:

C	F^eE = U
        U(1,1) = stretchNew(k,1)
        U(2,2) = stretchNew(k,2)
        U(3,3) = stretchNew(k,3)
        U(1,2) = stretchNew(k,4)
        U(2,3) = stretchNew(k,5)
        U(1,3) = stretchNew(k,6)
        U(2,1) = U(1,2)
        U(3,1) = U(1,3)
        U(3,2) = U(2,3)
C
C       B^ee = F^eE * F^eE^{T} = U^2
        B = zero
        B = matmul(U,U)

Is there any reason to choose one over the other?

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.inc'
c
c 3D Orthotropic Elasticity with Hashin 3d Failure criterion
c
c The state variables are stored as:
c state(*,1) = material point status
c state(*,2:7) = damping stresses
c
c User defined material properties are stored as
c * First line:
c props(1) --> Young's modulus in 1-direction, E1
c props(2) --> Young's modulus in 2-direction, E2
c props(3) --> Young's modulus in 3-direction, E3
c props(4) --> Poisson's ratio, nu12
c props(5) --> Poisson's ratio, nu13
c props(6) --> Poisson's ratio, nu23
c props(7) --> Shear modulus, G12
c props(8) --> Shear modulus, G13
c
c * Second line:
c props(9) --> Shear modulus, G23
c props(10) --> beta damping parameter
c props(11) --> "not used"
c props(12) --> "not used"
c props(13) --> "not used"
c props(14) --> "not used"
c props(15) --> "not used"
c props(16) --> "not used"
c
c * Third line:
c props(17) --> Ultimate tens stress in 1-direction, sigu1t
c props(18) --> Ultimate comp stress in 1-direction, sigu1c
c props(19) --> Ultimate tens stress in 2-direction, sigu2t
c props(20) --> Ultimate comp stress in 2-direction, sigu2c
c props(21) --> Ultimate tens stress in 2-direction, sigu3t
c props(22) --> Ultimate comp stress in 2-direction, sigu3c
c props(23) --> "not used"
c props(24) --> "not used"
c
c * Fourth line:
c props(25) --> Ultimate shear stress, sigu12
c props(26) --> Ultimate shear stress, sigu13
c props(27) --> Ultimate shear stress, sigu23
c props(28) --> "not used"
c props(29) --> "not used"
c props(30) --> "not used"
c props(31) --> "not used"
c props(32) --> "not used"
c

dimension props(nprops), density(nblock),
1 coordMp(nblock,*),
2 charLength(*), strainInc(nblock,ndir+nshr),
3 relSpinInc(nblock,nshr), tempOld(nblock),
4 stretchOld(nblock,ndir+nshr), defgradOld(nblock,ndir+nshr+nshr),
5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
6 stateOld(nblock,nstatev), enerInternOld(nblock),
7 enerInelasOld(nblock), tempNew(*),
8 stretchNew(nblock,ndir+nshr), defgradNew(nblock,ndir+nshr+nshr),
9 fieldNew(nblock,nfieldv), stressNew(nblock,ndir+nshr),
1 stateNew(nblock,nstatev),
2 enerInternNew(nblock), enerInelasNew(nblock)
*
character*80 cmname
*
parameter( zero = 0.d0, one = 1.d0, two = 2.d0, half = .5d0 )
*
parameter(
* i_svd_DmgFiberT = 1,
* i_svd_DmgFiberC = 2,
* i_svd_DmgMatrixT = 3,
* i_svd_DmgMatrixC = 4,
* i_svd_statusMp = 5,
* i_svd_dampStress = 6,
c * i_svd_dampStressXx = 6,
c * i_svd_dampStressYy = 7,
c * i_svd_dampStressZz = 8,
c * i_svd_dampStressXy = 9,
c * i_svd_dampStressYz = 10,
c * i_svd_dampStressZx = 11,
* i_svd_Strain = 12,
c * i_svd_StrainXx = 12,
c * i_svd_StrainYy = 13,
c * i_svd_StrainZz = 14,
c * i_svd_StrainXy = 15,
c * i_svd_StrainYz = 16,
c * i_svd_StrainZx = 17,
* n_svd_required = 17 )
*
parameter(
* i_s33_Xx = 1,
* i_s33_Yy = 2,
* i_s33_Zz = 3,
* i_s33_Xy = 4,
* i_s33_Yz = 5,
* i_s33_Zx = 6 )
*
* Structure of property array
parameter (
* i_pro_E1 = 1,
* i_pro_E2 = 2,
* i_pro_E3 = 3,
* i_pro_nu12 = 4,
* i_pro_nu13 = 5,
* i_pro_nu23 = 6,
* i_pro_G12 = 7,
* i_pro_G13 = 8,
* i_pro_G23 = 9,
*
* i_pro_beta = 10,
*
* i_pro_sigu1t = 17,
* i_pro_sigu1c = 18,
* i_pro_sigu2t = 19,
* i_pro_sigu2c = 20,
* i_pro_sigu3t = 21,
* i_pro_sigu3c = 22,
* i_pro_sigu12 = 25,
* i_pro_sigu13 = 26,
* i_pro_sigu23 = 27 )
* Temporary arrays
dimension eigen(maxblk*3)
*
* Read material properties
*
E1 = props(i_pro_E1)
E2 = props(i_pro_E2)
E3 = props(i_pro_E3)
xnu12 = props(i_pro_nu12)
xnu13 = props(i_pro_nu13)
xnu23 = props(i_pro_nu23)
G12 = props(i_pro_G12)
G13 = props(i_pro_G13)
G23 = props(i_pro_G23)
*
xnu21 = xnu12 * E2 / E1
xnu31 = xnu13 * E3 / E1
xnu32 = xnu23 * E3 / E2
*
*
* Compute terms of stiffness matrix
gg = one / ( one - xnu12*xnu21 - xnu23*xnu32 - xnu31*xnu13
* - two*xnu21*xnu32*xnu13 )
C11 = E1 * ( one - xnu23*xnu32 ) * gg
C22 = E2 * ( one - xnu13*xnu31 ) * gg
C33 = E3 * ( one - xnu12*xnu21 ) * gg
C12 = E1 * ( xnu21 + xnu31*xnu23 ) * gg
C13 = E1 * ( xnu31 + xnu21*xnu32 ) * gg
C23 = E2 * ( xnu32 + xnu12*xnu31 ) * gg
*
f1t = props(i_pro_sigu1t)
f1c = props(i_pro_sigu1c)
f2t = props(i_pro_sigu2t)
f2c = props(i_pro_sigu2c)
f3t = props(i_pro_sigu3t)
f3c = props(i_pro_sigu3c)
f12 = props(i_pro_sigu12)
f13 = props(i_pro_sigu13)
f23 = props(i_pro_sigu23)
*
beta = props(i_pro_beta)
*
* Assume purely elastic material at the beginning of the analysis
*
if ( totalTime .eq. zero ) then
if (nstatev .lt. n_svd_Required) then
call xplb_abqerr(-2,'Subroutine VUMAT requires the '//
* 'specification of %I state variables. Check the '//
* 'definition of *DEPVAR in the input file.',
* n_svd_Required,zero,' ')
call xplb_exit
end if
call OrthoEla3dExp ( nblock,
* stateOld(1,i_svd_DmgFiberT),
* stateOld(1,i_svd_DmgFiberC),
* stateOld(1,i_svd_DmgMatrixT),
* stateOld(1,i_svd_DmgMatrixC),
* C11, C22, C33, C12, C23, C13, G12, G23, G13,
* strainInc,
* stressNew )
return
end if
*
* Update total elastic strain
call strainUpdate ( nblock, strainInc,
* stateOld(1,i_svd_strain), stateNew(1,i_svd_strain) )
*
* Stress update
call OrthoEla3dExp ( nblock,
* stateOld(1,i_svd_DmgFiberT),
* stateOld(1,i_svd_DmgFiberC),
* stateOld(1,i_svd_DmgMatrixT),
* stateOld(1,i_svd_DmgMatrixC),
* C11, C22, C33, C12, C23, C13, G12, G23, G13,
* stateNew(1,i_svd_strain),
* stressNew )
*
* Failure evaluation
*
call copyr ( nblock,
* stateOld(1,i_svd_DmgFiberT), stateNew(1,i_svd_DmgFiberT) )
call copyr ( nblock,
* stateOld(1,i_svd_DmgFiberC), stateNew(1,i_svd_DmgFiberC) )
call copyr ( nblock,
* stateOld(1,i_svd_DmgMatrixT), stateNew(1,i_svd_DmgMatrixT) )
call copyr ( nblock,
* stateOld(1,i_svd_DmgMatrixC), stateNew(1,i_svd_DmgMatrixC) )
nDmg = 0
call eig33Anal ( nblock, stretchNew, eigen )
call Hashin3d ( nblock, nDmg,
* f1t, f2t, f3t, f1c, f2c, f3c, f12, f23, f13,
* stateNew(1,i_svd_DmgFiberT),
* stateNew(1,i_svd_DmgFiberC),
* stateNew(1,i_svd_DmgMatrixT),
* stateNew(1,i_svd_DmgMatrixC),
* stateNew(1,i_svd_statusMp),
* stressNew, eigen )
* -- Recompute stresses if new Damage is occurring
if ( nDmg .gt. 0 ) then
call OrthoEla3dExp ( nblock,
* stateNew(1,i_svd_DmgFiberT),
* stateNew(1,i_svd_DmgFiberC),
* stateNew(1,i_svd_DmgMatrixT),
* stateNew(1,i_svd_DmgMatrixC),
* C11, C22, C33, C12, C23, C13, G12, G23, G13,
* stateNew(1,i_svd_strain),
* stressNew )
end if
*
* Beta damping
if ( beta .gt. zero ) then
call betaDamping3d ( nblock,
* beta, dt, strainInc,
* stressOld, stressNew,
* stateNew(1,i_svd_statusMp),
* stateOld(1,i_svd_dampStress),
* stateNew(1,i_svd_dampStress) )
end if
*
* Integrate the internal specific energy (per unit mass)
*
call EnergyInternal3d ( nblock, stressOld, stressNew,
* strainInc, density, enerInternOld, enerInternNew )
*
return
end

************************************************************
* OrthoEla3dExp: Orthotropic elasticity - 3d *
************************************************************
subroutine OrthoEla3dExp ( nblock,
* dmgFiberT, dmgFiberC, dmgMatrixT, dmgMatrixC,
* C11, C22, C33, C12, C23, C13, G12, G23, G13,
* strain, stress )
*
include 'vaba_param.inc'

* Orthotropic elasticity, 3D case -
*
parameter( zero = 0.d0, one = 1.d0, two = 2.d0)
parameter(
* i_s33_Xx = 1,
* i_s33_Yy = 2,
* i_s33_Zz = 3,
* i_s33_Xy = 4,
* i_s33_Yz = 5,
* i_s33_Zx = 6,
* n_s33_Car = 6 )
*
dimension strain(nblock,n_s33_Car),
* dmgFiberT(nblock), dmgFiberC(nblock),
* dmgMatrixT(nblock), dmgMatrixC(nblock),
* stress(nblock,n_s33_Car)
* -- shear fraction in matrix tension and compression mode
parameter ( smt = 0.9d0, smc = 0.5d0 )
*
do k = 1, nblock
* -- Compute damaged stiffness
dft = dmgFiberT(k)
dfc = dmgFiberC(k)
dmt = dmgMatrixT(k)
dmc = dmgMatrixC(k)
df = one - ( one - dft ) * ( one - dfc )
*
dC11 = ( one - df ) * C11
dC22 = ( one - df ) * ( one - dmt ) * ( one - dmc ) * C22
dC33 = ( one - df ) * ( one - dmt ) * ( one - dmc ) * C33
dC12 = ( one - df ) * ( one - dmt ) * ( one - dmc ) * C12
dC23 = ( one - df ) * ( one - dmt ) * ( one - dmc ) * C23
dC13 = ( one - df ) * ( one - dmt ) * ( one - dmc ) * C13
dG12 = ( one - df )
* * ( one - smt*dmt ) * ( one - smc*dmc ) * G12
dG23 = ( one - df )
* * ( one - smt*dmt ) * ( one - smc*dmc ) * G23
dG13 = ( one - df )
* * ( one - smt*dmt ) * ( one - smc*dmc ) * G13
* -- Stress update
stress(k,i_s33_Xx) = dC11 * strain(k,i_s33_Xx)
* + dC12 * strain(k,i_s33_Yy)
* + dC13 * strain(k,i_s33_Zz)
stress(k,i_s33_Yy) = dC12 * strain(k,i_s33_Xx)
* + dC22 * strain(k,i_s33_Yy)
* + dC23 * strain(k,i_s33_Zz)
stress(k,i_s33_Zz) = dC13 * strain(k,i_s33_Xx)
* + dC23 * strain(k,i_s33_Yy)
* + dC33 * strain(k,i_s33_Zz)
stress(k,i_s33_Xy) = two * dG12 * strain(k,i_s33_Xy)
stress(k,i_s33_Yz) = two * dG23 * strain(k,i_s33_Yz)
stress(k,i_s33_Zx) = two * dG13 * strain(k,i_s33_Zx)
end do
*
return
end

************************************************************
* strainUpdate: Update total strain *
************************************************************
subroutine strainUpdate ( nblock,
* strainInc, strainOld, strainNew )
*
include 'vaba_param.inc'
*
parameter(
* i_s33_Xx = 1,
* i_s33_Yy = 2,
* i_s33_Zz = 3,
* i_s33_Xy = 4,
* i_s33_Yz = 5,
* i_s33_Zx = 6,
* n_s33_Car = 6 )
*
dimension strainInc(nblock,n_s33_Car),
* strainOld(nblock,n_s33_Car),
* strainNew(nblock,n_s33_Car)
*
do k = 1, nblock
strainNew(k,i_s33_Xx)= strainOld(k,i_s33_Xx)
* + strainInc(k,i_s33_Xx)
strainNew(k,i_s33_Yy)= strainOld(k,i_s33_Yy)
* + strainInc(k,i_s33_Yy)
strainNew(k,i_s33_Zz)= strainOld(k,i_s33_Zz)
* + strainInc(k,i_s33_Zz)
strainNew(k,i_s33_Xy)= strainOld(k,i_s33_Xy)
* + strainInc(k,i_s33_Xy)
strainNew(k,i_s33_Yz)= strainOld(k,i_s33_Yz)
* + strainInc(k,i_s33_Yz)
strainNew(k,i_s33_Zx)= strainOld(k,i_s33_Zx)
* + strainInc(k,i_s33_Zx)
end do
*
return
end

************************************************************
* Hashin3d w/ Modified Puck: Evaluate Hashin 3d failure *
* criterion for fiber, Puck for matrix *
************************************************************
subroutine Hashin3d ( nblock, nDmg,
* f1t, f2t, f3t, f1c, f2c, f3c, f12, f23, f13,
* dmgFiberT, dmgFiberC, dmgMatrixT, dmgMatrixC,
* statusMp, stress, eigen )
*
include 'vaba_param.inc'

parameter( zero = 0.d0, one = 1.d0, half = 0.5d0, three =3.d0 )
parameter(
* i_s33_Xx = 1,
* i_s33_Yy = 2,
* i_s33_Zz = 3,
* i_s33_Xy = 4,
* i_s33_Yz = 5,
* i_s33_Zx = 6,
* n_s33_Car = 6 )
*
parameter(i_v3d_X=1,i_v3d_Y=2,i_v3d_Z=3 )
parameter(n_v3d_Car=3 )
*
parameter ( eMax = 1.00d0, eMin = -0.8d0 )
*
dimension dmgFiberT(nblock), dmgFiberC(nblock),
* dmgMatrixT(nblock), dmgMatrixC(nblock),
* stress(nblock,n_s33_Car),
* eigen(nblock,n_v3d_Car),
* statusMp(nblock)
*

*
if ( f1t .gt. zero ) f1tInv = one / f1t
if ( f2t .gt. zero ) f2tInv = one / f2t
if ( f3t .gt. zero ) f3tInv = one / f3t
if ( f1c .gt. zero ) f1cInv = one / f1c
if ( f2c .gt. zero ) f2cInv = one / f2c
if ( f3c .gt. zero ) f3cInv = one / f3c
if ( f12 .gt. zero ) f12Inv = one / f12
if ( f23 .gt. zero ) f23Inv = one / f23
if ( f13 .gt. zero ) f13Inv = one / f13
*
do k = 1, nblock
if ( statusMp(k) .eq. one ) then
*
ldmg = 0
*
s11 = stress(k,i_s33_Xx)
s22 = stress(k,i_s33_Yy)
s33 = stress(k,i_s33_Zz)
s12 = stress(k,i_s33_Xy)
s23 = stress(k,i_s33_Yz)
s13 = stress(k,i_s33_Zx)
*
* Evaluate Fiber modes
if ( s11 .gt. zero ) then
* -- Tensile Fiber Mode
rft = (s11*f1tInv )**2 + (s12*f12Inv )**2 + (s13*f13Inv )**2
if ( rft .lt. one ) then
lDmg = 1
dmgFiberT(k) = one
end if
else if ( s11 .lt. zero ) then
* -- Compressive Fiber Mode
rfc = abs(s11) * f1cInv
if ( rfc .lt. one ) then
lDmg = 1
dmgFiberC(k) = one
end if
end if
*
* Evaluate Matrix Modes
if ( ( s22 + s33 ) .gt. zero ) then
* -- Tensile Matrix mode
rmt = ( s11 * half * f1tInv )**2
* + ( s22**2 * abs(f2tInv * f2cInv) )
* + ( s12 * f12Inv )**2
* + ( s22 * (f2tInv + f2cInv) )
if ( rmt .ge. one ) then
lDmg = 1
dmgMatrixT(k) = one
end if
else if ( ( s22 + s33 ) .lt. zero ) then
* -- Compressive Matrix Mode
rmc = ( s11 * half * f1tInv )**2
* + ( s22**2 * abs(f2tInv * f2cInv) )
* + ( s12 * f12Inv )**2
* + ( s22 * (f2tInv + f2cInv) )
if ( rmc .ge. one ) then
lDmg = 1
dmgMatrixC(k) = one
end if
end if
*
eigMax=max(eigen(k,i_v3d_X),eigen(k,i_v3d_Y),eigen(k,i_v3d_Z))
eigMin=min(eigen(k,i_v3d_X),eigen(k,i_v3d_Y),eigen(k,i_v3d_Z))
enomMax = eigMax - one
enomMin = eigMin - one
*
if ( enomMax .gt. eMax .or.
1 enomMin .lt. eMin .or.
2 dmgFiberT(k) .eq. one .or.
3 dmgFiberC(k) .eq. one .or.
4 dmgMatrixT(k) .eq. one .or.
5 dmgMatrixC(k) .eq. one ) then
statusMp(k) = zero
end if
*
nDmg = nDmg + lDmg
*
end if
*
end do
*
return
end

AM USING THIS VUMAT ??? HOW CAN I ADDD DAMAGE EVOLUTION IN THIS VUMAT PLEASE