工程仿真软件中的网格3：单元质量

有限元等仿真软件中网格的质量对于计算非常重要。高质量的单元可以提升计算速度，而低质量的单元会影响计算精度，甚至导致计算无法完成。常见的衡量单元质量的纲量有：纵横比，雅可比比值，最大内角，平行偏差，偏斜度，特征长度，正交质量，和翘曲因子等等。本文从实际应用与开发的角度，讨论常见的网格（单元）质量参数。 纵横比Aspect Ratio 纵横比描述单元的细长程度。用单元最长边/最短边的比值，表征单元是否过度细长；细长单元沿短边方向梯度捕捉会失效，使得弯曲、边界层应力误差巨大。理想值为1，数值越大表明单元质量越差。纵横比可以用于边界层、薄壁结构、拉伸网格的检查。 Image 对于通用的模型，整体的纵横比需要小于20，在应力集中区小于10。当纵横比大于20时，单元严重拉伸，有限元计算精度差、迭代难收敛。对于CFD中的边界层，法向加密区域允许 30~100的纵横比。 单元纵横比的数值方法直接。三角形通过计算最长边 / 等效内切圆直径的比值得到纵横比。四边形取两组对边平均长度比值max(L1, L2) / min(L1, L2)得到纵横比。四面体与三角形相似，基于内切球半径 r与外接球半径 R的比值得到纵横比。六面体取12条棱长度的最大 / 最小比值得到纵横比。通过计算方法可见，纵横比是一种纯几何长度计算，不使用高斯点、不使用雅可比。 雅可比比率Jacobian Ratio 雅可比比率是单元内最小与最大雅可比行列式的比值。雅可比变化的均匀性能反映扭曲梯度 ，可以衡量单元内部映射均匀性，专门捕捉高阶单元中节点扭曲、局部收缩膨胀，是二阶单元核心质控指标。值的范围是[-1, 1]，比值为1时，表明单元均匀无扭曲；比值越接近零，单元扭曲大。经验标准如下： 雅可比比值 0.3～1.0 0.1～0.3 0 ～0.1 小于等于0 单元质量 优质单元，求解稳定 畸变，应力精度下降，易震荡 严重畸变，不建议用于高精度计算 负雅可比 / 退化，求解报错、发散 雅可比比率的计算有两套标准，分别基于高斯积分点和顶点节点。由于高斯点的计算方式匹配求解器刚度矩阵积分位置，是有限元方法中可靠的网格指标。顶点节点计算方式仅几何顶点校验，计算效率高，但可能会出现 “顶点全正、积分点负” 漏判，常用于快速几何筛查。 Image 对于线性单元，雅可比行列式在单元内部为常数，雅可比比率恒为1，无法表现出线性单元的质量，此时可以对比理想单元的雅可比行列式，得到非恒为1的雅可比比率。同时由于单元内部雅可比行列式恒定，每个线性单元只需计算1 次，无需遍历所有节点或高斯点。 对于非线性单元，如二阶偏置或高阶单元，其形函数是多次的，此时单元内不同积分点的雅可比矩阵是不同的，需要在每个积分点上动态代入坐标计算。 最大内角度 Max Corner Angle 最大内角是一种直观的衡量单元质量的纲量，可以识别钝角极大、近乎扁平退化的单元。常用角度值作为单位，输出单元内最大角度。此数值越大，单元畸变越严重；越接近理想等角越好。四边形最优角度是90°，即直角；三角形最优角度是 60°；超过 150° 为劣质单元。超过165度的单元不建议使用。 Image 角度制的内角可以通过以下公式计算：θ=arccos(u⋅v /∥u∥⋅∥v∥ )⋅180/ π ，u 和v表示内角两条边的向量。三角形单元3 顶点的角度求得θ0, θ1, θ2，取其最大值 max { θ0, θ1, θ2 }，即为最大内角。同理，四边形单元的最大内角是max { θ0, θ1, θ2, θ3 }。 对于三维单元，取每一个面上的最大内角，再比较每个面上的最大内角得到单元最大内角。如四面体单元，4 个三角形面，每个面 3 个内角，取全部 12 个内角最大值。六面体有 6 个四边形面，每个面 4 个顶角，遍历所有24个面内角，即可得到该六面体的最大内角。 偏斜度 Skewness 偏斜度描述单元的倾斜扭曲，衡量四边形/六面体偏离矩形/立方体的程度。它是归一化的综合指标，同时惩罚过大角、过小角，是有限元求解器看重的单元质量，因此是最常用的通用单元指标之一。数值范围是[0, 1]，0表示单元理想正交，是标准规整单元，1表示单元完全歪斜，无法用于有限元计算。通常情况下，小于0.2表示优质单元；小于0.5的单元都可接受；大于0.8表示单元质量差，数值精度下降，迭代计算难以收敛、可能会导致应力奇异。 Image 偏斜度的通用计算方法如下， S=max{ (θmax−θideal)/(180−θideal), (θideal−θmin)/θideal } 其中，θmax是最大内角，θmax是最小内角。θideal是理想内角。三角形理想内角为60°；四边形是90°；四面体的理想二面角约 70.5°，六面体理想二面角是90°。对于三维单元如四面体和六面体单元，使用的是二面角的最大或最小值来计算偏斜度。此方法也叫等角法，支持所有的单元类型，包括棱柱与金字塔单元。 对于三角形和四面体单元，我们还可以使用体积法来计算偏斜度。公式如下：S= (Videal – Vactual)/ Videal，理想体积Videal是同外接球等边标准单元体积。 对于金字塔和楔形单元，可以分别计算所有四边形面的等角偏斜度，和所有三角形面的等体积偏斜度，取最大值即得到偏斜度。 单元综合质量 单元质量是一项综合网格质量评判指标，数值范围为 0～1。数值为 1 代表完美正方形（二维）或完美正方体（三维）；数值为 0 代表单元体积为零或出现负体积（网格畸变失效）。 Image 该指标的计算方式如下： 对于二维四边形 / 三角形单元，采用单元面积与各边长平方和的比值计算；Q = C* (A/∑L²)。三角形单元的面积可以通过向量叉乘模长的一半 得到。四边形单元可以划分为两个相邻的三角形单元，通过这个方法可以计算四边形单元的面积。 对于三维单元，则采用单元体积与各边长平方和立方的平方根的比值计算。Q = C (V/(∑L²)^1.5)。三维单元体积的可以通过计算单元的雅可比矩阵行列式实现。 常数C的大小取决于单元类型，具体数值如下： 单元 三角形 四边形 四面体 六面体 楔形 金字塔形 常数值 6.928 4.0 124.70 41.57 62.354 96 在单元质量评价中，质量值 >0.7为合格单元；应力集中区的单元值期望大于0.8。当单元值小于0.05时，存在严重坏单元，必须修复。 总结 本文讨论了常见单元类型的单元质量，计算并显示单元质量是现代通用有限元的一个必备功能。计算这些纲量不仅用于有限元网格划分，也用于网格的交互检查，云图显示单元质量，用户可以快速直观的了解网格的质量，并决定是否使用当前网格用于后续的计算。 本文未对三棱柱单元与楔形单元作详细讨论，但其单元质量的计算方式与四面体和六面体相似，可以用相同的方法计算得出。本文未对平行偏差，正交质量，翘曲因子等单元质量纲量做讨论。