有限元法讲义.doc

有限单元法

图3-1 简单的梁和桁架结构
现考虑对图3-1中结构的分析。在位移分析法中我们是把该结构看成是两个梁单元、一个桁架单元和一个弹簧单元的分割体，第一步是计算对应于结构总体自由度的单元刚度矩阵。在这种情况下，对于梁单元、弹簧单元和桁架单元，我们分别有
；
；
；
；
整个分割体的刚度矩阵可以由各个单元刚度矩阵通过直接刚度法有效地求得。在这个过程中，结构刚度矩阵K是通过各单元刚度矩阵直接相加而算得，即

而系统的平衡方程为

式中是系统的总体位移向量，R是作用在结构总体位移方向上的外力向量：
，
在求解结构的位移之前，我们需要利用边界条件和。这意味着我们可以只考虑含有五个未知位移的五个方程，即，式中是从K中删去第一和第七行以及第一和第七列后得到的，而
，
上面的讨论说明了有限元法的一些重要特点．基本的处理过程是先把整个结构看成为各个结构单元的分割体，计算对应于结构分割体总体自由度的各单元刚度矩阵，然后通过将各单元刚度矩阵叠加的方法形成结构刚度矩阵，求解单元分割体的平衡方程组就得到单元的位移，然后利用它们来计算单元的应力。虽然原来并不认为用位移法分析梁和衍架单元的分割体就是有限元分析，但以后我们将会看到，实际上我们可以把桁架和梁单元看作两种特殊的有限元。在这个定义上，上述的分析就是梁和衍架结构的有限元分析。而用另一种方法，我们可不通过求解平衡微分方程，而是用虚功原理计算其刚度系数。
导出表示弹性体平衡的相应方程的一个等效方法是利用虚位移原理，这个原理表示物体处于平衡的要求是：对于强加在该物体上的任意相容的微小的虚位移，总的内虚功应等于总的外虚功，即
（3.1）
式中
， ， （3.2）
是作用在弹性体上的体力、表面力和集中力。从未受荷载时的位置开始的弹性体的位移以表示，其中
（3.3）
相应于的应变为
（3.4）
相应的应力为
（3.5）
和表示虚应变和虚位移。
分析的基本步骤和上述桁架和梁结构的分析步骤一样。

图3-2 有限元平面分析
该问题的求解可按下列的步骤进行：
(1)假设每一单元节点i的两个未知位移为和，而和用简单多项式函数来表示，其中的末知参数是单元节点位移。对于图3-3中的单元，未知位移是。
(2)利用虚位移原理计算每个单元对应于节点自由度的刚度矩阵。
(3)将各单元刚度矩阵叠加得到结构的刚度矩阵，利用边界条件解平衡方程组求出节点位移，然后算出单元的应力。

图3-3 局部坐标系统中典型的二维四节点有限元

