1、第四节 有限单元法,Finite Element Method (FEM),2,内容简介有限元法是结构分析的一种数值计算方法。它在20世纪50年代初期随着计算机的发展应运而生。这一方法的理论基础牢靠，物理概念清晰，解题效率高，适应性强，目前已成为机械产品动、静、热特性分析的重要手段，它的程序包是机械产品计算机辅助设计方法库中不可缺少的内容之一。,第四节 有限单元法,本节介绍了如下内容:+有限元法的基本思想平面问题有限元分析原理及步骤有限元法的设计应用及计算实例,2018/10/15,3,一、概述,在工程分析和科学研究中，常常会遇到大量的由常微分方程、偏微分方程及相应的边界条件描述的场问题，如位
2、移场、应力场和温度场等问题。目前求解这类场问题的方法主要有两种： 用解析法求得精确解； 用数值解法求其近似解。其中， 能用解析法求出精确解的只能是方程性质比较简单且几何边界相当规则的少数问题。而对于绝大多数问题，则很少能得出解析解。这就需要研究它的数值解法，以求出近似解。,2018/10/15,4,目前，工程中实用的数值解法主要有三种： 有限差分法 有限元法 边界元法其中，以有限元法通用性最好，解题效率高，工程应用最广。目前它已成为机械产品动、静、热特性分析的重要手段，它的程序包是机械产品计算机辅助设计方法库中不可缺少的内容之一。,有限差分法：微分方程和积分微分方程数值解的方法。基本思想是把连
3、续的定解区域用有限个离散点构成的网格来代替， 这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似， 积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组 ， 解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。,边界元法：边界元法是一种继有限元法之后发展起来的一种新数值方法,与有限元法在连续体域内划分单元的基本思想不同,边界元法是公在定义域的边界上划分单元,用满足控制议程的函数去逼近边界条件.所以边界元法与有限元相比具有
4、单元的未知数少,数据准备简单等优点.但用边界元法解非线性问题时,遇到同非线性项相对应的区域积分,这种积分在奇异点附近有强烈的奇异性,使求解遇到困难.边界元法是在有限元法之后发展起来的一种较精确的、有效的工程数值分析方法 。 又称边界积分方程-边界元法。它以定义在边界上的边界积分方程为控制方程，通过对边界分元插值离散，化为代数方程组求解。它与基于偏微分方程的区域解法相比，由于降低了问题的维数，而显著降低了自由度数，边界的离散也比区域的离散方便得多，可用较简单的单元准确地模拟边界形状，最终得到阶数较低的线性代数方程组。,2018/10/15,7,“ 有限元法 ” 的基本思想早在20世纪40年代初期
5、就有人提出，但 真正用于工程中则是电子计算机出现以后。“ 有限元法 ” 这一名称是1960年美国的克拉夫（Clough，R.W.） 在一篇题为 “平面应力分析的有限元法” 论文中首先使用。此后，有 限元法的应用得到蓬勃发展。到20世纪80年代初期国际上较大型的结构分析有限元通用程序多 达几百种，从而为工程应用提供了方便条件。由于有限元通用程序 使用方便，计算精度高，其计算结果已成为各类工业产品设计和性 能分析的可靠依据。,2018/10/15,有限元法的分析过程可概括如下：,连续体离散化单元分析整体分析确定约束条件有限元方程求解结果分析与讨论,Guidelines,9,二、平面问题有限单元法,
6、有限元法分析问题的主要步骤：,连续体离散化,单元分析,整体分析,第三章 用常应变三角形,1、连续体离散化,常用单元,假定将连续的弹性体分割成由单元所组成的离散体，单元间由结点连接。结构离散是有限元单元分析的基础。,第三章 用常应变三角形单元解弹性力学平面问题,2、单元分析,1)、单元位移函数（或单元位移模式）,多项式 ：,按照有限元法的基本思想：首先需设定一种函数来近似表达单元内部的实际位移分布，称为位移函数，或位移模式。,三节点三角形单元有6个自由度，可以确定 6个待定系数，故三角形单元的广义坐标位移函数为,三角形单元,位移函数：,矩阵形式：,A为三角形单元的面积,代入水平位移分量和结点坐标
7、:,可逆矩阵,伴随矩阵。将其行列式中各元素的代数余子式按行列式中各元素的顺序排列成方阵，再转置后得的方阵。,行列式。其值为面积的2倍。,第三章 用常应变三角形单元解弹性力学平面问题,T的伴随矩阵,第三章 用常应变三角形单元解弹性力学平面问题,将垂直位移分量和结点坐标代入,第三章 用常应变三角形单元解弹性力学平面问题,第三章 用常应变三角形单元解弹性力学平面问题,第三章 用常应变三角形单元解弹性力学平面问题,1、在单元结点上形态函数的值为1或为0。,2、在单元中的任意一点上，三个形态函数之和等于1。,3、 三角形单元在单元边界上的形函数与第三个顶点的坐标无关,( i, j, m 轮换 ),形态函
8、数性质,为保证解答的收敛性，单元位移模式应满足以下条件,完备性条件:反映单元的刚体位移与常量应变。,2）协调性条件:相邻单元在公共边界上的位移连续，单元之间不能重叠，也不能脱离。即位移函数在单元之间连续。,第三章 用常应变三角形单元解弹性力学平面问题,2）、单元载荷移置,移植载荷遵循的原则:,非结点载荷移植到结点上,虚功等效原则,是指原载荷与结点载荷在任何虚位移上所做的虚功二者相等,第三章 用常应变三角形单元解弹性力学平面问题,单元的虚位移表示方法(线位移),结点载荷,实移位,虚位移,第三章 用常应变三角形单元解弹性力学平面问题,(1)集中载荷移植,由虚功等效原则,结点力作功,外力作功,第三章
9、 用常应变三角形单元解弹性力学平面问题,移植到结点上等效结点力,集中力,i,例题1：在均质、等厚的三角形单元ijm的任意一点o（0.4a，0.4a）上作用有集中载荷P=100N ，与水平方向成 =45，求单元的等效结点载荷。,1）求形函数矩阵,解：,等腰直角三角形的面积A为：,2）求单元等效结点载荷,(2) 体力的移植,令单元所受的均匀分布力为,由虚功等效原则,结点力作功,体力作功,第三章 用常应变三角形单元解弹性力学平面问题,（3）分布面力的移植,结点力作功,面力作功,由虚功等效原则,例:均质、等厚的三角形单元ijm的结点坐标如图所示，ij边上作用有沿y轴负方向呈三角形分布的载荷，载荷密度最
10、大值为q ,单 元的厚度为t，试求单元的等效结点载荷。,S,（ i ，j ，m轮换）,将i ，j ，m的坐标代入得：,（1分）,形函数矩阵为：,解：,(1）、计算形函数：,(2）、计算等效节点载荷：,在边界mj和mi上的面力为零，所以上式第二项和第三项积分应等于零。 在边界ij上的面力为:,qy,因为积分沿逆时针方向，所以有ds=dx,3）、由结点位移求单元的应变,根据单元的位移函数,由几何方程可以得到单元的应变表达式:,B矩阵称为 几何矩阵,( i, j, m 轮换 ),第三章 用常应变三角形单元解弹性力学平面问题,B矩阵可以表示为分块矩阵的形式,B矩阵称为几何矩阵 或应变转换矩阵。,( i
11、, j, m 轮换 ),称为应变矩阵,由于线性位移函数，应变矩阵为常数矩阵。因而单元中的应力与应变为常数，称这种单元为常应变单元。,4）、由结点位移求单元应力,由物理方程得：,D称为弹性矩阵,平面应力问题,称为弹性矩阵,第三章 用常应变三角形单元解弹性力学平面问题,应力矩阵,同理可以得到平面应变问题的应力矩阵,将应力矩阵分块表示为，,5）、由结点位移求单元结点力,外力作用下处于平衡状态的弹性体，如果发生虚位移，则所有外力在虚位移上做的虚功等于内应力在虚应变上做的虚功。,单元的结点力记为：,单元的虚应变为：,单元的外力虚功为 ：,单元的内力虚功为：,虚功原理：,由虚功原理得：,外力虚功,内力虚功
12、,*,第三章 用常应变三角形单元解弹性力学平面问题,定义为单元刚度矩阵。,在3结点等厚三角形单元中B和D均为常量，则单元刚度矩阵可以表示为：,6）、单元刚度矩阵,第三章 用常应变三角形单元解弹性力学平面问题,单元刚度矩阵表示为分块矩阵,r=i,j,m s=i,j,m,单元刚度矩阵的性质:,（1）对称性,（2）奇异性,第三章 用常应变三角形单元解弹性力学平面问题,3、整体分析,刚度集成法：由单元刚度矩阵中的元素累加得到整体刚度矩阵中的元素，即整体刚度矩阵是单元刚度矩阵的集成。,整体分析步骤： 1）建立整体刚度矩阵； 2）根据支承条件修改整体刚度矩阵； 3）解方程组，由已知的载荷，移置为节点力，求
13、出结点的位移； 4）根据结点位移，求出单元的应变和应力。,第三章 用常应变三角形单元解弹性力学平面问题,3.1 刚度集成法的物理意义,刚体集成法即结构中的结点力是相关单元结点力的叠加，整体刚度矩阵的元素是相关单元的单元刚度矩阵元素的集成。结点3在整体刚度矩阵的对应系数，应该是单元（1）、（3）、（4）中对应系数的集成。,第三章 用常应变三角形单元解弹性力学平面问题,3.2 刚度矩阵集成的规则,单元刚度矩阵元素取决于单元结点的局部编号顺序，必须知道单元结点的局部编号与该结点在整体结构中的总体编号之间的关系，才能得到单元刚度矩阵中的每个分块在整体刚度矩阵中的位置。将单元刚度矩阵中的每个分块按总体编
14、码顺序重新排列后，可以得到单元的扩大矩阵。,1）单元刚度矩阵中的每个分块放到在整体刚度矩阵中的对应 位置上，得到单元的扩大刚度矩阵。,单元结点的局部编号与整体的对应关系如下：,单元（2）的单元扩大矩阵的分块矩阵：,K= k(1) +k(2) + k(3) +k(4) 整体刚度矩阵：,2）将全部单元的扩大矩阵相加得到整体刚度矩阵。,4、约束条件的处理如图所示结构的约束和载荷情况，结点1、4上有水平方向的位移约束，结点4、6上有垂直方向的位移约束，结点3上作用有集中力（Px，Py）。,n结点水平方向对应的平衡方程为：,根据支承情况，方程应该换成下面的方程：,（3-32）,整体刚度在修改后可以得到以
15、下的形式，,可以采用同样的方法修改整体刚度矩阵。,第三章 用常应变三角形单元解弹性力学平面问题,修改矩阵：,整体刚度矩阵的特点与存储方法 整体刚度矩阵具有以下几个显著的特点：对称性，稀疏性， 非零系数带形分布。1）对称性 由单元刚度矩阵的对称性和整体刚度矩阵的集成规则，可知 整体刚度矩阵必为对称矩阵。利用对称性，只保存整体矩阵 上三角部分的系数即可。,单元刚度矩阵的多数元素为零，非零元素的个数只占较小的部分。如图所示的结构，结点2只和通过单元联接的1、3、4、5结点相关，结点5只和通过单元联接的2、3、4、6、8、9结点相关。由单元刚度矩阵的物理意义和整体刚度矩阵的形成方式可知，相关结点2、3
16、、4、6、8、9及结点5本身产生位移时，才使结点5产生结点力，其余结点产生位移时不在该结点处引起结点力。在用分块形式表示的整体矩阵中，与相关结点对应的分块矩阵具有非零的元素，其它位置上的分块矩阵的元素为零。,2）稀疏性,图3.11,图3.12整体刚度矩阵的分块矩阵示意,3）非零元素带形分布整体刚度矩阵的非零元素分布在以对角线为中心的带形区域内，这种矩阵称为带形矩阵。,半带宽：在包括对角线元素的半个带形区域内，每行具有的元素个数，用d表示。 d=(相邻节点的最大差值+1）2,图3.13(a),图3.13(b),把元素在K矩阵中的行、列编码记为r、s，在矩阵K*中的行、列编码记为r*、s*，对应关
17、系如下：,r*=r s*=s-r+1,二维等带宽存储,作业:,1、均质、等厚的三角形单元ijm的结点坐标如图所示，jm边上作用有沿y轴负方向按三角形分布的载荷，单元的厚度为1，求单元的等效结点载荷。,2018/10/15,将i ，j ，m的坐标代入得：,1、三角形面积:,2、计算形函数：,解:,3、计算等效节点载荷：,在边界ij和mi上的面力为零，所以上式第一项和第三项积分应等于零。,在边界jm上的面力为：,因为积分沿逆时针方向，所以有ds= -dx,2、如图所示三角形单元的结点坐标，单元的厚度为t,材料的弹性模量为E，泊松比 ，试求该单元的刚度矩阵.,3、一平面三角形薄板构件，离散为2个单元4个节点，如图所示。已知单元的编码顺序为（1，2，3），单元的编码顺序为（3，4， 1）。试分别写出：（1）单元的分块刚度矩阵；（2）单元的分块刚度矩阵；（3）总刚矩阵的分块矩阵表达式.,
