冲压工艺数值模拟预测成型极限图的方法和意义

5.2.1 概述

冲压工艺数值模拟发展较早，相关力学理论比较单一，是纯粹的板壳弹塑性变形理论，因此模拟商业软件以及核心计算方法相对成熟，在航空、汽车等领域获得广泛应用。大量的工业应用促进了商业软件的发展，Pam-stamp、Dynaform、Autoform、Optris都是国际上广泛流行的商业软件。近年来，这些商业软件根据用户需求，在核心计算和辅助功能方面都取得了很大进展。现在的成型模拟模块可以对多种工艺全过程进行模拟，包括坯料在重力作用下的变形、压边圈闭合过程、拉延过程、切边回弹、回弹补偿、翻边、胀形、液压成型、弯管成型，还可以仿真超塑性成型过程。通过成型仿真模块，可以预测起皱、开裂、回弹等多种成型缺陷；可以预测成型力、压边力、液压胀形的压力曲线、材料性能评估等。在辅助功能上，Dynaform开发了冲压模具设计模块DFE、下料排样模块BSE和模具强度寿命分析模块DSA。

5.2.2 薄壳单元类型选择

薄板冲压成型过程中，薄板坯料逐渐变形成为复杂形状的三维零件。由于工件厚度远小于零件另外两个自由度的尺寸，因此它们属于薄板—薄壳类零件，可以用其中面形状和板壳厚度描述零件的几何特征。在板料冲压成型模拟中通常使用壳单元，这种单元的应力应变分布满足板壳力学理论的基本假设。普通三维实体分析中位移、应变、应力都有确切的对应关系。但是板壳力学分析使用了在薄壁条件下对应力应变做出某些假设，由此建立的广义位移和广义力可能不是严格的一一对应。例如在薄壳变形直法线假设下，承认截面上存在切应力，但是假设对应的切应变为零；承认壳体表面法向载荷，但假设厚度方向的应变和应力为零。使用不同假设的壳单元，所定义的广义力、广义位移和广义应变可能不完全相同^[88]，对应的本构关系也不相同。在板料成型数值模拟软件中包含了各种各样的单元，需要认真阅读其理论文本，掌握各类单元的力学特征和使用范围，然后再根据需要正确选择单元类型。

图4-47所示为退化的4节点壳体单元，它是由三维8节点六面体单元退化而来^[89]，单元几何由壳中面的4个节点坐标和壳体厚度确定。用ψ^α表示形函数，用X^α_i表示壳中面节点坐标，记V_i^α=X^αup_i-X^αlower_i，则单元内任意点的坐标x_i可以表示为

式中 ξ、η、ζ——单元某点的局部坐标。

图4-47 退化的4节点壳体单元

这个单元的每个节点位移有5个自由度（u_i），其中i=1，2，3对应3个平移，i=4，5对应两个转动的角位移。用U^α_k表示4个节点的5个位移分量（k=1，2，3，4，5，α=1，2，3，4），那么单元内任意点的位移u_i（ξ，η，ζ）可以由单元节点位移U^α_k表示。应用小转角、直法线假设，即每个节点在一个小的时间增量步中，认为向量V^α_i保持直线，且转角很小，壳内任意一点的位移可以表示为

这类壳单元有5个广义位移分量可以同时处理板料拉伸和弯曲变形过程。对于纯粹的胀形过程，不存在弯曲变形，因此可以用膜单元建立有限元计算公式，在膜单元中只有3个平动位移分量，它们对应板料中面的两个轴向变形分量和与之对偶的两个膜力，这里的转角位移和弯矩等于零。

5.2.3 材料本构关系与FLD

冷冲压工艺模拟需要使用室温的弹塑性本构关系，软件中各种单元都附带某种弹塑性本构方程。需要注意的是，大多数坯料都是轧制的金属板材，这类板材多是正交各向异性的，应根据板材实际情况选择合适的本构关系。板料的成型极限图（FLD）是冲压工艺数值模拟预测拉伸破裂和起皱的判断依据，它表征抵抗拉应力不破裂、抵抗压应力不失稳的能力。材料本构方程的参数、FLD，需要由板材厂家提供或者自己实验测试。对于热冲压，则需要使用专门的单元和专用的本构关系（例如超塑性本构关系）。大多数情况下，可以简化为各向同性材料，忽略弹性变形，考虑应力对应变速率敏感特性。材料参数随温度变化，热状态下板料更容易起皱但延伸率增加，目前很少见到热状态下板材的FLD，使用者需要自己实验测试。

5.2.4 有限元核心算法

求解式（4-2）有动态显示方法和静态隐式方法，当前主流冲压模拟大都采用动态显示方法，用静态隐式方法求解式（4-2）的软件较少。对这两类算法的比较表明^[90-92]，隐式算法是非条件稳定的，它在解决低速接触问题中更有优势^[88]，但是对于形状复杂的工件，因为节点太多和时间步长太小，内存消耗和计算时间会急剧增大，甚至造成收敛困难^[91]。显式算法虽然无此缺点^[93]，但是其惯性项和阻尼项不符合板料成型的静态实际，需要控制上模速度不要太高，并使用低阶四边形板壳单元。