3.1 利用虚位移原理建立有限元法的公式
（一）平面应力分析的位移和应变－位移的变换矩阵
为了便于说明，再考虑一个平面内荷载作用下悬臂板分析的例子(图3-2)。该结构是处于平面应力的状态，因此可以将平桓理想化为二维平面应力有限元的分割体，如图3-2所示。
所需要的基本矩阵是对于分割体的每个单元的位移变换矩阵和应变-位移变换矩阵。为了推导这些矩阵，我们着重研究如图3-3所示的典型单元，并假设局部单元位移和是以局部坐标变量和的多项式的形式给出：
（3.6）
（3.7）
未知系数也称为广义坐标，它们可以用未知的单元节点位移和来表示。定义
（3.8）
我们可以把式（3.6）和（3.7）写成矩阵形式
（3.9）
其中
（3.10）
（3.11）
（3.12）
对于单元的各个节点，方程(3.9)都一定是成立的。因此，对于所有四个节点，利用式(3.9)，我们有
（3.13）
其中
（3.14）
因而广义坐标是
（3.15）
其中
（3.16）
将式（3.9）代入，可以得到
（3.17）
其中
（3.18）
我们用虚功原理来得到单元刚度矩阵。若是平面应力状态，单元应变为
（3.19）
其中
（3.20）
利用式（3.17），我们得
（3.21）
式中
（3.22）
（3.23）
单元应力是
（3.24）
假定是各向同性的线弹性材料，对于平面应力状态，有
（3.25）
其中
（3.26）
而和是材料的弹性常数，为杨氏模量，为泊松比。
要进行实际的有限元分析，就需要求出式(3.18)和式(3.22)中分别给出的分割体中每个单元的位移和应变-位移的变换矩阵，而算出每个单元的矩阵A，才能求出这些矩阵。
至今，我们所考虑的是通过单元局部节点位移来决定的每个有限元的位移。然而，在推导整个单元分割体的平衡方程的时候，通过整个单元分割体节点位移来表示单元的位移是方便的，即对单元m可写出
（3.27）
而对于单元应变和应力类似地可写成
（3.28）
（3.29）
式(3.27)至式(3.29)中的量与向量有关，而存储整个有限元分割体的总体坐标系下的全部节点位移。作为一个例子，考虑，即确定图3.2中悬臂板理想化后的单元2上位移矩阵。利用图3.2给出的整个单元分割体和图3.3给出的单元的节点位移的定义，我们有

（3.30）
式中是式（3.17）中考虑单元2时对应的单元。
3.2 建立有限元步骤归纳
有限元法与结构力学中的位移法相似。首先将连续体转化为离散化结构，即将连续体代之以仅在节点互相连结的许多单元组成的结构。这种离散化结构类似于结构力学中的桁架或刚架，但其中的单元不一定是杆件，可以是平面块体或空间块体等。然后对此离散化结构按类似结构力学的位移法进行分析，步骤如下：
1． 取每个单元的节点位移作为基本未知数。
2． 在单元内建立位移模式，。
3． 根据几何关系建立应变矩阵，。
4． 根据物理方程建立应力矩阵，。其中称为应力矩阵。
5． 根据虚位移原理建立单元刚度矩阵，。
6． 建立单元节点荷载向量，。
7． 建立平衡方程，，其中，，是全部节点位移向量。
8． 代入约束条件求解平衡方程，。
位移模式中的N称为形态函数。在有限单元法中，各种计算公式都依赖于位移模式，位移模式选择的恰当与否，与有限单元法的计算精度和收敛性有关。为了使位移模式尽可能地反映物体中的真实位移形态，它应满足下列条件：
(1)位移模式必须能反映单元的刚体位移。
(2)位移模式必须能反映单元的常量应变。
(3)位移模式应尽可能地反映位移的连续性。
在单元之间，除了节点处有共同的节点位移值外，还应尽可能反映在单元之间边界上位移的连续性。
第一章、 杆件结构的有限单元法
4.1 平面刚架结构
平面刚架结构的有限单元法分析，通常采用两端固接的平面固结单元(见图4-1)，现针对这种单元讨论有限元分析计算公式的建立。

图4-1平面刚架单元
（a）坐标系，单元节点位移与节点力向量
杆件结构的有限元分析需要建立单元局部坐标系与结构整体坐标系。对平面固结单元(如图4-1所示)，结构矩阵位移法取结构的节点位移为基本未知量。单元节点位移与基本未知量之间满足变形协调条件，在局部坐标和整体坐标下单元节点位移、节点力分别记为
（4.1）
（4.2）
（4.3）
（4.4）
图4-1中各量值所示方向均为正方向。设单元的弹性模量为E，惯性矩为I，截面积为A，杆长为L。
（b）位移模式，应力与应变矩阵
由材料力学知识可知，位移模式可以选择
（4.5）
（4.6）
注意，可得单元的位移模式
（4.7）
其中，形函数

