有限元期末大作业答案

1、第39 页 共 39 页. 一 等参单元及其应用11 概述11.1 等参单元的概念11.2 等参单元的原理11.3 工程应用的意义22 等参单元的数值积分方法32.1 数值积分方法32.2 确定积分阶的原理32.3 全积分单元与减缩积分单元讨论43 线性等参单元和非协调元54 等参单元的应用6二 分析与计算61 计算题一62 计算题二73 计算题三85 计算题四9三 上机实验133.1 实验一133.1.1 实验题目133.1.2 实验目的143.1.3 建模概述143.1.4 计算结果分析与结论153.1.5 实验体会与总结233.2实验二243.2.1实验题目243.2.2实验目的243.

2、2.3 建模概述253.2.4 计算结果分析与结论263.2.5 实验体会与总结273.3实验三273.3.1实验题目273.3.2实验目的273.3.3建模概述283.3.4计算结果分析与结论283.3.5实验体会与总结39一、课程论文：等参单元原理及应用1、 等参单元概述1.1概念 在有限元的学习中，我们在书本上经常可以看到三角形单元及四边形单元的应用，其边界都是直线和平面，对于结构复杂的曲边和曲面外形，只能通过减小单元尺寸，增加单元数量进行逐渐进行逼近，这样自由度的数目随之增加，并且使得计算时间长，工作量大。另外这些单元的位移模式是线性模式，是实际位移模式的最低级逼近模式，问题的求解精度

3、收到了很大的限制。从本学期书本的学习中了解到，之前介绍的各种2、3维单元主要受到两个方面的约束：第一是单元的精度，显然单元的节点数越多，单元精度越高。因此在这一点上，矩形单元优于3节点三角形单元，六面体单元优于四面体单元；第二是单元几何上的限制，上述矩形和六面体单元都不能模拟任意形状几何体，所有几种单元都是直线边界，处理曲边界几何体误差较大。解决上述矛盾的出路就是突破矩形单元和六面体单元几何方面的限制，使其成为任意四边形和任意六面体单元，如果再增加边中间节点，还可以成为曲边四边形和曲面六面体高精度实用单元。但是这类单元位移模式和形函数的构造和单元列式的导出不能沿用前面构造简单单元的方法，必须引

4、入所谓的等参变换，采用相同的插值函数对单元的节点坐标和节点位移在单元上进行插值。这种单元称为等参单元。等参单元的提出对于有限元法在工程实践中的应用具有重要意义。1.2原理 下图为一个4节点任意四边形单元，单元有8个自由度。将矩形单元放松为4节点任意四边形单元将带来许多好处。但在建立单元位移模式时产生了新的问题：单元上没有一个如矩形单元中的简单直接的局部坐标系，而又不能直接用x，y坐标系下的双线性位移模式（位移沿边界二次变化，不协调）。 须建立一种新的局部坐标系-（如图），使得4条边的坐标为常数（±1），则在-平面内，单元是一个边长为2的正方形。 同时，该局部坐标系的建立在x-y平面上

5、的任意四边形单元与-平面上的正方形之间形成了一个1-1对应的映射关系。 称-平面内的正方形单元为基本单元或母单元。x-y平面内的任意四边形单元称为实际单元。显然，母单元的节点相应于不同的x， y坐标就得到不同的任意四边形单元。 图1 4节点任意四边形单元及其母单元 建立了局部坐标系或映射后，我们只需要在-平面上的母单元中描述实际单元的位移模式和力学特性。任意四边形单元在母单元中的位移模式（或者称为-坐标系下的位移模式）与矩形单元相同： 当然，该位移模式在x，y坐标系下不是双线性位移模式，位移沿单元边界线性变化，能保证单元的协调性。 为了得到上述映射的数学表达，引入对母单元节点上x，y坐标进行插

6、值的思想，将母单元上每一点对应的x，y坐标看成是对节点坐标的插值，插值函数与位移插值中的形函数相同： 这样就得到了一个事实上的映射，该映射是用母单元描述实际单元力学特性的桥梁。由于该几何变换式中采用了与位移模式相同的参数（插值函数），因此称为等参变换。而所有采用等参变换的单元都称为等参单元。1.3工程应用的意义 由于等参变换的采用使等参单元的刚度、质量、阻尼、载荷等特性矩阵的计算仍在原单元的规则域内进行，因此不管各个积分式的被积函数如何复杂，都可以方便地采用标准化的数值积分方法进行计算，从而使各类不同工程实际问题的有限元分析纳入了统一的通用化的程序。借助于等参元可以对一般的任意几何形状的工程问

