一个有意思的材料本构模型设计方案，拉伸变形采用von Mises屈服，压缩侧 cap屈服本构模型设计。 原创

原始文献：《Mechanical modelling of indentation-induced densification in amorphous silica》

  该文章为了模拟非晶态二氧化硅的压缩力学性能，把拉伸与压缩分开处理：拉伸侧采用熟悉的 von Mises 屈服，压缩侧则切换到 cap 屈服面。这样的设计，正好对应了非晶二氧化硅在压痕加载下“既会发生剪切塑性，又会发生永久致密化”的真实特征。

  分享这个代码的主要原因：一方面，它很适合做玻璃、非晶材料、压痕问题中的压力敏感塑性分析；另一方面，它也是学习 cap 模型、致密化硬化和隐式本构积分的一个很好的范例。论文结果表明，这一模型能够较好复现实验载荷—位移曲线以及压痕致密化分布，不过需要明确指出的是，当前模型暂时还没有考虑剪切硬化，因此更适合用于理解“压痕致密化”这一核心机制，而不是直接覆盖所有复杂失效问题。作为一份用于科研复现和二次开发的代码，我觉得它很有参考价值。

原始程序如下：

    SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,

   1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,

   2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS,

   3 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,

   4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,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(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)

C 

    DIMENSION EELAS(NTENS),EPLAS(NTENS),FLOW(NTENS),DEPLAS(NTENS),

   1 PSTRESS(NTENS),DSTRESS(NTENS),DDSDDEEQ(NTENS,NTENS),

   2 FLOW2(NTENS,NTENS),EQSYS(NTENS+1,NTENS+1), BVECT(NTENS+1),

   3 PVECT(NTENS),VINDX(NTENS+1)

    PARAMETER (ONE=1.0,TWO=2.0,THREE=3.0,SIX=6.0, HALF =0.5) 

    DATA NEWTON,TOLER/40,1.D-6/ 

C 

C ----------------------------------------------------------- 

C     Material properties

C ----------------------------------------------------------- 

C     PROPS(1) - Young's modulus 

C     PROPS(2) - Poisson ratio 

C     PROPS(3) - Yield (von Mises) stress at zero pressur

C     PROPS(4) - Initial yield hydro at compression (negative)

C     PROPS(5) - Hardening coefficient

C ----------------------------------------------------------- 

C

C Elastic properties

C

    EMOD=PROPS(1)

    ENU=PROPS(2)

    YIELDS=PROPS(3)

    YHYDRO=PROPS(4)

    HPAR=PROPS(5)

C

    EG=EMOD/(TWO*(ONE+ENU))

    EG2=EG*TWO

    ELAM=EG2*ENU/(ONE-TWO*ENU)

C

C ----------------------------------------------------------- 

C Elastic stiffness tensor

C ----------------------------------------------------------- 

C

C Lirear elastic

C

    DO K1=1, NTENS

    DO K2=1, NTENS

     DDSDDE(K2, K1)=0.0

    END DO

    END DO

C

    DO K1=1, NDI

    DO K2=1, NDI

     DDSDDE(K2, K1)=ELAM

    END DO

    DDSDDE(K1, K1)=EG2+ELAM

    END DO 

C

    DO K1=NDI+1, NTENS

    DDSDDE(K1, K1)=EG

    END DO

C

C ----------------------------------------------------------- 

C  Elastic prediction

C ----------------------------------------------------------- 

C

    DO K2=1, NTENS

    DSTRESS(K2)=0.0

    END DO

    DO K1=1, NTENS

    DO K2=1, NTENS

     DSTRESS(K2)=DSTRESS(K2)+DDSDDE(K2, K1)*DSTRAN(K1)

    END DO

    END DO

C

    DO K2=1, NTENS

    STRESS(K2)=STRESS(K2)+DSTRESS(K2)

    END DO

C

C RECOVER ELASTIC AND PLASTIC STRAINS FROM STATEV

C 

    DO K1=1, NTENS 

    EELAS(K1)=STATEV(K1)+DSTRAN(K1)

    EPLAS(K1)=STATEV(K1+NTENS) 

    END DO

    EQPLAS=STATEV(1+2*NTENS)

C 

