采用试错接触算法的板料成形有限元模拟技术

1、采用试错接触算法的板料成形有限元模拟技术摘要：以金属板料成形过程为对象，采用基于板壳理论的8节点壳单元，以Kirchhoff应力张量和Green应变张量作为应力与应变的度量，采用有限变形的Updated Lagrangian列式，研究有限元方法模拟成形过程的技术。在此基础上，提出了一种直接试错的接触算法，将非线性的接触边界条件线性化处理，接触处理与有限元计算相对独立。通过试验结果与计算结果的比较，说明采用的模拟算法是合理可行的。    叙词：金属板料  成形  有限元  模拟  接触0 、前言板料冲压成形是机械工业

2、中一种重要的加工方法，在航空、宇航、汽车等制造领域均有着广泛的应用。长期以来，保证成形零件的成形质量、降低废品率，一直是板料成形研究的目标1，2。从力学角度而言，板料成形是一个同时涵盖几何非线性、材料非线性、边界非线性的复杂的力学过程。以往，分析成形问题多采用基于塑性理论的解析方法，由于成形过程的物理复杂性，不得不作出较多的简化和假设，这就使得解析方法只能进行很简单的成形问题的分析。近些年发展起来的板料成形的有限元模拟技术，使得对复杂成形问题的模拟分析成为可能，在优化成形工艺、提高成形质量、降低产品开发成本等方面将发挥潜在的重要作用。然而，目前的板料成形模拟技术仍存在诸多问题亟待解决，因此仍是

3、国际上成形领域的研究热点35。本文以金属板料成形过程为对象，研究利用弹塑性有限变形增量有限元方法模拟成形过程的技术。采用基于板壳理论的8节点壳单元，以Kirchhoff应力张量和Green应变张量作为应力与应变的度量，建立了有限变形的Updated Lagrangian列式。在此基础上，提出了一种直接试错的接触算法，在每一增量步进行接触搜索，之后交替进行几何协调处理与接触力的协调处理，直至几何协调与接触力的协调同时满足。通过试验结果与计算结果的比较，说明采用的模拟算法是合理可行的。1、8节点等参壳单元由于板料实际上是一种板壳结构，因此采用根据板壳理论建立的壳单元进行板料成形模拟，显然是合适的。

4、作者在板料成形模拟中采用了一种8节点等参壳单元，该单元满足板壳的两个假设，其一是中面法线的物质线元在板壳变形后仍为直线，其二是忽略与中面垂直的应力分量所产生的应变能，即认为垂直于中面的应力分量为零。在该单元中定义了4种不同的坐标系：总体坐标系(直角坐标系)，节点坐标、总刚度矩阵、位移以及外载力矢均定义于该坐标系，坐标用xi(i=1，2，3)表示，位移用ui表示，ei为3个坐标轴方向的单位矢量；随体的曲线坐标系，如图1所示，和是壳单元中面的两个曲线坐标，是厚度方向的线性坐标；中面节点处的节点坐标系vik，其中k指明此节点为单元的第k个节点；最后一种坐标系为单元内某一点处的局部坐标系xi，x1取该

5、点处方向的切矢，由该点处方向的切矢与方向的切矢叉乘得到x3，x2由x3与x1叉乘得到，xi单位矢量化可得到局部坐标系的坐标基矢量ei。壳单元内任意点的位置坐标，可以由节点坐标插值得到(1)式中m表示中面n壳单元的节点数k节点k处壳的厚度Nk(，)插值函数图1  8节点壳单元极坐标系在总体坐标系中，壳单元的位移场由节点处法线的5个自由度来描述，即法线中点的3个位移umik和法线绕V1k和V2k的2个转角k和k，如图2所示。这样，可由单元节点处的位移插值得到单元的位移场图2  壳单元的位移模式(2)为了引入厚度方向应力分量为零的假设，单元中的应变在局部坐标系xi中定义=x，y，

6、xy，yz，xzT(3)2、有限元求解列式物体的变形是连续的，从初始时刻到t+t时刻之间任一时刻的物体构成都可以作为参考构形建立有限元列式。实际应用中常采用两种参考构形，一种以未变形时的初始构形为参考构形，称作Total Lagrangian描述，另一种是以t时刻的平衡构形为参考构形来描述t+t时刻的变形，称作Updated Lagrangian描述(或UL方法)。考虑到成形的特点，采用UL有限元列式。2.1  有限变形的有限元列式当质点在邻域内做刚性运动时，Kirchhoff应力张量的各分量在固定于空间的坐标系中保持不变。而Green应变率张量也是与刚性转动无关的量，并且Green

7、应变ij与Kirchhoff应力ij在能量上是共轭的。因此，在本构方程中采用Kirchhoff应力和Green应变作为应力和应变的合理度量。以t时刻的构形为参考构形，t+t时刻的虚功方程为(4)式中V，At时刻构形的体积区和受载表面区p,kt+t时刻定义在t时刻构形上的面力b,kt+t时刻定义在t时刻构形上的体力变分运算符由于以t时刻的构形为参考构形，t时刻位移和应变均为已知，则式(4)可以写成VijT(ij+ij)dV=AukTp,kdA+VukTb,kdV(5)(k1+kn)a=-Rs(6)式中k1小位移刚度矩阵kn由于大位移而引出的矩阵a节点位移增量t+t时刻外载荷的等效节点力矢量Rs应