7、题和物理问题方便地进行有限元离散。因此，等参元的提出为有限元法成为现代工程实际领域最有效的数值分析方法迈出了重要的一步。2、等参单元数值积分方法2.1、数值积分方法 等参单元刚度矩阵的被积函数形式往往比较复杂，一般都不能进行显示积分而需求助于数值积分。计算单元特性矩阵的方法有Newton-Cotes积分、高斯积分、Irons积分、Hammer积分等。为了减少积分点的数目和便于程序编制，一般采用高斯数值积分方法。 高斯积分：一维高斯数值积分公式： 积分点 积分点数目，积分阶 权重系数结论：n阶高斯积分公式对2n-1次多项式被积函数可求得精确积分！ 同理，对二维高斯积分： 积分公式对，方向最高方次

8、为2n-1的多项式可求得精确值。2.2、确定积分阶的原理 （1）保证刚度矩阵积分精度的积分阶选择 完全积分方案选择积分阶的基本考虑是保证被积函数所有项次精确积分，这种积分方案称为完全积分方案。以二维单元刚度矩阵高斯数值积分为例：单元刚度矩阵为根据这个原理，具有规则形状的单元（|J|=常数）完全积分方案如下：二维4节点、三维8节点等参元分别是2×2、2×2×2积分；二维8节点、三维20节点等参元分别是3×3、3×3×3积分。对于|J| 常数，需要增加积分阶数。2）减缩积分方案实际应用中选取的积分阶往往可以低于被积函数所有项次精确积分所需

9、要的阶数，这种积分方案称为减缩积分。对二、三维连续体单元，通常按形函数中完全多项式阶数 p 来选取积分阶，即取 n=p。根据这个原理，具有规则形状的单元（|J|=常数）减缩积分方案如下：二维4节点、三维8节点等参元分别是1×1、1×1×1积分；二维8节点、三维20节点等参元分别是2×2、2×2×2积分。实际计算表明，采用减缩积分往往可以取得较完全精确积分更好的精度。（2）保证结构总刚度矩阵非奇异的积分阶选择 数值积分条件下刚度矩阵的秩单元刚度矩阵数值计算公式为：根据矩阵秩的基本规则有，是高斯积分点数秩 ，为保证正确求解，数值积分的秩不

10、能小于，即，是系统的独立自由度数，也就是系数矩阵的阶数。因此系数矩阵非奇异的必要条件是。 减缩积分情况下刚度矩阵的奇异性减缩积分得到的往往缺秩，即是奇异的，即存在单元发生非刚体位移模式情况下，单元应变能为零的情况零能模式（沙漏模式）。因此要注意检查奇异性的必要条件，即。2.3、全积分单元与减缩积分单元讨论 虽然减缩积分会出现零能模式，但在实际计算中减缩积分应用的十分广泛。实际计算表明减缩积分往往可以取得较完全精确积分更好的精度。这是主要由于：(1) 精确积分常常是由插值函数中非完全项的最高方次所要求，而决定有限元精度的，通常是完全多项式的方次。这些非完全的最高方次项往往并不能提高精度，反而可能

11、带来不好的影响。取较低阶的高斯积分，使积分精度正好保证完全多项式方次的要求，而不包括更高次的非完全多项式的要求，其实质是相当于用一种新的插值函数替代原来的插值函数，从而一定情况下改善了单元的精度。(2) 在最小位能原理基础上建立的位移有限元，其解答具有下限性质。即有限元的计算模型具有较实际结构偏大的整体刚度。选取减缩积分方案将使有限元计算模型的刚度有所降低，因此可能有助于提高计算精度。(3) 减缩积分方案对于泛函中包含罚函数的情况也常常是必须的，用以保证和罚函数相应的矩阵的奇异性，否则将可能导致完全歪曲了的结果。由于减缩积分单元的特性，因此在使用时须注意一些问题，在网格扭曲严重的情况下，优先使

