Hashin damage theory

classic Classic list List threaded Threaded
11 messages Options
Reply | Threaded
Open this post in threaded view
|

Hashin damage theory

Rubén Garrorena
Hi all,

I have some questions.

With the Hashin damage model, Can I remove the damage element? How could I do it?

The analysis is dinamic/explicit.

If I make a routine in fortran, how can I usen in Abaqus?


Thanks in advance.

Ruben

[Non-text portions of this message have been removed]




Community email addresses:
  Post message: [hidden email]
  Subscribe:    [hidden email]
  Unsubscribe:  [hidden email]
  List owner:   [hidden email]

Shortcut URL to this page:
  http://groups.yahoo.com/group/abaqus 
Yahoo! Groups Links

<*> To visit your group on the web, go to:
    http://groups.yahoo.com/group/ABAQUS/

<*> Your email settings:
    Individual Email | Traditional

<*> To change settings online go to:
    http://groups.yahoo.com/group/ABAQUS/join
    (Yahoo! ID required)

<*> To change settings via email:
    mailto:[hidden email]
    mailto:[hidden email]

<*> To unsubscribe from this group, send an email to:
    [hidden email]

<*> Your use of Yahoo! Groups is subject to:
    http://docs.yahoo.com/info/terms/
 

Reply | Threaded
Open this post in threaded view
|

Re: Hashin damage theory

Benz-2
Here is a VUMAT provided by Abaqus on Hashin's damage in fiber
reinforced composites :

       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)
*
      f1tInv = zero
      f2tInv = zero
      f3tInv = zero
      f1cInv = zero
      f2cInv = zero
      f3cInv = zero
      f12Inv = zero
      f23Inv = zero
      f13Inv = zero
*
      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
*    
         lFail = 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 .ge. one ) then
               lDmg = 1
               dmgFiberT(k) = one
            end if
         else if ( s11 .lt. zero ) then
*     -- Compressive Fiber Mode
            rfc = abs(s11) * f1cInv
            if ( rfc .ge. 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.
     *        enomMin .lt. eMin .or.
     *        dmgFiberT(k) .eq. one ) then
            statusMp(k) = zero
         end if
*
         nDmg = nDmk + lDmg
*
         end if
*
      end do
*
      return
      end


************************************************************
*   betaDamping: Add beta damping                          *
************************************************************
      subroutine betaDamping3d ( nblock,
     *     beta, dt, strainInc, sigOld, sigNew,
     *     statusMp, sigDampOld, sigDampNew )
*
      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  sigOld(nblock,n_s33_Car),
     *     sigNew(nblock,n_s33_Car),
     *     strainInc(nblock,n_s33_Car),
     *     statusMp(nblock),
     *     sigDampOld(nblock,n_s33_Car),
     *     sigDampNew(nblock,n_s33_Car)      
*
      parameter ( zero = 0.d0, one = 1.d0, two=2.0d0,
     *     half = 0.5d0, third = 1.d0/3.d0 )
      parameter ( asmall = 1.d-16 )
*    
      betaddt =  beta / dt
*
      do k =1 , nblock
         sigDampNew(k,i_s33_Xx) = betaddt * statusMp(k) *
     *        ( sigNew(k,i_s33_Xx)
     *        - ( sigOld(k,i_s33_Xx) - sigDampOld(k,i_s33_Xx) ) )
         sigDampNew(k,i_s33_Yy) = betaddt * statusMp(k) *
     *        ( sigNew(k,i_s33_Yy)
     *        - ( sigOld(k,i_s33_Yy) - sigDampOld(k,i_s33_Yy) ) )
         sigDampNew(k,i_s33_Zz) = betaddt * statusMp(k) *
     *        ( sigNew(k,i_s33_Zz)
     *        - ( sigOld(k,i_s33_Zz) - sigDampOld(k,i_s33_Zz) ) )
         sigDampNew(k,i_s33_Xy) = betaddt * statusMp(k) *
     *        ( sigNew(k,i_s33_Xy)
     *        - ( sigOld(k,i_s33_Xy) - sigDampOld(k,i_s33_Xy) ) )
         sigDampNew(k,i_s33_Yz) = betaddt * statusMp(k) *
     *        ( sigNew(k,i_s33_Yz)
     *        - ( sigOld(k,i_s33_Yz) - sigDampOld(k,i_s33_Yz) ) )
         sigDampNew(k,i_s33_Zx) = betaddt * statusMp(k) *
     *        ( sigNew(k,i_s33_Zx)
     *        - ( sigOld(k,i_s33_Zx) - sigDampOld(k,i_s33_Zx) ) )
*
         sigNew(k,i_s33_Xx) = sigNew(k,i_s33_Xx)+sigDampNew
(k,i_s33_Xx)
         sigNew(k,i_s33_Yy) = sigNew(k,i_s33_Yy)+sigDampNew
(k,i_s33_Yy)
         sigNew(k,i_s33_Zz) = sigNew(k,i_s33_Zz)+sigDampNew
(k,i_s33_Zz)
         sigNew(k,i_s33_Xy) = sigNew(k,i_s33_Xy)+sigDampNew
(k,i_s33_Xy)
         sigNew(k,i_s33_Yz) = sigNew(k,i_s33_Yz)+sigDampNew
(k,i_s33_Yz)
         sigNew(k,i_s33_Zx) = sigNew(k,i_s33_Zx)+sigDampNew
(k,i_s33_Zx)
*
      end do
*    
      return
      end


************************************************************
*   EnergyInternal3d: Compute internal energy for 3d case  *
************************************************************
      subroutine EnergyInternal3d(nblock, sigOld, sigNew ,
     *   strainInc, curDensity, enerInternOld, enerInternNew)
*
      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 )
*
      parameter( two = 2.d0, half = .5d0 )
*
      dimension sigOld (nblock,n_s33_Car), sigNew (nblock,n_s33_Car),
     *     strainInc (nblock,n_s33_Car), curDensity (nblock),
     *     enerInternOld(nblock), enerInternNew(nblock)