单元的线应变分为拉压应变和弯曲应变两部分(略去剪切变形)，于是有
（4.8）
其中

对应的应力为
（4.9）
（c）单元刚度矩阵
根据虚位移原理可知，在局部坐标系中
（4.10）
具体表达式为

（d）坐标转换矩阵，整体坐标下单元刚度矩阵
从图4-1可知，局部坐标系相对于整体坐标系逆时针旋转了θ角。根据单元的节点位移或节点力在两个坐标系中的投影关系，可得坐标转换关系
（4.11）
于是，单元节点位移或节点力的坐标转换矩阵(从整体坐标转到局部坐标)为
（4.12）
根据正交矩阵的特性和矩阵的相关知识，可知在整体坐标中的单元刚度矩阵为
（4.13）
（e）节点荷载向量
单元上的荷载按虚位移原理移置到节点上、得到等效节点荷载向量
（4.14）
其中是作用在杆上的均布荷载，是集中荷载。
整体坐标中的节点荷载同局部坐标中的关系为
（4.15）
（f）结构的平衡方程
根据虚位移原理，结构的平衡方程表示为
（4.16）
其中，，分别代表整体总刚度矩阵，总节点位移向量和总节点力向量。
4.2 平面桁架结构
平面桁架结构的有限元分析通常采用两端铰接的平面铰结单元，其坐标系如图4-2所示。
（a）坐标系，单元节点位移与节点力向量
杆件结构的有限元分析需要建立单元局部坐标系与结构整体坐标系。对平面固结单元(如图4-2所示)，在局部坐标和整体坐标下单元节点位移、节点力分别记为
（4.17）
（4.18）

图4-2 平面桁架单元

（4.18）
（4.19）
（b）局部坐标中的单元刚度矩阵
（4.20）
（c）坐标转换矩阵
（4.21）
（d）整体坐标中的单元刚度矩阵
（4.22）
其中，。
第二章、 平面问题的有限单元法
连续体的离散化与三角形常应变单元
（a）三角形单元的节点编码

图5-1 平面应力问题

我们用一组网格线把平板划分成若干三角形区域，每个三角形就算作是一个单元，网格线的交点就算是节点。把这些单元和节点按一定的顺序编号，并把各个节点沿坐标轴的位移取为基本未知量：
对于节点1
对于节点2
对于节点i
设有n个节点，则整体的总节点位移向量为

对每个单元来说，它与三个节点发生联系，设其整体序号是，依次标记为，则此单元的节点位移向量可以表达为（参看图5-2）

图5-2 三角形常应变单元
单元的节点位移向量是整体的总节点位移向量的一部分。

（b）位移模式
和桁架不同，这里没有现成的用节点位移向量来表示单元内任一点位移的现成公式，因而必须采用近似假设。我们知道，任一函数常常可以用有限项的幂级数来迫近。所以不妨假设单元上任何一点的位移为
（5.1）
也就是把高阶的项全都略去，只保留线性项，这自然会带来截尾误差。从几何意义上来解释：函数可认为是定义在单元上的某种复杂曲面，而现在用一平面来代表它，这样做当然是有误差的。可是我们从直觉上可以想像得出来，若把单元划分得很小，那末在足够小的范围内平面和曲面是相差不大的。所以可以用式（5.1）来代表真正的位移函数而不致造成很大的误差。
现在要用节点位移来表示式（5.1）中的系数，将节点坐标代入该式，位移应为节点位移，即
（5.2）
由上式可以得到用节点位移表示待定系数的式子，之后再代回到式（5.1）之第一式，就有
（5.3）
其中
（5.4）
，， （5.5）
为单元的面积。
其他的式子可用下述方法得到：将上边四个式子小的下标按“循环交换”的规则依次更改。所谓循环交换就是的交换方式。
同理可得
（5.6）
综合起来，就是
（5.7）
其中
（5.8）
通过以上处理方法，连续体上任一点处的位移皆可用若干选定的节点位移来表示。一个寻求未知函数的问题已经转化成一个寻求若干个未知数的问题。这样，我们就完成了将连续体离散化的关健一步。由此可见形状函数在有限单元法中占据的重要地位。
式（5.4）和式（5.6）表明，这种单元的形状函数矩阵是由三个节点处的形状函数组合而成的。每个形状函数都是线性函数，在几何上代表一个平面。在所代表的节点处，函数的值为单位1，在另外两个节点处，函数的值为零。图5-3是这类函数的图形表示(以为例)。

