有限元第9章有限元法程序设计.docx
《有限元第9章有限元法程序设计.docx》由会员分享，可在线阅读，更多相关《有限元第9章有限元法程序设计.docx（75页珍藏版）》请在冰豆网上搜索。
有限元第9章有限元法程序设计
第9章 有限元法程序设计
9．1 引言
在用有限元法进行结构分析时，将会遇到大量的数值计算，因而在实用上是一定要借助于计算机和有限元程序，才能完成这些复杂而繁重的数值计算工作。
事实上，有限元程序的设计是有限元研究的一个很重要的部分。
它是理论和方法的载体，是理论用于实际必不可少的桥梁，是有限元学术研究与实际应用水平的代表。
有好的、高深的理论和算法并不等于有好的程序，还必须有实际的程序开发经验的多年积累、丰富的计算机知识、大量的资金和人力的投入，多年的开发修正与改进才能编制出好的程序来。
一些著名的有限元程序开发的发展历史也体现出了这一规律。
设计一个用于结构分析的有限元法程序，要求设计者至少应该掌握下列知识：
（１）掌握一种程序开发工具，如VC（VisualC++），CB（C++Buildel），Delphi，VB（VisualBasic）或VF（VisualFortran）等。
在本书中所有程序均用VC写出。
（2）数值方法，如线性和非线性代数方程的求解，矩阵特征值的求解以及数值积分等。
（3）结构分析的基本理论，特别是用有限元法对结构进行分析的原理、方法和步骤。
由于一般的软件工程师不懂结构分析原理，因此，结构分析程序的开发任务主要应由结构工程师来承担。
掌握结构分析程序设计方法，是以计算机辅助设计为主要标志的现代工程设计方法对结构工程师的要求。
作为结构工程师，应该具有对结构分析程序的使用、阅读、修改和编制的基础知识和技术素质。
有限元程序的总体组成可分为三个部分：
前处理部分，有限元分析本体部分以及后处理部分。
有限元分析本体部分是有限元分析程序的核心。
它根据离散模型的数据文件进行有限元分析，有限元分析的原理和采用的数值方法集中于此。
因此，这一部分程序是有限元分析是否准确可靠的关键部分。
有限元分析所使用的离散模型的数据文件主要包括：
模型的节点数、节点坐标与节点编码，单元数据与单元编码；材料和载荷信息等。
实际工程问题的离散模型数据文件十分庞大。
一般情况下，用人工方法来生成工作量太大，并且容易出错，有时甚至是不可能的。
为解决这一问题，有限元程序必须有前处理程序。
前处理程序根据使用者提供的对计算模型外形及网格要求的简单数据描述，能自动地或半自动地生成离散模型的数据文件，并能绘制结构计算简图和网格图，供用户检查修改。
前处理程序的功能在很大程度上决定了有限元程序使用的方便性。
有限元分析程序的计算结果是由离散模型而得到的，输出的数据量往往很大，不易整理，也不易获得分析对象的全貌。
所以，一个使用方便的有限元分析程序还应具有较强的后处理功能。
能够按照用户的需要提供应力分级图、等值线图，结构变形图或振型图等图形显示功能，以及按照一定的要求对计算成果进行列表显示或打印。
因此，这部分程序设计的好坏，对整个有限元程序使用起来是否方便，具有举足轻重的作用。
程序设计工作经历了纯技巧阶段，已经形成了一门被称为软件工程的学科。
对于程序的质量评价也逐渐形成了一套客观标准。
一个质量较高的程序应该具有较好的可管理性和较高的运行可靠性。
可管理性要求程序的可读性好，易于调试、修改和发展，使用方便且效率高等。
可靠性要求程序能正确无误地完成规定的功能，当出现不正常情况时，能中止无价值的运行并输出有关的信息。
程序开发的过程大致可分为下述三个阶段：
（1）程序功能的规定；
（2）程序结构的设计，源程序及其说明的编写；（3）调试和纠错。
目前在实际的程序开发中，流行着两种截然不同的方法，即面向过程的方法和面向对象的方法。
大量的资料说明，在开发大型应用软件时，面向对象的方法与传统的过程化程序设计方法相比，显示出很大的优越性。
然而在开发一些规模不大的中小型程序时，面向过程的方法仍然有一定的优势。
本章将以平面杆系结构的静力分析为例，介绍用面向过程的方法进行有限元主体程序设计的方法。
9．1．1 结构化程序设计概述
结构化程序设计方法是一种传统的软件设计方法。
其基本要点是，自顶向下，逐步求精，以及模块化设计。
其基本思想是，把一个复杂问题的求解过程划分成若干阶段来进行。
每一个阶段所要解决的问题都控制在人们容易理解和处理的范围内，直到把原来的问题变换成若干个易于编写成程序的子问题（即模块）为止。
这种程序的逐步分解和精化是从抽象的做什么到具体的怎么做的发展过程。
程序展开的基本结构有下列三种：
（1）顺序结构。
把一个功能块展开成若干个顺序连接的语句块。
用元语言（即程序设计语言和自然语言的一种混合文体，也称伪语言）表示如下：
语句块1；//入口
语句块2；
…………
语句块n；
//出口
（2）选择结构。
把一个功能块展开成两个可供选择的语句块。
用元语言表示如下：
if（exp）//入口
语句块1；
else
语句块2；
//出口
在C/C++中，Switch语句组成的结构也属于这种结构。
（3）循环结构。
把一个功能块展开成需要重复执行的语句块。
用元语言表示如下：
while（exp）//入口
语句块；
//出口
在C/C++中类似的结构还可由dowhile语句以及for语句构成。
上述程序展开的三种基本结构表达了程序内部语句块的执行次序。
一个语句块可以是空的，可以是只含有一条简单语句，也可以是由一系列的语句组成的子问题块。
一个语句块只有一个入口点和一个出口点。
在一般情况下尽可能少使用goto之类的语句，最好不使用它。
这是因为这类跳转语句在程序中的频繁出现，会极大地破坏程序的逻辑清晰度，使程序难以阅读和理解。
结构化程序设计的另一个重要技术就是模块化编程。
逐步求精的结果是以子功能块为单位的算法描述。
以子功能块为单位进行程序设计，实现其求解算法的方式称为模块化设计。
模块化方法可以有效地简化问题的复杂性和提高程序的正确性。
模块化的准则是高独立性、高内聚性、低耦合性。
自顶向下法和自底向上法是两种最基本的模块化编程方法。
前者是逐步分解和求精的过程，而后者则是逐步组合的过程。
自顶向下法是从整个程序的功能描述开始的。
即首先把总功能分解为若干个较小的，然而相对来说较为简单的功能块。
连续不断地这种分解，将会发现有些功能块已经足够简单了。
于是就可以把它改写成程序。
如果用元语言的程序文本来描述这一设计过程的话，则此时程序文本中的一部分已可以用程序语言表示了，而另一部分则仍然用元语言表示。
随着不断地精化，最后将会把用元语言表示的整个设计文本全部由程序设计语言来表示，从而完成程序的编写工作。
上述用元语言表示程序流程的方法便于计算机存贮与操作，而且对表示功能的自然语言（做什么），在编写完相应的程序段（怎么做）之后，加上注解标记就可以成为程序中的注解行。
元语言表示正日益广泛地被程序设计人员采用。
最终的程序结构可以用一树状图来表示。
如图9.1所示，功能A被分解成B、C、D三个模块，然后再依次分解为E、F、G、H、I、J。
此时模块E、F、G、H可以用程序语言编写，H、J分解为K、L。
接着K、L可以用程序语言来编写。
于是逐层修改完善程序文本，使说明程序功能的自然语言成为程序的注解行。
图9.1自顶向下的程序设计
自底而上的编程方法同自顶而下法正好相反。
自顶而下法中的最后一个功能分解正是自底而上法中的第一个功能组合。
不断进行这种由多个简单的功能块逐步组合成一个较为复杂的功能块的过程，直至完成全部编程工作。
自顶而下法同自底而上法相比，具有效率高、适应性强，容易控制模块的复杂程序，以及能增加编程工作中的能见度等特点。
因而通常被认为是更好的设计方法。
但是当遇到比较复杂的实际问题时，为使设计更有把握，应该优先考虑关键部分。
另外，当使用自顶而下法逐步分解功能块的同时，应该考虑在较低层次的功能块中充分利用自己已有的程序经验和知识，来组合所需要的功能块的可能性。
因此，在实践中灵活交替地使用这两种方法可以提高工作效率。
9．1．2 程序的可读性和程序风格
几乎所有的程序在使用过程中都要经受各种各样的维护性修改。
据统计，程序在研制完成后的这种维护性的工作量大大超过了它在研制阶段的工作量。
因而程序除了能正确地完成一系列的计算之外，还应该让维护人员容易看清其算法的逻辑联系，即应该具有良好的可读性。
程序的可读性好了，维护和修改起来就会很方便。
提高程序的可读性是以良好的程序设计风格来体现的。
所谓程序风格是指不同的程序员，在其编写的程序中，在代码文件，语句构造以及变量命名等方面所表现出的特点。
具体来讲程序风格体现在以下几个方面：
1．程序的模块化、结构化。
如前所述，程序的模块化是按照功能划分模块的设计方法，目的是使程序的层次和作用更明确，也易于修改和调试。
而程序的结构化是指按照结构化方法进行程序设计。
结构化程序应具有以下特点：
（1）全部程序均由上节所讲述的三种基本结构组成，不包含其它类型的结构；
（2）只有一个入口和一个出口；
（3）程序的执行是有限的，无死循环；
（4）无死语句，即每条语句均有被执行的机会。
由上可见，结构化程序实际上是由许多相对完整独立的程序段构成，每一段均只有一个入口和一个出口。
条理清楚，层次分明，即使程序长一些，也仍然易读可靠，容易验证，而且修改某一段一般不会影响其它段。
2．程序内部的文档。
程序内部的文档通常包括以下几方面：
（1）必要的注解。
在程序的关键部分，程序、过程、函数和段落的开头，一些公式之前等处，应加上简明扼要的注解和说明语句。
这些说明性语句应能提供某些对理解程序有用的东西，而不仅仅是代码的语言解释。
请看下面的一个注解：
for（i=0；i｛
…
｝
该注解只是把语句解释了一下，丝毫没有增加有利于读者理解的新信息，因此，它是多余的。
换成下面的一种写法才是有意义的：
for（i=0；i｛
…
｝
（2）正确地命名和使用标记符。
这体现在正确地使用变量的类型和变量名两点上，要充分利用程序设计语言本身提供的特点，根据不同需要，分别选用整、实、字符、枚举、结构（或记录）和数组等各种变量类型。
例如，为了方便地表示一周中的七天，可以这样定义一个枚举类型变量：
enumDay
｛Sunday,Monday,Tuesday,Wednesday,Tursday,Friday,Saturday｝；
这使人看起来很习惯，也很自然。
又如平面刚架的节点信息，可以选用结构体这种数据类型，用如下方式定义：
structNode
｛
doubledX;//横坐标
doubledY;//纵坐标
intiaDOFIndex[3];//节点的三个自由度编号
｝；
以上两类问题如果仅用整、实型变量表示就不直观，不容易记忆。
为了使程序易读，要尽量采用容易记忆且能反映自身特征的变量名，如力可用F表示，面积用A表示，惯性矩用I表示，角度用Alpha，Beta表示等等。
诸如此类的问题，要尽量采用人们熟悉的，接近自然语言的，并为有关学科广泛采用的符号来表示。
在Windows程序设计中，匈牙利命名法是一个广泛采用的命名法。
这种命名法是以下面两条规则为基础的：
规则1：
标识符的名字是以一个或者多个小写字母开头（前缀），用这些字母来指定数据类型。
例如以小写字母c作前缀表示字符型，而以p作前缀表示指针类型，n或i开头表示整型等等。
经常使用的前缀如下表所示。
表9．1 匈牙利命名法的前缀
前 缀
数据类型
a
数组
b
布尔型
c
字符型
d
双精度型
f
单精度型
i
整型
n
整型
p
指针
s
字符串
规则２：
在标识符内，前缀以后就是一个或多个首字母大写的单词，这些单词清楚地指出了源代码内对象的意义和用途。
例如，用nＬoad表示载荷数，用sFileName表示文件名等。
用匈牙利命名法命名的变量往往可以提供比较丰富的变量信息，但有时显得较罗嗦，因此本书对于某些局部变量和循环变量等并未采用以上规则。
（3）充分利用分隔符。
充分利用括号、空格、空白行和连接符等的作用，以便明显隔开不同的语句成分。
（4）正确采用缩进规则。
在程序中各种嵌套之处采用适当的缩进表示。
如内层语句较外层语句缩进若干空格；不同的函数段与函数段之间用空行或注解行隔开等，都可增强程序的层次感和逻辑清晰度。
（5）明确程序间的调用关系。
在程序的首部对该程序调用的各种函数、过程应作说明。
3．数据说明。
数据的说明应该标准化，以便查阅、测试、调试和维护。
每一个变量都应作说明，当多个变量在一个说明语句时，应该按字母顺序排列这些变量。
有人建议一行只写一个变量声明。
不同层次最好不要使用相同名字的变量。
4．语句的构造。
语句应当简单而直接，不要为了节省空间或所谓的提高效率而使语句排列拥挤，或使程序显得复杂别扭，晦涩难懂。
5．效率。
合理解决运行时间和内存的矛盾。
不要为提高效率而牺牲程序的清晰性和可读性。
为了提高程序的可读性和整体质量，类似以上所讲的这些最基本的“条条框框”应当加以遵守。
当然，一个良好的程序设计风格，需要在长期的编程实践中有意识地培养和训练才能形成。
9．1．3 程序调试
为了验证程序的正确性，有两种截然不同的极端方法可供采用。
一种是深入每一个模块内部，使“所有的工作都过一遍筛子”，通过对程序的结构和语句的系统研究来检验它的正确性，这种方法也称作白箱法；另一种验证程序正确性的方法是不需要对模块作深入了解，只须通过一系列试算题的实际计算结果来判断该模块是否正确，这种方法也称作黑箱法。
在理论上前者可以不通过计算机的运行来检验程序的正确性，而后者则不必了解程序的结构即能考核程序的正确性。
然而只有研究了程序的结构后，才能提出一个用少量的考试题，即能检验程序正确性的方案。
最小量考试便是一例，此时只要用一组考试题便能保证程序中的每一个语句至少被执行一次。
程序的调试工作是一项细致而繁重的工作，需要认真而耐心地去做。
调试一般分为下以几个步骤：
（1）分调（亦称单元调试）——分别调试每个模块。
（2）联调（亦称集成调试）——将模块组合起来调试。
（3）验收——确认程序正确。
模块调试（分调）分为渐增式和非渐增式两种方式。
渐增式是在调试好N个模块时，再调试N+1个模块。
非渐增式是分别独立调试每一个模块，然后将它们联合成一个系统。
调试顺序可以自顶而下，也可以自底而上。
哪一种更好呢？
采用自顶而下的顺序，能较早地调好系统整体，特别是软件接口上的问题能较早暴露，但是调试起来比较困难。
采用自底而上的顺序实施起来比较容易，但是软件接口问题暴露得较晚。
因此，这两种调试顺序各有优缺点，而且其中一个的优点正好是另一个的缺点，最好结合起来使用，根据程序结构的具体情况，灵活地加以确定。
调试中若发现错误要认真分析。
纠错需要一定的经验和技巧，必须勤思考，不可蛮干。
正确的做法是采用原因排除或回溯两种方法。
原因排除法是先经过思考，提出出错的可能原因，然后进行必要的测试，一个个排除，最后找出出错的真正原因。
回溯是从出错的地方一步一步向前查，直到查出出错的地方为止。
有时通俗地称这种方法叫“顺藤摸瓜”。
回溯法一般适用于小程序。
VC的集成开发环境（IDE）为程序员提供了一系列优秀的调试工具，熟练地使用这些调试工具对于提高工作效率是十分有利的。
9．2 平面杆系结构静力分析程序设计
9.2.1程序功能的确定和模块划分
本章将要介绍的平面杆系结构静力分析程序，是一个示例性的小程序，目的是用来说明有限元法程序设计方法，因此，并不追求程序功能的完整和强大，将程序的功能限定为：
可以用来对一般平面杆系结构作静力分析。
结构可以是刚架，桁架，或由刚架单元和桁架单元组成的组合结构。
载荷除一般载荷外，还可对温度改变和支座位移进行计算。
程序的解题规模只决定于计算机内存大小。
根据程序的功能规定，平面杆系静力分析程序划分为一个主模块和15个子模块。
此外还设计了若干个辅助模块，供有关模块调用，如矩阵运算模块和二维数组的动态内存分配和释放模块等。
模块层次图如图9．2所示，程序流程图则如图9．3所示。
图9．2 平面杆系静力分析模块图
图9.3 流程图
9．2．2 结构的基本属性描述
我们将大量使用结构类型（struct）变量对平面杆系结构的基本属性进行描述。
之所以选择结构类型是基于以下两点考虑：
首先，结构类可以把类型不同但在逻辑上联系较紧密的数据组合成一个有机整体，具有较强的表现力，可以提高程序的清晰度；其次，可以比较容易地把这种程序改写为面向对象的形式。
事实上，在C++中关键字“struct”和“class”都可被用来声明一个类，区别仅仅是，在没有显式说明访问属性时，前者的所有成员默认为是公有的，而后者则是私有的。
另外，使用结构变量可以减少函数的参数表中的参数个数。
结构的属性包括以下内容：
（1）材料属性。
我们只考虑各向同性的线弹性材料，它应该提供以下数据：
弹性模量，泊松比和线膨胀系数。
用结构变量Material表示如下：
structMaterial
｛
doubledE;//弹性模量
doubledMu;//泊松比
doubledAlpha;//线膨胀系数
｝；
（2）杆载面几何特性用结构变量Section描述如下：
structSection
｛
doubledA;//横截面面积
doubledIz;//截面性矩
doubledH;//截面高
｝；
（3）节点属性用结构变量Node描述如下：
structNode
｛
intiType;//节点自由度类型
doubledX,dY;//节点的X，Y坐标值
intiaDOFIndex[3];//节点的X向，Y向及转角向位移编号
｝；
（4）由于桁架和刚架单元有很多的属性是相同的，因此，对这两类单元统一用结构变量Element描述：
structElement
｛
intiType;//单元类型
intiaNode[2];//杆单元两端节点号
intiSecion;;//截面几何特性索引号
intiMaterial[a]//材料索引号
doubltdLength;//杆长
doubledSin,dCos;//杆在整体坐标中的方向余弦
doubledaEndInterForce[6];//杆端力数组
｝；
（5）载荷属性用结构变量Load描述如下：
structLoad
｛
intiType;//载荷类型号
intiDirect;//载荷作用方向
doubledValue;//载荷值
intiLoadedElem;//载荷作用的单元号
intiLoadedNode;//载荷作用的节点号
doubledPosition;//载荷作用位置或分布长度
doubledT0,dT1;;//杆上、下表面温度变化值
｝；
结构变量Load中的载荷类型如表9．2所示。
表9.2 载荷类型
类型号
（iType）
载荷类型
载荷值
（dValue）
载荷方向
（iDirect）
载荷位置或分布
长度（dPosition）
局部坐标系方向
0
q，M
0（x向），1（y向），
2（转角向）
——
1
q
a
2
q
——
a
3
q
——
a
4
q
——
a
5
q
——
a
6
q
——
a
7
q
——
a
8
t0
t1
to，t1
——
——
9
支座位移
0（x向），1（y向），
2（转角向）
——
（6）受约束节点用结构变量ConstrainedNode描述如下：
structConstrainedNode
｛
intiNode;//受约束节点号
intiaConstrainedDOF;//受约束自由度特征数，“0”表示未知
//独立位移，“-1”表示已知位移包括零位移，“10000+节点号”表示该位移分量所从属的主节点号。
｝；
9．2．3 总刚度矩阵的存储
结构的总刚度矩阵是一个大型的对称带状稀疏矩阵，通常采用压缩的方法将它存储在计算机中。
这种压缩存储法有利于提高计算机的存储效率和处理大型系数矩阵的能力，并可使大量的零元素不参与运算，从而节省机时。
下面介绍实际计算中经常采用的两种存储方法。
1．半带存储法
因为总刚度矩阵是一个带状阵，非零元素对称地分布在主对角线的两旁，而且已经证明，对于对称带状系数矩阵的线性方程组，无论采用高斯消去法还是三角分解法对它进行求解，都不会使这个带形区以外的零元素变成非零元素，也就是说在方程求解过程中，不会使系数矩阵的半带宽发生变化。
这里所谓半宽带是指系数矩阵的上三角（或下三角）部分，从主对角元素开始，到与主对角元位于同一行的最右边（与下三角阵对应的是最左边）非零元素为止的元素个数。
如图9．4所示。
所以，我们可以只存储上三角矩阵最大半带宽以内的元素。
图9.4半带存储
由单元集成法可知，最大半带宽决定于单元内两节点编号的最大差值，其一般计算公式为
Mb=（各单元两节点编号之差最大值+1）×nDOF（9．1）
式中Mb为最大半带宽，nDOF为单元节点的自由度个数。
设原来的总刚度矩律[K]为n×n阶对称方阵，最大半带宽为MB。
当采用半带存储时，可将[K]的带形区域内的上三角部分元素，存储在n×MB阶的矩阵[K]*中。
在矩阵[K]和[K]*中，元素的行列号有下面的对应关系：
（1）矩阵[K]的第r行元素在矩阵[K]*中仍为第r行，即元素的行号不变。
（2）矩阵[K]的主对角元素在矩阵[K]*中位于第1列。
因此，矩阵[K]的上三角部分第r行各元素在矩阵[K]*中的列号，都要比同一元素在[K]中的列号小（r-1）。
设元素在矩阵[K]中的行号为r，列号为s，又设同一元素在[K]*中的行号为i，列号为j，则上述对应关系可用公式表示为：
（9.2）
2．一维压缩存储法
考虑到总刚度矩阵[K]中各行的半带宽并不相等，而且有时由于结构的几何形状的原因，使总刚度矩阵某些行的半带宽特别大，而另外一些行的半带宽又特别小，在这种情况下，如仍然采用半带存储法就会把许多零元素也包含了进去，这对节省存储量是很不利的。
一维压缩存储法是将总刚度矩阵[K]的下三角矩阵的每一行，从第一个非零元素开始按行将元素排成一序列，存放于一维数组GK[L]中，其中L是一维数组元素的个数。
但是为了确定GK中的元素在[K]中的行列号，还需将[K]中各行对角线上的元素在一维数组中的序号（称作对角元地址），存放在另一辅助数组KD[n]中，其中n是总刚度矩阵[K]的阶数。
现举例说明这一存储法。
设有一系数矩阵
它的下三角矩阵各行半带内的元素总共有13个，它们在一维数组GK[B]中的存放次序是
；
而辅助数组KD[6]中存放的对角元地址是
；
将一结构离散化后，对节点和节点位移进行编号，就可以计算出总刚度矩阵[K]各行的半带宽，由它依次累加，就可以得出其对角线元素在一维数组中的序号（地址），即得到KD数组中各元素的值。
显然，形成了KD数组后，就确定了[K]中被存储的元素的分布情况，以及GK和[K]中元素位置的对应关系。
例如用KD数组可求出[K]中第i行的半带宽为
Bi=KD[i]-KD[i-1] （9.3）
而[K]中第i行第j列的元素在一维数组GK中的地址Pij借助于KD数组可用下式求得
Pij=KD[i]-（i-j） （9.
升级会员 