弹塑性有限元方法

第三章 弹塑性有限元方法的实施

§3.1 增量平衡方程和切线刚度矩阵 1、 分段线性化的求解思想

塑性变形的特点决定了塑性本构关系的非线性和多值性，上面由塑性增量理论给

出了塑性应力—应变关系{}{}ep d D d σε=????

其中 [][]

{}{}[]{}[]{}

T

ep T

F

F

D D

D D F

F

A D

σσ

σ

σ

????=-

??+

??????

说明当前应力状态不仅与当前应变有关，而且和达到这一变形状态的路径（加载历史）有关。这里包含了屈服准则、强化条件和加卸载准则。

由此，对物理非线性问题，通常采用分段线性化的纯增量法和逐次迭代的方法求解。即将加载过程分成若干个增量步，选择其中任意一个增量步建立它的增量平衡方程并求解，对整个过程的求解有普遍意义。

2、 增量平衡方程和切线刚度矩阵

设t 时刻（加载至i -1步终），结构（单元）在当前载荷（广义体力{}v f 和表面力{}s f ）

的作用下处于平衡状态，此时物体内一点的应力、应变状态为{}{}σε、。在此基础上，施加一个载荷增量{}{}v s f f ??和，即从t t t →+?时刻，则在体内必然引起一个位移增量{}u ?和相应的{}σ?、{}ε?，只要{}{}v s f f ??和足够小，就有{}{}ep D σε?=?????。

倘若初始状态{}σ已知，加载过程已知，则ep D ????可以确定（即p

ij d ε?可以确定，然后

可在硬化曲线上得到1p ε所对应的硬化系数）于是上面的方程成为线性的。在t t t →+?这一增量过程中，应用于虚功原理可得到如下虚功方程：

()()()0e

e

T T T V V s s V S f f u dV f f u dS σσδεδδ??+?-+??-+??=??

?? （1）

根据小变形几何关系u N q B q ε?=??=?和，再由虚位移()q δ?的任意性，并设

()()e

e

T T v v s s V

S P P N f f dV N f f dS +?=

+?+

+??

?

，展开后，其中单元在t 时刻载荷等效节点

力：e

e

T T v s V

S P N f dV N f dS =

+

?

?

；t ?内增量载荷的等效力e

e

T T v s V

S P N f dV N f dS ?=

?+

??

?

。

这样，由方程（1）可得平衡方程：

[]{}{

}e

T

V B d V P P

σσ+?=+?? （2） 即： ()0e

e

T T t t V V F B dV B dV P P σσ+?=+?-+?=??

因为t 时刻（第i 步终）结构处于平衡状态 0e

T

t V F P B dV σ=-=? （3）

这样（2）式变为: e

T V P B d V σ?=

?? 即：0e

T V F B dV P σ?=?-?=? （4）

将{}{}ep d D d σε??=??和{}B q ε?=?代入上式得增量平衡方程：

e

T

ep V B D BdV q P F ??-?=?? （5）

对增量位移求导：

()

()e

T e ep t V d F B D BdV K d q ?==?? （6）

于是（5）式成为 e

t K q P ?=?? （7）

e

t K 为单元切向刚度矩阵。集合所有单刚后得到结构总的增量平衡方程

T K q P ?=?? （8）

方程（8）是线性的，可以直接求解。 3.2 硬化系数H '的数值表示

根据单一曲线定理，对于一般稳定性硬化材料，在其简单加载过程中，σ和ε之间存在着一一对应的确定的函数关系()εσΦ=，这一关系可用单向拉伸实验来确定。 例如，对于Mises 各向同性硬化材料

p d d H εσ/=' （8）

在有限元分析中，作为初始参数应把这一曲线输入（用函数或数字的形式），在加载过程中弹塑性矩阵不断地修改，根据当前的应力或应变来确定。目前，硬化曲线的输入格式有两种： 1） 解析表达式

根据单一曲线定理，由单向拉伸试验曲线直接得出硬化曲线的解析式。例如： （a ）Mises 各向同性线性硬化材料

单向拉伸曲线有： 当 s σσ≤ εσE =

当 s σσ≥ ()t s s E εεσσ-+= （9）

则有 111111

11t p p e

t

t EE d d d H d d d d E E E E

σσσεεεε'=

====--- （10） -----------------------------------------------------------------------------------------------------------------附：对于一般材料的硬化曲线的求法（求H '）

如单拉曲线 则硬化曲线