*
      do k = 1, nblock
         stressPower  = half * (
     *        ( sigOld(k,i_s33_Xx) + sigNew(k,i_s33_Xx) )
     *        * ( strainInc(k,i_s33_Xx) )
     *        +       ( sigOld(k,i_s33_Yy) + sigNew(k,i_s33_Yy) )
     *        * ( strainInc(k,i_s33_Yy))
     *        +       ( sigOld(k,i_s33_Zz) + sigNew(k,i_s33_Zz) )
     *        * ( strainInc(k,i_s33_Zz))
     *        + two * ( sigOld(k,i_s33_Xy) + sigNew(k,i_s33_Xy) )
     *        * strainInc(k,i_s33_Xy)
     *        + two * ( sigOld(k,i_s33_Yz) + sigNew(k,i_s33_Yz) )
     *        * strainInc(k,i_s33_Yz)
     *        + two * ( sigOld(k,i_s33_Zx) + sigNew(k,i_s33_Zx) )
     *        * strainInc(k,i_s33_Zx) )
*    
         enerInternNew(k) = enerInternOld(k) + stressPower/curDensity
(k)
      end do
*    
      return  
      end  

************************************************************
*   CopyR: Copy from one array to another                  *
************************************************************
      subroutine CopyR(nCopy, from, to )
*
      include 'vaba_param.inc'
*
      dimension from(nCopy), to(nCopy)
*
      do k = 1, nCopy
         to(k) = from(k)
      end do
*
      return
      end

*********************************************************************
*****
* eig33Anal: Compute eigen values of a 3x3 symmetric matrix
analytically *
*********************************************************************
*****
      subroutine eig33Anal( nblock, sMat, eigVal )
*
      include 'vaba_param.inc'
*
      parameter(i_s33_Xx=1,i_s33_Yy=2,i_s33_Zz=3 )
      parameter(i_s33_Xy=4,i_s33_Yz=5,i_s33_Zx=6 )
      parameter(i_s33_Yx=i_s33_Xy )
      parameter(i_s33_Zy=i_s33_Yz )
      parameter(i_s33_Xz=i_s33_Zx,n_s33_Car=6 )
*
      parameter(i_v3d_X=1,i_v3d_Y=2,i_v3d_Z=3 )
      parameter(n_v3d_Car=3 )
*
      parameter ( zero = 0.d0, one = 1.d0, two = 2.d0,
     *     three = 3.d0, half = 0.5d0, third = one / three,
     *     pi23 = 2.094395102393195d0,
     *     fuzz = 1.d-8,
     *     preciz = fuzz * 1.d4 )
*
      dimension eigVal(nblock,n_v3d_Car), sMat(nblock,n_s33_Car)
*
      do k = 1, nblock
        sh  = third*(sMat(k,i_s33_Xx)+sMat(k,i_s33_Yy)+sMat
(k,i_s33_Zz))
        s11 = sMat(k,i_s33_Xx) - sh
        s22 = sMat(k,i_s33_Yy) - sh
        s33 = sMat(k,i_s33_Zz) - sh
        s12 = sMat(k,i_s33_Xy)
        s13 = sMat(k,i_s33_Xz)
        s23 = sMat(k,i_s33_Yz)
*
        fac  = max(abs(s11), abs(s22), abs(s33))
        facs = max(abs(s12), abs(s13), abs(s23))
        if( facs .lt. (preciz*fac) ) then
          eigVal(k,i_v3d_X) = sMat(k,i_s33_Xx)
          eigVal(k,i_v3d_Y) = sMat(k,i_s33_Yy)
          eigVal(k,i_v3d_Z) = sMat(k,i_s33_Zz)
        else
          q = third*((s12**2+s13**2+s23**2)+half*
(s11**2+s22**2+s33**2))
          fac = two * sqrt(q)
          if( fac .gt. fuzz ) then
            ofac = two/fac
          else
            ofac = zero
          end if
          s11 = ofac*s11
          s22 = ofac*s22
          s33 = ofac*s33
          s12 = ofac*s12
          s13 = ofac*s13
          s23 = ofac*s23
          r = s12*s13*s23
     *         + half*(s11*s22*s33-s11*s23**2-s22*s13**2-s33*s12**2)
          if( r .ge. one-fuzz ) then
            cos1 = -half
            cos2 = -half
            cos3 = one
          else if( r .le. fuzz-one ) then
            cos1 = -one
            cos2 = half
            cos3 = half
          else
            ang = third * acos(r)
            cos1 = cos(ang)
            cos2 = cos(ang+pi23)
            cos3 =-cos1-cos2
          end if
          eigVal(k,i_v3d_X) = sh + fac*cos1
          eigVal(k,i_v3d_Y) = sh + fac*cos2
          eigVal(k,i_v3d_Z) = sh + fac*cos3
        end if
      end do
*
      return
      end





BenZ.








--- In [hidden email], Rubén Garrorena <garrorena@...> wrote:
>
> Hi all,
>
> I have some questions.
>
> With the Hashin damage model, Can I remove the damage element? How
could I do it?

>
> The analysis is dinamic/explicit.
>
> If I make a routine in fortran, how can I usen in Abaqus?
>
>
> Thanks in advance.
>
> Ruben
>
> [Non-text portions of this message have been removed]
>





Community email addresses:
  Post message: [hidden email]
  Subscribe:    [hidden email]
  Unsubscribe:  [hidden email]
  List owner:   [hidden email]

Shortcut URL to this page:
  http://groups.yahoo.com/group/abaqus 
Yahoo! Groups Links

<*> To visit your group on the web, go to:
    http://groups.yahoo.com/group/ABAQUS/

<*> Your email settings:
    Individual Email | Traditional

<*> To change settings online go to:
    http://groups.yahoo.com/group/ABAQUS/join
    (Yahoo! ID required)

<*> To change settings via email:
    mailto:[hidden email]
    mailto:[hidden email]

<*> To unsubscribe from this group, send an email to:
    [hidden email]

<*> Your use of Yahoo! Groups is subject to:
    http://docs.yahoo.com/info/terms/
 

Reply | Threaded
Open this post in threaded view
|

Re: Re: Hashin damage theory

Rubén Garrorena
Thanks BenZ.

This subroutine is great.

