角形单元的有限元法.ppt

第五章 三角形单元的有限元法,5.1 基本思想,把整体结构离散为有限个单元，研究单元的平衡和变形协调；再把这有限个离散单元集合还原成结构，研究离散结构的平衡和变形协调。划分的单元大小和数目根据计算精度和计算机能力来确定。,弹性悬臂板剖分与集合,单元、节点需编号,（1-1）,2、单元内任意点的体积力列阵qV,（1-2）,1、单元表面或边界上任意点的表面力列阵qs,5.2 基本力学量矩阵表示,3、单元内任意点的位移列阵f,（1-3）,4、单元内任意点的应变列阵 ,（1-4）,5、单元内任意点的应力列阵,（1-5）,6、几何方程,（1-6）,将上式代入式(1-4），,7、物理方程矩阵式,（1-7）,式中 E、弹性模量、泊松比。,上式可简写为,（1-8）,其中,对于弹性力学的平面应力问题，物理方程的矩阵形式可表示为：,(1-9）,矩阵D称为弹性矩阵。对于平面应变问题，将式(1-9）中的E换为 ，换为 。,5.3 位移函数和形函数,1、位移函数概念 由于有限元法采用能量原理进行单元分析，因而必须事先设定位移函数。 “位移函数”也称 “位移模式”，是单元内部位移变化的数学表达式，设为坐标的函数。 一般而论，位移函数选取会影响甚至严重影响计算结果的精度。在弹性力学中，恰当选取位移函数不是一件容易的事情；但在有限元中，当单元划分得足够小时，把位移函数设定为简单的多项式就可以获得相当好的精确度。这正是有限单元法具有的重要优势之一。,不同类型结构会有不同的位移函数。这里，仍以平面问题三角形单元（图1-2）为例，说明设定位移函数的有关问题。,图1-2是一个三节点三角形单元，其节点i、j、m按逆时针方向排列。每个节点位移在单元平面内有两个分量：,（1-10）,一个三角形单元有3个节点（以 i、j、m为 序），共有6个节点位移分量。其单元位移或单元节点位移列阵为：,2、位移函数设定,本问题选位移函数（单元中任意一点的位移与节点位移的关系）为简单多项式：,(1-12）,式中：a1、a2、a6待定常数，由单元位移的 6个分量确定。a1、a4代表刚体位移，a2、 a3 、 a5 、 a6 代表单元中的常应变，而且，位移函数是连续函数。,(1-11）,选取位移函数应考虑的问题,（1）位移函数的个数 等于单元中任意一点的位移分量个数。本单元中有u和v，与此相应，有2个位移函数；,（3）位移函数中待定常数个数 待定常数个数应等于单元节点自由度总数，以便用单元节点位移确定位移函数中的待定常数。本单元有6个节点自由度，两个位移函数中共包含6个待定常数。,（2）位移函数是坐标的函数 本单元的坐标系为：x、y；,（4）位移函数中必须包含单元的刚体位移。,（5）位移函数中必须包含单元的常应变。,（6）位移函数在单元内要连续。相邻单元间要尽 量协调。,条件（4）、（5）构成单元的完备性准则。 条件（6）是单元的位移协调性条件。 理论和实践都已证明，完备性准则是有限元解收敛于真实解的必要条件。单元的位移协调条件构成有限元解收敛于真实解的充分条件。 容易证明，三角形三节点常应变单元满足以上必要与充分条件。,（7）位移函数的形式 一般选为完全多项式。为实现（4）（6）的要求，根据Pascal三角形由低阶到高阶按顺序、对称地选取；多项式的项数等于（或稍大于）单元节点自由度数。,例：平面应力矩形板被划分为若干三角形单元。,对任一单元，如单元，取位移函数：,、单元的位移函数都是,可以看出： 位移函数在单元内是连续的；,以、的边界26为例,两条直线上有两个点重合，此两条直线必全重合。,位移函数在单元之间的边界上也连续吗？是。,3、形函数,形函数是用单元节点位移分量来描述位移函数的插值函数。,(1-13）,（1）形函数确定,现在，通过单元节点位移确定位移函数中的待定常数a1、a2、a6 。设节点i、j、m的坐标分别为（xi、yi）、（ xj、yj ）、（ xm、ym ），节点位移分别为（ui、vi）、 （uj、vj） 、 （um、vm）。将它们代入式（1-12），有,从式(1-13）左边3个方程中解出待定系数a1、a2、a3为,(1-14）,式中， A为三角形单元的面积，有,(1-15）,特别指出：为使求得面积的值为正值，本单元节点号的次序必须是逆时针转向，如图所示。至于将哪个节点作为起始节点i，则没有关系。,将式(1-14）代入式(1-12）的第一式，整理后得,同理,(1-16）,式中,(1-16）,令,(1-18）,位移模式(1-16）可以简写为,(1-19）,式(1-19）中的Ni、Nj、Nm是坐标的函数，反应了单元的位移形态，称为单元位移函数的形函数。数学上它反应了节点位移对单元内任一点位移的插值，又称插值函数。,(1-16）,用形函数把式(1-16）写成矩阵，有,缩写为,(1-20）,形函数是有限单元法中的一个重要函数，它具有以下性质：,N为形函数矩阵，写成分块形式：,(1-21）,其中子矩阵,(1-22）,I是22的单位矩阵。,（2）形函数性质,性质1 形函数Ni在节点i上的值等于1，在其它节点 上的值等于0。对于本单元，有,（i、j、m）,性质2 在单元中任一点，所有形函数之和等于1。对 于本单元，有,图1-3,？,图1-4,也可利用行列式代数余子式与某行或列元素乘积的性质（等于行列式值或0）证明。,性质3 在三角形单元的边界ij上任一点（x，y），有,证,图1-5,（1）,性质4 形函数在单元上的面积分和在边界上的线积分公式为,(1-23）,式中 为 边的长度。,5.4 单元应变和应力,根据几何方程(1-6）和位移函数(1-16）可以求得单元应变。,1、单元应变,对位移函数（式(1-16）,(1-24）,(1-16）,求导后代入式(1-6），得到应变和节点位移的关系式。,上式简写一般式：,(1-25）,式中， B单元应变矩阵。,对本问题，维数为36。它的分块形式为：,子矩阵,(1-26）,由于 与x、y无关，都是常量，因此B矩阵也是常量。单元中任一点的应变分量是B矩阵与单元位移的乘积，因而也都是常量。因此，这种单元被称为常应变单元。,2、单元应力,将式(1-25）代入物理方程式(1-8），得 单元应力,(1-27）,也可写为,(1-28）,其中：S称为单元应力矩阵，并有,(1-29）,这里，D是33矩阵，B是36矩阵，因此S也是36矩阵。它可写为分块形式,(1-30）,将弹性矩阵（式(1-9） 和应变矩阵（式(1-26）代入，得子矩阵Si,由式(1-29）,(1-31）,式(1-31）是平面应力的结果。对于平面应变问题，只要将上式中的E换成 ，换成 即得。,(1-32）,由于同一单元中的D、B矩阵都是常数矩阵，所以S矩阵也是常数矩阵。也就是说，三角形三节点单元内的应力分量也是常量。 当然，相邻单元的bi、ci(i，j，m)一般不完全相同，因而具有不同的应力，这就造成在相邻单元的公共边上存在着应力突变现象。但是随着网格的细分，这种突变将会迅速减小，收敛于平衡被满足。,5.5 单元平衡方程,1、 单元应变能,对于平面应力问题中的三角形单元，设单元厚度为h 。,将式(1-25）和(1-8）代入上式进行矩阵运算，并注意到弹性矩阵D的对称性，有,应变能 U为,由于和T是常量，提到积分号外，上式可写成,引入矩阵符号k，且有,(1-33a）,式(1-33a）是针对平面问题三角形单元推出的。注意到其中hdxdy的实质是任意的微体积dv，于是得计算k的一般式。,(1-33）,式(1-33）不仅适合于平面问题三角形单元，也是计算各种类型单元k的一般式。,5.6节中将明确k的力学意义是单元刚度矩阵。式(1-33）便是计算单元刚度矩阵的基本矩阵式。它适合于各种类型的单元。,单元应变能写成,(1-34）,2、 单元外力势能,单元受到的外力一般包括体积力、表面力和集中力。自重属于体积力范畴。表面力指作用在单元表面的分布载荷，如风力、压力，以及相邻单元互相作用的内力等。,(1-33）,（1） 体积力势能,单位体积中的体积力如式(1-35）所示。,单元上体积力具有的势能Vv为,注意到式(1-20）,有,（2） 表面力势能,面积力虽然包括单元之间公共边上互相作用的分布力，但它们属于结构内力，成对出现，集合时互相抵消，在结构整体分析时可以不加考虑，因此单元分析时也就不予考虑。,现在，只考虑弹性体边界上的表面力，它只在部分单元上形成表面力（右下图）。设边界面上单位面积受到的表面力如下式：,l单元边界长度 h单元厚度 A表面力作用面积,qs,qs沿厚度均匀分布，则单元表面力的势能Vs为,（3） 集中力势能,当结构受到集中力时，通常在划分单元网格时就把集中力的作用点设置为节点。于是单元集中力Pc的势能Vc为,（4）总势能,把(1-35）式中原括符内的部分用列阵Fd代替，,综合以上诸式，单元外力的总势能V为,(1-35）,Fd具有和相同的行、列数。则,(1-36）,由单元的应变能U(1-34）和外力势能V(1-36），可得单元的总势能,(1-37）,将式(1-37）代入，,根据弹性力学最小势能原理：结构处于稳定平衡的必要和充分条件是总势能有极小值。,3、单元平衡方程,于是有，,式(1-38）是从能量原理导出的单元平衡方程。这个方程表达了单元力与单元位移之间的关系。其中，Fd和单元节点力F具有相同的意义。,(1-38）,即得单元平衡方程,5.6 单元刚度矩阵,平衡方程(1-38）中的矩阵k是单元力和单元位移关系间的系数矩阵，代表了单元的刚度特性，称为单元刚度矩阵。单元刚度矩阵的体积为nj nj， nj 是单元位移总数。其一般计算公式为：,1、一般计算公式,它与单元应变矩阵B和弹性矩阵D有关。,对于平面应力三角形单元，应变矩阵B是常数矩阵，同时弹性矩阵D也是常数矩阵，于是式(1-33）可以化简为,式中A表示三角形单元的面积。h是单元厚度。,2、平面问题三角形单元刚度矩阵,（1）平面应力三角形单元,将式(1-9）和(1-26）代入上式，,即得平面应力三角形单元刚度矩阵。写成分块形式，有,(1-40）,式(1-40)中子矩阵krs为22矩阵，有,(1-41）,（2）平面应变三角形单元,对于平面应变问题，须将上式中的E换为 ， 换为 ，于是有,其中，bi(j,m)、ci(j,m)是形函数式(1-16）中的系数(式2-17)。,(1-42）,平面问题的单元刚度矩阵k不随单元（或坐标轴）的平行移动或作n角度（n为整数）的转动而改变。 由公式(1-41）、(1-42）知，krs矩阵和其中的br、cr 、 bs、cs （r、s=i、j、m）有关。 单元平移时， bi、ci不变。,(3) 三角形单元刚度矩阵与坐标系无关, 单元转动时， bi、ci不变。 当单元旋转时，各节点的编号保持不变。如图1-7所示，图a所示的单元旋转时，到达图b所示位置。,可以证明，这两种情形的k是相同的。 其实，推演公式(1-40）、(1-41）、(1-42）时并没有规定坐标系的方位，当坐标系旋转任意角度时，也不影响刚度矩阵的结果。因此，平面问题的单元刚度矩阵可以认为是结构坐标系中的单元刚度矩阵，没有坐标变换问题。,(1-38）,（1）单元刚度矩阵中每个元素有明确的物理意义例如，kij表示单元第j个自由度产生单位位移(j=1)，其他自由度固定（=0）时，在第i个自由度产生的节点力Fi。,主对角线上元素kii(i=1,nj)恒为正值。,3、单元刚度矩阵性质,（2）k的每一行或每一列元素之和为零,以上式中第i行为例，,当所有节点沿x向或y向 都产生单位位移时，,单元作平动运动，无应变，也无应力。则有：,即：k的每一行元素之和为零。根据对称性，每一列元素之和也为零。,（3）k是对称矩阵 由k各元素的表达式，可知k具有对称性。,njnj,对于主对角线元素对称。对称表达式：,kij = kji,证明, kij表示当单元位移中第j个元素为1(j=1)其余元素为零时，引起的单元力中的第i个节点力Fi, kji表示当单元位移中第i个元素为1(i=1)其余元素为零时，引起的单元力中的第j个节点力Fj,由虚功原理，得,kij = kji,（4）单元刚度矩阵是奇异矩阵 即k的行列式为零（由行列式性质） 。 单元刚度矩阵是在单元处于平衡状态的前提下得出的。单元作为分离体看待，作用在它上面的外力（单元力）必定是平衡力系。然而，研究单元平衡时没有引入约束。承受平衡力系作用的无约束单元，其变形是确定的，但位移不是确定的。所以出现性质（3）中的“平动问题”，即单元可以发生任意的刚体运动。从数学上讲，方程(1-28）的解不是唯一的或不能确定的。由此，单元刚度矩阵一定是奇异的。,（5）单元刚度矩阵是常量矩阵,单元力和单元位移成线性关系是基于弹性理论的结果。,4、例：平面应力直角三角形单元刚度矩阵,图1-8示出一平面应力直角三角形单元，直角边长分别为a、b，厚度为h，弹性模量为E，泊松比为，计算单元刚度矩阵。,图1-8,第一步：计算bi、ci和单元 面积A。,图1-8,（1-17）,表2-1 单元节点坐标和bi、ci值（i、j、m）,参数,节点,单元面积: A=ab/2, 计算步骤,第二步：求子矩阵 由式(1-41），算得,其他从略。,第三步：形成k 将kii等按式(1-40）组集成k 。,(1-43a）,2i-1 2i 2j-1 2j 2m-1 2m,红色号码是单元位移（1、 2、）在结构中对应的节点位移的序号。,i,j,m,i,j,m,i、j、m表示单元中3个节点在结构系统中的编号。,当a=b时，即等腰直角三角形单元，有,(1-43b）,i j m,i,j,m,5.7 等价节点力,从前面单元分析可以看出：单元平衡所用到的的量均要属于节点的量，如单元位移、单元力。载荷亦应如此，必须将体积力、表面力转化到节点上去，成为等价节点力（载荷）。在第2.5节中已经得到了公式(1-35）和(1-36） 。,这里，Fd就是体积力、表面力和集中力之和的总等价节点力。,(1-44）,把总等价节点力 Fd 分解成体积力、表面力和集中力的等价节点力之和，有,FV单元上体积力的等价节点力 FS单元上表面力的等价节点力 pC单元上节点上的集中力,注意到式(1-35），得体积力等价节点力计算公式：,表面力的等价节点力计算公式：,(1-45）,(1-46）,1、体积力的等价节点力,2、表面力的等价节点力,3、等价节点力计算举例,（1）单元自重,图1-9所示平面应力三角形单元，单元厚度为h。单元单位体积自重为，自重指向y轴的负方向。,(1-45）, 计算式,注意到形函数的性质4：,(1-23）,得自重荷载的等价节点力,根据体积力和式(1-45）、(1-21）、(1-22），得,(1-47）,上式表明：自重载荷的等价节点力为单元重量的1/3。,（2）均布面力,单元边界上作用了均匀的分布力，如图1-10所示，其集度为qs。,(1-46）,(1-21）,根据式(1-46）、(1-21）和(1-22）, 计算式,注意到形函数性质4 ：,(1-23）,得,(1-48）,(1-22）,均匀分布力的等价节点力为,式(1-48）表明：在ij边上受均布面力的平面问题三角形单元，其等价节点力等于将均布面力合力之半简单地简化到i、j节点上，方向与分布力方向相同。m节点上为零。,(1-48）,（3）线性分布面力,表面力集度在i点为qsx qsyT，而在j点为0。设坐标轴s的原点取在j点，沿ji为正向， 。,ij边上任一点的面力集度qs,在ij边上有：,将qs和上式代入式(1-46），有,由形函数的性质3：,(1-49）,式(1-49）表明：ij边受线性分布面力： i点为qsx, qsyT，j点为0 时，其等价节点力可将总载荷的2/3分配给i点，1/3分配给j点，m点为零得出。,体积力和表面力向节点的移置符合静力等效原理的前提条件是：线性位移模式。,5.7 系统分析,5.7.1 坐标系,研究各离散单元集合成整体结构，集合整体结构的平衡和变形协调，建立整体结构平衡方程。,单元分析时采用的坐标系成为局部坐标或单元坐标（单元刚度矩阵的通用性）。而结构系统分析时，必须在统一的坐标系内进行（各力学量才能叠加），称为“结构坐标”或“整体坐标”，如图1-13所示。,单元坐标系下，单元位移、单元力、单元刚度矩阵表示为：,整体坐标系下，单元位移、单元力、单元刚度矩阵表示为：,如何从单元坐标转化为结构坐标将在第4章中讨论。,5.7.2 整体刚度矩阵,假设整体结构被划分为ne个单元和n个节点，在整体坐标系下，对于每个单元均有：,将上述这些方程集合起来（整体坐标下叠加），便可得到整个结构的平衡方程。为此，需要将k、F体积膨胀，分别扩大为n1n1、n11和n11的矩阵才能相加。膨胀后，原有节点号对应位置的元素不变，而其它元素均为零。,组装方法：建立一个体积为n1n1的方阵，按单元序号依次把结构坐标单元刚度矩阵的元素放入该方阵中。 放入方法：（1）按单元节点编码对号入座； （2）同位置元素累加。,式中：K为整体刚度矩阵，为整体节点位移列阵；P为整体等价节点荷载列阵。如下：,（1-50）,i,j,m,i,j,m,例：平面三角单元,双行双列,5.7.3 结构刚度矩阵特性,1、结构刚度矩阵元素的力学意义,把方程(1-50）写开，,=1,(1-51）,2、结构刚度矩阵是对称矩阵 已知单元刚度矩阵是对称矩阵（5.7节），用单元刚度矩阵组集 结构刚度矩阵的 过程又没有破坏 其对称性，结构 刚度矩阵必然也 是对称的。当然， 对称性也可以通 过虚功原理得到证明。,结构刚度矩阵中的任一元素kij是j为单位位移（ j =1），其它位移为零时的Pi。,3、结构刚度矩阵主对角线上的元素恒为正值 由性质(1)可知，任一主对角线上元素kii是使节点位移i为一单位位移，其它节点位移为零时必须在第i号位移方向施加的力Pi。它的方向自然应与位移方向相同，因而是正值。,4、结构刚度矩阵是一个稀疏矩阵,稀疏矩阵指：存在大量零元素。非零元素稀疏排列。,矩阵的每一列都有很多零元素。考察矩阵中第j列。,再分析图（1-14）。设节点b发生单位位移j=1，其它位移为零时， j只能在与点节b有直接联系的 q 、 r节点引起节点力，不能 在其它节点引起节点 力。所以式（1-52） 中，只有和q、p、r、 b节点位移的相关元素 才不为零，其余的元素 都是零元素。,任一元素kij是j=1（其它=0）引起的Pi（i=1、2）,（1-52）,b,其它各列的情况也是类似的。 结构的节点总数通常都比直接环绕于任何一个节点的节点数大得多，因而，结构刚度矩阵中很大一部分元素是零，即所谓的稀疏矩阵。,5、结构刚度矩阵是一个奇异矩阵,从单元刚度矩阵的奇异性讨论中知，处于静力平衡状态的无约束单元可以发生任意的刚体位移。与单元刚度矩阵是奇异矩阵的理由一样，无约束结构的结构刚度矩阵K也是奇异矩阵，即K的行列式为零。,5.7.6 引入支承约束的结构节点平衡方程,6、结构刚度矩阵是常量矩阵,结构刚度矩阵是常量矩阵。结构的节点力和节点位移成线性关系都是基于弹性理论的结果。,（1-53）,用平衡方程（1-53）是解不出结构的节点位移的，因为结构刚度矩阵是奇异矩阵。因此，必须引入约束，排除任何刚体位移，使结构为几何不变体系。,方程（1-53）中的刚度矩阵K和节点荷载向量列阵P可分割为约束和自由两部分：,式中，Pr是支承反力，约束位移,自由,约束,(1-55）,(1-56）,展开（154），有：,Kff引入约束后的结构刚度矩阵。它通对K引入约束后获得，具体方法: 从无约束的结构刚度矩阵K中删去与受约束位移号对应的行和列，再将矩阵压缩排列成nn阶方阵，即为约化后的结构刚度矩阵Kff 。 Kff这是一个非奇异矩阵，它存在逆矩阵。,方程(1-55）是引入约束后的结构节点平衡方程，用于计算结构所有非刚性约束节点的节点位移。而方程(1-60）可以用来计算结构所有受刚性约束节点的反力。,(1-61）,由式(1-55）即可解出全部未知的节点位移：,5.7.7 节点位移和单元力的解答,1、结构节点位移,2、支座反力,把解出的f代入(1-56），即得支座反力Pr：,关于方程 (1-61）的解算方法，当Kff采用本章中上述方法组集时，可直接采用结构力学中的高斯（Gauss）法求解。,(1-56）,至此，我们可以看出：系统分析的主要任务是： （1）组集引入约束后的结构刚度矩阵Kff； （2）求解式(1-55）给出的线性代数方程组。算出全部未知的节点位移。,至于支座反力的计算，实际计算时，根本不去组集式(1-56）中的矩阵Krf，不用式(1-62）而是直接对支承点使用节点平衡方程计算。,3、单元力,(1-62）,（1）全解公式,全解是指结构在未经简化的实际荷载作用下（全解系统）的解答。,余解是结构在节点荷载作用下（余解系统）的解答。节点荷载包括原来给定的节点荷载和由分布荷载简化而来的等价节点荷载。,特解是单元从实际分布荷载到等价节点力过程中（特解系统）导致的单元节点力。,通过式(1-55）解出获得全部未知节点位移，据此找出任何一个单元的单元位移（单元节点位移），然后可以求得单元力。,图(1-15）,=,+,（b）特解系统,（c）余解系统,（2）余解,（3）特解,把单元等价节点力反向作用在单元的相应节点，即为单元力的特解。,4、结论,（1）全解的节点位移等于余解的节点位移 因为，特解给出的节点位移总是为零。,（2）对不承受分布荷载的单元，内力由余解直接给出。,（3）对承受分布荷载作用的单元，其内力由余解和特解的迭加给出。,F全e = F余e Fqe (1-64）,等价节点力,3.9 举例,如图1-17a所示两端固支的矩形深梁，跨度为2a，梁高为a，截面宽度为h，已知， =0、E，承受均布压力q。试用有限元法解此平面应力问题。,图1-17,(a),利用对称性，可取梁的一半分析，共分2个单元，4个结点，支座约束如(a)图。,解 1、划分单元。,（C)自由度,等腰直角三角形单元dde k,2、单元刚度矩阵,由书公式（1-46b）,3 4 5 6 1 2,1,3,5 6 3 4 7 8,3,3,3、结构刚度矩阵,引入约束前，结构刚度矩阵是88阶的。受约束的自由度号为1，3，4，5，7，8。删除相应行和列后，结构刚度矩阵剩下2行2列，相应的自由度号是：2，6。,3/2,-1,1,-1,1/2,+,4、结构节点荷载向量,一般情况节点荷载向量有8个元素。引入约束后剩下2个，相应序号为：2，6。,5、结构无约束节点平衡方程,解得,根据结点位移找单元位移 ，代入（1-38）式可求各单元的单元力F 。,（1-38）,利用单元力F计算支点反力。,将单元位移代入式（1-28），其中S可由式（1-30）、（1-31）求得，从而可求出单元应力。,（1-28）,习题：算出式（1-28）的结果.,有限元分析之大腕版 一定要选最变态的题目，什么材料非线性啊，几何非线性啊，接触非线性啊，多物理场耦合啊，都给他弄进去。是个模型就几百万个单元，上千万个节点，画个剖面图就要十几个小时。再整一并行机群，TOP500的，张口就是 High Performance Computation，一口地道的伦敦腔，倍儿有面子。题目一扔进去就跑个把月，你要是一个星期以内出结果，你都不好意思和别人打招呼。你说这样一趟算下来要发多少Paper? 10篇？10篇？就1篇！你还别嫌少，说不定人家还发在会议上。你得琢磨牛人的心理啊，有能耐算这样题目的人，根本就不在乎多发一篇两篇文章。什么叫大牛知道么？大牛就是不求灌水，但求经典。,