普通金属板料冲压成型的上模速度较低，属于准静态变形范畴，可以忽略式（4-2）中的惯性项和阻尼项，这样式（4-2）就变成了静态方程，即

KU=F （4-5）

求解式（4-5）即静态隐式方法。其“隐式”是指U非显式表达，需要求解式（4-5）得到。因为材料非线性和几何非线性，必须用非线性迭代方法求解方程。由于冲压过程时间较长，零件几何尺寸很大，时间步长和单元尺寸很小，所以模拟计算量很大，还有可能迭代不收敛，这些问题造成静态隐式算法应用有限。但是因为这种方法更符合弹塑性理论，更符合冲压工艺实际，因此模拟结果更可靠。也正是因此仍然有不少人研究如何改善静态隐式算法。

动态显示算法本来是针对高速冲击而发展起来的计算方法，显然普通板料冲压不符合高速冲击的条件。但是这种算法不用迭代，模拟计算用时少，因此受到用户欢迎，使得这种算法成为板料成型软件的主流算法。获得位移矢量显式表达的关键是在时间维度离散时使用了中心差分方法。改写式（4-2）得到

方程右边第一项F为外力，第二相P=KU可理解为变形能对位移的微分，是内力。

应用中心差分法得

将式（4-7）代入式（4-6），得到

使用集中质量矩阵，则M为对角矩阵；设C=2ωM，ω为系统最低圆频率，C也是对角矩阵。这样U^t+Δt就可以显式地表达为U^t、U^t−Δt、外力F与内力P的函数。其中内力P可以由上一步的应力算出。为了保证计算稳定，时间步长Δt与单元特征尺寸L之比应足够小：

(www.xing528.com)

式中 ——弹性体积波波速。

因为高阶单元特征长度小，时间步长小，因此动态显示算法最好使用4节点退化壳单元，单元的质量被认为集中在4个节点的上表面和下表面上。

使用动态显式算法时，因为时间步长很小，冲压时间过长，为了提高计算速度通常采用提高上模速度方法，多数情况下实际上模速度在0.1～0.5m/s之间，但是模拟时选择压力机速度为4～8m/s。模拟与实际情况的差别增加了惯性力影响，这种影响将显著增加模拟结果的误差，因此上模速度不可过大，这就是动态显式算法的缺点。

逆算法（一步模拟法）。法国Batoz和Guo等人^[93]开发了基于形变理论的板料成型一步模拟方法（one step simulation）或逆算法（inverse approach）：这是一种粗糙的但是快速的板料成型模拟方法。这种方法的基本思路是用塑性变形的全量理论来分析从坯料到零件的变形过程。塑性全量理论的核心是在简单加载条件下，应力应变之间存在单值对应关系。据此建立有限元方法。通过有限元方法确定在满足一定的边界条件下工件中任意节点在初始平板毛坯中的位置，由此得到毛坯节点到工件节点的位移，进而得到工件的应力、应变的分布。Chung^[94]等人在希尔的最小塑性功路径理论的基础上，提出了理想成型（ideal forming）理论，并用于板料成型的分析，可以得到按最小塑性功路径变形的工件的应力、应变。这种算法可用于预测坯料展开图、零件开裂与起皱的可成型性分析或毛坯材料选择。

5.2.5 回弹预测与模具补偿

回弹是板料冲压不可避免的现象，它严重影响最终零件的形状尺寸，因此回弹预测和模具补偿一直是板料冲压工艺模拟的研究重点。但是非常遗憾，经过几十年的努力，回弹预测的精度一直不高。多届numisheet学术会议报告显示，对于相同的成型实验，几乎所有数值模拟给出的零件壁厚预测与实验符合都很好，但是回弹的预测结果与实验结果却差别很大，并且不同模拟给出的结果彼此也有很大差别。通常的板料冲压工艺模拟软件中，普遍采用动态显式算法计算成型过程，然后再用静态隐式算法计算回弹。在这种计算方法中，回弹变形的大小取决于成型过程完成时成型零件内的应力状态和应力水平。回弹预测误差与很多因素有关，但是更与现行冲压工艺模拟的基本方程和基本算法密切相关。现行所有冲压工艺模拟算法，无论是静态隐式算法还是动态显示算法，其所解方程[式（4-2）]的基本变量是节点位移。改进方程解法固然可以提高位移的精度，但是作为位移微分的应变变量，其精度比位移精度低一个量级，由于应力与应变具有相同量级的精度，因此应力的精度比位移精度低一级，自然应力卸载引起的弹性应变与应力具有相同的精度级别，所以回弹变形的计算精度必然低于位移的预报精度。基于这个认识，如果根据最小余能原理^[95]推导出以应力为基本变量的方程，并用它作为有限元的基本方程来模拟成型过程，然后再由应力算出回弹，这样有可能显著提高回弹的预测精度。

