切线刚度法修正

      SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,     1 RPL,DDSDDT,DRPLDE,DRPLDT,     2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,     3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,     4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,JSTEP,KINC)C      INCLUDE 'ABA_PARAM.INC'C      CHARACTER*80 CMNAME      DIMENSION STRESS(NTENS),STATEV(NSTATV),     1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),     2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),     3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3),     4 JSTEP(4)      double precision E,V,yield0,H,lamdba,G,STATEV,M,W,R      double precision hydroStaticPressure,misesEqualStress      double precision hydroStaticPressureTry      double precision dEqualPlasticStrain,misesEqualStressTry      double precision, DIMENSION(NTENS)::deviatoricStress,stressTry      double precision, DIMENSION(NTENS)::deviatoricStressTry      INTEGER i,jC     Elastic properties      E = PROPS(1)      V = PROPS(2)      yield0 = PROPS(3)      H = (PROPS(5)-PROPS(3))/(PROPS(6)-PROPS(4))
      lamdba = E*V/((1+V)*(1-2*V))      G = E/(2*(1+V)) 
C     Calculate the material jacobian matrix      do i=1, NTENS        do j=1, NTENS          DDSDDE(i,j)=0        end do      end do
      do i=1,NDI        do j=1,NDI          DDSDDE(i,j)=lamdba        end do         DDSDDE(i,i)=lamdba+2*G        DDSDDE(i+NDI,i+NDI)=G      end do  C     Calculate stress-try      do i=1,ntens        stressTry(i)=stress(i)      end do       do i=1,NTENS        do j=1,NTENS          stressTry(i)=stressTry(i)+DDSDDE(i,j)*DSTRAN(j)        end do       end do 
      call GetMisesEqualStress(stressTry,NTENS,NDI,misesEqualStressTry     1   ,deviatoricStressTry,hydroStaticPressureTry)
C     get the nth equalPlasticStrain from STATEV      equalPlasticStrain=STATEV(1)
C     Calculate the nth yield      yieldn=yield0+H*equalPlasticStrain
      if (misesEqualStressTry>yieldn) then
      call GetMisesEqualStress(stress,NTENS,NDI,misesEqualStress     1   ,deviatoricStress,hydroStaticPressure)
      if (yieldn==yield0) then C         if (NOEL==400) then C           if(NPT==1) then C             write(*,*) "misesEqualStress=",misesEqualStressC             write(*,*) "misesEqualStressTry=",misesEqualStressTryC             write(*,*) "yieldn=",yieldnC             write(*,*) "yield0=",yield0C           end ifC         end if          R=(yieldn-misesEqualStress)/(misesEqualStressTry-misesEqualStress)          do i=1,NTENS            do j=1,NTENS              stress(i)=stress(i)+R*DDSDDE(i,j)*dstran(j)            end do          end do          call GetMisesEqualStress(stress,NTENS,NDI,misesEqualStress     1   ,deviatoricStress,hydroStaticPressure)          do i=1,NTENS            dstran(i)=(1-R)*dstran(i)          end do      end if            dEqualPlasticStrain=0      M=3*G/(misesEqualStress*(H+3*G))      do i=1,NTENS        dEqualPlasticStrain=dEqualPlasticStrain+M*deviatoricStress(i)*dstran(i)      end do
C    update jacobian matrix      w = 9*G**2/(misesEqualStress**2*(H+3*G))      do i = 1,NTENS        do j=1,ntens          DDSDDE(i,j)=DDSDDE(i,j)-w*deviatoricStress(i)*     &deviatoricStress(j)         end do       end do 
C     Calculate stress-(n+1)      do i=1,NTENS        do j=1,NTENS          stress(i)=stress(i)+DDSDDE(i,j)*dstran(j)        end do       end do              else      do i=1,NTENS        do j=1,NTENS          stress(i)=stress(i)+DDSDDE(i,j)*dstran(j)        end do       end do
      dEqualPlasticStrain=0      end if 
      STATEV(1)=equalPlasticStrain+dEqualPlasticStrain
      RETURN      END
      SUBROUTINE GetMisesEqualStress(STRESS,NTENS,NDI,misesEqualStress     1   ,deviatoricStress,hydroStaticPressure)      INCLUDE 'ABA_PARAM.INC'      DIMENSION STRESS(NTENS)      double precision, DIMENSION(NTENS)::deviatoricStress      double precision misesEqualStress,hydroStaticPressure      INTEGER NTENS,NDI      hydroStaticPressure=0      do i=1,NDI         hydroStaticPressure=hydroStaticPressure+STRESS(i)      enddo      hydroStaticPressure=hydroStaticPressure/3      do i=1,NDI         deviatoricStress(i)=STRESS(i)-hydroStaticPressure      enddo      do i=1+NDI,NTENS         deviatoricStress(i)=STRESS(i)      enddo      misesEqualStress=0      do i=1,NDI         misesEqualStress=misesEqualStress+deviatoricStress(i)**2      enddo      do i=1+NDI,NTENS         misesEqualStress=misesEqualStress+2*deviatoricStress(i)**2      enddo      misesEqualStress=sqrt(1.5*misesEqualStress)      END