原始文献：《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 propertiesC ----------------------------------------------------------- C PROPS(1) - Young's modulus C PROPS(2) - Poisson ratio C PROPS(3) - Yield (von Mises) stress at zero pressurC PROPS(4) - Initial yield hydro at compression (negative)C PROPS(5) - Hardening coefficientC ----------------------------------------------------------- CC Elastic propertiesC 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)CC ----------------------------------------------------------- C Elastic stiffness tensorC ----------------------------------------------------------- CC Lirear elasticC DO K1=1, NTENS DO K2=1, NTENS DDSDDE(K2, K1)=0.0 END DO END DOC 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 DOCC ----------------------------------------------------------- C Elastic predictionC ----------------------------------------------------------- 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 DOC DO K2=1, NTENS STRESS(K2)=STRESS(K2)+DSTRESS(K2) END DOCC RECOVER ELASTIC AND PLASTIC STRAINS FROM STATEVC 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))/THREEC 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)CC ----------------------------------------------------------- C Calculating hardening pressure resistanceC ----------------------------------------------------------- C HM=YHYDRO+HPAR*EQPLAS HP=-1.0*HMCC ----------------------------------------------------------- C ***********************************************************C Testing yied criterionC ***********************************************************C ----------------------------------------------------------- CC 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 DOC DO K1=1,NTENS FLOW(K1)=FLOW(K1)*HALF/SMISES END DOC FFUN=SMISES-YIELDSC 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/VNEVEZC 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 DOC DO K2=1, NTENS DO K1=1, NTENS PSTRESS(K2)=PSTRESS(K2)+DDSDDE(K2, K1)*DEPLAS(K1) END DO END DOC DO K1=1, NTENS STRESS(K1)=STRESS(K1)-PSTRESS(K1) END DOC DO K2=1, NTENS DO K4=1, NTENS DDSDDEEQ(K2,K4)= 0.0 END DO END DOC 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 DOC DO K2=1, NTENS DO K1=1, NTENS DDSDDE(K2,K1)=DDSDDE(K2,K1)-DDSDDEEQ(K2,K1)/VNEVEZ END DO END DOC ENDIF ELSECC If the system is in compression: Cap (HYDRO =< 0) FFUN=((2*HYDRO-HP-HM)/(HP-HM))**2.0+(SMISES/YIELDS)**2.0-1.0C IF (FFUN.GT.TOLER) THENC ----------------------------------------------------------- C *** *** *** *** *** *** *** *** *** *** *** *** *** *** ***C Initiating inside teration - Newton-RhapsonC *** *** *** *** *** *** *** *** *** *** *** *** *** *** ***C ----------------------------------------------------------- C DLAMB=0.0 CNT=0C 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 DOC DO K1=1,NTENS FLOW(K1)=FLOW(K1)/YIELDS**2.0 END DOC DO K1=1,NDI FLOW(K1)=FLOW(K1)+4.0/3.0*(TWO*HYDRO-HP-HM)/(HP-HM)**2.0 END DOCC ----------------------------------------------------------- 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 DOC 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 DOC DO K1=1, NDI FLOW2(K1,K1)=FLOW2(K1,K1)+3.0/YIELDS**2.0 END DOC DO K1=NDI+1,NTENS FLOW2(K1,K1)=6.0/YIELDS**2.0 END DOCC ----------------------------------------------------------- C Linear equation systemC ----------------------------------------------------------- C DO K1=1, NTENS+1 DO K2=1, NTENS+1 EQSYS(K1,K2)=0.0 END DO END DOC DO K1=1, NTENS EQSYS(K1,K1)=ONE END DOC 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 DOC DO K1=1, NTENS EQSYS(NTENS+1,K1)=FLOW(K1) END DOC DO K1=1, NTENS DO K2=1, NTENS EQSYS(K1,NTENS+1)=EQSYS(K1,NTENS+1)+DDSDDE(K1,K2)*FLOW(K2) END DO END DOCC ----------------------------------------------------------- C Result vectorC ----------------------------------------------------------- C DO K1=1, NTENS PVECT(K1)=0.0 BVECT(K1)=0.0 END DOC 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*FFUNC HMOD=-2.0*HPAR*HYDRO**2.0/HM**3.0*(FLOW(1)+FLOW(2)+FLOW(3)) EQSYS(NTENS+1,NTENS+1)=HMODC C Solving EQ systemC CALL LUDCMP(EQSYS,NTENS+1,VINDX,D,CODE) CALL LUBKSB(EQSYS,NTENS+1,VINDX,BVECT)CC Extracting plastic multiplierC DLAMB=DLAMB+BVECT(NTENS+1)CC Calculating Stress stepC DO K1=1, NTENS DSTRESS(K1)=DSTRESS(K1)+BVECT(K1) END DOC C Calculating stress stateC DO K1=1, NTENS STRESS(K1)=STRESS(K1)+BVECT(K1) END DOCC Recalculating stress invariantsC 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)CC Recalculating plastic strain and hardening strainC 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 DOC EQPLAS=EQPLAS+(DEPLAS(1)+DEPLAS(2)+DEPLAS(3))CC Recalculating yield strengthC HM=YHYDRO+HPAR*EQPLAS HP=-1*HMC C Recalculating yield function C FFUN=((2.0*HYDRO-HP-HM)/(HP-HM))**2.0+(SMISES/YIELDS)**2.0-1.0C CNT=CNT+1 IF (CNT.GT.100) THEN CALL XIT ENDIFC END DOC ----------------------------------------------------------- C ---- END --- END --- End of RN --- END --- END --- END ----C ----------------------------------------------------------- CC ----------------------------------------------------------- C Recalculating FLOW and HardeningC ----------------------------------------------------------- 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 DOC DO K1=1,NTENS FLOW(K1)=FLOW(K1)/YIELDS**2.0 END DOC DO K1=1,NDI FLOW(K1)=FLOW(K1)+4.0/3.0*(2.0*HYDRO-HP-HM)/(HP-HM)**2.0 END DOC HMOD=-2.0*HPAR*HYDRO**2.0/HM**3.0*(FLOW(1)+FLOW(2)+FLOW(3))CC ----------------------------------------------------------- C Calculating tangent stiffness matrixC ----------------------------------------------------------- C VNEVEZ=0.0 DO K1=1, NTENS DO K2=1, NTENS VNEVEZ=VNEVEZ+FLOW(K2)*DDSDDE(K2, K1)*FLOW(K1) END DO END DOC VNEVEZ=VNEVEZ-HMODC DO K1=1, NTENS PSTRESS(K1)=0.0 END DOC DO K2=1, NTENS DO K4=1, NTENS DDSDDEEQ(K2,K4)= 0.0 END DO END DOC 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 DOC DO K2=1, NTENS DO K1=1, NTENS DDSDDE(K2,K1)=DDSDDE(K2,K1)-DDSDDEEQ(K2,K1)/VNEVEZ END DO END DOC ENDIF ENDIFCC ----------------------------------------------------------- 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)=EQPLASC RETURN ENDCC ------------------------------------------------------------ C ---------------------- END OF UMAT FILE --------------------C ------------------------------------------------------------ C SUBROUTINE LUDCMP(A,N,VINDX,D,CODE)C INCLUDE 'ABA_PARAM.INC' CHARACTER*80 CMNAMEC 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 DOC 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 DOC 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 IFC VINDX(J) = VIMAX IF(ABS(A(J,J)) < TINY) A(J,J) = TINYC 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 DOC RETURN ENDCC ******************************************************************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 = 0C 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 DOC 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 DOC RETURN END来源：我的博士日记
