有限元分析法

有限元分析法 麻省理工学院 材料科学与工程系

2001 年 2 月 28 日引言 有限元分析法（FEA）近年来已应用得非常广泛，现已成为年创收达数十亿美元的相关产业的基础。即使是很复杂的应力问题的数值解，现在用有限元分析的常规方法就能得到。此方法是如此的重要，以至于即便像这些只对材料力学作入门性论述的模块，也应该略述其主要特点。不管有限元法是如何的卓有成效，当你应用此法及类似的方法时，计算机解的缺点必须牢记在心头：这些解不一定能揭示诸如材料性能、几何特征等重要的变量是如何影响应力的。一旦输入数据有误，结果就会大相径庭，而分析者却难以觉察。所以理论建模最重要的作用可能是使设计者的直觉变得敏锐。有限元程序的用户应该为此目标部署设计策略，以尽可能多的封闭解和实验分析作为计算机仿真的补充。与现代微机上许多字处理和电子制表软件包相比，有限元的程序不那么复杂。然而，这些程序的复杂程度依然使大部分用户无法有效地编写自己所需的程序。可以买到一些预先编好的商用程序1，其价格范围宽，从微机到超级计算机都可兼容。但有特定需求的用户也不

必对程序的开发望而生畏，你会发现，从诸如齐凯维奇（Zienkiewicz2）等的教材中提供的

程序资源可作为有用的起点。大部分有限元软件是用Fortran语言编写的，但诸如felt等某些更新的程序用的是C语言或其它更时新的程序语言。 在实践中，有限元分析法通常由三个主要步骤组成：1、预处理：用户需建立物体待分析部分的模型，在此模型中，该部分的几何形状被分割成若干个离散的子区域——或称为“单元”。各单元在一些称为“结点”的离散点上相互连接。这些结点中有的有固定的位移，而其余的有给定的载荷。准备这样的模型可能极其耗费时间，所以商用程序之间的相互竞争就在于：如何用最友好的图形化界面的“预处理模块”，来帮助用户完成这项繁琐乏味的工作。有些预处理模块作为计算机化的画图和设计过程的组成部分，可在先前存在的CAD文件中覆盖网格，因而可以方便地完成有限元分析。2、分析：把预处理模块准备好的数据输入到有限元程序中，从而构成并求解用线性或非线性代数方程表示的系统

式中，u和f分别为各结点的位移和作用的外力。矩阵K的形式取决于求解问题的类

3、分析的早期，用户需仔细地研读程序运算后产生的大量数字，即型，本模块将概述桁架与线弹性体应力分析的方法。商用程序可能带有非常大的单元库，不同类型的单元适用于范围广泛的各类问题。有限元法的主要优点之一就是：许多不同类型的问题都可用相同的程序来处理，区别仅在于从单元库中指定适合于不同问题的单元类型。后处理：在有限元

1 C.A.Brebbia, ed.,有限元系统(Finite Element System), A Handbook, Springer-Verlag, Berlin, 1982.

2 O. C. Zienkiewicz and R.L. Taylor, 有限元法(The Finite Element Method), McGraw-Hill Co., London, 1989.

1列出的模型内各离散位置处的位移和应力。这种方法容易漏掉重要的趋向与热点，而最新的程序则利用图形显示来帮助用户直接观察运算结果。典型的后处理模块能显示遍布于模型上的彩色等应力线图，以表示不同的应力水平，显示的整个应力场的图像类似于光弹性法或云纹法的实验结果。 于特定的程序，其运行方法通常会在软件对的附件中详述，那些较为昂贵的程序的销售商经桁架的矩阵分析法 这为引入有限元分析的概念提供了很好的途径的矩阵分析法逐一考虑每根杆件的刚度，然后用这些刚度来确定由于节点（在有限元分常会提供实习场所或举办培训班，以帮助用户弄清楚程序运行中错综复杂的情况。但用户即使经过这种培训后，仍可能遇到问题：程序往往是“黑匣子”，其内部的运行机理是难以明了的。在本模块中，我们将略述大部分当今使用的有限元应力分析程序的基本原理，但目前的讨论仅限于线弹性范围内的分析。掌握这一理论有助于消除黑匣子综合症，同时也是对固体力学的基本分析方法作一总结。

在模块5中，比较充分地讨论过铰接的桁架，。对桁架精确地进行静力学分析，即使是复杂的桁架，其方程也能集成为矩阵形式，以便于用数值方法求解。此法有时称为“矩阵分析法”，它为早期有限元分析法的发展奠定了基础。 桁架析中通常称为结点）位移而产生的杆件力。注意到每根杆件对结点作用力的矢量和必须等于作用在结点处的外力，我们可以列出一系列线性代数方程。在这些方程中，结点位移是未知量而作用在结点上的力是已知量。这些方程可以方便地写成如下的矩阵形式（该方法亦由此得名）：

