材料本构模型设计方案，拉伸变形采用von Mises屈服，压缩侧 cap屈服本构模型设计。_ACT

材料本构模型设计方案，拉伸变形采用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)

INCLUDE 'ABA_PARAM.INC'

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)

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 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 Elastic properties

EMOD=PROPS(1)

ENU=PROPS(2)

YIELDS=PROPS(3)

YHYDRO=PROPS(4)

HPAR=PROPS(5)

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

EG2=EG*TWO

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

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

C Elastic stiffness tensor

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

C Lirear elastic

DO K1=1, NTENS

DO K2=1, NTENS

DDSDDE(K2, K1)=0.0

END DO

DO K1=1, NDI

DO K2=1, NDI

DDSDDE(K2, K1)=ELAM

END DO

DDSDDE(K1, K1)=EG2+ELAM

END DO

DO K1=NDI+1, NTENS

DDSDDE(K1, K1)=EG

END DO

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

C Elastic prediction

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

DO K2=1, NTENS

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

END DO

C RECOVER ELASTIC AND PLASTIC STRAINS FROM STATEV

DO K1=1, NTENS

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

EPLAS(K1)=STATEV(K1+NTENS)

END DO

EQPLAS=STATEV(1+2*NTENS)

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

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

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 Calculating hardening pressure resistance

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

HM=YHYDRO+HPAR*EQPLAS

HP=-1.0*HM

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

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

C Testing yied criterion

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

DO K1=1,NTENS

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

END DO

FFUN=SMISES-YIELDS

VNEVEZ=0.0

DO K1=1, NTENS

DO K2=1, NTENS

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

END DO

DLAMB=FFUN/VNEVEZ

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))

DO K1=1, NTENS

PSTRESS(K1)=0.0

END DO

DO K2=1, NTENS

DO K1=1, NTENS

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

END DO

DO K1=1, NTENS

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

END DO

DO K2=1, NTENS

DO K4=1, NTENS

DDSDDEEQ(K2,K4)= 0.0

END DO

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

DO K2=1, NTENS

DO K1=1, NTENS

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

END DO

ENDIF

ELSE

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

IF (FFUN.GT.TOLER) THEN

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

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

C Initiating inside teration - Newton-Rhapson

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

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

DLAMB=0.0

CNT=0

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

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

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

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

DO K1=1,NTENS

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

END DO

DO K1=1,NDI

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

END DO

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

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

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

DO K1=1, NTENS

DO K2=1, NTENS

FLOW2(K1,K2)=0.0

END DO

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

DO K1=1, NDI

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

END DO

DO K1=NDI+1,NTENS

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

END DO

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

C Linear equation system

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

DO K1=1, NTENS+1

DO K2=1, NTENS+1

EQSYS(K1,K2)=0.0

END DO

DO K1=1, NTENS

EQSYS(K1,K1)=ONE

END DO

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

DO K1=1, NTENS

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

END DO

DO K1=1, NTENS

DO K2=1, NTENS

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

END DO

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

C Result vector

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

DO K1=1, NTENS

PVECT(K1)=0.0

BVECT(K1)=0.0

END DO

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

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

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

C Solving EQ system

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

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

C Extracting plastic multiplier

DLAMB=DLAMB+BVECT(NTENS+1)

C Calculating Stress step

DO K1=1, NTENS

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

END DO

C Calculating stress state

DO K1=1, NTENS

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

END DO

C Recalculating stress invariants

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 Recalculating plastic strain and hardening strain

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

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

C Recalculating yield strength

HM=YHYDRO+HPAR*EQPLAS

HP=-1*HM

C Recalculating yield function

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

CNT=CNT+1

IF (CNT.GT.100) THEN

CALL XIT

ENDIF

END DO

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

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

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

C Recalculating FLOW and Hardening

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

DO K1=1,NTENS

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

END DO

DO K1=1,NDI

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

END DO

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

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

C Calculating tangent stiffness matrix

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

VNEVEZ=0.0

DO K1=1, NTENS

DO K2=1, NTENS

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

END DO

VNEVEZ=VNEVEZ-HMOD

DO K1=1, NTENS

PSTRESS(K1)=0.0

END DO

DO K2=1, NTENS

DO K4=1, NTENS

DDSDDEEQ(K2,K4)= 0.0

END DO

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

DO K2=1, NTENS

DO K1=1, NTENS

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

END DO

ENDIF

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

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

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

DO K1=1,NTENS

STATEV(K1)=EELAS(K1)

STATEV(K1+NTENS)=EPLAS(K1)

END DO

STATEV(1+2*NTENS)=EQPLAS

RETURN

END

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

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

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

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

