面向对象有限元程序设计及其VC_与Matlab混合编程实现.pdf

《面向对象有限元程序设计及其VC_与Matlab混合编程实现.pdf》由会员分享，可在线阅读，更多相关《面向对象有限元程序设计及其VC_与Matlab混合编程实现.pdf（5页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

1、第 26 卷第 12 期 岩 土 力 学 Vol.26 No.12 2005 年 12 月 Rock and Soil Mechanics Dec.2005 收稿日期：2004-11-02 修改稿收到日期：2005-01-18 作者简介：史贵才，男，1979 年生，博士，主要从事岩土工程数值模拟方向的研究工作。E-mail: 文章编号：文章编号：10007598(2005)12200505 面向对象有限元程序设计及其面向对象有限元程序设计及其 VC+与与 Matlab 混合编程实现混合编程实现 史贵才史贵才1,2，葛修润，葛修润1(1.中国科学院 岩土力学重点实验室，武汉 430071；2.常

2、州工学院 土木建筑工程学院，江苏 常州 213002)摘摘 要：要：应用面向对象方法来研究有限元，是对有限元新方法的有益尝试和创新性发展。通过对比面向过程和面向对象的程序设计方法，讨论了面向对象方法与有限元程序设计相结合的优点，并简要回顾了国内外面向对象的程序设计方法的研究进展。应用面向对象的程序分析方法，建立了三维脆塑性有限元分析类库。采用 VC+和 Matlab 混合编程的手段，设计了基于Windows98/2000/NT 操作平台的面向对象的三维脆塑性有限元分析软件，成功地分析了国内某大型水电站地下硐室群围岩稳定性，验证了该面向对象有限元分析程序的有效性和实用性。关关 键键 词：词：面向

3、对象；脆塑性；有限元；程序设计 中图分类号：中图分类号：O 241.82；TB 115 文献标志码：文献标志码：A Object-oriented finite element method and programming by combining VC+with Matlab SHI Gui-cai1,2,GE Xiu-run1 (1.Key Laboratory of Rock and soil Mechanics,Institute of Rock and Soil Mechanics,Chinese Academy of Sciences,Wuhan 430071,China;2.Sc

4、hool of Civil Engineering and Architecture,Changzhou Institute of Technology,Changzhou 213002,China)Abstract:Applying object-oriented programming to researching finite element method is a beneficial effort and creativity development.By comparing the conventional procedural programming and object-ori

5、ented programming,the advantages of combining the FEM program and object-oriented programming method are discussed;and the progress of object-oriented programming at home and abroad is reviewed.Based on a description of the major characters of object-oriented programming,a 3D brittle-plastic finite

6、element analysis class library is constructed.By combining VC+with Matlab,a geotechnical FEM software based on Windows98/2000/NT is designed to deal with 3D brittle-plastic problems.The successful application to the stability analysis of an underground excavation of a certain hydropower project prov

7、es the validity and practicability of this object-oriented program,which can be of beneficial reference to analogous projects.Key words:object-oriented;brittle-plasticity;finite element method;programming 1 引 言 有限单元法是计算力学在过去 50 余年中的最大进展，在固体力学、热传导、计算流体动力学、渗流、电磁场等领域得到广泛应用，已成为数学物理中几乎无所不能的数值计算方法1。传统的有限元

8、程序设计一般采用结构化的程序设计方法，其特点是代码和数据的独立性，即数据结构和其操作过程彼此分离。这就造成当系统需求改变时，这些程序将不能被重用和移植，甚至十分简单的改动会产生整个程序的代码修改，特别是数据结构的变化将贯穿整个程序代码的更动2。其程序的扩展能力有限，代码的重利用率低，调试复杂。研究表明3,4，这种以算法为核心，过程和数据分离的结构化程序设计方法，已越来越不适应有限元软件的发展要求。面向对象程序设计方法集数据抽象、抽象数据类型和类型继承为一体，使软件设计中人们普遍遵循的模块化、信息隐藏、抽象和代码共享等思想在面向对象机制下得以充分实现，从而成为一种强有力的软件设计模式；另外，面向

9、对象的程序设计，岩 土 力 学 2005 年 由于程序具有封装性、继承性和多态性等优点，使得程序设计概念清楚，调试容易，代码的重利用率高，已成为现代程序设计的主要方法之一。另一方面，作者认为，尽管面向对象技术在有限元分析中有巨大的优势和广阔的应用前景，但是其并没有从根本上改变有限元方法的核心，并不能突破有限元本身发展的瓶颈。将 OOP 应用于有限元分析源于 1990 年，Forde首次将 OOP 应用于有限元程序设计，提示了 OOP将有限元知识结构在程序空间更清晰地加以表达的可能4。受 Forde 工作的启示，在国际上，De Oliveira Sousa 和 Ingraffea 用 OOFEM

