有限元素法课程学习

FEM的应用领域50年代，有限元法首先从结构力学中得到应用，而后推广到流体力学、热传导等领域,现在已被广泛用于物理和工程设计计算的很多领域。该方法逐渐成为求解许多数学物理方程边值问题的高效能的方法。第1页/共109页第一页，共110页。1FEM的基本思想及其特点利用FEM方法求解问题包括单元划分、假定单元分布函数、列出单元方程、联立求解等几个步骤。例如解决一个超静定梁的问题，每一个杆件就可以作为一个单元来处理,可以列出自己的载荷与位移的关系方程(单元方程),然后联立求解。

第2页/共109页第二页，共110页。1-2FEM的基本思想-里兹法用插值（或其它）方法用节点处的函数值构建单元分布函数；将代入原来的微分方程，将得到有关节点变量值的代数方程组。求解代数方程组，获得节点函数值。离散化方程仍受未知函数的微分方程支配。第3页/共109页第三页，共110页。单元划分的意义在整个计算域内寻找试探函数不太容易。如果将计算域划分成一定数量的单元，单元分布函数容易建立。因为选取分段分布的结果，一定的离散方程只与少数几个节点有关。一方面分布函数容易构造，而且也容易求解。第4页/共109页第四页，共110页。分布函数的随意性由于分布假设以及推导方法的不同，一个微分方程的离散化方程不是惟一的。因为只需要满足本质性边界条件，而不必考虑自然边界条件（第二、第三类边界条件自动满足），试探函数的选取是比较容易的。试探函数阶次提高，解的精度也提高。第5页/共109页第五页，共110页。当网格特别细密时，相邻节点之间的变化就很小，因此单元内分布假设的实际细节变得不再重要。离散化方程的解将趋近于相应微分方程的精确解。第6页/共109页第六页，共110页。单元形状单元的形状没有限制。例如椭圆孔应力集中集中问题，图中将它划分为三角形网络。把原来的连续体简化为由有限个三角形单元组成的离散体。其中三角形单元之间只在节点处用铰链相连，把载荷按照静力等效原则也转移到节点处。第7页/共109页第七页，共110页。1-2FEM的特点物理概念清晰，特别是对于力学问题。灵活性与通用性。由于单元形状灵活，易于处理复杂区域、复杂边界条件。而对于具有规则的几何特性和均匀的材料特性问题,差分法的程序设计比较简单,收敛性也比有限元法好。有限元法同时具有里兹法与差分法的优点，使变分问题的直接解法变成了工程计算中的现实。

第8页/共109页第八页，共110页。FEM的特点有限元素方法是物理量的矩阵分析方法在连续体中的有效推广。每个元素都采用有限个参数来描述它的物理特性。有限元素方法是基于虚功原理,或者说是变分原理。它不象差分法那样直接去解场方程,而是求解一个虚功取极小值的变分问题。FEM是解决复杂区域、边界条件数学物理方程边值问题的一种比较完美的离散化方法。

第9页/共109页第九页，共110页。本章提要有限元素法的原理概要举例说明如何运用有限元素法弹性力学平面问题热传导问题比较有限元素法和有限差分法。第10页/共109页第十页，共110页。2弹性力学基础知识弹性力学研究宏观均匀、各向同性固体的弹性变形，例如刃型位错应力场计算。严格地说，任何弹性体总是处在空间应力状态，因而实际问题都是三维空间问题。但是，有些弹性力学问题可以简化为平面问题。第11页/共109页第十一页，共110页。平面应力问题

例如平面薄板的拉伸变形问题，由于厚度很小，而载荷又平行板面且沿厚度方向均匀分布，因此可以近似认为沿厚度方向的应力分量等于零。第12页/共109页第十二页，共110页。平面应变问题

水库大坝的长度比高度和宽度要大得多，而载荷又都与横断面平行且沿长度方向均匀分布，可以认为沿长度方向的应变分量等于零，这种问题称平面应变问题。

第13页/共109页第十三页，共110页。2-1变形参数单元体弹性变形参数包括位移沿x,y,z轴的三个分量u,v,w力或载荷沿x,y,z轴的三个分量X,Y,Z应变ij（作用面垂直于i轴，指向j方向）i,j=x,y,z应力ij(规定同应变) 第14页/共109页第十四页，共110页。一般规定当i=j,ij为正应变,表示线段伸长或缩短,可简化为i

