线性与非线性有限元及应用-电子书

第 1 章  绪 论

1.1   有限元及其应用

在 20 世纪 50 年代，有限单元法是作为处理固体力学问题的方法出现的。追溯历史，早在 1943 年，Courant 已应用了单元的概念。在 1945 年至 1955 年，Argyris 等人在结构矩阵分析方面有了很大的发展。1956 年，Turner 和Clough 等人将刚架位移法的思路推广应用于弹性力学平面问题，将连续体划分为若干个三角形和矩形单元，单元中的位移函数首先采用了近似表达式，推导了单元刚度矩阵，建立了节点力和节点位移之间的单元刚度方程。在 1960 年，Clough 首先将这种方法命名为有限单元法。几乎与此同时， 我国的冯康也独立提出了类似的方法。随后，伴随着计算机硬件的高速发展，计算速度的不断提高，有限元方法也迅速发展，日新月异，日臻完善。

线性和非线性有限元分析是计算机辅助设计的基本组成部分。由于它提供了更快捷和低成本的方式评估设计的概念和细节，因此，人们越来越多地应用有限元仿真的方法代替样品原型的试验。例如，在汽车设计领域中，对初期设计概念和最终设计细节的评估，碰撞的仿真代替了整车的试验，如布置判定气囊释放的加速计，内部的缓冲装置，以及选择材料和满足碰撞准则的构件截面。在土木工程领域中，对桥梁结构的频率和周期计算，地震或者风雨激振的动态响应评估，以及在桥梁建造过程中变形和应力的模拟，以此指导桥梁设计和施工。在许多制造领域中，可以进行加工过程的仿真，从而加速了设计过程，例如金属薄板成型、挤压和铸造。在电子工业中，为了评估产品的耐久性，如手机和便携式计算机的设计，利用仿真分析代替了跌落试验。所有这一切均将有限元分析和仿真模拟作为必要的环节，从而使有限元技术的发展成为现代工业与工程技术领域中必不可少的组成部分。

对于有限元分析的初学者，以及应用和开发有限元程序的人员必须理解有限元分析的基本概念。若不理解基本概念，有限元程序就只是一个提供仿真的黑匣子，进而导致有限元分析者常常面对许多选择和困惑。若不理解所包含的内容和这些选择和困惑的内涵，分析者将处于非常不利的状态，甚至使分析完全偏离实际结果，给出错误结论。本书的目的是描述固体力学线性和非线性有限元分析的理论和方法。我们的意图是提供一种全面的阐述，使读者能够获得对基本方法的理解，对不同近似计算用途的比较，正确处理在应用有限元时所遇到的困难。同时，也足够详细地给出了各种列式和算法的实施过程以便读者能够编程。

有限元分析包含下列步骤：

1 ． 建立模型

2 ． 推导方程的公式

3 ． 建立离散方程

4 ． 求解方程

5 ． 表述结果

对于初学者必须掌握这 5 方面的内容。对于有限元软件的使用者，第 2 至 4 项已在典型的分析程序中实现，而分析者的工作体现在第 1 和第 5 项。

直到二十世纪八十年代，有限元建模注重提取反映力学性能的基本单元，目的是使这些基本单元能够与所复 制的研究力学性能的最简单模型一致。而现在，建模的目标已发展为建立一个单一的详细的设计模型并应用它检验所有的必要的工业准则，这在工业界已经成为非常普遍的方法。必须指出，这种建模形式的改变和提高是有限元应用的一大发展，它的动力在于，对于一种工业产品生成几种网格的成本远高于生成对每种应用都适用的特殊网格的成本。例如，同样一个便携式计算机的有限元模型可以用来进行跌落仿真动力分析、线性静力分析和热应力分析。通过使用同一个模型进行所有这些分析，节省了大量的工程研制时间。

在最近的将来，有限元模型可能成为“虚拟”的样品原型，可以用来检验设计性能的许多方面。计算机机时的节省和计算机速度的提高使得这一方式更加有效。然而，对于特殊的分析，有限元软件的使用者还必须能够评估有限元模型的适用性和限制条件。今天，计算公式的推导和模型离散主要掌握在软件开发者手中。然而，由于工业产品及工程问题的千姿百态，常用的某些方法和软件可能不合适， 一位不理解软件基本内容的分析者会面对许多风险。而且，为了将试验数据转换为输入文件，分析者必须清楚在程序中所应用的和由实验人员提供的材料应力和应变数据的度量。分析者必须理解和知道如何评估数据响应的敏感程度，必须清楚容易产生误差的来源，如何检查这些误差和评价误差的量级，以及各种算法的限制和误差影响量。如要正确理解和解决这些问题，除了必须的相应工程知识外，还必须对有限元的基本理论和方法有透彻的理解。

