有限元法基本原理及应用第3章重庆大学龙雪峰

有限元法作为一种数值分析的方法现已被广泛应用于固体力学、流体力学、传热学、电磁学等多个领域。这种方法最早萌芽于弹性力学问题的求解。以位移、力或者力和位移的组合为未知量分别把有限元方法分为有限元力法、有限元位移法和混合有限元法。有限元位移法由于计算过程的系统性和通用性而被广泛应用，目前大多数商业有限元软件在对弹性力学问题的分析中大多采用位移法。本章介绍求解弹性力学问题的有限元位移法。第3章弹性力学有限元法3.1有限元法求解问题的基本步骤有限元法作为一种分析方法，即可以分析比较简单结构的问题如一根轴、一块板等，也可以分析比较复杂结构的如潜艇、体育场馆等，但其分析问题的基本步骤大同小异，可归纳如下：1．建立力学模型首先，在尽可能反映工程实际的基础上对工程结构进行几何简化，一般分为一维问题、二维问题和三维问题。其中二维问题要分清是平面应变问题还是平面应力问题，如果结构具有对称性可利用来简化计算。其次，要把结构上作用的载荷根据作用区域的大小简化为集中力、线载荷、面载荷或体力等。第三，根据边界约束和运动情况提出边界条件。2．连续体离散化根据所建立的力学模型选用合适的单元将连续体划分为有限个具有规则形状的单元集合。单元的选取应视所分析问题的性质、规模和精度要求而定，如杆单元有二维单元、三维单元，二维问题可选三角形单元、四边形单元，三维问题有四面体单元、六面体单元等。3．单元分析单元分析包括位移模式选择，单元力学分析两个内容。位移模式也称位移函数或插值函数，在有限元位移法中是以节点位移为基本未知量，再由这些节点位移插值得到单元内任意一点的位移值。单元的位移模式一般采用多项式，因为多项式计算简便，并且随着项数的增加，可以逼近任何一段光滑的函数曲线。单元力学分析根据所选单元的节点数和单元材料性质，应用弹性力学几何方程和物理方程得到单元刚度矩阵。由于连续体离散化后假定力是通过节点在单元间传递的，因此要利用插值函数把作用在单元上的体积力、面积力和集中力按静力等效原则移到节点上。4．整体分析和有限元方程求解由已知的单元刚度矩阵和等效节点载荷列阵组装成整个结构的整体刚度矩阵和载荷列阵，得到一个由总体刚度矩阵[K]、总载荷向量{F}和整体节点位移向量{δ}表示的平衡方程式：[K]{δ}＝{F}。引进位移边界条件后求解得到整体节点位移向量。有限元离散方程是一个代数方程组，代入边界条件处理以后的刚度矩阵是一个正定的对称稀疏方阵，这样一个代数方程组可以用高斯消元法、三角分解法、波前法和雅可比迭代法等多种方法求解。5．结果后处理和分析求解线性方程组得到位移矢量后，由几何和物理关系可以得到应变和应力。由于应变（应力）来自位移的微分可能导致单元间应力不连续，这会使应力计算误差较大，要在节点附近进行平均化处理。通过后处理还可得到位移、应变和应力的最大最小值及其所在位置以及主应力、主应变或其它定义的等效应力。结果的输出可以应用图表、动画等各种方式。最后还要对这些结果进行分析以指导工程设计、产品开发等等。3.2连续体离散化一个结构经过力学简化以后通常变为梁杆、板壳、实体或它们的组合，要进行有限元分析必须进行离散化，即建立有限元模型，就是把结构用不同的单元划分为只在结点上彼此连接的有限个单元组成的离散体。为了使有限元模型能够准确地体现实际结构，必须选择适当的单元类型，下面简要介绍几种常用的单元类型。