12、用网格细化的线性减缩积分单元，对于接触问题，采用细网格的线性减缩积分单元，对于规则单元等选用完全积分单元。3、线性等参单元在ANSYS软件中采用等截面悬臂梁端部受弯矩，求端部位移的平面应力问题来说明全积分、减缩积分线性等参单元在计算精度、剪力自锁、零能模式与总刚的奇异性等方面的差异。悬臂梁左端固定，右端受一大小为M6 N·m的弯矩。梁的长度l=0.3m,横截面高度h=0.03m，宽度b=0.01m。材料为钢，弹性模量E=210GPa，泊松比u=0.3，速度的理论解：（1）采用ANSYS14.0中的PLANE182完全积分双线性单元，计算结果如下：单元1×102×2

13、04×208×40PLANE1820.386e-40.521e-40.520e-40.561e-4（2）采用ANSYS14.0中的PLANE182减缩积分单元，计算结果如下：单元1×102×204×208×40PLANE1820.1590.756e-40.614e-40.588e-4从计算结果发现：虽然与理论计算结果有差异，但是在网格划分密度相当的情况下，完全积分的计算结果较为精确 这是由于完全积分的线性单元出现了剪力自锁，剪力自锁引起单元在弯曲时过于刚硬，造成应变偏小，同时随网格的加密，剪力自锁效应可以得到改善，但与实际结果偏差仍较

14、大。而减缩积分单元存在沙漏问题（零能模式），使得结构过于柔软，变形偏大，甚至在在粗划网格情况下，产生无意义的结果，但在加密网格后，结果能够改善，但依然偏大。非协调元在一次单元中引入了非节点的附加自由度，克服了剪力自锁和零能模式的不良影响，可以产生与二次单元相当的结果，而计算量却明显的降低，但其对于单元的过度变形很敏感，可能影响结果，所以在网格划分上，对于使用者有着较高的要求。4、等参单元的应用 等参单元通常用来模拟二维、三维实体单元。应用时可选择线性位移模式或二次位移模式，完全积分或减缩积分。二维四边形和三维六面体线性单元有非协调位移模式可选择。 平面问题中的单元使用一般应力分析，优先使用8节

15、点四边形二次单元（减缩积分、完全积分），其次使用四节点非协调元、6节点三角形二次单元，但非协调元要避免扭曲的单元形状。在网格扭曲严重的情况下，优先使用网格细化的线性减缩积分单元。有弯曲变形情况下，避免使用4节点双线性位移模式完全积分单元。因为该单元有剪力自锁，会引起单元弯曲时过于刚硬。对于接触问题，采用细网格的线性减缩积分单元或者非协调模式单元。 三维问题中的单元选择对于中小规模的一般结构分析，优先采用20节点六面体减缩积分、完全积分单元，应力集中区域，应采用完全积分单元。对大规模分析或接触问题，优先采用8节点六面体非协调模式单元。网格扭曲严重的情况下，应使用细网格8节点六面体减缩积分单元。要

16、求快速建模情况下，可考虑采用10节点四面体二次单元，但得到的模型规模（节点数）大于等效网格的六面体单元，否则结果不精确。二、分析与计算1、证明平面4节点等参元满足收敛的协调性准则。（10分）证明：任意四边形单元在母单元中的位移模式（或者称为-坐标系下的位移模式）与矩形单元相同： （1） （i=1,2,3,4） 该位移模式在x，y坐标系下不是双线性位移模式，位移沿单元边界线性变化，能保证单元的协调性。在实际单元中取四个节点的坐标（）为参数，用上式的形函数构造插值函数，建立用坐标变量确定实际单元中任一点坐标（x,y）的关系式： （2） 这样可以将局部坐标系中矩形单元的节点变为整体坐标系中任意四边形

17、单元的对应节点。位移模式（1）在单元边界上线性变化，并且根据单元公共边界上两个节点位移插值得到，因此单元的协调性得到满足。同理，坐标变换式（2）也能保证变换后的实际图形仍然保持单元公共边界上的坐标值处处相同，满足收敛的协调性准则。 2、计算图示平面等参单元在2-6-3边作用均布水平载荷q时的等效节点力。单元厚度为t。（10分）由题意可知，建立如图2.1所示直角坐标系图2.1直角坐标下的平面等参单元由于只有2-6-3边受均布载荷作用，故只有2、6、3这三个节点有等效节点力，由8节点平面等参单元的形函数可得： 由x-y坐标系与坐标系重合，O点为两坐标系的原点，则有如下的坐标变换又因为2-6-3边对

