根据论坛最大应力准则改写VUMAT，计算结果不对，大佬们帮忙看看

zhuhao3802

************************************************************
*   define the material of matrix            *
************************************************************
      subroutine VUMAT(
     1  nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
     2  stepTime, totalTime, dt, cmname, coordMp, charLength,
     3  props, density, strainInc, relSpinInc,
     4  tempOld, stretchOld, defgradOld, fieldOld,
     5  stressOld, stateOld, enerInternOld, enerInelasOld,
     6  tempNew, stretchNew, defgradNew, fieldNew,
c Write only -
     7  stressNew, stateNew, enerInternNew, enerInelasNew )
c
      include 'vaba_param.inc'
      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)   
c
       character*80 cmname
   
      parameter( zero = 0.d0, one = 1.d0, two = 2.d0, half = .5d0 )
*
      parameter(
     *     i_svd_DmgMatrixT       = 1,
     *     i_svd_DmgMatrixC       = 2,
     *     i_svd_DmgMatrixS       = 3,
     *     i_svd_statusMp    = 4,
     *     i_svd_dampStress  = 5,
c     *    i_svd_dampStressXx = 5,
c     *    i_svd_dampStressYy = 6,
c     *    i_svd_dampStressZz = 7,
c     *    i_svd_dampStressXy = 8,
c     *    i_svd_dampStressYz = 9,
c     *    i_svd_dampStressZx = 10,
     *     i_svd_Strain   = 11,
c     *    i_svd_StrainXx = 11,
c     *    i_svd_StrainYy = 12,
c     *    i_svd_StrainZz = 13,
c     *    i_svd_StrainXy = 14,
c     *    i_svd_StrainYz = 15,
c     *    i_svd_StrainZx = 16,
     *     n_svd_required = 16 )
*
      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 )
*
*  the metarial of matrix is ragard as isotropic
*
      parameter (
     *     i_pro_E         = 1,
     *     i_pro_nu        = 2,
     *     i_pro_sigu1t    = 3,            
     *     i_pro_sigu1c    = 4,            
     *     i_pro_sigus     = 5,
*
     *     i_pro_beta  = 6)
*
*     *     i_pro_rate   = 7,               
*     *     i_pro_c      = 8 )               

* Temporary arrays
*
* Read material properties
*
      E     = props(i_pro_E)
      xnu   = props(i_pro_nu)
*
      beta = props(i_pro_beta)
*      
      G=half*E/(one+xnu)
      gg = one / ( one - xnu*xnu - xnu*xnu - xnu*xnu
     *     - two*xnu*xnu*xnu )
      C11  = E * ( one - xnu*xnu ) * gg
      C22  = E * ( one - xnu*xnu ) * gg
      C33  = E * ( one - xnu*xnu ) * gg
      C12  = E * ( xnu + xnu*xnu ) * gg
      C13  = E * ( xnu + xnu*xnu ) * gg
      C23  = E * ( xnu + xnu*xnu ) * gg
*
      f1t = props(i_pro_sigu1t)
      f1c = props(i_pro_sigu1c)
      f1s = props(i_pro_sigus)

*
* 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 isoelastic ( nblock,
     *        stateOld(1,i_svd_DmgMatrixT),
     *        stateOld(1,i_svd_DmgMatrixC),stateOld(1,i_svd_DmgMatrixS),
     *        strainInc, stressNew,
     *        C11, C22, C33, C12, C23, C13, G)
         return
      end if
*
*  应变更新
      call strainUpdate ( nblock,strainInc,
     *     stateOld(1,i_svd_strain), stateNew(1,i_svd_strain) )
*
* 应力更新
      call isoelastic ( nblock,     
     *     stateOld(1,i_svd_DmgMatrixT),
     *     stateOld(1,i_svd_DmgMatrixC),
     *     stateOld(1,i_svd_DmgMatrixS),
     *     stateNew(1,i_svd_strain), stressNew,
     *     C11, C22, C33, C12, C23, C13, G)
