产品
Abaqus2021的接口程序；Johnson-Cook本构
看了很多本构模型的代码，也走了很多弯路，“error”千奇百怪
最后总算能跑完程序了
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,JSTEP,KINC)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 CMNAME
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),DFGRD0(3,3),DFGRD1(3,3),
4 JSTEP(4)
C 数值及矩阵尺寸声明
PARAMETER(ZERO=0.0D0,ONE=0.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0)
DIMENSION EELAS(6),EPLAS(6),FLOW(6)
DATA NEWTON,TOLER/20,1.D-6/
C 读取材料参数
E=PROPS(1)
ENU=PROPS(2)
A=PROPS(3)
B=PROPS(4)
EN=PROPS(5)
C 计算材料参数
EBULK3=E/(ONE-TWO*ENU)
EG2=E/(ONE+ENU)
EG=EG2/TWO
EG3=THREE*EG
ELAM=(EBULK3-EG2)/THREE
C 计算材料弹性刚度矩阵DDSDDE
DO 20 K1=1,NTENS
DO 10 K2=1,NTENS
DDSDDE(K2,K1)=0.0
10 CONTINUE
20 CONTINUE
DO 40 K1=1,NDI
DO 30 K2=1,NDI
DDSDDE(K2,K1)=ELAM
30 CONTINUE
DDSDDE(K1,K1)=ELAM+EG2
40 CONTINUE
DO 50 K1=NDI+1,NTENS
DDSDDE(K1,K1)=EG
50 CONTINUE
C 读取等效塑性应变、弹性应变、塑性应变
EQPLAS=STATEV(1+2*NTENS)
DO 60 K1=1,NTENS
EELAS(K1)=STATEV(K1)+DSTRAN(K1)
EPLAS(K1)=STATEV(K1+NTENS)
60 CONTINUE
C 计算试应力
DO 80 K1=1,NTENS
DO 70 K2=1,NTENS
STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)
70 CONTINUE
80 CONTINUE
IF(NPROPS.GT.2.AND.PROPS(3).GT.0.0) THEN
C 计算米塞斯等效试应力SMISES
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))
SMISES=SMISES+SIX*(STRESS(4)*STRESS(4)+STRESS(5)*STRESS(5)
1 +STRESS(6)*STRESS(6))
SMISES=SQRT(SMISES/TWO)
C 计算硬化率和当前屈服应力
CALL USERHARD(SYIEL0,HARD,EQPLAS,PROPS,NPROPS)
C 判断材料是否处于屈服阶段
IF (SMISES.GT.(1.0+TOLER)*SYIEL0) THEN
C 计算塑性流动方向
C 计算静水压力SHYDRO
SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE
ONESY=ONE/SMISES
DO 90 K1=1,NDI
FLOW(K1)=(STRESS(K1)-SHYDRO)*ONESY
90 CONTINUE
DO 100 K1=NDI+1,NTENS
FLOW(K1)=STRESS(K1)*ONESY
100 CONTINUE
C 牛顿迭代
SYIELD=SYIEL0
DEQPL=(SMISES-SYIELD)/EG3
DO 110 KEWTON=1,NEWTON
CALL USERHARD(SYIELD,HARD,EQPLAS+DEQPL,PROPS,NPROPS)
RHS=SMISES-EG3*DEQPL-SYIELD
DEQPL=DEQPL+RHS/(EG3+HARD)
IF (ABS(RHS).LE.TOLER*SYIEL0) GOTO 10
110 CONTINUE
C 将警告消息写入.dat文件
WRITE(6,2) NEWTON
2 FORMAT(//,30X,'***WARNING - PLASTICITY ALGORITHM DID NOT ',
1 'CONVERGE AFTER ',I3,' ITERATIONS')
120 CONTINUE
C 计算此增量步硬化段斜率
EFFHRD=EG3*HARD/(EG3+HARD)
C 计算应力和更新应变
DO 130 K1=1,NDI
STRESS(K1)=FLOW(K1)*SYIELD+SHYDRO
EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL/TWO
EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL/TWO
130 CONTINUE
DO 140 K1=DNI+1,NTENS
STRESS(K1)=FLOW(K1)*SYIELD
EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL
EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL
140 CONTINUE
EQPLAS=EQPLAS+DEQPL
C 计算塑性变形下的雅克比矩阵
EFFG=EG*SYIELD/SMISES
EFFG2=TWO*EFFG
EFFG3=THREE*EFFG2/TWO
EFFLAM=(EBULK3-EFFG2)/THREE
DO 160 K1=1,NDI
DO 150 K2=1,NDI
DDSDDE(K2,K1)=EFFLAM
150 CONTINUE
DDSDDE(K1,K1)=EFFG2+EFFLAM
160 CONTINUE
DO 170 K1=NDI+1,NTENS
DDSDDE(K1,K1)=EFFG
170 CONTINUE
DO 190 K1=1,NTENS
DO 180 K2=1,NTENS
DDSDDE(K2,K1)=DDSDDE(K2,K1)+FLOW(K2)*FLOW(K1)
1 *(EFFHRD-EFFG3)
180 CONTINUE
190 CONTINUE
END IF
END IF
C 将弹性应变，塑性应变分量和等效塑性应变保存到状态变量中，并
C 传递到下一个增量步
DO 200 K1=1,NTENS
STATEV(K1)=EELAS(K1)
STATEV(K1+NTENS)=EPLAS(K1)
200 CONTINUE
STATEV(1+2*NTENS)=EQPLAS
STATEV(2+2*NTENS)=SYIELD
STATEV(3+2*NTENS)=DEQPL
STATEV(4+2*NTENS)=HARD
STATEV(5+2*NTENS)=EQPLAS**EN
STATEV(6+2*NTENS)=B*(EQPLAS**(EN-1))
RETURN
END
SUBROUTINE USERHARD(SYIELD,HARD,EQPLAS,PROPS,NPROPS)
C
INCLUDE 'ABA_PARAM.INC'
DIMENSION PROPS(NPROPS)
A=PROPS(3)
B=PROPS(4)
EN=PROPS(5)
HARD=0.0
IF(EQPLAS.EQ.0.0) THEN
SYIELD=A
ELSE
HARD=EN*B*EQPLAS**(EN-1)
SYIELD=A+B*EQPLAS**EN
END IF
RETURN
END注意事项（每项都是血的教训）：
1、Fortran语言遵循IN原则，即i、j、k、l、m、n（包括大写）都是整型变量，说白了就是整数，所以参数赋值的时候避开以这些开头。
2、最让我不能理解的是里面的“continue”如果换成“end do”就会出错，到现在我也不知道为什么，反正就老老实实打上编号用“continue”就完事了。
3、如果程序老出错可以多弄点中间变量，在后处理数据时好查错。
未完...
后续有时间可能会把验证子程序用的有限元模型和子程序里面的部分公式附上（可以用来对比Abaqus和UMAT子程序编的J-C），因为加上模型出错后，会分不清到底是模型错了，还是程序错了。
155-2731-8020