Ruben.

  ----- Original Message -----
  From: BenZ
  To: [hidden email]
  Sent: Tuesday, October 24, 2006 8:40 AM
  Subject: [ABAQUS] Re: Hashin damage theory


  Here is a VUMAT provided by Abaqus on Hashin's damage in fiber
  reinforced composites :

  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)
  *
  f1tInv = zero
  f2tInv = zero
  f3tInv = zero
  f1cInv = zero
  f2cInv = zero
  f3cInv = zero
  f12Inv = zero
  f23Inv = zero
  f13Inv = zero
  *
  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
  *
  lFail = 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 .ge. one ) then
  lDmg = 1
  dmgFiberT(k) = one
  end if
  else if ( s11 .lt. zero ) then
  * -- Compressive Fiber Mode
  rfc = abs(s11) * f1cInv
  if ( rfc .ge. 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.
  * enomMin .lt. eMin .or.
  * dmgFiberT(k) .eq. one ) then
  statusMp(k) = zero
  end if
  *
  nDmg = nDmk + lDmg
  *
  end if
  *
  end do
  *
  return
  end

  ************************************************************
  * betaDamping: Add beta damping *
  ************************************************************
  subroutine betaDamping3d ( nblock,
  * beta, dt, strainInc, sigOld, sigNew,
  * statusMp, sigDampOld, sigDampNew )
  *
  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 sigOld(nblock,n_s33_Car),
  * sigNew(nblock,n_s33_Car),
  * strainInc(nblock,n_s33_Car),
  * statusMp(nblock),
  * sigDampOld(nblock,n_s33_Car),
  * sigDampNew(nblock,n_s33_Car)
  *
  parameter ( zero = 0.d0, one = 1.d0, two=2.0d0,
  * half = 0.5d0, third = 1.d0/3.d0 )
  parameter ( asmall = 1.d-16 )
  *
  betaddt = beta / dt
  *
  do k =1 , nblock
  sigDampNew(k,i_s33_Xx) = betaddt * statusMp(k) *
  * ( sigNew(k,i_s33_Xx)
  * - ( sigOld(k,i_s33_Xx) - sigDampOld(k,i_s33_Xx) ) )
  sigDampNew(k,i_s33_Yy) = betaddt * statusMp(k) *
  * ( sigNew(k,i_s33_Yy)
  * - ( sigOld(k,i_s33_Yy) - sigDampOld(k,i_s33_Yy) ) )
  sigDampNew(k,i_s33_Zz) = betaddt * statusMp(k) *
  * ( sigNew(k,i_s33_Zz)
  * - ( sigOld(k,i_s33_Zz) - sigDampOld(k,i_s33_Zz) ) )
  sigDampNew(k,i_s33_Xy) = betaddt * statusMp(k) *
  * ( sigNew(k,i_s33_Xy)
  * - ( sigOld(k,i_s33_Xy) - sigDampOld(k,i_s33_Xy) ) )
  sigDampNew(k,i_s33_Yz) = betaddt * statusMp(k) *
  * ( sigNew(k,i_s33_Yz)
  * - ( sigOld(k,i_s33_Yz) - sigDampOld(k,i_s33_Yz) ) )
  sigDampNew(k,i_s33_Zx) = betaddt * statusMp(k) *
  * ( sigNew(k,i_s33_Zx)
  * - ( sigOld(k,i_s33_Zx) - sigDampOld(k,i_s33_Zx) ) )
  *
  sigNew(k,i_s33_Xx) = sigNew(k,i_s33_Xx)+sigDampNew
  (k,i_s33_Xx)
  sigNew(k,i_s33_Yy) = sigNew(k,i_s33_Yy)+sigDampNew
  (k,i_s33_Yy)
  sigNew(k,i_s33_Zz) = sigNew(k,i_s33_Zz)+sigDampNew
  (k,i_s33_Zz)
  sigNew(k,i_s33_Xy) = sigNew(k,i_s33_Xy)+sigDampNew
  (k,i_s33_Xy)
  sigNew(k,i_s33_Yz) = sigNew(k,i_s33_Yz)+sigDampNew
  (k,i_s33_Yz)
  sigNew(k,i_s33_Zx) = sigNew(k,i_s33_Zx)+sigDampNew
  (k,i_s33_Zx)
  *
  end do
  *
  return
  end

  ************************************************************
  * EnergyInternal3d: Compute internal energy for 3d case *
  ************************************************************
  subroutine EnergyInternal3d(nblock, sigOld, sigNew ,
  * strainInc, curDensity, enerInternOld, enerInternNew)
  *
  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 )
  *
  parameter( two = 2.d0, half = .5d0 )
  *
  dimension sigOld (nblock,n_s33_Car), sigNew (nblock,n_s33_Car),
  * strainInc (nblock,n_s33_Car), curDensity (nblock),
  * enerInternOld(nblock), enerInternNew(nblock)
  *
  do k = 1, nblock
  stressPower = half * (
  * ( sigOld(k,i_s33_Xx) + sigNew(k,i_s33_Xx) )
  * * ( strainInc(k,i_s33_Xx) )
  * + ( sigOld(k,i_s33_Yy) + sigNew(k,i_s33_Yy) )
  * * ( strainInc(k,i_s33_Yy))
  * + ( sigOld(k,i_s33_Zz) + sigNew(k,i_s33_Zz) )
  * * ( strainInc(k,i_s33_Zz))
  * + two * ( sigOld(k,i_s33_Xy) + sigNew(k,i_s33_Xy) )
  * * strainInc(k,i_s33_Xy)
  * + two * ( sigOld(k,i_s33_Yz) + sigNew(k,i_s33_Yz) )
  * * strainInc(k,i_s33_Yz)
  * + two * ( sigOld(k,i_s33_Zx) + sigNew(k,i_s33_Zx) )
  * * strainInc(k,i_s33_Zx) )
  *
  enerInternNew(k) = enerInternOld(k) + stressPower/curDensity
  (k)
  end do
  *
  return
  end

  ************************************************************
  * CopyR: Copy from one array to another *
  ************************************************************
  subroutine CopyR(nCopy, from, to )
  *
  include 'vaba_param.inc'
  *
  dimension from(nCopy), to(nCopy)
  *
  do k = 1, nCopy
  to(k) = from(k)
  end do
  *
  return
  end

  *********************************************************************
  *****
  * eig33Anal: Compute eigen values of a 3x3 symmetric matrix
  analytically *
  *********************************************************************
  *****
  subroutine eig33Anal( nblock, sMat, eigVal )
  *
  include 'vaba_param.inc'
  *
  parameter(i_s33_Xx=1,i_s33_Yy=2,i_s33_Zz=3 )
  parameter(i_s33_Xy=4,i_s33_Yz=5,i_s33_Zx=6 )
  parameter(i_s33_Yx=i_s33_Xy )
  parameter(i_s33_Zy=i_s33_Yz )
  parameter(i_s33_Xz=i_s33_Zx,n_s33_Car=6 )
  *
  parameter(i_v3d_X=1,i_v3d_Y=2,i_v3d_Z=3 )
  parameter(n_v3d_Car=3 )
  *
  parameter ( zero = 0.d0, one = 1.d0, two = 2.d0,
  * three = 3.d0, half = 0.5d0, third = one / three,
  * pi23 = 2.094395102393195d0,
  * fuzz = 1.d-8,
  * preciz = fuzz * 1.d4 )
  *
  dimension eigVal(nblock,n_v3d_Car), sMat(nblock,n_s33_Car)
  *
  do k = 1, nblock
  sh = third*(sMat(k,i_s33_Xx)+sMat(k,i_s33_Yy)+sMat
  (k,i_s33_Zz))
  s11 = sMat(k,i_s33_Xx) - sh
  s22 = sMat(k,i_s33_Yy) - sh
  s33 = sMat(k,i_s33_Zz) - sh
  s12 = sMat(k,i_s33_Xy)
  s13 = sMat(k,i_s33_Xz)
  s23 = sMat(k,i_s33_Yz)
  *
  fac = max(abs(s11), abs(s22), abs(s33))
  facs = max(abs(s12), abs(s13), abs(s23))
  if( facs .lt. (preciz*fac) ) then
  eigVal(k,i_v3d_X) = sMat(k,i_s33_Xx)
  eigVal(k,i_v3d_Y) = sMat(k,i_s33_Yy)
  eigVal(k,i_v3d_Z) = sMat(k,i_s33_Zz)
  else
  q = third*((s12**2+s13**2+s23**2)+half*
  (s11**2+s22**2+s33**2))
  fac = two * sqrt(q)
  if( fac .gt. fuzz ) then
  ofac = two/fac
  else
  ofac = zero
  end if
  s11 = ofac*s11
  s22 = ofac*s22
  s33 = ofac*s33
  s12 = ofac*s12
  s13 = ofac*s13
  s23 = ofac*s23
  r = s12*s13*s23
  * + half*(s11*s22*s33-s11*s23**2-s22*s13**2-s33*s12**2)
  if( r .ge. one-fuzz ) then
  cos1 = -half
  cos2 = -half
  cos3 = one
  else if( r .le. fuzz-one ) then
  cos1 = -one
  cos2 = half
  cos3 = half
  else
  ang = third * acos(r)
  cos1 = cos(ang)
  cos2 = cos(ang+pi23)
  cos3 =-cos1-cos2
  end if
  eigVal(k,i_v3d_X) = sh + fac*cos1
  eigVal(k,i_v3d_Y) = sh + fac*cos2
  eigVal(k,i_v3d_Z) = sh + fac*cos3
  end if
  end do
  *
  return
  end

  BenZ.

  --- In [hidden email], Rubén Garrorena <garrorena@...> wrote:
  >
  > Hi all,
  >
  > I have some questions.
  >
  > With the Hashin damage model, Can I remove the damage element? How
  could I do it?
  >
  > The analysis is dinamic/explicit.
  >
  > If I make a routine in fortran, how can I usen in Abaqus?
  >
  >
  > Thanks in advance.
  >
  > Ruben
  >
  > [Non-text portions of this message have been removed]
  >



   