18、应=-1，进一步化简得作用在边界上的面积力为T=q，单元等效节点力计算公式为：=,其中 N= 将N和T代入上式得3、证明平面问题三节点三角形单元发生刚体位移（小位移平动和转动）时，单元中将不产生应力。（10分）证明：设3节点单元在x方向有刚性位移，y方向有刚性位移，则有位移列阵为 应变系数行列式 得 所以 可得 ，同理 应变矩阵 应力矩阵 ，即单元中不产生应力。5、如图示一根直杆，长度2L，截面积A，弹性模量E。杆受到沿轴向的线分布力：q=cx。试用2个3节点一维杆单元求解其位移、应力。要求推导详细的有限元求解列式，设置合理的参数值将求解结果绘制成曲线，并与精确解进行对比分析。（10分）解：将

19、直杆分成两个单元，在每个单元内假设简单的位移场。单元1中根据单元内3个节点上的位移未知量进行插值得到单元上任意点的位移，设 ，由此可得当 时，当 时，当 时，联立方程解得 对单元2也同样可得到插值函数形式的假定线性位移场。单元2： 在两个单元分别进行总势能计算有单元内有 计算可得 所以一个单元内的应变能为：其中 单元2的与单元1相同。两个单元杆上的载荷分别为：代入外力功积分式，对两个单元计算外力功。单元1： 单元2： 位移 系统总势能表达为 应用驻值条件：，得到节点平衡方程：即可得 考虑到 ，并划去第一个方程，解出其余四个方程得到：，对于一维杆的精确解为：位移精确解： 应力的精确解：用MATL

20、AB绘制位移和应力曲线图如下： 图 5.1 从上图可以看出，单元1,2的理论解与精确解几乎重合，说明3节点的杆单元的解有很高的精度。 图 5.2 从上图可以看出，单元1,2的理论解逼近精确解，说明3节点杆单元解应力也很准确。综上所述得到的位移场在一般位置上均为近似解。单元应力由公式：得到。三、上机实验3.1实验一3.1.1实验题目1、图示一个简支梁平面应力模型。梁截面为矩形，高度h=160mm，长度L=1000mm，厚度t=10mm。上边承受均布压力q =1N/mm2，材料E=206GPa，=0.29。X方向正应力弹性力学理论解为：分别应用3节点三角形单元、4节点线性等参元（完全积分、减缩积分

21、、非协调模式）、8节点二次等参元进行下列数值实验：1）用较粗单元网格求解梁中部应力分量的最大值和上下边法向应力分量，并对计算精度进行比较；2）用密网格进行上述计算并比较计算精度；3）对粗网格下梁中部最大位移进行对比和分析。总结出研究结论，撰写实验报告。3.1.2实验目的（1）掌握平面应力问题下的有限元解法，学习如何使用ANSYS软件对三节点三角形单元和四节点线性等参单元（完全积分、减缩积分、非协调模式）及八节点二次等参单元进行数值实验。并会使用软件进行有限元的几何建模、网格划分、载荷施加、运算求解以及后处理。通过实验比较单元网格粗细和算法的选择对求解精度的影响；（2）通过实验比较3节点三角形单

22、元、4节点线性等参元和理论解的差异，学习网格密度对系统求解精度的影响，考察有限元解的收敛性。3.1.3建模概述（1） 选用单元类型对于3节点三角形单元和4节点线性等参元选用solid182单元，并对solid182单元的算法进行相应设置。对于8节点二次等参元则采用solid182单元对应的高阶单元solid183单元。（2） 设定材料属性 创建材料，创建截面属性，并给部件赋予截面属性设定材料参数E=206GPa，=0.29。（3） 建模由于问题具有对称性，建模时只需以其中的一个平面来建模。建立长度为1m，高度为0.16m的平面四边形区域。（4） 划分网格为单元分配属性并划分成单元。（5） 施加

23、载荷选择图元左侧中点并限制全部自由度以及右侧中点Y方向的自由度。选择图元上表面并施加压强q=106N/m。（6） 提交计算（7） 查看结果3.1.4计算结果分析与结论1）理论解X方向正应力由下式计算：已知q=1N/mm2,h=160mm,L=1000mm，代入上式得2） 粗网格（）情况下应力分析结果如下：图1粗网格三角形单元 图2粗网格四边形单元图3三角形单元应力云图图4三角形单元上边法向应力 图5三角形单元下边法向图6四边形全积分单元应力云图 图7四边形全积分单元上边法向应力 图8四边形全积分单元下边法向应力 图9四边形减缩积分单元应力云图 图10四边形减缩积分单元上边法向应力 图11四边形

