有限元思路框图.ppt

《有限元思路框图.ppt》由会员分享，可在线阅读，更多相关《有限元思路框图.ppt（86页珍藏版）》请在三一办公上搜索。

1、有限元思路框图,（1）剖分结构时应对单元、节点分别用连续正整数编号。,（2）从结构中取出单元，进行单元分析,杆件单元,板单元,第二章 单元分析 平面问题常应变单元,在用矩阵描述单元各种力学量时，不同性质单元的同一力学量可采用相同的矩阵符号，不同的仅仅是矩阵体积和矩阵元素。,本章主要讲单元分析的一般理论、方法。但为了便于理解，以平面问题常应变三角形单元为对象进行说明、演引。必须指出：尽管说明、演引中具有明显的针对性（平面问题三角形单元），但原理、方法和主要矩阵公式都具有普遍性。,单 元 分 析 的 内 容,结点位移,(1),单元内部各点位移,单元应变,单元应力,(2),(3),结点力,(4),位

2、移协调模式,几何方程,物理方程,平衡方程边界条件,单元分析,单元刚度矩阵,（2-1）,2、单元内任意点的体积力列阵qV,（2-2）,1、单元表面或边界上任意点的表面力列阵qs,2.1 基本力学量矩阵,3、单元内任意点的位移列阵f,（2-3）,4、单元内任意点的应变列阵,（2-4）,5、单元内任意点的应力列阵,（2-5）,6、几何方程列阵,（2-6）,将上式代入式（2-4）,7、物理方程矩阵式,（2-7）,式中 E、弹性模量、泊松比。,上式可简写为,（2-8）,对于弹性力学的平面应力问题，物理方程的矩阵形式可表示为：,（2-9）,矩阵D称为弹性矩阵。式（2-9）给出的弹性矩阵D的矩阵元素是按照平

3、面应力问题的物理方程得出的；对于平面应变问题，需将式（2-9）中的 E 换为，换为。,各种类型结构的弹性物理方程都可用式（2-8）描述。但结构类型不同，力学性态(应力分量、应变分量)有区别，弹性矩阵D的体积和元素是不同的。,其中:,单 元 分 析 的 内 容,结点位移,(1),单元内部各点位移,单元应变,单元应力,(2),(3),结点力,(4),位移协调模式,几何方程,物理方程,平衡方程边界条件,单元分析,单元刚度矩阵,？,2.2 位移函数和形函数,1、位移函数概念“位移函数”也称“位移模式”，是单元内部位移变化的数学表达式，是坐标的函数。有限元法采用能量原理进行单元分析，因而必须事先给出（设

4、定）位移函数。一般而论，位移函数选取会影响甚至严重影响计算结果的精度。弹性力学中，恰当选取位移函数不是一件容易的事情。有限单元法中当单元划分得足够小时，把位移函数设定为简单的多项式也可得到相当精确的结果。这正是有限单元法具有的重要优势之一。,不同类型结构会有不同的位移函数。这里，仍以平面问题三角形单元（图2-2）为例，说明设定位移函数的有关问题。,图2-2是一个三节点三角形单元，其节点i、j、m按逆时针方向排列。每个节点位移在单元平面内有两个分量：,（2-10）,2、位移函数设定举例,一个三角形单元有3个节点（以 i、j、m为 序），共有6个节点位移分量。其单元位移或单元节点位移列阵为：,（2

5、-11）,本问题选位移函数为：,（2-12）,式中:,a1、a2、a6待定常数，由单元位移的6个分量确定。,式（2-12）位移函数中，a1、a4代表刚体位移，a2、a3、a5、a6 代表单元中有常应变，且位移函数是连续函数。,3、选取位移函数应考虑的问题,（1）单元有几个位移函数 单元中任意一点有几个位移分量就有几个位移函数。本单元中有u和v，与此相应，有2个位移函数；,（3）位移函数中待定常数个数 待定常数个数应等于单元位移列阵中的位移分量数。以便用单元位移确定位移函数中的待定常数。本单元位移列阵中有6个分量，为了能把2个位移函数（u、v）和单元位移的6个分量联系起来，两个位移函数中包含的待

