1、收稿日期 :2003 - 10 - 29第 22 卷 第 3 期 计 算 机 仿 真 2004 年 5 月 文章编号 :1006 - 9348 (2005) 03 - 0088 - 04有限元数值模拟的若干关键技术朱伟 ,张质良(上海交通大学国家模具 CAD 工程研究中心 , 上海 200030)摘要 :该文针对当前在材料加工领域中进行有限元数值模拟各种加工样式下材料塑性成形过程时 ,会经常遇到的、几个棘手的关键技术处理难题 ,简要地介绍了当中关于如何选择适当的求解算法和迭代算法 ,如何处理表面接触问题 ,如何定义网络划分与网络重划分的原则问题以及其它一些应注意的问题的解决途径。关键词 :求解
2、算法 ;迭代算法 ;接触问题 ;网格划分 ;网格再划分中图分类号 :TP391. 9 文献标识码 :BSeveral Key Techniques About Using FEM SimulationZHU Wei ,ZHANG Zhi - liang(National Die Iterative arithmetic ; Contact problem ; Mesh generation ; Mesh regeneration1 引言目前 ,不管是在挤压、轧制等体积成形还是在板料成形过程中 ,以有限元为代表的数值模拟方法得到了较为广泛的应用。它不仅可以很方便地、且形象地描述金属在塑性成形过程
3、中的金属流动行为 ,还能提供工件变形体及模具在塑性成形过程中各种物理学场量的分布及其变化规律 ,是现阶段分析和解决日益复杂多样化的塑性成形加工工艺的经济型有力工具。但在这利用有限元进行数值模拟过程中 ,经常会遇到几个棘手的关键技术问题 ,如求解算法和迭代算法的选择 ,接触问题的处理以及网络如何划分和再划分的判据定义。这些问题如果处理不好 ,不光会直接影响到计算结果的终精度 ,而且还可能导致有限元数值模拟的正常进行。本文即对这几个关键技术问题的解决途径进行简要的阐述。2 求解算法为确保数值模拟计算的顺利进行 ,并得到正确的结果 ,第一步则是算法的收敛性及稳定性是必须保证的。目前来说 ,非线性有限
4、元计算从根本上说有两大类算法 :即所谓的显式法和隐式法。一般来说 ,静力或准静力问题多采用隐式算法求解。因为已证明在适当的参数选择下 ,它是无条件稳定。而对于动态接触问题 ,由于接触力与边界条件的不断变动 ,若再采用隐式算法求解 ,收敛性就难以保证。比如说 ,在求解板料成形问题时 ,多趋向于采用中心差分法列式的显式算法来求解板料成形问题。动态显式积分算法通常采用基于动力学的时间中心差分格式 ,有限元方程显式化 ,避免了准静态隐式积分算法中所需的迭代计算。由于中心差分法是条件稳定的 ,其时间步长不应超过临界步长。文献 1 中作者用商用模拟软件 PAMSTAMP 分析了板料的回弹现象 ,所采用的计
5、算方法是混合法 ,即以动态显式算法模拟其成形过程 ,以静态隐式算法来计算回弹 ,整个算法同时兼有动态显式算法和静态隐式算法的优点。计算结果证明 1 :由于回弹过程中卸载是占主导地位 ,因此变形是弹性占优的。这一特点便使隐式算法容易收敛 ,因此用静态隐式算法模拟回弹问题是可靠的。3 迭代算法第二步必须考虑的则是迭代算法的选择 。迭代算法选取的好坏直接影响到计算的精度和速度。如在模拟软件 De2form 中 ,就有 Newton - Raphson法、直接迭代 (Direct Iteration)法和 BFGS(Quasi - Newton)法等三种迭代算法可以选择。其中 ,Newton - Ra
6、phson法相比其它两种迭代方法来说 ,因其迭代次88 1995-2006 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved.数可更少 ,因而可适合用于大多数有限元问题的求解 ,但是该 方法的收敛性不如其它两种算法。直接迭代法则在求解刚塑性有限元的情况下 ,在其线性迭代阶段期间 ,向真实速度解的收敛速度比较快 ,因而开始时容易收敛 ,但是当迭代值很接近真实速度解时 ,收敛就会变得很慢 ,这样就使最终总迭代次数相对较多。 BFGS法比 Newton - Raphson法更容易收敛 ,常应用于弹塑性问题的求解 ,但是其精确度差
7、一些 ,一般情况下不采用这种迭代算法。其实一般来说 ,最好的数值计算效果是将前二种算法综合起来考虑 ,这样就可以取其各自的优点而避开其缺点 ,从而大大提高收敛速度和计算精度。具体综合方法是 :1)通常情况下直接使用 Newton - Raphson方法 ;2)当 Newton - Raphson 方法不能收敛时 ,改用直接迭代法求解。4 接触问题接触问题的处理一直是有限元法数值模拟金属塑性成型过程的一个主要难点 。它其实包括多个方面 ,如最重要的便是摩擦问题以及边界接触判定问题。摩擦问题一直是金属加工工作者研究的一个中心问题之一。摩擦条件的好坏直接与产品质量、设备效能、生产效率、工具磨损密切相
8、关。但在实际上 ,压力加工中的摩擦状态受很多因素的影响 ,尤其是在润滑条件下 ,摩擦状态更趋于复杂化。目前 ,在理论分析上处理金属塑性加工过程中的摩擦应用较为广泛的有以下两种方法 2 :1)滑动摩擦理论 :即库仑摩擦理论 ,认为接触面的摩擦力与接触面上的正压力成正比 ,即 f = N。式中 :f 为摩擦力 ,N 为正压力 , 为摩擦系数。2)剪切摩擦理论 :认为接触面间的摩擦力与金属的剪切屈服极限之比为常数 ,即 f/ k = m。式中 : k 为金属剪切屈服极限 ,m为小于或等于 1 的正数 ,称为剪切摩擦系数。一般对于弹性接触问题来说 ,通常采用库伦摩擦理论。但当正压力较大时 ,摩擦力可能
9、会大于剪切屈服极限 ,导致较大的误差。因此对于大多数塑性加工问题 ,通常可采用剪切摩擦理论。在应用有限元软件模拟加工时 ,为了简化模拟模型 ,一般简单地假定摩擦条件在整个加工过程中为一个不变量。但在实际当中 ,摩擦的状态在不断发生变化。因此 ,摩擦问题总是最终影响计算结果偏差较大的一个关键原因之一 2 。所以如何定义摩擦问题一直是众多科学研究者的热门课题。由于在模拟塑性变形过程中 ,工件的边界节点与模具表面的接触状态在不断地发生变化。这就需要在模拟计算中不断修改变形体构形的同时 ,同时也要对边界节点的接触状态进行动态自动识别 ,不断修改已产生变化的边界接触条件。这主要包括以下三个方面 :1)触
10、模的判断 ,即工件边界节点与模具表面接触的判断。这是由于变形体与模具之间存在着相对的运动 ,所以在每一加载步中 ,某些自由边界节点可能会与模具表面相接触。2)接触边界节点的位置修正。3)原接触边界节点脱离模具表面的判断 (即脱模的判断 ) 。如果某一接触节点有脱离模具的趋势 ,则它势必要沿模具型腔表面的法线方向向型腔内部运动。但是 ,在当前加载步内该节点是作为接触约束节点来处理的 ,其相对于模具型腔表面的法向速度应为零。因此 ,该节点处的约束反力 F在模具型腔表面的法线方向分量 Fn 必为负值 ,同时该边界节点处的法向应力 n 也相应地变为拉应力。所以当某一接触节点满足上述条件 ,即 n 0
11、或 Fn 0 则可判定该节点将脱离模具 ,相应地解除其边界约束条件 3 。5 网格划分有限元网格划分是有限元数值模拟前处理中的一部分重要工作 ,网格划分的质量和优劣对整个有限元计算过程及结果产生相当大的影响。不同的网格划分有着不同的计算时间复杂度和计算精度。如文献 4 中 ,作者通过研究一种网格划分边界与炉四壁重合的常规网格划分方法 ,以及另一种网格划分选取与射流方向近似、锅炉边界与网格边界不重合的方法对四周切向煤粉炉流场进行了数值模拟比较。结果明显表明 :采用后一种网格划分方法 ,可明显改善虚假扩散的影响 ,使得计算结果比较符合实验。又如在文献 5 中 ,作者针对二维矩形结构稳态传热为例 ,
12、分析了计算网格尺寸的不同对模拟稳态计算的结果及模拟计算时间的影响 ,进而提出了一种采用组合网格划分的计算新方法 (即先采用高迭代规定限额 ,利用大网格计算出稳态温度场 ,然后插值得到相应的小网格节点值 ,再取适当迭代规定限额 ,利用小网格重复上述运算 ,以便校正各插值点 ) 。最后结果表明 :此组合划分网格法与单一采用小网格尺寸模拟计算相比 ,计算时间可减少 80 % ,计算结果的准确程度可以提高 20 %。近年来随着计算机技术和各种数学算法理论的不断深入研究 ,各种智能式的划分网格算法层出不穷。如大型通用有限元分析软件 ANSYS就提供了智能网格划分功能。智能网格划分是一种比较高效的自由网格
13、划分方法。它考虑了几何图形的曲率以及线与线之间的接近程度 ,从而自动进行网格划分。它共分 10 个级别 ,级别越高网格划分就越粗 ,级别越低网格划分就越细。一般系统默认值为第 6 级 ,一般将级别设置在 4 8 之间。智能网格划分应用非常简便 ,用户不必考虑在划分区域内网格划分的具体细节 ,只须将划分设定为一定的级数即可 ,省掉了手工划分的麻烦。但采用智能网格划分同时也使用户失去了对自己感兴趣的区域网格划分的控制。文献 6 中 ,作者给出了采用 ANSYS 自适应网格划分的功能 ,对任一 CAD 模型进行了网格划分和数据重构 ,建立了可用于 C - MOLD 软件分析的 CAE模型 ,并将它与
14、由CAD 模型直接转换为 C - MOLD软件模拟的模型最终模拟结98 1995-2006 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved.果进行比较 ,验证了采用 ANSYS自适应网格划分的方法的可行性。又如在文献 7 中 ,作者对四边形网格自动生成方法进行了研究 ,并且对目前网格划分效果较好的铺路法进行了改进 ,提高了此算法的可靠性和稳定性。也有部分人 8 ,9 研究如何运用神经网络技术去划分网格 ,通过输入基于几何形状特征和外加载荷以及边界条件特征的各项参数 ,利用神经网络系统的自学习和自指导功能 ,输出网格单元的
15、大小或者反映单元稀疏的定性表示语言。文献 10 中 ,作者研究了基于特征技术的有限元网格划分问题 ,通过结合冲压产品的特征设计 ,定性描述了特征对有限元网格划分的意义 ,提出了基于特征变化过程的有限元网格划分方法 ,并讨论了如何实现特征对有限元网格划分的定量指导。6 网格重新划分在一般金属塑性成形问题中 ,塑性变形区有些部位的有效应变值可能达到或超过 2 ,尤其是在模拟金属体积成形时 ,坯料经过一系列的中间变形后 ,其形状简单的坯料会转变成形状复杂的锻件。在这之间 ,坯料的变形通常是很大的。因此 ,在用有限元模拟塑性成形过程中 ,网格重新划分将是不可避免的 ,不然受到畸变的网格不仅难以真实反应
16、塑变形区的真实情况 ,同时会使雅可比矩阵 J 为负 ,导致计算不能顺利进行。对于如何判断何时进行网格划分 ,通常可采用射线法。即在程序中设计一段程序 ,它可任取工件上某一质点 ,向任意方向作多条射线 ,如有其中一条与模具相交于两点 ,即认为质点进入了模具内部 ,此时就需进行网格再划分。如目前生成四边形网格的算法就有 11 :映射网格生成法、节点生成法、子域分离的网格自动生成方法、叉树法、扫描法以及栅格法等。其中节点生成法又可分为 Delaunay 三角形剖分网格自动生成方法和四边形网格自动生成 (三角形单元转换法和 Paving 单元法 (即所谓的铺路法 ) 、堆砌法 ) 。在文献 11 中
17、,作者对以上所有方法进行了横向比较 ,结果见表 111 。网格再划分一般可分为两个步骤 :1)先按照上述任一方法重新划分网格 ;2)将原来网格系统计算出的信息全部转换到新划分的网格系统上来。在这转换信息过程时 ,应先将旧网格系统中的有效应变等场变量插值到结点上 ,然后再将这些信息插值到新划分的网格系统上。表 1 几种典型网格生成算法对比分类方法 三角形 四边形 六面体 单元形状 密度控制 自动程度 时间复杂度映射法 也可 是 是 好 能 低 (N)Delaunay 法 是 否 否 好 好 高 (N)三角合成法 是 是 好 好 高 (N)Paving 法 也可 是 是 好 一般 一般 (N)子域
18、分离法 是 也可 否 一般 能局部控制 一般叉树法 是 是 是 内部好 能局部控制 高 (Nlog(N) )栅格法 也可 是 是 内部好 能局部控制 高 (Nlog(N) )7 其它问题另外还有其它一些因素 ,诸如粘性系数等。用粘塑性有限元法进行模拟研究时 ,必须要考虑材料的粘性系数。粘性系数的取值对最后的计算结果有一定的影响。但实际材料粘性系数的大小和多方面的因素有关 ,难以给出精确的解析表达式。文献 12 中 ,作者以塑性元件和粘性元件并联的模型为基础 ,研究了材料粘性系数的获取方法 ,给出了一些材料的粘性系数和影响因素之间的关系曲线 ,为模拟研究中粘性系数的获取提供了一种方便有效的方法。
19、当前如何把神经网络和人工智能学科等新兴学科有机地结合到传统的塑性加工当中去 ,使经典的数值模拟计算能够有效地减少误差和工作量 ,是当前的一个热点问题。如文献 13 、 14 即在这方面进行了尝试 ,实验结果均证明此种尝试是可行、切实有效的 13 ,14 。8 结束语本文针对当前有限元数值模拟中若干个对计算结果影响因素较大的关键技术问题处理进行了简要阐述 。相信随着成形塑性理论、有限元软件的不断完善以及计算机技术的飞速发展 ,以上关键技术问题都将逐步得到解决 ,使得数值模拟技术越来越完备 ,精度越来越高。参考文献 :1 张立力 ,齐恬 ,戴映荣 . 数值模拟参数和工艺参数对板材成形回弹影响的研究
20、 J .锻压技术 , 2002 ,27(6) : 18 - 21.2 李虎兴 . 压力加工过程的摩擦与润滑 M .北京 :冶金工业出版社 ,1993.3 单德彬 , 吕炎 , 王真 . 金属体积成形过程三维刚塑性有限元模拟技术的研究 J .塑性工程学报 ,1997 , 4 (3) : 98 - 102.4 汪箭 , 刘团结 ,范维澄 . 网格划分对四周切向煤粉炉流场计算的影响 J .计算物理 ,1997 , 14(4 ,5) : 595 - 596.5 刘晓燕 ,吴国忠 ,庞丽萍 ,高德生 . 网格划分对传热模拟计算的影响 J .大庆石油学院学报 ,1999 , 23(4) : 58 - 60
21、.6 黄英勇 ,辛勇 ,杨国泰 . 网格生成方法对注塑流动模拟的影响J .锻压技术 ,2002 ,27(2) :34 - 37.09 1995-2006 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved.7 张清萍 ,赵国群 ,王鹏程 . 金属成形有限元模拟四边形网格自动划 分方法 J .锻压技术 ,2003 , 27(2) : 32 - 35.8 D N Dyck , etal. Determining an Approximate Finite Element MeshDensity Using NeuralNetwo
22、rk Techniques J . IEEE Transactions onMagenetics ,1992 ,28 (2) :5.9 D A Low ther ,etal. A Density Driven Mesh Generator Guidedbya Neu2ralNetworkJ . IEEE Transactions on Magenetics , 1993 ,29(2) .10 虞春 ,周雄辉 ,张永清 . 基于特征的有限元网格划分研究 J .机械科学与技术 ,1999 , 18(2) : 210 - 212.11 黄志超 ,包忠诩 ,周天瑞 ,陈泽中 . 有限元网格划分技术研究
23、J .南昌大学学报 ,2001 , 23(4) : 25 - 31.12 周明智 ,洪深泽 ,许莹 .刚塑性有限数值模拟中粘性系数的获取方法 J .金属成形工艺 ,2002 , 20(6) : 25 - 26.13 李帆 ,翟福宝 ,张质良 . 使用神经网络对强力旋压有限元模拟进行参数预选 J .锻压技术 ,2002 , 27(1) : 32 - 35.14 李少平 ,郑静风 ,何丹农 . 利用神经网络及数值模拟获取变压边力控制曲线 J .金属成形工艺 ,2002 , 20(3) : 43 - 45.作者简介 朱 伟 (1976 - ) ,男 (汉族 ) ,江西萍乡人 ,博士研究生 ,研究方向
24、为材料成形过程的数值模拟仿真与工艺优化 ;张质良 (1935 - ) ,男 (汉 ) ,上海人 ,上海交通大学国家模具 CAD 工程研究中心博士生导师 ,教授 ,研究方向为材料加工过程的模拟和材料成形理论与新技术。(上接第 70 页 )图 4 仿真结果由仿真图 4b、 c 的结果可以看出 ,电机的转矩随负载的变化而变化 ;c 图的仿真结果显示 ,当供电频率达到 50Hz (即额定频率时 ) ,电机带负载的转速达到 2946rpm ,接近于额定的转速 2950rpm ,可见系统具有良好的动特性。4 结束语本文依据转差频率控制方式对电机进行数学建模并在Simulink 下进行了仿真 ,输入不同的电
25、源频率 f 就能实现对电机转速的改变 ,它具有以下特点 :1)建立在传递函数基础上的仿真模型 ,结构简单 ;2) Simulink 中的仿真系统 ,可以很容易的进行扩展 ,如加入泵的模型以及管道模型可以较容易地构建出一个完整的恒压供水系统仿真模型 ;3)根据自己的需要改变相应的系统参数 ,就可以得到不同情况下的仿真模型 ,通过观测器观察各项输出值 ,可以检验是否与真实情况相符。参考文献 :1 杜坤梅 , 李铁才 . 电机控制技术 M .哈尔滨工业大学出版社 ,2002. 2 张志涌 . 精通 MATLAB :5. 3 版 M .北京航空航天大学出版社 ,2000.3 贾建强 ,韩如成 ,左龙
26、. 基于 MATLAB/ SIMULINK的交流电机调速系统建模与仿真 J .电机与控制学报 , 2000 - 6.作者简介 程 伟 (1978 - ) ,男 (汉族 ) ,安徽淮北人 ,硕士研究生 ,研究方向 :计算机仿真 ,机电一体化技术 ;马树元 (1960 - ) ,男 (汉族 ) ,河北人 ,教授 ,博士生导师 ,研究方向 :测量技术及仪器 ,精密机械 ;吴平东 (1942 - ) ,男 (汉族 ) ,江苏无锡人 ,工学博士、教授、博士生导师 ,享受国务院有突出贡献津贴 ,研究方向 :机电一体化技术 ,远程控制与操作等 ;陈之龙 (1952 - ) ,男 (汉族 ) ,河北香河人 ,副研究员 ,研究方向 :机电一体化 ,远程操作与控制。19 1995-2006 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved.
