ABAQUS子程序UMAT里弹塑本构的实现

TAO_YZ= STRESS(6)

CEGMA_CP=(CEGMA_X+CEGMA_Y+CEGMA_Z)/THREE CEGMA_SX=CEGMA_X-CEGMA_CP CEGMA_SY=CEGMA_Y-CEGMA_CP CEGMA_SZ=CEGMA_Z-CEGMA_CP TAO_SXY=TAO_XY TAO_SZX=TAO_ZX TAO_SYZ=TAO_YZ

CEGMA_EQ=SQRT(THREE/TWO*(CEGMA_SX**2+CEGMA_SY**2+CEGMA_SZ**2+ 1TWO*(TAO_SXY**2+TAO_SYZ**2+TAO_SZX**2)))

DO k1=1,NTENS DO k2=1,NTENS DDSDDE(K1,K2)=ZERO END DO END DO

C 平均应力 C 6个应力偏量

c Mises等效应力

C 雅可比矩阵,初始化默认为

c 计算弹性矩阵 DO K1=1,NDI

DO K2=1,NDI

DDSDDE(K1,k2)=E*(ONE-Mu)/(ONE+Mu)/(ONE-TWO*Mu)*(Mu/(ONE-Mu)) END DO

DDSDDE(K1,K1)=E*(ONE-Mu)/(ONE+Mu)/(ONE-TWO*Mu) END DO

DO K1=NDI+1,NTENS

DDSDDE(K1,K1)=E*(ONE-Mu)/(ONE+Mu)/(ONE-TWO*Mu)*(ONE-TWO*Mu)/TWO/ END DO

1(ONE-Mu)

CCCCCCCCCCCCCCCCCCCCCCCCCC4：根据计算的Mises等效应力和读取的屈服应力Yield0比较，

CCCCCCCCCCCCCCCCCCCCCCCCCC如果Mises等效应力小于屈服应力，表明此时材料未屈服，那么转到，否则转到

C 按照弹性理论更新应力

DO K1=1,NTENS DO K2=1,NTENS

STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1) END DO

END DO ELSE

IF(CEGMA_EQ.LT.Yield0) THEN

CCCCCCCCCCCCCCCCCCCCCCCCCC5：没有屈服，按照弹性理论计算

CCCCCCCCCCCCCCCCCCCCCCCCCC6：屈服发生，按照塑性理论计算

31

C 切线模量H，根据本构关系求导 IF(EQPLAS.EQ.0) THEN

H=E ELSE

H=A*B*EQPLAS**(B-ONE) END IF

C 求W

W=9.D0*G/TWO/CEGMA_EQ**2/(H+THREE*G) C 等效塑性应变增量

DEQPLAS=(6.D0*G*CEGMA_EQ/(TWO*CEGMA_EQ**2*H+9.D0*G*(CEGMA_SX 1**2+CEGMA_SY**2+CEGMA_SZ**2+TWO*TAO_XY**2+TWO*TAO_ZX**2+TWO*TAO_YZ 2**2)))*(CEGMA_SX*DSTRAN(1)+CEGMA_SY*DSTRAN(2)+CEGMA_SZ* C

3DSTRAN(3)+TAO_XY*DSTRAN(4)+TAO_ZX*DSTRAN(5)+TAO_YZ*DSTRAN(6)) write(10,*) \

C 更新状态变量

STATEV(1)=EQPLAS+DEQPLAS C 计算塑性应变增量

DO K1=1,NTENS

DPLAS(1)= DEQPLAS*THREE*CEGMA_SX/TWO/CEGMA_EQ DPLAS(2)= DEQPLAS*THREE*CEGMA_SY/TWO/CEGMA_EQ DPLAS(3)= DEQPLAS*THREE*CEGMA_SZ/TWO/CEGMA_EQ DPLAS(4)= DEQPLAS*THREE*TAO_XY/CEGMA_EQ DPLAS(5)= DEQPLAS*THREE*TAO_ZX/CEGMA_EQ DPLAS(6)= DEQPLAS*THREE*TAO_ZY/CEGMA_EQ

C 按照塑性理论更新应力 DO K2=1,NTENS

STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*(DSTRAN(K1)-DPLAS(K1)) END DO

END DO END IF RETURN END

CCCCCCCCCCCCCCCCCCCCCCCCCC计算完成

5.4. 切线刚度法程序设计

算法设计

1：定义程序需要用到的常数和变量

2：读取ABAQUS定义的材料常数和状态变量(这里只定义了一个状态变量)，材料常数为，弹性模量E，泊松比Mu，屈服应力Yield0，参数A,B,C,并且计算出剪切模量G,状态变量为等效塑性应变EQPLAS

32

3：读取应力分量，计算平均应力，应力偏量以及Mises等效应力 平均应力：

?cp?(?x??y??z)/3

'??x??x??cp?'??y??y??cp?'??z??z??cp?'??xy??xy?'??yz??yz??'??zx应力偏量：?zx

??Mises等效应力：

2'2'2'2'2'2(?x??y??z'2?2(?xy??yz??zx))3

4：根据3计算的Mises等效应力和2读取的屈服应力Yield0比较，如果Mises等效应力小于屈服应力，表明此时材料未屈服，那么转到5，否则转到6 5：雅可比矩阵,初始化为0，然后计算弹性矩阵，按照弹性理论更新应力

?1????1??????1??E(1??)?De?(1??)(1?2?)?0???0???0??

d????Ded???1对10001?2?2(1??)001?2?2(1??)0称?1??000??????????????1?2??2(1??)??

6：雅可比矩阵,初始化为0

1）计算切线模量H

H'?(A(?P)B?C)'?A*B*?PB?1

33

注意到当等效塑性应?P?0时对应于本构关系的屈服点，此时的H不能通过上式计算，可以取此时的H为弹性模量 2）计算w

根据前一章推导的公式

??9G2?2(H'?3G)

3）计算等效塑性应变增量DEQPLAS，并更新状态变量 根据前一章推导的公式：

??T}[D]e?{?}d?p?d{?}??T??'H?{}[D]e?{?}?{?}

{[D]e??33G''''''?[D]e[?x,?y,?z,2?xy,2?yz,2?zx]T?[?x,?y,?z,?xy,?yz,?zx]T????2??T

??T??????([D]e)???[D]e????????????

带入后可得等效塑性应变增量DEQPLAS

d?p?'''''6*G*?*??x,?y,?z',?xy,?yz,?zx?2*?*H?9*G*(??????2??2??2?)2''2x'2y'2z'2xy'2yz'2zx*d???

然后EQPLAS=DEQPLAS+EQPLAS更新状态变量 4）计算雅可比矩阵

首先初始化默认为0，然后用下式计算雅可比矩阵

34