求解离散方程也面临许多选择。一种不恰当的选择将导致冗长的计算时间消耗，从而使分析者在规定的时间内无法获得正确结果，甚至出现振荡和发散。为了实现建立一个合理的模型和选择最佳的求解过程的良好策略，了解各种求解过程的优势和劣势以及所需要的大致计算机机时是非常必要的。

分析者最重要的任务是表述结果。除固有的近似之外，即便是线性有限元模型，分析对于许多参数常常是敏感的，它们可能将一个模拟引入歧途。非线性固体可能经历非稳态过程，其对缺陷的反应可能是非常敏感的，它们的结果可能主要取决于材料的参数。除非分析者对于这些现象非常清醒，否则很有可能错误地描述模拟的结果。

尽管有这些困惑，但是对于有限元分析的用途和前景，人们持非常乐观的态度。在许多工业领域中，非线性有限元分析已经缩短了设计周期，大大减少了原型试验的成本。由于仿真产生了各种变化的输出，很容易实现各种尝试，因此仿真能够极大地改进工程师对其产品在各种环境下的基本物理性能的理解程度。而从试验中只能得到产品是否能经受一种确定环境的大体的并且是重要的结果，如果产品不能满足要求，当产品重新设计时，这些试验通常几乎无法提供所需要依据的性能细节。另一方面，计算机仿真给出了详细的应力和应变以及其它状态变量的历史，这些数据掌握在一名优秀的有洞察力的工程师手中，将能为他提供如何重新设计产品的有价值的依据。

本书展示了对于工程和科学问题有限元求解的各种方法和技巧。然而，为了保持适合于教学的特点，我们在书中也紧密地结合了我们认为是有限元分析的核心问题的几个主题，它们是：

1 ． 如何构造各种类型的单元，以及它们的优劣；

2 ． 对于要解决的问题，如何选择近似的方法；

3 ． 对于给定的问题，如何选择合适的网格，以及动力学和运动学的描述；

4 ． 如何检验结果和求解过程的稳定性；

5 ． 如何认识模型的平滑响应和所隐含的求解质量和困难；

6 ． 如何判断主要假设的作用和误差的来源。

对于许多包含过程仿真的大变形问题和破坏分析，选择合适的网格描述是重要的，例如，是否应用 Lagrangian（拉格朗日）、Eulerian（欧拉）或者 Arbitrary Lagrangian Eulerian （任意的拉格朗日－欧拉）网格，需要认识网格畸变的影响，在选择网格时必须牢牢记住不同类型网格描述的优点。

在仿真非线性过程中，普遍存在着稳定性的问题。在数值模拟中，很可能获得物理上不稳定的因而相对无意义的解答。对于不完备的材料和荷载参数，许多解答是敏感的；在某些情况下，甚至是敏感于所采用的网格。 一位睿智的有限元应用者必须清醒这些特性和所遇到的困难，否则，由计算机仿真精心制作的结果可能是相当错误的和导致不正确的设计精度。

在非线性有限元分析中，平滑性也是普遍存在的问题。缺乏平滑性会降低大多数算法的强劲功能，并可能在结果中引入不期望的波动。目前已经发展了改进响应平滑性的技术；它们称为调整过程。然而，调整过程常常并不基于物理现象，在许多情况下难以确定与调整相关的常数。因此，分析者常常面临进退两难的窘境，是否选择导致平滑求解的方法或者处理不连续的响应，例如在接触－碰撞中的罚函数方法，以及正确评价这些方法的益处。

在分析中，结果的精度和稳定性是重要的问题。这些问题以多种方式出现。例如，在选择单元的过程中，分析者必须清楚稳定性和各种单元的锁定特性。单元的明智选择包括多种因素，例如，对于求解问题的单元的稳定性，和结果的期望平滑性，以及期望变形的量级。此外，分析者必须清楚非线性分析的复杂性。对于出现物理和数值非稳定性的可能性，必须给予密切关注和在求解过程中进行检查。

因此，精通有限元软件的应用，必定要求分析者重视对于有限元方法的理解，提供这种理解和使读者能够清楚在有限元分析中尚存在许多有兴趣的挑战和机遇，这也是本书的目的之一。

1.2  有限元的有关著作和简要历史

