有限元2-弹性力学平面问题(23程序设计)

1、土木工程学院有限单元法有限单元法P-1/57第2章 弹性力学平面问题有限单元法2.1 三角形单元三角形单元 2.2 三角形单元中几个问题的讨论三角形单元中几个问题的讨论 2.3 平面问题有限元程序设计平面问题有限元程序设计 2.4 矩形单元矩形单元 2.5 六结点三角形单元六结点三角形单元2.6 四结点四边形单元四结点四边形单元 2.7 八结点曲线四边形等参元八结点曲线四边形等参元 2.8 几个问题的补充几个问题的补充 土木工程学院有限单元法有限单元法P-2/572.3 平面问题有限元程序设计平面问题有限元程序设计 土木工程学院有限单元法有限单元法P-3/57一、程序设计方法与结构分析程序的特

2、点一、程序设计方法与结构分析程序的特点 1程序设计方法论简述程序设计方法论简述 借助计算机来完成某项工作，通常都要先编写相应的计算机程序，或叫程序设计。完成一个结构分析或结构CAD系统也必然要经过程序设计才能实现。程序设计要使用专门的程序语言。我国结构程序设计中所采用的语言，在60年代和70年代初以ALGOL语言为主。此后逐步广泛使用的主要是BASIC语言和FORTRAN语言，随着CAD和人工智能技术的发展，PASCAL、 C、LISP、 PROLOG等有着各自特长的程序语言也逐步进入土木工程领域的计算机程序设计中。土木工程学院有限单元法有限单元法P-4/57过去人们通常认为，程序设计的中心问

3、题就是学会使用一种程序语言，用以编写程序。然而学会用程序语言编程只是整个程序设计中的一部分。据有关资料介绍，编写程序在整个系统的研据有关资料介绍，编写程序在整个系统的研制过程中仅占制过程中仅占15的工作量。在一个大型程序的工作量。在一个大型程序系统的整个存在阶段的工作量中，在系统投入使系统的整个存在阶段的工作量中，在系统投入使用后的维护工作量为原来研制工作量总和的两倍用后的维护工作量为原来研制工作量总和的两倍。维护工作量是如此之高，这就使我们必须注意到，在程序研制阶段便即应当考虑为以后的维护工作提供方便，哪怕是为此要增加一些额外的工作量也是值得的。 土木工程学院有限单元法有限单元法P-5/57

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

5、动的费用。总的来说，程序设计已从强调灵活总的来说，程序设计已从强调灵活的技巧和局部效率向着强调程序结构化和整体功的技巧和局部效率向着强调程序结构化和整体功能的方向发展。能的方向发展。这实际上就是逐步发展起来的关于程序设计编写与调试的一套方法论，其要点可归纳为以下几方面： 土木工程学院有限单元法有限单元法P-7/57(1) 编程结构化编程结构化为了使程序设计者能按照一定的结构形式，而不是随心所欲随心所欲地设计和编写程序，使编制的程序易读、易修改，以提高程序设计和维护工作的效率，荷兰学者Dijkstra提出了“结构化程序设结构化程序设计方法计方法”。结构化程序规定了三种基本的结构形式，他们是顺序结

6、构、分支选择结构和循环结构。编程结构化又称结构化程序设计，他可使编写的程序层次分明，逻辑清楚，容易阅读。 土木工程学院有限单元法有限单元法P-8/57 (2) 分层处理技术分层处理技术为了解决现实世界中的许多复杂问题，人们往往需要根据问题的内在联系将其分割成有层次的一系列问题来分别求解。对于一个大型程序系统设计来说，也需采用分层的办法来处理，在每一层里集中解决一个问题，并为下一层的执行作好准备。分层处理技术的主要内容是将程序划分为多个层次的若干模块，每个模块完成一个（或几个）预定的功能。为了保证模块的独立性，各模块之间只能通过接口接口与其他模块连接。另外，对于一个较大的软件系统要由多人合作才能

7、完成，模块化也为此提供较好的合作条件。模块化的一种典型结构形式就是子程序结构。土木工程学院有限单元法有限单元法P-9/57(3) 避免过多使用GOTO语句，特别是逆转的GOTO语句。这是结构化程序设计的基本要求之一。（4）可移植性 对于应用软件，特别是大型的CAD软件，可移植性高低同样是衡量软件的质量的重要指标。因为研制一个大型应用软件系统不但要耗费大量的人力物力，而且还要花费相当长的时间才能研制成功。在其研制和应用期间，不可避免地发生运行环境的变化计算机硬件设备的换代和系统软件的更新。可移植性强即等于延长了软件的生存期，从而可节省新的开发投资。土木工程学院有限单元法有限单元法P-10/57可

