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

0 、前言

　　板料冲压成形是机械工业中一种重要的加工方法，在航空、宇航、汽车等制造领域均有着广
泛的应用。长期以来，保证成形零件的成形质量、降低废品率，一直是板料成形研究的目标^［1，2］。
从力学角度而言，板料成形是一个同时涵盖几何非线性、材料非线性、边界非线性的复杂的力学过程。以往，分析成形问题多采用基于塑性理论的解析方法，由于成形过程的物理复杂性，不得不作出较多的简化和假设，这就使得解析方法只能进行很简单的成形问题的分析。近些年发展起来的板料成形的有限元模拟技术，使得对复杂成形问题的模拟分析成为可能，在优化成形工艺、提高成形质量、降低产品开发成本等方面将发挥潜在的重要作用。然而，目前的板料成形模拟技术仍存在诸多问题亟待解决，因此仍是国际上成形领域的研究热点^［3～5］。
　　本文以金属板料成形过程为对象，研究利用弹塑性有限变形增量有限元方法模拟成形过程的技术。采用基于板壳理论的8节点壳单元，以Kirchhoff应力张量和Green应变张量作为应力与应变的度量，建立了有限变形的Updated Lagrangian列式。在此基础上，提出了一种直接试错的接触算法，在每一增量步进行接触搜索，之后交替进行几何协调处理与接触力的协调处理，直至几何协调与接触力的协调同时满足。通过试验结果与计算结果的比较，说明采用的模拟算法是合理可行的。

1、8节点等参壳单元

　　由于板料实际上是一种板壳结构，因此采用根据板壳理论建立的壳单元进行板料成形模拟，显然是合适的。作者在板料成形模拟中采用了一种8节点等参壳单元，该单元满足板壳的两个假设，其一是中面法线的物质线元在板壳变形后仍为直线，其二是忽略与中面垂直的应力分量所产生的应变能，即认为垂直于中面的应力分量为零。
　　在该单元中定义了4种不同的坐标系：总体坐标系(直角坐标系)，节点坐标、总刚度矩阵、位移以及外载力矢均定义于该坐标系，坐标用x_i(i=1，2，3)表示，位移用u_i表示，e_i为3个坐标轴方向的单位矢量；随体的曲线坐标系，如图1所示，ξ和η是壳单元中面的两个曲线坐标，ζ是厚度方向的线性坐标；中面节点处的节点坐标系v_ik，其中k指明此节点为单元的第k个节点；最后一种坐标系为单元内某一点处的局部坐标系x′_i，x′₁取该点处ξ方向的切矢，由该点处ζ方向的切矢与η方向的切矢叉乘得到x′₃，x′₂由x′₃与x′₁叉乘得到，x′_i单位矢量化可得到局部坐标系的坐标基矢量e′_i。
　　壳单元内任意点的位置坐标，可以由节点坐标插值得到

　　　(1)

式中　m——表示中面
　　　n——壳单元的节点数
　　　δ_k——节点k处壳的厚度
　　　N_k(ξ，η)——插值函数

图1  8节点壳单元极坐标系

　　在总体坐标系中，壳单元的位移场由节点处法线的5个自由度来描述，即法线中点的3个位移u^m_ik和法线绕V_1k和V_2k的2个转角α_k和β_k，如图2所示。这样，可由单元节点处的位移插值得到单元的位移场

图2  壳单元的位移模式

(2)

　　为了引入厚度方向应力分量为零的假设，单元中的应变在局部坐标系x′_i中定义

ε=｛ε_x′，ε_y′，γ_x′y′，γ_y′z′，γ_x′z′｝^T　　　(3)

2、有限元求解列式

　　物体的变形是连续的，从初始时刻到t+Δt时刻之间任一时刻的物体构成都可以作为参考构形建立有限元列式。实际应用中常采用两种参考构形，一种以未变形时的初始构形为参考构形，称作Total Lagrangian描述，另一种是以t时刻的平衡构形为参考构形来描述t+Δt时刻的变形，称作Updated Lagrangian描述(或U_L方法)。考虑到成形的特点，采用U_L有限元列式。
2.1  有限变形的有限元列式
　　当质点在邻域内做刚性运动时，Kirchhoff应力张量的各分量在固定于空间的坐标系中保持不变。而Green应变率张量也是与刚性转动无关的量，并且Green应变ε_ij与Kirchhoff应力σ_ij在能量上是共轭的。因此，在本构方程中采用Kirchhoff应力和Green应变作为应力和应变的合理度量。
　　以t时刻的构形为参考构形，t+Δt时刻的虚功方程为

