5.本构模型-应力更新专题-UMAT和VUMAT

1、计算固体力学计算固体力学应力更新专题应力更新专题 柳占立柳占立 庄茁庄茁 2022年年4月月13日日导出有限元方程（完全，更新、共旋）不同用户要求： 知道有限元基本流程； 编写本构关系、节点内力 开发有限元程序第4-6章求解有限元方程（隐式和显式）有限元基本知识有限元基本知识有限元高级知识有限元高级知识单元技术、结构单元、接触等第7-10章应力更新应力更新1.1. 什么是应力更新什么是应力更新2.2. 几种客观率几种客观率3.3. 前推后拉及前推后拉及LieLie导数导数4.4. 常用客观率应力更新流程常用客观率应力更新流程5.5. ABAQUSABAQUS用户自定义材料用户自定义材料理论理论

2、应用应用有限元程序有限元程序赋予力学以活力赋予力学以活力0()( )ttttdt dBTITIfint1.1.应力更新应力更新应力增量应力增量xyttt , FRDD0t, FRTJFFS1TTTJJJFFFFFFS111 FSFTJsL JsJsLT sT1JFDDtF1JsFFT svsL ssLT tT tL ttLT FDDtF1tFTFTLvt2.2.几种客观率几种客观率TruesdellTruesdell率：率：前推LieLie导数导数tt Js s:SESCE2.2.几种客观率几种客观率Green-Green-NaghdiNaghdi率：率：前推LieLie导数导数TTTTGDt

3、DRRRRRRRR TR RGTDDtRRGT :DsCDTRR共轴旋转应力率共轴旋转应力率2.2.几种客观率几种客观率JTW W1112TWRU UUURGT JaumannJaumann率：率： = W Rigid body rotation or The principal directions of the strain remaining constantJTW W如何得到如何得到2.2.几种客观率几种客观率TTL LvGT TruesdellTruesdell率：率：Green-Green-NaghdiNaghdi率：率：JTW WJaumannJaumann率：率：0()( )t

4、tttdt 应力更新应力更新: :本构关系本构关系例例4.1 考虑处于剪切状态的一个单元，如图所示。对于次弹性各向同考虑处于剪切状态的一个单元，如图所示。对于次弹性各向同性材料，应用性材料，应用Jaumann，Truesdell和和GreenNaghdi率求出剪切率求出剪切应力应力 x(t) = X + k t Y , y(t) = Y , z(t) = Z :Ds CD2.2.几种客观率几种客观率CCCJTGDssssC 对于不同客观率采用了对于不同客观率采用了相同的材料常数相同的材料常数，其差别是非常大的。，其差别是非常大的。事实上，这是误用了材料模型。材料模型必须根据不同的率转换。事实上

