基于三角形膜单元的钣金成形有限元逆算法

基于三角形膜单元的钣金成形有限元逆算法
本文以锥形盆为例，验证了根据有限元逆算法设计的毛料展开程序在毛料展开中的有效性和精确性，实例分析表明在工程允许精度范围内，有限元逆算法可对板料成形工艺方案进行快速评价，优化冲压工艺参数，是一种有效的钣金展开算法。

0 引言
钣金成形是一种生产效率高、成形质量好的金属加工方法，被广泛地应用于航空航天、汽车、轻工业等领域，在钣金成形加工中，合理的毛料形状和尺寸能够改善材料流动和应力分布状态，使钣金件厚度分布比较均匀，提高成形极限；还能降低突耳高度，减少后续工序的切边余量，提高材料的利用率。

目前毛料展开算法主要有经验法、势场模拟法、几何映射法、有限元增量法和有限元逆算法，钣金展开有限元逆算法的基本思想是：采用形变理论即全量理论，通过最小化塑性位能函数建立毛料形状与零件形状之间的关系，求出毛料形状和零件最终构形的厚度分布。

在零件的最终构形上采用三角形膜单元进行离散，三角形膜单元具有处理简单，计算量少，计算速度快的优点，适合逆算法快速高效的要求。

主要研究了基于三角形膜单元的有限元逆算法，并用商用软件PAM-STAMP模拟验证了逆算法的准确性和可行性。

1 有限元逆算法
在板料展开有限元逆算法中采用了以下假设：
①只考虑板料变形过程中的薄膜应变效应，不考虑板料的弯曲效应；
②材料不可压缩；
③变形过程是比例加载，基于形变理论即全量理论来分析变形；
④凸模、拉延筋、压边圈对板料的接触作用简化为法向作用力和切向摩擦力；
⑤采用对数应变来描述变形。

1.1 几何方程
逆算法只考虑冲压件的初始形态和最终形态，忽略中间变形过程，应变由两种构形间的几何差别定义。

如图1所示，当物质点从C o运动到C时，有：
点击图片查看大图
点击图片查看大图
图1 变形中的几何关系
式(1)中[F]-1为：
点击图片查看大图
式(2)中u，v，w是三角形膜单元的位移函数，[I]是单位矩阵。

[F]为线性变换梯度矩阵，描述了质点从初始构形映射到现时构形中的变形。

在钣金成形加工中，板料要经历较大的弹性、塑性变形，所以采用左Cauchy-Creen应变张量[B]来描述板料的相对变形。

[B]可由变形梯度矩阵求得：
[B]=[F][F]T，
板料上一点的变形可以用该点的三个方向的主伸长λi；来表示，各点的面内主伸长分量λ1，λ2可由该点[B]-1的特征值求得：
点击图片查看大图
式(3)中的θ是面内主伸长λ与最终构形的局部坐标系X轴之间的夹角，B ij为B-1的分量。

根据塑性变形中体积不变的假设，板料上各点的厚向主伸长λ3可由式(4)得，由式(5)可求得对数应变εx，εy，εxy。

点击图片查看大图
1.2 单元类型
在零件的最终构形上采用三角形膜单元进行离散，典型的三节点三角形单元如图2所示，节点编码为i，j，m以逆时针方向为正。

每个节点有2个位移分量u，v。

点击图片查看大图
图2 三节点三角形单元
三角形单元的位移函数取一次多项式，由三个节点的位移分量表示为：
点击图片查看大图
式中N i为形函数，A为三角形膜单元的面积。

1.3 屈服准则和本构方程
在平面应力状态下，采用Hill各项异性屈服准则和全量理论，得比例加载条件下的弹塑性板材的本构方程为：
点击图片查看大图
点击图片查看大图
式(7)中{σ}为Cauchy应力张量，{ε}为应变矩阵；σ和ε分别为等效应力和等效应变，E s为等效应力应变的割线模量，矩阵[P]可由式(8)求得。

其中R0，R45，R90分别为板料的各项异性参数：
点击图片查看大图
1.4 有限元控制方程及求解
根据Hill假设，变形体在平衡状态下塑性势能最小，由于板料成形可以看作是准静力问题，体力影响可以忽略。

塑性势能可以表示为：
点击图片查看大图
式(9)中右边第一项为塑性变形能的变化率，第二项为外力的功率，可近似的表示为变形能与等效外力功的差值，根据虚功原理和力平衡方程可算出等效外力功，所以近似的塑性势能表示为：
点击图片查看大图
式(10)中W p(x)为板料的塑性变形能，W f(x)为模具与板料间的摩擦功，W d(x)为压边圈与板料间的摩擦功，W b(x)为拉延筋与板料间的摩擦功，这些量都是毛料初始位置X的函数。

用共轭梯度法方法求解时，需要给出一个初始猜测值x o，如果猜测值比较合理，可以大大提高收敛的速度。

一般采用线性映射法来获得板料的初始猜测值，其基本思想是假设两种构形之间的变形是弹性变形，给定一个毛料形状的猜测值，由该猜测值跟最终构形之间的弹性变形来求得初始猜测值。

