当板料成形数值模拟的有限元模型及应用

板料成形数值模拟的有限元模型及应用

摘要：本文回顾了板料成形数值模拟的二维截面和三维实体单元模型研究发展概况，阐述和比较了用于板料成形数值模拟几个典型有限元模型的性能和特点，指出了当前国际上该领域的研究热点。

关键词：有限元模型；自锁；沙漏；隐式算法；显式算法

1 引言

板料成形数值模拟技术在经历上世纪七、八十年代的起步探索阶段后，自八十年代末至今进入快速发展阶段，在发达国家的某些生产领域如汽车制造业已进入实用阶段，如今正朝着CAD/CAE/CAM一体化方向发展，形成所谓的“虚拟生产系统”。每三年一度的板料成形数值模拟国际会议“NUMISHEET”可以很好的说明板料成形数值模拟技术的发展水平。从NUMISHEET’2002发表的文献来看，板料成形的成形缺陷预测如起皱、拉裂已比较成熟，但回弹预测离实际应用尚有一段距离。特别是有限元模型及接触力计算是影响回弹预测精度的两个关键因素[1】，因此对合理的有限元模型及接触算法的研究仍然是当今世界提高数值模拟精度的前沿课题。

大变形弹塑性有限元法分析板料成形问题，不仅能计算工件的变形和应力、应变分布，而且还能计算工件的回弹和残余应力、残余应变分布。实践证明，大变形弹塑性有限元法能够有效模拟板料成形，是板料成形模拟的主要方法。在有限变形理论中，板料发生的大转动和大变形的描述常采用欧拉描述或拉格朗日描述。拉格朗日描述又有“完全的拉格朗日描述”(T．L格式)和“更新的拉格朗日描述”(U．L格式)之分。T.L格式以初始构形为参考构形。由于它是基于一个固定的构形来描述运动，从而使运动学变量的处理和积分都比较简单。但是，因为参考构形建立在初始构形上，根据T.L格式所得到的应力只是数学上的一种定义，这对荷载和边界条件随变形而变化的问题造成了应用上的困难。U.L格式采用前一个增量步结束时的构形来定义当前增量步的各变量值，格在不断刷新，参考构形也在不断变化，适用于格产生畸变时的情况。在板料成形数值模拟的实际应用中，经常采用近似于欧拉描述的有限元模型，采用率形式的本构关系，欧拉应力可以在增量步计算出来，便于实际应用。Belytschko等提出了随动坐标系的概念，在大转动、小应变的条件下，可以认为随动坐标系和材料的物质坐标系重合，本构关系可进一步简化为Green-Naghdi或Jaumann率形式的本构关系，并可在随动坐标系上施加材料的各向异性假设，大大简化了计算。由于Mindlin壳体理论比Kirchhoff壳体理论更适合于板料成形的实际情况，基于 Mindlin 壳体理论的退化壳元已成为应用的主流。针对板料成形回弹分析尚不成熟的实际情况，最近几年国际计算力学界普遍关注两种单元模型：一是采用修正的 Mindin壳体理论（考虑厚度方向的应力应变）固壳元[]和七参数壳元[]；一是采用EAS方法或B-bar方法构造的三维实体单元[]。它们都应用于板料成形的模具小园角处，更能合理地描述板料在此处的实际应力应变状态。下面的论文里论述了用于板料成形数值模拟的二维截面和三维实体单元模型，论文Ⅱ部分将详细地论述用于板料成形的壳单元模型。

2 板料成形二维截面单元模型

板料成形的某些特殊况可以简化为在其截面上为平面应力或平面应变状态。在这种情况下，可以用二维单元模拟其成形过程或进行回弹分析。Flanagan 和Belytschko[13]在八十年 代初提出了用于显式算法的一点积分带正交沙漏控制(γ稳定向量正交于线性速度场)的二维有限元模型，需要人工输入摄动沙漏控制参数保证其计算稳定性。剪切自锁问题被巧妙地避开，但这种有限元模型难于处理体积自锁问题。Belytschko等【14】随后又指出这种正交沙漏控制技术可由正确估价线性速度场的梯度从而保证有限元方程的一致性来获得，可通过分片试验因而保证收敛性。基于Hu-Washizu混合变分原理讨论了沙漏控制参数的性质及选择方法，但这篇文章只涉及了体积自锁问题。1986年Belytschko等基于Hu-Washizu广义变分原理提出了一个四节点四边形单元【15】，从位移场导出的应变场被投射到一个小空间以消除体积自锁和剪切自锁，沙漏模式被物理稳定因而避免了人工输入参数。以后Belytschko 等又发表论文【16】进一步阐述了采用应变投射方法消除体积自锁及剪切自锁的原理（假设应变法即B-bar方法），并根据简明的Hu-Washizu广义变分原理建立了一个四节点四边形单元模型。对于在非线性领域的应用，节点参数采用速度，基于随动坐标系的有限变形理论，建立了一个适用于板料成形的二维单元模型，并成功地消除了二维板料成形中的体积自锁及剪切自锁问题。

u等提出了沙漏控制的另外一种途径【17】(URI方法)，对二维单元应变及应力场在单元中心被展开成一阶泰勒级数，γ稳定向量可由应变场对自然坐标求偏导数获得。根据URI方法，沙漏控制无需任何人工参数，单元刚度阵及内力向量可被显式表达，但他们没有提及如何消除剪切自锁问题。以后u等基于应变场的一阶泰勒展开式，采用B-bar方法及选择多点积分技术（SI），在随动坐标系下建立了一个适用于板料成形的二维有限元模型【18】，并成功地消除了体积及剪切自锁行为。u等采用假设应变法，但B-bar阵由B矩阵的解析表达式修正而来，建立了可用于板料成形二维有限元模型，并模拟了杯型件的拉深问题【19】。

3 板料成形三维实体单元模型

由于板料的厚度较薄，三维实体单元用于板料成形过程模拟格划分不宜太疏，否则易造成刚度矩阵的奇异性，从而导致求解失败。但格划分太密无疑会大大增加计算工作量，降低了计算效率，所以在板料成形的一般区域宜采用壳单元。在板料成形的特殊区域（如模具的小园角处）应力应变状态特别复杂，而采用三维实体单元能较好地描述板料在该区域的实际应力应变状态，具有独特的优势,特别是用于板料成形的回弹分析。所以自上世纪八十年代至今适用于弹塑性大变形的三维实体单元仍然是研究的热点。

Belytschko等在上世纪八十年代初采用正交沙漏控制技术建立了适用于显式算法的三维实体单元模型【13，14】，但它们需要输入人工摄动沙漏控制参数，且难于处理体积自锁问题。以后又在九十年代采用B-bar方法提出了一个八节点三维实体单元模型【9】，不需要输入人工参数，消除了体积及剪切自锁行为，可用于隐式算法及显式算法。u等采用URI方法，应变及应力场在单元中心被展开成二阶泰勒级数，提出了一个三维实体单元模型【17】，但遗憾的是此单元模型没有通过分片试验，不能保证收敛性。以后u等基于应变场的二阶泰勒展开式，采用B-bar方法及选择多点积分技术（SI），在随动坐标系下建立了一个适用于板料成形的三维有限元模型【18】，并成功地消除了体积及剪切自锁行为。由于该单元模型将常量应变改为拟应变，可使单元通过分片试验，从而保证了收敛性。九十年代末，u等对该模型进一步研究，在每一个选择积分点建立随动坐标系从而保证了单元的坐标不变性，采用切向预测径向返回算法对弹塑性本构关系进行积分，可较好地计算板料成形的弹塑性大变形问题【10】。u等基于B矩阵的解析表达式，采用假设应变方法，分别按URI和SI积分技术建立了一个三维实体单元模型，并模拟了园板的方杯拉深问题【11】。等基于B矩阵的解析表达式，采用B-bar方法提出了一个三维实体单元模型，并模拟了圆形杯及方形杯的拉深问题[12]。在该文献里，等提出了一个自动确定的系数以消除剪切自锁，可适用于任意厚度的板料。

图1 园杯的拉深成形

如图1所示为一圆形杯的拉深成形数值模拟实例【12】，为验证等提出的三维实体单元性能，采用了他们所提出的单元模型。园板的几何及载荷数据如图所示，冲头直径D=80mm,模具园角R=8mm,材料数据如下：

杨氏模量：E＝206Pa 泊松比： ν＝0.3

摩擦系数：μ＝0.15 初始d 试样上应力最大处的直径屈服应力：σS=167.14Kpa

采用分段线性强化弹塑性模型和Von Mises屈服准则，由于只考虑板料的拉深成形问题，沿厚向设置一层单元，并根据对称性只对板料的1/4进行计算。图2给出了冲程为60mm时的板料的变形格，图3给出了板料的厚向应变分布。

图2 冲头行程为60mm时板料的变形格

图3 冲头行程为60mm时板料的厚向应变分布

美国Lawrence Livermore国家实验室的so最近采用大变形分析的非协调模式（EAS: The Enhanced Assumed Strain）构造了适用于板壳结构分析的三维实体单元模型[8]。在这篇文献里，so分别采用假设应变法及EAS方法消除了单元的剪切自锁及体积自锁行为，提出了增广应变从参数空间向物理空间转换时矩阵的简化计算方法，对各向同性材料可简化为对角形式，获得了可观的计算精度及计算效率。

EAS方法用于非线性有限元分析也存在诸多问题。一个普遍存在的问题是独立假设的增广应变参数在单元水平上消去。对板料成形中的弹塑性大变形问题，增广应变参数的消去需要在每一单元、每一迭代步进行，大大提高了计算成本，降低了计算效率，这已成为EAS方法用于非线性有限元分析的发展瓶颈。具体表现在以下两个方面：

（1）本构关系的限制。EAS方法用于非线性有限元分析，比较适合于各向同性弹性或拟弹性的本构关系。Puso在其论文里提出了的一种简化计算方法，对于各向同性材料，可简化成对角形式，大大实验机指针零点调剂对检测结果的影响实验机在使用之前降低了计算成本。对于弹塑性或各向异性本构关系，计算工作量会显著增大。

（2）求解方法的限制。EAS方法用于大变形分析，比较适合于用T.L.格式求解。无论是U.L.格式还是采用近似的欧拉描述方法都要求在每一个时间步计算一次，这显然会大大增加了计算工作量。而为了解决工程中的实际问题，用T.L.格式计算得到的Kirchhoff应力及Green应变还需要转换成真实应力及真实应变，无疑也会增大计算工作量，这也是T.L.格式用于大变形分析的一个缺陷。

参考文献

[1] Wang C T. An industrial outlook for springback predictability,measurement realibility, and compensation technology. In: NUMISHEET2002,2002:

[2] Parisch H. A continuum-based shell theory for non-linear applications. Int J Numer Meth Engng 1995,38:

[3] Hauptmann R,Schweizerhof K. A systematic development of ‘solid-shell’ element formulations for linear and non-linear analyses employing only displacement degrees of J Numer Meth Engng,1998,42:

[4] Miehe C. A theoretical and computational model for isotropic elastoplastic stress analysis in shells at large strains. Comput Meth Appl Mech Engng,1998,155:

[5] Büchter N,et al. Three dimensional extension of non-linear shell formulation based on the enhanced assumed strain concept. Int J Numer Meth Engng,1994,37:

[6] Sansour C. A theory and finite element fo从1上任就提出民族复兴rmulation of shells at finite deformation involving thickness change. Arch Appl Mech,1995,65:

[7] Bischoff M,Ramm E. Shear deformable shell elements for large strains and J Numer Meth Engng, 1997,40:

[8] Puso M A. A highiy efficient enhanced assumed strain physically stabilized hexahedral element. Int J Numer Meth Engng,2000,49:1029－1064

[9] Belytschko T,et al. Assumed strain stabilization of the eight node hexahedral element. Comput Meth Appl Mech Engng,1993,105:

[10] Liu W K,et al. A multiple-quadrature eight-node hexahedral finite element for large deformation elastoplastic analysis. Comput Meth Appl Mech Engng,1998,154:

[11] Zhu Y Y,et al. Unified and mixed formulation of the 8-node hexahedral elements by assumed strain method. Comput Meth Appl Mech En11月11日gng,1996,129:

[12] Li K P,et al. An 8-node brick element with mixed formulation for large deformation analyses. Comput Meth Appl Mech Engng,1997,141:

[13] Flanagan D P,Belytschko T. A uniform strain hexahedron and quadrilateral with orthogonal hourglass control. Int J Numer Meth Engng,1981,17:679－706

[14] Belytschko T,et al. Hourglass control in linear and nonlinear problems. Comput Meth Appl Mech Engng,1984,43:

[15] Belytschko T,et al. Efficient implementation of quadrilaterals with high coarse-mesh accuracy. Comput Meth Appl Mech Engng,1986,54:

[16] Belytschko T,et al. Assumed strain stabilization of the 4-node quadrilateral with 1-point quadature for nonlinear problems. Comput Meth Appl Mech Engng,1991,88:

[17] Liu W K,et al. Finite element stabilization matrices-a unification approach. Comput Meth Appl Mech Engng,1985,53:

[18] Liu W K,et al. Multiple quadrature underintegrated finite elements. Int J Numer Meth Engng,1994,37:3263－3289

[19] Zhu Y Y,et al. Unified and mixed formulation of the 4-node quadrilateral elements by assumed strain method: application to thermomechanical problems. Int J Numer Meth Engng,1995,38:685－716

王金彦，男，1964年生，博士生，副教授，现从事板料成形数值模拟理论及应用研究工作 联系：（021），e-mail: gerald_wang@(end)