3.2.1杆状单元平面杆单元b)空间杆单元图3.1杆单元平面梁单元b)空间梁单元图3.2梁单元平面桁架、空间网架是工程中常见的结构，这种结构的主要受力构件是杆件，其截面尺寸远小于轴向尺寸。这类单元主要有平面杆单元和空间杆单元，如图3.1所示。杆单元有两个节点，只能承受轴向的拉压载荷，平面杆单元每个节点有两个自由度，空间杆单元每个节点有三个自由度。杆状单元还有用于求解平面刚架和空间刚架问题的梁单元，如图3.2所示。常用梁单元有两个节点，平面梁单元每个节点有三个自由度，两个线位移和一个角位移，可承受平面内的体力、集中力、分布力和垂直平面的弯矩的作用，只考虑单元的拉压和一个方向的弯曲变形。空间梁单元每个节点有六个自由度，三个线位移和三个角位移，一般可承受各个方向的集中力、分布力和弯矩，能够考虑单元的弯曲、拉压、扭转变形。当切应变对挠度的影响比较大时，必须修正单元刚度矩阵以考虑切应变。3.2.2平面单元平面问题是三维问题根据结构和外力特点，近似简化的结果。这种简化不仅为结构的理论分析带来了很大的便利，也为结构的有限元分析带来了便利。弹性力学平面问题根据其应力与应变特点分为平面应力问题和平面应变问题，在进行有限元分析时他们所采用的基本单元是相同的，区别仅在于弹性矩阵不同。平面单元属于二维单元，只能承受单元平面内的分布力和集中力，不能承受面外载荷。常见的平面单元有三角形单元和四边形单元，矩形单元是经常采用的特殊的四边形单元，如图3.3所示。常用的三角形单元有三个节点，每个节点有两个自由度，因此只能采用线性模式，导致应变值为一常量，因此也称为常应变单元或常应力单元。三角形单元计算精度较差，但可适用于复杂不规则形状的结构。为了提高精度也发展了六节点三角形单元，其位移模式是二次的，即可适用于复杂形状也有比较高的精度，但增加了未知数数量。四边形单元有四个节点，每个节点也有两个自由度，采用双线性位移模式，计算精度较高。3.2.3薄板弯曲单元和薄板单元板壳结构是工程上经常采用的一类结构形式，其特点是在一个方向上的尺度远小于另外两个方向的，通过一定的弹性力学假设，简化为特殊的二维结构，即便于解析方法求解也给有限元分析带来了很大的方便。这类结构通常有压力容器、舰船外壳，体育馆屋顶，建筑物楼板等等。薄板弯曲单元通常也有三角形单元和四边形单元两种，矩形单元为后者的特殊形式，通常三角形单元有三个节点，四边形单元有四个节点。主要承受横向载荷和绕水平轴的弯矩。如果挠度与板厚相比是小量时，板的中面应变可以忽略不及，如图3.4所示单元的每个节点有三个自由度，这样的单元一般称为薄板弯曲单元。如果挠度与板厚相比不再为小量，如金属板，当挠度w与板厚t的关系在范围内，板的中面应变就不能忽略，如图3.5所示，面内的两个自由度也要一并考虑，导致单元的每个节点上就要有五个自由度，此类单元一般称为薄板单元。四边形弯曲单元b)三角形弯曲单元图3.4薄板弯曲单元四边形薄板单元b)三角形薄板单元图3.5薄板单元3.2.4多面体单元对于实体结构就要用三维多面体单元进行分析，如机床工作台、机械基础件等等。常用的三维多面体单元有四节点四面体单元和八节点六面体单元，六面体单元有规则六面体和不规则六面体，如图3.6所示。为了提高精度也有八节点四面体单元和二十节点六面体单元。四面体单元b)六面体单元图3.6多面体单元3.2.5等参单元一般而言，随着单元节点数目的增加，位移模式的次数增加，相应的计算精度也会得到提高，同时计算规模相应增大，计算时间增加。在求解实际问题时人们总希望用最少的单元实现比较高的计算精度，而且所选用的单元对复杂结构也有比较好的适应性，因此将形状规则的单元变换为边界为曲线或曲面的单元是比较好的选择。最简便易行的方法是把单元形状的变化和单元内位移函数的变化用相同数目的节点参数和相同的插值函数进行变换，即所谓等参变换，采用等参变换的单元称为等参单元。下面以四节点四边形单元为例说明等参单元。如图3.7所示，一个以ξOη为局部坐标的矩形单元，单元边界上的ξ、η值为±1，如果能够通过坐标变换使其与整体坐标系中的一个单元建立映射关系，即矩形单元的节点和内部任意一点都与总体坐标系中单元的节点和内部点形成一一对应的关系，用数学表达式可写为任意四边形单元b)对应的矩形单元图3.7四边形等参元的原理3.2.6轴对称单元在一些实际问题中，如飞轮、转轴、活塞、汽缸套等都是回转结构，如果约束条件和载荷都对称于回转轴，其应力、应变和位移也都对称于回转轴线，这类应力应变问题称为轴对称问题。进行轴对称问题有限元分析时，通常采用柱坐标系来描述轴对称结构的应力、应变和位移。轴向、径向和周向坐标分别用z、r和θ表示，单元类型为实心圆环体。圆环体截面可以是三节点三角形，如图3.8a所示，也可以是四节点矩形，如图3.8b所示，单元之间用圆环形的铰相互连接，每个节点有两个自由度。此类单元的分析可看作是位移限制在rz平面的二维问题。当回转体为薄壳结构时，采用回转薄壁壳单元，如图3.9所示，该单元有两个节点，每个节点有三个自由度，两个位移和一个转角自由度，此单元的分析类似于平面梁单元。三角形环单元b)四边形环单元图3.8轴对称参元图3.9回转圆锥薄壳单元3.3单元分析单元分析是有限元分析中最为关键的一环。单元分析包括：

