模型介绍:
弹性模量20GPa,泊松比0.3,延续上一节Umat的模型,硬化曲线方程:
本次代码中用到了13个状态变量,1~6个为弹性应变,7~12个为塑性应变,13为等效塑性应变。
代码讲解:
部分代码如下(篇幅有限):
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)
double precision, DIMENSION(NTENS)::deviatoricStress
double precision misesEqualStress,hydroStaticPressure,W
1 ,ONESY,DEQPL,SYIELD
C
DIMENSION EELAS(6),EPLAS(6),FLOW(6)
PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0)
DATA NEWTON,TOLER/1000,1.D-6/
INTEGER K1,K2,cn
C PROPS(1) - E
C PROPS(2) - NU
C PROPS(3) - SYIELD
C CALLS AHARD FOR CURVE OF SYIELD VS. PEEQ
C ELASTIC PROPERTIES
C
EMOD=PROPS(1)
EMU=PROPS(2)
EG2=EMOD/(ONE+EMU)
EG=EG2/TWO
EG3=THREE*EG
ELAM=EMOD*EMU/((ONE+EMU)*(ONE-TWO*EMU))
C
C ELASTIC STIFFNESS
C
DO K1=1,NTENS
DO K2=1,NTENS
DDSDDE(K2,K1)=0.0
enddo
enddo
C
DO K1=1,NDI
DO K2=1,NDI
DDSDDE(K2,K1)=ELAM
enddo
DDSDDE(K1,K1)=EG2+ELAM
enddo
DO K1=NDI+1,NTENS
DDSDDE(K1,K1)=EG
enddo
......
107~119行:(其中112~119行的循环结构要看懂)
123~133行:
138~145行:调用等效Mises应力,更新弹塑性雅克比矩阵;
150~154行:将弹性应变、塑性应变、等效塑性应变储存在状态变量中。
160~185行:编写计算非线性强化系数子程序;
187~217行:编写等效Mises应力子程序。
本次的代码原本要为大家验证一下的,打开Abaqus的时候才发现子程序关联失常,就把代码先写出来吧,有兴趣大家可以自行验证,然后一起交流。
来源:易木木响叮当