已经发表了一些成功的试验和专题文章，它们完全地或者部分地对有限元分析作出了贡献。仅论述非线性有限元分析的著作， 国际上有 Oden（1972），Crisfield（1991），Kleiber （1989）和 Zhong（1993）。特别值得注意的是 Oden  的书，因为它是固体和结构非线性有限元分析的先驱著作。最近的著作有 Simo 和 Hughes（1998），Bonet 和 Wood（1997）。某些著作还部分地对非线性分析作出了贡献，它们是 Belytschko 和 Hughes（1983），Zienkiewicz和 Taylor（1991），Bathe（1996），以及 Cook 、Malkus 和 Plesha（1989）。对于非线性有限元分析，这些书提供了有益的入门指南。作为姐妹篇，线性有限元分析的论述更是不胜枚举，内容最全面的著作是 Hughes（1987），Zienkiewicz 和 Taylor（1991）。

下面我们回顾有限元方法和相关软件的简单历史。在这个信息－计算机时代， 像许多其它方面的进步一样，在有限元分析中，软件常常比文献更好地代表了最新的进展。

有限元方法有多种溯源。通过波音研究组的工作和 Turner、Clough 、Martin  和 Topp（1956）的著名文章，使线性有限元分析得以闻名，不久之后，在许多大学和研究所里，工程师们开始将其扩展至非线性、小位移的静态问题。对于许多必须涉及工程问题的工程师们，他们非常清楚有限元方法的前途：它提供了一种处理复杂形状的真实问题的可能性。

在二十世纪六十年代，由于 Ed Wilson 发布了他的第一个程序，这种激 情终于被点燃了。这些程序的第一代没有名字。在遍布世界的许多实验室里，通过改进和扩展这些早期在 Berkeley 开发的软件，工程师们扩展了新的用途；这些带来了对工程分析的巨大冲击和有限元软件的随之发展。在 Berkeley 开发的第二代线性程序称之为 SAP（Structural Analysis Program）。由 Berkeley 的工作发展起来的第一个非线性程序是 NONSAP，它具有隐式积分进行平衡求解和瞬时问题求解的功能。

第一批非线性有限元方法文章的主要贡献者有 Argyris（1965）以及 Marcal  和 King （1967）。不久，一大批文章激增，而且软件随之诞生。当时在 Brown  大学任教的 Pedro Marcal，为了第一个非线性商业有限元程序进入市场，于 1969 年建立了一个公司；程序命名为 MARC，目前它仍然是主要软件。大约在同期， John  Swanson   为了核能应用在Westinghouse 发展了一个非线性有限元程序，后命名为 ANSYS。

在早期的商用软件舞台上，另外两个主要人物是 David Hibbitt 和 Klaus- Jürgen Bathe。 Hibbitt 与 Pedro Marcal 合作到了 1972 年，后来与其它人合作建立了 HKS 公司，使 ABAQUS商用软件进入市场。因为该程序是能够引导研究人员增加用户单元和材料模型的早期有限元程序之一，它对软件行业带来了实质性的冲击。Jürgen Bathe 是在 Ed Wilson 的指导下在Berkeley 获得博士学位的，不久之后开始在 MIT 任教，这期间他便发布了他的程序。这是NONSAP 软件的派生产品，称为 ADINA。

直到大约 1990  年，商用有限元程序集中在静态解答和隐式方法的动态解答。当代非线性软件的另一支血脉是显式有限元程序。DOE 实验室的工作强烈地影响了早期的显式有限元方法，Wilkins （1964），特别是命名为 hydro-codes 的软件。在 1964 年，Costantino 在芝加哥的 IIT 研究院发展了可能是第一个显式有限元程序（Costantino,  1967）。它局限于线性材料和小变形，由带状刚度矩阵乘以节点位移计算内部的节点力。

在 1969  年，研究人员开发了著名的单元乘单元的技术；节点力的计算不必应用刚度矩阵。因此，发展了名为 SAMSON 的二维有限元程序。在 1972 年，该程序功能扩展至结构的完全非线性三维瞬态分析，称为 WRECKER。这一工作使人们在七十年代初期就预言汽车的碰撞试验可能被仿真所代替。而在 Argonne，由 Belytschko 发展的显式程序被移植应用在核安全工业上，其程序命名为 SADCAT 和 WHAMS。

在 DOE 国家实验室，开始了平行的研究工作。1975 年，工作在 Sandia 的 Sam Key 完成了 HONDO，它也是具有单元乘单元功能的显式方法。程序可以处理材料非线性和几何非线性问题，并且有精心编辑的文件。得益于 Northwestern 大学的研究生 Dennis  Flanagan 的工作，这些程序得到进一步的发展，他将程序命名为 PRONTO。