[Non-text portions of this message have been removed]




Community email addresses:
  Post message: [hidden email]
  Subscribe:    [hidden email]
  Unsubscribe:  [hidden email]
  List owner:   [hidden email]

Shortcut URL to this page:
  http://groups.yahoo.com/group/abaqus 
Yahoo! Groups Links

<*> To visit your group on the web, go to:
    http://groups.yahoo.com/group/ABAQUS/

<*> Your email settings:
    Individual Email | Traditional

<*> To change settings online go to:
    http://groups.yahoo.com/group/ABAQUS/join
    (Yahoo! ID required)

<*> To change settings via email:
    mailto:[hidden email]
    mailto:[hidden email]

<*> To unsubscribe from this group, send an email to:
    [hidden email]

<*> Your use of Yahoo! Groups is subject to:
    http://docs.yahoo.com/info/terms/
 

Reply | Threaded
Open this post in threaded view
|

Re: Hashin damage theory

zgwei2002
In reply to this post by Benz-2
Hi,BenZ,

This user-subroutine looks great. Could you post more or give a link
to more information about 3D/Hashin Damage model/Abaqus/Explicit?
such as papers, documents. Thanks!

Regards,

ZGWEI

--- In [hidden email], "BenZ" <benjamin_hagege@...> wrote:

