这是对编写的UMAT用于三维实体单元的一个验证。UMAT虽然是从一维Johnson-Cook模型中建立起来,但是在UMAT中率相关塑性公式将它扩展到了三维情况,所以能用于三维实体单元也是意料之中。
虽然取得的结果是一致的,但是计算时间上却有很大差别,在个人计算机上,使用轴对称二维网格完成一次计算只需要10分钟,然后完成一次三维网格的计算需要3个小时。相比而言使用轴对称二维模型要经济得多。
18.5 UMAT的Fortran程序
18.5.1 UMAT
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 MATERL
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), 4 DFGRD0(3,3),DFGRD1(3,3) C
DIMENSION EELAS(6),EPLAS(6),FLOW(6)
PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0, HALF =0.5d0)
DATA NEWTON,TOLER/40,1.D-6/ C
C ----------------------------------------------------------- C UMAT FOR JOHNSON-COOK MODEL
C ----------------------------------------------------------- C PROPS(1) - YANG'S MODULUS C PROPS(2) - POISSON RATIO
C PROPS(3) - INELASTIC HEAT FRACTION
19-31
C PARAMETERS OF JOHNSON-COOK MODEL: C PROPS(4) - A C PROPS(5) - B C PROPS(6) - n C PROPS(7) - C C PROPS(8) - m
C ----------------------------------------------------------- C
IF (NDI.NE.3) THEN WRITE(6,1)
1 FORMAT(//,30X,'***ERROR - THIS UMAT MAY ONLY BE USED FOR ',
1 'ELEMENTS WITH THREE DIRECT STRESS COMPONENTS') ENDIF C
C ELASTIC PROPERTIES C
EMOD=PROPS(1) ENU=PROPS(2)
IF(ENU.GT.0.4999.AND.ENU.LT.0.5001) ENU=0.499 EBULK3=EMOD/(ONE-TWO*ENU) EG2=EMOD/(ONE+ENU) EG=EG2/TWO EG3=THREE*EG
ELAM=(EBULK3-EG2)/THREE C
C ELASTIC STIFFNESS C
DO 20 K1=1,NTENS DO 10 K2=1,NTENS
DDSDDE(K2,K1)=0.0 10 CONTINUE 20 CONTINUE C
DO 40 K1=1,NDI DO 30 K2=1,NDI
DDSDDE(K2,K1)=ELAM 30 CONTINUE
DDSDDE(K1,K1)=EG2+ELAM 40 CONTINUE
DO 50 K1=NDI+1,NTENS DDSDDE(K1,K1)=EG
19-32
50 CONTINUE C
C CALCULATE STRESS FROM ELASTIC STRAINS C
DO 70 K1=1,NTENS DO 60 K2=1,NTENS
STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1) 60 CONTINUE 70 CONTINUE C
C RECOVER ELASTIC AND PLASTIC STRAINS C
DO 80 K1=1,NTENS
EELAS(K1)=STATEV(K1)+DSTRAN(K1) EPLAS(K1)=STATEV(K1+NTENS) 80 CONTINUE
EQPLAS=STATEV(1+2*NTENS) C
C CALCULATE MISES STRESS C
IF(NPROPS.GT.5.AND.PROPS(4).GT.0.0) THEN
SMISES=(STRESS(1)-STRESS(2))*(STRESS(1)-STRESS(2)) + 1 (STRESS(2)-STRESS(3))*(STRESS(2)-STRESS(3)) + 1 (STRESS(3)-STRESS(1))*(STRESS(3)-STRESS(1)) DO 90 K1=NDI+1,NTENS
SMISES=SMISES+SIX*STRESS(K1)*STRESS(K1) 90 CONTINUE
SMISES=SQRT(SMISES/TWO) C
C CALL USERHARD SUBROUTINE, GET HARDENING RATE AND YIELD STRESS C C
CALL USERHARD(SYIEL0,HARD,EQPLAS,PROPS(4)) C DETERMINE IF ACTIVELY YIELDING C
IF (SMISES.GT.(1.0+TOLER)*SYIEL0) THEN C
C MATERIAL RESPONSE IS PLASTIC, DETERMINE FLOW DIRECTION C
SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE ONESY=ONE/SMISES DO 110 K1=1,NDI
19-33
FLOW(K1)=ONESY*(STRESS(K1)-SHYDRO) 110 CONTINUE
DO 120 K1=NDI+1,NTENS
FLOW(K1)=STRESS(K1)*ONESY 120 CONTINUE C
C READ PARAMETERS OF JOHNSON-COOK MODEL C
A=PROPS(4) B=PROPS(5) EN=PROPS(6) C=PROPS(7) EM=PROPS(8) C
C NEWTON ITERATION C
SYIELD=SYIEL0
DEQPL=(SMISES-SYIELD)/EG3 DSTRES=TOLER*SYIEL0/EG3
DEQMIN=HALF*DTIME*EXP(1.0D-4/C) DO 130 KEWTON=1,NEWTON
DEQPL=MAX(DEQPL,DEQMIN)
CALL USERHARD(SYIELD,HARD,EQPLAS+DEQPL,PROPS(4)) TVP=C*LOG(DEQPL/DTIME) TVP1=TVP+ONE
HARD1=HARD*TVP1+SYIELD*C/DEQPL SYIELD=SYIELD*TVP1
RHS=SMISES-EG3*DEQPL-SYIELD DEQPL=DEQPL+RHS/(EG3+HARD1)
IF(ABS(RHS/EG3) .LE. DSTRES ) GOTO 140 130 CONTINUE
WRITE(6,2) NEWTON
2 FORMAT(//,30X,'***WARNING - PLASTICITY ALGORITHM DID NOT ',
1 'CONVERGE AFTER ',I3,' ITERATIONS') 140 CONTINUE
EFFHRD=EG3*HARD1/(EG3+HARD1) C
C CALCULATE STRESS AND UPDATE STRAINS C
DO 150 K1=1,NDI
STRESS(K1)=FLOW(K1)*SYIELD+SHYDRO
19-34
EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL/TWO EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL/TWO 150 CONTINUE
DO 160 K1=NDI+1,NTENS
STRESS(K1)=FLOW(K1)*SYIELD
EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL 160 CONTINUE
EQPLAS=EQPLAS+DEQPL
SPD=DEQPL*(SYIEL0+SYIELD)/TWO RPL = PROPS(3)*SPD/DTIME C
C JACOBIAN C
EFFG=EG*SYIELD/SMISES EFFG2=TWO*EFFG
EFFG3=THREE*EFFG2/TWO
EFFLAM=(EBULK3-EFFG2)/THREE DO 220 K1=1,NDI DO 210 K2=1,NDI
DDSDDE(K2,K1)=EFFLAM 210 CONTINUE
DDSDDE(K1,K1)=EFFG2+EFFLAM 220 CONTINUE
DO 230 K1=NDI+1,NTENS DDSDDE(K1,K1)=EFFG 230 CONTINUE
DO 250 K1=1,NTENS DO 240 K2=1,NTENS
DDSDDE(K2,K1)=DDSDDE(K2,K1)+FLOW(K2)*FLOW(K1) 1 *(EFFHRD-EFFG3) 240 CONTINUE 250 CONTINUE ENDIF ENDIF C
C STORE STRAINS IN STATE VARIABLE ARRAY C
DO 310 K1=1,NTENS
STATEV(K1)=EELAS(K1)
STATEV(K1+NTENS)=EPLAS(K1) 310 CONTINUE
19-35