10、 求解了一个岩石裂隙和流体的有限元问题；Udo Meissner等人用OOFEM对一个三维岩土工程问题进行了分析5；Wemer 等人研究了在隧道工程中如何用面向对象的模型来设计与分析问题6。在国内，公开报导的有：孔详安、张向、魏泳涛、曹本清、马其罕、李会平、项阳等。面向对象的有限元程序设计方法已经成为一个研究热点，而且取得了相当进展。近年来，面向对象程序设计方法在线性静态有限元分析领域的有效性已得到了论证，而在非线性领域的探索还很少见报道。本文将给出自己编制的三维弹脆塑性有限元的程序结构及其实现。2 面向对象的基本概念 2.1 对象对象(Object)对象是客观世界中的一个实体，是其自身所具有

11、的状态特征及对这些状态施加的操作结合在一起所构成的独立实体。它把客观世界的实体和计算机系统运行实体有机的结合在一起。要设计一个面向对象的有限元程序，对象的确定和划分非常重要，对象的确定和划分是否得当直接影响所编制的程序的质量。目前，关于如何合理划分有限元分析中的对象，共同性的结论还不多。2.2 数 据 抽 象 和 封 装（数 据 抽 象 和 封 装（Data Abstraction and Encapsulation）数据抽象是指从较特殊的类或对象中抽出一般属性以建立一个超类的过程7。封装又称数据隐藏，是指将一个数据和与这个数据有关的操作集合放在一起，形成一个能动的实体对象。这一特性大大地降低

12、了模块间的耦合性，从而提高了程序的可靠性，尽可能地排除了对数据进行任意访问造成的隐患。2.3 继承性（继承性（Inheritance）继承性就是指一个类可以继承其父类的所有数据和成员函数，同时又可以定义自己的数据和成员函数。这样做一方面可以减少代码冗余，另一方面可以通过协调性来减少相互之间的接口和界面。图1 给出了面向对象有限元程序中单元类的简单的继承关系。图图 1 单元继承关系图单元继承关系图 Fig.1 The heritance diagram of elements 2.4 多态性（多态性（Polymorphism）多态性是指同一消息可以根据发送消息的对象的不同采用不同的行为方式。由于

13、这种特性，不同的对象接收到用户统一发送的消息就可完成不同的工作。C+语言支持两种多态，即函数的重载和虚函数。例如：不同类型的单元有不同的计算单元刚度矩阵的函数，可以通过虚函数的重载来实现。3 面向对象的有限元程序分析 有限元方法作为求解数学表述的连续体的一般离散化方法，是一种求解连续体问题的近似方法。其一般方法是：（1）把连续体分成有限个部分，每个部分的性态由有限个参数所规定；（2）求解作为单元的集合体的整个系统，此时其单元所遵循的是标准的离散体问题。根据有限元分析的具体过程，可大致划分为 7个步骤，即：（1）离散化，（2）插值，（3）单元刚度计算，（4）单元集成，（5）引入边界条件，（6）求

14、解有限元控制方程（7）进行辅助计算。按照有限元的分析方法，有限元分析的主要数据包括：（1）描述有限元分析的整体数据，如单元总数、节点总数、问题的维数、材料种类数，问题类型指示数，屈服准则指示数，收敛容差等；（2）单元数据，包括单元类型、单元材料号、单元包含的节点号，刚度矩阵和质量矩阵，高斯积分点维数等；（3）节点数据，包括节点坐标、节点自由度、节点力、节点位移、节点约束等；（4）高斯积分点数据，包括高斯积分点坐标、高斯积分点应力，高斯积分点应变等。2006 第 12 期 史贵才等：面向对象有限元程序设计及其 VC+与 Matlab 混合编程实现 4 面向对象有限元程序设计实例化 按照上面的面向

15、对象有限元程序分析，一个完整的面向对象有限元的分析程序应该设计的类主要有：（1）有限元整体类和相关的方法；（2）单元数据类和相关的方法；（3）节点数据和相关的方法；（4）高斯积分点数据相关方法；（5）材料数据类和方法等。下面给出了作者最近编制的三维弹脆塑性有限元的程序结构及其实现（限于篇幅，仅仅给出了部分代码）。4.1 有限元整体类和实现方法有限元整体类和实现方法 有限元分析的整体数据，如单元总数、节点总数、问题的维数、材料种类数、问题类型指示数、屈服准则指示数、收敛容差等。其具体定义如下：class CFEM:public:/控制数据 BOOL m_bTwoDimesion;/问题的维数 i

