有限元 2-弹性力学平面问题有限单元法(2.3程序设计,2.4矩形单元,2.5六节点三角形单元)

1、有限元讲义 1 2.3 平面问题有限元程序设计 2.3 平面问题有限元程序设计 一、程序设计方法与结构分析程序的特点 一、程序设计方法与结构分析程序的特点 1程序设计方法论简述程序设计方法论简述 借助计算机来完成某项工作，通常都要先编写相应的计算机程序，或叫程序设计。完 成一个结构分析或结构CAD系统也必然要经过程序设计才能实现。 程序设计要使用专门的程序语言。我国结构程序设计中所采用的语言，在60年代和70 年代初以ALGOL语言ALGOL语言为主。此后逐步广泛使用的主要是BASICBASIC语言和FORTRANFORTRAN语言，随着CAD 和人工智能技术的发展，PASCAL、 C、LIS

2、P、 PROLOG等有着各自特长的程序语言也逐步 进入土木工程领域的计算机程序设计中。 过去人们通常认为，程序设计的中心问题就是学会使用一种程序语言，用以编写程序。然 而学会用程序语言编程只是整个程序设计中的一部分。据有关资料介绍，编写程序在整个 系统的研制过程中仅占1515的工作量。在一个大型程序系统的整个存在阶段的工作量中， 在系统投入使用后的维护工作量为原来研制工作量总和的两倍（这一点在作者所从事的软 件开发工作中也得到充分的证明）。 维护工作量是如此之高，这就使我们必须注意到，在程序研制阶段便即应当考虑为以后的 维护工作提供方便，哪怕是为此要增加一些额外的工作量也是值得的。 要编制一个

3、好的程序系统并没有一种绝对的规则，就象是工程设计没有一种绝对规则 一样。但对于程序设计的好坏现在已逐渐形成了一套评价的客观标准。这些标准大致分为 以下几个主要方面： (1) 程序的可读性; (2) 正确性与可靠性; (3) 使用方便且效率高; (4) 软件的可移置性; (5) 易于调试与维护。 直到1970年代中期人们才认识到软件的维护是软件研究的一个关键领域。造成软件维 护工作量大的原因之一是与程序研制过程中所采用的设计方法不够科学化有关。为了解决 这一问题，人们开展了对于程序设计方法论的研究与实践，其目标是使软件正确、可靠和 降低整个软件研制活动的费用。总的来说，程序设计已从强调灵活的技巧

4、和局部效率向着 强调程序结构化和整体功能的方向发展。这实际上就是逐步发展起来的关于程序设计编写 与调试的一套方法论，其要点可归纳为以下几方面： (1) 编程结构化。为了使程序设计者能按照一定的结构形式，而不是随心所欲地设计 和编写程序，使编制的程序易读、易修改，以提高程序设计和维护工作的效率，荷兰学者 Dijkstra提出了“结构化程序设计方法”。结构化程序规定了三种基本的结构形式，他们 是顺序结构、分支选择结构和循环结构。编程结构化又称结构化程序设计，他可使编写的 程序层次分明，逻辑清楚，容易阅读。 (2) 分层处理技术。为了解决现实世界中的许多复杂问题，人们往往需要根据问题的 内在联系将其

5、分割成有层次的一系列问题来分别求解。对于一个大型程序系统设计来说， 也需采用分层的办法来处理， 在每一层里集中解决一个问题， 并为下一层的执行作好准备。 分层处理技术的主要内容是将程序划分为多个层次的若干模块，每个模块完成一个（或几 个）预定的功能。 为了保证模块的独立性，各模块之间只能通过接口与其他模块连接。另外，对于一个 较大的软件系统要由多人合作才能完成，模块化也为此提供较好的合作条件。模块化的一 种典型结构形式就是子程序结构。 (3) 避免过多使用GOTO语句，特别是逆转的GOTO语句。 这是结构化程序设计的基本要求之一。 有限元讲义 2 （4）可移植性 对于应用软件， 特别是大型的C