* Failure evaluation
      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) )
      call copyr ( nblock,
     *     stateOld(1,i_svd_DmgMatrixS), stateNew(1,i_svd_DmgMatrixS) )
      nDmg = 0
      call Failure3d  ( nblock, nDmg,
     *     f1t, f1c, f1s,
     *     stateOld(1,i_svd_DmgMatrixT),
     *     stateOld(1,i_svd_DmgMatrixC),
     *     stateOld(1,i_svd_DmgMatrixS),      
     *     stateOld(1,i_svd_statusMp),
     *     stressNew)
*     -- Recompute stresses if new Damage is occurring
      if ( nDmg .gt. 0 ) then
         call isoelastic ( nblock,
     *        stateNew(1,i_svd_DmgMatrixT),
     *        stateNew(1,i_svd_DmgMatrixC),
     *        stateNew(1,i_svd_DmgMatrixS),
     *        stateNew(1,i_svd_strain), stressNew,
     *        C11, C22, C33, C12, C23, C13, G)
      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 isoelastic ( nblock, strain, stress,
     *     DmgMatrixC, DmgMatrixS,DmgMatrixT,
     *     C11, C22, C33, C12, C23, C13, G)
*
      include 'vaba_param.inc'
      DOUBLE PRECISION :AMAGE,MaxDam
*  isotropic elasticity, 3D case - （刚度退化计算）
*
      parameter( zero = 0.d0, one = 1.d0, two = 2.d0,ten=10.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), stress(nblock,n_s33_Car),
     *     DmgMatrixT(nblock), DmgMatrixC(nblock), DmgMatrixS(nblock)
*     
      DAMAGE = ZERO
      MaxDam = ZERO
      do k = 1, nblock
*     -- Stress update
          MaxDam=max(DmgMatrixT(k),DmgMatrixT(k),DmgMatrixT(k))
          if (MaxDam.eq.one)then
              DAMAGE=one/ten  
          else
              DAMAGE=ZERO
          endif
*        
         dC11 = ( one - DAMAGE ) * C11
         dC22 = ( one - DAMAGE ) * C22
         dC33 = ( one - DAMAGE ) * C33
         dC12 = ( one - DAMAGE ) * C12
         dC23 = ( one - DAMAGE ) * C23
         dC13 = ( one - DAMAGE ) * C13
         dG = ( one - DAMAGE ) * G
         dG = ( one - DAMAGE ) * G
         dG = ( one - DAMAGE ) * G
*     -- 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 * dG * strain(k,i_s33_Xy)
         stress(k,i_s33_Yz) = two * dG * strain(k,i_s33_Yz)
         stress(k,i_s33_Zx) = two * dG * 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      
************************************************************
*                 （失效判据）
************************************************************
      subroutine Failure3d ( nblock, nDmg,
     *     f1t, f1c, f1s, DmgMatrixT, DmgMatrixC,DmgMatrixS,  
     *     statusMp, stress)
*
      include 'vaba_param.inc'
      PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,NINE=9.0D0,FOUR=4.D0)
      PARAMETER (PI=3.141592653589,TOL=0.0001)   
      DOUBLE PRECISION :: I1,I2,I3,A,B,C,E1,E2,Fai,SP1,SP2,SP3,SSP1,
     1                    SSP2,SSP3,S1,S2,ST
      parameter ( eMax = 1.d0, eMin = -0.5d0 )
*
      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 )
*     
      dimension  DmgMatrixS(nblock),
     *     DmgMatrixT(nblock), DmgMatrixC(nblock),
     *     stress(nblock,n_s33_Car),
     *     statusMp(nblock)      
