1引言

　　板料成形是材料加工成形的一个重要分支，它广泛应用于汽车、航天、航空、家电等各个部门。随着现代工业的发展，板料成形件越来越复杂，人们对板料成形的质量和速度的要求也越来越高。传统的板料成形模具的设计依赖的是经验和直觉，并且通过反复试验调试来保证成形的质量。这不仅需要消耗大量的人力物力，而且周期长，效率低，不能适应社会发展的需要。上世纪七十年代以来，人们逐渐以数值模拟技术为辅助设计手段，大大降低了生产制造的成本。然而，由于板料成形是一种复杂的力学过程，其中包含几何非线性、材料非线性、接触非线性等强非线性问题，影响的参数非常多，这对数值模拟技术造成了极大的挑战。虽然目前板料成形的数值模拟软件已经商业化，但板料成形的模拟技术还不够完善，仍然是国内国外研究的热点。本文将主要介绍金属板料成形数值模拟的研究现状。

　　2板料成形有限元算法

　　用于板料成形的有限元算法大体可以分为弹-(粘)塑性和刚-(粘)塑性。粘塑性有限元法主要应用于热加工，而刚塑性有限元法在板料成形中应用有限；目前，弹塑性有限元法在板料成形数值模拟中应用较广。用弹塑性有限元法分析板料成形问题，不仅能计算工件的变形和应力、应变分布；而且还能计算工件的回弹和残余应力、残余应变等。由于板料成形过程中板料与模具具有相对滑动、粘着和脱落，所以必须控制增量步的大小从而尽量反映真实情况。

　　根据对时间积分方法的不同，板料成形有限元算法可以分为静力隐式、静力显式和动力显式。隐式算法是非条件稳定的，它在解决低速接触问题中更有优势[1]，而在解决复杂的三维模型时将会遇到困难[2]：当时间步长减小时，内存消耗会急剧增大，甚至造成收敛问题；局部的不稳定性很难达到力的平衡，这也不符合静态隐式的先决条件。显式算法克服了隐式算法的缺点[3]，然而它的不足之处在于，在解决像板料成形这样的条件稳定问题时，必须尽量消除惯性力的影响。对此一般有两个办法可采用：一个就是将运动能限制在应变能的5%以下；另外就是限制元素类型，一般只采用四节点的四边形或者八节点的实体型。由于回弹是一个准静态问题宜采用隐式算法。板料成形中常常先用显式算法模拟成形阶段，而用隐式算法模拟回弹。

　　Belytschco[4]等人提出并行混合时间算法，即将两种方法——显式方法和隐式方法综合起来，并行处理算法的一大难点就是进程间的交互。研究表明[5]并行混合时间算法可以用来进行板料成形的模拟，但是当计算大的工件时，时间步长增大了，计算精度也不够高；当时间步长变小时，并行程度就会下降。

　　Kwansoo Chung[6，7]等将形变理论与极值功路径相结合，提出了理想形变的思想，兰箭等人[8]编写了一步法有限元程序，并将其成功地应用于板料成形的数值模拟中。

　　3板料成形接触处理

　　板料成形过程是板料与模具不断的接触、作用与分离的过程，所以接触问题是板料成形数值模拟中关键技术之一。

　　要解决接触问题首先必须对模具进行描述，模具的一般表示方法有：解析函数法、参数曲面法、网格法等。解析函数法只能用来表示类似于圆筒件等简单的模型，而参数曲面法虽然能较准确的表达模具曲面，但也存在算法复杂效率低的问题。网格法在一定程度上克服了以上两者的缺点，加上有限元研究对象也用网格法表示，网格法成为大多数研究者采用的方法。陈中奎[9]提出了接触特征空间的概念，并建立了接触特征空间和模具曲面网格的对应关系，成功地解决了三维板料冲压成形过程有限元分析中的边界条件非线性问题。

　　对于几何形状复杂的模具往往采用曲面单元描述法[10]，也就是将参数曲面描述法和有限元网格描述法相结合，即首先利用CAD等造型软件生成模具的参数曲面，再通过IGES读入FEA系统，并沿参数方向将参数曲面划分为一系列曲面片，从而得到由曲面单元描述的模具。王晓林等[10]在此基础上采用逐级更新的搜索算法，推导出求迭代初值优化选择方案及满足无穿透原则的判断准则，并建立了板料成形的数值模拟接触搜索模型，这种模型不但改善了数值模拟的收敛性，而且提高了计算效率和稳定性。

