有限单元法

弹塑性有限单元法

汽车车体冲压件通常使用弹塑性材料，在冲压过程中，这种材料的变形成形过程非常复杂，一般用刚塑性FEM与弹塑性FEM二种方法来评价整个冲压成形过程。在刚塑性FEM中，忽略弹性变形，仅将塑性应变作为计算指标。因此，在冲压成形过程中，当材料放置到模具上因自重产生的弯曲挠度，从模具中取出冲压件厚产生的弹性恢复等材料变形不能进行计算。因此，有人提出了根据刚塑性FEM的计算结果，再用弹性FEM计算其卸载过程，但是，刚塑性FEM很难正确地预测在冲压过程中产生的缺陷。弹塑性FEM可以再空间上时间上交替考虑弹性变形与塑性变形，从理论上讲可以正确地描述整个冲压过程，所以弹塑性FEM可以说是评价冲压过程的最好解析方法。 现有的弹塑性FEM，根据其时间积分方法的不同，可分为“静态显函数法”“静态隐函数法”和“动态显函数法”。讲加速度项加入平衡方程式求解的称为动态，反之，平衡方程式中不包含加速度项的解法称为静态。隐函数与显函数是常微分方程数值计算方法中的数学用语。显函数求方程的解不需要反复计算，而隐函数常微分方程求解时需要迭代多次逼近其解。显函数解法要求增分补偿不能取得太大，解析冲压成形过程需要较多的计算解析次数。隐函数解法通常可以保证应力平衡方程式成立，因而增分步长可以取得较大些，以减少解析计算次数。

各种弹塑性FEM的优缺点如下：

动态显函数法：该方法求解各节点的独立性运动方程以获得节点变

形，因而不需要组成刚度矩阵，即使单元划分得再细，节点再多也占用的计算内存较少，并且每一模拟步骤的计算速度也比其他方法快，因此可以计算对象的单元分割得很细。但是这种方法是用动态的冲击求解变形问题，时间增量需控制在10-6秒以下，要模拟一秒钟的冲压过程，就需要计算10^6次，实际上这种计算方法十分耗时，为减少运算时间，常常将物理意义不十分清楚的衰减项加入到运动过程，人为地将质量附以加权常数以减少模拟计算次数。此外，即使在方程式中加入了衰减项，应力值还是会发生振动，增加了弹性恢复计算的难度。

静态隐函数法，本方法的优点是考虑静态变形，模拟计算的每一步都能保证应力平衡方程式的成立，但是，在冲压件与模具相接触的极为复杂的冲压成形过程中。对应于每一个增分的反复计算，其接触状态部分发生变化，因而有可能造成计算结果的不收敛。所以静态隐函数解法存在的解不收敛的危险性使得计算结果的可信度大为降低。 静态显函数法 该方法将增分步长取得足够小，在保证每一步计算中应力平衡方程式中的线性关系条件下直接求解。因而不需要像隐函数数值计算法中反复迭代求解。也不会产生由于解的发散而得不到最终计算结果。但是，该方法必须将增分步长取得较小，因而需要较多的计算次数，模拟计算所需的时间较长，此外，该方法需根据静态应力平衡方程式组建整体刚度矩阵，对于单元多的大规模解析模型而言，需要较大的计算机内存。

综上分析，本方法采用静态显函数法数值模拟整个冲压过程。

二，弹塑性有限单元法的基本方程

在冲压成形数值模拟过程中，必须考虑以下三个非线性问题： 1材料非线性问题：板材在变形过程中，其应力与应变关系由弹性范围向塑性范围急剧变化，在塑性变形的范围内应力与应变的关系随着其变形会伴随产生所谓的冷作硬化现象。

几何学非线性问题：冲压成形过程中板材的变形很大，推到基本方程时以哪一个位置作为基准，会对基本方程的形式产生很大的影响。 接触非线性的问题：冲压成形过程中，板材与模具的接触位置，接触力，摩擦力等所谓的解除状态随着变形而变化。

为考虑上述非线性问题，考虑材料非线性，根据微小弹性有限塑性构建弹塑性本构式，考虑几何学非线性，采用常态更新基准位置的Updated Lagrange式构建刚性方程式，这样可以保证在增分步长较小时，求解增发量内的线性方程式，因而可以用显函数法求变分增量，本方法中，推到刚性方程式时采用不考虑加速度项的静力平衡方程式，因而称之为静态显函数法。实际的冲压成形过程中，可以说是静态或准静态变形加速度量的影响极其微小，因而本方法是切实可行的。 弹塑性本构式