显式有限元程序发展的里程碑来自于Lawrence Livermore 实验室的John Hallquist 的工作。于 1975 年，John 开始他的工作，在 1976 年，他首先发布 DYNA 程序。他慧眼吸取了前面许多人的成果，并且与 Berkeley 的研究人员紧密交流合作，包括 Jerry Goudreau ，Bob Taylor ，Tom Hughes 和 Juan Simo。他之所以成功的部分关键因素是与 Dave Benson 合作发展了接触－碰撞相互作用，和他的令人敬畏的编程效率， 以及计算程序 DYNA-2D  和DYNA-3D 的广泛传播。因此，像 Wilson 的程序和 John 的程序，不久后在全世界的大学、政府和工业实验室里到处可见。

Hallquist 关于有效接触－碰撞算法的发展（与今天的有效算法相比，是原始的第一批算法，但是仍然常常采用它们），采用一点积分单元和高阶矢量使得工程仿真得以有显著性突破的可能。矢量似乎已经与新的一代计算机无关，但是，在八十年代以 Cray 机为主的计算机上运行大型问题，矢量是至关重要的。 一点积分单元与沙漏控制的一致性，通过几乎是一阶量值的完全积分三维单元，可以提高三维分析的速度。

在八十年代，DYNA  程序首先被法国 ESI  公司商品化，命名为 PAMCRASH，它与WHAMS 也有许多相关的子程序。在 1989 年，John Hallquist 扩展LSDYNA－商业版的 DYNA程序。

在过去的十年，计算机成本的迅速下降和显式程序功能的强健带来了设计的革命。第一个主要的应用领域是有价值的汽车碰撞，然而，应用的领域迅速地扩展。在越来越多的工业领域，非线性有限元仿真正在代替样品原型的试验。产品设计依赖于仿真正常工作、跌落－试验和其它极端加载情况的帮助，例如蜂窝电话、手提电脑、洗衣机、链锯和许多其它产品。制造过程也应用有限元进行仿真，例如锻压、薄金属板成型和挤压。对于某些仿真问题，隐式方法的功能也变得越来越强，很明显，两种方法的功能都是必要的。例如，显式方法可能最适合仿真薄金属板成型的加工过程，在回弹过程模拟中，隐式方法是更合适的。

今天，隐式方法比显式方法的功能增加得更加迅速。对于处理非线性约束， 例如接触和摩擦，隐式方法已经有了明显的改进。稀疏迭代求解器也已经成为更加有效的工具。今天，具有强健功能的有限元软件需要两种方法的结合。

非线性有限元分析代表了三个领域的关联：（1）线性有限元方法，它是变分原理极值问题的一种表达形式，也可以认为是结构分析矩阵方法的扩展；（2）非线性连续介质力学；

（3）数学，包括数值分析，线性代数和泛函（Hughes ，1996）。

第 2 章  线性有限元法的一般原理和基本方程

在众多的数值方法中，有限单元法（Finite Element Method，缩写为 FEM）具有适用性强，处理均质、非均质、非线性、复杂边界问题方便， 且有计算精度高等优点。因此作为一种最适用的数值分析方法广泛地应用于工程分析中。

本章我们将通过弹性力学变分原理建立弹性力学问题有限元法（线性有限元）的表达格式。基于最小势能原理建立的有限元基本未知量是位移。以单元的节点位移作基本未知量是有限元中最常用的单元，称为位移元。为了正确阐述线性有限元法的一般原理，建立有限元方程和求解过程，将通过最为简单的平面问题 3 节点三角形常应变单元为实例加以系统讨论。

作为一种数值方法，有限元解的收敛性和精度估计必然是一个十分重要的问题。本章将简单地加以讨论，这里所阐述的一般原则将可供在以后各章中作参考或进一步具体化。

2.1  单元和形函数

有限元法的基本思想是将一个连续的求解区域离散(剖分)为一组有限个单元，并按一定方式相互联结在一起的集 合体，即用一个有限个单元的集 合体来替代几何形状相同或近似相同的连续求解区域。由于单元能按不同的联结方式进行组合， 且单元本身又可以用不同的形状描述，因此它能模拟几何形状十分复杂的连续体。同时， 它作为数值方法的一个特点是在单元内假设的一个近似函数来分片地表示求解域上待求的函数场。这样就将求解域内连续的未知的函数场（无限个自由度）转化为有限个单元离散节点的新未知量（即自由度）的函数场（有限个自由度）。如果单元尺度趋向无穷小，则有限个自由度离散的未知函数场向无限个自由度的连续未知函数场转化，近似解也就收敛于精确解。

下面我们将按有限元分析的过程先对单元剖分和形函数加以讨论。

2.1.1 离散化