6、定常数一共应有6个。,（2）位移函数是坐标的函数 本单元的坐标系为：X、Y；,（4）位移函数中必须包含单元的刚体位移。,（5）位移函数中必须包含单元的常应变。,（6）位移函数在单元内要连续;相邻单元间要尽量协调。,条件（4）、（5）构成单元的完备性准则，条件（6）是单元的协调性条件。理论和实践都已证明:完备性准则是有限元解收敛于真实解的必要条件，再加上位移协调条件(充分条件)才构成有限元解的充要条件。容易证明，三角形三节点常应变单元满足以上必要与充分条件。,例：平面应力矩形板被划分为若干三角形单元。,位移函数中包含了单元的常应变。,(a2,a6,a3+a5),位移函数中包含了单元的刚体位移:,

7、对任一单元，如单元，取位移函数：,、单元的位移函数都是,可以看出：位移函数在单元内是连续的；位移函数在单元之间的边界上也连续吗？是。,以、的边界26为例:,两条直线上有两个点重合，此两条直线必全重合。,4、形函数,形函数是用单元位移分量来描述位移函数的插值函数。,（1）形函数定义,现在，通过单元位移确定位移函数中的待定常数a1、a2、a6。设节点i、j、m的坐标分别为（xi、yi）、（xj、yj）、（xm、ym），节点位移分别为（ui、vi）、（uj、vj）、（um、vm）。将它们代入式（2-12），有:,（2-13）,从式（2-13）左边3个方程中解出待定系数a1、a2、a3为:,（2-14

8、）,式中，A为三角形单元的面积，有:,（2-15）,特别指出：为使求得面积的值为正值，本单元节点号的次序必须是逆时针转向，如图所示。至于将哪个结点作为起始结点i，则没有关系。,将式（2-14）代入式（2-12）的第一式，整理后得,同理:,（2-16）,式中:,式（2-17）中（i、j、m）意指：按i、j、m依次轮换下标，可得到aj、bj、cjam、bm、cm。后面出现类似情况时，照此推理。式（2-17）表明：ai、bi、ciam、bm、cm是单元三个节点坐标的函数。,（2-16）,令,（2-18）,位移模式（2-16）可以简写为,（2-19）,式（2-19）中的Ni、Nj、Nm是坐标的函数，反

9、应了单元的位移形态，称为单元位移函数的形函数。数学上它反应了节点位移对单元内任一点位移的插值，又称插值函数。,（2-16）,用形函数把式（2-16）写成矩阵，有,缩写为,（2-20）,形函数是用单元位移分量来描述位移函数的插值函数。,N为形函数矩阵，进一步写成分块形式：,（2-21）,其中子矩阵,（2-22）,I是22的单位矩阵。,下面将会看到，形函数是有限单元法中的一个重要函数。了解它的一些基本性质是有益的。,（2）形函数性质,图2-3,性质1 形函数Ni在节点i上的值等于1，在其它节点上的值等于0。对于本单元，有:,（i、j、m）,性质2 在单元中任一点，所有形函数之和等于1。对于本单元，