5、，这是误用了材料模型。材料模型必须根据不同的率转换。这是变形体，若是刚体转动这是变形体，若是刚体转动，Jaumann率与率与Truesdell率是一致率是一致的的。例例4.12.2.几种客观率几种客观率Green-Green-NaghdiNaghdi率率JaumannJaumann率率发生有限剪切时发生有限剪切时慎用慎用JaumannJaumann率率切线模量之间的关系切线模量之间的关系 对于各向同性材料对于各向同性材料Jaumann率的切线模量为率的切线模量为 次次弹性本构关系共同应用的形式为弹性本构关系共同应用的形式为DC:JJsDC:TTsDC:GGs)(jkiljlikklijJijk

6、lCs对于同一种材料，切线模量不同，材料反应的率形式不同，如对于同一种材料，切线模量不同，材料反应的率形式不同，如 *CCICCCJJTssJsCTsC如果如果是常数，是常数，不是常数。不是常数。 切线模量切线模量证明见第证明见第5.4.5节，推导复杂节，推导复杂GT*spinCCCC2.2.几种客观率的关系几种客观率的关系如何得到正确的结果？如何得到正确的结果？2.2.几种客观率的关系几种客观率的关系如何得到正确的结果？如何得到正确的结果？ECS:SE :TSETTTT F S FFFCF FFFE: or :TSESEijklimjnkplqmnpqCF F F F CttFFFFCDCD

7、假设假设C CSE已知：已知：假设假设 已知：已知：DsC:():GDsRRRRCD如何处理各向异性材料？:C():GDss RRRRC2.2.几种客观率的关系几种客观率的关系如何得到正确的结果？如何得到正确的结果？如何处理各向异性材料？:C():GDss RRRRC:TSECt FFFFCxyPK2PK2和共轴旋转应力和共轴旋转应力横观各向同性横观各向同性:SESCE:DsCD2.2.几种客观率的关系几种客观率的关系什么时候可以不区分几种客观率？什么时候可以不区分几种客观率？TTsssssvLLGT JTW WAssume F R, ignoring the stretch componen

8、t of FAssume = W退化退化简化简化Comparison of different objective stress rateTruesdellDifficult to implementNot used in commercial softwareGreen-NaghdiKinematically consistent with the rate of Cauchy stressMust accurately determine the rotation tensor RJaumannRelatively easy to implementProduces symmetric t

9、angent moduliGT TTvL LJTW W2.2.几种客观率的关系几种客观率的关系总结一：总结一：2.2.几种客观率的关系几种客观率的关系总结二：总结二：:C():GDss RRRRC:TSECt FFFFCxyu当参考构型和当前构型重合时：当参考构型和当前构型重合时：F=Iu当旋转构型和当前构型重合时：当旋转构型和当前构型重合时：R=I0, , d=dTSECtCS = C, GDssC使使TLTL和和ULUL的互相转换变得非常方便！的互相转换变得非常方便！SolverElement TypeObjective RateAbaqus/StandardSolid (Continuu

10、m)JaumannStructural (Shells, Membranes, Beams, Trusses)Green-NaghdiAbaqus/ExplicitSolid (Continuum)Jaumann&Green-NaghdiStructural (Shells, Membranes, Beams, Trusses)Green-Naghdi2.2.几种客观率的关系几种客观率的关系总结三：总结三：后拉、前后拉、前推推Euler张量张量 Lagrangian张量张量 两点张量两点张量 Green应变张量应变张量E PK2应力应力 线单元线单元dX 线单元线单元dxCauchy应

11、力应力 速度梯度速度梯度 LDW Lagrangian矢量矢量dX和和Eulerian矢量矢量dx定义的二阶张量定义的二阶张量 可以由后拉和前推运算给出可以由后拉和前推运算给出E-LE-L张量之间映射的统一描述。张量之间映射的统一描述。例如，例如，L L矢量矢量dXdX由由F F前推到当前构形给出前推到当前构形给出E E矢量矢量dxdxXXFxddd*E E矢量矢量dxdx由由1F后拉到参考构形给出后拉到参考构形给出L L矢量矢量dXdX xxFXddd*1LagrangianEulerian Eulerian-Lagrangian 前推运算前推运算 后拉运算后拉运算 通过在拓扑空间的分析，通

12、过在拓扑空间的分析，获得各种获得各种张量之间的关系和映射张量之间的关系和映射 3.3.前推后拉及前推后拉及LieLie导数导数后拉、前后拉、前推推 二阶张量的后拉和前推运算给出了在变形和未变形构形情况下二阶张量的后拉和前推运算给出了在变形和未变形构形情况下张量之间的关系，例如张量之间的关系，例如GreenGreen应变率和应变率和变形率的关系，变形率的关系，PK2PK2应力与应力与CauchyCauchy应力的关系。应力的关系。 这些定义取决于是否一个张量是这些定义取决于是否一个张量是动力学动力学还是还是运动学运动学的，区别在的，区别在于由这些张量所观察到于由这些张量所观察到的功的的功的共轭性

13、：如功共轭的运动学和动力学共轭性：如功共轭的运动学和动力学张量被后拉或前推，张量被后拉或前推，则功必须则功必须保持不变。保持不变。 许多关系来自于框许多关系来自于框3.23.2，这些概念，这些概念能够使我们发现那些能够使我们发现那些不容易显不容易显示的关系。一些重要的二阶张量的后拉和前推在框示的关系。一些重要的二阶张量的后拉和前推在框5.165.16给出。给出。ijklTikijTFDFE FDFE或11ljklTikijTFEFD FEFD或前推运算前推运算后拉运算后拉运算TljklikijTFJFSJs11 或FFSTljklikijTFSFJJ11 s或FSF后拉后拉 前推前推 3.3.

14、前推后拉及前推后拉及LieLie导数导数LieLie导数导数 后拉和前推的概念为定义张量的时间导数提供了数学上的一致性后拉和前推的概念为定义张量的时间导数提供了数学上的一致性LieLie导数导数。如框。如框5.175.17，KirchhoffKirchhoff应力的应力的LieLie导数是其应力的后拉的时导数是其应力的后拉的时间导数的前推间导数的前推。) F(F) (1*DtDDtDLv 不严格地说，在不严格地说，在LieLie导数中，导数中，是在固定的参考构形中对时间求导，是在固定的参考构形中对时间求导，再前推到当前构形。再前推到当前构形。 在框在框5.175.17中给出了用势共轭方式定义的

15、运动学张量的中给出了用势共轭方式定义的运动学张量的LieLie导数。导数。 以上计算是将应力后拉到参考构型上，对时间求普通导数，再前推回以上计算是将应力后拉到参考构型上，对时间求普通导数，再前推回到当前构型到当前构型。否则是偏导数。否则是偏导数。 ,tXx 3.3.前推后拉及前推后拉及LieLie导数导数TTTvDtDDtDDtDLFFFFFSF)()(1*后拉、前推和后拉、前推和LieLie导数导数 证明证明KirchhoffKirchhoff应力的对流率对应于它的应力的对流率对应于它的LieLie导数导数TFSFt(框框3.2) 材料时间导数的计算材料时间导数的计算 TTTTTTTTTTv

16、DtDLFFFFFFFFFFFFFFFFFFFF)()()(1111111*111FFFF1FFLTcvLLL应用应用 得到得到 LieLie导数等价于在公式导数等价于在公式(5.4.22)(5.4.22)中定义的中定义的TruesdellTruesdell应力的对流率应力的对流率 ) F(F) (1*DtDDtDLv3.3.前推后拉及前推后拉及LieLie导数导数4.4.常用客观率应力更新流程常用客观率应力更新流程JTW WJaumannJaumann率：率：0()( )ttttdt 应力更新应力更新: :0()( )tJTtttdtW Wxyttt, , FRWDDD速度梯度速度梯度jii

17、jTTxvL grad或vvxvL速度梯度张量可以分解为对称部分和偏对称部分为速度梯度张量可以分解为对称部分和偏对称部分为 2121TTLLLLLijjiijTxvxvD21 21或LLDijjiijTxvxv21 W 21或LLW令令变形率变形率转动率转动率任何一个二阶张量都可以表示为它的对称部分和偏对称部分的和任何一个二阶张量都可以表示为它的对称部分和偏对称部分的和 ijijjiijTWDvL, 或WDvL回顾第回顾第3 3章章 4.4.常用客观率应力更新流程常用客观率应力更新流程4.4.常用客观率应力更新流程常用客观率应力更新流程LS-LS-DynaDyna AnsysAnsys等商业软

18、件采用此等商业软件采用此算法算法0()( )tJTtttdtW Wxyttt , , FRWDDD2tt4.4.常用客观率应力更新流程常用客观率应力更新流程AbaqusAbaqus采采用此算法用此算法DDDDTtttR R+ +=鬃+=鬃+xyttt , , FRWDDD0()( )tJTtttdtW WHughes-Winget近似近似 WFRUD= DD11122RII-骣骣鼢珑D=-D +D鼢珑鼢珑桫桫刚体转动求解刚体转动求解法一：通过QR分解计算法二：通过法二：通过Hughes-Hughes-WingetWinget更新更新，近似求解，近似求解4.4.常用客观率应力更新流程常用客观率应

19、力更新流程JaumannJaumann率率xyttt , , FRWDDD100tttttttttxxxFFFxxx-+ D+ D+ D骣骣抖鼢珑鼢D=珑鼢珑鼢珑抖桫桫()ttFL FtD= DHughes & Winget, 1980 (见网络学堂课外阅读）见网络学堂课外阅读）OR W4.4.常用客观率应力更新流程常用客观率应力更新流程0()( )ttttdt 应力更新应力更新: :GT Green-Green-NaghdiNaghdi率：率：0()( )tGTtttdt 4.4.常用客观率应力更新流程常用客观率应力更新流程0()( )ttttdt 应力更新应力更新: :TDRD R

20、:DsCD把课本知识和实际程序应用结合起来！AbaqusAbaqus，LS-LS-DynaDyna4.4.常用客观率应力更新流程常用客观率应力更新流程11122RII-骣骣鼢珑D=-D +D鼢珑鼢珑桫桫刚体转动求解（商业软件显式程序中采用算法）刚体转动求解（商业软件显式程序中采用算法）J. K. Dienes, On the Analysis of Rotation and Stress Rate In Deforming Bodies, Acta Mechanica, Vol. 32, pp 217-232, 1979.LS-LS-DynaDyna AnsysAnsys等商业软件采用等商业软

21、件采用4.4.常用客观率应力更新流程常用客观率应力更新流程0()( )tJTtttdtW WDDDDTtttR R+ +=鬃+=鬃+11TnnR R+=鬃=鬃小结：小结：xyttt , FRDD0t, FR1 ABAQUS的的UMAT 用户材料子程序用户材料子程序(User-defined Material Mechanical (User-defined Material Mechanical BehaviorBehavior，简称，简称UMAT)UMAT)通过与通过与ABAQUSABAQUS主求解程序的接口实现与主求解程序的接口实现与ABAQUSABAQUS的数据交流。的数据交流。在输入文

22、件中，使用关键字在输入文件中，使用关键字“* *USER USER MATERIALMATERIAL”表示定义用户材料属性。表示定义用户材料属性。uUMATUMAT子程序具有强大的功能，使用子程序具有强大的功能，使用UMATUMAT子程序：子程序：u可以定义材料的本构关系，使用可以定义材料的本构关系，使用ABAQUSABAQUS材料库中没有包含的材料库中没有包含的材料进行计算，扩充程序功能；材料进行计算，扩充程序功能；u几乎可以用于力学行为分析的任何分析过程，可以把用户材几乎可以用于力学行为分析的任何分析过程，可以把用户材料属性赋予料属性赋予ABAQUSABAQUS中的任何单元；中的任何单元；

23、u必须在必须在UMATUMAT中提供材料本构模型的雅可比（中提供材料本构模型的雅可比（JacobianJacobian）矩阵，）矩阵，即应力增量对应变增量的变化率；即应力增量对应变增量的变化率；u可以和用户子程序可以和用户子程序“USDFLDUSDFLD”联合使用，通过联合使用，通过“USDFLDUSDFLD”重重新定义单元每一物质点上传递到新定义单元每一物质点上传递到UMATUMAT中场变量的数值。中场变量的数值。5.5. ABAQUS ABAQUS用户自定义材料用户自定义材料ABAQUSABAQUS用户子程序用户子程序l使用方法使用方法 要在模型中包含用户子程序，可以利用要在模型中包含用户

24、子程序，可以利用ABAQUSABAQUS执行程执行程序，在执行程序中应用序，在执行程序中应用useruser选项指明包含这些子程序的选项指明包含这些子程序的FORTRANFORTRAN源程序或者目标程序的名字。源程序或者目标程序的名字。 ABAQUSABAQUS的输入文件除了可以通过的输入文件除了可以通过ABAQUS/CAEABAQUS/CAE的作业的作业模块提交运行外，还可以在模块提交运行外，还可以在ABAQUS CommandABAQUS Command窗口中输入窗口中输入ABAQUSABAQUS执行程序直接运行：执行程序直接运行：ABAQUS job=ABAQUS job=输入文件名输入

25、文件名 user=user=用户子程序的用户子程序的FortranFortran文文件名件名注注：ABAQUS/StandardABAQUS/Standard和和ABAQUS/ExplicitABAQUS/Explicit都支持用户子都支持用户子程序功能，但是它们所支持的用户子程序种类不尽相同。程序功能，但是它们所支持的用户子程序种类不尽相同。 由于主程序与由于主程序与UMATUMAT之间存在数据传递，甚至共用一些变量，之间存在数据传递，甚至共用一些变量，因此必须遵守有关因此必须遵守有关UMATUMAT的书写格式，的书写格式，UMATUMAT中常用的变量在文件开中常用的变量在文件开头予以定义，

26、通常格式为：头予以定义，通常格式为： SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,DRPLDE,SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,RPL,DDSDDT,DRPLDE, 1 DRPLDT,STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME, 1 DRPLDT,STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME, 2 NDI,NSHR,NTENS

27、,NSTATV,PROPS,NPROPS,COORDS, 2 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROTDROT,PNEWDT,PNEWDT, 3 CELENT, 3 CELENT,DFGRD0DFGRD0, ,DFGRD1DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC),NOEL,NPT,LAYER,KSPT,KSTEP,KINC)C C INCLUDE ABA_PARAM.INC INCLUDE ABA_PARAM.INCC C CHARACTER CHARACTER* *80 CMNAME80 CMNAME D