假设材料各向同性，变形限定在微小弹性应变及有限塑性应变范围内，本构式可以表示为：

根据Hooke微小变形线性弹性理论：

ijCijklkl

e

e

Cijkl2G(ik

e

jl



12

ijkl)

„„„„„„„„„„

（1）

这里,G为杨氏模量，为泊松比，为Kronecker Delta。

ij

根据Hill各向异性屈服条件，塑性本构式表示为：

32(FGH)

1

{F(yz)G(zx)H(xy)2L

222

2yz

2M

2zx

2N}2y0

2xy

„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„（2）

这里F、G、H、L、M、N为各向异性材料店物理参数，假设这些参数不随塑性变形而变化的常数。用各向异性板的lankford参数R值的关系式如（4）所示。为柯西应力。

ij



'

ij

为偏差应力，



ij

为当量应力，为比例常数，为冷作硬化参数，其关系式为

p

p

()，

FH1ry

为相当塑性应变。其中函数为单向拉伸实验值。

GH1rx

， ，

NH



(rxry)(2rxy1)

2rxry

„„„„„„（4）

弹塑性应变关系式

DijDijDij

e

p

„„„„„„（5）

其中Dije为弹性应变速度，Dijp为塑性应变速度，由（1）-（5）式可以推



导

ep

出速度弹性本构式为：

ijCijklDkl

„„„„„„„„„„„„„„„„

„„（6）

C

ep

ijkl

2Gik



jl





12

ijkl



2'

2(H3G)

9Gijkl

''

„„„„„„„

„（7）

其中ij为Cauchy应力的Jaumann导函数，H'dH



d

p

为冷作硬化曲

0

线的斜率。弹塑性本构式中，在弹性状态或者卸掉外载荷时，本

ep

构式中矩阵Cijkl的杨氏弹性模量与泊松比均有常数决定。塑性状态时

1，弹性状态Cijkl

e中的偏差应力，相当应力，冷作硬化曲线斜率等

组成右项。因而，随着板材的变形，应力与应变的关系由线性逐渐变为非线性。

3、平衡方程式与边界条件

物体有变形初期状态进过时间T变形后，到达另一个新的平衡状态，如图1所示。考虑从时间t开始发生变形时，根据Updated lagrange式进行定式化。

在Updated lagrange式中，需适时更新基准位置，为与速度形的弹塑性本构式（6）直接关联定式，平衡方程式，边界条件等方程式也需要勇速度形来描述。

如图1所示，密度为p的物体，在单位体积力bi的作用下，单位面积上作用的表面力ti的作用下产生准静态变形。边界面S表面力作用部分st与位移i产生部分Su组成，bi、ti、ui均为指定值。 变形条件，用名义应力可以将total lagrange 平衡方程式表示为：

ijxi

bj

=0 „„„„„„（8）

上式中将时间t的基准位置x用现在位置x修正，以x作为基准位置，时刻t的速度平衡方程式可以描述为：



ijxi



bj0

„„„„„„（9）

上式即为update lagrange式的平衡方程式。

其次，考虑边界条件时，边界si上的表面力速度必须与内部应力速度相平衡。





ijnitj

„„„„„„„„（10）

由边界上位移速度的连续性条件





ii

即 i

i

„„„„„„„„（11）

4、虚功原理

由平衡方程式（9）边界条件式（10）、（11）根据Gauss发散原理，在时间t的基准位置时，updated lagrange式的虚功原理。





v

ji

(vi)dvxj





st

tividS

bvdv „„„„„„„„„„（12）

v

i

i



根据虚位移在v内一阶连续，并且不破坏Su上的几何学边界条件



avi

=0 „„„„„„„

„„„（13）

反之，假设虚功原理（12）式及满足虚位移全部条件的（13）式成立





v

ji

(vi)dvxj









st

tjvjdS

b

v

j

vjdv



ij

st(ijnitj)vjdSvxbj

i



vjdv0

„„„„„„„„„„

（14）

Su

边界以外



Si不受任何约束，因而式（13）的恒等式成立

ijxi

bj0



（v内） „„„„„„„„„„

（15）





ijnitj

（

Su

上） „„„„„„„„„„„„„（16） 式（13）以下式的成立为前提条件

vivi

（Su上） „„„„„„„„„„„„

（11）