10、有,为什么？,图2-4,形 函 数,i,j,m,在单元任一点上三个形函数之和等于1,第一列与它的代数余子式之和,第一列与第二列的代数余子式之和,第一列与第三列的代数余子式之和,1.三个形函数只有两个是独立的,2 当三角形单元的三个结点的位移相等,图2-4,上图说明：形函数是线性的（1）单元位移场是线性的；（2）单元位移场与结点位移是协调的；（3）结点位移将影响位移场的数值大小,性质3 在三角形单元的边界ij上任一点（x，y），有:,证,图2-5,性质4 形函数在单元上的面积分和边界上的线积分公式为,（2-23）,式中 为 边的长度。,相邻单元的位移在公共边上是连续的,i,j,p,m,i(xi，

11、yi),j(xj，yj),m(xm，ym),2.3 单元应变矩阵和应力矩阵,式（2-6）给出了单元内任一点的应变和位移之间关系。,（2-6）,1、单元应变矩阵,对位移函数（式（2-16）,（2-24）,（2-16）,求导后代入式（2-6），得到应变和节点位移的关系式。,（2-25）,式中,B单元应变矩阵。,对本问题，维数为36。它的分块形式为：,子矩阵:,（2-26）,由于 与x、y无关，都是常量，因此B矩阵也是常量。单元中任一点的应变分量是B矩阵与单元节点位移的乘积，因而也都是常量。因此，这种单元被称为常应变单元。,2、单元应力矩阵,将式（2-25）代入物理方程式（2-8），得,（2-8）,

12、（2-27）,上式也可写为：,（2-28）,这是单元内任一点应力与单元位移的关系式。其中S称为单元应力矩阵，并有：,（2-29）,D是33 弹性矩阵，B是36应变矩阵，因此S也是36 矩阵。它可写为分块形式,（2-30）,将弹性矩阵（式（2-9）和应变矩阵（式（2-26）代入，得子矩阵Si,由式（2-29）得：,（2-31）,式（2-31）是平面应力问题的结果。对于平面应变问题，只要将上式中的E换成，换成 即得。,（2-32）,由于同一单元中的D、B矩阵都是常数矩阵，所以S矩阵也是常数矩阵。也就是说，三角形三节点单元内的应力分量也是常量。当然，相邻单元的E，A和bi、ci(i，j，m)一般不完

13、全相同，因而具有不同的应力，这就造成在相邻单元的公共边上存在着应力突变现象。但是随着网格的细分，这种突变将会迅速减小，平衡被满足。,几何关系位移函数,几何关系,?,平衡关系,单元刚度矩阵,2.4 单元应变能和外力势能的矩阵表达,1、单元应变能,仍以平面应力问题中的三角形单元说明，设单元厚度为h,将式（2-25）和（2-27）代入上式进行矩阵运算，并注意到弹性矩阵D的对称性，有,应变能 U为：,由于和T是常量，提到积分号外，上式可写成,引入矩阵符号k，且有：,（2-33a）,式（2-33a）是针对平面问题三角形单元推出的。注意到其中hdxdy的实质是任意的微体积dv，于是得 k的一般式。,（2-

14、33）,式（2-33）不仅适合于平面问题三角形单元，也是计算各种类型单元k的一般式。,2.6节中将明确k的力学意义是单元刚度矩阵。式（2-33）便是计算单元刚度矩阵的基本矩阵式。它适合于各种类型的单元。,单元应变能写成,（2-34）,2、单元外力势能,单元受到的外力一般包括体积力、表面力和集中力。自重属于体积力范畴。表面力指作用在单元表面的分布载荷，如风力、压力，以及相邻单元互相作用的内力等。,（2-33）,（1）体积力势能,单位体积中的体积力如式（2-2）所示。,单元上体积力具有的势能Vv为,注意到式（2-20）,有,（2）表面力势能,面积力虽然包括单元之间公共边上互相作用的分布力，但它们属

15、于结构内力，成对出现，集合时互相抵消，在结构整体分析时可以不加考虑，因此单元分析时也就不予考虑。,现在，只考虑弹性体边界上的表面力，它只在部分单元上形成表面力（右下图）。设边界单位长度上受到的表面力如式（2-1）。,l单元边界长度h单元厚度A表面力作用面积,qs,则单元表面力的势能Vs为,（3）集中力势能,当结构受到集中力时，通常在划分单元网格时就把集中力的作用点设置为节点。于是单元集中力Pc的势能Vc为,（4）外力总势能,如果把（2-35）式中原括号内的部分用列阵Fd代替，,综合以上诸式，单元外力的总势能V为,Fd具有和相同的行、列数。则：,由单元的应变能U（2-34）和外力势能V（2-36

