第三章、有限元基本理论

1、第3章 有限元基本理论摘要：从一般的边值问题数值解理论出发，讲解了有限元法的基本过程和基本理论。有限元法基本过程包括问题几何区域的离散、近似解待定参数的确定、方程的建立等；基本理论包括单元的分类、单元形函数的性质、等参单元、单元积分和节点等。本章讲述的内容不受应用领域的限制。有限元法是为了解决结构分析而发展起来的一种新的数值方法。经过近50 年发展，它不但是结构分析强有力的工具，而且，在结构分析获得重大成功后，其理论也已日趋成熟，商务化软件系统也已有一定规模和数量，在其它领域边值问题的数值计算方面同样获得巨大成功。设由边界，其基本解为未知函数 u 的某一连续介质边值问题。在第一章中我们将此围成

2、区域问题转化成等效积分形式，并用加权残数法进行数值解； 第二章中对具有泛函极值形式的问题采用Litz进行数值解。但是以上两章并没有解决数值解中的试探函数（有限元中称形函数）的选取问题。有限元方法的关键是待定参数和形函数的选取及计算，那么采用有限元数值解法，需要经过哪些基本理论和过程呢§ 3.1有限元法概述区域的离散化将区域近似地离散成有限数量的，基本形状有一定限制的，尺寸远小于和的子区域集Elems ， Elems 称为有限单元（ Element ）集，它的元素称单元，记为e 或 ei ，对每个单元给予编号，即 Elems ei | i 1, 2, , M e (3.1.1)单元边界

3、 1节点e12e2单元 ee1e3e4 2图单元位置与形状由结点控制图单元协调性图区域离散成单元单元的基本形状可根据的几何维数选择，例如一维几何区域为线单元；二维区域可选择三角形或四边形单元； 而三维区域选择四面体、五面体和六面体单元等。图的平面区域被离散成有限个三角形单元，详细的单元分类和性质请见的讨论。控制单元形状和位置的点称为单元节点（element node,也有称结点或接点），简称节点（Node），例如图。节点的集合记为Nodes ，称节点集，并给予编号，即Nodes ni | i1,2, M n 围成单元的几何元素称为单元边界，例如图中四边形单元的四条边（edge）、四个顶点节点和

4、四个中间节点都属于单元边界。单元边界比之单元在几何维数上要低，根据几何维数不同，单元边界又可以分单元面、单元边、单元节点。在离散区域时，为了保证问题解的唯一性和连续性，两相邻单元的边界必须保持完全重合，即单元边界的节点被相邻单元完全共享。例如图中的节点1 被单元 e2 与 e4 共享，而节点 2 被 e1、e2 、e3 与 e4 四个单元共享。如何保证单元之间的问题解的连续性将在、节中讨论。确定待定参数集在第一章中已经指出，边值问题的数值解u ( x )是待定参数矢量集 的线性组合iu ( x ) u ( x )N (3.1.3)设节点 ni 的问题解的值为ui， u i组成的集合记为U ，即

5、U u i| i1, 2, M n (3.1.4)虽然 U 还不能完全等同近似解的待定参数集，但如果试探函数看成是对U 插值函数， 从矢量运算角度考虑， (3.2.3) 可以改写成和u ( x ) u ( x )NU(3.1.5)其中试探函数（插值函数）在有限元中称为形函数（shape function ），所以在得到U 后，就获得了问题的近似解，只是选定合适的形函数。例如，图由四个四边形单元组成固体力学平面应力应变问题，则U 由所有单元节点位移矢量所组成，简称位移矢量。所以问题的单元集、节点集和位移矢量分别为Elems e1, e2 , e3 , e4 Nodes n1, n2 , n9Uu