>
> Here is a VUMAT provided by Abaqus on Hashin's damage in fiber
> reinforced composites :
>
>        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)
> *
>       f1tInv = zero
>       f2tInv = zero
>       f3tInv = zero
>       f1cInv = zero
>       f2cInv = zero
>       f3cInv = zero
>       f12Inv = zero
>       f23Inv = zero
>       f13Inv = zero
> *
>       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
> *    
>          lFail = 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 .ge. one ) then
>                lDmg = 1
>                dmgFiberT(k) = one
>             end if
>          else if ( s11 .lt. zero ) then
> *     -- Compressive Fiber Mode
>             rfc = abs(s11) * f1cInv
>             if ( rfc .ge. 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.
>      *        enomMin .lt. eMin .or.
>      *        dmgFiberT(k) .eq. one ) then
>             statusMp(k) = zero
>          end if
> *
>          nDmg = nDmk + lDmg
> *
>          end if
> *
>       end do
> *
>       return
>       end
>
>
> ************************************************************
> *   betaDamping: Add beta damping                          *
> ************************************************************
>       subroutine betaDamping3d ( nblock,
>      *     beta, dt, strainInc, sigOld, sigNew,
>      *     statusMp, sigDampOld, sigDampNew )
> *
>       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  sigOld(nblock,n_s33_Car),
>      *     sigNew(nblock,n_s33_Car),
>      *     strainInc(nblock,n_s33_Car),
>      *     statusMp(nblock),
>      *     sigDampOld(nblock,n_s33_Car),
>      *     sigDampNew(nblock,n_s33_Car)      
> *
>       parameter ( zero = 0.d0, one = 1.d0, two=2.0d0,
>      *     half = 0.5d0, third = 1.d0/3.d0 )
>       parameter ( asmall = 1.d-16 )
> *    
>       betaddt =  beta / dt
> *
>       do k =1 , nblock
>          sigDampNew(k,i_s33_Xx) = betaddt * statusMp(k) *
>      *        ( sigNew(k,i_s33_Xx)
>      *        - ( sigOld(k,i_s33_Xx) - sigDampOld(k,i_s33_Xx) ) )
>          sigDampNew(k,i_s33_Yy) = betaddt * statusMp(k) *
>      *        ( sigNew(k,i_s33_Yy)
>      *        - ( sigOld(k,i_s33_Yy) - sigDampOld(k,i_s33_Yy) ) )
>          sigDampNew(k,i_s33_Zz) = betaddt * statusMp(k) *
>      *        ( sigNew(k,i_s33_Zz)
>      *        - ( sigOld(k,i_s33_Zz) - sigDampOld(k,i_s33_Zz) ) )
>          sigDampNew(k,i_s33_Xy) = betaddt * statusMp(k) *
>      *        ( sigNew(k,i_s33_Xy)
>      *        - ( sigOld(k,i_s33_Xy) - sigDampOld(k,i_s33_Xy) ) )
>          sigDampNew(k,i_s33_Yz) = betaddt * statusMp(k) *
>      *        ( sigNew(k,i_s33_Yz)
>      *        - ( sigOld(k,i_s33_Yz) - sigDampOld(k,i_s33_Yz) ) )
>          sigDampNew(k,i_s33_Zx) = betaddt * statusMp(k) *
>      *        ( sigNew(k,i_s33_Zx)
>      *        - ( sigOld(k,i_s33_Zx) - sigDampOld(k,i_s33_Zx) ) )
> *
>          sigNew(k,i_s33_Xx) = sigNew(k,i_s33_Xx)+sigDampNew
> (k,i_s33_Xx)
>          sigNew(k,i_s33_Yy) = sigNew(k,i_s33_Yy)+sigDampNew
> (k,i_s33_Yy)
>          sigNew(k,i_s33_Zz) = sigNew(k,i_s33_Zz)+sigDampNew
> (k,i_s33_Zz)
>          sigNew(k,i_s33_Xy) = sigNew(k,i_s33_Xy)+sigDampNew
> (k,i_s33_Xy)
>          sigNew(k,i_s33_Yz) = sigNew(k,i_s33_Yz)+sigDampNew
> (k,i_s33_Yz)
>          sigNew(k,i_s33_Zx) = sigNew(k,i_s33_Zx)+sigDampNew
> (k,i_s33_Zx)
> *
>       end do
> *    
>       return
>       end
>
>
> ************************************************************
> *   EnergyInternal3d: Compute internal energy for 3d case  *
> ************************************************************
>       subroutine EnergyInternal3d(nblock, sigOld, sigNew ,
>      *   strainInc, curDensity, enerInternOld, enerInternNew)
> *
>       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 )
> *
>       parameter( two = 2.d0, half = .5d0 )
> *
>       dimension sigOld (nblock,n_s33_Car), sigNew
(nblock,n_s33_Car),

>      *     strainInc (nblock,n_s33_Car), curDensity (nblock),
>      *     enerInternOld(nblock), enerInternNew(nblock)
> *
>       do k = 1, nblock
>          stressPower  = half * (
>      *        ( sigOld(k,i_s33_Xx) + sigNew(k,i_s33_Xx) )
>      *        * ( strainInc(k,i_s33_Xx) )
>      *        +       ( sigOld(k,i_s33_Yy) + sigNew(k,i_s33_Yy) )
>      *        * ( strainInc(k,i_s33_Yy))
>      *        +       ( sigOld(k,i_s33_Zz) + sigNew(k,i_s33_Zz) )
>      *        * ( strainInc(k,i_s33_Zz))
>      *        + two * ( sigOld(k,i_s33_Xy) + sigNew(k,i_s33_Xy) )
>      *        * strainInc(k,i_s33_Xy)
>      *        + two * ( sigOld(k,i_s33_Yz) + sigNew(k,i_s33_Yz) )
>      *        * strainInc(k,i_s33_Yz)
>      *        + two * ( sigOld(k,i_s33_Zx) + sigNew(k,i_s33_Zx) )
>      *        * strainInc(k,i_s33_Zx) )
> *    
>          enerInternNew(k) = enerInternOld(k) +
stressPower/curDensity

> (k)
>       end do
> *    
>       return  
>       end  
>
> ************************************************************
> *   CopyR: Copy from one array to another                  *
> ************************************************************
>       subroutine CopyR(nCopy, from, to )
> *
>       include 'vaba_param.inc'
> *
>       dimension from(nCopy), to(nCopy)
> *
>       do k = 1, nCopy
>          to(k) = from(k)
>       end do
> *
>       return
>       end
>
>
*********************************************************************
> *****
> * eig33Anal: Compute eigen values of a 3x3 symmetric matrix
> analytically *
>
*********************************************************************