INCLUDE 'ABA_PARAM.INC'

CHARACTER*80 CMNAME

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

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

D=1

CODE=0

DO I=1,N

VINDX(I)=0

VV(I)=0.0

END DO

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

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

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

VINDX(J) = VIMAX

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

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

RETURN

END

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)

INCLUDE 'ABA_PARAM.INC'

CHARACTER*80 CMNAME

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

II = 0

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

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

RETURN

END

来源：我的博士日记

HILL48 +各向同性Voce硬化umat子程序

!=====================================================================! HILL48 PLASTICITY UMAT WITH VOCE ISOTROPIC HARDENING! IMPLEMENTATION: GENERALIZED RADIAL-RETURN IN EIGENSPACE!===============================================================================!! AUTHOR: ! Mohammad Hasaninia! Computational Advanced Manufacturing and Materials Laboratory (CAMML)! Department of Mechanical Engineering! University of Wyoming!! REFERENCE:! Versino, D. and Bennett, K.C. (2018). Generalized radial-return mapping ! algorithm for anisotropic von Mises plasticity framed in material ! eigenspace. Int. J. Numer. Meth. Engng, 116: 202-222.! https://doi.org/10.1002/nme.5921!!-------------------------------------------------------------------------------! INPUT PROPERTIES (PROPS) DEFINITION:! PROPS(1) : Young's Modulus (E)! PROPS(2) : Poisson's Ratio (nu)! PROPS(3) : Initial Yield Stress (Sigma_Y0)! PROPS(4) : Saturation Stress (Q_inf) - Voce Parameter! PROPS(5) : Hardening Rate (b) - Voce Parameter! PROPS(6) : Hill F! PROPS(7) : Hill G! PROPS(8) : Hill H! PROPS(9) : Hill L! PROPS(10) : Hill M! PROPS(11) : Hill N!! STATE VARIABLES (STATEV) OUTPUT:! STATEV(1) : Equivalent Plastic Strain (alpha)! STATEV(2) : Current Yield Stress! STATEV(3-8) : Plastic Strain Tensor (11, 22, 33, 12, 13, 23)! STATEV(9) : Convergence Flag (0=Converged, 1=Failed)!=====================================================================SUBROUTINE UMAT(STRESS, STATEV, DDSDDE, SSE, SPD, SCD, RPL, & DDSDDT, DRPLDE, DRPLDT, STRAN, DSTRAN, TIME, DTIME, & TEMP, DTEMP, PREDEF, DPRED, CMNAME, NDI, NSHR, NTENS, & NSTATV, PROPS, NPROPS, COORDS, DROT, PNEWDT, & CELENT, DFGRD0, DFGRD1, NOEL, NPT, LAYER, KSPT, & KSTEP, KINC) IMPLICIT NONE ! Abaqus UMAT interface variables CHARACTER*80 CMNAME INTEGER NDI, NSHR, NTENS, NSTATV, NPROPS INTEGER NOEL, NPT, LAYER, KSPT, KSTEP, KINC REAL(8) STRESS(NTENS), STATEV(NSTATV), DDSDDE(NTENS,NTENS) REAL(8) SSE, SPD, SCD, RPL, DRPLDT, DTIME, TEMP, DTEMP, PNEWDT, CELENT REAL(8) STRAN(NTENS), DSTRAN(NTENS), TIME(2), PREDEF(1), DPRED(1) REAL(8) PROPS(NPROPS), COORDS(3), DROT(3,3), DFGRD0(3,3), DFGRD1(3,3) REAL(8) DDSDDT(NTENS), DRPLDE(NTENS) ! Local variables for plasticity computation REAL(8) :: alpha_tr, alpha, StressYield, s_y, s_y_D, s_y_tr, f_gamma_tr REAL(8) :: Sigma_Y0, Q_inf, b_voce REAL(8) :: MatM(6,6), MatD(6,6), MatI(6,6) REAL(8) :: MatM_copy(6,6), Mat_inv(6,6) REAL(8) :: W(6), Q(6,6), Gamma(6,6), K(6,6), dK_dDelta_gamma_matrix(6,6) REAL(8) :: s_tr(6), strain_p(6), Stress_New(6), s_tilde_tr(6) REAL(8) :: Vec(6), dev_corrected(6) REAL(8) :: E, XNUE, EBULK3, EG2, EG, ELAM REAL(8) :: HillF, HillG, HillH, HillL, HillM, HillN REAL(8) :: pressure, d_ga, omega, d_al, f_gamma_delta, delta_alpha REAL(8) :: Df_gamma_d_d_gamma, d_omega_d_d_gamma, d_sy_d_d_gamma REAL(8) :: TOLERANCE, exp_term INTEGER :: iteration, K3, K5, K6, K7, K8, INFO LOGICAL :: converged ! Extract material properties E = PROPS(1) ! Young's modulus XNUE = PROPS(2) ! Poisson's ratio Sigma_Y0 = PROPS(3) ! Initial yield stress Q_inf = PROPS(4) ! Saturation stress (Q∞) b_voce = PROPS(5) ! Voce hardening rate parameter HillF = PROPS(6) ! Hill parameter F HillG = PROPS(7) ! Hill parameter G HillH = PROPS(8) ! Hill parameter H HillL = PROPS(9) ! Hill parameter L HillM = PROPS(10) ! Hill parameter M HillN = PROPS(11) ! Hill parameter N ! Initialize state variables IF (TIME(2) < 1.0D-6) THEN alpha_tr = 0.0D0 ! Voce: σ_y0 = σ_Y0 + Q∞(1 - e^0) = σ_Y0 StressYield = Sigma_Y0 strain_p = 0.0D0 STATEV(1:9) = 0.0D0 ELSE alpha_tr = STATEV(1) StressYield = STATEV(2) strain_p(1:6) = STATEV(3:8) END IF ! Create elasticity matrix EBULK3 = E/(1.0D0-2.0D0*XNUE) EG2 = E/(1.0D0+XNUE) EG = EG2/2.0D0 ELAM = (EBULK3-EG2)/3.0D0 MatD = 0.0D0 MatD(1:3,1:3) = ELAM DO K3 = 1, 3 MatD(K3,K3) = EG2 + ELAM END DO DO K3 = 4, 6 MatD(K3, K3) = EG END DO ! Create Hill anisotropy matrix MatM = 0.0D0 MatM(1,1) = HillG + HillH MatM(2,2) = HillF + HillH MatM(3,3) = HillF + HillG MatM(1,2) = -HillH MatM(1,3) = -HillG MatM(2,3) = -HillF MatM(2,1) = MatM(1,2) MatM(3,1) = MatM(1,3) MatM(3,2) = MatM(2,3) MatM(4,4) = 2.0D0*HillN MatM(5,5) = 2.0D0*HillM MatM(6,6) = 2.0D0*HillL ! Identity matrix MatI = 0.0D0 DO K5 = 1, 6 MatI(K5,K5) = 1.0D0 END DO ! Compute trial stress (elastic predictor) Vec = MATMUL(MatD, DSTRAN) Stress_New = Vec + STRESS s_tr = Stress_New ! Extract deviatoric part pressure = (Stress_New(1)+Stress_New(2)+Stress_New(3))/3.0D0 s_tr(1:3) = s_tr(1:3) - pressure ! Eigendecomposition of anisotropy tensor MatM_copy = MatM CALL compute_eigendecomposition(MatM_copy, W, Q, Gamma, INFO) IF (INFO /= 0) THEN WRITE(*,*) 'Error in eigendecomposition, INFO = ', INFO RETURN END IF ! Transform trial stress to eigenspace s_tilde_tr = MATMUL(TRANSPOSE(Q), s_tr) ! Initialize trial values alpha_tr = STATEV(1) s_y_tr = StressYield ! Compute trial yield function f_gamma_tr = 0.5D0 * DOT_PRODUCT(s_tilde_tr, MATMUL(Gamma, s_tilde_tr)) - s_y_tr**2 ! Check yield condition IF (f_gamma_tr <= 0.0D0) THEN ! Elastic step STRESS = Stress_New d_ga = 0.0D0 alpha = alpha_tr s_y = StressYield ELSE ! Plastic step - Newton-Raphson iteration tolerance = 5.0D-4 d_ga = 0.0D0 alpha = alpha_tr converged = .FALSE. DO iteration = 1, 250 ! Build K matrix (diagonal) K = 0.0D0 DO K6 = 1, 6 K(K6,K6) = Gamma(K6,K6) / (1.0D0 + 2.0D0*EG*d_ga*Gamma(K6,K6))**2 END DO ! Compute auxiliary variable omega omega = SQRT(0.5D0 * DOT_PRODUCT(s_tilde_tr, MATMUL(K, s_tilde_tr))) IF (omega <= 1.0D-14) THEN WRITE(*,*) 'Warning: omega is too s mall' EXIT END IF ! Update alpha delta_alpha = 2.0D0 * d_ga * omega alpha = alpha_tr + delta_alpha ! ======================================================== ! VOCE HARDENING MODEL ! σ_y = σ_Y0 + Q∞(1 - exp(-b·α)) ! ======================================================== exp_term = EXP(-b_voce * alpha) s_y = Sigma_Y0 + Q_inf * (1.0D0 - exp_term) ! Derivative: ∂σ_y/∂α = Q∞·b·exp(-b·α) s_y_D = Q_inf * b_voce * exp_term ! ======================================================== ! Compute residual function f_gamma_delta = (s_y / omega) - 1.0D0 ! Check convergence IF (ABS(f_gamma_delta) <= tolerance) THEN StressYield = s_y converged = .TRUE. EXIT END IF ! Compute derivatives for Newton-Raphson dK_dDelta_gamma_matrix = 0.0D0 DO K7 = 1, 6 dK_dDelta_gamma_matrix(K7,K7) = -4.0D0*EG*Gamma(K7,K7)**2 & / (1.0D0 + 2.0D0*EG*d_ga*Gamma(K7,K7))**3 END DO d_omega_d_d_gamma = (1.0D0 / (4.0D0 * omega)) & * DOT_PRODUCT(s_tilde_tr, MATMUL(dK_dDelta_gamma_matrix, s_tilde_tr)) ! Chain rule: ∂σ_y/∂Δγ = (∂σ_y/∂α)·(∂α/∂Δγ) d_sy_d_d_gamma = 2.0D0 * s_y_D * (omega + d_ga * d_omega_d_d_gamma) ! Derivative of residual Df_gamma_d_d_gamma = (1.0D0 / omega) * d_sy_d_d_gamma & - (s_y / (omega**2)) * d_omega_d_d_gamma ! Newton-Raphson update d_ga = d_ga - f_gamma_delta / Df_gamma_d_d_gamma END DO ! Update stress after Newton loop Mat_inv = 0.0D0 DO K8 = 1, 6 Mat_inv(K8,K8) = 1.0D0 / (1.0D0 + 2.0D0*EG*d_ga*Gamma(K8,K8)) END DO dev_corrected = MATMUL(Q, MATMUL(Mat_inv, MATMUL(TRANSPOSE(Q), s_tr))) ! Assemble full Cauchy stress STRESS(1:3) = dev_corrected(1:3) + pressure STRESS(4:6) = dev_corrected(4:6) IF (.NOT. converged) THEN WRITE(*,*) 'WARNING: Newton-Raphson did not converge!' StressYield = s_y STATEV(9) = 1.0D0 END IF END IF ! Consistent tangent operator DDSDDE = MatD ! Update state variables STATEV(1) = alpha ! Equivalent plastic strain STATEV(2) = s_y ! Current yield stress strain_p = strain_p + d_ga * MATMUL(MatM, dev_corrected) STATEV(3:8) = strain_p(1:6) ! Plastic strain components RETURN CONTAINS !********************************************************************** SUBROUTINE compute_eigendecomposition(A, W, Q, Gamma, INFO) IMPLICIT NONE INTEGER, INTENT(OUT) :: INFO REAL(8), DIMENSION(6,6), INTENT(INOUT) :: A REAL(8), DIMENSION(6), INTENT(OUT) :: W REAL(8), DIMENSION(6,6), INTENT(OUT) :: Q REAL(8), DIMENSION(6,6), INTENT(OUT) :: Gamma INTEGER :: LWORK, i, j REAL(8), DIMENSION(:), ALLOCATABLE :: WORK LOGICAL :: is_symmetric ! Check symmetry is_symmetric = .TRUE. DO i = 1, 6 DO j = i + 1, 6 IF (ABS(A(i,j) - A(j,i)) > 1.0D-12) THEN is_symmetric = .FALSE. EXIT END IF END DO IF (.NOT. is_symmetric) EXIT END DO IF (.NOT. is_symmetric) THEN WRITE(*,*) "Error: Matrix not symmetric" INFO = -1 RETURN END IF ! LAPACK eigendecomposition LWORK = -1 ALLOCATE(WORK(1)) CALL DSYEV('V', 'U', 6, A, 6, W, WORK, LWORK, INFO) LWORK = INT(WORK(1)) DEALLOCATE(WORK) ALLOCATE(WORK(LWORK)) CALL DSYEV('V', 'U', 6, A, 6, W, WORK, LWORK, INFO) IF (INFO == 0) THEN Q = A Gamma = 0.0D0 DO i = 1, 6 Gamma(i,i) = W(i) END DO ELSE WRITE(*,*) "Eigendecomposition failed, INFO:", INFO Gamma = 0.0D0 END IF DEALLOCATE(WORK) END SUBROUTINE compute_eigendecompositionEND SUBROUTINE UMAT!*****************************BOTTOM*****************************原始链接：https://github.com/mhasaninia/Hill48-Voce-UMAT/tree/main来源：我的博士日记