、有限元基本理论

1、Elems , Elems称为有限单元（Hement）集,它的元素称单元，记为 即Elems e | i1,2,|,Me单元边界单元e图1.2单元位置与形状由结点控制第3章有限元基本理论摘要：从一般的边值问题数值解理论出发，讲解了有限元法的基本过程和基本理论。有限元法基本过程包括问题几何区域的离散、近似解待定参数的确定、方程的建立等；基本理论包括单元的分类、单 元形函数的性质、等参单元、单元积分和节点等。本章讲述的内容不受应用领域的限制。有限元法是为了解决结构分析而发展起来的一种新的数值方法。经过近50年发展，它不但是结构分析强有力的工具，而且，在结构分析获得重大成功后，其理论也已日趋成熟，商

2、务化软件系统也已有 一定规模和数量，在其它领域边值问题的数值计算方面同样获得巨大成功。设由边界r围成区域q,其基本解为未知函数 u的某一连续介质边值问题。在第一章中我们将此问题转化成等效积分形式， 并用加权残数法进行数值解； 第二章中对具有泛函极值形式的问题采用Utz进行数值解。但是以上两章并没有解决数值解中的试探函数（有限元中称形函数）的选取问题。有限元方法的关键是待定参数和形函数的选取及计算，那么采用有限元数值解法，需要经过哪些基本理论和过程呢？§ 3.1 元法概述3.1.1 区域的离散化将区域 Q近似地离散成有限数量的，基本形状有一定限制的，尺寸远小于Q和r的子区域集e或e ,

3、对每个单元给予编号，单元的基本形状可根据 Q的几何维数选择，例如一维几何区域为线单元;边形单元；而三维区域选择四面体、五面体和六面体单元等。图1.1的平面区域被离散成有限个三角形单元，详细的单元分类和性质请见3.3的讨论。控制单元形状和位置的点称为单元节点（element node,也有称结点或接点），简称节点（Node）,例如图1.2。节点的集合记为 Nodes,称节点集，并给予编 号，即Nodesni | i 1,2,|,Mn|（3.1.2）围成单元的几何元素称为单元边界，例如图1.2中四边形单元的四条边（edge）、四个顶点节点和四个中间节点都属于单元边界。单元边界比之单元在几何维数上要

4、低，根据几何维数不同，单元边界又可以分单元面、单元边、单元节点。在离散区域时，为了保证问题解的唯一性和连续性，两相邻单元的边界必须保持完全重合，即单元边界的节点被相邻单元完全共享。例如图1.3中的节点1被单元金与e4共享，而节点2被e1、e2、e3与e4四个单元共享。如何保证单元之间的问题解的连续性将在3.3、3.4节中讨论。3.1.2 确定待定参数集在第一章中已经指出，边值问题的数值解u (x)是待定参数矢量集 四的线性组合u (x) u (x) N a(3.1.3)设节点n的问题解的值为Ui , 5组成的集合记为 U ,即U uJi 1,2,|(|,M n(3.1.4)虽然U还不能完全等同

5、近似解的待定参数集，但如果试探函数看成是对U插值函数，从矢量运算角度考虑，(3.2.3)可以改写成和u (x) u (x) NU(3.1.5)其中试探函数(插值函数)在有限元事称为形函数( shape function ),所以在得到 U后，就获得了问 题的近似解，只是选定合适的形函数。例如，图1.4由四个四边形单元组成固体力学平面应力应变问题，则 U由所有单元节点位移矢量所组成，简称位移矢量。所以问题的单元集、节点集和位移矢量分别为Elems ee2,e3,e4Nodes nn2, 飞U u1,u2,|,u9 T （ Ui，Vi）,（U2”2）J|,（U9，V9）T2图1.4单元集与结点集在

6、U中，并不是所有参数是待定的。在本质边界上，节点的u值是确定，在混合边界上，节点的 u受到边界条件方程的约束。例如固 体力学问题位移解法中，位移边界上节点的位移值属于已知，混合边 界上节点的位移受混合边界条件方程约束。但是不管节点的u值如何获得，(3.1.5)的近似式仍然成立。所以在有限元方法对单元讨论， 暂时把U看成待定参数集，只是在后面求界待定参数方程组时，把 已知的参数和约束方程代入方程组，从而减少方程组的数量，详细讨 论见下章讨论。3.1.3 单元形函数的基本要求在单元e中，设有me个节点。为了分析方便，节点的编号仍然采用1至me,称之为局部编号以区别节点的整体编号。记第i个单元局部节

7、点的问题解在有限元法中采用以下假定：1)单元内的问题解近似值只是该单元节点问题解的值所决定，与其他单元问题解的值无关。2)问题解u的每个分量都采用相同的形函数。所以单元e内近似解的插值的矢量形式和分量形式为m em eu e (x)N：(x)u：u：(x)N：(x)u；(3.1.6)j 1j 1其中ue(x)为单元e内问题解u(x)的第i个分量，N：(x)为单元e的第j个节点的形函数，u：是第j个 节点的u的第i个分量。以上插值显然是Langrange插值法，只保证了近似解的 C0阶连续。如果要提高问题解连续性阶数，则需采用Hermite插值法，这时以上第一条假定得取消。为了保证问题解的唯一性