，如x。一般规定，伸长应变为正值。i≠j，ij表示切应变（角应变），表示两线段之间夹角的变化。一般规定，直角变成锐角切应变大于零。应变分量有9个,一般有6个独立分量。应力分量与应变分量类似，也有9个。

第15页/共109页第十五页，共110页。2-2几何方程

第16页/共109页第十六页，共110页。2-3广义虎克定律

第17页/共109页第十七页，共110页。广义虎克定律第18页/共109页第十八页，共110页。拉梅方程的其它形式第19页/共109页第十九页，共110页。广义虎克定律第20页/共109页第二十页，共110页。平面应力问题的虎克定律第21页/共109页第二十一页，共110页。平面应变问题的虎克定律第22页/共109页第二十二页，共110页。2-4力学平衡方程

应力与体积力，如重力之间的平衡其中fi分别是体积力在i方向上的分量。第23页/共109页第二十三页，共110页。力学平衡方程应力与表面力之间的平衡式中X、Y、Z是表面力的三个分量，l,m,n是表面外法线的方向余弦。第24页/共109页第二十四页，共110页。3变分方法与虚功原理有限元法建立离散方程的方法有三类。直接法例如，超静定桁架问题，每个组件就是一个元素。易于理解,但只能用于较简单的问题，实际用途不大。变分法把有限元法归结为求泛函的极值问题。使有限元法建立在更加坚实的数学基础上,扩大了有限元法的应用范围。加权余量法直接从基本微分方程出发,求得近似解。对于根本不存在泛函的工程领域都可采用,从而扩大了有限元法的应用范围。第25页/共109页第二十五页，共110页。3-1变分原理变量与变量间的关系称为函数。如果对于某一类函数y(x)中的每一个函数y(x),I都有一值与其对应，则变量I叫做依赖于函数y(x)的泛函记为I[y(x)]。泛函是函数集合的函数,也就是函数的函数。

第26页/共109页第二十六页，共110页。泛函举例-曲线长度在直角坐标系中两定点M和N,连接两点曲线的长度L与曲线的形状y=y(x)有关。曲线长度是由曲线的方程y=y(x)所确定因此，曲线的长度L就是曲线函数的一个泛函，可以记为L[y(x)]。

第27页/共109页第二十七页，共110页。泛函举例-应变能弹性体在各种力作用下引起弹性应变（x）,物体具有不同的位移函数（或应变分布函数）具有不同的应变能，那么应变能就是一个依赖位移函数的泛函。

第28页/共109页第二十八页，共110页。其它例子外力势能F也是位移函数的泛函f为体积力，q为面积力。此外总势能P=V-F也是位移函数的泛函。第29页/共109页第二十九页，共110页。变分的定义如果对应于函数y(x)的微小改变，有泛函的微小改变与之对应，称泛函是连续的。函数的变分y是指两个非常接近的函数y(x)与y1(x)的差y=y(x)—y1(x)而泛函的变分是第30页/共109页第三十页，共110页。变分问题变分问题就是研究泛函的极值问题。这个问题的解法类似于函数极值的求法。函数f(x)的极值可以根据df=0这个条件去寻找判断。类似于普通函数取极值的条件,若具有变分的泛函在y=y0(x)处取得极值,那末泛函在该处的变分应为零。用数学公式写出为

即可以根据变分等于零这个条件去寻找y0(x)。

第31页/共109页第三十一页，共110页。一维欧拉方程可以证明，泛函Q实现极值的必要条件是函数y(x)满足一维欧拉方程第32页/共109页第三十二页，共110页。二维欧拉方程若f(x,y)在R内二阶可微，区域R的边界分为B、C两部分，在B边界满足f=f1(x,y)，C边界上为未知函数。泛函

函数G（f）沿边界C取值。假定泛函有极值，必须满足欧拉方程

在区R内，

在边界上，

式中l,m是边界外法线的方向余弦。第33页/共109页第三十三页，共110页。变分问题举例第34页/共109页第三十四页，共110页。变分原理由于使泛函实现极值的函数必然满足相应的欧拉方程，这样可以把求解某些微分方程的问题转化为泛函极值问题。如果求得了使泛函实现极值的函数，也就是微分方程的解函数。变分方法的主要难点就是寻找适当形式的泛函（许多微分方程不存在对应的泛函）。泛函的寻找一般要根据物理概念，如虚功原理、最小能量原理等来进行。第35页/共109页第三十五页，共110页。弹性力学问题的变分原理最小总势能原理指出，在满足边界条件的所有位移中，以保持力的平衡状态的位移造成的总势能最小。因此，求平衡态的位移函数就等同于求应变能（位移函数的泛函）的极小值，即P(u0(x))=0的解。这就是弹性力学问题有限元方法的变分原理。

