-
前
言
有限元
法是工程中广泛使用的一种数值计算方法。它是力学、计算方法和计算机
技术相结合的产
物。在工程应用中，有限元法比其它数值分析方法更流行的一个重要
原因在于：相对与其
它数值分析方法，有限元法对边界的模拟更灵活，近似程度更高。
所以，伴随着有限元理
论以及计算机技术的发展，大有限元软件的应用证变得越来越
普及。
ABAQUS
软件一直以非线性有限元分析软件而闻名，
这也是它和
ANSYS,Nastran
等
软件的区别所在。非线性有限元分析的用处越来越大，因为在所用材料非常复杂很多
情况下，用线性分析来近似已不再有效。比方说，一个复合材料就不能用传统的线性
分析
软件包进行分析。任何与时间有关联，有较大位移量的情况都不能用线性分析法
来处理。
多年前，虽然非线性分析能更适合、更准确的处理问题，但是由于当时计算
设备的能力不
够强大、非线性分析软件包线性分析功能不够健全，所以通常采用线性
处理的方法。
这种情况已经得到了极大的改善，计算设备的能力变得更加强大
、类似
ABAQUS
这
样的产品功能日
臻完善，应用日益广泛。
非线性有限元分析在各个制造行业得
到了广泛应用，有不少大型用户。航空航天
业一直是非线性有限元分析的大客户，一个重
要原因是大量使用复合材料。新一代波
音
787
客机将全部采用复合材料。只有像
ABAQUS
这样的软件，才能分析包括多个子
系统的产品耐久性能。在汽车业，用
线性有限元分析来做四轮耐久性分析不可能得到
足够准确的结果。分析汽车的整体和各个
子系统的性能要求（如悬挂系统等）需要进
行非线性分析。
在土
木工程业，
ABAQUS
能处理包括混凝土静动力开裂分析以
及沥青混
凝土方面的静动力分析，还能处理高度复杂非线性材料的损伤和断裂问题，这对
于大
型桥梁结构，高层建筑的结构分析非常有效。
瞬态、大变形、高级材料的碰撞问题必须用非线性有限元分析
来计算。线性分析
在这种情况下是不适用的。以往有一些专门的软件来分析碰撞问题，但
现在
ABAQUS
在
通用有限元软件包
就能解决这些问题。所以，
ABAQUS
可以在一个软件完成线
性和非线
性分析。
ABAQUS<
/p>
给用户提供了强大二次开发接口，尤其是在材料本构方面，给用户开发符
< br>合实际工程的材料本构模型提供了强大帮助，
本文将针对其用户材料子程序展开研
究，
总结常用材料模型的开发方法。
目
录
摘
要
..........................
..................................................
..................................................
...........
I
ABSTRACT ..............................
..................................................
............................................
II
1.
绪论
.........................
..................................................
..................................................
......
1
1.1.
课题的研究背景
.
< br>............................................... .................................................. ...
1
1.2.
本文的研究内容和方法
.
............................................ ............................................
2
2.
<
/p>
基于
ABAQUS
软件的二次开发
.
..................
..................................................
................
3
2.1.
ABAQUS
介绍
........
..................................................
.............................................
3
2.2.
ABAQUS
各模块简介
.....
..................................................
....................................
3
2.3.
ABAQUS
的二次开发平台
.
..............................
..................................................
...
5
2.4.
ABAQUS
的二次开发语言
.
..............................
..................................................
...
6
3.
用户材料子程序
UMAT
.....
..................................................
............................................
8
3.1.
UMAT
开发环境设置
......
..................................................
....................................
8
3.2.
UMAT
注意事项
........
..................................................
..........................................
9
3.3.
UMAT
接口的原理
..............................................
.................................................
10
3.4.
UMAT
的使用方法
..............................................
.................................................
12
4.
材料非线性问题
....................
..................................................
.......................................
14
4.1.
材料的弹塑性本构关系
.
............................................ ..........................................
14
4.2.
非线性有限元算法理论
.
............................................ ..........................................
17
4.3.
增量理论常刚度法公式推导
.
..........................................
....................................
20
4.4.
增量理论切线刚度法公式推导
.
.........................................
.................................
21
5.
UMAT
程序设计和编码
............................................ .................................................. ...
25
5.1.
本构关系描述
.
................................................ .................................................. ....
25
5.2.
常刚度法程序设计
.
..............................................
................................................
27
5.3.
常刚度法程序编码
.
..............................................
................................................
29
5.4.
切线刚度法程序设计
.
.............................................
.............................................
32
5.5.
切线刚度法程序编码
.
.............................................
.............................................
36
5.6.
程序的调试
.
.................................................
..................................................
.......
39
6.
程序验证
.......................
..................................................
................................................
40
1
6.1.
问题描述
.
..................................................
..................................................
..........
41
6.2.
本构关系
.
..................................................
..................................................
..........
42
6.3.
ABAQUS
自带材料模型计算
..
..................................................
.........................
42
6.4.
常刚度法的
UMAT
验证
.................
..................................................
..................
44
6.5.
切线刚度法的
UMAT
验证
................
..................................................
...............
46
6.6.
两种算法的比较分析
.
.............................................
.............................................
48
7.
结论与展望
......................
..................................................
.............................................
52
7.1.
结论
.
..
..................................................
..................................................
................
52
7.2.
展望
.
..
..................................................
..................................................
................
52
致
谢
..
..................................................
..................................................
................................
54
参考文献
.......................
..................................................
..................................................
.....
55
附
1
：
ABAQUS
自带弹塑性
材料验证的
INP
文件
......
..................................................
...
56
附
2
：用于算法验证的
INP
文件<
/p>
.
................
..................................................
......................
62
2
摘
要
ABAQUS
软件功能强大，特别是
能够模拟复杂的非线性问题，它包括了多种材料
本构关系及失效准则模型，并具有良好的
开放性，提供了若干个用户子程序接口，允
许用户以代码的形式来扩展主程序的功能。<
/p>
本文主要研究了
ABAQUS
用户子程序
UMAT
的开发方法，采用
FORTRAN
语言
编制了各向同性硬化材料
模型的接口程序，研究该类材料的弹塑性本构关系极其实现
方法。
本文紧紧围绕
UMAT
的二次开发技术
，首先对其接口原理做了详细介绍，然后针
对非线性有限元增量理论中的常刚度法和切线
刚度法的算法理论做了深入的剖析，推
导出了常刚度法和切线刚度法的算法理论的具体表
达式，然后分别编制了两种算法的
UMAT
程序，
最后建立了一个具体的验算模型，
通过与
ABAQU
S
自带弹塑性本构关系
的计算结果相比较，验证两者的正确性。
本文还对常刚度法和切线刚度法得算法效率做了对比，
得出了在非线
性程度较高时
切线刚度法效率高于常刚度法的结论。
关键字：
ABAQUS
、
UMAT
、有限元、材料非线性、
FORTRAN
、切线刚度
I
ABSTRACT
ABAQUS
software
powerful,
especially
to
simulate
complex
non-linear
problem,
which
includes
a
wide
range
of
material
constitutive
model
and
failure
criteria,
and
has
a
good open, providing a
number of user subroutine interface that allows
users to code form to
expand the
functions of the main program.
This paper studies the user
subroutine UMAT of ABAQUS development methods, the
use
of
FORTRAN
language
isotropic
hardening
material
model
of
the
interface
program,
studied the effects of such material is
extremely elastic-plastic constitutive relation
method.
This article UMAT tightly around the
secondary development of technology, the first
principle
of
its
interface
detail,
and
then
for
the
theory
of
nonlinear
finite
element
incremental
stiffness of the regular tangent stiffness method
and the theory of algorithms to
do an
in-depth analysis of deduced a regular tangent
stiffness and rigidity of the law of the
specific
expression
of
algorithm
theory,
and
then
the
preparation
of
the
two
algorithms,
respectively,
of
the
UMAT
program,
and
finally
the
establishment
of
a
specific
model
checking, bringing with ABAQUS elasto-
plastic constitutive relation of the calculated
results
compared to verify the
correctness of the two.
This article
also often stiffness and tangent stiffness method
was to do a comparison of
algorithm
efficiency is obtained when a higher degree in the
non-linear tangent stiffness
method
more efficient than the conclusions of law often
stiffness.
KEY
WORDS
：
ABAQUS
、
UMAT
、
Finit
e
element
、
Material
nonlinearity
、
FORT
RAN
、
Tangent
stiffness
II
1.
绪论
1.1.
课题的研究背景
有限单元法基本思想的提出，可以追溯到克劳夫
()
在
1943
年的工作
[1]<
/p>
，
他第一次尝试应用定义在三角形区域上的分片连续函数和最小位
能原理相结合，来求
解
St. Venant
< br>扭转问题。
1960
年克劳夫进一步处理了平面弹性问题
，
并第一次提出了
“有
限单元法”的名
称，使人们开始认识了有限单元法的功效。
四十多
年来，随着电子计算机的广泛应用和发展，有限单元法的理论和应用都得
到迅速的，持续
不断的发展，其应用己由弹性力学平面问题扩展到空间问题、板壳问
题，由静力学问题扩
展到稳定问题、动力问题和波动问题。分析的对象从弹性材料扩
展到塑性、粘弹性、粘塑
性和复合材料等，从固体力学扩展到流体力学、传热学等连
续介质力学领域。在工程分析
中的作用已从分析和校核扩展到优化设计并和计算机辅
助设计。
利用有限元软件解决工程和科学问题，是有限元理论应用于工程设计和科学研究
实践的主要形式。由于工程设计的巨大市场需要，有限元软件的发展是很迅速的，目
< br>前常用的大型有限元软件常见的有
Sap2000
，
ADINA
，
MSC/NASTRAN
，
MSC Marc
，
ANSYS
，
ABAQUS
等，这
些软件的共同特点是具有丰富的单元库和求解器，强大而可
靠的分析功能，
人们利用这些软件解决了很多工程建设和工业产品设计中遇到的问题，
取得了
巨大的经济技术效益。
由于工程问题的千差万别，不同的
用户有不同的专业背景和发展方向，通用软件
不免在具体的专业方面有所欠缺，针对这些
不足，大部分的通用软件都提供了二次开
发功能，以帮助用户减少重复性的编程工作、提
高开发起点、缩短研发周期、降低开
发成本，并能简化后期维护工作，给用户带来很多方
便。基于通用软件平台进行开发，
是目前研究的一个重要发展方向。
ABAQUS
也提供了若干用户子程序（
User Subroutines
）接口，它是一个功能非
1
常强
大且适用的分析工具，与命令行的程序格式相比，用户子程序的限制少得多，
从而使用更加灵活方便。针对
ABAQUS
所提
供的本构关系模型种类有限，无法满足
工程应用需要的问题，
用户子程序中的用户材料子程序（
User-defined Materia
Mechanical Behavior
，简称
UMAT
）接口可以帮助用户定义自己的材料本构模型和算
< br>
法，这是
ABAQUS
的独到
之处。由于其操作方便，能被灵活地应用于各个领域中，
尤其受到用户的青睐。
1.2.
本文的研究内容和方法
ABA
QUS
中用户材料子程序
UMAT
的开
发主要解决两方面的问题：
本构模型的建
立和积分算法的选择。
本文主要研究非线性材料的<
/p>
UMAT
实现方法，并重点研究其迭代算法部分，目前，
用户材料子程序
UMAT
的迭代算法主要是常刚
度法，
常刚度法的优点在于算法原理较简
单，程序编写较方便，
缺点是当遇到复杂非线性材料时，其迭代次数较多，收敛速度
也较慢，在这个情况下，本
文采取的是一种迭代次数较少且收敛速度较快的切线刚度
法，具体就是采用
FORTRAN
语言编制了基于
V
< br>on-Mises
模型的接口程序
,
并采用切线
刚度算法，通过与
ABAQUS
自带本构关系计算的结果相比较，验证其正确性。
本文的
研究工作紧紧围绕
UMAT
的二次开发技术，首先根据有限元方
法推导材料
非线性问题算法的公式，然后参考
UMAT
接口规范设计程序的算法流程，继而编写出
该程序，最后建立一个具体的
本构和具体的模型做测试，验证程序的正确性，在这一
过程中，调试是一个非常重要的过
程，占用了大量的时间，在调试程序时采用了将中
间变量输出到文本的方式，这样能明确
跟进迭代过程，发现算法或程序的缺陷。
本文采用的本构关系是经过归纳和抽象的，
也就是说本文的程序并不仅仅是只
针对
某个具体模型和问题，而是针对所有符合抽象出的各向同性硬化材料，这样做的好处
是能保证程序的通用性和复用性，避免以后的重复劳动，当然，这也是符合
ABAQUS
软件设计
UMAT
接口的宗旨的。
2
2.
基于
A
BAQUS
软件的二次开发
2.1. ABAQUS
介绍
ABAQUS
是一套功能强大的基于有限元法的工程模拟软件
[2]
，
其解决问题的
范围
从相对简单的线性分析到最富有挑战性的非线性模拟问
题。
ABAQUS
具备十分丰富的、
可
模拟任意实际形状的单元库。并与之对应拥有各种类型的材料模型库，可以模拟大
多数典
型工程材料的性能，其中包括金属、橡胶、高分子材料、复合材料、钢筋混凝
土、
可压缩弹性的泡沫材料以及岩石和土这样的地质材料。
作为通用的模拟分析工具，
ABAQUS
不仅能解决结构分析中的问题，
还能模拟和研究各种领域中的问题，
如热传
导、
质量扩散、
电子元
器件的热控制
（热一电耦合分析）
、
声
学分析、
土壤力学分析
（渗
流——应力
耦合分析）和压电介质力学分析。
ABA
QUS
为用户提供了广泛的功能，且使用起来又十分简明。最复杂的问题也可
以很容易地建立模型
[3]
。
例如复杂的多部件问题可以通过对每个部件定义材料模型和几
何形状，然后再把它们
组装起来而构成。在大部分模拟分析问题中，甚至在高度非线
性问题中，用户也只需要提
供结构的几何形状、材料性能、边界条件和荷载工况这样
的工程数据就可以进行分析。在
非线性分析中，
ABAQUS
能自动选择合适的荷载增量
和收敛精度。不仅能选择这些参数值，而且能在分析过程中不断地调整参数来保证有
效地得到高精度的解，很少需用户去定义这些参数。
2.2.
ABAQUS
各模块简介
ABAQUS
有两个主要的分析模块：
ABAQUS/Standard
和
ABAQUS/Explicit
。
ABAQUS/Standard
还有两个特殊用途的附加分
析模块：
ABAQUS/Aqua
和
A
BAQUS/Design
。
另外，还有
ABAQUS
分别与
ADAMS/Flex
，
C-MOLD
和
Mold
flow
的接口模块：
ABAQUS/ADAMS
，
ABAQUS/C-MOLD
和
< br>ABAQUS/ MOLDFLOW
。
ABAQUS/C
AE
是完全的
ABAQUS
3
工作环境模块，它包括了
ABAQUS
模型的构造，交互式提交作业、监控作业过程以及评
价结果的能力。
ABAQUS/Viewer
是
ABAQUS/CAE
的子集，它具有后处理功能，这些模
块之间的关系见图
2- 1
图
2-1
ABAQUS/Standard
ABAQUS/Standard
是一个通用分析模块，在数值方法上采用有限元方法常
用的
隐式积分。它能够求解广泛的线性和非线性问题，包括结
构的静态、动态问题、热
力学场和电磁场问题等。对于通常同
时发生作用的几何、材料和接触非线性可以采
用自动控制技术处理，也可以由用户自己控制。
ABAQUS/Explicit
ABAQUS/Explicit
是一个在数值方法上采
用有限元显式积分的特殊模块，
它利用对
时间的显式积分求解动
态有限元方程。它适合于分析诸如冲击和爆炸这样短暂、瞬时
的动态问题，同时对高度非
线性问题如模拟加工成型过程中接触条件的改变等也非常
有效。
ABAQUS/CAE
AB
AQUS/CAE
是
ABAQUS
进行
有限元分析的前后处理模块，
也是建模、
分析和后处理
的人机交互平台。该模块根据结构的几何图形生成网格，将材料和截面的特性分配到
网格上，并施加载荷和边界条件。该模块可以进一步将生成的模型投入到分析模块中
< br>进行高效率的后台运行，
并对运行情况进行监测，
对计算
结果进行后处理。
ABAQUS/CAE
的后处理支持
ABAQUS
分析模块的所有功能，
并且对计算
结果的描述和解释提供了范围
4
很广的选择，除了通常的云图，等值线和动画显示之外，还可以用列表，曲线
(
包括部
分常用运算
)<
/p>
等其他常用工具来完成对结果数据的处理。该模块的许多独特功能与特
点，例如
CAD
特征化建模、参数化建模、适应设计者要求
的数据管理系统等极大的方
便了
ABAQUS
< br>的使用者。
ABAQUS/Aqua
ABAQUS/Aqua
的
一系列功能可以附加在
ABAQUS/Standard
中应用
。
它偏向于模拟海
上结构，如海洋石油平台。它的功能包括模拟
波浪，风载荷及浮力的
影响。在本指南
中不讨论
ABAQUS/Aqua
。
ABAQUS/ADAMS
ABAQUS/ADAMS
允许
ABAQUS
< br>有限元模型作为柔性部件进入到
MDIADAMS
产品族
中去
进行分析。
ABAQUS/C-MOLD
ABAQUS/C-MOLD
把注模分析软件
C-MOLD
中有限元网格、材料性质和初始应力数据
转换成为
< br>ABAQUS
输入文件。
ABAQUS/Design
ABAQUS/Design
的一系列功能可附加在
ABAQUS/Standard
中进行设计敏度计算。
ABAQUS/MOLDFLOW
ABAQUS/MOLDFLOW
模块把
MOLDFLOW
分析软件中的有限元模型信息
转换成
ABAQUVS
输入文件的一部分。
2.3.
A
BAQUS
的二次开发平台
ABAQUS
的脚本语言接口非
常友好，其自嵌的脚本语言是
Python
[4]
，系国际上广泛
使用、功能强大、具有良好开放性的一种面向对象程序设计语
言。所以，应用
Python
在
ABA
QUS
中进行二次开发也比较方便，且可移植性强。
ABAQU
S
以基于
Python
的语法
规则向二次开发者提供了许多库函数，这些库函数主要是用来增强
ABA
QUS
的交互式
(GUI)
操作功能。
用户可以通过
ABAQUS
的交互式
(
GUI)
界面实现分析对象的特征造型、
指定材料属性、完成网
格剖分和控制、提交并监控分析作业，也可以使用
ABAQUS
脚本
5
语言越过
ABAQUS
的交互式
(GU
I)
界面直接高效地向
ABAQUS
内
核提交任务。
使用
Python
可以进
行参数化建模，修改交互式建立的模型，还可以一次提交多个作业。
出了脚本语言接口，
ABAQ
US
还为用户提供了功能强大的用户子程序接口（
Abaqus
User Subroutines
）
，以帮助用户开发基于
ABAQUS
内核的程序，常用的用户
子程序包
括
UEL(User subroutine to
define an element
，
用户单元子程序
)
，
UMAT(User
subroutine
to define a material's
mechanical
behavior
，
用户材料子程序
)
[5]
，
其
中
UMAT
的使用最为
广泛，它主要用
于用户开发自己的材料模型，以弥补
ABAQUS
自带材料模型
的不足，
帮助用户完成各种材料分析，功能极为强大。
在国外，众多的有限元分析和研究者热衷于使用
ABAQUS
，一个很重要的原因就在
于
ABAQUS
给用户提供了功能强大，使用方便的二次开发工具和接
口，使得用户可以方
便的进行富含个性化的有限元建模、分析和后处理，满足特定工程问
题的需要。通过
用户材料子程序接口，用户可定义任何补充的材料模型，不但任意数量的
材料常数都
可以作为资料被读取，
而且
ABAQUS
对于任何数量的与解相关的状态变量在每一材料计
数点都提供了存储功能，以便在这些子程序中应用。
2.4.
A
BAQUS
的二次开发语言
ABAQUS
的二次开发语言主要有
3
种：
P
ython
，
< br>FORTRAN
，
C++
Py
thon
语言主要用于
GUI
开发，<
/p>
FORTRAN
语言主要用于用户子程序开发，而
c++
语
言主要专注于其他高级开发部分。
本文主要是针对用户子程序的开发，所
以采用
FORTRAN
语言，下面简要介绍一
< br>下该语言极其特点：
FORTRAN
语言是世界上第一个被正式推广使用的高级语言
[6]
。它是
1954
年被提
出来的，
1956
年开始正式使用，至今已有三
十多年的历史，但仍历久不衰，它始终是
数值计算领域所使用的主要语言。
FORTRAN
语言是
Formula
Translation
的缩写，意为“公式翻译”
。它是为科学、
工程问题或企事业管理中的那些能够用数学公式表达的问题
而设计的，其数值计算的
功能较强。
FORTRAN
语言问世以来，根据需要几经发展，先后推出了不同
的版本，主要版
本有
FORTRAN
77
，
FORTRAN
90
，
FORTRAN
95
，
ABAQUS
采用
FORTRAN
77
，通
常用固定格式编写代码。
6
FO
RTRAN
77
语言同
C
语言一样，是一种结构化编程语言
结构化程序设计
方法规定，在结构化的程序中，只能有三种基本结构：
(1)
顺序结构
这是一种最简单的基本结构形式，它的特点是，在这个结构内的各个功能模块或
语句
序列，是按其出现的先后顺序执行的，如赋值语句、输入
/
输出
语句等。它有一个
入口和一个出口，并在入口和出口之间包含着若干个功能块，其中每一
个功能块可以
是一个非转移语句。因此，顺序基本结构块是由一系列的顺序执行语句组成
的。
(2)
分支选择结构
在给定的条件下，分支选择结构判断选择哪一条路径执行，不同路径完成的功能
是不同的。实现分支选择结构主要由块
IF
语句、
ELSE
语句、
END
IF
语句以及
ELSE
IF
语句组成的
IF-THEN-
ELSE
结构。
(3)
循环结构
循环结构也称重复处理结构，即重复执行某一功能块，直到满足（或不满足）某
一条
件为止。实现循环结构的
FORTRAN90
语句主要是
DO
语句、块
IF
语
句和逻辑
IF
语
句的结合。
以上三种基本结构，是组成结构化程序的基本结构形式。这里有两层意思
：一是
结构化的程序中，各个模块均由这三种基本结构组成；二是结构化程序本身，从宏
观
上也是这三种基本结构形式之一。
7
3.
用户材料子程序
UMAT
3.1. UMAT
开发环境设置
由于
UMAT
是采用
FORTRAN
语言编写，
那么要运行
UMAT
就需要安装
FORTRAN
的开发环境，
同时还需要
ABAQU
S
的支持，
本文采用的
ABAQUS<
/p>
版本为
6.81
，
支持
INTEL
Fortran9.1-10.1
，
Intel
Fortran
安装时又需要安装
Microsoft
Visual Studio
的相应版
本，
经过比较，
本文选用
ABAQUS6.81+Intel
Fortran10.1+Microsoft VisualC++ 2005
，
相
对于
ABAQUS
来说，
UMAT
开发环境的设置较为繁琐，
< br>这给子程序的使用带来诸多不
便，为了解决这一问题，我用
C#
语言编制了
ABAQUS
子程序
编译环境设置工具，只
需要将安装文件解压到
ABAQUS
的安装目录，
运行安装程序就可以了，
整个
过程不需
要人工干预，也不需要安装庞大的
VisualC++
2005
，如图
3-1
所示
图
3-1
8
3.2. UMAT
注意事项
ABAQUS
的用户子程序是根据<
/p>
ABAQUS
提供的相应接口，按照
Fo
rtran
语法，
用户自己编写的代
码。它是一个独立的程序单元，可以独立的被存储和编译，也能
被其它程序单元引用，因此，利用它可带回大量数据供引用程序使用，也可以用它
来完成各种特殊的功能。它的一般结构形式是：
SUBROUTINE
S(x1,x2,
……
,xn)
IN
CLUDE
‘
ABA_
’
（用于
ABAQUS/Standard
用户子程序
中）
OR INCLUDE
‘
V
ABA_
’
)
（用于
ABAQUS/Explicit
用
户子程序中）
……
RETURN
END
x1
，
x2
，……，
xn
是
ABAQUS
提供的用户子程序的接口参数，有些参数是
ABAQUS
传到用
户子程序中的，例如
SUBROUTINE DLOAD
中的<
/p>
KSTEP
、
KINC
< br>、
COORDS
，
有些是需要用
户自己定义的，例如
F
，文件
aba_
和
vaba_
随着
ABAQUS
软件的安装而包含在操作系统中，它们含有重要的参数，帮助
ABAQUS
主
求解程序对用户子程序进行编译
和链接。
当控制遇到
RETURN
语句
时便返回到引用程
序单元中去，
END
语句是用户子程序结束的标志。
在一个算例中，
< br>用户可以用到多个用户子程序，
但必须把它们放在一个以
.for
为扩
展名的文件中。运行带有用户子程序的算例同时有
两种方法：一是在
CAE
中运行，在
E
DIT JOB
菜单中的
GENERAL
子菜单的
USER SUBROUTINE FILE
对话框
中选择用户
子程序所在的文件即可；另外是在
D
中运行，语法如下：
abaqus
job=job-name user={source-file|object-
file}
编制用户子程序时应注意
(1)
用户子程序相互之间不能调用，但可以调用用户自己编写
的
Fortran
子程序和
ABAQU
S
应用程序。
ABAQUS
应用程序必
须由用户子程序调用。当用户编写
Fortran
子程序时，建
议子程序名以
K
开头，以免和
ABAQ
US
内部程序冲突。
(2)
当用户在用户子程序中利用
OPEN
打开外部文
件时，要注意以下两点：一是设备号
的选择是有限制的，只能取
15
～
18
和大于
100
的设备号，其余的都已被
ABAQUS
占
用；二是用户需提供外部文件的绝对路径而不是相对路径。
9
< br>(3)
对于不同的用户子程序
ABAQUS
调用的时间是不同的，有的是在每个
STEP
的开
始，有的是
STEP
的结尾，有的是在每个
INCREMENT
的开始等等。当
ABAQUS
调用用户子程序时，都会把当前的
STEP
和
INCREMENT
利用用户子程序的两个实
< br>
参
KSTEP
和
KINC
传给用户子程序，用户可把它们输出到外部文件中，这样就可清
楚的知道
ABAQUS
< br>何时调用该用户子程序。
为保证用户子程序的正确执行
，子程序的书写必须遵循
ABAQUS
的相关规定，
下面以用户材料子程序为例详细说明。
3.3.
U
MAT
接口的原理
用户材料子程序（
User-defined
Material Mechanical Behavior
，简称
< br>UMAT
）是
ABAQUS<
/p>
提供给用户定义自己的材料属性的
Fortran
程序接口
[7][8]
，使用户能使用
ABAQUS
材料库中没有定义的材料模型。用户材
料子程序
UMAT
通过与
ABAQUS
主
求解程序的接口实现与
ABAQUS
的资料交流。在输入文件中，使用关键词“
*USER
MATERIAL
”表示定义用户材料属性。
UMAT
子程序具有强大的功能，使用
UMAT
子程序：
（
1
）可以定义材料的本构关系，使用
ABAQUS
材料库中没有包含的材料进行
计算，扩充程序功能。
（
2
）几乎可以用于力学行为分析的任何分析过程，几乎可以把用户材料属性
赋
予
ABAQUS
< br>中的任何单元。
（
3
）必须在
UMAT
中提供材料本构的雅可比（<
/p>
Jacobian
）矩阵，即应力增量对
应变增量的变化率。
由于主程序与<
/p>
UMAT
之间存在数据传递，甚至共享一些变量，因此必须遵守有
关
UMAT
的书写格式，
UMAT
中常用的变量在文件开头予以定义，通常
格式
SUBROUTINE UMA
T(STRESS,STATEV
,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,
2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PR
EDEF,DPRED,CMNAME,
3 NDI,NSHR,NTENS,NST
A
TV
,PROPS,NPROPS,COORDS,DROT
,PNEWDT,
4 CELENT,DFGRD0,DFGRD1,NOEL,NP
T,LAYER,KSPT,KSTEP,KINC)
IN
CLUDE
‘
ABA_
’
10
CHARACTER*80 CMNAME
DIMENSION STRESS(NTENS),STATEV(NSTATV),
1
DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PR
EDEF(1),DPRED(1),
3 PROPS(NPROPS),COORD
S(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)
user coding to define
DDSDDE,STRESS,STATEV
,SSE,SPD,SCD
and,if necessary,RPL,DDSDDT,DRPLDE,DRPL
DT,PNEWDT
RETURN
END
UMAT
中的应力矩阵、应变矩阵以及矩阵
DDSDDE
、
DDSDDT
、
DRPLDE
< br>等，都是直接分量存
储在前，剪切分量存储在后。直接分量有
NDI
个，剪切分量有
NSHR
个
。各分量之间的
顺序根据单元自由度的不同有一些差异，
所以编
写
UMAT
时要考虑到所使用单元的类别。
下面对
UMAT
中用到的一些变量进行说明：
DDSDDE(NTENS
NTENS)
：一个
NTENS
×
NTENS
的矩阵，称作
Jacobian
矩阵
?
?
/
?
?
?
?
是应力的增量，
?
?
< br>是应变的增量，
DDSDDE(i
，
j)
表示增量步结束时第
j
个应变
分
量的改变引起的第
i
个应力分量的变
化。通常
Jacobian
矩阵是一个对称矩阵，除非在
“
*USER MATERIAL
”语句中加
入了“
UNSYMM
”参数。
STRESS(NTENS)
：应力张量数组，对应
< br>NDI
个直接分量和
NSHR
个
剪切分量。在增量步
的开始，应力张量矩阵中的数值通过
UMA
T
和主程序之间的接口传递到
UMAT
中，在增
量步的结束
UMAT
将对应力
张量矩阵更新。对于包含刚体转动的有限应变问题，一个增
量步调用
UMAT
之前就已经对应力张量进行了刚体转动，因此
UM
AT
中只需处理应力张
量的共旋部分。
UMAT
中应力张量的度量为柯西应力。
STATEV(NSTATEV)
：用于存储与解有关的状态变量的数组，在增量步
开始时将数值传递
到
UMAT
中，也可
在子程序
USDFLD
或
UEXPAN
中先更新数据，然后增量步开始时将更新
后的资料传递到
UMAT
中。在增量步的结束必须更新状态变量矩阵中的数据。和应力
张
量矩阵不同的是：对于有限应变问题，除了材料本构行为引起的资料更新以外，与解<
/p>
有关的状态变量矩阵中的任何向量或者张量都必须通过旋转来考虑材料的刚体运动。
11
状态变量
矩阵的维数通过
ABAQUS
输入文件中的关键词“
*DEPVAR
”定义，关键词下面
数据行的数值
即为状态变量矩阵的维数。
PROPS(NPROPS)
：材料常数数组。材料常数的个数，等于关键词“
*USER
MATERIAL
”中
“
< br>CONSTANTS
”常数设定的值。矩阵中元素的数值对应于关键词“
USER
MATERIAL
”下
面的数据行。
SSE
，
SPD
，
SCD
：分别定义每一增量步的弹性应变能，塑性耗散和蠕
变耗散。
它们对计算结果没有影响，仅仅作为能量输出。
STRAN(NTENS)
：应变数组。
DSTRAN(NTENS)
：应变增量数组。
DTIME
：增量步的时间增量。
NDI
：直接应力分量的个数。
NSHR
：剪切应力分量的个数。
<
/p>
NTENS
：总应力分量的个数，
NTE
NS
＝
NDI
＋
NSHR
。
由于
< br>UMAT
子程序在单元的积分点上调用，增量步开始时，主程序路径将通过
UMAT
的接口进入
UMAT
，单元当前积分点必要变量的初始值将随之传递给
U
MAT
的相应变量。在
UMAT
结束
时，变量的更新值将通过接口返回主程序。
3.4.
U
MAT
的使用方法
我们知
道，有限元计算（增量方法）的基本问题
[7]
是：已知第
n
步的结果（应力，
应变等）
?
n
，
?
n
，然后给出一个应变增量
d
?
n
?
1
，计
算
?
n
?
1<
/p>
，
UMAT
要完成这一计算，并
要计算
DDSDDE(I,J)=
?
?
/
?
?
。
?
?
是应力增量矩阵，
?
?
是应变增量矩阵，
DDSDDE(I,J)
定义了第
J
个应变分量的微小变化对第
I
个应力分量带来的变化。
该矩阵只影
响收
敛速度，不影响计算结果的准确性
（当然，不收敛自然得不到结果）
。
有限元计算的中心问题就是求得节点的位移
（进而应变、应力）
，以使内力和
外力达到平衡：
int
ext
F
(
d
)
?
F
（
3-1
）
d
是节点位移矩阵，黑体字表示矩阵或矢量。除了小变形、线
弹性问题，方程
2-1
是
非性的，要用
迭代的方法解出：
12
ext
int
i
K
?
d
?
F
?
F
(
d
T
n
?
1
n
?
1
)
（
3-2
）
i
?
1
i
d
?
d
n
?
1
n
?
1
?
?
< br>d
(3-3)
i
表示一个增量步内的第
i
次迭代，
n
表示第
n
个增量步。
K
T
是切线刚度，由材
料的
Jacobian
< br>矩阵结合单元计算组装而得。刚度矩阵其实就是力对位移的梯度。要想
快速收敛，
位移增量应沿该梯度方向变化，也就是说，如果
Jacobian
矩阵不是那么准
确，自然
K
T
也不怎么准确，那么满足
3-1
式的位移被找到的速度也就变慢，甚至发
散，
根本找不到。但收敛速度无论慢快，
3-1
式才是判断结果准确
与否的唯一标准。所
以
Jacobian
矩阵不影响结果的准确性，只影响收敛速度的快慢。
以大变形、非性材料为例，整个计算步骤是这样的：
整个外力不是一次加上，而是一点点加上的，不然会发散得不
到结果的。所以，
每一个增量步开始时就是在原来的外力上加上一点点，
得到
ext
F
int
。
根据
3-2
得到位移
增
量
?
d
，此
时要知道力对位移的梯度
K
T
，以尽快
找到满足平衡条件的位移，由材料的
Jacobian
矩阵和
单元结合起来组装得到
(
此处使用
UM
AT
提供的
Jacobian
矩阵
)
。
然后
可计
算应变增量
i
d
?
n
?
1
，调用
UMAT
，得到新的应力，进而得到新的内力，所以，程序
不在乎新的应力是由增量方法得到，还是全量方法得到，而只在乎新应力是否准确。
然后
回到
3-2
，如此循环，直至
3-2<
/p>
右端为
0
，也即满足
3-1
。这样第
n+1
步就完成了，
然后开始第
n+2
步，
即
外力
加上一点点，
按同样的方法求解新的位移。
直至整个外力
全部施加并得到满足
3-1
的位移。
13
4.
材料非线性问题
弹性力学作为精确理论，
从本质上都是非线性的，
< br>早期
Cauchy, Green,Kirchhoff
和
Kelvin
在这些方面都作出了重要贡献。
后来又提出超弹性
(
即具有弹性势的
< br>)
有限变形理
论，由于理论方程的冗长而复杂，且工程应
用也没有提出这方面要求而被搁置。
20
世
纪
40
年代以后，由于橡胶材料、高分子合成材料的迅速发
展和工业领域的大量应用，
非线性弹性与超弹性的研究再次引起科学和工程界的重视，除
了一般理论研究有了新
的发展以外，工程应用计算方法也得到长足的发展。
4.1.
材料的弹塑性本构关系
弹塑性材料进
入塑性的特征是当荷载卸去后存在不可恢复的永久变形。所以，在
卸载情况下，应力应变
之间不再是唯一的对应关系。这是区别于非线性弹性材料的基
本属性。只以加载时应力应
变关系成非线性，还不足以判断材料是非线性弹性还是弹
塑性。但是一经卸载就可以看出
两者的区别。非线性弹性材料沿原路径返回，而弹塑
性材料将依据不同的加载历史卸载后
产生不同的永久变形。
对大多数材料来说，在单调加载的情况
下，存在一个明显的极限应力
?
s
，当
应力低于
?
s
时，材料保持线弹性。而
当应力达到
?
s
以后，则材料开始进入
弹塑性状态。如继
续加载，然后在卸载，材料始终保持永久的塑性变形。如果应力达到<
/p>
?
s
后，应力不再
增加，而材料变形可以继续增加，及变形处于不定的流动状态，则称材料为理想弹塑
性
的。
反之如果应力达到
?
s
后，
再增加变形，
应力也必须增加，
则材料是应变硬化的，
这时应力
?
< br>s
是塑性应变
?
p
的函数，可解析为：
?
s
?
?
s
(
?
p
)
(4-1)
本构关系反应着应力应变
之间的关系。对于弹性材料变形是可以恢复的；而塑性
材料变形是不可以恢复的。典型的
弹塑性应变在卸载后要保持一个永久的变形。如图
3-2
14
图
4-1
塑性应变有下列特性：
（
1
）总应变分为弹性和塑性两部分，即
e
?
ij
?<
/p>
?
ij
?
?
ij
p
(4-2)
或者：
e
<
/p>
d
?
ij
?
d
?
ij
?
d
?
ij
p
(4-3)
（<
/p>
2
）塑性变形取决于加载路径，而应力应变之间没有一一对应的关
系。所以必须确定
二则之间的本构关系，这种本构关系可以用偏微分方程或者增量形式来
描述。
总之，弹塑性理论主要包括以下几个方面：
< br>（
1
）应变张量的分解；
（
2
）应力空间的屈服条件；
（
3
）流动法则；
（
4
）强化法则；
（
5
）协调
性条件。
1
：本构模型
塑性力学的应力
-
应变曲线通常有
5
种简化模型
[8]
：
< br>
（
1
）理想弹塑性模型，用于
低碳钢或强化性质不明显的材料。
（
2
）线性强化弹塑性模型，用于有显著强化性质的材料。
（
3
）理想刚塑性模型，用于弹性应变比塑
性应变小得多且强化性质不明显的材料。
（
< br>4
）线性强化刚塑性模型，用于弹性应变比塑性应变小得多且强化性质明显的材料
。
15
（
5
）幂强化模型，为简化计算中的解
析式，可将应力
-
应变关系的解析式写为
σ＝σ
y
（ε／ε
y
）
n
，式中σ
y
为屈服应力，ε
y
为与σ
y
相对应的应变，
n
< br>为材料常
数。
图
4-2
2
：屈服条件
在复杂应力状态下，
判断物体屈服状态的准则称为屈服条件
[9]
。
屈服条件是各应力
分量组合应满足的条件。
对于金属材料，
最常用的屈服条件
为最大剪应力屈服条件
（又
称
Tres
ca
屈服条件）和弹性形变比能屈服条件（又称
V
on Mises
条件）
。对于岩土材料
则常用
Tresca
屈服条件、
< br>Drucker-Prager
屈服条件和
Mohr-C
oulomb
屈服条件。对于强
化或软化材料，屈服条件将随塑
性变形的增长而变化，改变后的屈服条件称为后继屈
服条件。当已知主应力的大小次序时
，使用
Tresca
屈服条件较为方便；若不知道主应
力的大小次序，
则使用
V
on Mises
屈服条件较为方便。
对于韧性较好的材料，
V
on
Mises
屈服条件与试验数据符合较好。
V
on Mises
屈服准则具体形式
是，
对于各项同性材料，
应力偏量第二不变量
< br>J
2
等于
某一定值时，材料开始
进入了塑性状态。
F
?
3
J
2
?
< br>?
y
(
k
)
?
0
(4-4)
3
：强化法则
对理想的弹塑性材料而言，因无强化作用，所以，整个塑性变形过程中，屈服函
数值保
持一个常量，强化定义了屈服面在应力空间的演化准则。
f<
/p>
(
?
ij
,
k
)
?
0
(4-5)
16
其
中，
k
是强化参数。
通常采用的强化法则有以下几种：
(1)
各向同性强化
此法则规定材料进入塑性变形以后，加载曲面在各方向均匀的向外扩张，没有畸
变。而其形状、中心及其在应力空间的方位均保持不变
[10]
。需要指出的是：各向同性
强化法则主要适用于单调加载情况。如果用于卸载情况，它
只适合反向屈服应力等于
应力反转点的材料，而通常材料不具备这种性质，因此在塑性力
学中还发展了其它强
化准则。
(2)
随动强化
此法则规定材料进入塑性状
态以后，
加载曲面在应力空间作刚体移动而没有转动，
因此初始
屈服面的形状、大小和方向仍然保持不变。
(3)
混合强化
把各向同性强化模型和随动
强化模型加以组合，得到混合强化模型。它假定在塑
性变形过程中，加载曲面不但作刚性
平移，还同时在各个方向作均匀扩大。
在以上几种强化模型中
，各向同性强化模型应用最为广泛。本文也是采用该硬化
法则，这一方面是由于它便于进
行数学处理；另一方面，如果在加载过程中应力方向
（或各个应力分量的比值）变化不大
，采用各向同性强化模型的计算结果与实际情况
也比要符合。随动强化模型可以考虑材料
的包兴格（
Bauschinger
）效应，在循环加载
或可能出现反向屈服的问题中，需要采用这种模型。
由于塑性变形与变形历史有关，
因此
反映塑性应力
-
应变关系的本构关系用应变增
< br>量形式给出比较方便。用应变增量形式表示塑性本构关系的理论称为塑性增量理论。
增量理论的本构关系在理论上是合理的，但应用比较麻烦，因为要积分整个变形路径
才
能得到最后结果。因此，又发展出塑性全量理论，即采用全量应力和全量应变表示
塑性本
构关系的理论。在比例变形的条件下，可通过积分增量理论的本构关系获得全
量理论的本
构关系。当偏离比例变形条件不多时，全量理论的计算结果和实险结果比
较接近。本文的
程序都是基于增量理论。
4.2.
非线性有限元算法理论
对于非
线性问题，在有限元求解该问题时，对一个自由度总可以表达成
H
?
?
f
?
0
，式中，
?
为基本未矢量。如果是线性
问题，
H
与
?
无关，而
?
是一次项，
17
显然这是一个线性方程。如果
H
与
?
相关，
则方程的
?
出现非一次项，变成非线性问
题，在实际工程中，特别是塑性成型问题，材料的几何方程，本构方程以及边界条件
往
往是非线性的也体现在
K
中出现了
?<
/p>
，所以变为了非线性问题，要得到最基本的未
知量，就必须求解非
线性方程组
1
：直接迭代法
又称常刚度法
[11]
，这是种最简
.
单的求解方法，在每次求解前，利用上次的
?
解来求出
这一次的
H
值，然
后利用
f
和
H
的倒数的乘积求出
?
的当前值
?
1
?
?
?
H
(
?
)
f
(4-6)
表达为迭代形式
(
< br>n
?
1)
?
?
H
(
?
n
)
?
1
f<
/p>
(4-7)
?
上式可
以看出，这种方法首先需要有一个初始的
?
值，以便开始迭代。
另外，每一次
求解都需要对
H
求倒数，
如果求解方程组，就是对刚度矩阵求逆，这种方法在求解中
控制两次求解
?
之差，当其值很小时，就认为接近真实值了，迭代结束
图
4-3
2
：
Newton-
Raphson
方法
18
Newton-Raphso
n
方法的算法与常刚度法不同
[12]
，如果
?
得近似表达式
H
?
?
f
?
< br>0
是不
成立的，存在着残余值，即
H
?
?
f
?
0
，此式也可以作为近似值与真实值的差值量度，
实际上在具体计算时，也可以控制其值，当
?
?
H
?
?
f
?
0
极小时，就认为
?
接近真实
值了，当第
n
?
1
次迭代的
?
值
?
n
< br>?
1
?
?
n
?
?
?
n
是真实解，则可以按照
Taylor
级
数展开得到
?
(
?
n
?
1
)
?
?
(
?<
/p>
n
)
?
(
d
?
n
n
)
?
?
?
0
d
?
d
?
n
H
(
?
n
)
?<
/p>
n
n
n
?
?
?
?
?
(
?
)
/(
)
?
?
?
< br>(
?
)
/(
)
?
?
?
(
?
n
)
/(
H
(
?
n
)
?
n
?
H
(
?
n
))
d
?
d
?
n
n
图
4-4
3
：切线刚度法
在复杂非线性问题求解中，刚度
H
与
?
的大小是有一定关系的，在用增量法来求解这
种问题时，<
/p>
H
就等于结构任一点处力与位移的曲线的局部梯度，称为切线刚度
[13]
，刚
度矩阵的倒数很难用自变
量显示表达，通过增量方式求解，在每一步荷载增量范围内
把问题线性化，求解方法与<
/p>
Newton-Raphson
方法相同
19
总结以上可以得到：
以上几
种算法中，通过比较，不难发现，直接迭代法采用了固定的刚度，适合解
决非线性程度不
高的本构关系，而切线刚度法采用了变化的刚度，在每一步上都做了
实时的修正，对非线
性程度较高本构关系任然有效，在效率和迭代精度方面，切线刚
度法采用的修正更符合非
线性材料的应力应变关系，具有较大的优势，这也是本文采
用切线刚度法计算的原因，当
然，非线性有限元算法还有很多，切线刚度法也不见得
就是最好的能解决所有问题的算法
，但是它是在程序开发难度不高和精度方面较高的
条件下相对来说最好的
本文采用的本构关系是同性硬化弹塑性模型
< br>[14]
，采用
Mises
屈服
准则，下面将根
据
J2
理论
[15]
，分别推导常刚度法和切线刚度法计算该问题的的算法公式
4.3.
增量理论常刚度法公式推导
由应力应变关系得：
d
{
?
}
< br>?
d
{
?
}
e
?
d
{
?
}
p
(4-8)
(4-9)
d
{
?
}
?
[
D
]
e
(
d
{
?
}
?
d
{
?
}
p
)
< br>其中
[
D
]
e
是弹性矩阵，它的表达式为
?
1
?
?
?
?
1
?
?
?
?
?
?
1
?
?
E
< br>(1
?
?
)
?
D
e
?
(1
?
?
)(1
?
2
?
)
?
0
?
?
?
0
?
?
?
0
?
?
d
{
?
}
p
< br>1
对
1
0
0
0
1
?
2
?
2(1
?
?
)
0
0
1
?
2
?
2(1
?
?
)
0
称
?
1
?
?
0
0
0
< br>?
?
?
?
?
?
?
?
?
?
?
?
?
?
1
?
2
?
?
2(1
?
?
)
?
?
是塑性应变增量，它的表达式为
d
{
?
}
p
?
d
?
p
?
?
?
{<
/p>
?
}
20
其中
d<
/p>
?
p
为等效塑性应变增量，它的表达式为
{
?
?
T
}
[
D
]
e
?
{
?
}
d
?
< br>p
?
d
{
?
}
?
?
T
?
?
'
H
?
{
}
[
D
]
e
?
{
?
}
?
< br>{
?
}
[
D
]
e
?
?
3
3
G
'
'
'
'
'
?
[
D
]
e
[
?
< br>x
,
?
y
,
?
z
'
,
2
?
xy
,<
/p>
2
?
yz
,
2
?
zx
]
T
?
[
?
x
,
?
y
,
?
z
,
?
xy
,
?
yz
,
?
zx
]
T
?
?
?
?
2
?
?
([
D
]
?
?
T
T
?
?
?
?
?
?
e
?
?
?
?
)
?
?
?
?
?
?<
/p>
?
?
?
?
[
D
]
?
e
H
'
为切线模量
对于
3
维空间问题，流动方向：
?
< br>?
?
{
?
}
?
3
2
?
[
?
'
'
'
'
'
T
x
,
?
y
,
?
'
z
< br>,
2
?
xy
,
2
?
yz
,
2
?
zx
]
等效应力
?
?
1
2
{
[(
?
x
?
?
y
)
2
?
(
?
2
y
?
?
z
)
?
(
?
z
< br>?
?
x
)
2
]}
应力偏量
?
'
x
?
?
x<
/p>
?
(
?
x
?
?
y
?
?
z
)
/
3
?
'
y
?
?
y
?
(
?
x
?
?<
/p>
y
?
?
z
)
/
3
?
'
z
?
?
z
?
(
?
x
?
?
y
?
?
z
)
/<
/p>
3
?
'
xy
?
?
xy
?
'
yz
?
?
yz
?
'
zx
?
?
zx
4.4.
增量理论切线刚度法公式推导
本文采用的是一种切线刚度法，其应力应变关系为
d
{
?
}
?
d
{
?
}
e
?
d
{
?
}
p
d
{
?
}
?
[
D
]
e
(
d
{
< br>?
}
?
d
{
?
}
p
)
21
(4-10)
(4-11)
{
?
?
}
T
3-14
两端
左右同乘
?
{
?
}
，得到：
{
?
?
?
{
?
}
}<
/p>
T
d
{
?
}
?
{
?
?
T
?
{
?
}
}
[
D
]
e
(
d
{
?
}
?<
/p>
d
{
?
}
p
)
对于强化材料
，
Mises
准则
< br>d
?
?
{
?
?
?
{
?
}
}
T
d
{
?
}
?
H
'
d
?
p
代入可得
H
'
d
?
?
?
p
?
{
?
{
?
}
}
T
[
D<
/p>
]
?
?
T
?
?
e
d
?
?
{
?
{
?
}
}
[
D
]
e
?
{
?
}
d<
/p>
?
p
其中，<
/p>
d
?
p
为等效塑
性应变增量，它的表达式为
{
?
?
}
T
d
?
?
{
?
}
[
D
]
< br>e
p
?
d
{
H
'
?
{
?
?
?
{
?
}
}
T
[
D
]
?
?
?
}
e
< br>?
{
?
}
由于：
[
D
]
?
?
?
[
D
]
3
'
'
3
G
[
?
'
'
'
T
e
?
< br>?
?
?
e
2
?
[
?
x
,
?
y
,
?
'
z
,
2
?
T
xy
,
2
?
yz
,
2
?
zx
]
?
?
x
,
?
y
,
?
z
,
?
xy
,
?
yz
,<
/p>
?
zx
]
T
([
D
]
?
?
T
?
?
?
?
?
?
e
?
?
?
?
)
?
?
?
?
?
?
?<
/p>
?
?
?
[
D
]
?
e
由流动法则可知
d<
/p>
{
?
}
p
?
d
?
?
?
p
?
{
?
}
应用上式得
?
d
{
?
}
?
?
[
D
]
?
?
e
?
?
{
?
}
{
?
?
< br>?
{
?
}
}
T
[
D
]
?
e
?
?
[
D
]
e
?
?
d
{
?
?
H
'
< br>?
{
?
?
?
{
?
}
}
T
[
D
]
?
?
e
?
{
?
}
?
?
}
?
?
< br>
所以得到
d
{
?
}
?
[
D
]
ep
d
{
?
}
这就是切线刚度法的矩阵表达式
22
(4-12)
其中
[
D<
/p>
]
ep
为弹塑性矩阵，在
ABAQUS
里面称为雅可比矩阵，它的表达式为
[<
/p>
D
]
ep
?
[
D
]
e
?
[
D
]
p
其中
[
D
]
e
是弹性矩阵，它的表达
式为
?
?
1
?
?
?
1
对
?
1
?
?
?
?
?
1
?
?
1
< br>称
D
E
(1
?
?
)
?
e
?
?
1
?<
/p>
?
(1
?
?
)(1
?
2
?
)
?
?
0
0
0
1
?
2
?
?
2(1
?
?
)
?
< br>?
0
0
0
0
1
?
2
?
?
2(1
?
?
)
?
?
?
0
0
0
0
0
D
p
为塑性矩阵，它
的表达式为
[
D
]
?
?
e
{
?
?
}<
/p>
T
[
D
]
?
{
?
}
?
{
?
}
[
D
]
e
p
?
H
'
?
{
?
?
T<
/p>
?
?
?
{
?
}
}
[
D
]
e
?
{
?
}
对于
3
维空
间问题
流动方向
?
?
?
{
?
}<
/p>
?
3
2
?
[
?
'
'
,
?
'
'
'
'
T
x
,
?
y
z
,
2
?
xy
,
2
?
yz
,<
/p>
2
?
zx
]
所以
?
?
'2
?
x
?
?
?
'
'
x
?
< br>y
?
'2
y
对
?
?
D
]
9
G
2
[<
/p>
?
?
'
'
x
?
y
?
'
?
'
'2
y
z
?
z
< br>称
?
p
?
(
H
'
?
3
G
)
?
2
?
?
?
'
'
x
?
?
'
'
?
'
< br>'2
?
?
xy
< br>y
?
xy
?
'
z
xy
?
xy
?
?
'
?
'
?
x
yz
?
'
'
'
'
'
?
'2
y
?
yz
?
'
z
?
yz
?
xy
?
yz
yz
?
?
?
?
'
x
?
< br>'
zx
?
'
?
'
'
'
'
?
'
'
y<
/p>
zx
?
'
z
?
zx
?
xy
?
zx
yz
?
zx
?
'2
?
zx
?
?
最后推导可得，弹塑性矩阵的表达式
23
?
?
?
?
?
?
?
?
?
?
?
?
?
1
?
2
?
?
< br>2(1
?
?
)
< br>?
?
?
?
1
?
?
'2
?
??
x
?
1
?
2
?
?
?
?
'
?
??
?
y
x
?
1
?
2
?
?
?
1
?
?<
/p>
'
?
??
?
z
x
E
?
1
?
2
?
[
D
]
ep
?
?
1
?
?
?
'
'
?
??
?
xy
z
?
?
?
?<
/p>
??
'
?
'
z
yz
?
?
'
?
?
??
z
'
?
zx
?
其中
1
?
?
'2
?
??
y
1
< br>?
2
?
对
1
?
?
?
?
?
z
'2
1
?
2
?
'
?
??
z
'
?
xy
'
?
??
z
'
?
yz
'
?
??
z
'
?
zx
?
1
?
2
?
< br>'
?
??
y
?
z
'
称
1
'2
?
??
xy
2
'
'
?
??
xy
?
y
z
'
'
?
??
xy
?
zx
'
'
?
??
y<
/p>
?
xy
'
'
?
??
y
?
yz
'
'
?
??
y
?
zx
1
'2
?
??
yz
2
'
'
?
??
xy
?
zx
?
?
?
?
?
?
?
< br>?
?
?
?
?
?
?
?
1
'2
?
?
??
zx
2
?
?<
/p>
?
9
G
E
G
?
2
?
2
(
H
'
?
3
G
)
2(1
?
?
)
H
'
为切线模量，对本构关系求导得到
H
'
?
总结推导过程：上面的推导过程看
似复杂，其实核心的问题只有一个，即两者对于塑
性阶段应力更新的算法不同。
常刚度法采用的是：
d
{
?
}
?
[
D
]
e
(
d
{
?
}
?
d
{
?
}
p
)
< br>d
?
d
?
p
而切线刚度法采用的是：
d
{
?
}
?
[
D
]
ep
d
{
?
}
把握好了两者的本质上的区别，对于两者的算法设
计和程序开发问题便迎刃而解
24
5.
UMAT
程序设计和编码
本章将严格按照前一章推导的公式展开程序设计和编码，为了
便于编程，本文将
本构关系做了抽象化处理，即将其描述成一个含参数的表达式，改变参
数即可应用于
不同的模型，
这样做的好处是能保证程序的复用性
，
这也是本文反复强调的使用
UMAT
的原则。
5.1.
本构关系描述
本文采用各向同性硬化弹塑性材料，材料参数如下：
?
?
?
E
?
e
(
?
?
?
< br>T
)
?
B
?
?
A
(
?
)
?
C
(
?
?
?
T
)
?
P
(5-1)
弹性部分：
?
?
E
?
e
，
弹性模量
E=200000Mpa
，泊松比
Mu=0.3
图
5-1
弹性部分本构关系
B
?
?
A
(
?
)
?
C
，为了研究方便，取
A=700
，<
/p>
B=0.5
，
C=400
P
塑性部分：
25
图
5-2
塑性部分本构关系
将
2
个曲线统一到同一个坐标系（为方便显示，
x
轴标注时扩大了
1000
< br>倍）
图
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
：
读取
AB
AQUS
定义的材料常数和状态变量
(
这里只定义了一个状态变量
)
，
材料常
数
为，弹性模量
E
，泊松比
Mu
，屈服应力
Yield0
，参数
A,B,C,
并且计算出剪切模量
< br>G
,
状态变量为等效塑性应变
E
QPLAS
3
：读取应力分量，计算平均应力，应力偏量以及
Mises
等效应力
平均应力：
?
cp
?
(
?
x
?
?
y
?
?
< br>z
)
/
3
'
?
?
x
?
?
x
?
?
cp
?
'
?
?
y
?
?
y
?
?
cp
?
'
< br>?
?
z
?
?
z
?
?
c
p
?
'
?
?<
/p>
xy
?
?
xy<
/p>
?
'
?
?
yz
?
?
yz
?
?
'
?
?
zx
应力偏量：
?
zx
?
?
Mises
等效应力：
2
'2
'2
'2
< br>'2
'2
(
?
< br>x
?
?
y
?
?
z
'2
?
2(
?
xy
?
?
yz
?
?
zx
)
3
<
/p>
4
：根据
3
计算
的
Mises
等效应力和
2
读取的屈服应力
Yield0
比较，如果
Mises
等效应
力小于屈服应力，表明此时材
料未屈服，那么转到
5
，否则转到
6
5
：雅可比矩阵
,
初始化为
0
，计算弹性矩阵，按照弹性理论更新应力
27
?
1
?
?
?
?
1
?
?
?
?
?
?
1
?
?
E
(1
?
?
)
?
D
e
?
(1
?
?
)(1<
/p>
?
2
?
)
?
0
?
?
?
0
?
?
?
0
?
?
d
?
?
?
?
D
e
d<
/p>
?
?
?
1
对
1
0
0
0
1
?
2
?
2(1
?
?
)
0
0
1
?
2
?
2(1
?
?
)
0
称
?
1
?
?<
/p>
0
0
0
?
?
?
?
?
?
?
?
?
?
?
?
?
?
1
?
2
?
?
2(1
?
?
)
?
?
6
：雅可比矩阵
,
初始化为
0
1
）：计算切线模量
H'
H
'
?
A
*
B
*(
?
< br>P
?
Yield0/E)
B
?
1
,
注意到当等
效塑性应
?
P
?
0
时对应于本构关系的屈服点，
此
时
的
H
不能通过上式计算，可以取此时的
H
为弹性模量
2
）：计算等效塑性应变增量并更新
'2
'2
'2
'2
'2
6*
G
*
?
*
?
?
x
,
?
y
,
?
z
'2
,
?
xy
,
?
yz
,
?
zx
?
d
?
p
?
2*
?
*
H
?
9*
< br>G
*(
?
?
?
?
?
?
2
?
?
2
?<
/p>
?
2
?
)
2
'
'2
x
'2
y
'2
z
'2
xy
'2
yz
'2
zx
*
d
?
?
?
EQPLAS=DEQPLAS+EQPLAS
更新状态变量
3
）计算流动方向
< br>?
?
3
'
'
?
[
?
x
,
?
y
,
?
z
'
,
2
?
xy
,
2
?
yz
,
2
?
zx
]
T
?
?
?
?
2
?
4
）计算塑性应变增量
d
{
?
}
< br>p
?
d
?
p
?
?
?
{
?
}
5
）更新应力
d
{
?
}
?
[
D
]
e
(
d
{
?
}
?
d
{
?
}
p
)
< br>
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,CEL
ENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
29
DIMENSION
STRESS(NTENS),STATEV(NSTATV),DDSDDE(NTENS,NTENS),
1
DDSDDT(NTE
NS),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)
DOUBLE PRECISION
E,Mu,Yield0,A,B,C
DOUBLE
PRECISION
H,W,CEGMA_EQ,CEGMA_X,CEGMA_Y,CEGMA_Z,TAO_XY,
1
TAO_YZ,TAO_ZX,CEGMA_CP,CEGMA_SX,
CEGMA_SY,CEGMA_SZ,TAO_SXY,TAO_SYZ
C
定义材料常数
C
定义中间变量
常刚度法
INCLUDE
'ABA_'
CHARACTER
*8 CMNAME
2
,TAO_SZX
C
定义状态变量
C
DOUBLE
PRECISION
EQPLAS
定义更新变量
DOUBLE PRECISION
DEQPLAS
CCCCCCCCCCCCCCCCCCCCCCCCCC2
：读取变量
C
C
C
弹性模量
泊松比
Mu=PROPS(2)
G=E/TWO/(ONE+Mu)
屈服应力
Yield0=PROPS(3)
A=PROPS(4)
B=PROPS(5)
C=PROPS(6)
EQPLAS=STATEV(1)
Yield=A*(EQPLAS+Yield0/E)**B+C
E=PROPS(1)
C
剪切模量
C
参数
C
读取状态参数
CCCCCCCCCCCCCCCCCCCCCCCCCC3
：读取应力分量，
计算平均应力，应力偏量以及
Mises
等效应力
C 6
个应力分量
CEGMA_X= STRESS(1)
CEGMA_Y= STRESS(2)
TAO_XY=
STRESS(4)
TAO_ZX= STRESS(5)
CEGMA_Z= STRESS(3)
30