C++面向对象的有限元程序设计

【实例简介】
C++面向对象的有限元程序设计,空间8结点等参元的分析计算
第二章有限元分析计算的理论基础(空间8结点等参元) 有限元法进行的前提是满足四个假定:连续性假定、完全弹性假定、均匀与各 向同性假定、小变形和小位移假定。有限元法是建立在弹性力学的基础之上的。它 是以平衡微分方程(数学上)、变形协调方程(几何方程)、本构方程(物理方程) 作为基本的理论方程,同时又有圣维南原理、基于能量形式的虚位移原理作为解决 问题的手段。有限元法是依照弹性力学的基本解法进行求解的,不论有限元法的哪 种单元类型,在求解的过程和步骤上都是一样的,只不过求解的具体方法和细节 处理有所不同。本文就以空间8结点等参元的求解的过程为例来进行编程,别的单 元类型可以自己开发。而且,有了面向对象的于段,许多单元类型可以用继承的方 式加以实现。 2.1自然坐标系与位移模式 自然坐标系与直角坐标系,如图 (-1,-1,1)5 8(-1,1,1) X (1,-1,1)6 7(1,1,1 4(-1 s 3(1.1,1) 局部自然坐标系与整体直角坐标系之间的关系 (2-1) 其屮N是关于三个自然坐标系变量;.,t的矩阵,x,y1,21,x2,…,zs是六面体的八 个顶点(结点)在直角坐标系下的位置坐标,它们是已知的,xy,z与;s,t就有了对 应的关系。 N被称为形函数矩阵, N100N 0 N=0N100 NNI 小(22) 其中N1=N(r,s,t)=(+m)1+.+t11≤i≤8 (2-3) 整体直角坐标与局部自然坐标的关系又可以表达成: x ∑Ny ∑N 偎定结点间的位移变化是线性的,则位移模式可以表达为: (2-5) 其中1,V1,w1,u2,…,ws是八个顶点(结点)的位移。 2.2应变与位移间的关系 根据变形协调方程,有: Oulo z A E ou Ov (2-6) 6×1 6×24v Ow cu 其中[B=[B1B B (2-7) ON ON 0 aN B 1≤i<8 (2-8) 6×3 ON. ON 0 ay ax ON. aN aN aN 要求出[B,需要按照以下步骤推导求出: 由于 aN aN, ax an ay ON az ON Ox ar ay ar az or ax AN aN aN ay ON ON (2-9) OS oV N aN Ox aN, ay aN a aN at Ox ot ay az at 其中[为雅克比矩阵, OX ON ∑ ON ∑ OI OX aN aN x∑ ON ∑ ∑ (2-10) asas a ax av a ON N ∑ ∑|y∑ ON at at at X at aN r(+s:1+m1) N s、(+r)1 (1+r)1+)(2 Or 8 aN ON X CI aM 由式(2-9)可知 ON 小]1{∞,由此可求,即可得到B aN ON 2.3应力与应变间的关系 σ=[D]6,其中[D是弹性矩阵,它是由弹性模量E和泊松比A确定的 6×1 6×6 2.4单元刚度矩阵的求解 K(o]=SIB]IDIBI (2-12) 由于[B是关于自然坐标系;s,t的矩阵,上式的积分变量就是;s,t,所以 dv=dxdydz= det[J drdsdt (2-13) 这样= [B]LD[]det[ ] drdsdt (2-14) 利用数值积分中两点的高斯积分,可得 kx"]∑∑∑(Dle (2-15) 其中;;,t分别只能取得对应的两个值,即 R1=S1=T1=-0.5735,R2=S2=T2 2=0.5735 2.5整体刚度矩阵组装 整体刚度矩阵是由所有的单元刚度矩阵组装而成的,在组装过程中采取人家熟 知的“对号入座”的原则和方法。 对于空间8结点的单元刚度矩阵为24行24列的矩阵,可以分为8×8的子块, 每个子块是3行3列的矩阵;找到每个结点在本单元中的局部编号(从1至8之间) 和对应的整体编号(从1至n之间,n为整体的单元总数),整体刚度矩阵是一个3n 行3n列的矩阵,同样可以分为nxn的子块,每个子块同样是3行3列的矩阼。设 在单元刚度矩阵中有子块k,其中动为局部编号,其对应的整体编号为pq,则应 该填入整体刚度矩阵的k。子块中。在整体刚度矩阵中的同一个子块中可能填入多 个单元的单元刚度矩阵子块,这些子块应该先求和再填入整体刚度矩阵。将所有的 单元刚度矩阵如此“对号入座”,便组装成整体刚度矩阵。 整体刚度矩阵是个稀疏矩阵(0元素占绝人多数)、带状对称矩阵(非0元素只 对称分布于主对角线周围的狭窄的区域内)、奇异矩阵。 2.6通过边界条件求解 对于满足胡克定律的线弹性物体来说,有 K 3n×13n×3n3n×1 其中[为载荷向量,[为位移向量。 =[F 2 丁(2-17) [F为所有结点的各个方向的载荷的结合,它是其中所有元素是已知的; ]为所有结点的在各个方向的位移,它是部分已知、部分未知的,我们要求解 的就是那些未知量,已知的部分就是所谓的约束 由于整体刚度矩阵[K是奇异矩阵,就要通过这些约束将其修改为非奇异矩阵, 再来求解。修改方法如下 将[K]中主对角线上与已知约束同行的元素乘上一个大数,如103,同时将[F 同行的元素换成已知约束的位移值与该大数的成积。如,给定了w1,就将[K中的对 角元素k3乘以105,同时将F1×k3×103。 至此,便可以求解方程(2-16)了,通常用到稀疏矩阵和带状方程的求解器。 2.7求出单元的应力应变 通过以上的步骤已经求得了每个结点的位移,将每个单元中对应结点的位移代 入(2-6)屮,再通过数值高斯积分就可以得出应变,有了应变,根据2.3节屮的关 系算出应力。 根据本章对空间8结点等参元的分析计算的理论叙述,我们就可以利用面向对 象的方法进行编程了。其中还运用了大量数值分析的算法,在此不再详细说明。 第三章面向对象的程序设计 前面讲述了结构有限元法的基本原理,本章就是讲述简单的面向对象的有限元 程序设计的原理。面向对象的程序设计主要包括面向对象的分析(OOA)和面向对象 的设计(OOD) 3.1面向对象分析 面向对象的分析是重要的建模过程,它包括对象的识别、对象关系的确定、对 象的属性和方法的确定。 3.1.1对象的识别 对象的识别考虑的是在一个面向对象的程序系统中识别出对象。 我们很容易了解到在有限元分析计算中的模型是由若干个单元组成的,在本例 中是空间8结点等参元,每一个单元又是由若干个节点组成的。这就使我们至少抽 象出三个对象: Entity(实体类)、 Element(单元类)、Node(结点类) 对有限元的分析过程进行归纳整理,识别出了一个通用的数学对象: Matrix(矩阵类) 由于本文只是简单的对模型进行分析,有些构造大型有限元分析的类并没有加 以认定,如 ElemMng(单元管理类)、 Nodemng(结点管理类)、Load(载荷类)、 Constraint(约束类)、 LoadIng(载荷管理类)、 ConstrMng(约束管理类)等。未 加认定的对象在程序中均未实现。而单元的管理操作交给了实体类管理,结点的管 理操作也分给了实体类和单元类来管理 3.1.2对象关系的确定 对于有限元这个复杂的系统来说,类/对象之间的关系的确定最为重要 类之间的关系主要有两种:绊承关系和关联关系。 (1)继承关系 继承关系就是特殊与一般的关系,我们又形象地称为is-a的关系。继承关 系在本文这个简单的模型中最名显的体现是单元类,因为有限元分析中的单元 类型有很多,但它们都拥有一般单元的性质和操作方法,所以无论哪一种单元 类型都可以从通用的单元类中继承。本文的空间8结点等参元就是从单元类中 继承出来的。 以UML图表示的继承关系如下图所示 Element Hexahedron 图3-1空间8结点等参元是通用单元的继承 (2)关联关系 关联关系是指类之间通过某种方式发生联系。最主要的有聚合和委托关 系。聚合关系是指一个聚合类由其他若干个被聚合类组成,这些被聚合类作为 聚合的属性而存在,我们通常称为has-a关系。本例中的实体是由若干个单元 和结点组成的,而单元又是由8个结点组成的,这构成了最为明显的聚合关系。 以UML图表示的聚合关系如下图: Entity Element N ode 图3-2实体包含单元和结点,单元也包含结点 委托关系是指两个类木无任何联系,一个类通过一定的方式委托另一个类 执行部分的工作,最后完成自己的工作。有限元分析计算的过程中用到大量的 数学数值运算和关于矩阵、向量的运算,它需要委托这些数学对象进行运算, 帮助其解决问题。在此例中,由于实体的部分计算操作和单元的计算操作都是 由关于矩阵的运算操作,这些类就当委托类的角色,矩阵就是被委托的类。 以UML图表示的委托关系如下图 Element Matrix Entity 图3-3单元类和实体类委托矩阵类进行运算 本例中的矩阵类实现也委托了一个链表类。 对于本例中的所有类,总体的关系如下图: Entit Node Matrix Element Linklist H exahearor Listnode 图3-4总体对象间关系 3.1.3对象的属性和方法的确定 在确定了对象的关系以后就可以确定对象的属性和方法了。 下面的表只列举重要的类的属性和方法: 【实例截图】
【核心代码】