面向对象的有限元程序设计.pdf
《面向对象的有限元程序设计.pdf》由会员分享,可在线阅读,更多相关《面向对象的有限元程序设计.pdf(60页珍藏版)》请在淘文阁 - 分享文档赚钱的网站上搜索。
1、 1面向对象的有限元程序设计面向对象的有限元程序设计 空间空间 8 结点等参元的分析计算结点等参元的分析计算 第一章第一章 概述概述 有限元法是用来分析各种结构问题的强有力工具。在经典的固体力学理论工作中,尽管人们已经进行了几百年的努力,但所能解决的实际问题为数不多。然而,有限元法却能成功地解决各种各样的固体力学问题,如杆系、板与壳、复杂的三维物体和大变形问题等。不论结构的几何形状和边界条件多么复杂,不论材料性质和外加载荷如何多变,使用有限元法均可获得满意的答案。有限元解决实际问题的能力远远超过了经典的方法,并且已经取得了很大的成就,因而受到普遍重视。无论大型飞机、大型舰船还是高层建筑、水利大
2、坝,有限元法均可方便地进行结构的分析。现在,掌握有限元法的原理和应用,对于一个从事结构分析与设计的工程师来说,已经是必不可少的了。有限元法的实施,离不开程序设计。如果一个结构工程师,不仅懂得有限元法的原理,而且还会应用有限元程序设计方法,无疑会在工作中如虎添翼。有限元的程序设计内容包括三部分:前处理、分析计算、后处理。其中,前处理主要负责读入数据、生成模型、网格划分,为有限元的计算做好准备;分析计算主要是进行有限元矩阵的计算、组装、求解;后处理主要对计算的结果进行各个方式、各个角度的输出。本文主要进行的是分析计算部分的程序设计,不考虑前处理和后处理的程序设计,数据的输入和结果的输出都以固定格式
3、的文件进行处理。分析计算阶段的主要内容有:(1)计算单元的刚度矩阵;(2)组装总体刚度矩阵;(3)根据已知的位移边界条件,进行约束处理,消除总体刚度矩阵的奇异性;(4)求解整体刚度矩阵,得出节点的位移向量;(5)根据节点的位移向量求出每个单元的应力和应变状态。面向对象的程序设计方法是一种新兴的程序设计方法,或者说是一种新的程序设计思想,这种基本思想是运用对象、类、继承、封装、聚合、消息传递、多态性 2等概念来构造系统进行程序设计。该方法强调直接以问题域中的事物为中心来思考问题、认识问题,并根据这些事物的本质特征,把它们表示为系统中的对象,作为系统构成的基本单位。它对描述事物的属性、特点及事物之
4、间的关系具有一种先进的软件开发方法,较传统的程序设计思想更接近于人的思维。它把待求解问题变简单,而且是一整套关于如何看待软件系统与现实关系,以什么观点来研究问题并进行求解,以及进行系统构造的软件方法学。与以往的结构化程序设计相比,面向对象的程序设计更结构化、更模块化、更抽象。面向对象方法的主要特点是:(1)从问题域中客观存在的事物出发来构造软件系统,用对象作为这些事物的抽象表示,并以此作为系统的基本构成单位。(2)事物的静态特征由对象的属性表示,动态特征由对象的方法表示。(3)对象的属性与方法结合成一个整体。(4)对事物进行分类。(5)运用抽象的原则。(6)复杂对象可以由简单对象聚合。(7)对
5、象间通过消息进行通信。(8)通过关联表达对象之间的静态关系。面向对象方法的优点有:(1)维护简单。面向对象支持、鼓励软件工程实践中的信息隐藏、数据抽象和封装,在一个对象内的修改被局部隔离。(2)可扩充性。可以根据需要在各个类中增加内容,而不影响其他类。(3)代码重用。面向对象开发鼓励重用,不仅软件的重用,还包括分析、设计的模型重用。所以面向对象方法是当今计算机领域的主流技术,对促进计算机科学技术发展具有十分重要的意义。第二章第二章 有限元分析计算的理论基础有限元分析计算的理论基础(空间空间 8 结点等参元结点等参元)有限元法进行的前提是满足四个假定:连续性假定、完全弹性假定、均匀与各向同性假定
6、、小变形和小位移假定。有限元法是建立在弹性力学的基础之上的。它是以平衡微分方程(数学上)、变形协调方程(几何方程)、本构方程(物理方程)作为基本的理论方程,同时又有圣维南原理、基于能量形式的虚位移原理作为解决问题的手段。有限元法是依照弹性力学的基本解法进行求解的,不论有限元法的哪一种单元类型,在求解的过程和步骤上都是一样的,只不过求解的具体方法和细节处理有所不同。本文就以空间 8 结点等参元的求解的过程为例来进行编程,别的单元类型可以自己开发。而且,有了面向对象的手段,许多单元类型可以用继承的方式加以实现。21 自然坐标系与位移模式自然坐标系与位移模式 自然坐标系与直角坐标系,如图:局部自然坐
7、标系与整体直角坐标系之间的关系:=82111zxzyxNzyx?(2-1)其中 N 是关于三个自然坐标系变量 r,s,t 的矩阵,是六面体的八个顶点(结点)在直角坐标系下的位置坐标,它们是已知的,x,y,z 与 r,s,t 就有了对应的关系。82111,zxzyx?N 被称为形函数矩阵,3INININNNNNNN821811210000000000?=(2-2)其中()()()()81,11181,+=ittssrrtsrNNiiiii (2-3)整体直角坐标与局部自然坐标的关系又可以表达成:=818181iiiiiiiiizNyNxNzyx (2-4)假定结点间的位移变化是线性的,则位移模式
8、可以表达为:=82111wuwvuNwvu?(2-5)其中是八个顶点(结点)的位移。82111,wuwvu?22 应变与位移间的关系应变与位移间的关系 根据变形协调方程,有:=+=8211124616wvwvuBzuxwywzvxvyuzwyvxuzxyzxyzzyyxx?(2-6)其中 (2-7)821246BBBB?=46 3000000,18000iiiiiiiiiiNxNyNzBiNNyxNNzyNNzx=(2-8)要求出B,需要按照以下步骤推导求出:由于 =+=zNyNxNJtzzNtyyNtxxNszzNsyyNsxxNrzzNryyNrxxNtNsNrNiiiiiiiiiiiii
9、ii (2-9)其中J为雅克比矩阵,=81818181818181818133iiiiiiiiiiiiiiiiiiiiiiiiiiiztNytNxtNzsNysNxsNzrNyrNxrNtztytxszsysxrzryrxJ (2-10)5()(iiiittssr)rN+=1181,()(iiiittrrssN+=)1181,()(iiiissrrttN+=)1181(2-11)由式(2-9)可知=tNsNrNJzNyNxNiiiiii1,由此 iB可求,即可得到B。23 应力与应变间的关系应力与应变间的关系 166616=D,其中D是弹性矩阵,它是由弹性模量E和泊松比确定的。24 单元刚度矩
10、阵的求解单元刚度矩阵的求解 =)()(evTedVBDBK (2-12)由于B 是关于自然坐标系r,s,t的矩阵,上式的积分变量就是r,s,t,所以 drdsdtJdxdydzdVdet=(2-13)这样 =111111)(detdrdsdtJBDBKTe (2-14)利用数值积分中两点的高斯积分,可得 ()()=212121,)(|detRRrSSsTTttsrTeJBDBK (2-15)其中r,s,t分别只能取得对应的两个值,即 5735.0,5735.0222111=TSRTSR 25 整体刚度矩阵组装整体刚度矩阵组装 整体刚度矩阵是由所有的单元刚度矩阵组装而成的,在组装过程中采取大家熟
11、知的“对号入座”的原则和方法。对于空间 8 结点的单元刚度矩阵为 24 行 24 列的矩阵,可以分为 88 的子块,每个子块是 3 行 3 列的矩阵;找到每个结点在本单元中的局部编号(从 1 至 8 之间)和对应的整体编号(从 1 至 n 之间,n 为整体的单元总数),整体刚度矩阵是一个 3n行 3n 列的矩阵,同样可以分为 n n 的子块,每个子块同样是 3 行 3 列的矩阵。设在单元刚度矩阵中有子块,其中i,j为局部编号,其对应的整体编号为p,q,则应该填入整体刚度矩阵的子块中。在整体刚度矩阵中的同一个子块中可能填入多个单元的单元刚度矩阵子块,这些子块应该先求和再填入整体刚度矩阵。将所有的
12、单元刚度矩阵如此“对号入座”,便组装成整体刚度矩阵。)(,ejikqpk,整体刚度矩阵是个稀疏矩阵(0 元素占绝大多数)、带状对称矩阵(非 0 元素只对称分布于主对角线周围的狭窄的区域内)、奇异矩阵。26 通过边界条件求解通过边界条件求解 对于满足胡克定律的线弹性物体来说,有 6 133313=nnnnKF (2-16)其中F为载荷向量,为位移向量。TnTznxzyxwuwvuFFFFFF?21112111,=(2-17)F为所有结点的各个方向的载荷的结合,它是其中所有元素是已知的;为所有结点的在各个方向的位移,它是部分已知、部分未知的,我们要求解的就是那些未知量,已知的部分就是所谓的约束。由
13、于整体刚度矩阵K是奇异矩阵,就要通过这些约束将其修改为非奇异矩阵,再来求解。修改方法如下:将K中主对角线上与已知约束同行的元素乘上一个大数,如 1015,同时将F同行的元素换成已知约束的位移值与该大数的成积。如,给定了,就将K中的对角元素乘以 101w3,3k15,同时将。153,3110kFz至此,便可以求解方程(2-16)了,通常用到稀疏矩阵和带状方程的求解器。27 求出单元的应力应变求出单元的应力应变 通过以上的步骤已经求得了每个结点的位移,将每个单元中对应结点的位移代入(2-6)中,再通过数值高斯积分就可以得出应变,有了应变,根据 2.3 节中的关系算出应力。根据本章对空间 8 结点等
14、参元的分析计算的理论叙述,我们就可以利用面向对象的方法进行编程了。其中还运用了大量数值分析的算法,在此不再详细说明。7 8第三章 面向对象的程序设计第三章 面向对象的程序设计 前面讲述了结构有限元法的基本原理,本章就是讲述简单的面向对象的有限元程序设计的原理。面向对象的程序设计主要包括面向对象的分析(OOA)和面向对象的设计(OOD)。31 面向对象分析面向对象分析 面向对象的分析是重要的建模过程,它包括对象的识别、对象关系的确定、对象的属性和方法的确定。311 对象的识别对象的识别 对象的识别考虑的是在一个面向对象的程序系统中识别出对象。我们很容易了解到在有限元分析计算中的模型是由若干个单元
15、组成的,在本例中是空间 8 结点等参元,每一个单元又是由若干个节点组成的。这就使我们至少抽象出三个对象:Entity(实体类)、Element(单元类)、Node(结点类)对有限元的分析过程进行归纳整理,识别出了一个通用的数学对象:Matrix(矩阵类)由于本文只是简单的对模型进行分析,有些构造大型有限元分析的类并没有加以认定,如 ElemMng(单元管理类)、NodeMng(结点管理类)、Load(载荷类)、Constraint(约束类)、LoadMng(载荷管理类)、ConstrMng(约束管理类)等。未加认定的对象在程序中均未实现。而单元的管理操作交给了实体类管理,结点的管理操作也分给了
16、实体类和单元类来管理。312 对象关系的确定对象关系的确定 对于有限元这个复杂的系统来说,类/对象之间的关系的确定最为重要。类之间的关系主要有两种:继承关系和关联关系。(1)继承关系 继承关系就是特殊与一般的关系,我们又形象地称为 is-a 的关系。继承关系在本文这个简单的模型中最名显的体现是单元类,因为有限元分析中的单元类型有很多,但它们都拥有一般单元的性质和操作方法,所以无论哪一种单元类型都可以从通用的单元类中继承。本文的空间 8 结点等参元就是从单元类中继承出来的。以 UML 图表示的继承关系如下图所示:图图 3-1 空间空间 8 结点等参元是通用单元的继承结点等参元是通用单元的继承(2
17、)关联关系 关联关系是指类之间通过某种方式发生联系。最主要的有聚合和委托关系。聚合关系是指一个聚合类由其他若干个被聚合类组成,这些被聚合类作为聚合的属性而存在,我们通常称为 has-a 关系。本例中的实体是由若干个单元和结点组成的,而单元又是由 8 个结点组成的,这构成了最为明显的聚合关系。以 UML 图表示的聚合关系如下图:图图 3-2 实体包含单元和结点,单元也包含结点实体包含单元和结点,单元也包含结点 委托关系是指两个类本无任何联系,一个类通过一定的方式委托另一个类执行部分的工作,最后完成自己的工作。有限元分析计算的过程中用到大量的数学数值运算和关于矩阵、向量的运算,它需要委托这些数学对
18、象进行运算,帮助其解决问题。在此例中,由于实体的部分计算操作和单元的计算操作都是由关于矩阵的运算操作,这些类就当委托类的角色,矩阵就是被委托的类。以 UML 图表示的委托关系如下图:9 图图 3-3 单元类和实体类委托矩阵类进行运算单元类和实体类委托矩阵类进行运算 本例中的矩阵类实现也委托了一个链表类。对于本例中的所有类,总体的关系如下图:图图 3-4 总体对象间关系总体对象间关系 313 对象的属性和方法的确定对象的属性和方法的确定 在确定了对象的关系以后就可以确定对象的属性和方法了。下面的表只列举重要的类的属性和方法:10 11Matrix 类 分类分类 描述描述 说明说明 unsigne
19、d int Row 矩阵行数 unsigned int Col 矩阵列数 属性 LinkList Link 矩阵中元素 void setDimension(const unsigned int newLength)设置向量长度 void setDimension(const unsigned int newRow,const unsigned int newCol)设置矩阵大小 Matrix&operator=(const Matrix&rhs)矩阵赋值 unsigned int getRow()const 取得行数 unsigned int getCol()const 取得列数 T&oper
20、ator()(const unsigned int index)得到向量元素 T&operator()(const unsigned int x,const unsigned int y)得到矩阵元素 Matrix&changerow(const unsigned int x,const unsigned int y)交换两行 Matrix&changecol(const unsigned int x,const unsigned int y)交换两列 Matrix trans()矩阵转置 T det()求方阵的行列式 Matrix inv()矩阵求逆 Matrix copny()求方阵的伴随
21、矩阵 Matrix operator+(Matrix&rhs)两矩阵相加 Matrix operator-(Matrix&rhs)两矩阵相减 Matrix operator*(Matrix&rhs)两矩阵相乘 Matrix operator*(T&rhs)矩阵乘以变量 Matrix operator*(const T&rhs)矩阵乘以常数 方法 void Display()矩阵的显示 12Node 类 分类分类 描述描述 说明说明 unsigned int Id 结点编号 NodeType Kind 结点类型 double x,y,z 结点的空间位置 double fx,fy,fz 结点所受载
22、荷 double dispx,dispy,dispz 结点位移 NodeProperty Prop 结点属性 属性 LinkList Adj 结点的邻接结点 unsigned int getId()const 取得结点编号 void setId(const unsigned int newId)设置结点编号 NodeType getKind()const 取得结点类型 void setKind(NodeType newKind)设置结点类型 void setx(double newx)void sety(double newy)void setz(double newz)设置结点的空间位置 v
23、oid setfx(double newfx)void setfy(double newfy)void setfz(double newfz)设置结点的载荷 void setdispx(double newdispx)void setdispy(double newdispy)void setdispz(double newdispz)设置结点的位移 double getx()const double gety()const double getz()const 取得结点的空间位置 double getfx()const double getfy()const double getfz()co
24、nst 取得结点所受载荷 double getdispx()const double getdispy()const double getdispz()const 取得结点的位移 NodeProperty getProp()const 取得结点属性 方法 void setProp(NodeProperty newProp)设置结点属性 13void addAdjNode(unsigned int NodeId)添加一个邻接结点 unsigned int getAdjNodeId(const unsigned long int index)取得一个邻接结点的编号 Element 类 分类分类 描
25、述描述 说明说明 unsigned int Id 单元编号 double E,mu 单元弹性模量与柏松比 ElementType Kind 单元类型 unsigned int NodeNum 单元包含的结点数量 Matrix K 单元刚度矩阵 Matrix D 单元弹性矩阵 Matrix Stress 单元应力 Matrix Strain 单元应变 Node*NodeList 单元所包含的结点 属性 Matrix B 单元 B 矩阵 unsigned int getId()const 取得单元编号 void setId(const unsigned int newId)设置单元编号 doubl
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 面向 对象 有限元 程序设计
限制150内