结构分析的有限元法 与MATLAB程序设计 徐荣桥 编著 （浙江大学建筑工程学院）内容简介本书以有限单元基本理论为重点，以 MATLAB 程序为平台，以工程实例为背景，介绍了有限元法及其程序设计方法。本书讲述结构分析中有限元法的基本原理，单元类型包括平面杆系，空间杆系，平面等参元，空间等参元，薄板壳单元和厚板壳单元等。内容涉及杆系结构、平面问题、空间问题和板壳问题。以结构线弹性静力分析为主，同时也讲述了结构的振动、稳定和动力响应分析。本书介绍了 MATLAB 编程环境下编写有限元程序的方法和技巧，并附有若干算例的有限元程序。还有若干推导有限元列式的 MATLAB 符号运算程序示例。这些都为读者深入理解有限元理论和掌握其实施技巧提供了极好的手段。本书可以作为高等院校土木工程专业高年级本科生或研究生有限元法结构分析的教材，也适合于其他如工程力学、机械工程等相关专业的科研人员在学习和研究工作中参考。前 言有限元法经过半个世纪的发展，在理论和实践上均取得了引人瞩目的成就，事实上它已经发展成为工程领域中一门不可或缺的技术，同时也是科技工作者进行科学研究的有力工具。因此，有限元法逐渐成为高等院校理工科专业的必修课程。本书是编者在为浙江大学桥梁与隧道专业硕士生所讲授的《结构分析的有限元法》课程讲义的基础上编写而成。本书的一大特色是采用 MATLAB 作为编程平台，利用 MATLAB强大的科学计算和符号运算功能，帮助读者轻松跨越繁琐的公式推导和复杂的编程技巧，获得最佳的学习效率。在写作上，采用理论和程序紧密结合的方法，以增加读者的感性认识，更好地理解有限元理论。不仅每章后面有极强工程背景的数值算例和源程序，而且还在有限元理论的叙述过程中，插入相应的 MATLAB 程序段以帮助理解。本书的内容共分八章，前面两章主要讲述有限元和 MATLAB 的发展概况及其基础知识。第 3 章介绍了在工程中有极大应用价值的杆系结构的各类单元，包括杆单元、平面和空间梁单元。第 4 章讲述经典的平面问题三角形单元和四边形单元。第5 章把第四章的平面单元推广到空间单元，还讨论了空间轴对称单元。第 6 章介绍了应用最广的平面和空间等参数单元以及轴对称等参数单元。第 7 章在简要介绍板理论的基础上，讨论了三角形和矩形薄板单元，也给出了基于 Mindlin 板理论的四变形板单元。最后介绍了用平板单元和平面应力单元组合而成的壳单元。第 8 章讲述了结构的动力和稳定问题，并介绍了几种常用的特征值问题和动力响应问题的算法。在写作过程中，本书得到了浙江大学丁皓江教授和徐兴教授的大力指导，详细审阅了初稿，指出并修改了若干错误，在此表示衷心的感谢。亦得到了浙江大学项贻强教授、陈伟球教授、叶贵如教授、蔡金标副教授、赵阳讲师和张治成博士的大力指导和帮助。我要特别感谢我的家人，感谢他们在我写作此书时所做的牺牲。由于编者水平有限，时间仓卒，书中难免有许多缺点甚至错误，热情欢迎专家和读者批评指正。徐荣桥2005 年 12 月 第一章 绪论许多工程问题，我们虽然已经得到了他们的基本方程和边界条件及初始条件，但是能用解析方法求解的只是少数方程简单、边界规则的问题，而绝大多数只能通过其他途径解决。随着计算机硬件和软件的发展，数值方法越来越受到人们的青睐，其中有限元法以其方法的统一性和理论的普遍性而独领风 骚，已经成为处理各类科学和工程问题的有效方法之一，亦有众多的有限元商用软件流行，如 ANSYS，NASTRAN，ADINA，ABAQUS 等，它们包含众多单元类型，能求解各类问题。有限元法的基本思想是“化整为零”，其实这个思想很早就在各个领域被人们采用。老子就说“道生一，一生二，二生三，三生万物”，他把世间万物全部分解成最基本的“道”。古代数学家将圆用多边形近似，并以此估算圆周率 π 值，达到了 40 位数字的精度。有限元法把一个复杂的结构分解成相对简单的“单元”，各单元之间通过结点相互连接。单元内的物理量由单元结点上的物理量按一定的假设内插得到，这样就把一个复杂结构从无限多个自由度简化为有限个单元组成的结构。我们只要分析每个单元的力学特性，然后按照有限元法的规则把这些单元“拼装”成整体，就能够得到整体结构的力学特性。可见有限元“化整为零”的思想十分简单明了。但更重要的是，有限元可以建立在严格的数学基础上[ 1]，成为求解微分方程的标准方法之一，从而它不仅能够解决结构分析问题，也能解决工程中如电磁场、流体力学、热传导、渗流等领域的诸多问题[2]，因为它们在数学上都能用微分方程来描述。有限元法是一种数值方法，应用性很强，因此在学习的过程中，必要的程序设计、编写和调试能加深理解有限元理论，增加从实际工程出发建立有限元模型的感性认识，极大地提高学习的兴趣和效率。传统的有限元书籍，凡给出具体程序代码的[2,3,4] ，基本上都是以FORTRAN 为程序设计语言。由于一个完整的有限元计算程序，要包括模型的输入，稀疏矩阵的存储，求解，分析结果的后处理等，导致程序十分庞大难懂，尤其对于初学者。因此有的就只用一个求解单元刚度矩阵的子程序来说明程序的编写，而这样一段程序无法独立编译连接成一个可执行文件，读者就无法验证其算例获得感性认识，示例效果大大降低。本书的特色就是采用目前最流行的科学分析计算软件 MATLAB 作为编程环境。 MATLAB 的优势在于采用矩阵作为它的基本数据类型，并提供大量实用的矩阵运算函数，可以把有限元理论中用矩阵表示的公式直接写成形式上一样的程序代码，它还有专门的极易使用的大型稀疏矩阵存储和求解的功能，从而使得一个完整的有限元程序十分简洁明了，方便读者理解和验证，大大提高了示例程序的效果。目前虽也已经有采用 MATLAB 作为程序设计语言的有限元专著[5] 出版，并被翻译成中文[6]。但该书以介绍 MATLAB 编写的有限元程序为主，几乎不涉及有限元理论。另一本结合 MATLAB 编写的有限元专著[7] 的情况有所不同，它更倾向于把有限元作为一个求解偏微分方程的数值方法来介绍，目前还未见中译本。1.1 有限元法简介有关有限元的起源，很多关于有限元的专著[2,3,8,9, 10, 11]都会提到 Hrenikoff[ 12]，Courant[ 13]， McHenry[ 14] ，Newmark[ 15] ，Turner et al[ 16]和 Clough[17] 的工作。但是，我们同样也应该提到Ritz[ 18]，他在 1909 年曾提出一个求解连续体力学问题的非常有用的近似方法，后人称为 Ritz法。该方法把待求函数用一组已知的试函数的加权和来表示，通过最小势能原理可求解每个试函数的权系数。 Ritz 法最基本的缺点就是每一个试函数必须满足指定的边界条件。 Courant[ 13]对 Ritz 法作了重要的推广，他把整个求解区域分成很多三角形的子区域，然后在子区域内假定待求函数为线性函数，把三角形顶点处的函数值作为未知数，并且把边界条件放宽，只要求在边界的有限点上满足，这样，应用 Ritz 法的一个很大的困难被克服了。更为重要的是，原先的 Ritz 法解的精度很到程度上取决于试函数的选取，这要求应用人员具有丰富的工程经验。后来人们发现，Courant[ 13]应用的 Ritz 法其实与若干年后由 Clough[ 17]提出的有限元法是一致的。当然，在 1960 年有限元获得迅速发展的原因在于该法中大量的数值运算能够由当时正发展的电子计算机来实现，而在 1943 年，还没有这个工具为 Courant[ 13]所用。随着计算机硬件和软件的发展，有限元迅速发展，并渐趋成熟，目前它已经被推广至三维问题，非线性问题，时变问题，甚至已经超越了结构分析领域，如流体流动，热传导和电磁场分析等。应用有限元法进行结构分析时，应把所分析结构物离散成有限个单元，并在每个单元上指定有限个结点，单元通过这些结点连接构成整个有限元模型，用来模拟所分析的实际结构。同时选定所求物理量的结点值，例如结点位移作为基本未知量。然后对每个单元假设一个简单的插值函数（称为形函数），近似地表示未知量在单元内的分布规律，再利用变分原理或其他方法，建立单元结点力和结点位移之间的关系，得到一组以结点位移为未知量的代数方程组，从而求解结点的位移分量。一经解出，就可以利用插值函数确定单元内任意一点的位移值。显然，如果单元满足问题的收敛性要求，那么随着缩小单元的尺寸，增加求解区域内单元的数量，解的近似程度将不断改进，最终收敛于精确解。当然，有限单元法中并不一定要求取结点位移作为基本未知量，也可以取结点内力为未知量，因而，随着所取未知量的不同，有所谓位移法、力法、杂交法和混合法之分。本书采用最为普遍的位移法，介绍结构分析中的限单元法基本理论和方法。有限元法概念浅显，容易掌握，可以在不同的水平上建立起对该法的理解，即可以通过非常直观的物理解释理解，也可以建立基于严格的数学分析的理论[ 1, 10, 19,20]。它不仅对结构物的复杂几何形状有很强的适应性，也能应用于结构物的各种物理问题，如静力问题，动力问题，非线性问题，热应力问题等。还能处理非均质材料、各向异性材料、以及复杂边界条件等难题。因而，有限元法已经被公认为工程分析的有效工具，受到普遍重视。有限元法还有一个特点是，它的理论采用矩阵形式表达。这并不利于一般的计算机语言编制计算机程序，因为传统的计算机语言处理的对象是标量，使用矩阵形式的有限元理论时，必须把矩阵形式的公式转换成标量表示的公式。而如果采用 MATLAB，这个特点就变成了有限元法的优点。这也是本书采用 MATLAB 作为编程环境的一个重要原因。1.2 MATLAB 简介MATLAB 是当今国际科学界最具影响力和活力的软件。它起源于矩阵运算，并已经发展成一种高度集成的计算机语言。它提供了强大的科学计算，灵活的程序设计流程，高质量的图形可视化与界面设计，便捷的与其他程序和语言接口的功能。MATLAB 在各国高校与研究单位起着重大的作用。MATLAB 语言的首创者 Cleve Moler 教授在数值分析，特别是在数值线性代数的领域中很有影响，他参与编写了数值分析领域两个重要的 Fortran 程序包 EISPACK 和 LINPACK 。他曾在密西根大学、斯坦福大学和新墨西哥大学任数学与计算机科学教授。1980 年前后，当时的新墨西哥大学计算机系主任 Moler 教授在讲授线性代数课程时，发现了用其他高级语言编程极为不便，便构思并开发了 MATLAB(MATrix LABoratory，即矩阵实验室)，这一软件利用了当时数值线性代数领域最高水平的 EISPACK 和 LINPACK 两大软件包中可靠的子程序，用 Fortran 语言编写了集命令翻译、科学计算于一身的一套交互式软件系统。所谓交互式语言，是指人们给出一条命令，立即就可以得出该命令的结果。该语言无需像 C 和 Fortran 语言那样，首先要求使用者去编写源程序，然后对之进行编译、连接，最终形成可执行文件。这无疑会给使用者带来极大的方便。早期的 MATLAB 是用 Fortran 语言编写的，只能作矩阵运算；绘图也只能用极其原始的方法，即用星号描点的形式画图；内部函数也只提供了几十个。但即使其当时的功能十分简单，当它作为免费软件出现以来，还是吸引了大批的使用者。后来，Cleve Moler 和 John Little 等人成立了 MathWorks 公司，Cleve Moler 一直任该公司的首席科学家。该公司于 1984 年推出了第一个 MATLAB 的商业版本。当时的 MATLAB版本已经用 C 语言作了完全的改写，其后又增添了丰富多彩的图形图像处理功能、多媒体功能、符号运算和与其它流行软件的接口功能，使得 MATLAB 的功能越来越强大。MathWorks 公司于 1992 年推出了具有划时代意义的 MATLAB 4.0 版本，并于 1993 年推出了其微机版，可以配合 Microsoft Windows 一起使用，使之应用范围越来越广。1994 年推出的 4.2 版本扩充了4.0 版本的功能，尤其在图形界面设计方面更提供了新的方法。1997 年推出的 MATLAB 5.0 版允许了更多的数据结构，如单元数据、数据结构体、多维矩阵、对象与类等，使其成为一种更方便编程的语言。1999 年初推出的 MATLAB 5.3 版在很多方面又进一步改进了 MATLAB 语言的功能。2000 年 10 月底推出了其全新的 MATLAB 6.0 正式版(Release 12)，在核心数值算法、界面设计、外部接口、应用桌面等诸多方面有了极大的改进。最近又推出了 MATLAB7.0。虽然 MATLAB 语言是计算数学专家倡导并开发的，但其普及和发展离不开自动控制领域学者的贡献。甚至可以说，MATLAB 语言是自动控制领域学者和工程技术人员捧红的，因为在 MATLAB 语言的发展进程中，许多有代表性的成就和控制界的要求与贡献是分不开的。迄今为止，大多数工具箱也都是控制方面的。MATLAB 具有强大的数学运算能力、方便实用的绘图功能及语言的高度集成性，它在其他科学与工程领域的应用也是越来越广，并且有着更广阔的应用前景和无穷无尽的潜能。“工欲善其事，必先利其器”。如果有一种十分有效的工具能解决在教学与研究中遇到的问题，那么 MATLAB 语言正是这样的一种工具。它可以将使用者从繁琐、无谓的底层编程中解放出来，把有限的宝贵时间更多地花在解决问题中，这样无疑会提高工作效率。目前，MATLAB 已经成为国际上最流行的科学与工程计算的软件工具，现在的MATLAB已经不仅仅是一个“矩阵实验室”了，它已经成为了一种具有广泛应用前景的全新的计算机高级编程语言了，有人称它为“第四代”计算机语言，它在国内外高校和研究部门正扮演着重要的角色。MATLAB 语言的功能也越来越强大，不断适应新的要求提出新的解决方法。可以预见，在科学运算、自动控制与科学绘图领域 MATLAB 语言将长期保持其独一无二的地位。1.3 本书的目的和内容本书的目的是把结构分析有限元法的基本概念、基本理论和基本的程序设计方法介绍给土木工程和交通工程专业的高年级本科生和研究生。然而， 由于各种概念是以非常简单的形式给出的，本书对其他背景的学生和实际工程人员也是有帮助的。本书也适合于那些想把有限元应用于解决实际工程问题的各类人员。内容涉及杆系单元，平面单元，三维空间单元，等参数单元，板壳单元，包括静力问题，动力问题，稳定问题。并给出了丰富的以土木工程，尤其是桥梁与隧道工程为背景的算例，可以帮助读者增加如何从工程实际抽象到有限元模型的感性认识。参考文献[1] Brenner S. C., Scott L. R. 1993 The Mathematical Theory of Finite Element Methods, Springer-Verlag.[2] Zienkiewicz O. C. 1977, The Finite Element Method, Third Edition. McGraw-Hill Inc.[3] 王勖成 2003 有限单元法 清华大学出版社[4] 凌道盛 徐兴 2004 非线性有限元及程序 浙江大学出版社[5] Kattan P. I. 2003 MATLAB Guide to Finite Element Springer-Verlag[6] Kattan P. I. (韩来彬 译) 2003 MATLAB 有限元分析与应用 清华大学出版社[7] Kwon Y. W., Bang H. 1997 The Finite Element Method using MATLAB, CRC Press, Inc. Boca Raton, Florida[8] Bathe K. J., Wilson E. L. 1976 Numerical Methods in Finite Element An alysis, Prentice-Hall, Inc.[9] Cook R. D. 1982, Concepts and Applications of Finite Element Ana lysis. John Wiley & Sons Inc.[10] Logan D. L. 1993, A First Course in the Finite Element Method, Second Edition, PWS Publishing Company, ITP[11] 龙驭球 龙志飞 岑松 2004 新型有限元论 清华大学出版社[12] Hrenikoff A. 1941, Solution of problems in elasticity by the framework method, J. Appl. Mech. A8:169-175.[13] Courant R. 1943, Variational methods for the solution of problems of equilibrium and vibration. Bull. Am. Math. Soc., 49:1-23[14] McHenry D. 1943, A lattice ana logy for the solution of plane stress problems, J. Inst. Civ.Eng. 21: 59-82.[15] Newmark N. M. 1949, Numerical methods of an alysis in bars plates and elastic bodies, Numerical Methods in Ana lysis in Engineering (ed. L. E. Grinter), Macmiilan[16] Turner M. J., Clough R. W., Matrin H. C. and Topp L. J., 1956, Stiffness and deflection ana lysis of complex structures, J. Aero. Sci., 23: 805-23[17] Clough R. W. 1960, The finite element in plane stress an alysis, Proc. 2nd ASCE Conf. on Electronic Computation, Pitts burgh, Pa., Sept.[18] Ritz W. 1909 Über eine neue Methode zur Lösung gewissen Variations-Probleme der mathematischen Physik, J. Reine Angew Math., 135:1-61[19] 丁皓江 何福保 谢贻权 徐兴 1989 弹性和塑性力学中的有限单元法 机械工业出版社[20] Reddy J. N. 1993 An Introduction to the Finite Element Method, Second Edition, McGraw-Hill.[21] Bathe K. J. 1982 Finite Element Procedures in Engineering An alysis, Prentice-Hall, Inc.第二章 有限元法的预备知识本章主要介绍有限元法的一些预备知识，包括矩阵和线性代数的基本概念，MATLAB编程基础，弹性力学中的变分原理和能量原理基础以及有限元法的基本过程。如果对这些内容十分熟悉，可跳过本章，直接阅读第三章。2.1 矩阵、线性代数和 MATLAB在数学上，有限元法通过变分法或加权残数法把微分方程化为近似的代数方程组进行求解。如果采用矩阵形式进行中间过程的推导，就能使有限元公式简洁紧凑。为此， 我们必须掌握矩阵和线性代数的基本知识，了解矩阵的详细操作，这非常有助于阅读和理解有限元法的基本理论并组织有效的算法。这一节的目标是扼要地介绍以后要用到的矩阵和线性代数的基本知识。我们可以发现，其实这些都是非常基础的。如果希望了解更详细的内容，可以参考有关矩阵和线性代数的专著[ 1,2,3]。在叙述矩阵、线性代数基本知识的同时，将介绍使用 MATLAB 实现矩阵定义，运算等基本操作，这有助于增加感性认识，加深理解。当然，这里给出的 MATLAB 语言的用法都是最基本的，详细的说明可以阅读 MATLAB 软件随带的帮助文档，如图 2-1 所示。图 2-1 MATLAB 的帮助窗口另外也有许多关于介绍 MATLAB 的专著和网站，后面的参考文献中列出一些比较著名的，读者也可以用网络搜索工具比如 Google 等在互联网上查找更新更多的有关 MATLAB 的介绍和说明资料[4-20]。2.1.1 矩阵简介在实际计算中矩阵的表达效率可以从求解线性方程组的过程中看出，例如 式中，未知数为δ1 ， δ2 ， δ3 和δ4 。使用矩阵符号，这组方程可写成 在这里，我们把未知数的系数按一定的顺序排成一个阵列，左边的未知数( δ1 ， δ2 ， δ3 和δ4 ）和右边的已知数也各自被排成一列。虽然写法不同，式(2.1)和式(2.2)表示了同样的内容。可以进一步把式(2.2)写成Aδ = b (2.3)式中，A 是线性方程组的系数矩阵， δ 是未知数矩阵，b 是已知数矩阵，即 显然，式(2.3)比式(2.1)或(2.2)简洁得多。而且，可以任意指定矩阵A 和未知数δ 及已知数b的阶数，用来表示任意阶数的线性方程组。现在矩阵的正式定义看上起非常明显了，即矩阵是按一定顺序排列的数据阵列。在 MATLAB 中，定义矩阵变量十分方便和直观，例如式(2.4)中的矩阵 A，可以这样定义 >>A=[ 1.4 0.0 -1.0 -0.2 0.0 1.4 -0.4 0.4 -1.0 -0.4 1.8 0.6 -0.2 0.4 0.6 2.4 ] 这里的符号“>> ”是 MATLAB 的命令提示符，表示等待输入命令状态。输入上述命令后，在 MATLAB的命令窗口立刻显示 我们说这个矩阵是n m×阶矩阵，在不引起歧义的情况下，可简写为A。当矩阵只有一行（1=m）或一列（1=n）时，我们称它为行向量或列向量。本书统一用大写的黑体字母表示矩阵或行向量，用小写的黑体字母表示列向量。另外，用ija表示矩阵A的第i行第j列的元素。2.1.2 特殊矩阵1、零矩阵元素全部为零的m × n 阶矩阵，称为零矩阵，记为0m×n 或者简写为0 。MATLAB 中定义零矩阵十分方便，比如定义一个 4×3 的零矩阵的命令是命令 zeros 接受可接受任意个参数，一般我们使用的都是 2 维的矩阵，前一个参数表示行数，后一个是列数。2、对称矩阵我们把m × n 阶矩阵A 的转置，记为 AT ，它由A 中的行和列互换得到。如果m = n ，则称A 为n 阶方阵。如果n 阶方阵 A 的元素aij = aji ，则称 A 为对称矩阵。因此对称矩阵必须是方阵。如果 A 是对称矩阵，则有 A = AT 。显然，式(2.4)中的矩阵 A 是对称矩阵。例如设这里需要指出的是，函数 transpose 与单引号“’”对于实矩阵的操作是完全一样的，而对于复矩阵则不同。因为复矩阵在数学上的转置定义不仅包括行列互换，还要执行共轭操作。函数 transpose 只执行行列互换的操作，而单引号则是数学意义上的转置。3、单位矩阵另一个特殊矩阵是单位矩阵I n 。它是n 阶方阵，除了对角线上的元素为1外，其余元素都为零。例如，一个 3 阶的单位阵为在实际计算时，单位矩阵的阶数常常是隐含的。与单位矩阵相似，我们定义一个n 阶的单位列向量ei ，式中i 表示该向量是单位矩阵中的第i 列。MATLAB 中用 eye 来定义单位矩阵。也许读者可能奇怪为什么不用 I，这是因为小写的 i 已经被用来定义复数单位，而大写的 I 容易跟它混淆，而且 MATLAB 中内部的变量和函数名全部是用小写的。下面的命令定义了一个三阶的单位矩阵4、带状稀疏矩阵在有限元中，我们将接触到大量的对称带状稀疏矩阵，因为一般情况下，刚度矩阵和质量矩阵都是带状稀疏矩阵。带状矩阵是指矩阵中位于带宽以外的元素全为零。如果带状矩阵A 是对称的，我们可以把这种情况用算式表示为式中2HBW+1是矩阵A 的带宽。作为一个例子，下述矩阵是一个 6 阶的对称带状矩阵，它的半带宽HBW是 2：如果一个矩阵的半带宽为零，该矩阵只有位于对角元的元素非零，此时称它为对角矩阵。例如单位矩阵就是对角矩阵。如果我们用普通的矩阵变量来存储这种矩阵，将浪费大量的空间。因为稀疏矩阵中的很多零元素不参加运算，也不用存储。由于计算机的内存有限，因此我们必须合理地安排存储空间，使得在计算中要用到的数据都能放入内存，否则频繁地读取外存，将使计算速度大大降低。对于这种带状稀疏矩阵，我们应该组织一种合理的存储方法，使得内存中能容下高阶稀疏矩阵。图 2-2 所示为一对称带状稀疏矩阵。在轮廓线以外的元素为零，而且在以后的运算中永远为零，因此我们不必存储。而在轮廓线以内的元素必须存储，即使某一元素可能暂时为零，但是经过运算后就会变成非零，因此一定要存储。它的一维存储方法如图 2-2 所示。在有限元程序中，整体刚度矩阵就属于这类矩阵，它们阶数很大，但是大部分元素为零。在 MATLAB 中，有专门处理这类矩阵的工具。这样就极大地节省了我们处理稀疏矩阵的时间，使我们可以绕开很多烦人的细节而专心于有限元程序设计的主要方面，这对于初学者尤其重要。因为这些问题往往是他们编写有限元程序的主要障碍，而不是有限元的基本理论。但是目前的 MATLAB 中，稀疏矩阵还有一个缺点，就是不能指定它是对称的，因此将浪费不少内存。更多内容见附件免责声明：本页面/内容部分素材来源于互联网公 开 信 息，旨在传递更多信息，不代表本平台立场。版权归原作者或机构所有，如涉及侵权，请通过平台联系我们，我们将在核实后第一时间处理。本平台对转载内容的真实性、准确性不作任何保证，用户需自行判断并承担使用风险。
