《有限元方法理论及其应用》综合考试.doc

有限元方法理论及其应用综合考试1 等参单元及其应用1.1 概述通常用的一些三角形、矩形、六面体单元都是形状很规则的单元。对于形状规则的连续体，用这些单元来离散可以获得比较好的结果。但是，对于一些几何形状比较复杂的连续体，再用这些单元离散就比较困难。因为单元的节点数越多，单元精度越高，因此在这一点上，矩形单元优于3节点三角形单元，六面体单元优于四面体单元；其次是单元几何上的限制，上述矩形和六面体单元都不能模拟任意形状几何体，所有几种单元都是直线边界，处理曲边界几何体误差较大。于是，有人就想出了坐标变换的方法来解决这个问题。通过一一对应的坐标变换，把规则的单元转变成形状不规则的单元，就可以用它们来离散几何形状复杂的连续体。 解决上述矛盾的出路就是突破矩形单元和六面体单元几何方面的限制，使其成为任意四边形和任意六面体单元，如果再增加边中间节点，还可以成为曲边四边形和曲面六面体高精度实用单元。但是这类单元位移模式和形函数的构造和单元列式的导出不能沿用前面构造简单单元的方法，必须引入所谓的等参变换，采用相同的插值函数对单元的节点坐标和节点位移在单元上进行插值。这种单元称为等参单元。图1为一个4节点任意四边形单元，单元有8个自由度。将矩形单元放松为4节点任意四边形单元将带来许多好处。但在建立单元位移模式时产生了新的问题：单元上没有一个如矩形单元中的简单直接的局部坐标系，而又不能直接用x，y坐标系下的双线性位移模式。须建立一种新的局部坐标系-（如图），使得4条边的坐标为常数（1），则在-平面内，单元是一个边长为2的正方形。同时，该局部坐标系的建立在x-y平面上的任意四边形单元与-平面上的正方形之间形成了一个1-1对应的映射关系。图1 4节点任意四边形单元及其母单元称-平面内的正方形单元为基本单元或母单元。x-y平面内的任意四边形单元称为实际单元。显然，母单元的节点相应于不同的x， y坐标就得到不同的任意四边形单元。建立了局部坐标系或映射后，我们只需要在-平面上的母单元中描述实际单元的位移模式和力学特性。任意四边形单元在母单元中的位移模式（或者称为-坐标系下的位移模式）与矩形单元相同： (i=1,2,3,4)当然，该位移模式在x，y坐标系下不是双线性位移模式，位移沿单元边界线性变化，能保证单元的协调性。为了得到上述映射的数学表达，引入对母单元节点上x，y坐标进行插值的思想，将母单元上每一点对应的x，y坐标看成是对节点坐标的插值，插值函数与位移插值中的形函数相同：这样就得到了一个事实上的映射，只要证明该映射确实把母单元映射成为实际单元，就是所需要的映射。该映射是用母单元描述实际单元力学特性的桥梁。由于该几何变换式中采用了与位移模式相同的参数（插值函数），因此称为等参变换。而所有采用等参变换的单元都称为等参单元。为了使实际问题物理坐标系内的单元刚度、质量、阻尼、载荷等特性矩阵的计算也能在局部的自然坐标表示的规则域内进行计算，还需要研究这些矩阵积分式内被积函数中所涉及的导数、体积微元、面积微元、线段微元的变换以及积分限的置换。实现了这种变换和置换，则不管各个积分式中的被积函数如何复杂，都可以方便地采用标准化的数值积分方法进行计算，从而使各类不同工程实际问题的有限元分析纳入统一的通用化的程序。借助于等参元可以对于一般的任意几何形状的工程问题和物理问题方便地进行有限元离散。因此，等参元的提出为有限元法成为现代工程实际领域最有效的数值分析方法迈出了重要的一部。1.2 等参单元的数值积分1.2.1 等参单元刚度矩阵的数值积分方法通常都用数值积分代替函数积分，即在单元内选出某些点，称为积分点，算出被积函数在这些积分点处的函数值，然后用对应的加权系数乘亡这些函数值，再求出总和，将其作为近似的积分值。1.2.1.1 一维数值积分首先构造一个多项式,使在（i=1，2，3，n）上有=，然后用近似函数的积分来近似原被积函数的积分。称为积分点或取样点。积分点的数目和位置决定了近似的程度，因而也就决定了数值积分的精度。对于n个积分点，按照积分点位置的不同选择，通常采用两种不同的数值积分方案，即Newton-Cotes积分方案和高斯积分方案。 Newton-Cotes积分对于n个积分点，根据积分点上的被积函数值可以构造一个近似多项式，使在积分点上有= (i=1,2,n)上式用拉格朗日多项式表示=其中是n-1阶拉格朗日插值函数。由于拉格朗日插值函数有如下性质： 的积分为并令则可得=式中为积分的权系数。Newton-Cotes积分中，积分点的位置按等间距分布，即用近似，可以写出+式中为余项。引入一些变量后 高斯积分在此积分方案中，积分点不是等间距分布。定义n次多项式由下列条件确定积分点的位置（i=1,2,n-1）有以下性质：a. 在积分点上。b.多项式与在域内正交被积函数可由2n-1次多项式来近似，即用近似，可的结果：1.2.1.2 二维和三维高斯积分对于二维问题的积分首先令为常数，进行内层积分用同样的方法进行外层积分 类似地，对于三维数值积分，则有1.2.1.3 三维六面体单元的Irons积分关于三维积分Irons给出的公式是 1.2.1.4 二维三角形单元和三维四面锥单元的Hammer积分对于二维积分，有对于三角形积分，有积分限中包含了变量自身。1.2.2 确定积分阶的原理1.2.2.1 保证积分的精度以一维问题刚度矩阵的积分为例，如果插值函数N中的多项式阶数为P，微分算子L中导数的阶次是m，则有限元得到的被积函数是2（p-m）次多项式。为保证原积分的精度，应选择高斯积分的阶次n=p-m+1,这时可以精确积分至2（p-m）+1次多项式，可以达到精确积分刚度矩阵的要求。还需要指出，由于位移有限元所根据的最小位能原理，所以当单元尺寸h不断减小时，有限元解将单调地收敛于精确解。对于二维、三维单元刚度矩阵进行精确积分的条件下，将保证当单元尺寸h不断减小时，有限元解单调地收敛于精确解。采用减缩积分往往可以取得较完全精确积分更好的精度。因为精确积分常常是由插值函数中非完全项的最高次所要求，而决定有限元精度的，通常是完全多项式的方次。这些非完全的最高方次项往往并不能提高精度，反而可能带来不好的影响。另外，在最小位能原理基础上建立的位移有限元，其解答具有下限性质。即有限元的计算模型具有较实际结构偏大的整体刚度。选取减缩积分方案将使有限元计算模型的刚度有所降低，因此可能有助于提高计算精度。1.2.2.2 保证结构总刚度矩阵K是非奇异的求解系统方程K a=P，要求方程有解则必须系数矩阵的逆矩阵是存在的，即在引入强迫边界条件后K必须是非奇异的。系数矩阵K非奇异的条件是，或称K是满秩的。在实际计算中，只在采用减缩积分方案计算矩阵K时，才需要检查矩阵K非奇异的必要条件是否得到满足因为采用精确积分方案计算矩阵K时，不仅矩阵K非奇异的必要条件而且它的充分条件都是恒被满足的。如果采用精确积分方案进行矩阵K的计算，实际上就是对给定离散方案下的应变能进行精确的计算，而总是大于零的，所以K必然是正定的，也即K是非奇异的。如果采用减缩积分方案，情况及不同了，此时K 非奇异性的必要条件不是恒被满足的。1.3 线性等参单元和非协调元1.3.1 线性等参单元对于二维、三维单元刚度矩阵的数值计算，则需要对被积函数作进步的分析。例如二维4结点双线形单元，它的插值函数中包含项，在假设单元的是常数的情况下，刚皮矩阵的被积函数中包含项。出于被积函数在和方向的最高方次为2，所以要达到精确积分，应采用阶高斯积分，如果单元的常数，则需要选取更多的积分点。对于二维8结点单元也可作类似的分析、结论是：为精确积分单元刚度矩阵，在常数条件下，应采用阶高斯积分。如果常数，则需要采用更高阶的高斯积分。正如前面已指出的，在对单元刚度短阵进行精确积分的条件下，将保证当单元尺寸h小断减小时，有限元解单调地收敛于精确解。 但是在很多情况下，实际选取的高斯积分点数低于精确积分的要求。例如按单元插值函数中完全多项式的阶数产来选取，仍以上述二维4结点和8结点单元为例，它们的插值函数中完全多项式阶数p分别等于l和2。保证这部分被积函数积分的精度，只需要分别采用和的高斯积分。这种高斯积分阶数低于被积函数所有项次精确积分所需要阶数的积分方案称之为减缩积分。实际汁算表明；采用减缩积分往往可以取得较完全精确积分更好的精度。1.3.2 非协调元E. Wilson提出了二维和三维非协调元，改善单元的性质，提高精度。在二维线性单元的位移模式中，增加与节点位移无关的位移附加项：这些附加项与原位移模式共同构成完全二次式，改善了单元的模拟能力，但引起边界位移的不协调，称为非协调项。位移模式含非协调项的单元称为非协调元。非协调元破坏了收敛条件的协调性要求。然而，可以证明，对型问题，当单元尺寸不断缩小（常应变）的极限情况下，如果单元之间的位移协调性能得到恢复，则非协调元的解仍然收敛。平面4节点和三维8节点非协调元在计算中可以达到与平面8节点和三维20节点完全二次位移模式协调元同级的计算精度，而计算量却大大减小。然而，非协调元对于网格扭曲很敏感。应用时，尽量采用规则的单元形状。实体单元应用时可选择线性位移模式或二次位移模式，完全积分或减缩积分。二维四边形和三维六面体线性单元有非协调位移模式可选择。结构分析时应尽量避免使用线性位移模式的平面三角形单元和三维四面体单元。除非网格密度非常大，对于一般实际问题尤其是有明显应力梯度的问题这类单元求解精度很差。只有在特殊场合采用。一般应力分析，优先使用8节点四边形二次单元（减缩积分、完全积分），其次使用四节点非协调元、6节点三角形二次单元。但非协调元要避免扭曲的单元形状。在网格扭曲严重的情况下，优先使用网格细化的线性、减缩积分单元。有弯曲变形情况下，避免使用4节点双线性位移模式完全积分单元。因为该单元有剪力自锁（shear locking),会引起单元弯曲时过于刚硬。对于接触问题，采用细网格的线性、减缩积分单元或者非协调模式单元。三维问题应用中的单元选择:对于中小规模的一般结构分析，优先采用20节点六面体减缩积分、完全积分单元。应力集中区域，应采用完全积分单元。对大规模分析或接触问题，优先采用8节点六面体非协调模式单元。网格扭曲严重的情况下，应使用细网格8节点六面体减缩积分单元。要求快速建模情况下，可考虑采用10节点四面体二次单元。但得到的模型规模（节点数）大于等效网格的六面体单元，否则结果不精确。零能模式的存在将使解答失真，甚至使求解无法进行。因此在实际分析中，必须防止零能模式的出现。亦即在采用减缩积分方案时，必须注意检查K的非奇异性条件是否得到满足。1.4 等参单元的应用与结论等参单元形状、方位任意，容易构造高阶单元，适应性好，精度高。其优点是有较大的选择单元的自由、能很好地模拟曲线边界，计算精度高，输人数据少，这一点对复杂区域的求解时特别突出等参单元列式具有统一的形式，规律性强，采用数值积分计算，程序处理方便。由于等参单元涉及单元几何形状的变换，对实际单元的形态有一定要求。单元形态好坏影响计算结果的精度。单元形态应满足：单元各方向的尺寸尽量接近；单元边界不能过于曲折，不能有拐点和折点，尽量接近直线或抛物线；边之间夹角接近直角。构造等参数单元是以局部坐标为出发点，并使整个讨论立足于局部坐标上。即在局部坐标系中的规则单元上构造相容的插值函数导出相容的坐标变换式(得到整体坐标系中单元形状)由局部到整体坐标系的坐标变换形成有限元计算格式(其中用到数值积分)。整个讨论和计算都是在局部坐标系中规则单元内进行的。等参单元在有限元法的发展中占有重要的位置，由于它能使局部坐标系内的形状规则的单元变换为总体坐标系内形状为扭曲的单元，从而为求解域是任意形状的实际问题的求解提供了有效的单元形式。要注意等参单元在实际应用中的环节，首先是由于等参单元的特性矩阵是建立在单元局部坐标系的，因此必须进行导数、体积、面积、长度等的变换。等参单元中重要的问题是数值积分方法和阶次的选择，虽然由于等参元的特性矩阵是建立于局部坐标系的规则域内，为数值积分方法的采用提供了很大的方便，同时高斯积分也已被证实是最方便而有效地方法而被广泛采用，但积分阶次的选择仍是等参单元应用中需要认真对待的问题。2 分析与计算2.1 2个节点一维杆求解直杆的位移、应力解：采用有限元法求解根据图示结构，应用里兹法分两个区域假设位移场： 代入势能泛函： 对单元插值得；。对单元插值也有同样的位移场，同上。整个杆上，由各单元上假设的位移场拼接而成的位移试探函数是连续的，只要我们记住，得到的就是全域可能位移场。这样的位移场已经把节点位移自由度作为广义坐标。在两个单元上分片进行总势能计算：首先计算应变：单元内有：， 所以一个单元内的应变能为：或者，则载荷在两个单元内分别表达为：， 带入外力功积分式，对各个单元分别计算外力功。单元外力功为：单元外力功为： 总势能是两个单元总应变能减去两个单元外力功。为了能够以矩阵形式相加，将单元势能矩阵表达式中的单元节点位移列阵用结构整体位移向量代替，单元刚度矩阵等与单元有关的矩（列）阵扩展成结构整体规模（55），则相加后系统总势能表达为：上述表达形式上的变换，不改变总势能的大小！上式简写为：应用驻值条件：，得到节点平衡方程，即：考虑到，并划去第一个方程，解出其余四个方程得到：由此得到的位移场在一般位置上均为近似值（小于精确解）。单元应力由公式： 得到： 根据材料力学知识可求得精确解如下：位移和应力的计算结果与精确解的比较如下图所示：应力值与精确解对比位移值与精确解对比有限元解 由看出，位移的近似度比应力更好。2.2 证明平面三节点三角形单元发生位移时，不产生应力证明：平面问题三节点三角形单元发生刚体位移时有两种情况: 刚体的平动;刚体的转动。三节点三角形单元的应力为：当刚体发生平动时因此：同理：故，当刚体发生平动时单元中不产生应力。当刚体发生转动时，假设刚体绕点的转过的角度为，三角形单元的边长分别为，则：； ；由应力应变公式可得：将代入上式可得：当刚体转动时也不产生应力由综合可得，平面问题三节点三角形单元发生刚体位移时，单元中将不产生应力。2.3证明对四节点矩形单元，存在剪力自锁和不存在的情况解：局部节点编号1，2，3，4。在单元中心（x0,y0）建立一个局部坐标系-。坐标轴平行于矩形的两边。则，与x，y之间有简单的变换关系：x=x0+ay=y0+b矩形单元,在单元4个节点上的值分别为1，因此看作自然坐标。其位移模式：其中为插值函数（形函数）。具体表达式为：每个节点的自然坐标用符号i(i,i)（i=1,2,3,4）表示,则上述形函数可写成通式： 利用双一次线性拉格朗日多项式构造插值函数可得：，。 ,单元应变矩阵：，式中：应力矩阵：，其中，其中给单元施加如图所示对称载荷，单元内将产生纯弯曲。由于在单元变形变形过程中，单元边界不能弯曲，直边只能变成直边。可设,则：于是得到；而精确解中，可见产生了虚加的剪应力，存在剪力自锁。引入Wilson非协调项后， ，解得，如果取适合的值，当，则，因此剪力自锁可消除。3 上机实验报告3.1 建立减缩积分模型，网格密度、位移约束对奇异性、求解精度的影响实验报告3.1.1 实验题目采用线性减缩积分等参元的有限元模型，通过数值试验验证网格密度、位移约束条件对系统奇异性、求解精度的影响。3.1.2 实验目的通过数值试验验证网格密度、位移约束条件对系统奇异性、求解精度的影响。3.1.3 建模概述薄板长宽均为20mm,正左下角上挖掉半径为5mm的圆。右端面受50Mpa的外压力，弹性模量E=200Gpa,泊松比为0.3。切出中间典型的一片可按平面应变问题进行计算，由于对称，取其1/4圆筒来建模。采用矩形单元，其网格、约束、载荷和位移如下图细网格，双约束 应力云图 细网格，单约束 应力云图 粗网格，单约束 应力云图 粗网格，双约束 应力云图3.1.4 计算结果分析与结论由材料力学知识可得，精确解是150MPa左右。比较应力云图可知，单约束时，细网格的应力最大值147MPa, 粗网格的应力最大值172MPa，故细网格的精度比粗网格的更高。随着网格的加密，线性缩减积分所计算的结果也就越精确。最大应力是随着网格的加密而逐渐由小到大而接近精确解，而最小应力则相反。当采用不同的约束时，对有限元的分析结果影响很大。如果约束条件较少或约束较多时，则会发生奇异现象，不能得到正确的结果。网格细化对提高求解精度非常重要，网格越细，求解精度越高。采用减缩积分方案，此时K非奇异性的必要条件不是恒满足的，3.1.5 实验体会与总结个人觉得，刚开始看书时，有些郁闷,比较迷茫,通过慢慢的操作ANSYS软件，熟悉了一些命令，看了一些实例子，做题过程中遇到了很多困难与挫折,多次向身边的朋友请教了一下，通过操作本题目，对ANSYS软件和有限元法有了更进一步的了解。用ANSYS软件分析有限元时，首先要画好网格，这是一个很关键的环节，模型的好坏主要取决于网格的划分，同时还要注意选取单元的约束条件，约束要适当选取，不能过多或过少。操作软件时一定要细心，不能马马虎虎，刚开始操作时出现了很多低级错误，通过自己看书和请教同学，这些问题都很好的解决了。3.2 管接头静态响应分析报告3.2.1 实验题目 管接头静态响应分析报告3.2.2 实验目的利用20节点六面体单元对其建模并求解3.2.3 建模概述管接头内壁受均匀压力作用，此管接头为对称结构，分析时可利用对称性，从而可以减少计算分析的工作量。管接头为钢材，材料的弹性模量为，泊松比0.3。管接头内部具有1个大气压的压力载荷。管接头三维模型如下图2-1所示：图2-1 管接头三维模型计算中管接头采用简支约束，左端约束、三个平动自由度，右端约束、两个平动自由度，对称面上约束方向平动自由度,载荷为压力载荷，约束与载荷与下图2-2,2-3所示。 图2-2 20节点六面体单元约束和网格 图2-3载荷 计算结束后的图3.2.4 计算结果分析与结论由计算分析得，管接头应力云图如下图2-5所示，最大应力值位置在角点处，该处由于结构突变而产生应力集中。由于两端面平动自由度被约束，故应力值较小。中间部位受压力作用，结构向外膨胀，符合实际应力情况。图2-5 应力云图位移云图下图2-6所示图2-6位移云图3.2.5 实验体会与总结管道三通，这是机械结构应力分析中典型的三维结构，通常利用对称性，取结构的1/4建立有限元模型，进行各种载荷作用下的结构应力分析。在此试验过程中，熟练掌握了用建模软件构建三通管，体会到用云图的直观性，并学会在云图上观察每一点的应力和位移值。进一步提高了运用ANSYS软件的能力。通过本题目，更加熟悉了ANSYS 软件的操作，操作技巧更流畅了。3.3 矩形平板动态响应分析报告3.3.1 实验题目矩形平板动态响应分析3.3.2 实验目的用粗网格的线性和二次六面体全积分等参单元计算其前六阶自由振动频率和振型；用隐式直接积分和振型叠加法对该平板进行瞬态响应分析。3.3.3 自由振动频率和振型的建模概述3.3.3.1 自由振动频率与振型的计算用线性和二次六面体的等参单元分析计算前六阶的自由振动：1） 先建立800*100*5mm矩形平板。2） 定义单元类型：线性单元8节点的六面体单元选择Brick-8node-45；二次单元20结点的六面体单元选择Brick-20node-95。3） 添加材料定义：取弹性模量EX=，泊松比，密度。4） 划分网格：对平板进行网格划分。5） 分析类型选择模态分析，并在板的两个侧边施加全约束。6） 开始计算结果并输出，计算完成后可以输出各阶自振频率和各阶振型图。本计算中计算两种不同网格单元下前六阶自由振动频率和振型，并对比结果。网格为线性和二次六面体全积分等参元。矩形长板三维实体网格如下图3-1和图3-2所示。 图3-1线性网格 图3-2二次六面体网格3.3.4 自由振动频率和振型的计算结果分析与结论3.3.4.1 线性与二次六面体计算结果比较线性与二次六面体计算出的结构各阶振型对比如下，图中左边为线性单元计算结果，右边为二次单元计算结果。 一阶振型对比如下图3-3 图3-3 一阶振型对比 二阶振型对比如下图3-4 图3-4 二阶振型对比 三阶振型对比如下图3-5 图3-5 三阶振型对比 四阶振型对比如下图3-6 图3-6 四阶振型对比 五阶振型对比如下图3-7 图3-7 五阶振型对比 六阶振型对比如下图3-8 图3-8 六阶振型对比线性和二次单元各阶频率值如下表1，表中单位。由表中数据可以看出二次单元计算出的各阶固有振型和固有频率均较线性单元大，虽然二次单元较一次单元计算出的结果更准确，但由于实际结构均含有阻尼，有阻尼的固有频率和无阻尼的固有频率的关系是，而结构频率分析中为考虑阻尼，所以在实际使用中可能一阶单元计算出的值更符合实际。 * INDEX OF DATA SETS ON RESULTS FILE * SET TIME/FREQ LOAD STEP SUBSTEP CUMULATIVE 1 41.551 1 1 1 2 115.03 1 2 2 3 207.93 1 3 3 4 227.23 1 4 4 5 379.64 1 5 5 6 422.20 1 6 6 线性自由振动的各阶频率 * INDEX OF DATA SETS ON RESULTS FILE * SET TIME/FREQ LOAD STEP SUBSTEP CUMULATIVE 1 41.399 1 1 1 2 114.19 1 2 2 3 206.82 1 3 3 4 224.37 1 4 4 5 372.14 1 5 5 6 419.23 1 6 6 二次自由振动的各阶频率表1 线性与二次六面体单元各阶频率值二次线性一阶41.39941.551二阶114.19115.03三阶206.82207.93四阶224.37227.23五阶372.14379.64六阶419.23422.20 从线性和二次六面体的自由振动的各阶频率的数值来看，各阶相对应的频率相差很小，只有略微的变化，二次等参单元计算的自由振动频率要比线性等参单元的自振频率略小，更接近理论值。3.3.5 隐式积分分析动态响应建模概述建立800*100*5mm矩形平板。定义单元类型 添加材料定义：取弹性模量EX=，泊松比，密度。划分网格施加边界约束：两边简支约束。施加载荷时间分别设为0.2s、0.4s、0.8、1.2s的时间进行计算。首先要先设置一半的时间的加载三角形脉冲均布载荷，最大值为0.05MPa，这个区间设为5步计算就可以了，记为第一步；然后再设置后一半时间的卸载斜坡载荷到0，这个区间也设为5步计算，记为第二步，最后开始分析计算。平板模型网格3.3.6 隐式积分分析动态响应的计算结果分析与结论一个脉宽的情况下，都是在载荷最大的时刻同时出现应力的最大值，因此，着重分析在载荷最大的时刻的应力云图。图3.17 脉宽为0.2s的0.1s时的应力云图脉宽为0.2s的0.1s时，板中部的动态应力响应最大为9137Pa。图3.18 脉宽为0.4s的0.2s时的应力云图脉宽为0.4s