90-网格划分密度与有限元求解精度研究

第三届中国CAE工程分析技术年会论文集497网格划分密度与有限元求解精度研究高素荷姚河省（太原重型机械（集团）公司技术中心太原科技大学应用科学学院030024)摘要摘要摘要摘要：：：：本文通过对不同网格密度、不同单元类型的有限元力学模型计算结果与精确解的分析比较，探索研究单元网格划分与有限元求解精度的内在联系，为在保证有限元解满足工程实际精度要求的前提下，确定合理的网格密度，提高有限元分析效率进行了有益的探索。关键词关键词关键词关键词：：：：有限元网格划分密度求解精度有限单元法的基本思想是把一个连续体人为地分割成有限个单元，对通过节点连接的单元进行单元分析，然后再把这些单元组合起来代表原来的结构。从数学的角度来看，有限元法是将一个偏微分方程化成一个代数方程组，并利用计算机求解的一种数值分析方法。它的分析过程可以分为建立力学模型（前处理）、计算及后处理三个阶段。其中，根据结构实际形状和实际工况条件建立有限元分析的计算模型，为有限元数值计算提供必要的原始输入数据，是整个有限分析过程的关键。由于工程结构形状和工况条件的复杂性，要建立一个符合实际的有限元模型不仅需要考虑多种因素，而且输入数据的误差也将直接决定计算结果的精度。所以，其力学模型的正确性和求解精度就成为衡量有限元分析结果精确与否的重要指标。对于有限元这样一种数值分析方法，在单元形状确定之后，当单元网格划分越来越细时，位移近似解将收敛于精确解。增加网格数量和密度，计算精度一般也会随之提高。但是，如果盲目地增加网格数量，将会大大增加单元网格划分时间及求解方程时间。有时还会因计算的累积误差反而会降低计算精度。所以，在实际工作中，如何划分网格才能既保证计算结果有较高的精度，又不致使计算量太大，一直困扰着许多分析人员。本文将通过对不同网格密度、单元类型的分析比较，确定出合理的网格密度，期望能为提高有限元求解精度提供参考依据。1111有限元网格划分的指导思想有限元网格划分的指导思想有限元网格划分的指导思想有限元网格划分的指导思想众所周知,建立有限元力学模型的中心任务是进行单元网格划分和处理许多与之相关的工作，如结构形式处理、力学模型建立、单元特性定义、单元质量检查以及模型边界条件的定义等。有限元网格生成就是将实体模型离散成由节点和单元组成的有限元力学模型。它的指导思想是进行总体模型规划，包括物理模型的构造、单元类型的选择、网格密度的确定等多方面的内容。有限元分析的精度和效率均与单元的密度和几何形状有着密切的关系。在网格剖分时，单元应满足以下要求：(i)一个单元的节点不能落入其他单元内部，单元边界上的节点均应作为单元的节点。(ⅱ)单元必须落在实体模型内部，不能落入外部，且单元集合边界应逼近实体模型的边界。(ⅲ)单元应具有良好的形状，如正多边形或正多面体。（ⅳ）单元之间过渡相对平稳。由此可见，网格划分是建立有限元模型的重要环节，有限元模型的合理性很大程度上可以通过所划分的网格形式反映出来。目前，有限元网格剖分方法可分为拓扑分解法、结点连元法、网格模板法、映射法和几何分解法五种。在通用有限元分析软件中，广泛采用自动或半自动网格划分方法。但由于结构和网格生成过程的复杂性，划分出来的网格有时存在一些问题，如网格形状较差，单元和节点编号顺序不合理等，这些都将影响有限元分析的计算精度和计算时间，有时还需进行人工改进。第三届中国CAE工程分析技术年会论文集4982单元单元单元单元网格网格网格网格划分划分划分划分对求解精度的影响对求解精度的影响对求解精度的影响对求解精度的影响分析分析分析分析有限元解的误差主要来自离散误差、插值误差（即逼近误差和边界误差）以及数值误差几个方面。其中，逼近误差指的是用有限尺寸的单元及单元插值函数代替精确解后产生的误差；边界误差则是在结构边界以直代曲引起的误差。另外，输入数据不正确，也会导致较大的计算误差。就输入数据而言，单元的质量和数量对求解结果和求解过程影响较大，如果结构单元全部由等边三角形、正方形、正四面体、立方六面体等单元构成，则求解精度可接近实际值，但这种理想情况在实际工程结构中很难做到。必须根据模型的不同特征，设计不同形状种类的网格，才能保证有限元分析的正确性，有效地改善网格的质量和求解精度。因此，评价一个有限元网格的优劣，应当从选择单元类型、单元网格划分的数量、网格密度（即疏密程度）、单元阶次及网格质量几方面考虑。2.12.12.12.1单元位移模式对有限元分析精度的影响单元位移模式对有限元分析精度的影响单元位移模式对有限元分析精度的影响单元位移模式对有限元分析精度的影响有限元分析软件的精度主要取决于单元。在进行有限元分析时首先要选择合适的单元类型。在广泛应用的位移法中，4节点四面体单元每个节点的多项式位移模式为zayaxaau4321+++=zbybxbbv4321+++=zcycxccw4321+++=形函数为)(611zdycxbaVNiiii+++=其中，i=1,…，4为常数，V为单元体积。10节点四面体单元的多项式位移模式为基于自然坐标的形函数为8节点六面体单元每个节点的位移模式为xyzaxzayzaxyazayaxaau87654321+++++++=xyzbxzbyzbxybzbybxbbv87654321+++++++=xyzcxzcyzcxyczcycxccw87654321+++++++=21029287654321zayaxaxzayzaxyazayaxaau+++++++++=21029287654321zbybxbxzbyzbxybzbybxbbv+++++++++=21029287654321zcycxcxzcyzcxyczcycxccw+++++++++=34102491483173262154,4,44,4,4)4,3,2,1()12(LLNLLNLLNLLNLLNLLNiLLNiii=======-=yx0324z14节点四面体常应变单元yx0324z1610798510节点四面体高阶单元ζξo1234576η88节点六面体等参数单第三届中国CAE工程分析技术年会论文集499形函数为从以上分析可见，高阶单元的曲线或曲面边界能更好地逼近结构的边界曲线或曲面，高次位移函数能更好地逼近结构复杂的位移分布。单元的多项式位移插值函数阶次越高，精度越高，收敛速度也越快。因此，可选用高阶单元以提高计算精度。但是高阶单元节点较多，使用时应权衡计算精度和模型规模两个因素，处理好单元阶次和节点数量的关系。2.2低阶单元网格与高阶单元网格的比较低阶单元网格与高阶单元网格的比较低阶单元网格与高阶单元网格的比较低阶单元网格与高阶单元网格的比较表1、表2给出的是应用Unigraphics程序采用全自动网格划分方法进行圆轴受纯扭转时的有限元计算结果。单元类型为4节点和10节点四面体三维实体单元。从表1可以看出，选用4节点四面体单元时，随着单元边长从55㎜逐渐加密至30㎜，有限元分析计算结果与精确解的相对误差逐渐从7.15%减小到2.51%，求解精度明显提高。但单元网格继续加密时，有限元力学模型的节点和单元数量增长很快，而求解精度增加很少。从提高工作效率的角度看，过密的网格密度，对提高求解精度意义并不大，换来的只是过长的求解时间。这是很不经济的。表1不同密度4节点四面体单元网格时圆轴受纯扭转的有限元计算结果单元网格参数计算分析结果及相对误差单元类型单元大小(㎜)节点总数单元总数位移(㎜)相对误差剪应力(MPa)相对误差553818158031.2807.15%63.722.43%504253181391.2975.92%64.211.68%455453236501.3095.03%67.012.61%406628291171.3164.54%66.772.24%358660387711.3283.70%65.750.68%3012599576741.3442.51%66.281.49%20296431423301.3680.78%64.081.88%18378491838501.3710.55%64.411.38%16495382430051.3780.07%64.021.97%4节点四面体单元151621289347271.3710.06%64.281.57%注：1）圆轴直径D=340㎜，长度L=1705㎜，扭矩Mn=504×106N.㎜2）精确解为最大剪应力τmax=65.31MPa3）应用软件为Unigraphics，全自动网格划分从表2计算结果可知，如果选择10节点四面体高阶单元，则较稀疏的单元网格也可获得较高的计算精度。而且，同等网格密度下，选用10节点四面体高阶单元的有限元计算结果要比选用4节点四面体单元有限元计算结果更精确。这说明10节点四面体高阶单元的收敛速度明显比4节点四面体常应变单元快得多。8,,2,1,,,)1)(1)(1(81000000L====+++=iNiiiiζζζηηηξξξζηξ其中：第三届中国CAE工程分析技术年会论文集500表2不同密度10节点四面体单元网格时圆轴受纯扭转的有限元计算结果单元网格参数计算分析结果及相对误差单元类型单元大小(㎜)节点总数单元总数位移(㎜)相对误差剪应力(MPa)相对误差120760143091.3920.97%66.331.56%1001007658081.3920.97%66.111.23%7017356103631.3910.90%66.351.60%6021244129241.3910.90%66.331.56%5525701158031.3920.97%66.331.56%5029017181391.4021.70%66.421.70%3088721576741.3931.04%65.850.83%10节点四面体单元202134111423301.3941.12%65.950.98%注：1）圆轴直径D=340㎜，长度L=1705㎜，扭矩Mn=504×106N.㎜2）精确解为最大剪应力τmax=65.31MPa3）应用软件为Unigraphics，全自动网格划分2.32.32.32.3四面体单元与六面体单元网格的比较四面体单元与六面体单元网格的比较四面体单元与六面体单元网格的比较四面体单元与六面体单元网格的比较表3为应用MSC.MARC程序采用相同节点排序进行圆轴受纯扭转时半自动网格划分时的有限元计算结果。单元类型以8节点六面体单元为主体。从表中可以看出，当单元网格逐渐加密，节点数和单元数不断增加时，与精确解的相对误差可从3.71%减小到0.36%。它虽然提高了求解精度，求解时间也迅速增加。与表1、表2的计算结果相比较，相同的单元网格密度条件下，8节点六面体实体单元的求解精度在4节点四面体单元与10节点四面体单元之间。表3不同密度8节点六面体实体单元为主体网格圆轴受纯扭转时的有限元计算结果单元网格参数计算分析结果及相对误差单元类型单元大小(㎜)节点总数单元总数位移(㎜)相对误差剪应力(MPa)相对误差求解时间(秒)90157013441.42983.71%67.613.52%270181615601.42903.66%67.603.51%350432337281.42053.04%67.303.05%6359168103571.39551.23%66.311.53%248节点6面体实体单元为主体2050836468721.38350.36%66.511.84%1006注：1）圆轴直径D=340㎜，长度L=1705㎜，扭矩Mn=504×106N.㎜2）精确解为最大剪应力τmax=65.31MPa3）应用软件为MSC.MARC，半自动网格划分3333单元单元单元单元网格网格网格网格质量质量质量质量的的的的评价评价评价评价一个高质量的单元网格划分，应当综合计算机容量、计算时间及工程实际中的精度要求等因素，来决定单元的大小及阶次，设定合理的网格数量和相对网格密度。选择高阶单元，网格密度可以相对稀疏一些。换句话说，在应力集中区域应采用较密的网格，而对于其它应力变化平缓的区域，则