3 弹性力学平面问题的有限元法 本章包括以下的内容： 3.1弹性力学平面问题的基本方程 3.2单元位移函数 3.3单元载荷移置 3.4单元刚度矩阵 3.5单元刚度矩阵的性质与物理意义 3.6整体分析 3.7约束条件的处理 3.8整体刚度矩阵的特点与存储方法 3.9方程组解法 3.1弹性力学平面问题的基本方程 弹性力学是研究弹性体在约束和外载荷作用下应力和变形分布规律的一门学科。在弹性力学中针对微小的单元体建立基本方程，把复杂形状弹性体的受力和变形分析问题归结为偏微分方程组的边值问题。弹性力学的基本方程包括平衡方程、几何方程、物理方程。 弹性力学的基本假定如下： 1）完全弹性，2）连续，3）均匀，4）各向同性，5）小变形。 3.1.1基本变量 弹性力学中的基本变量为体力、面力、应力、位移、应变，各自的定义如下。 体力 体力是分布在物体体积内的力，例如重力和惯性力。 面力 面力是分布在物体表面上的力，例如接触压力、流体压力。 应力 物体受到约束和外力作用，其内部将产生内力。物体内某一点的内力就是应力。  图3.1 如图3.1假想用通过物体内任意一点p的一个截面mn将物理分为Ⅰ、Ⅱ两部分。将部分Ⅱ撇开，根据力的平衡原则，部分Ⅱ将在截面mn上作用一定的内力。在mn截面上取包含p点的微小面积，作用于面积上的内力为。 令无限减小而趋于p点时，的极限S就是物体在p点的应力。  应力S在其作用截面上的法向分量称为正应力，用σ表示；在作用截面上的切向分量称为剪应力，用τ表示。 显然，点p在不同截面上的应力是不同的。为分析点p的应力状态，即通过p点的各个截面上的应力的大小和方向，在p点取出的一个平行六面体，六面体的各楞边平行于坐标轴。  图3.2 将每个上的应力分解为一个正应力和两个剪应力，分别与三个坐标轴平行。用六面体表面的应力分量来表示p点的应力状态。应力分量的下标约定如下： 第一个下标表示应力的作用面，第二个下标表示应力的作用方向。 ，第一个下标x表示剪应力作用在垂直于X轴的面上，第二个下标y表示剪应力指向Y轴方向。 正应力由于作用表面与作用方向垂直，用一个下标。表示正应力作用于垂直于X轴的面上，指向X轴方向。 应力分量的方向定义如下： 如果某截面上的外法线是沿坐标轴的正方向，这个截面上的应力分量以沿坐标轴正方向为正； 如果某截面上的外法线是沿坐标轴的负方向，这个截面上的应力分量以沿坐标轴负方向为正。 剪应力互等： 物体内任意一点的应力状态可以用六个独立的应力分量、、、、、来表示。 位移 位移就是位置的移动。物体内任意一点的位移，用位移在x，y，z坐标轴上的投影u、v、w表示。 应变 物体的形状改变可以归结为长度和角度的改变。 各线段的单位长度的伸缩，称为正应变，用ε表示。 两个垂直线段之间的直角的改变，用弧度表示，称为剪应变，用γ表示。 物体内任意一点的变形，可以用六个应变分量表示。 3.1.2平面应力和平面应变问题 弹性体在满足一定条件时，其变形和应力的分布规律可以用在某一平面内的变形和应力的分布规律来代替，这类问题称为平面问题。平面问题分为平面应力问题和平面应变问题。 平面应力问题 设有很薄的等厚薄板，只在板边上受到平行于板面并且不沿厚度变化的面力，体力也平行于板面且不沿厚度变化。  图3.3 设板的厚度为t，在板面上： ， ，  由于平板很薄，外力不沿厚度变化，因此在整块板上有， ， ，  剩下平行于XY平面的三个应力分量未知。 2）平面应变问题 设有很长的柱形体，支承情况不沿长度变化，在柱面上受到平行于横截面而且不沿长度变化的面力，体力也如此分布。  图3.4 以柱体的任一横截面为XY平面，任一纵线为Z轴。假定该柱体为无限长，则任一截面都可以看作对称面。由对称性， ，， 由于没有Z方向的位移，Z方向的应变。 未知量为平行于XY平面的三个应力分量，物体在Z方向处于自平衡状态。 3.1.3平衡方程 弹性力学中，在物体中取出一个微小单元体建立平衡方程。平衡方程代表了力的平衡关系，建立了应力分量和体力分量之间的关系。对于平面问题，在物体内的任意一点有，  （3-1） 3.1.4几何方程 由几何方程可以得到位移和变形之间的关系。对于平面问题，在物体内的任意一点有，  （3-2） 刚体位移 由位移u=0，v=0可以得到应变分量为零，反过来，应变分量为零则位移分量不为零。应变分量为零时的位移称为刚体位移。刚体位移代表了物体在平面内的移动和转动。 由  可以得到刚体位移为以下形式，  由可得，  将代入可得，  积分后得到，  得到位移分量，  当时，物体内任意一点都沿x方向移动相同的距离，可见代表物体在x方向上的刚体平移。 当时，物体内任意一点都沿y方向移动相同的距离，可见代表物体在y方向上的刚体平移。 当时，可以假定，此时的物体内任意一点P（x，y）的位移分量为，  P点位移与y轴的夹角为α，   P点合成位移为，  r为P点到原点的距离，可见ω代表物体绕z轴的刚体转动。 3.1.5物理方程 弹性力学平面问题的物理方程由广义虎克定律得到。 1）平面应力问题的物理方程   （3-3）  平面应力问题有，   2）平面应变问题的物理方程   （3-4）  平面应变问题有，   在平面应力问题的物理方程中，将E替换为、替换为，可以得到平面应变问题的物理方程；在平面应变问题的物理方程中，将E替换为、替换为，可以得到平面应力问题的物理方程。  图3.5 求解弹性力学平面问题，可以归结为在任意形状的平面区域内已知控制方程、在位移边界上约束已知、在应力边界上受力条件已知的边值问题。然后以应力分量为基本未知量求解，或以位移作为基本未知量求解。 如果以位移作为未知量求解，求出位移后，由几何方程可以计算出应变分量，得到物体的变形情况；再由物理方程计算出应力分量，得到物体的内力分布，就完成了对弹性力学平面问题的分析。 3.2单元位移函数 根据有限元法的基本思路，将弹性体离散成有限个单元体的组合，以结点的位移作为未知量。弹性体内实际的位移分布可以用单元内的位移分布函数来分块近似地表示。在单元内的位移变化可以假定一个函数来表示，这个函数称为单元位移函数、或单元位移模式。 对于弹性力学平面问题，单元位移函数可以用多项式表示，   （3-5） 多项式中包含的项数越多，就越接近实际的位移分布，越精确。具体取多项，由单元形式来确定。即以结点位移来确定位移函数中的待定系数。  图3.6 如图3.6所示的3结点三角形单元，结点I、J、M的坐标分别为、、，结点位移分别为、、、、、。六个节点位移只能确定六个多项式的系数，所以3结点三角形单元的位移函数如下，  （3-6） 将3个结点上的坐标和位移分量代入公式（3-6）就可以将六个待定系数用结点坐标和位移分量表示出来。 将水平位移分量和结点坐标代入（3-6）中的第一式，  写成矩阵形式，  （3-7） 令 ， 则有  （3-8）  ，A为三角形单元的面积。 [T]的伴随矩阵为，  （3-9） 令  （3-10） 则  （3-11） 同样，将垂直位移分量与结点坐标代入公式（3-6）中的第二式，可得，  （3-12） 将（3-11）、（3-12）代回（3-6）整理后可得，   令 （下标i，j，m轮换） 可得  （3-13） 单元内的位移记为  单元的结点位移记为  单元内的位移函数可以简写成，  （3-14） 把[N]称为形态矩阵，Ni称为形态函数。 选择单元位移函数应满足以下条件： 反映单元的刚体位移与常量应变。 2）相邻单元在公共边界上的位移连续，即单元之间不能重叠，也不能脱离。 由（3-6）可以将单元位移表示成以下的形式，   反映了刚体位移和常应变。 单元位移函数是线性插值函数，因此单元边界上各点的位移可以由两个结点的位移完全确定。两个单元的边界共用两个结点，所以边界上的位移连续。 形态函数Ni具有以下性质： 在单元结点上形态函数的值为1或为0。 2）在单元中的任意一点上，三个形态函数之和等于1。 用来计算三角形面积时，要注意单元结点的排列顺序，当三个结点i，j，m取逆时针顺序时，；当三个结点i，j，m取顺时针顺序时，。 例题：如图3.7所示等腰三角形单元，求其形态矩阵[N]。  解： 由    在公式中轮换下标可以计算得 ，，  ，，  ，  三角形积为 形态函数为    形态矩阵为  三角形面积的计算公式可得，  如果把三个结点按顺时针方向排列，即i（a，0），j（0，0），m（0，a）  3.3单元载荷移置 有限元法的求解对象是单元的组合体，因此作用在弹性体上的外力，需要移置到相应的结点上成为结点载荷。载荷移置要满足静力等效原则。静力等效是指原载荷与结点载荷在任意虚位移上做的虚功相等。 单元的虚位移可以用结点的虚位移表示为，  （3-15） 令结点载荷为 1）集中力的移置 如图3.7所示，在单元内任意一点作用集中力  图3.8 由虚功相等可得，  由于虚位移是任意的，则  （3-16） 例题1：在均质、等厚的三角形单元ijm的任意一点p（xp，yp）上作用有集中载荷。  体力的移置 令单元所受的均匀分布体力为 由虚功相等可得，   （3-17） 分布面力的移置 设在单元的边上分布有面力，同样可以得到结点载荷，  （3-18） 例题2：设有均质、等厚的三角形单元ijm，受到沿y方向的重力载荷qy的作用。求均布体力移置到各结点的载荷。      同理， 例题3：在均质、等厚的三角形单元ijm的ij边上作用有沿x方向按三角形分布的载荷，求移置后的结点载荷。  取局部坐标s，在i点s=0，在j点s=l，L为ij边的长度。在ij边上，以局部坐标表示的插值函数为， ，， 载荷为   3.4单元刚度矩阵 根据单元的位移函数，  由几何方程可以得到单元的应变表达式，  （3-19） 记为 ，[B]矩阵称为几何矩阵。 [B]矩阵可以表示为分块矩阵的形式  （3-20） 由物理方程，可以得到单元的应力表达式，  （3-21） [D]称为弹性矩阵，对于平面应力问题，  定义为应力矩阵。 将应力矩阵分块表示为，   （3-22） 应用虚功原理可以建立单元结点位移与结点力的关系矩阵，单元刚度矩阵。 虚功原理：在外力作用下处于平衡状态的弹性体，如果发生了虚位移，则所有外力在虚位移上做的虚功等于内应力在虚应变上做的虚功。 单元的结点力记为 单元的虚应变为 单元的外力虚功为，  单元的内力虚功为，  由虚功原理可得，  （3-23）     （3-24） 定义为单元刚度矩阵。 在3结点等厚三角形单元中[B]和[D]的分量均为常量，则单元刚度矩阵可以表示为，  （3-25） 单元刚度矩阵表示为分块矩阵：   （3-26）  （3-27） 3.5单元刚度矩阵的性质与物理意义 （一）单元刚度矩阵的物理意义 假设单元的结点位移如下： 由，得到结点力如下：  （3-28） 表示i结点在水平方向产生单位位移时，在结点i的水平方向上需要施加的结点力。 表示i结点在水平方向产生单位位移时，在结点i的垂直方向上需要施加的结点力。 选择不同的单元结点位移，可以得到单元刚度矩阵中每个元素的物理含义： 表示s结点在水平方向产生单位位移时，在结点r的水平方向上需要施加的结点力。 表示s结点在水平方向产生单位位移时，在结点r的垂直方向上需要施加的结点力。 表示s结点在垂直方向产生单位位移时，在结点r的水平方向上需要施加的结点力。 表示s结点在垂直方向产生单位位移时，在结点r的垂直方向上需要施加的结点力。 因此单元刚度矩阵中每个元素都可以理解为刚度系数，即在结点产生单位位移时需要施加的力。 （二）单元刚度矩阵的性质 对称性 利用分块矩阵的性质证明如下：    即 奇异性 即单元刚度矩阵的行列式为零，。 将定单元产生了x方向的刚体移动，，此时对应的单元结点力为零。  可以得到，在单元刚度矩阵中1，3，5列中对应行的系数相加为零，由行列式的性质可知，。 同样如果假定单元产生了y方向上的刚体位移，可以得到，在单元刚度矩阵中2，4，6列中对应行的系数相加为零。 3.6整体分析 得到了单元刚度矩阵后，要将单元组成一个整体结构，根据结点载荷平衡的原则进行分析，即整体分析。 整体分析包括以下4个步骤： 建立整体刚度矩阵， 根据支承条件修改整体刚度矩阵， 解方程组，求出结点的位移， 根据结点位移，求出单元的应变和应力。 在这里把结点位移作为基本未知量求解。 如何得到整体刚度矩阵？基本方法是刚度集成法，即整体刚度矩阵是单元刚度矩阵的集成。  图3.9 如图3.9所示，一个划分为6个结点、4个单元的结构。得到了每个单元的单元刚度矩阵后，要集成为整体刚度矩阵。 3.6.1刚度集成法的物理意义 由单元刚度矩阵的物理意义可知，单元刚度矩阵的系数是由单元结点产生单位位移时引起的单元结点力。 在如图3.9所示的结构中，使结点3产生单位位移时，在单元（1）中的结点2上引起结点力。由于结点2、3同时属于单元（1）、（3），在单元（2）中的结点2上同样也引起结点力，因此，在整体结构中当结点3产生位移时，结点2上的结点力应该是单元（1）、（2）在结点2上的结点力的叠加。 刚体集成法即结构中的结点力是相关单元结点力的叠加，整体刚度矩阵的系数是相关单元的单元刚度矩阵系数的集成。结点3在整体刚度矩阵的对应系数，应该是单元（1）、（3）、（4）中对应系数的集成。 3.6.2刚度矩阵集成的规则 将单元刚度矩阵中的每个分块放到在整体刚度矩阵中的对应位置上，得到单元的扩大刚度矩阵。 单元刚度矩阵系数取决于单元结点的局部编号顺序，必须知道单元结点的局部编号与该结点在整体结构中的总体编号之间的关系，才能得到单元刚度矩阵中的每个分块在整体刚度矩阵中的位置。将单元刚度矩阵中的每个分块按总体编码顺序重新排列后，可以得到单元的扩大矩阵。 假定单元结点的局部编号与整体的对应关系如下： 单元编号 单元结点局部编号 单元结点整体编号  1 i 3  1 j 1  1 m 2  2 i 5  2 j 2  2 m 4  3 i 5  3 j 3  3 m 2  4 i 3  4 j 5  4 m 6   单元（2）的单元扩大矩阵的分块矩阵形式如下，只列出非零的分块： 局部 编号          整体 编号 1 2 3 4 5 6   1         2         3         4         5         6         2）将全部单元的扩大矩阵相加得到整体刚度矩阵。  整体刚度矩阵如下所示： 整体 编号 1 2 3 4 5 6  1        2   + +  +   +   3     + +   +   4        5   +  +   + +   6         3.7约束条件的处理 图3.9所示的结构的约束和载荷情况，如图3.10所示。结点1、4上有水平方向的位移约束，结点4、6上有垂直方向的约束，结点3上作用有集中力（Px，Py）。  图3.10 整体刚度矩阵[K]求出后，结构上的结点力可以表示为：  根据力的平衡，结点上的结点力与结点载荷或约束反力平衡。用表示结点载荷和支杆反力，则可以得到结点的平衡方程：  （3-29） 这样构成的结点平衡方程组，在右边向量{P}中存在未知量，因此在求解平衡方程之前，要根据结点的位移约束情况修改方程(3-29)。先考虑结点n有水平方向位移约束，与n结点水平方向对应的平衡方程为：  （3-30） 根据支承情况，方程（3-30）应该换成下面的方程：  （3-31） 对比公式（3-30）和（3-31），在式（3-29）中应该做如下修正： 在[K]矩阵中，第2n-1行的对角线元素改为1，该行中全部非对角线元素改为0；在{P}中，第2n-1个元素改为0。为了保持[K]矩阵的对称性，将第2n-1列的全部非对角元素也改为0。 同理，如果结点n在垂直方向有位移约束，则（3-29）中的第2n个方程修改为，  在[K]矩阵中，第2n行的对角线元素改为1，该行中全部非对角线元素改为0；在{P}中，第2n个元素改为0。为了保持[K]矩阵的对称性，将第2n列的全部非对角元素也改为0。  （3-32） 对图3.9所示结构的整体刚度在修改后可以得到以下的形式，  （3-33） 如果结点n处存在一个已知非零的水平方向位移，这时的约束条件为，  （3-34） 在[K]矩阵中，第2n-1行的对角线元素乘上一个大数A，向量{P}中的对应换成，其余的系数保持不变。 方程改为，  （3-35） A的取值要足够大，例如取1010。只有这样，方程（3-35）才能与方程（3-34）等价。 如果结点n处存在一个已知非零的垂直方向位移，这时的约束条件为， 。 也可以采用同样的方法修改整体刚度矩阵。 3.8整体刚度矩阵的特点与存储方法 用有限元方法分析复杂工程问题时，结点的数目比较多，整体刚度矩阵的阶数通常也是很高的。那么，是否在进行计算时要保存整体刚度矩阵的全部元素？能否根据整体刚度矩阵的特点提高计算效率？ 整体刚度矩阵具有以下几个显著的特点：对称性，稀疏性，非零系数带形分布。 对称性 由单元刚度矩阵的对称性和整体刚度矩阵的集成规则，可知整体刚度矩阵必为对称矩阵。利用对称性，只保存整体矩阵上三角部分的系数即可。 稀疏性 单元刚度矩阵的多数元素为零，非零元素的个数只占较小的部分。如图3.11所示的结构，结点2只和通过单元联接的1、3、4、5结点相关，结点5只和通过单元联接的2、3、4、6、8、9结点相关。由单元刚度矩阵的物理意义和整体刚度矩阵的形成方式可知，相关结点2、3、4、6、8、9及结点5本身产生位移时，才使结点5产生结点力，其余结点产生位移时不在该结点处引起结点力。在用分块形式表示的整体矩阵中，与相关结点对应的分块矩阵具有非零的元素，其它位置上的分块矩阵的元素为零，如图3.12所示。   图3.11 图3.12整体刚度矩阵的分块矩阵示意 3）非零元素带形分布 在图3.12中，明显可以看出，整体刚度矩阵的非零元素分布在以对角线为中心的带形区域内，这种矩阵称为带形矩阵。 在包括对角线元素的半个带形区域内，每行具有的元素个数叫做半带宽，用hbd表示。  图3.11所示结构的相邻结点编码的最大差值为4，所以半带宽为10。 二维等带宽存储 设整体刚度矩阵[K]为一个的矩阵，最大半带宽为d。利用带形矩阵的特点和对称性，只需要保存以d为固定带宽的上半带的元素，称为二维等带宽存储。进行存储时，把整体刚度矩阵[K]每行中的上半带元素取出，保存在另一个矩阵[K*]的对应行中，得到一个矩阵[K*]。 把元素在[K]矩阵中的行、列编码记为r、s，在矩阵[K*]中的行、列编码记为r*、s*，对应关系如下： r*=r s*=s-r+1   图3.13(a) 图3.13(b) 如图3.13(a)所示的最大半带宽为d的整体刚度矩阵[K]，采用二维等带宽存储后得到如图3.13(b)所示的矩阵[K*]。用新的方法存储后，[K]矩阵中的对角线元素保存在新矩阵中的第1列中，[K]矩阵中的r行元素仍然保存在新矩阵的r行中，[K]矩阵中的s列元素则按照新的列编码保存在新矩阵的不同列中。 采用二维等带宽存储，需要保存的元素数量与[K]矩阵中的总元素数量之比为。所存储的元素数量取决于最大半带宽d的值，d的值则由单元结点的编码方式决定。 虽然在采用二维等带宽存储时，仍然会保存一些零元素，但是采用这种方法时元素寻址很方便。  图3.14(a) 图3.14(b) 对于同样的有限元单元网格，按照图3.14(a)的结点编码，最大的半带宽为14；按照图3.14(b)的结点编码，最大的半带宽为18；按照图3.11的结点编码，最大的半带宽为10。 3.9线性方程组解法 由于有限元分析需要使用较多的单元，线性方程组的阶数很高，有限元求解的效率很大程度上取决于线性方程组的解法。利用矩阵的对称、稀疏、带状分布等特点提高方程求解效率是关键。 线性方程组的解法分为两大类： 直接解法， 迭代解法。 直接解法以高斯消去法为基础，包括高斯消去法、等带宽高斯消去法、三角分解法，以及适用于大型方程组求解的分块算法和波前法等。 迭代算法有高斯-赛德尔迭代、超松弛迭代和共轭梯度法等。 在方程组的阶数不是特别高时，通常采用直接解法。当方程组的阶数过高时，为避免舍入误差和消元时有效数损失等对计算精度的影响，可以选择迭代方法。 这里给出了用Fortran语言编写的等带宽高斯消去法的代码，其中NROW为矩阵行的数目，NHBW为最大半带宽。 SUBROUTINE SOLVERB(NROW,NHBW,STIFF,DISPL) C Band elimination method C .................................... C : SOLVE EQUATIONS [KS]*{H}=(Q) : C :..................................: C [KS] is stored with the half band width method C Input: C NROW - the quantity of rows C NHBW - half band width C STIFF[NROW,NHBW] C - coefficient matrix stored with the half width method C DISPL[NROW] - the right hand vector C Output: C DISPL - results of variables C Warning: C array STIFF is changed. C***ADD:DPR*** IMPLICIT DOUBLE PRECISION ( A-H,O-Z ) C***END:DPR*** C C INTEGER NROW,NHBW DIMENSION STIFF(NROW,NHBW),DISPL(NROW) IN = NROW KD = NHBW DO K = 1, IN-1 IF((K+KD-1).GE.IN) THEN IM = IN ELSE IM = K + KD - 1 ENDIF DO I = K + 1, IM L = I - K + 1 CC = STIFF(K,L)/STIFF(K,1) DO J = 1,KD-L+1 M = J+I-K STIFF(I,J) = STIFF(I,J)-CC*STIFF(K,M) ENDDO DISPL(I) = DISPL(I)-CC*DISPL(K) ENDDO ENDDO DISPL(IN) = DISPL(IN)/STIFF(IN,1) DO I = IN-1,1,-1 IF(KD.GT.(IN-I+1)) THEN JM = IN-I+1 ELSE JM = KD ENDIF DO J = 2,JM IH1 = J+I-1 DISPL(I) = DISPL(I)-STIFF(I,J) * DISPL(IH1) ENDDO DISPL(I) = DISPL(I) / STIFF(I,1) ENDDO RETURN END ANSYS提供了多种求解器供选择，如图3.15所示，分为直接解法和迭代解法。 直接解法包括： 波前法（Frontal Solver），稀疏法（Sparse Direct Sovler）。 迭代解法包括： 雅可比共轭梯度法（Jacobi Conjugate Gradient Solver， JCG）， 不完全共轭梯度法（Incomplete Cholesky Conjugate Gradient Solver， ICCG） 预处理共轭梯度法（Preconditioned Conjugate Gradient Solver， PCG） 代数多格法（Algebraic Multigrid Solver，AMG） 区域分割法（Distributed Domain Solver， DDS）  图3.15 ANSYS提供的方程组求解方法