28、IMENSION STRESS(NTENS),STATEV(NSTATV), DIMENSION STRESS(NTENS),STATEV(NSTATV), 1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS), 1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS), 2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1), 2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1), 3 PRO

29、PS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3) 3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3) user coding to define DDSDDE, STRESS, STATEV, SSE, SPD, SCD user coding to define DDSDDE, STRESS, STATEV, SSE, SPD, SCD and, if necessary, RPL, DDSDDT, DRPLDE, DRPLDT, PNEWDT and, if ne

30、cessary, RPL, DDSDDT, DRPLDE, DRPLDT, PNEWDT RETURN RETURN END ENDl变量介绍变量介绍 STRAN(NTENS)：应变矩阵：应变矩阵 DSTRAN(NTENS)：应变增量矩阵：应变增量矩阵 DTIME：增量步的时间增量：增量步的时间增量 NDI：直接应力分量的个数：直接应力分量的个数 NSHR：剪切应力分量的个数：剪切应力分量的个数 NTENS：总应力分量的个数：总应力分量的个数SSE,SPD,SCD 分别定义每一增量步的弹性应变能，塑性耗散和蠕变分别定义每一增量步的弹性应变能，塑性耗散和蠕变耗散。它们对计算结果没有影响，仅仅作为

31、能量输出。耗散。它们对计算结果没有影响，仅仅作为能量输出。空间坐标系下的分量空间坐标系下的分量 UMAT UMAT中的应力矩阵、应变矩阵以及矩阵中的应力矩阵、应变矩阵以及矩阵 DDSDDEDDSDDE，DDSDDTDDSDDT，DRPLDE DRPLDE 等，都是直接分量存储在前，剪切分量存储在后。直接等，都是直接分量存储在前，剪切分量存储在后。直接分量有分量有 NDINDI个，剪切分量有个，剪切分量有 NSHR NSHR 个。各分量之间的顺序根据单个。各分量之间的顺序根据单元自由度的不同有一些差异，所以编写元自由度的不同有一些差异，所以编写UMATUMAT时要考虑到所使用单时要考虑到所使用单