6、AD软件， 可移植性高低同样是衡量软件的质量的重要指标。 因为研制一个大型应用软件系统不但要耗费大量的人力物力，而且还要花费相当长的时间 才能研制成功。在其研制和应用期间，不可避免地发生运行环境的变化计算机硬件设 备的换代和系统软件的更新。可移植性强即等于延长了软件的生存期，从而可节省新的开 发投资。可移植性主要体现在软件对支撑环境的独立性和软件本身的封闭性。例如，用 FORTRAN、BASIC 、C等语言编写，并尽量避免采用非标准语句和函数。此外，在软件中采 用统一的I/O 模块，也是提高可移植性的手段之一。 应当指出，程序设计方法论仍在发展探索之中，千万不能把上述有关内容当成一成不 变的教

7、条套用，而应当通过实践来发展和丰富其内容。然而程序设计发展到今天，已经奠 定了很多必要的理论基础。我们正在达到一个可以谈论程序设计是一门科学而不仅仅是一 种技巧的阶段。 2、结构化程序设计 2、结构化程序设计 结构化程序设计，又称结构程序设计(Structured Programming)是 荷兰学者E. W. dijkstra首先提出来的，简称SP。 人们对SP有各种定义和解释: 有人说它是: 指导程序员编程的一般方法； 有人说它是: 不使用goto语句的程序设计; 有人说它是: 自顶向下的程序设计。 有人把层次结构、顺序结构、选择结构和重复结构定义为SP的全部内容，而把一个程 序结构模块定

8、义为上述四种基础结构的某种集合，即： 结构模块(程序)(层次结构, 顺序结构, 选择结构, 重复结构) 认为SP就是定义为这些基本结构在软件开发和维护中的严格应用。 3、程序的正确性及其验证 3、程序的正确性及其验证 程序正确性(Program Correctness)被定义为一个程序和它打算要实现的功能之间的 一次符合(a correspondness) Gries说，在程序设计的初期阶段，人们很少看到程序正确性方面的问题，那时人们 往往着重于调试(debugging),但是后来发现调试好的程序并不能表示错误不存在。于是人 们就不得不考虑程序的正确性证明问题。 什么是程序正确性证明呢？从词义

9、上讲，“证明”就是提供一种有力的证据使人们在 思想上不得不接受一条真理或一件事实。但是在程序设计中要做到这一点并不容易，因为 现在还没有找到一种象数学一样非常完整的形式。目前主要还是通过选择较好的算法，大 量考题，与已有结果比较等办法来证明一个程序的正确性。但不管怎样，程序正确性理论 已经被提出来，并有了一些实际想法。 程序的验证(Program Verification)实际上就是检验程序的正确性。 一个程序如果有错误， 主要是两方面的：一方面是语法错误，这部分比较好解决，一般是在调试阶段（编辑阶段） 完成。但是一个语法完全没问题的程序并不一定是正确的。因为程序中的许多部分往往是 靠逻辑关系

10、来达到所要实现的目的。目前应用程序的验证主要靠针对程序每一功能，每一 逻辑分支进行各种类型的考题，包括考题的规模。 4、结构分析程序的特点 4、结构分析程序的特点 1. 规律性、 通用性好。 有限元讲义 2. 计算工作量大、运算时间长（形成单、总刚、解方程、动力、非线性）。 3. 矩阵阶数高，与微机的存贮量发生矛盾。 二、平面问题有限元程序设计（三角形单元） 二、平面问题有限元程序设计（三角形单元） 三角形单元是平面问题中最简单的单元形式，但其他单元形式的有限元程序与三角形 单元程序的结构基本相同，并无本质区别。故本节通过一个FORTRAN 源程序，详细介绍三 角形单元的程序设计。 1.1.

11、程序总框图 程序总框图 输入结构控制参数 输入其它数据 形成整体刚度阵 引入支承条件 解方程，输出位移 求应力，输出应力 形成节点荷载向量 开始 结束 1 单元面积 求弹性矩阵 单元刚度矩阵 位移-应变矩阵 2 3 4 5 6 7 8 9 10 上述总框图左边主线对应着源程序中的主程序流程（见源程序）。 现按照总框图中的10个子框图，结合详细的子框图设计，介绍附录中给出的“平面问题三 角形单元源程序”。 2子框图1（SDATA2子框图1（SDATA） 输入其他数据并计算半带宽 输入5个控制参数后，程序运行中所需的其他原始数据均放在该子程序中输入。其框 图如下： 3 有限元讲义 4 有限元讲义