C ----------------------------------------------------------- 

C Stress invarient calculation (HYDRO = -1*pressure, von Mises)

C ----------------------------------------------------------- 

C

    HYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE

C

    SMISES=(STRESS(1)-STRESS(2))*(STRESS(1)-STRESS(2)) + 

   1 (STRESS(2)-STRESS(3))*(STRESS(2)-STRESS(3)) + 

   2 (STRESS(3)-STRESS(1))*(STRESS(3)-STRESS(1)) 

    DO K1=NDI+1,NTENS 

    SMISES=SMISES+SIX*STRESS(K1)*STRESS(K1) 

    END DO

    SMISES=SQRT(SMISES/TWO)

C

C ----------------------------------------------------------- 

C Calculating hardening pressure resistance

C ----------------------------------------------------------- 

C

    HM=YHYDRO+HPAR*EQPLAS

    HP=-1.0*HM

C

C ----------------------------------------------------------- 

C ***********************************************************

C Testing yied criterion

C ***********************************************************

C ----------------------------------------------------------- 

C

C If the system is in extension: von MISES (HYDRO > 0)

    IF (HYDRO.GT.ZERO) THEN

    FFUN=SMISES-YIELDS

    IF (FFUN.GT.TOLER) THEN

    FLOW(1)=2.0*STRESS(1)-STRESS(2)-STRESS(3)

    FLOW(2)=2.0*STRESS(2)-STRESS(1)-STRESS(3)

    FLOW(3)=2.0*STRESS(3)-STRESS(2)-STRESS(1)

    DO K1=NDI+1,NTENS

    FLOW(K1)=6.0*STRESS(K1)

    END DO

C

    DO K1=1,NTENS

    FLOW(K1)=FLOW(K1)*HALF/SMISES

    END DO

C

    FFUN=SMISES-YIELDS

C

    VNEVEZ=0.0

    DO K1=1, NTENS

    DO K2=1, NTENS

     VNEVEZ=VNEVEZ+FLOW(K2)*DDSDDE(K2, K1)*FLOW(K1)

    END DO

    END DO 

C    

    DLAMB=FFUN/VNEVEZ

C

    DO K1=1, NTENS

    DEPLAS(K1)=DLAMB*FLOW(K1)

    EPLAS(K1)=EPLAS(K1)+DEPLAS(K1)

    EELAS(K1)=EELAS(K1)+DSTRAN(K1)-DEPLAS(K1)

    END DO

    EQPLAS=EQPLAS+(DEPLAS(1)+DEPLAS(2)+DEPLAS(3))

C

    DO K1=1, NTENS

    PSTRESS(K1)=0.0

    END DO

C

   DO K2=1, NTENS

    DO K1=1, NTENS

     PSTRESS(K2)=PSTRESS(K2)+DDSDDE(K2, K1)*DEPLAS(K1)

    END DO

    END DO

C

    DO K1=1, NTENS

    STRESS(K1)=STRESS(K1)-PSTRESS(K1)

    END DO

C

    DO K2=1, NTENS

    DO K4=1, NTENS

     DDSDDEEQ(K2,K4)= 0.0

    END DO

    END DO

C

    DO K1=1, NTENS

    DO K2=1, NTENS

     DO K3=1, NTENS

     DO K4=1, NTENS

      DDSDDEEQ(K2,K4)= DDSDDEEQ(K2,K4)+DDSDDE(K2, K1)*FLOW(K1)*

   1  FLOW(K3)*DDSDDE(K3, K4)

     END DO

     END DO

    END DO

    END DO

C    

    DO K2=1, NTENS

    DO K1=1, NTENS

     DDSDDE(K2,K1)=DDSDDE(K2,K1)-DDSDDEEQ(K2,K1)/VNEVEZ

    END DO

    END DO

C    

    ENDIF

    ELSE

C

C If the system is in compression: Cap (HYDRO =< 0)

    FFUN=((2*HYDRO-HP-HM)/(HP-HM))**2.0+(SMISES/YIELDS)**2.0-1.0

C

    IF (FFUN.GT.TOLER) THEN

C ----------------------------------------------------------- 

C *** *** *** *** *** *** *** *** *** *** *** *** *** *** ***

