ABAQUS在反复加载卸载应力应变计算说明

ABAQUS 反复加载卸载应力应变计算说明
1、 塑性模型
采用Armstrong and Frederick model (AF 模型)[1] （1） 屈服准则（Mises 屈服）
屈服函数为2)~
~(:)~~(2
3Y s s F ---=
αα 式中：s ~
为应力偏量，α~为总背应力，Y 为屈服极限 （2） 流动准则（the associated flow rule ）
σ
λε~~∂∂=F
p 式中：p
ε
~为塑性应变对时间的微分，λ 为待定量[2]，σ~为应力张量 （3） 硬化准则
<1> p h i i p i i )()()()(~~3
2~αζεα
-= 式中：)
(~i α
为背应力分量对时间的微分；)
(i h ，)(i ζ为材料常数为已知量，p 为等效塑性应变对时间的微分。

<2> )1()
(1
)(0p
M
i i i e r Y Y ζ
-=-+
=∑
式中：Y 为对应于等效塑性应变p 的屈服极限，0Y 为初始屈服极限为已知量，)
(i r 为材料常
数为已知量。

2、 反复加载卸载应力应变计算过程说明
假设在受载前，物体的初始应力、应变以及背应力均为零
2.1加载过程计算
<1> 外力不足以使得物体中的任何一点的Mises 应力值大于屈服极限
此时：塑性应变、背应力均保持为零，屈服极限保持不变。

应力由e e D εσ~:~~~=计算，
总应变值等于弹性应变。

<2> 外力使得物体中的任何一点的Mises 应力值大于屈服极限
为了说明ABAQUS 是如何确定应变增量1~+∆n ε，有必要对ABAQUS 求解材料非线性问
题进行简单介绍[3]。

ABAQUS 首先将载荷分为若干个微小增量，如图1所示。

当结构收到一个微小增量P ∆时，ABAQUS 用与初始结构位移相对应的初始刚度矩阵K 和载荷增量P ∆计算出结构在这
一增量后的位移修正a c 、修正后的位移值a u 和相应的新的刚度矩阵a K 。

ABAQUS 用新的刚度矩阵计算结构的内力a I ，载荷P 和a I 的差值为迭代的残余力a R 。

如果a R 在模型内的每一个自由度上的值都为零或小于一个给定的容差，如图1所示的a 点，则结构处于平衡状态即ABAQUS 计算到的内外力是平衡的。

若假设整个物体只有一个单元，则位移修正a c 就为应变增量1
~+∆n ε，刚度矩阵就是UMAT 程序中的雅可比矩阵。

根据塑性变形时的雅可比矩阵计算公式：
d e p I L G D D ~~:~~4~~~~12--=其中e D ~~为弹性状态下的雅可比矩阵，可知p D ~~较e D ~~小。

UMAT 计
算出的应力便是结构的内力a I 。

在UMA T 程序中通过计算：
2111)~~(:)~~(2
3Y s s F n tri n n tri n tri n ---=
+++αα
式中：11~:~~~~++∆+=n e n tri n D εσσ 并判断若tri n F 1+是否大于零来决定节点是否达到屈服。

若tri n F 1+小于零则节点未达到屈服，若tri n F 1+大于等于零节点达到屈服，需要进行塑性迭代与重新计算雅可比矩阵。

图1 增量法迭代原理
ABAQUS 处理塑性加载过程问题可以描述为：首先读入上一状态的应力、应变、背应力及屈服强度，使用下标n 表示（对应时间为n t ）。

再根据上一步的刚度矩阵n K 和载荷增
量P ∆计算出给定的应变增量1~+∆n ε，然后在UMA T 程序中求出1+n t （1
1++∆+=n n n t t t ）时刻满足所用塑性模型给出的屈服准则、流动准则、硬化准则下的应力、应变、背应力及屈服强度。

UMAT 程序利用已知的n σ~，)(~i n α，n ε~，p n
ε~，n p ，n Y 和给定的1~+∆n ε首先求出可满足所用塑性模型给出的屈服准则、流动准则、硬化准则的∆p ，再使用公式依次求出其余量。

等效塑性应变：11++∆+=n n n p p p 屈服极限：)1(1
)
(1
)(01+-=+-+
=∑n i p M
i i n e r Y Y ζ
塑性流动：)
~~(:)~~()
~~(~1
)()(111
)()(111
)
()(111
∑∑∑=++=++=+++---=M
i i n i n tri
n M
i i n i n tri n M
i i n i n tri n n s s s n αθαθαθ
式中：1
)()
(111
++∆-=
n i i n p ζθ
塑性应变增量：1
11
~2
3~+++∆=∆n n p
n n p ε 背应力：)~3
2~(~1)()
()
(1)
(1p
n i i n
i n i n h +++∆+=εαθα 应力：tri
m
n n i n n I n Y σασ~~3
2~~11)
(1
1++=++++ 式中：)(31332211tri tri tri tri
m σσσσ++=
再计算出：n ε
~，p n
ε~ 加载过程按照此方式进行计算。

2.2当外力卸载时
同样已知：n σ~，)(~i n α，n ε~，p n
ε~，n p ，n Y 并且给定：1~+∆n ε。

这里考虑卸载过程的第一步迭代。

由于外力开始卸载，计算得到的)~~(:)~~()~~(:)~~(11n n n n n tri n n tri n s s s s αααα--<--++。

由于上一步迭代材料处于屈服状态，所以在UMAT 程序运行结束时有下面关系成立：
0)~~(:)~~(232=---=n n n n n Y s s F αα所以此时的0)~~(:)~~(2
32111<---=+++n
n tri n n tri n n Y s s F αα根据程序判断不进行塑性计算，只计算弹性雅可比矩阵即按照弹性方式进行卸载。

各个参数有关系： )()(1~~i n i n αα=+， p n p n εε~~1=+，n n p p =+1，n n Y Y =+1 弹性应变：1
1~~~++∆+=n e
n e
n εεε 总应变：p
n e
n n 1
1~~~+++=εεε 按照此方式，直到卸载完成。

这里不讨论反向加载。

2.3 再次加载时
由于)(~i n α，p n ε~，n Y 在卸载过程中不变化，所以在第二次加载过程中，读入的n σ~，)(~i n α，n ε~，p n ε~，n p ，n Y 中)(~i n α，p n
ε~一定是不为零的，并且0Y Y n >。

对于n σ~是否为零，与残余应变是否引起了残余应力有关。

对于一个单元残余塑性应变是无法引起残余应力的；对于多
个单元构成的物体，当卸载完成后，一些单元的残余塑性应变使得它们与周围的单元之间产
生作用力从而产生残余应力。

此时读入的n
σ~是不为零的。

这是与第一次加载初始应力状态有区别的。

根据UMAT 弹性试应力的计算公式11~:~~~~++∆+=n e n tri n D εσσ，若物体中含有残余
应力则在计算过程中直接代入程序计算。

由给定的应变增量1
~+∆n ε计算出下一时刻的2111)~~(:)~~(2
3Y s s F n tri n n tri n tri
n ---=+++αα。

若
tri n F 1+小于零则节点未达到屈服按照2.1中的弹性加载过程进行计算；若tri n F 1+大于等于零按照2.1中的塑性加载过程进行计算。

分析：若要此模型出现图2所示的应力应变曲线必须进行双向加载。

在单向反复加载时，若所受最大值是不变的则，在后续的加载过程中均为弹性变形。

图2 应力应变曲线
问题：我这样的分析合理吗？
[1] 2002-implementation of cyclic plasticity models based on a general form of kinematic hardening
[2] 有限单元法 王勛成 清华大学
[3] 2005-ABAQUS 用户材料子程序开发及应用（硕士论文-华中科技大学）。