有限元收敛性

本文讨论了有限元网格的重要概念 包括单元的分类 有限元误差的分类与影响因素 并讨 论分析结果的收敛性控制方法 并由实例说明了网格质量及收敛性对取得准确分析结果的 重要性 同时讨论了一些重要网格控制的建 议及其他网格设定的说明 一 基本有限元网格概念一 基本有限元网格概念 1 单元概述 几何体划分网格之前需要确定单元类型 单元类型的选择应该根据分析类型 形状特 征 计算数据特点 精度要求和计算的硬件条件等因素综合考虑 为适应特殊的分析对象 和边界条件 一些问题需要采用多种单元进行组合建模 2 单元分类 选择单元首先需要明确单元的类型 在结构有限元分析中主要有以下一些单元类型 平面应力单元 平面应变单元 轴对称实体单元 空间实体单元 板单元 壳单元 轴对 称壳单元 杆单元 梁单元 弹簧单元 间隙单元 质量单元 摩擦单元 刚体单元和约 束单元等 根据不同的分类方法 上述单元可以分成以下不同的形式 3 按照维度进行单元分类 根据单元的维数特征 单元可以分为一维单元 二维单元和三维单元 一维单元的网格为一条直线或者曲线 直线表示由两个节点确定的线性单元 曲线代 表由两个以上的节点确定的高次单元 或者由具有确定形状的线性单元 杆单元 梁单元 和轴对称壳单元属于一维单元 如图 1 图 3 所示 二维单元的网格是一个平面或者曲面 它没有厚度方向的尺寸 这类单元包括平面单 元 轴对称实体单元 板单元 壳单元和复合材料壳单元等 如图 4 所示 二维单元的形 状通常具有三角形和四边形两种 在使用自动网格剖分时 这类单元要求的几何形状是表 面模型或者实体模型的边界面 采用薄壳单元通常具有相当好的计算效率 三维单元的网格具有空间三个方向的尺寸 其形状具有四面体 五面体和六面体 这 类单元包括空间实体单元和厚壳单元 如图 5 所示 在自动网格划分时 它要求的是几何 模型是实体模型 厚壳单元是曲面也可以 4 按照插值函数进行单元分类 根据单元插值函数多项式的最高阶数多少 单元可以分为线性单元 二次单元 三次 单元和更高次的单元 线性单元具有线性形式的插值函数 其网格通常只具有角节点而无边节点 网格边界 为直线或者平面 这类单元的优点是节点数量少 在精度要求不高或者结果数据梯度不太 大的情况下 采用线性单元可以得到较小的模型规模 但是由于单元位移函数是线性的 单元内的位移呈线性变化 而应力是常数 因此会造成单元间的应力不连续 单元边界上 存在着应力突变 如图 6 所示 二次单元的插值函数是二次多项式 其网格不仅在每个顶点处有角节点 而且在棱边 上还存在一个边节点 因此网格边界可以是二次曲线或者曲面 这类单元的优点是几何和 物理离散精度较高 单元内的位移呈二次变化 应力呈线性变化 因此单元边界上的应力 是连续的 但是在单元数量相同的条件下二次单元的节点数比线性单元的节点数多 模型 的规模较大 如图 7 和图 8 所示 三次单元的插值函数是三次多项式 其网格的每条边上存在两个节点 有些三次单元 还具有内部节点 这类单元的离散精度更高 但是由于单元节点数较多 网格划分较为困 难 模型规模很大 一般用于具有特殊精度要求的场合 如图 9 所示 对于一阶和二阶单元 我们通常也称其为 H 单元 三阶及以上的单元 我们也称其为 P 单元 高阶次的 P 单元可以更好地拟合变形形状 特别对于曲率或者应力梯度变化较大 的区域会较为真实的模拟 但会比 H 单元有较多的运算量 如图 10 所示 5 结构单元与非结构单元 根据单元能否离散成实际结构 可以将单元分为结构单元和非结构单元 能离散成实 际结构的称为结构单元 如轴对称单元离散轴对称结构 杆 梁单元用于离散杆件结构 实体单元用于离散空间结构等 这些单元都属于结构单元 除此之外 还有一类单元并不 用于实际结构的离散 而是在模型中模拟一些特殊的结构和边界条件 如质量单元用于实 际的物体质量效应 弹簧和阻尼单元用于模拟结构的弹性支承和减振吸能部件 间隙和接 触单元用于结构之间的相互接触作用 螺栓预紧力单元用于模拟螺栓的预紧力 刚体单元 用于模拟节点之间的刚性连接等 这些单元称为非结构单元 由于非结构单元非常抽象 使用起来有一定的难度 在设计仿真一体化分析里面通常 会将其工程化 帮助使用者淡化其力学概念 6 节点和单元的重要力学概念 针对前述的单元分类 此处要澄清关于节点和单元的一些非常有用的总结性概念 有限元分析首先计算节点的位移量 接着再推算其对应单元的应变值 再计算积分 点的应力 因此位移的准确性高于应变 应变高于应力 当结构静力平衡时计算变形的单元是求得准确有限元分析结果的关键 因此线性计 算中单元不可以变形过大 否则会造成求解失败 网格质量概括来说 初始网格必须可呈现初始模型的几何形状 而且要足够 弹性 以符合静力平衡后的变形几何形状 在预计会有应力梯度变化剧烈的位置上 为预测其准确变形情况 细小特征几何必 须要更精确符合 以利于准确计算这些位置上的应力值 在理想曲率边线与网格曲率边线之间的差距称之为离散误差 二 有限元误差分析二 有限元误差分析 1 有限元误差 有限元的误差主要来自两个方面 一是模型误差 一是计算误差 模型误差是指将实际工程问题抽象为适合计算机求解的有限元模型时所产生的误差 即有限元模型和实际问题之间的差异 它包括有限元离散处理所固有的原理性误差 也可 能包括几何模型处理 实际工况转化为模型边界条件时所带来的偶然性误差 计算误差是指采用数值方法对有限元模型进行计算所产生的误差 误差的性质是舍入 误差和截断误差 模型误差包含离散误差 边界条件误差和单元形状误差 离散误差包含物理离散误差 和几何离散误差 2 离散误差 物理离散误差是插值函数和真实函数之间的差异 其大小与单元尺寸和插值多项式的 阶次有关 单元尺寸减小也就是网格划分越密 插值函数的阶次增加 将使有限元的解收 敛于精确解 几何离散误差是指离散后的几何体与原有几何形状上的差异 对于由直线或者平面边 界构成的规则结构 这类误差较小 对于具有复杂曲线或者曲面边界构成的结构 离散后 会产生较大的形状误差 本文下面通过 SolidWorks Simulation 来详细讨论物理离散误差与几何离散误差的具体 操作细节 三 收敛性及自动收敛方法三 收敛性及自动收敛方法 一般而言 网格拥有较多的单元 可得到较准确的结果 会有更多的节点可供计算 所以结果会较准确 较多的单元也就表示单元大小较小 所以物理离散误差可减小 实际 分析上也有极限 在收敛性分析过程中网格尺寸一再缩减也不一定会对精确结果有帮助 对一给定的几何而言 要达到收敛性的网格会与外部负载条件及边界约束条件有关 见 图 11 在线性静态分析中 载荷大小不是收敛性的系数 下面以 SolidWorks Simulation 的收敛性为例简单介绍收敛性的处理方法和技巧 SolidWorks Simulation 提供三种收敛性的技术 包括有手动控制收敛性和软件自动控 制收敛性技术 其中自动控制收敛性的方法我们也称之为自适应方法 如图 12 所示 包 括自动 H 自适应方法 H adaptive 自动 P 自适应方法 P adaptive 1 自动 H 自适应方法 H adaptive H 方法的本质就是根据应力梯度的变化情况自动在应力梯度大的地方 根据预先规定 的收敛准则 重新自动剖分网格 进行自动加密 见图 13 原始网格与 H 自适应网格结果 SolidWorks Simulation 的 H 方法具有以下特性 1 适用于实体零件及装配体 仅支持实体单元 的静态分析研究 2 在应变能误差较高的区域使用较小网格尺寸 3 可以在应变能误差较低区域网格粗化 加大网格尺寸 便于在后面的优化计算中降 低计算规模 大大提高优化效率 4 目标精度定义应变能量密度范数的精度等级默认值是 98 此处可以调整能量密度 范数的精度等级 一般情况下默 认的精度可以达到分析的要求 5 精度偏差设置 见图 13 精度偏差设置有局部 本地 和全局 整体 滑动杆朝局部移动 指示程序以较少的单元取得精确 的峰值应力结果 滑动杆朝全局移动 指示程序取得整体零件刚度精确 的结果 而不是应力结果 6 若不确定 保持默认值即可 7 网格粗糙化的目的是对应力梯度变化不大的区域 加大此处的单元尺寸 可以使用 较小的网格得到较好的结果 同时也便于后继的优化求解 图 14 展示了某机械零件 采用一阶单元不同单元大小 并采用自适应方法进行分析 方案 1 采用平均单元大小 13 6mm 进行网格划分 然后采用 H 自适应网格划分 此时得 到的最大应力点的应力是 44MPa 方案 2 采用平均单元大小为 3 4mm 进行网格划分 然 后采用 H 自适应网格划分 此时得到的最大应力点的应力是 75 6MPa 两者之间有 42 的差异 说明采用线性单元 使用 H 方法得不到准确的结果 同样的模型 不做任何修改采用二阶单元进行网格划分 如图 15 所示 然后采用 H 方法 得到的应力误差小于 3 说明采用较大的全局二阶单元 然后采用 H 方法可以得 到相当准确的应力结果 2 自动 P 自适应方法 P adaptive P 方法的本质就是根据约束条件 如应变能 的变化情况自动在约束条件大的地方 根 据预先规定的收敛准则 调整该处的单元形函数的阶次 在单元大小不变的情况下提高单 元内部应力的精确性 见图 16 SolidWorks Simulation 的 P 方法具有以下特性 1 适用于实体零件及装配体的静态分析研究 但装配体仅支持结合方式 不可以有其 他接触存在 2 收敛准则有总应变能 均方根合位移 均方根 von Mises 应力 3 默认收敛准则是总应变能 均方根合位移及 von Mises 应力准则并不常用 4 默认的设定通常就足够 由于系统通常会提前满足设 定精度 因此最大 p order 及 最大循环数很少用到 5 开始 p 阶序起始于 2 设为 1 会报错 6 必须使用二阶单元为初始网格 一般而言初始的网格尺寸影响很小 见图 17 Jacobian 雅可比检查 对误差有较大的影响 在局部位置无法反馈结果 建议将 Jacobian 检查设定在节点 见图 18 四 手动收敛性检查四 手动收敛性检查 1 相对收敛性检查 在大多数复杂情况下很难通过自适应方法得到好的结果 必须通过相对收敛性检查得 到收敛的结果 见图 19 其操作步骤如下 1 执行多个分析研究 逐步调整加密网格 检查应力值的变化情形 2 每次以 2 1 比例调整加细网格尺寸 3 如果局部网格尺寸远小于整体网格尺寸 要留意扭曲 失真的情况 2 等值线质量检查 见图 20 1 应力等值线应该和连续几何体一样连续 使用不连续选项可以更清楚看到不连续的 结果 如果几何体光滑连续而结果呈锯齿状 表明此处结果不好 需要加密网格或提高网 格质量 2 没有一个收敛的绝对测试标准 但是可以显示区域检查 3 显示带有结果的单元结果可以显示他们如何影响结果 3 误差估算方法 1 能量范数值 见图 21 1 ERR 能量范数误差绘图可以显示出相邻元素之应力值差异 理论上要愈小愈好 2 并非绝对是收敛性测试方法 但会显示出要特别注意的高误差区域 3 最佳的应用方式是在比较不同的设计过程 如果在某个区域的误差估算已经调整至 最佳状况并维持一致 其分析结果的比较应该较具有实际意义 4 误差估算方法 2 节点和单元应力值比较 见图 22 1 节点解是临近单元的节点应力的平均值 2 单元解是每个单元所有节点应力的平均值 3 评定标准 理论上节点和单元应力值应该有较小的差异 4 一般情况下 节点应力和单元应力的误差不允许超过 5 五 手动收敛性控制和自动收敛性