12、半带宽的计算： 如图所示结构，9个单元10个结点。将单刚按结点分块，采用直接刚度法以子块对号入座 的方式可形成图示总刚。 图中整体刚度矩阵K的非零元素分布在以主对角线为中心的斜带形区域内(图中用 粗线标明),这种矩阵称作带形矩阵。在半个斜带形区域中(包括主对角线元素在内),每行 具有的元素个数称为半带宽,用d表示。由图中看出,在半带中,每行有五个子块,即十个元 素,因此半带宽d=10。半带宽d的一般公式是： 半带宽d=(相邻结点码的最大差值+1)结点位移未知量数 半带宽d=(相邻结点码的最大差值+1)结点位移未知量数 三角形单元： 三角形单元： 半带宽d=(相邻结点码的最大差值+1)2 半带宽

13、d=(相邻结点码的最大差值+1)2 图中相邻结点码的最大差值是4,故d=(4+1)2=10 根据带形矩阵的特点,并利用矩阵的对称性,则在计算机中可只存贮上半带的元素。这种存 贮方式称为半带存贮。 如下面左图总刚，半带存贮时,只从K中取出上半部斜带中的元素,存贮在右图的矩 阵K *的竖带中,矩阵K *为n行d列,只有nd个元素。因此,实际存贮量与K中元素总数 之比为d/n。由此可见,d值愈小,则存贮量愈省。 由K改成K *时，元素的行码不变，新的列码改为： 新列码=原列码-行码+1 （存上三角） 新列码= （存下三角）？ 5 有限元讲义 3子框图2（STE）3子框图2（STE） 计算单元刚度矩阵

14、KE 6 有限元讲义 4子框图3（AE） 计算单元面积AE 4子框图3（AE） 计算单元面积AE 5子框图4（DTE） 计算弹性参数矩阵D 5子框图4（DTE） 计算弹性参数矩阵D 7 有限元讲义 6子框图5（BTE） 计算几何矩阵B 6子框图5（BTE） 计算几何矩阵B 8 有限元讲义 7子框图6（STIFF） 形成总刚度矩阵KS 7子框图6（STIFF） 形成总刚度矩阵KS (1) 结点局部码与总码的对应关系（结点换码） 图示为单元3个结点的两种编码，某个结点的局部码如果是i(1),则它在总刚中的对应总码 是LND(IE,1)。 9 有限元讲义 （2）同一子块在单刚KE和总刚KS的对应位置

15、（子块搬家） （3）同一元素在KE、总刚KS和带状总刚KS *中的对应位置（元素搬家） 现在讨论形成总刚中某个元素的搬家问题 设该元素属于单刚矩阵KE中的I行J列子块,该元素在这个子块中的位置是II行JJ列。则 此元素在单元刚度矩阵KE的位置应为： 单元行码： IH=2*(I-1)+II 单元列码: L =2*(J-1)+JJ 此元素在整体刚度矩阵KS中的位置为： 整体行码： ID=2*(LND(IE,I)-1)+II 整体列码: IL=2*(LND(IE,J)-1)+JJ 此元素在半带存贮的整体刚度矩阵KS *中的 位置为 半带行码： IDH=ID 半带列码: IDL=IL-IDH+1 8子

16、框图7（EQUPE） 形成结点荷载向量P 8子框图7（EQUPE） 形成结点荷载向量P 第一步：将作用在结点上的荷载送入P中 10 有限元讲义 第二步：将作用在单元上的荷载转化为等效结点力后送入P 单元自重引起的等效结点力(竖向力)： PE=-V(材料容重)*AE(单元面积)*T(厚度)/3 11 有限元讲义 9子框图8（INSCD） 引入支承条件修改总刚 9子框图8（INSCD） 引入支承条件修改总刚 子框图8对应于子程序INSCD(子程序8)。其功能是引入支承条件, 引入支承条件就是 对带形矩阵KS *和向量P进行修改。 设第I个支承结点的结点码为IR,即IR =JR(I,1),又设支承结

