0 前言
目前，板材冲压成形过程的数值技术已取得了很大的进展。根据有限元模拟技术预示的成形载荷、板材的几何变形、应力应变分布和加工条件，调整模具的几何形状、材料等级或边界条件，从而改进模具设计。绝大多数冲压过程应用三维模型来分析，然而这种总体分析将面临许多技术障碍，如大量的计算时间，繁杂的工具数据前处理以及未知的边界条件等。然而，复杂冲压件的许多局部变形常常可以近似地采用平面应变的截面分析方法对成形过程加以模拟，而且这种截面分析方法所需的数据量极少，在有效的精度范围内能更快、更稳定地对模具设计的合理性给出定量的评价。早在1980年，hughes和liu就给出了一个总的二维问题的非线性有限元公式。经过多年的研究，1991年keum、wang、saran和wagoner编写的有限元程序应用线单元模拟了任意形状的平面应变截面深拉或胀形的成形过程。
本文基于有限变形虚功率增量型原理的弹塑性大变形有限元理论，建立了比线单元精度更高的八节点四边形单元平面应变状态的截面分析模型，用三次b一样条曲线描述凸、凹模模具截面的几何形状，界面滑动应用库仑摩擦定律，进而实现模具表面与板料之间的接触判断。该模型可在很短的时间内，经过少量的数据处理，得到较理想的板料成形的模拟结果。文末用该模型数值模拟了金属板料成形过程中的塑性流动规律，并将计算结果与试验结果或三维有限元程序的计算结果进行了比较，证实了这种平面应变的截面分析技术的正确性、可行性和高效性。
1 有限元模型的建立
采用逐级更新lagrange法，在xi坐标下以t时刻构形为参考构形的虚功率增率原理为
(1)
式中 v，a——参考构形的体积和表面积
——参考构形的面积力率和体积力率
——第一类kirchhoff应力率
假定塑性变形体积不可压缩，则有
(2)
式中 σij,j——xi坐标下cauchy应力σij的jaumann速率
σij,j与cauchy应变速率
的本构关系可一般性地写成下式
(3)
对于平面应变模型，不妨假定材料为各向同性，则dijkl可以写成
(4)
式中 σ′ij——cauchy应力偏量
h′*，μ*,e*——材料参数
取不同的表达方法，将对应于不同的塑性本构理论，本文采用的是“拟流动角点理论”。将式(2)和式(3)引入式(1)，得
(5)
式中
(6)
将式(6)代入式(5)并离散化，可得单元刚度方程。
2 接触判断与摩擦约束
模具的截面曲线是由非均匀有理b样条(nurbs)曲线来描述的。给定控制顶点位置矢量di，i＝0，1，…，n，次数k及确定节点的参数矢量u＝［u0,u1,…，un+k+1］，就定义了一条k次b样条曲线。如若给出曲线定义域内一参数值u∈［ui,ui+1］
［uk,un+1］，欲计算该b样条曲线上对应一点位置矢量p(u)，采用德布尔算法
(7)
模具的截面线是由若干段非均匀有理b样条曲线来描述的。每段b样条曲线采用疏密不同的点，即可以更好的描述模具的危险截面，又可以提高接触判断的效率。
所谓接触判断就是求出板料上的节点与模具表面的接触点。根据模具不可穿透原则，对于计算后进入模具的点，必须拉回到模具表面上来。截面分析采用八节点四边形单元，节点的拉回方向为板料的法线方向，并求出平均外法线与b样条曲线的交点。首先要确定外法线与哪一段b样条曲线相交，称这一段b样条曲线为目标曲线。在每一步接触判断中，具体的工作就集中在求节点的外法线与目标曲线的交点。再将节点坐标调整到交点处，就完成了接触判断的几何调整。
对于已经接触到模具的节点，增量步内的位移不再是自由位移，该节点必须沿着模具表面滑动，在模具表面的截面上，板料节点i的x方向的位移增量δuxi与y方向的位移增量δuyi存在下列约束关系
δuyi＝kiδuxi+δy(8)
式中 δy——该增量步内模具的位移
ki——模具截面线与节点i的接触点的斜率
如图1a所示。
图1
在冲压问题中，板料与模具接触并发生了界面滑动，必然存在界面滑动摩擦。现建立简化的力学模型，如图1b所示。t为板料节点n在模具的接触点与模具相切的方向矢量，即为b样条曲线的一阶导矢，由于节点n沿冲头滑动，则滑动切向方向矢量t与节点n的滑动速度无关。采用库仑摩擦定律，有
ff＝μfn(9)
则在笛卡尔坐标系下，节点n的切向摩擦力为
ft＝fft(10)
且可表示为
ft＝(-fn2，fn1)(11)
式中 fn2，fn1——ft在x轴和z轴上的投影
3 截面分析数值算例
3.1 方冲头圆坯料成形过程的模拟
方冲头冲压圆形坯料(如图2)，现截取截面a进行截面分析。材料的物理参数见表1。
图2 方冲头成形的结构与尺寸和板料平面视图
表1
| 库仑摩擦系数μ | 弹性模量e/gpa | 泊松比γ | 屈服点σs/mpa | 硬化准则 | 硬化指数n | 板料厚度δ/mm |
| 0.15 | 200 | 0.3 | 147.9 | σ＝cεn | 0.2 | 1.0 |
对于胀形问题，板料四周为固定约束，取单元数为120，节点数为603。计算结果与文献［6］采用三维厚曲壳单元的模拟结果几乎重合，而前者模拟的计算时间仅为后者的7%，计算得到的应变分布如图3a所示。对于深拉延问题，板料四周受拉深筋阻力约束。采用平面应变等效拉伸筋阻力模型［7］，取单元数为160，节点数为803。计算结果与试验结果［2］相比较，其中应变分布图如图3b所示。计算结果表明：截面分析方法不仅具有高效率，而且可以获得令人满意的模拟结果。
 