> *****
>       subroutine eig33Anal( nblock, sMat, eigVal )
> *
>       include 'vaba_param.inc'
> *
>       parameter(i_s33_Xx=1,i_s33_Yy=2,i_s33_Zz=3 )
>       parameter(i_s33_Xy=4,i_s33_Yz=5,i_s33_Zx=6 )
>       parameter(i_s33_Yx=i_s33_Xy )
>       parameter(i_s33_Zy=i_s33_Yz )
>       parameter(i_s33_Xz=i_s33_Zx,n_s33_Car=6 )
> *
>       parameter(i_v3d_X=1,i_v3d_Y=2,i_v3d_Z=3 )
>       parameter(n_v3d_Car=3 )
> *
>       parameter ( zero = 0.d0, one = 1.d0, two = 2.d0,
>      *     three = 3.d0, half = 0.5d0, third = one / three,
>      *     pi23 = 2.094395102393195d0,
>      *     fuzz = 1.d-8,
>      *     preciz = fuzz * 1.d4 )
> *
>       dimension eigVal(nblock,n_v3d_Car), sMat(nblock,n_s33_Car)
> *
>       do k = 1, nblock
>         sh  = third*(sMat(k,i_s33_Xx)+sMat(k,i_s33_Yy)+sMat
> (k,i_s33_Zz))
>         s11 = sMat(k,i_s33_Xx) - sh
>         s22 = sMat(k,i_s33_Yy) - sh
>         s33 = sMat(k,i_s33_Zz) - sh
>         s12 = sMat(k,i_s33_Xy)
>         s13 = sMat(k,i_s33_Xz)
>         s23 = sMat(k,i_s33_Yz)
> *
>         fac  = max(abs(s11), abs(s22), abs(s33))
>         facs = max(abs(s12), abs(s13), abs(s23))
>         if( facs .lt. (preciz*fac) ) then
>           eigVal(k,i_v3d_X) = sMat(k,i_s33_Xx)
>           eigVal(k,i_v3d_Y) = sMat(k,i_s33_Yy)
>           eigVal(k,i_v3d_Z) = sMat(k,i_s33_Zz)
>         else
>           q = third*((s12**2+s13**2+s23**2)+half*
> (s11**2+s22**2+s33**2))
>           fac = two * sqrt(q)
>           if( fac .gt. fuzz ) then
>             ofac = two/fac
>           else
>             ofac = zero
>           end if
>           s11 = ofac*s11
>           s22 = ofac*s22
>           s33 = ofac*s33
>           s12 = ofac*s12
>           s13 = ofac*s13
>           s23 = ofac*s23
>           r = s12*s13*s23
>      *         + half*(s11*s22*s33-s11*s23**2-s22*s13**2-
s33*s12**2)

>           if( r .ge. one-fuzz ) then
>             cos1 = -half
>             cos2 = -half
>             cos3 = one
>           else if( r .le. fuzz-one ) then
>             cos1 = -one
>             cos2 = half
>             cos3 = half
>           else
>             ang = third * acos(r)
>             cos1 = cos(ang)
>             cos2 = cos(ang+pi23)
>             cos3 =-cos1-cos2
>           end if
>           eigVal(k,i_v3d_X) = sh + fac*cos1
>           eigVal(k,i_v3d_Y) = sh + fac*cos2
>           eigVal(k,i_v3d_Z) = sh + fac*cos3
>         end if
>       end do
> *
>       return
>       end
>
>
>
>
>
> BenZ.
>
>
>
>
>
>
>
>
> --- In [hidden email], Rubén Garrorena <garrorena@> wrote:
> >
> > Hi all,
> >
> > I have some questions.
> >
> > With the Hashin damage model, Can I remove the damage element?
How

> could I do it?
> >
> > The analysis is dinamic/explicit.
> >
> > If I make a routine in fortran, how can I usen in Abaqus?
> >
> >
> > Thanks in advance.
> >
> > Ruben
> >
> > [Non-text portions of this message have been removed]
> >
>


Reply | Threaded
Open this post in threaded view
|

Re: Hashin damage theory

Benz-2
Yes :
1/The Abaqus doc of version 6.6 ! The VUMAT does exactly the same as
the built-in Hashin damage of Abaqus.
2/The Hashin's paper.
3/I do researches on 3D damage, but I have just started so I have no
other links at the moment. But I can say that you must know the
papers of Ladevèze and his PhD candidate on "damage meso-models".
They are from the well-known LMT, Cachan, laboratory.


BenZ.

--- In [hidden email], "zgwei2002" <zgwei2002@...> wrote:
>
> Hi,BenZ,
>
> This user-subroutine looks great. Could you post more or give a
link

