1、深入剖析 UMAT线形粘弹性体liuxiongsimwe用户名：LXCAD摘要本文介绍了线形粘弹体的 UMAT 程序说明，与 ZHANG chunyu 一起学习 UMAT 为姐妹篇。UMAT 其实并不难学，你要把握几点即可：1， 必须提供准确的 雅可比距阵，程序收敛速度快2， 必须用增量法更新应力 must update the stresses and solution-dependent state variables to their values at the end of the increment for which it is called; must provide the m
2、aterial Jacobian matrix3，与本构方程相关的状态变量必须更新把我分析的子程序过程多读几遍，相关的弹性力学，朔性力学概念弄懂，可能就理解更为清晰！我的目的就是让大家学的轻松！请大家鼓励！三单元体的固体模型(线形粘弹性)下图为虎克体和开尔文体的串联一维方向的应力与应变的行为推广到实体其中PROPS(1) PROPS(2) PROPS(3) PROPS(4) PROPS(5)vSUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,1 RPL,DDSDDT,DRPLDE,DRPLDT,2 STRAN,DSTRAN,TIME,DTIME,
3、TEMP,DTEMP,PREDEF,DPRED,CMNAME,3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)CINCLUDE ABA_PARAM.INCCCHARACTER*80 CMNAMEDIMENSION STRESS(NTENS),STATEV(NSTATV),1 DDSDDE(NTENS,NTENS),2 DDSDDT(NTENS),DRPLDE(NTENS),3 STRAN(NTENS),DSTRAN(NTEN
4、S),TIME(2),PREDEF(1),DPRED(1),4 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)DIMENSION DSTRES(6),D(3,3)CC EVALUATE NEW STRESS TENSORCEV = 0.DEV = 0.DO K1=1,NDIEV = EV + STRAN(K1) ! 直接应力和 xyxvDEV = DEV + DSTRAN(K1) !直接增量应力和 xyEND DOCTERM1 = .5*DTIME + PROPS(5) ! vt2TERM1I = 1./TERM1 ! t1T
5、ERM2 = (.5*DTIME*PROPS(1)+PROPS(3)*TERM1I*DEV ! vt22TERM3 = (DTIME*PROPS(2)+2.*PROPS(4)*TERM1I ! vtC更新正应力DO K1=1,NDIDSTRES(K1) = TERM2+TERM3*DSTRAN(K1) 1 +DTIME*TERM1I*(PROPS(1)*EV2 +2.*PROPS(2)*STRAN(K1)-STRESS(K1)本构方程程序执行：+ +vt22t2)*2(vtSTRESS(K1) = STRESS(K1) + DSTRES(K1)n1END DOC 更新剪应力TERM2= vt2
6、TERM2 = (.5*DTIME*PROPS(2) + PROPS(4)*TERM1II1 = NDIDO K1=1,NSHRI1 = I1+1本构方程程序执行：* + * *xyvt2xytv21)(xyDSTRES(I1) = TERM2*DSTRAN(I1)+1 DTIME*TERM1I*(PROPS(2)*STRAN(I1)-STRESS(I1)STRESS(I1) = STRESS(I1)+DSTRES(I1)END DOC雅可比距阵开始C CREATE NEW JACOBIANCTERM2= *vt21)2)(tTERM2 = (DTIME*(.5*PROPS(1)+PROPS(
7、2)+PROPS(3)+1 2.*PROPS(4)*TERM1ITERM3 = *vt2)(tTERM3 = (.5*DTIME*PROPS(1)+PROPS(3)*TERM1IDO K1=1,NTENSDO K2=1,NTENSDDSDDE(K2,K1) = 0.END DOEND DOCDO K1=1,NDIDDSDDE(K1,K1) = TERM2END DOC填充距阵如下位置232termttertttDO K1=2,NDIN2 = K11DO K2=1,N2DDSDDE(K2,K1) = TERM3DDSDDE(K1,K2) = TERM3END DOEND DOTERM2 = *v
8、t21)(tTERM2 = (.5*DTIME*PROPS(2)+PROPS(4)*TERM1II1 = NDI !l1=3按新 TERM2填充距阵如下位置 22termtterDO K1=1,NSHR !按剪切力运行 3次I1 = I1+1DDSDDE(I1,I1) = TERM2END DO小结：雅可比距阵填充完毕C关于物质能量的增长C TOTAL CHANGE IN SPECIFIC ENERGYC )2(ETDE = 0.DO K1=1,NTENSTDE = TDE + (STRESS(K1)+.5*DSTRES(K1)*DSTRAN(K1)END DOC C关于朔性能量的增长C CH
9、ANGE IN SPECIFIC ELASTIC STRAIN ENERGYCTERM1 = PROPS(1) + 2.*PROPS(2) !TERM1= 2填充距阵如下位置 11termtterDO K1=1,NDID(K1,K1) = TERM1END DO填充距阵如下位置 DO K1=2,NDIN2 = K1-1DO K2=1,N2D(K1,K2) = PROPS(1)D(K2,K1) = PROPS(1)END DOEND DO直接应力部分外循环数内循环数K1K2TERM1 TERM2 DEE1 1 1 )1(*,D)1(*,(D2 1 2+2)(, +2,)(13 1 3+3*1+)
10、2(,D+3*,1+)2(D,)1()21131DD1 2 1 )1(*,)1(*(2 2 2+2)(,D+2,)(D23 2 3+3*+)2(,1+3*,+)2(1,)()21131DD)()1223211 3 1 )(*,D)(*3(D2 3 2+23)1(, +2,)1(33 3 3+*,D+)2(1,+3*,(D+)21(,)()21131DD)()122321)3(*)213313DD分析以下循环，实质为能量变化，力*位移DEE = 0.DO K1=1,NDITERM1 = 0.TERM2 = 0.DO K2=1,NDITERM1 = TERM1 + D(K1,K2)*STRAN(K
11、2)TERM2 = TERM2 + D(K1,K2)*DSTRAN(K2)END DODEE = DEE + (TERM1+.5*TERM2)*DSTRAN(K1)END DO剪切应力部分L1 K1 DDE4 1 44)2(*5 2 5544 )()(6 3 665544 )2(*)2(*)2(* I1 = NDIDO K1=1,NSHRI1 = I1+1DEE = DEE + PROPS(2)*(STRAN(I1)+.5*DSTRAN(I1)*DSTRAN(I1)END DO最终的 DDE为直接应力与剪切应力发生能量改变之和最终能量消散所改变的值 SSE,SCDSSE = SSE + DEE !朔性消散SCD = SCD + TDE DEE !徐变消散RETURNEND
