作者介绍
长江学渣
博士
浙江大学在读博士，已在国际知名SCI期刊上发表论文多篇，并在美国加州
理工大学、新加坡南洋理工大学等世界顶级高校举办的国际会议上报告，拥有数年UMAT编写，本构理论研究及
有限元数值模拟经验。
本科毕业于华中科技大学土木工程专业，随后进入浙江大学岩土工程专业直接攻读博士学位。
现已通过国家留学基金委评审，即将公派留学至美国西北大学进行交流学习。
一．Abaqus中UMAT的使用背景
ABAQUS中虽然提供了很多材料本构模型，但是有时候我们的分析的实际情况是比较复杂的，比如当我们要计算一种特定的橡胶材料的时候，已有的橡胶本构比如Mooney等都有着这样、那样的不足，很难满足我们的要求。
又比如我们要计算土壤应力的时候，土壤本身受压不受拉，有着很强的非线性因素，这时候就需要我们建立足够精确的土壤本构模型。如果你正在处理的问题，ABAQUS的常规本构不能满足你的要求，那么UMAT将是你的首选。
(a)试验滞回曲线(b)有限元模拟滞回曲线
(c)锈蚀前混凝土应力分布图(d)锈蚀后混凝土应力分布图
图1.UMAT应用实例一：将锈蚀钢筋本构通过UMAT编入ABAQUS分析不均匀锈蚀钢筋混凝土剪力墙力学性能
(a)水平沉积地基破坏形态及承载力
(b)竖向沉积地基破坏形态及承载力 图2.UMAT应用实例二：将能够考虑各向异性的土体本构通过UMAT编入ABAQUS从而分析各向异性对地基承载力及破坏形态的影响 二．Abaqus中UMAT编写需要面对的主要难点问题 难点一：本构模型的高度非线性
由于许多本构模型都是高度非线性的，在对这些模型计算得到的增量进行积分时，想要取得准确的积分结果是十分困难的，因此所采用的积分方法是十分关键的。以最简单的欧拉法为例，已知t时刻的应力σt，求t+△t时刻的应力σt+△t，采用的方法是：
(1)计算t时刻的刚度矩阵D(σt);
(2)通过求得的刚度矩阵计算应力增量dσ= D(σt)×dε；
(3)通过计算得到的应力增量更新应力，得到t+△t时刻的应力σt+△t=σt+dσ=σt+ D(σt)×dε。
然而对于高度非线性的本构模型，由于在每个增量步内由t到t+△t的时间里其刚度矩阵D的变化是非常大的，因此如果像欧拉法一样将初始t时刻的刚度矩阵作为整个过程的刚度矩阵，必然会产生较大的积分误差。
在有限元计算中，每个增量步的误差如果不能控制在合理的范围内，那这一误差就将在计算过程中逐渐累积，导致严重失真的模拟结果。下图比较了不同积分方法计算得到的应力-应变曲线，实线表示精确值。从图中可以看出，若不采用具有高阶精度的积分方法，计算得到的结果会严重偏离真实值；而随着积分方法精度的提高，得到的结果也会逐渐收敛于精确值。
图3 不同积分方法得到的应力应变曲线 然而，即使采用高阶的积分方法，在面对非线性程度较高的本构模型以及积分过程中误差的累计效应时，都会显得力不从心。此外，积分精度的提高必然会增大积分方法的复杂性，从而导致计算效率的下降，因此如何在积分精度和计算效率之间平衡，又是UMAT中本构积分的一个需要考虑的问题。 难点二：应力穿越屈服面
众所周知，当应力在屈服面内移动时，材料的响应是纯弹性的；而当应力由屈服面上移动到屈服面外时，屈服面要随之发生扩张，同时材料的响应是弹塑性的。所以当应力在一个增量步内由屈服面内穿越到屈服面外时，在屈服面内部分需要用弹性刚度矩阵更新，而屈服面外部分需要用弹塑性刚度矩阵更新，如图4(a)所示。
如何准确计算这两部分各占多少比例，是需要十分注意的，因为如果计算误差过大，同样将导致有限元模拟结果失真。同时，要注意应力穿越屈服面有一些不同的形式，图4(a)表示的是一种最一般的穿越；但同时图4(b)也是很可能出现的一种穿越情况，即应力穿越进入屈服面再穿越出屈服面。因此一个完备的UMAT中，需要先判断应力出现了哪种形式的穿越，再用合理的方法去计算屈服面内外部分所占的比例。
(a)
(b) 图4 应力穿越屈服面 难点三：由于积分误差导致的应力偏离屈服面
当应力穿越屈服面时，屈服面也会随之扩大，因此更新后的应力理论上是应该处在屈服面上的。但是在实际计算中，由于积分误差，往往会导致更新后的应力偏离屈服面。如果这时不对偏离的应力进行修正，有可能因为误差累积，导致应力严重偏离屈服面，即超出屈服面很多。因此，在每次积分后，都要检测更新后的应力是否处于屈服面上，若其偏离屈服面，则要对其进行修正。图5展示了模拟滑坡问题时，在UMAT中不将应力修正回屈服面(图5(a))以及将应力修正回屈服面(图5(b))的结果。
从图5(a)中可以看出，若不将应力修正回屈服面，在计算过程中由于误差累积应力将越来越偏离屈服面，并最终严重超出屈服面，导致严重失真的模拟结果。而从图5(b)中可以看出，采用应力拉回屈服面处理后，每次更新后的应力都会落在屈服面上，最终模拟出比较好的滑坡结果。
因此，将应力拉回屈服面是一个完备的UMAT必须拥有的步骤。但是，如何将应力拉回屈服面，又是一个值得研究的问题，不合理的拉回方式同样会造成计算结果的失真。
(a)不将应力拉回屈服面的模拟结果(上：位移云图；下：应变云图)
(b)将应力拉回屈服面的模拟结果(上：位移云图；下：应变云图)
图5 模拟滑坡问题时不将应力拉回屈服面以及将应力拉回屈服面模拟结果对比
三．Abaqus中基于UMAT的本构二次开发教程由入门到精通：
通过详细的教程讲解编写一个UMAT的完整步骤，解决上文所说的几个UMAT编写难题。
本课程将由浅入深介绍UMAT的编写方法，包含了一些数值处理技巧如自动误差控制积分方法，弹性试探应力，应力拉回屈服面等。同时附件还会有课程里的课件，参考资料，论文，UMAT代码等。 通过本课程你可以学到：
(1)了解UMAT的基本使用方法，书写格式； (2)掌握如何利用高阶精度的积分方法对高度非线性的本构方程进行积分； (3)掌握弹塑性模型的积分步骤，例如弹性试探应力及其位置的判断； (4)掌握如何对超出屈服面的应力进行修正； (5)通过对UMAT实例细致深入的解释使读者掌握如何将这些数值方法编写成代码； (6)介绍一些UMAT编写的技巧及注意事项； 课程目录：
(1)UMAT基本书写格式及简单例子
(2)自动控制误差的应力显式积分法
(3)弹塑性模型积分方法_理论(弹性试探应力，应力回拉)
(4)UMAT实例详解(上)：弹性试探应力及其位置的判断
(5)UMAT实例详解(中)：应变增量在屈服面内外的比例计算
(6)UMAT实例详解(下)：高度非线性本构方程的积分方法
(7)UMAT实例详解(补充)：如何将应力拉回屈服面
(8)UMAT编写的一些技巧及注意事项
