编辑代码

      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,KSTEP,KINC)
C
      INCLUDE 'ABA_PARAM.INC'
C 
      CHARACTER*80 CMNAME
      DIMENSION STRESS(NTENS),STATEV(NSTATV),
     1 DDSDDE(NTENS,NTENS),
     2 DDSDDT(NTENS),DRPLDE(NTENS),
     3 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
     4 PROPS(NPROPS),COORDS(2),DROT(2,2),DFGRD0(2,2),DFGRD1(2,2)
      double precision E,V,yield0,H,lamdba,G,STATEV
      double precision hydroStaticPressure,misesEqualStress
      double precision dEqualPlasticStrain
      double precision,  dimension(NTENS) : :deviatoricStress
      integer i,j
C    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 (NTENS表示应力分量/应变分量的数量，三维=6，二维=3；NDI表示独立的正应变数量，三维=3，二维=2）
    do=1, NTENS
        do j=1, NTENS
            DDSDDE(i,j) = lamdba
    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
        do j=1,NTENS
            stress(i)=stress(i)+DDSDDE(i,j)*DSTRAN(j)
        end do
    end do
        
C       Calculate hydroStaticPressure
    hydroStaticPressure=0
    do i=1,NDI
        hydroStaticPressure=hydroStaticPressure+stress(i)
    end do
        hydroStaticPressure=hydroStaticPressure/3
        
C       Calculate the deviatoricStress-try        
    do i=1,NDI
        deviatoricStress(i)=stress(i)-hydroStaticPressure
    end do
    do i=1+NDI,NTENS
        deviatoricStress(i)=stress(i)
    end do

C       Calculate misesEqualStress
    misesEqualStress=0
    do i=1,NDI
        misesEqualStress=miseEqualStress+deviatoricStress(i)**2
    end do
    misesEqualStress=SQRT(1.5*misesEqualStress)
        
C       get the nth equalPlasticStrain form STATEV(值依赖于材料模型中所需跟踪的状态变量数量)
    equalPlasticStrain=STATEV(1)
        
C       Calculate the nth yield
    yieldn=yield0+H*eqaulPlasticStrain
        
    if (misesEqualStress>yieldn) then
        dEqualPlasticStrain=（misesEqualStress-yieldn)/(H+3*G)
           
        do i=1,NTENS
            deviatoricStress(i)=deviatoricStress(i)-3*G*dEqualPlasticStrain*deviatoricStress
        end do
           
        do i=1,NDI
            stress(i)=deviatoricStress(i)+hydroStaticPressure
        end do
        do i=1+NDI,NTENS
            stress(i)=deviatoricStress(i)
        end do
    else
        dEqualPlasticStrain=0
    end if
        
    STATEV(1)=equalPlasticStrain+dEqualPlasticStrain
                
      RETURN
      END