第36页/共109页第三十六页，共110页。当单元内任意一点的函数值根据插值法用节点函数值表示时，单元积分也成为节点函数值的函数。由于求解域R已经划分成若干单元，泛函变成这些子域上积分的和。显然，泛函也成为节点函数值（i=1,2,3……n）的函数。由泛函极值的必要条件，得方程组，i=1,2,3,……n第37页/共109页第三十七页，共110页。3-2虚功原理设变形体处于平衡受力状态，体积力为{f},在自由边界上的表面力为{q}。设变形体产生任一虚位移{u*}，相应的虚应变为{*}=[B]{u*}，则体积力和表面力所作的虚功，恒等于应力在虚应变上所作的虚变形功（应变能增量）

第38页/共109页第三十八页，共110页。如果外力是集中力P1,P2，P3……Pn,而各力相应的虚位移分量*i,以上结果与变分原理的结果相同。

第39页/共109页第三十九页，共110页。3-3加权余量法用变分法求解微分方程，首先要找到相应的泛函。对于有些问题相应的泛函尚未找到，或者根本不存在相应的泛函。在这种情况下，就无法用变分法求解。加权余量法（也称加权余值法）是求微分方程近似解的一种有效方法。

第40页/共109页第四十页，共110页。设有微分方程假设有一个满足边界条件和具有一定连续程度的试探函数（其中含有若干待定系数）使余量R与所选的试探函数有关，希望R在某种意义上比较小。通过把余量与加权函数Wi(x)正交化的途径，化为代数方程组，即令第41页/共109页第四十一页，共110页。加权函数比较常用的加权函数有狄拉克函数(x),配置法；幂函数,i=0,1,2,……矩量法；里兹基函数，Galerkin法；精度高，常用余数R,最小二乘法。有限元法中积分在单元内进行。第42页/共109页第四十二页，共110页。4弹性力学问题的FEM单元分析

总体分析边界条件处理方程求解第43页/共109页第四十三页，共110页。4-1单元分析单元分析的任务是推导基本未知量节点位移{}与节点力{F}之间的关系。一般可从节点位移、单元内各点位移(插值)、应变（几何方程）、应力（虎克定律）入手，最终求得节点位移{}与节点力{F}之间的关系第44页/共109页第四十四页，共110页。平面三角形单元图示一个三角形单元三个节点按逆时针方向的顺序编码为i,j,m。节点坐标分别为（xk,yk）k=i,j,m。第45页/共109页第四十五页，共110页。在弹性力学平面问题中，每个节点有两个位移分量，因此三角形单元的位移向量可写成与此对应的物理量是六个节点力分量单元分析的任务就是推导以上两向量之间的关系，即其中转换矩阵[K]e称单元刚度矩阵。第46页/共109页第四十六页，共110页。位移插值函数

单元（元素）分析的第一步是由单元节点位移分量推算单元内部任意一点的位移。在任一单元内可以通过线性插值构造三角形内部任意一点的位移u和v，即可设在第47页/共109页第四十七页，共110页。代入节点位移可得到方程组，解方程组得到系数值第48页/共109页第四十八页，共110页。其中按ijmi的顺序轮换可得到其它参数。三组常数都是与节点坐标有关的常数.A是三角形元素的面积ak,bk,ck是行列式第k行对应元素的代数余子式。

第49页/共109页第四十九页，共110页。位移插值函数将系数i代入u,v表达式，整理可得式中N是插值基函数,也称形函数

l=i,j,m。第50页/共109页第五十页，共110页。形函数的性质形函数描述了单元内部研究对象例如位移函数的分布情况。在单元任意一点，三个形函数的和等于一（设三节点位移相等，就可以证明）。在单元的边上（如ij边），有一个形函数（如Nm）等于零。在单元节点上第51页/共109页第五十一页，共110页。单元刚度矩阵

（1）由位移求应变几何方程用矩阵表示，B称几何矩阵。由于插值函数是一次函数，在单元内的应变分量都是常量。