16、）,可得单元的总势能,（2-37）,以节点位移为未知量，对总势能取极值问题变成了一个多元函数的极值问题。有极值条件,2.5 能量原理和单元平衡方程,（2-36）,式（2-38）是从能量原理导出的单元平衡方程。这个方程表达了单元力与单元位移之间的关系。其中，Fd和单元节点力F具有相同的意义。,（2-38）,于是，将式（2-37）代入，即得单元平衡方程:,根据弹性力学能量原理：结构处于稳定平衡的必要和充分条件是总势能有极小值。,2.6 单元刚度矩阵,平衡方程（2-38）中的矩阵k是单元力和单元位移关系间的系数矩阵，代表了单元的刚度特性，称为单元刚度矩阵。单元刚度矩阵的体积为nj nj，nj 是单元

17、位移总数。,1、计算单元刚度矩阵的一般公式,计算各类单元的单元刚度矩阵可用式（2-33）执行。它与单元应变矩阵B和弹性矩阵D有关。,对于平面应力三角形单元，应变矩阵B是常数矩阵，同时弹性矩阵D也是常数矩阵，于是式（2-33）可以化简为,式中A表示三角形单元的面积。,2、平面问题三角形单元刚度矩阵,（1）平面应力三角形单元,将式（2-9）和（2-26）代入上式，,即得平面应力三角形单元刚度矩阵。写成分块形式，有,（2-40）,式中子矩阵为22矩阵，有,（2-41）,（2）平面应变三角形单元,对于平面应变问题，须将上式中的E换为，换为，于是有：,其中，bi(j,m)、ci(j,m)是形函数式（2-

18、16）中的系数。,（2-42）,平面问题的单元刚度矩阵ke不随单元（或坐标轴）的平行移动或作n角度（n为整数）的转动而改变。这可由kij及bi、bj、ci、cj（i、j、m）计算公式（2-41）、（2-42）、（2-17）的分析中得到。?,图2-7,应当注意的是：当单元旋转时，各节点的编号保持不变。如图2-7所示，图a所示的单元旋转时，到达图b所示位置。,这两种情形的k是相同的。,当坐标系旋转任意角度时，由于公式（2-40）、（2-41）、（2-42）中已经包含坐标系的影响，它们仍然照用。,由此得出重要结论：平面问题的单元刚度矩阵k计算式不因坐标系不同而变。因此，这组公式是规格化的式子。,（3

19、）示例 平面应力直角三角形单元刚度矩阵,图2-8示出一平面应力直角三角形单元，直角边长分别为a、b，厚度为h，弹性模量为E，泊松比为，计算单元刚度矩阵。,第一步：计算bi、ci和单元面积A,图2-8,（1-17）,表2-1 单元节点坐标和bi、ci值（i、j、m）,参数,节点,单元面积:A=ab/2,计算步骤,第二步：求子矩阵(由式（2-41），算得),其他从略。,第三步：形成k 将kii等按式（2-40）组集成k。,（2-40）,式中上边和右边的i、j、m表示子矩阵对应的节点号。,当a=b时，即等腰直角三角形单元，有,子程序框图,子程序,见附页 SUBROUTING SME3(),3、单元刚

20、度矩阵性质,njnj,单元刚度矩阵k的详细内容为（i、j是行列号）：,（2-38）,单元刚度矩阵具有以下的性质：,（1）单元刚度矩阵中每个元素有明确的物理意义例如，kij表示当单元位移中第j个元素为1(j=1)其余元素为零时，引起的单元力中的第i个节点力Fi。,把平衡方程（2-38）写开,主对角线上元素kii(i=1,nj)恒为正值。,（2）k的每一行或每一列元素之和为零,以上式中第i行为例，,当所有节点沿x向或y向都产生单位位移时，单元作平动运动，无应变，也无应力，因而单元结点力为零（不含初应力）。所以有,即，k的每一行元素之和为零。由于对称性，每一列元素之和也为零。,（3）k是对称矩阵 由

