有限元学习xyt.doc

有限元分析学习心得

有限单元法是20世纪50年代以来随着电子计算机的广泛应用而发展起来的有一种数值解法。有限元分析（FEA，Finite‌Element‌Analysis）的基本概念是用较简单的问题代替复杂问题‌有限元分析‌后再求解。‌它将求解域看成是由许多称为有限元的小的互连子域组成，‌对每一单元假定一个‌合适的(较简单的）近似解，然后推导求解这个域总的满足条件(如结构的平衡条件），从而‌得到问题的解。这个解不是准确解，而是近似解，因为实际问题被较简单的问题所代替。
有限单元法[1]是应用局部的近似解来建立整个定义域的解的一种方法。先把注意力集中在单个单元上，进行上述所谓的单元分析。基本前提是每一单元要尽可能小，以致其边界值在整个边界上的变化也是小的。这样，边界条件就能取某一在结点间插值的光滑函数来近似，在单元内也容易建立简单的近似解。因此，比起经典的近似法，有限元法具有明显的优越性。比如经典的Ritz法，要求选取一个函数来近似描述整个求解区域中的位移，并同时满足边界条件，这是相当困难的。而有限元法采用分块近似，只需对一个单元选择一个近似位移函数，且不必考虑位移边界条件，只须考虑单元之间位移的连续性即可。对于具有复杂几何形状或材料、荷载有突变的实际结构，不仅处理简单，而且合理适宜。
有限单元法，是一种有效解决数学问题的解题方法。在有限元方法中，把计算域离散剖分为有限个互不重叠且相互连接的单元，在每个单元内选择基函数，用单元基函数的线形组合来逼近单元中的真解，整个计算域上总体的基函数可以看为由每个单元基函数组成的，则整个计算域内的解可以看作是由所有单元上的近似解构成。在河道数值模拟中[2]，常见的有限元计算方法是由变分法和加权余量法发展而来的里兹法和伽辽金法、最小二乘法等。根据所采用的权函数和插值函数的不同，有限元方法也分为多种计算格式。从权函数的选择来说，有配置法、矩量法、最小二乘法和伽辽金法，从计算单元网格的形状来划分，有三角形网格、四边形网格和多边形网格，从插值函数的精度来划分，又分为线性插值函数和高次插值函数等。不同的组合 同样构成不同的有限元计算格式。对于权函数，伽辽金(Galerkin)法是将权函数取为逼近函数中的基函数 ；最小二乘法是令权函数等于余量本身，而内积的极小值则为对代求系数的平方误差最小；在配置法中，先在计算域 内选取N个配置点 。令近似解在选定的N个配置点上严格满足微分方程，即在配置点上令方程余量为0。
插值函数一般由不同次幂的多项式组成，但也有采用三角函数或指数函数组成的乘积表示，但最常用的多项式插值函数。有限元插值函数分为两大类，一类只要求插值多项式本身在插值点取已知值，称为拉格朗日(Lagrange)多项式插值；另一种不仅要求插值多项式本身，还要求它的导数值在插值点取已知值，称为哈密特(Hermite)多项式插值。单元坐标有笛卡尔直角坐标系和无因次自然坐标，有对称和不对称等。常采用的无因次坐标是一种局部坐标系，它的定义取决于单元的几何形状，一维看作长度比，二维看作面积比，三维看作体积比。在二维有限元中，三角形单元应用的最早，近来四边形等参元的应用也越来越广。对于二维三角形和四边形电源单元，常采用的插值函数有Lagrange插值直角坐标系中的线性插值函数及二阶或更高阶插值函数、面积坐标系中的线性插值函数、二阶或更高阶插值函数等。
在工程计算过程中，对于许多力学问题，人们可以给出他们的数学模型，即基本方程和定解条件。但能用解析方法求出精确解的只是少数方程性质比较简单，并且几何形状要非常规则。对于大多数的问题，由于几何形状的不规则等原因，只能采用数值分析的方法。随着计算机的广泛应用，有限单元法已经成为求解复杂问题的一条很适用的方法。
已经发展的偏微分方程数值分析方法可以分为两类。一类以有限差分法为代表，主要应用在流体问题的分析。而另一类即是有限单元法。有限单元法区别与传统的加权余量法和求解泛函驻值法，该法不是在整个求解域上假设近似函数，而是在各个单元上分片假设近似函数。这样就克服了在全域上假设近似函数所遇到的困难，是近代工程仿真分析方法领域的重大突破。
有限单元法一开始就对一个连续体用有限个（然而是大量的）坐标或自由度来近似地（然而是系统的）加以描绘。一个离散化的结构可由许多结构单元组成，这些单元仅在有限个结点上彼此铰结。每一单元所受的已知体力和面力都按静力等效原则移置到结点上，成为结点荷载。计算通常采用位移法，取结点的未知位移分量{δ}e为基本未知量。为了在求得结点位移后可求得应力，必须建立单元中应力与结点位移的关系，由应力转换矩阵[S]表达。
有限单元法基本方程的推导有很多途径，被广泛接受的是变分法，即结合最小势能原理推导有限单元法的过程。由最小势能原理可以推导下列方程式：

故

可得

利用弹性力学的几何方程写出单元应变与结点位移的关系矩阵，称应变矩阵[B]，即
再由材料的本构关系（即物理方程），得到单元弹性矩阵[D]，从而推出用结点位移表示单元应力表达式

其中，[S] = [D][B]。
然后考虑结点平衡求得单元结点力与结点位移的关系，由矩阵[k]e表示，称单元刚度矩阵。根据虚功原理或最小势能原理（平衡条件），也可导出用结点位移表示结点力的表达式
其中，单元刚度矩阵
利用虚功原理（或变分原理）可同时导出单元等效结点力{F}e。
在经逐个单元（逐个结点）叠加其贡献予以集合（整体分析）后，生成结构刚度矩阵[K]（也称总刚）、荷载列阵{F}和结构结点位移列阵{δ}，并利用平衡条件建立表达结构的力-位移的关系式，即所谓结构刚度方程：

考虑几何边界条件作适当修改后，求解上式所示的高阶线性代数方程组，得到结构所有的未知结点位移（同矩阵位移法）。最后利用已求出的结点位移计算各个单元的应力，并经后处理软件整理、显示计算结果。
从上面的理论推导过程可以总结出有限单元法分析问题的步骤的几个部分：对结构进行离散、生成单元的刚度矩阵和等效节点荷载矩阵、集成结构的刚度矩阵和等效节点荷载列阵、引入强制的边界条件、求解有限元求解方程，得到节点位移、计算单元应变和应力。有限元中要解决的问题也就是在这几个方面。
复杂结构的离散是有限元分析的基础，也决定着计算结果的精确度。一个复杂的结构总可以离散为一维、二维、三维的小单元。当然对二维和三维单元，其离散后的形状可以为任意的，但是为了计算的方便性和精确性的结合，二维单元一般采用三角形和四边形，而三维单元则采用四面体和六面体。
简单的说，复杂结构的离散就是网格的划分。有限元网格[3]的划分有很多原则，一是网格数量，网格数量直接影响计算精度和计算时耗, 网格数量增加会提高计算精度, 但同时计算时耗也会增加。当网格数量较少时增加网格计算精度可明显提高, 但计算时耗不会有明显增加; 当网格数量增加到一定程度后, 再继续增加网格时精度提高就很小, 而计算时耗却大幅度增加。所以在确定网格数量时应权衡这两个因素综合考虑。二是网格密度，为了适应应力等计算数据的分布特点, 在结构不同部位需要采用大小不同的网格。如在孔的附近有集中应力，因此网格需要加密，周边应力梯度相对较小，网格划分较稀。该网格反映了疏密不同的网格划分原则：在计算数据变化梯度较大的部位。为了较好地反映数据变化规律，需要采用比较密集的网格； 而在计算数据变化梯度较小的部位，为减小模型规模，网格则应相对稀疏。三是单元阶次，单元阶次与有限元的计算精度有着密切的关联，单元一般具有线性、二次和三次等形式，其中二次和三次形式的单元称为高阶单元。高阶单元的曲线或曲面边界能够更好地逼近结构的曲线和曲面边界，且高次插值函数可更高精度地逼近复杂场函数，所以增加单元阶次可提高计算精度。但增加单元阶次的同时网格的节点数也会随之增加，在网格数量相同的情况下由高阶单元组成的模型规模相对较大,，因此在使用时应权衡考虑计算精度和时耗。四是网格形状，网格单元形状的好坏对计算精度有着很大的影响，单元形状太差的网格甚至会中止计算。在网格划分时应保证合理的单元形状，即使只有一个单元形状很差或畸形时,也可能给计算结果带来很大的误差，甚至使得计算无法进行下去。
单元形状评价一般有以下几个指标：( 1) 单元的边长比、面积比或体积比以正三角形、正四面体、正六面体为参考基准，理想单元的边长比为一, 线性单元可接受的边长比小于三, 二次单元小于十。( 2) 扭曲度: 单元面内的扭转和面外的翘曲程度。( 3) 节点编号: 节点编号对于求解过程中总刚矩阵的带宽和波前因数有较大的影响, 从而影响计算时耗和存储容量的大小。因此合理的节点编号有利于刚度矩阵对称、带状分布等求解效率, 从而提高计算速度。
我们对各维模型的单元划分做简要的介绍。一维单元可分为两种。一类是单元的节点参数中只包含场函数的节点值C0型，另一类是单元的节点参数中，除场函数的结点值外，还包含场函数导数的节点值的C1型单元。这分别是拉格朗日单元和Hermite单元。也就是说拉格朗日是一次插值单元，而后者是二次插值，这样就能保证导数的连续性，也就是能保证在连接处除了位移连续，连接的交点也是光滑的。
对二维单元，可以采用三角形和四边形单元。对三角形单元，如同一维单元的情形，可以利用总体笛卡尔坐标，也可以利用无量纲的局部自然坐标以构造三角形单元的插值函数。利用总体笛卡尔坐标构造三结点三角形单元的差值函数较复杂，更普遍采用的是局部自然坐标来直接构造一般三角行单元的差值函数，这时运算比较简单。三角形单元的插值一般采用面积坐标，把一个三角形用线段分成等分块，由插值函数的性质等可以推导出差值函数。
通常情况下，采用矩形单元比三角形单元更为方便而有效。其差值函数的推导和一维情况也很相似，也可以构造二维的拉格朗日矩形单元和Hermite矩形单元。此时后者的精度同样比拉格朗日单元的精度要高。由于有时四边形单元的节点在矩形内部，所以一个偶然的发现，Serendipity四边形单元被发现，这个单元有很多优势，一方面由于在实际用用中优势希望统一单元的不同边界有不同数目的节点，这样可以实现不同阶次单元之间的过渡，从而可能在求解的不同区域采用不同精度的单元，另一方面通过它阐述构造单元插值函数的一般方法。
三维单元可能有的几何形状要比二维单元多得多，在应用中只讨论几种常用的形状，又因为构造其插值函数的方法只是二维的推广，所以其形式是很容易构造出来的。其四面体单元也可以用体积坐标，同时也存在Serendipity单元。
单元刚度矩阵的一般表达式为

其中，B是应变矩阵，D是材料弹性矩阵，V是单元体积。
考虑单元存在的初应力和初应变情况，单元等效节点荷载列阵的一般表达式为

其中，等式后面的四项分别是和作用与单元的体积力f，边界分布力T，单元内的初应力、出应变等效的节点荷载列阵。
在形成单元矩阵和等效荷载矩阵的过程中，由于划分后的单元形状并不规则，所以要用到等参元的概念，即单元的几何形状和单元内的场函数采用相同数目的节点参数及相同的插值函数进行变换。需要注意，由于等参元的位移场是以斜角或曲线坐标表达的，而母单元的位移场是以正则坐标表达的。因此，等参元的位移分布与母单元的位移分布即使在位移相同的情况下也可能是不同的。由于单元位移是在局部坐标系（ξ,η）下描述的，而以位移求应变所用的公式是以整体坐标系（x,y）表达的。要将形函数利用坐标变换式写成整体坐标的显式一般十分困难，需要将对直角坐标的求导运算变换成对斜角或曲线坐标的求导运算，因此引入雅可比变换矩阵J（Jacobi matrix）。等参元当=0时将失效，离散时要注意避免。
建立了单元位移场后，可按统一的单元分析步骤推导出单元刚度矩阵。在等参元的单元刚度矩阵和等效结点荷载的计算公式中，被积函数十分复杂，很难精确积分得到解析结果，要采用数值积分方法，有限元分析中通常采用高斯积分法，只用较少的积分点就能达到较高的精度，从而节省计算时间。
等参元的过程可以通过可以通过四边形的变换来了解：
首先建立规整形状的母单元，如取边长为2的正方形单元（如图5.3.4a所示），其形心处设局部坐标ξOη，得到如同（5.3.19）式的形函数，作坐标变换：

(a) 母单元 (b) 等参元
图1 四结点平面等参元
使得图5.3.4a中的ξη平面上的4个角点分别映射成图1b中xy平面上的4个点，其坐标为xi，yi（i=1，2，3，4）。由于形函数Ni是双线性的，ξη平面上的正方形被映射到xy平面上以xi、yi为角点的四边形。所以上面的坐标变换式起了把xy平面上的所有四边形单元（称子单元）都映射到ξη平面上的正方形单元（称母单元）的作用。同时，还可以把ξOη坐标看成为子单元的局部坐标,该局部坐标系是用一组不超过1的无量纲数来定出单元中的点，单元各边的方程分别是ξ=±1和η=±1。
假如四边形单元（子单元）的位移模式也采用别的形函数，即
(5.3.21）
可以证明,收敛准则中的完备性和协调性要求能够得到满足。由于上述单元的位移模式和坐标变换式采用相同的形函数（即母单元形函数Ni），故称之为等参数单元（或等参元）。
通过等参元的变换，任何形状的单元刚度都可以计算出来，对于六面体单元原来的刚度计算式如下：

此时应该变化为：

在生成单元的刚度矩阵和结点等效荷载矩阵后，应集成结构的总体刚度矩阵和总体荷载矩阵，以方便方程的求解。其集成方程式如下：

其中是直接作用于结点的集中力。
一般边界条件有三种形式，分为本质边界条件(狄里克雷边界条件 )、自然边界条件(黎曼边界条件)、混合边界条件(柯西边界条件)。对于自然边界条件， 一般在积分表达式中可自动得到满足。对于本质边界条件和混合边界条件，需按一定法 则对总体有限元方程进行修正满足。
在实际建模过程中, 合理地划分单元网格很重要, 合理地从计算对象中提取出各类边界条件同样重要。如果认为有限元建模的主要技术在于如何把握单元质量, 而忽视边界条件的处理技术,其有限元模型不仅有可能成为“残次品”, 还有可能带来不易发现的危险。
在建立的方程中，如果不引入边界条件，将无法求解，也就是说将有无数的结果，引入边界条件就使方程有唯一解，也就让问题的求解成为可能。
有限元方程是一组方程组，方程在经历平衡问题中就是以节点位移为基本未知量的系统结点平衡方程。有限元求解的效率及计算结果的精确很大成都上取决于线性代数方程组的解法。特别是随着研究对象的更加复杂，有限元分析需要采用更多单元的离散模型来近似实际结构或力学问题的几何构形时，线性代数方程租的阶数就愈来愈高。因而，线性方程组采用何种有限的方法求解，以保证求解的效率和精度就称为更加重要的问题。
不仅在线性静力分析中，求解代数方程组的时间在整个解题时间中站友很大比重，而且在动力分析和非线性分析中这部分比重也是相当大的。若不采用适当的求解方法，不仅计算费用大量增加，严重时可能导致求解过程的不稳定和求解的失败。
线性代数方程组的解法可以分为两大类，即直接解法和迭代解法。直接解法的特点是，选定某种形式的直接解法以后，对于一个给定的线性代数方程组，事先可以按规定的算法步骤计算出它所需要的算数运算操作数，直接给出最后的结果。
迭代解法的特点是，对于一个给定的线性代数方程租，首先假设一个初始解，然后按一定的算法共识进行迭代。在每次迭代过程中对解的误差进行检查，并通过增加迭代次数不断降低解的误差，直至满足解的精度要求，并输出最后的解答。迭代解法的优点之一是，它不要求保存洗漱矩阵中高度轮廓线以下的零元素，并且不对它们进行运算，即它们保持为零不变。这样一来，计算机只需存储洗漱矩阵的非零元素以及记录它们位置的辅助数组。这不仅可以最大限度地节约了存储空间，而且提高了计算效率。另一方面，迭代解法在计算过程中可以对解的误差进行检查，并通过增加迭代次数来降低误差，直至满足解的精度要求。其不足之处是，每一种迭代算法可能只适合某一类问题，常缺乏通用的有效性，如使用不当，可能会出现迭代收敛很慢，甚至不收敛的情况。
在求解过程中，计算机数据的存储方法也是很重要的，在有限单元法中，线性代数方程组的系数矩阵是对称的，因此可以只存储一个上三角矩阵。但是由于矩阵的稀疏性，仍然会发生零元素占绝大多数的情况。考虑到非零元素的分布呈带状特点，在计算机中系数矩阵的存储一般采用二维等带宽存储或一维变带宽存储。
对于n阶的系数矩阵，若取最大的半带宽D为带宽，则上三角阵中的全部非零元素都将包括在这条以主对角为一边的一条等带宽中。二维等带宽存储就是将这样一条带中的元素，以二维数组形式存储在计算机中。采用二维等带宽存储，消除了最大带宽以外的全部零元素，较之于存全部上三角阵大大节省了内存。但是由于取最大带宽为存储范围，因此它不能排除在带宽范围内的零元素。当系数矩阵的带宽变化不大时，采用二维等带宽存储是合适的，求解也是方便的。但当出现局部带宽特别大的情况时，采用二维等带宽存储时将由于带宽过大而使整个系数矩阵的存储大大增加，此时可以采取一维变带宽存储。
一维变带宽存储就是将变化的带宽内的元素按一定的顺序存储在一维数组中。由于它不按最大带宽存储，因此较二维等带宽存储更能节省内存。按照解法可分为按行一维变带宽存储及按列一维变带宽存储。一维变带宽存储是比二维等带宽存储更节省内存的一种存储方法，但由于寻找元素较二维等带宽存储复杂，因而编写程序亦较麻烦，并且计算机耗时可能也比二维等带宽存储时要多。因此在选用存储方式上要权衡两者的利弊，统盘考虑。通常，当带宽变化不大而计算机内存又允许时，采用二维等带宽存储是合适的。
方程的解最终得到的是结点位移，而在实际工程中人们还关心着应变和应力，这些值都可以用弹性力学[4]的方法从位移中推导出来。其方程式如下：

得到应力与应变后，有限元的求解也就进入末段。
有限单元法已经成为现代力学领域分析问题的一个最重要的途径，为了方便用户的使用和适应问题复杂性的要求，目前有限单元法发展方向主要集中在以下几个方面：
当今有限元分析软件的一个发展趋势是与通用CAD软件的集成使用，即在用CAD软件完成部件和零件的造型设计后，能直接将模型传送到CAE软件中进行有限元网格划分并进行分析计算，如果分析的结果不满足设计要求则重新进行设计和分析，直到满意为止，从而极大地提高了设计水平和效率。为了满足工程师快捷地解决复杂工程问题的要求，许多商业化有限元分析软件都开发了和著名的CAD软件的接口。
有限元法求解问题的基本过程主要包括：分析对象的离散化、有限元求解、计算结果的后处理三部分。由于结构离散后的网格质量直接影响到求解时间及求解结果的正确性与否，近年来各软件开发商都加大了其在网格处理方面的投入，使网格生成的质量和效率都有了很大的提高，但在有些方面却一直没有得到改进，如对三维实体模型进行自动六面体网格划分和根据求解结果对模型进行自适应网格划分，除了个别商业软件做得较好外，大多数分析软件仍然没有此功能。自动六面体网格划分是指对三维实体模型程序能自动的划分出六面体网格单元，现在大多数软件都能采用映射、拖拉、扫略等功能生成六面体单元，但这些功能都只能对简单规则模型适用，对于复杂的三维模型则只能采用自动四面体网格划分技术生成四面体单元。对于四面体单元，如果不使用中间节点，在很多问题中将会产生不正确的结果，如果使用中间节点将会引起求解时间、收敛速度等方面的一系列问题，因此人们迫切的希望自动六面体网格功能的出现。自适应性网格划分是指在现有网格基础上，根据有限元计算结果估计计算误差、重新划分网格和再计算的一个循环过程。对于许多工程实际问题，在整个求解过程中，模型的某些区域将会产生很大的应变，引起单元畸变，从而导致求解不能进行下去或求解结果不正确，因此必须进行网格自动重划分。自适应网格往往是许多工程问题如裂纹扩展、薄板成形等大应变分析的必要条件。
随着科学技术的发展，线性理论已经远远不能满足设计的要求，许多工程问题如材料的破坏与失效、裂纹扩展等仅靠线性理论根本不能解决，必须进行非线性分析求解，例如薄板成形就要求同时考虑结构的大位移、大应变（几何非线性）和塑性（材料非线性）；而对塑料、橡胶、陶瓷、混凝土及岩土等材料进行分析或需考虑材料的塑性、蠕变效应时则必须考虑材料非线性。众所周知，非线性问题的求解是很复杂的，它不仅涉及到很多专门的数学问题，还必须掌握一定的理论知识和求解技巧，学习起来也较为困难。为此国外一些公司花费了大量的人力和物力开发非线性求解分析软件，如ADINA、ABAQUS等。它们的共同特点是具有高效的非线性求解器、丰富而实用的非线性材料库，ADINA还同时具有隐式和显式两种时间积分方法。
有限元分析方法最早应用于航空航天领域，主要用来求解线性结构问题[5]，实践证明这是一种非常有效的数值分析方法。而且从理论上也已经证明，只要用于离散求解对象的单元足够小，所得的解就可足够逼近于精确值。现在用于求解结构线性问题的有限元方法和软件已经比较成熟，发展方向是结构非线性、流体动力学和耦合场问题的求解。例如由于摩擦接触而产生的热问题，金属成形时由于塑性功而产生的热问题,需要结构场和温度场的有限元分析结果交叉迭代求解，即"热力耦合"的问题。当流体在弯管中流动时，流体压力会使弯管产生变形，而管的变形又反过来影响到流体的流动，这就需要对结构场和流场的有限元分析结果交叉迭代求解，即所谓"流固耦合"的问题。由于有限元的应用越来越深入，人们关注的问题越来越复杂，耦合场的求解必定成为有限元的发展方向。
随着商业化的提高，各软件开发商为了扩大自己的市场份额，满足用户的需求，在软件的功能、易用性等方面花费了大量的投资，但由于用户的要求千差万别，不管他们怎样努力也不可能满足所有用户的要求，因此必须给用户一个开放的环境，允许用户根据自己的实际情况对软件进行扩充，包括用户自定义单元特性、用户自定义材料本构、用户自定义边界条件、用户自定义结构断裂判据和裂纹扩展规律等等。
有限元法综合了力学理论、数值计算、计算程序、工程结构等多方面知识，运用力学分析方法解决各种实际工程问题和进行科学研究。其强大的功能和适用性使对它的研究充满着挑战性和重要性。
三维实体结构分析属于弹性力学空间问题，空间弹性力学问题一般情况下很难得到解析解，有限单元法是解决这类问题的有效手段。三维实体结构分析常见于各类机械零部件的力学分析，建筑结构构建节点的分析中也有时用到。
对于弹性力学的空间问题，有三个唯一分量，六个应力（或应变）分量，但其分析方法与平面问题或轴对称问题完全雷同。也要经过结构离散化为单元，找单刚，组总刚，引入边界条件，求解以及结果分析等环节。下面将以一个实例加以说明。材料特性：弹性模量为30*e6,泊松比为0.
单元类型：选择Structural Solid和Brick 20node 95.
网格划分：激活智能网格划分器，对模型进行四面体单元的自由网格划分.
理论和实践都已证明：为了有限单元法的解答在单元的尺寸逐步取小时能够收敛于正确解答，反映刚体位移和常量应变是必要条件，加上反映相邻单元的位移连续性，就是充分条件。
有限单元法在将来的工作生活中有着重要的作用，它的功能如此强大，‌前景是很美好的，是值得我们好好用心学习和研究的。

有限元学习心得

土木10-17班
20103956
徐洋特