5. UMAT程序设计和编码
本章将严格按照前一章推导的公式展开程序设计和编码，为了便于编程，本文将本构关系做了抽象化处理，即将其描述成一个含参数的表达式，改变参数即可应用于不同的模型，这样做的好处是能保证程序的复用性，这也是本文反复强调的使用UMAT的原则。
5.1. 本构关系描述
本文采用各向同性硬化弹塑性材料，材料参数如下：
???E?e (??B??A(?)?C (??P
????TT)) (5-1)
弹性部分：??E?e，弹性模量E=200000Mpa，泊松比Mu=0.3
图5-1 弹性部分本构关系
塑性部分：
??A(?P)?CB，为了研究方便，取A=700，B=0.5，C=400
25
图5-2 塑性部分本构关系
将2个曲线统一到同一个坐标系（为方便显示，x轴标注时扩大了1000倍）
图5-3 本构关系
由此可以求出两条线的交点即初始屈服点的应力Yield0=400Mpa（注意两条曲线相差了一个屈服应变，因为两者其实不是一个坐标系）
综上定义的材料常数见表5-1：
表5-1
弹性模量E 泊松比Mu
200000 0.3
屈服应力Yield0 400 A
700
26
B C
0.5 400
注意：上面A,B,C的取值只是为了便于理解和分析本文的材料模型，为了保证程序的 通用性，本文的参数在程序中的A,B,C一律用变量表示。
5.2. 常刚度法程序设计
算法设计
1：定义程序需要用到的常数和变量
2：读取ABAQUS定义的材料常数和状态变量(这里只定义了一个状态变量)，材料常数为，弹性模量E，泊松比Mu，屈服应力Yield0，参数A,B,C,并且计算出剪切模量G,状态变量为等效塑性应变EQPLAS
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，计算弹性矩阵，按照弹性理论更新应力
27
?1????1??????1??E(1??)?De?(1??)(1?2?)?0???0???0??1对10001?2?2(1??)001?2?2(1??)0称?1??000??????????????1?2??2(1??)??
d????Ded???
6：雅可比矩阵,初始化为0
1）：计算切线模量H'
H?A*B*(?P?Yield0/E)'B?1,注意到当等效塑性应?P?0时对应于本构关系的屈服点，此
时的H不能通过上式计算，可以取此时的H为弹性模量 2）：计算等效塑性应变增量并更新
d?p?6*G*?*??x,?y,?z,?xy,?yz,?zx?'2'2'2'2'2'22*?*H?9*G*(?2''2x??'2y??'2z?2?'2xy?2?'2yz?2?)'2zx*d???
EQPLAS=DEQPLAS+EQPLAS更新状态变量 3）计算流动方向
???????32?''
'T[?x,?y,?z,2?xy,2?yz,2?zx]
4）计算塑性应变增量
d{?}p?d?p???{?}
5）更新应力
d{?}?[D]e(d{?}?d{?}p)
28
算法流程图
图5-4 常刚度法算法流程图
5.3. 常刚度法程序编码
根据算法流程，用FORTRAN77 固定格式编制了常刚度法的计算程序如下：
SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,
1 DRPLDE,DRPLDT,STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,
2 CMNAME,NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT, 3 PNEWDT,CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
29