> to more information about 3D/Hashin Damage model/Abaqus/Explicit?
> such as papers, documents. Thanks!
>
> Regards,
>
> ZGWEI
>
> --- In [hidden email], "BenZ" <benjamin_hagege@> wrote:
> >
> > Here is a VUMAT provided by Abaqus on Hashin's damage in fiber
> > reinforced composites :
> >
> >        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)
> > *
> >       f1tInv = zero
> >       f2tInv = zero
> >       f3tInv = zero
> >       f1cInv = zero
> >       f2cInv = zero
> >       f3cInv = zero
> >       f12Inv = zero
> >       f23Inv = zero
> >       f13Inv = zero
> > *
> >       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
> > *    
> >          lFail = 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 .ge. one ) then
> >                lDmg = 1
> >                dmgFiberT(k) = one
> >             end if
> >          else if ( s11 .lt. zero ) then
> > *     -- Compressive Fiber Mode
> >             rfc = abs(s11) * f1cInv
> >             if ( rfc .ge. 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.
> >      *        enomMin .lt. eMin .or.
> >      *        dmgFiberT(k) .eq. one ) then
> >             statusMp(k) = zero
> >          end if
> > *
> >          nDmg = nDmk + lDmg
> > *
> >          end if
> > *
> >       end do
> > *
> >       return
> >       end
> >
> >
> > ************************************************************
> > *   betaDamping: Add beta damping                          *
> > ************************************************************
> >       subroutine betaDamping3d ( nblock,
> >      *     beta, dt, strainInc, sigOld, sigNew,
> >      *     statusMp, sigDampOld, sigDampNew )
> > *
> >       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  sigOld(nblock,n_s33_Car),
> >      *     sigNew(nblock,n_s33_Car),
> >      *     strainInc(nblock,n_s33_Car),
> >      *     statusMp(nblock),
> >      *     sigDampOld(nblock,n_s33_Car),
> >      *     sigDampNew(nblock,n_s33_Car)      
> > *
> >       parameter ( zero = 0.d0, one = 1.d0, two=2.0d0,
> >      *     half = 0.5d0, third = 1.d0/3.d0 )
> >       parameter ( asmall = 1.d-16 )
> > *    
> >       betaddt =  beta / dt
> > *
> >       do k =1 , nblock
> >          sigDampNew(k,i_s33_Xx) = betaddt * statusMp(k) *
> >      *        ( sigNew(k,i_s33_Xx)
> >      *        - ( sigOld(k,i_s33_Xx) - sigDampOld(k,i_s33_Xx) ) )
> >          sigDampNew(k,i_s33_Yy) = betaddt * statusMp(k) *
> >      *        ( sigNew(k,i_s33_Yy)
> >      *        - ( sigOld(k,i_s33_Yy) - sigDampOld(k,i_s33_Yy) ) )
> >          sigDampNew(k,i_s33_Zz) = betaddt * statusMp(k) *
> >      *        ( sigNew(k,i_s33_Zz)
> >      *        - ( sigOld(k,i_s33_Zz) - sigDampOld(k,i_s33_Zz) ) )
> >          sigDampNew(k,i_s33_Xy) = betaddt * statusMp(k) *
> >      *        ( sigNew(k,i_s33_Xy)
> >      *        - ( sigOld(k,i_s33_Xy) - sigDampOld(k,i_s33_Xy) ) )
> >          sigDampNew(k,i_s33_Yz) = betaddt * statusMp(k) *
> >      *        ( sigNew(k,i_s33_Yz)
> >      *        - ( sigOld(k,i_s33_Yz) - sigDampOld(k,i_s33_Yz) ) )
> >          sigDampNew(k,i_s33_Zx) = betaddt * statusMp(k) *
> >      *        ( sigNew(k,i_s33_Zx)
> >      *        - ( sigOld(k,i_s33_Zx) - sigDampOld(k,i_s33_Zx) ) )
> > *
> >          sigNew(k,i_s33_Xx) = sigNew(k,i_s33_Xx)+sigDampNew
> > (k,i_s33_Xx)
> >          sigNew(k,i_s33_Yy) = sigNew(k,i_s33_Yy)+sigDampNew
> > (k,i_s33_Yy)
> >          sigNew(k,i_s33_Zz) = sigNew(k,i_s33_Zz)+sigDampNew
> > (k,i_s33_Zz)
> >          sigNew(k,i_s33_Xy) = sigNew(k,i_s33_Xy)+sigDampNew
> > (k,i_s33_Xy)
> >          sigNew(k,i_s33_Yz) = sigNew(k,i_s33_Yz)+sigDampNew
> > (k,i_s33_Yz)
> >          sigNew(k,i_s33_Zx) = sigNew(k,i_s33_Zx)+sigDampNew
> > (k,i_s33_Zx)
> > *
> >       end do
> > *    
> >       return
> >       end
> >
> >
> > ************************************************************
> > *   EnergyInternal3d: Compute internal energy for 3d case  *
> > ************************************************************
> >       subroutine EnergyInternal3d(nblock, sigOld, sigNew ,
> >      *   strainInc, curDensity, enerInternOld, enerInternNew)
> > *
> >       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 )
> > *
> >       parameter( two = 2.d0, half = .5d0 )
> > *
> >       dimension sigOld (nblock,n_s33_Car), sigNew
> (nblock,n_s33_Car),
> >      *     strainInc (nblock,n_s33_Car), curDensity (nblock),
> >      *     enerInternOld(nblock), enerInternNew(nblock)
> > *
> >       do k = 1, nblock
> >          stressPower  = half * (
> >      *        ( sigOld(k,i_s33_Xx) + sigNew(k,i_s33_Xx) )
> >      *        * ( strainInc(k,i_s33_Xx) )
> >      *        +       ( sigOld(k,i_s33_Yy) + sigNew(k,i_s33_Yy) )
> >      *        * ( strainInc(k,i_s33_Yy))
> >      *        +       ( sigOld(k,i_s33_Zz) + sigNew(k,i_s33_Zz) )
> >      *        * ( strainInc(k,i_s33_Zz))
> >      *        + two * ( sigOld(k,i_s33_Xy) + sigNew(k,i_s33_Xy) )
> >      *        * strainInc(k,i_s33_Xy)
> >      *        + two * ( sigOld(k,i_s33_Yz) + sigNew(k,i_s33_Yz) )
> >      *        * strainInc(k,i_s33_Yz)
> >      *        + two * ( sigOld(k,i_s33_Zx) + sigNew(k,i_s33_Zx) )
> >      *        * strainInc(k,i_s33_Zx) )
> > *    
> >          enerInternNew(k) = enerInternOld(k) +
> stressPower/curDensity
> > (k)
> >       end do
> > *    
> >       return  
> >       end  
> >
> > ************************************************************
> > *   CopyR: Copy from one array to another                  *
> > ************************************************************
> >       subroutine CopyR(nCopy, from, to )
> > *
> >       include 'vaba_param.inc'
> > *
> >       dimension from(nCopy), to(nCopy)
> > *
> >       do k = 1, nCopy
> >          to(k) = from(k)
> >       end do
> > *
> >       return
> >       end
> >
> >
>
*********************************************************************
> > *****
> > * eig33Anal: Compute eigen values of a 3x3 symmetric matrix
> > analytically *
> >
>
*********************************************************************

