首页/文章/ 详情

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

26天前浏览249

有限元等仿真软件中网格的质量对于计算非常重要。高质量的单元可以提升计算速度,而低质量的单元会影响计算精度,甚至导致计算无法完成。常见的衡量单元质量的纲量有:纵横比,雅可比比值,最大内角,平行偏差,偏斜度,特征长度,正交质量,和翘曲因子等等。本文从实际应用与开发的角度,讨论常见的网格(单元)质量参数。

纵横比Aspect Ratio

纵横比描述单元的细长程度。用单元最长边/最短边的比值,表征单元是否过度细长;细长单元沿短边方向梯度捕捉会失效,使得弯曲、边界层应力误差巨大。理想值为1,数值越大表明单元质量越差。纵横比可以用于边界层、薄壁结构、拉伸网格的检查。

对于通用的模型,整体的纵横比需要小于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

单元质量

优质单元,求解稳定

畸变,应力精度下降,易震荡

严重畸变,不建议用于高精度计算

负雅可比 / 退化,求解报错、发散

雅可比比率的计算有两套标准,分别基于高斯积分点和顶点节点。由于高斯点的计算方式匹配求解器刚度矩阵积分位置,是有限元方法中可靠的网格指标。顶点节点计算方式仅几何顶点校验,计算效率高,但可能会出现 “顶点全正、积分点负” 漏判,常用于快速几何筛查。

对于线性单元,雅可比行列式在单元内部为常数,雅可比比率恒为1,无法表现出线性单元的质量,此时可以对比理想单元的雅可比行列式,得到非恒为1的雅可比比率。同时由于单元内部雅可比行列式恒定,每个线性单元只需计算1 次,无需遍历所有节点或高斯点。

对于非线性单元,如二阶偏置或高阶单元,其形函数是多次的,此时单元内不同积分点的雅可比矩阵是不同的,需要在每个积分点上动态代入坐标计算。

最大内角度 Max Corner Angle

最大内角是一种直观的衡量单元质量的纲量,可以识别钝角极大、近乎扁平退化的单元。常用角度值作为单位,输出单元内最大角度。此数值越大,单元畸变越严重;越接近理想等角越好。四边形最优角度是90°,即直角;三角形最优角度是 60°;超过 150° 为劣质单元。超过165度的单元不建议使用。

角度制的内角可以通过以下公式计算:θ=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表示单元质量差,数值精度下降,迭代计算难以收敛、可能会导致应力奇异。

偏斜度的通用计算方法如下,

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 代表单元体积为零或出现负体积(网格畸变失效)。

该指标的计算方式如下: 对于二维四边形 / 三角形单元,采用单元面积与各边长平方和的比值计算;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时,存在严重坏单元,必须修复。

总结

本文讨论了常见单元类型的单元质量,计算并显示单元质量是现代通用有限元的一个必备功能。计算这些纲量不仅用于有限元网格划分,也用于网格的交互检查,云图显示单元质量,用户可以快速直观的了解网格的质量,并决定是否使用当前网格用于后续的计算。

本文未对三棱柱单元与楔形单元作详细讨论,但其单元质量的计算方式与四面体和六面体相似,可以用相同的方法计算得出。本文未对平行偏差,正交质量,翘曲因子等单元质量纲量做讨论。

来源:WELSIM
ACT非线性通用WELSIM
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2026-07-30
最近编辑:26天前
WELSIM
一枚搞仿真的老员工
获赞 27粉丝 87文章 271课程 0
点赞
收藏
作者推荐

生成CalculiX求解器的接触与边界条件

