1有限元法简介1.1有限单法的形成在工程技术领域内，经常会遇到两类典型的问题。
其中的第一类问题，可以归结为有限个已知单元体的组合。
例如，材料力学中的连续梁、建筑结构框架和桁架结构。
我们把这类问题，称为离散系统。
如图1-1所示平面桁架结构，是由6个承受轴向力的“杆单元”组成。
尽管离散系统是可解的，但是求解图1-2所示这类复杂的离散系统，要依靠计算机技术。
图1-1 平面桁架系统图1-2 大型编钟“中华和钟”的振动分析及优化设计（曾攀教授）第二类问题，通常可以建立它们应遵循的基本方程，即微分方程和相应的边界条件。
例如弹性力学问题，热传导问题，电磁场问题等。
由于建立基本方程所研究的对象通常是无限小的单元，这类问题称为连续系统。
图1-3 V6引擎的局部下面是热传导问题的控制方程与换热边界条件：t Tc Q z T z y T y x T x ∂∂=+⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂ρλλλ （1- 1）初始温度场也可以是不均匀的，但各点温度值是已知的：() 00x,y,z T Tt ==（1- 2）通常的热边界有三种，第三类边界条件如下形式：()f T-T h nTλ=∂∂-（1- 3）尽管我们已经建立了连续系统的基本方程，由于边界条件的限制，通常只能得到少数简单问题的精确解答。
对于许多实际的工程问题，还无法给出精确的解答，例如，图1-3所示V6引擎在工作中的温度分布。
这为解决这个困难，工程师们和数学家们提出了许多近似方法。
在寻找连续系统求解方法的过程中，工程师和数学家从两个不同的路线得到了相同的结果，即有限元法。
有限元法的形成可以回顾到二十世纪50年代，来源于固体力学中矩阵结构法的发展和工程师对结构相似性的直觉判断。
从固体力学的角度来看，桁架结构等标准离散系统与人为地分割成有限个分区后的连续系统在结构上存在相似性。
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”的论文中首次提出了有限元（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法，导出标准的有限元过程来求解非结构问题。
我国的力学工作者为有限元方法的初期发展做出了许多贡献，其中比较著名的有：陈伯屏（结构矩阵方法），钱令希（余能原理），钱伟长（广义变分原理），胡海昌（广义变分原理），冯康（有限单元法理论）。
遗憾的是，从1966年开始的近十年期间，我国的研究工作受到阻碍。
1.2 有限元法的基本思路有限元法的基本思路可以归结为：将连续系统分割成有限个分区或单元，对每个单元提出一个近似解，再将所有单元按标准方法组合成一个与原有系统近似的系统。
下面用在自重作用下的等截面直杆来说明有限元法的思路。
等截面直杆在自重作用下的材料力学解答图1-4 受自重作用的等截面直杆 图1-5 离散后的直杆受自重作用的等截面直杆如图所示，杆的长度为L ，截面积为A ，弹性模量为E ，单位长度的重量为q ，杆的内力为N 。
试求：杆的位移分布，杆的应变和应力。
)()(x L q x N -= EAdxx L q EA dx x N x dL )()()(-==⎰-==x x Lx EA q EA dx x N x u 02)2()()(（1- 4）)(x L EAq dx du x -==ε )(x L AqE x x -==εσ等截面直杆在自重作用下的有限元法解答1）离散化如图1-5所示，将直杆划分成n 个有限段，有限段之间通过一个铰接点连接。
称两段之间的连接点为结点，称每个有限段为单元。
第i 个单元的长度为L i ，包含第i ，i+1个结点。
2）用单元节点位移表示单元内部位移第i 个单元中的位移用所包含的结点位移来表示，)()(1i iii i x x L u u u x u --+=+ （1- 5）其中i u 为第i 结点的位移，i x 为第i 结点的坐标。
第i 个单元的应变为i ε，应力为i σ，内力为i N ：iii i L u u dx du -==+1ε （1- 6）ii i i i L u u E E )(1-==+εσ（1- 7）ii i i i L u u EA A N )(1-==+σ（1- 8）3）把外载荷集中到节点上把第i 单元和第i+1单元重量的一半2)(1++i i L L q ，集中到第i+1结点上。
图1-6 集中单元重量4）建立结点的力平衡方程对于第i+1结点，由力的平衡方程可得：2)(11+++=-i i i i L L q N N（1- 9）令1+=i ii L L λ，并将（1- 8）代入得：221)11(2)1(i ii i i i i L EA q u u u λλλ+=-++-++（1-10）根据约束条件，01=u 。
对于第n+1个结点，2nn qL N =EAqL u u n n n 221=+-+ （1-11）建立所有结点的力平衡方程，可以得到由n+1个方程构成的方程组，可解出n+1个未知的接点位移。
1.3 有限元法的计算步骤有限元法的计算步骤归纳为以下三个基本步骤：网格划分，单元分析，整体分析。
1.3.1网格划分有限元法的基础是用有限个单元体的集合来代替原有的连续体。
因此首先要对弹性体进行必要的简化，再将弹性体划分为有限个单元组成的离散体。
单元之间通过单元节点相连接。
由单元、结点、结点连线构成的集合称为网格。
通常把三维实体划分成4面体或6面体单元的网格，平面问题划分成三角形或四边形单元的网格。
图1-7四面体四节点单元图1-8 六面体8节点单元图1-9 三维实体的四面体单元划分图1-10 三维实体的六面体单元划分图1-11 三角形3节点单元图1-12 四边形4节点单元图1-13 平面问题的三角形单元划分图1-14 平面问题的四边形单元划分1.3.2单元分析对于弹性力学问题，单元分析，就是建立各个单元的节点位移和节点力之间的关系式。
由于将单元的节点位移作为基本变量，进行单元分析首先要为单元内部的位移确定一个近似表达式，然后计算单元的应变、应力，再建立单元中节点力与节点位移的关系式。
以平面问题的三角形3结点单元为例。
如图1-15所示，单元有三个结点I、J、M，每个结点有两个位移u、v和两个结点力U、V。
图1-15 三角形3结点单元单元的所有结点位移、结点力，可以表示为结点位移向量（vector ）：结点位移{}⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=m m j j i i ev u v u v u δ结点力{}⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=m m j j i i eV U V U V U F 单元的结点位移和结点力之间的关系用张量（tensor ）来表示，{}[]{}e e e K F δ=（1-12）1.3.3整体分析对由各个单元组成的整体进行分析，建立节点外载荷与结点位移的关系，以解出结点位移，这个过程为整体分析。
再以弹性力学的平面问题为例，如图1-16所示，在边界结点i 上受到集中力iy ix P P ,作用。
结点i 是三个单元的结合点，因此要把这三个单元在同一结点上的结点力汇集在一起建立平衡方程。
图1-16 整体分析i 结点的结点力： ∑=++ee i i i i U U U U )()3()2()1(∑=++ee i i i i V V V V )()3()2()1(i 结点的平衡方程：⎪⎭⎪⎬⎫=∑∑=i y ee ie i x e iP V P U )()( （1-13）1.4有限元法的进展与应用有限元法不仅能应用于结构分析，还能解决归结为场问题的工程问题，从二十世纪六十年代中期以来，有限元法得到了巨大的发展，为工程设计和优化提供了有力的工具。
1.4.1算法与有限元软件从二十世纪60年代中期以来，进行了大量的理论研究，不但拓展了有限元法的应用领域，还开发了许多通用或专用的有限元分析软件。
理论研究的一个重要领域是计算方法的研究，主要有： 大型线性方程组的解法， 非线性问题的解法，动力问题计算方法。
目前应用较多的通用有限元软件如下表所列：另外还有许多针对某类问题的专用有限元软件，例如金属成形分析软件Deform、Autoform，焊接与热处理分析软件SysWeld等。
1.4.2应用实例有限元法已经成功地应用在以下一些领域：固体力学，包括强度、稳定性、震动和瞬态问题的分析；传热学；电磁场；流体力学。
转向机构支架的强度分析（刘道勇，东风汽车工程研究院动，用MSC/Nastran 完成）图1-17 转向机构支架的强度分析金属成形过程的分析（用Deform软件完成）分析金属成形过程中的各种缺陷。
图1-18 型材挤压成形的分析。
型材在挤压成形的初期，容易产生形状扭曲。
图1-19 螺旋齿轮成形过程的分析图1-20 T形锻件的成形分析焊接残余应力分析（用Sysweld完成）图1-21 结构与焊缝布置图1-22 焊接过程的温度分布与轴向残余应力热处理过程的分析BMW曲轴的感应淬火（Induction quenching of crankshafts at BMW，用SysWeld软件完成）在曲轴表面获得压应力，可以提高曲轴的疲劳寿命。
曲轴的有限元模型有限元模型的局部沿网格线52的残余应力分布，红线为预测的轴向应力与径向应力之差，黑点为实测值。
复杂形状工件的组织转变预测（石伟，用NSHT3D完成）预测工件的组织分布和机械性能二分之一工件的有限元模型淬火3.06 min 时的温度分布淬火3.06 min 时的马氏体分布参考文献S. I. Oh, W. T. Wu, K. Arimoto. Recent developments in process simulation for bulk forming processes. Journal of Materials Processing Technology III (2001) 2-9郭和德编. 有限单元法概论，清华大学，1998思考题：1）什么是有限元法？2）简述有限元法的基本思路。