17、点IR的约束位移所对应的 位移分量码为II。现分别说明对KS、KS *和P的修改内容。 1.对矩阵KS来说,第II行的主对角线元素应改为1,该行的其它元素改为0.此外,II列的 全部非主对角线元素也改为0,如图所示。 2.对半带矩阵KS *(程序中仍用KS表示)来说,第II行中第一个元素应改为1,该行其它元 素都改为0。此外,以该行第1个元素为起点,向右上方画45度斜线,则此斜线上的元素也都 改为0,如图所示。斜线上的元素如果列码是JJ,则其行码为II一JJ+1。 为了确定斜线元素的最大列码JO,需分两种情况考虑： 在图a中,行码II半带宽NW,故JO =NW, 在图b中,行码II半带宽NW,

18、 故JO =II。 总括起来,最大列码JO应等于NW和II两个数中的较小值 3.将P中对应行的元素改为零。 其框图如下： 12 有限元讲义 10子框图9 10子框图9 消元法解方程(略) 11子框图10（SIGME）求应力 11子框图10（SIGME）求应力 子框图10对应于子程序SIGME（子程序10）。其功能是计算单元应力。 上面解方程时，得出结构的节点位移，并存于P中。首先从P中取出单元e的节点位移 a e（程序中用DE表示），然后根据单元的节点位移ae求单元应力。求应力的基本公式 是 eee aBDaS= 求矩阵D，B时，要分别调研子框图4，5 关于主应力和主平面角的计算公式，可以参考

19、下图的应力圆得出： 平均应力： 2/ )( yx ASG+= 应力圆半径： 2 2 4 )( xy yx RSG + = 最大主应力：SGMA=ASG+RSG 最小主应力： SGMI=ASG-RSG 13 有限元讲义 主平面角： = = SGMI arcctg SGMI arcctg y xy y xy 29578.5790 180 0 0 下面写出详细框图。 14 有限元讲义 算例： 1.输入数据 1）主程序中输入结构控制参数 结点总数 NJ=10； 单元总数 NE=9 支承节点数 NR=2; 结点荷载数 NPJ=1 15 有限元讲义 问题类型码 IPS=0 （2）子程序1中输入数据 弹性模

20、量 E=1.0; 泊松比 PR=0.25 单元厚度 T=1.0； 容重 V=0.0； = 1098 896 976 865 562 673 743 632 521 39 LND 支承节点数组JR 节点荷载数组PJ 16 = 114 111 32 JR 0 . 100 0 . 10 31 = PJ 2. 计算半带宽NW 相邻节点码最大差值为4，NW=2*(4+1)=10 3.程序标识符说明 (1)整型变量 (1)整型变量 NJ-结点总数； NE-单元总数； NS-支承结点数 NPJ-有荷载作用的结点数, IPS-平面问题类型码(平面应力问题IPS=0，平面应变问题IPS=1) NJ2-位移分量总

21、数,NJ2=NJ*2, NW-半带宽；IE-单元序号。 (2)实型变量 (2)实型变量 E-弹性模量； PR泊松比； T-单元厚度； V-材料容重 有限元讲义 AE-单元面积 SGX,SGY,TXY-应力分量 ASG,RSG-平均应力,应力圆半径, SGMA,SGMI-最大与最小主应力, CETA-主平面角。 （3）数组 （3）数组 X(NJ)，Y(NJ) -结点坐标数组, LND(NE,3)单元的节点码数组 PJ(NPJ,3) - 结点荷载数组, JR(NS,3)支承节点信息数组 D(3，3)- 弹性矩阵D B(3，6) - 位移-应变转换矩阵（几何矩阵B） S(3.6) - 位移-应力转换

22、矩阵（应力矩阵S） ST(3) - 单元应力数组; DE(6) - 单元结点位移向量, KE(6,6) - 单元刚度矩阵K e KS(NJ2,NW) - 整体刚度矩阵半带存宽 P(NJ2)- 荷载向量、解方程后存放结点位移。 2.4 矩形单元 2.4 矩形单元 如图示深梁,设划分成具有八个矩形单元的网格，从中取出一典形单元e，单元有个 结点(i,j,m,p)。 一、局部坐标 一、局部坐标 设单元的边界平行于结构整体坐标x, y， 边长分别为2a, 2b, 结点编号从左下角开始逆时 钟编成i，j, m, p。 单元局部坐标的原点取在单元的形心(x0, y0)上,并采用无因次坐标,。,轴分别平 行

23、于x, y轴(如图),局部坐标,与结构坐标x, y的关系是: () () 0 0 yy xx 1 1 b a = = ()142 x0 , y0 为坐标原点。由此得到个结点的坐标分别是（-1, -1)、（1,-1）,(1,1),(-1,1)。 17 有限元讲义 二、结点位移列阵和结点力列阵 二、结点位移列阵和结点力列阵 每个结点2个位移分量，共个位移分量，设结点位移和结点力列阵分别为： T p ()242 pmmjjii e VUVUVUVUd= T p Y()342 pmmjjii e XYXYXYXF= 三、单元位移函数和形函数 三、单元位移函数和形函数 单元共有个位移分量， 将结点位移分