式中，和分别表示第个结点的位移和作用在第i

u

jfij

个节点上的力（实际上这些量都是

矢量，其分量沿各坐标轴）。系数数组称为整体刚度矩阵，其第分量的物理意义是第

ij

Kijj

个位移对第个力的影响。该矩阵方程可简写为 i

式中，用下标或黑体字表示矢量和矩阵。 位移这两者中有一个是已知的，但对一个给定的结单根桁架杆件的刚度矩阵 为了得出描述桁架系统的矩阵方程，第一步是需要建立单根桁架杆件每一端的力与位开始时，每个结点所受的外力或结点的点，同时指定任意的位移和任意的外力是不可能的。这些预先给出的结点力和结点位移是问题的边界条件。矩阵分析的任务就是要确定与给定的位移相对应的力、以及结点处有已知外力作用时的位移。

2移之间的关系。考虑如图1所示在yx−平面内的杆件，该杆件连接第i个和第j个结点，与水平线的夹角为θ。

图1 单根桁架杆件 沿杆件轴向及垂直于轴线的方向研究杆件待求的伸长量 δ，δ可用杆件两端点的位移

之差 来表示：

式中的水平及垂直分量（图1中画出的第个结点的两个位移分量都是 ，u和v分别为位移i负的）。上式可写成矩阵形式：

式中，c= cosθ、s= sinθ。

图2 结点力的分量 对线弹性体，与伸长量 δ对应的轴向力P可由胡克定律得出，即δ）（LAEP/=。水平

和垂直的结点力如图2所示。这些关系可写成总轴向力的形式：

进行矩阵乘法运算后得： 3 方括号中的量乘以后，称为“单元刚度矩阵”。其每一项均有物理意义，表示某个位移对某个力的贡献。对每根桁架杆件依次列出其单元刚度矩阵，将这些单元刚度矩阵组合后，就能得到整个系统的方程，因此单元刚度矩阵的计算是矩阵结构分析法的核心。桁架的矩阵分析法与一般的有限元法的主要区别在于如何形成单元刚度矩阵，而大部分其他的计算机运算是相同的。

LAE/ijK

多个杆件力的集成

图3 各杆件对总结点力的贡献 下一步是考虑各节点连接的多个桁架杆件的集成。与某节点（或结点）相连的各杆件都对该结点施加了作用力，此力由各杆件两端结点的位移来确定（见图3）。对给定的结点，

为保持其静力平衡，所有杆件力的合力必须等于作用在该结点处的外力：

elemifext

if

每个单元的刚度矩阵均加在总的或“整体”刚度矩阵的适当位置。整体刚度矩阵建

立了桁架所有位移与所有外力之间的关系。这个过程称为“集成”。上式中下标的数字必须是在整个桁架结构中的“整体”编号。但通常用局部的受力图计算单个单元的刚度矩阵较为方便，然后在集成单个矩阵时由计算机将其转换为整体编号。

elemijkijK

例1 集成过程是有限元法的核心，因此用手工求解一个简单的实例、以了解该方法的工作原理是值得的。考虑图4所示的两杆件桁架问题，其中结点从1到3的整体编号是任意编排的。因为每个结点通常能沿两个方向运动，故此问题中总的自由度为3╳2=6个。相应的

整体刚度矩阵是一个6╳6的方阵，它建立了六个位移与六个外力之间的关系。此问题中只有一个位移是未知的，因为除结点2的垂直位移（即第4个自由度）外，其余的位移均限定为零。图4所示为对每根杆件列出的整体编号和“局部”编号。

4 图4 两杆件桁架的整体与局部编号 应用局部编号，根据式（2），可算出两杆件中每一根杆件的4╳4单元刚度矩阵。倾角可

按结点坐标计算：

得出杆件1的刚度矩阵为：

杆件2的刚度矩阵为： （单位的统一十分重要。此处长度的单位为英寸、力的单位为磅、弹性模量的单位为磅/英寸（psi）。两杆的弹性模量均为2=E10Mpsi、横截面面积均为=0.1in。）这些矩阵行和列的编号均为从1到4，与杆件的局部自由度数相对应。而每个局部的自由度又可以与总系统的整体自由度中的一个相对应。通过对图4的观察，我们可列出下表，以表示局部与整体的自由度编号之间的关系：

A2

局部 整体 整体 杆件1 杆件2

例如，由此表可看出，杆件2的第2个自由度即整体编号系统中的第4个自由度，其第3个局部自由度对应于第5个整体自由度。因此，杆件2的刚度矩阵中第2行第3列元素的值

（记为），应该放在6╳6整体刚度矩阵中第4行第5列的位置。我们将此过程记为 )2(23k

两根杆件中每一杆件的刚度矩阵中都有16个元素，这些元素都要根据上表给出的对应关系放到整体刚度矩阵中。由此可得如下结果

5