改进型逆有限元法在某轿车翼子板快速成形模拟中的应用

改进型逆有限元法在某轿车翼子板 快速成形模拟中的应用那景新陆善彬 李 威 胡 平（吉林大学汽车工程学院，长春，130022） 摘 要 本文为进一步提高一步成形逆有限元法的模拟精度，将四节点等参膜单元应用于一步成形计算，并对空间四边形单元的四个节点可能不在同一平面引起的单元翘曲进行了修正，编写了计算机程序，对某轿车翼子板进行了实际模拟，通过一步模拟结果与增量法模拟结果和实验结果的对比，发现改进后的一步成形模拟结果与增量模拟结果已非常接近，而计算时间不到增量法的三十分之一，进而验证了本文方法的有效性。

关键词 一步成形 逆有限元法 等参膜单元 翘曲修正1 引言翼子板是轿车上比较典型外覆盖件件之一，质量要求高、成形难度较大。

工业发达国家普遍采用数值模拟技术，对其成形工艺进行分析与优化，以最大限度地缩短模具调试周期，降低制造成本。

目前国际上流行的成形模拟软件，大多采用基于动力学显式积分的增量有限元算法，与隐式有限元方法相比，计算效率已经有明显改善，但要完成翼子板类拉延件的成形模拟，一般也需要10～12小时的计算时间。

一步成形模拟方法是基于全量理论的一种有限元逆算法[1]，J.L. Batoz、 Y.Q. Guo、C.H. Lee 等人在这方面作了大量研究工作，并在坯料形状预示方面取得了很大成功[2-3]，但在成形性模拟方面，效果一直不够理想。

为进一步完善一步成形理论，作者自2001年开始，作了一系列的相关研究，先后提出基于截面线展开的初始场预示算法[4]、松弛因子优化算法[5]及处理摩擦边界条件的等效刚度法[6]等。

研究中作者发现，现有研究大多采用三角形膜单元[7,8]，而精度更高的四节点等参单元却很少使用，其主要原因是板成形过程非常复杂，很难保证四节点单元的四个节点始终在同一个平面上。

对四边形单元而言，如果四个节点不在同一个平面上，将会引起单元翘曲，此时的刚度矩阵难以琢磨。

为进一步提高一步成形模拟的计算精度，作者对翘曲四边形单元的刚度矩阵进行了修正，将四边形等参膜单元应用于一步成形模拟计算，并通过某汽车翼子板的一步成形模拟结果与增量法模拟结果及实验结果的比较，验证了本文方法的有效性。

2 一步成形基本理论一步模拟将板料成形过程简化为一个简单加载的变形过程，采用全量理论进行分析。

计算中可以只考虑初始状态和变形终了构形，而忽略中间状态和构形的变化[9]。

对板料终 了构形建立虚功方程如下： {}{}{}{}()0int =−=−=∑∫∫∑∑∗∗ev A TT eeext ee eedA f u dv W W W σε (1)其中eW int ，e ext W 分别表示单元内力虚功和外力虚功；{∗u }，{∗ε}分别表示虚位移和虚应变，{}f 是外力，在有限元数值计算过程中，（1）式往往不能严格满足，需要采用Newton-Raphson方法进行迭代求解[8]。

