《有限元素法》PPT课件.ppt

计算机在材料中的应用,主讲教师 武建军,第四章 有限元素法,数值方法的本质就是把计算域内有限位置（节点）上的因变量值当作基本未知量来处理，任务是提供一组关于这些未知量的代数方程并规定求解这些方程的算法。 有限元法是以变分原理为基础吸取差分格式的思想而发展起来的一种有效的数值解法。 有限元素法（FEM）是求解微分方程的一个常用数值解法。,FEM的应用领域,50年代，有限元法首先从结构力学中得到应用，而后推广到流体力学、热传导等领域,现在已被广泛用于物理和工程设计计算的很多领域。该方法逐渐成为求解许多数学物理方程边值问题的高效能的方法。,1 FEM的基本思想及其特点,利用FEM方法求解问题包括单元划分、假定单元分布函数、列出单元方程、联立求解等几个步骤。 例如解决一个超静定梁的问题，每一个杆件就可以作为一个单元来处理,可以列出自己的载荷与位移的关系方程(单元方程),然后联立求解。,1-2 FEM的基本思想-里兹法,用插值（或其它）方法用节点处的函数值构建单元分布函数； 将代入原来的微分方程，将得到有关节点变量值的代数方程组。 求解代数方程组，获得节点函数值。 离散化方程仍受未知函数的微分方程支配。,单元划分的意义,在整个计算域内寻找试探函数不太容易。如果将计算域划分成一定数量的单元，单元分布函数容易建立。 因为选取分段分布的结果，一定的离散方程只与少数几个节点有关。一方面分布函数容易构造，而且也容易求解。,分布函数的随意性,由于分布假设以及推导方法的不同，一个微分方程的离散化方程不是惟一的。 因为只需要满足本质性边界条件，而不必考虑自然边界条件（第二、第三类边界条件自动满足），试探函数的选取是比较容易的。 试探函数阶次提高，解的精度也提高。,当网格特别细密时，相邻节点之间的变化就很小，因此单元内分布假设的实际细节变得不再重要。离散化方程的解将趋近于相应微分方程的精确解。,单元形状,单元的形状没有限制。例如椭圆孔应力集中集中问题，图中将它划分为三角形网络。把原来的连续体简化为由有限个三角形单元组成的离散体。其中三角形单元之间只在节点处用铰链相连，把载荷按照静力等效原则也转移到节点处。,1-2 FEM的特点,物理概念清晰，特别是对于力学问题。 灵活性与通用性。由于单元形状灵活，易于处理复杂区域、复杂边界条件。 而对于具有规则的几何特性和均匀的材料特性问题,差分法的程序设计比较简单,收敛性也比有限元法好。 有限元法同时具有里兹法与差分法的优点，使变分问题的直接解法变成了工程计算中的现实。,FEM的特点,有限元素方法是物理量的矩阵分析方法在连续体中的有效推广。每个元素都采用有限个参数来描述它的物理特性。 有限元素方法是基于虚功原理,或者说是变分原理。它不象差分法那样直接去解场方程,而是求解一个虚功取极小值的变分问题。 FEM是解决复杂区域、边界条件数学物理方程边值问题的一种比较完美的离散化方法。,本章提要,有限元素法的原理概要 举例说明如何运用有限元素法 弹性力学平面问题 热传导问题 比较有限元素法和有限差分法。,2 弹性力学基础知识,弹性力学研究宏观均匀、各向同性固体的弹性变形，例如刃型位错应力场计算。 严格地说，任何弹性体总是处在空间应力状态，因而实际问题都是三维空间问题。但是，有些弹性力学问题可以简化为平面问题。,平面应力问题,例如平面薄板的拉伸变形问题，由于厚度很小，而载荷又平行板面且沿厚度方向均匀分布，因此可以近似认为沿厚度方向的应力分量等于零。,平面应变问题,水库大坝的长度比高度和宽度要大得多，而载荷又都与横断面平行且沿长度方向均匀分布，可以认为沿长度方向的应变分量等于零，这种问题称平面应变问题。,2-1 变形参数,单元体弹性变形参数包括 位移沿x,y,z轴的三个分量 u,v,w 力或载荷沿x,y,z轴的三个分量 X,Y,Z 应变ij（作用面垂直于i轴，指向j方向）i,j=x,y,z 应力ij(规定同应变 ),一般规定,当i=j, ij为正应变,表示线段伸长或缩短,可简化为i ，如x。一般规定，伸长应变为正值。 ij， ij表示切应变（角应变），表示两线段之间夹角的变化。一般规定，直角变成锐角切应变大于零。 应变分量有9个,一般有6个独立分量。 应力分量与应变分量类似，也有9个。,2-2 几何方程,2-3 广义虎克定律,广义虎克定律,拉梅方程的其它形式,广义虎克定律,平面应力问题 的虎克定律,平面应变问题的虎克定律,2-4 力学平衡方程,应力与体积力，如重力之间的平衡 其中fi分别是体积力在i方向上的分量。,力学平衡方程,应力与表面力之间的平衡 式中X、Y、Z是表面力的三个分量，l,m,n是表面外法线的方向余弦。,3 变分方法与虚功原理,有限元法建立离散方程的方法有三类。 直接法 例如，超静定桁架问题，每个组件就是一个元素。易于理解,但只能用于较简单的问题，实际用途不大。 变分法 把有限元法归结为求泛函的极值问题。使有限元法建立在更加坚实的数学基础上,扩大了有限元法的应用范围。 加权余量法 直接从基本微分方程出发,求得近似解。对于根本不存在泛函的工程领域都可采用,从而扩大了有限元法的应用范围。,3-1 变分原理,变量与变量间的关系称为函数。 如果对于某一类函数y(x)中的每一个函数y(x),I都有一值与其对应，则变量I叫做依赖于函数y(x)的泛函记为Iy(x)。 泛函是函数集合的函数,也就是函数的函数。,泛函举例-曲线长度,在直角坐标系中两定点M和N,连接两点曲线的长度L与曲线的形状y=y(x)有关。曲线长度是由曲线的方程y=y(x)所确定,因此，曲线的长度L就是曲线函数的一个泛函，可以记为Ly(x)。,泛函举例-应变能,弹性体在各种力作用下引起弹性应变（x）, 物体具有不同的位移函数（或应变分布函数）具有不同的应变能，那么应变能就是一个依赖位移函数的泛函。,其它例子,外力势能F也是位移函数的泛函 f为体积力，q为面积力。 此外总势能P=V-F也是位移函数的泛函。,变分的定义,如果对应于函数y(x)的微小改变，有泛函的微小改变与之对应，称泛函是连续的。 函数的变分y是指两个非常接近的函数y(x)与y1(x)的差 y= y(x)y1(x) 而泛函的变分是,变分问题,变分问题就是研究泛函的极值问题。这个问题的解法类似于函数极值的求法。 函数f(x)的极值可以根据df=0这个条件去寻找判断。类似于普通函数取极值的条件,若具有变分的泛函在y=y0(x)处取得极值,那末泛函在该处的变分应为零。用数学公式写出为 即可以根据变分等于零这个条件去寻找y0(x)。,一维欧拉方程,可以证明，泛函Q 实现极值的必要条件是函数y(x)满足一维欧拉方程,二维欧拉方程,若f(x,y)在R内二阶可微，区域R的边界分为B、C两部分，在B边界满足f=f1(x,y)，C边界上为未知函数。泛函 函数G（f）沿边界C取值。 假定泛函有极值，必须满足欧拉方程 在区R内， 在边界上， 式中l,m是边界外法线的方向余弦。,变分问题举例,变分原理,由于使泛函实现极值的函数必然满足相应的欧拉方程，这样可以把求解某些微分方程的问题转化为泛函极值问题。如果求得了使泛函实现极值的函数，也就是微分方程的解函数。 变分方法的主要难点就是寻找适当形式的泛函（许多微分方程不存在对应的泛函）。泛函的寻找一般要根据物理概念，如虚功原理、最小能量原理等来进行。,弹性力学问题的变分原理,最小总势能原理指出，在满足边界条件的所有位移中，以保持力的平衡状态的位移造成的总势能最小。 因此，求平衡态的位移函数就等同于求应变能（位移函数的泛函）的极小值，即P(u0(x)=0的解。这就是弹性力学问题有限元方法的变分原理。,当单元内任意一点的函数值根据插值法用节点函数值表示时，单元积分也成为节点函数值的函数。 由于求解域R已经划分成若干单元，泛函变成这些子域上积分的和。显然，泛函也成为节点函数值 （i=1,2,3n）的函数。 由泛函极值的必要条件，得方程组 ，i=1,2,3,n,3-2 虚功原理,设变形体处于平衡受力状态，体积力为f,在自由边界上的表面力为q。设变形体产生任一虚位移u*，相应的虚应变为*=Bu*，则体积力和表面力所作的虚功，恒等于应力在虚应变上所作的虚变形功（应变能增量）,如果外力是集中力P1,P2，P3Pn,而各力相应的虚位移分量*i, 以上结果与变分原理的结果相同。,3-3 加权余量法,用变分法求解微分方程，首先要找到相应的泛函。 对于有些问题相应的泛函尚未找到，或者根本不存在相应的泛函。在这种情况下，就无法用变分法求解。 加权余量法（也称加权余值法）是求微分方程近似解的一种有效方法。,设有微分方程 假设有一个满足边界条件和具有一定连续程度的试探函数 （其中含有若干待定系数）使 余量R与所选的试探函数有关，希望R在某种意义上比较小。通过把余量与加权函数Wi(x)正交化的途径，化为代数方程组，即令,加权函数,比较常用的加权函数有 狄拉克函数(x),配置法； 幂函数 ,i=0,1,2,矩量法； 里兹基函数，Galerkin法；精度高，常用 余数R,最小二乘法。 有限元法中积分在单元内进行。,4 弹性力学问题的FEM,单元分析 总体分析 边界条件处理 方程求解,4-1 单元分析,单元分析的任务是推导基本未知量节点位移与节点力F之间的关系。 一般可从节点位移、单元内各点位移(插值)、应变（几何方程）、应力（虎克定律）入手，最终求得节点位移与节点力F之间的关系,平面三角形单元,图示一个三角形单元 三个节点按逆时针方向的顺序编码为i,j,m。 节点坐标分别为（xk,yk）k=i,j,m。,在弹性力学平面问题中，每个节点有两个位移分量，因此三角形单元的位移向量可写成 与此对应的物理量是六个节点力分量 单元分析的任务就是推导以上两向量之间的关系，即 其中转换矩阵Ke称单元刚度矩阵。,位移插值函数,单元（元素）分析的第一步是由单元节点位移分量推算单元内部任意一点的位移。在任一单元内可以通过线性插值构造三角形内部任意一点的位移u和v，即可设在,代入节点位移可得到方程组，解方程组得到系数值,其中 按ijmi的顺序轮换可得到其它参数。三组常数都是与节点坐标有关的常数. A是三角形元素的面积 ak,bk,ck是行列式第k行对应元素的代数余子式。,位移插值函数,将系数i代入u,v表达式，整理可得 式中N是插值基函数,也称形函数 l=i,j,m。,形函数的性质,形函数描述了单元内部研究对象例如位移函数的分布情况。 在单元任意一点，三个形函数的和等于一（设三节点位移相等，就可以证明）。 在单元的边上（如ij边），有一个形函数（如Nm）等于零。 在单元节点上,单元刚度矩阵,（1）由位移求应变 几何方程用矩阵表示，B称几何矩阵。由于插值函数是一次函数，在单元内的应变分量都是常量。,（2）位移与应力 根据虎克定律， 其中转换矩阵S=DB (3)由节点位移求节点力 根据虚功原理 ，,单元刚度矩阵的性质,1 单元刚度矩阵是66对称矩阵,符号规定,刚度矩阵中kpq,p代表节点力编号和方向， q代表位移的编号与方向 带代表水平方向，不带代表垂直方向。 例如， 代表m节点力垂直方向分量Vm与i节点水平位移分量ui之间的作用系数。,2 K是奇异矩阵，对应的行列式等于零，不存在逆矩阵。实际无法独立求解 从几何关系也可得到这一结论。如果给定节点位移，按照这个关系可以计算惟一的节点力。但是反过来，根据节点力却不能得出节点位移的唯一解。因为单元可以产生任意的刚体位移-不影响节点力。,4-2 整体分析,平面弹性力学有限元法计算的最终结果是要求出区域R中的位移分布、应变分布及应力分布。 由于已经把解域R划分成有限个三角形单元及若干个节点，并把整个场离散到这些节点上，因此，我们的任务就是把这些节点上的位移、应变和应力求出来。,整体分析的任务,由于一个节点不止属于一个单元，因而要将单元方程加以合并组合，得到反映所有节点节点力与位移关系的总方程组。即根据单元刚度矩阵建立连续体的总体刚度矩阵、根据边界条件建立节点载荷向量。 整体分析的任务就是建立整体刚度矩阵和节点载荷向量。,合成方法与步骤,在刚度矩阵合成时，存在两个问题。一是单元矩阵与总体矩阵的阶数不同；二是单元节点是局部号码。 刚度矩阵的集成分两步进行。第一步把单元刚度矩阵扩大到nn阶，然后把单元刚度矩阵中Ke的9个2 2子块搬家，变成按总体编码的9个子块,其它子块用零来填充。例如某单元局部编码i,j,m对应于总码（节点编号）2,4,5,则该单元矩阵在扩大、搬家之后将变成（设n=7）,扩维,叠加,第二步，把扩充以后的单元矩阵迭加，即得到整体刚度矩阵。为了节省存储容量，以上两步实际上是穿插进行的，即按对号入座、边搬家、边累加的方法进行集成。 节点载荷向量的构建方法与刚度矩阵相似。由作用力与反作用力的关系，大多数节点的节点力的合力都等于零，只有那些第一类边界节点或受到体积力的节点才不等于零。,4-3 第一类边界条件的处理,在有限元法中，自然边界条件（包括第二 、第三类边界条件）已经包含在变分公式当中，因而是自动满足的。 第一类边界条件必须强制性地引入到有限元方程当中，即对有关边界节点的方程进行修改。为了修改方便，总体节点编号一般先编内部节点，最后编第一类边界节点。,例如节点9的水平位移 ,只需要将K中第17行主对角元素改为1,其余全为零，对应的载荷向量也改成零。 例如节点7的垂直位移 ,则可进行如下修改。将K中第14行主对角元素乘以一个较大的数M（例如108），对应的载荷向量也改成 ,其余各项不变。这时，这个方程中除包含大数M的两项外，其余各项相对较小，可以忽略不计，可以得到 。,4-4 方程特点与求解,1.整体刚度矩阵的特点 （1）整体刚度矩阵是一个对称矩阵。 （2）整体刚度矩阵是一个稀疏矩阵，它的大多数元素都是零。 （3）整体刚度矩阵中的非零元素分布在以主对角线为中心的斜带型区域内。 利用矩阵的对称性和带型分布规律，在计算机中可以只存储上半带元素。,2.方程的求解 可以用带消去法或高斯赛德尔迭代法等求解方程组，得到节点位移。 然后再用单元分析公式计算每个单元的应变或应力。,5 弹性力学问题有限元解法,下面以一个简单的模型，示范弹性力学平面问题有限元解法的基本程序。 如图所示的试件，受到一组集中力。,5-1 单元分析,第一单元，节点坐标分别是i=1(0,20),j=2(0,10),m=3(10,10)（逆时针方向编号），面积A=50,弹性系数,设为平面应力问题，泊松比0.28, E=21000，板厚t=10,1，2，4单元刚度矩阵,对于第二单元，节点坐标分别是i=2(0,10),j=4(0,0),m=5(10,0),A=50 B及K与第一、四单元相同,第三单元，节点坐标分别是 i=2(0,10), j=5(10,0), m=3(10,10) A=50,单元3刚度矩阵,5-2 整体分析,第三单元刚度矩阵扩容后的形式,整体刚度矩阵,强加第一类边界条件后刚度矩阵,5-3 平衡方程及其解,其中，,节点位移,解方程，可以得到节点位移。,在单元内 利用位移与应变的关系（几何方程）可得到节点应变； 根据虎克定律可以计算节点应力。 例如对于第一单元，根据节点位移，计算1 2 3三节点的应变和应力,节点应变,节点应力,6二维热传导问题FEM解法,插值函数中的系数,A为三角形单元的面积,因为传热问题没有对应的泛函，以下我们采用加权余量法构造单元方程。 应用加辽金法，令,因为所取的插值函数为一次函数，首先应用分部积分公式，把二阶导数项化为一阶导数 S是单元e的周界。,如果单元周界外法线n与x轴的夹角余弦为nx, 与y轴的夹角余弦为ny,上式中沿单元周界的积分项可以写成，,此时，微分方程变成如下形式 如果温度插值函数使得温度在每个单元之间连续的话，则经过单元集合，相邻边界的积分必将互相抵消，上式中的单元积分仅在边界单元才有作用。,单元温度分布,单元方程,对于内部单元，不必计算曲线积分 对于热流边界单元,换热边界单元方程,规定记号,6-2 单元分析,对于常用的三角形线性单元，由于 l=i,j,m 系数都是节点坐标的函数。一阶导数 均为常数，可以提到积分号之前,单元热传导矩阵,单元热交换矩阵,根据里兹基函数的性质,单元热容矩阵,向量,单元方程,对于内部单元 对于热流边界单元 对于换热边界,整体传热方程,把所有单元集合可求得总体方程 其中,6-3 时间域的离散,不稳定温度场计算含有 是一个常系数微分方程，可以用差分方法求得数值解。,向后差分格式,将 代入 得到,C-N格式,加辽金格式,两点格式通式,三点向后差分公式,式中a是变步长系数。如果a=1，步长固定。 经过时间域的离散化，最后可将问题归结为求解代数方程组问题。,7 FEM与FDM的比较-处理方法,有限元素法是基于变分原理,同时吸收了差分法中的区域划分的思想。它将所要求解的物理问题化为泛函求极值的问题,再通