32、元的类别。元的类别。,DDSDDE NTENS NTENS是一个是一个NTENSNTENS维方阵，称为雅克比矩阵，即维方阵，称为雅克比矩阵，即 ，切线模量。，切线模量。/,DDSDDE I J 是表示增量步结束时第是表示增量步结束时第 J J 个应变分量的改变引起的个应变分量的改变引起的第第 I 个应力分量的变化。通常雅可比是一个对称矩阵。个应力分量的变化。通常雅可比是一个对称矩阵。STRESS NTENS在增量步开始，应力张量矩阵中的数值通过在增量步开始，应力张量矩阵中的数值通过UMATUMAT和主程序之间的和主程序之间的接口传递到接口传递到UMATUMAT中；在增量步结束，中；在增量步结束

33、，UMATUMAT将对应力张量矩阵更新。将对应力张量矩阵更新。对于包含刚体转动的有限应变问题，一个增量步调用对于包含刚体转动的有限应变问题，一个增量步调用UMATUMAT之前就之前就已经对应力张量进行了刚体转动，因此在已经对应力张量进行了刚体转动，因此在UMATUMAT中只需处理应力张中只需处理应力张量的共旋量的共旋( (co-rotational)co-rotational)变形部分。变形部分。UMATUMAT中应力张量的度量为中应力张量的度量为CauchyCauchy应力。应力。DDDDTtttR R+ +=鬃+=鬃+:JCsDD= =用于存储状态变量的矩阵，在增量步开始时将数值传递到用于