2 有限元逆算法的模拟验证
基于以上算法，研究开发了毛料展开设计程序。

钣金件的数值建模可在CAD造型软件中完成，用三角形单元将该模型离散后，导入到毛料展开程序中，就可以得到由上述算法计算得到的零件的毛料轮廓、应力、应变和厚度变化率等信息。

2.1 锥形盆的展开和厚度变薄率
根据不可展曲面近似展开方法的合理性判断标准，由有限元逆算法展开的平面板料通过PAM-STAMP加工仿真、逆向建模得到的零件的表面积与原零件的表面积比较，相对误差在10%以内，且重建零件的变薄率小于40%，最大应变值不超过成形极限时，就说明展开的平面图形的工艺性好。

图3(a)所示的是一个半径为50mm，深度为20mm，侧壁倾角为90°，各处倒角半径为10mm的一个锥形盆的CAD模型，锥形盆的厚度为1mm，表面积为15707.96mm2。

在I-Deas中采用三角形膜单元均匀离散CAD模型，网格大小为10mm，离散后的模型如图3(b)所示。

本例计算是在Windows XP上完成的，CPU是Intel P4/1.7GHz，内存512MB。

盆的材料为AlMg3，板料的各项异性指数为1.87，材料的硬化曲线如下式所示：
点击图片查看大图
经逆算法98s运算得到的锥形盆的毛料展开轮廓如图3(b)所示。

点击图片查看大图
图3 锥形盆的离散模型和毛料展开图
把展开后的网格文件，导入到PAMSTAMP中并提取网格节点坐标信息，利用提取的信息和CATIA造型软件的的逆向建模功能重建平面展开图形，测量展开后的表面积为24132.855mm2。

计算展开前和展开后表面积的相对误差为1.08%，小于10%，可以近似的认为展开前后两者的表面积是一致的。

图4是锥形盆的厚度变薄率分布云图，图5是对称平面上的厚向变形率分布曲线。

从图4中可以看出锥形盆边缘部分的材料变厚，存在起皱的危险，而盆底变薄，最大的变薄处发生在侧壁和盆底过渡底圆角部分，最大变形率为-0.3804E，该处最可能被拉裂，与实际情况相符。

点击图片查看大图
图4 锥形盆的应变分布图
点击图片查看大图
图5 锥形盆在对称平面上的厚向变形率分布曲线
2.2 模拟结果的验证
将有限元逆算法数值模拟结果与增量有限元计算结果进行分析，可以验证逆算法有限元数值模拟结果的可靠性和准确性。

为了验证有限元逆算法的结果准确与否，根据逆算法展开的毛料轮廓的边缘向外加一个偏移量为10mm的加工余量，然后导入到有限元增量模拟软件PamStamp-2C中，并划分网格，锥形盆凸模的几何及材料参数不变，由于锥形盆是轴对称零件，为了节省计算机资源，减少计算时间，只模拟锥形盆的1/4。

经多次模拟验证得知当冲压力为3.8kN，模型的成形性能最好，经过16min20s计算后，锥形盆的厚度变薄率如图6所示。

从图中可以看出锥形盆的最大变薄处在盆底过渡底圆角部分，最大变形率为-0.3838E，与有限元逆算法的模拟结果基本一致。

图7是压边力为3.8kN时，锥形盆的FLC曲线，从图7可以发现，所有单元的应变都在FLC曲线以下，根据PAMSTAMP的判别规则可知，所有的单元都在安全区域，都没有被拉裂，这说明冲压参数的选取是合理的，符合加工工艺的要求。

点击图片查看大图
图6 压边力为3.8kN时锥形盆的变薄率云图
点击图片查看大图
图7 压边力为3.8kN时的锥形盆的FLC曲线图
3 结论
根据最小化塑性位能函数结合板料成形的情况开发了基于三角形膜单元的有限元逆算法钣金展开程序，利用编写的程序对锥形盆进行了展开和模拟计算，得到了毛坯形状和厚度变薄率云图。

将得到的结果与增量有限元模拟软件PAMSTAMP进行了比较分析，得到如下结论：
(1)使用逆算法展开后的锥形盆的表面积和展开前的面积基本一致，表明逆算法具有较高的计算精度，展开的平面图形的工艺性较好。

(2)通过逆算法模拟还可以得到冲压工件最终构形上的厚度变化。

通过比较厚度变薄率，发现逆算法的模拟结果与PAMSTAMP的模拟结果基本一致，说明有限元逆算法可快速预测板料在最终构形上的变化情况，从而可以分析潜在的成形缺陷，预测板料的成形性能。

(3)有限元逆算法采用全量理论的本构关系，只考虑初始构形和最终构形，从而提高了计算速度，比使用PAMSTAMP速度更快，时间更短，节省了大量的计算时间，
总之，有限元逆算法是一种有效的钣金展开算法，能应用到钣金零件的早期设计阶段以及工艺参数的优化设计中，快速简单，在冲压工艺及模具设计中具有良好的发展前景。