当前位置：文档库 › UMAT-Neo-Hookean

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)

INCLUDE ’ABA_PARAM.INC’

CHARACTER*8 CMNAME

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 (AN

D EELAS)

C DISTGR - DEVIATORIC DEFORMATION GRADIENT (DISTORTION TENSOR)

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 USE

D FOR PLAN

E 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 AN

D DISTORTION TENSOR C

1.程序里面的DET=Z

z Y z X z Z

y Y y X y

Z x Y x X x F J ??????????????????=

=)det(,当平面应变问题时,Z y

Z x ????,=0 2. DISTGR=F dev =J -1/3.F 为偏量变形梯度, 体积变形梯度 F V0L =J 1/3.I, 且F VOL .F dev =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

C CALCULATE DEVIATORIC LEFT CAUCHY-GREEN DEFORMATION TENSOR C

3.计算左Cauchy-Green 变形张量B=F .F T , 那么左Cauchy-Green 变形偏量

()

T T

dev

dev DISTGR DISTGR F F BBAR B ?=?==,且B J B 3/2-=)4(),1(B B 表示平

面应变问题, 加上)6(),5(B B 表示空间问题,排列按照运用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. 计算应力

上式为应变能函数, 21,I I 分别为第一第二应变不变量.对应左Cauchy-Green 变形张量B ,且

)

(),(21B B trace I B trace I ?==)1(2

,0,1

2101-=??=??=??J D J U I U C I U , 应力偏量()

V B B J C B B I U B I U

I I U D E V J S 01022

22-=??

???????-???? ????+??=

,21,I I ,对应左Cauchy-Green

变形偏量B ,且13/23

/21)()(I J B J trace B trace I --===,因此,

3/21I U

J I U ??=??,又因为pI S +=σ,)1(2)(311321--=??-=++-=J D J U p σσσ(静水

压力)

TRBBAR=(BBAR(1)+BBAR(2)+BBAR(3))/THREE !!!VOL

B TRBBAR = 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(ijkl C )

()kl

ij mm kl ij kl ij kl ij jk il jk il jl ik jl ik J ijkl J D B B B B B B B C C δδδδδδδδδδ)12(2923232211

102

-+??? ??+--+++=[]????

??????=??????????=12121222

1211

22122222

2211

111211221111

3332

2322

131211

C C C C C C C C C C C C C C C C C C C (平面问题)

按照V oigt 规则

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

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 (AN

D EELAS)

=, SCALE=DET**(-ONE/THREE)

B3/2-

!!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 STORE ELASTIC STRAINS IN STATE VARIABLE ARRAY

7. 把弹性应变存为状态变量

DO K1=1, NTENS

STATEV(K1)=EELAS(K1)

END DO

RETURN

END