24、量全部作为已知边界条件， 则位移函数可取为： yxv yxu 765 321 xy xy 8 4 += += ()442 它是单值连续的。在平行于x轴的直线上,位移分量是x的线性函数；在平行于y轴的直 线上它是y的线性函数，故 2-4-4称为坐标x, y的双线性函数。 将各结点坐标式代入上式 （如u(-a， -b)=ui ） ， 通过解方程组便可求得待定常数 1 ， 2 8 ，将这些参数代回式（2-4-4），经整理得： ()()()() p U pmmjjii NUNUNUNu,+= ()()()() p V pmmjjii NVNVNVNv,+= (2-4-5) 以矩阵形式表示： (2-4-

25、6) e pmji pmji d NNNN NNNN v u f 0000 0000 = = e dN= 形函数矩阵N中的元素： ()=1 4 1 i N()1 ()+=1 4 1 j N()1 ()+=1 4 1 m N()+1 ()=1 4 1 p N()+1 18 如果引进参数: i = 0 , i = 0i (i=i, j, m, p), ( , j )是矩形单元4个结点的 局部坐标。 结点i( i , i )的坐标值分别是 i(-1,-1), j(1,-1),m(1,1), p(-1,-1)。 )， 代 入上式，则可将上式简记成： ()()( iii N+=11 4 1 ,) ()7

26、32()pm, ji, 或() 0 11 4 1 +=() 0 + 矩形单元的形函数，具有与三角形单元类似的特点，即： 在结点上的形函数值： 有限元讲义 19 i Ni( i )=1 Ni j ( j )=Ni( m m )=Ni p ( p )=0 即, Ni + Nj + Nm + Np =1 在单元形心上, Ni(0, 0)=1/4 (i,j,m,p) 四、几何矩阵B,应力矩阵S 四、几何矩阵B,应力矩阵S 在三角形单元中，我们已经推知： 应变： e dB= 应力: e dS e dBD= 矩形单元自然具有相同关系，只是B, S的内容有所区别。 仿照三角形单元： pm B 00 00 0

27、 0 ji pi pi BBB NN NN xy y x NHB= = L L (2-4-8) 式中: i i N 0 = i N ba a b ab B 0 0 0 1 () () ()() ji ba a b ab iiii ii ii + + + = 11 10 01 4 1 ()pm 单元应力矩阵： p S mji SSSBDS= 式中： () ()() ()() ()() i ba ab ab ab E S iiii iiii iii i , 1 2 1 1 2 1 11 11 14 2 + + + + = ()pmj, dVB 五、单元刚度矩阵 五、单元刚度矩阵 将上述公式代入单刚

28、的一般表达式 DBK T v = 对于矩形单元, 有限元讲义 20 dd BDBabtdxdyBDBtK T T a a b b = 1 1 1 1 计算可得按结点分块后的子矩阵为: () sr 3 1 + + + + + + + + = srsrsr srsr srsrsrsrsrsr rs a b b a b a a b Et K 1 2 1 3 1 1 2 1 2 1 3 1 1 2 1 3 1 1 14 2 ()pmjisr, = 2.5 六结点三角形单元 2.5 六结点三角形单元 一、面积坐标 一、面积坐标 在开始单元推导之前，先介绍有关面积坐标的概念。 1面积坐标的定义 在第二节中