16、nt ProbTyp;/问题类型指示数 int Crition;/屈服准则指示数 int Tolerance;/收敛容差 int NumOfNodes;/节点总数 int NumOfElements;/单元总数 int NumOfMaterials;/材料总数 public:/实现方法 void FormGeneralRigidityMatrix();/形成整体刚度矩阵 void FormGeneralForceMatrix();/形成整体载荷列阵 void DealWithCondition();/处理边界约束条件 void Solution();/根据问题类型指示数选择/求解器求解控制方程

17、 void Residual();/计算残余节点力 void JudgeConvergence();/判断收敛 4.2 单元类和实现方法单元类和实现方法 在岩土工程有限元三维分析中，普遍使用的单元类型主要四面体、三棱柱、六面体等，另外还有一些特殊的单元，如锚杆单元、节理单元、无限元等。如图 1 所示，对不同的单元类型进行数据抽象，形成单元基类，包括单元类型、单元材料号、单元包含的节点号，单元刚度矩阵，高斯积分点数等。具体定义如下：class CElement :public CObject public:/单元相关数据 static CObArray*m_pNodes;/节点指针 static

18、 CObArray*m_Arr_pMaterials;/材料指针 CGaussPoint*m_gaussp;/高斯点指针 int m_iElementIDNumber;/单元序号 int m_iMaterialNumber;/单元材料号 int m_nNodeAmount;/单元节点数 int*m_pNodesNumber;/包含的节点号 int m_nGaussAmount;/高斯积分点数 int m_nDamageType;/单元的破坏类型 static mwArray*m_mwA_pGenRigMatrix;/总刚指针 mwArray m_mwA_RigidityMatrix;/单刚 p

19、ublic:/实现方法 virtual void ComputeRigidityMatrix();/虚函数，计算单元刚度矩阵 void AddElementtoGeneralRigidityMatrix();/虚函数，组装总刚矩阵 virtual void Computer_Stress();/虚函数，计算单元应力 virtual double Ritem(int,int,int);/虚函数，计算右端项 virtual void JudgeElementDamage();/虚函数，判断单元破坏与否 不同的单元类型将继承基类的所有数据和成员函数，同时又可以定义自己的数据和成员函数。例如，不同单元

20、的单元刚度矩阵计算方法不同，可以通过改写虚函数的方法实现自己的方法；但是，不同单元的单刚向总刚的组装方法是可以相同的，可以通过继承基类的方法来实现。4.3 点类和实现方法点类和实现方法 节点类的数据，主要包括节点坐标、节点自由度、节点力、节点位移、节点约束等。节点类描述如下：class CNode :public CObject public:int m_iNodeIDNumber;/节点编号 double m_pCoor3;/节点坐标 double m_dDisplacement3;/节点位移 double m_dForce3;/节点力 mwArray m_dStrain;/节点应变 mwA

21、rray m_dStress;/节点应力 int m_iDOFConstraints;/节点约束 节点约束状态用一个整型数表示，若 x 向自由度被约束其值为 1，y 向自由度被约束其值为 2，x2007 岩 土 力 学 2005 年 和 y 向自由度均被约束其值为 3，z 向自由度被约束其值为 4，依此类推。4.4 高斯积分点类和实现方法高斯积分点类和实现方法 高斯积分点数据，包括高斯积分点坐标、高斯积分点应力、高斯积分点应变等。高斯点类描述如下：class CGaussPoint public:mwArray m_mwA_gaussB;/高斯点的 B 矩阵 int m_nStatus;/高斯

22、点的弹塑性状态 int m_nDamageHistory;/高斯点的破坏历史 mwArray m_dGaussStrain;/高斯点应变 mwArray m_dGaussStress;/高斯点应力 int m_nGaussStress;/高斯点应力分量数 public:virtual mwArray ComputerJacobi();/计算高斯点的 Jacobi 矩阵 virtual mwArray ComputerBMaxtrix();/计算高斯点的 B 矩阵 virtual void ComputerGaussStress();/计算高斯点的应力 virtual double Comput

23、erInvar();/计算等效应力和单向屈服应力 virtual void ComputerFlowVector();/计算流动矢量 virtual mwArray ComputeResidualNodeForce();/计算残余节点力 对于一般的弹塑性问题，高斯点的弹塑性状态可以用一个 BOOL 型的变量标识。但是，对如图 2所示的岩土工程中常见的脆塑性问题分析，两个状态还不够，故本文用一个整型数表示，其值为 0，表示该高斯点处于弹性状态，且从未发生过应力跌落；1 表示该高斯点发生过应力跌落，但目前处于弹性状态；2 表示该高斯点处于塑性状态。图图 2 脆塑性材料示意脆塑性材料示意 Fig.2