图3
注：■——上表面值 ×——下表面值
— ——平均值 □——试验值
3.2 发动机油底壳横截面成形过程的模拟
油底壳是汽车覆盖件冲压成形中典型的深拉延件。取油底壳的一个危险截面，位置如图4a箭头所示，截面的几何形状如图4所示。坯料的物理参数如表2。
 
图4
表2 |
| 库仑摩擦系数μ | 弹性模量e/gpa | 泊松比γ | 屈服点σs/mpa | 硬化准则 | 硬化指数n | 板料厚度δ/mm |
| 0.1 | 207 | 0.3 | 147.9 | σ＝kεn | 0.2 | 1.0 |
油底壳的成形过程为拉深过程，边界条件为拉伸筋阻力。仍采用平面应变等效拉伸筋阻力模型［7］，将坯料划分成500单元，节点数为2503。若采用三维曲壳单元，需划分整个油底壳成形坯料为16570个节点。在dec/433工作站上分别进行计算。两者的计算时间相差近20倍，前者仅需十几min即可完成全部模拟过程。计算得到的等效应变分布如图5a所示。板料厚度分布如图5b所示。将试验的油底壳冲压件采用激光切割方法按计算选取的截面切割开来，图中的.号为实测得到的变形终了的厚度。
 
图5
——截面分析方法的模拟结果
——三维壳单元的模拟结果.——试验点
4 结论
基于有限变形虚功率增量型原理建立的八节点四边形弹塑性大变形二维截面分析有限元模型，模拟了方盒典型件胀形和深拉延及油底壳的深拉延成形过程。通过计算结果与试验结果相比较，可以得出以下结论：
(1)对边界条件明显且近似满足平面应变的截面，本方法可以比较精确的模拟其成形过程；对不完全满足平面应变条件的截面，也可以得到较严格的结果。
(2)本模型有良好的求解效率，可在相当短的时间内完成油底壳横截面的分析过程，从而大大提高了板料成形性的分析效率。
(3)采用八节点四边形单元，可以考虑板料上下表面的不同应力、应变情况，在保证求解效率的前提下，其计算精度要高于膜单元的二维截面分析。 |
