弹性力学平面问题有限元课件

1、第8章弹性力学平面问题的有限元分析8.1概述在有限元法中，把单元与单元之间设置的相互连接点，称为结点。一般用号码1、2、进行结点编号。有了结点，才可将实际连续体看成是仅在结点处相互连接的单元集合组成的离散型结构。由单元、结点、结点连线构成的集合称为有限元模型。 8.1.1单元划分类型单元类型：三角形、四边形单元数目：根据计算精度要求来确定结点设置：使单元的的结点编号尽量靠近有限元模型：由单元、结点、结点连线构成的集合8.1.2位移函数在选择多项式时，为了使有限单元法的计算精度和收敛性得到保障，还需要满足完备性和连续性的要求。为了使位移模式尽可能地反映物体中的真实位移形态，它应满足下列条件：（1

2、）位移模式必须能反映单元的刚体位移；（2）位移模式必须能反映单元的常量应变；（3）位移模式应尽可能地反映位移的连续性。弹性力学平面问题一般选择多项式函数作为位移函数。 8.2平面三角形单元三角形单元是一种简单方便、对边界适应性强的单元，由于以三角形的三个顶点作为结点，因此又成为三结点三角形单元。这种单元的计算精度较低，使用的时候必须进行精细的网格划分，但他仍然是一种常用的单元 (a) 均匀受力板力学模型 (b) 力学模型离散化 平面问题有限单元法的计算力学模型根据完备性和连续性的要求，选取3结点三角形单元的位移场函数如下： （8-1a）将3个结点上的坐标和位移分别代入式（8-1a）就可以将六个

3、待定系数用结点坐标和结点位移分量表示出来。将水平位移分量和结点坐标分别代入（8-1a）中的第一式，得到 （8-1b）8.2.1位移函数的选取如图所示的3结点平面三角形单元，结点的坐标分别为 、 、 ，结点位移分别为 、 、 、。 记单元的结点位移向量 和结点力向量 为： （8-1a）写成矩阵形式，有： （8-2）（8-3）为了避免出现 的情况，三个结点 按逆时针顺序排列。 其中A，为三角形单元的面积。将（8-3）代回（8-1） 整理后可得：（8-1）（8-3） （下标 轮换）（8-4） 把称为形函数矩阵，称为形函数。 单元内的位移场函数可以简写成：写成矩阵形式（8-4）（8-5b）形函数 是单

4、元内任意一点坐标的线性函数。单元内位移场的插值函数。也就是说，单元的节点位移通过N(x,y)控制着单元的位移场的形态。所以N(x,y)称为单元的形态函数或形函数。插值函数具有如下性质：（1） 在节点 i 上值为1，在节点 j 和 m 上值为0，即有：也即有由上式可见，当即在节点i，应有其他两个形函数也具有同样的性质。因此也必然要求（2）在三角形单元的一条边上，如i j 上，形函数与第三个顶点m的坐标无关（3）在单元中任一点各插值函数之和应等于1，即：因为若单元发生刚体位移，如x方向有刚体位移 ，则单元内（包括结点上到处应有位移 ，即 ，又有：（4）对于现在的单元，插值函数是线性的，在单元内部及

5、单元的边界上位移也是线性的，可由结点上的位移值唯一地确定。由于相邻单元公共结点的结点位移是相等的，因此保证了相邻单元在公共边界上位移的连续性。8.2.2单元的应变场由几何方程知（8-6）根据单元的位移场函数式（8-5），由几何方程可以得到单元的应变场表达式：（8-5）求偏导数简写为： 由于矩阵B是常量，单元内各点应变分量也都是常量，这是由于采用了线性位移函数的缘故，这种单元称为常应变三角形单元。8.2.3单元的应力场由物理方程,可以得到单元的应力场表达式： 其中 为应力矩阵，称为弹性矩阵，对于平面应力问题，(8-9) (8-8) 将应力矩阵表示为分块矩阵的形式，有： 其中：(8-10) 可以看