8、力场的等效节点力矢量2.2  本构方程采用Euler描述时，屈服函数写成一般形式为f(ij，E)=(ij)-s(E)=0(7)采用Hill的正交各向异性屈服准则，并考虑到壳单元的假设，即3=0，于是得到2=a112+2a1212+a222+a3122+a4232+a5132(8)6个各向异性参数a1，a12，a2，a3，a4及a5可通过试验确定。在t时刻，根据塑性流动法则，可导出本构方程ij=Cijkl,tkl(9)上式写成矩阵形式，则可得到正交异性的弹塑性物理阵(10)式中塑性流矢量弹性物理阵H塑性模量(单向拉伸试验中应力塑性应变曲线的切线斜率)3、成形模拟中的接触算法接触问题一直

9、是成形模拟中的难点和研究热点。经过十几年的发展，借助于有限元方法，产生了众多用于板料成形模拟的接触算法。本文中采用的是一种直接试错的接触算法，其主要特点是，有限元求解过程与接触处理过程相对独立，接触条件以边界条件的形式作用于有限元求解方程中。接触处理过程如图3所示，模具产生一个小的行程后，首先进行接触搜索，找到穿透模具表面的板料节点；然后，进行几何协调处理，根据这些节点的当前接触状态，以相应的方式将它们拉回到模具表面，之后，检查这些节点上所受的接触力是否协调，若不协调则进行接触力调整。这个过程重复进行，直至同时达到几何协调和接触力协调。以上的接触求解过程实际上包含两层迭代，外层是将非线性接触边

10、界条件线性化的迭代，内层则是根据给定的边界条件进行的弹塑性有限变形有限元求解迭代。图3  接触处理过程3.1  接触搜索方法由三角平面片组成的离散网格处理起来相对简单，并可用于描述复杂形状的模具，因此在接触算法中采用离散网格法描述模具的几何表面。接触处理时，认为模具是刚性的不发生变形，只须判别板料节点是否穿透模具表面，因此接触搜索称作单向接触搜索。模具表面称作主表面，板料节点称作从节点。接触状态在从节点处考察，因此要在每个从节点x(s)处定义间隙函数g(x(s)=(x(s)-x(m).n(x(m)(11)式中x(m)x(s)在主表面上最近的投影点n主表面上x(m)处的单位外

11、法矢接触搜索时，在每个从节点处检查函数g的值，若g0则该节点与主表面是分离的未发生接触，而g0则说明该节点穿透主表面进入了模具内部，应采取相应的几何协调处理。3.2  接触中的几何协调根据板料与模具之间的分离、附着和滑动这3种接触状态，将板料节点分为3类，即分离节点、附着节点和滑动节点。成形过程中，板料节点的状态不断变化，但必须满足接触条件：板料节点不能穿透模具表面。与模具接触的板料节点同模具之间不能作用拉力。图4  间隙函数的定义相应于板料节点的分离、附着和滑动3种状态，进行不同的几何协调处理：(1)若从节点x(s)在当前几何协调处理前处于分离状态，则将其拉回到它在主表面

12、上的投影位置x(m)，并认为该节点当前为滑动状态。(2)若从节点x(s)在当前几何协调处理前处于附着状态，则沿模具行程将其移动一个距离，并认为该节点当前仍为附着状态。(3)若从节点x(s)在当前几何协调处理前处于滑动状态，则将其拉回到它在主表面上的投影位置x(m)，并认为该节点当前仍为滑动状态。3.3  接触力协调处理由于将从节点拉回到模具表面的位置是近似的，因此得到的接触力不一定是协调合理的。接触力(模具作用于板料的力)存在两种形式的不协调：法向接触力为拉力，在这样的从节点处应将接触力释放掉，并令从节点成为分离节点。切向接触力不协调，处于滑动状态的从节点所受切向接触力小于摩擦力，或

13、处于附着状态的从节点所受切向接触力大于摩擦力。切向接触力不协调又分为两种情况：(1)t时刻从节点的状态为附着状态，如图5a所示。t+t时刻，若Fc,tFc，n，则应在该从节点处施加调整力Fc=-(Fc,t-Fc，n)t(12)式中Fc,tt+t时刻的切向接触力矢Fc,nt+t时刻的法向接触力矢tFc,t的方向矢量松弛因子(01)同时，将该节点由附着状态变为滑动状态。(2)t时刻从节点的状态为滑动状态，如图5b所示。t+t时刻，若Fc.u0则显然违背了能量守衡，从节点不能滑动而应将状态转变为附着状态，并对从节点施加调整力Fc=-Fc(13)图5  从节点所受的接触力的变化4、试验及模拟

14、算例4.1  方盒拉深零件的拉裂失稳分析凹模尺寸为39mm×39mm，凸模尺寸为36.5mm×36.5mm，凸凹模圆角半径均为5mm。毛料为91mm×91mm的正方形板料，厚度为1.0mm，材料为LY12M。试验中，凸模行程13mm时零件靠近底部圆角处发生破裂，如图6a。通过成形模拟，得到了成形深度为13mm1和2投放到成形极限图上，如图7所示，显然其中某些点已超出了成形极限，这说明该模拟算法能够用于预测拉裂失稳。图6  拉深成形的方盒零件图7  应变在FLD上的分布4.2  球头形拉深件成形模拟凹模直径为30mm，凸模球头半径为13.75mm，凹模圆角半径为5mm。采用56mm×56mm的正方形毛料，厚度为0.8mm，材料为LY12M，凸模行程为15mm。分别用本文的成形模拟程序和非线形有限元软件ABAQUS进行了成形模拟，对得到的成形力曲线进行了对比，图8为成形力的模拟结果，显然两条曲线的趋势是大体相同的。图8  模