8、和单元之间问题解的C0连续，（3.1.6 ）式形函数必须满足以下性质:1）唯一性：在每个节点上插值函数的值有n e(Xi)(3.1.7)图1.5六面体20结点单元2）连续性：单元边界（或是单元面，或是单元边，或是单元节点）上的形函数值，除了此边界上节点的形函数外，其他节点的形函数必须为 0,即N：(x)0 x F,n： F(3.1.8)1、2、6和5围成一个单元面，此面上的形函数值除了1、2、5、6、9、10、13和17节点的形函数外，其他节点的形函数必须等于0;节点1、2组成的单元边，此边上的形函数值除了 1、2和13节点外，其他节点的形函数必须等于把此原则推广到单元的节点上便得到以上第一条

9、性质。满足 满足了单元与单元的交接边界上，问题解的插值是连续的。3)所以有:0。(3.1.8) 式也就常数性：如果单元上每个节点的问题解值相同，则此单元内每个坐标的问题解值也相同,即:u e(x)mN：(x)u：n e(x)C1c n e(x) j 1e(x )(3.1.9)单元边界F可以是。例如图1.5三维20节点的六面体单元，节点常数性在固体力学中可解释为保证单元的插值能反映刚体位移。以上三点性质是选择单元形函数的必要条件。例1.研究图1.6对于任意一点坐标所示平面问题3节点三角形单元的形函数。p(x, y),节点 1、2、3 （按逆时针方向）的形函数分别取为Ni(x, y)AiA2A2A

10、x2x3yy2V3x1x2x3y1V2V3(3.1.10-1)N2(x, y)A2A2A22AN3(x, y)A3A2A32Axy/1xiy1x3V31x2y2xiy/1x3y3xy/1xYixiy1y2X2y2/1x3Y3111(3.1.10-2)(3.1.10-3)其中A为三角形面积,A,4,A3含义见图1.6。此3个节点的形函数满足了以上提出的三点性质要求,显然它们都是坐标（x, y）的线性插值函数。3.1.4 建立待定参数计算方程因为有限元法中单元内的问题解近似值只取决所在单元节点的问题解值，与其他节点的解值无关, 所以对于迦辽金法，其区域Q和边界的等效积分(1.4.6)可以变成M e