假定外力与内力的差值（残余力）为)}({u R ，则： 0)}({)}({)}({int ≠−=u F u F u R ext (2))}({u F ext 为节点外力向量，)}({int u F 为节点内力向量，该有限元方程的Newton-Raphson迭代格式为[10]：)}({})]{([iiiT u R u u K =∆ (3) }{}{}{1i i i u u u∆+=+ω (4)式中，ω为迭代收敛因子，)]([iT u K 为第i 步迭代的切线刚度矩阵，在全局坐标下可以写为：}{}{)}({)]([i u u T u u R u K =⎥⎦⎤⎢⎣⎡∂∂−= (5) 3 对四节点等参单元的改进由于板成形分析处理的是空间问题，为最大限度地减少单元翘曲对模拟结果的影响，必须对普通四节点等参单元模型进行修正[11],图1示出翘曲的四边形1234，在作修正性调整中，第一步是建立局部坐标轴xyz 以便四个节点交替地高出和低于xy 平面H 个单位。

单元在xy 面上的投影是ABCD 。

按照通常的方法，我们形成作用在八个自由度{D B A A v u v u ⋅⋅⋅}上的单元ABCD 的刚度阵]'[k 。

下一步需要变换矩阵][T ，以使]'[k 扩展为作用于{42111w u w v u ⋅⋅⋅}的12×12刚度矩阵。

为了求得][T ，我们考虑一个力的变换，这个变换满足三个不在一个平面的平衡方程，并对节点编号的次序是不敏感的。

在每一个标有字母的节点，平行于x 和平行于y 的节点力分解成与它们的边相平行的分量，如同图2中平行于AB 边的12f 和21f 一样，当12f 和21f 从标有字母的节点转换到标有数字的节点时，要求力矩平衡方程成立，因此可求得z 方向的“分离力”。

从图3可得：图1 局部坐标系xyz 中的翘曲四边形图2 单元在xy 面上的投影图3 与xy 面平行且垂直于AB 边的视图)(2112f f L Hf f ABzB zA +== (6) 把所有四个边的结果汇总起来，则有zD zC z zC zB z zB zA z zA zD z f f f f f f f f f f f f −=−=−=−=4321,,, (7)力的变换}'{][}{f T f T=是TyD yB xB yA xA TTz x z y x f f f f f T f f f f f }{][}{42111⋅⋅⋅=⋅⋅⋅ (8)}{f 为修正后的单元节点力向量, }'{f 为修正前的单元节点力向量, T T 的3，6，9和12行表示由方程(7)归纳的变换。

这些行的系数[8]是H 和加标字母节点的局部xy 坐标的函数。

T T 的其余各行各包含一个1和七个0，并表示xA x f f =1， ⋅⋅⋅=yA y f f 1等等条件。

于是翘曲单元的刚度阵是]]['[][][T k T k T=。

4 计算实例为了验证本文算法，采用了一个实际汽车覆盖件——某汽车翼子板（图4），通过逆向工程方法得到其几何模型，然后进行有限元离散化，最后运用本文算法进行了实际模拟。

图5显示了本文算法采用的有限元模型及拉延筋布置图，单元数为11218，节点数为11422；所用材料参数：毛坯厚度0.8mm ，各项异性系数62.1=γ，应力应变的关系：238.09.557εσ=MPa ，屈服函数使用Hill 各项异性屈服函数[12]，工艺参数：摩擦系数0.1，压边力600KN ，拉延筋水平阻力和举力如表1。

图6是本文算法预示的毛坯形状，图7是本文算法预示厚度分布，图8是利用我们自行开发的KMAS 显式算法软件，在在相同的材料、工艺条件下预示的厚度分布；增量模拟所用的板料即是本文算法预示的毛坯形状。

表2给出若干典型点两种算法所得厚度的比较结果，图9显示了截面'A A −(略去拉延筋处影响)的厚度分布曲线。

表1拉延筋工艺参数 (单位：KN/mm)R1R2R3R4R5R6R7图4 某汽车前翼子板照片图5 前翼子板计算模型及等效拉延筋布置图1R 2R 3R 4R 5R 6R 7R 'A A水平阻力 0.165 0.1650.1650.1650.08250.165 0.140 举力0.1300.1300.1300.1300.1300.1300.110表2 典型点处两种算法所得厚度与实验结果测得厚度的比较结果 1234 5 6 789实验测试 0.743 0.755 0.6130.7410.7070.7290.727 0.767 0.670增量法模拟 0.695 0.713 0.6180.7350.6880.6890.693 0.764 0.704一步成形模拟0.690 0.723 0.5840.7130.6480.6840.671 0.717 0.677从上面结果可以看出，改进后的一步成形模拟结果与增量法模拟结果已经非常接近，但模拟时间不到增量法的三十分之一。

不仅可以用于设计阶段做成形性快速评估，更适合在制造阶段进行工艺条件优化。

与增量法相比，一步模拟不仅效率高、数据准备也简化很多，很有可能经过进一步的努力，将其与CAD 系统集成，开发出可以自动进行成形性分析和成形工艺优化的CAD 知识系统软件产品。

5 结论一步成形模拟理论是建立在比例加载假定的基础上，对于大多数以拉延和胀形为主要变形特征的问题，由于其成形过程接近比例加载，因而，可以得到比较理想的模拟结果。

对于拉延过程中有明显的局部反成形，成形过程不符合比例加载假定的问题，会产生一定的计算误差。

C.H. Lee 和H. Huh 等人作过这方面的研究[13]，研究认为这种情况可以将整个过程分解为两个或多个一步成形问题来处理，这样可以使模拟精度有明显改善，当然计算时间会相应增加。

但整个计算过程与增量法相比，其计算效率仍有很大优势。

一步成形模拟的毛坯轮实际毛坯轮廓线图6 模拟预示坯料形状与实际坯料形状对比图7 一步成形方法预示厚度分布图8 增量法模拟预示厚度分布图9 截面'A A −厚度分布图567894123参考文献（略）。