34、存储状态变量的矩阵，在增量步开始时将数值传递到UMATUMAT中。中。也可在子程序也可在子程序USDFLDUSDFLD或或UEXPANUEXPAN中先更新数据，然后在增量步开始中先更新数据，然后在增量步开始时将更新后的数据传递到时将更新后的数据传递到UMATUMAT中。在增量步结束时必须更新状态中。在增量步结束时必须更新状态变量矩阵中的数据。变量矩阵中的数据。 和应力张量矩阵不同的是：对于有限应变问题，除了材料本和应力张量矩阵不同的是：对于有限应变问题，除了材料本构行为引起的数据更新以外，构行为引起的数据更新以外，状态变量矩阵状态变量矩阵NSTATEVNSTATEV中的任何矢中的任何矢量或者张

35、量都必须通过旋转来考虑材料的刚体运动。量或者张量都必须通过旋转来考虑材料的刚体运动。STATEV NSTATEVPROPS NPROPS材料常数矩阵，矩阵中元素的数值对应于输入文件关键字材料常数矩阵，矩阵中元素的数值对应于输入文件关键字 * *USER USER MATERIALMATERIAL下面的数据行下面的数据行。DDDDTtttR R+ +=鬃+=鬃+1 ABAQUS的的UMAT2 ABAQUS的的VUMAT5.5. ABAQUS ABAQUS用户自定义材料用户自定义材料ubroutine vumat(1 nblock, ndir, nshr, nstatev, nfieldv, np

36、rops, lanneal, 2 stepTime, totalTime, dt, cmname, coordMp, charLength, 3 props, density, strainInc, relSpinInc, 4 tempOld, stretchOld, defgradOld, fieldOld, 5 stressOld, stateOld, enerInternOld, enerInelasOld, 6 tempNew, stretchNew, defgradNew, fieldNew, 7 stressNew, stateNew, enerInternNew, enerIne

37、lasNew ) dimension props(nprops), density(nblock), coordMp(nblock,*), 1 charLength(nblock), strainInc(nblock,ndir+nshr), 2 relSpinInc(nblock,nshr), tempOld(nblock), 3 stretchOld(nblock,ndir+nshr), 4 defgradOld(nblock,ndir+nshr+nshr), 5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr), 6 stateOl

38、d(nblock,nstatev), enerInternOld(nblock), 7 enerInelasOld(nblock), tempNew(nblock), 8 stretchNew(nblock,ndir+nshr), 8 defgradNew(nblock,ndir+nshr+nshr), 9 fieldNew(nblock,nfieldv), 1 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev), 2 enerInternNew(nblock), enerInelasNew(nblock), 2 ABAQUS的的VUMA

39、TVariables to be definedstressNew (nblock, ndir+nshr)Stress tensor at each material point at the end of the increment.stateNew (nblock, nstatev)State variables at each material point at the end of the increment.Variables passed in for informationStrainInc (nblock, ndir+nshr)Strain increment tensor at each material point.