UMAT-Neo-Hookean

UMAT 分析Neo-Hookean材料 !由于本人学ABA不久, 到现在为止, 对umat还是处于半懂状态, 从ABA论坛下载了大大们提供的资料和读了一些大大的经验之谈, 受益非浅! 以下这段是我关于neo-hookean材料子程序的部分阐释,不当之处还望各为大大指正.(由于本人水平有限,但本着大家互相帮助的出发点,所以抛点砖头出来,^_^) Neo-Hookean材料模型是各向同性线性定律(Hookean定律)至大变形的扩展. 在关键词*ELASTIC中定义弹性模量不能有效的解决有限变形问题,必须定义一个合适的势能函数 SUBROUTINE UMAT(STRESS, STATEV, DDSDDE, SSE, SPD, SCD, RPL, 1 DDSDDT, DRPLDE, DRPLDT, STRAN, DSTRAN, TIME, DTIME, TEMP, DTEMP, 2 PREDEF, DPRED, CMNAME, NDI, NSHR, NTENS, NSTATV, PROPS, NPROPS, 3 COORDS, DROT, PNEWDT, CELENT, DFGRD0, DFGRD1, NOEL, NPT, LAYER, 4 KSPT, KSTEP, KINC) C INCLUDE ’ABA_PARAM.INC’ C CHARACTER*8 CMNAME C DIMENSION STRESS(NTENS), STATEV(NSTATV), DDSDDE(NTENS, NTENS), 1 DDSDDT(NTENS), DRPLDE(NTENS), STRAN(NTENS), DSTRAN(NTENS), 2 PREDEF(1), DPRED(1), PROPS(NPROPS), COORDS(3), DROT(3, 3), 3 DFGRD0(3, 3), DFGRD1(3, 3) C LOCAL ARRAYS C ---------------------------------------------------------------- C EELAS - LOGARITHMIC ELASTIC STRAINS C EELASP - PRINCIPAL ELASTIC STRAINS C BBAR - DEVIATORIC RIGHT CAUCHY-GREEN TENSOR RIGHT 应为LEFT?? C BBARP - PRINCIPAL VALUES OF BBAR C BBARN - PRINCIPAL DIRECTION OF BBAR (AND EELAS) C DISTGR - DEVIATORIC DEFORMATION GRADIENT (DISTORTION TENSOR) C ---------------------------------------------------------------- C DIMENSION EELAS(6), EELASP(3), BBAR(6), BBARP(3), BBARN(3, 3), 1 DISTGR(3,3) C PARAMETER(ZERO=0.D0, ONE=1.D0, TWO=2.D0, THREE=3.D0, FOUR=4.D0, 1 SIX=6.D0) C C ---------------------------------------------------------------- C UMAT FOR COMPRESSIBLE NEO-HOOKEAN HYPERELASTICITY C CANNOT BE USED FOR PLANE STRESS C ---------------------------------------------------------------- C PROPS(1) - E C PROPS(2) – NU C ---------------------------------------------------------------- C C ELASTIC PROPERTIES C EMOD=PROPS(1) ENU=PROPS(2) C10=EMOD/(FOUR*(ONE+ENU)) D1=SIX*(ONE-TWO*ENU)/EMOD C C JACOBIAN AND DISTORTION TENSOR C 1.程序里面的DET= ,当平面应变问题时, =0 2. DISTGR=Fdev=J-1/3.F 为偏量变形梯度, 体积变形梯度 FV0L=J1/3.I, 且FVOL.Fdev=F DET=DFGRD1(1, 1)*DFGRD1(2, 2)*DFGRD1(3, 3) 1 -DFGRD1(1, 2)*DFGRD1(2, 1)*DFGRD1(3, 3) IF(NSHR.EQ.3) THEN DET=DET+DFGRD1(1, 2)*DFGRD1(2, 3)*DFGRD1(3, 1) 1 +DFGRD1(1, 3)*DFGRD1(3, 2)*DFGRD1(2, 1) 2 -DFGRD1(1, 3)*DFGRD1(3,1)*DFGRD1(2, 2) 3 -DFGRD1(2, 3)*DFGRD1(3, 2)*DFGRD1(1, 1) END IF SCALE=DET**(-ONE/THREE) DO K1=1, 3 DO K2=1, 3 DISTGR(K2, K1)=SCALE*DFGRD1(K2, K1) END DO END DO C CALCULATE DEVIATORIC LEFT CAUCHY-GREEN DEFORMATION TENSOR C 3.计算左Cauchy-Green变形张量B=F.FT, 那么左Cauchy-Green变形偏量 ,且 EMBED Equation.3 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 示平面应变问题, 加上 表示空间问题,排列按照运用Voigt规则. BBAR(1)=DISTGR(1, 1)**2+DISTGR(1, 2)**2+DISTGR(1, 3)**2 BBAR(2)=DISTGR(2, 1)**2+DISTGR(2, 2)**2+DISTGR(2, 3)**2 BBAR(3)=DISTGR(3, 3)**2+DISTGR(3, 1)**2+DISTGR(3, 2)**2 BBAR(4)=DISTGR(1, 1)*DISTGR(2, 1)+DISTGR(1, 2)*DISTGR(2, 2) 1 +DISTGR(1, 3)*DISTGR(2, 3) IF(NSHR.EQ.3) THEN BBAR(5)=DISTGR(1, 1)*DISTGR(3, 1)+DISTGR(1, 2)*DISTGR(3, 2) 1 +DISTGR(1, 3)*DISTGR(3, 3) BBAR(6)=DISTGR(2, 1)*DISTGR(3, 1)+DISTGR(2, 2)*DISTGR(3, 2) 1 +DISTGR(2, 3)*DISTGR(3, 3) END IF C C CALCULATE THE STRESS 4. 计算应力 上式为应变能函数, 分别为第一第二应变不变量.对应左Cauchy-Green变形张量 ,且 , 应力偏量 , ,对应左Cauchy-Green变形偏量 ,且 ,因此, ,又因为 , (静水压力) TRBBAR=(BBAR(1)+BBAR(2)+BBAR(3))/THREE !!! EG=TWO*C10/DET EK=TWO/D1*(TWO*DET-ONE) PR=TWO/D1*(DET-ONE) !!! PR=-P DO K1=1,NDI STRESS(K1)=EG*(BBAR(K1)-TRBBAR)+PR !!!(我下载两个pdf的文档,其中一个名为umat_all中的程序与这不一样,我个人认为这个是对的.) END DO DO K1=NDI+1,NDI+NSHR STRESS(K1)=EG*BBAR(K1) END DO C CALCULATE THE STIFFNESS C 5. 计算DDSDDE( ) 平面问题 按照Voigt规则 EG23=EG*TWO/THREE DDSDDE(1, 1)= EG23*(BBAR(1)+TRBBAR)+EK DDSDDE(2, 2)= EG23*(BBAR(2)+TRBBAR)+EK DDSDDE(3, 3)= EG23*(BBAR(3)+TRBBAR)+EK DDSDDE(1, 2)=-EG23*(BBAR(1)+BBAR(2)-TRBBAR)+EK DDSDDE(1, 3)=-EG23*(BBAR(1)+BBAR(3)-TRBBAR)+EK DDSDDE(2, 3)=-EG23*(BBAR(2)+BBAR(3)-TRBBAR)+EK DDSDDE(1, 4)= EG23*BBAR(4)/TWO DDSDDE(2, 4)= EG23*BBAR(4)/TWO DDSDDE(3, 4)=-EG23*BBAR(4) DDSDDE(4, 4)= EG*(BBAR(1)+BBAR(2))/TWO IF(NSHR.EQ.3) THEN DDSDDE(1, 5)= EG23*BBAR(5)/TWO DDSDDE(2, 5)=-EG23*BBAR(5) DDSDDE(3, 5)= EG23*BBAR(5)/TWO DDSDDE(1, 6)=-EG23*BBAR(6) DDSDDE(2, 6)= EG23*BBAR(6)/TWO DDSDDE(3, 6)= EG23*BBAR(6)/TWO DDSDDE(5, 5)= EG*(BBAR(1)+BBAR(3))/TWO DDSDDE(6, 6)= EG*(BBAR(2)+BBAR(3))/TWO DDSDDE(4,5)= EG*BBAR(6)/TWO DDSDDE(4,6)= EG*BBAR(5)/TWO DDSDDE(5,6)= EG*BBAR(4)/TWO END IF DO K1=1, NTENS DO K2=1, K1-1 DDSDDE(K1, K2)=DDSDDE(K2, K1) END DO END DO C C CALCULATE LOGARITHMIC ELASTIC STRAINS (OPTIONAL) 6.计算对数应变 CALL SPRIND(BBAR, BBARP, BBARN, 1, NDI, NSHR) C EELAS - LOGARITHMIC ELASTIC STRAINS C EELASP - PRINCIPAL ELASTIC STRAINS C BBAR - DEVIATORIC RIGHT CAUCHY-GREEN TENSOR !!RIGHT 应为LEFT?? C BBARP - PRINCIPAL VALUES OF BBAR C BBARN - PRINCIPAL DIRECTION OF BBAR (AND EELAS) , SCALE=DET**(-ONE/THREE) !!EELASP(1~3)为左变形张量B主值的对数值 然后计算弹性应变的各个分量,公式的推导很多弹性力学书上有! EELASP(1)=LOG(SQRT(BBARP(1))/SCALE) EELASP(2)=LOG(SQRT(BBARP(2))/SCALE) EELASP(3)=LOG(SQRT(BBARP(3))/SCALE) EELAS(1)=EELASP(1)*BBARN(1,1)**2+EELASP(2)*BBARN(2, 1)**2 1 +EELASP(3)*BBARN(3, 1)**2 EELAS(2)=EELASP(1)*BBARN(1, 2)**2+EELASP(2)*BBARN(2, 2)**2 1 +EELASP(3)*BBARN(3, 2)**2 EELAS(3)=EELASP(1)*BBARN(1, 3)**2+EELASP(2)*BBARN(2, 3)**2 1 +EELASP(3)*BBARN(3, 3)**2 EELAS(4)=TWO*(EELASP(1)*BBARN(1, 1)*BBARN(1, 2) 1 +EELASP(2)*BBARN(2, 1)*BBARN(2, 2) 2 +EELASP(3)*BBARN(3, 1)*BBARN(3, 2)) IF(NSHR.EQ.3) THEN EELAS(5)=TWO*(EELASP(1)*BBARN(1, 1)*BBARN(1, 3) 1 +EELASP(2)*BBARN(2, 1)*BBARN(2, 3) 2 +EELASP(3)*BBARN(3, 1)*BBARN(3, 3)) EELAS(6)=TWO*(EELASP(1)*BBARN(1, 2)*BBARN(1, 3) 1 +EELASP(2)*BBARN(2, 2)*BBARN(2, 3) 2 +EELASP(3)*BBARN(3, 2)*BBARN(3, 3)) END IF C C STORE ELASTIC STRAINS IN STATE VARIABLE ARRAY 7. 把弹性应变存为状态变量 DO K1=1, NTENS STATEV(K1)=EELAS(K1) END DO C RETURN END _1258273840.unknown _1258274117.unknown _1258350509.unknown _1258352567.unknown _1258353070.unknown _1258353439.unknown _1258352101.unknown _1258350256.unknown _1258350385.unknown _1258274624.unknown _1258274363.unknown _1258273956.unknown _1258274023.unknown _1258273893.unknown _1258268949.unknown _1258270422.unknown _1258272309.unknown _1258273728.unknown _1258270954.unknown _1258270275.unknown _1258268873.unknown