29、，我们讨论过三角形单元形函数的几何意义，如果从三角形单元e内任意 点p向三个顶点i,j,m 引连线,使其分割成三个小三角形便可得到这些小三角形面积与形函 数的关系。 A A N i i = A A N j j = A Am =Nm 反过来说，只要给定了三角形面积 Ai , Aj , Am， p点即可确定。因为 Ai, Aj ,Am 中 只有两个是独立的(AiAjAm=A)，这与P点坐标(x, y)对应。因此，我们可以用三角形面 积来定义三角形内的任意一点P的位置，并称它为面积坐标。P点的面积坐标定义为以下三 个面积比: A A L i i = A Aj Lj= A A L m m = 由此定义

30、,必有：Li +Lj +Lm =1。由于面积坐标Li , Lj , Lm 只限于用来确定三角形单 元上(内)点的位置（只有当L(i,j,m)1才有定义）,并且只是在单元分析中使用,所以L(i,j,m) 属于局部坐标。 （直角坐标与坐标原点、平移转动等有关，可因人而异，故称它为人为坐标，而面积坐标 具有不变性，故称其为自然坐标。） 2面积坐标Li和整体坐标x,y之间的关系 面积坐标Li和整体坐标x,y之间的关系为: () ()xxyyyxyx A yx yx yx AA A L mmjjmmj mm jj i i += 2 1 1 1 1 2 1 ()yxj ()mj, i 有限元讲义 21 或

31、: ()ycxba A L iii += 2 1 i ()mji, ()152 将上式与(2-1-3) Ni =1/2A(ai+bix+ci y)比较得知,三角形 常应变单元的形函数Ni , Nj , Nm 分别等于三角形面积坐 标Ni , Nj , Nm ,因此在三角形的三个顶点上亦具有与形函 数相同的结点坐标: () = =yxL ssr 0 1 , sr sr 三角形内部点和边界点上的面积坐标值可以用右图 中的一些等值线表示（面积坐标的等值线）。从图中可看 出：在平行于jm边的任何一条直线上各点都具有相同的Li 值(可从几何作图上看出)其值等于该直线到对应jm边界 的距离与结点i到jm边

32、界距离的比值 其值等于该直线到对应jm边界 的距离与结点i到jm边界距离的比值。其它亦然。 3面积坐标与直角坐标的关系 将式(2-5-1)写成矩阵形式得到用x、y表示的面积坐标: y x 1 = cba cba cba A L L L mmm jjj iii m j i 2 1 ()252 其反变换式为: m j i L L L = mji mji yyy xxx y x 1111 ()352 上式的数学意义是:三角形ijm内某点的坐标x,y,可以用三角形顶点(3个结点)的坐标 (xi , yi ), (xj , yj ), (xm , ym )进行函数插值而获得,其插值多项式是三角形的面积坐

33、 标Li , Lj , Lm， 也就是三角形的形函数Ni , Nj , Nm 。这种坐标（x, y）表达式的插值 函数和位移场 （u, v） 表达式的插值函数采用同一函数的单元称为等参单元 坐标（x, y）表达式的插值 函数和位移场 （u, v） 表达式的插值函数采用同一函数的单元称为等参单元(简称等参元)。 （又称：坐标变换函数和位移插值函数中的形函数相同）。 “等参元”的另一种定义是：（参考书：湖南科技社有限元法概论 ） 设 位移变量u的近似式: = = s i u 1 iiu ( ) i x u是用s-1次插值函数近似表示的。坐标变换式： r i i x = = 1 通常,坐标变换式的次

