刚塑性有限元数值模拟中产生误差的原因及改进方法

刚塑性有限元数值模拟中产生误差的原因及改进方法

作者：陈军邓清张卫刚

1 引言

塑性加工过程的有限元数值模拟，可以获得金属变形的详细规律，如网格变形、速度场、应力和应变场的分布规律，以及载荷-行程曲线。通过对模拟结果的可视化分析，可以在现有的模具设计上预测金属的流动规律，包括缺陷的产生（如角部充不满、折叠、回流和断裂等）。利用得到的力边界条件对模具进行结构分析，从而改进模具设计，提高模具设计的合理性和模具的使用寿命，减少模具重新试制的次数。在制造技术飞速发展、市场竞争日益加剧的今天，塑性加工过程的计算机模拟可在模具虚拟设计、制造阶段就能充分检验模具设计的合理性，减少新产品模具的开发研制时间，对用户需求做出快速响应，提高市场竞争能力。由此可见，金属成型过程的有限元模拟已是模具计算机集成制造系统中必不可少的模具设计检验环节。

金属成形工艺分体积成形和板料成形两大类，相应地，用于分析其流动规律的有限元法也分为两类，即：刚塑性、刚粘塑性有限元和弹塑性有限元。体积成形中的挤压成形和锻造成形在实际生产中应用很广，中外学者在这方面进行了很多研究，其中二维模拟技术已相当成熟，三维模拟是目前的世界研究热点。刚塑性、刚粘塑性有限元模拟能否对模具设计的合理性做出可靠校验，取决于模拟的精度和效率。作者结合从事二维塑性有限元模拟的经验和当前的三维塑性有限元模拟系统开发的实践，对刚塑性、刚粘塑性有限元模拟过程中产生误差的原因进行了全面的详细分析，并提出相应的解决方法，同时以具体实例说明。

2 刚塑性、刚粘塑性有限元模拟中产生误差的原因及改进方法

2.1 刚塑性有限元法求解的数学基础

刚塑性有限元法是假设材料具有刚塑性的特点，把实际的加工过程定义为边值问题，从刚塑性材料的变分原理或上界定理出发，接有限元模式把能耗率表示为节点速度的非线性函数，利用数学上的最优化原理，在给定变形体某些表面的力边界条件和速度边界条件的情况下，求满足平衡方程、本构方程和体积不变条件的速度场和应力场。速度场的真实解使以动可容速度场建立的能量泛函取极小值。但所得到的塑性力学的微分方程组一般不能用解析法求

解，常采用数值解近似，而采用数值解，则会出现各种误差。误差取决于所用的数值方法。下述处理方式易引起系统误差。

2.1.1时间和空间的离散化

刚塑性有限元分析的对象是一个非线性变化过程，即材料应力-应变关系的非线性和几何边界条件的非线性。解决这一问题可以采用线性小变形拟合非线性大变形，如图1所示。每一个小变形过程的选取须足够小，同时兼顾逼近的精度和效率。对于刚塑性材料来说，每个加载步长△S，即（△t·Vdz）应小于某一规定值（坯料当前高度的1.0%）[1]。作者认为，三维模拟的位移加载步长不应超过边界单元最小边长的1/4，以减缓接触边界非平面性的程度，更好的模拟金属的流动规律。

另一方面，所分析的变形材料是一个空间连续体，而有限元法的思想是把无限的连续用有限的连续近似，即用网格离散变形体。变形场量在单元内连接，这时将产生离散误差，一般地，有限元网格划分得越细，引起的离散误差越小。但是刚塑性有限元分析必须同时兼顾精度和效率，所以单元不可能过于细化，采用局部网格细分可以满足要求。如图2所示是作者对方坯反挤工艺三维刚塑性有限元模拟时采用局部网格细分后的网格变形。文献[2～4]中采用自适应网格离散变形体。

但对于一定的变形过程，即使网格划分得再细，误差仍然存在，这就是形函数误差，也叫做插值误差，是指单元的形函数反映真实变形规律的能力。所以高阶单元反映材料的变形能力较好，但计算效率偏低；如果同时考虑计算精度、效率以及可视化问题，二维问题宜采用四节点四面体等参元，三维问题宜采用八节点六面体等参元。三维问题考虑到网格再划分的方便，可采用十节点四面体等参元[5]，但可视性较差，文献[6]中采用多种单元类型的处理方法，可以减小一定的插值误差，但系统实现繁杂。

