有限元法分析结果的误差影响，四类误差您了解吗？
本文指出了有限元法分析结果的误差影响存在于其每一操作步骤，并对这些误差进行了归类分析。随后，结合工程实例，通过改变单元类型（形状和精度）、调整单元尺寸大小和应用多种分网方式，显示理想化误差和离散化误差对计算结果的影响。最后，提出建议和今后的研究方向。
有限元法分析起源于50年代初杆系结构矩阵的分析。随后，Clough于1960年第一次提出了“有限元法”的概念。其基本思想是利用结构离散化的概念，将连续介质体或复杂结构体划分成许多有限大小的子区域的集合体，每一个子区域称为单元（或元素），单元的集合称为网格，实际的连续介质体（或结构体）可以看成是这些单元在它们的节点上相互连接而组成的等效集合体；通过对每个单元力学特性的分析，再将各个单元的特性矩阵组集成可以建立整体结构的力学方程式，即力学计算模型；按照所选用计算程序的要求，输入所需的数据和信息，运用计算机进行求解。
当前，有限元方法/理论已经发展的相当成熟和完善，而计算机技术的不断革新，又在很大程度上推进了有限元法分析在工程技术领域的应用。然而，如此快速地推广和应用使得人们很容易忽视一个前提，即有限元分析软件提供的计算结果是否可靠、满足使用精度的前提，是合理地使用软件和专业的工程分析。只有这两者很好地结合，我们才能得到工程上切实可信的计算结果，否则只会在工程上造成极大的浪费，甚至带来严重的工程事故。
有限元法分析一般包括四个步骤：物理模型的简化、数学模型的程序化、计算模型的数值化和计算结果的分析。每一个步骤在操作过程中都或多或少地引入了误差，这些误差的累积最终可能会对计算结果造成灾难性的影响，进而蒙蔽我们的认识和判断。
第一步，物理模型的简化，主要有几何实体、连接/装配关系、环境边界条件和材料特性的简化，进而构建数学模型。这些简化或者说假设，是必要的，也是必须的，但是也由此在模型中引入了理想化误差（idealization error）。
有些理想化误差是非良性奇异的，比如几何实体简化时细节部位上忽略小的圆/倒角，连接/装配关系简化时忽略焊缝和螺栓连接等，往往导致模型发生结构方面（诸如L形截面的角点）的奇异，即结构奇异（奇异的数学定义是在某一点处导数无穷）；
有些理想化误差是良性奇异的，比如边界条件简化时添加集中载荷和孤立点约束，导致模型发生边界条件的奇异，即边界奇异；
其它理想化误差，比如几何实体简化时三维壳/面体简化为二维壳/面、三维梁简化为一维梁，边界条件简化时非均匀温度场和压力场简化为均匀温度场和压力场等，只会影响计算结果的准确度，不会引发计算结果方面的数值奇异，即应力奇异和位移奇异等。
理想化误差是在有限元法分析开始之前引入的，因此我们不可能通过改进有限元分析技术来达到消除其的目的，而只能通过修改数学模型本身来实现消除其的目的。
第二步，数学模型的程序化，主要有几何实体的单元离散、单元网格的装配连接、模型环境边界条件的添加，进而构建计算模型。几何实体的离散，和单元类型（形状和精度）、单元尺寸以及分网方式的选择有关，不可避免地会引入离散化误差（discretization error）。
离散化误差，是根植于有限元法分析本身的，因此只能通过改进有限元分析技术或者技巧来尽力消除/减小这方面的误差，比如采用规则化的单元形状避免单元在形状上产生奇异（即单元奇异）、提高单元精度和增加网格密度减小计算方面的误差等方法。
单元网格的装配连接一般采用MPC多点约束法，因而会引入人为误差（artificial error），这方面误差的消除更多是需要长期计算经验的积累。模型环境边界条件的添加，其误差影响依赖于第一步的理想化简化。
第三步，计算模型的数值化，主要是用数值计算方法（程序求解器）求解、逼近真实的解析值，因而必然存在数值化误差（numerical error）。
数值计算方法的精度（非人为可控）越高，计算结果的误差就越小，但计算的工作量也越大。实际考虑到计算精度和计算资源的利用，必然要做一个适当的统一。
第四步，计算结果的分析，主要是利用数值计算结果来分析、评判，或预知真实的物理模型，由此也存在着认知误差（recognized error）。
认知误差的消除，一方面需要真实物理试验的指导，另一方面依赖于分析人员的工程经验和认知能力。同时，不要忘记了我们的前提假设，即第一步物理模型的简化，或假设。
下文，将通过一个简单的例子来说明理想化误差和离散化误差对有限元法分析结果的影响。计算时，采用有限元数值分析软件ANSYS11.0版本，32位操作系统软件WindowsXP版本，HP xw4200服务器硬件平台，保证了程序求解器及其运行环境的统一，以消除数值化误差。
图1 起竖支耳模型
图1中所示为工程上最常见的起竖支耳模型，其包括两个部分：横板和竖耳。工艺生产上，即可以将横板和竖耳做成一个整体铸件（如图1.a），也可以将二者作为两个单独的部件焊接而成（如图1.b）。
有限元分析时，一般会忽略铸件上小的过渡圆角，也经常忽略焊接件上的焊缝，即而简化为适于分析的几何实体（如图1.c）。工程使用中，起竖机构通过销轴作用于竖耳销孔以推动横板连接机构完成起竖，具体压力载荷P的数学形式可表示为：
上式(1)中，P0为最大压力载荷幅值，θ为载荷作用面上某点的周向角度，F为真实的载荷力，R为销孔半径，L为销孔纵向长度。数学建模时，仅为满足工程简单计算的需求，材料模型通常取为线性弹性模型。
为了显示计算结果的误差影响，这里选用三种评估方法：智能化自由网格划分、规则化网格划分和自适应P改进，见表1。
表1 三种评估方法
智能化自由网格划分，利用SMRTSIZE选项控制单元尺寸大小，分割几何实体为四面体单元，单元形状较差，计算效率低；规则化网格划分，利用线段分割数参数NSize控制网格单元大小，保证几何实体规则划分，以避免单元奇异，进而消除离散化带来的误差影响；自适应P改进在规则化网格划分的基础上，通过提高单元精度以展示所关心位置处节点的应力收敛过程，进而显示理想化误差带来的影响。
关心位置包括：几何简化引起结构奇异的位置，竖耳根部左下角LUGDL、右下角LUGDR、左上角LUGUL和右上角LUGRL，这些位置的节点应力通常不会收敛，即应力奇异；边界约束引起边界奇异的位置，横板约束根部左下角PLATEDL、右下角PLATEDR、左上角PLATEUL和右上角PLATEUR，这些位置受刚性约束限制了其横向位移（泊松比效应）从而导致应力不收敛；实体离散可能引起单元奇异的位置，销孔左下角PINDL、右下角PINRL、左上角PINUL和右上角PINUL，这些位置若网格规则其应力集中，若网格不规则其应力奇异。
1. 智能化自由网格划分
按照有限元分析软件前处理步骤，首先定义了单元类型（三维十节点四面体结构单元SOLID187）、实常数和材料模型参数；然后，采用智能化自由网格划分技术对几何实体图1.c进行分网（见图2所示），单元尺寸大小由SMRTSIZE参数控制；最后，对模型进行检查，确保无误后退出前处理模块。
图2 有限单元模型
进入有限元分析软件计算模块，首先给计算模型添加环境边界条件：横板左、右面固定约束，竖耳销孔位置施加式(1)形式的载荷，其合力沿水平正Z轴方向；随后，设置求解环境参数，静态分析；最后，对模型进行检查，确保无误后求解退出计算模块。
图3 销孔上下边缘应力变化曲线
图4 竖耳根部角点应力变化曲线
图5 横板约束角点应力变化曲线
由图3-图5所示，可以发现：
1) 竖耳根部、销孔上下缘和横板约束根部的节点应力具有比较明显的网格敏感特性，当网格密度达到一定程度时其整体趋势是随着网格密度的增加（SMRTSIZE越小）应力值无限增加，即不收敛；
2) 由于网格不太规则，使得载荷和位移边界条件偏离对称特性，因而对称位置的应力数值有明显偏差。
众所周知，特定工况下结构的真实应力，其数值只能有一个，具有唯一性。考虑到数值计算的精度问题，其应力计算值可能是一个围绕真实应力在一定精度范围内波动的数值，而随着数值精度的不断提高，应力计算值将会逐渐趋于一个有限数值，即真实应力。由此看来，图3~图5所示关心位置处应力并非真实的应力，而是计算得到的虚假应力。
2. 规则化网格划分
采用三维二十节点六面体实体单元SOLID186，利用线段分割数参数NSize控制网格单元大小，扫略划分图1.c几何实体，见图6所示。
图6 有限单元模型
由图7~图8所示，可以看出：
1) 由于网格规则，使得载荷和位移边界条件满足对称特性，因而对称位置的应力数值基本完全重合；
2) 竖耳根部和横板约束根部的节点应力仍具有比较明显的网格敏感特性，其整体趋势是随着网格密度的增加（NSize越大）应力值无限增加，即不收敛；
3) 销孔上下边缘的节点应力在4.45~4.50MPa范围内波动，且有收敛态势，因而其应力计算值为真实应力。
以上分析表明，销孔上下边缘的应力结果误差为离散化误差，而竖耳根部角点和横板约束角点的应力结果误差为理想化误差。
图3 销孔上下边缘应力变化曲线
图8 竖耳根部角点应力变化曲线
图9 横板约束角点应力变化曲线
3. 自适应P改进
采用三维二十节点六面体实体P单元SOLID147（其位移模式中形函数的阶次可以在2~8之间自适应选择），设置线段分割数参数NSize=3控制网格单元大小，扫略划分图1.c几何实体，见图10所示。
图10 有限单元模型
由图11~图13所示，可以看出：
1) 竖耳根部和横板约束根部的节点应力随单元精度的提高无限增加，不收敛；
2) 销孔上下边缘的节点应力随单元精度的提高逐渐收敛于4.39MPa左右。
由此再次表明，销孔上下边缘的应力结果误差为离散化误差，而竖耳根部角点和横板约束角点的应力结果误差为理想化误差。
图11 销孔上下边缘应力变化曲线
图12 竖耳根部角点应力变化曲线
图13 横板约束角点应力变化曲线
本文转载自e-works网，原文来自于安亚世太（peraglobal.com），作者：商霖。
声明：本微信转载文章出于非商业性的教育和科研目的，并不意味着赞同其观点或证实其内容的真实性。版权归原作者所有，如转载稿涉及版权等问题，请立即联系我们，我们会予以更改或删除相关文章，保证您的权利！