6、出，应力分量也是一个常量。在一个三角形单元中各点应力相同，一般用形心一点表示。其应变也可同样表示。平面应变问题：8.2.4单元刚度矩阵最小势能原理（1）单元的应变能 (8-11) 单元的刚度矩阵：结点力和结点位移间的关系。 8.2.4单元刚度矩阵式中 分别表示单位体积的体积力、单元上的表面力、单元结点上的结点荷载。(8-12) （2）单元上外力的势能单元的总势能为：利用最小势能原理，取结点位移的变分，得到： 由的任意性，有：(8-13) 考虑到的对称性，对式（8-13）求偏导得到：记：(8-14) (8-15) (8-13) 则式（8-14）可写为：这就是描述单元结点力和结点位移向量之间关系的

7、平衡方程。其中 称为单元刚度矩阵。在3结点等厚三角形单元中 和 的分量均为常量，则单元刚度矩阵可以表示为：其中t、分别为单元的厚度和面积。单元刚度矩阵 可以表示为分块矩阵的形式：对于平面应力问题，其中：单元刚度矩阵仅与单元的几何特征和材料性质有关。 单元的横截面积A、惯性矩I、单元长度L、单元的弹性模量E有关。b. 单元刚度矩阵是一个对称阵。在单元刚度矩阵对角线两侧对称位置上的两个元素数值相等，即，根据是反力互等定理。c. 单元刚度矩阵是一个奇异阵。d. 主元恒正。e. 带状f. 单元刚度矩阵可以分块矩阵的形式表示。具有确定的物理意义。单元的刚度矩阵的性质 当单元的第j个结点位移为单位位移而其

8、他结点位移为零时，需要在第i个结点位移方向上施加的结点力的大小。 单元的刚性大，则使结点产生单位位移所需施加的结点力就大。因此，单元刚度矩阵中的每个元素反映了单元的刚性大小，称为刚度系数。 对于单元刚度矩阵的每一列（行）元素应有：单元的刚度矩阵的性质 用虚位移原理来建立结点力和结点位移间的关系式，从而得出三角形单元的刚度矩阵。 (a) 实际力系 (b) 虚设位移弹性体虚位移原理的应用8.2.5单元刚度矩阵位移变分方程外力虚功为 结点虚位移列向量和虚应变列向量为用虚位移原理建立三角形单元的虚功方程为 由于虚位移是任意的，等号两边可左乘三角形单元的刚度矩阵可写成 在式中括号中的前两个量分别为体积力

9、、表面力移置到结点上的等效结点力，依次定义为、，即：得到了单元刚度矩阵后，需要将一系列的将单元组成一个整体结构，由单元刚度矩阵得到整体刚度矩阵的基本方法是刚度集成法，即整体刚度矩阵是单元刚度矩阵的集成。8.2.6 整体刚度矩阵得到整体刚度矩阵。整体分析包括以下4个步骤：（1）建立整体刚度矩阵，（2）根据支承条件修改整体刚度矩阵，（3）解方程组，求出结点的位移，（4）根据结点位移，求出单元的应变和应力。 结构的平衡条件可用所有结点的平衡条件表示。假定i 结点为结构中的任一公共结点，则该结点平衡条件为：i 结点的结点力列向量 i结点的载荷列向量。 单元编号单元结点局部编号单元结点整体编号（1）31

10、2（2）524（3）532（4）356整体编号123456局部编号00010000002000300000040005000600000000000000000由上表可以看出，整体刚度矩阵: 具有对称性;稀疏矩阵，即其中有大量的零元素，并且非零元素都集中于主对角线附近呈带状。和单元刚度矩阵一样，由于位移函数中包含刚体位移，所以整体刚度矩阵也是一个奇异矩阵。必须要排除刚体位移后，才能变为正定矩阵。 8.2.7约束条件的处理将方程分组（1）对角线元素改1，同行同列其他元素改0（2）对角线元素乘以一个大数（1）对角线元素改1，同行同列其他元素改0若第r个位移分量 为已知值 ,则将整体刚度矩阵中主对角