21、k单元的表达式，可见，由此可知k具有对称性。,njnj,对于主对角线元素对称。对称表达式：,kij=kji,证明：,kij表示当单元位移中第j个元素为1(j=1)其余元素为零时，引起的单元力中的第i个节点力Fi,kji表示当单元位移中第i个元素为1(i=1)其余元素为零时，引起的单元力中的第j个节点力Fj,由虚功原理，得,kij=kji,（4）单元刚度矩阵是奇异矩阵（即k的行列式为零）单元刚度矩阵是在单元处于平衡状态的前提下得出的。单元作为分离体看待，作用在它上面的外力（单元力）必定是平衡力系。然而，研究单元平衡时没有引入约束。承受平衡力系作用的无约束单元，其变形是确定的，但位移不是确定的。所

22、以出现性质（3）中的“平动问题”，即单元可以发生任意的刚体运动。从数学上讲，方程（2-38）的解不是唯一的或不能确定的。由此，单元刚度矩阵一定是奇异的。,（5）单元刚度矩阵是常量矩阵,单元力和单元位移成线性关系是基于弹性理论的结果。,2.7 等价节点力,从前面单元分析可以看出：单元平衡所用到的量均属于节点的量，如单元位移、单元力。载荷亦应如此，必须将体积力、表面力转化到节点上去，成为等价节点力（载荷）。在第2.5节中已经得到了公式（2-35）和（2-36）。,（2-35）,（2-36）,这里，Fd就是体积力、表面力和集中力之和的总等价节点力。,（2-44）,把总等价节点力 Fd 分解成体积力、

23、表面力和集中力的等价节点力之和，有,FV单元上体积力的等价节点力 FS单元上表面力的等价节点力 pC单元上节点上的集中力,注意到式（2-35），得体积力等价节点力计算公式：P39,表面力的等价节点力计算公式：,（2-45）,（2-46）,1、体积力的等价节点力,2、表面力的等价节点力,3、等价节点力计算举例,（1）单元自重,图2-9所示平面应力三角形单元，单元厚度为h。单元单位体积自重为，自重指向y轴的负方向。,（2-45）,计算式,注意到形函数的性质4：,（2-23）,得自重荷载的等价节点力,根据体积力和式（2-45）、（2-21）、（2-22），得,（2-47）,上式表明：自重载荷的等价节

24、点力为单元重量的1/3。,子程序,见 SUBJPE.FOR,（2）均布面力,单元边界上作用了均匀的分布力，如图2-10所示，其集度为qs。,（2-46）,（2-21）,根据式（2-46）、（2-21）和（2-22）,计算式,注意到形函数性质4：,（2-23）,得,（2-48）,（2-22）,均匀分布力的等价节点力为,式（2-48）表明：在ij边上受均布面力的平面问题三角形单元，其等价节点力等于将均布面力合力之半简单地简化到i、j节点上，方向与分布力方向相同。m节点上为零。,子程序,见 SUBJPE.FOR,（3）线性分布面力,表面力集度在i点为qsx,qsyT，而在j点为0。设坐标轴s的原点取

25、在i点，沿ij为正向，。,ij边上任一点的面力集度qs,计算式,在ij边上有：,将qs和上式代入式（2-46），有,由形函数的性质3,用坐标s表示：,（2-49）,式（2-49）表明：ij边受线性分布面力：i点为qsx,qsyT，j点为0时，其等价节点力可将总载荷的2/3分配给i点，1/3分配给j点，m点为零得出。,子程序,见 SUBJPE.FOR,单元上的体积力和表面力向结点的移置都是符合直观的静力等效原理的，并与工程中的简单的处理方法相一致。应当指出，这种移置方法是线性位移模式三结点三角形单元的必然结果。对于非线性位移模式的单元，上述这种简单的载荷移置方法一般是不成立的，而应当采用公式（2-35）进行计算。,（2-35）,本章小结：单元分析的主要任务是：一、组集单元刚度矩阵；二、组集单元等价节点力矩阵。,下接第3章,