位移插值函数的选择，

单元应力应变分析，作用在单元上的载荷向单元节点的移置等3.3.1单元的插值函数由弹性力学理论可知，如果弹性体内的位移分量已知，则通过弹性力学几何方程和物理方程，应力分量和应变分量也都唯一的确定。但是如果仅知道每个单元节点的位移，既不能确定单元内任意一点的位移，更不能直接求得单元上任意一点的应力应变。只是在假定了一个可用节点位移表示单元内任意一点位移的插值函数以后，单元内任意一点的应力应变求解才成为可能。用一个位移函数表示整个结构的变形，不仅对于简单结构比较困难，而且对于复杂结构更不可能。怎么办呢？？？把结构划分成足够多的单元，对于每个单元来说，却是一件简单易行的事。每个单元的局部可以采用比较简单的函数来近似表达单元的真实位移。把各单元的位移连接起来，就可以近似表达整个区域的真实位移。在有限元分析中一般都采用多项式作为插值函数，多项式的项数由所选取的单元和单元的节点数决定，如对于平面三节点三角形单元有如下插值函数式中的上标e表示单元，而对于图3.10所示六节点三角形单元则有如下插值函数3.13.2由上式可以看到在选择多项式时，多项式系数的个数与每个单元的总自由度数是相等的。通常以形函数的方式表示插值函数。以平面矩形单元为例，单元如图3-3b所示。四节点矩形单元的插值多项式为在已知四个节点位移的前提下任意一点位移可表示为3.33.4式中的N就是形函数，其表达式为上式如果令，形函数可写如下成无量纲形式

3.53.6则(3.4)式就可以用矩阵表示为式中{δ}e

为节点位移列向量。

3.7形函数应有如下特征：(1)本节点上为1，其它节点上为0，即(2)在单元内任一点各形函数之和等于1，即这一性质反映单元的刚体位移。(3)单元任意一条边上的形函数，仅与该边两端节点的坐标有关，而与其他节点无关，这一性质可以保证相邻单元在公共边界上位移的连续性。

3.83.9为了保证有限元解答的正确性，下列条件：(1)位移模式必须包括单元的刚体位移和结构变形引起的位移(2)位移模式还必须包括单元的常应变状态和与坐标有关的应变。(3)位移模式在单元内连续，在单元的公共边界处协调，即单元的变形在单元之间不能开裂或重叠。在这三个条件中，前（1）、（2）条是必要条件，第（3）条是充分条件。3.3.2单元分析单元位移确定以后，就可以利用几何方程和物理方程求得单元的应力和应变。下面仍以平面四节点矩形单元为例推导单元刚度矩阵。由弹性力学平面问题的几何方程可知，单元内任意一点的应变为式中[B]称为应变矩阵，其分块子矩阵为3.113.12代入无量纲插值函数上式可写为式中，。再由弹性力学平面问题物理方程单元内任意一点的应力可表示为式中[S]为应力矩阵，[D]为弹性矩阵其表达式为3.123.13E是杨氏模量，μ是泊松比。应力矩阵[S]的分块子矩阵为3.143.15代入无量纲形函数上式变为对于平面应变问题，只需将上式中的和分别用和代替即可。利用虚功原理可求得单元节点力与单元节点位移之间的关系式令，上式简化为[K]e

称为单元刚度矩阵，t为单元厚度3.163.173.18单元刚度矩阵[K]e可分块表示为式中3.193.20对于平面应力问题每一个子块为r和p遍历i、j、l和m得到单元刚度矩阵。对于平面应变问题，只需将上式中的E和μ分别用和代替即可。单元刚度矩阵具有如下性质：