图5-3 形状函数平面

（c）应力应变矩阵
利用前面的知识，可知应变
（5.9）
其中应变矩阵
（5.10）
； （5.11）
可见由于采用了线性的形状函数，单元上的应变(以及应力)为常数。所以这种单元也可称之为常应变(或常应力)单元。这当然是近似的，这种应力实际上只是真实应力的某种平均值。对于平面应力问题，应力为
（5.12）

其中为材料的泊松比，为材料的弹性模量。对于平面应变问题，只要将弹性矩阵中的用、用代替即可。将应变代入，可得
（5.13）
其中称为应力矩阵（针对平面应力情形）
（5.14）
； （5.15）

（d）单元刚度矩阵
根据虚位移原理，由式（3.39）可知，单元的刚度矩阵为

当平面问题的板厚度为时，上式改写为

此处的标记均表示第个单元。将式（5.9）到式（5.15）代入上式，可得
（5.16）
式中每个子矩阵都是阶的，其下标与单元的节点编号相对应，具体的计算公式如下：
；

（e）节点荷载向量
对于连续体来说，在单元的边界上有着分布作用的边界力，它们是其他单元或是结构的外部环境对该单元的作用力。而在节点上，一般却没有集中的节点力，这和杆系结构不一样。不过出于习慨在物理意义上我们常常可以设想把这些边界力向节点处简化，形成一组等价的虚拟的节点力。显然当单元尺寸很小时用这种等价节点力代表真实的边界力不致引起很大的误差。
由第三章的知识，我们知道三节点三角形单元的节点荷载向量表示为

下面给出几种情况下的节点荷载向量的表达式：
（1） 设在三角形单元的点受到集中荷载作用，则
（5.17a）
（2） 设单元受方向的体力作用，例如受重力作用。其合力，为容重。
（5.17b）
（3） 设单元的边上受有沿方向的集中力，其作用点距离节点的距离分别为，，，则
（5.17c）
（4） 设单元的边上受有沿方向的三角形分布荷载，荷载在1点的集度为，在2点的集度为零，边长为，则
（5.17d）
（5） 设单元的边上受有均匀分布的侧压力，则
（5.17e）
5.1 面积坐标
在采用较精密的三角形单元时，利用面积坐标，可以大大简化荷载向量、应力矩阵和刚度矩阵矩阵的推导工作。在如图5-4所示的三角形单元中，任意一点P的位置，可以用如下的三个比值来确定

图5-4 面积坐标
，， （5.18）
其中A为三角形的面积，、、分别为三角形、、的面积。这三个比值就称为P点的面积坐标。注意，三个面积坐标并不是互相独立的，即
（5.19）
这里所引用的面积坐标，只限于用在一个三角形单元之内，在该三角形之外并没有定义，因而是一种所谓局部坐标。与此相反，以前所用的直角坐标和，则是一种所谓整体坐标，它通用于所有单元的总体也就是通用于整个结构物。
根据面积坐标的定义，在图5.4中不难看出，在平行于边的一根直线上的所有各点，都具有相同的坐标而且这个坐标就等于“该直线至边的距离”与“节点i至边的距离”的比值。图中示出的一些等值线。
根据面积的计算公式，可知
， （5.20）
同式（5.4）相比，可知三角形单元的形状函数，就是对应的。用矩阵表示为
（5.21）
将上三式分别乘以然后相加，可得