　　目前求解接触问题常用的算法有拉格朗日乘子法和罚函数法，它们的共同点在于将有约束条件的变分问题转化为无约束条件的变分问题。拉格朗日乘子法模拟出的结果更为准确，但是计算效率较差，而且对于变形大的单元容易造成收敛困难。罚函数控制方程的阶数和带宽都小于拉格朗日乘子法，但是罚函数法的因子的取值对计算结果的精度影响很大。Shimizu和Sano[11，12]将罚函数发展为可以同时处理接触问题和摩擦问题，其主要思想就是将由于接触压力造成的罚因子和由于摩擦造成的罚因子加起来作为总的罚因子来进行搜索计算。

　　4摩擦处理

　　板料成形过程是板料相对于模具运动变形的过程，所以板料成形的结果必定受到板料与模具之间摩擦力的影响。研究表明摩擦力与以下一些参数有关：接触压力、滑动速度、钢板和模具的材料特性、表面粗糙度、润滑剂和同向变形。尤其是板料厚度/板料面积很小时，摩擦力对板料的流动和最后的应变分布有重要影响。除此之外[13]，研究还表明摩擦力的大小还与模具的几何特性有关。因此对摩擦力描述的正确与否直接影响到数值模拟的精确性。传统的用于有限元模拟的摩擦力模型有两种：一是Amonton准则，也就是库仑准则τ=μq，q是正应力；另一个是固定摩擦(constant friction)τ=mk，k是剪应力。一般情况下，正应力比剪应力大得多，所以第二种模型运用得比较多，但是它们只有在主应力值不是很大时模拟效果比较准确。1976年Wanheim和Bay[14]提出了一个通用的模型用来描述摩擦力：在正应力比较低(q/σ0<1.5，)时摩擦力正比于正应力，而当正应力比较高(q/σ0>3)时摩擦力为一常数，中间为过渡区间。后来陆续有人对此模型进行了大量的研究，绝大多数的研究集中在小变形的情况。S.B.Petersen[15]等人将其运用于大塑性变形的情况，综合运用模拟和实验的方法对三种模型进行了比较，研究结果表明Wanheim-Bay模型能够更加精确地反映实际情况。也有研究[16]表明基于局部接触条件的摩擦模型，特别是在板料高速滑动时，比传统的库仑模型更有优势。赵振铎等人[17]从微观角度出发，指出在变形过程中，真实接触面积只占名义接触面积很小的一部分，并对板料与模具相对滑动时的摩擦力学模型进行了初步的探讨。也有人[18]运用神经网络的原理对摩擦行为研究。有关数值模拟之摩擦模型的研究还有待进一步深入。

　　5成形极限图FLD(forming limit diagram)

　　在上世纪六十年代，Keeler首先提出成形极限图的思想，成形极限图的两个数据轴分别代表主应变和次应变，连接材料发生缩颈或断裂时对应的应变状态所得的点，就得到成形极限图，起初它只局限于次应变大于零的情况，后来由Goowin发展了次应变小于零的情况。虽然FLD的原理比较简单，但是其制作方法却很复杂，比较常用的是网格应变分析法，即将试样首先划分为一定的网格，根据样品失效后网格的形变求出应变值，从而作出成形极限图。这种成形极限图的缺陷在于：它依赖于材料的应变轨迹，对工业生产中材料的非线性应变轨迹，则不能用来预测成形。R.Arrieux[19]提出了基于应力而作的成形极限图即FLSD(Forming limit stress diagram)，它不受材料应变轨迹的影响。R.Arrieux分别将其应用于研究各向同性材料和各向异性材料，都得到不错的效果。它的优点在于，它可以用来预测复杂及多道次工序零件的缩颈的发生，可以分段表示出应变状况，这是经典的成形极限图所做不到的。Z.Zimniak[20]将FLD与MARC结合起来，数值模拟和物理实验证明基于摄动理论形成的FLSD能够比较精确的预测板料成形过程的开裂问题。

　　FLD可以由实验得方法测得，但是测量的周期比较长，而且无论哪种极限图都必须受到实验条件的影响，人们一方面寻找尽可能方便简洁的实验方法来制作成形极限图，另一方面也在积极寻求从理论上来对FLD进行预测，这就必须知道FLD的影响因素。实验研究结果表明，FLD受润滑剂、板料弯曲率，厚度以及落料方向的影响。材料性质也影响着FLD的位置，曲线在坐标轴上的截距和它的斜率取决于n和r值；不仅如此，预应变和板料厚度也能改变FLD图。D.W.A.Rees[21]对塑性镀锌钢的FLD进行保守的估算，其中考虑了板料落料方向、预应变、厚度以及n和r的值，研究表明：n值必须相对较高，预应变要尽量消除，落料方向为45度较好。为了判别材料成形状态，人们又先后发展了几种不同的屈服准则，其中具有代表意义的有von Mises模型、Hill(1948、1979、1993)模型、Gotoh模型以及Barlat-Lian模型等等[22]。