2.1.2解析式的数值计算

对于刚塑性有限元分析，单元的刚度矩阵由下式求出：

式中[B]——形函数矩阵

上式无法直接积分求出，必须进行插值，最有效的方法是高斯插值。二维问题可采取形心和另外四个高斯积分点插值，三维问题采取形心和八个高斯积分点插值，误差极小。例如三维问题采用8个高斯积分点插值，上模速度若设为1.0，与27个高斯积分点相比，节点速度的绝对误差小于10-5。

与单刚的求法相似，摩擦力边界条件对单刚和载荷列阵的贡献可采用辛普森插值的方法，二维问题的线性边界取5个积分点，三维问题的类似平面边界中的两个局部坐标方向各取5个积分点，即可满足精度。

在采用计算机进行数值计算时，还会产生截断误差和舍入误差，误差的积累可能产生病态的

线性方程，无法获得精确解。应根据采用的计算机硬件和编程语言，合理选择数值精度。

2.2 刚塑性有限元法求解的力学基础

刚塑性有限元求解时，假定材料各向同性，而且体积不变，即忽略弹性变形。实际上，各种原材料由于生产方式的不同和材料成分的不均匀，多呈现一定程度的各相异性，一般在分析轧制工艺和板材成形时考虑材料的各向异性，而对于其它体积成形问题多不考虑，这种做法仅产生较小的偶然误差。另一方面，因为在刚塑性有限元分析时，处处体积不变的条件不易满足，常采用三种方法实现近似满足，即Lagrangian乘子法、罚函数法和泊松比接近于0.5法。第一种方法对力的求解最精确，但线性方程组的求解量大，后一种方法的求解的精度取决于材料的体积可压缩率与实际的接近程度，常用于分析特殊成型过程，如粉末烧结成形。而罚函数法由于其求解效率高而应用最广。但这种方法是用惩罚因子与单元平均等效应变速率的乘积近似作为单元的静水应力，这也是产生误差的主要原因。

当模具的塑性有限元法是根据变形能量的泛函进行变形场量的求解，其实质上是一种稍微精确的上限法，上限法所引起的误差是存在的。

刚塑性有限元模拟时，一般不考虑体积力和惯性力，前者引起的误差足够小；对于惯性力，当材料高速成形时，能量泛函中需考虑惯性力做功，此时材料在模腔内的充填规律与低速成形时不大相同，尤其对于刚粘塑性材料[8]。

温度条件的简化也是产生误差的原因。对于冷挤压问题，从一开始材料接触模具边界的成分相当多，剧烈的摩擦容易引起局部温度的升高，变形载荷和应力分布也发生相应的改变。因此耦合模拟是必要的。对于热锻成形过程，尤其对于多工位成形，温度变化相当大，也需要耦合模拟。

刚塑性有限元模拟时，材料的应力-应变函数关系（多采用各相同性强化模式，而极少采用随动强化模式）是一项极其重要的初始条件。这一具体的函数关系式必须要由专门的实验确定，它的精确性决定了对实际成形过程模拟的近似程度。文献[9]中的实验表明，变形力和应力的值具有与应力-应变曲线同等的精度，流动特性对应力-应变曲线的偏差则不敏感，该文还尤其强调了应力-应变曲线斜率的精确性要求。

2.3 物理模型在模拟过程中的技术处理

刚塑性有限元理论在具体的系统实现时有诸多技术问题。这些问题的处理直接决定了模拟系统的精度和效率。若处理不恰当，则会产生误差。产生误差的原因主要有以下几方面：

1）初始速度场的生成刚塑性有限元模拟开始要生成精度合理的初始速度场。生成方法有多种[10]，实践证明，适用于二维和三维任意边界以及速率敏感材料的应属有限元法线性化本构关系和直接迭代生成初始速度场。三维问题中由触节点的局部坐标方向和坯料放置不当而生成的速度场可能无法用于加载迭代运算，尤其在多工步成形时

2）摩擦力边界条件的施加在模具的作用下发生塑性变形的金属与模具表面之间存在着剧烈的摩擦，这在数值模拟过程中，是一项很重要的边界条件。该边界条件的简化直接影响模拟系统的可靠性。目前有多种摩擦力的数学模型，应用最为广泛的是C.C.Chen和S.Kobayashi 提出的反正切函数模型[11]，可以有效地处理各种情况的边界，尤其是具有分流点的属流动，其模型表达式为：