　　　(4)

式中　V，A_σ——t时刻构形的体积区和受载表面区
　　　_p,k——t+Δt时刻定义在t时刻构形上的面力
　　　_b,k——t+Δt时刻定义在t时刻构形上的体力
　　　δ——变分运算符
　　由于以t时刻的构形为参考构形，t时刻位移和应变均为已知，则式(4)可以写成

∫_VδΔε_ij^T(σ_ij+Δσ_ij)dV=　　
∫_AσδΔu_k^T_p,kdA+∫_VδΔu_k^T_b,kdV　　(5)

　　对上式进一步整理，并采用矩阵形式，可得到线性化的U.L.求解方程

(k₁+k_n)Δa=-R_s　　　(6)

式中　k₁——小位移刚度矩阵
　　　k_n——由于大位移而引出的矩阵
　　　Δa——节点位移增量
　　　——t+Δt时刻外载荷的等效节点力矢量
　　　R_s——应力场的等效节点力矢量
2.2  本构方程
　　采用Euler描述时，屈服函数写成一般形式为

f(σ_ij，E)=(σ_ij)-σ_s(E)=0　　(7)

　　采用Hill的正交各向异性屈服准则，并考虑到壳单元的假设，即σ₃=0，于是得到

²=a₁σ₁²+2a₁₂σ₁σ₂+a₂σ₂²+a₃τ₁₂²+a₄τ₂₃²+a₅τ₁₃²　　　(8)

　　6个各向异性参数a₁，a₁₂，a₂，a₃，a₄及a₅可通过试验确定。
　　在t时刻，根据塑性流动法则，可导出本构方程

σ_ij=C_ijkl,tkl　　(9)

　　上式写成矩阵形式，则可得到正交异性的弹塑性物理阵

　　　(10)

式中　——塑性流矢量
　　　——弹性物理阵
　　　H′——塑性模量(单向拉伸试验中应力塑性应变曲线的切线斜率)

3、成形模拟中的接触算法

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

图3  接触处理过程

3.1  接触搜索方法
　　由三角平面片组成的离散网格处理起来相对简单，并可用于描述复杂形状的模具，因此在接触算法中采用离散网格法描述模具的几何表面。接触处理时，认为模具是刚性的不发生变形，只须判别板料节点是否穿透模具表面，因此接触搜索称作单向接触搜索。模具表面称作主表面，板料节点称作从节点。接触状态在从节点处考察，因此要在每个从节点x^(s)处定义间隙函数

g(x^(s))=(x^(s)-x^(m))^.n(x^(m))　　　(11)

式中　x^(m)——x^(s)在主表面上最近的投影点
　　　n——主表面上x^(m)处的单位外法矢
　　接触搜索时，在每个从节点处检查函数g的值，若g＞0则该节点与主表面是分离的未发生接触，而g＜0则说明该节点穿透主表面进入了模具内部，应采取相应的几何协调处理。
3.2  接触中的几何协调
　　根据板料与模具之间的分离、附着和滑动这3种接触状态，将板料节点分为3类，即分离节点、附着节点和滑动节点。成形过程中，板料节点的状态不断变化，但必须满足接触条件：①板料节点不能穿透模具表面。②与模具接触的板料节点同模具之间不能作用拉力。

图4  间隙函数的定义