CalculiX 是一款著名的开源有限元分析(FEA)软件,主打结构力学(静力、动力、非线性)与热分析,适用于科研、教学和中小型工程场景。CalculiX的核心优势有:开源免费。个人、企业、科研机构均可自由使用、修改和分发代码,遵循 GNU GPL 许可证。兼容支持 Abaqus 的 INP 输入文件格式,多数 Abaqus 模型也可直接导入 CalculiX 求解,减少学习曲线。求解器代码精简,对硬件资源要求低于商业软件,适合中小型模型的快速计算。对第三方前后端软件支持也非常好,易于前后处理器与求解器的整合。笔者曾在《使用WELSIM生成CalculiX求解器文件》一文中,介绍了如何快速地生成CalculiX输入文件。目前,WELSIM已经可以很好的作为CalculiX的前处理器生成输入文本,极大的方便用户使用CalculiX进行有限元分析。本文以更详细的方式介绍生成接触与边界条件的命令。接触在多体有限元结构分析中,接触是必要的设置之一。本章节描述生成两种常见的接触类型:绑定接触和可分离接触。绑定接触绑定接触是一种用于模拟多体之间完全刚性连接的接触条件,适用于部件间无相对滑移、无分离、无间隙的工况(如焊接、胶粘、过盈配合后固化、螺栓紧固的刚性连接)。绑定接触会强制接触界面的节点位移完全一致,等同于将多个独立部件 “融合” 为一个整体结构。 绑定接触设置如下图所示,转换成CalculiX命令时,使用*TIE 命令。生成的输入命令如下*Tie, Name=ID_15, Position tolerance=0.01, Adjust=NoTarget_Surface_15, Master_Surface_15在网格文件中,定义的表面单元如下:** Box_to_Box_1*SURFACE, TYPE=ELEMENT, NAME=Master_Surface_15 70, S2 63, S2 74, S2 61, S2 66, S1 69, S2** Box_to_Box_1*SURFACE, TYPE=ELEMENT, NAME=Target_Surface_15 166, S3 185, S2 160, S3 173, S2 171, S4 169, S3可分离接触可分离接触是一种用于模拟多体接触界面允许接触、滑移、分离三种状态的非线性接触条件,核心特点是接触界面仅传递法向压力,无拉力作用,且切向可自由滑移或黏着,适用于部件间存在动态接触与分离的工况(如齿轮啮合、冲压成型、轴承滚动、机械碰撞)。根据摩擦力假设,又分为无摩擦接触和有摩擦接触。WELSIM设置接触如下图所示,转换成CalculiX命令时,使用CONTACT命令。这里值得注意的是,CalculiX规定,主接触面又称作独立接触面,必须是以单元定义的表面。被动接触面又称作是依赖接触面,可以是单元定义的表面,也可以是节点定义的表面。生成对应的CalculiX命令为** Box_to_Box 1*Surface interaction, Name=SurfaceInteraction_15*Surface behavior, Pressure-overclosure=Hard*Friction0.15*Contact Pair, Interaction=SurfaceInteraction_15, Type=Surface to surfaceTarget_Surface_15, Master_Surface_15结构分析中的边界条件设置边界条件是有限元分析中重要的步骤,也是前处理软件生成求解器命令的重要部分。WELSIM支持大量CalculiX的边界条件。本章节演示前处理器生成结构分析中的边界条件,并显示对应的求解器命令。1. 固定边界或位移位移边界条件是用于约束模型节点位移自由度的本质边界条件,其作用是限制结构的刚体运动、模拟实际支撑约束(如固定支座、铰支座),是结构静力、动力分析的基础条件。 3维模型中,对于实体单元,可以施加X/Y/Z方向的平移(Ux/Uy/Uz),对于壳单元,可以施加额外的绕 X/Y/Z 轴转动的转动(Rx/Ry/Rz)。对应生成的CalculiX求解器命令是,** BC Name: Displacement *Boundary ID_22, 1, 1, 1 *Boundary ID_22, 2, 2, 2 *Boundary ID_22, 3, 3, 3 *Boundary ID_22, 6, 6, 0.5235992. 压力压力边界条件是一种作用于结构表面的面载荷,属于分布载荷,其本质是通过指定单位面积上的法向力来模拟流体压力、接触压力、气体载荷等工程场景(如压力容器内壁压力、风载荷、土壤侧压力)。压力是结构有限元分析中常见的边界条件之一。压力的符号规则决定载荷方向:正压力:载荷方向指向表面内侧,对结构产生压缩效应(如压力容器内壁压力);负压力:载荷方向背离表面外侧,对结构产生拉伸效应(如真空容器外壁的大气压力)。压力边界条件如图所示,对应生成的求解器命令是,** BC Name: Pressure *DLOAD ID_20, P, 1233. 力力边界条件是作用于几何关键点、线、或者面的载荷,用于模拟工程中局部位置的受力工况(如螺栓预紧力、吊环拉力、销钉作用力)。 力是结构力学中最常见的边界条件,然而, 对于有限元数值方法中,几何面上的力边界条件的施加并不简单,通常是将多个节点的自由度耦合到一个参考节点,常用于将集中载荷分散到节点组(如把力加在参考节点,自动分配到耦合组)。如图所示,在前处理器的力设置中,用户可以选定几何体的面,线,或点,并给定三个方向上的力大小。得到力的求解器命令如下,*COUPLING, REF NODE=46, SURFACE=ID_17, CONSTRAINT NAME=c17 *KINEMATIC 1 3** BC Name: Force *CLOAD, OP=NEW 46, 1, 100 *CLOAD, OP=NEW 46, 2, 200 *CLOAD, OP=NEW 46, 3, 3004. 节点力节点力边界条件是直接作用于单个或多个节点的集中力载荷,是有限元离散化后最基础的载荷施加形式,与之前提到的集中力边界条件高度关联,但更偏向于网格离散后的直接操作,常用于精细化载荷控制或特殊工况模拟。 节点力是直接施加在网格节点上的力,根据节点数量大小,综合施加的力会叠加。生成的求解器命令也非常简单明了,根据所选的节点,施加CLOAD载荷。 ** BC Name: Nodal Force *CLOAD ID_18, 1, 10 *CLOAD ID_18, 2, 20 *CLOAD ID_18, 3, 30** Nodal_Force*NSET, NSET=ID_18 82, 84, 86, 88, 95, 96, 103, 104, 107, 108, 111, 112, 145, 146, 147, 148, 149, 1505. 重力重力边界条件属于体积力(Body Force) 范畴,用于模拟地球引力场对结构的作用,是通过给模型的每个单元分配与质量成正比的惯性力,适用于几乎所有包含自重影响的结构分析(如建筑梁板、机械部件、航天器地面工况)。 重力的本质是加速度场作用下的体积力,其大小与单元的质量直接相关。前处理器设置如下,生成的求解器命令如下图所示。** DC Name: Earth Gravity *Dload Eall, Grav, 9.807, 0, 0, -16. 旋转角速度旋转角速度边界条件用于模拟结构绕某一轴做匀速转动的工况,会在结构内部产生离心力和科里奥利力(高速转动时需考虑),属于惯性载荷范畴,广泛应用于旋转机械分析(如飞轮、叶轮、离心机、涡轮转子)。输入的参数有旋转轴,旋转原点,和角速度的大小。输入界面如下图所示。生成的求解器命令如下所示。** DC Name: Rotational Velocity *Dload ID_21, CENTRIF, 100, 5, 5, 0, 0, 0, 1热分析中的常见边界与初始条件1. 温度在有限元热分析中,温度边界条件是直接定义模型特定区域温度值的约束条件,属于本质边界条件(Dirichlet 边界条件),适用于已知表面或节点温度的场景(如恒温水箱壁面、与大热源接触的部件、温控设备的设定温度面)。 在热传导控制方程中,温度边界条件是强制约束,会直接固定边界节点的温度自由度,求解时优先满足该约束,再计算内部温度场分布。 前处理器设置如下,生成的CalculiX求解器命令如下,** BC Name: Temperature *Boundary, op=New ID_23, 11, 11, 02. 热流密度热流密度(Heat Flux)边界条件是直接定义单位面积上通过固体表面的热流量的边界设置,属于热载荷的一种,适用于已知表面热流输入 / 输出的场景(如电加热器表面、太阳能集热板、高温燃气冲刷壁面)。热流密度是核心参数,正方向为热量流入固体,负方向为热量流出。输入界面如下:生成的CalculiX求解器命令如下:** BC Name: Heat Flux *Dflux ID_24, S, 0.53. 体热流密度体积热流是一种体载荷,用于模拟物体内部热源产生或消耗热量的过程,区别于作用在表面的热流密度边界条件,它直接作用于模型的体单元,适用于内热源驱动的温度场分析。体积热流密度,SI单位为 W/m3(正号表示生成热量,负号表示消耗热量),用户输入界面如下:生成的CalculiX求解器命令如下, *DFLUXID_31,BF,10.4. 对流对流(Convection)边界条件又称膜(Film)边界条件。是用于模拟固体表面与相邻流体(气体或液体)之间对流换热的边界设置,是工程热分析中最常用的边界条件之一,广泛应用于电子设备散热、管道换热、汽车发动机冷却等场景。 对流换热系数是表征对流换热强度的核心参数,环境温度参数是远离固体表面的主流流体温度。前处理器中的输入界面如下,生成的CalculiX求解器命令如下:** BC Name: Convection *Film, op=New ID_29_S1, F1, 80, 123 ID_29_S2, F2, 80, 123 ID_29_S3, F3, 80, 123 ID_29_S4, F4, 80, 1235. 热辐射边界条件辐射边界条件是用于模拟物体表面与周围环境(或其他物体表面)之间通过热辐射方式传递热量的边界设置,其核心遵循斯蒂芬–玻尔兹曼定律,适用于真空、气体等介质的热交换场景(如航天器热控、高温设备散热)。辐射边界条件,是一种非线性边界条件,热流与温度的四次方成正比,求解时需要迭代计算。用户需要输入的参数为发射率(Emissivity),介于0-1之间的数值,和环境温度。输入界面如下,生成的热辐射边界条件命令如下,** BC Name: Radiation *Radiate, op=New ID_30_S1, R1, 30, 0.9 ID_30_S2, R2, 30, 0.9 ID_30_S3, R3, 30, 0.9 ID_30_S4, R4, 30, 0.9总结有限元分析中,边界条件的类型较多,涉及了条件的参数本身,也涉及到所选择的单元,或节点。对于前处理器来说,生成各种边界条件的命令,是一项复杂的工程。WELSIM已经可以很好的设置各种边界条件,能够快速生成CalculiX求解器所需要的输入文件,并可以直接用于计算。用户可以使用WELSIM作为CalculiX的前处理器。CalculiX还有少量边界条件,如Constraint和Equation等, 由于常规分析中较少使用,本文未作讨论。CalculiX与Abaqus的输入文件格式有着很大的相似性,因此,本文所描述的内容,也可以用于Abaqus求解器。WelSim与作者不隶属于CalculiX, Abaqus。和CalculiX与Abaqus开发团队与机构没有直接关系。这里引用开源软件的名称和图片仅用作技术博客文章与软件使用的参考。 来源:WELSIM

未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习计划 福利任务
下载APP
联系我们
帮助与反馈