1有限元法简介 有限元法作为一个有效的数值分析工具,在许多科学领域当中有成功的应 用。本课程为机械类专业的高年级本科生开设,课程目标如下: 1)了解什么是有限元法、有限元方法的基本思想 2)了解有限元软件的发展水平,学会用有限元软件分析工程问题的方法 3)学习有限元法的原理,主要结合弹性力学问题来介绍有限元法的基本方法 包括单元分析、整体分析、载荷与约束处理、等参单元的概念等。 1.1有限单法的形成 在工程技术领域内,经常会遇到两类典型的问题。其中的第一类问题,可以 归结为有限个已知单元体的组合。例如,材料力学中的连续梁、建筑结构框架和 桁架结构。我们把这类问题,称为离散系统。如图1-1所示平面桁架结构,是由 6个只承受轴向力的“杆单元”组成,其中每根杆的受力状况相似。尽管离散系 统是可解的,但是求解图1-2所示的复杂离散系统,要依靠计算机技术。 P ua.u4Ho0Doo 图1-1平面桁架系统 图1-2大型编钟“中华和钟”的振动分 析及优化设计(曾攀教授) 第二类问题是针对连续介质,通常可以建立它们应遵循的基本方程,即微分 方程和相应的边界条件。例如弹性力学问题,热传导问题,电磁场问题等。由于 建立基本方程所研究的对象通常是无限小的单元,这类问题称为连续系统。连续 介质中的热传导问题可以归结为以下的控制方程、初始条件与换热边界条件:
1 有限元法简介 有限元法作为一个有效的数值分析工具, 在许多科学领域当中有成功的应 用。本课程为机械类专业的高年级本科生开设，课程目标如下： 1)了解什么是有限元法、有限元方法的基本思想。 2)了解有限元软件的发展水平，学会用有限元软件分析工程问题的方法。 3)学习有限元法的原理，主要结合弹性力学问题来介绍有限元法的基本方法， 包括单元分析、整体分析、载荷与约束处理、等参单元的概念等。 1.1 有限单法的形成 在工程技术领域内，经常会遇到两类典型的问题。其中的第一类问题，可以 归结为有限个已知单元体的组合。例如，材料力学中的连续梁、建筑结构框架和 桁架结构。我们把这类问题，称为离散系统。如图 1-1 所示平面桁架结构，是由 6 个只承受轴向力的“杆单元”组成，其中每根杆的受力状况相似。尽管离散系 统是可解的，但是求解图 1-2 所示的复杂离散系统，要依靠计算机技术。 图 1-1 平面桁架系统 图 1-2 大型编钟“中华和钟”的振动分 析及优化设计（曾攀教授） 第二类问题是针对连续介质，通常可以建立它们应遵循的基本方程，即微分 方程和相应的边界条件。例如弹性力学问题，热传导问题，电磁场问题等。由于 建立基本方程所研究的对象通常是无限小的单元，这类问题称为连续系统。连续 介质中的热传导问题可以归结为以下的控制方程、初始条件与换热边界条件： t T Q c z T y z T x y T x + = + + （1- 1）
初始温度场也可以是不均匀的,但各点温度值是已知的 7=0=7(xy:) (1-2) 通常的热边界有三种,第三类边界条件如下形式 (1-3) of 图1-3V6引擎的局部 尽管我们已经建立了连续系统的基本方程,由于边界条件的限制,通常只能 得到少数简单问题的精确解答。对于许多实际的工程问题,还无法给出精确的解 答,例如,图1-3所示V6引擎在工作中的温度分布。这为解决这个困难,工程 师们和数学家们提出了许多近似方法。 在寻找连续系统求解方法的过程中,工程师和数学家从两个不同的路线得到 了相同的结果,即有限元法( Finite element method)。有限元法的形成可以回顾到 二十世纪的50年代甚至更早些时间,基本思路来源于固体力学中矩阵结构法的 发展和工程师对结构相似性的直觉判断。对不同结构的杆系、不同的载荷,用矩 阵结构法求解都可以得到统一的矩阵公式]。从固体力学的角度来看,桁架结 构等标准离散系统与人为地分割成有限个分区后的连续系统在结构上存在相似 性,可以把矩阵结构法推广到非杆系结构的求解。 1956年 M.J. Turner, R W. Clough, H C. Martin, LJTopp在纽约举行的航空学会 年会上介绍了一种新的计算方法,将矩阵位移法推广到求解平面应力问题。他们 把结构划分成一个个三角形和矩形的“单元”,利用单元中近似位移函数,求得 单元节点力与节点位移关系的单元刚度矩阵 1954-1955年, J H. Argyris在航空工程杂志上发表了一组能量原理和结构分 析论文 1960年, Clough在他的名为“ The finite element in plane stress analysis”的论 文[]中首次提出了有限元( Finite element)这一术语 数学家们则发展了微分方程的近似解法,包括有限差分方法,变分原理和加 权余量法。在1963年前后,经过 J F Besseling, RJ.Melosh, REJones, RH Gallaher, THH Pian(卞学磺)等许多人的工作,认识到有限元法就是变分原理中Rz近
初始温度场也可以是不均匀的，但各点温度值是已知的： ( ) 0 0 T T x,y,z t= = （1- 2） 通常的热边界有三种，第三类边界条件如下形式： ( ) h T-Tf n T λ = − （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和 YK. Cheung(张佑启)发现能写成变分形式的所 有场问题,都可以用与固体力学有限元法的相同步骤求解。1969年B. A Szabo 和 G. CLee指出可以用加权余量法特别是 Galerkin法,导出标准的有限元过程来 求解非结构问题[] 我国的力学工作者为有限元方法的初期发展做出了许多贡献,其中比较著名 的有:陈伯屏(结构矩阵方法),钱令希(余能原理),钱伟长(广义变分原理) 胡海昌(广义变分原理),冯康(有限单元法理论)。遗憾的是,从1966年开始 的近十年期间,我国的研究工作受到阻碍。 12有限元法的基本思路 有限元法是求解数学物理问题的一种数值计算近似方法。它发源于固体力 学,以后迅速扩展到流体力学、传热学、电磁学、声学等其它物理领域山。有限 元法的基本思路可以归结为:将连续系统分割成有限个分区或单元,对每个单元 提出一个近似解,再将所有单元按标准方法组合成一个与原有系统近似的系统 下面用在自重作用下的等截面直杆来说明有限元法的思路 (一)等截面直杆在自重作用下的材料力学解答 L L计1 n+1 图1-4受自重作用的等截面直杆 图1-5离散后的直杆 受自重作用的等截面直杆如图1-4所示,杆的长度为L,截面积为A,弹性 模量为E,单位长度的重量为q,杆的内力为N。试求:杆的位移分布,杆的应 变和应力 N(x=qlL-x
似法的一种变形，发展了用各种不同变分原理导出的有限元计算公式。 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。试求：杆的位移分布，杆的应 变和应力。 N(x) = q(L − x)
du(x) N(x)dx q(L-x)dx l(x)= (1-4) EA q (二)等截面直杆在自重作用下的有限元法解答 1)离散化 如图1-5所示,将直杆划分成n个有限段,有限段之间通过一个铰接点连接 称两段之间的铰接点为结点,称每个有限段为单元。其中,第i个单元的长度为 L;,包含第i,计+1个结点 2)用单元节点位移表示单元内部位移 第i个单元中的位移用所包含的结点位移来表示, (x)=l1+ (x-x) L 其中u,为第i结点的位移,x,为第i结点的坐标。 第i个单元的应变为E,应力为σ,内力 为N, du l1+1-l Ea E =+1 (1-7) q(L+l+1) N=Ag EA(ui+ 图1-6集中单元重量 3)把外载荷集中到节点上 把第i单元和第计+1单元重量的一半9(L+L) ,集中到第计1结点上
EA q L x dx EA N x dx du x ( ) ( ) ( ) − = = = = − x x Lx EA q EA N x dx u x 0 2 ) 2 ( ( ) ( ) （1- 4） (L x) EA q dx du x = = − (L x) A q x = E x = − (二)等截面直杆在自重作用下的有限元法解答 1） 离散化 如图 1-5 所示，将直杆划分成 n 个有限段，有限段之间通过一个铰接点连接。 称两段之间的铰接点为结点，称每个有限段为单元。其中，第 i 个单元的长度为 Li，包含第 i，i+1 个结点。 2） 用单元节点位移表示单元内部位移 第 i 个单元中的位移用所包含的结点位移来表示， ( ) ( ) 1 i i i i i x x L u u u x u − − = + + （1- 5） 其中 i u 为第 i 结点的位移， i x 为第 i 结点的坐标。 第 i 个单元的应变为 i ，应力为 i ，内力 为 Ni ： i i i i L u u dx du − = = +1 （1- 6） i i i i i L E u u E ( ) 1 − = = + （1- 7） i i i i i L EA u u N A ( ) 1 − = = + （1- 8） 图 1-6 集中单元重量 3）把外载荷集中到节点上 把第 i 单元和第 i+1 单元重量的一半 2 ( ) q Li + Li+1 ，集中到第 i+1 结点上
4)建立结点的力平衡方程 对于第计1结点,由力的平衡方程可得: q(L1+L1) (19) 2 L,并将(1-8)代入得: 令M1=L14 1+(+,)u1-Ln2=9(1+)L2 (1-10) 2EA 根据约束条件,1=0。 对于第n+1个结点, -u+un+l 2EA (1-11) 建立所有结点的力平衡方程,可以得到由n+1个方程构成的方程组,可解出 n+1个未知的接点位移 例:将受自重作用的等截面直杆划分成3个等长的单元,试按有限元法的思路求 解。 定义单元的长度为a L L 对于结点1 对于结点2,由(1-10)式可得, L: l1+2l2 同样,对于结点3有, L3 2+2a23-u,=y EA 对于结点4,可以有两种处理方法 1)直接用第3个单元的内力与结点4上的载荷建立平衡方程 EA(4-l3)
4）建立结点的力平衡方程 对于第 i+1 结点，由力的平衡方程可得： 2 ( ) 1 1 + + + − = i i i i q L L N N （1- 9） 令 +1 = i i i L L ，并将（1- 8）代入得： 2 1 2 ) 1 (1 2 (1 ) i i i i i i i L EA q u u u − + + + − + = + （1-10） 根据约束条件， u1 = 0。 对于第 n+1 个结点， 2 n n qL N = EA qL u u n n n 2 2 − + +1 = （1-11） 建立所有结点的力平衡方程，可以得到由 n+1 个方程构成的方程组，可解出 n+1 个未知的接点位移。 例：将受自重作用的等截面直杆划分成 3 个等长的单元，试按有限元法的思路求 解。 定义单元的长度为 3 L a = 对于结点 1， u1 = 0 对于结点 2，由（1-10）式可得， EA qa u u u 2 − 1 + 2 2 − 3 = 同样，对于结点 3 有， EA qa u u u 2 − 2 + 2 3 − 4 = 对于结点 4，可以有两种处理方法。 1）直接用第 3 个单元的内力与结点 4 上的载荷建立平衡方程 2 3 qa N = ， a EA u u N ( ) 4 3 3 − =
13T4 2EA 2)假定存在一个虚拟结点5,与结点4构成了虚拟单元4 0, u =u 13 L 3 L 在结点4上应用(1-10), l3+(+3)u4-35-2E+2 q u, +u 整理后得到线性方程组, 2-101 EA 12 EA 5q 2EA 解得, iqa 2EA 13有限元法的计算步骤 有限元法的计算步骤归纳为以下三个基本步骤:网格划分,单元分析,整体 分析 13.1网格划分 有限元法的基础是用有限个单元体的集合来代替原有的连续体。因此首先要对弹 性体进行必要的简化,再将弹性体划分为有限个单元组成的离散体。单元之间通 过单元节点相连接。由单元、结点、结点连线构成的集合称为网格。通常把三维 实体划分成4面体或6面体单元的网格,如图1-9、1-10所示;平面问题划分成 三角形或四边形单元的网格,如图1-13、1-14所示
EA qa u u 2 2 − 3 + 4 = 2）假定存在一个虚拟结点 5，与结点 4 构成了虚拟单元 4 L4 = 0 ，u5 = u4， = → 4 3 3 L L 在结点 4 上应用（1-10）， 2 3 3 3 4 3 5 ) 1 (1 2 (1 ) a EA q u u u − + + − = + EA qa u u 2 2 − 3 + 4 = 整理后得到线性方程组， = − − − − EA qa EA qa EA qa u u u 2 1 1 1 2 1 2 1 0 2 2 2 4 3 2 解得， = = = EA qa u EA qa u EA qa u 2 9 4 2 5 2 4 2 3 2 2 1.3 有限元法的计算步骤 有限元法的计算步骤归纳为以下三个基本步骤：网格划分，单元分析，整体 分析。 1.3.1 网格划分 有限元法的基础是用有限个单元体的集合来代替原有的连续体。因此首先要对弹 性体进行必要的简化，再将弹性体划分为有限个单元组成的离散体。单元之间通 过单元节点相连接。由单元、结点、结点连线构成的集合称为网格。通常把三维 实体划分成 4 面体或 6 面体单元的网格，如图 1-9、1-10 所示；平面问题划分成 三角形或四边形单元的网格，如图 1-13、1-14 所示
tetrahed ron with 4 nodes brick with 8 nodes 图1-7四面体四节点单元 图1-8六面体8节点单元 图1-9三维实体的四面体单元划分图1-10三维实体的六面体单元划分 triangle with 3 nodes quadrilateral with 4 nodes 图1-11三角形3节点单元 图1-12四边形4节点单元 132单元分析 对于弹性力学问题,单元分析,就是建立各个单元的节点位移和节点力之间 的关系式。由于将单元的节点位移作为基本变量,进行单元分析首先要为单元内
图 1-7 四面体四节点单元 图 1-8 六面体 8 节点单元 图 1-9 三维实体的四面体单元划分 图 1-10 三维实体的六面体单元划分 图 1-11 三角形 3 节点单元 图 1-12 四边形 4 节点单元 1.3.2 单元分析 对于弹性力学问题，单元分析，就是建立各个单元的节点位移和节点力之间 的关系式。由于将单元的节点位移作为基本变量，进行单元分析首先要为单元内
部的位移确定一个近似表达式,然后计算单元的应变、应力,再建立单元中节点 力与节点位移的关系式。 以平面问题的三角形3结点单元为例。如图1-15所示,单元有三个结点I M,每个结点有两个位移u、v和两个结点力U、V。 大大个 △天X 图1-13平面问题的三角形单元划分图1-14平面问题的四边形单元划分 U 图1-15三角形3结点单元 单元的所有结点位移、结点力,可以表示为结点位移向量( vector):
部的位移确定一个近似表达式，然后计算单元的应变、应力，再建立单元中节点 力与节点位移的关系式。 以平面问题的三角形 3 结点单元为例。如图 1-15 所示，单元有三个结点 I、 J、M，每个结点有两个位移 u、v 和两个结点力 U、V。 图 1-13 平面问题的三角形单元划分 图 1-14 平面问题的四边形单元划分 图 1-15 三角形 3 结点单元 单元的所有结点位移、结点力，可以表示为结点位移向量（vector）：
uvu 结点位移{6} 结点力y=P 单元的结点位移和结点力之间的关系用张量( tensor)来表示, F}=[]{y (1-12) 133整体分析 对由各个单元组成的整体进行分析,建立节点外载荷与结点位移的关系,以 解出结点位移,这个过程为整体分析。再以弹性力学的平面问题为例,如图1-16 所示,在边界结点i上受到集中力P,P作用。结点i是三个单元的结合点,因 此要把这三个单元在同一结点上的结点力汇集在一起建立平衡方程。 图1-16整体分析 i结点的结点力 U(+U1(2+U0=∑U r(+(2)+V(3)=∑ i结点的平衡方程:
结点位移 = m m j j i i e v u v u v u 结点力 = m m j j i i e V U V U V U F 单元的结点位移和结点力之间的关系用张量（tensor）来表示， e e e F = K （1-12） 1.3.3 整体分析 对由各个单元组成的整体进行分析，建立节点外载荷与结点位移的关系，以 解出结点位移，这个过程为整体分析。再以弹性力学的平面问题为例，如图 1-16 所示，在边界结点 i 上受到集中力 i y i Px ,P 作用。结点 i 是三个单元的结合点，因 此要把这三个单元在同一结点上的结点力汇集在一起建立平衡方程。 图 1-16 整体分析 i 结点的结点力： + + = e e Ui Ui Ui Ui (1) (2) (3) ( ) + + = e e Vi Vi Vi Vi (1) (2) (3) ( ) i 结点的平衡方程：
∑Ue=Pr 14有限元法的进展与应用 有限元法不仅能应用于结构分析,还能解决归结为场问题的工程问题,从二 十世纪六十年代中期以来,有限元法得到了巨大的发展,为工程设计和优化提供 了有力的工具。 141算法与有限元软件 从二十世纪60年代中期以来,进行了大量的理论研究,不但拓展了有限元 法的应用领域,还开发了许多通用或专用的有限元分析软件。 理论研究的一个重要领域是计算方法的研究,主要有 大型线性方程组的解法 √非线性问题的解法, √动力问题计算方法 目前应用较多的通用有限元软件如表1-1所列 表1-1常用的有限元分析软件 简介 MSC/Nastran 著名结构分析程序,最初由NASA研制 MSC/Dytran 动力学分析程序 MSC/Marc 线性分析软件 ANSYS 通用结构分析软件 ADINA 非线性分析软件 ABAQUS 非线性分析软件 另外还有许多针对某类问题的专用有限元软件,例如金属成形分析软件 Deform、 Autoform,焊接与热处理分析软件 Sys Weld等。下面列出了一些有限元 软件开发商的网址: MSC中国, http://www.mscsoftware.comcn/ ANSYS中国http://www.ansys.comcn/ ABAQUS, http://www.abaqus.com ADINA http://www.adina.com DEFORM http://www.deformcomcn/ AUTOFORM http://www.autoform.del YSWELD http://www.esi-group.com/products/sysweld/
= = i y e e i e i x e i V P U P ( ) ( ) （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/