第52页/共109页第五十二页，共110页。（2）位移与应力根据虎克定律，其中转换矩阵[S]=[D][B](3)由节点位移求节点力

根据虚功原理，第53页/共109页第五十三页，共110页。单元刚度矩阵的性质

1单元刚度矩阵是6×6对称矩阵第54页/共109页第五十四页，共110页。符号规定刚度矩阵中kpq,p代表节点力编号和方向，q代表位移的编号与方向带—代表水平方向，不带—代表垂直方向。例如，代表m节点力垂直方向分量Vm与i节点水平位移分量ui之间的作用系数。第55页/共109页第五十五页，共110页。2K是奇异矩阵，对应的行列式等于零，不存在逆矩阵。实际无法独立求解从几何关系也可得到这一结论。如果给定节点位移，按照这个关系可以计算惟一的节点力。但是反过来，根据节点力却不能得出节点位移的唯一解。因为单元可以产生任意的刚体位移-不影响节点力。第56页/共109页第五十六页，共110页。4-2整体分析平面弹性力学有限元法计算的最终结果是要求出区域R中的位移分布、应变分布及应力分布。由于已经把解域R划分成有限个三角形单元及若干个节点，并把整个场离散到这些节点上，因此，我们的任务就是把这些节点上的位移、应变和应力求出来。第57页/共109页第五十七页，共110页。整体分析的任务由于一个节点不止属于一个单元，因而要将单元方程加以合并组合，得到反映所有节点节点力与位移关系的总方程组。即根据单元刚度矩阵建立连续体的总体刚度矩阵、根据边界条件建立节点载荷向量。整体分析的任务就是建立整体刚度矩阵和节点载荷向量。第58页/共109页第五十八页，共110页。合成方法与步骤在刚度矩阵合成时，存在两个问题。一是单元矩阵与总体矩阵的阶数不同；二是单元节点是局部号码。刚度矩阵的集成分两步进行。第一步把单元刚度矩阵扩大到n×n阶，然后把单元刚度矩阵中[K]e的9个2×2子块搬家，变成按总体编码的9个子块,其它子块用零来填充。例如某单元局部编码i,j,m对应于总码（节点编号）2,4,5,则该单元矩阵在扩大、搬家之后将变成（设n=7）第59页/共109页第五十九页，共110页。扩维第60页/共109页第六十页，共110页。叠加第二步，把扩充以后的单元矩阵迭加，即得到整体刚度矩阵。为了节省存储容量，以上两步实际上是穿插进行的，即按对号入座、边搬家、边累加的方法进行集成。节点载荷向量的构建方法与刚度矩阵相似。由作用力与反作用力的关系，大多数节点的节点力的合力都等于零，只有那些第一类边界节点或受到体积力的节点才不等于零。第61页/共109页第六十一页，共110页。4-3第一类边界条件的处理

在有限元法中，自然边界条件（包括第二、第三类边界条件）已经包含在变分公式当中，因而是自动满足的。第一类边界条件必须强制性地引入到有限元方程当中，即对有关边界节点的方程进行修改。为了修改方便，总体节点编号一般先编内部节点，最后编第一类边界节点。第62页/共109页第六十二页，共110页。例如节点9的水平位移,只需要将K中第17行主对角元素改为1,其余全为零，对应的载荷向量也改成零。例如节点7的垂直位移,则可进行如下修改。将K中第14行主对角元素乘以一个较大的数M（例如108），对应的载荷向量也改成,其余各项不变。这时，这个方程中除包含大数M的两项外，其余各项相对较小，可以忽略不计，可以得到。第63页/共109页第六十三页，共110页。4-4方程特点与求解

1.整体刚度矩阵的特点（1）整体刚度矩阵是一个对称矩阵。（2）整体刚度矩阵是一个稀疏矩阵，它的大多数元素都是零。（3）整体刚度矩阵中的非零元素分布在以主对角线为中心的斜带型区域内。利用矩阵的对称性和带型分布规律，在计算机中可以只存储上半带元素。

第64页/共109页第六十四页，共110页。2.方程的求解可以用带消去法或高斯—赛德尔迭代法等求解方程组，得到节点位移。然后再用单元分析公式计算每个单元的应变或应力。第65页/共109页第六十五页，共110页。5弹性力学问题有限元解法下面以一个简单的模型，示范弹性力学平面问题有限元解法的基本程序。如图所示的试件，受到一组集中力。