根据 11~εσ ===》 p εσ~===》??==p ij p p d d εεε1 其中单拉时等效应变为 ()113

2

32ευεεε+==

ij ij 因为 εε=1，υεεε-==32，平均应变为 ()()132********

ευεεεε-=++=

所以10εεε∴=+ ，当 2

1

=υ时 εε=1

（b ）Mises 各向同性幂硬化材料

单向拉伸曲线有： 当 s σσ≤ εσE =

当 s σσ≥ m A σε= （11）

由屈服点条件：m

S S

E A εε=得 1m

S

A E ε-= 据（8）式得 111

11111

1111111p e m e d d EB H d d d E B d d Am E d d σσεεεσσεεε---'=

====--????-- ? ?????

（12）

其中： 1111)/(--==m S m Em Am B εεε

2） 根据离散的单拉实验数据，采用样条插值计算H '

（参看清华大学孟凡中教材：弹塑性有限变形理论和有限元方法） 3.3 过渡单元弹塑性矩阵的确定

1． 三种变形状态

弹塑性变形体中，在一个载荷增量步内可能有三种变形状态： 1）弹性区：加载前后均处于弹性状态，故采用弹性阵不变。

2）塑性区：加载前后均处于塑性状态，其弹塑性矩阵ep D ????由塑性增量理论确定（与当前应力水平和塑性变形增量的总量有关）

3）过渡区：加载前处于弹性状态，加载后进入塑性状态，所以，在这一过程中采用弹

性矩阵[]D 或最终的ep D ????都不合适，必须寻找一个合适的弹塑性矩阵ep D ????。

2． 加权平均的弹塑性矩阵ep D ????

1）过渡单元在加载后的应力计算（以单拉状态为例）

在t ?时间步内施加一个增量载荷后，讨论某单元的应力应变状态。

设某单元加（卸）载前的应力状态0σ，相应的应变0ε（A 点）处于弹性状态（弹性区间O ’C ）。加载后，按弹性计算得到应变增量ε?，到达B 点。显然B 点不是实际的应力状态，因为已经超过了C 点，进入了塑性变形阶段，假设实际应到达D 点。

该增量步的弹塑性矩阵ep D ????是未知的，

它的大小应该和该增量步内弹塑性应变所占比例有关，只能经过迭代试算得出。

因为: p e 111εεε?+?=? （13） 且设 11εε??=?m p （14） 则 11)1(εε??-=?m e （15） 显然，10≤≤m ，且m =0时是全弹性，m =1时是全塑性。 实际应力增量为

1111111

11

(1)(1)e p AC CD ep ep ep D m D m m D mD D σσσσσεεεε?=?+?=?+?=-?+???=-+?=??? （16）

推广到一般的应力状态{}σ?，{}ε?为

()[]1ep ep m D m D D σεε???????=-+?=??????? （17） ()[]1ep ep D m D m D ????=-+????

（18） ep D ????

－加权平均弹塑性矩阵；m －比例系数。 2）比例系数m 的迭代公式

已知A 点的0σ和0ε，同时由到达A 点的路径确定s

ε'和s σ'

s

p

εεεεε

'-?+=?=?0 由定义：

001p s s m εεεεεεεεε

'

'+?--?===+??? （19）

3． 过渡单元m 的确定

1）确定是过渡单元。即在第（i -1）个增量步终（求解结束时）某单元是弹性的应力、

应变状态0σ和0ε，且0s

σσ'≤（或 0s εε'≤），进入第i 个增量步（t ?内载荷增量 i P ?），按弹性计算到达B 点，其应力0s σσ'≥，应变0s εε'≥，可以确定该单元在第i 个增量步内是过渡单元。 2）关于m 的迭代过程：

按弹性矩阵[]D 计算该单元的切线刚度矩阵t k ，然后和其他单元集合成总刚T K ，代入结构的增量平衡方程并求解得总位移向量{})

1(i q ?，从{}1

i q ?中提取该单元的

{})1(i q ?，并求出{})1(i ε?，{})1(i σ?及{})1(i ε?。代入（19）式，计算出)

1(m ，再将)

1(m 代入（18）得()

1ep D ????， ep D ?

???中的ep D ????与当前应力和应变状态有关。当前应力为： {}{}{}()11i i σσσ?+=-，{}{}{}()11i i εεε-=+?

（a ）按第1次迭代的计算值()

1ep D ????代入该单元计算切线刚度阵，并与其他单元集

合组装求解总的增量方程得{}()

2i

q ?及相应的{}()

2i

ε?，{})