11、TT _ej tejN i A(u)d Q N i B(u)dr(Ni j) A(Uij(x)dQQrj i ejMe e T-e ,(Ni j) A (Uij(x)dT 0i (1,2,., n) (3.1.11)j i巧等效“弱”积分形式(1.4.7)变成MeCT(Ni)D(U)dQET(Ni)F(U)dr e CT (N ej )D(ujx)d Q,j 1 jMe t e；e；(3.1.12)(3.1.13)E (Ni j)F(Ui j(x)dr 0 i (1,2,., n)j 1 ej同理，对于最小势能原理 (2.3.12)变成p(Ui) q(WfiUi)dQr PiUidrM eM

12、fe(W fi u i) d Qpiu i d rj 1 ejj 1 efj其他几种变分也可以化成对区域单元的积分和对边界上单元边界面的积分之和。对于导数不超过两阶的物理问题，不管采用哪种形式建立的数值计算方程，最终可以得到KU P(3.1.14)这样形式的方程。如果是对于固体力学问题，K为刚度矩阵，其中的系数由单元的材料参数D (例如弹性力学中的E与)、物理参数T (例如板结构中的板厚度)、形函数和形函数导数组合而成的代数 式对单元的积分，而且是若干单元的积分和，即kjeR(D,T,N ,N')dQ(3.1.15)而P矢量的分量由单元的体力与形函数组合代数式的积分，应力边界上面力与形

13、函数组合在单元面上的积分所合成，即PieS(f ,N )d Qe Q(p, N)dr(3.1.16)具体如何获得，每个量代表什么物理意义将在下章讨论。虽然形函数为一般的多项式，但是由于形函数形式很多，对不同的类型、不同的结构和材料(3.1.12 )和(3.1.13 )表达形式都有所不同，所以对它们的单元积分也是采用数值积分方法，并以高斯积分方法为多数。详细刚度矩阵、力矢量、高斯积分、 位移与约束条件的解除、方程的求解等内容将在后面章节中讨论。§ 3.2 点与节点密切相关的一个重要的概念是自由度，所谓自由度是问题解 u的维数。自由度的多少也同时决定了边界条件维数。在固体力学中，最多自由

14、度可达6个，三个线位移u,v,w和三个角位移 x, y, z,对应的应力边界条件是线力 X,Y,Z和力矩M x,My,Mz, 一般结构是以上 6这个自 由度的子集。例如平面应力应变结构为u,v；平板结构为w, x, y；三维实体结构为u,v,w；平面框架结构为u,v, z；三维框架结构为全部 6个等。当然结构不同建立的基本微分方程也不同，从 而导致对应单元的计算方法不同，例如梁单元、平面应力单元、三维实体单元等等。节点在有限元法中承载许多模型方面信息：1）表达位置的坐标；2）连接单元；3）在它之上施加边界条件；4）放置数值分析的计算结果。前两点表达了节点组成有限元网格的几何信息，而后两点表 达

15、了节点模型中的物理信息。有限元法在求解方程组例如方程式（3.1.13 ）前，必须把任何边界条件简化成节点边界条件。例如固体力学中，设某边界属于某边界条件的几何元素，则必须由有限而且完整的单元边界去离散。如果n是位移边界，则离散后n上的根据位移边界条件方程u u（x）在上获得n上任一节点nj的位移值u（nj） u（Xj）nj 口，Xj 是节点 1 坐标如果是应力边界条件，则离散后把作用在口上的作用力，按照静力等效原则分配到节点上。如果是混合边界条件，一些自由度方向（不一定是坐标方向）按照位移边界条件处理，留下的自由度方向按照 应力边界处理。例2.图2.1a为一平面应力问题的力学模型，几何形状为正

16、方形。假如单元全部采用四节点的四边形单元，并有限元网格时在垂直方向均匀划分，则此问题的有限元模型为图2.1b。节点除了需要接受处理后的边界条件外，还需要保存计算后的结果。例如在固体力学问题， 对于位移是未知的节点保存位移（假如是位移解），位移已知的节点保存支座反力。而混合边界条件上部分自由度方向位移，其余保存支座反力。b.有限元模型图2.1有限元法的边界条件处理所以总结地讲，对于固体力学，所有节点（不仅仅是区域边界上的节点）需要保存的信息，作为已知条件，需要保存已知或已知作用力；作为计算结果保存位移或支座反力。这种思想理解为本质边界条件和自然边界条件，也可以应用到其他领域问题的处理上。

17、7; 3.3 元及其几何分类3.3.1 单元单元是有限单元的简称，单元是对问题区域的几何离散。在有限元计算过程中，在结构（结构决定 了基本方程和边界条件方程的形式）确定的情况下单元还需要包含几何、材料参数、物理参数三方面的信息。在几何上按照求解问题物体形状的几何测度（几何维），有限单元可分为一维、二维和三维单元。一维单元是对问题可以抽象为一维几何形状物体的离散，例如工程中的杆件结构、 弦等；二维单元是对问题可以抽象为二维几何形状物体的离散，例如平面问题、薄板壳结构等；同样的三维单元是对三维几何形状物体的离散。其中三维问题最有广泛意义。除此之外，还有一些应用领域特殊的单元，例如在固a.一维线性单

18、元b.一维抛物线单元c. 一维三次抛物线单元图3.1. 一维单元(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)a,线性三角形单元b.抛物线三角形单元 c.含中间结点的抛物线三角形单元d.三阶抛物线三角形单元e.含中间结点三阶抛物线三角形单元f.线性四边形单元g,抛物线四边形单元h.含中间结点的抛物线四边形单元i,三阶抛物线四边形单元j.含中间结点三阶抛物线四边形单元图3.2.二维单元a.线性四面体单元d.线性五面体单元g.线性六面体单元b.抛物线四面体单元 e.抛物线五面体单元 h.抛物线六面体单元图3.4,三维单元c.三次抛物线四面体单元 f.三次抛物线五面体单元 i.三次抛物

19、线六面体单元体力学有限元方法中，存在质点单元、刚体单元、弹簧单元、阻尼单元、粘弹性单元和伪单元等一些特 殊的单元。基本的有限单元除了按照几何测度分类外，根据单元的插值函数多项式阶数的需求，在单元的边界线（见图3.1 ）上，可以有两个节点、三个节点甚至四个节点，分别称线性单元、抛物线单元和三次抛 物线单元。边界上的节点的数量越多，插值函数多项式的阶数也越高，问题求解的精度也越高，但是求 解问题的未知数数量也随之增加。对于特殊情况，除了单元边界上存在节点外，单元内部也可能存在节点（见图3.2 ）。每一个单元必须选择一种材料（一种材料可以有多个单元），在固体力学中，材料参数是根据材料本构关系需要而确

20、定需要什么参数，与问题结构无关。材料性质可以分线弹性材料、弹塑性材料、蠕变 材料等。不同材料有不同的材料选择模式。对于各向异性材料需要输入不同方向的材料参数。材料性质是由材料参数表描述，材料的参数可以独立与单元存在，可以在单元生成之前建立。物理参数是对单元几何特性的补充，例如二维单元的厚度、梁单元横截面的性质等。单元厚度是二维单元向第三个几何方向的几何补充，梁横截面是一维单元向第二、第三个几何方向的几何补充。与材料特性一样，物理参数也是单元计算中需要的参数，可以在网格生成前建立。 但并不是所有单元都需要物理参数，是否需要取决求解的问题结构，对于平面应变、板壳单元，需要参考单元厚度物理参数，对

21、于梁单元，需要参考梁横截面物理参数，而对于平面应力、轴对称、三维等问题，则不需要物理参数。3.3.2 单元几何分类一、一维基本有限单元在固体力学有限单元方法中，一维单元主要用来解决杆件与绳结构问题，像桁架、框架、网架和悬 索等结构，所以一维单元在土建工程中有着广泛的应用。一维线性单元、抛物线单元和三次抛物线单元（见图3.1 ）。二、二维基本有限单元二维基本有限单元分三角形和四边形两种基本几何形状，平面单元适应平面问题和空间曲面的几何区域离散。在固体力学中，平面单元用在平面应力、平面应变、轴对称、板壳等结构等的有限元方法中。 三角形和四边形型单元又分线性单元、抛物线单元和三次抛物线单元（见图3.

22、2）。三、三维基本有限单元三维基本有限单元分四面体、五面体和六面体三种基本几何形状，原则上讲，三维有限基本单元内部也可以有中间节点， 但是使用情况比较少， 在图3.3中没有绘出。图3.3中把单元边界都绘制成线性 边界，实际拟合是可以用曲线和曲面。三维单元适应任何能够适用有限元方法三维问题的几何区域离散问题。3.3.3 单元的几何协调条件单元是对求解问题几何区域的离散，离散后，问题的几何区域被单元的集合所替代。为了保持所求解状态结果的连续性和一定的精度，单元的划分并不是随意的，必须满足一定的几何协调条件和形状要求。显然，单元的几何性质是由节点控制的，节点的状态解构成了有限元的最终解。在生成的网格

23、中， 单元与节点必须保证协调性：1 .对于连续的区域，在离散后，单元与单元之间不能重叠，非边界单元的单元边界与另外单元的边界公享，公享部分的单元边界包括公享单元面、单元边和单元节点；2 .单元的形状不能太奇形，理想的形状是等边和等角度的几何形状。§ 3.4 一维Lagrange插值法与 Hermite插值法插值函数，即单元的形函数，常用的有 Lagrange插值法和Hermite插值法两种。Lagrange插值 法只考虑问题解的节点值，只保证了近似解的C0阶连续；而Hermite插值法在考虑问题解的节点值时， 同时考虑n阶导数值，所以达到了近似解的Cn阶连续，称为n阶Hermite插

24、值。显然Lagrange插值是属于0阶Hermite插值。在有限元的一般系统中，为防止计算规模的急剧增加和插值函数过于复杂，大都采用Lagrange插值法。如果问题解的本身需要考虑的导数连续性，常常其导数也作为问题解的维。固体力学中的梁板结构，把转角也作为问题解的维。例如空间梁结构同时考虑了u,v,w, x，y, z和三个角位移；平板结构考虑了w, x, y等。图4.2 维单元Hermite插值图4.1 一 维单元Lagngrange插值下面仅以一维单元图示说明两种插值方法的区别，详细讨论请 参考有关有限元书籍。对于一维单元，从图 4.1可以看出，Lagrange插值法就是简单 的线性插值和抛

25、物线，单元内问题的近似解可以写成(3.1.6)形式。Lagrange插值法的特点是形函数的数量与节点的数量相同，问题解整个离散区域保持C。阶连续，即在单元连接处只是连续但不可导。*.M-101一阶 Hermite插值法可以写成(4.4.1)u(x)(0) i(x)UiH i(x)()i xa.线单元b.面单元u(x)N i(x)Qii 1(4.4.2)Hermite多项式具有以下性质也3)。Hi(0)(x)xiHi(1)(x)0Hi(x)1c.体单元xi(4.4.3)一阶Hermite插值形函数具有三阶多项式。数具有五阶多项式，形式可以写成而二阶Hermite插值形函图5.1等参单元局部坐标H

26、(x)(i 12 u、，、T)i (4.4.4) x22 ，八 uu(x)Hi(0)(x)uHi(x)(Tii 1i 1xHermite插值法不但形函数的阶数高，而且形函数的数量与节点数量也不同。§ 3.5 参单元所谓等参单元就是单元形函数的数量与节点数量一致，而且同一类单元(例如面单元的8节点四边形单元)，则采用相同的形函数，单元的几何形状采用等参变换。例如上面例举的三节点三角形单元 的形函数是以面积比(L A/A)为参数的等参单元，与具体的节点坐标无关。 等参单元属于 Lagrange 插值法，所以它符合 3.3.1小节提出的形函数三点基本要求。3.5.1 坐标变换对单元的形函数

27、，不是取整体坐标(x, y,z)的函数，而是统一取局部等参数坐标(，)的函数,即Ni W(,)(3.5.1)避免了因坐标不同而形函数不同的困难。局部参数坐标更多的是采用正则坐标，其含义可见图5.1 ,面积参数坐标一般用在三角形单元和四面体单元上。在单元的局部参数坐标(，)上，问题解可表达为 meu( , , ) Ni( , , )Uii 1(3.5.2)整体坐标(x, y, z)与局部等参数坐标,)的变换关系为mex Ni(i 1me)xi y Ni( i 1me)yi x Ni( , , )zii 1(3.5.3)3.5.2导数变换Ni(,)对可表示成N j N i x N i y N i

28、z xyz(3.5.4)根据对称性，可以写出其他两局部参数坐标的导数，并写成矩阵形式Ni_x_y_z_NiNixxNixyzNiJNiyyN ixyzN iN i(3.5.5)式中J为Jacobi矩阵，记为J (x, y, z)/ ()，利用(3.5.3)可以计算出Jacobi矩阵。求(3.5.5)逆得(3.5.6)3.5.3积分变换在形成近似解的计算方程时,常常用到对单元体的体积分、如固体力学中的单元刚度矩阵系数是由单元体积分形成(见单元面的面积分和单元边的边积分。例3.1.4式)，右端的力矢量元素由单元体力的体积分、单元面上面分布力的面积分和单元边上线分布力的线积分形成(见 对于单元体局部

29、参数坐标的体积微元3.1.5dVdE(dr d?)(3.5.7)而局部参数坐标矢量微分Ed-dd t)Ed-d(3.5.8)d?Xd-d所以dV(3.5.9)对单元面积分,取某参数坐标为常数，例如在1面上面微元为dAd” d?(3.5.10)c Ad d(3.5.11)对于单元边积分，取两参数坐标为常数，例如在CI1,c21 边上dSd?C1 ,C2Sd(3.5.12)把(3.5.8)的第三子式代入可计算得(3.5.13)有了以上的单元体、面和边积分变换，可以把对整体坐标的积分变换成局部参数坐标的积分eG(x, y,z)dV1111 1G (x( , , ), y( , , ),z( , ,

30、) J d d d(3.5.14)1c 面 ±H(x,y,z)dA 111H (x(c, , ), y(c, , ),z(c, , )Ad d(3.5.15)d, =c2边上 I(x,y,z)dA11I(x(q,c2, ),y(c1,c2, ),z(Ci,C2, )Sd(3.5.16)以上推:建立在三维坐标系之上,如果是二维、一维问题，以上公式需要进行退化处理，例如二维问题Jacobi矩阵为(3.5.17)(3.5.15) 、(3.5.16)和(3.5.17)这样的数值积分，虽然许多领域问题的被积函数也可常常写出 其解析形式，但是因为单元形式和领域问题的多样性，加上数值积分仍然能保持

31、极高的精度，所以对这些积分也采用数值解，一般采用高斯积分法。具体计算方法请参考有关有限元书籍。§ 3.6 参单元形函数的构造技巧前面讲了等参单元的坐标变换、形函数导数求解和积分计算，但是还没有构造出等参单元的形函 数。等参单元的形函数需要一定技巧，非常有规律，很容易掌握。3.6.1 三节点三角形等参单元的形函数图6.1所示，对于节点1,有单元边10没有通过它，所以取Ni( , )(1)显然满足了N1( j, j)0 j 1利用 N1(0,0)1 得 N1( , ) (1)(3.6.1-1)同理得N2(,)(3.6.1-2)(3.6.1-3)图5.3.八节点四边形等参单元3.6.2 四节点四边形等参单元的形函数图5.2所示，对于节点1 ,有两条单元边10没有通过它，所以取N1( , )(1)(1)显然满足了N1( j, j)0利用N1( 1, 1) 1得1N1(,)4(1)(1结合1节点的坐标值，并用同样方法得各节点的形函数1Ni( , )(1 i )(14i