C Initiating inside teration - Newton-Rhapson

C *** *** *** *** *** *** *** *** *** *** *** *** *** *** ***

C ----------------------------------------------------------- 

C    

    DLAMB=0.0

    CNT=0

C    

    DO WHILE (FFUN.GT.1e-4)

C    

C ----------------------------------------------------------- 

C Determining the FLOW direction: df/dsig (first gradient)

C ----------------------------------------------------------- 

C    

    FLOW(1)=2.0*STRESS(1)-STRESS(2)-STRESS(3)

    FLOW(2)=2.0*STRESS(2)-STRESS(1)-STRESS(3)

    FLOW(3)=2.0*STRESS(3)-STRESS(2)-STRESS(1)

    DO K1=NDI+1,NTENS

    FLOW(K1)=6.0*STRESS(K1)

    END DO

C

    DO K1=1,NTENS

    FLOW(K1)=FLOW(K1)/YIELDS**2.0

    END DO

C    

    DO K1=1,NDI

    FLOW(K1)=FLOW(K1)+4.0/3.0*(TWO*HYDRO-HP-HM)/(HP-HM)**2.0

    END DO

C

C ----------------------------------------------------------- 

C Determining the FLOW2 curvature: d2f/dsig2 (second gradient)

C ----------------------------------------------------------- 

C

   DO K1=1, NTENS

    DO K2=1, NTENS

     FLOW2(K1,K2)=0.0

    END DO

    END DO

C

   DO K1=1, NDI

    DO K2=1, NDI

     FLOW2(K1,K2)=8.0/(9.0*(HP-HM)**2.0)-1.0/YIELDS**2.0

    END DO

    END DO

C

   DO K1=1, NDI

    FLOW2(K1,K1)=FLOW2(K1,K1)+3.0/YIELDS**2.0

   END DO

C

   DO K1=NDI+1,NTENS

    FLOW2(K1,K1)=6.0/YIELDS**2.0

   END DO

C

C ----------------------------------------------------------- 

C Linear equation system

C ----------------------------------------------------------- 

C

   DO K1=1, NTENS+1

    DO K2=1, NTENS+1

     EQSYS(K1,K2)=0.0

    END DO

    END DO

C 

    DO K1=1, NTENS

    EQSYS(K1,K1)=ONE

    END DO

C

    DO K1=1, NTENS

    DO K2=1, NTENS

     DO K3=1, NTENS

     EQSYS(K1,K3)=EQSYS(K1,K3)+DDSDDE(K1,K2)*FLOW2(K2,K3)*DLAMB

     END DO

    END DO

    END DO

C    

    DO K1=1, NTENS

    EQSYS(NTENS+1,K1)=FLOW(K1)

    END DO

C

    DO K1=1, NTENS

    DO K2=1, NTENS

     EQSYS(K1,NTENS+1)=EQSYS(K1,NTENS+1)+DDSDDE(K1,K2)*FLOW(K2)

    END DO

    END DO

C

C ----------------------------------------------------------- 

C Result vector

C ----------------------------------------------------------- 

C

    DO K1=1, NTENS

     PVECT(K1)=0.0

     BVECT(K1)=0.0

    END DO

C

    DO K1=1, NTENS

    PVECT(K1)=DSTRESS(K1)

    DO K2=1, NTENS

     PVECT(K1)=PVECT(K1)-DDSDDE(K1,K2)*DSTRAN(K2)+DDSDDE(K1,K2)*

   1 FLOW(K2)*DLAMB

    END DO

     BVECT(K1)=-PVECT(K1)

    END DO

    BVECT(NTENS+1)=-1.0*FFUN

C

    HMOD=-2.0*HPAR*HYDRO**2.0/HM**3.0*(FLOW(1)+FLOW(2)+FLOW(3))

    EQSYS(NTENS+1,NTENS+1)=HMOD

C    

C Solving EQ system

C

    CALL LUDCMP(EQSYS,NTENS+1,VINDX,D,CODE) 

    CALL LUBKSB(EQSYS,NTENS+1,VINDX,BVECT)

C

C Extracting plastic multiplier

C

    DLAMB=DLAMB+BVECT(NTENS+1)