　　相应于板料节点的分离、附着和滑动3种状态，进行不同的几何协调处理：
　　(1)若从节点x^(s)在当前几何协调处理前处于分离状态，则将其拉回到它在主表面上的投影位置x^(m)，并认为该节点当前为滑动状态。
　　(2)若从节点x^(s)在当前几何协调处理前处于附着状态，则沿模具行程将其移动一个距离，并认为该节点当前仍为附着状态。
　　(3)若从节点x^(s)在当前几何协调处理前处于滑动状态，则将其拉回到它在主表面上的投影位置x^(m)，并认为该节点当前仍为滑动状态。
3.3  接触力协调处理
　　由于将从节点拉回到模具表面的位置是近似的，因此得到的接触力不一定是协调合理的。接触力(模具作用于板料的力)存在两种形式的不协调：①法向接触力为拉力，在这样的从节点处应将接触力释放掉，并令从节点成为分离节点。②切向接触力不协调，处于滑动状态的从节点所受切向接触力小于摩擦力，或处于附着状态的从节点所受切向接触力大于摩擦力。
　　切向接触力不协调又分为两种情况：
　　(1)t时刻从节点的状态为附着状态，如图5a所示。t+Δt时刻，若｜F′_c,t｜＞μ｜F′_c，n｜，则应在该从节点处施加调整力

ΔF_c=-β(│F′_c,t│-μ│F′_c，n│)t　　　(12)

式中　F′_c,t——t+Δt时刻的切向接触力矢
　　　F′_c,n——t+Δt时刻的法向接触力矢
　　　t——F′_c,t的方向矢量
　　　β——松弛因子(0＜β＜1)
　　同时，将该节点由附着状态变为滑动状态。
　　(2)t时刻从节点的状态为滑动状态，如图5b所示。t+Δt时刻，若F′_c^.Δu＞0则显然违背了能量守衡，从节点不能滑动而应将状态转变为附着状态，并对从节点施加调整力

ΔF_c=-βF′_c　　(13)

图5  从节点所受的接触力的变化

4、试验及模拟算例

4.1  方盒拉深零件的拉裂失稳分析
凹模尺寸为39mm×39mm，凸模尺寸为36.5mm×36.5mm，凸凹模圆角半径均为5mm。毛料为91mm×91mm的正方形板料，厚度为1.0mm，材料为LY12M。试验中，凸模行程13mm时零件靠近底部圆角处发生破裂，如图6a。通过成形模拟，得到了成形深度为13mm的零件构形(如图6b)及其应力应变分布(由于采用U.L.方法，每一增量步结束时都将应力应变转移到最新的构形上，因此最终得到的是真实的应力应变)。将模拟得到的零件内抽样点处的两个主应变ε₁和ε₂投放到成形极限图上，如图7所示，显然其中某些点已超出了成形极限，这说明该模拟算法能够用于预测拉裂失稳。

图6  拉深成形的方盒零件

图7  应变在FLD上的分布

4.2  球头形拉深件成形模拟
　　凹模直径为30mm，凸模球头半径为13.75mm，凹模圆角半径为5mm。采用56mm×56mm的正方形毛料，厚度为0.8mm，材料为LY12M，凸模行程为15mm。分别用本文的成形模拟程序和非线形有限元软件ABAQUS进行了成形模拟，对得到的成形力曲线进行了对比，图8为成形力的模拟结果，显然两条曲线的趋势是大体相同的。

图8  模拟得到的成形力曲线的对比

5  结论

　　在金属板料成形中，材料的变形量很大，此时再使用小应变假设就不能正确地反映变形的物理本质，必须采用有限变形理论。本文的直接试错接触算法，使有限元求解过程与接触处理过程相对独立，接触条件以边界条件的形式应用到有限元求解方程中。算法的稳定性好计算效率也较高，可适用于板料状零件成形时的接触处理。通过试验验证及与ABAQUS的对比，说明本算法是合理有效的。

塑性成形模拟及模具技术国家重点实验室开放课题资助项目(00-4)。

参 考 文 献

1，Zhou D,Wagoner R H.Development and application of sheet?forming simulation.J.Mater.Process.Technol.,1995,50?(1～4)∶?1～16
2，Taylor L,Cao J.Numerical simulations of sheet?metal forming.J.Mater.Process.Technol.,1995,50(1～4)∶168～179
3，曹简板料成形数值模拟研究进展.海内外青年科学家论坛，1998
4，Kanto Y,Yagawa G.A dynamic contact bulkling analysis by the penalty finite element method.Int J.Numer.Methods Eng.1990,29(4)∶755～774
5，Galbraith P C,Hallquist J O.Shell?element formulations in LS?DYNA3D:Their use in the modelling of sheet?metal forming.J.Mater.Process.Technol.,1995,50(1～4)∶158～167