同样，用坐标，可得

上二式表明了直角坐标与面积坐标之间的关系，写成矩阵形式
（5.22）
对面积坐标函数求导数时，可以利用下面的公式
（5.23）
求面积坐标的幂函数在三角形单元上的积分值时，可应用积分公式
（5.24）
求面积坐标的幂函数在三角形莱一边上的积分值时，可应用积分公式
（5.25）
其中为边的长度。
（a）位移模式
假定位移模式按照二阶多项式变化，即
（5.26）

图5-5 六节点三角形单元

上式中包含着十二个待定系数，为了把它们用节点位移表示出来，每个单元应当有六个节点。通常的作法是在三个顶点之外再取三个边的个点，如图5-5所示。对应的形状函数矩阵是
（5.27）
其中，并且
（5.28）
式中六个系数与三角形的几何尺寸有关

式中的下标仍然按照循环。

（b）应变矩阵
对应的应变矩阵为
（5.29）
式中，可见应变是线性变化的，并可以表示为。

（c）单元刚度矩阵
同理，单元刚度矩阵表达为
（5.30）
式中任意一个子矩阵可以写成

其中

此外，计算时应注意到，分别为三角形单元的体积、静矩、惯性矩和惯性积等等。
5.2 矩形单元
图5-8 矩形单元
三角形单元的优点之一是它的“适应性”， 任何复杂边界的弹性体，总是可以划分成三角形。可是在规则边界的情况下，显然划分成矩形更加方便。另外，许多计算实例已经证明，矩形单元的计算精度也比三角形单元好。这是因为在每个小矩形范围内，矩形单元有连续变化的应力场，而对应的两个三角形单元却只有阶梯变化的应力场。所以，在复杂边界的情况下，同时使用三角形和矩形单元将是可取的计算方案。
与三角形单元不问，不可能采用完全的多项式作为位移函数。因为一次的完全多项式只有3个待定系数而矩形单元却有4个角点。为使二者的数目一致应在二次的多项式中选取补充项。设我们补充的项，则在矩形的上、下边界上位移曲线是二次的，它不可能由两角点处的位移来唯一地确定，因而是不相容的。同样的道理，也不可以取的项。所以，应取
（5.41）
这个函数虽然含有二次项，但在单元的任何一条边界上都是线性的，所以只要两点就可以唯一的确定位移，因而满足相容条件。
根据上式可以写出：
；；；
若改用自然坐标表示（见图5-8），，可得
；
；
或者概括成
； （5.42）
单元的应变为

式中为单元节点位移向量。

； （5.43）
对于平面应力问题，单元的刚度矩阵为

式中

；（5.44）
其中，h为单元厚度。对于平面应变问题按照前述同样处理。
从式（5.43）可以看出，在单元内部，应力（应变）不是常数。它虽然不是完全的一次多项式表达的，但却是一种线性变化的规律。
6.1 等参单元
等参单元是在1966年公开引入的（由Irons），它们使得非矩形四边形单元的存在成为可能。等参单元最明显的特点是单元可以有曲边，并且它们有特殊的坐标系（称为自然坐标系）。等参单元在模拟具有曲线边界的结构和从粗网格到细网格作的单元过渡时是很有用的。等参单元具有多方面的适用性；在二维和三维弹性分析、壳体分析以及非结构的应用中已被证明是有效的。等参单元的构造过程初看起来似乎生疏，但并不复杂。一旦掌握了一种类型单元的构造过程，就很容易应用到大多数其它类型的等参单元上。
单元节点确定两件事情
1． 节点位移向量（亦称节点自由度）表示单元内部的位移，具有；
2． 节点坐标向量定义了单元内部任意一点的坐标（用向量表示），具有。
矩阵和都是自然坐标的函数，如果上述两项中的节点集相同，且和也相同，则这个单元是等参的。
6.2 平面线性等参单元
首先，需要再次说明，不同的等参单元终究是类似的。本节的列式可直接推广到更复杂的单元：主要的变化是增加节点和利用不同的形函致。
矩形单元，最大的不足在于其对边界的要求高，即只能是规则的、平行于整体坐标轴的。如果能有任意四边形的单元，这样的问题就解决了，但任意的四边形单元的形函数的构造却不容易，等参单元能够很好的解决这样的困难。