C

C Calculating Stress step

C

    DO K1=1, NTENS

     DSTRESS(K1)=DSTRESS(K1)+BVECT(K1)

    END DO

C  

C Calculating stress state

C

    DO K1=1, NTENS

     STRESS(K1)=STRESS(K1)+BVECT(K1)

    END DO

C

C Recalculating stress invariants

C

    HYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE

    SMISES=(STRESS(1)-STRESS(2))*(STRESS(1)-STRESS(2)) + 

   1 (STRESS(2)-STRESS(3))*(STRESS(2)-STRESS(3)) + 

   2 (STRESS(3)-STRESS(1))*(STRESS(3)-STRESS(1)) 

    DO K1=NDI+1,NTENS 

    SMISES=SMISES+SIX*STRESS(K1)*STRESS(K1) 

    END DO

    SMISES=SQRT(SMISES/TWO)

C

C Recalculating plastic strain and hardening strain

C

    DO K1=1, NTENS

    DEPLAS(K1)=BVECT(NTENS+1)*FLOW(K1)

    EPLAS(K1)=EPLAS(K1)+DEPLAS(K1)

    EELAS(K1)=EELAS(K1)+DSTRAN(K1)-DEPLAS(K1)

    END DO

C

    EQPLAS=EQPLAS+(DEPLAS(1)+DEPLAS(2)+DEPLAS(3))

C

C Recalculating yield strength

C

    HM=YHYDRO+HPAR*EQPLAS

    HP=-1*HM

C    

C Recalculating yield function  

C    

    FFUN=((2.0*HYDRO-HP-HM)/(HP-HM))**2.0+(SMISES/YIELDS)**2.0-1.0

C

    CNT=CNT+1

    IF (CNT.GT.100) THEN

    CALL XIT

    ENDIF

C

    END DO

C ----------------------------------------------------------- 

C ---- END --- END --- End of RN --- END --- END --- END ----

C ----------------------------------------------------------- 

C

C ----------------------------------------------------------- 

C Recalculating FLOW and Hardening

C ----------------------------------------------------------- 

C    

    FLOW(1)=2.0*STRESS(1)-STRESS(2)-STRESS(3)

    FLOW(2)=2.0*STRESS(2)-STRESS(1)-STRESS(3)

    FLOW(3)=2.0*STRESS(3)-STRESS(2)-STRESS(1)

    DO K1=NDI+1,NTENS

    FLOW(K1)=6.0*STRESS(K1)

    END DO

C

    DO K1=1,NTENS

    FLOW(K1)=FLOW(K1)/YIELDS**2.0

    END DO

C    

    DO K1=1,NDI

    FLOW(K1)=FLOW(K1)+4.0/3.0*(2.0*HYDRO-HP-HM)/(HP-HM)**2.0

    END DO

C

    HMOD=-2.0*HPAR*HYDRO**2.0/HM**3.0*(FLOW(1)+FLOW(2)+FLOW(3))

C

C ----------------------------------------------------------- 

C Calculating tangent stiffness matrix

C ----------------------------------------------------------- 

C

    VNEVEZ=0.0

    DO K1=1, NTENS

    DO K2=1, NTENS

     VNEVEZ=VNEVEZ+FLOW(K2)*DDSDDE(K2, K1)*FLOW(K1)

    END DO

    END DO

C

    VNEVEZ=VNEVEZ-HMOD

C

    DO K1=1, NTENS

    PSTRESS(K1)=0.0

    END DO

C

    DO K2=1, NTENS

    DO K4=1, NTENS

     DDSDDEEQ(K2,K4)= 0.0

    END DO

    END DO

C

    DO K1=1, NTENS

    DO K2=1, NTENS

     DO K3=1, NTENS

     DO K4=1, NTENS

      DDSDDEEQ(K2,K4)= DDSDDEEQ(K2,K4)+DDSDDE(K2, K1)*FLOW(K1)*

   1  FLOW(K3)*DDSDDE(K3, K4)

     END DO

     END DO

    END DO

    END DO

C    

    DO K2=1, NTENS

    DO K1=1, NTENS

     DDSDDE(K2,K1)=DDSDDE(K2,K1)-DDSDDEEQ(K2,K1)/VNEVEZ

    END DO

    END DO