　　随着计算机技术的发展，人们寄希望于利用数值模拟的方法来制作FLD，Marciniak和Kuczynski在上世纪60年代第一次对FLD进行模拟，即M-K模型，实验证明由这种模型作出的FLD比较精确[23]。为了得到更好的效果Choi等人和Pishbin以及Gillis提出了另外一种数值模型称为JG方法，但是它的用户界面友好性不如M-K模型[23]。建立更好的数值模型仍然是FLD的研究方向之一。

　　用数值模拟的方法结合实验数据研究FLD的影响因素以及优化屈服准则也是FLD的热点之一。Antonio F.Ávila[23]等人在前人的基础上基于M-K模型来预测FLD的算法，其中用到了五种不同的屈服准则：von Mises、Hill(1948)、Hill(1979)、Logan和Hosford、Hill(1993)；FLD随着判据的不同而有显著的不同，在此与实验结果最为符合的是用Hill(1993)准则。

　　6缺陷的研究

　　板料成形的缺陷有起皱、开裂和回弹。

　　6.1起皱

　　在薄板冲压成形过程中当切向压应力达到或超过板料的临界应力时，就可能发生起皱。实际冲压成形时，一般常采用调节压边力和设置拉伸筋的方法控制拉深过程中冲压件突缘与侧壁起皱，但由此又会导致板材成形极限下降，使冲压件过早地发生破裂。特别是对于形状复杂或深度大的冲压件来说，既要保证成形形状或深度又要抑制起皱发生，在实际加工工艺中常常是比较困难的。当起皱与破裂两者发生矛盾时，起皱的抑制必须以板材不破裂为基本条件，所以防止起皱的难点在于抑制起皱成形规律的强有力工具。由于起皱经常发生在汽车的覆盖件上，形成波纹或折叠，极大地影响了表面质量，因此近年来对起皱的研究，特别是理论研究越来越多。

　　起皱的形成可以由Hill的分叉准则(Hutchinson将其特别用于薄壁薄壳的情况)来预测，可以假设一个局部的分歧位移区，在这个位移区内忽略接触条件和连续性条件，就可以推导出只决定于材料性能和几何性质的临界应力。以主应力为坐标轴的起皱极限图就是依据这种方法制得的。研究表明[24]应力状态的各向异性以及进入侧壁的板料的曲率对起皱的影响最大，各向异性指数增大时，板料发生起皱时拉伸距离也变大。起皱还与坯料的长宽比有关，也与冲头冲压的距离有关[25]。

　　用于起皱的有限元方法有两种：一种是具有完美结构的分叉分析，另外一种是具有初始非理想度的非分叉分析。对于第一种方法，它能够非常精确地预测起皱现象，但是分析过程复杂，而且需要一定的起皱判据；第二种方法的精确度受到初始非理想度的影响。Xu Weili[26]采用了一种简单的模拟起皱现象的方法：在接触算法中采用节点直接投影以及使用BHG(blankholder gap)，这种方法避免使用起皱判别准则，不需要复杂的塑性分叉理论以及坯料的初始非理想度。谢晖等人[27]从16个单元区域的起皱分析入手，用能量法计算，初步得到各单元的起皱临界因子，然后进一步进行区域搜索或选定区域计算，找到包含某单元最可能发生起皱的区域，得到能反映各单元起皱发生难易程度的起皱临界因子，并用云图将它描述出来，从而较精确、直观地预测起皱的发生。

　　对不均匀拉应力下的起皱机理、起皱规律、影响起皱高度的因素、消皱措施等各国都进行了大量研究工作。但对剪应力下的起皱的研究还很少。崔令江[28]研究表明：(1)剪应力起皱实验中试样的起皱区长度越大，其抗剪应力起皱能力就越差；(2)板材的硬化指数n值和厚向异性系数r值越大，抗剪应力起皱的能力越强。

　　6.2开裂

　　开裂是由于材料的强度或塑性不足，当拉应力超过临界值时便会发生。开裂发生的位置主要在：凸模端部、侧壁、凸模圆角部位、法兰部分和拉伸筋部分[29]。由于开裂的影响因素也很多，因此到目前为止还没有十分精确的判断开裂的准则。对于开裂的判断现在用的比较多的主要还是成形极限图[30]：包括FLCN(forming limit curves at neck)以及FLCF(forming limit curves at fracture)；最大拉深率也是主要判据之一，Min Wan等人[31]结合了圆柱和圆锥杯形的内在联系，从理论上给出了圆锥杯形的极限拉深系数的判断公式。