2(i

σ?及

{}{}{}{}()212i i i

σσσ?+=-。此时0ε和s ε没有改变，再代入（19）式计算)2(m 。将(2)m 和{}{}

2σ（若是硬化材料，还要根据当时塑性应变总量确定H '的值）代入

计算()

2ep D ????。

(b) 依次类推，求出 ()

3ep D ????；()

4ep D ????…………直至前后两次的m 值十分接近（到达给定的允许误差范围）停止迭代。

(c ）将迭代终止时的ep D ????作为该单元的弹塑性矩阵，求单刚集合，解方程，求出

{}n i q ?及相应的{}n i ε?，{}n i σ?将其累加到上一步的终值上作为下一步的初值。

总位移 {}{}{}n

i i q q q ?+=-1

{}{}{}n

i i εεε?+=-1

{}{}{}n

i i σσσ?+=-1

并记下每个单元的{}c σ和s ε' {}{}{}N

i i i σσσ?+=-1，以此作为（i +1）步的初始

状态，继续加载。

3) 讨论

上面采用的是最简单的纯增量法，并取其中一个增量步（i 步）内对m 值的迭代，最终确定加权平均的弹塑性矩阵ep D ????

（1） 采用加权平均的弹塑性矩阵ep D ????，在同一增量步内，对过渡单元的m 值往往要迭

代若干次，每次迭代都要重新计算单元的切线刚度阵，并重新组装总刚和解方程。显然求解过程比较复杂繁琐，由此增加了许多工作量，但从提高精度和加速收敛两方面是大有好处的。

实践证明：即使加载步长比较大，势必在这一步内新进入屈服的单元（过渡元）比较多，然而采用对m 迭代计算出比较准确的ep D ????后，仍能获得比较满意的结果。 如果不采用对过渡单元迭代的办法，则为了保证解的精度，必须控制每个增量步的大小，以保证每一步内新增加的塑性单元较少，否则将越来越偏离正确解，使求解失真，甚至发散。

（2） 迭代收敛准则：m 值的迭代从理论上讲要求两次值非常接近时方可结束。而大量实

践说明，一般迭代2~3次就可以，所以往往用迭代次数来控制即可，这是对一个过渡单元而言，而从整个结构来看，还要求在前后两次迭代中不再有新的塑性单元产生来决定是否可进入下一增量步的计算。

§3.4 采用纯增量法作弹塑性有限元分析的步骤

以下仅限于简单加载过程（无反复加卸载过程）和Mises 各向同性强化材料 1. 开始，输入初始参数（几何；材料性质，0s σ，H '；边界条件；外载荷P ） 2. 将外载荷一次加上作线弹性分析 {}{}{}max q εσσ→→→（Mi.条件） 如果 0max s σσ≤ 不存在塑性区则为弹性问题→直接输出结果 结束！ 否则 0max s σσ> 作弹塑性分析 3. 计算弹性极限{}e P

设 0max /s ασσ= ， 则 {}{}e P P α

=

并可输出弹性极限载荷{}e P 下的结果{}{}{}e e e q εσ、、。 4. 对剩余载荷{}{}{}r e P P P =-作弹塑性分析

如果采用等增量步格式，则将{}r P 等分为N 个增量步，即每一增量步载荷为：

{}r

P P N

?=

。下面5.中是对N 个增量步循环。

5. 在i 步上施加一个增量载荷i P ?。已知当前状态下（i -1步终），各单元的（or 高斯点）σ，

ε，s σ。判断三种类型的单元：1）弹性 2）塑性 3）过渡单元。对本增量步内所有过渡单元经过2~3次迭代得到合适的ep D ????，计算各单元的t k ，并集合所有单元，形成总刚T K ，求解{}[]T K q P ?=?得{}i q ? 得到第i 步的解。

{}{}{}1i i i q q q -=+? 和 {}{}{}1i i i εεε-=+?； {}{}{}1

i i i

σσσ-=+? 同时记录下各单元的当前状态。,s s σε'' 6. 直至全部载荷施加完毕，输出结果，结束

讨论:

1）以上步骤只考虑简单加载情况，变形与加载路径无关，因此没有采用加卸载准则。 2）在()t P ??增量步内方程是线性的，因为只要当前状态,σε已知，[]ep D 是可确定对

于每一步{}{}[]T K q P ?=?是线性的，但对于整个加载过程{}{}[()]K q q P =仍是非线性的，这就是分段线性化的思想。