第5章 有限元法-1

null第5章 有限元法(1)第5章 有限元法(1) Ⅴ Finite Element Method (FEM)null第5章 有限元法内容简介 　 有限元法是结构分析的一种数值计算方法。它在20世纪50年代初期随着计算机的发展应运而生。 　　这一方法的理论基础牢靠，物理概念清晰，解题效率高，适应性强，目前已成为机械产品动、静、热特性分析的重要手段，它的程序包是机械产品计算机辅助设计方法库中不可缺少的内容之一。 本章介绍了如下内容: ■ 有限元法的基本思想及应用 ■ 平面问题有限元分析原理及步骤 ■ 有限元分析的前后处理 ■ 有限元法的设计应用及计算实例5.1 概述5.1 概述 在工程分析和科学研究中，常常会遇到大量的由常微分方程、偏微分方程及相应的边界条件描述的场问题，如位移场、应力场和温度场等问题。目前求解这类场问题的方法主要有两种： 　　　　　　　　● 用解析法求得精确解； 　　　　　　　　● 用数值解法求其近似解。 其中， 能用解析法求出精确解的只能是方程性质比较简单且几何边界相当规则的少数问题。 　　 　　而对于绝大多数问题，则很少能得出解析解。这就需要研究它的数值解法，以求出近似解。null 目前，工程中实用的数值解法主要有三种： 　　　　  有限差分法 　　　　  有限元法 　　　　  边界元法 其中，以有限元法通用性最好，解题效率高，工程应用最广。目前它已成为机械产品动、静、热特性分析的重要手段，它的程序包是机械产品计算机辅助设计方法库中不可缺少的内容之一。null　　“ 有限元法 ” 的基本思想早在20世纪40年代初期就有人提出，但 真正用于工程中则是电子计算机出现以后。 　　 “ 有限元法 ” 这一名称是1960年美国的克拉夫（Clough，R.W.） 在一篇题为 “平面应力分析的有限元法” 论文中首先使用。此后，有 限元法的应用得到蓬勃发展。 到20世纪80年代初期国际上较大型的结构分析有限元通用程序多 达几百种，从而为工程应用提供了方便条件。由于有限元通用程序 使用方便，计算精度高，其计算结果已成为各类工业产品设计和性 能分析的可靠依据。 null有限元法的分析过程可概括如下：● 连续体离散化 ● 单元分析 ● 整体分析 ● 确定约束条件 ● 有限元方程求解 ● 结果分析与讨论null　　1. 连续体离散化 　　连续体：是指所求解的对象（如物体或结构）。 　　离散化（划分网格或网络化）：是将所求解的对象划分为有限 个具有规则形状的微小块体，把每个微小块体称为单元，相邻两个 单元之间只通过若干点互相连接，每个连接点称为节点。 　　相邻单元只在节点处连接，载荷也只通过节点在各单元之间传 递，这些有限个单元的集合体，即原来的连续体。 　　  单元划分后，给每个单元及节点进行编号； 　　  选定坐标系，计算各个节点坐标； 　　  确定各个单元的形态和性态参数以及边界条件等。null 图5-1 所示为将一悬臂梁建立有限元分析模型的例子。 图中将该悬臂梁划分为许多三角形单元； 三角形单元的三个顶点都是节点。图5-1　 悬臂梁及其有限元模型 null2. 单元分析 连续体离散化后，即可对单元体进行特性分析，简称为单元分析。 　　单元分析工作主要有两项： 　　(1)选择单元位移模式(位移函数) 　　 用节点位移来表示单元体内任一点的位移、应变和应力，就需 搞清各单元中的位移分布。  一般是假定单元位移是坐标的某种简单函数，用其模拟内位移的分布规律，这种函数就称为位移模式或位移函数。通常采用的函数形式多为多项式。 　　 根据所选定的位移模式，就可以导出用节点位移来表示单元体内任一点位移的关系式。null　　 进行单元力学特性分析，将作用在单元上的所有力（表面 力、体积力、集中力）等效地移置为节点载荷； 　　 采用有关的力学原理建立单元的平衡方程，求得单元内节 点位移与节点力之间的关系矩阵　　　 单元刚度矩阵。 (2) 分析单元的特性，建立单元刚度矩阵 3. 整体分析  把各个单元的刚度矩阵集成为总体刚度矩阵，以及将各单 元的节点力向量集成总的力向量，求得整体平衡方程。  集成过程所依据的原理是节点变形协调条件和平衡条件。 null4. 确定约束条件 　　由上述所形成的整体平衡方程是一组线性代数方程，在求解之前，必修根据具体情况分析，确定求解对象问题的边界约束条件，并对这些方程进行适当修正。 5. 有限元方程求解 　　通过求解整体平衡方程，即可求得各节点的位移， 进而根据位移可计算单元的应力及应变。 6. 结果分析与讨论null 应用有限元法求解机械结构应力类问题时，根据未知量和分析 方法的不同，有三种基本解法： 位移法 力法 混合法null (1)位移法 　　此法是以节点位移作为基本未知量，通过选择适当的位移函数，进行单元的力学特性分析。在节点处建立单元刚度方程，再组合成整体刚度矩阵，求解出节点位移后，进而由节点位移求解出应力。 位移法优点是比较简单，规律性强，易于编写计算机程序。所以得到广泛应用，其缺点是精度稍低。 　　(2)力法 　　该法是以节点力作为基本未知量，在节点处建立位移连续方 程，求解出节点力后，再求解节点位移和单元应力。 　　力法的特点是计算精度高。 　(3)混合法 　　此法是取一部分节点位移和一部分节点力作为基本未知量，建立平衡方程进行求解。null 有限元法最初用于飞机结构的强度设计，由于它在理论上的通用性，因而它可用于解决工程中的许多问题。 　　目前，它可以解决几乎所有的连续介质和场的问题，包括热传导、电磁场、流体动力学、地质力学、原子工程和生物医学等方面的问题。 机械设计中，从齿轮、轴、轴承等通用零部件到机床、汽车、飞机等复杂结构的应力和变形分析（包括热应力和热变形分析）。 　　有限元法不仅可以解决工程中的线性问题、非线性问题，而且对于各种不同性质的固体材料，如各向同性和各向异性材料，粘弹性和粘塑性材料以及流体均能求解； 　　对于工程中最有普遍意义的非稳态问题也能求解。null　　到20世纪80年代初期，国际上已开发出了多种用于结构分析的有 限元通用程序，其中著名的有NASTRAN、ANSYS、ASKA、ADINA、 SAP等。 表5-1列出了几种国际上流行的商用有限元程序的应用范围。 null表5-1 几种有限元程序及其应用范围null表5-1(续) 几种有限元程序及其应用范围null5.2 单元特性的推导方法5.2.1 单元划分方法及原则 　　连续体离散化时，要根据设计对象的具体情况（结构物的形状、 载荷特性、边界条件等），确定单元（网格）的大小和形状、单元的数目以及划分。 图5-2　杆状单元图5-2 所示为杆状单元。可有一维、二维和三维梁单元。null图5-3 平面单元和多面体单元常用的平面单元和多面体单元，如图5-3所示。null　　单元的划分基本上是任意的，一个结构体可以有多种划分结果。但应遵循以下划分原则： 　　(1) 分析清楚所讨论对象的性质，例如，是桁架结构还是结构物，是平面问题还是空间问题等等。 　　(2) 单元的几何形状取决于结构特点和受力情况，单元的几何尺寸(大小)要按照要求确定。一般来说，单元几何形体各边的长度比不能相差太大。 　　(3) 有限元模型的网格划分越密，其计算结果越精确，但计算工作量就越大。 因此，在保证计算精度的前提下，单元网格数量应尽量少。 　　(4) 在进行网格疏密布局时，应力集中或变形较大的部位，单元网格应取小一些，网格应划分得密一些，而其他部分则可疏一些。null　　(5) 在设计对象的厚度或者弹性系数有突变的情况下，应该取相应的突变线作为网格的边界线； 　　(6) 相邻单元的边界必须相容，不能从一单元的边或者面的内部产生另一个单元的顶点。 　　(7) 网格划分后，要将全部单元和节点按顺序编号，不允许有错漏或者重复；　　(8) 划分的单元集合成整体后，应精确逼近原设计对象。原设计对象的各个顶点都应该取成单元的顶点。 所有网格的表面顶点都应该在原设计对象的表面上。所有原设计对象的边和面都应被单元的边和面所逼近。 null5.2.2 单元特性的推导方法 　　单元刚度矩阵的推导是有限元分析的基本步骤之一。 　　目前，建立单元刚度矩阵的方法主要有以下四种： 直接刚度法  虚功原理法  能量变分法  加权残数法null　 1. 直接刚度法 直接刚度法是直接应用物理概念来建立单元的有限元方程和分析单元特性的一种方法。这一方法仅能适用于简单形状的单元，如梁单元。但它可以帮助理解有限元法的物理概念。 　　图5-4所示是xoy平面中的一简支梁简图，现以它为例，来说明用直接刚度法建立单元刚度矩阵的思想和过程。图5-4　平面简支梁元及其计算模型null　　 梁在横向外载荷（可以是集中力或分布力或力矩等）作用下产生弯曲变形，在水平载荷作用下产生线位移。 　　 对于该平面简支梁问题： 　　梁上任一点受有三个力的作用： 水平力Fx， 　　　　　　　　　　　　　　　　 剪切力Fy ， 　　　　　　　　　　　　　　　　 和弯矩Mz。　　 　　相应的位移为： 水平线位移u， 　　　　　　　　 　挠度v ， 　　　　　　　　　 和转角 z 。 由上图可见： 水平线位移和水平力向右为正， 挠度和剪切力向上为正， 转角和弯矩逆时针方向为正。 通常：null　　为使问题简化，可把图示的梁看作是一个梁单元。 　　如图5-4所示，当令左支承点为节点 i ，右支承点为节点 j 时， 　　则该单元的节点位移和节点力可以分别表示为：称为单元的节点位移列阵。称为单元的节点力列阵；若 {F} 为外载荷，则称为载荷列阵。 （5-1）（5-2）写成矩阵形式为null　　显然，梁的节点力和节点位移是有联系的。在弹性小变形范围 内，这种关系是线性的，可用下式表示 或（5-3b）（5-3a）null　　上式(5-3b)称为单元有限元方程，或称为单元刚度方程，它代表了单元的载荷与位移之间（或力与变形之间）的联系； 　　式中，[K](e)称为单元刚度矩阵，它是单元的特性矩阵。 　 对于图5-4所示的平面梁单元问题，利用材料力学中的杆件受力与变形间的关系及叠加原理，可以直接计算出单元刚度矩阵[K](e)中的各系数 kst( s, t = i, j ) 的数值，具体方法如下： 　　(1) 假设 ui=１，其余位移分量均为零，即 vi=iz=uj=vj=zj=0，此时梁单元如图5-5(a)所示，由梁的变形公式得 图5-5(a)null伸缩：挠度：转角：由上述三式可以解得 根据静力平衡条件 null由式(5-3a)解得 　　(2) 同理，设vi=１，其余位移分量均为零，即ui=iz=uj=vj=zj= 0, 此时梁单元如图5-5(b)所示，由梁的变形公式得伸缩：挠度：转角：图5-5(b)null由上述三式可以解得利用静力平衡条件由式(5-3a) 解得 null　　(3)同理，设 ，其余位移分量均为零，即 此时梁单元如图5-5(c)所示，由梁的变形公式得 伸缩：挠度：转角：由上述三式可以解得图5-5(c)null利用静力平衡条件，可得 由式(5-3a)解得 　　剩余三种情况，仿此可推出。最后可以得到平面弯曲梁元的单元刚度矩阵为 这样，便求得式（5-3a）单元刚度矩阵中前三列刚度系数。null可以看出， [K](e)为对称矩阵。（5-4）2. 虚功原理法2. 虚功原理法　　下面以平面问题中的三角形单元为例，说明利用虚功原理法来建立单元刚度矩阵的步骤。 　　如前所述，将一个连续的弹性体分割为一定形状和数量的单元，从而使连续体转换为有限个单元组成的组合体。单元与单元之间仅通过节点连结，除此之外再无其他连结。也就是说，一个单元上的只能通过节点传递到相邻单元。 　　从分析对象的组合体中任取一个三角形单元： 　　设其编号为 e ， 　　三个节点的编号为i、j、m， 在定义的坐标系 xoy 中，节点坐 标分别为(x j , y j)、(xi , y i)、(xm, ym)，如图5-6所示。图5-6　三节点三角形单元null　　由弹性力学平面问题的特点可知，单元每个节点有两个位移分 量，即每个单元有6个自由度，相应有6个节点载荷，写成矩阵形式， 即 单元节点载荷矩阵：单元节点位移矩阵：图5-6　三节点三角形单元null(1)设定位移函数 按照有限元法的基本思想：首先需设定一种函数来近似表达单元内部的实际位移分布，称为位移函数，或位移模式。 三节点三角形单元有6个自由度，可以确定 6个待定系数，故三角形单元的位移函数为 （5-5）式(5-5)为线性多项式，称为线性位移函数，相应的单元称为线性单元。null上式(5-5)也可用矩阵形式表示，即 式中，{d}为单元内任意点的位移列阵。 （5-6）null　　由于节点 i、j、m 在单元上，它们的位移自然也就满足位移函数式(5-5)。 　　设三个节点的位移值分别为( ui, vi)、( uj, vj )、( um, vm )，将节 点位移和节点坐标代入式(5-5)，得 null（5-7）式中（5-8）　　由上可知，共有6个方程，可以求出6个待定系数。解方程，求得 各待定系数和节点位移之间的表达式为 为三角形单元的面积。其中： null（5-9）将式(5-7)及式(5-8)、式(5-9)代入式(5-6)中，得到 null式中，矩阵[N] 称为单元的形函数矩阵； 为单元节点位移列阵。 其中， 为单元的形函数，它们反映单元内位移的 分布形态，是x, y 坐标的连续函数，且有 （5-12）式(5-10)又可以写成（5-13）上式清楚地表示了单元内任意点位移可由节点位移插值求出。 null(2) 利用几何方程由位移函数求应变　　根据弹性力学的几何方程 ， 　　线应变 　　　 　　剪切应变 　　则应变列阵可以写成 式中，[B]称为单元应变矩阵，它是仅与单元几何尺寸有关的常量矩阵，即 （5-14）null（5-15） 上述方程(5-14)称为单元应变方程，它的意义在于： 单元内任意点的应变分量亦可用基本未知量即节点位移分量来表示。 null(3)利用广义虎克定律求出单元应力方程根据广义虎克定律，对于平面应力问题上式（5-16）也可写成（5-16）（5-17）式中， 为应力列阵； [D] 称为弹性力学平面问题的弹性矩阵，并有null则有如下单元应力方程由式(5-19)可求单元内任意点的应力分量，它也可用基本未知量即 节点位移分量来表示。（5-18）（5-19）null(4)由虚功原理求单元刚度矩阵 　　根据虚功原理，当弹性结构受到外载荷作用处于平衡状态时， 在任意给出的微小的虚位移上，外力在虚位移上所做的虚功 AF等 于结构内应力在虚应变上所存储的虚变形势能 A ，即　　设处于平衡状态的弹性结构内任一单元发生一个微小的虚位移， 则单元各节点的虚位移 为 （5-21）（5-22）（5-20）则单元内部必定产生相应的虚应变，故单元内任一点的虚应变 为null显然，虚应变和虚位移之间关系为设节点力为 则外力虚功为 （5-25）（5-23）（5-24）单元内的虚变形势能为 null根据虚功原理 因为（5-27）（5-26）代人式(5-26)，则有式中，　　 ，　　均与坐标 x, y 无关，故可以从积分符号中提出，可得：null其中，单元刚度矩阵 （5-28）式(5-28)称为单元有限元方程，或称单元刚度方程，其中 是单元刚度矩阵。（5-29）　　因为三角形单元是常应变单元，其应变矩阵[B] 、弹性矩阵[D] 均为常量，而 ，所以式(5-29)可以写成 （5-30）null式中，t 为三角形单元的厚度；  为三角形单元的面积。 对于图5-6所示的三角形单元，将[D] 及[B]代入式(5-29)，可以得到单元刚度为 （5-31）式中: [K]为6×6阶矩阵，其中每个子矩阵为2×2阶矩阵，由下式给出 （5-32）null　　按照力学的一般说法，任何一个实际状态的弹性体的总位能是 这个系统从实际状态运动到某一参考状态(通常取弹性体外载荷为 零时状态为参考状态)时它的所有作用力所做的功。弹性体的总位 能 是一个函数的函数，即泛函，位移是泛函的容许函数。　　从能量原理考虑，变形弹性体受外力作用处于平衡状态时，在 很多可能的变形状态中，使总位能最小的就是弹性体的真正变形， 这就是最小位能原理。用变分法求能量泛函的极值方法就是能量变 分原理。能量变分原理除了可解机械结构位移场问题以外，还扩展 到求解热传导、电磁场、流体力学等连续性问题。3. 能量变分原理法null　该方法是将假设的场变量的函数(称为试函数)引入问题的控制方程式及边界条件，利用最小二乘法等方法使残差最小，便得到近似的场变量函数形式。 　该方法的优点是不需要建立要解决问题的泛函式，所以，即使没有泛函表达式也能解题。 4. 加权残数法