1有限元法简介 有限元法作为一个有效的数值分析工具, 在许多科学领域当中有成功的应用。本课程为机械类专业的高年级本科生开设,课程目标如下: 1)了解什么是有限元法、有限元方法的基本思想。 2)了解有限元软件的发展水平,学会用有限元软件分析工程问题的方法。 3)学习有限元法的原理,主要结合弹性力学问题来介绍有限元法的基本方法,包括单元分析、整体分析、载荷与约束处理、等参单元的概念等。 1.1有限单法的形成 在工程技术领域内,经常会遇到两类典型的问题。其中的第一类问题,可以归结为有限个已知单元体的组合。例如,材料力学中的连续梁、建筑结构框架和桁架结构。我们把这类问题,称为离散系统。如图1-1所示平面桁架结构,是由6个只承受轴向力的“杆单元”组成,其中每根杆的受力状况相似。尽管离散系统是可解的,但是求解图1-2所示的复杂离散系统,要依靠计算机技术。    图1-1 平面桁架系统 图1-2 大型编钟“中华和钟”的振动分析及优化设计(曾攀教授)   第二类问题是针对连续介质,通常可以建立它们应遵循的基本方程,即微分方程和相应的边界条件。例如弹性力学问题,热传导问题,电磁场问题等。由于建立基本方程所研究的对象通常是无限小的单元,这类问题称为连续系统。连续介质中的热传导问题可以归结为以下的控制方程、初始条件与换热边界条件:  (1- 1) 初始温度场也可以是不均匀的,但各点温度值是已知的:  (1- 2) 通常的热边界有三种,第三类边界条件如下形式:  (1- 3)  图1-3 V6引擎的局部 尽管我们已经建立了连续系统的基本方程,由于边界条件的限制,通常只能得到少数简单问题的精确解答。对于许多实际的工程问题,还无法给出精确的解答,例如,图1-3所示V6引擎在工作中的温度分布。这为解决这个困难,工程师们和数学家们提出了许多近似方法。 在寻找连续系统求解方法的过程中,工程师和数学家从两个不同的路线得到了相同的结果,即有限元法(Finite Element Method)。有限元法的形成可以回顾到二十世纪的50年代甚至更早些时间,基本思路来源于固体力学中矩阵结构法的发展和工程师对结构相似性的直觉判断。对不同结构的杆系、不同的载荷,用矩阵结构法求解都可以得到统一的矩阵公式[1]。从固体力学的角度来看,桁架结构等标准离散系统与人为地分割成有限个分区后的连续系统在结构上存在相似性,可以把矩阵结构法推广到非杆系结构的求解。 1956年M..J.Turner, R.W.Clough, H.C.Martin, L.J.Topp在纽约举行的航空学会年会上介绍了一种新的计算方法,将矩阵位移法推广到求解平面应力问题。他们把结构划分成一个个三角形和矩形的“单元”,利用单元中近似位移函数,求得单元节点力与节点位移关系的单元刚度矩阵。 1954-1955年,J.H.Argyris在航空工程杂志上发表了一组能量原理和结构分析论文。 1960年,Clough在他的名为“The finite element in plane stress analysis”的论文[1]中首次提出了有限元(Finite Element)这一术语。 数学家们则发展了微分方程的近似解法,包括有限差分方法,变分原理和加权余量法。在1963年前后,经过J.F.Besseling, R.J.Melosh, R.E.Jones, R.H.Gallaher, T.H.H.Pian(卞学磺)等许多人的工作,认识到有限元法就是变分原理中Ritz近似法的一种变形,发展了用各种不同变分原理导出的有限元计算公式。 1965年O.C.Zienkiewicz和Y.K.Cheung(张佑启)发现能写成变分形式的所有场问题,都可以用与固体力学有限元法的相同步骤求解。1969年B.A.Szabo和G.C.Lee指出可以用加权余量法特别是Galerkin法,导出标准的有限元过程来求解非结构问题[1]。 我国的力学工作者为有限元方法的初期发展做出了许多贡献,其中比较著名的有:陈伯屏(结构矩阵方法),钱令希(余能原理),钱伟长(广义变分原理),胡海昌(广义变分原理),冯康(有限单元法理论)。遗憾的是,从1966年开始的近十年期间,我国的研究工作受到阻碍。 1.2 有限元法的基本思路 有限元法是求解数学物理问题的一种数值计算近似方法。它发源于固体力学,以后迅速扩展到流体力学、传热学、电磁学、声学等其它物理领域[1]。有限元法的基本思路可以归结为:将连续系统分割成有限个分区或单元,对每个单元提出一个近似解,再将所有单元按标准方法组合成一个与原有系统近似的系统。 下面用在自重作用下的等截面直杆来说明有限元法的思路。 等截面直杆在自重作用下的材料力学解答    图1-4 受自重作用的等截面直杆 图1-5 离散后的直杆   受自重作用的等截面直杆如图1-4所示,杆的长度为L,截面积为A,弹性模量为E,单位长度的重量为q,杆的内力为N。试求:杆的位移分布,杆的应变和应力。    (1- 4)   (二)等截面直杆在自重作用下的有限元法解答 离散化 如图1-5所示,将直杆划分成n个有限段,有限段之间通过一个铰接点连接。称两段之间的铰接点为结点,称每个有限段为单元。其中,第i个单元的长度为Li,包含第i,i+1个结点。 用单元节点位移表示单元内部位移 第i个单元中的位移用所包含的结点位移来表示,  (1- 5) 其中为第i结点的位移,为第i结点的坐标。 第i个单元的应变为,应力为,内力为:  (1- 6)  (1- 7)  (1- 8) 图1-6 集中单元重量 3)把外载荷集中到节点上 把第i单元和第i+1单元重量的一半,集中到第i+1结点上。 4)建立结点的力平衡方程 对于第i+1结点,由力的平衡方程可得:  (1- 9) 令,并将(1- 8)代入得:  (1-10) 根据约束条件,。 对于第n+1个结点,   (1-11) 建立所有结点的力平衡方程,可以得到由n+1个方程构成的方程组,可解出n+1个未知的接点位移。 例:将受自重作用的等截面直杆划分成3个等长的单元,试按有限元法的思路求解。 定义单元的长度为 对于结点1, 对于结点2,由(1-10)式可得,  同样,对于结点3有,  对于结点4,可以有两种处理方法。 直接用第3个单元的内力与结点4上的载荷建立平衡方程 ,  假定存在一个虚拟结点5,与结点4构成了虚拟单元4  ,,  在结点4上应用(1-10),   整理后得到线性方程组,  解得,  1.3 有限元法的计算步骤 有限元法的计算步骤归纳为以下三个基本步骤:网格划分,单元分析,整体分析。 1.3.1网格划分 有限元法的基础是用有限个单元体的集合来代替原有的连续体。因此首先要对弹性体进行必要的简化,再将弹性体划分为有限个单元组成的离散体。单元之间通过单元节点相连接。由单元、结点、结点连线构成的集合称为网格。通常把三维实体划分成4面体或6面体单元的网格,如图1-9、1-10所示;平面问题划分成三角形或四边形单元的网格,如图1-13、1-14所示。    图1-7四面体四节点单元 图1-8 六面体8节点单元      图1-9 三维实体的四面体单元划分 图1-10 三维实体的六面体单元划分      图1-11 三角形3节点单元 图1-12 四边形4节点单元   1.3.2单元分析 对于弹性力学问题,单元分析,就是建立各个单元的节点位移和节点力之间的关系式。由于将单元的节点位移作为基本变量,进行单元分析首先要为单元内部的位移确定一个近似表达式,然后计算单元的应变、应力,再建立单元中节点力与节点位移的关系式。 以平面问题的三角形3结点单元为例。如图1-15所示,单元有三个结点I、J、M,每个结点有两个位移u、v和两个结点力U、V。    图1-13 平面问题的三角形单元划分 图1-14 平面问题的四边形单元划分    图1-15 三角形3结点单元 单元的所有结点位移、结点力,可以表示为结点位移向量(vector): 结点位移 结点力 单元的结点位移和结点力之间的关系用张量(tensor)来表示,  (1-12) 1.3.3整体分析 对由各个单元组成的整体进行分析,建立节点外载荷与结点位移的关系,以解出结点位移,这个过程为整体分析。再以弹性力学的平面问题为例,如图1-16所示,在边界结点i上受到集中力作用。结点i是三个单元的结合点,因此要把这三个单元在同一结点上的结点力汇集在一起建立平衡方程。  图1-16 整体分析 i结点的结点力:   i结点的平衡方程:  (1-13) 1.4有限元法的进展与应用 有限元法不仅能应用于结构分析,还能解决归结为场问题的工程问题,从二十世纪六十年代中期以来,有限元法得到了巨大的发展,为工程设计和优化提供了有力的工具。 1.4.1算法与有限元软件 从二十世纪60年代中期以来,进行了大量的理论研究,不但拓展了有限元法的应用领域,还开发了许多通用或专用的有限元分析软件。 理论研究的一个重要领域是计算方法的研究,主要有: 大型线性方程组的解法, 非线性问题的解法, 动力问题计算方法。 目前应用较多的通用有限元软件如表1-1所列: 表1-1 常用的有限元分析软件 软件名称 简介  MSC/Nastran 著名结构分析程序,最初由NASA研制  MSC/Dytran 动力学分析程序  MSC/Marc 非线性分析软件  ANSYS 通用结构分析软件  ADINA 非线性分析软件  ABAQUS 非线性分析软件   另外还有许多针对某类问题的专用有限元软件,例如金属成形分析软件Deform、Autoform,焊接与热处理分析软件SysWeld等。下面列出了一些有限元软件开发商的网址: MSC中国, http://www.mscsoftware.com.cn/; ANSYS中国, http://www.ansys.com.cn/; ABAQUS, http://www.abaqus.com/; ADINA, http://www.adina.com/; DEFORM, http://www.deform.com.cn/; AUTOFORM, http://www.autoform.de/; SYSWELD, http://www.esi-group.com/products/sysweld/。 1.4.2应用实例 有限元法已经成功地应用在以下一些领域: 固体力学,包括强度、稳定性、震动和瞬态问题的分析; 传热学; 电磁场; 流体力学。 下面介绍一些有限元法应用的实例。 1)转向机构支架的强度分析(刘道勇,东风汽车工程研究院,用MSC/Nastran完成)  图1-17 转向机构支架的强度分析 2)金属成形过程的分析(用Deform软件完成) 分析金属成形过程中的各种缺陷。如图1-18所示,型材在挤压成形的初期,容易产生形状扭曲。  图1-18 型材挤压成形的分析。  图1-19 螺旋齿轮成形过程的分析  (a)模拟结果 (b)实物 图1-20 T形锻件的成形分析 3)焊接残余应力分析(用Sysweld完成)  图1-21 结构与焊缝布置  图1-22 焊接过程的温度分布与轴向残余应力 4)热处理过程的分析 BMW曲轴的感应淬火 (Induction quenching of crankshafts at BMW,用SysWeld软件完成) 在曲轴表面获得压应力,可以提高曲轴的疲劳寿命。  图1-23曲轴的有限元模型  图1-24有限元模型的局部  图1-25沿网格线52的残余应力分布, 红线为预测的轴向应力与径向应力之差,黑点为实测值。 复杂形状工件的组织转变预测(石伟,用NSHT3D完成) 预测工件的组织分布和机械性能  图1-26 四分之一工件的有限元模型  图1-27 淬火3.06 min 时的温度分布  图1-28 淬火3.06 min 时的马氏体分布 参考书目 王勖成,邵敏编著. 有限单元法基本原理和数值方法 . 北京 : 清华大学出版社, 1997 Saeed Moaveni. Finite Element Method Theory and Application with ANSYS. New Jersey: Prentice-Hall, Inc., 1999 郭和德编. 有限单元法概论(校内教材),1998 Zienkiewicz, R.L. Taylor. The finite element method( 5th ed). Oxford ; Boston : Butterworth-Heinemann, 2000 朱伯芳著. 有限单元法原理与应用(第2版). 北京: 中国水利水电出版社, 1998 思考题: 1、什么是有限元法? 2、简述有限元法的基本思路。 参考文献 [1] 蒋维城编著. 固体力学有限元分析. 北京:北京理工大学出版社,1989 [2] R.W.Clough. The finite element in plane stress analysis. Proc. 2nd ASCE Conf. On Electronic Computation. Pittsburg, Pa., Sept. 1960 [3] S. I. Oh, W. T. Wu, K. Arimoto. Recent developments in process simulation for bulk forming processes. Journal of Materials Processing Technology 111 (2001) 2-9