11、线元素 改为 1 ，第r行、第r列的其他元素改为 0 ，荷载矩阵中第r行的元素改为 ，其他元素都减去结点位移的已知值和原中这行相应列元素的乘积。 （3）对角线元素乘以一个大数若第r个位移分量 为已知值 则将整体刚度矩阵中主对角线元素 乘以一个大数（如 ） ，第r行、第r列的其他元素不变 ，荷载矩阵中第r行的元素改为 ，其他元素不变。 8.3平面矩形单元矩形单元也是常用的单元之一，由于采用了比常应变三角形单元更高次数的位移模式，故可以更好地反映弹性体的位移状态和应力状态。如下图所示四结点矩形单元，记单元的结点位移向量 和结点力向量为： 为了能推导出简洁的结果，在这里引入无量纲坐标： 8.3.1单

12、元位移场结点条件共有8个，因此，x和y方向的位移场可以各有4个待定系数，可以取以下多项式作为单元的位移场模式：它们是具有完全一次项的非完全二次项，其中以上两式中右端的第四项是考虑到x和y方向的对称性而取的。 (8-17)由结点条件，在 处，有 将(8-18)代入(8-17)中，可以求解出待定系数，然后再代回(8-17)中经整理后，有(8-18)(8-19)其中，N为单元的形函数矩阵，如以无量纲坐标系来表达，则(8-20)式可以写成其中：(8-20)(8-21)8.3.2单元应变场根据单元的位移场函数式（8-19，由几何方程可以得到单元的应变场表达式，记为：(8-22)(8-23)这里，B矩阵称

13、为几何矩阵。B矩阵可以表示为分块矩阵的形式其中 (8-25)(8-24)8.3.3单元应力场由物理方程，可以得到单元的应力场表达式，其中为应力矩阵，D称为弹性矩阵，对于平面应力问题，(8-26)将应力矩阵表示为分块矩阵的形式其中：对于平面应变问题，只需将E换为，换为。 （8-28）（8-27）8.3.4单元刚度矩阵和三角形单元一样，可以根据最小势能原理导出结点位移向量和结点力向量之间关系，即单元的刚度矩阵，可以将其写成分块的形式。其中（8-30）（8-29）对于平面应力问题，如果单元厚度t为常数，则得到（8-56）的显式形式： （8-31）8.4 广义坐标有限单元法选择单元位移函数的一般原则

14、单元中的位移模式一般采用以广义坐标为待定参数的有限多项式作为近似函数，如3结点三角形单元的格式，有限项多项式选取的原则应考虑以下几点：1、广义坐标为待定参数，由结点场变量确定，因此它的个数与结点自由度应相等。3结点三角形单元 4结点矩形单元2、选取多项式时，常数项和坐标的一次项必须完备。常数项和坐标的一次项反映了单元常应变及刚体位移特性。3、多项式选取，由低到高，多项式具有对坐标的对称性。几何各向同性8.5 有限元解的性质和收敛性8.5.1 有限元解的收敛准则 有限单元法作为一种数值方法可以认为是里兹法的一种特殊形式，不同在于有限单元法的试探函数是定义于单元而不是全域。因此，有限元解的收敛性可

15、以与里兹法的收敛性对比进行讨论。 以一个待求的标量场为例，其微分方程为相应的泛函为在有限单元法中，场函数的总体泛函是由单元泛函集成的。如果采用完全多项式作为单元的插值函数，则有限元解在一个有限尺寸的单元内可以精确地和真正解一致。但是实际上有限元的试探函数只能取有限项多项式，因此有限元解只能是真正解的一个近似解答，我们需要研究当单元尺寸趋于零时，有限元解趋于真正解。常应变要求完备性要求：在单元内，协调性要求：在单元交界面，有限元解收敛反映常应变及刚体位移8.5.2 位移元解的下限性质8.6平面问题程序设计有限元方法在分析工程问题时的优越性，在于与计算机的结合。通过编写的计算程序，确定所研究问题的

16、信息和数据，问题就可以得到解决。前面两节讲述了平面问题的有限元计算原理。本节主要介绍采用三角形单元，计算在给定荷载作用下的弹性力学平面静力问题时，有限元法的程序编制与框图设计。通常有限元法的程序设计和编写分为以下几步： （2）框图设计计算方案选定以后，需先画出总框图。平面问题的有限元法的总框图如所示。 （1）问题分析使用计算机解决具体问题时，首先要根据所要求解的问题，找出已知的数据和条件，建立数学力学模型，选择计算方案，准备所需要的原始数据。 （3）设计程序根据总框图的规划和有限元法的基本原理，用自己所使用的程序语言把这个解决方案严格地描述出来，也就是编写出程序代码。它就是用于求解问题的源程序

