1、如文档存 在问 题，请发 邮件 至：n uaal i zhe n 1 26 .com 版 权所有 ，翻版必 究 适用于金属材料的 UMAT 子程序公式详解 * 注：本文如未特殊说明，标注的公式和章节均来自 Dunne 的Introduction to computational plasticity 。 参考文献 ： Dunne, Fionn, and Nik Petrinic. Introduction to computational plasticity. New York: Oxford University Press, 2005. * P2 ：JOHNSON-COOK 模型 Joh
2、nson-Cook 模型是 一种经典的本构模型， 适用于模拟高应变率 （如冲击、 弹道侵彻等） 下的金属材 料，经过适 当修正后同 样适用于陶 瓷等非金属 材料。Johnson-Cook 强化 模型假 设 材料是各向同性硬化(Isotropic hardening) ，表示为三项的乘积，分别反映了应变硬化、应变率 硬化和温度软化，表达式为： room 0r o o m () 1 l n ( ) 1 () nm m TT AB c TT 其中 为等效流动应力； 为等效塑性应变； 为等效塑性应变率 （ 1 s ） ； 0 为参考应变率， 一般取 1 1s ； room T 为室温； m T 为材料
3、熔点温度。 J-C 模型中共含有 5 个经 验参数，各参数的物理意义为：屈服强度(MPa) ；B 为硬化 模量 (MPa) ；c 为经验性应变率敏感系数；n 为硬化系 数；m 为温 度软化效应。 本文程序不考虑温度软化效应，即温度始终保持室温，即 room TT ，此时 J-C 模型变为： 0 () 1 l n ( ) n AB c 并且，本文 UMAT 程序不考虑塑性耗散。 参考文献 ： 1. Johnson, Gordon R., and William H. Cook. “Fracture characteristics of three metals subjected to vari
4、ous strains, strain rates, temperatures and pressures.“Engineering fracture mechanics 21.1 (1985): 31-48. 2.Umbrello, D., R. Msaoubi, and J. C. Outeiro. “The influence of JohnsonCook material constants on finite element simulation of machining of AISI 316L steel.“ International Journal of Machine To
5、ols and Manufacture 47.3 (2007): 462-470. 3. 范亚夫,段祝 平. “Johnson-Cook 材料 模型参数的实验测定.“ 力学与实践 25.5 (2003): 40-43. 4. 庄茁, 由小 川, 廖剑辉. 基于 ABAQUS 的有限元分析和应用. 清华大学出版社, 2009. 5. 常列珍， 潘玉田， 张治民等. 一种调质 50SiMnVB 钢 Johnson-Cook 本构模 型的建立J. 兵器材 料科学与工程,2010,33(4):68-72. * P3 ：材料参数 NDI: 直接应力分量个数。本文的 UMAT 只适用 于单元有三个直接应力
6、分量。 ENU: 泊松比。判断泊松比大小，避免计算得到的 EBULK3 值无穷大。 K ：体积模量，elastic bulk moduls ，详 见 Chapter 5.2.1 。 如文档存 在问 题，请发 邮件 至：n uaal i zhe n 1 26 .com 版 权所有 ，翻版必 究 EG3 ：在应力更新部分用到，只是简便。 ELAM ：The Lame constant 11 (3 2 ) ( ) 3 312 1 ( 12) ( 1 ) EE E KG * P4 ：计算刚度矩阵 2 2 2 G G G C G G G 详见 Chapter 5.3.1 注：此处刚度阵 C 暂时 存储在
7、 DDSDDE 数组中。 * P4 ：计算试探应力 试探应力为公式(5.7) 中的 Elastic predictor ，表达 式由公式(6.11) 可得： 2G ( ) : tr t t I I C 详见 Chapter 5.3.1 和 Chapter 6.4.2 注：试探应力存储在 STRESS 数组中。 * P5 ：读取状态变量 读取状态变量矩阵中弹性应变值、塑性应变值和等效塑性应变。 EELAS: 弹性应 变 ee tt t EPLAS： 塑性应 变 pp tt t EQPLAS: 等 效塑性应变 如文档存 在问 题，请发 邮件 至：n uaal i zhe n 1 26 .com 版
8、 权所有 ，翻版必 究 tt t p p * P5 ：计算 V on Mises 应力 V on Mises 应力表达式为 2222 2 2 1 () () () 6 ( ) 2 tr ex yy zz xx y y z z x 注：此应力存储在 SMISES 数组中。 * P5 ：判断是否屈服 如果 V on Mises 应力 tr e 大于 初始屈服应力 0 y ，则屈服发生，定义流动法则。 SHYDRO ：为静水压应力（hydrostatic stress ） ，为 1 () 3 tr Tr I 。 ONESY ：为 1 tr e FLOW ：为 tr tr e ，其中 tr 为试探应力
9、 tr 的偏应力量， 1 () 3 tr tr tr Tr I 。 * P6 ：Newton 迭代（隐式应力更新径向返回方法） SYIELD: 屈服应力，表示为 0 yy DEQPL ：初 始等效塑性应变增量，表示为 () 3 tr ey k p G Newton 迭代 结束后，需使得屈服方程满足 () 30 tr k eye y fG p DSTRES ：表 示为 0 3 y tol G 其中 tol 为容 差。 DEQMIN ：为预设等效塑性应变增量值，保证迭代时增量不至于过小。 如文档存 在问 题，请发 邮件 至：n uaal i zhe n 1 26 .com 版 权所有 ，翻版必 究
10、 %Newton iteration% 调用 USERHARD 程序，求得当前等效塑性应变下的 () 0 () k y pp 和应变硬化系数 () 0 () k Hpp ： () 0 () k yy pp TVP ：表示为 () ln( ) k p c t TVP1 ：表示 为 () 1l n () k p c t HARD1 ：屈服面硬化系数表示为 () 0 10 () 1 ln ( ) k yy y k c p HH c ppt p SYIELD: 屈服应力更新为 () 1 ln( ) k yy p c t RHS ：表示为 () 3 tr k ey Gp DEQPL ：更 新为 (1 )
11、 () kk p pd p 其中： () 1 3 3 tr k ey Gp dp GH 如 ABS(RHS/EG3) 小于 DSTRES 则退出 Newton 迭 代，即 () 0 3 33 tr k eyy Gp t o l GG %Newton iteration% EFFHRD: 用于 Jacobian 矩阵计算， 表示为 1 1 3 3 GH GH 如文档存 在问 题，请发 邮件 至：n uaal i zhe n 1 26 .com 版 权所有 ，翻版必 究 注：DEQPL 的推导过程（参考公式（5.18 ） 、 （5.19 ）和（5.20 ） ） ： 屈服方程为： () 30 tr
12、k eye y fG p 屈服方程进行泰勒展开： 0 f fd p p 则可得到： () 0 () 0 () 3 3 1 l n ( ) 0 k y tr k ey k c p Gp GH c dp tp () 1 3 3 tr k ey Gp dp GH * P7: 更新 应力 和应变 由公式（6.14 ）可知： 3 2 2: tr p tr e ep ee p GII 由公式（6.15 ）可知： t t p pp STRESS ：应 力 1 () 3 tr tr tt y tr e Tr I EPLAS： 塑性应 变 ppp tt tt EELAS: 弹性应 变 eep tt tt EQP
13、LAS: 等 效塑性应变 tt tt p pp RPL ：塑性功转化为热量的热量值 如文档存 在问 题，请发 邮件 至：n uaal i zhe n 1 26 .com 版 权所有 ，翻版必 究 注： 1. 应力更新的 推导过程（参考公式（5.18 ） 、 （5.19）和 （5.20 ） ） ： 由： 1 () 3 tr tr tr Tr I 又由于 Newton 迭代后： 3 tr ye Gp 则可得： 1 () 3 (3) 3 tr tr tt y tr e tr tr tr tr e tr e tr tr tr e Tr I Gp Gp 将公式（6.11 ）和（6.14 ） 代入： 2
14、2: 3 3 2: p tt t e t GIIG GII 由于： :() 0 pp II T r 则： 2: ee tt t t GII 2. 在 UMAT 中，剪切应变为工程应变。 * P7: 计算 Jacobian 矩阵 EFFG ：表示 为 y e tr tr ee GGG R EFFG2 ：表示为 2GR EFFG3 ：表示为 3GR EFFLAM ：表示为 12 (3 2 ) ( ) 33 KG RKG R EFFHRD-EFFG3 ：表 示为 1 11 33 1 32( ) 2 32 1 ( 3 ) e tr e GH GR G GQ GH GH 如文档存 在问 题，请发 邮件
15、至：n uaal i zhe n 1 26 .com 版 权所有 ，翻版必 究 Jacobian 矩阵 表达式( 参考公式（5.22）至 （5.31 ） ： 2 22 () 3 tr tr tr tr ee JG QG R I K G R I I 注： R 和 Q 的定义参考公式（5.27 ） * P8: 计算 硬化 率和屈服应力 本文程序中的 J-C 模型表达式为： 0 ( )1 ln( ) ( )1 ln( ) n yy p AB p cp p c t 其中， 0 () n y p AB p 。 当 p=0 时， 0 () y p A ； 当 p0 时， 0 () n y p AB p ，
16、应变硬化系数 0 1 0 () y n p Hn B p p 。 注： 1.USERHARD 子程序计 算的只是 0 () y p 和 0 H ，在初始时 ，将 0 () y p 认为是当前的屈 服应力 y 。而在 Newton 迭 代中会产生应变率，才开始考虑公式第二项的应变率硬化。 2. 由 Chapter 2.2.3 结论可 知， 对于单轴拉伸情况， p 。 试件受 到高塑性应变作用时 ， ep 且 p 。 3. FORTRAN 中是地址调用（CALL BY ADDRESS/REFERENCE ） 。主程序 和子程序对应 的参数名可以不同，但都使用同一个地址。 主程序： CALL USE
17、RHARD(SYIEL0,HARD,EQPLAS,PROPS(4) CALL USERHARD(SYIELD,HARD,EQPLAS+DEQPL,PROPS(4) 子程序： SUBROUTINE USERHARD(SYIELD,HARD,EQPLAS,TABLE) 则传递参数的对应关系为： 主程序 子程序 A PROPS(4) TABLE(1) B PROPS(5) TABLE(2) n PROPS(6) TABLE(3) 参考文献： 彭国伦. Fortran 95 程序设 计M. 中国 电力出版社, 2002. * 如文档存 在问 题，请发 邮件 至：n uaal i zhe n 1 26 .com 版 权所有 ，翻版必 究 P0: UMAT 子 程序流程 图 UMAT 计算流程可以参考 Chapter 6.4.2 UMAT 材料常数和状态变量矩阵如下所示： 参考文献 ： 庄茁, 由小川, 廖剑辉. 基于 ABAQUS 的有限元分析和应用. 清 华大学出版社, 2009. *