(1)单元刚度矩阵与所选单元的位移模式、几何形状、大小及单元的材料性质有关。不同位移模式、不同形状和大小的单元，单元刚度矩阵不同，计算精度也不同。

(2)单元刚度矩阵具有对称性。

(3)单元刚度矩阵是主元恒为正的奇异矩阵，即单元矩阵没有逆矩阵且Krr>

0。

3.3.3载荷移置根据有限元原理，所有作用在结构上的载荷必须等效的移置到节点上。因此，需将结构上的载荷按照静力等效的原则向节点移置，化为等效节点载荷。仍以平面问题四节点单元为例，假设作用在矩形单元上的体力为，分布面力是，集中力为。(1)集中力。等效载荷为由虚功原理得到：或3.223.21(2)体力。等效载荷为由虚功原理有或：3.233.243.25(3)分布面力。等效载荷为由虚功原理有其分量形式为s为单元上作用有外载荷的边。3.263.273.283.4整体分析仍以平面问题为例，单元为矩形四节点单元。假设弹性体被划分为N个单元，整个结构有n个节点。对于每一个单元都有一组如式(3.18)表示的方程组，将这些方程集合起来就得到表征整个结构弹性力学平衡的代数方程组。假定整个结构的节点位移列阵为{δ}2n×1，它是把节点位移按总体编码由小到大排列起来得到的，即式中节点i的位移分量表示为3.29同理节点载荷列阵为{F}2n×1，它是把每个移置到节点上的载荷相加并按总体编码由小到大排列起来得到的，即有

式中是节点i的等效节点载荷。3.30将单元刚度矩阵根据单元内节点的局部编码和节点整体编码的对应关系扩充为一个2n×2n的矩阵。对于每一个单元都可以根据编码关系得到以下一个矩阵3.31上式要遍历每一个单元，其中的i,j,l,m也是针对每一个单元的节点编号，显然由于每个单元中i,j,l,m的编号对应的整体节点编号不尽相同，因此由每个单元扩充得到的上述矩阵的i,j,l,m的位置也是不同的。把上式得到的N个单元进行叠加求和，就得到整体刚度矩阵，即其分块形式为3.333.32式中krp是的子矩阵，按下式计算上式是单元刚度矩阵扩充为2n×

2n阶之后，在同一位置上的子矩阵之和。至此就得到关于节点位移的所有2n×

2n个线性代数方程组，即式中[K]称为整体刚度矩阵。3.343.35整体刚度矩阵[K]具有以下性质：(1)整体刚度矩阵是对称的稀疏矩阵，矩阵中各个元素都集中分布于对角线附近，形成“带状”，其余元素均为零。(2)由于单元刚度矩阵对角线上的元素均大于零，由整体刚度形成的方式可知，整体刚度矩阵的主对角线元素必然大于零。(3)未经约束条件处理的刚度矩阵是奇异矩阵。故在求解有限元方程时，需要根据约束条件，修正结构刚度矩阵以消除奇异性。3.5边界条件处理整体刚度矩阵的奇异性是由于刚体位移的存在，代入边界条件可以消除刚体位移，从而消除刚体矩阵的奇异性，最终使方程可解。通常的位移边界条件有两种，一种是位移为零，另一种是指定位移。常用下列方法处理边界条件。

3.5.1划行划列法当某一位移为零时，如（r为整体编号），可将整体刚度矩阵中的第r行和第r列划去，同时划掉第r行载荷列阵和位移列阵元素，对于多个零位移条件依次处理。3.5.2对角线元素置1法当给定位移为零时，如（r为整体编号），在整体刚度矩阵除了让主对角元素外，整体刚度矩阵中的第r行和第r列元素均改为零，同时在整体载荷列阵中让

。这样修正以后，解方程时

。对于多个零位移边界条件依次处理，全部修正完毕再进行求解。这种方法不改变原来方程的阶数和节点未知量的顺序，便于计算机编程实现。3.5.3对角线元素乘大数法当节点位移为定值，如，类似对角线元素置1法，不同之处在于对角元素Krr

乘以大数a（a可取1010

左右数量级），并将Fr用aKrr

uc代替既可。解方程可得

，对于多个位移依次作上述修正。也可以看到对角线元素