8、移植性主要体现在软件对支撑环境的独立性和软件本身的封闭性。例如，用FORTRAN、BASIC 、C等语言编写，并尽量避免采用非标准语句和函数。此外，在软件中采用统一的I/O 模块，也是提高可移植性的手段之一。应当指出，程序设计方法论仍在发展探索之中，千万不能把上述有关内容当成一成不变的教条套用，而应当通过实践来发展和丰富其内容。然而程序设计发展到今天，已经奠定了很多必要的理论基础。我们正在达到一个可以谈论程序设计是一门科学而不仅仅是一种技巧的阶段。 土木工程学院有限单元法有限单元法P-11/572、结构化程序设计、结构化程序设计结构化程序设计，又称结构程序设计(Structured Progr

9、amming)是 荷兰学者E. W. dijkstra首先提出来的，简称SP。人们对SP有各种定义和解释:有人说它是: 指导程序员编程的一般方法；有人说它是: 不使用goto语句的程序设计;有人说它是: 自顶向下的程序设计。有人把层次结构、顺序结构、选择结构和重复结构定义为SP的全部内容，而把一个程序结构模块定义为上述四种基础结构的某种集合，即：结构模块(程序)(层次结构, 顺序结构, 选择结构, 重复结构)，认为SP就是定义为这些基本结构在软件开发和维护中的严格应用。土木工程学院有限单元法有限单元法P-12/573、程序的正确性及其验证、程序的正确性及其验证程序正确性(Program Cor

10、rectness)被定义为一个程序和它打算要实现的功能之间的一次符合(a correspondness) Gries说，在程序设计的初期阶段，人们很少看到程序正确性方面的问题，那时人们往往着重于调试(debugging),但是后来发现调试好的程序并不能表示错误不存在。于是人们就不得不考虑程序的正确性证明问题正确性证明问题。 土木工程学院有限单元法有限单元法P-13/57什么是程序正确性证明呢？从词义上讲，“证明”就是提供一种有力的证据使人们在思想上不得不接受一条真理或一件事实。但是在程序设计中要做到这一点并不容易，因为现在还没有找到一种象数学一样非常完整的形式。目前主要还是通过选择较好的算法，

11、大量考题，与已有结果比较等办法来证明一个程序的正确性。但不管怎样，程序正确性理论已经被提出来，并有了一些实际想法。 土木工程学院有限单元法有限单元法P-14/57程序的验证(Program Verification)实际上就是检验程序的正确性。一个程序如果有错误，主要是两方面的：一方面是语法错误，这部分比较好解决，一般是在调试阶段（编辑阶段）完成。但是一个语法完全没问题的程序并不一定是正确的。因为程序中的许多部分往往是靠逻辑关系来达到所要实现的目的。目前应用程序的验证主要靠针对程序每一功能，每一逻辑分支进行各种类型的考题，包括考题的规模。土木工程学院有限单元法有限单元法P-15/574、结构分

12、析程序的特点、结构分析程序的特点1. 规律性、 通用性好。2. 计算工作量大、运算时间长（形成单、总刚、解方程、动力、非线性）。3. 矩阵阶数高，与微机的存贮量发生矛盾。土木工程学院有限单元法有限单元法P-16/57二、平面问题有限元程序设计（三角形单元）二、平面问题有限元程序设计（三角形单元） 三角形单元是平面问题中最简单的单元形式，但其他单元形式的有限元程序与三角形单元程序的结构基本相同，并无本质区别。故本节通过一个FORTRAN 源程序，详细介绍三角形单元的程序设计。下面总框图左边主线对应着源程序中的主程序流程。现按照总框图中的10个子框图，结合详细的子框图设计，介绍 “平面问题三角形单

13、元源程序”。土木工程学院有限单元法有限单元法P-17/57输入结构控制参数输入其它数据 形成整体刚度阵 引入支承条件 解方程，输出位移 求应力，输出应力形成节点荷载向量 开始结束1单元面积 求弹性矩阵 单元刚度矩阵位移-应变矩阵2345678910土木工程学院有限单元法有限单元法P-18/572子框图子框图1（SDATA）输入其他数据并计算半带宽输入5个控制参数(结点数、单元数、支承数、荷载数、类型)后，程序运行中所需的其他原始数据均放在该子程序中输入。其框图如下：土木工程学院有限单元法有限单元法P-19/57土木工程学院有限单元法有限单元法P-20/57土木工程学院有限单元法有限单元法P-2

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

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