采用有限元法求解问题的第一步是将一个连续的求解区域剖分为有限个单元，这一个过程习惯上称为结构离散化。以最简单的单元 3 节点三角形平面常应变单元为例，将弹性平面结构体分割为互不重叠、也没有间隙的有限个三角形， 并确定三角形的三个顶点为节点，以此来代替原来的结构。图 2-1 所示为承受内压的无限长的厚壁圆管，按弹性力学理论可以简化为一个平面问题，且利用其对称性可仅分析其四分之一。此时弹性体的曲线边界为若干段直线所代替，随着单元的增加，这种拟合将会更精确。然后， 对所有的单元和节点依照一定规则的次序从 1 开始不间断地加以编号。

2.1.2 位移模式和形函数

线性有限元法可以说是建立在最小势能原理基础上的一种近似数值方法，它以位移作为基本未知量。为此，我们必须首先建立以单元节点位移表示单元内任一点位移地近似表达式。

图 2-2 中表示一个典型的平面 3 节点三角形单元，其节点i ，j ，m 按逆时针方向排列。每个节点位移在单元平面内有两个位移分量(u, v ) ，整个单元将有六个节点位移分量（以后也称单元有六个自由度），可记为：

其中子向量：

由于考察对象是二维的小变形弹性体，单元内任意一点的位移分量是x，的函数。有限元分析中以单元节点位移作为基本未知量。为了能使单元内任意一点的应变和应力可以由单元的节点位移来表示，需要假定一个单元位移模式（单元的位移插值函数）。作为三角形单元，可以选择最为简单的线性多项式表示，即：

式中u ，v 表示单元内任一点的位移， α1 …α6 为待定常数，亦称广义坐标，且六个常数可以由单元节点位移六个分量来完全确定。

设节点i ，j ，m 的坐标分别为(xk , yk) ，其中k = i, j, m ，将它们代入（2-2）式可得：

上式第一式由k = i, i, m 可得联立的三个方程。由这三个联立方程，即可解得 3 个待定常数：

式中：

由解析几何可知， Δ 等于三角形ijm 的面积。为了使面积不为负， 节点i ，j ，m 的次序必须是逆时针转向（如图 2-2 中所示），或者对 Δ 取绝对值。

将式（2-4）代入式（2-3）第一式，稍加整理可得 

这里出现的记号表示按下标进行轮换来得到另两组表达式，以后相同。 ( , , ) i jm , , i jm

同理由式（2-3）第二式的三个联立方程可得

从以上讨论可以发现，单元内任一点位移形态完全决定于（2-2）式假定的位移模式，而位移模式中多项式阶数的选择又依据于单元节点的自由度数。从这个角度考察 3 节点三角形平面单元是最简单的一种位移模式，也是形函数中最简单的一种形式。

下面讨论形函数的几个重要性质:

1) 形函数Ni 在节点i 上的值

而在其余节点j ，m 上的值

也就是说在i 节点上Ni  = 1，在 j、m 节点上Ni  = 0 。如此类推Nj  (xi , yi) = 0,  Nj  (xj , yj) = 1,  Nj  (xm , ym) = 0

2）单元任一点上三个形函数之和等于 1，即Ni  + Nj  + Nm   = 1。它可以证明如下：

由ai  ，bi  和ci （ i ，j ，m ）定义的行列式可知，上式右端第一项圆括号内的值为2Δ ,其余两项均为零，故上式等于 1。由此可见Ni (i, j, m) 中只有两个是独立的。

同时，如果单元发生刚体 位移，假如x 方向有刚体 位移u0 ，则单元内任一点（包括节点）到处都应有位移u0 ，即 ui  = uj  = um  = u0 ，又由式（2-7）有

因此必然要求Ni  + Nj  + Nm   = 1，若形函数不满足此要求，则形函数不能反映刚体 位移要求，以此求解不能得到正确结果。

3) 三角形单元ijm的某边界上一点的形函数，与该点所在边界相对的第三个节点的坐标无关。如ij边上某一点p(x, y )有：

由此可见，在公共边 上的位移，v是完全由两个公共节点ij的位移所确定，所以相邻单元的位移是连续的，这 一点在以后讨论单元协调性是很重要的。但这里提出的性质仅对该单元的线性位移模式而 言，是否可以推广应用，应对具体 位移模式进行证明。

更多内容见附件

免责声明：

本页面/内容部分素材来源于互联网公 开 信 息，旨在传递更多信息，不代表本平台立场。

版权归原作者或机构所有，如涉及侵权，请通过平台联系我们，我们将在核实后第一时间处理。

本平台对转载内容的真实性、准确性不作任何保证，用户需自行判断并承担使用风险。