图6-2 四边形等参单元映射关系

图6-2中的坐标轴和通过四边形对边中点，它们不需要正交，也不必平行于x或y轴。坐标的方位由节点号确定。就是说，对于如下定义的形函数N，节点1在处，节点2在处等等。
整体坐标和位移为
； （6.9）
其中
（6.10a）
（6.10b）
（6.10c）
有时将式（6.9）写成如下形式更为方便
和 （6.11）
其中，各个形函数为
； （6.12）
同矩形单元的完全一样。
通常，u和v分别是平行于x和y方向的位移，在一般情形下它们不平行于和。一个特殊情形是其边平行于整体坐标xy的矩形，此时和变成无量纲的形心坐标。
下述推导类似上节的推导。为了求得单元刚度矩阵K，我们需要B矩阵，我们不能用x和y来表示B矩阵，而用和来写出它，这就要求引用下列导数的坐标变换。设是某一x和y的函数。那么由复合求导法则得到
或者 （6.13）
其中J为雅可比矩阵
（6.14）
从式（6.13）可以得到逆关系
（6.15）
上述关系是一般的情况。对于本章的平面单元，就是u或v，雅可比阵J的系数的数值与单元的尺寸、形状和单元所处的方位（节点坐标和节点位置编号）有关。对于目前讨论的4节点单元。从式（6.11）有
（6.16）
其它的具有类似的表达式。如果和，则。
现在，可以计算B矩阵了。
（6.17）
（6.18）
（6.19）
B矩阵就是式（6.17）、（6.18）和（6.19）中三个矩阵的乘积。单元刚度矩阵为
（6.20）
其中是雅可比矩阵的行列式值，h为单元厚度。在平面问题中
（6.21）
雅可比行列式是单元内位置的函数，并且等于从面积变到的放大倍数。矩阵B中的一个典型系数是与节点坐标有关的，并且在分子和分母上都是变量的多项式。因此，式（6.20）中的积分必须以数值方法求出。对于单元荷载向量中的体力和初始应力部分（见式（3.36）的部分）都需要体积积分也必须采用数值积分（见下节）
（6.22）
我们可以根据伴随矩阵求得式（6.18）中雅可比矩阵的逆矩阵
（6.23）
式（6.20）和（6.22）中的积分微元的变换，是来源于如下关系。在直角坐标中

图6-3 坐标变换关系
由于为非正交坐标系（见图6-3），所以
（6.24）
其中

代入到式（6.24），可得
（6.25）
6.3 8节点平面等参单元
将单元的边节点加到线性单元上，可以得到等参族较高次的单元。图6-4所示8节点单元，边缘的形状和位移按照或变化。
在二次单元中的形状和位移的插值多项式是以不完全三次多项式为基础的，但是，在这里我们不打算用广义位移替换节点位移的方法建立形函数N，而是用一种更直接、更有明显的物理意义的方法来构造形函数。我们来考察各种模式并把它们组合起来，在图6-5a中，节点位移按照二次方、按照线性变化。u除掉在节点5之处为1外，在其它节点处皆为零，因而形函数，类似的，形函数。在图6-5c中，位移u按和线性变化，而在节点5和8分别取。

图6-4 8节点平面等参单元

图6-5 二次单元位移模式

这样在节点1将节点5和8的位移的一般减去，便可以得到节点1的形函数 （6.29）
其中

表6-2 平面二次等参单元的形函数
4个线性边
3，2，1或0个线性边，其余为二次边
包括1～4节点
加节点5
加节点6
加节点7
加节点8

28