因此，updated lagrange虚功原理方程式（12）（13）是平衡方程式（15）基准位置边界条件式（16）几何学边界条件式（11）的等价方程式。 虚功原理方程式（12）中，以公称应力速度ji表示应力，因而不能直接带入本构式（6），根据公称应力与Kirchhoff应力的速度关系。









LDDL

T

12

„„„„„„„„„„„（17）

T

其





中

L

vx

，

D(LL)

，

W

12

(LL)

T

，

WW

„„„„„„（18）

将式（17）代入式（13），并考虑对称变形，速度形虚功原理式（12）可以表述为下式：



ij2ikDkjDij

LLdvjkikij







st

tividS

bvdv „„„„„„

v

i

i

„（19）

Updated lagrange形式的虚功原理方程（19）的成立于应力应变无关，

因此，不仅是弹性问题，对于塑性问题也可以直接应用方程式（19） 刚性方程式

为将虚功原理方程式（19）应用于有限单元法中，需将其离散化，速度导数L和变形导数D为

LEv

，

DBv „„„„„„„„„„„„（20）

考虑在有限的时间增量t内可以保证线性关系时。

LijLijtwijwijt

，

dijDijt

，

ijit

，

„„„„„„„„„„（21）

根据式（21）将虚功原理式（19）离散化，可得出下述刚性方程式：

BC

v

E

ep

FBEGEdv

T





TT

tdS „„NbdvNvS

E

E

（22）

式（19）中以应力速度，弹塑性变形时，材料的体积变形很小，因而可以假定，将本构式代入式（19）

Kf „„„„„„„

E

E

E

„„„（23）

KB(CF)BEGEdv

E

T

ep

T

v

fNbdvNtdS „„„„„„„„„

E

T

v

e

„（24）

其中，N为形函数矩阵，当单元形状与节点数确定之后，可以求得E和B。

从式（23）单元的刚度方程，可以得到整体刚度矩阵式：

Kf

E

E

E

E1

E1

NENE

„„„„„„„„„

„„„（25）

求解上式，即可得到节点的变形位移增量。

图XX为本方法所使用的三维板壳单元，单元有四个节点，8个积分点，每个节点5个自由度，即三个移动自由度与两个转动自由度，单元内的变形位移速度表示为：



p

p

p

iNiN

h2

p2

p

pri



(1i2i1)Nr

p

(i1,3 r4,5) „„„

（26）

N

p



14



(111)(122)

D

D



p4

1

p

p





p5

2

p

NirN,r1,2,3

vr

p



根据上述形函数，可以导出速度倒数标量式为：



pkr

Lklk1l(N

)lN

rppkr1l

r



p

„„„„„„„„

„„„„（29）

为防止剪切形变产生过度约束而发生锁死现象，速度导数标量式中的横向剪切应变成分的相关项作如下近似：

L13(,,)(1)L13(0,1,0)(1)L13(0,1,0) 22

„„„„„„„„„

„（30）

L23(,,)(1)L23(1,0,0)(1)L23(1,0,0)

22

„„„„„„„„„

„（31）

最大步长

应用非线性弹塑性FEM静态显函数解法，刚性方程式成立的前提

条件为从时间t到tt的增量过程中，必须得保证其为线性关系。但是，冲压成形过程是塑性加压过程，是非线性关系的变形，要满足这个条件并不容易。在本方法中，对于材料非线性与几何学非线性，以保证在每一个增量步长范围内，节点力与节点位移见保证线性关系，为前提。决定每一次的增量步长。最大步长的决定方法如下： 1、 根据给定的增量步长值求解刚性方程式，求得位移增量，载荷

增量。

2、 对所有增量作如下检验： 从弹性状态向塑性状态过渡

弹性状态与塑性状态中的应力与应力关系大不相同，要保证刚性方程式中的线性关系，在一次增量过程中，不能发生从弹性状态转变为塑性状态的现象。

针对所有积分点，查看其是处于弹性状态还是塑性状态，求其应力正好达到屈服值时的增量比率Rep。

+如图XX所示，时间t时点A的应力状态位于弹性区域内，

超越屈服曲面，到达C点时，这时AC与屈服曲面交点B，

RAB

AC

点A的相当应力为：



32(FGH)

F(

y

z)G(

2

z

x)H(

2

x

y)2L

2

2

yz

2M

2zx

2N

2xy



12

相当应力值达到屈服曲面上点B时

3

2(FGH)F(yz)G(2zx)H(2xy)2L22yz2M2zx2N2xy12