模具回弹补偿的含义是，根据预测出的零件回弹变形补偿模具型腔形状，再用修改后的模具进行成型和回弹计算，看回弹后的零件形状是否和零件的设计形状一致。计算实例表明，这种模具补偿方法是有效的，有些商业软件已经增加了回弹计算和模具补偿的模块。

5.2.6 拉延筋模型

在板料成型中，在模具型腔外围设置拉延筋是控制材料流动、防止起皱的主要技术措施。准确模拟拉延筋阻止材料流动的作用对预测拉延成型起皱、开裂缺陷非常重要。因为拉延筋尺寸较小、形状复杂，要精确模拟板料与拉延筋的接触作用，就需要把单元网格尺寸划分得非常细小，这会显著增加计算量。为了不过多增加计算量，目前多数模拟软件的做法是将复杂形状截面的拉延筋抽象为一条能够对板料流动施加一定阻力的拉延筋曲线。计算时这个拉延筋线也被离散为多个线单元。应用二维有限元模拟方法可以建立拉延筋截面几何参数与其力学特征的关系，软件根据这种关系设置板料流动的阻力。冲压工艺模拟之前，用户要根据实际拉延筋截面积、长度形状和在模具上的位置在模具CAD图上设置拉延筋。

5.2.7 板料热成型的工艺模拟

板料热成型包括热冲压成型、液压、气压热胀性、超塑性成型以及近年来发展起来的高强钢热冲压成型。材料在高温状态下的力学性能与室温状态下有很大区别，主要表现在：流动应力显著降低，弹性应变显著降低，流动应力表现出明显的应变速率敏感性，合适温度与细晶条件下还会出现超塑性，即使不是超塑性状态材料延伸率也大幅度提高。温度不同材料的力学性能也不同。根据上述分析，板料热成型的模拟需要考虑两个问题，首先是使用热状态下的材料本构关系，因为热状态下弹性变形很小，因此成型过程模拟可使用刚塑性本构关系，可不考虑回弹对零件形状的影响。但是需要考虑零件出模重力引起的零件塑性变形以及出模后降温引起的塑性变形。因为材料流动应力低，变形幅度大，在主要成型过程中弯曲的作用远远低于板料拉伸的影响，因此热冲压成型模拟可以忽略弯曲作用，使用膜单元可以大幅度减少计算量。对于非均温的成型情况，工艺模拟需要增加传热计算模块。因为零件壁很薄，成型过程中工件向空气和模具传热很快。因为传统的膜单元和壳单元都不考虑厚度方向的变化，因此需要编制专用的传热模块处理薄壁零件的传热分析。对于高强钢热冲压成型，成型过程模拟不仅要进行零件变形模拟，还要模拟零件成型后的淬火过程以及零件材料从奥氏体到马氏体的相变过程。正确预测零件各部位温度以及温度变化速率是准确预测材料相变和微观组织的关键。因为存在零件急冷和淬火过程，零件强度随温度降低急剧上升，因此弹性变形、弯曲效应、回弹以及冷缩热胀对零件最终形状的影响不可忽略，为此应使用随温度变化的弹塑形本构关系，使用包括弯曲和膜力效应的壳单元。这是一个最复杂的板料冲压工艺模拟，它包括冷热结合、宏微观耦合、弯曲和拉伸结合。PaulÅkerström在他的博士学位论文中建立了一个热力模型并插入到显式有限元代码dyna3d中，用来模拟薄壁件成型和淬火。它用一种交错的方法处理热力耦合分析，其中每个分析使用不同的时间步长。所使用的热壳单元沿厚度方向和平面方向假设线性温度分布，包括接触热传导。单元使用温度相关的非线性各向同性弹塑形本构关系。

5.2.8 算例

（1）应用逆算法，计算出的冲压件毛坯展开图如图4-48所示。

（2）使用动态显示算法，割草机罩拉延成型的数值模拟。图4-49显示出了零件壁厚分布的模拟结果，其中蓝色预示可能出现开裂缺陷，左下边台阶拐角处显示出了一些折叠缺陷。

（3）使用动态显示算法，汽车车身结构件的冲压成型数值模拟。如图4-50所示，工艺模拟的目的是优化模具设计，消除拐弯处的起皱缺陷，成型设备为双动压力机。

（4）高强钢B柱热冲压成型工艺模拟，图4-51所示为壁厚分布。

（5）高强钢零件热冲压成型工艺模拟。材料为高强钢Boron02，图4-52所示为铁素体体积分数在零件内的分布图。

图4-48 精确预测翻边成型冲压件的坯料形状

图4-49 拉延成型后的割草机罩壁厚分布

图4-50 结构件冲压模拟壁厚分布

图4-51 热冲压高强钢B柱壁厚分布

图4-52 热冲压高强钢Boron02零件铁素体体积分数分布^[96]