静力隐式时间积分算法假设板料冲压整个过程均处于类似平衡状态，其每一时间积分步均保持平衡条件.把物体离散为有限元单元后，将所有的节点方程集合起来就可以得到：
Ku=F
式中K为总刚度矩阵；u为位移向量；F为节点载荷向量。
在隐式算法中，对于第i个给定的加载增量，用Newton-Raphson迭代法，需要求解下面的方程：
板料冲压成形及回弹有限元模拟分析
摘要
回弹是板材冲压成形过程中不可避免的普遍现象，直接影响到冲压件的尺寸精度和零件最终形状。本文利用ANSYS/LS-DYNA有限元软件中的非线性动力的显式、隐式连续求解功能，模拟了板料冲压成形过程与卸载后板料回弹变形的全过程，得到了成形过程中任一时刻各处Von-Mises应力云图和应变值及卸载后板料的回弹结果，帮助我们更好的认识分析板料冲压成形以及回弹过程中物质内部的变化。
ANSYS／LS-DYNA是一个通用的显示非线性动力分析有限元程序，近年来开发的板料冲压成形数值模拟新功能，取得了很大成功。通过计算，可以观察板料冲压成形过程中的变形状态、应力状态和壁厚变化，了解可能出现的起皱和开裂现象。此外，ANSYS／LS-DYNA程序具有显式、隐式求解功能，用显式求解模拟动态成形过程，然后用隐式求解模拟线性回弹变形。这不仅能够模拟材料变形过程，而且也能较好地计算回弹现象，比较准确地得到材料最终成形状态。
关键词：板材冲压，回弹，非线性有限元分析，数值模拟
Sheet metal stamping and rebound finite element simulation analysis
Abstract
The rebound is inevitable common phenomenon in sheet metal forming process, a direct impact on the final shape to the dimensional accuracy of the stampings and parts. In this paper, the nonlinear dynamic finite element software ANSYS / LS-DYNA explicit, implicit sequential solution function to simulate thespringbackdeformation of the sheet after sheet metal stamping process and uninstall the whole process, forming process at any time throughout the Von-Misesstress cloud and strain and after unloading sheetspringbackresults, help us to a better understanding analysis sheet metal stamping and rebound process material internal changes.
2.2 有限元法
在拉格朗日坐标系中，对应于某一瞬时t，具有总体积V，内作用面 ，位移边界 与外力边界 的相互作用系统应满足虚功原理，即
式中 是密度； 是位置坐标； 是加速度； 和 分别为虚位移及其导数的增量； 和 是体积力与面力。
将体积离散化为n个单元，为了保持4面体单元的抗畸变能力，同时又使计算精度得以保证在4节点4面体单元的基础上，建立了8节点6面体超级单元。每个单元的节点数为8，位置坐标 可用形函数 与节点位置 的插值表示为
4.1 板料变形分析
图4-1板料变形图和图4-2位移矢量的等值线图为时间t=0.01s板料的变形程度,从图中可以清楚的知道在冲压过程结束后，板料的下部变形程度最大。
在板料上可能与模具曲率变化较大的区域相接触的部位也应细化板料网格。则板料网格划分中存在网格密度分布问题。在板料冲压成形及回弹过程中，板料会发生大位移、大转动和大变形，在冲压成形前预先确定扳料上的网格密度分布情况是非常困难的。目前有限元软件中采用自适应网格划分技术解决此问题。所谓自适应网格划分。是指板料在接触过程中遇到变形较剧烈的部位自动进行局部区域的网格细化，以提高有限元在这些部位的计算准确度．网格划分如图3-4所示。
式中 为当前的切向刚度矩阵； 为位移纠正向量； 为外力向量； 为内力向量。
经Newton-Raphson迭代计算后的位移增越为：
算法中的要点如下：
l）对于每一个加载增量都要进行迭代，直至收敛为止。
2）每一节点都要达到力平衡和力矩平衡。
3）每一节点都要满足接触条件。
4）相对于位移增量来说，位移校正量是很小的。
图3-2SHELL163 单元几何图示
SOLID164是用于3维的显式结构实体单元，如图3-3所示，有8节点构成。节点在x、y、z方向的平移、速度和加速度的自由度。这个单元只用在动力显式分析，它支持所有的许可的非线性特性。
图3-3SOLID164 实体单元几何特征
3.3 网格划分
有限元网格的划分是数值模拟技术中的重要环节，网格划分是否合理在一定程度上会影响到计算精度和计算时间。在网格划分过程中，网格的形状、尺寸和位置都是需要考虑的问题。对于板料冲压成形回弹问题，坯料是典型的大变形构件，又是成形分析过程中的从接触面，因此必须采用比摸具单元较精细的网格模型，单元的形状尽量采用四边形或者六面体。而对于模具，因在分析过程中认为是刚性体，若形状较简单，可不划分单元直接采用其解析刚性面进行分析，若模具形状较为复杂，因对其进行网格划分，在模具圆角及形状变化较大的区域采用较密的网格，以保证精确的描述模具型面的几何形状．
2 基本原理
2.1 弹塑性材料的应力应变关系
板料冲压成形属于弹塑性大变形行为，弹塑性材料进入塑性后应力应变不再存在唯一的线性对应关系每一时刻的应力~应变状态都是此前的加载路径和应变状态综合作用的结果应引入增量理论。本文采用Mises屈服准则即
式中 是材料单轴应力状态时的屈服应力； 是应力偏张量， 。假设材料发生的是各向同性强化后续屈服函数应为 ，即
形函数 定义为
式中 ， 与 是单元局部坐标，t是时间变量。
完整的回弹数值模拟分为2步：冲头落下板料成形的过程模拟和去除冲头等工具板料回弹的过程模拟。前1步计算是回弹过程模拟的基础，为其提供应力、应变等数据，该步计算的准确与否将直接影响回弹计算的精度，2步计算可分别采用不同的算法。
数值模拟主要有2种算法：静力隐式算法和动力显式算法。2种算法的区别在于静力隐式算法需要构造和求解刚度矩阵，每一步迭代都要进行接触判断，往往会出现迭代不收敛的情况，计算速度慢，但计算精度高；而动力显示算法不需要构造和求解刚度矩阵，因此该算法的最大优势是计算速度快。实践证明：动力显式算法计算的应力是不准确的，而残余应力恰恰是造成回弹现象的关键所在，因此认为在回弹模拟的第1步采用动力显式算法是不合适的。虽然动力显式算法在计算以接触为主的问题时速度相当快，但是当用于不存在接触并且以低频率方式变形的柔性结构时（最典型的是回弹模拟的第2步计算），为保证算法的稳定性，分析时间常常比模拟成形过程所要的时间还要长。质量加大技术虽然可以加速成形模拟过程，但它对回弹第2步计算是毫无帮助的，因为加大质量的同时也加大了达到静平衡的时间。为预测回弹后零件的形状，大部分学者建议采用静力隐式方法。
本次分析选用单元类型是：板料模型采用显式壳单元SHELL163，且使用能容忍弯曲能力的Belytschko-Wong-Chiang算法；凸、凹摸则采用SOLID164体单元。
SHELL163 是一个4节点显示结构薄壳单元,如图3-2所示，有弯曲和膜特征。可加平面和法向载荷。单元在每个节点上有12个自由度: 在节点x,y,和z方向的平动,加速度,速度和绕x,y,和z轴的转动。该单元只用于显示动力学分析所有非线性特性.
Keywords: sheet metal stamping, rebound, nonlinear finite element analysis, numerical simulation
1 引言
金属板料冲压成形是压力加工的重要组成部分，薄板金属在压力作用下由模具引导成形的过程是一个十分复杂的物理过程，由于模具几何结构尺寸、接触摩擦和压边力等因素的影响，在金属的成形过程中常发生起皱现象。为了得到高质量的钣金产品，必须克服和抑制成形过程中起皱的发生和发展。过去人们常常采用试错的方法，但这种单凭设计人员经验的方法在很大程度上带有盲目性，在面对复杂零件或新产品时往往无法解决。
图 3-4网格划分
3.4施加边界条件
凸模和凹模与板材的接触都采用面—面接触方式，并选择自动接触类型，摩擦系数取0.1。 有限元数值模拟时要尽量反映成型过程中的真实工艺情况，有了真实有效的工艺参数才能保证数值模拟结果的正确，按照实际加工时的参数来设定凸模的速度，求解时间为0.01秒，输出步数为100。
应力连续增长到超过弹性极限 ，其相应的应变为 ，然后下降到 ，其相应的应变为 。如在卸载的过程中，应力改变为 ，应变改变量为 ，则卸载后的应力与应变应等于卸载前的应力与应变分别减去它们的改变量，即
因为卸载是服从弹性规律的，所以上面式中的 和 是呈线性关系的，即 。式中，E为弹性模量。由此可以看出，在全部卸载后，即外载荷等于零的情况下，在物体内不仅会留下残余变形，而且有残余应力。
式中 为现时的后续屈服应力。增量理论中，塑性流动法则为
式中 是塑性应变增量； 为瞬时的非负比例系数，具体数值与材料硬化法则有关，与Mises屈服条件相关联的流动法则为 ，于是有
在金属成形过程中，回弹被定义为变形部件(板料)离开模具后的空间变化。板料冲压成形后的全应变包括弹性应变和塑性应变两部分，即
式中， 为弹性应变； 为塑性应变。在卸载过程中，只有弹性变形得到恢复，而塑性变形保持不变，并且卸载过程服从线性关系。
4 数值模拟结果与分析
计算完成后通过ANSYS后处理器POST1与POST26对结果进行后处理。首先先检验结果的有效性，一般来说，沙漏能量（hourglass energy）不超过内能（internal energy）的10%即被认为是可以接受的分析结果。本次分析的沙漏能为零，因此数值模拟分析的过程是可靠的。