24、 Description of brittle-plastic material 4.5 材料类和实现方法材料类和实现方法 材料类数据主要包括材料号、材料的力学参数（如弹性模量、泊松比、密度等），实现方法主要是计算弹性矩阵等。本文限于篇幅，不再罗列。另外，有限元程序中大量使用整型数组，经常遇到矩阵代数运算且总体刚度矩阵具有对称性、稀疏性，呈带形分布等特点。查阅所有公开发表的关于面向对象有限元程序设计文献，几乎无一例外的选择了作者自己编制矩阵类，但是花费了大量的时间和精力设计出的矩阵类，总体而言，并不完美。MATLAB 提供了 C/C+数学库，其中的 C+数学库功能很强，使用它可以节省大量的程序

25、编制工作量，而且可以使程序十分简洁。在 VC+中包含Matlab.hpp 后形成的单元刚度矩阵的程序段如下：CRockMaterial*pmater;/单元的材料指针 pmater=(CRockMaterial*)m_Arr_pMaterials-GetAt(m_iMaterialNumber);m_mwA_RigidityMatrix=zeros(24,24);/为单元刚度矩阵分配空间 mwArray j,J,BTD,BTDB,BTDBJ,Volume;/定义变量为 mwArray 类型 for(intiGuass=0;iGuassElasticMatrixD()；/计算BTD BTDB=m

26、times(BTD,m_gausspiGuass.m_mwA_gaussB);/计算BTDB Volume=det(J)*m_gausspiGuass.m_gweight;BTDBJ=times(BTDB,Volume.ExtractScalar(1);m_mwA_RigidityMatrix=m_mwA_RigidityMatrix+BTDBJ;本文在 VC+中嵌入 MATLAB 进行混合编程，以替代矩阵类。实践证明，虽然混合编程编制的程序在速度上仍然比纯粹的 VC+程序慢，但是由此换来的高效的开发效率和可靠性是值得肯定的。MPB2008 第 12 期 史贵才等：面向对象有限元程序设计及其

27、VC+与 Matlab 混合编程实现 作者使用本程序成功地对国内某大型水电站地下硐室群进行了围岩稳定性分析，取得了良好的效果，验证了在 VC+中嵌入 MATLAB 进行混合编程所生成的有限元分析程序的有效性、实用性和先进性。图 3 给出了使用本程序（EBPFEM）进行该地下硐室群的计算的用户界面。图图 3 EBPFEM 的用户界面的用户界面 Fig.3 Interface of EBPFEM 5 结 论(1)随着有限元技术处理问题的复杂性和规模的增加，传统的面向过程编程在大型有限元分析中暴露出了越来越多的缺点和问题，而面向对象编程技术可有效解决这些问题。(2)面向对象编程方法的多态性和重载机制

28、使得整个问题域的信息响应变得越来越简单。同时直接的指针操作和内存使用的动态分配，使得大型有限元分析中的速度和效率都有很大的提高。(3)采用在 VC+中嵌入 MATLAB 进行混合编程易于进行可视化开发工作，其高效的开发效率和可靠性是值得肯定的。(4)如前所述，面向对象技术应用于有限元有诸多优越性，但是，其作用也不应该被过分夸大。面向对象技术仅仅是一种程序设计理念和方法的改进，它并不能突破有限元技术本身的瓶颈而使有限元变得“无所不能”。参参 考考 文文 献献 1 魏泳涛,于建华,陈君楷.面向对象有限元程序设计基本数据类J.四川大学学报（工科版）,2001,33(2):1721.WEI Yong-

29、tao,YU Jian-hua,Chen Jun-kai.Object-oriented approach to finite element programming:basic data classesJ.Journal of Sichuan University(Engineering Science Edition),2001,33(2):1721.2 马永其,陈罕,李斯特.面向对象有限元程序的研究J.计算机工程与应用,2001,9:120122.MA Yong-qi,CHEN Han,LI Si-te.Object-oriented programming for finite ele

30、ment analysisJ.Computer Engineering and Application,2001,9:120122.3 Sanal Z.Finite element programming and CJ.Computer&Structure,1994,51(6):671686.4 Forde B W R,Foschi R O,Stiemer S F.Object-oriented finite element analysisJ.Computer&Structure,1990,34(3):355374.5 Udo Meissner,Joaquin Diaz,et al.Obje

31、ct-oriented analysis of three dimensional geotechnical engineering systemsA.Computing in Civil and Building EngineeringC.Rootterdam:Balkema A A,1995.6165.6 Wemer H,Mackert M.Object oriented models and tools in tunnel design and analysisA.Computing in Civil and Building EngineeringC.Rootterdam:Balkema A A,1995.107112.7 项阳,平扬,葛修润.面向对象有限元方法在岩土工程中的应用J.岩土力学,2000,21(4):346349.XIANG Yang,PING Yang,GE Xiu-run.The application of object oriented FEM(OOFEM)to geotechnical engineeringJ.Rock and Soil Mechanics,2000,21(4):346349.2009