24、减缩积分单元下边法向应力图12非协调元应力云图 图13非协调元上边法向应力 图14非协调元下边法向应力图15二次等参元应力云图 图16二次等参元上边法向应力 图17二次等参元下边法向应力现将粗网格下梁中部应力分量x的最大最小值列于下表：单元类型三角形单元全积分单元减缩积分非协调元二次等参元xmax(MPa)12.725.519.128.729.6xmin(MPa)-12.6-25.5-19.1-28.7-29.6从计算结果发现：在粗网格划分情况下，二次等参元计算精度最高，其次是非协调元，其它单元计算误差均非常大，尤其是减缩积分单元，产生了严重的沙漏，影响了计算的精确性。粗网格下梁中部最大位移如

25、下表所示单元类型三角形单元全积分单元减缩积分非协调元二次等参元最大位移0.131mm0.165mm1.055mm0.193mm0.197mm由表可知二次等参元和非协调元的计算结果相当，而减缩积分单元存在严重的沙漏，变形严重失真，完全不能反映实际情况。其次还可发现全积分单元计算精度比三角形单元高。通过网格加密后，在密网格（）的情况下的应力分析结果如下图所示，密网格三角形单元：图1密网格三角形单元图 2密网格三角形单元图3三角形单元应力云图 图4三角形单元上边法向应力 图5三角形单元下边法向应力图6四边形全积分单元应力云图 图7四边形全积分单元上边法向应力 图8四边形全积分单元下边法向应力图9四边

26、形减缩积分单元应力云图 图10四边形减缩积分单元上边法向应力 图11四边形减缩积分单元下边法向应力图12非协调元应力云图 图13非协调元上边法向应力 图14非协调元下边法向应力图15二次等参元应力云图 图16二次等参元上边法向应力 图17二次等参元下边法向应力现将密网格下梁中部应力分量x的最大最小值列于下表：单元类型三角形单元全积分单元减缩积分非协调元二次等参元xmax(MPa)27.829.726.129.529.5xmin(MPa)-27.8-29.7-26.1-29.5-29.5从计算结果发现：在细网格划分情况下，各单元均取得了较理想的计算结果，尤其是二次等参元和非协调元，两者取得了和理

27、论解相一致的解。3.1.5实验体会与总结通过此次实验，并且结合了书本上的知识以及仿真软件的运用，使得我对三角形单元、四边形完全积分单元、四边形减缩积分单元、非协调元和二次等参元有了深刻认识，尤其是减缩积分单元。减缩积分单元在计算时会产生沙漏，因此对这类单元应使用足够密的单元网格。另外还可以看出在结构分析计算中进行数值积分时，积分单元、积分方法、积分阶次、网格密度的不同选择会导致计算结果的差异。如果单元算法选择不当，除影响到计算精度外，还有可能导致计算的失败。3.2实验二3.2.1 实验题目图示一管接头，内壁受均匀压力。自行建立其三维几何模型，运用二次六面体单元对其建模并求解。要求利用对称性，自

28、行设置载荷大小和位移约束条件,并撰写实验报告。（10分）3.2.2 实验目的（1）掌握平面应力问题下的有限元解法，学习如何使用软件对二次六面体等参单元进行数值实验。并会使用软件进行有限元的几何建模、网格划分、载荷施加、运算求解以及后处理。3.2.3 建模概述 （1）首先建立管接头实体模型，建模从Solidworks建好，尺寸如下：大管道外径D1=80mm，内径d1=70mm,小管道外径D2=40mm,内径=30mm,长度为100mm，且管接头参数如下：内壁承受均布压力q=20N/mm2,材料E=200Gpa, =0.3,模型如下： 图 1 管接头实体模型（2）设置有限单元类型。选择类型为C3D20。（3）定义材料属性。 （4）进行网格划分。 图 2 网格划分图3.2.4 计算结果与分析 进行仿真之后所得的应力云图如下图： 图 3 管接头应力云图3.2.5 实验总结