1、第四章第四章 平面三角形单元平面三角形单元 41 41 有限元法的基本思想有限元法的基本思想 42 42 三角形常应变单元三角形常应变单元 43 43 形函数的性质形函数的性质 44 44 刚度矩阵刚度矩阵 45 45 等效节点力载荷列阵等效节点力载荷列阵 46 46 有限元分析的实施步骤有限元分析的实施步骤 47 47 计算实例计算实例第四章第四章 平面三角形单元平面三角形单元一、有限元法的基本思想一、有限元法的基本思想 假想的把一连续体分割成数目有限的小体（单元），彼此间只在数目有限的指定点（结点）出相互连结只在数目有限的指定点（结点）出相互连结，组成一个单元的集合体以代替原来的连续体，再
2、在结点上在结点上引进等效力以代替实际作用于单元上的外力引进等效力以代替实际作用于单元上的外力。选择一个简单的函数来近似地表示位移分量的分布规律，建立位移和节点力之间的关系。有限元法的实质是：把有无限个自由度的连续体，把有无限个自由度的连续体，理想化为只有有限个自由度的单元集合体理想化为只有有限个自由度的单元集合体，使问题简化为适合于数值解法的结构型问题。二、经典解与有限元解的区别：二、经典解与有限元解的区别：微分 数目增到 建立一个描述连续体经 典 解 法（解析法）大小趋于 0 性质的偏微分方程 有限单元 离散化 集合 总体分析解有限元法连续体单元代替原连续体 （近似法）（单元分析）线性方程组
3、xy为平面应力问题，由于结构的对称性可取结构的1/4来研究，故所取的力学模型三、有限元法算题的基本步骤三、有限元法算题的基本步骤1.力学模型的选取力学模型的选取（平面问题，平面应变问题，平面应力问题，轴对称问题，空间问题，板，梁，杆或组合体等，对称或反对称等）例如：根据题目的要求，可选择适当的单元把结构离散化。对于平面问题可用三角元，四边元等。2.单元的选取、结构的离散化单元的选取、结构的离散化例如：结构离散化后，要用单元内结点的位移通过插值来获得单元内各点的位移。在有限元法中，通常都是假定单元的位移模式是多项式，一般来说，单元位移多项式的项数应与单单元位移多项式的项数应与单元的自由度数相等元
4、的自由度数相等。它的阶数至少包含常数项常数项和一次项一次项。至于高次项要选取多少项，则应视单元的类型而定。eNf3.选择单元的位移模式选择单元的位移模式（4-1）f单元内任一点的位移列阵；e单元的结点位移列阵；N单元的形函数矩阵；（它的元素是任一点位置坐标的函数）eB eBD4.单元的力学特性分析单元的力学特性分析 把（4-1）式代入几何方程可推导出用单元结点位移表示的单元应变表达式：（4-2）式中：单元内任一点应变列阵；B单元的应变矩阵；（它的元素仍为位置坐标的 函数）再把（4-）式代入物理方程，可导出用单元结点位移列阵表示的单元应力表达式：（4-3）最后利用弹性体的虚功方程建立单元结点力阵
5、与结点位移列阵之间的关系，即形成单元的刚度方程式：eeekR vTedxdydzBDBk式中：单元内任一点的应力列阵；D单元的弹性矩阵，（它与材料的特性有关）式中：单元刚度矩阵（4-4）（4-5）考虑整体结构的约束情况，修改整体刚度方程之后，（4-6）式就变成以结点位移为未知数的代数方程组。解此方程组可求出结点位移。用直接刚度法将单刚组集成总纲，并将组集成总载荷列阵，形成总体结构的刚度方程：ek K eR R（4-6）解出整体结构的结点位移列阵后，再根据单元结点的编号找出对应于单元的位移列阵，将代入（4-3）式就可求出各单元的应力分量值。e e RK5.建立整体结构的刚度方程建立整体结构的刚度
6、方程6.求解修改后的整体结构刚度方程求解修改后的整体结构刚度方程7.由单元的结点位移列阵计算单元应力由单元的结点位移列阵计算单元应力 求解出整体结构的位移和应力后，可有选择地整理输出某些关键点的位移值和应力值，特别要输出结构的 变形图、应力图、应变图、结构仿真变形过程动画图及整体结构的弯矩、剪力图等等。8.计算结果输出计算结果输出一、离散化一、离散化 在运用有限单元法分析弹性力学平面问题时，第一步就是要对弹性体进行离散化，把一个连续的弹性体变换为一个离散的结构物。对于平面问题，三角形单元是最简单、也是三角形单元是最简单、也是最常用的单元最常用的单元，在平面应力问题中，单元为三角形板，而在平面应
7、变问题中，则是三棱柱。假设采用三角形单元，把弹性体划分为有限个互不重叠的三角形。这些三角形在其顶点（即节点）处互相连接，组这些三角形在其顶点（即节点）处互相连接，组成一个单元集合体，以替代原来的弹性体。成一个单元集合体，以替代原来的弹性体。同时，将所有作用在单元上的载荷（包括集中载荷、表面载荷和体积载荷），都按虚功等效的原则移置到节点上，成为等效节点载荷等效节点载荷。由此便得到了平面问题的有限元计算模型，如图4-1所示。图4-1 弹性体和有限元计算模型 图4-2 平面三角形单元 ui(Ui)um(Um)uj(Uj)vj(Vj)vi(Vi)um(Um)j i m x y o 二、位二、位 移移
8、TmmjjiiTTmTjTievuvuvu Tiiivu 首先，我们来分析一下三角形单元的力学特性，即建立以单元节点位移表示单元内各点位移的关系式。设单元e的节点编号为i、j、m，如图4-2所示。由弹性力学平面问题可知，每个节点在其单元平面内的位移可以有两个分量，所以整个三角形单元将有六个节点位移分量，即六个自由度。用列阵可表示为：其中的子矩阵（i，j，m 轮换）(a)式中 ui、vi 是节点i在x轴和y轴方向的位移。(4-7)从弹性力学平面问题的解析解法中可知，如果弹性如果弹性体内的位移分量函数已知，则应变分量和应力分量也就体内的位移分量函数已知，则应变分量和应力分量也就确定了确定了。但是，
9、如果只知道弹性体中某几个点的位移分量的值，那么就不能直接求得应变分量和应力分量。因此，在进行有限元分析时，必须先假定一个位移模式位移模式。由于在弹性体内，各点的位移变化情况非常复杂，很难在整个弹性体内选取一个恰当的位移函数来表示位移的复杂变化，但是如果将整个区域分割成许多小单元，那么在每个单元的局部范围内就可以采用比较简单的函数在每个单元的局部范围内就可以采用比较简单的函数来近似地表示单元的真实位移来近似地表示单元的真实位移，将各单元的位移式连接 在有限单元法中，虽然是用离散化模型来代替原来的连续体，但每一个单元体仍是一个弹性体，所以在其每一个单元体仍是一个弹性体，所以在其内部依然是符合弹性力
10、学基本假设的内部依然是符合弹性力学基本假设的，弹性力学的基本方程在每个单元内部同样适用。uxyvxy123456起来，便可近似地表示整个区域的真实位移函数。这种化繁为简、联合局部逼近整体的思想，正是有限单元法的绝妙之处。基于上述思想，我们可以选择一个单元位移模式，单元内各点的位移可按此位移模式由单元节点位移通过插值而获得。线性函数是一种最简单的单元位移模式线性函数是一种最简单的单元位移模式，故设(b)式中 1、2、6是待定常数。因三角形单元共有六个自由度，且位移函数u、v在三个节点处的数值应该等于这些点处的位移分量的数值。假设节点i、j、m的坐标分别为（xi ,yi）、（xj ,yj）、（xm
11、 ,ym），代入(b)式，得：uxyvxyuxyvxyuxyvxyiiijiijjjjjjmmmmmm123456123456123456 ,mmjjiimmjjiimmmjjjiiiuxuxuxyuyuyuyxuyxuyxu11121,11121,213212111 xyxyxyiijjmm(c)由(c)式左边的三个方程可以求得 (d)其中(4-8)从解析几何可知，式中的 就是三角形i、j、m的面积。为保证求得的面积为正值，节点i、j、m的编排次序必须是逆逆时针方向时针方向，如图4-2所示。图4-2 平面三角形单元 ui(Ui)um(Um)uj(Uj)vj(Vj)vi(Vi)um(Um)j
12、i m x y o 将(d)式代入(b)式的第一式，经整理后得到 uab xc y uab xc y uab xc y uiiiijjjjmmmm12(e)mjmjimjmjijmmjmmjjixxxxcyyyybyxyxyxyxa1111vab xc y vab xc y vab xc y viiiijjjjmmmm12Nab xc yiiii12其中同理可得若令这样，位移模式(e)和(f)就可以写为（i,j,m轮换）(4-10)（i,j,m轮换）(4-9)(f)式中 I是二阶单位矩阵；Ni 、Nj 、Nm 是坐标的函数，它们反映了单元的位移状态，所以一般称之为形状函数，简称形函数。矩阵矩阵
13、 N 叫做形函数矩阵叫做形函数矩阵。三节点三角形单元的形函数是坐标的线性函数。单元中任一条直线发生位移后仍为一条直线，即只要两单元在公共节点处保持位移相等。则公共边线变形后仍为密合。uN uN uN uvN vN vN viijjmmiijjmm eemjiNINININvuf(4-11)也可写成矩阵形式(4-12)三、应三、应 变变 xyxyuxvyuyvx 12000000bbbccccbcbcbijmijmiijjmme有了单元的位移模式，就可以利用平面问题的几何方程求得应变分量。将(e)、(f)两式代入上式，即得：(g)uab xc y uab xc y uab xc y uiiiij
14、jjjmmmm12(e)vab xc y vab xc y vab xc y viiiijjjjmmmm12(f)Be BBBBijmBbccbiiiii1200可简写成 其中 B 矩阵叫做单元应变矩阵，可写成分块形式而子矩阵由于和bi、bj、bm、ci、cj、cm 等都是常量，所以矩阵B中的诸元素都是常量，因而单元中各点的应变分量也都是常量，通常称这种单元为常应变单元。（i,j,m轮换）(4-15)(4-14)(4-13)四、应四、应 力力 D Be SD B Se D 求得应变之后，再将（4-13）式代入物理方程 ，便可推导出以节点位移表示的应力。即(4-16)(h)(4-17)令则 SD
15、 BBBSSSijmijm DE11100122对称 SD BEbcbccbiiiiiiii2 112122其中 S叫做应力矩阵，若写成分块形式，有对于平面应力问题，弹性矩阵D为(4-18)(i)所以，S的子矩阵可记为（i,j,m轮换）(4-19)DE111211100122 1对称 SD BEbcbccbiiiiiiii12 11211122 1122 1 对于平面应变问题，只要将(i)式中的E换成E/1-2，换成 /1-，即得到其弹性矩阵(j)（i,j,m轮换）(4-20)SSSiijjmm注意到（4-7）式，则有(4-21)由（4-19）、（4-20）式不难看出，S中的诸元素都是常量，所
16、以每个单元中的应力分量也是常量应力分量也是常量。可见，对于常应变单元，由于所选取的位移模式是线性的，因而其相邻单元将具有不同的应力和应变相邻单元将具有不同的应力和应变，即在单在单元的公共边界上应力和应变的值将会有突变元的公共边界上应力和应变的值将会有突变，但位移却是连续的。Nab xc yiiii122111 xyxyxyiijjmm在上节中，提出了形函数的概念，即其中（i,j,m轮换）现在我们来讨论一下形函数所具有的一些性质。根据行列式的性质：行列式的任一行（或列）的元素与其相应的代数余子式的乘积之和等于行列式的值，而任一行（或列）的元素与其他行（或列）对应元素的代数余子式乘积之和为零，并注
17、意到（4-9）式中的常数ai、bi、ci，aj、bj、Nxyab xc yiiiiiiim,121Nxyab xc yijjiijij,120Nxyab xc yimmiimim,120cj 和am、bm、cm 分别是行列式的第一行、第二行和第三行各元素的代数余子式，我们有 形函数在各单元节点上的值，具有形函数在各单元节点上的值，具有“本点是本点是1、它、它点为零点为零”的性质的性质，即在节点i上，在节点j、m上，(a)(b)(c)NxyNxyNxyNxyNxyNxyjiijjjjmmmiimjjmmm,101000 12121,ycccxbbbaaaycxbaycxbaycxbayxNyxN
18、yxNmjimjimjimmmjjjiiimji类似地有(d)在单元的任一节点上，三个形函数之和等于在单元的任一节点上，三个形函数之和等于1，即，即(e)NNNijm 10,1,yxNxxxxyxNxxxxyxNmijijijii简记为(4-22)这说明，三个形函数中只有二个是独立的。三角形单元任意一条边上的形函数，仅与该边的两端节三角形单元任意一条边上的形函数，仅与该边的两端节点坐标有关、而与其它节点坐标无关。点坐标有关、而与其它节点坐标无关。例如，在i j 边上，有(4-23)uN uN uvN vN viijjiijj 例如，对图4-3所示的单元ijm和ijn，具有公共边ij。这样，不论
19、按哪个单元来计算，根据（4-11）式，公共边ij上的位移均由下式表示jinmxyo图4-3由(4-23)式可知，在ij边上式中 Ni ,Nj 的表达形式如（4-23）式所示。(i)利用形函数的这一性质可以证明，相邻单元的位移分别进行线性插值之后，在其公共边上将是连续的。mmjjiiLLL由此可见，在公共边上的位移u、v 将完全由公共边上的两个节点i、j 的位移所确定，因而相邻单元的位移是保持连续的。为了在以后讨论问题中能够比较方便地确定单元中任意一点处的形函数数值，这里引入面积坐标的概念。在图4-4所示的三角形单元ijm中，任意一点P(x,y)的位置可 以用 以下三个比值来确定oyxLi =0
20、Li =1/4Li =1/2Li =3/4Li =1Pjim图4-4 式中 为三角形单元ijm的面积，i、j、m 分别是三角形Pjm、Pmi、Pij的面积。这三个比值就叫做P点的面积坐标。(4-24)mji1mjiLLLycxbayxyxyxiiimmiii2111121ycxbaLiiiii21显然这三个面积坐标并不是完全独立的，由于所以有：而三角形pjm的面积为：故有：Lab xc yjjjjj12Lab xc ymmmmm12类似地有(4-25)(4-26)由此可见，前述的三角形常应变单元中的形函数Ni、Nj、Nm 就是面积坐标Li、Lj、Lm 。根据面积坐标的定义，我们不难发现，在平行
21、jm边的直线上的所有各点，都有相同的坐标Li，并且该坐标就等于“该直线至jm边的距离”与“节点i至jm边的距离”之比，图4-4中给出了Li 的一些等值线。xx Lx Lx Lyy Ly Ly LLLLiijjmmiijjmmijm 1容易看出，单元三个节点的面积坐标分别为节点 i：Li =1 Lj =0 Lm =0节点 j：Li =0 Lj =1 Lm =0 节点m：Li =0 Lj =0 Lm =1不难验证，面积坐标与直角坐标之间存在以下变换关系：(4-27)xLxLLxLLxLbLbLbLiijjmmiijjmm222yLyLLyLLyLcLcLcLiijjmmiijjmm222当面积坐标
22、的函数对直角坐标求导时，可利用下列公式：(4-28)一一.单元刚度矩阵单元刚度矩阵 RRRRUVUVUVeiTjTmTTiijjmmT eiijjmmTuvuvuv 为了推导单元的节点力和节点位移之间的关系，可应用虚位移原理对图4-2中的单元e进行分析。单元e是在等效节点力的作用下处于平衡的，而这种节点力可采用列阵表示为(a)假设在单元e中发生有虚位移，则相应的三个节点i、j、m 的虚位移为且假设单元内各点的虚位移为f*，并具有与真实位移相同的位移模式。fNe Be()eTeR Ttdxdy故有(c)参照（4-13）式，单元内的虚应变*为于是，作用在单元体上的外力在虚位移上所做的功可写为(d)
23、(f)而单元内的应力在虚应变上所做的功为(g)tdxdyBDBeTTe)()()eTeeTTeRBD Btdxdy RBD B tdxdyeTe这里我们假定单元的厚度t为常量。把（d）式及（4-16）式代入上式，并将虚位移提到积分号的前面，则有根据虚位移原理，由（f）和（h）式可得到单元的虚功方程，即注意到虚位移是任意的，所以等式两边与相乘的项应该相等，即得 tdxdyBDBkTe eeekR记(4-32)则有(4-33)上式就是表征单元的节点力和节点位移之间关系的刚节点力和节点位移之间关系的刚度方程度方程，ke就是单元刚度矩阵。如果单元的材料是均质的，那么矩阵D 中的元素就是常量，并且对于三
24、角形常应变单元，B矩阵中的元素也是常量。当单元的厚度也是常量时，所以（4-32）式可以简化为ke=BT DBt (4-34)二二 整体刚度矩阵整体刚度矩阵 2112nTTnTT iiiTuv 讨论了单元的力学特性之后，就可转入结构的整体分析。假设弹性体被划分为N个单元和n个节点，对每个单元按前述方法进行分析计算，便可得到N组形如（4-33）式的方程。将这些方程集合起来，就可得到表征整个弹性体的平衡关系式。为此，我们先引入整个弹性体的节点位移列阵 2n1，它是由各节点位移按节点号码以从小到大的顺序排列组成，即其中子矩阵(j)(i=1,2,n)(k)是节点节点i的位移分量的位移分量。RRRRnTT
25、nTT2112 RXYUViiiTieeNieeNT11 继而再引入整个弹性体的载荷列阵R2n1，它是移置到节点上的等效节点载荷依节点号码从小到大的顺序排列组成，即(l)其中子矩阵(i=1,2,n)(m)是节点i上的等效节点载荷等效节点载荷。nmjinmjikkkkkkkkkkmmmjmijmjjjiimijiinn1 122(q)同样，将六阶方阵k加以扩充，使之成为2n阶的方阵 123421q图 4-5组装总刚组装总刚k的一般规则：的一般规则：1.当krs中r=s时，该点被哪几个单元所共有，则总刚子矩阵krs就是这几个单元的刚度矩阵子矩阵krse的相加。2.当krs中r s时，若rs边是组合
26、体的内边，则总体刚度矩阵krs就是共用该边的两相邻单元单刚子矩阵krse的相加。3.当krs中r和s不同属于任何单元时，则总体刚度矩阵krs=0。下面，我们考查一个组装总刚的实例：下面，我们考查一个组装总刚的实例：1.整体刚度矩阵及载荷列阵的整体刚度矩阵及载荷列阵的 组集组集 根据叠加原理，整体结构的各个刚度矩阵的元素显然是由有关单元的单元刚度矩阵的元素组集而成的，为了便于理解，现结合图4-5说明组集过程。子块33333231232221131211KKKKKKKKKKe 图中有两种编码：一是节点总码：1、2、3、4；二是节点局部码，是每个单元的三个节点按逆时针方向的顺序各自编码为1，2，3。
27、图中两个单元的局部码与总码的对应关系为：单元 1 ：1，2，3 1，2，3 单元 2 ：1，2，3 3，4，1或：单元 1 ：1，2，3 1，2，3 单元 2 ：1，2，3 1，3，4单元e的刚度矩阵分块形式为：子块4444434241343332312423222114131211KKKKKKKKKKKKKKKKK整体刚度矩阵分块形式为：其中每个子块是按照节点总码排列的。通常，采用刚度集成法或直接刚度法来组集整体结构刚度矩阵。刚度集成法分两步进行。eK eK 第一步，把单元刚度矩阵 扩大成单元的贡献矩阵 ，使单元刚度矩阵的四个子块按总体编号排列，空白处作零子块填充。2K 第二步，以单元 2
28、为例，局部码1，2，3 对应于总码3，4，1，按照这个对应关系扩充后，可得出单元 2 的贡献矩阵 。213444342413433323124232221141312112KKKKKKKKKKKKKKKKK 总码 1 2 3 4 1 2 3 4 3 1 2 局部码 1K用同样的方法可得单元 1 的贡献矩阵 。K 第二步，把各单元的贡献矩阵对应行和列的子块相叠加，即可得出整体结构的刚度矩阵 ，如（4-42）式。K22nn22 在这里应该指出，整体刚度矩阵 中每个子块为 阶矩阵，所以若整体结构分为n个节点，则整体刚度矩阵的阶数是 。TeynxnyxyxeFFFFFFFF 212211 总码 1 2
29、 3 4 1 2 3 （4-42）3 1 2 局部码 21332100432123323223122322213313222113112312212121321211311221111121KKKKKKKKKKKKKKKKKKKKee F eF 至于整体结构的节点载荷列阵 的组集，只需将各单元的等效节点力列阵 扩大成2n行的列阵，然后按各单元的节点位移分量的编号，对应相叠加即可三三 整体刚度矩阵的性质整体刚度矩阵的性质 nynxyxyxnnnnnnnnRRRRRRvuvuvukkkkkkkkk221122112,22,21,22,222212,11211 由总刚度方程可知：欲 使 弹 性体的某
30、一节点在坐标轴方向发生单位位移，而其它节点都保持为零的变形状态，在各节点上所需要施加的节点力。刚度矩阵刚度矩阵K中每一列元素的物理意义为中每一列元素的物理意义为：RRRRRRkkkkkkxyxynxnyTnnT1122112131412121 1 ijk 令节点1在坐标轴x方向的位移u1=1，而其余的节点位移v1=u2=v2=u3=v3=u2n =v2n =0，这样就可得到节点载荷列阵等于K的第一列元素组成的列阵，即即表示：是在j节点有单位位移时，而在I节点所需施加的力。(s)KkBD BtBD BtkKsrTsrTeNsTrTeNrTseNrseNrs1111 刚度矩阵刚度矩阵K中主对角元素
31、总是正的。中主对角元素总是正的。例如，刚度矩阵K中的元素k33 是表示节点2在x方向产生单位位移，而其它位移均为零时，在节点2的x方向上必须施加的力，很显然，力的方向应该与位移方向一致，故应为正号。刚度矩阵刚度矩阵K是一个对称矩阵，即是一个对称矩阵，即Krs=Ksr T。由（4-32）、（4-36）式得所以，可以只存储上三角或下三角矩阵。(t)刚度矩阵刚度矩阵K是一个稀疏矩阵。是一个稀疏矩阵。如果遵守一定的节点编号规则，就可使矩阵的非零元素都集中在主对角线附近呈带状。前面在讨论总刚子矩阵的计算时曾指出，总刚中第r双行的子矩阵Krs，有很多位置上的元素都等于零，只有当第二个下标s等于r或者s与r
32、同属于一个单元的节点号码时才不为零，这就说明，在第r双行中非零子矩阵的块数，应该等于节点r周围直接相邻的节点数目加一。可见，K的元素一般都不是填满的，而是呈稀疏状（带状）。以图4-6a所示的单元网格为例，其整体刚度矩阵中的非零子块（每个子块为2行2列）的分布情况如图4-6b所示。123456789101112131415图 4-6 a 1 2 3 4 5 6 7 8 9 10 11 12 13 1415123456789101112131415图 4-6 b 半带宽B=（相邻节点号的最大差值D+1）*2 若第r双行的第一个非零元素子矩阵是Krl，则从Krl 到Krr 共有（r-l+1）个子矩阵
33、，于是K的第2r行从第一个非零元素到对角元共有2(r-l+1)个元素。显然，带状刚度矩带状刚度矩阵的带宽取决于单元网格中相邻节点号码的最大差值阵的带宽取决于单元网格中相邻节点号码的最大差值D。我们把半个斜带形区域中各行所具有的非零元素的最大个数叫做刚度矩阵的半带宽（包括主对角元），用B表示，即 B=2(D+1)。通常的有限元程序，一般都利用刚度矩阵的对称和稀疏带状的特点，在计算求解中，只存储上半带的元素，即所谓的半带存储半带存储。因此，在划分完有限元网格进行节点编号时，要采用合理的编码方式，使同一单元中相邻两节点的号码差尽可能小，以便节省存储空间、提高计算效率。刚度矩阵刚度矩阵K是一个奇异矩阵
34、，在排除刚体位移后，是一个奇异矩阵，在排除刚体位移后，它是正定阵。它是正定阵。弹性体在R的作用下处于平衡，R的分量应该满足三个静力平衡方程。这反映在整体刚度矩阵K中就意味着存在三个线性相关的行或列，所以K是个奇异阵，不存在逆矩阵。()eTeTTTRfGfq tdsfp tdxdy 在上节讨论整体刚度矩阵时已经指出载荷列阵R 是由弹性体的全部单元的等效节点力集合而成，而其中单元的等效节点力单元的等效节点力Re 则是由作用在单元上的集中力、则是由作用在单元上的集中力、表面力和体积力分别移置到节点上，再逐点加以合成求表面力和体积力分别移置到节点上，再逐点加以合成求得。得。根据虚位移原理，等效节点力的
35、大小，应按其所做等效节点力的大小，应按其所做的功与作用在单元上的三种力在任何虚位移上所做的功的功与作用在单元上的三种力在任何虚位移上所做的功相等这一原则来确定相等这一原则来确定。即 上式中等号的左边表示单元的等效节点力单元的等效节点力Re 所做所做的虚功的虚功；等号右边的第一项是集中力第一项是集中力G所做的虚功所做的虚功、第二项的积分是沿着单元的边界进行，表示面力q所做的虚功、第三项的积分则是遍及整个单元，表示体积力p所做的虚功；t为单元的厚度，假定为常量。(a)一、体力的移置：一、体力的移置：303030WWWFey0 xijmcwuP 01010130303032cTuuuextPPPPc
36、x容量tW eF如果任意三角形单元ijk的重心c上受有自重 则按刚体静力等效原理可把W直接移置到i，j，m三个节点上而组成 ：2cuXtP如果单元的重心c受有惯性力Pu作用，且 ，则Pu移置到i，j，m节点上的等效结点力为：式中：旋转角速度 是单元重心处x坐标。二、面力的移置：二、面力的移置：y0 xijmlq2ql2ql已知在ij边受有面力q，则移置到i、j结点上的等效节点力为：TeqltQ0010102y0 xijml0q20lq3l eQ TeltqQ0320310020当某一边上有三角形分布的面力时，可由刚体静力等效直接写出三、集中力的移置：三、集中力的移置：如集中力G做用于其一边界上
37、如图：先将G分解为 ，后，分别按线段的比例把 和 分别移置到i，j两点上。即：xPyPxPyP 001122llPllPllPllPFyxyxev 即按静力平衡方法分配。即按静力平衡方法分配。0 xijmG2lxPyP1lly根据前面的讨论，现以三角形常应变单元为例来说明应用有限元法求解弹性力学平面问题的具体步骤。确定力学模型，确定力学模型，根据工程实际情况确定问题的力学模型，并按一定比例绘制结构图、注明尺寸、载荷和约束情况等。将计算对象进行离散化将计算对象进行离散化，即弹性体划分为许多三角形单元，并对节点进行编号。确定全部节点的坐标值，对单元进行编号，并列出各单元三个节点的节点号。计算载荷的
38、等效节点力计算载荷的等效节点力（要求的输入信息）。由各单元的常数bi、ci、bj、cj、bm、cm 及行列式2，计算单元刚度矩阵计算单元刚度矩阵。组集整体刚度矩阵组集整体刚度矩阵，即形成总刚的非零子矩阵。处理约束处理约束，消除刚体位移。求解线性方程组求解线性方程组，得到节点位移。计算应力矩阵，求得单元应力计算应力矩阵，求得单元应力，并根据需要计算主应力和主方向。整理计算结果整理计算结果（后处理部分）。为了提高有限元分析计算的效率、达到一定的精度，应该注意以下几个方面。一一.对称性的利用对称性的利用在划分单元之前，有必要先研究一下计算对象的对称或反对称的情况，以便确定是取整个物体，还是部分物体作
39、为计算模型。例如，图4-11(a)所示受纯弯曲的梁，其结构对于x、y轴都是几何对称的，而所受的载荷则是对于y轴对称，对于x轴反对称。可知，梁的应力和变形也将具有同样的对称特性，所以只需取1/4梁进行计算即可。取分离体如图4-11(b)所示，对于其它部分结构对此分离体的影响，可以作相应的处理，即对处于y轴对称面内各节点的x方向位移都设置为零，而对于在x轴反对称面上的各节点的x方向位移也都设置为零。这些条件就等价于在图4-11(b)中相应节点位置处施加约束，图中o点y方向施加的约束是为了消除刚体位移。RRRRyxooyxR(a)(b)图4-11节点的布置是与单元的划分互相联系的。通常，集中载荷的作
40、用点、分布载荷强度的突变点，分布载荷与自由边界的分界点、支承点等都应该取为节点。并且，当物体是由不同的材料组成时，厚度不同或材料不同的部分，也应该划分为不同的单元。节点的多少及其分布的疏密程度（即单元的大小），一节点的多少及其分布的疏密程度（即单元的大小），一般要根据所要求的计算精度等方面来综合考虑。般要根据所要求的计算精度等方面来综合考虑。从计算结果的精度上讲，当然是单元越小越好，但计算所需要的时间也要大大增加。另外，在微机上进行有限元分析时，还要考虑计算机的容量。因此，在保证计算精度的前提下，应力求采用较少的单元。为了减少单元，在划分单元时，对于应力变化梯度较大的部位单元可小一些，而在应力
41、变化比较平缓的区域可以划分得粗一些。二二.节点的选择及单元的划分节点的选择及单元的划分(a)(b)图4-12还有一点值得注意的是，单元各边的长度不要相差太大，以免出现过大的计算误差或出现病态矩阵。例如，图4-12所示的(a)、(b)两种单元划分，虽然都是同样的四个节点，但(a)的划分方式显然要比(b)的方式好。在进行节点编号时，应该注意要尽量使同一单元的相邻节点的号码差尽可能地小，以便最大限度地缩小刚度矩阵的带宽，节省存储、提高计算效率。如前所述，平面问题的半带宽为B=2(d+1)三三.节点的编号节点的编号若采取带宽压缩存储，则整体刚度矩阵的存储量N 最多为 N=2nB=4n(d+1)其中 d
42、为相邻节点的最大差值，n为节点总数。1 2 3 4 5 6 71 3 5 7 9 11 138 9 10 11 12 13 142 4 6 8 10 12 14例如在图3-13中，(a)与(b)的单元划分相同，且节点总数都等于14，但两者的节点编号方式却完全不同。(a)是按长边进行编号，d=7，N=488；而(b)是按短边进行编号，d=2，N=168。显然(b)的编号方式可比(a)的编号方式节省280个存储单元。(a)(b)图4-13四四.单元节点单元节点i、j、m的次序的次序在前面章节中，我们曾指出，为了在计算中保证单元的面积 不会出现负值，节点i、j、m的编号次序必须是逆必须是逆时针方向时
43、针方向。事实上，节点i、j、m的编号次序是可以任意安排的，只要在计算刚度矩阵的各元素时，对取绝对值，即可得到正确的计算结果。在实际计算时，应该注意所选有限元分析软件的使用要求。五五.边界条件的处理及整体刚度矩阵的修正边界条件的处理及整体刚度矩阵的修正在前面讨论整体刚度矩阵时，已经提到，整体刚度矩阵的奇异性可以采用提高考虑边界约束条件的方法来排除弹性体的刚体位移，以达到求解的目的。一般情况下，求解的问题，其边界往往已有一点的位移位移约束条件约束条件，本身已排除了刚体运动的可能性。否则的话，就必须适当指定某些节点的位移值适当指定某些节点的位移值，以避免出现刚体位移。这里介绍两种比较简单的引入已知节
44、点位移的方法，这两种方法都可保持原K矩阵的稀疏、带状和对称等特性。下面我们来实际考察一个只有四个方程的简单例子。1.划划0置置1法。法。保持方程组为2n2n系统，仅对K和R进行修正。例如，若指定节点i在方向y的位移为vi ，则令K中的元素k2i,2i 为1，而第2i行和第2i列的其余元素都为零。R中的第2i个元素则用位移vi 的已知值代入，R中的其它各行元素均减去已知节点位移的指定值和原来K中该行的相应列元素的乘积。KKKKKKKKKKKKKKKKuvuvRRRR111213142122232431323334414243441122123410000000100022244244112212
45、21123334411433KKKKuvuvRKKRKK假定该系统中节点位移u1 和u2分别被指定为当引入这些节点的已知位移之后，方程(a)就变成然后，就用这组维数不变的方程来求解所有的节点位移。显然，其解答仍为原方程(a)的解答。u1=1 ，u2=2KKKKKKKKKKKKKKKKuvuvKRKR111512131421222324313233153441424344112211115233315410101010KuK vK uK vK11151121132142111151010 乘大数法。乘大数法。将K中与指定的节点位移有关的主对角元素乘上一个大数，如1015，同时将R中的对应元素换成
46、指定的节点位移值与该大数的乘积。实际上，这种方法就是使K中相应行的修正项远大于非修正项。若把此方法用于上面的例子，则方程(a)就变成事实上，该方程组的第一个方程为 图4-16所示为一平面应力问题离散化以后的结构图，其中图（a）为离散化后的总体结构，图（b）为单元1，2，3，4的结构，图（c）为单元3的结构。用有限元法计算节点位移、单元应变及单元应力（为简便起见，取泊松比 ，单元厚度t=1）。0 xyp1234651234a3ijmivaiXiujXjuiYjvjYmvmYmXmuaa1,2,4ijmiXiviujXjuiYjvjYmvmYmXmu 图 4-16 计算实例2的结构图首先求确定各单
47、元刚度所需的系数 及面积A，对于单元1，2，4有：mjimjicccbbb,2/0,02aAcacacababbmjimji解：解：对于单元3有:2/,0,0,2aAacaccabbabmjimji 其次，求出各单元的单元刚度矩阵。对于1，2，4单元，其单元刚度矩阵为：10110102020010312112130100202010110144,2,1Ekijmijm 各单元的节点编号与总体结构的总编号之间的对应关系见表4-2。31211013011220200011011011011002000243Ek 对于单元3，其单元刚度矩阵为：ijmijm 各单元节点号与总体节点号对应表各单元节点号
48、与总体节点号对应表 单元号 1 2 3 4 节点号 节 点 总 编 号 I 1 2 2 3 j 2 4 5 5 m 3 5 3 6表表4-2 将各单元刚度矩阵按节点总数及相应的节点号关系扩充成12*12矩阵,分别如下:654321000000000000000000000000000000000000000000000000000000000000000000000000000000101101000000020200000000103121000000121301000000002020000000101101465432111212Ek654321000000000000000000000
49、000001011000100000202000000001031002100001213000100000000000000000000000000000020002000001011000100000000000000000000000000465432121212Ek654321000000000000000000000000002000200100000100111000000000000000000000000000002100311000000100131200000100111000000000020200000000000000000000000000465432131212E
50、k654321101100010000020200000000103100210000121300010000000000000000000000000000002000200000101100010000000000000000000000000000000000000000000000000000465432141212Ek将扩充后的各单元刚度矩阵相加,得总体刚度矩阵K，即：6543211011000100000202000000001061114101001216021210000010310021000012130001000041006121011012001614000001202