第66页/共109页第六十六页，共110页。5-1单元分析第一单元，节点坐标分别是i=1(0,20),j=2(0,10),m=3(10,10)（逆时针方向编号），面积A=50

第67页/共109页第六十七页，共110页。弹性系数设为平面应力问题，泊松比0.28,E=21000，板厚t=10第68页/共109页第六十八页，共110页。1，2，4单元刚度矩阵第69页/共109页第六十九页，共110页。对于第二单元，节点坐标分别是i=2(0,10),j=4(0,0),m=5(10,0),A=50B及K与第一、四单元相同第70页/共109页第七十页，共110页。第三单元，节点坐标分别是i=2(0,10),j=5(10,0),m=3(10,10)A=50

第71页/共109页第七十一页，共110页。单元3刚度矩阵第72页/共109页第七十二页，共110页。5-2整体分析第三单元刚度矩阵扩容后的形式第73页/共109页第七十三页，共110页。整体刚度矩阵第74页/共109页第七十四页，共110页。强加第一类边界条件后刚度矩阵第75页/共109页第七十五页，共110页。5-3平衡方程及其解其中，第76页/共109页第七十六页，共110页。节点位移解方程，可以得到节点位移。第77页/共109页第七十七页，共110页。在单元内利用位移与应变的关系（几何方程）可得到节点应变；根据虎克定律可以计算节点应力。例如对于第一单元，根据节点位移，计算123三节点的应变和应力

第78页/共109页第七十八页，共110页。节点应变第79页/共109页第七十九页，共110页。节点应力第80页/共109页第八十页，共110页。6二维热传导问题FEM解法

第81页/共109页第八十一页，共110页。插值函数中的系数A为三角形单元的面积

第82页/共109页第八十二页，共110页。因为传热问题没有对应的泛函，以下我们采用加权余量法构造单元方程。

应用加辽金法，令

第83页/共109页第八十三页，共110页。因为所取的插值函数为一次函数，首先应用分部积分公式，把二阶导数项化为一阶导数S是单元e的周界。

第84页/共109页第八十四页，共110页。如果单元周界外法线n与x轴的夹角余弦为nx,与y轴的夹角余弦为ny,上式中沿单元周界的积分项可以写成，第85页/共109页第八十五页，共110页。此时，微分方程变成如下形式如果温度插值函数使得温度在每个单元之间连续的话，则经过单元集合，相邻边界的积分必将互相抵消，上式中的单元积分仅在边界单元才有作用。第86页/共109页第八十六页，共110页。单元温度分布第87页/共109页第八十七页，共110页。单元方程对于内部单元，不必计算曲线积分对于热流边界单元第88页/共109页第八十八页，共110页。换热边界单元方程第89页/共109页第八十九页，共110页。规定记号第90页/共109页第九十页，共110页。6-2单元分析对于常用的三角形线性单元，由于

l=i,j,m

系数都是节点坐标的函数。一阶导数

均为常数，可以提到积分号之前第91页/共109页第九十一页，共110页。单元热传导矩阵第92页/共109页第九十二页，共110页。单元热交换矩阵根据里兹基函数的性质第93页/共109页第九十三页，共110页。单元热容矩阵第94页/共109页第九十四页，共110页。向量第95页/共109页第九十五页，共110页。单元方程对于内部单元对于热流边界单元对于换热边界第96页/共109页第九十六页，共110页。整体传热方程把所有单元集合可求得总体方程其中第97页/共109页第九十七页，共110页。6-3时间域的离散不稳定温度场计算含有是一个常系数微分方程，可以用差分方法求得数值解。

第98页/共109页第九十八页，共110页。

向后差分格式

将代入得到第99页/共109页第九十九页，共110页。C-N格式

第100页/共109页第一百页，共110页。加辽金格式

第101页/共109页第一百零一页，共110页。两点格式通式第102页/共109页第一百零二页，共110页。三点向后差分公式

式中a是变步长系数。如果a=1，步长固定。经过时间域的离散化，最后可将问题归结为求解代数方程组问题。

第103页/共109页第一百零三页，共110页。7FEM与FDM的比较-处理方法有限元素法是基于变分原理,同时吸收了差分法中的区域划分的思想。它将所要求解的物理问题化为泛函求极值的问题,再通过有限元素划分来构造插值函数。将连续问题离散化而建立起计算格式。因而有限元素法在处理物理问题