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),DPLAS(6) C----------------------------
CCCCCCCCCCCCCCCCCCCCCCCCCC1：定义变量 c 定义常数
PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0)
常刚度法
C 定义材料常数
DOUBLE PRECISION E,Mu,Yield0,A,B,C
DOUBLE PRECISION H,W,CEGMA_EQ,CEGMA_X,CEGMA_Y,CEGMA_Z,TAO_XY, 1TAO_YZ,TAO_ZX,CEGMA_CP,CEGMA_SX,CEGMA_SY,CEGMA_SZ,TAO_SXY,TAO_SYZ
C 定义中间变量
2,TAO_SZX C 定义状态变量 C
DOUBLE PRECISION EQPLAS 定义更新变量
DOUBLE PRECISION DEQPLAS CCCCCCCCCCCCCCCCCCCCCCCCCC2：读取变量 C C
弹性模量 泊松比 Mu=PROPS(2)
E=PROPS(1)
C 剪切模量 C
G=E/TWO/(ONE+Mu) 屈服应力 Yield0=PROPS(3)
C 参数
A=PROPS(4) B=PROPS(5) C=PROPS(6) EQPLAS=STATEV(1)
Yield=A*(EQPLAS+Yield0/E)**B+C
C 读取状态参数
CCCCCCCCCCCCCCCCCCCCCCCCCC3：读取应力分量，计算平均应力，应力偏量以及Mises等效应力 C 6个应力分量
CEGMA_X= STRESS(1) CEGMA_Y= STRESS(2)
CEGMA_Z= STRESS(3)
TAO_XY= STRESS(4) TAO_ZX= STRESS(5)
30
TAO_YZ= STRESS(6)
C 平均应力
CEGMA_CP=(CEGMA_X+CEGMA_Y+CEGMA_Z)/THREE
C 6个应力偏量
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 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
IF(CEGMA_EQ.LT.Yield0) THEN
CCCCCCCCCCCCCCCCCCCCCCCCCC5：没有屈服，按照弹性理论计算
DO K2=1,NTENS
STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1) END DO
END DO ELSE
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
CCCCCCCCCCCCCCCCCCCCCCCCCC计算完成
RETURN END
5.4. 切线刚度法程序设计
算法设计
1：定义程序需要用到的常数和变量
2：读取ABAQUS定义的材料常数和状态变量(这里只定义了一个状态变量)，材料常数为，弹性模量E，泊松比Mu，屈服应力Yield0，参数A,B,C,并且计算出剪切模量G,状态变量为等效塑性应变EQPLAS
32
3：读取应力分量，计算平均应力，应力偏量以及Mises等效应力 平均应力：
?cp?(?x??y??z)/3??x??cp??y??cp??z??cp??xy??yz??zx23
'??x?'??y?'??z?'??xy?'??yz?'应力偏量：??zx
(?x??'2'2y??Mises等效应力：
??z?2(?xy??yz??zx))'2'2'2'2
4：根据3计算的Mises等效应力和2读取的屈服应力Yield0比较，如果Mises等效应力小于屈服应力，表明此时材料未屈服，那么转到5，否则转到6 5：雅可比矩阵,初始化为0，然后计算弹性矩阵，按照弹性理论更新应力
?1????1??????1??E(1??)?De?(1??)(1?2?)?0???0???0????????????????1?2??2(1??)??
1对10001?2?2(1??)001?2?2(1??)0称?1??000
d????Ded???
6：雅可比矩阵,初始化为0
1）计算切线模量H
H?(A(?P)?C)?A*B*?P'B'B?1
33
注意到当等效塑性应?P?0时对应于本构关系的屈服点，此时的H不能通过上式计算，可以取此时的H为弹性模量 2）计算w
根据前一章推导的公式
??9G2?(H?3G)2'
3）计算等效塑性应变增量DEQPLAS，并更新状态变量 根据前一章推导的公式：
{??}[D]eTd?p??{?}d{?}??T??'H?{}[D]e?{?}?{?}
[D]e???????[D]e32?[?x,?y,?z,2?xy,2?T''',2?zx]yzT?3G?[?x,?y,?z,?xy,?'''yz,?zx]T
???T?????([D]e)???[D]e????????????
带入后可得等效塑性应变增量DEQPLAS
d?p?6*G*?*??x,?y,?z,?xy,?yz,?zx?''''''2*?*H?9*G*(?x??2''2'2y??z?2?xy?2?yz?2?zx)'2'2'2'2*d???
然后EQPLAS=DEQPLAS+EQPLAS更新状态变量 4）计算雅可比矩阵
首先初始化默认为0，然后用下式计算雅可比矩阵
34