C

    ENDIF

    ENDIF

C

C ----------------------------------------------------------- 

C ******* Updating used variables ********

C ----------------------------------------------------------- 

C

    DO K1=1,NTENS 

    STATEV(K1)=EELAS(K1) 

    STATEV(K1+NTENS)=EPLAS(K1) 

    END DO

    STATEV(1+2*NTENS)=EQPLAS

C

    RETURN

    END

C

C ------------------------------------------------------------    

C ---------------------- END OF UMAT FILE --------------------

C ------------------------------------------------------------    

C    

    SUBROUTINE LUDCMP(A,N,VINDX,D,CODE)

C    

    INCLUDE 'ABA_PARAM.INC'

    CHARACTER*80 CMNAME

C

    PARAMETER(NMAX=100,TINY=1.5D-16)

    DIMENSION VV(N),A(N,N),VINDX(N)

C

    D=1

    CODE=0

    DO I=1,N

    VINDX(I)=0

    VV(I)=0.0

    END DO

C

    DO I=1,N

     AMAX=0.0

     DO J=1,N

      IF (ABS(A(I,J)).GT.AMAX) THEN 

      AMAX=ABS(A(I,J))

      END IF

     END DO

     IF(AMAX.LT.TINY) THEN

      CODE = 1

      RETURN

     END IF

     VV(I) = 1.0 / AMAX

    END DO

C

    DO J=1,N

     DO I=1,J-1

      SUMM = A(I,J)

      DO K=1,I-1

       SUMM = SUMM - A(I,K)*A(K,J) 

      END DO

      A(I,J) = SUMM

     END DO

     AMAX = 0.0

     DO I=J,N

      SUMM = A(I,J)

      DO K=1,J-1

       SUMM = SUMM - A(I,K)*A(K,J) 

      END DO

      A(I,J) = SUMM

      DUM = VV(I)*ABS(SUMM)

      IF(DUM.GE.AMAX) THEN

       VIMAX = I

       AMAX = DUM

      END IF

     END DO 

C  

     IF(J.NE.VIMAX) THEN

      DO K=1,N

       DUM = A(VIMAX,K)

       A(VIMAX,K) = A(J,K)

       A(J,K) = DUM

      END DO

      D = -D

      VV(VIMAX) = VV(J)

     END IF

C

     VINDX(J) = VIMAX

     IF(ABS(A(J,J)) < TINY) A(J,J) = TINY

C

     IF(J.NE.N) THEN

      DUM = 1.0 / A(J,J)

      DO I=J+1,N

       A(I,J) = A(I,J)*DUM

      END DO

     END IF 

    END DO

C

    RETURN

    END

C

C ******************************************************************

C * Solves the set of N linear equations A . X = B. Here A is  *

C * input, not as the matrix A but rather as its LU decomposition, *

C * determined by the routine LUDCMP. VINDX is input as the permuta-*

C * tion vector returned by LUDCMP. B is input as the right-hand *

C * side vector B, and returns with the solution vector X. A, N and*

C * VINDX are not modified by this routine and can be used for suc- *

C * cessive calls with different right-hand sides. This routine is *

C * also efficient for plain matrix inversion.          *

C ******************************************************************

    SUBROUTINE LUBKSB(A,N,VINDX,B)

C    

    INCLUDE 'ABA_PARAM.INC'

    CHARACTER*80 CMNAME

    DIMENSION A(N,N),VINDX(N),B(N)

C

    II = 0

C

    DO I=1,N

     LL = VINDX(I)

     SUMM = B(LL)

     B(LL) = B(I)

     IF(II.NE.0) THEN

      DO J=II,I-1

       SUMM = SUMM - A(I,J)*B(J)

      END DO

     ELSE IF(SUMM.NE.0.0) THEN

      II = I

     END IF

     B(I) = SUMM

    END DO

C

    DO I=N,1,-1

     SUMM = B(I)

     IF(I < N) THEN

      DO J=I+1,N

       SUMM = SUMM - A(I,J)*B(J)

      END DO

     END IF

     B(I) = SUMM / A(I,I)

    END DO

C

    RETURN

    END