34、数与u 近似式的次数并不相等，由此可将单元分成三类。 有限元讲义 、子参元 rs 又名 亚参元（Sub-parametric element） 、等参元 r =s (isoperimetric element) 、超参元 r (Super-parametric element) 亚参元概念： 有时为了减少结点坐标数据的输入量，用于位移描述的结点数和用于坐标变换的结点 数可以不一样。如图示四边形单元，可将位移取为8个结点,即位移函数为 22 iiu N iiv N iix f iiy f T v3 T y3 2 6y 2 12y i u = = 8 1 i v = = 8 1 而坐标变换式则取为

35、: ijmp x = ijmp y = 我们称这种坐标结点少,而位移结点多的单元为亚参元 二、六结点三角形单元 二、六结点三角形单元 在三角形单元三条边界的中点上各增设一个结点，使三角 形单元共有12个位移自由度。 1、结点位移列阵： (2-5-4） mmii e vuvuvud 11 LL= 结点力阵: （2-5-5） mmii e yxyxyxF 11 LL= 、单元位移函数和形函数 设为次多项式： 2 54321 xxyyxu+= 2 1110987 xxyyxv+= 为了计算简单,采用面积坐标Li, Lj ,Lm(Li=1/2A(ai+bix+ciy) (i,j,m) ,单元上六个结点

36、的 面积坐标计算如图示。 将六个结点的12个位移分量代入， 并采用面积坐标， 按前面类似的推导过程，便可得到用形函数表达的单元 有限元讲义 位移函数： 23 N 33v NvNvNv ijii33u uNuNu jjii +=L =+L e dN （2-5-6） 写成矩阵形式: i ji ji v u NNN NNN v u f= = = 3 3 3 000 000 L L 式中形函数: ()12= ii LN i Lmji, mjL L4()imj, () N1= , 3 , 2 , 1 、几何矩阵B, 应力矩阵S 应变=Bd x u x = y v y = x v y u xy + = 面

37、积坐标函数的运算: 设函数 Z= f(Li Lj Lm )， 由于面积坐标Li, Lj Lm 都是x,y的函数。 因此将面积坐标对直角坐标求导时，应按复合函数求导法则： x Lm m L z x L L z x L L z x z j j i i + + = 同理: y Lm m L z y L L z y L L z y z j j i i + + = 由(2-5-1) (Li=(ai+bix+ciy)/2A) (i j m) A b x Li 2 = ()mji i 所以，面积坐标函数对直角坐标的偏导数为: m m L z + + = j j i i b L z b L z b Ax z

38、2 1 ()752 + + = m m j j i i L z c L z c L z c Ay z 2 1 同理: () () i L A c y N L A b L z A b x N i ii i i i ii 14 2 14 22 = = = ()mj 有限元讲义 () () ()imj A LcLc y N A LbLb x N jmmj jmmj ,321 2 4 2 4 1 1 + = + = 由此得： ()()() m u1 mmjjjiiix LbuLbuLb Ax u 41414 2 1 += =+ +()()() 3 uLi 21 444bLbuLbLbuLbLb jji

39、miimjmmj + 24 同理可得y和 xy 的类似表达式，于是得到六结点三角形单元的几何方程: d()852BdBBBBBB mji = 321 = 式中： () () () () = 414 40 014 2 1 iiii ii ii i LbLc Lc Lb A B 1 1 ()mji () () () j j L L () + + + = mmjjmmj mmj jmmj bLbLcLc cLc LbLb A B 44 40 04 2 1 1 ()i m j 3 2 1 dSdBDD xy y x = = 单元应力方程: ()952 应力矩阵： 3 S 21 SSSSSS mji =

40、 式中: () () ()() i i b = i ii i i i c cb cb A LE S 11 22 22 14 14 2 ()mji () ()() ) () jm j j L L ()imj321 ()( ()()() + + + = mjjmmj mmjjmmj mmjjmmj cLbLcLc LcLcLbLb cLcLbLb A E S 1414 88 88 14 2 1 从（2-5-8、9）两式可见，六结点三角形单元的应力和应变分量均为整体坐标 x,y或 面积坐标L的一次函数,因此它比常应变单元精度好。而且，它们沿任何方向又都是线性变 化的，能满足边界的变形协调条件，即单元的两条收敛准则,完备性和协调性都能得到满 足。 、单元刚度矩阵 有限元讲义 单刚为1212阶，由一般公式K