c      
c  
      do k = 1, nblock
        if ( statusMp(k) .eq. one ) then
            
         I1 = stress(k,1)+stress(k,2)+stress(k,3)
         I2 = stress(k,1)*stress(k,2)
     1        +stress(k,2)*stress(k,3)
     2        +stress(k,1)*stress(k,3)
     3        -stress(k,4)**TWO-stress(k,5)**TWO
     4        -stress(k,6)**TWO
         I3 = stress(k,1)*stress(k,2)*stress(k,3)
     1        -stress(k,1)*stress(k,5)**TWO
     2        -stress(k,2)*stress(k,6)**TWO
     3        -stress(k,3)*stress(k,4)**TWO
     4        +TWO*stress(k,4)*stress(k,5)*stress(k,6)
C
C calculate principal stress
         A = I1/THREE
         E1 = I1**TWO-THREE*I2
         IF(E1.GT.TOL)THEN
            B=SQRT(E1)
         ELSE
C         WRITE(*,*) 'Warning: B is less than 0.01,set to 0'
            B=ZERO
         ENDIF
C      
         IF(B.GT.0)THEN
            E2=(TWO*I1**THREE-NINE*I1*I2+NINE*THREE*I3)/
     1          (B**(THREE/TWO)*TWO)
           IF(E2.GT.ONE)THEN
C          WRITE(*,*) 'Warning:cos(Fai)>1,E2= ', E2, 'adjust to 1'
              E2=ONE
           ELSEIF(E2.LT.-ONE)THEN
C           WRITE(*,*) 'Warning:cos(Fai)<-1,E2= ', E2, 'adjust to -1'
              E2=-ONE
           ENDIF
           Fai=ACOS(E2)/THREE
           SP1=A+TWO*B*COS(Fai)/THREE
           SP2=A+TWO*B*COS(Fai+TWO*PI/THREE)/THREE
           SP3=A+TWO*B*COS(Fai+FOUR*PI/THREE)/THREE
         ELSE
           SP1=A
           SP2=A
           SP3=A
         ENDIF
         S1=SP1+SP2+SP3
         S2=SP1*SP2+SP2*SP3+SP3*SP1
         ST=SQRT(ABS(TWO*S1**TWO/NINE-TWO*S2/THREE))
C
C  make sure SP1>=SP2>=SP3   
         IF((SP1.GT.SP2).and.(SP2.GT.SP3))THEN
            SSP1=SP1
            SSP2=SP2
            SSP3=SP3
         ELSEIF ((SP1.GT.SP3).and.(SP3.GT.SP2))THEN
            SSP1=SP1
            SSP2=SP3
            SSP3=SP2
         ELSEIF ((SP2.GT.SP1).and.(SP1.GT.SP3))THEN
            SSP1=SP2
            SSP2=SP1
            SSP3=SP3
         ELSEIF ((SP2.GT.SP3).and.(SP3.GT.SP1))THEN
            SSP1=SP2
            SSP2=SP3
            SSP3=SP1
         ELSEIF ((SP3.GT.SP1).and.(SP1.GT.SP2))THEN
            SSP1=SP3
            SSP2=SP1
            SSP3=SP2
         ELSEIF ((SP3.GT.SP2).and.(SP2.GT.SP1))THEN
            SSP1=SP3
            SSP2=SP2
            SSP3=SP1
         ENDIF
C         write(*,*),'A=',SSP1,SSP2,SSP3
c   
         IF (SSP3.GT.ZERO.AND.SSP1.GT.f1t) THEN
             matrixdmg = one
             DmgMatrixT(k) = one
         END IF
C
         IF (SSP1.LT.ZERO.AND.ABS(SSP3).GT.f1c) THEN
             matrixdmg = one
             DmgMatrixC(k) = one
         END IF
c      
         IF (ABS(ST).GT.f1s) THEN
             matrixdmg = one
             DmgMatrixS(k) = one
         END IF
*
         nDmg =  matrixdmg
*            
         if (DmgMatrixT(k).eq.one.or.
     *       DmgMatrixT(k).eq.one.or.
     *       DmgMatrixT(k).eq.one)then
             statusMp(k) = zero
         end if
       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

zhuhao3802

有大佬帮忙看下么