　　在有限元模拟上，人们一方面在探索新的模型，另一方面将现有的判断准则与有限元模拟结合起来进行研究。K.Komori[32]提出了一种节点分离的模型：将服从开裂准则的单元体和与它邻近的单元体共有的节点分离开来，也就是假设开裂沿着单元体的边界扩展；而一般的有限元模型中常常将服从开裂准则的单元直接删除，在一定程度上改变了材料本来的特性，导致模拟结果误差相对较大。Z.H.Chen等人[33]应用了混合有限元法(Mixed finite element method)，混合了位移和压力的有限元法，可以对板料成形中的不可压缩的塑性变形进行比较精确的模拟。而Ridha Hambli等人[34]用了一种反向技术的模型，可以识别开裂准则的临界值，从而将其应用于开裂模拟中，对研究开裂准则和开裂模拟都有很大的帮助。

　　另外，由于摩擦力在板料成形过程常常充当为拉应力的一部分，所以摩擦和润滑也是研究开裂问题的主要方向之一。

　　6.3回弹

　　回弹是板料成形后不可避免的现象，回弹现象主要表现为整体卸载回弹、切边回弹和局部卸载回弹，当回弹量超过允许容差后，就成为成形缺陷，影响零件的几何精度。因此，回弹一直是影响、制约模具和产品质量的重要因素。早期的工作主要基于解析法对一些简单零件纯弯曲或拉弯成形的回弹进行分析[35]。Gardiner最早基于弯曲工程理论模型对理想弹塑性板弯曲的回弹问题进行了研究[36]。Hill[37]首先建立了板弯曲的精确数学理论，并对刚塑性宽板纯弯曲问题进行了研究。现在人们逐渐寻找一些新的研究回弹问题的模型，这样能够研究更多种因素。Lumin Geng[38]采用了一种新的各向不等硬化模型，并且将这种硬化模型和四个屈服公式用在了ABQUS中：von Mises、Hill二次方程式、Barlat三参量和Barlat1996。

　　由于应用于板料成形的有限元方法的不断发展，人们逐渐对有限元模拟回弹进行了研究。由于回弹是成形的最后一步，成形过程模拟中产生的任何误差都会积累到回弹计算阶段；因此，回弹模拟结果的准确性很大程度上取决于成形过程的模拟精度。回弹问题的有限元模拟受到许多因素的影响：厚度方向积分点的数量、元素类型、网格大小、与凹模台肩处每个元素接触角度、惯性和接触算法的影响。不仅如此回弹还受到很多物理参量的影响：材料性质、硬化规则、摩擦系数、压边力(BHF)和卸载过程。K.P.Li[39]等人以及Luc Papeleux等人[40]利用有限元模拟的方法对影响板料回弹的物理参数和数值因素进行了详细的研究。Jenn-Terng Gau等人[41]在研究回弹问题时对包申格效应进行了不少研究。研究包申格效应普遍用单轴拉伸/压缩实验和循环扭曲实验。Jenn-Terng Gau等人研究发现两种相似的样品有相同的总的应变，却有完全不同的回弹效果。这是因为在应变区它们变形的历程不一样，包申格效应不仅与材料有关还与变形的历程有关。包申格效应对钢类材料影响不大对铝材影响较大。P.Xue等人[42]基于膜壳理论和能量方法研究发现弯曲半径越大回弹越大。

　　在有限元算法上普遍采用显式算法计算成形过程，再用隐式算法计算回弹。由于影响因素较多，过程复杂，今后回弹仍将是板料成形研究的主攻方向之一。

　　7展望

　　板料成形的数值模拟经过20多年的研究已经有了很大的发展。在实际生产中往往由于零件所要求的变形程度超过了材料的一次成形所允许的最大变形程度，必须采用多道次拉深成形的方法。如何处理多道次拉深变形的分配，以及如何精确确定材料首次拉深后的性能参数是板料成形方面亟待解决的问题。另外由于板料成形模拟技术在汽车覆盖件方面的应用，如何精确模拟复杂曲面的成形过程也是大家正在研究的热点之一。今后板料成形的数值模拟研究还要注意以下几方面：

　　1)进一步优化时间积分方案，能够将几种时间积分方案综合起来运用，各取所长。2)继续对摩擦行为进行研究，建立更合适的接触算法。3)优化各种本构关系，能够将成形极限图和有限元分析更好地结合起来。4)对缺陷问题还要深入研究，优化数值模拟模型。5)将有限元模拟系统和神经网络、自动控制等结合起来，形成大的分析系统。