> > *****
> >       subroutine eig33Anal( nblock, sMat, eigVal )
> > *
> >       include 'vaba_param.inc'
> > *
> >       parameter(i_s33_Xx=1,i_s33_Yy=2,i_s33_Zz=3 )
> >       parameter(i_s33_Xy=4,i_s33_Yz=5,i_s33_Zx=6 )
> >       parameter(i_s33_Yx=i_s33_Xy )
> >       parameter(i_s33_Zy=i_s33_Yz )
> >       parameter(i_s33_Xz=i_s33_Zx,n_s33_Car=6 )
> > *
> >       parameter(i_v3d_X=1,i_v3d_Y=2,i_v3d_Z=3 )
> >       parameter(n_v3d_Car=3 )
> > *
> >       parameter ( zero = 0.d0, one = 1.d0, two = 2.d0,
> >      *     three = 3.d0, half = 0.5d0, third = one / three,
> >      *     pi23 = 2.094395102393195d0,
> >      *     fuzz = 1.d-8,
> >      *     preciz = fuzz * 1.d4 )
> > *
> >       dimension eigVal(nblock,n_v3d_Car), sMat(nblock,n_s33_Car)
> > *
> >       do k = 1, nblock
> >         sh  = third*(sMat(k,i_s33_Xx)+sMat(k,i_s33_Yy)+sMat
> > (k,i_s33_Zz))
> >         s11 = sMat(k,i_s33_Xx) - sh
> >         s22 = sMat(k,i_s33_Yy) - sh
> >         s33 = sMat(k,i_s33_Zz) - sh
> >         s12 = sMat(k,i_s33_Xy)
> >         s13 = sMat(k,i_s33_Xz)
> >         s23 = sMat(k,i_s33_Yz)
> > *
> >         fac  = max(abs(s11), abs(s22), abs(s33))
> >         facs = max(abs(s12), abs(s13), abs(s23))
> >         if( facs .lt. (preciz*fac) ) then
> >           eigVal(k,i_v3d_X) = sMat(k,i_s33_Xx)
> >           eigVal(k,i_v3d_Y) = sMat(k,i_s33_Yy)
> >           eigVal(k,i_v3d_Z) = sMat(k,i_s33_Zz)
> >         else
> >           q = third*((s12**2+s13**2+s23**2)+half*
> > (s11**2+s22**2+s33**2))
> >           fac = two * sqrt(q)
> >           if( fac .gt. fuzz ) then
> >             ofac = two/fac
> >           else
> >             ofac = zero
> >           end if
> >           s11 = ofac*s11
> >           s22 = ofac*s22
> >           s33 = ofac*s33
> >           s12 = ofac*s12
> >           s13 = ofac*s13
> >           s23 = ofac*s23
> >           r = s12*s13*s23
> >      *         + half*(s11*s22*s33-s11*s23**2-s22*s13**2-
> s33*s12**2)
> >           if( r .ge. one-fuzz ) then
> >             cos1 = -half
> >             cos2 = -half
> >             cos3 = one
> >           else if( r .le. fuzz-one ) then
> >             cos1 = -one
> >             cos2 = half
> >             cos3 = half
> >           else
> >             ang = third * acos(r)
> >             cos1 = cos(ang)
> >             cos2 = cos(ang+pi23)
> >             cos3 =-cos1-cos2
> >           end if
> >           eigVal(k,i_v3d_X) = sh + fac*cos1
> >           eigVal(k,i_v3d_Y) = sh + fac*cos2
> >           eigVal(k,i_v3d_Z) = sh + fac*cos3
> >         end if
> >       end do
> > *
> >       return
> >       end
> >
> >
> >
> >
> >
> > BenZ.
> >
> >
> >
> >
> >
> >
> >
> >
> > --- In [hidden email], Rubén Garrorena <garrorena@>
wrote:

> > >
> > > Hi all,
> > >
> > > I have some questions.
> > >
> > > With the Hashin damage model, Can I remove the damage element?
> How
> > could I do it?
> > >
> > > The analysis is dinamic/explicit.
> > >
> > > If I make a routine in fortran, how can I usen in Abaqus?
> > >
> > >
> > > Thanks in advance.
> > >
> > > Ruben
> > >
> > > [Non-text portions of this message have been removed]
> > >
> >
>


Reply | Threaded
Open this post in threaded view
|

Re: Hashin damage theory

MacGyver1980
In reply to this post by Benz-2
Hello there.

This VUMAt is exactly what I need. Unfortunately it seems that this VUMAT doesn't work with Abaqus 6.7 or 6.8. Can anybody tell me what causes this problem? I think this VUMAT was created for 6.6. Is this correct? Where can i get informations about changes between 6.6 and 6.7/6.8? I looked in the Abaqus Documentation but didn't find what i looked for.

Many thanks for your HELP.
Reply | Threaded
Open this post in threaded view
|

Re: Hashin damage theory

franpesca
This post was updated on .
MacGyver1980 wrote
Hello there.

This VUMAt is exactly what I need. Unfortunately it seems that this VUMAT doesn't work with Abaqus 6.7 or 6.8. Can anybody tell me what causes this problem? I think this VUMAT was created for 6.6. Is this correct? Where can i get informations about changes between 6.6 and 6.7/6.8? I looked in the Abaqus Documentation but didn't find what i looked for.

Many thanks for your HELP.

Hi everybody

I'm trying to use the Hashin damage model to simulate the crush of a crfp specimen with abaqus version 6.8
Hashin's damage model is supposed to be implemented inside the program and no more as an external VUMAT, but I need something as the VUMAT posted earlier referred to a 2D case instead of 3D... can anybody help?

thanks a lot
Reply | Threaded
Open this post in threaded view
|

Re: Hashin damage theory

amin1075
This post has NOT been accepted by the mailing list yet.
In reply to this post by Rubén Garrorena
hi
What should I do if I want to change the model of progressive hashin subroutin?
What's to be learned?
vumat, umat, usdfld, ...?
Reply | Threaded
Open this post in threaded view
|

Re: Hashin damage theory

kinvi dossou
This post has NOT been accepted by the mailing list yet.
In reply to this post by Benz-2
Hi,
I tested this subroutine on a model with one element but I always get this error:
 "Bad Material definition in element number 1 instance PART-1-1: zero or negative initial dilatational modulus caused by bad material data. Please check your material input and any initial conditions if necessary."
I think that my materials properties are not well defined.
Could you send me some properties that should work on the model in order to compare them with mine?
Thanks

Reply | Threaded
Open this post in threaded view
|

Re: Hashin damage theory

kinvi dossou
This post has NOT been accepted by the mailing list yet.
In reply to this post by amin1075
Hi, I tested this subroutine on a model with one element but I always get this error: "Bad Material definition in element number 1 instance PART-1-1: zero or negative initial dilatational modulus caused by bad material data. Please check your material input and any initial conditions if necessary." I think that my materials properties are not well defined. Could you send me some properties that should work on the model in order to compare them with mine? Thanks
Reply | Threaded
Open this post in threaded view
|

Re: Hashin damage theory

prakash_chakrapani
This post has NOT been accepted by the mailing list yet.
can you please send me the updated hashin fortran file