17、。 （4）程序的编译和调试将编写好的源程序输入计算机，进行编译，改正其中的语法等错误。并用实际的输入数据对编好的程序进行调试，分析所得到的运行结果，进行程序的测试和调整，直至获得预期的结果。 8.7 ANSYS平面结构计算示例 8.7.1问题描述 如图8-7所示长方形板ABCD，板厚0.04m，孔半径r=0.2m，E=210GPa，泊松比=0.3，约束条件：在长方形底边AD约束全部自由度, BC边施加垂直向下均布载荷g=10000000N/m。图8-7 长方形板结构 8.7.2 ANSYS求解操作过程 打开Ansys软件，在Ansys环境下做如下操作。 设置相关工作目录和工作文件名，注意：工作

18、目录不要用中文 定义分析类型 选择 Main Menu Preferences 命令，弹出如下优选对话框（1）设置相关工作目录和工作文件名选择 Utility Menu File Change Title 命令，弹出如下对话框 输入分析标题： Plane strs Analysis （2）定义分析标题图2-8 单元类型对话框 （3）选择单元类型 运行PreprocessorElement TypeAdd/Edit/Delete，弹出Element Types对话框，如图2-8所示。单击Add，弹出Library of Element Types窗口，如图2-9所示，依次选择Structural

19、 Solid，Quad 8 node 82，单击OK。图2-9 单元类型库对话框 在Element Types对话框中，如图8-10所示，单击Options，弹出如图8-11所示对话框，设置K3选项栏为Plane strs w/thk，设置K5选项栏为Nodal stress，设置K6选项栏为No extra output。表示单元是应用于平面应力问题，且单元是有厚度的。图8-10 单元类型对话框 图8-11 PLANE82 单元选项设置对话框 （2）定义实常数 运行PreprocessorReal ConstantsAdd /Edit/Delete，弹出如图8-12所示对话框，点击Add，弹

20、出如图8-13所示对话框，点击OK，弹出如图8-14所示对话框，在THK选项栏中设置板厚度为0.04m。设置完毕单击OK按钮。图8-12 实常数对话框 图8-13 选择要设置实常数的单元类型 图8-14 PLANE82实常数设置 （3）设置材料属性 运行PreprocessorMaterial PropsMaterial Models，弹出如图8-15所示对话框，依次双击Structural，Linear，Elastic，Isotropic，弹出图8-16所示对话框，在EX选项栏中设置数值2.1e11，在PRXY选项栏中设置数值0.3。设置完毕单击OK按钮。图8-15 选择材料属性对话框 图8

21、-16 设置材料属性对话框 （4）建立模型 运行PreprocessorModelingCreateAreaRectangleBy 2 Corners，弹出如图8-17所示对话框，设置参数，WP X选项栏中填写0，WP Y选项栏中填写0，Width选项栏中填写1.5，Height选项栏中填写1，单击OK。继续运行PreprocessorModelingCreateAreaCircleSolid Circle，得到如图8-18所示对话框，在WP X选项栏中填写0.75，WP Y选项栏中填写0.5，在Radius选项栏中填写0.2，设置完毕点击OK按钮。图8-17 建立矩形对话框图8-18 创建实

22、 心圆对话框 进行布尔运算：PreprocessorModelingOperateBooleansSubtractAreas，先选矩形面单击OK，再单击圆面，单击OK。得到如图8-19所示图形。 图8-19 长方形板模型 （5）划分网格 运行MeshingSize CntrlsManual SizeAreasAll Areas，弹出如图8-20所示对话框，在SIZE选项栏中填写0.05，点击OK按钮。图8-20 设置网格尺寸对话框 66 运行MeshMesh Tool，弹出如图8-21所示对话框，在Shape选项栏后面，选择Tri和Free，单击Mesh.划分网格，网格划分如图8-22所示。图8-21 网格划分对话框 图8-22