有限单元法综述

有限单元法开展与趋势综述摘要：本文简述了有限单元法的提出过程，开展历史及根本思想，重点论述了有限单元法的根本解题步骤及有限单元法的开展趋势。文章最后介绍了一个基于ANSYS软件运用有限单元法分析某隧道的实例关键词：有限单元法；开展趋势；ANSYS；隧道ReviewoftheFiniteElementMethod’sDevelopmentsandTrendsAbstract：Thispaperbrieflydiscussestheproposedprocessofthefiniteelement(FEM),thedevelopmentoftheFEMandthebasicideaoftheFEM,withemphasisonthebasicproblem-solvingstepsoftheFEMandthedevelopmenttrendoftheFEM.Finally,thepaperanalyzesatunnelproblembyusingtheFEM,whichisbasedontheANSYSsoftwareplatform.Keywords：thefiniteelement；developmenttrends；ANSYS；tunnel引言随着现代工业、生产技术的开展，不断要求设计高质量、高水平的大型、复杂和精密的机械及工程结构。为此目的，人们必须预先通过有效的计算手段，确切的预测即将诞生的机械和工程结构，在未来工作时所发生的应力、应变和位移。但是传统的一些方法往往难以完成对工程实际问题的有效分析。弹性力学的经典理论，由于求解偏微分方程边值问题的困难，只能解决结构形状和承受载荷较简单的问题，对于几何形状复杂、不规那么边界、有裂缝或厚度突变，以及几何非线性、材料非线性等问题往往遇到很多麻烦，试图按经典的弹性力学方法获得解析解是十分困难的，甚至是不可能的。因此，需要寻求一种简单而又精确的数值分析方法。有限单元法正是适应这种要求而产生和开展起来的一种十分有效的数值计算方法。有限单元法概述有限单元法早在40年代初期就已有人提出，但当时由于没有计算工具而搁置，一直到50年代中期，高速数字电子计算机的出现和开展为有限单元法的应用提供了重要的物质条件，才使有限单元法得以迅速开展。，R.W.Clough,H.C.Martin和L.J.Topp等人于1956年发表了一篇题为《复杂结构的刚度和挠度分析》一文，此文提出了计算复杂结构刚度影响系数的方法，说明了如何利用计算机进行分析。美国教授R.W.Clough于1960年在一篇介绍平面应力分析的论文中，首次提出了有限单元法的名字。1965年英国的Zienliewice教授及其合作者解决了将有限元应用于所有场的问题，使有限单元法的应用范围更加广泛。有限单元法的优点很多，其中最突出的优点是应用范围广。开展至今，不仅能解决静态的、平面的、最简单的杆系结构，而且还可以解决空间问题、板壳问题、结构的稳定性问题、动力学问题、弹塑性问题和粘弹性问题、疲劳和脆性断裂问题以及结构的优化设计问题。而且不管物体的结构形式和边界条件如何复杂，也不管材料的性质和外载荷的情况如何，原那么上都能应用。有限单元法的根本思想有限单元法的根本思想，是在力学模型上将一个原来连续的物体离散成为有限个具有一定大小的的单元，这些单元仅在有限个节点上相连接，并在节点上引进等效力以代替实际作用于单元上的外力。对于每个单元，根据分块近似的思想，选择一种简单的函数来表示单元内位移的分布规律，并按弹性理论中的能量原理〔或用变分原理〕建立单元节点力和节点位移之间的关系。最后，把所有的单元的这种关系式集合起来，就得到一组以节点位移为未知量的代数方程组，解这些方程组就可以求出物体上有限个离散节点上的位移。有限单元法的根本步骤有限单元法分析计算的根本步骤可归纳为以下五点：1.结构的离散化结构的离散化是有限单元法分析的第一步，它是有限元法的根底。将某个机械结构划分为由各种单元组成的计算模型，这一步称作单元划分。离散后单元与单元之间利用单元的节点相互连接起来，将求解区域变成为用点、线或面划分的有限数目的单元组合成的集合体。单元的形状原那么上是任意的。例如，在平面问题中通常采用三角形单元，有时也采用矩形或任意四边形单元。在空间问题中，可以采用四面体、长方体或任意六面体单元。可见，不管单元取什么样的形状，在一般情况下，单元的便捷总不可能与求解区域的真实边界完全吻合，这就带来了有限元法的一个根本近似性——几何近似。在一个具体的结构中，确定单元的类型和数目以及哪些部位的单元可以取得大一些，哪些部位单元应该取得小一些，需要由经验来做出判断。单元划分越细那么描述变形情况越精确，即越接近实际变形，但计算量越大。所以有限元法中分析的结构已不是原有的物体或结构物，而是同样材料的众多单元以一定方式连接成的离散物体。这样，用有限元分析计算所得的结果只是近似的。如果划分单元数目非常多而又合理，那么所获得的计算结果就越逼近实际情况。1)网格划分根底与划分原那么复杂结构的离散是有限元分析的根底，也决定着计算结果的精确度。一个复杂的结构总可以离散为一维、二维、三维的小单元。当然对二维和三维单元，其离散后的形状可以为任意的，但是为了计算的方便性和精确性的结合，二维单元一般采用三角形和四边形，而三维单元那么采用四面体和六面体。简单的说，复杂结构的离散就是网格的划分。有限元网格[3]的划分有很多原那么，一是网格数量，网格数量直接影响计算精度和计算时耗,网格数量增加会提高计算精度,但同时计算时耗也会增加。当网格数量较少时增加网格计算精度可明显提高,但计算时耗不会有明显增加;当网格数量增加到一定程度后,再继续增加网格时精度提高就很小,而计算时耗却大幅度增加。所以在确定网格数量时应权衡这两个因素综合考虑。二是网格密度，为了适应应力等计算数据的分布特点,在结构不同部位需要采用大小不同的网格。如在孔的附近有集中应力，因此网格需要加密，周边应力梯度相对较小，网格划分较稀。该网格反映了疏密不同的网格划分原那么：在计算数据变化梯度较大的部位。为了较好地反映数据变化规律，需要采用比拟密集的网格；而在计算数据变化梯度较小的部位，为减小模型规模，网格那么应相对稀疏。三是单元阶次，单元阶次与有限元的计算精度有着密切的关联，单元一般具有线性、二次和三次等形式，其中二次和三次形式的单元称为高阶单元。高阶单元的曲线或曲面边界能够更好地逼近结构的曲线和曲面边界，且高次插值函数可更高精度地逼近复杂场函数，所以增加单元阶次可提高计算精度。但增加单元阶次的同时网格的节点数也会随之增加，在网格数量相同的情况下由高阶单元组成的模型规模相对较大,，因此在使用时应权衡考虑计算精度和时耗。四是网格形状，网格单元形状的好坏对计算精度有着很大的影响，单元形状太差的网格甚至会中止计算。在网格划分时应保证合理的单元形状，即使只有一个单元形状很差或畸形时,也可能给计算结果带来很大的误差，甚至使得计算无法进行下去。2〕对网格的评价单元形状评价一般有以下几个指标：(1)单元的边长比、面积比或体积比以正三角形、正四面体、正六面体为参考基准，理想单元的边长比为一,线性单元可接受的边长比小于三,二次单元小于十。(2)扭曲度:单元面内的扭转和面外的翘曲程度。(3)节点编号:节点编号对于求解过程中总刚矩阵的带宽和波前因数有较大的影响,从而影响计算时耗和存储容量的大小。因此合理的节点编号有利于刚度矩阵对称、带状分布等求解效率,从而提高计算速度。3〕不同维数模型划分介绍我们对各维模型的单元划分做简要的介绍。一维单元可分为两种。一类是单元的节点参数中只包含场函数的节点值C0型，另一类是单元的节点参数中，除场函数的结点值外，还包含场函数导数的节点值的C1型单元。这分别是拉格朗日单元和Hermite单元。也就是说拉格朗日是一次插值单元，而后者是二次插值，这样就能保证导数的连续性，也就是能保证在连接处除了位移连续，连接的交点也是光滑的。对二维单元，可以采用三角形和四边形单元。对三角形单元，如同一维单元的情形，可以利用总体笛卡尔坐标，也可以利用无量纲的局部自然坐标以构造三角形单元的插值函数。利用总体笛卡尔坐标构造三结点三角形单元的差值函数较复杂，更普遍采用的是局部自然坐标来直接构造一般三角行单元的差值函数，这时运算比拟简单。三角形单元的插值一般采用面积坐标，把一个三角形用线段分成等分块，由插值函数的性质等可以推导出差值函数。通常情况下，采用矩形单元比三角形单元更为方便而有效。其差值函数的推导和一维情况也很相似，也可以构造二维的拉格朗日矩形单元和Hermite矩形单元。此时后者的精度同样比拉格朗日单元的精度要高。由于有时四边形单元的节点在矩形内部，所以一个偶然的发现，Serendipity四边形单元被发现，这个单元有很多优势，一方面由于在实际用用中优势希望统一单元的不同边界有不同数目的节点，这样可以实现不同阶次单元之间的过渡，从而可能在求解的不同区域采用不同精度的单元，另一方面通过它阐述构造单元插值函数的一般方法。三维单元可能有的几何形状要比二维单元多得多，在应用中只讨论几种常用的形状，又因为构造其插值函数的方法只是二维的推广，所以其形式是很容易构造出来的。其四面体单元也可以用体积坐标，同时也存在Serendipity单元。2.单元分析1〕选择位移模式位移模式是表示单元内任意点的位移随位置变化的函数式，由于所采用的函数是一种近似的试函数，一般不能精确地反映单元中真实的位移分布，这就带来了有限元法的另一种根本近似性。采用位移法时，物体或结构物离散化之后，就可把单元中的一些物理量如位移、应变和应力等由节点位移来表示。这时可以对单元中位移的分布采用一些能逼近原函数的近似函数予以描述。通常，有限单元法中我们将位移表示为坐标变量的简单函数，这种函数称为位移模式或位移函数，如式中，是待定系数；是与坐标有关的某种函数。2〕建立单元刚度矩阵选定单元的类型和位移模式以后，就可按虚功原理或最小势能原理建立单元刚度方程，它实际上是单元各个节点的平衡方程，其系数矩阵称为单元刚度矩阵式中，为单元编号；为单元的节点位移向量；为单元的节点力向量；为单元刚度矩阵，它的每一个元素都反映了一定的刚度特性。根据单元的材料性质、形状、尺寸、节点数目、位置及其含义等，找出单元节点力和节点位移的关系式，这是单元分析中的关键一步。此时需要应用弹性力学中的几何方程和物理方程来建立力和位移的方程式，从而导出单元刚度矩阵，这是有限元法的根本步骤之一。3〕计算等效节点力物体离散化后，假定力是通过节点从一个单元传递到另一个单元。但是，对于实际的连续体，力是从单元的公共边界传递到另一个单元中去的。因而，这种作用在单元边界的外表力、体积力或集中力都需要等效的移到节点上去，也就是用等效的节点力来代替所有作用在单元上的力。3.整体分析有限元法的分析过程就是先分后合，即先进行单元分析，在建立了单元刚度方程以后，再进行整体分析，把这些方程集成起来，形成求解区域的刚度方程，称为有限元位移法根本方程。集成所遵循的原那么是各相邻单元在共同节点处具有相同的位移。利用结构力的平衡条件和边界条件把各个单元按原来的结构重新连接起来；形成整体的有限元方程式中，K为整体结构的刚度矩阵；为整体节点位移向量；F为整体载荷向量。4.求解方程，得出节点位移解有限元方程式得出位移、这里，可以根据方程组的具体特点来选择适宜的计算方法。有限元方程是一组方程组，方程在经历平衡问题中就是以节点位移为根本未知量的系统结点平衡方程。有限元求解的效率及计算结果的精确很大成都上取决于线性代数方程组的解法。特别是随着研究对象的更加复杂，有限元分析需要采用更多单元的离散模型来近似实际结构或力学问题的几何构形时，线性代数方程租的阶数就愈来愈高。因而，线性方程组采用何种有限的方法求解，以保证求解的效率和精度就称为更加重要的问题。不仅在线性静力分析中，求解代数方程组的时间在整个解题时间中占有很大比重，而且在动力分析和非线性分析中这局部比重也是相当大的。假设不采用适当的求解方法，不仅计算费用大量增加，严重时可能导致求解过程的不稳定和求解的失败。线性代数方程组的解法可以分为两大类，即直接解法和迭代解法。直接解法的特点是，选定某种形式的直接解法以后，对于一个给定的线性代数方程组，事先可以按规定的算法步骤计算出它所需要的算数运算操作数，直接给出最后的结果。迭代解法的特点是，对于一个给定的线性代数方程租，首先假设一个初始解，然后按一定的算法共识进行迭代。在每次迭代过程中对解的误差进行检查，并通过增加迭代次数不断降低解的误差，直至满足解的精度要求，并输出最后的解答。迭代解法的优点之一是，它不要求保存洗漱矩阵中高度轮廓线以下的零元素，并且不对它们进行运算，即它们保持为零不变。这样一来，计算机只需存储洗漱矩阵的非零元素以及记录它们位置的辅助数组。这不仅可以最大限度地节约了存储空间，而且提高了计算效率。另一方面，迭代解法在计算过程中可以对解的误差进行检查，并通过增加迭代次数来降低误差，直至满足解的精度要求。其缺乏之处是，每一种迭代算法可能只适合某一类问题，常缺乏通用的有效性，如使用不当，可能会出现迭代收敛很慢，甚至不收敛的情况。5.由节点位移计算单元的应变与应力解出节点位移以后，根据需要，可由弹性力学的几何方程和物理方程来计算应变和应力。通过上述分析，可以看出，有限单元法的根本思想是“一分一合”，化整为零，集零为整，把复杂的结构看成由有限个单元组成的整体。有限单元法的开展趋势有限单元法已经成为现代力学领域分析问题的一个最重要的途径，为了方便用户的使用和适应问题复杂性的要求，目前有限单元法开展方向主要集中在以下几个方面：更为强大的网格处理能力有限元法求解问题的根本过程主要包括：分析对象的离散化、有限元求解、计算结果的后处理三局部。由于结构离散后的网格质量直接影响到求解时间及求解结果的正确性与否，近年来各软件开发商都加大了其在网格处理方面的投入，使网格生成的质量和效率都有了很大的提高，但在有些方面却一直没有得到改良，如对三维实体模型进行自动六面体网格划分和根据求解结果对模型进行自适应网格划分，除了个别商业软件做得较好外，大多数分析软件仍然没有此功能。自动六面体网格划分是指对三维实体模型程序能自动的划分出六面体网格单元，现在大多数软件都能采用映射、拖拉、扫略等功能生成六面体单元，但这些功能都只能对简单规那么模型适用，对于复杂的三维模型那么只能采用自动四面体网格划分技术生成四面体单元。对于四面体单元，如果不使用中间节点，在很多问题中将会产生不正确的结果，如果使用中间节点将会引起求解时间、收敛速度等方面的一系列问题，因此人们迫切的希望自动六面体网格功能的出现。自适应性网格划分是指在现有网格根底上，根据有限元计算结果估计计算误差、重新划分网格和再计算的一个循环过程。对于许多工程实际问题，在整个求解过程中，模型的某些区域将会产生很大的应变，引起单元畸变，从而导致求解不能进行下去或求解结果不正确，因此必须进行网格自动重划分。自适应网格往往是许多工程问题如裂纹扩展、薄板成形等大应变分析的必要条件。由单一结构场求解开展到耦合场问题的求解有限元分析方法最早应用于航空航天领域，主要用来求解线性结构问题[5]，实践证明这是一种非常有效的数值分析方法。而且从理论上也已经证明，只要用于离散求解对象的单元足够小，所得的解就可足够逼近于精确值。现在用于求解结构线性问题的有限元方法和软件已经比拟成熟，开展方向是结构非线性、流体动力学和耦合场问题的求解。例如由于摩擦接触而产生的热问题，金属成形时由于塑性功而产生的热问题,需要结构场和温度场的有限元分析结果交叉迭代求解，即"热力耦合"的问题。当流体在弯管中流动时，流体压力会使弯管产生变形，而管的变形又反过来影响到流体的流动，这就需要对结构场和流场的有限元分析结果交叉迭代求解，即所谓"流固耦合"的问题。由于有限元的应用越来越深入，人们关注的问题越来越复杂，耦合场的求解必定成为有限元的开展方向。由求解线性问题开展到求解非线性问题随着科学技术的开展，线性理论已经远远不能满足设计的要求，许多工程问题如材料的破坏与失效、裂纹扩展等仅靠线性理论根本不能解决，必须进行非线性分析求解，例如薄板成形就要求同时考虑结构的大位移、大应变〔几何非线性〕和塑性〔材料非线性〕；而对塑料、橡胶、陶瓷、混凝土及岩土等材料进行分析或需考虑材料的塑性、蠕变效应时那么必须考虑材料非线性。众所周知，非线性问题的求解是很复杂的，它不仅涉及到很多专门的数学问题，还必须掌握一定的理论知识和求解技巧，学习起来也较为困难。为此国外一些公司花费了大量的人力和物力开发非线性求解分析软件，如ADINA、ABAQUS等。它们的共同特点是具有高效的非线性求解器、丰富而实用的非线性材料库，ADINA还同时具有隐式和显式两种时间积分方法。程序面向用户的开放性随着商业化的提高，各软件开发商为了扩大自己的市场份额，满足用户的需求，在软件的功能、易用性等方面花费了大量的投资，但由于用户的要求千差万别，不管他们怎样努力也不可能满足所有用户的要求，因此必须给用户一个开放的环境，允许用户根据自己的实际情况对软件进行扩充，包括用户自定义单元特性、用户自定义材料本构、用户自定义边界条件、用户自定义结构断裂判据和裂纹扩展规律等等。与图形软件如CAD软件的无缝集成当今有限元分析软件的一个开展趋势是与通用CAD软件的集成使用，即在用CAD软件完成部件和零件的造型设计后，能直接将模型传送到CAE软件中进行有限元网格划分并进行分析计算，如果分析的结果不满足设计要求那么重新进行设计和分析，直到满意为止，从而极大地提高了设计水平和效率。为了满足工程师快捷地解决复杂工程问题的要求，许多商业化有限元分析软件都开发了和著名的CAD软件的接口。基于ANSYS的计算实例分析题目:某双线铁路隧道，长700米，除进出口洞段〔约2100m〕岩体为碎石状松散结构体，岩体级别为V级外，其余500m洞段主要为块碎石状镶嵌结构体，岩体级别为IV级，埋深50m~150m不等。现根据《铁路隧道设计标准》(TB10003—2005)规定，IV级岩体段隧道衬砌采用复合衬砌，其中的初期支护及二次衬砌采用工程类比法进行设计，并通过数值方法进行验算。初衬和二衬设计参数见表1，围岩及支护结构物理力学参数见表2。表1铁路双线隧道IV级围岩复合式衬砌设计参数围岩级别初期支护二次衬砌(C30)喷射混凝土(C20)厚度(cm)锚杆(直径25mm)钢筋网钢架拱、墙(cm)仰拱(cm)拱、墙仰拱位置长度(m)间距(m)IV2020拱、墙3.01.0拱、墙、仰拱@2525—4045表2围岩及支护结构根本物理力学参数工程密度/kg/m3弹性(变形)模量/GPa泊松比内摩擦角/粘聚力/MPaIV级围岩22003.20.32350.5C20喷砼2500210.2——C30混凝土2500300.2——锚杆78002000.3——图1隧道开挖轮廓为两段圆弧(见图1)，顶拱及边墙段圆弧圆心在隧道轴线上，半径5.5m，圆心角约223；仰拱段圆弧由顶拱及边墙段圆弧的两个角点及仰拱顶点组成，仰拱顶点位于顶拱及边墙段圆弧圆心以下3.8m处。施工采用CRD工法。计算当隧道埋深分别为50m、100m和150m〔至隧道轴线〕时，初衬内力和围岩的应力分布特征，并进行比拟。计算方法:对于双线铁路隧道，采用CRD工法施工，其衬砌包括初期支护和二次衬砌。本次计算中采用开挖数值模拟〔地层结构模型〕来分别计算隧道埋深为50m，100m，150m时初衬内力和围岩的应力分布。地层结构模型需要用到ansys软件中平面单元〔Plane42〕模拟围岩和杆单元〔Link〕模拟锚杆，采用Beam3单元模拟喷射混凝土和钢筋网。锚杆布置于边墙轴线以上的局部，喷混凝土为封闭结构，不考虑钢筋网对喷砼的加强作用。计算过程:荷载计算：隧道埋深分别为50m，100m，150m时，模拟面均取为50m埋深。100m，150m时模拟面以上的岩体以荷载形式施加，顶面垂直荷载为上部岩体〔未模拟局部〕自重，侧边水平荷载，按弹性力学方法计算，等于垂直荷载乘侧压力系数μ/(1-μ)，此题中侧压力系数为0.4706。由于该铁路隧道地下水不丰富，按水土合算侧压力。模拟的地层范围，横向两端延伸至4倍洞径即44m，高度为埋深50m，隧道底下3倍洞高即31.7m。整个模拟地层为99m×81.7m。〔1〕．前处理主要包括材料、实常数和单元类型定义，建立几何模型和单元网格划分。具体见命令流。生成的关键点，生成的线