6、 1, u 2 , , u 9T, v),( u, v ) T( u , v ),( u229119789在 U 中，并不是所有参数是待定的。在本质边界上， 节点的 u 值34是确定，在混合边界上，节点的u 受到边界条件方程的约束。例如固456体力学问题位移解法中，位移边界上节点的位移值属于已知，混合边1223界上节点的位移受混合边界条件方程约束。但是不管节点的 u 值如何1图1.4单元集与结点集获得，的近似式仍然成立。所以在有限元方法对单元讨论，暂时把U 看成待定参数集，只是在后面求界待定参数方程组时，把已知的参数和约束方程代入方程组，从而减少方程组的数量，详细讨论见下章讨论。单元形函数的基

7、本要求在单元 e 中，设有 me 个节点。为了分析方便，节点的编号仍然采用1 至 me ，称之为局部编号以区别节点的整体编号。记第i 个单元局部节点的问题解在有限元法中采用以下假定：1） 单元内的问题解近似值只是该单元节点问题解的值所决定，与其他单元问题解的值无关。2） 问题解 u 的每个分量都采用相同的形函数。所以单元 e 内近似解的插值的矢量形式和分量形式为m em e(3.1.6)u e( x )N ej ( x ) u eju ie( x )N ej( x ) u ijej 1j 1其中 uie (x) 为单元 e 内问题解 u (x ) 的第 i 个分量， N ej (x) 为单元

8、e 的第 j 个节点的形函数，uije 是第 j 个节点的 u 的第 i 个分量。以上插值显然是 Langrange 插值法， 只保证了近似解的C0 阶连续。 如果要提高问题解连续性阶数，则需采用 Hermite插值法，这时以上第一条假定得取消。为了保证问题解的唯一性和单元之间问题解的C0 连续，（ 3.1.6）式形函数必须满足以下性质：1）唯一性： 在每个节点上插值函数的值有208e1ij5(3.1.7)1719N j ( x i )ij0ij91812761162）连续性： 单元边界 F （或是单元面，或是单元边，或是单元节点）411101315上的形函数值， 除了此边界上节点的形函数外，

9、其他节点的形函数必须为0，2143图六面体 20即结点单元3N ej (x )0xF , n ejF(3.1.8)单元边界 F 可以是。例如图三维20 节点的六面体单元，节点1、2、6 和 5A 2pA 1围成一个单元面，此面上的形函数值除了1、 2、 5、 6、 9、10、 13和 17节(x,y)A 320；节点 1、 2 组成的单元边，此1点的形函数外，其他节点的形函数必须等于图1.6 三角形单元1、 2 和 13 节点外，其他节点的形函数必须等于0。把边上的形函数值除了此原则推广到单元的节点上便得到以上第一条性质。满足(3.1.8) 式也就满足了单元与单元的交接边界上，问题解的插值是连

10、续的。3） 常数性： 如果单元上每个节点的问题解值相同，则此单元内每个坐标的问题解值也相同，即：m em em eu e ( x )CN ej ( x ) u ejN ej ( x ) CCN ej ( x )j1j1j1m e所以有：N ej ( x )1(3.1.9)j1常数性在固体力学中可解释为保证单元的插值能反映刚体位移。以上三点性质是选择单元形函数的必要条件。例 1. 研究图所示平面问题节点三角形单元的形函数。对于任意一点坐标p(x, y) ，节点 1、2、 3（按逆时针方向）的形函数分别取为A12A11xy1x1y1N11 x2y21x2y2( x, y)2 AA1x3y31x3y

11、3A22 A21xy1x1y1N21 x3y31x2y2( x, y)2 AA1x1y11x3y3A32A31xy1x1y1N3 (x, y)1x1y11x2y2(3.1.10-3)A2 A1x2y21x3y3其中 A 为三角形面积，A1, A2 , A3 含义见图。此3 个节点的形函数满足了以上提出的三点性质要求，显然它们都是坐标( x, y) 的线性插值函数。建立待定参数计算方程因为有限元法中单元内的问题解近似值只取决所在单元节点的问题解值，与其他节点的解值无关，所以对于迦辽金法，其区域和边界的等效积分 可以变成M eTA ( u ) d TB ( u ) de jTN iN i( N i

12、) A ( uje j1M eeje j)T( N iA (u i（ x ) d 0e jj 1e ji( x ) d i(1, 2,., n )(3.1.11)等效“弱”积分形式 (1.4.7)变成TTM eTeejjC (N i )D ( u )dE (N i )F (u )dC(N i)D(ui（ x ) dje j1M eTeje jE(N i)F (ui( x ) d 0 i (1, 2,., n )(3.1.12)ejj 1同理，对于最小势能原理(2.3.12)变成p ( u i )(Wf i u i ) d p i u i d fM eM fe(Wf iu i ) d p i u

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

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

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

16、而后两点表达了节点模型中的物理信息。有限元法在求解方程组 例如方程式（） 前，必须把任何边界条件简化成节点边界条件。例如固体力学中， 设某边界 i 属于某边界条件的几何元素，则 i 必须由有限而且完整的单元边界去离散。如果 i是位移边界，则离散后i 上的根据位移边界条件方程u u (x)在 i上获得 i上任一节点 n j 的位移值u (n j ) u (x j )n ji ， x j 是节点 n j 坐标如果 i是应力边界条件，则离散后把作用在i 上的作用力，按照静力等效原则分配到节点上。如果是混合边界条件，一些自由度方向（不一定是坐标方向）按照位移边界条件处理，留下的自由度方向按照应力边界处

17、理。例 2图 2.1a 为一平面应力问题的力学模型，几何形状为正方形。假如单元全部采用四节点的四边形单元，并有限元网格时在垂直方向均匀划分，则此问题的有限元模型为图。节点除了需要接受处理后的边界条件外， 还需要保存计算后的结果。 例如在固体力学问题， 对于位移是未知的节点保存位移（假如是位移解） ，位移已知的节点保存支座反力。而混合边界条件上部分自由度方向位移，其余保存支座反力。qPPPPP/2P/2P=qL/5La. 力学模型b. 有限元模型图 2.1有限元法的边界条件处理所以总结地讲，对于固体力学，所有节点（不仅仅是区域边界上的节点）需要保存的信息，作为已知条件， 需要保存已知或已知作用力

18、；作为计算结果保存位移或支座反力。这种思想理解为本质边界条件和自然边界条件，也可以应用到其他领域问题的处理上。§3.3 单元及其几何分类单元单元是有限单元的简称，单元是对问题区域的几何离散。在有限元计算过程中，在结构（结构决定了基本方程和边界条件方程的形式）确定的情况下单元还需要包含几何、材料参数、 物理参数三方面的信息。在几何上按照求解问题物体形状的几何测度（几何维），有限单元可分为一维、二维和三维单元。一维单元是对问题可以抽象为一维几何形状物体的离散，例如工程中的杆件结构、弦等； 二维单元是对问题可以抽象为二维几何形状物体的离散，例如平面问题、 薄板壳结构等； 同样的三维单元是对

19、三维几何形状物体的离散。其中三维问题最有广泛意义。除此之外，还有一些应用领域特殊的单元，例如在固a. 一 维 线 性 单b. 一维抛物线单元c. 一维三次抛物线单图一维单元(a)(b)(c)(d)(f)(g)(h)(i)a 线性三角形单元b抛物线三角形单元c含中间结点的抛物线三角形单元d三阶抛物线三角形单元e含中间结点三阶抛物线三角形单元f 线性 四边形单元g 抛物线四边形单元h 含中间结点的抛物线四边形单元i 三阶抛物线四边形单元(a)(b)(c)(d)(e)(f)(g)(h)(i)a 线性四面体单元d线性五面体单元g线性六面体单元b抛物线四面体单元e抛物线五面体单元h抛物线六面体单元c三次

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

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

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

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

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

25、 Hermite插值。显然 Lagrange 插值是属于 0 阶 Hermite插值。在有限元的一般系统中， 为防止计算规模的急剧增加和插值函数过于复杂，大都采用 Lagrange插值法。如果问题解的本身需要考虑的导数连续性，常常其导数也作为问题解的维。固体力学中的梁板结构，把转角也作为问题解的维。例如空间梁结构同时考虑了 u, v, w, x , y,z 和三个角位移；平板结构考虑了 w, x , y 等。uuue1e2e1e2e1e2a. 线 性 插 值b. 抛 物 线 插 值图 4.1 一 维 单 元 Lagngrange 插 值图 4.2一维 单 元 Hermite 插 值下面仅以一维

26、单元图示说明两种插值方法的区别，详细讨论请参考有关有限元书籍。对于一维单元，从图可以看出，Lagrange插值法就是简单的线性插值和抛物线，单元内问题的近似解可以写成形式。Lagrange 插值法的特点是形函数的数量与节点的数量相同，问题解整-101a. 线单元(-1,1)(1,1)个离散区域保持C0 阶连续，即在单元连接处只是连续但不可导。一阶 Hermite插值法可以写成(1,-1)22u i ) i(-1,-1)(0)(1)u ( x )H i( x ) u iH i( x )(b. 面单元i 1i 1x(4.4.1)1114或u ( x )N i ( x ) Q ii 1(4.4.2)

27、1Hermite多项式具有以下性质11H i(0)( xi )ijH i(0)( x)0c. 体单元xxi图5.1等参单元局部坐标(1) ( x)H i(1) ( xi )0H iijxxi(4.4.3)一阶 Hermite插值形函数具有三阶多项式。而二阶 Hermite插值形函数具有五阶多项式，形式可以写成u(x)2(0)2(1)ui)i2(2)2ui)i (4.4.4)Hi (x)uiHi(x)(xHi (x)(2i 1i 1i 1xHermite插值法不但形函数的阶数高，而且形函数的数量与节点数量也不同。§ 3.5等参单元所谓等参单元就是单元形函数的数量与节点数量一致，而且同一

28、类单元（例如面单元的8 节点四边形单元），则采用相同的形函数，单元的几何形状采用等参变换。例如上面例举的三节点三角形单元的形函数是以面积比 （ Li Ai / A）为参数的等参单元， 与具体的节点坐标无关。 等参单元属于 Lagrange 插值法，所以它符合小节提出的形函数三点基本要求。坐标变换对单元的形函数，不是取整体坐标(x, y, z) 的函数，而是统一取局部等参数坐标( , ) 的函数，即N iN i (,)避免了因坐标不同而形函数不同的困难。局部参数坐标更多的是采用正则坐标，其含义可见图， 面积参数坐标一般用在三角形单元和四面体单元上。在单元的局部参数坐标( , ) 上，问题解可表达

29、为m eu (,)N i (,) u ii1整体坐标 ( x, y, z) 与局部等参数坐标( , ) 的变换关系为mememexN i (,) xiyN i (,) yixN i (,) zii1i1i1导数变换Ni ( , ) 对可表示成N iN ixN iyN izxyz根据对称性，可以写出其他两局部参数坐标的导数，并写成矩阵形式NNNiiixyzNxxyzNyxyzNziiiJNxNyNziii式中 J 为 Jacobi 矩阵，记为 J( x, y, z)( ,) ，利用 可以计算出Jacobi矩阵。求逆得TTN iN iN iJ 1N iN iN ixyz积分变换在形成近似解的计算方

30、程时，常常用到对单元体的体积分、单元面的面积分和单元边的边积分。例如固体力学中的单元刚度矩阵系数是由单元体积分形成（见3.1.4 式），右端的力矢量元素由单元体力的体积分、单元面上面分布力的面积分和单元边上线分布力的线积分形成（见式）。对于单元体局部参数坐标的体积微元dVd(d d?)(3.5.7)而局部参数坐标矢量微分d x diy djz dkd x diydjz dk(3.5.8)d ?x diy djz dk所以xyzdVxyzdd dJd d d(3.5.9)xyz对单元面积分，取某参数坐标为常数，例如在c1 面上面微元为dAd d?cAd d把的第二和第三子式代入可计算得A(yzy

31、z ) 2 +(zxzx ) 2 +(xyxy )2对于单元边积分，取两参数坐标为常数，例如在c11,c21 边上dSd?c1 ,c2Sd把的第三子式代入可计算得S( x )2 +( y ) 2 + ( z ) 2有了以上的单元体、面和边积分变换，可以把对整体坐标的积分变换成局部参数坐标的积分( ,)111e GzdV1 11 G ( x( , , ), y( , , ), z( , , ) J d d dx y11c面上 H ( x, y, z) dA1 1 H ( x(c, , ), y(c, , ), z(c, , ) Ad dI ( x, y, z)dA1I ( x(c1 ,c2 ,

32、), y(c1 ,c2 , ), z(c1, c2 , )Sdc1 , =c2 边上1以上推导是建立在三维坐标系之上，如果是二维、 一维问题， 以上公式需要进行退化处理，例如二维问题 Jacobi矩阵为J 2xyxy对、 和这样的数值积分，虽然许多领域问题的被积函数也可常常写出其解析形式，但是因为单元形式和领域问题的多样性， 加上数值积分仍然能保持极高的精度， 所以对这些积分也采用数值解，一般采用高斯积分法。具体计算方法请参考有关有限元书籍。§3.6 等参单元形函数的构造技巧前面讲了等参单元的坐标变换、形函数导数求解和积分计算，但是还没有构造出等参单元的形函数。等参单元的形函数需要一

33、定技巧，非常有规律，很容易掌握。三节点三角形等参单元的形函数图所示，对于节点1，有单元边 10 没有通过它，所以取N1(,)(1)显然满足了N1 ( j ,j )0j1利用 N1 (0,0)1 得N1 (,)(1)(3.6.1-1)同理得N 2 (,)(3.6.1-2)3.6.2N 3 (,)3(0,1)3(0,1)(0,0.5)6(0.5,0.5)512142(0,0)(1,0)(0,0)(0.5,0)(1,0)图6.1. 三节点三角形等参单元图6.2.六节点三角形等参单元4103473(-1,1)(1,1)(-1,1)(1,1)108610(-1,-1)10(1,-1)(-1,-1)5(1,-1)1212图5.2.四节点四边形等参单元图5.3.八节点四边形等参单元四节点四边形等参单元的形函数图所示，对于节点1，有两条单元边 10 和 10 没有通过