由于上式中α的取值极小，所以一般边界自由节点一旦与模具边界接触，摩擦力几乎是最大值（m·k）。这可能与实际存在误差，较完善的方法是附加一项考虑相对滑动速度的修正系数。另外，三维模拟时若采用六面体八节点等参单元，单元与模具的接触部分是四个节点组成的直纹面（即四个节点可能不在同一个平面内），构造与该直纹面对应的虚拟接触平面（关键是其法线方向）至关重要。

3）收敛准则速度场迭代是否收敛常采用的判定准则有三种：节点相对速度误差泛数收敛，节点相对力误差泛数收敛和一阶功率泛函收敛。如果衰减因子β选取合理（收敛因子的自动调整是决定计算效率最重要的因素），则第一个条件最先满足。应当指出的是，收敛精度不能取的太低，否则容易引起以后加载步速度场迭代的发散。同时，对边界接触模具节点的脱模法向力的判断也不准确。

4）刚性区的处理刚塑性有限元法将变形体的弹性变形区视为刚性区，这种区分的准则是极限应变速率ε0。当单元的平均等效速率小于极限应变速率时，则认为该单元为刚塑性区。刚性区变形功率的处理有以下两种方法：

对于不同的锻压工艺过程的模拟，极限应变速率的取值有所不同。一般对于变形成形过程，初始变形时刚性区较大，边界约束较小，极限应变速率的取值稍大一点，为10-2～10-3，稍后可逐渐减小极限应变速率的值；而对于挤压变形过程，边界约束较多，极限应变速率的取值可以稍小一点，为10-3～10-4。如果极限应变速率的取值不合适，则速度场迭代难以收敛，甚至得到完全不同的速度场，模拟的经验就显得特别重要。

文献[12]中提出一种新的不必区分刚塑性区的处理方法，可以减少刚塑性区的区分而引起的误差。

5）动态边界的自动处理由于刚塑性有限元模拟的是金属的整个变形过程（非稳态成形过程），因此网格形状的动态刷新是模拟得以进行下去的保证，动态边界的自动处理正是实现这一功能的前提条件。这一处理技术主要包括：①边界自由节点接触模具时间步长增一的确定；②时间加载步长增量确定后的位置刷新；③边界触模节点位置的调整；④边界触模节点的脱模判断。

时间加载步长增量若以下面的式子确定：

式中△ti——边界自由节点接触模具的时间加载步长增量

△t1——保证速度场迭代不发散的时间加载步长增量

△t2——保证体积损失不超过1.0%的时间加载步长增量

则可以保证较高的位置刷新精度，但计算效率低，有时在模拟模具的局部形状比较特殊的成型过程时，时间加载步长增量经常为零，出现迭代过程中的死点。为了避免这种情况，可以采用下式确定时间加载步长增量：

这样即可以保证刷新的精度，又可以提高模拟的效率。作者在进行连杆锻造过程滚挤预成型工步的三维刚塑性有限元模拟时，在初始的加载步内，时间加载步长增量的确定先采用式（5），随后再采用式（6），实践证明效率较高。图3为连杆多工位锻造成型中的预成形滚挤工艺三维仿真的变形网格图。

6）网格再划分刚塑性有限元模拟的是金属的大变形过程，由于摩擦力的作用和模具表面形状的变化，容易引起不均匀的网格变形，最终导致网格畸变。为了使模拟得以继续进行下去，必须重新划分网格，然后采用新的网格系统继续进行模拟。二维和三维自动重新划分的网格应尽可能减小离散误差，同时对可能出现缺陷的位置采用自适应网格加密。另外，由于历史量的插值在边界处误差较大，可采用外推法进行插值[13]。刚塑性有限元模拟应尽量减少网格再划分的次数，因为每一次网格再划分都要降低一次精度。

3 结论

本文系统全面地分析了刚塑性有限元数值模拟中产生误差的原因，并提出了相应的改进方法。应当指出，有些处理方法引进技术起的误差只能分析其趋势，目前还无法进行量的预测。

同时，物理模型的有些机制还未真正揭示出来，简化方法还有待于改进。全面分析刚塑性有限元模拟产生误差的原因，提出有效的改进方法，可以提高模拟结果的可靠性，使其真正起到在模具制造之前检验模具设计合理性的功能。(end)