16、元法有限单元法P-24/57土木工程学院有限单元法有限单元法P-25/573子框图子框图2（STE） 计算单元刚度矩阵KE两个矩阵相乘土木工程学院有限单元法有限单元法P-26/57须调用子程序3，4，5，输出参数KE土木工程学院有限单元法有限单元法P-27/574子框图子框图3（AE） 计算单元面积AE对应子程序ATE，被子框图2，7，10调用土木工程学院有限单元法有限单元法P-28/575子框图子框图4（DTE） 计算弹性参数矩阵D土木工程学院有限单元法有限单元法P-29/576子框图子框图5（BTE） 计算几何矩阵B土木工程学院有限单元法有限单元法P-30/57土木工程学院有限单元法有限单

17、元法P-31/577子框图子框图6（STIFF） 形成总刚度矩阵KS节点循环自由度循环土木工程学院有限单元法有限单元法P-32/57单元列码整体列码半带列码土木工程学院有限单元法有限单元法P-33/57(1) 结点局部码与总码的对应关系（结点换码） 图示为单元3个结点的两种编码，某个结点的局部码如果是i(1),则它在总刚中的对应总码是LND(IE,1)。 土木工程学院有限单元法有限单元法P-34/57（2）同一子块在单刚KE和总刚KS的对应位置（子块搬家）(a)矩阵KE土木工程学院有限单元法有限单元法P-35/57（3）同一元素在KE、总刚KS和带状总刚KS*中的对应位置（元素搬家） 形成总刚

18、中某个元素的搬家问题： 设该元素属于单刚矩阵KE中的I行J列子块,该元素在这个子块中的位置是II行JJ列。则此元素在单元刚度矩阵KE的位置应为：单元行码：IH=2*(I-1)+II单元列码: L =2*(J-1)+JJI行J列子块KE土木工程学院有限单元法有限单元法P-36/57此元素在整体刚度矩阵KS中的位置为：整体行码: ID=2*(LND(IE,I)-1)+II整体列码: IL=2*(LND(IE,J)-1)+JJ土木工程学院有限单元法有限单元法P-37/57此元素在半带存贮的整体刚度矩阵KS*中的位置为半带行码：IDH=ID半带列码: IDL=IL-IDH+1土木工程学院有限单元法有限

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

20、设支承结点IR的约束位移所对应的位移分量码为II。现分别说明对KS、KS*和P的修改内容。 土木工程学院有限单元法有限单元法P-41/57（1）对矩阵KS来说，第II行的主对角线元素应改为1，该行的其它元素改为0。此外，II列的全部非主对角线元素也改为0，如图所示。土木工程学院有限单元法有限单元法P-42/57（2）对半带矩阵KS*(程序中仍用KS表示)来说，第II行中第一个元素应改为1，该行其它元素都改为0。此外，以该行第1个元素为起点，向右上方画45度斜线，则此斜线上的元素也都改为0，如图所示。斜线上的元素如果列码是JJ，则其行码为II-JJ+1。为了确定斜线元素的最大列码JO，需分两种情

21、况考虑：在图a中,行码II半带宽NW，故JO =NW,在图b中,行码II半带宽NW， 故JO =II。总括起来，最大列码JO应等于NW和II两个数中的较小值。土木工程学院有限单元法有限单元法P-43/57图a：行码II半带宽NW,故JO =NW图b：行码II半带宽NW, 故JO =II土木工程学院有限单元法有限单元法P-44/57（3） 将P中对应行的元素改为零。其框图如下：II=2*(IR-1)+J-1土木工程学院有限单元法有限单元法P-45/57土木工程学院有限单元法有限单元法P-46/5710子框图子框图9 消元法解方程(略)11子框图子框图10（SIGME）求应力 子框图10对应于子程

22、序SIGME（子程序10）。其功能是计算单元应力。上面解方程时，得出结构的节点位移，并存于P中。首先从P中取出单元e的节点位移ae（程序中用DE表示），然后根据单元的节点位移ae求单元应力。求应力的基本公式是 eeeaBDaS求矩阵D，B时，要分别调研子框图4，5土木工程学院有限单元法有限单元法P-47/57关于主应力和主平面角的计算公式，可以参考下图的应力圆得出：平均应力：应力圆半径：最大主应力：SGMA=ASG+RSG最小主应力： SGMI=ASG-RSG主平面角：2/ )(yxASG224)(xyyxRSGSGMIarcctgSGMIarcctgyxyyxy29578.5790 18000土木工程学院有限单元法有限单元法P-48/57详细框图 土木工程学院有限单元法有限单元法P-49/57土木工程学院有限单元法有限单元法P-50/57土木工程学院有限单元法有限单元法P-51/57